Line data Source code
1 : /*
2 : * Copyright (c) 2018 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 "MatrixElementCache.hpp"
21 :
22 : #include "Constants.hpp"
23 : #include "QuantumDefect.hpp"
24 : #include "SQLite.hpp"
25 : #include "filesystem.hpp"
26 : #include "utils.hpp"
27 : #include "version.hpp"
28 :
29 : #include <boost/tokenizer.hpp>
30 :
31 : #include <cctype>
32 : #include <iostream>
33 : #include <limits>
34 : #include <memory>
35 : #include <sstream>
36 : #include <stdexcept>
37 : #include <string>
38 :
39 87158 : bool selectionRulesMomentumNew(StateOne const &state1, StateOne const &state2, int q) {
40 87158 : bool validL = state1.getL() == state2.getL();
41 87158 : bool validJ = std::fabs(state1.getJ() - state2.getJ()) <= 1;
42 87158 : bool validM = state1.getM() == state2.getM() + q;
43 87158 : bool validQ = abs(q) <= 1;
44 87158 : return validL && validJ && validM && validQ;
45 : }
46 :
47 13482 : bool selectionRulesMomentumNew(StateOne const &state1, StateOne const &state2) {
48 13482 : bool validL = state1.getL() == state2.getL();
49 13482 : bool validJ = std::fabs(state1.getJ() - state2.getJ()) <= 1;
50 13482 : bool validM = (std::fabs(state1.getM() - state2.getM()) <= 1);
51 13482 : return validL && validJ && validM;
52 : }
53 :
54 251558 : bool selectionRulesMultipoleNew(StateOne const &state1, StateOne const &state2, int kappa, int q) {
55 439906 : bool validL = (abs(state1.getL() - state2.getL()) <= kappa) &&
56 188348 : (kappa % 2 == abs(state1.getL() - state2.getL()) % 2);
57 435979 : bool validJ = (std::fabs(state1.getJ() - state2.getJ()) <= kappa) &&
58 184421 : (state1.getJ() + state2.getJ() >= kappa);
59 251558 : bool validM = state1.getM() == state2.getM() + q;
60 251558 : bool validQ = abs(q) <= kappa;
61 : bool noZero =
62 258868 : !(kappa == 2 && state1.getJ() == state2.getJ() && state2.getJ() == 1.5 &&
63 7310 : state1.getM() == -state2.getM() && std::fabs(state1.getM() - state2.getM()) == 1);
64 251558 : return validL && validJ && validM && validQ && noZero;
65 : }
66 :
67 1813870 : bool selectionRulesMultipoleNew(StateOne const &state1, StateOne const &state2, int kappa) {
68 3498658 : bool validL = (abs(state1.getL() - state2.getL()) <= kappa) &&
69 1684788 : (kappa % 2 == abs(state1.getL() - state2.getL()) % 2);
70 3493800 : bool validJ = (std::fabs(state1.getJ() - state2.getJ()) <= kappa) &&
71 1679930 : (state1.getJ() + state2.getJ() >= kappa);
72 1813870 : bool validM = (std::fabs(state1.getM() - state2.getM()) <= kappa);
73 : bool noZero =
74 1852148 : !(kappa == 2 && state1.getJ() == state2.getJ() && state2.getJ() == 1.5 &&
75 38278 : state1.getM() == -state2.getM() && std::fabs(state1.getM() - state2.getM()) == 1);
76 1813870 : return validL && validJ && validM && noZero;
77 : }
78 :
79 : ////////////////////////////////////////////////////////////////////
80 : /// Constructors ///////////////////////////////////////////////////
81 : ////////////////////////////////////////////////////////////////////
82 :
83 33 : MatrixElementCache::MatrixElementCache()
84 99 : : defectdbname(""), dbname(""), pid_which_created_db(utils::get_pid()) {}
85 :
86 25 : MatrixElementCache::MatrixElementCache(std::string const &cachedir)
87 : : defectdbname(""),
88 50 : dbname((fs::absolute(cachedir) / ("cache_elements_" + version::cache() + ".db")).string()),
89 50 : db(new sqlite::handle(dbname)), stmt(new sqlite::statement(*db)),
90 125 : pid_which_created_db(utils::get_pid()) {
91 :
92 : // Speed up database access
93 25 : stmt->exec("PRAGMA synchronous = OFF"); // do not wait on write, hand off to OS and carry on
94 50 : stmt->exec(
95 50 : "PRAGMA journal_mode = MEMORY"); // keep rollback journal in memory during transaction
96 :
97 : // Create cache tables (reduced_momentum_s and reduced_momentum_l need not to be cached since
98 : // they are trivial)
99 50 : stmt->exec("create table if not exists cache_radial ("
100 : "method int, species text, k integer, n1 integer, l1 integer, j1 double,"
101 : "n2 integer, l2 integer, j2 double, value double, primary key (method, species, k, "
102 50 : "n1, l1, j1, n2, l2, j2)) without rowid;");
103 :
104 50 : stmt->exec(
105 : "create table if not exists cache_angular ("
106 : "k integer, j1 double, m1 double,"
107 50 : "j2 double, m2 double, value double, primary key (k, j1, m1, j2, m2)) without rowid;");
108 :
109 50 : stmt->exec(
110 : "create table if not exists cache_reduced_commutes_s ("
111 : "s double, k integer, l1 integer, j1 double,"
112 50 : "l2 integer, j2 double, value double, primary key (s, k, l1, j1, l2, j2)) without rowid;");
113 :
114 50 : stmt->exec(
115 : "create table if not exists cache_reduced_commutes_l ("
116 : "s double, k integer, l1 integer, j1 double,"
117 50 : "l2 integer, j2 double, value double, primary key (s, k, l1, j1, l2, j2)) without rowid;");
118 :
119 50 : stmt->exec("create table if not exists cache_reduced_multipole ("
120 : "k integer, l1 integer,"
121 50 : "l2 integer, value double, primary key (k, l1, l2)) without rowid;");
122 25 : }
123 :
124 1 : void MatrixElementCache::setDefectDB(std::string const &path) {
125 : // Set path to user defined database for quantum defects and model potential parameters
126 1 : defectdbname = path;
127 :
128 : // Prevent using the cache file to avoid inconsistencies through the user defined database
129 1 : dbname = "";
130 1 : }
131 :
132 12 : void MatrixElementCache::setMethod(method_t const &m) {
133 : // Set method for calculating the radial matrix element
134 12 : method = m;
135 12 : }
136 :
137 1 : void MatrixElementCache::loadElectricDipoleDB(std::string const &path, std::string const &species) {
138 1 : int kappa_angular = 1;
139 1 : float s = 0.5;
140 1 : if (std::isdigit(species.back()) != 0) {
141 0 : s = ((species.back() - '0') - 1) / 2.;
142 : }
143 :
144 : // Get the radial matrix elements from the CSV file
145 5 : boost::escaped_list_separator<char> els("", ";", "\"\'");
146 2 : std::vector<std::pair<CacheKey_cache_radial, double>> radial_raw_data;
147 2 : std::ifstream ifs(path.c_str());
148 2 : std::string line;
149 1 : bool firstline = true;
150 73 : while (std::getline(ifs, line)) {
151 144 : boost::tokenizer<boost::escaped_list_separator<char>> tk(line, els);
152 144 : std::vector<std::string> tk_vec(tk.begin(), tk.end());
153 :
154 : try {
155 72 : int n1 = std::stoi(tk_vec[0]);
156 71 : int l1 = std::stoi(tk_vec[1]);
157 71 : float j1 = std::stof(tk_vec[2]);
158 71 : int n2 = std::stoi(tk_vec[3]);
159 71 : int l2 = std::stoi(tk_vec[4]);
160 71 : float j2 = std::stof(tk_vec[5]);
161 71 : double val = std::stod(tk_vec[6]);
162 :
163 : // Convert the radial matrix element to um
164 71 : val *= au2um;
165 :
166 : // Store the radial matrix element
167 : auto key1 =
168 142 : CacheKey_cache_radial(method, species, kappa_angular, n1, n2, l1, l2, j1, j2);
169 71 : radial_raw_data.emplace_back(key1, val);
170 :
171 : // Store which expressions are needed for converting the radial matrix element to a
172 : // reduced radial matrix element
173 71 : auto key3 = CacheKey_cache_reduced_commutes(s, kappa_angular, l1, l2, j1, j2);
174 71 : cache_reduced_commutes_s_missing.insert(key3);
175 71 : auto key4 = CacheKey_cache_reduced_multipole(kappa_angular, l1, l2);
176 71 : cache_reduced_multipole_missing.insert(key4);
177 :
178 2 : } catch (std::invalid_argument &e) {
179 1 : if (!firstline) { // Skip header
180 : std::cerr << "WARNING: During loading the electric dipole database, the error '"
181 0 : << e.what() << "' occured." << std::endl;
182 : }
183 : }
184 :
185 72 : firstline = false;
186 : }
187 1 : ifs.close();
188 :
189 : // Precaculate the expressions
190 1 : this->update();
191 :
192 : // Calculate the reduced radial matrix element and save it to the in-memory cache (it's not
193 : // saved to the sqlite database)
194 72 : for (auto const &pair : radial_raw_data) {
195 71 : auto key1 = pair.first;
196 142 : auto key3 = CacheKey_cache_reduced_commutes(s, kappa_angular, key1.l[0], key1.l[1],
197 71 : key1.j[0], key1.j[1]);
198 71 : auto key4 = CacheKey_cache_reduced_multipole(kappa_angular, key1.l[0], key1.l[1]);
199 71 : double val = pair.second /
200 71 : (key3.sgn * cache_reduced_commutes_s[key3] * key4.sgn * cache_reduced_multipole[key4]);
201 :
202 71 : cache_radial[key1] = val;
203 : }
204 1 : }
205 :
206 : ////////////////////////////////////////////////////////////////////
207 : /// Keys ///////////////////////////////////////////////////////////
208 : ////////////////////////////////////////////////////////////////////
209 :
210 554189 : MatrixElementCache::CacheKey_cache_radial::CacheKey_cache_radial(method_t method,
211 : std::string species, int kappa,
212 : int n1, int n2, int l1, int l2,
213 554189 : float j1, float j2)
214 554189 : : species(std::move(species)), method(method), kappa(kappa) {
215 554189 : if ((n1 < n2) || ((n1 == n2) && ((l1 < l2) || ((l1 == l2) && (j1 <= j2))))) {
216 231070 : n = {{n1, n2}};
217 231070 : l = {{l1, l2}};
218 231070 : j = {{j1, j2}};
219 : } else {
220 323119 : n = {{n2, n1}};
221 323119 : l = {{l2, l1}};
222 323119 : j = {{j2, j1}};
223 : }
224 554189 : }
225 :
226 553878 : MatrixElementCache::CacheKey_cache_angular::CacheKey_cache_angular(int kappa, float j1, float j2,
227 553878 : float m1, float m2)
228 553878 : : kappa(kappa) {
229 553878 : if ((j1 < j2) || ((j1 == j2) && (m1 <= m2))) {
230 328946 : j = {{j1, j2}};
231 328946 : m = {{m1, m2}};
232 328946 : sgn = 1;
233 : } else {
234 224932 : j = {{j2, j1}};
235 224932 : m = {{m2, m1}};
236 224932 : sgn = std::pow(-1, int(j1 - m1 + j2 - m2));
237 : }
238 553878 : }
239 :
240 562840 : MatrixElementCache::CacheKey_cache_reduced_commutes::CacheKey_cache_reduced_commutes(
241 562840 : float s, int kappa, int l1, int l2, float j1, float j2)
242 562840 : : s(s), kappa(kappa) {
243 562840 : if ((l1 < l2) || ((l1 == l2) && (j1 <= j2))) {
244 299616 : l = {{l1, l2}};
245 299616 : j = {{j1, j2}};
246 299616 : sgn = 1;
247 : } else {
248 263224 : l = {{l2, l1}};
249 263224 : j = {{j2, j1}};
250 263224 : sgn = std::pow(-1, int(l1 + j1 + l2 + j2 + 2 * s)); // TODO is this formula always correct?
251 : }
252 562840 : }
253 :
254 545200 : MatrixElementCache::CacheKey_cache_reduced_multipole::CacheKey_cache_reduced_multipole(int kappa,
255 : int l1,
256 545200 : int l2)
257 545200 : : kappa(kappa) {
258 545200 : if (l1 <= l2) {
259 292639 : l = {{l1, l2}};
260 292639 : sgn = 1;
261 : } else {
262 252561 : l = {{l2, l1}};
263 252561 : sgn = std::pow(-1, kappa);
264 : }
265 545200 : }
266 :
267 570 : MatrixElementCache::CacheKey_cache_radial::CacheKey_cache_radial() =
268 : default; // TODO the default constructors seem to be needed for serialization, can one somehow
269 : // circumvent the need of default constructors?
270 :
271 920 : MatrixElementCache::CacheKey_cache_angular::CacheKey_cache_angular() =
272 : default; // TODO the default constructors seem to be needed for serialization, can one somehow
273 : // circumvent the need of default constructors?
274 :
275 120 : MatrixElementCache::CacheKey_cache_reduced_commutes::CacheKey_cache_reduced_commutes() =
276 : default; // TODO the default constructors seem to be needed for serialization, can one somehow
277 : // circumvent the need of default constructors?
278 :
279 40 : MatrixElementCache::CacheKey_cache_reduced_multipole::CacheKey_cache_reduced_multipole() =
280 : default; // TODO the default constructors seem to be needed for serialization, can one somehow
281 : // circumvent the need of default constructors?
282 :
283 551479 : bool MatrixElementCache::CacheKey_cache_radial::operator==(const CacheKey_cache_radial &rhs) const {
284 551479 : return (method == rhs.method) && (species == rhs.species) && (kappa == rhs.kappa) &&
285 1102958 : (n == rhs.n) && (l == rhs.l) && (j == rhs.j);
286 : }
287 :
288 551678 : bool MatrixElementCache::CacheKey_cache_angular::operator==(
289 : const CacheKey_cache_angular &rhs) const {
290 551678 : return (kappa == rhs.kappa) && (j == rhs.j) && (m == rhs.m);
291 : }
292 :
293 562249 : bool MatrixElementCache::CacheKey_cache_reduced_commutes::operator==(
294 : const CacheKey_cache_reduced_commutes &rhs) const {
295 562249 : return (s == rhs.s) && (kappa == rhs.kappa) && (l == rhs.l) && (j == rhs.j);
296 : }
297 :
298 546506 : bool MatrixElementCache::CacheKey_cache_reduced_multipole::operator==(
299 : const CacheKey_cache_reduced_multipole &rhs) const {
300 546506 : return (kappa == rhs.kappa) && (l == rhs.l);
301 : }
302 :
303 : ////////////////////////////////////////////////////////////////////
304 : /// Hasher /////////////////////////////////////////////////////////
305 : ////////////////////////////////////////////////////////////////////
306 :
307 : std::size_t
308 577645 : MatrixElementCache::CacheKeyHasher_cache_radial::operator()(const CacheKey_cache_radial &c) const {
309 577645 : size_t seed = 0;
310 577645 : utils::hash_combine(seed, c.method);
311 577645 : utils::hash_combine(seed, c.species);
312 577645 : utils::hash_combine(seed, c.kappa);
313 577645 : utils::hash_combine(seed, c.n);
314 577645 : utils::hash_combine(seed, c.l);
315 577645 : utils::hash_combine(seed, c.j);
316 577645 : return seed;
317 : }
318 :
319 577603 : std::size_t MatrixElementCache::CacheKeyHasher_cache_angular::operator()(
320 : const CacheKey_cache_angular &c) const {
321 577603 : size_t seed = 0;
322 577603 : utils::hash_combine(seed, c.kappa);
323 577603 : utils::hash_combine(seed, c.j);
324 577603 : utils::hash_combine(seed, c.m);
325 577603 : return seed;
326 : }
327 :
328 586019 : std::size_t MatrixElementCache::CacheKeyHasher_cache_reduced_commutes::operator()(
329 : const CacheKey_cache_reduced_commutes &c) const {
330 586019 : size_t seed = 0;
331 586019 : utils::hash_combine(seed, c.s);
332 586019 : utils::hash_combine(seed, c.kappa);
333 586019 : utils::hash_combine(seed, c.l);
334 586019 : utils::hash_combine(seed, c.j);
335 586019 : return seed;
336 : }
337 :
338 560774 : std::size_t MatrixElementCache::CacheKeyHasher_cache_reduced_multipole::operator()(
339 : const CacheKey_cache_reduced_multipole &c) const {
340 560774 : size_t seed = 0;
341 560774 : utils::hash_combine(seed, c.kappa);
342 560774 : utils::hash_combine(seed, c.l);
343 560774 : return seed;
344 : }
345 :
346 : ////////////////////////////////////////////////////////////////////
347 : /// Get matrix elements ////////////////////////////////////////////
348 : ////////////////////////////////////////////////////////////////////
349 :
350 410070 : double MatrixElementCache::getElectricDipole(StateOne const &state_row, StateOne const &state_col) {
351 410070 : return getElectricMultipole(state_row, state_col, 1, 1);
352 : }
353 :
354 85098 : double MatrixElementCache::getElectricMultipole(StateOne const &state_row,
355 : StateOne const &state_col, int k) {
356 85098 : return getElectricMultipole(state_row, state_col, k, k);
357 : }
358 :
359 9400 : double MatrixElementCache::getDiamagnetism(StateOne const &state_row, StateOne const &state_col,
360 : int k) {
361 9400 : return 2. / 3. * elementary_charge * getElectricMultipole(state_row, state_col, 2, k);
362 : }
363 :
364 4410 : double MatrixElementCache::getMagneticDipole(StateOne const &state_row, StateOne const &state_col) {
365 4410 : if (state_row.getSpecies() != state_col.getSpecies()) {
366 0 : throw std::runtime_error("The species must be the same for the final and initial state.");
367 : }
368 :
369 4410 : const std::string &species = state_row.getSpecies();
370 4410 : const float &s = state_row.getS();
371 :
372 : // Search cache for constituents of the matrix element
373 13230 : auto key1 = CacheKey_cache_radial(method, species, 0, state_row.getN(), state_col.getN(),
374 17640 : state_row.getL(), state_col.getL(), state_row.getJ(),
375 4410 : state_col.getJ());
376 4410 : auto iter1 = cache_radial.find(key1);
377 4410 : if (iter1 == cache_radial.end()) {
378 12 : cache_radial_missing.insert(key1);
379 : }
380 :
381 17640 : auto key2 = CacheKey_cache_angular(1, state_row.getJ(), state_col.getJ(), state_row.getM(),
382 4410 : state_col.getM());
383 4410 : auto iter2 = cache_angular.find(key2);
384 4410 : if (iter2 == cache_angular.end()) {
385 15 : cache_angular_missing.insert(key2);
386 : }
387 :
388 13230 : auto key3 = CacheKey_cache_reduced_commutes(s, 1, state_row.getL(), state_col.getL(),
389 4410 : state_row.getJ(), state_col.getJ());
390 4410 : auto iter3 = cache_reduced_commutes_s.find(key3);
391 4410 : if (iter3 == cache_reduced_commutes_s.end()) {
392 12 : cache_reduced_commutes_s_missing.insert(key3);
393 : }
394 :
395 13230 : auto key4 = CacheKey_cache_reduced_commutes(s, 1, state_row.getL(), state_col.getL(),
396 4410 : state_row.getJ(), state_col.getJ());
397 4410 : auto iter4 = cache_reduced_commutes_l.find(key4);
398 4410 : if (iter4 == cache_reduced_commutes_l.end()) {
399 12 : cache_reduced_commutes_l_missing.insert(key4);
400 : }
401 :
402 : // Update cache by calculate missing constituents
403 4410 : if (this->update() != 0) {
404 17 : iter1 = cache_radial.find(key1);
405 17 : iter2 = cache_angular.find(key2);
406 17 : iter3 = cache_reduced_commutes_s.find(key3);
407 17 : iter4 = cache_reduced_commutes_l.find(key4);
408 : }
409 :
410 4410 : return -bohr_magneton * iter1->second * key2.sgn * iter2->second *
411 4410 : (gL * key3.sgn * iter3->second *
412 4410 : std::sqrt(state_row.getL() * (state_row.getL() + 1) * (2 * state_row.getL() + 1)) +
413 8820 : gS * key4.sgn * iter4->second * std::sqrt(s * (s + 1) * (2 * s + 1)));
414 : }
415 :
416 504568 : double MatrixElementCache::getElectricMultipole(StateOne const &state_row,
417 : StateOne const &state_col, int kappa_radial,
418 : int kappa_angular) {
419 504568 : if (state_row.getSpecies() != state_col.getSpecies()) {
420 0 : throw std::runtime_error("The species must be the same for the final and initial state.");
421 : }
422 :
423 504568 : const std::string &species = state_row.getSpecies();
424 504568 : const float &s = state_row.getS();
425 :
426 : // Search cache for constituents of the matrix element
427 1009136 : auto key1 = CacheKey_cache_radial(method, species, kappa_radial, state_row.getN(),
428 2018272 : state_col.getN(), state_row.getL(), state_col.getL(),
429 504568 : state_row.getJ(), state_col.getJ());
430 504568 : auto iter1 = cache_radial.find(key1);
431 504568 : if (iter1 == cache_radial.end()) {
432 97 : cache_radial_missing.insert(key1);
433 : }
434 :
435 1513704 : auto key2 = CacheKey_cache_angular(kappa_angular, state_row.getJ(), state_col.getJ(),
436 504568 : state_row.getM(), state_col.getM());
437 504568 : auto iter2 = cache_angular.find(key2);
438 504568 : if (iter2 == cache_angular.end()) {
439 48 : cache_angular_missing.insert(key2);
440 : }
441 :
442 : auto key3 = CacheKey_cache_reduced_commutes(
443 504568 : s, kappa_angular, state_row.getL(), state_col.getL(), state_row.getJ(), state_col.getJ());
444 504568 : auto iter3 = cache_reduced_commutes_s.find(key3);
445 504568 : if (iter3 == cache_reduced_commutes_s.end()) {
446 38 : cache_reduced_commutes_s_missing.insert(key3);
447 : }
448 :
449 504568 : auto key4 = CacheKey_cache_reduced_multipole(kappa_angular, state_row.getL(), state_col.getL());
450 504568 : auto iter4 = cache_reduced_multipole.find(key4);
451 504568 : if (iter4 == cache_reduced_multipole.end()) {
452 34 : cache_reduced_multipole_missing.insert(key4);
453 : }
454 :
455 : // Update cache by calculate missing constituents
456 504568 : if (this->update() != 0) {
457 118 : iter1 = cache_radial.find(key1);
458 118 : iter2 = cache_angular.find(key2);
459 118 : iter3 = cache_reduced_commutes_s.find(key3);
460 118 : iter4 = cache_reduced_multipole.find(key4);
461 : }
462 :
463 504568 : return elementary_charge * iter1->second * key2.sgn * iter2->second * key3.sgn * iter3->second *
464 1009136 : key4.sgn * iter4->second;
465 : }
466 :
467 240 : double MatrixElementCache::getRadial(StateOne const &state_row, StateOne const &state_col,
468 : int kappa) {
469 240 : if (state_row.getSpecies() != state_col.getSpecies()) {
470 0 : throw std::runtime_error("The species must be the same for the final and initial state.");
471 : }
472 :
473 240 : const std::string &species = state_row.getSpecies();
474 :
475 : // Search cache for constituents of the matrix element
476 720 : auto key1 = CacheKey_cache_radial(method, species, kappa, state_row.getN(), state_col.getN(),
477 960 : state_row.getL(), state_col.getL(), state_row.getJ(),
478 240 : state_col.getJ());
479 240 : auto iter1 = cache_radial.find(key1);
480 240 : if (iter1 == cache_radial.end()) {
481 149 : cache_radial_missing.insert(key1);
482 : }
483 :
484 : // Update cache by calculate missing constituents
485 240 : if (this->update() != 0) {
486 149 : if (iter1 == cache_radial.end()) {
487 149 : iter1 = cache_radial.find(key1);
488 : }
489 : }
490 :
491 480 : return iter1->second;
492 : }
493 :
494 158169 : const std::string &MatrixElementCache::getDefectDB() const { return defectdbname; }
495 :
496 : ////////////////////////////////////////////////////////////////////
497 : /// Precalculate matrix elements ///////////////////////////////////
498 : ////////////////////////////////////////////////////////////////////
499 :
500 : // TODO improve methods for precalculation, for every method getSomething create the method
501 : // precalculateSomething
502 :
503 287 : void MatrixElementCache::precalculateMultipole(const std::vector<StateOne> &basis_one, int k) {
504 287 : int q = std::numeric_limits<int>::max();
505 287 : precalculate(basis_one, k, q, k, true, false, false);
506 287 : }
507 :
508 0 : void MatrixElementCache::precalculateRadial(const std::vector<StateOne> &basis_one, int k) {
509 0 : int q = std::numeric_limits<int>::max();
510 0 : precalculate(basis_one, k, q, k, false, false, true);
511 0 : }
512 :
513 70 : void MatrixElementCache::precalculateElectricMomentum(const std::vector<StateOne> &basis_one,
514 : int q) {
515 70 : precalculate(basis_one, 1, q, 1, true, false, false);
516 70 : }
517 :
518 90 : void MatrixElementCache::precalculateMagneticMomentum(const std::vector<StateOne> &basis_one,
519 : int q) {
520 90 : precalculate(basis_one, 1, q, 0, false, true, false);
521 90 : }
522 :
523 170 : void MatrixElementCache::precalculateDiamagnetism(const std::vector<StateOne> &basis_one, int k,
524 : int q) {
525 170 : precalculate(basis_one, k, q, 2, true, false, false);
526 170 : }
527 :
528 : ////////////////////////////////////////////////////////////////////
529 : /// Utility methods ////////////////////////////////////////////////
530 : ////////////////////////////////////////////////////////////////////
531 :
532 2464 : double MatrixElementCache::calcRadialElement(const QuantumDefect &qd1, int power,
533 : const QuantumDefect &qd2) {
534 2464 : if (method == NUMEROV) {
535 2454 : return std::pow(au2um, power) * IntegrateRadialElement<Numerov>(qd1, power, qd2);
536 : }
537 10 : if (method == WHITTAKER) {
538 10 : return std::pow(au2um, power) * IntegrateRadialElement<Whittaker>(qd1, power, qd2);
539 : }
540 : std::string msg("You have to provide all radial matrix elements on your own because you have "
541 0 : "deactivated the calculation of missing radial matrix elements!");
542 0 : std::cout << msg << std::endl;
543 0 : throw std::runtime_error(msg);
544 : }
545 :
546 617 : void MatrixElementCache::precalculate(const std::vector<StateOne> &basis_one, int kappa_angular,
547 : int q, int kappa_radial, bool calcElectricMultipole,
548 : bool calcMagneticMomentum, bool calcRadial) {
549 :
550 1234 : std::string species;
551 617 : float s = std::numeric_limits<float>::max();
552 :
553 : // --- Determine elements ---
554 :
555 16880 : for (size_t idx_col = 0; idx_col < basis_one.size(); ++idx_col) {
556 16263 : const auto &state_col = basis_one[idx_col];
557 16263 : if (state_col.isArtificial()) {
558 0 : continue;
559 : }
560 :
561 16263 : if (species.empty()) {
562 617 : species = state_col.getSpecies();
563 617 : s = state_col.getS();
564 : }
565 :
566 485076 : for (size_t idx_row = 0; idx_row <= idx_col; ++idx_row) {
567 468813 : const auto &state_row = basis_one[idx_row];
568 468813 : if (state_row.isArtificial()) {
569 0 : continue;
570 : }
571 :
572 468813 : if (q != std::numeric_limits<int>::max() && state_row.getM() - state_col.getM() != q) {
573 294420 : continue;
574 : }
575 :
576 335304 : if ((calcElectricMultipole &&
577 174393 : selectionRulesMultipoleNew(state_row, state_col, kappa_angular)) ||
578 348786 : (calcMagneticMomentum && selectionRulesMomentumNew(state_row, state_col)) ||
579 : calcRadial) {
580 44900 : if (calcElectricMultipole || calcMagneticMomentum || calcRadial) {
581 : auto key = CacheKey_cache_radial(
582 134700 : method, species, kappa_radial, state_row.getN(), state_col.getN(),
583 89800 : state_row.getL(), state_col.getL(), state_row.getJ(), state_col.getJ());
584 44900 : if (cache_radial.find(key) == cache_radial.end()) {
585 19403 : cache_radial_missing.insert(key);
586 : }
587 : }
588 :
589 44900 : if (calcElectricMultipole || calcMagneticMomentum) {
590 : auto key =
591 134700 : CacheKey_cache_angular(kappa_angular, state_row.getJ(), state_col.getJ(),
592 44900 : state_row.getM(), state_col.getM());
593 44900 : if (cache_angular.find(key) == cache_angular.end()) {
594 20272 : cache_angular_missing.insert(key);
595 : }
596 : }
597 :
598 44900 : if (calcElectricMultipole || calcMagneticMomentum) {
599 89800 : auto key = CacheKey_cache_reduced_commutes(s, kappa_angular, state_row.getL(),
600 134700 : state_col.getL(), state_row.getJ(),
601 44900 : state_col.getJ());
602 44900 : if (cache_reduced_commutes_s.find(key) == cache_reduced_commutes_s.end()) {
603 18651 : cache_reduced_commutes_s_missing.insert(key);
604 : }
605 : }
606 :
607 44900 : if (calcMagneticMomentum) {
608 8820 : auto key = CacheKey_cache_reduced_commutes(s, kappa_angular, state_row.getL(),
609 13230 : state_col.getL(), state_row.getJ(),
610 4410 : state_col.getJ());
611 4410 : if (cache_reduced_commutes_l.find(key) == cache_reduced_commutes_l.end()) {
612 3428 : cache_reduced_commutes_l_missing.insert(key);
613 : }
614 : }
615 :
616 44900 : if (calcElectricMultipole) {
617 80980 : auto key = CacheKey_cache_reduced_multipole(kappa_angular, state_row.getL(),
618 40490 : state_col.getL());
619 40490 : if (cache_reduced_multipole.find(key) == cache_reduced_multipole.end()) {
620 15169 : cache_reduced_multipole_missing.insert(key);
621 : }
622 : }
623 : }
624 : }
625 : }
626 617 : }
627 :
628 509219 : int MatrixElementCache::update() {
629 :
630 : // --- Return if the cache is already up-to-date ---
631 1527112 : if (cache_radial_missing.empty() && cache_angular_missing.empty() &&
632 2036046 : cache_reduced_commutes_s_missing.empty() && cache_reduced_commutes_l_missing.empty() &&
633 508934 : cache_reduced_multipole_missing.empty()) {
634 508934 : return 0;
635 : }
636 :
637 : // --- Load from database ---
638 285 : if (!dbname.empty()) {
639 :
640 : // Reopen the connection to the database if it was opend by a different process (without
641 : // doing this, there can be problems with Python multiprocessing)
642 230 : if (pid_which_created_db != static_cast<long>(utils::get_pid())) {
643 0 : db = std::make_unique<sqlite::handle>(dbname);
644 0 : stmt = std::make_unique<sqlite::statement>(*db);
645 0 : pid_which_created_db = utils::get_pid();
646 :
647 : // Speed up database access
648 0 : stmt->exec(
649 0 : "PRAGMA synchronous = OFF"); // do not wait on write, hand off to OS and carry on
650 0 : stmt->exec("PRAGMA journal_mode = MEMORY"); // keep rollback journal in memory during
651 : // transaction
652 : }
653 :
654 230 : if (!cache_radial_missing.empty()) {
655 420 : stmt->set("select value from cache_radial where `method` = ?1 and `species` = ?2 and "
656 : "`k` = ?3 and `n1` = ?4 and `l1` = ?5 and `j1` = ?6 and `n2` = ?7 and `l2` = "
657 420 : "?8 and `j2` = ?9;");
658 210 : stmt->prepare();
659 :
660 1782 : for (auto cached = cache_radial_missing.begin();
661 1782 : cached != cache_radial_missing.end();) {
662 1572 : stmt->bind(1, cached->method);
663 1572 : stmt->bind(2, cached->species);
664 1572 : stmt->bind(3, cached->kappa);
665 1572 : stmt->bind(4, cached->n[0]);
666 1572 : stmt->bind(5, cached->l[0]);
667 1572 : stmt->bind(6, cached->j[0]);
668 1572 : stmt->bind(7, cached->n[1]);
669 1572 : stmt->bind(8, cached->l[1]);
670 1572 : stmt->bind(9, cached->j[1]);
671 1572 : if (stmt->step()) {
672 477 : cache_radial.insert({*cached, stmt->get<double>(0)});
673 477 : cached = cache_radial_missing.erase(cached);
674 : } else {
675 1095 : ++cached;
676 : }
677 1572 : stmt->reset();
678 : }
679 : }
680 :
681 230 : if (!cache_angular_missing.empty()) {
682 100 : stmt->set("select value from cache_angular where `k` = ?1 and `j1` = ?2 and `m1` = ?3 "
683 100 : "and `j2` = ?4 and `m2` = ?5;");
684 50 : stmt->prepare();
685 :
686 1284 : for (auto cached = cache_angular_missing.begin();
687 1284 : cached != cache_angular_missing.end();) {
688 1234 : stmt->bind(1, cached->kappa);
689 1234 : stmt->bind(2, cached->j[0]);
690 1234 : stmt->bind(3, cached->m[0]);
691 1234 : stmt->bind(4, cached->j[1]);
692 1234 : stmt->bind(5, cached->m[1]);
693 1234 : if (stmt->step()) {
694 404 : cache_angular.insert({*cached, stmt->get<double>(0)});
695 404 : cached = cache_angular_missing.erase(cached);
696 : } else {
697 830 : ++cached;
698 : }
699 1234 : stmt->reset();
700 : }
701 : }
702 :
703 230 : if (!cache_reduced_commutes_s_missing.empty()) {
704 68 : stmt->set("select value from cache_reduced_commutes_s where `s` = ?1 and `k` = ?2 and "
705 68 : "`l1` = ?3 and `j1` = ?4 and `l2` = ?5 and `j2` = ?6;");
706 34 : stmt->prepare();
707 :
708 293 : for (auto cached = cache_reduced_commutes_s_missing.begin();
709 293 : cached != cache_reduced_commutes_s_missing.end();) {
710 259 : stmt->bind(1, cached->s);
711 259 : stmt->bind(2, cached->kappa);
712 259 : stmt->bind(3, cached->l[0]);
713 259 : stmt->bind(4, cached->j[0]);
714 259 : stmt->bind(5, cached->l[1]);
715 259 : stmt->bind(6, cached->j[1]);
716 259 : if (stmt->step()) {
717 60 : cache_reduced_commutes_s.insert({*cached, stmt->get<double>(0)});
718 60 : cached = cache_reduced_commutes_s_missing.erase(cached);
719 : } else {
720 199 : ++cached;
721 : }
722 259 : stmt->reset();
723 : }
724 : }
725 :
726 230 : if (!cache_reduced_commutes_l_missing.empty()) {
727 20 : stmt->set("select value from cache_reduced_commutes_l where `s` = ?1 and `k` = ?2 and "
728 20 : "`l1` = ?3 and `j1` = ?4 and `l2` = ?5 and `j2` = ?6;");
729 10 : stmt->prepare();
730 :
731 79 : for (auto cached = cache_reduced_commutes_l_missing.begin();
732 79 : cached != cache_reduced_commutes_l_missing.end();) {
733 69 : stmt->bind(1, cached->s);
734 69 : stmt->bind(2, cached->kappa);
735 69 : stmt->bind(3, cached->l[0]);
736 69 : stmt->bind(4, cached->j[0]);
737 69 : stmt->bind(5, cached->l[1]);
738 69 : stmt->bind(6, cached->j[1]);
739 69 : if (stmt->step()) {
740 13 : cache_reduced_commutes_l.insert({*cached, stmt->get<double>(0)});
741 13 : cached = cache_reduced_commutes_l_missing.erase(cached);
742 : } else {
743 56 : ++cached;
744 : }
745 69 : stmt->reset();
746 : }
747 : }
748 :
749 230 : if (!cache_reduced_multipole_missing.empty()) {
750 62 : stmt->set("select value from cache_reduced_multipole where `k` = ?1 and `l1` = ?2 and "
751 62 : "`l2` = ?3;");
752 31 : stmt->prepare();
753 :
754 129 : for (auto cached = cache_reduced_multipole_missing.begin();
755 129 : cached != cache_reduced_multipole_missing.end();) {
756 98 : stmt->bind(1, cached->kappa);
757 98 : stmt->bind(2, cached->l[0]);
758 98 : stmt->bind(3, cached->l[1]);
759 98 : if (stmt->step()) {
760 20 : cache_reduced_multipole.insert({*cached, stmt->get<double>(0)});
761 20 : cached = cache_reduced_multipole_missing.erase(cached);
762 : } else {
763 78 : ++cached;
764 : }
765 98 : stmt->reset();
766 : }
767 : }
768 : }
769 :
770 : // --- Calculate missing elements and write them to the database --- // TODO parallelize
771 :
772 285 : if (!dbname.empty()) {
773 230 : stmt->exec("begin transaction;");
774 : }
775 :
776 285 : if (!cache_radial_missing.empty()) {
777 251 : if (!dbname.empty()) {
778 400 : stmt->set("insert or ignore into cache_radial (method, species, k, n1, l1, j1, n2, l2, "
779 400 : "j2, value) values (?1, ?2, ?3, ?4, ?5, ?6, ?7, ?8, ?9, ?10);");
780 200 : stmt->prepare();
781 : }
782 :
783 2715 : for (auto &cached : cache_radial_missing) {
784 2464 : double val = 0;
785 4928 : QuantumDefect qd1(cached.species, cached.n[0], cached.l[0], cached.j[0], defectdbname);
786 4928 : QuantumDefect qd2(cached.species, cached.n[1], cached.l[1], cached.j[1], defectdbname);
787 2464 : val = calcRadialElement(qd1, cached.kappa, qd2);
788 :
789 2464 : cache_radial.insert({cached, val});
790 :
791 2464 : if (!dbname.empty()) {
792 1095 : stmt->bind(1, cached.method);
793 1095 : stmt->bind(2, cached.species);
794 1095 : stmt->bind(3, cached.kappa);
795 1095 : stmt->bind(4, cached.n[0]);
796 1095 : stmt->bind(5, cached.l[0]);
797 1095 : stmt->bind(6, cached.j[0]);
798 1095 : stmt->bind(7, cached.n[1]);
799 1095 : stmt->bind(8, cached.l[1]);
800 1095 : stmt->bind(9, cached.j[1]);
801 1095 : stmt->bind(10, val);
802 1095 : stmt->step();
803 1095 : stmt->reset();
804 : }
805 : }
806 :
807 251 : cache_radial_missing.clear();
808 : }
809 :
810 285 : if (!cache_angular_missing.empty()) {
811 67 : if (!dbname.empty()) {
812 90 : stmt->set("insert or ignore into cache_angular (k, j1, m1, j2, m2, value) values (?1, "
813 90 : "?2, ?3, ?4, ?5, ?6);");
814 45 : stmt->prepare();
815 : }
816 :
817 1998 : for (auto &cached : cache_angular_missing) {
818 1931 : float q = cached.m[0] - cached.m[1];
819 :
820 1931 : double val = std::pow(-1, int(cached.j[0] - cached.m[0])) *
821 1931 : WignerSymbols::wigner3j(cached.j[0], cached.kappa, cached.j[1], -cached.m[0], q,
822 1931 : cached.m[1]);
823 1931 : cache_angular.insert({cached, val});
824 :
825 1931 : if (!dbname.empty()) {
826 830 : stmt->bind(1, cached.kappa);
827 830 : stmt->bind(2, cached.j[0]);
828 830 : stmt->bind(3, cached.m[0]);
829 830 : stmt->bind(4, cached.j[1]);
830 830 : stmt->bind(5, cached.m[1]);
831 830 : stmt->bind(6, val);
832 830 : stmt->step();
833 830 : stmt->reset();
834 : }
835 : }
836 :
837 67 : cache_angular_missing.clear();
838 : }
839 :
840 285 : if (!cache_reduced_commutes_s_missing.empty()) {
841 49 : if (!dbname.empty()) {
842 58 : stmt->set(
843 : "insert or ignore into cache_reduced_commutes_s (s, k, l1, j1, l2, j2, value) "
844 58 : "values (?1, ?2, ?3, ?4, ?5, ?6, ?7);");
845 29 : stmt->prepare();
846 : }
847 :
848 615 : for (auto &cached : cache_reduced_commutes_s_missing) {
849 :
850 : // Check triangle conditions (violated e.g. for WignerSymbols::wigner6j(0, 0.5, 0.5,
851 : // 0.5, 0, 1))
852 566 : double val = 0;
853 566 : if (cached.l[0] >= std::abs(cached.j[0] - cached.s) &&
854 566 : cached.l[0] <= cached.j[0] + cached.s &&
855 566 : cached.l[0] >= std::abs(cached.l[1] - cached.kappa) &&
856 1102 : cached.l[0] <= cached.l[1] + cached.kappa &&
857 551 : cached.j[1] >= std::abs(cached.j[0] - cached.kappa) &&
858 1102 : cached.j[1] <= cached.j[0] + cached.kappa &&
859 1683 : cached.j[1] >= std::abs(cached.l[1] - cached.s) &&
860 551 : cached.j[1] <=
861 551 : cached.l[1] + cached.s) { // TODO implement this in the wignerSymbols library
862 551 : val = std::pow(-1, int(cached.l[0] + cached.s + cached.j[1] + cached.kappa)) *
863 551 : std::sqrt((2 * cached.j[0] + 1) * (2 * cached.j[1] + 1)) *
864 551 : WignerSymbols::wigner6j(cached.l[0], cached.j[0], cached.s, cached.j[1],
865 551 : cached.l[1], cached.kappa);
866 : }
867 566 : cache_reduced_commutes_s.insert({cached, val});
868 :
869 566 : if (!dbname.empty()) {
870 199 : stmt->bind(1, cached.s);
871 199 : stmt->bind(2, cached.kappa);
872 199 : stmt->bind(3, cached.l[0]);
873 199 : stmt->bind(4, cached.j[0]);
874 199 : stmt->bind(5, cached.l[1]);
875 199 : stmt->bind(6, cached.j[1]);
876 199 : stmt->bind(7, val);
877 199 : stmt->step();
878 199 : stmt->reset();
879 : }
880 : }
881 :
882 49 : cache_reduced_commutes_s_missing.clear();
883 : }
884 :
885 285 : if (!cache_reduced_commutes_l_missing.empty()) {
886 18 : if (!dbname.empty()) {
887 18 : stmt->set(
888 : "insert or ignore into cache_reduced_commutes_l (s, k, l1, j1, l2, j2, value) "
889 18 : "values (?1, ?2, ?3, ?4, ?5, ?6, ?7);");
890 9 : stmt->prepare();
891 : }
892 :
893 145 : for (auto &cached : cache_reduced_commutes_l_missing) {
894 :
895 : // Check triangle conditions (violated e.g. for WignerSymbols::wigner6j(0, 0.5, 0.5,
896 : // 0.5, 0, 1))
897 127 : double val = 0;
898 127 : if (cached.s >= std::abs(cached.j[0] - cached.l[0]) &&
899 254 : cached.s <= cached.j[0] + cached.l[0] &&
900 127 : cached.s >= std::abs(cached.s - cached.kappa) &&
901 254 : cached.s <= cached.s + cached.kappa &&
902 127 : cached.j[1] >= std::abs(cached.j[0] - cached.kappa) &&
903 254 : cached.j[1] <= cached.j[0] + cached.kappa &&
904 381 : cached.j[1] >= std::abs(cached.s - cached.l[0]) &&
905 127 : cached.j[1] <=
906 127 : cached.s + cached.l[0]) { // TODO implement this in the wignerSymbols library
907 127 : val = std::pow(-1, int(cached.l[0] + cached.s + cached.j[0] + cached.kappa)) *
908 127 : std::sqrt((2 * cached.j[0] + 1) * (2 * cached.j[1] + 1)) *
909 127 : WignerSymbols::wigner6j(cached.s, cached.j[0], cached.l[0], cached.j[1],
910 127 : cached.s, cached.kappa);
911 : }
912 127 : cache_reduced_commutes_l.insert({cached, val});
913 :
914 127 : if (!dbname.empty()) {
915 56 : stmt->bind(1, cached.s);
916 56 : stmt->bind(2, cached.kappa);
917 56 : stmt->bind(3, cached.l[0]);
918 56 : stmt->bind(4, cached.j[0]);
919 56 : stmt->bind(5, cached.l[1]);
920 56 : stmt->bind(6, cached.j[1]);
921 56 : stmt->bind(7, val);
922 56 : stmt->step();
923 56 : stmt->reset();
924 : }
925 : }
926 :
927 18 : cache_reduced_commutes_l_missing.clear();
928 : }
929 :
930 285 : if (!cache_reduced_multipole_missing.empty()) {
931 45 : if (!dbname.empty()) {
932 52 : stmt->set(
933 : "insert or ignore into cache_reduced_multipole (k, l1, l2, value) values (?1, "
934 52 : "?2, ?3, ?4);");
935 26 : stmt->prepare();
936 : }
937 :
938 238 : for (auto &cached : cache_reduced_multipole_missing) {
939 :
940 193 : double val = std::pow(-1, cached.l[0]) *
941 193 : std::sqrt((2 * cached.l[0] + 1) * (2 * cached.l[1] + 1)) *
942 193 : WignerSymbols::wigner3j(
943 193 : cached.l[0], cached.kappa, cached.l[1], 0, 0,
944 193 : 0); // TODO call WignerSymbols::wigner3j(cached.kappa,
945 : // cached.l[1], 0, 0, 0) and loop over the resulting vector
946 193 : cache_reduced_multipole.insert({cached, val});
947 :
948 193 : if (!dbname.empty()) {
949 78 : stmt->bind(1, cached.kappa);
950 78 : stmt->bind(2, cached.l[0]);
951 78 : stmt->bind(3, cached.l[1]);
952 78 : stmt->bind(4, val);
953 78 : stmt->step();
954 78 : stmt->reset();
955 : }
956 : }
957 :
958 45 : cache_reduced_multipole_missing.clear();
959 : }
960 :
961 285 : if (!dbname.empty()) {
962 230 : stmt->exec("commit transaction;");
963 : }
964 :
965 285 : return 1;
966 : }
967 :
968 6 : size_t MatrixElementCache::size() {
969 6 : return cache_radial.size() + cache_angular.size() + cache_reduced_commutes_s.size() +
970 6 : cache_reduced_commutes_l.size() + cache_reduced_multipole.size();
971 : }
|
