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