Line data Source code
1 : /*
2 : * Copyright (c) 2016 Sebastian Weber, Henri Menke. All rights reserved.
3 : *
4 : * This file is part of the pairinteraction library.
5 : *
6 : * The pairinteraction library is free software: you can redistribute it and/or modify
7 : * it under the terms of the GNU Lesser General Public License as published by
8 : * the Free Software Foundation, either version 3 of the License, or
9 : * (at your option) any later version.
10 : *
11 : * The pairinteraction library is distributed in the hope that it will be useful,
12 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : * GNU Lesser General Public License for more details.
15 : *
16 : * You should have received a copy of the GNU Lesser General Public License
17 : * along with the pairinteraction library. If not, see <http://www.gnu.org/licenses/>.
18 : */
19 :
20 : #include "MatrixElements.hpp"
21 :
22 : #include "Constants.hpp"
23 : #include "QuantumDefect.hpp"
24 : #include "SQLite.hpp"
25 :
26 : #include <fmt/format.h>
27 :
28 : #include <cctype>
29 : #include <iostream>
30 : #include <limits>
31 : #include <sstream>
32 : #include <stdexcept>
33 : #include <string>
34 : #include <utility>
35 :
36 65178 : bool selectionRulesMomentum(StateOneOld const &state1, StateOneOld const &state2, int q) {
37 65178 : bool validL = state1.l == state2.l;
38 65178 : bool validJ = std::fabs(state1.j - state2.j) <= 1;
39 65178 : bool validM = state1.m == state2.m + q;
40 65178 : bool validQ = abs(q) <= 1;
41 65178 : return validL && validJ && validM && validQ;
42 : }
43 :
44 8604 : bool selectionRulesMomentum(StateOneOld const &state1, StateOneOld const &state2) {
45 8604 : bool validL = state1.l == state2.l;
46 8604 : bool validJ = std::fabs(state1.j - state2.j) <= 1;
47 8604 : bool validM = (std::fabs(state1.m - state2.m) <= 1);
48 8604 : return validL && validJ && validM;
49 : }
50 :
51 187080 : bool selectionRulesMultipole(StateOneOld const &state1, StateOneOld const &state2, int kappa,
52 : int q) {
53 187080 : bool validL =
54 187080 : (abs(state1.l - state2.l) <= kappa) && (kappa % 2 == abs(state1.l - state2.l) % 2);
55 187080 : bool validJ = (std::fabs(state1.j - state2.j) <= kappa) && (state1.j + state2.j >= kappa);
56 187080 : bool validM = state1.m == state2.m + q;
57 187080 : bool validQ = abs(q) <= kappa;
58 189360 : bool noZero = !(kappa == 2 && state1.j == state2.j && state2.j == 1.5 &&
59 2280 : state1.m == -state2.m && std::fabs(state1.m - state2.m) == 1);
60 187080 : return validL && validJ && validM && validQ && noZero;
61 : }
62 :
63 141090 : bool selectionRulesMultipole(StateOneOld const &state1, StateOneOld const &state2, int kappa) {
64 141090 : bool validL =
65 141090 : (abs(state1.l - state2.l) <= kappa) && (kappa % 2 == abs(state1.l - state2.l) % 2);
66 141090 : bool validJ = (std::fabs(state1.j - state2.j) <= kappa) && (state1.j + state2.j >= kappa);
67 141090 : bool validM = (std::fabs(state1.m - state2.m) <= kappa);
68 141594 : bool noZero = !(kappa == 2 && state1.j == state2.j && state2.j == 1.5 &&
69 504 : state1.m == -state2.m && std::fabs(state1.m - state2.m) == 1);
70 141090 : return validL && validJ && validM && noZero;
71 : }
72 :
73 7 : MatrixElements::MatrixElements(std::string const &species, std::string dbname)
74 7 : : species(species), dbname(std::move(dbname)) {
75 7 : method = "Modelpotentials";
76 7 : muB = 0.5;
77 7 : gS = 2.0023192;
78 7 : gL = 1;
79 :
80 7 : s = 0.5;
81 7 : if (std::isdigit(species.back()) != 0) {
82 0 : s = ((species.back() - '0') - 1) / 2.; // TODO think of a better solution
83 : }
84 7 : }
85 :
86 7 : MatrixElements::MatrixElements(const Configuration &config, std::string const &species,
87 7 : std::string const &dbname)
88 7 : : MatrixElements(species, dbname) {
89 7 : if (config["missingCalc"].str() == "true") {
90 7 : method = "Modelpotentials";
91 0 : } else if (config["missingWhittaker"].str() == "true") {
92 0 : method = "Whittaker";
93 : } else {
94 0 : method = "Error";
95 : }
96 7 : }
97 :
98 4 : void MatrixElements::precalculateMultipole(std::shared_ptr<const BasisnamesOne> const &basis_one,
99 : int k) {
100 4 : int q = std::numeric_limits<int>::max();
101 4 : precalculate(basis_one, k, q, k, true, false, false);
102 4 : }
103 :
104 0 : void MatrixElements::precalculateRadial(std::shared_ptr<const BasisnamesOne> const &basis_one,
105 : int k) {
106 0 : int q = std::numeric_limits<int>::max();
107 0 : precalculate(basis_one, k, q, k, false, false, true);
108 0 : }
109 :
110 3 : void MatrixElements::precalculateElectricMomentum(
111 : std::shared_ptr<const BasisnamesOne> const &basis_one, int q) {
112 3 : precalculate(basis_one, 1, q, 1, true, false, false);
113 3 : }
114 :
115 3 : void MatrixElements::precalculateMagneticMomentum(
116 : std::shared_ptr<const BasisnamesOne> const &basis_one, int q) {
117 3 : precalculate(basis_one, 1, q, 0, false, true, false);
118 3 : }
119 :
120 6 : void MatrixElements::precalculateDiamagnetism(std::shared_ptr<const BasisnamesOne> const &basis_one,
121 : int k, int q) {
122 6 : precalculate(basis_one, k, q, 2, true, false, false);
123 6 : }
124 :
125 1332 : double MatrixElements::getElectricMomentum(
126 : StateOneOld const &state_row,
127 : StateOneOld const &
128 : state_col) { // TODO remove this method and use directly getMultipole(state_row,state_col,1)
129 1332 : return getMultipole(state_row, state_col, 1);
130 : }
131 :
132 1209 : double MatrixElements::getMagneticMomentum(StateOneOld const &state_row,
133 : StateOneOld const &state_col) {
134 2418 : double val = au2GHz / au2G * muB *
135 1209 : cache_radial[0][StateTwoOld({{state_row.n, state_col.n}}, {{state_row.l, state_col.l}},
136 1209 : {{state_row.j, state_col.j}}, {{0, 0}})
137 2418 : .order()] *
138 2418 : cache_angular[1][StateTwoOld({{0, 0}}, {{0, 0}}, {{state_row.j, state_col.j}},
139 1209 : {{state_row.m, state_col.m}})] *
140 2418 : (gL *
141 2418 : cache_reduced_commutes_s[1][StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}},
142 1209 : {{state_row.j, state_col.j}}, {{0, 0}})] *
143 1209 : std::sqrt(state_row.l * (state_row.l + 1) * (2 * state_row.l + 1)) +
144 2418 : gS *
145 2418 : cache_reduced_commutes_l[1][StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}},
146 1209 : {{state_row.j, state_col.j}}, {{0, 0}})] *
147 : std::sqrt(
148 : 0.5 * (0.5 + 1) *
149 : (2 * 0.5 +
150 1209 : 1))); // TODO give the radial matrix element the unit au2GHz / au2G * muB (for the
151 : // electric momentum, the radial matrix element already carries the unit)
152 1209 : return val;
153 : }
154 :
155 3045 : double MatrixElements::getDiamagnetism(StateOneOld const &state_row, StateOneOld const &state_col,
156 : int k) {
157 3045 : double val = elementary_charge * inverse_electron_rest_mass * 1. / 12. *
158 3045 : cache_radial[2][StateTwoOld({{state_row.n, state_col.n}}, {{state_row.l, state_col.l}},
159 3045 : {{state_row.j, state_col.j}}, {{0, 0}})
160 6090 : .order()] *
161 6090 : cache_angular[k][StateTwoOld({{0, 0}}, {{0, 0}}, {{state_row.j, state_col.j}},
162 3045 : {{state_row.m, state_col.m}})] *
163 6090 : cache_reduced_commutes_s[k][StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}},
164 3045 : {{state_row.j, state_col.j}}, {{0, 0}})] *
165 6090 : cache_reduced_multipole[k][StateTwoOld(
166 3045 : {{0, 0}}, {{state_row.l, state_col.l}}, {{0, 0}},
167 3045 : {{0, 0}})]; // TODO do the multiplication with inverse_electron_rest_mass outside this
168 : // class // TODO remove this method and use directly
169 : // getMultipole(state_row,state_col,kappa_angular = ..., kappa_radial = 2)
170 3045 : return val;
171 : }
172 :
173 1588 : double MatrixElements::getMultipole(StateOneOld const &state_row, StateOneOld const &state_col,
174 : int k) {
175 : double val =
176 1588 : cache_radial[k][StateTwoOld({{state_row.n, state_col.n}}, {{state_row.l, state_col.l}},
177 1588 : {{state_row.j, state_col.j}}, {{0, 0}})
178 3176 : .order()] *
179 3176 : cache_angular[k][StateTwoOld({{0, 0}}, {{0, 0}}, {{state_row.j, state_col.j}},
180 1588 : {{state_row.m, state_col.m}})] *
181 3176 : cache_reduced_commutes_s[k][StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}},
182 1588 : {{state_row.j, state_col.j}}, {{0, 0}})] *
183 3176 : cache_reduced_multipole[k][StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}}, {{0, 0}},
184 1588 : {{0, 0}})];
185 1588 : return val;
186 : }
187 :
188 0 : double MatrixElements::getRadial(StateOneOld const &state_row, StateOneOld const &state_col,
189 : int k) {
190 : double val =
191 0 : cache_radial[k][StateTwoOld({{state_row.n, state_col.n}}, {{state_row.l, state_col.l}},
192 0 : {{state_row.j, state_col.j}}, {{0, 0}})
193 0 : .order()];
194 0 : return val;
195 : }
196 :
197 16 : void MatrixElements::precalculate(const std::shared_ptr<const BasisnamesOne> &basis_one, int kappa,
198 : int q, int kappar, bool calcElectricMultipole,
199 : bool calcMagneticMomentum, bool calcRadial) {
200 32 : sqlite::handle db(dbname);
201 32 : sqlite::statement stmt(db);
202 :
203 : // --- create cache tables if necessary (reduced_moemntumS and reduced_moemntumL need not to be
204 : // cached since they are trivial) ---
205 :
206 16 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
207 16 : stmt.exec("CREATE TABLE IF NOT EXISTS cache_radial ("
208 : "method text, species text, k integer, n1 integer, l1 integer, j1 double,"
209 : "n2 integer, l2 integer, j2 double, value double, UNIQUE (method, species, k, "
210 32 : "n1, l1, j1, n2, l2, j2));");
211 16 : stmt.exec("CREATE TEMPORARY TABLE tmp_radial ("
212 : "n1 integer, l1 integer, j1 double,"
213 32 : "n2 integer, l2 integer, j2 double);");
214 : }
215 :
216 16 : if (calcElectricMultipole || calcMagneticMomentum) {
217 16 : stmt.exec("CREATE TABLE IF NOT EXISTS cache_angular ("
218 : "k integer, j1 double, m1 double,"
219 32 : "j2 double, m2 double, value double, UNIQUE (k, j1, m1, j2, m2));");
220 16 : stmt.exec("CREATE TEMPORARY TABLE tmp_angular ("
221 : "j1 double, m1 double,"
222 32 : "j2 double, m2 double);");
223 : }
224 :
225 16 : if (calcElectricMultipole || calcMagneticMomentum) {
226 16 : stmt.exec("CREATE TABLE IF NOT EXISTS cache_reduced_commutes_s ("
227 : "k integer, l1 integer, j1 double,"
228 32 : "l2 integer, j2 double, value double, UNIQUE (k, l1, j1, l2, j2));");
229 16 : stmt.exec("CREATE TEMPORARY TABLE tmp_reduced_commutes_s ("
230 : "l1 integer, j1 double,"
231 32 : "l2 integer, j2 double);");
232 : }
233 :
234 16 : if (calcMagneticMomentum) {
235 3 : stmt.exec("CREATE TABLE IF NOT EXISTS cache_reduced_commutes_l ("
236 : "k integer, l1 integer, j1 double,"
237 6 : "l2 integer, j2 double, value double, UNIQUE (k, l1, j1, l2, j2));");
238 3 : stmt.exec("CREATE TEMPORARY TABLE tmp_reduced_commutes_l ("
239 : "l1 integer, j1 double,"
240 6 : "l2 integer, j2 double);");
241 : }
242 :
243 16 : if (calcElectricMultipole) {
244 13 : stmt.exec("CREATE TABLE IF NOT EXISTS cache_reduced_multipole ("
245 : "k integer, l1 integer,"
246 26 : "l2 integer, value double, UNIQUE (k, l1, l2));");
247 13 : stmt.exec("CREATE TEMPORARY TABLE tmp_reduced_multipole ("
248 : "l1 integer,"
249 26 : "l2 integer);");
250 : }
251 :
252 : // --- determine elements ---
253 :
254 32 : std::stringstream ss;
255 :
256 16 : stmt.exec("begin transaction;");
257 :
258 2264 : for (const auto &state_col : *basis_one) {
259 372600 : for (const auto &state_row : *basis_one) {
260 370352 : if (q != std::numeric_limits<int>::max() && state_row.m - state_col.m != q) {
261 236394 : continue;
262 : }
263 :
264 : // if (state_row.idx < state_col.idx) { // TODO
265 : // continue;
266 : //}
267 :
268 133958 : if ((calcElectricMultipole && selectionRulesMultipole(state_row, state_col, kappa)) ||
269 267916 : (calcMagneticMomentum && selectionRulesMomentum(state_row, state_col)) ||
270 : calcRadial) {
271 23622 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
272 : StateTwoOld state_nlj =
273 0 : StateTwoOld({{state_row.n, state_col.n}}, {{state_row.l, state_col.l}},
274 23622 : {{state_row.j, state_col.j}}, {{0, 0}})
275 47244 : .order();
276 23622 : auto missing_cache_radial = cache_radial[kappar].insert(
277 23622 : {state_nlj, std::numeric_limits<double>::max()});
278 23622 : if (missing_cache_radial.second) {
279 867 : ss.str(std::string());
280 867 : ss << "insert into tmp_radial (n1,l1,j1,n2,l2,j2) values ("
281 867 : << state_nlj.n[0] << "," << state_nlj.l[0] << "," << state_nlj.j[0]
282 867 : << "," << state_nlj.n[1] << "," << state_nlj.l[1] << ","
283 867 : << state_nlj.j[1] << ");";
284 867 : stmt.exec(ss.str());
285 : }
286 : }
287 :
288 23622 : if (calcElectricMultipole || calcMagneticMomentum) {
289 : StateTwoOld state_jm =
290 23622 : StateTwoOld({{0, 0}}, {{0, 0}}, {{state_row.j, state_col.j}},
291 47244 : {{state_row.m, state_col.m}});
292 : auto missing_cache_angular =
293 23622 : cache_angular[kappa].insert({state_jm, std::numeric_limits<double>::max()});
294 23622 : if (missing_cache_angular.second) {
295 992 : ss.str(std::string());
296 992 : ss << "insert into tmp_angular (j1,m1,j2,m2) values (" << state_jm.j[0]
297 992 : << "," << state_jm.m[0] << "," << state_jm.j[1] << "," << state_jm.m[1]
298 992 : << ");";
299 992 : stmt.exec(ss.str());
300 : }
301 : }
302 :
303 23622 : if (calcElectricMultipole || calcMagneticMomentum) {
304 23622 : StateTwoOld state_lj = StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}},
305 47244 : {{state_row.j, state_col.j}}, {{0, 0}});
306 23622 : auto missing_cache_reduced_commutes_s = cache_reduced_commutes_s[kappa].insert(
307 23622 : {state_lj, std::numeric_limits<double>::max()});
308 23622 : if (missing_cache_reduced_commutes_s.second) {
309 123 : ss.str(std::string());
310 123 : ss << "insert into tmp_reduced_commutes_s (l1,j1,l2,j2) values ("
311 123 : << state_lj.l[0] << "," << state_lj.j[0] << "," << state_lj.l[1] << ","
312 123 : << state_lj.j[1] << ");";
313 123 : stmt.exec(ss.str());
314 : }
315 : }
316 :
317 23622 : if (calcMagneticMomentum) {
318 2268 : StateTwoOld state_lj = StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}},
319 4536 : {{state_row.j, state_col.j}}, {{0, 0}});
320 2268 : auto missing_cache_reduced_commutes_l = cache_reduced_commutes_l[kappa].insert(
321 2268 : {state_lj, std::numeric_limits<double>::max()});
322 2268 : if (missing_cache_reduced_commutes_l.second) {
323 17 : ss.str(std::string());
324 17 : ss << "insert into tmp_reduced_commutes_l (l1,j1,l2,j2) values ("
325 17 : << state_lj.l[0] << "," << state_lj.j[0] << "," << state_lj.l[1] << ","
326 17 : << state_lj.j[1] << ");";
327 17 : stmt.exec(ss.str());
328 : }
329 : }
330 :
331 23622 : if (calcElectricMultipole) {
332 : StateTwoOld state_l =
333 42708 : StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}}, {{0, 0}}, {{0, 0}});
334 21354 : auto missing_cache_reduced_multipole = cache_reduced_multipole[kappa].insert(
335 21354 : {state_l, std::numeric_limits<double>::max()});
336 21354 : if (missing_cache_reduced_multipole.second) {
337 39 : ss.str(std::string());
338 39 : ss << "insert into tmp_reduced_multipole (l1,l2) values (" << state_l.l[0]
339 39 : << "," << state_l.l[1] << ");";
340 39 : stmt.exec(ss.str());
341 : }
342 : }
343 : }
344 : }
345 : }
346 :
347 16 : stmt.exec("end transaction;");
348 :
349 : // --- load from database ---
350 :
351 : int n1, n2, l1, l2;
352 : float j1, j2, m1, m2;
353 : double value;
354 :
355 16 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
356 16 : ss.str(std::string());
357 : ss << "SELECT c.n1, c.l1, c.j1, c.n2, c.l2, c.j2, c.value FROM cache_radial c INNER JOIN "
358 : "tmp_radial t ON ("
359 : << "c.n1 = t.n1 AND c.l1 = t.l1 AND c.j1 = t.j1 AND c.n2 = t.n2 AND c.l2 = t.l2 AND "
360 : "c.j2 = t.j2) "
361 32 : << "WHERE c.method= '" << method << "' AND c.species = '" << species
362 16 : << "' AND c.k = " << kappar << ";";
363 32 : sqlite::statement stmt(db, ss.str());
364 16 : stmt.prepare();
365 274 : while (stmt.step()) {
366 258 : n1 = stmt.get<decltype(n1)>(0);
367 258 : l1 = stmt.get<decltype(l1)>(1);
368 258 : j1 = stmt.get<decltype(j1)>(2);
369 258 : n2 = stmt.get<decltype(n2)>(3);
370 258 : l2 = stmt.get<decltype(l2)>(4);
371 258 : j2 = stmt.get<decltype(j2)>(5);
372 258 : value = stmt.get<decltype(value)>(6);
373 258 : cache_radial[kappar][StateTwoOld({{n1, n2}}, {{l1, l2}}, {{j1, j2}}, {{0, 0}})] = value;
374 : }
375 : }
376 :
377 16 : if (calcElectricMultipole || calcMagneticMomentum) {
378 16 : ss.str(std::string());
379 : ss << "SELECT c.j1, c.m1, c.j2, c.m2, c.value FROM cache_angular c INNER JOIN tmp_angular "
380 : "t ON ("
381 : << "c.j1 = t.j1 AND c.m1 = t.m1 AND c.j2 = t.j2 AND c.m2 = t.m2) "
382 16 : << "WHERE c.k = " << kappa << ";";
383 32 : sqlite::statement stmt(db, ss.str());
384 16 : stmt.prepare();
385 188 : while (stmt.step()) {
386 172 : j1 = stmt.get<decltype(j1)>(0);
387 172 : m1 = stmt.get<decltype(m1)>(1);
388 172 : j2 = stmt.get<decltype(j2)>(2);
389 172 : m2 = stmt.get<decltype(m2)>(3);
390 172 : value = stmt.get<decltype(value)>(4);
391 172 : cache_angular[kappa][StateTwoOld({{0, 0}}, {{0, 0}}, {{j1, j2}}, {{m1, m2}})] = value;
392 : }
393 : }
394 :
395 16 : if (calcElectricMultipole || calcMagneticMomentum) {
396 16 : ss.str(std::string());
397 : ss << "SELECT c.l1, c.j1, c.l2, c.j2, c.value FROM cache_reduced_commutes_s c INNER JOIN "
398 : "tmp_reduced_commutes_s t ON ("
399 : << "c.l1 = t.l1 AND c.j1 = t.j1 AND c.l2 = t.l2 AND c.j2 = t.j2) "
400 16 : << "WHERE c.k = " << kappa << ";";
401 32 : sqlite::statement stmt(db, ss.str());
402 16 : stmt.prepare();
403 38 : while (stmt.step()) {
404 22 : l1 = stmt.get<decltype(l1)>(0);
405 22 : j1 = stmt.get<decltype(j1)>(1);
406 22 : l2 = stmt.get<decltype(l2)>(2);
407 22 : j2 = stmt.get<decltype(j2)>(3);
408 22 : value = stmt.get<decltype(value)>(4);
409 22 : cache_reduced_commutes_s[kappa][StateTwoOld({{0, 0}}, {{l1, l2}}, {{j1, j2}},
410 22 : {{0, 0}})] = value;
411 : }
412 : }
413 :
414 16 : if (calcMagneticMomentum) {
415 3 : ss.str(std::string());
416 : ss << "SELECT c.l1, c.j1, c.l2, c.j2, c.value FROM cache_reduced_commutes_l c INNER JOIN "
417 : "tmp_reduced_commutes_l t ON ("
418 : << "c.l1 = t.l1 AND c.j1 = t.j1 AND c.l2 = t.l2 AND c.j2 = t.j2) "
419 3 : << "WHERE c.k = " << kappa << ";";
420 6 : sqlite::statement stmt(db, ss.str());
421 3 : stmt.prepare();
422 3 : while (stmt.step()) {
423 0 : l1 = stmt.get<decltype(l1)>(0);
424 0 : j1 = stmt.get<decltype(j1)>(1);
425 0 : l2 = stmt.get<decltype(l2)>(2);
426 0 : j2 = stmt.get<decltype(j2)>(3);
427 0 : value = stmt.get<decltype(value)>(4);
428 0 : cache_reduced_commutes_l[kappa][StateTwoOld({{0, 0}}, {{l1, l2}}, {{j1, j2}},
429 0 : {{0, 0}})] = value;
430 : }
431 : }
432 :
433 16 : if (calcElectricMultipole) {
434 13 : ss.str(std::string());
435 : ss << "SELECT c.l1, c.l2, c.value FROM cache_reduced_multipole c INNER JOIN "
436 : "tmp_reduced_multipole t ON ("
437 : << "c.l1 = t.l1 AND c.l2 = t.l2) "
438 13 : << "WHERE c.k = " << kappa << ";";
439 26 : sqlite::statement stmt(db, ss.str());
440 13 : stmt.prepare();
441 21 : while (stmt.step()) {
442 8 : l1 = stmt.get<decltype(l1)>(0);
443 8 : l2 = stmt.get<decltype(l2)>(1);
444 8 : value = stmt.get<decltype(value)>(2);
445 8 : cache_reduced_multipole[kappa][StateTwoOld({{0, 0}}, {{l1, l2}}, {{0, 0}}, {{0, 0}})] =
446 : value;
447 : }
448 : }
449 :
450 : // --- calculate missing elements and write them to the database ---
451 :
452 16 : stmt.exec("begin transaction;");
453 :
454 16 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
455 1963 : for (auto &cache : cache_radial[kappar]) {
456 1947 : if (cache.second == std::numeric_limits<double>::max()) {
457 1218 : auto state = cache.first;
458 :
459 1218 : QuantumDefect qd1(species, state.n[0], state.l[0], state.j[0]);
460 609 : QuantumDefect qd2(species, state.n[1], state.l[1], state.j[1]);
461 609 : cache.second = calcRadialElement(qd1, kappar, qd2);
462 :
463 609 : ss.str(std::string());
464 : ss << "insert into cache_radial (method, species, k, n1, l1, j1, n2, l2, j2, "
465 : "value) values ("
466 609 : << "'" << method << "'"
467 : << ","
468 609 : << "'" << species << "'"
469 609 : << "," << kappar << "," << state.n[0] << "," << state.l[0] << "," << state.j[0]
470 609 : << "," << state.n[1] << "," << state.l[1] << "," << state.j[1] << ","
471 609 : << cache.second << ");";
472 609 : stmt.exec(ss.str());
473 : }
474 : }
475 : }
476 :
477 16 : if (calcElectricMultipole || calcMagneticMomentum) {
478 2608 : for (auto &cache : cache_angular[kappa]) {
479 2592 : if (cache.second == std::numeric_limits<double>::max()) {
480 820 : auto state = cache.first;
481 820 : float q = state.m[0] - state.m[1];
482 :
483 820 : cache.second = std::pow(-1, state.j[0] - state.m[0]) *
484 820 : WignerSymbols::wigner3j(state.j[0], kappa, state.j[1], -state.m[0], q,
485 820 : state.m[1]);
486 :
487 820 : ss.str(std::string());
488 820 : ss << "insert into cache_angular (k, j1, m1, j2, m2, value) values (" << kappa
489 820 : << "," << state.j[0] << "," << state.m[0] << "," << state.j[1] << ","
490 820 : << state.m[1] << "," << cache.second << ");";
491 820 : stmt.exec(ss.str());
492 : }
493 : }
494 : }
495 :
496 16 : if (calcElectricMultipole || calcMagneticMomentum) {
497 407 : for (auto &cache : cache_reduced_commutes_s[kappa]) {
498 391 : if (cache.second == std::numeric_limits<double>::max()) {
499 101 : auto state = cache.first;
500 :
501 : // Check triangle conditions (violated e.g. for WignerSymbols::wigner6j(0, 0.5, 0.5,
502 : // 0.5, 0, 1))
503 202 : if (state.l[0] >= std::abs(state.j[0] - 0.5) && state.l[0] <= state.j[0] + 0.5 &&
504 101 : state.l[0] >= std::abs(state.l[1] - kappa) &&
505 200 : state.l[0] <= state.l[1] + kappa &&
506 100 : state.j[1] >= std::abs(state.j[0] - kappa) &&
507 302 : state.j[1] <= state.j[0] + kappa && state.j[1] >= std::abs(state.l[1] - 0.5) &&
508 100 : state.j[1] <= state.l[1] + 0.5) {
509 100 : cache.second = std::pow(-1, state.l[0] + 0.5 + state.j[1] + kappa) *
510 100 : std::sqrt((2 * state.j[0] + 1) * (2 * state.j[1] + 1)) *
511 100 : WignerSymbols::wigner6j(state.l[0], state.j[0], 0.5, state.j[1], state.l[1],
512 : kappa);
513 : } else {
514 1 : cache.second = 0;
515 : }
516 :
517 101 : ss.str(std::string());
518 101 : ss << "insert into cache_reduced_commutes_s (k, l1, j1, l2, j2, value) values ("
519 101 : << kappa << "," << state.l[0] << "," << state.j[0] << "," << state.l[1] << ","
520 101 : << state.j[1] << "," << cache.second << ");";
521 101 : stmt.exec(ss.str());
522 : }
523 : }
524 : }
525 :
526 16 : if (calcMagneticMomentum) {
527 54 : for (auto &cache : cache_reduced_commutes_l[kappa]) {
528 51 : if (cache.second == std::numeric_limits<double>::max()) {
529 17 : auto state = cache.first;
530 :
531 17 : cache.second = std::pow(-1, state.l[0] + 0.5 + state.j[0] + kappa) *
532 17 : std::sqrt((2 * state.j[0] + 1) * (2 * state.j[1] + 1)) *
533 17 : WignerSymbols::wigner6j(0.5, state.j[0], state.l[0], state.j[1], 0.5, kappa);
534 :
535 17 : ss.str(std::string());
536 17 : ss << "insert into cache_reduced_commutes_l (k, l1, j1, l2, j2, value) values ("
537 17 : << kappa << "," << state.l[0] << "," << state.j[0] << "," << state.l[1] << ","
538 17 : << state.j[1] << "," << cache.second << ");";
539 17 : stmt.exec(ss.str());
540 : }
541 : }
542 : }
543 :
544 16 : if (calcElectricMultipole) {
545 108 : for (auto &cache : cache_reduced_multipole[kappa]) {
546 95 : if (cache.second == std::numeric_limits<double>::max()) {
547 31 : auto state = cache.first;
548 :
549 31 : cache.second = std::pow(-1, state.l[0]) *
550 31 : std::sqrt((2 * state.l[0] + 1) * (2 * state.l[1] + 1)) *
551 31 : WignerSymbols::wigner3j(state.l[0], kappa, state.l[1], 0, 0, 0);
552 :
553 31 : ss.str(std::string());
554 31 : ss << "insert into cache_reduced_multipole (k, l1, l2, value) values (" << kappa
555 31 : << "," << state.l[0] << "," << state.l[1] << "," << cache.second << ");";
556 31 : stmt.exec(ss.str());
557 : }
558 : }
559 : }
560 :
561 16 : stmt.exec("end transaction;");
562 16 : }
563 :
564 609 : double MatrixElements::calcRadialElement(const QuantumDefect &qd1, int power,
565 : const QuantumDefect &qd2) {
566 :
567 609 : double converter = (power == 0) ? 1 : au2GHz / au2Vcm * std::pow(au2um, power - 1);
568 :
569 609 : if (method == "Modelpotentials") {
570 609 : return converter * IntegrateRadialElement<Numerov>(qd1, power, qd2);
571 : }
572 0 : if (method == "Whittaker") {
573 0 : return converter * IntegrateRadialElement<Whittaker>(qd1, power, qd2);
574 : }
575 : std::string msg("You have to provide all radial matrix elements on your own because you have "
576 0 : "deactivated the calculation of missing radial matrix elements!");
577 0 : std::cout << msg << std::endl;
578 0 : throw std::runtime_error(msg);
579 : }
580 :
581 0 : void MatrixElements::precalculateMultipole(const std::vector<StateOneOld> &basis_one, int k) {
582 0 : int q = std::numeric_limits<int>::max();
583 0 : precalculate(basis_one, k, q, k, true, false, false);
584 0 : }
585 :
586 0 : void MatrixElements::precalculateRadial(const std::vector<StateOneOld> &basis_one, int k) {
587 0 : int q = std::numeric_limits<int>::max();
588 0 : precalculate(basis_one, k, q, k, false, false, true);
589 0 : }
590 :
591 0 : void MatrixElements::precalculateElectricMomentum(const std::vector<StateOneOld> &basis_one,
592 : int q) {
593 0 : precalculate(basis_one, 1, q, 1, true, false, false);
594 0 : }
595 :
596 0 : void MatrixElements::precalculateMagneticMomentum(const std::vector<StateOneOld> &basis_one,
597 : int q) {
598 0 : precalculate(basis_one, 1, q, 0, false, true, false);
599 0 : }
600 :
601 0 : void MatrixElements::precalculateDiamagnetism(const std::vector<StateOneOld> &basis_one, int k,
602 : int q) {
603 0 : precalculate(basis_one, k, q, 2, true, false, false);
604 0 : }
605 :
606 0 : void MatrixElements::precalculate(const std::vector<StateOneOld> &basis_one, int kappa, int q,
607 : int kappar, bool calcElectricMultipole, bool calcMagneticMomentum,
608 : bool calcRadial) {
609 0 : sqlite::handle db(dbname);
610 0 : sqlite::statement stmt(db);
611 :
612 : // --- Speed up database access ---
613 :
614 0 : stmt.exec("PRAGMA synchronous = OFF"); // do not wait on write, hand off to OS and carry on
615 0 : stmt.exec("PRAGMA journal_mode = MEMORY"); // keep rollback journal in memory during transaction
616 :
617 : // --- Create cache tables if necessary (reduced_moemntum_s and reduced_moemntum_l need not to
618 : // be cached since they are trivial) --- // TODO put this into the constructor of the
619 : // prospective cache object
620 :
621 0 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
622 0 : stmt.exec("create table if not exists cache_radial ("
623 : "method text, species text, k integer, n1 integer, l1 integer, j1 double,"
624 : "n2 integer, l2 integer, j2 double, value double, primary key (method, species, "
625 0 : "k, n1, l1, j1, n2, l2, j2)) without rowid;");
626 : }
627 :
628 0 : if (calcElectricMultipole || calcMagneticMomentum) {
629 0 : stmt.exec(
630 : "create table if not exists cache_angular ("
631 : "k integer, j1 double, m1 double,"
632 0 : "j2 double, m2 double, value double, primary key (k, j1, m1, j2, m2)) without rowid;");
633 : }
634 :
635 0 : if (calcElectricMultipole || calcMagneticMomentum) {
636 0 : stmt.exec("create table if not exists cache_reduced_commutes_s ("
637 : "s double, k integer, l1 integer, j1 double,"
638 : "l2 integer, j2 double, value double, primary key (s, k, l1, j1, l2, j2)) "
639 0 : "without rowid;");
640 : }
641 :
642 0 : if (calcMagneticMomentum) {
643 0 : stmt.exec("create table if not exists cache_reduced_commutes_l ("
644 : "s double, k integer, l1 integer, j1 double,"
645 : "l2 integer, j2 double, value double, primary key (s, k, l1, j1, l2, j2)) "
646 0 : "without rowid;");
647 : }
648 :
649 0 : if (calcElectricMultipole) {
650 0 : stmt.exec("create table if not exists cache_reduced_multipole ("
651 : "k integer, l1 integer,"
652 0 : "l2 integer, value double, primary key (k, l1, l2)) without rowid;");
653 : }
654 :
655 : // --- Determine elements ---
656 :
657 0 : for (const auto &state_col : basis_one) {
658 0 : for (const auto &state_row : basis_one) {
659 0 : if (q != std::numeric_limits<int>::max() && state_row.m - state_col.m != q) {
660 0 : continue;
661 : }
662 :
663 0 : if (state_row.species.empty() || state_col.species.empty()) {
664 0 : continue; // TODO artifical states !!!
665 : }
666 :
667 : // if (state_row.idx < state_col.idx) { // TODO
668 : // continue;
669 : //}
670 :
671 0 : if ((calcElectricMultipole && selectionRulesMultipole(state_row, state_col, kappa)) ||
672 0 : (calcMagneticMomentum && selectionRulesMomentum(state_row, state_col)) ||
673 : calcRadial) {
674 0 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
675 : StateTwoOld state_nlj =
676 0 : StateTwoOld({{state_row.n, state_col.n}}, {{state_row.l, state_col.l}},
677 0 : {{state_row.j, state_col.j}}, {{0, 0}})
678 0 : .order();
679 0 : cache_radial[kappar].insert(
680 : {state_nlj,
681 0 : std::numeric_limits<double>::max()}); // TODO can this be speed up?
682 : }
683 :
684 0 : if (calcElectricMultipole || calcMagneticMomentum) {
685 : StateTwoOld state_jm =
686 0 : StateTwoOld({{0, 0}}, {{0, 0}}, {{state_row.j, state_col.j}},
687 0 : {{state_row.m, state_col.m}});
688 0 : cache_angular[kappa].insert({state_jm, std::numeric_limits<double>::max()});
689 : }
690 :
691 0 : if (calcElectricMultipole || calcMagneticMomentum) {
692 0 : StateTwoOld state_lj = StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}},
693 0 : {{state_row.j, state_col.j}}, {{0, 0}});
694 0 : cache_reduced_commutes_s[kappa].insert(
695 0 : {state_lj, std::numeric_limits<double>::max()});
696 : }
697 :
698 0 : if (calcMagneticMomentum) {
699 0 : StateTwoOld state_lj = StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}},
700 0 : {{state_row.j, state_col.j}}, {{0, 0}});
701 0 : cache_reduced_commutes_l[kappa].insert(
702 0 : {state_lj, std::numeric_limits<double>::max()});
703 : }
704 :
705 0 : if (calcElectricMultipole) {
706 : StateTwoOld state_l =
707 0 : StateTwoOld({{0, 0}}, {{state_row.l, state_col.l}}, {{0, 0}}, {{0, 0}});
708 0 : cache_reduced_multipole[kappa].insert(
709 0 : {state_l, std::numeric_limits<double>::max()});
710 : }
711 : }
712 : }
713 : }
714 :
715 : // --- Load from database ---
716 :
717 0 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
718 0 : stmt.set(
719 : "select value from cache_radial where `method` = ?1 and `species` = ?2 and `k` = ?3 "
720 0 : "and `n1` = ?4 and `l1` = ?5 and `j1` = ?6 and `n2` = ?7 and `l2` = ?8 and `j2` = ?9;");
721 :
722 0 : stmt.prepare();
723 0 : for (auto &cache : cache_radial[kappar]) {
724 0 : auto state = cache.first;
725 0 : stmt.bind(1, method);
726 0 : stmt.bind(2, species);
727 0 : stmt.bind(3, kappar);
728 0 : stmt.bind(4, state.n[0]);
729 0 : stmt.bind(5, state.l[0]);
730 0 : stmt.bind(6, state.j[0]);
731 0 : stmt.bind(7, state.n[1]);
732 0 : stmt.bind(8, state.l[1]);
733 0 : stmt.bind(9, state.j[1]);
734 0 : if (stmt.step()) {
735 0 : cache.second = stmt.get<double>(0);
736 : }
737 0 : stmt.reset();
738 : }
739 : }
740 :
741 0 : if (calcElectricMultipole || calcMagneticMomentum) {
742 0 : stmt.set("select value from cache_angular where `k` = ?1 and `j1` = ?2 and `m1` = ?3 and "
743 0 : "`j2` = ?4 and `m2` = ?5;");
744 :
745 0 : stmt.prepare();
746 0 : for (auto &cache : cache_angular[kappa]) {
747 0 : auto state = cache.first;
748 0 : stmt.bind(1, kappa);
749 0 : stmt.bind(2, state.j[0]);
750 0 : stmt.bind(3, state.m[0]);
751 0 : stmt.bind(4, state.j[1]);
752 0 : stmt.bind(5, state.m[1]);
753 0 : if (stmt.step()) {
754 0 : cache.second = stmt.get<double>(0);
755 : }
756 0 : stmt.reset();
757 : }
758 : }
759 :
760 0 : if (calcElectricMultipole || calcMagneticMomentum) {
761 0 : stmt.set("select value from cache_reduced_commutes_s where `s` = ?1 and `k` = ?2 and `l1` "
762 0 : "= ?3 and `j1` = ?4 and `l2` = ?5 and `j2` = ?6;");
763 :
764 0 : stmt.prepare();
765 0 : for (auto &cache : cache_reduced_commutes_s[kappa]) {
766 0 : auto state = cache.first;
767 0 : stmt.bind(1, s);
768 0 : stmt.bind(2, kappa);
769 0 : stmt.bind(3, state.l[0]);
770 0 : stmt.bind(4, state.j[0]);
771 0 : stmt.bind(5, state.l[1]);
772 0 : stmt.bind(6, state.j[1]);
773 0 : if (stmt.step()) {
774 0 : cache.second = stmt.get<double>(0);
775 : }
776 0 : stmt.reset();
777 : }
778 : }
779 :
780 0 : if (calcMagneticMomentum) {
781 0 : stmt.set("select value from cache_reduced_commutes_l where `s` = ?1 and `k` = ?2 and `l1` "
782 0 : "= ?3 and `j1` = ?4 and `l2` = ?5 and `j2` = ?6;");
783 :
784 0 : stmt.prepare();
785 0 : for (auto &cache : cache_reduced_commutes_l[kappa]) {
786 0 : auto state = cache.first;
787 0 : stmt.bind(1, kappa);
788 0 : stmt.bind(2, s);
789 0 : stmt.bind(3, state.l[0]);
790 0 : stmt.bind(4, state.j[0]);
791 0 : stmt.bind(5, state.l[1]);
792 0 : stmt.bind(6, state.j[1]);
793 0 : if (stmt.step()) {
794 0 : cache.second = stmt.get<double>(0);
795 : }
796 0 : stmt.reset();
797 : }
798 : }
799 :
800 0 : if (calcElectricMultipole) {
801 0 : stmt.set("select value from cache_reduced_multipole where `k` = ?1 and `l1` = ?2 and `l2` "
802 0 : "= ?3;");
803 :
804 0 : stmt.prepare();
805 0 : for (auto &cache : cache_reduced_multipole[kappa]) {
806 0 : auto state = cache.first;
807 0 : stmt.bind(1, kappa);
808 0 : stmt.bind(2, state.l[0]);
809 0 : stmt.bind(3, state.l[1]);
810 0 : if (stmt.step()) {
811 0 : cache.second = stmt.get<double>(0);
812 : }
813 0 : stmt.reset();
814 : }
815 : }
816 :
817 : // --- Calculate missing elements and write them to the database ---
818 :
819 0 : stmt.exec("begin transaction;");
820 :
821 0 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
822 0 : stmt.set("insert or ignore into cache_radial (method, species, k, n1, l1, j1, n2, l2, j2, "
823 0 : "value) values (?1, ?2, ?3, ?4, ?5, ?6, ?7, ?8, ?9, ?10);");
824 :
825 0 : stmt.prepare();
826 0 : for (auto &cache : cache_radial[kappar]) {
827 0 : if (cache.second == std::numeric_limits<double>::max()) {
828 0 : auto state = cache.first;
829 :
830 0 : QuantumDefect qd1(species, state.n[0], state.l[0], state.j[0]);
831 0 : QuantumDefect qd2(species, state.n[1], state.l[1], state.j[1]);
832 0 : cache.second = calcRadialElement(qd1, kappar, qd2);
833 :
834 0 : stmt.bind(1, method);
835 0 : stmt.bind(2, species);
836 0 : stmt.bind(3, kappar);
837 0 : stmt.bind(4, state.n[0]);
838 0 : stmt.bind(5, state.l[0]);
839 0 : stmt.bind(6, state.j[0]);
840 0 : stmt.bind(7, state.n[1]);
841 0 : stmt.bind(8, state.l[1]);
842 0 : stmt.bind(9, state.j[1]);
843 0 : stmt.bind(10, cache.second);
844 :
845 0 : stmt.step();
846 0 : stmt.reset();
847 : }
848 : }
849 : }
850 :
851 0 : if (calcElectricMultipole || calcMagneticMomentum) {
852 0 : stmt.set("insert or ignore into cache_angular (k, j1, m1, j2, m2, value) values (?1, ?2, "
853 0 : "?3, ?4, ?5, ?6);");
854 :
855 0 : for (auto &cache : cache_angular[kappa]) {
856 0 : if (cache.second == std::numeric_limits<double>::max()) {
857 0 : auto state = cache.first;
858 0 : float q = state.m[0] - state.m[1];
859 :
860 0 : cache.second = std::pow(-1, state.j[0] - state.m[0]) *
861 0 : WignerSymbols::wigner3j(state.j[0], kappa, state.j[1], -state.m[0], q,
862 0 : state.m[1]);
863 :
864 0 : if (cache.second > 1e9) {
865 : std::cout << "Warning: Error in the calculation of the Wigner-3j-symbol "
866 0 : "detected and resolved."
867 0 : << std::endl;
868 0 : cache.second = 0;
869 : }
870 :
871 0 : stmt.prepare(); // @henri: can this be moved outside the loop?
872 0 : stmt.bind(1, kappa);
873 0 : stmt.bind(2, state.j[0]);
874 0 : stmt.bind(3, state.m[0]);
875 0 : stmt.bind(4, state.j[1]);
876 0 : stmt.bind(5, state.m[1]);
877 0 : stmt.bind(6, cache.second);
878 0 : stmt.step();
879 0 : stmt.reset(); // @henri: can this be moved outside the loop?
880 : }
881 : }
882 : }
883 :
884 0 : if (calcElectricMultipole || calcMagneticMomentum) {
885 0 : stmt.set("insert or ignore into cache_reduced_commutes_s (s, k, l1, j1, l2, j2, value) "
886 0 : "values (?1, ?2, ?3, ?4, ?5, ?6, ?7);");
887 :
888 0 : stmt.prepare();
889 0 : for (auto &cache : cache_reduced_commutes_s[kappa]) {
890 0 : if (cache.second == std::numeric_limits<double>::max()) {
891 0 : auto state = cache.first;
892 :
893 : // Check triangle conditions (violated e.g. for WignerSymbols::wigner6j(0, 0.5, 0.5,
894 : // 0.5, 0, 1))
895 0 : if (state.l[0] >= std::abs(state.j[0] - s) && state.l[0] <= state.j[0] + s &&
896 0 : state.l[0] >= std::abs(state.l[1] - kappa) &&
897 0 : state.l[0] <= state.l[1] + kappa &&
898 0 : state.j[1] >= std::abs(state.j[0] - kappa) &&
899 0 : state.j[1] <= state.j[0] + kappa && state.j[1] >= std::abs(state.l[1] - s) &&
900 0 : state.j[1] <= state.l[1] + s) {
901 0 : cache.second = std::pow(-1, state.l[0] + s + state.j[1] + kappa) *
902 0 : std::sqrt((2 * state.j[0] + 1) * (2 * state.j[1] + 1)) *
903 0 : WignerSymbols::wigner6j(state.l[0], state.j[0], s, state.j[1], state.l[1],
904 : kappa);
905 : } else {
906 0 : cache.second = 0;
907 : }
908 :
909 0 : stmt.bind(1, s);
910 0 : stmt.bind(2, kappa);
911 0 : stmt.bind(3, state.l[0]);
912 0 : stmt.bind(4, state.j[0]);
913 0 : stmt.bind(5, state.l[1]);
914 0 : stmt.bind(6, state.j[1]);
915 0 : stmt.bind(7, cache.second);
916 0 : stmt.step();
917 0 : stmt.reset();
918 : }
919 : }
920 : }
921 :
922 0 : if (calcMagneticMomentum) {
923 0 : stmt.set("insert or ignore into cache_reduced_commutes_l (s, k, l1, j1, l2, j2, value) "
924 0 : "values (?1, ?2, ?3, ?4, ?5, ?6, ?7);");
925 :
926 0 : stmt.prepare();
927 0 : for (auto &cache : cache_reduced_commutes_l[kappa]) {
928 0 : if (cache.second == std::numeric_limits<double>::max()) {
929 0 : auto state = cache.first;
930 :
931 : // Check triangle conditions (violated e.g. for WignerSymbols::wigner6j(0, 0.5, 0.5,
932 : // 0.5, 0, 1))
933 0 : if (s >= std::abs(state.j[0] - state.l[0]) && s <= state.j[0] + state.l[0] &&
934 0 : s >= std::abs(s - kappa) && s <= s + kappa &&
935 0 : state.j[1] >= std::abs(state.j[0] - kappa) &&
936 0 : state.j[1] <= state.j[0] + kappa && state.j[1] >= std::abs(s - state.l[0]) &&
937 0 : state.j[1] <= s + state.l[0]) {
938 0 : cache.second = std::pow(-1, state.l[0] + s + state.j[0] + kappa) *
939 0 : std::sqrt((2 * state.j[0] + 1) * (2 * state.j[1] + 1)) *
940 0 : WignerSymbols::wigner6j(s, state.j[0], state.l[0], state.j[1], s, kappa);
941 : } else {
942 0 : cache.second = 0;
943 : }
944 :
945 0 : stmt.bind(1, s);
946 0 : stmt.bind(2, kappa);
947 0 : stmt.bind(3, state.l[0]);
948 0 : stmt.bind(4, state.j[0]);
949 0 : stmt.bind(5, state.l[1]);
950 0 : stmt.bind(6, state.j[1]);
951 0 : stmt.bind(7, cache.second);
952 0 : stmt.step();
953 0 : stmt.reset();
954 : }
955 : }
956 : }
957 :
958 0 : if (calcElectricMultipole) {
959 0 : stmt.set("insert or ignore into cache_reduced_multipole (k, l1, l2, value) values (?1, ?2, "
960 0 : "?3, ?4);");
961 :
962 0 : stmt.prepare();
963 0 : for (auto &cache : cache_reduced_multipole[kappa]) {
964 0 : if (cache.second == std::numeric_limits<double>::max()) {
965 0 : auto state = cache.first;
966 :
967 0 : cache.second = std::pow(-1, state.l[0]) *
968 0 : std::sqrt((2 * state.l[0] + 1) * (2 * state.l[1] + 1)) *
969 0 : WignerSymbols::wigner3j(state.l[0], kappa, state.l[1], 0, 0, 0);
970 :
971 0 : if (cache.second > 1e9) {
972 : std::cout << "Warning: Error in the calculation of the Wigner-3j-symbol "
973 0 : "detected and resolved."
974 0 : << std::endl;
975 0 : cache.second = 0;
976 : }
977 :
978 0 : stmt.bind(1, kappa);
979 0 : stmt.bind(2, state.l[0]);
980 0 : stmt.bind(3, state.l[1]);
981 0 : stmt.bind(4, cache.second);
982 0 : stmt.step();
983 0 : stmt.reset();
984 : }
985 : }
986 : }
987 :
988 0 : stmt.exec("commit transaction;");
989 0 : }
|
