ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
QkTable.ipp
1#pragma once
2#include "CoulombIntegrals.hpp"
3#include "IO/ChronoTimer.hpp"
4#include "IO/FRW_fileReadWrite.hpp"
5#include "QkTable.hpp"
6#include <algorithm>
7#include <cassert>
8#include <cstring> // for memcpy
9#include <string_view>
10
11namespace Coulomb {
12
13//==============================================================================
14template <Symmetry S>
15void CoulombTable<S>::summary() const {
16 std::cout << "Summary: \n";
17 int k = 0;
18 auto total = 0ul;
19 for (auto &qk : m_data) {
20 std::cout << "k=" << k << ": " << qk.size() << " [" << qk.bucket_count()
21 << "]\n";
22 total += qk.size();
23 ++k;
24 }
25 std::cout << "Total: " << total << " non-zero integrals\n";
26}
27
28template <Symmetry S>
29std::size_t CoulombTable<S>::count() const {
30 auto total = 0ul;
31 for (auto &qk : m_data) {
32 total += qk.size();
33 }
34 return total;
35}
36
37//==============================================================================
38template <Symmetry S>
39void CoulombTable<S>::add(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
40 Real value) {
41 add(k, NormalOrder(a, b, c, d), value);
42}
43//------------------------------------------------------------------------------
44template <Symmetry S>
45void CoulombTable<S>::add(int k, const DiracSpinor &a, const DiracSpinor &b,
46 const DiracSpinor &c, const DiracSpinor &d,
47 Real value) {
48 add(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index(), value);
49}
50//------------------------------------------------------------------------------
51template <Symmetry S>
52void CoulombTable<S>::add(int k, nk4Index index, Real value) {
53 const auto sk = std::size_t(k);
54 if (sk >= m_data.size()) {
55 m_data.resize(sk + 1);
56 }
57 m_data.at(sk).insert({index, value});
58}
59
60//==============================================================================
61template <Symmetry S>
62// Updates Q in table. If not present, adds new Q
63void CoulombTable<S>::update(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
64 Real value) {
65 update(k, NormalOrder(a, b, c, d), value);
66}
67
68//------------------------------------------------------------------------------
69template <Symmetry S>
70// Updates Q in table. If not present, adds new Q
71void CoulombTable<S>::update(int k, const DiracSpinor &a, const DiracSpinor &b,
72 const DiracSpinor &c, const DiracSpinor &d,
73 Real value) {
74 update(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index(), value);
75}
76
77//------------------------------------------------------------------------------
78template <Symmetry S>
79void CoulombTable<S>::update(int k, nk4Index index, Real value) {
80 const auto sk = std::size_t(k);
81 if (sk >= m_data.size()) {
82 m_data.resize(sk + 1);
83 }
84 m_data.at(sk).insert_or_assign(index, value);
85}
86//==============================================================================
87template <Symmetry S>
88bool CoulombTable<S>::contains(int k, nkIndex a, nkIndex b, nkIndex c,
89 nkIndex d) const {
90 const auto sk = std::size_t(k);
91 if (sk >= m_data.size())
92 return false;
93 return m_data.at(sk).find(NormalOrder(a, b, c, d)) != m_data.at(sk).cend();
94}
95//------------------------------------------------------------------------------
96template <Symmetry S>
97bool CoulombTable<S>::contains(int k, const DiracSpinor &a,
98 const DiracSpinor &b, const DiracSpinor &c,
99 const DiracSpinor &d) const {
100 return contains(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
101}
102//------------------------------------------------------------------------------
103template <Symmetry S>
104bool CoulombTable<S>::contains(int k, nk4Index index) const {
105 const auto sk = std::size_t(k);
106 if (sk >= m_data.size())
107 return false;
108 return m_data.at(sk).find(index) != m_data.at(sk).cend();
109}
110
111//==============================================================================
112//==============================================================================
113
114template <Symmetry S>
115double *CoulombTable<S>::get(int k, nk4Index index) {
116 const auto sk = std::size_t(k);
117 if (sk >= m_data.size())
118 return nullptr;
119 // check valid k? Probably faster to lookup in table
120 const auto map_it = m_data.at(sk).find(index);
121 if (map_it == m_data.at(sk).cend())
122 return nullptr;
123 return &(map_it->second);
124}
125
126//==============================================================================
127template <Symmetry S>
128double CoulombTable<S>::Q(int k, nkIndex a, nkIndex b, nkIndex c,
129 nkIndex d) const {
130 const auto sk = std::size_t(k);
131 if (sk >= m_data.size())
132 return 0.0;
133 // check valid k? Probably faster to lookup in table
134 const auto map_it = m_data.at(sk).find(NormalOrder(a, b, c, d));
135 if (map_it == m_data.at(sk).cend())
136 return 0.0;
137 return map_it->second;
138}
139//------------------------------------------------------------------------------
140template <Symmetry S>
141double CoulombTable<S>::Q(int k, const DiracSpinor &a, const DiracSpinor &b,
142 const DiracSpinor &c, const DiracSpinor &d) const {
143 return Q(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
144}
145//------------------------------------------------------------------------------
146template <Symmetry S>
147double CoulombTable<S>::Q(int k, nk4Index index) const {
148 const auto sk = std::size_t(k);
149 if (sk >= m_data.size())
150 return 0.0;
151 // check valid k? Probably faster to lookup in table
152 const auto map_it = m_data.at(sk).find(index);
153 if (map_it == m_data.at(sk).cend())
154 return 0.0;
155 return map_it->second;
156}
157
158//==============================================================================
159template <Symmetry S>
160double CoulombTable<S>::R(int k, nkIndex a, nkIndex b, nkIndex c,
161 nkIndex d) const {
162 const auto tQk = Q(k, a, b, c, d);
163 if (tQk == 0.0)
164 return 0.0;
165 const auto s = Angular::neg1pow(k);
166
167 const auto [na, ka] = Angular::index_to_nk(a);
168 const auto [nb, kb] = Angular::index_to_nk(b);
169 const auto [nc, kc] = Angular::index_to_nk(c);
170 const auto [nd, kd] = Angular::index_to_nk(d);
171
172 const auto tCkac = Angular::tildeCk_kk(k, ka, kc);
173 const auto tCkbd = Angular::tildeCk_kk(k, kb, kd);
174 return tQk / (s * tCkac * tCkbd);
175}
176//------------------------------------------------------------------------------
177template <Symmetry S>
178double CoulombTable<S>::R(int k, const DiracSpinor &a, const DiracSpinor &b,
179 const DiracSpinor &c, const DiracSpinor &d) const {
180 return R(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
181}
182
183//==============================================================================
184template <Symmetry S>
185double CoulombTable<S>::P(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
186 const Angular::SixJTable *const sj) const {
187 double Pk_abcd{0.0};
188
189 const auto ka = Angular::nkindex_to_kappa(a);
190 const auto kb = Angular::nkindex_to_kappa(b);
191 const auto kc = Angular::nkindex_to_kappa(c);
192 const auto kd = Angular::nkindex_to_kappa(d);
193
194 const auto tja = Angular::twoj_k(ka);
195 const auto tjb = Angular::twoj_k(kb);
196 const auto tjc = Angular::twoj_k(kc);
197 const auto tjd = Angular::twoj_k(kd);
198
199 // 6j(s) Triads: {a,c,k}, {k,b,d}, {c,b,l}, {d,a,l}
200 if (Angular::triangle(tja, tjc, 2 * k) == 0 ||
201 Angular::triangle(2 * k, tjb, tjd) == 0)
202 return 0.0;
203
204 const auto [lmin, lmax] = k_minmax_Q(ka, kb, kd, kc); // exchange
205 for (int l = lmin; l <= lmax; l += 2) {
206 const auto ql = this->Q(l, a, b, d, c); // exchange
207 if (ql == 0.0)
208 continue;
209 const auto sixj = sj ? sj->get_2(tja, tjc, 2 * k, tjb, tjd, 2 * l) :
210 Angular::sixj_2(tja, tjc, 2 * k, tjb, tjd, 2 * l);
211 Pk_abcd += sixj * ql;
212 }
213 Pk_abcd *= double(2 * k + 1);
214 return Pk_abcd;
215}
216
217//------------------------------------------------------------------------------
218
219template <Symmetry S>
220double CoulombTable<S>::P(int k, const DiracSpinor &a, const DiracSpinor &b,
221 const DiracSpinor &c, const DiracSpinor &d,
222 const Angular::SixJTable *const sj) const {
223
224 return P(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index(), sj);
225}
226
227//==============================================================================
228template <Symmetry S>
229double CoulombTable<S>::P2(int k, const DiracSpinor &a, const DiracSpinor &b,
230 const DiracSpinor &c, const DiracSpinor &d,
231 const Angular::SixJTable &sj,
232 const std::vector<double> &fk) const {
233 double Pk_abcd{0.0};
234
235 // 6j(s) Triads: {a,c,k}, {k,b,d}, {c,b,l}, {d,a,l}
236 if (Coulomb::triangle(a, c, k) == 0 || Coulomb::triangle(k, b, d) == 0)
237 return 0.0;
238
239 const auto [lmin, lmax] = k_minmax_Q(a, b, d, c); // exchange
240 for (int l = lmin; l <= lmax; l += 2) {
241 const auto ql = this->Q(l, a, b, d, c); // exchange
242 if (ql == 0.0)
243 continue;
244 const auto sixj = sj.get(a, c, k, b, d, l);
245
246 // include effective Coulomb screening:
247 const auto f_scr_l = (l < (int)fk.size()) ? fk[std::size_t(l)] : 1.0;
248
249 Pk_abcd += f_scr_l * sixj * ql;
250 }
251 Pk_abcd *= double(2 * k + 1);
252 return Pk_abcd;
253}
254
255//==============================================================================
256template <Symmetry S>
257double CoulombTable<S>::W(int k, const DiracSpinor &a, const DiracSpinor &b,
258 const DiracSpinor &c, const DiracSpinor &d,
259 const Angular::SixJTable *const sj) const {
260 return W(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index(), sj);
261}
262
263template <Symmetry S>
264double CoulombTable<S>::W(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
265 const Angular::SixJTable *const sj) const {
266 return Q(k, a, b, c, d) + P(k, a, b, c, d, sj);
267}
268
269//==============================================================================
270template <Symmetry S>
271double CoulombTable<S>::g(const DiracSpinor &a, const DiracSpinor &b,
272 const DiracSpinor &c, const DiracSpinor &d, int tma,
273 int tmb, int tmc, int tmd) const {
274
275 if (tmc - tma != tmb - tmd)
276 return 0.0;
277 const int twoq = tmc - tma;
278
279 //
280 double g = 0.0;
281 const auto [k0, ki] = k_minmax_Q(a, b, c, d);
282 for (int k = k0; k <= ki; k += 2) {
283 if (std::abs(twoq) > 2 * k)
284 continue;
285
286 // const auto A =
287 // Angular::neg1pow_2(twoq) *
288 // Angular::Ck_kk_mmq(k, a.kappa(), c.kappa(), tma, tmc, -twoq) *
289 // Angular::Ck_kk_mmq(k, b.kappa(), d.kappa(), tmb, tmd, twoq);
290 // if (A != 0.0) {
291 // g += A * R(k, a, b, c, d);
292 // }
293
294 const auto tjs1 =
295 Angular::threej_2(a.twoj(), 2 * k, c.twoj(), -tma, -twoq, tmc);
296 const auto tjs2 =
297 Angular::threej_2(b.twoj(), 2 * k, d.twoj(), -tmb, twoq, tmd);
298 const auto s = Angular::neg1pow_2(2 * k + tmb - tma + twoq);
299 g += s * tjs1 * tjs2 * Q(k, a, b, c, d);
300 }
301 return g;
302}
303
304//==============================================================================
305//==============================================================================
306
307//==============================================================================
308template <>
309inline nk4Index
310CoulombTable<Symmetry::Qk>::NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
311 nkIndex d) {
312
313 // abcd -> ijkl, with i = min(a,b,c,d)
314
315 // abcd = adcb
316 // badc = bcda
317 // cbad = cdab
318 // dabc = dcba
319 // nb: there must be a more efficient way of doing this!
320
321 // const auto tmp1 = FormIndex(a, b, c, d);
322 // const auto tmp2 = FormIndex(a, d, c, b);
323 // const auto tmp3 = FormIndex(b, a, d, c);
324 // const auto tmp4 = FormIndex(b, c, d, a);
325 // const auto tmp5 = FormIndex(c, b, a, d);
326 // const auto tmp6 = FormIndex(c, d, a, b);
327 // const auto tmp7 = FormIndex(d, a, b, c);
328 // const auto tmp8 = FormIndex(d, c, b, a);
329 // return std::min({tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8});
330
331 nk4Index out = FormIndex(a, b, c, d);
332 const auto min = std::min({a, b, c, d});
333 if (a == min) {
334 const auto tmp = FormIndex(a, d, c, b);
335 if (tmp < out)
336 out = tmp;
337 }
338 if (b == min) {
339 const auto tmp = std::min(FormIndex(b, a, d, c), FormIndex(b, c, d, a));
340 if (tmp < out)
341 out = tmp;
342 }
343 if (c == min) {
344 const auto tmp = std::min(FormIndex(c, b, a, d), FormIndex(c, d, a, b));
345 if (tmp < out)
346 out = tmp;
347 }
348 if (d == min) {
349 const auto tmp = std::min(FormIndex(d, a, b, c), FormIndex(d, c, b, a));
350 if (tmp < out)
351 out = tmp;
352 }
353 return out;
354}
355//------------------------------------------------------------------------------
356template <>
357inline nk4Index
358CoulombTable<Symmetry::Wk>::NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
359 nkIndex d) {
360 const auto tmp1 = FormIndex(a, b, c, d);
361 const auto tmp2 = FormIndex(b, a, d, c);
362 const auto tmp3 = FormIndex(c, d, a, b);
363 const auto tmp4 = FormIndex(d, c, b, a);
364 return std::min({tmp1, tmp2, tmp3, tmp4});
365}
366
367//------------------------------------------------------------------------------
368template <>
369inline nk4Index
370CoulombTable<Symmetry::Lk>::NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
371 nkIndex d) {
372 const auto tmp1 = FormIndex(a, b, c, d);
373 const auto tmp2 = FormIndex(b, a, d, c);
374 return std::min({tmp1, tmp2});
375}
376
377//------------------------------------------------------------------------------
378template <>
379inline nk4Index
380CoulombTable<Symmetry::none>::NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
381 nkIndex d) {
382 return FormIndex(a, b, c, d);
383}
384
385//==============================================================================
386template <Symmetry S>
387nk4Index CoulombTable<S>::NormalOrder(nkIndex a, nkIndex b, nkIndex c,
388 nkIndex d) const {
389 return NormalOrder_impl(a, b, c, d);
390}
391//------------------------------------------------------------------------------
392template <Symmetry S>
393nk4Index
394CoulombTable<S>::NormalOrder(const DiracSpinor &a, const DiracSpinor &b,
395 const DiracSpinor &c, const DiracSpinor &d) const {
396 return NormalOrder(a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
397}
398
399//==============================================================================
400template <Symmetry S>
401nk4Index CoulombTable<S>::CurrentOrder(nkIndex a, nkIndex b, nkIndex c,
402 nkIndex d) const {
403 return FormIndex(a, b, c, d);
404}
405//------------------------------------------------------------------------------
406template <Symmetry S>
407nk4Index CoulombTable<S>::CurrentOrder(const DiracSpinor &a,
408 const DiracSpinor &b,
409 const DiracSpinor &c,
410 const DiracSpinor &d) const {
411 return CurrentOrder(a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
412}
413
414//==============================================================================
415template <Symmetry S>
416nk4Index CoulombTable<S>::FormIndex(nkIndex a, nkIndex b, nkIndex c,
417 nkIndex d) {
418 // create a nk4Index from four small nkIndex, by laying out in memory
419
420 // static_assert(sizeof(nk4Index) == 4 * sizeof(nkIndex));
421 // nk4Index big_index;
422 // const auto p_index = reinterpret_cast<nkIndex *>(&big_index);
423 // std::memcpy(p_index, &d, sizeof(nkIndex));
424 // std::memcpy(p_index + 1, &c, sizeof(nkIndex));
425 // std::memcpy(p_index + 2, &b, sizeof(nkIndex));
426 // std::memcpy(p_index + 3, &a, sizeof(nkIndex));
427 // return big_index;
428
429 // this seems slightly faster...though variance large
430 static_assert(sizeof(nk4Index) == 4 * sizeof(nkIndex));
431 static_assert(sizeof(nkIndex) * 8 == 16);
432 return ((nk4Index)a << 48) + ((nk4Index)b << 32) + ((nk4Index)c << 16) +
433 (nk4Index)d;
434}
435
436//==============================================================================
437template <Symmetry S>
438std::array<nkIndex, 4>
439CoulombTable<S>::UnFormIndex(const nk4Index &index) const {
440 static_assert(sizeof(nk4Index) == sizeof(std::array<nkIndex, 4>));
441 std::array<nkIndex, 4> set;
442 std::memcpy(&set, &index, sizeof(nk4Index));
443 // note: this returns {d, c, b, a}, not {a, b, c, d}
444 // A warning: this might depend on "endiness"
445 // It's OK, this function is never used???
446 // If there's an issue, a unit test will fail, so it's "safe"
447 std::swap(set[0], set[3]);
448 std::swap(set[1], set[2]);
449 return set;
450}
451
452//==============================================================================
453//==============================================================================
454template <Symmetry S>
455void CoulombTable<S>::fill(const std::vector<DiracSpinor> &basis,
456 const YkTable &yk, int k_cut, bool print) {
457 static_assert(S == Symmetry::Qk);
458 IO::ChronoTimer t("fill");
459
460 /*
461 In order to make best use of CPU and Memory, we "fill" the QK table in a
462 strange 4-step manner.
463 1. Count the number of non-zero Q's for each k
464 2. Use info from aboe to 'reserve()' space in the map
465 3. Fill all the non-zero Q entries in the map with 0. Note that this adds new
466 elements into the map, so cannot be done in parallel. Since we store a
467 different map for each k, we may parallelise over k (though this is not super
468 efficient)
469 4. Now that we have a fully sized map, we can update each of its values in
470 parallel
471 Note that we make use of the symmetry to ensure we do not access any element
472 from more than 1 thread (also, allowing us to not do any more calculations
473 than necisary)
474 */
475
476 const auto tmp_max_k =
477 std::size_t(std::max(DiracSpinor::max_tj(basis), 1) - 1);
478
479 const auto max_k =
480 (k_cut <= 0) ? tmp_max_k : std::min(tmp_max_k, std::size_t(k_cut));
481
482 if (m_data.size() < max_k + 1)
483 m_data.resize(max_k + 1);
484
485 // 1) Count non-zero Q integrals (each k). Use this to 'reserve' map space
486 t.start();
487 std::vector<std::size_t> count_non_zero_k(max_k + 1);
488#pragma omp parallel for
489 for (auto k = 0ul; k <= max_k; ++k) {
490 const auto ik = static_cast<int>(k);
491 for (const auto &a : basis) {
492 for (const auto &b : basis) {
493 for (const auto &c : basis) {
494 if (!Angular::Ck_kk_SR(ik, a.kappa(), c.kappa()))
495 continue;
496 for (const auto &d : basis) {
497 // due to symmetry, only calculate each 'unique' integral once
498 if (NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d)) {
499 if (Angular::Ck_kk_SR(ik, b.kappa(), d.kappa())) {
500 ++count_non_zero_k[k];
501 }
502 }
503 }
504 }
505 }
506 }
507 }
508 if (print)
509 std::cout << "Count non-zero: " << t.lap_reading_str() << std::endl;
510
511 // 2) Reserve space in each sub-map
512 t.start();
513#pragma omp parallel for
514 for (auto ik = 0ul; ik <= max_k; ++ik) {
515 m_data[ik].reserve(count_non_zero_k[ik]);
516 }
517 if (print)
518 std::cout << "Reserve: " << t.lap_reading_str() << std::endl;
519
520 // 3) Create space in map (set each element to zero).
521 t.start();
522#pragma omp parallel for
523 for (auto k = 0ul; k <= max_k; ++k) {
524 for (const auto &a : basis) {
525 for (const auto &c : basis) {
526 if (!Angular::Ck_kk_SR(int(k), a.kappa(), c.kappa()))
527 continue;
528 for (const auto &b : basis) {
529 for (const auto &d : basis) {
530 if (!Angular::Ck_kk_SR(int(k), b.kappa(), d.kappa()))
531 continue;
532 const auto normal_index = NormalOrder(a, b, c, d);
533 if (normal_index == CurrentOrder(a, b, c, d)) {
534 add(int(k), normal_index, 0.0);
535 }
536 }
537 }
538 }
539 }
540 }
541 if (print)
542 std::cout << "Fill w/ zeros: " << t.lap_reading_str() << std::endl;
543
544 // 4) Fill the pre-constructed map with values, in parallel. Since we are
545 // not adding any new elements to map, and since we are guarenteed to only
546 // access each map element once, we can do this part in parallel. nb: This
547 // //isation is not very efficient, though in theory it can be 100%
548 if (print)
549 std::cout << "Fill w/ values: " << std::flush;
550 t.start();
551#pragma omp parallel for collapse(2)
552 for (auto ia = 0ul; ia < basis.size(); ++ia) {
553 for (auto ib = 0ul; ib < basis.size(); ++ib) {
554 const auto &a = basis[ia];
555 const auto &b = basis[ib];
556 for (const auto &c : basis) {
557 for (const auto &d : basis) {
558 const auto normal_index = NormalOrder(a, b, c, d);
559 if (normal_index == CurrentOrder(a, b, c, d)) {
560 const auto [kmin, kmax] = k_minmax_Q(a, b, c, d);
561 for (int k = kmin; k <= kmax && k <= int(max_k); k += 2) {
562 double *ptr = get(k, normal_index);
563 assert(ptr != nullptr);
564 if (*ptr == 0.0) {
565 // only calculate if not already in table
566 // This saves some time, but surprisingly little!
567 *ptr = yk.Q(k, a, b, c, d);
568 }
569 // update(k, a, b, c, d, yk.Q(k, a, b, c, d));
570 }
571 }
572 }
573 }
574 }
575 }
576
577 if (print)
578 std::cout << t.lap_reading_str() << std::endl;
579
580 if (print)
581 summary();
582}
583
584//==============================================================================
585template <Symmetry S>
586void CoulombTable<S>::fill_if(const std::vector<DiracSpinor> &basis,
587 const YkTable &yk,
588 const SelectionRules &SelectionFunction,
589 int k_cut, bool print) {
590 static_assert(S == Symmetry::Qk);
591 IO::ChronoTimer t("fill");
592
593 /*
594 Same as above, but only calculate integrals satisfying the SelectionFunction
595 */
596
597 const auto tmp_max_k =
598 std::size_t(std::max(DiracSpinor::max_tj(basis), 1) - 1);
599
600 const auto max_k =
601 (k_cut <= 0) ? tmp_max_k : std::min(tmp_max_k, std::size_t(k_cut));
602
603 if (m_data.size() < max_k + 1)
604 m_data.resize(max_k + 1);
605
606 // 1) Count non-zero Q integrals (each k). Use this to 'reserve' map space
607 t.start();
608 std::vector<std::size_t> count_non_zero_k(max_k + 1);
609#pragma omp parallel for
610 for (auto k = 0ul; k <= max_k; ++k) {
611 const auto ik = static_cast<int>(k);
612 for (const auto &a : basis) {
613 for (const auto &b : basis) {
614 for (const auto &c : basis) {
615 if (!Angular::Ck_kk_SR(ik, a.kappa(), c.kappa()))
616 continue;
617 for (const auto &d : basis) {
618 // due to symmetry, only calculate each 'unique' integral once
619 if (NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d)) {
620 if (Angular::Ck_kk_SR(ik, b.kappa(), d.kappa()) &&
621 SelectionFunction(ik, a, b, c, d)) {
622 ++count_non_zero_k[k];
623 }
624 }
625 }
626 }
627 }
628 }
629 }
630 if (print)
631 std::cout << "Count non-zero: " << t.lap_reading_str() << std::endl;
632
633 // 2) Reserve space in each sub-map
634 t.start();
635#pragma omp parallel for
636 for (auto ik = 0ul; ik <= max_k; ++ik) {
637 m_data[ik].reserve(count_non_zero_k[ik]);
638 }
639 if (print)
640 std::cout << "Reserve: " << t.lap_reading_str() << std::endl;
641
642 // 3) Create space in map (set each element to zero).
643 t.start();
644#pragma omp parallel for
645 for (auto k = 0ul; k <= max_k; ++k) {
646 for (const auto &a : basis) {
647 for (const auto &c : basis) {
648 if (!Angular::Ck_kk_SR(int(k), a.kappa(), c.kappa()))
649 continue;
650 for (const auto &b : basis) {
651 for (const auto &d : basis) {
652 if (!Angular::Ck_kk_SR(int(k), b.kappa(), d.kappa()))
653 continue;
654 if (!SelectionFunction(int(k), a, b, c, d))
655 continue;
656 const auto normal_index = NormalOrder(a, b, c, d);
657 if (normal_index == CurrentOrder(a, b, c, d)) {
658 add(int(k), normal_index, 0.0);
659 }
660 }
661 }
662 }
663 }
664 }
665 if (print)
666 std::cout << "Fill w/ zeros: " << t.lap_reading_str() << std::endl;
667
668 // 4) Fill the pre-constructed map with values, in parallel. Since we are
669 // not adding any new elements to map, and since we are guarenteed to only
670 // access each map element once, we can do this part in parallel. nb: This
671 // //isation is not very efficient, though in theory it can be 100%
672 if (print)
673 std::cout << "Fill w/ values: " << std::flush;
674 t.start();
675#pragma omp parallel for collapse(2)
676 for (auto ia = 0ul; ia < basis.size(); ++ia) {
677 for (auto ib = 0ul; ib < basis.size(); ++ib) {
678 const auto &a = basis[ia];
679 const auto &b = basis[ib];
680 for (const auto &c : basis) {
681 for (const auto &d : basis) {
682 const auto normal_index = NormalOrder(a, b, c, d);
683 if (normal_index == CurrentOrder(a, b, c, d)) {
684 const auto [kmin, kmax] = k_minmax_Q(a, b, c, d);
685 for (int k = kmin; k <= kmax && k <= int(max_k); k += 2) {
686 if (!SelectionFunction(k, a, b, c, d))
687 continue;
688 double *ptr = get(k, normal_index);
689 assert(ptr != nullptr);
690 if (*ptr == 0.0) {
691 // only calculate if not already in table
692 // This saves some time, but surprisingly little!
693 *ptr = yk.Q(k, a, b, c, d);
694 }
695 }
696 }
697 }
698 }
699 }
700 }
701
702 if (print)
703 std::cout << t.lap_reading_str() << std::endl;
704
705 if (print)
706 summary();
707}
708
709//==============================================================================
710template <Symmetry S>
711void CoulombTable<S>::fill(const std::vector<DiracSpinor> &basis,
712 const CoulombFunction &Fk,
713 const SelectionRules &Fk_SR, int k_cut, bool print) {
714 IO::ChronoTimer t("fill");
715
716 /*
717 In order to make best use of CPU and Memory, we "fill" the QK table in a
718 strange 4-step manner.
719 1. Count the number of non-zero Q's for each k
720 2. Use info from aboe to 'reserve()' space in the map
721 3. Fill all the non-zero Q entries in the map with 0. Note that this adds new
722 elements into the map, so cannot be done in parallel. Since we store a
723 different map for each k, we may parallelise over k (though this is not super
724 efficient)
725 4. Now that we have a fully sized map, we can update each of its values in
726 parallel
727 Note that we make use of the symmetry to ensure we do not access any element
728 from more than 1 thread (also, allowing us to not do any more calculations
729 than necisary)
730 */
731
732 const auto tmp_max_k = std::size_t(DiracSpinor::max_tj(basis));
733 // May have different parity rule, so don't -1 here
734
735 const auto max_k =
736 (k_cut <= 0) ? tmp_max_k : std::min(tmp_max_k, std::size_t(k_cut));
737
738 if (m_data.size() < max_k + 1)
739 m_data.resize(max_k + 1);
740
741 // 1) Count non-zero Q integrals (each k). Use this to 'reserve' map space
742 t.start();
743 std::vector<std::size_t> count_non_zero_k(max_k + 1);
744#pragma omp parallel for
745 for (auto k = 0ul; k <= max_k; ++k) {
746 const auto ik = static_cast<int>(k);
747 for (const auto &a : basis) {
748 for (const auto &b : basis) {
749 for (const auto &c : basis) {
750 for (const auto &d : basis) {
751 // due to symmetry, only calculate each 'unique' integral once
752 if (NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d)) {
753 if (Fk_SR(ik, a, b, c, d)) {
754 ++count_non_zero_k[k];
755 }
756 }
757 }
758 }
759 }
760 }
761 }
762 if (print)
763 std::cout << "Count non-zero: " << t.lap_reading_str() << std::endl;
764
765 // 2) Reserve space in each sub-map
766 t.start();
767#pragma omp parallel for
768 for (auto ik = 0ul; ik <= max_k; ++ik) {
769 m_data[ik].reserve(count_non_zero_k[ik]);
770 }
771 if (print)
772 std::cout << "Reserve: " << t.lap_reading_str() << std::endl;
773
774 // 3) Create space in map (set each element to zero).
775 t.start();
776#pragma omp parallel for
777 for (auto k = 0ul; k <= max_k; ++k) {
778 for (const auto &a : basis) {
779 for (const auto &b : basis) {
780 for (const auto &c : basis) {
781 for (const auto &d : basis) {
782 if (Fk_SR(int(k), a, b, c, d)) {
783 const auto normal_index = NormalOrder(a, b, c, d);
784 if (normal_index == CurrentOrder(a, b, c, d)) {
785 add(int(k), normal_index, 0.0);
786 }
787 }
788 }
789 }
790 }
791 }
792 }
793 if (print)
794 std::cout << "Fill w/ zeros: " << t.lap_reading_str() << std::endl;
795
796 // 4) Fill the pre-constructed map with values, in parallel. Since we are
797 // not adding any new elements to map, and since we are guarenteed to only
798 // access each map element once, we can do this part in parallel. nb: This
799 // //isation is not very efficient, though in theory it can be 100%
800 if (print)
801 std::cout << "Fill w/ values: " << std::flush;
802 t.start();
803
804#pragma omp parallel for collapse(2)
805 for (std::size_t ia = 0; ia < basis.size(); ++ia) {
806 for (std::size_t ib = 0; ib < basis.size(); ++ib) {
807 const auto &a = basis[ia];
808 const auto &b = basis[ib];
809 for (const auto &c : basis) {
810 for (const auto &d : basis) {
811 const auto normal_index = NormalOrder(a, b, c, d);
812 if (normal_index == CurrentOrder(a, b, c, d)) {
813 for (int k = 0; k <= int(max_k); ++k) {
814 if (Fk_SR(k, a, b, c, d)) {
815 double *ptr = get(k, normal_index);
816 assert(ptr != nullptr);
817 if (*ptr == 0.0) {
818 // only calculate if not already in table
819 *ptr = Fk(k, a, b, c, d);
820 }
821 // update(k, a, b, c, d, Fk(k, a, b, c, d));
822 }
823 }
824 }
825 }
826 }
827 }
828 }
829 if (print)
830 std::cout << t.lap_reading_str() << std::endl;
831
832 if (print)
833 summary();
834}
835
836//==============================================================================
837template <Symmetry S>
838void CoulombTable<S>::write(const std::string &fname) const {
839 if (fname == "false")
840 return;
841 std::cout << "Writing " << count() << " integrals to file: " << fname << ".."
842 << std::flush;
843 std::fstream f;
844 const auto rw = IO::FRW::write;
845 IO::FRW::open_binary(f, fname, rw);
846
847 auto size = m_data.size();
848 rw_binary(f, rw, size);
849 for (const auto &Q_k : m_data) {
850 auto size_k = Q_k.size();
851 rw_binary(f, rw, size_k);
852 for (auto [key, value] : Q_k) {
853 auto key_copy = key; // have no pass non-const reference!
854 rw_binary(f, rw, key_copy, value);
855 }
856 }
857 std::cout << "\n" << std::flush;
858}
859
860//==============================================================================
861template <Symmetry S>
862bool CoulombTable<S>::read(const std::string &fname) {
863 if (fname == "false")
864 return false;
865 std::fstream f;
866 const auto rw = IO::FRW::read;
867 IO::FRW::open_binary(f, fname, rw);
868
869 if (!f.good())
870 return false;
871
872 std::size_t size{0};
873 rw_binary(f, rw, size);
874 m_data.resize(size);
875 if (m_data.size() == 0)
876 return false;
877 for (std::size_t i = 0; i < m_data.size(); ++i) {
878 auto &Q_k = m_data[i];
879 std::size_t size_k{0};
880 rw_binary(f, rw, size_k);
881 Q_k.reserve(size_k);
882 for (std::size_t ik = 0; ik < size_k; ++ik) {
883 nk4Index key;
884 Real value;
885 rw_binary(f, rw, key, value);
886 Q_k[key] = value;
887 }
888 }
889 std::cout << "Read " << count() << " integrals from file: " << fname << "\n"
890 << std::flush;
891 return true;
892}
893
894} // namespace Coulomb
Angular::SixJTable
Lookup table for Wigner 6J symbols.
Definition SixJTable.hpp:24
Angular::SixJTable::get_2
double get_2(int a, int b, int c, int d, int e, int f) const
Return 6j symbol {a/2,b/2,c/2,d/2,e/2,f/2}. Note: takes in 2*j/2*k as int.
Definition SixJTable.hpp:53
Angular::SixJTable::get
double get(const A &a, const B &b, const C &c, const D &d, const E &e, const F &f) const
Return "magic" 6j. Pass in either integer (for k), or DiracSpinor, F. e.g.: (F,F,k,...
Definition SixJTable.hpp:65
Coulomb::CoulombTable
Base (pure virtual) class to store Coulomb integrals, and similar. 3 derived classes,...
Definition QkTable.hpp:57
Coulomb::CoulombTable::P2
Real P2(int k, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d, const Angular::SixJTable &sj, const std::vector< double > &fk) const
'Exchange-only', defined via W = Q + P. Optionally, takes pointer to 6J table (faster eval of 6J symb...
Definition QkTable.ipp:229
Coulomb::CoulombTable::get
double * get(int k, nk4Index index)
Directly gets one of the stored elements, given normal-ordered nk4Index.
Definition QkTable.ipp:115
Coulomb::CoulombTable::fill_if
void fill_if(const std::vector< DiracSpinor > &basis, const YkTable &yk, const SelectionRules &SelectionFunction, int k_cut=-1, bool print=true)
Takes a constructed YkTable, and fills Coulomb table with all possible non-zero Qk Coulomb integrals ...
Definition QkTable.ipp:586
Coulomb::CoulombTable::P
Real P(int k, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d, const Angular::SixJTable *const sj=nullptr) const
'Exchange-only', defined via W = Q + P. Optionally, takes pointer to 6J table (faster eval of 6J symb...
Definition QkTable.ipp:220
Coulomb::CoulombTable::read
bool read(const std::string &fname)
Reads coulomb integrals to disk. Returns false if none read in.
Definition QkTable.ipp:862
Coulomb::CoulombTable::update
void update(int k, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d, Real value)
Updates value in table. If not present, adds new value.
Definition QkTable.ipp:71
Coulomb::CoulombTable::contains
bool contains(int k, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d) const
Checks if given {k,a,b,c,d} is in the table.
Definition QkTable.ipp:97
Coulomb::CoulombTable::NormalOrder
nk4Index NormalOrder(const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d) const
Creates single 'nk4Index' corresponding to 'NormalOrder' symmetry of {a,b,c,d}.
Definition QkTable.ipp:394
Coulomb::CoulombTable::CurrentOrder
nk4Index CurrentOrder(const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d) const
Creates single 'nk4Index', WITHOUT accounting for 'NormalOrder'. Can be used to check if {a,...
Definition QkTable.ipp:407
Coulomb::CoulombTable::FormIndex
static nk4Index FormIndex(nkIndex a, nkIndex b, nkIndex c, nkIndex d)
Converts given set of nkIndex's (in any order) to nk4Index.
Definition QkTable.ipp:416
Coulomb::CoulombTable::R
Real R(int k, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d) const
Returns 'R', defined via: R := Q / (angular_coef)
Definition QkTable.ipp:178
Coulomb::CoulombTable::UnFormIndex
std::array< nkIndex, 4 > UnFormIndex(const nk4Index &index) const
Breaks nk4Index back into {ia,ib,ic,id}. Not used.
Definition QkTable.ipp:439
Coulomb::CoulombTable::fill
void fill(const std::vector< DiracSpinor > &basis, const YkTable &yk, int k_cut=-1, bool print=true)
Takes a constructed YkTable, and fills Coulomb table with all possible non-zero Qk Coulomb integrals,...
Definition QkTable.ipp:455
Coulomb::CoulombTable::write
void write(const std::string &fname) const
Writes coulomb integrals to disk.
Definition QkTable.ipp:838
Coulomb::CoulombTable::summary
void summary() const
For testing: prints details of coulomb integrals stored.
Definition QkTable.ipp:15
Coulomb::CoulombTable::g
Real g(const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d, int tma, int tmb, int tmc, int tmd) const
Returns 'g_abcd'.
Definition QkTable.ipp:271
Coulomb::CoulombTable::add
void add(int k, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d, Real value)
adds a new value. Note: does nothing if {k,a,b,c,d} already exists
Definition QkTable.ipp:45
Coulomb::CoulombTable::Q
Real Q(int k, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d) const
Retrieve a stored Q. If not present, returns 0. (Returns exactly as stored in table....
Definition QkTable.ipp:141
Coulomb::CoulombTable::W
Real W(int k, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d, const Angular::SixJTable *const sj=nullptr) const
W^k_abcd = Q^k_abcd + \sum_l [k] 6j * Q^l_abdc. Optionally, takes pointer to 6J table (faster eval of...
Definition QkTable.ipp:257
Coulomb::CoulombTable::count
std::size_t count() const
Returns number of stored integrals.
Definition QkTable.ipp:29
Coulomb::YkTable
Calculates + stores Hartree Y functions + Angular (w/ look-up), taking advantage of symmetry.
Definition YkTable.hpp:26
Coulomb::YkTable::Q
double Q(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb, const DiracSpinor &Fc, const DiracSpinor &Fd) const
Calculates Qk using the existing yk integrals. Note: Yk and Ck tables must include all required value...
Definition YkTable.cpp:127
DiracSpinor
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:41
DiracSpinor::max_tj
static int max_tj(const std::vector< DiracSpinor > &orbs)
Returns maximum (2j) found in {orbs}.
Definition DiracSpinor.cpp:266
DiracSpinor::nk_index
Index nk_index() const
(n,kappa) index (see AtomData)
Definition DiracSpinor.hpp:97
Angular::nkindex_to_kappa
int nkindex_to_kappa(int index)
Returns kappa, given nk_index.
Definition Wigner369j.hpp:87
Angular::threej_2
double threej_2(int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3)
Calculates wigner 3j symbol. Takes in 2*j (or 2*l) - intput is integer.
Definition Wigner369j.hpp:146
Angular::sixj_2
double sixj_2(int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6)
6j symbol {j1 j2 j3 \ j4 j5 j6} - [takes 2*j as int]
Definition Wigner369j.hpp:236
Angular::neg1pow
constexpr int neg1pow(int a)
Evaluates (-1)^{a} (for integer a)
Definition Wigner369j.hpp:119
Angular::tildeCk_kk
double tildeCk_kk(int k, int ka, int kb)
tildeCk_kk = (-1)^{ja+1/2}*Ck_kk
Definition Wigner369j.hpp:302
Angular::index_to_nk
std::pair< int, int > index_to_nk(int index)
return {n, kappa} given nk_index:
Definition Wigner369j.hpp:78
Angular::Ck_kk_SR
bool Ck_kk_SR(int k, int ka, int kb)
Ck selection rule only. Returns false if C^k=0, true if non-zero.
Definition Wigner369j.hpp:293
Angular::neg1pow_2
constexpr int neg1pow_2(int two_a)
Evaluates (-1)^{two_a/2} (for integer a; two_a is even)
Definition Wigner369j.hpp:121
Angular::triangle
constexpr int triangle(int j1, int j2, int J)
Returns 1 if triangle rule is satisfied. nb: works with j OR twoj!
Definition Wigner369j.hpp:131
Angular::twoj_k
constexpr int twoj_k(int ka)
returns 2j given kappa
Definition Wigner369j.hpp:14
Coulomb
Functions (+classes) for computing Coulomb integrals.
Definition Coulomb.hpp:8
Coulomb::Real
double Real
Data type used to store integrals.
Definition QkTable.hpp:15
Coulomb::SelectionRules
std::function< bool(int, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d)> SelectionRules
Function type for determining Coulomb(-like) integral selection rules. Takes k and 4 DiracSpinors,...
Definition QkTable.hpp:33
Coulomb::nk4Index
uint64_t nk4Index
index type for set of 4 orbitals {nk,nk,nk,nk} -> nk4Index
Definition QkTable.hpp:17
Coulomb::Pk_abcd
double Pk_abcd(const int k, const DiracSpinor &Fa, const DiracSpinor &Fb, const DiracSpinor &Fc, const DiracSpinor &Fd)
Exchange only version of W (W-Q): W = Q + P [see Qk above].
Definition CoulombIntegrals.cpp:428
Coulomb::nkIndex
uint16_t nkIndex
index type for each {nk} (orbital)
Definition QkTable.hpp:19
Coulomb::k_minmax_Q
std::pair< int, int > k_minmax_Q(const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d)
Returns min and max k (multipolarity) allowed for Q^k_abcd. Parity rule is included,...
Definition CoulombIntegrals.cpp:556
Coulomb::CoulombFunction
std::function< double(int, const DiracSpinor &a, const DiracSpinor &b, const DiracSpinor &c, const DiracSpinor &d)> CoulombFunction
Function type for calculating Coulomb(-like) integrals. Takes k and 4 DiracSpinors,...
Definition QkTable.hpp:27
qip::add
void add(std::vector< T > *first, const std::vector< T > &second, const Args &...rest)
Adds any number of vectors, in place (modifies first vector). Must be of same type....
Definition Vector.hpp:106
Generated on Wed Jan 29 2025 05:05:22 for ampsci by doxygen 1.9.8