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 "HamiltonianTwo.hpp"
21 :
22 : #include "Constants.hpp"
23 :
24 : #include <fmt/format.h>
25 :
26 : #include <stdexcept>
27 : #include <utility>
28 :
29 : template <typename Scalar>
30 2 : HamiltonianTwo<Scalar>::HamiltonianTwo(
31 : const Configuration &config, fs::path &path_cache,
32 : const std::shared_ptr<HamiltonianOne<Scalar>> &hamiltonian_one)
33 : : hamiltonian_one1(hamiltonian_one), hamiltonian_one2(hamiltonian_one),
34 2 : path_cache(path_cache) { // TODO
35 :
36 2 : samebasis = true;
37 :
38 2 : calculate(config);
39 2 : }
40 :
41 : template <typename Scalar>
42 0 : HamiltonianTwo<Scalar>::HamiltonianTwo(const Configuration &config, fs::path &path_cache,
43 : std::shared_ptr<HamiltonianOne<Scalar>> hamiltonian_one1,
44 : std::shared_ptr<HamiltonianOne<Scalar>> hamiltonian_one2)
45 0 : : hamiltonian_one1(std::move(hamiltonian_one1)), hamiltonian_one2(std::move(hamiltonian_one2)),
46 0 : path_cache(path_cache) {
47 :
48 0 : samebasis = false;
49 :
50 0 : calculate(config);
51 0 : }
52 :
53 : template <typename Scalar>
54 2 : void HamiltonianTwo<Scalar>::calculate(const Configuration &conf_tot) {
55 4 : fs::path path_cache_mat;
56 : if (utils::is_complex<Scalar>::value) {
57 0 : path_cache_mat = path_cache / "cache_matrix_complex";
58 : } else {
59 2 : path_cache_mat = path_cache / "cache_matrix_real";
60 : }
61 2 : if (!fs::exists(path_cache_mat)) {
62 0 : fs::create_directory(path_cache_mat);
63 : }
64 :
65 2 : double tol = 1e-32;
66 :
67 2 : if (hamiltonian_one1->size() != hamiltonian_one2->size()) {
68 0 : throw std::runtime_error(
69 : "The number of single atom Hamiltonians must be the same for both atoms.");
70 : }
71 :
72 2 : size_t nSteps_one = hamiltonian_one1->size();
73 :
74 : // --- generate configuration ---
75 4 : std::vector<Configuration> conf_mat;
76 2 : conf_mat.reserve(nSteps_one);
77 :
78 : // new, pair hamiltonian specific configuration
79 :
80 2 : if (samebasis) {
81 2 : this->basis = std::make_shared<BasisnamesTwo>(hamiltonian_one1->names()); // TODO remove
82 : } else {
83 0 : this->basis = std::make_shared<BasisnamesTwo>(hamiltonian_one1->names(),
84 0 : hamiltonian_one2->names()); // TODO remove
85 : }
86 4 : Configuration conf_matpair = this->basis->getConf();
87 2 : conf_matpair["deltaEPair"] = conf_tot["deltaEPair"];
88 2 : conf_matpair["deltaNPair"] = conf_tot["deltaNPair"];
89 2 : conf_matpair["deltaLPair"] = conf_tot["deltaLPair"];
90 2 : conf_matpair["deltaJPair"] = conf_tot["deltaJPair"];
91 2 : conf_matpair["deltaMPair"] = conf_tot["deltaMPair"];
92 : // conf_matpair["conserveM"] = conf_tot["conserveM"];
93 : // conf_matpair["conserveParityL"] = conf_tot["conserveParityL"];
94 2 : conf_matpair["exponent"] = conf_tot["exponent"];
95 :
96 4 : for (size_t i = 0; i < nSteps_one; ++i) {
97 : // old, single atom hamiltonian specific configuration
98 4 : Configuration conf_matsingle = *hamiltonian_one1->getParams(i); // TODO
99 2 : conf_matsingle += conf_matpair;
100 2 : conf_mat.push_back(conf_matsingle);
101 : // conf_mat.push_back(conf_matsingle + conf_matpair); // conf_matpair overwrites settings in
102 : // conf_matsingle // TODO
103 : }
104 :
105 : // setup variables
106 4 : conf_mat.back()["species1"] >>
107 2 : species1; // TODO order state inside cinfiguration object config.order()
108 2 : conf_mat.back()["species2"] >> species2; // TODO order state inside cinfiguration object
109 2 : conf_mat.back()["deltaEPair"] >> deltaE;
110 2 : conf_mat.back()["deltaNPair"] >> deltaN;
111 2 : conf_mat.back()["deltaLPair"] >> deltaL;
112 2 : conf_mat.back()["deltaJPair"] >> deltaJ;
113 2 : conf_mat.back()["deltaMPair"] >> deltaM;
114 : // conserveM = conf_tot["conserveM"].str() == "true";
115 : // conserveParityL = conf_tot["conserveParityL"].str() == "true";
116 2 : conf_tot["steps"] >> nSteps_two;
117 2 : conf_tot["minR"] >> min_R;
118 2 : conf_tot["maxR"] >> max_R;
119 2 : conf_tot["exponent"] >> multipoleexponent;
120 :
121 : double minEx, minEy, minEz, maxEx, maxEy, maxEz, minBx, minBy, minBz, maxBx, maxBy, maxBz;
122 2 : conf_tot["minEx"] >> minEx;
123 2 : conf_tot["minEy"] >> minEy;
124 2 : conf_tot["minEz"] >> minEz;
125 2 : conf_tot["maxEx"] >> maxEx;
126 2 : conf_tot["maxEy"] >> maxEy;
127 2 : conf_tot["maxEz"] >> maxEz;
128 2 : conf_tot["minBx"] >> minBx;
129 2 : conf_tot["minBy"] >> minBy;
130 2 : conf_tot["minBz"] >> minBz;
131 2 : conf_tot["maxBx"] >> maxBx;
132 2 : conf_tot["maxBy"] >> maxBy;
133 2 : conf_tot["maxBz"] >> maxBz;
134 : // bool fields_change_m = (minEx != 0) || (minEy != 0) || (maxEx != 0) || (maxEy != 0) || (minBx
135 : // != 0) || (minBy != 0) || (maxBx != 0) || (maxBy != 0); // TODO wie richtig? so ist es eine
136 : // variable, die von mehreren matrizen abhaengt bool fields_change_l = (minEx != 0) || (minEy !=
137 : // 0) || (minEz != 0) || (maxEx != 0) || (maxEy != 0) || (maxEz != 0); // TODO wie richtig? so
138 : // ist es eine variable, die von mehreren matrizen abhaengt
139 :
140 : // fields_change_m = true; // TODO
141 : // fields_change_l = true; // TODO
142 :
143 4 : std::vector<parity_t> sym_inversion;
144 2 : if (conf_tot["invE"].str() == "true") {
145 2 : sym_inversion.push_back(EVEN);
146 : }
147 2 : if (conf_tot["invO"].str() == "true") {
148 2 : sym_inversion.push_back(ODD);
149 : }
150 2 : if (sym_inversion.empty()) {
151 0 : sym_inversion.push_back(NA);
152 : }
153 :
154 4 : std::vector<parity_t> sym_permutation;
155 2 : if (conf_tot["perE"].str() == "true") {
156 2 : sym_permutation.push_back(EVEN);
157 : }
158 2 : if (conf_tot["perO"].str() == "true") {
159 2 : sym_permutation.push_back(ODD);
160 : }
161 2 : if (sym_permutation.empty()) {
162 0 : sym_permutation.push_back(NA);
163 : }
164 :
165 4 : std::vector<parity_t> sym_reflection;
166 2 : if (conf_tot["refE"].str() == "true") {
167 0 : sym_reflection.push_back(EVEN);
168 : }
169 2 : if (conf_tot["refO"].str() == "true") {
170 0 : sym_reflection.push_back(ODD);
171 : }
172 2 : if (sym_reflection.empty()) {
173 2 : sym_reflection.push_back(NA);
174 : }
175 :
176 2 : bool conserveM = conf_tot["conserveM"].str() == "true";
177 :
178 2 : bool sametrafo = conf_tot["sametrafo"].str() == "true";
179 :
180 2 : bool zerotheta = conf_tot["zerotheta"].str() == "true";
181 :
182 2 : if (min_R == max_R && nSteps_one == 1) {
183 0 : nSteps_two = 1;
184 : }
185 :
186 : ////////////////////////////////////////////////////////
187 : ////// Restrict single atom states /////////////////////
188 : ////////////////////////////////////////////////////////
189 :
190 : // === Restrict states of atom 1 ===
191 :
192 4 : auto basis_one1 = hamiltonian_one1->names();
193 4 : std::vector<StateOneOld> initial1 = basis_one1->initial();
194 4 : std::vector<bool> necessary1(basis_one1->size(), false);
195 :
196 226 : for (const auto &state : *basis_one1) {
197 224 : bool validN = false;
198 224 : bool validL = false;
199 224 : bool validJ = false;
200 224 : bool validM = false;
201 :
202 672 : for (const auto &initial : initial1) {
203 448 : if (deltaN < 0 || std::abs(state.n - initial.n) <= deltaN) {
204 448 : validN = true;
205 : }
206 448 : if (deltaL < 0 || std::abs(state.l - initial.l) <= deltaL) {
207 448 : validL = true;
208 : }
209 448 : if (deltaJ < 0 || std::abs(state.j - initial.j) <= deltaJ) {
210 448 : validJ = true;
211 : }
212 448 : if (deltaM < 0 || std::abs(state.m - initial.m) <= deltaM) {
213 448 : validM = true;
214 : }
215 : }
216 :
217 224 : if (validN && validL && validJ && validM) {
218 224 : necessary1[state.idx] = true;
219 : }
220 : }
221 :
222 2 : hamiltonian_one1->removeUnnecessaryStates(necessary1);
223 :
224 : // === Restrict states of atom 2 ===
225 :
226 2 : if (!samebasis) {
227 0 : auto basis_one2 = hamiltonian_one2->names();
228 0 : std::vector<StateOneOld> initial2 = basis_one2->initial();
229 0 : std::vector<bool> necessary2(basis_one2->size(), false);
230 :
231 0 : for (const auto &state : *basis_one2) {
232 0 : bool validN = false;
233 0 : bool validL = false;
234 0 : bool validJ = false;
235 0 : bool validM = false;
236 :
237 0 : for (const auto &initial : initial2) {
238 0 : if (deltaN < 0 || std::abs(state.n - initial.n) <= deltaN) {
239 0 : validN = true;
240 : }
241 0 : if (deltaL < 0 || std::abs(state.l - initial.l) <= deltaL) {
242 0 : validL = true;
243 : }
244 0 : if (deltaJ < 0 || std::abs(state.j - initial.j) <= deltaJ) {
245 0 : validJ = true;
246 : }
247 0 : if (deltaM < 0 || std::abs(state.m - initial.m) <= deltaM) {
248 0 : validM = true;
249 : }
250 : }
251 :
252 0 : if (validN && validL && validJ && validM) {
253 0 : necessary2[state.idx] = true;
254 : }
255 : }
256 :
257 0 : hamiltonian_one2->removeUnnecessaryStates(necessary2);
258 : }
259 :
260 : ////////////////////////////////////////////////////////
261 : ////// Build pair state basis //////////////////////////
262 : ////////////////////////////////////////////////////////
263 :
264 : // === Build pair state basis ===
265 :
266 2 : std::cout << "Two-atom Hamiltonian, build pair state basis" << std::endl;
267 :
268 2 : if (samebasis) {
269 2 : this->basis = std::make_shared<BasisnamesTwo>(hamiltonian_one1->names());
270 : } else {
271 0 : this->basis =
272 0 : std::make_shared<BasisnamesTwo>(hamiltonian_one1->names(), hamiltonian_one2->names());
273 : }
274 :
275 2 : std::cout << "Two-atom Hamiltonian, basis size without restrictions: " << this->basis->size()
276 2 : << std::endl;
277 :
278 : // === Determine necessary symmetries ===
279 2 : std::cout << "Two-atom Hamiltonian, determine symmetrized subspaces" << std::endl;
280 :
281 4 : StateTwoOld initial = this->basis->initial();
282 2 : parity_t initalParityL = (std::pow(-1, initial.l[0] + initial.l[1]) > 0) ? EVEN : ODD;
283 2 : int initalM = initial.m[0] + initial.m[1];
284 2 : int initalJ = initial.j[0] + initial.j[1];
285 2 : bool samestates = initial.first() == initial.second();
286 :
287 4 : std::vector<int> sym_rotation;
288 2 : if (conserveM) {
289 2 : if (!sametrafo) {
290 0 : for (const auto &state : *this->basis) {
291 0 : sym_rotation.push_back(state.m[0] + state.m[1]);
292 : }
293 2 : } else if (!zerotheta) {
294 0 : for (int M = -initalJ; M <= initalJ; ++M) {
295 0 : sym_rotation.push_back(M);
296 : }
297 : } else {
298 2 : sym_rotation.push_back(initalM);
299 : }
300 : } else {
301 0 : sym_rotation.push_back(NA);
302 : }
303 :
304 : Symmetry sym;
305 4 : std::set<Symmetry> symmetries_set;
306 6 : for (const parity_t &inv : sym_inversion) {
307 4 : sym.inversion = inv;
308 :
309 : // In case of even inversion symmetry and the same inital state for the first and second
310 : // atom: the inital state ist not contained in the corresponding block
311 4 : if (sym.inversion == EVEN && samestates && sametrafo) {
312 0 : continue;
313 : }
314 :
315 12 : for (const parity_t &per : sym_permutation) {
316 8 : sym.permutation = per;
317 :
318 : // In case of even permutation symmetry and the same inital state for the first and
319 : // second atom: the inital state ist not contained in the corresponding block
320 8 : if (sym.permutation == EVEN && samestates && sametrafo) {
321 0 : continue;
322 : }
323 :
324 : // In case of inversion and permutation symmetry: the orbital parity is conserved and we
325 : // just use the blocks with the same parity as the inital state
326 8 : if (sym.inversion != NA && sym.permutation != NA && sametrafo) {
327 8 : if (sym.inversion * sym.permutation != initalParityL) {
328 4 : continue;
329 : }
330 : }
331 :
332 8 : for (const parity_t &ref : sym_reflection) {
333 4 : sym.reflection = ref;
334 8 : for (const int &rot : sym_rotation) {
335 :
336 : // In case of reflection symmetry: do just use the absolute value of rot
337 4 : if (sym.reflection != NA) {
338 0 : sym.rotation = std::abs(rot);
339 : } else {
340 4 : sym.rotation = rot;
341 : }
342 :
343 4 : symmetries_set.insert(sym);
344 : }
345 : }
346 : }
347 : }
348 4 : std::vector<Symmetry> symmetries(symmetries_set.begin(), symmetries_set.end());
349 :
350 : // === Build up the list of necessary pair states ===
351 2 : std::cout << "Two-atom Hamiltonian, build up the list of necessary pair states" << std::endl;
352 :
353 : // Apply energy cutoff
354 4 : std::vector<bool> necessary_tmp(this->basis->size(), false);
355 :
356 2 : auto nSteps_one_i = static_cast<int>(nSteps_one);
357 :
358 2 : #pragma omp parallel for
359 : for (int i = 0; i < nSteps_one_i; ++i) {
360 : energycutoff(*(hamiltonian_one1->get(i)), *(hamiltonian_one2->get(i)), deltaE,
361 : necessary_tmp);
362 : }
363 :
364 : // Apply restrictions due to symmetries
365 4 : std::vector<bool> necessary(this->basis->size(), false);
366 :
367 50178 : for (const auto &state : *this->basis) {
368 150528 : for (Symmetry sym : symmetries) {
369 100352 : float M = state.m[0] + state.m[1];
370 100352 : int parityL = std::pow(-1, state.l[0] + state.l[1]);
371 :
372 : // In case of rotation symmetry: skip pair states with wrong total magnetic momentum
373 100352 : if (sym.rotation != NA && sym.rotation != M &&
374 87262 : (sym.reflection == NA || sym.rotation != -M)) {
375 87262 : continue;
376 : }
377 :
378 : // In case of inversion and permutation symmetry: skip pair states with wrong orbital
379 : // parity
380 13090 : if (sym.inversion != NA && sym.permutation != NA && parityL != initalParityL &&
381 : sametrafo) {
382 6690 : continue;
383 : }
384 :
385 6400 : necessary[state.idx] = necessary_tmp[state.idx];
386 : }
387 : }
388 :
389 2 : int numNecessary = std::count(necessary.begin(), necessary.end(), true);
390 2 : std::cout << "Two-atom Hamiltonian, basis size with restrictions: " << numNecessary
391 2 : << std::endl;
392 2 : std::cout << fmt::format(">>BAS{:7d}", numNecessary) << std::endl;
393 :
394 : // === Save pair state basis ===
395 2 : std::cout << "Two-atom Hamiltonian, save pair state basis" << std::endl;
396 :
397 : // initialize uuid generator
398 : boost::uuids::random_generator generator;
399 :
400 : // generate uuid
401 4 : std::string uuid;
402 2 : boost::uuids::uuid u = generator();
403 2 : boost::algorithm::hex(u.begin(), u.end(), std::back_inserter(uuid));
404 :
405 : // save pair state basis
406 4 : fs::path path_basis = fs::temp_directory_path();
407 2 : path_basis /= "basis_two_" + uuid + ".csv";
408 4 : this->basis->save(
409 : path_basis
410 4 : .string()); // TODO save only necessary entries, i.e. save pair state basis in sparse
411 : // format (possibility, remove basis states but keep their idx - this would
412 : // also make "if (necessary) continue" unneeded; then, "combine" has to
413 : // check existence of basis element and the python script has to be adapted)
414 :
415 4 : std::cout << fmt::format(">>STA {:s}", path_basis.string()) << std::endl;
416 :
417 : ////////////////////////////////////////////////////////
418 : ////// Construct atom-atom interaction /////////////////
419 : ////////////////////////////////////////////////////////
420 :
421 : // Construct pair Hamiltonians for all orders of the multipole expansion
422 :
423 4 : std::vector<int> exponent_multipole;
424 4 : std::vector<Hamiltonianmatrix<Scalar>> mat_multipole;
425 4 : MatrixElements matrixelements_atom1(conf_tot, species1,
426 4 : (path_cache / "cache_elements.db").string());
427 4 : MatrixElements matrixelements_atom2(conf_tot, species2,
428 4 : (path_cache / "cache_elements.db").string());
429 4 : std::vector<idx_t> size_mat_multipole;
430 :
431 2 : int idx_multipole_max = -1;
432 :
433 2 : if (multipoleexponent > 2) {
434 :
435 : // --- Initialize two-atom interaction Hamiltonians ---
436 2 : std::cout << "Two-atom Hamiltonian, initialize interaction Hamiltonians" << std::endl;
437 :
438 2 : int kappa_min = 1; // spherical dipole operators
439 2 : int kappa_max = multipoleexponent - kappa_min - 1;
440 2 : int sumOfKappas_min = kappa_min + kappa_min;
441 2 : int sumOfKappas_max = kappa_max + kappa_min;
442 2 : idx_multipole_max = sumOfKappas_max - sumOfKappas_min;
443 :
444 2 : exponent_multipole.reserve(idx_multipole_max + 1);
445 2 : mat_multipole.reserve(idx_multipole_max + 1);
446 2 : size_mat_multipole.resize(idx_multipole_max + 1);
447 :
448 : // --- Precalculate matrix elements --- // TODO parallelization
449 2 : std::cout << "Two-atom Hamiltonian, get one-atom states needed for the pair state basis"
450 2 : << std::endl;
451 :
452 4 : auto basis_one1_needed =
453 4 : std::make_shared<BasisnamesOne>(BasisnamesOne::fromFirst(this->basis));
454 4 : auto basis_one2_needed =
455 4 : std::make_shared<BasisnamesOne>(BasisnamesOne::fromSecond(this->basis));
456 :
457 4 : for (int kappa = kappa_min; kappa <= kappa_max; ++kappa) {
458 2 : std::cout << "Two-atom Hamiltonian, precalculate matrix elements for kappa = " << kappa
459 2 : << std::endl;
460 2 : matrixelements_atom1.precalculateMultipole(basis_one1_needed, kappa);
461 2 : matrixelements_atom2.precalculateMultipole(basis_one2_needed, kappa);
462 : }
463 :
464 : // TODO if (samebasis) ...
465 :
466 : // --- Count entries of two-atom interaction Hamiltonians ---
467 : std::cout
468 2 : << "Two-atom Hamiltonian, count number of entries within the interaction Hamiltonians"
469 2 : << std::endl;
470 :
471 4 : for (int sumOfKappas = sumOfKappas_min; sumOfKappas <= sumOfKappas_max; ++sumOfKappas) {
472 2 : int idx_multipole = sumOfKappas - sumOfKappas_min;
473 :
474 50178 : for (const auto &state_col : *this->basis) { // TODO parallelization
475 50176 : if (!necessary[state_col.idx]) {
476 50056 : continue;
477 : }
478 :
479 120 : int M_col = state_col.first().m + state_col.second().m;
480 :
481 6021240 : for (const auto &state_row : *this->basis) {
482 6021120 : if (!necessary[state_row.idx]) {
483 6006720 : continue;
484 : }
485 :
486 14400 : if (state_row.idx < state_col.idx) {
487 7140 : continue;
488 : }
489 7260 : int M_row = state_row.first().m + state_row.second().m;
490 7260 : if (M_col != M_row) {
491 0 : continue;
492 : }
493 :
494 : // multipole interaction with 1/R^(sumOfKappas+1) = 1/R^(idx_multipole+3) decay
495 14392 : for (int kappa1 = kappa_min; kappa1 <= sumOfKappas - 1; ++kappa1) {
496 7260 : int kappa2 = sumOfKappas - kappa1;
497 :
498 : // allowed deltaL, deltaJ, and deltaM?
499 15128 : if (selectionRulesMultipole(state_row.first(), state_col.first(), kappa1) &&
500 7868 : selectionRulesMultipole(state_row.second(), state_col.second(),
501 : kappa2)) {
502 128 : int q1 = state_row.first().m - state_col.first().m;
503 128 : int q2 = state_row.second().m - state_col.second().m;
504 :
505 : // total momentum preserved?
506 128 : if (q1 == -q2) {
507 128 : size_mat_multipole[idx_multipole]++;
508 128 : break;
509 : }
510 : }
511 : }
512 : }
513 : }
514 : }
515 :
516 : // --- Construct two-atom interaction Hamiltonians ---
517 2 : size_t size_basis = this->basis->size();
518 :
519 4 : for (int sumOfKappas = sumOfKappas_min; sumOfKappas <= sumOfKappas_max; ++sumOfKappas) {
520 : std::cout
521 2 : << "Two-atom Hamiltonian, construct interaction Hamiltonian that belongs to 1/R^"
522 2 : << sumOfKappas + 1 << std::endl;
523 :
524 2 : int idx_multipole = sumOfKappas - sumOfKappas_min;
525 :
526 2 : exponent_multipole.push_back(sumOfKappas + 1);
527 2 : mat_multipole.emplace_back(
528 : size_basis,
529 2 : 2 * size_mat_multipole[idx_multipole]); // factor of 2 because triangular matrix is
530 : // not sufficient
531 :
532 50178 : for (const auto &state_col : *this->basis) { // TODO parallelization
533 50176 : if (!necessary[state_col.idx]) {
534 50056 : continue;
535 : }
536 :
537 120 : int M_col = state_col.first().m + state_col.second().m;
538 :
539 6021240 : for (const auto &state_row : *this->basis) {
540 6021120 : if (!necessary[state_row.idx]) {
541 6006720 : continue;
542 : }
543 :
544 14400 : if (state_row.idx < state_col.idx) {
545 7140 : continue;
546 : }
547 7260 : int M_row = state_row.first().m + state_row.second().m;
548 7260 : if (M_col != M_row) {
549 0 : continue;
550 : }
551 :
552 : // construct basis
553 7260 : if (state_row.idx == state_col.idx) {
554 120 : mat_multipole[idx_multipole].addBasis(state_row.idx, state_col.idx, 1);
555 : }
556 :
557 : // multipole interaction with 1/R^(sumOfKappas+1) = 1/R^(idx_multipole+3) decay
558 7260 : double val = 0;
559 :
560 14520 : for (int kappa1 = kappa_min; kappa1 <= sumOfKappas - 1; ++kappa1) {
561 7260 : int kappa2 = sumOfKappas - kappa1;
562 :
563 : // allowed deltaL, deltaJ, and deltaM?
564 15128 : if (selectionRulesMultipole(state_row.first(), state_col.first(), kappa1) &&
565 7868 : selectionRulesMultipole(state_row.second(), state_col.second(),
566 : kappa2)) {
567 128 : int q1 = state_row.first().m - state_col.first().m;
568 128 : int q2 = state_row.second().m - state_col.second().m;
569 :
570 : // total momentum preserved?
571 128 : if (q1 == -q2) {
572 384 : double binomials = boost::math::binomial_coefficient<double>(
573 128 : kappa1 + kappa2, kappa1 + q1) *
574 256 : boost::math::binomial_coefficient<double>(kappa1 + kappa2,
575 128 : kappa2 - q2);
576 128 : val += inverse_electric_constant * std::pow(-1, kappa2) *
577 256 : std::sqrt(binomials) *
578 128 : matrixelements_atom1.getMultipole(state_row.first(),
579 256 : state_col.first(), kappa1) *
580 128 : matrixelements_atom2.getMultipole(state_row.second(),
581 256 : state_col.second(), kappa2);
582 : }
583 : }
584 : }
585 :
586 7260 : if (std::abs(val) > tol) {
587 128 : mat_multipole[idx_multipole].addEntries(state_row.idx, state_col.idx, val);
588 128 : if (state_row.idx != state_col.idx) {
589 128 : mat_multipole[idx_multipole].addEntries(
590 128 : state_col.idx, state_row.idx,
591 : val); // triangular matrix is not sufficient because of basis change
592 : }
593 : }
594 :
595 : // TODO state_two soll std::array<state_one, 2> sein! Dann geht auch die Abfrage
596 : // der selection rules eindeutiger
597 : }
598 : }
599 :
600 : std::cout
601 2 : << "Two-atom Hamiltonian, compress interaction Hamiltonian that belongs to 1/R^"
602 2 : << sumOfKappas + 1 << std::endl;
603 :
604 6 : mat_multipole[idx_multipole].compress(
605 2 : this->basis->dim(),
606 2 : this->basis->dim()); // TODO substitute dim() by size()
607 : }
608 : }
609 :
610 : ////////////////////////////////////////////////////////
611 : ////// Prepare processing of Hamiltonians //////////////
612 : ////////////////////////////////////////////////////////
613 :
614 : // TODO Put the logic in its own class
615 :
616 2 : std::cout << "Two-atom Hamiltonian, process Hamiltonians" << std::endl;
617 :
618 : // === Open database ===
619 4 : fs::path path_db;
620 :
621 : if (utils::is_complex<Scalar>::value) {
622 0 : path_db = path_cache / "cache_matrix_complex.db";
623 : } else {
624 2 : path_db = path_cache / "cache_matrix_real.db";
625 : }
626 4 : sqlite::handle db(path_db.string());
627 4 : sqlite::statement stmt(db);
628 :
629 : // === Initialize variables ===
630 2 : bool flag_perhapsmissingtable = true;
631 :
632 2 : this->matrix_path.resize(nSteps_two * symmetries.size());
633 :
634 2 : auto indices_symmetry_i = static_cast<int>(symmetries.size());
635 :
636 : // --- Determine combined single atom matrices ---
637 : // Construct pair Hamiltonian consistent of combined one-atom Hamiltonians (1 x Hamiltonian2 +
638 : // Hamiltonian1 x 1)
639 :
640 4 : std::vector<Hamiltonianmatrix<Scalar>> mat_single;
641 :
642 : // Check if one_atom Hamiltonians change with step_two
643 : // It is assumed that nSteps_one = 1 if nSteps_two != nSteps_one // TODO introduce variable
644 : // "is_mat_single_const" to improve readability
645 2 : if (nSteps_two != nSteps_one) {
646 : std::cout
647 2 : << "Two-atom Hamiltonian, construct contribution of combined one-atom Hamiltonians"
648 2 : << std::endl;
649 :
650 2 : mat_single.resize(symmetries.size());
651 :
652 2 : #pragma omp parallel for
653 : for (int idx_symmetry = 0; idx_symmetry < indices_symmetry_i; ++idx_symmetry) {
654 : Symmetry sym = symmetries[idx_symmetry];
655 :
656 : // Combine the Hamiltonians of the two atoms
657 : mat_single[idx_symmetry] = combine(
658 : *(hamiltonian_one1->get(0)), *(hamiltonian_one2->get(0)), deltaE, this->basis, sym);
659 :
660 : // Remove more or less empty basis vectors
661 : mat_single[idx_symmetry].removeUnnecessaryBasisvectors();
662 : }
663 : }
664 :
665 : // --- Determine transformed interaction matrices ---
666 4 : std::vector<Hamiltonianmatrix<Scalar>> mat_multipole_transformed;
667 :
668 : // Check if one_atom Hamiltonians change with step_two
669 2 : if (nSteps_two != nSteps_one) {
670 2 : std::cout << "Two-atom Hamiltonian, construct transformed interaction matrices"
671 2 : << std::endl;
672 :
673 2 : mat_multipole_transformed.resize(symmetries.size() * (idx_multipole_max + 1));
674 :
675 2 : #pragma omp parallel for
676 : for (int idx_symmetry = 0; idx_symmetry < indices_symmetry_i; ++idx_symmetry) {
677 : for (int idx_multipole = 0; idx_multipole <= idx_multipole_max; ++idx_multipole) {
678 : mat_multipole_transformed[idx_symmetry * (idx_multipole_max + 1) + idx_multipole] =
679 : mat_multipole[idx_multipole].changeBasis(mat_single[idx_symmetry].basis());
680 : }
681 : }
682 : }
683 :
684 : ////////////////////////////////////////////////////////
685 : ////// Loop through steps and symmetries ///////////////
686 : ////////////////////////////////////////////////////////
687 :
688 4 : std::cout << fmt::format(">>TOT{:7d}", nSteps_two * symmetries.size()) << std::endl;
689 :
690 2 : auto nSteps_two_i = static_cast<int>(nSteps_two);
691 :
692 2 : #pragma omp parallel for schedule(static, 1)
693 :
694 : // Loop through steps
695 : for (int step_two = 0; step_two < nSteps_two_i; ++step_two) {
696 :
697 : // Loop through symmetries
698 : for (size_t idx_symmetry = 0; idx_symmetry < symmetries.size(); ++idx_symmetry) {
699 : Symmetry sym = symmetries[idx_symmetry];
700 :
701 : size_t step = step_two * symmetries.size() + idx_symmetry;
702 :
703 : // === Get parameters for the current position inside the loop ===
704 : int single_idx = (nSteps_two == nSteps_one) ? step_two : 0;
705 :
706 : // Get interatomic distance
707 : double normalized_position = (nSteps_two > 1) ? step_two / (nSteps_two - 1.) : 0;
708 : double position = min_R + normalized_position * (max_R - min_R);
709 :
710 : // Get configuration and save postions and symmetries
711 : Configuration conf = conf_mat[single_idx];
712 : conf["R"] << position;
713 : conf["inversion"] << ((sym.inversion == NA) ? "" : std::to_string(sym.inversion));
714 : conf["permutation"] << ((sym.permutation == NA) ? "" : std::to_string(sym.permutation));
715 : conf["reflection"] << ((sym.reflection == NA) ? "" : std::to_string(sym.reflection));
716 : conf["rotation"] << ((sym.rotation == NA) ? "" : std::to_string(sym.rotation));
717 :
718 : // === Create table if necessary ===
719 : std::stringstream query;
720 : std::string spacer;
721 :
722 : if (flag_perhapsmissingtable) {
723 : query << "CREATE TABLE IF NOT EXISTS cache_two (uuid text NOT NULL PRIMARY KEY, "
724 : "created TIMESTAMP DEFAULT CURRENT_TIMESTAMP, "
725 : "accessed TIMESTAMP DEFAULT CURRENT_TIMESTAMP";
726 : for (auto p : conf) {
727 : query << ", " << p.first << " text";
728 : }
729 : query << ", UNIQUE (";
730 : for (auto p : conf) {
731 : query << spacer << p.first;
732 : spacer = ", ";
733 : }
734 : query << "));";
735 :
736 : #pragma omp critical(database)
737 : stmt.exec(query.str());
738 :
739 : flag_perhapsmissingtable = false;
740 : }
741 :
742 : // === Get uuid as filename ===
743 : std::string uuid;
744 : query.str(std::string());
745 : spacer = "";
746 : query << "SELECT uuid FROM cache_two WHERE ";
747 : for (auto p : conf) {
748 : query << spacer << p.first << "='" << p.second.str() << "'";
749 : spacer = " AND ";
750 : }
751 : query << ";";
752 :
753 : #pragma omp critical(database)
754 : {
755 : sqlite::statement stmt(db, query.str());
756 : stmt.prepare();
757 : if (stmt.step()) {
758 : uuid = stmt.get<std::string>(0);
759 : }
760 : }
761 :
762 : if (!uuid.empty()) {
763 : query.str(std::string());
764 : query << "UPDATE cache_two SET accessed = CURRENT_TIMESTAMP WHERE uuid = '" << uuid
765 : << "';";
766 : #pragma omp critical(database)
767 : stmt.exec(query.str()); // TODO check whether this slows down the program
768 :
769 : } else {
770 : boost::uuids::uuid u = generator();
771 : boost::algorithm::hex(u.begin(), u.end(), std::back_inserter(uuid));
772 :
773 : query.str(std::string());
774 : query << "INSERT INTO cache_two (uuid";
775 : for (auto p : conf) {
776 : query << ", " << p.first;
777 : }
778 : query << ") values ( '" << uuid << "'";
779 : for (auto p : conf) {
780 : query << ", "
781 : << "'" << p.second.str() << "'";
782 : }
783 : query << ");";
784 : #pragma omp critical(database)
785 : stmt.exec(query.str());
786 : }
787 :
788 : // === Check existence of files ===
789 :
790 : // Check whether .mat and .json file exists and compare settings in program with
791 : // settings in .json file
792 : fs::path path, path_mat, path_json;
793 :
794 : path = path_cache_mat / ("two_" + uuid);
795 : path_mat = path;
796 : path_mat.replace_extension(".mat");
797 : path_json = path;
798 : path_json.replace_extension(".json");
799 :
800 : bool is_existing = false;
801 : if (fs::exists(path_mat)) {
802 : if (fs::exists(path_json)) {
803 : Configuration params_loaded;
804 : params_loaded.load_from_json(path_json.string());
805 : if (conf == params_loaded) {
806 : is_existing = true;
807 : }
808 : }
809 : }
810 :
811 : // Create .json file if "is_existing" is false
812 : if (!is_existing) {
813 : conf.save_to_json(path_json.string());
814 : }
815 :
816 : // === Build and diagonalize total matrix if not existent ===
817 : Hamiltonianmatrix<Scalar> totalmatrix;
818 :
819 : if (!is_existing || !totalmatrix.load(path_mat.string())) {
820 :
821 : // --- Combine single atom matrices ---
822 : if (nSteps_two == nSteps_one) {
823 : totalmatrix =
824 : combine(*(hamiltonian_one1->get(step_two)),
825 : *(hamiltonian_one2->get(step_two)), deltaE, this->basis, sym);
826 : totalmatrix.removeUnnecessaryBasisvectors();
827 : } else {
828 : totalmatrix = mat_single[idx_symmetry];
829 : }
830 :
831 : // --- Add interaction ---
832 : for (int idx_multipole = 0; idx_multipole <= idx_multipole_max; ++idx_multipole) {
833 : double pos = 1. / std::pow(position, exponent_multipole[idx_multipole]);
834 : if (nSteps_two == nSteps_one) {
835 : totalmatrix +=
836 : mat_multipole[idx_multipole].changeBasis(totalmatrix.basis()) * pos;
837 : } else {
838 : totalmatrix +=
839 : mat_multipole_transformed[idx_symmetry * (idx_multipole_max + 1) +
840 : idx_multipole] *
841 : pos;
842 : }
843 : }
844 :
845 : // Stdout: Hamiltonian assembled
846 : #pragma omp critical(textoutput)
847 : {
848 : std::cout << fmt::format(">>DIM{:7d}", totalmatrix.num_basisvectors())
849 : << std::endl;
850 : std::cout << "Two-atom Hamiltonian, " << step + 1 << ". Hamiltonian assembled"
851 : << std::endl;
852 : }
853 :
854 : // --- Diagonalize matrix and save diagonalized matrix ---
855 : totalmatrix.diagonalize();
856 : totalmatrix.save(path_mat.string());
857 :
858 : // Stdout: Hamiltonian diagonalized
859 : #pragma omp critical(textoutput)
860 : {
861 : std::cout << fmt::format(">>OUT{:7d}{:7d}{:7d}{:7d} {:s}", (step + 1), step_two,
862 : symmetries.size(), idx_symmetry, path.string())
863 : << std::endl;
864 : std::cout << "Two-atom Hamiltonian, " << step + 1
865 : << ". Hamiltonian diagonalized" << std::endl;
866 : }
867 : } else {
868 :
869 : // Stdout: Hamiltonian loaded
870 : #pragma omp critical(textoutput)
871 : {
872 : std::cout << fmt::format(">>DIM{:7d}", totalmatrix.num_basisvectors())
873 : << std::endl;
874 : std::cout << fmt::format(">>OUT{:7d}{:7d}{:7d}{:7d} {:s}", (step + 1), step_two,
875 : symmetries.size(), idx_symmetry, path.string())
876 : << std::endl;
877 : std::cout << "Two-atom Hamiltonian, " << step + 1 << ". Hamiltonian loaded"
878 : << std::endl;
879 : }
880 : }
881 :
882 : // === Store path to configuration and diagonalized matrix ===
883 : this->matrix_path[step] = path.string();
884 : }
885 : }
886 :
887 2 : std::cout << "Two-atom Hamiltonian, all Hamiltonians processed" << std::endl;
888 2 : }
889 :
890 : template class HamiltonianTwo<std::complex<double>>;
891 : template class HamiltonianTwo<double>;
|
