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 "fmt/format.hpp"
8#include "qip/Widgets.hpp"
9#include <algorithm>
10#include <cassert>
11#include <cstring> // for memcpy
12#include <numeric>
13#include <string_view>
14
15namespace Coulomb {
16
17//==============================================================================
18template <Symmetry S>
19void CoulombTable<S>::summary() const {
20 // const auto symm = to_string(S);
21 std::cout << "Stored non-zero two-body integrals (" << to_string(S)
22 << " symmetry): \n";
23 int k = 0;
24 auto total = 0ul;
25 for (auto &qk : m_data) {
26 fmt::print("k = {:2}: {}\n", k, qk.size());
27 // << " [" << qk.bucket_count()<< "]\n";
28 total += qk.size();
29 ++k;
30 }
31 std::cout << "Total: " << total << "\n";
32}
33
34template <Symmetry S>
35std::size_t CoulombTable<S>::count() const {
36 auto total = 0ul;
37 for (auto &qk : m_data) {
38 total += qk.size();
39 }
40 return total;
41}
42
43//==============================================================================
44template <Symmetry S>
45void CoulombTable<S>::add(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
46 Real value) {
47 add(k, NormalOrder(a, b, c, d), value);
48}
49//------------------------------------------------------------------------------
50template <Symmetry S>
51void CoulombTable<S>::add(int k, const DiracSpinor &a, const DiracSpinor &b,
52 const DiracSpinor &c, const DiracSpinor &d,
53 Real value) {
54 add(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index(), value);
55}
56//------------------------------------------------------------------------------
57template <Symmetry S>
58void CoulombTable<S>::add(int k, nk4Index index, Real value) {
59 const auto sk = std::size_t(k);
60 if (sk >= m_data.size()) {
61 m_data.resize(sk + 1);
62 }
63 m_data.at(sk).insert({index, value});
64}
65
66//==============================================================================
67template <Symmetry S>
68// Updates Q in table. If not present, adds new Q
69void CoulombTable<S>::update(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
70 Real value) {
71 update(k, NormalOrder(a, b, c, d), value);
72}
73
74//------------------------------------------------------------------------------
75template <Symmetry S>
76// Updates Q in table. If not present, adds new Q
77void CoulombTable<S>::update(int k, const DiracSpinor &a, const DiracSpinor &b,
78 const DiracSpinor &c, const DiracSpinor &d,
79 Real value) {
80 update(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index(), value);
81}
82
83//------------------------------------------------------------------------------
84template <Symmetry S>
85void CoulombTable<S>::update(int k, nk4Index index, Real value) {
86 const auto sk = std::size_t(k);
87 if (sk >= m_data.size()) {
88 m_data.resize(sk + 1);
89 }
90 m_data.at(sk).insert_or_assign(index, value);
91}
92//==============================================================================
93template <Symmetry S>
94bool CoulombTable<S>::contains(int k, nkIndex a, nkIndex b, nkIndex c,
95 nkIndex d) const {
96 const auto sk = std::size_t(k);
97 if (sk >= m_data.size())
98 return false;
99 return m_data.at(sk).find(NormalOrder(a, b, c, d)) != m_data.at(sk).cend();
100}
101//------------------------------------------------------------------------------
102template <Symmetry S>
103bool CoulombTable<S>::contains(int k, const DiracSpinor &a,
104 const DiracSpinor &b, const DiracSpinor &c,
105 const DiracSpinor &d) const {
106 return contains(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
107}
108//------------------------------------------------------------------------------
109template <Symmetry S>
110bool CoulombTable<S>::contains(int k, nk4Index index) const {
111 const auto sk = std::size_t(k);
112 if (sk >= m_data.size())
113 return false;
114 return m_data.at(sk).find(index) != m_data.at(sk).cend();
115}
116
117//==============================================================================
118//==============================================================================
119
120template <Symmetry S>
121double *CoulombTable<S>::get(int k, nk4Index index) {
122 const auto sk = std::size_t(k);
123 if (sk >= m_data.size())
124 return nullptr;
125 // check valid k? Probably faster to lookup in table
126 const auto map_it = m_data.at(sk).find(index);
127 if (map_it == m_data.at(sk).cend())
128 return nullptr;
129 return &(map_it->second);
130}
131
132//==============================================================================
133template <Symmetry S>
134double CoulombTable<S>::Q(int k, nkIndex a, nkIndex b, nkIndex c,
135 nkIndex d) const {
136 const auto sk = std::size_t(k);
137 if (sk >= m_data.size())
138 return 0.0;
139 // check valid k? Probably faster to lookup in table
140 const auto map_it = m_data.at(sk).find(NormalOrder(a, b, c, d));
141 if (map_it == m_data.at(sk).cend())
142 return 0.0;
143 return map_it->second;
144}
145//------------------------------------------------------------------------------
146template <Symmetry S>
147double CoulombTable<S>::Q(int k, const DiracSpinor &a, const DiracSpinor &b,
148 const DiracSpinor &c, const DiracSpinor &d) const {
149 return Q(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
150}
151//------------------------------------------------------------------------------
152template <Symmetry S>
153double CoulombTable<S>::Q(int k, nk4Index index) const {
154 const auto sk = std::size_t(k);
155 if (sk >= m_data.size())
156 return 0.0;
157 // check valid k? Probably faster to lookup in table
158 const auto map_it = m_data.at(sk).find(index);
159 if (map_it == m_data.at(sk).cend())
160 return 0.0;
161 return map_it->second;
162}
163
164//==============================================================================
165template <Symmetry S>
166double CoulombTable<S>::R(int k, nkIndex a, nkIndex b, nkIndex c,
167 nkIndex d) const {
168 const auto tQk = Q(k, a, b, c, d);
169 if (tQk == 0.0)
170 return 0.0;
171 const auto s = Angular::neg1pow(k);
172
173 const auto [na, ka] = Angular::index_to_nk(a);
174 const auto [nb, kb] = Angular::index_to_nk(b);
175 const auto [nc, kc] = Angular::index_to_nk(c);
176 const auto [nd, kd] = Angular::index_to_nk(d);
177
178 const auto tCkac = Angular::tildeCk_kk(k, ka, kc);
179 const auto tCkbd = Angular::tildeCk_kk(k, kb, kd);
180 return tQk / (s * tCkac * tCkbd);
181}
182//------------------------------------------------------------------------------
183template <Symmetry S>
184double CoulombTable<S>::R(int k, const DiracSpinor &a, const DiracSpinor &b,
185 const DiracSpinor &c, const DiracSpinor &d) const {
186 return R(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
187}
188
189//==============================================================================
190template <Symmetry S>
191double CoulombTable<S>::P(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
192 const Angular::SixJTable *const sj) const {
193 double Pk_abcd{0.0};
194
195 const auto ka = Angular::nkindex_to_kappa(a);
196 const auto kb = Angular::nkindex_to_kappa(b);
197 const auto kc = Angular::nkindex_to_kappa(c);
198 const auto kd = Angular::nkindex_to_kappa(d);
199
200 const auto tja = Angular::twoj_k(ka);
201 const auto tjb = Angular::twoj_k(kb);
202 const auto tjc = Angular::twoj_k(kc);
203 const auto tjd = Angular::twoj_k(kd);
204
205 // 6j(s) Triads: {a,c,k}, {k,b,d}, {c,b,l}, {d,a,l}
206 if (Angular::triangle(tja, tjc, 2 * k) == 0 ||
207 Angular::triangle(2 * k, tjb, tjd) == 0)
208 return 0.0;
209
210 // We may store non-Coulomb integrals in this table (for example, Breit)
211 // Then, the parity selection rule doesn't apply,
212 // and we don't know what the actual rule is, so use whole range
213 const auto [lmin, lmax] = [&]() -> std::pair<int, int> {
214 if constexpr (S == Symmetry::Qk)
215 return k_minmax_Q(ka, kb, kd, kc);
216 else
217 return {0, max_k()};
218 }();
219 const int dl = (S == Symmetry::Qk) ? 2 : 1;
220
221 for (int l = lmin; l <= lmax; l += dl) {
222 const auto ql = this->Q(l, a, b, d, c); // exchange
223 if (ql == 0.0)
224 continue;
225 const auto sixj = sj ? sj->get_2(tja, tjc, 2 * k, tjb, tjd, 2 * l) :
226 Angular::sixj_2(tja, tjc, 2 * k, tjb, tjd, 2 * l);
227 Pk_abcd += sixj * ql;
228 }
229 Pk_abcd *= double(2 * k + 1);
230 return Pk_abcd;
231}
232
233//------------------------------------------------------------------------------
234
235template <Symmetry S>
236double CoulombTable<S>::P(int k, const DiracSpinor &a, const DiracSpinor &b,
237 const DiracSpinor &c, const DiracSpinor &d,
238 const Angular::SixJTable *const sj) const {
239
240 return P(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index(), sj);
241}
242
243//==============================================================================
244template <Symmetry S>
245double CoulombTable<S>::P2(int k, const DiracSpinor &a, const DiracSpinor &b,
246 const DiracSpinor &c, const DiracSpinor &d,
247 const Angular::SixJTable &sj,
248 const std::vector<double> &fk) const {
249 double Pk_abcd{0.0};
250
251 // 6j(s) Triads: {a,c,k}, {k,b,d}, {c,b,l}, {d,a,l}
252 if (Coulomb::triangle(a, c, k) == 0 || Coulomb::triangle(k, b, d) == 0)
253 return 0.0;
254
255 // We may store non-Coulomb integrals in this table (for example, Breit)
256 // Then, the parity selection rule doesn't apply,
257 // and we don't know what the actual rule is, so use whole range
258 const auto [lmin, lmax] = [&]() -> std::pair<int, int> {
259 if constexpr (S == Symmetry::Qk)
260 return k_minmax_Q(a, b, d, c);
261 else
262 return {0, max_k()};
263 }();
264 const int dl = (S == Symmetry::Qk) ? 2 : 1;
265
266 for (int l = lmin; l <= lmax; l += dl) {
267 const auto ql = this->Q(l, a, b, d, c); // exchange
268 if (ql == 0.0)
269 continue;
270 const auto sixj = sj.get(a, c, k, b, d, l);
271
272 // include effective Coulomb screening:
273 const auto f_scr_l = (l < (int)fk.size()) ? fk[std::size_t(l)] : 1.0;
274
275 Pk_abcd += f_scr_l * sixj * ql;
276 }
277 Pk_abcd *= double(2 * k + 1);
278 return Pk_abcd;
279}
280
281//==============================================================================
282template <Symmetry S>
283double CoulombTable<S>::W(int k, const DiracSpinor &a, const DiracSpinor &b,
284 const DiracSpinor &c, const DiracSpinor &d,
285 const Angular::SixJTable *const sj) const {
286 return W(k, a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index(), sj);
287}
288
289template <Symmetry S>
290double CoulombTable<S>::W(int k, nkIndex a, nkIndex b, nkIndex c, nkIndex d,
291 const Angular::SixJTable *const sj) const {
292 return Q(k, a, b, c, d) + P(k, a, b, c, d, sj);
293}
294
295//==============================================================================
296template <Symmetry S>
297double CoulombTable<S>::g(const DiracSpinor &a, const DiracSpinor &b,
298 const DiracSpinor &c, const DiracSpinor &d, int tma,
299 int tmb, int tmc, int tmd) const {
300
301 if (tmc - tma != tmb - tmd)
302 return 0.0;
303 const int twoq = tmc - tma;
304
305 //
306 double g = 0.0;
307 const auto [k0, ki] = [&]() -> std::pair<int, int> {
308 if constexpr (S == Symmetry::Qk)
309 return k_minmax_Q(a, b, c, d);
310 else
311 return {0, max_k()};
312 }();
313 const int dk = (S == Symmetry::Qk) ? 2 : 1;
314
315 for (int k = k0; k <= ki; k += dk) {
316 if (std::abs(twoq) > 2 * k)
317 continue;
318
319 // const auto A =
320 // Angular::neg1pow_2(twoq) *
321 // Angular::Ck_kk_mmq(k, a.kappa(), c.kappa(), tma, tmc, -twoq) *
322 // Angular::Ck_kk_mmq(k, b.kappa(), d.kappa(), tmb, tmd, twoq);
323 // if (A != 0.0) {
324 // g += A * R(k, a, b, c, d);
325 // }
326
327 const auto tjs1 =
328 Angular::threej_2(a.twoj(), 2 * k, c.twoj(), -tma, -twoq, tmc);
329 const auto tjs2 =
330 Angular::threej_2(b.twoj(), 2 * k, d.twoj(), -tmb, twoq, tmd);
331 const auto s = Angular::neg1pow_2(2 * k + tmb - tma + twoq);
332 g += s * tjs1 * tjs2 * Q(k, a, b, c, d);
333 }
334 return g;
335}
336
337//==============================================================================
338//==============================================================================
339
340//==============================================================================
341template <>
342inline nk4Index
343CoulombTable<Symmetry::Qk>::NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
344 nkIndex d) {
345
346 // abcd -> ijkl, with i = min(a,b,c,d)
347
348 // abcd = adcb
349 // badc = bcda
350 // cbad = cdab
351 // dabc = dcba
352 // nb: there must be a more efficient way of doing this!
353
354 // const auto tmp1 = FormIndex(a, b, c, d);
355 // const auto tmp2 = FormIndex(a, d, c, b);
356 // const auto tmp3 = FormIndex(b, a, d, c);
357 // const auto tmp4 = FormIndex(b, c, d, a);
358 // const auto tmp5 = FormIndex(c, b, a, d);
359 // const auto tmp6 = FormIndex(c, d, a, b);
360 // const auto tmp7 = FormIndex(d, a, b, c);
361 // const auto tmp8 = FormIndex(d, c, b, a);
362 // return std::min({tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8});
363
364 nk4Index out = FormIndex(a, b, c, d);
365 const auto min = std::min({a, b, c, d});
366 if (a == min) {
367 const auto tmp = FormIndex(a, d, c, b);
368 if (tmp < out)
369 out = tmp;
370 }
371 if (b == min) {
372 const auto tmp = std::min(FormIndex(b, a, d, c), FormIndex(b, c, d, a));
373 if (tmp < out)
374 out = tmp;
375 }
376 if (c == min) {
377 const auto tmp = std::min(FormIndex(c, b, a, d), FormIndex(c, d, a, b));
378 if (tmp < out)
379 out = tmp;
380 }
381 if (d == min) {
382 const auto tmp = std::min(FormIndex(d, a, b, c), FormIndex(d, c, b, a));
383 if (tmp < out)
384 out = tmp;
385 }
386 return out;
387}
388//------------------------------------------------------------------------------
389template <>
390inline nk4Index
391CoulombTable<Symmetry::Wk>::NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
392 nkIndex d) {
393 const auto tmp1 = FormIndex(a, b, c, d);
394 const auto tmp2 = FormIndex(b, a, d, c);
395 const auto tmp3 = FormIndex(c, d, a, b);
396 const auto tmp4 = FormIndex(d, c, b, a);
397 return std::min({tmp1, tmp2, tmp3, tmp4});
398}
399
400//------------------------------------------------------------------------------
401template <>
402inline nk4Index
403CoulombTable<Symmetry::Lk>::NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
404 nkIndex d) {
405 const auto tmp1 = FormIndex(a, b, c, d);
406 const auto tmp2 = FormIndex(b, a, d, c);
407 return std::min({tmp1, tmp2});
408}
409
410//------------------------------------------------------------------------------
411template <>
412inline nk4Index
413CoulombTable<Symmetry::none>::NormalOrder_impl(nkIndex a, nkIndex b, nkIndex c,
414 nkIndex d) {
415 return FormIndex(a, b, c, d);
416}
417
418//==============================================================================
419template <Symmetry S>
420nk4Index CoulombTable<S>::NormalOrder(nkIndex a, nkIndex b, nkIndex c,
421 nkIndex d) const {
422 return NormalOrder_impl(a, b, c, d);
423}
424//------------------------------------------------------------------------------
425template <Symmetry S>
426nk4Index
427CoulombTable<S>::NormalOrder(const DiracSpinor &a, const DiracSpinor &b,
428 const DiracSpinor &c, const DiracSpinor &d) const {
429 return NormalOrder(a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
430}
431
432//==============================================================================
433template <Symmetry S>
434nk4Index CoulombTable<S>::CurrentOrder(nkIndex a, nkIndex b, nkIndex c,
435 nkIndex d) const {
436 return FormIndex(a, b, c, d);
437}
438//------------------------------------------------------------------------------
439template <Symmetry S>
440nk4Index CoulombTable<S>::CurrentOrder(const DiracSpinor &a,
441 const DiracSpinor &b,
442 const DiracSpinor &c,
443 const DiracSpinor &d) const {
444 return CurrentOrder(a.nk_index(), b.nk_index(), c.nk_index(), d.nk_index());
445}
446
447//==============================================================================
448template <Symmetry S>
449nk4Index CoulombTable<S>::FormIndex(nkIndex a, nkIndex b, nkIndex c,
450 nkIndex d) {
451 // create a nk4Index from four small nkIndex, by laying out in memory
452
453 // static_assert(sizeof(nk4Index) == 4 * sizeof(nkIndex));
454 // nk4Index big_index;
455 // const auto p_index = reinterpret_cast<nkIndex *>(&big_index);
456 // std::memcpy(p_index, &d, sizeof(nkIndex));
457 // std::memcpy(p_index + 1, &c, sizeof(nkIndex));
458 // std::memcpy(p_index + 2, &b, sizeof(nkIndex));
459 // std::memcpy(p_index + 3, &a, sizeof(nkIndex));
460 // return big_index;
461
462 // this seems slightly faster...though variance large
463 static_assert(sizeof(nk4Index) == 4 * sizeof(nkIndex));
464 static_assert(sizeof(nkIndex) * 8 == 16);
465 return ((nk4Index)a << 48) + ((nk4Index)b << 32) + ((nk4Index)c << 16) +
466 (nk4Index)d;
467}
468
469//==============================================================================
470template <Symmetry S>
471std::array<nkIndex, 4>
472CoulombTable<S>::UnFormIndex(const nk4Index &index) const {
473 static_assert(sizeof(nk4Index) == sizeof(std::array<nkIndex, 4>));
474 std::array<nkIndex, 4> set;
475 std::memcpy(&set, &index, sizeof(nk4Index));
476 // note: this returns {d, c, b, a}, not {a, b, c, d}
477 // A warning: this might depend on "endiness"
478 // It's OK, this function is never used???
479 // If there's an issue, a unit test will fail, so it's "safe"
480 std::swap(set[0], set[3]);
481 std::swap(set[1], set[2]);
482 return set;
483}
484
485//==============================================================================
486template <Symmetry S>
487std::array<std::size_t, 4>
488CoulombTable<S>::count_non_zero_integrals(const std::vector<DiracSpinor> &basis,
489 std::size_t k_cutoff,
490 double eF) const {
491 static_assert(S == Symmetry::Qk);
492
493 // -1 account for parity selection rule
494 // (should be same as 2*max_l, without the -1)
495 const auto tmp_max_k =
496 std::size_t(std::max(DiracSpinor::max_tj(basis), 1) - 1);
497 const auto max_k =
498 (k_cutoff <= 0) ? tmp_max_k : std::min(tmp_max_k, std::size_t(k_cutoff));
499
500 std::vector<std::size_t> count_non_zero_k_tot(max_k + 1);
501 std::vector<std::size_t> count_non_zero_k_Excited(max_k + 1);
502 std::vector<std::size_t> count_non_zero_k_MBPT(max_k + 1);
503
504 const auto select_Q_Excited = [eF](int, const DiracSpinor &s,
505 const DiracSpinor &t, const DiracSpinor &u,
506 const DiracSpinor &v) {
507 auto num = Coulomb::number_below_Fermi(s, t, u, v, eF);
508 return num == 0;
509 };
510
511 const auto select_Q_MBPT = [eF](int, const DiracSpinor &s,
512 const DiracSpinor &t, const DiracSpinor &u,
513 const DiracSpinor &v) {
514 auto num = Coulomb::number_below_Fermi(s, t, u, v, eF);
515 return num == 1 || num == 2;
516 };
517
518#pragma omp parallel for schedule(dynamic, 1)
519 for (auto k = 0ul; k <= max_k; ++k) {
520 const auto ik = static_cast<int>(k);
521 for (const auto &a : basis) {
522 for (const auto &b : basis) {
523 for (const auto &c : basis) {
524 if (!Angular::Ck_kk_SR(ik, a.kappa(), c.kappa()))
525 continue;
526 for (const auto &d : basis) {
527 // due to symmetry, only calculate each 'unique' integral once
528 if (NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d)) {
529 if (Angular::Ck_kk_SR(ik, b.kappa(), d.kappa())) {
530 ++count_non_zero_k_tot[k];
531 if (eF != 0.0) {
532 if (select_Q_Excited(0, a, b, c, d)) {
533 ++count_non_zero_k_Excited[k];
534 }
535 if (select_Q_MBPT(0, a, b, c, d)) {
536 ++count_non_zero_k_MBPT[k];
537 }
538 }
539 }
540 }
541 }
542 }
543 }
544 }
545 }
546 // return count_non_zero_k;
547
548 const auto count_non_zero_tot = std::accumulate(
549 count_non_zero_k_tot.begin(), count_non_zero_k_tot.end(), 0ul);
550
551 const auto count_non_zero_Excited = std::accumulate(
552 count_non_zero_k_Excited.begin(), count_non_zero_k_Excited.end(), 0ul);
553
554 const auto count_non_zero_MBPT = std::accumulate(
555 count_non_zero_k_MBPT.begin(), count_non_zero_k_MBPT.end(), 0ul);
556
557 return {count_non_zero_tot, count_non_zero_Excited, count_non_zero_MBPT,
558 max_k};
559}
560
561//==============================================================================
562//==============================================================================
563template <Symmetry S>
564void CoulombTable<S>::fill(const std::vector<DiracSpinor> &basis,
565 const YkTable &yk, int k_cut, bool print) {
566 static_assert(S == Symmetry::Qk);
567 IO::ChronoTimer t(print ? "fill" : "");
568
569 /*
570 In order to make best use of CPU and Memory, we "fill" the QK table in a
571 strange 4-step manner.
572 1. Count the number of non-zero Q's for each k
573 2. Use info from aboe to 'reserve()' space in the map
574 3. Fill all the non-zero Q entries in the map with 0. Note that this adds new
575 elements into the map, so cannot be done in parallel. Since we store a
576 different map for each k, we may parallelise over k (though this is not super
577 efficient)
578 4. Now that we have a fully sized map, we can update each of its values in
579 parallel
580 Note that we make use of the symmetry to ensure we do not access any element
581 from more than 1 thread (also, allowing us to not do any more calculations
582 than necisary)
583 */
584
585 // -1 account for parity selection rule
586 // (should be same as 2*max_l, without the -1)
587 const auto tmp_max_k =
588 std::size_t(std::max(DiracSpinor::max_tj(basis), 1) - 1);
589
590 const auto max_k =
591 (k_cut <= 0) ? tmp_max_k : std::min(tmp_max_k, std::size_t(k_cut));
592
593 if (m_data.size() < max_k + 1)
594 m_data.resize(max_k + 1);
595
596 // 1) Count non-zero Q integrals (each k). Use this to 'reserve' map space=
597 std::vector<std::size_t> count_non_zero_k(max_k + 1);
598#pragma omp parallel for
599 for (auto k = 0ul; k <= max_k; ++k) {
600 const auto ik = static_cast<int>(k);
601 for (const auto &a : basis) {
602 for (const auto &b : basis) {
603 for (const auto &c : basis) {
604 if (!Angular::Ck_kk_SR(ik, a.kappa(), c.kappa()))
605 continue;
606 for (const auto &d : basis) {
607 // due to symmetry, only calculate each 'unique' integral once
608 if (NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d)) {
609 if (Angular::Ck_kk_SR(ik, b.kappa(), d.kappa())) {
610 ++count_non_zero_k[k];
611 }
612 }
613 }
614 }
615 }
616 }
617 }
618
619 // 2) Reserve space in each sub-map
620#pragma omp parallel for
621 for (auto ik = 0ul; ik <= max_k; ++ik) {
622 m_data[ik].reserve(count_non_zero_k[ik]);
623 }
624
625 // 3) Create space in map (set each element to zero).
626#pragma omp parallel for schedule(dynamic, 1)
627 for (auto k = 0ul; k <= max_k; ++k) {
628 for (const auto &a : basis) {
629 for (const auto &c : basis) {
630 if (!Angular::Ck_kk_SR(int(k), a.kappa(), c.kappa()))
631 continue;
632 for (const auto &b : basis) {
633 for (const auto &d : basis) {
634 if (!Angular::Ck_kk_SR(int(k), b.kappa(), d.kappa()))
635 continue;
636 const auto normal_index = NormalOrder(a, b, c, d);
637 if (normal_index == CurrentOrder(a, b, c, d)) {
638 add(int(k), normal_index, 0.0);
639 }
640 }
641 }
642 }
643 }
644 }
645 if (print)
646 std::cout << "Construct empty table: " << t.lap_reading_str() << std::endl;
647
648 // 4) Fill the pre-constructed map with values, in parallel. Since we are
649 // not adding any new elements to map, and since we are guarenteed to only
650 // access each map element once, we can do this part in parallel.
651 if (print)
652 std::cout << "Fill w/ values:\n" << std::flush;
653 t.start();
654 qip::ProgressBar prog(int(basis.size()), print);
655#pragma omp parallel for schedule(dynamic, 1)
656 for (auto ia = 0ul; ia < basis.size(); ++ia) {
657 const auto &a = basis[ia];
658 for (const auto &b : basis) {
659 for (const auto &c : basis) {
660 for (const auto &d : basis) {
661
662 const auto normal_index = NormalOrder(a, b, c, d);
663 if (normal_index != CurrentOrder(a, b, c, d))
664 continue;
665
666 const auto [kmin, kmax] = k_minmax_Q(a, b, c, d);
667 for (int k = kmin; k <= kmax && k <= int(max_k); k += 2) {
668 double *ptr = get(k, normal_index);
669 assert(ptr != nullptr);
670 if (*ptr == 0.0) {
671 // only calculate if not already in table
672 // This saves some time, but surprisingly little!
673 *ptr = yk.Q(k, a, b, c, d);
674 }
675 }
676 }
677 }
678 }
679 prog.update();
680 }
681
682 if (print)
683 std::cout << t.lap_reading_str() << std::endl;
684
685 if (print)
686 summary();
687}
688
689//==============================================================================
690template <Symmetry S>
691void CoulombTable<S>::fill_if(const std::vector<DiracSpinor> &basis,
692 const YkTable &yk,
693 const SelectionRules &SelectionFunction,
694 int k_cut, bool print) {
695 static_assert(S == Symmetry::Qk);
696 IO::ChronoTimer t(print ? "fill" : "");
697
698 /*
699 Same as above, but only calculate integrals satisfying the SelectionFunction
700 */
701
702 const auto tmp_max_k =
703 std::size_t(std::max(DiracSpinor::max_tj(basis), 1) - 1);
704
705 const auto max_k =
706 (k_cut <= 0) ? tmp_max_k : std::min(tmp_max_k, std::size_t(k_cut));
707
708 if (m_data.size() < max_k + 1)
709 m_data.resize(max_k + 1);
710
711 // 1) Count non-zero Q integrals (each k). Use this to 'reserve' map space
712 std::vector<std::size_t> count_non_zero_k(max_k + 1);
713#pragma omp parallel for
714 for (auto k = 0ul; k <= max_k; ++k) {
715 const auto ik = static_cast<int>(k);
716 for (const auto &a : basis) {
717 for (const auto &b : basis) {
718 for (const auto &c : basis) {
719 if (!Angular::Ck_kk_SR(ik, a.kappa(), c.kappa()))
720 continue;
721 for (const auto &d : basis) {
722 // due to symmetry, only calculate each 'unique' integral once
723 if (NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d)) {
724 if (Angular::Ck_kk_SR(ik, b.kappa(), d.kappa()) &&
725 SelectionFunction(ik, a, b, c, d)) {
726 ++count_non_zero_k[k];
727 }
728 }
729 }
730 }
731 }
732 }
733 }
734
735 // 2) Reserve space in each sub-map
736#pragma omp parallel for
737 for (auto ik = 0ul; ik <= max_k; ++ik) {
738 m_data[ik].reserve(count_non_zero_k[ik]);
739 }
740
741 // 3) Create space in map (set each element to zero).
742#pragma omp parallel for
743 for (auto k = 0ul; k <= max_k; ++k) {
744 for (const auto &a : basis) {
745 for (const auto &c : basis) {
746 if (!Angular::Ck_kk_SR(int(k), a.kappa(), c.kappa()))
747 continue;
748 for (const auto &b : basis) {
749 for (const auto &d : basis) {
750 if (!Angular::Ck_kk_SR(int(k), b.kappa(), d.kappa()))
751 continue;
752 if (!SelectionFunction(int(k), a, b, c, d))
753 continue;
754 const auto normal_index = NormalOrder(a, b, c, d);
755 if (normal_index == CurrentOrder(a, b, c, d)) {
756 add(int(k), normal_index, 0.0);
757 }
758 }
759 }
760 }
761 }
762 }
763 if (print)
764 std::cout << "Construct empty table: " << t.lap_reading_str() << std::endl;
765
766 // 4) Fill the pre-constructed map with values, in parallel. Since we are
767 // not adding any new elements to map, and since we are guarenteed to only
768 // access each map element once, we can do this part in parallel.
769 if (print)
770 std::cout << "Fill w/ values:\n" << std::flush;
771 t.start();
772 qip::ProgressBar prog(int(basis.size()), print);
773#pragma omp parallel for schedule(dynamic, 1)
774 for (auto ia = 0ul; ia < basis.size(); ++ia) {
775 const auto &a = basis[ia];
776 for (const auto &b : basis) {
777 for (const auto &c : basis) {
778 for (const auto &d : basis) {
779
780 const auto normal_index = NormalOrder(a, b, c, d);
781 if (normal_index != CurrentOrder(a, b, c, d))
782 continue;
783
784 const auto [kmin, kmax] = k_minmax_Q(a, b, c, d);
785 for (int k = kmin; k <= kmax && k <= int(max_k); k += 2) {
786 if (!SelectionFunction(k, a, b, c, d))
787 continue;
788 double *ptr = get(k, normal_index);
789 assert(ptr != nullptr);
790 if (*ptr == 0.0) {
791 // only calculate if not already in table
792 // This saves some time, but surprisingly little!
793 *ptr = yk.Q(k, a, b, c, d);
794 }
795 }
796 }
797 }
798 }
799 prog.update();
800 }
801
802 if (print)
803 std::cout << t.lap_reading_str() << std::endl;
804
805 if (print)
806 summary();
807}
808
809//==============================================================================
810template <Symmetry S>
811void CoulombTable<S>::fill(const std::vector<DiracSpinor> &basis,
812 const CoulombFunction &Fk,
813 const SelectionRules &Fk_SR, int k_cut, bool print) {
814 IO::ChronoTimer t(print ? "fill" : "");
815
816 /*
817 In order to make best use of CPU and Memory, we "fill" the QK table in a
818 strange 4-step manner.
819 1. Count the number of non-zero Q's for each k
820 2. Use info from aboe to 'reserve()' space in the map
821 3. Fill all the non-zero Q entries in the map with 0. Note that this adds new
822 elements into the map, so cannot be done in parallel. Since we store a
823 different map for each k, we may parallelise over k (though this is not super
824 efficient)
825 4. Now that we have a fully sized map, we can update each of its values in
826 parallel
827 Note that we make use of the symmetry to ensure we do not access any element
828 from more than 1 thread (also, allowing us to not do any more calculations
829 than necisary)
830 */
831
832 const auto tmp_max_k = std::size_t(DiracSpinor::max_tj(basis));
833 // May have different parity rule, so don't -1 here
834
835 const auto max_k =
836 (k_cut <= 0) ? tmp_max_k : std::min(tmp_max_k, std::size_t(k_cut));
837
838 if (m_data.size() < max_k + 1)
839 m_data.resize(max_k + 1);
840
841 // 1) Count non-zero Q integrals (each k). Use this to 'reserve' map space
842 std::vector<std::size_t> count_non_zero_k(max_k + 1);
843#pragma omp parallel for
844 for (auto k = 0ul; k <= max_k; ++k) {
845 const auto ik = static_cast<int>(k);
846 for (const auto &a : basis) {
847 for (const auto &b : basis) {
848 for (const auto &c : basis) {
849 for (const auto &d : basis) {
850 // due to symmetry, only calculate each 'unique' integral once
851 if (NormalOrder(a, b, c, d) == CurrentOrder(a, b, c, d)) {
852 if (Fk_SR(ik, a, b, c, d)) {
853 ++count_non_zero_k[k];
854 }
855 }
856 }
857 }
858 }
859 }
860 }
861
862 // 2) Reserve space in each sub-map
863#pragma omp parallel for
864 for (auto ik = 0ul; ik <= max_k; ++ik) {
865 m_data[ik].reserve(count_non_zero_k[ik]);
866 }
867
868 // 3) Create space in map (set each element to zero).
869#pragma omp parallel for
870 for (auto k = 0ul; k <= max_k; ++k) {
871 for (const auto &a : basis) {
872 for (const auto &b : basis) {
873 for (const auto &c : basis) {
874 for (const auto &d : basis) {
875 if (Fk_SR(int(k), a, b, c, d)) {
876 const auto normal_index = NormalOrder(a, b, c, d);
877 if (normal_index == CurrentOrder(a, b, c, d)) {
878 add(int(k), normal_index, 0.0);
879 }
880 }
881 }
882 }
883 }
884 }
885 }
886 if (print)
887 std::cout << "Construct empty table: " << t.lap_reading_str() << std::endl;
888
889 // 4) Fill the pre-constructed map with values, in parallel. Since we are
890 // not adding any new elements to map, and since we are guarenteed to only
891 // access each map element once, we can do this part in parallel.
892 if (print)
893 std::cout << "Fill w/ values:\n" << std::flush;
894 t.start();
895
896 qip::ProgressBar prog(int(basis.size()), print);
897#pragma omp parallel for schedule(dynamic, 1)
898 for (std::size_t ia = 0; ia < basis.size(); ++ia) {
899
900 const auto &a = basis[ia];
901 for (const auto &b : basis) {
902 for (const auto &c : basis) {
903 for (const auto &d : basis) {
904
905 const auto normal_index = NormalOrder(a, b, c, d);
906 if (normal_index != CurrentOrder(a, b, c, d))
907 continue;
908
909 for (int k = 0; k <= int(max_k); ++k) {
910
911 if (!Fk_SR(k, a, b, c, d))
912 continue;
913 double *ptr = get(k, normal_index);
914 assert(ptr != nullptr);
915 if (*ptr == 0.0) {
916 // only calculate if not already in table
917 *ptr = Fk(k, a, b, c, d);
918 }
919 }
920 }
921 }
922 }
923
924 prog.update();
925 }
926 if (print)
927 std::cout << t.lap_reading_str() << std::endl;
928
929 if (print)
930 summary();
931}
932
933//==============================================================================
934template <Symmetry S>
935void CoulombTable<S>::update(const std::vector<DiracSpinor> &basis,
936 const CoulombFunction &Fk, double damp,
937 bool print) {
938 IO::ChronoTimer t(print ? "update" : "");
939
940 const auto max_k = m_data.size() - 1;
941
942 // zero mean no damping (default)
943 // damp must be between [0, 1) - cannot equal 1 (meaningless)
944 assert(damp != 1.0); // also: catches accidentally passing 'true' (print)
945 const auto x = std::clamp(damp, 0.0, 1.0);
946 const auto y = 1.0 - damp;
947
948 // Fill the pre-constructed map with values, in parallel. Since we are
949 // not adding any new elements to map, and since we are guarenteed to only
950 // access each map element once, we can do this part in parallel.
951
952 qip::ProgressBar prog(int(basis.size()), print);
953#pragma omp parallel for schedule(dynamic, 1)
954 for (std::size_t ia = 0; ia < basis.size(); ++ia) {
955
956 const auto &a = basis[ia];
957 for (const auto &b : basis) {
958 for (const auto &c : basis) {
959 for (const auto &d : basis) {
960
961 // Now, we can check symmetries, selectrion rules etc.
962 // just be seeing if this one is in the current table!
963 const auto c_index = CurrentOrder(a, b, c, d);
964 const auto n_index = NormalOrder(a, b, c, d);
965 if (c_index != n_index)
966 continue;
967
968 for (int k = 0; k <= int(max_k); ++k) {
969 double *ptr = get(k, n_index);
970 // if not null, means already in table: so update it
971 if (ptr != nullptr) {
972 *ptr = (*ptr == 0.0 || x == 0.0) ?
973 Fk(k, a, b, c, d) :
974 y * Fk(k, a, b, c, d) + x * *ptr;
975 }
976 }
977 }
978 }
979 }
980
981 prog.update();
982 }
983}
984
985//==============================================================================
986template <Symmetry S>
987void CoulombTable<S>::write(const std::string &fname, bool verbose) const {
988 if (fname == "false" || fname == "")
989 return;
990 if (verbose)
991 std::cout << "Writing " << count() << " integrals to file: " << fname
992 << ".." << std::flush;
993 std::fstream f;
994 const auto rw = IO::FRW::write;
995 IO::FRW::open_binary(f, fname, rw);
996
997 auto size = m_data.size();
998 rw_binary(f, rw, size);
999 for (const auto &Q_k : m_data) {
1000 auto size_k = Q_k.size();
1001 rw_binary(f, rw, size_k);
1002 for (auto [key, value] : Q_k) {
1003 auto key_copy = key; // have to pass non-const reference!
1004 rw_binary(f, rw, key_copy, value);
1005 }
1006 }
1007 if (verbose)
1008 std::cout << "\n" << std::flush;
1009}
1010
1011//==============================================================================
1012template <Symmetry S>
1013bool CoulombTable<S>::read(const std::string &fname, bool verbose) {
1014 if (fname == "false" || fname.empty())
1015 return false;
1016
1017 std::fstream f;
1018 const auto rw = IO::FRW::read;
1019 IO::FRW::open_binary(f, fname, rw);
1020
1021 // store max_n for each l: allows us to reconstruct the basis string
1022 // This adds noticable overhead, but is useful
1023 std::vector<int> max_n_l;
1024 constexpr bool count_orbs = true;
1025
1026 if (!f.good())
1027 return false;
1028
1029 if (verbose)
1030 std::cout << "Reading: " << fname << "\n" << std::flush;
1031
1032 std::size_t size{0};
1033 rw_binary(f, rw, size);
1034 m_data.resize(size);
1035 if (m_data.size() == 0)
1036 return false;
1037 for (std::size_t i = 0; i < m_data.size(); ++i) {
1038 auto &Q_k = m_data[i];
1039 std::size_t size_k{0};
1040 rw_binary(f, rw, size_k);
1041 Q_k.reserve(size_k);
1042 for (std::size_t ik = 0; ik < size_k; ++ik) {
1043 nk4Index key;
1044 Real value;
1045 rw_binary(f, rw, key, value);
1046 Q_k[key] = value;
1047
1048 // determine basis string (max n for each l)
1049 if constexpr (count_orbs) {
1050 const auto indexes = UnFormIndex(key);
1051 for (const auto index : indexes) {
1052 const auto [n, ki] = Angular::nkindex_to_n_kindex(index);
1053 const auto l = std::size_t(Angular::kindex_to_l(ki));
1054 if (max_n_l.size() < l + 1)
1055 max_n_l.resize(l + 1);
1056 if (n > max_n_l.at(l))
1057 max_n_l.at(l) = n;
1058 }
1059 }
1060 }
1061 }
1062
1063 if (verbose) {
1064 // convert list to basis string (max n for each l)
1065 int last_n = 0;
1066 int l = 0;
1067 std::string basis_str;
1068 for (const auto &n : max_n_l) {
1069 if (n != last_n)
1070 basis_str += std::to_string(n);
1071 basis_str += AtomData::l_symbol(l);
1072 last_n = n;
1073 l++;
1074 }
1075 if (!basis_str.empty())
1076 basis_str = "[" + basis_str + "] ";
1077
1078 std::cout << "Read " << count() << " integrals " << basis_str
1079 << "from file: " << fname << "\n"
1080 << std::flush;
1081 }
1082
1083 return true;
1084}
1085
1086} // namespace Coulomb
Angular::SixJTable
Lookup table for Wigner 6j symbols.
Definition SixJTable.hpp:68
Angular::SixJTable::get_2
double get_2(int a, int b, int c, int d, int e, int f) const
Returns 6j symbol {a/2, b/2, c/2; d/2, e/2, f/2}.
Definition SixJTable.hpp:100
Angular::SixJTable::get
double get(const A &a, const B &b, const C &c, const D &d, const E &e, const F &f) const
"Magic" table lookup: pass integers (for k) or DiracSpinors (for j).
Definition SixJTable.hpp:116
Coulomb::CoulombTable
Base class template to store Coulomb integrals, and similar. 3 specific cases (by template instantiat...
Definition QkTable.hpp:102
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 integral P^k_abcd with effective Coulomb screening.
Definition QkTable.ipp:245
Coulomb::CoulombTable::get
double * get(int k, nk4Index index)
Directly gets one of the stored elements, given normal-ordered nk4Index.
Definition QkTable.ipp:121
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)
Fill table with Coulomb Q integrals satisfying a selection rule.
Definition QkTable.ipp:691
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 integral P^k_abcd = (2k+1) sum_l {6j} Q^l_abdc.
Definition QkTable.ipp:236
Coulomb::CoulombTable::update
void update(const std::vector< DiracSpinor > &basis, const CoulombFunction &Fk, double damp, bool print=true)
Re-calculate all existing table entries using a CoulombFunction.
Definition QkTable.ipp:935
Coulomb::CoulombTable::read
bool read(const std::string &fname, bool verbose=true)
Reads coulomb integrals to disk. Returns false if none read in.
Definition QkTable.ipp:1013
Coulomb::CoulombTable::count_non_zero_integrals
std::array< std::size_t, 4 > count_non_zero_integrals(const std::vector< DiracSpinor > &basis, std::size_t max_k=99, double eF=0.0) const
Definition QkTable.ipp:488
Coulomb::CoulombTable::write
void write(const std::string &fname, bool verbose=true) const
Writes coulomb integrals to disk.
Definition QkTable.ipp:987
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:103
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:427
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:440
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:449
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:184
Coulomb::CoulombTable::UnFormIndex
std::array< nkIndex, 4 > UnFormIndex(const nk4Index &index) const
Breaks nk4Index back into {ia,ib,ic,id}. Not often used.
Definition QkTable.ipp:472
Coulomb::CoulombTable::fill
void fill(const std::vector< DiracSpinor > &basis, const YkTable &yk, int k_cut=-1, bool print=true)
Fill table with all non-zero Coulomb Q integrals from a YkTable.
Definition QkTable.ipp:564
Coulomb::CoulombTable::summary
void summary() const
For testing: prints details of coulomb integrals stored.
Definition QkTable.ipp:19
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
Full matrix element g_abcd with explicit magnetic quantum numbers.
Definition QkTable.ipp:297
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:51
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:147
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
Antisymmetrised integral W^k_abcd = Q^k_abcd + P^k_abcd.
Definition QkTable.ipp:283
Coulomb::CoulombTable::count
std::size_t count() const
Returns number of stored integrals.
Definition QkTable.ipp:35
Coulomb::YkTable
Calculates + stores Hartree Y functions + Angular (w/ look-up), taking advantage of symmetry.
Definition YkTable.hpp:31
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: Ck tables must include all required values,...
Definition YkTable.cpp:129
DiracSpinor
Stores radial Dirac spinor: F_nk = (f, g)
Definition DiracSpinor.hpp:42
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:101
IO::ChronoTimer
Class that uses std::chrono to easily time code.
Definition ChronoTimer.hpp:35
IO::ChronoTimer::lap_reading_str
std::string lap_reading_str() const
Current lap elapsed time as a human-readable string.
Definition ChronoTimer.hpp:157
IO::ChronoTimer::start
void start()
Starts (or resumes) timing.
Definition ChronoTimer.hpp:93
qip::ProgressBar
Thread-safe progress bar for OpenMP parallel loops.
Definition Widgets.hpp:85
qip::ProgressBar::update
void update()
Atomically increment progress counter and print updated bar.
Definition Widgets.hpp:108
Angular::index_to_nk
std::pair< int, int > index_to_nk(std::uint64_t nk_index)
Returns {n, kappa} given nk_index.
Definition Wigner369j.hpp:191
Angular::threej_2
double threej_2(int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3)
Wigner 3j symbol. Inputs are 2*j (or 2*l) as integers.
Definition Wigner369j.hpp:277
Angular::sixj_2
double sixj_2(int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6)
Wigner 6j symbol {j1 j2 j3 | j4 j5 j6}. Inputs are 2*j as integers.
Definition Wigner369j.hpp:431
Angular::nkindex_to_kappa
int nkindex_to_kappa(std::uint64_t nk_index)
Returns kappa, given nk_index.
Definition Wigner369j.hpp:209
Angular::neg1pow
constexpr int neg1pow(int a)
Evaluates (-1)^{a} (for integer a)
Definition Wigner369j.hpp:237
Angular::nkindex_to_n_kindex
std::pair< int, std::uint64_t > nkindex_to_n_kindex(std::uint64_t nk_index)
Returns {n, kappa_index} given nk_index.
Definition Wigner369j.hpp:201
Angular::tildeCk_kk
double tildeCk_kk(int k, int ka, int kb)
Tilde variant of the C^k reduced ME: (-1)^{ja+1/2} * <ka||C^k||kb>
Definition Wigner369j.hpp:557
Angular::Ck_kk_SR
bool Ck_kk_SR(int k, int ka, int kb)
Selection rule check for C^k: returns true if <ka||C^k||kb> may be non-zero.
Definition Wigner369j.hpp:537
Angular::kindex_to_l
constexpr int kindex_to_l(std::uint64_t i)
returns l, given kappa_index
Definition Wigner369j.hpp:116
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:239
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:249
Angular::twoj_k
constexpr int twoj_k(int ka)
returns 2j given kappa
Definition Wigner369j.hpp:46
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 CoulombBreit.cpp:13
Coulomb::Real
double Real
Data type used to store integrals.
Definition QkTable.hpp:32
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:54
Coulomb::to_string
std::string_view to_string(Symmetry s)
Convert Symmetry to string.
Definition QkTable.hpp:17
Coulomb::number_below_Fermi
int number_below_Fermi(const DiracSpinor &i, const DiracSpinor &j, const DiracSpinor &k, const DiracSpinor &l, double eFermi)
Returns number of orbitals that are below Fermi level. Used for Qk selection.
Definition CoulombIntegrals.cpp:443
Coulomb::nk4Index
uint64_t nk4Index
index type for set of 4 orbitals {nk,nk,nk,nk} -> nk4Index [max n: 256]
Definition QkTable.hpp:36
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:352
Coulomb::nkIndex
uint16_t nkIndex
index type for each {nk} (orbital) [max n: 256]
Definition QkTable.hpp:40
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:496
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:48
IO::FRW::open_binary
void open_binary(std::fstream &stream, const std::string &fname, RoW row)
Opens a binary fstream for reading or writing according to row.
Definition FRW_fileReadWrite.hpp:108