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