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 "HamiltonianOne.hpp"
21 : #include "filesystem.hpp"
22 :
23 : #include <fmt/format.h>
24 :
25 : #include <stdexcept>
26 : #include <utility>
27 :
28 : template <typename Scalar>
29 3 : HamiltonianOne<Scalar>::HamiltonianOne(const Configuration &config, fs::path &path_cache,
30 : std::shared_ptr<BasisnamesOne> basis_one)
31 3 : : path_cache(path_cache) {
32 3 : this->basis = std::move(basis_one);
33 3 : configure(config);
34 3 : build();
35 3 : }
36 :
37 : template <typename Scalar>
38 0 : const Configuration &HamiltonianOne<Scalar>::getConf()
39 : const { // TODO in Configurable Klasse auslagern, von der geerbt werrden soll
40 0 : return basicconf;
41 : }
42 :
43 : template <typename Scalar>
44 8 : void HamiltonianOne<Scalar>::changeToSpherical(double val_x, double val_y, double val_z,
45 : double &val_p, double &val_m, double &val_0) {
46 8 : if (val_y != 0) {
47 0 : std::string msg("For fields with non-zero y-coordinates, a complex data type is needed.");
48 0 : std::cout << fmt::format(">>ERR{:s}", msg) << std::endl;
49 0 : throw std::runtime_error(msg);
50 : }
51 8 : val_p = -val_x / std::sqrt(2);
52 8 : val_m = val_x / std::sqrt(2);
53 8 : val_0 = val_z;
54 8 : }
55 :
56 : template <typename Scalar>
57 4 : void HamiltonianOne<Scalar>::changeToSpherical(double val_x, double val_y, double val_z,
58 : std::complex<double> &val_p,
59 : std::complex<double> &val_m,
60 : std::complex<double> &val_0) {
61 4 : val_p = std::complex<double>(-val_x / std::sqrt(2), -val_y / std::sqrt(2));
62 4 : val_m = std::complex<double>(val_x / std::sqrt(2), -val_y / std::sqrt(2));
63 4 : val_0 = std::complex<double>(val_z, 0);
64 4 : }
65 :
66 : template <typename Scalar>
67 3 : void HamiltonianOne<Scalar>::configure(const Configuration &config) {
68 3 : basicconf = this->basis->getConf();
69 3 : basicconf["deltaESingle"] << config["deltaESingle"];
70 3 : basicconf["diamagnetism"] << config["diamagnetism"];
71 :
72 3 : basicconf["deltaESingle"] >> deltaE;
73 3 : basicconf["species1"] >> species;
74 :
75 3 : diamagnetism = basicconf["diamagnetism"].str() == "true";
76 :
77 3 : config["minBx"] >> min_B_x;
78 3 : config["minBy"] >> min_B_y;
79 3 : config["minBz"] >> min_B_z;
80 3 : config["minEx"] >> min_E_x;
81 3 : config["minEy"] >> min_E_y;
82 3 : config["minEz"] >> min_E_z;
83 3 : config["maxBx"] >> max_B_x;
84 3 : config["maxBy"] >> max_B_y;
85 3 : config["maxBz"] >> max_B_z;
86 3 : config["maxEx"] >> max_E_x;
87 3 : config["maxEy"] >> max_E_y;
88 3 : config["maxEz"] >> max_E_z;
89 :
90 3 : if ((min_B_x == max_B_x) && (min_B_y == max_B_y) && (min_B_z == max_B_z) &&
91 3 : (min_E_x == max_E_x) && (min_E_y == max_E_y) && (min_E_z == max_E_z)) {
92 3 : nSteps = 1;
93 : } else {
94 0 : config["steps"] >> nSteps;
95 : }
96 3 : }
97 :
98 : template <typename Scalar>
99 3 : void HamiltonianOne<Scalar>::build() {
100 6 : fs::path path_cache_mat;
101 : if (utils::is_complex<Scalar>::value) {
102 1 : path_cache_mat = path_cache / "cache_matrix_complex";
103 : } else {
104 2 : path_cache_mat = path_cache / "cache_matrix_real";
105 : }
106 3 : if (!fs::exists(path_cache_mat)) {
107 3 : fs::create_directory(path_cache_mat);
108 : }
109 :
110 3 : double tol = 1e-32;
111 :
112 : ////////////////////////////////////////////////////////
113 : ////// Build single atom basis and Hamiltonian /////////
114 : ////////////////////////////////////////////////////////
115 :
116 : // === Calculate one-atom Hamiltonian ===
117 :
118 : // --- Count entries of one-atom Hamiltonian ---
119 3 : size_t size_basis = this->basis->size();
120 3 : size_t size_energy = this->basis->size();
121 :
122 : // --- Construct one-atom Hamiltonian and basis ---
123 3 : std::cout << "One-atom Hamiltonian, construct diagonal Hamiltonian" << std::endl;
124 :
125 6 : Hamiltonianmatrix<Scalar> hamiltonian_energy(size_basis, size_energy);
126 :
127 3 : double energy_initial = 0;
128 9 : for (const auto &state : this->basis->initial()) {
129 6 : energy_initial += energy_level(species, state.n, state.l, state.j);
130 : }
131 3 : energy_initial /= this->basis->initial().size(); // TODO save it to the json file
132 :
133 6 : std::vector<bool> is_necessary(this->basis->size(), false);
134 3 : idx_t idx = 0;
135 1149 : for (const auto &state : *this->basis) {
136 1146 : double val = energy_level(species, state.n, state.l, state.j) - energy_initial;
137 1146 : if (std::abs(val) < deltaE + 1e-11 || deltaE < 0) { // TODO
138 374 : is_necessary[state.idx] = true;
139 374 : hamiltonian_energy.addEntries(idx, idx, val);
140 374 : hamiltonian_energy.addBasis(idx, idx, 1);
141 374 : ++idx;
142 : }
143 : }
144 3 : std::cout << "One-atom Hamiltonian, basis size without restrictions: " << this->basis->size()
145 3 : << std::endl;
146 :
147 3 : this->basis->removeUnnecessaryStates(is_necessary);
148 :
149 3 : hamiltonian_energy.compress(this->basis->dim(), this->basis->dim());
150 :
151 3 : std::cout << "One-atom Hamiltonian, basis size with restrictions: " << this->basis->size()
152 3 : << std::endl;
153 6 : std::cout << fmt::format(">>BAS{:7d}", this->basis->size()) << std::endl;
154 :
155 : // === Save single atom basis ===
156 3 : std::cout << "One-atom Hamiltonian, save single atom basis" << std::endl;
157 :
158 : // initialize uuid generator
159 : boost::uuids::random_generator generator;
160 :
161 : // generate uuid
162 6 : std::string uuid;
163 3 : boost::uuids::uuid u = generator();
164 3 : boost::algorithm::hex(u.begin(), u.end(), std::back_inserter(uuid));
165 :
166 : // save basis
167 6 : fs::path path_basis = fs::temp_directory_path();
168 3 : path_basis /= "basis_one_" + uuid + ".csv";
169 3 : this->basis->save(path_basis.string());
170 :
171 6 : std::cout << fmt::format(">>STA {:s}", path_basis.string()) << std::endl;
172 :
173 : ////////////////////////////////////////////////////////
174 : ////// Construct atom-field interaction ////////////////
175 : ////////////////////////////////////////////////////////
176 :
177 : // --- Obtain existence of fields ---
178 1 : Scalar min_E_0, min_E_p, min_E_m, min_B_0, min_B_p, min_B_m, max_E_0, max_E_p, max_E_m, max_B_0,
179 1 : max_B_p, max_B_m;
180 3 : changeToSpherical(min_E_x, min_E_y, min_E_z, min_E_p, min_E_m, min_E_0);
181 3 : changeToSpherical(max_E_x, max_E_y, max_E_z, max_E_p, max_E_m, max_E_0);
182 3 : changeToSpherical(min_B_x, min_B_y, min_B_z, min_B_p, min_B_m, min_B_0);
183 3 : changeToSpherical(max_B_x, max_B_y, max_B_z, max_B_p, max_B_m, max_B_0);
184 :
185 3 : bool exist_E_0 = (std::abs(min_E_0) != 0 || std::abs(max_E_0) != 0);
186 3 : bool exist_E_1 = (std::abs(min_E_p) != 0 || std::abs(max_E_p) != 0);
187 3 : bool exist_B_0 = (std::abs(min_B_0) != 0 || std::abs(max_B_0) != 0);
188 3 : bool exist_B_1 = (std::abs(min_B_p) != 0 || std::abs(max_B_p) != 0);
189 :
190 : // --- Precalculate matrix elements --- // TODO parallelization
191 3 : std::cout << "One-atom Hamiltonian, precalculate matrix elements" << std::endl;
192 :
193 9 : MatrixElements matrix_elements(basicconf, species, (path_cache / "cache_elements.db").string());
194 :
195 3 : if (exist_E_0) {
196 1 : matrix_elements.precalculateElectricMomentum(this->basis, 0);
197 : }
198 3 : if (exist_E_1) {
199 1 : matrix_elements.precalculateElectricMomentum(this->basis, 1);
200 : }
201 3 : if (exist_E_1) {
202 1 : matrix_elements.precalculateElectricMomentum(this->basis, -1);
203 : }
204 :
205 3 : if (exist_B_0) {
206 1 : matrix_elements.precalculateMagneticMomentum(this->basis, 0);
207 : }
208 3 : if (exist_B_1) {
209 1 : matrix_elements.precalculateMagneticMomentum(this->basis, 1);
210 : }
211 3 : if (exist_B_1) {
212 1 : matrix_elements.precalculateMagneticMomentum(this->basis, -1);
213 : }
214 :
215 3 : if (diamagnetism && (exist_B_0 || exist_B_1)) {
216 1 : matrix_elements.precalculateDiamagnetism(this->basis, 0, 0);
217 : }
218 3 : if (diamagnetism && (exist_B_0 || exist_B_1)) {
219 1 : matrix_elements.precalculateDiamagnetism(this->basis, 2, 0);
220 : }
221 3 : if (diamagnetism && exist_B_0 && exist_B_1) {
222 1 : matrix_elements.precalculateDiamagnetism(this->basis, 2, 1);
223 : }
224 3 : if (diamagnetism && exist_B_0 && exist_B_1) {
225 1 : matrix_elements.precalculateDiamagnetism(this->basis, 2, -1);
226 : }
227 3 : if (diamagnetism && exist_B_1) {
228 1 : matrix_elements.precalculateDiamagnetism(this->basis, 2, 2);
229 : }
230 3 : if (diamagnetism && exist_B_1) {
231 1 : matrix_elements.precalculateDiamagnetism(this->basis, 2, -2);
232 : }
233 :
234 : // --- Count entries of atom-field Hamiltonian ---
235 3 : std::cout << "One-atom Hamiltonian, count number of entries within the field Hamiltonian"
236 3 : << std::endl;
237 :
238 3 : size_basis = this->basis->size();
239 3 : size_t size_electricMomentum_0 = 0;
240 3 : size_t size_electricMomentum_p = 0;
241 3 : size_t size_electricMomentum_m = 0;
242 :
243 3 : size_t size_magneticMomentum_0 = 0;
244 3 : size_t size_magneticMomentum_p = 0;
245 3 : size_t size_magneticMomentum_m = 0;
246 :
247 3 : size_t size_diamagnetism_00 = 0;
248 3 : size_t size_diamagnetism_20 = 0;
249 3 : size_t size_diamagnetism_2p = 0;
250 3 : size_t size_diamagnetism_2m = 0;
251 3 : size_t size_diamagnetism_2pp = 0;
252 3 : size_t size_diamagnetism_2mm = 0;
253 :
254 377 : for (const auto &state_col : *this->basis) { // TODO parallelization
255 73050 : for (const auto &state_row : *this->basis) {
256 72676 : if (state_row.idx < state_col.idx) { // lower triangle only
257 36151 : continue;
258 : }
259 :
260 36525 : if (exist_E_0 && selectionRulesMultipole(state_row, state_col, 1, 0)) {
261 468 : size_electricMomentum_0++;
262 36057 : } else if (exist_E_1 && selectionRulesMultipole(state_row, state_col, 1, 1)) {
263 432 : size_electricMomentum_p++;
264 35625 : } else if (exist_E_1 && selectionRulesMultipole(state_row, state_col, 1, -1)) {
265 432 : size_electricMomentum_m++;
266 : }
267 :
268 36525 : if (exist_B_0 && selectionRulesMomentum(state_row, state_col, 0)) {
269 480 : size_magneticMomentum_0++;
270 36045 : } else if (exist_B_1 && selectionRulesMomentum(state_row, state_col, 1)) {
271 426 : size_magneticMomentum_p++;
272 35619 : } else if (exist_B_1 && selectionRulesMomentum(state_row, state_col, -1)) {
273 303 : size_magneticMomentum_m++;
274 : }
275 :
276 47850 : if (diamagnetism && (exist_B_0 || exist_B_1) &&
277 11325 : selectionRulesMultipole(state_row, state_col, 0, 0)) {
278 300 : size_diamagnetism_00++;
279 47250 : } else if (diamagnetism && (exist_B_0 || exist_B_1) &&
280 11025 : selectionRulesMultipole(state_row, state_col, 2, 0)) {
281 450 : size_diamagnetism_20++;
282 46350 : } else if (diamagnetism && (exist_B_0 && exist_B_1) &&
283 10575 : selectionRulesMultipole(state_row, state_col, 2, 1)) {
284 672 : size_diamagnetism_2p++;
285 45006 : } else if (diamagnetism && (exist_B_0 && exist_B_1) &&
286 9903 : selectionRulesMultipole(state_row, state_col, 2, -1)) {
287 561 : size_diamagnetism_2m++;
288 43884 : } else if (diamagnetism && (exist_B_1) &&
289 9342 : selectionRulesMultipole(state_row, state_col, 2, 2)) {
290 579 : size_diamagnetism_2pp++;
291 42726 : } else if (diamagnetism && (exist_B_1) &&
292 8763 : selectionRulesMultipole(state_row, state_col, 2, -2)) {
293 483 : size_diamagnetism_2mm++;
294 : }
295 : }
296 : }
297 :
298 : // --- Construct atom-field Hamiltonian ---
299 3 : std::cout << "One-atom Hamiltonian, construct field Hamiltonian" << std::endl;
300 :
301 6 : Hamiltonianmatrix<Scalar> hamiltonian_electricMomentum_0(size_basis, size_electricMomentum_0);
302 6 : Hamiltonianmatrix<Scalar> hamiltonian_electricMomentum_p(size_basis, size_electricMomentum_p);
303 6 : Hamiltonianmatrix<Scalar> hamiltonian_electricMomentum_m(size_basis, size_electricMomentum_m);
304 :
305 6 : Hamiltonianmatrix<Scalar> hamiltonian_magneticMomentum_0(size_basis, size_magneticMomentum_0);
306 6 : Hamiltonianmatrix<Scalar> hamiltonian_magneticMomentum_p(size_basis, size_magneticMomentum_p);
307 6 : Hamiltonianmatrix<Scalar> hamiltonian_magneticMomentum_m(size_basis, size_magneticMomentum_m);
308 :
309 6 : Hamiltonianmatrix<Scalar> hamiltonian_diamagnetism_00(size_basis, size_diamagnetism_00);
310 6 : Hamiltonianmatrix<Scalar> hamiltonian_diamagnetism_20(size_basis, size_diamagnetism_20);
311 6 : Hamiltonianmatrix<Scalar> hamiltonian_diamagnetism_2p(size_basis, size_diamagnetism_2p);
312 6 : Hamiltonianmatrix<Scalar> hamiltonian_diamagnetism_2m(size_basis, size_diamagnetism_2m);
313 6 : Hamiltonianmatrix<Scalar> hamiltonian_diamagnetism_2pp(size_basis, size_diamagnetism_2pp);
314 6 : Hamiltonianmatrix<Scalar> hamiltonian_diamagnetism_2mm(size_basis, size_diamagnetism_2mm);
315 :
316 377 : for (const auto &state_col : *this->basis) { // TODO parallelization
317 73050 : for (const auto &state_row : *this->basis) {
318 72676 : if (state_row.idx < state_col.idx) {
319 36151 : continue;
320 : }
321 :
322 36525 : if (state_row.idx == state_col.idx) {
323 374 : hamiltonian_electricMomentum_0.addBasis(state_row.idx, state_col.idx, 1);
324 374 : hamiltonian_electricMomentum_p.addBasis(state_row.idx, state_col.idx, 1);
325 374 : hamiltonian_electricMomentum_m.addBasis(state_row.idx, state_col.idx, 1);
326 :
327 374 : hamiltonian_magneticMomentum_0.addBasis(state_row.idx, state_col.idx, 1);
328 374 : hamiltonian_magneticMomentum_p.addBasis(state_row.idx, state_col.idx, 1);
329 374 : hamiltonian_magneticMomentum_m.addBasis(state_row.idx, state_col.idx, 1);
330 :
331 374 : hamiltonian_diamagnetism_00.addBasis(state_row.idx, state_col.idx, 1);
332 374 : hamiltonian_diamagnetism_20.addBasis(state_row.idx, state_col.idx, 1);
333 374 : hamiltonian_diamagnetism_2p.addBasis(state_row.idx, state_col.idx, 1);
334 374 : hamiltonian_diamagnetism_2m.addBasis(state_row.idx, state_col.idx, 1);
335 374 : hamiltonian_diamagnetism_2pp.addBasis(state_row.idx, state_col.idx, 1);
336 374 : hamiltonian_diamagnetism_2mm.addBasis(state_row.idx, state_col.idx, 1);
337 : }
338 :
339 36525 : if (exist_E_0 && selectionRulesMultipole(state_row, state_col, 1, 0)) {
340 468 : double val = matrix_elements.getElectricMomentum(state_row, state_col);
341 468 : if (std::abs(val) > tol) {
342 468 : hamiltonian_electricMomentum_0.addEntries(state_row.idx, state_col.idx, val);
343 : }
344 36057 : } else if (exist_E_1 && selectionRulesMultipole(state_row, state_col, 1, 1)) {
345 432 : double val = matrix_elements.getElectricMomentum(state_row, state_col);
346 432 : if (std::abs(val) > tol) {
347 432 : hamiltonian_electricMomentum_p.addEntries(state_row.idx, state_col.idx, val);
348 : }
349 35625 : } else if (exist_E_1 && selectionRulesMultipole(state_row, state_col, 1, -1)) {
350 432 : double val = matrix_elements.getElectricMomentum(state_row, state_col);
351 432 : if (std::abs(val) > tol) {
352 432 : hamiltonian_electricMomentum_m.addEntries(state_row.idx, state_col.idx, val);
353 : }
354 : }
355 :
356 36525 : if (exist_B_0 && selectionRulesMomentum(state_row, state_col, 0)) {
357 480 : double val = matrix_elements.getMagneticMomentum(state_row, state_col);
358 480 : if (std::abs(val) > tol) {
359 480 : hamiltonian_magneticMomentum_0.addEntries(state_row.idx, state_col.idx, val);
360 : }
361 36045 : } else if (exist_B_1 && selectionRulesMomentum(state_row, state_col, 1)) {
362 426 : double val = matrix_elements.getMagneticMomentum(state_row, state_col);
363 426 : if (std::abs(val) > tol) {
364 426 : hamiltonian_magneticMomentum_p.addEntries(state_row.idx, state_col.idx, val);
365 : }
366 35619 : } else if (exist_B_1 && selectionRulesMomentum(state_row, state_col, -1)) {
367 303 : double val = matrix_elements.getMagneticMomentum(state_row, state_col);
368 303 : if (std::abs(val) > tol) {
369 303 : hamiltonian_magneticMomentum_m.addEntries(state_row.idx, state_col.idx, val);
370 : }
371 : }
372 :
373 47850 : if (diamagnetism && (exist_B_0 || exist_B_1) &&
374 11325 : selectionRulesMultipole(state_row, state_col, 0, 0)) {
375 300 : double val = matrix_elements.getDiamagnetism(state_row, state_col, 0);
376 300 : if (std::abs(val) > tol) {
377 300 : hamiltonian_diamagnetism_00.addEntries(state_row.idx, state_col.idx, val);
378 : }
379 47250 : } else if (diamagnetism && (exist_B_0 || exist_B_1) &&
380 11025 : selectionRulesMultipole(state_row, state_col, 2, 0)) {
381 450 : double val = matrix_elements.getDiamagnetism(state_row, state_col, 2);
382 450 : if (std::abs(val) > tol) {
383 450 : hamiltonian_diamagnetism_20.addEntries(state_row.idx, state_col.idx, val);
384 : }
385 46350 : } else if (diamagnetism && (exist_B_0 && exist_B_1) &&
386 10575 : selectionRulesMultipole(state_row, state_col, 2, 1)) {
387 672 : double val = matrix_elements.getDiamagnetism(state_row, state_col, 2);
388 672 : if (std::abs(val) > tol) {
389 660 : hamiltonian_diamagnetism_2p.addEntries(state_row.idx, state_col.idx, val);
390 : }
391 45006 : } else if (diamagnetism && (exist_B_0 && exist_B_1) &&
392 9903 : selectionRulesMultipole(state_row, state_col, 2, -1)) {
393 561 : double val = matrix_elements.getDiamagnetism(state_row, state_col, 2);
394 561 : if (std::abs(val) > tol) {
395 555 : hamiltonian_diamagnetism_2m.addEntries(state_row.idx, state_col.idx, val);
396 : }
397 43884 : } else if (diamagnetism && (exist_B_1) &&
398 9342 : selectionRulesMultipole(state_row, state_col, 2, 2)) {
399 579 : double val = matrix_elements.getDiamagnetism(state_row, state_col, 2);
400 579 : if (std::abs(val) > tol) {
401 579 : hamiltonian_diamagnetism_2pp.addEntries(state_row.idx, state_col.idx, val);
402 : }
403 42726 : } else if (diamagnetism && (exist_B_1) &&
404 8763 : selectionRulesMultipole(state_row, state_col, 2, -2)) {
405 483 : double val = matrix_elements.getDiamagnetism(state_row, state_col, 2);
406 483 : if (std::abs(val) > tol) {
407 483 : hamiltonian_diamagnetism_2mm.addEntries(state_row.idx, state_col.idx, val);
408 : }
409 : }
410 : }
411 : }
412 :
413 3 : std::cout << "One-atom Hamiltonian, compress field Hamiltonian" << std::endl;
414 :
415 3 : hamiltonian_electricMomentum_0.compress(this->basis->dim(), this->basis->dim());
416 3 : hamiltonian_electricMomentum_p.compress(this->basis->dim(), this->basis->dim());
417 3 : hamiltonian_electricMomentum_m.compress(this->basis->dim(), this->basis->dim());
418 :
419 3 : hamiltonian_magneticMomentum_0.compress(this->basis->dim(), this->basis->dim());
420 3 : hamiltonian_magneticMomentum_p.compress(this->basis->dim(), this->basis->dim());
421 3 : hamiltonian_magneticMomentum_m.compress(this->basis->dim(), this->basis->dim());
422 :
423 3 : hamiltonian_diamagnetism_00.compress(this->basis->dim(), this->basis->dim());
424 3 : hamiltonian_diamagnetism_20.compress(this->basis->dim(), this->basis->dim());
425 3 : hamiltonian_diamagnetism_2p.compress(this->basis->dim(), this->basis->dim());
426 3 : hamiltonian_diamagnetism_2m.compress(this->basis->dim(), this->basis->dim());
427 3 : hamiltonian_diamagnetism_2pp.compress(this->basis->dim(), this->basis->dim());
428 3 : hamiltonian_diamagnetism_2mm.compress(this->basis->dim(), this->basis->dim());
429 :
430 : ////////////////////////////////////////////////////////
431 : ////// Prepare processing of Hamiltonians //////////////
432 : ////////////////////////////////////////////////////////
433 :
434 : // TODO Put the logic in its own class
435 :
436 3 : std::cout << "One-atom Hamiltonian, processe Hamiltonians" << std::endl;
437 :
438 : // === Open database ===
439 6 : fs::path path_db;
440 :
441 : if (utils::is_complex<Scalar>::value) {
442 1 : path_db = path_cache / "cache_matrix_complex.db";
443 : } else {
444 2 : path_db = path_cache / "cache_matrix_real.db";
445 : }
446 6 : sqlite::handle db(path_db.string());
447 6 : sqlite::statement stmt(db);
448 :
449 : // === Initialize variables ===
450 3 : bool flag_perhapsmissingtable = true;
451 :
452 3 : this->matrix_path.resize(nSteps);
453 3 : this->matrix_diag.resize(nSteps); // TODO maybe remove
454 3 : this->params.resize(nSteps); // TODO maybe remove
455 :
456 : ////////////////////////////////////////////////////////
457 : ////// Loop through steps //////////////////////////////
458 : ////////////////////////////////////////////////////////
459 :
460 6 : std::cout << fmt::format(">>TOT{:7d}", nSteps) << std::endl;
461 :
462 3 : auto nSteps_i = static_cast<int>(nSteps);
463 :
464 3 : #pragma omp parallel for schedule(static, 1)
465 :
466 : // Loop through steps
467 : for (int step = 0; step < nSteps_i; ++step) {
468 :
469 : // === Get parameters for the current position inside the loop ===
470 :
471 : // Get fields
472 : double normalized_position = (nSteps > 1) ? step / (nSteps - 1.) : 0;
473 :
474 : double Ex = min_E_x + normalized_position * (max_E_x - min_E_x);
475 : double Ey = min_E_y + normalized_position * (max_E_y - min_E_y);
476 : double Ez = min_E_z + normalized_position * (max_E_z - min_E_z);
477 : double Bx = min_B_x + normalized_position * (max_B_x - min_B_x);
478 : double By = min_B_y + normalized_position * (max_B_y - min_B_y);
479 : double Bz = min_B_z + normalized_position * (max_B_z - min_B_z);
480 :
481 : Scalar E_0 = min_E_0 + normalized_position * (max_E_0 - min_E_0);
482 : Scalar E_p = min_E_p + normalized_position * (max_E_p - min_E_p);
483 : Scalar E_m = min_E_m + normalized_position * (max_E_m - min_E_m);
484 : Scalar B_0 = min_B_0 + normalized_position * (max_B_0 - min_B_0);
485 : Scalar B_p = min_B_p + normalized_position * (max_B_p - min_B_p);
486 : Scalar B_m = min_B_m + normalized_position * (max_B_m - min_B_m);
487 :
488 : // Get configuration and save fields
489 : Configuration conf = basicconf;
490 : conf["Ex"] << Ex;
491 : conf["Ey"] << Ey;
492 : conf["Ez"] << Ez;
493 : conf["Bx"] << Bx;
494 : conf["By"] << By;
495 : conf["Bz"] << Bz;
496 :
497 : // === Create table if necessary ===
498 : std::stringstream query;
499 : std::string spacer;
500 :
501 : if (flag_perhapsmissingtable) {
502 : query << "CREATE TABLE IF NOT EXISTS cache_one (uuid text NOT NULL PRIMARY KEY, "
503 : "created TIMESTAMP DEFAULT CURRENT_TIMESTAMP, "
504 : "accessed TIMESTAMP DEFAULT CURRENT_TIMESTAMP";
505 : for (auto p : conf) {
506 : query << ", " << p.first << " text";
507 : }
508 : query << ", UNIQUE (";
509 : for (auto p : conf) {
510 : query << spacer << p.first;
511 : spacer = ", ";
512 : }
513 : query << "));";
514 :
515 : flag_perhapsmissingtable = false;
516 :
517 : #pragma omp critical(database)
518 : stmt.exec(query.str());
519 : }
520 :
521 : // === Get uuid as filename === // TODO put code in its own method
522 : std::string uuid;
523 : spacer = "";
524 : query.str(std::string());
525 : query << "SELECT uuid FROM cache_one WHERE ";
526 : for (auto p : conf) {
527 : query << spacer << p.first << "='" << p.second.str() << "'";
528 : spacer = " AND ";
529 : }
530 : query << ";";
531 :
532 : #pragma omp critical(database)
533 : {
534 : sqlite::statement stmt(db, query.str());
535 : stmt.prepare();
536 : if (stmt.step()) {
537 : uuid = stmt.get<std::string>(0);
538 : }
539 : }
540 :
541 : if (!uuid.empty()) {
542 : query.str(std::string());
543 : query << "UPDATE cache_one SET accessed = CURRENT_TIMESTAMP WHERE uuid = '" << uuid
544 : << "';";
545 : #pragma omp critical(database)
546 : stmt.exec(query.str()); // TODO check whether this slows down the program
547 :
548 : } else {
549 : boost::uuids::uuid u = generator();
550 : boost::algorithm::hex(u.begin(), u.end(), std::back_inserter(uuid));
551 :
552 : query.str(std::string());
553 : query << "INSERT INTO cache_one (uuid";
554 : for (auto p : conf) {
555 : query << ", " << p.first;
556 : }
557 : query << ") values ( '" << uuid << "'";
558 : for (auto p : conf) {
559 : query << ", "
560 : << "'" << p.second.str() << "'";
561 : }
562 : query << ");";
563 : #pragma omp critical(database)
564 : stmt.exec(query.str());
565 : }
566 :
567 : // === Check existence of files === // TODO put code in its own method
568 :
569 : // Check whether .mat and .json file exists and compare settings in program with settings in
570 : // .json file
571 : fs::path path, path_mat, path_json;
572 :
573 : path = path_cache_mat / ("one_" + uuid);
574 : path_mat = path;
575 : path_mat.replace_extension(".mat");
576 : path_json = path;
577 : path_json.replace_extension(".json");
578 :
579 : bool is_existing = false;
580 : if (fs::exists(path_mat)) {
581 : if (fs::exists(path_json)) {
582 : Configuration params_loaded;
583 : params_loaded.load_from_json(path_json.string());
584 : if (conf == params_loaded) {
585 : is_existing = true;
586 : }
587 : }
588 : }
589 :
590 : // Create .json file if "is_existing" is false
591 : if (!is_existing) {
592 : conf.save_to_json(path_json.string());
593 : }
594 :
595 : // === Build and diagonalize total matrix if not existent ===
596 : Hamiltonianmatrix<Scalar> totalmatrix;
597 :
598 : // calculate Hamiltonian if "is_existing" is false
599 : std::shared_ptr<Hamiltonianmatrix<Scalar>> mat;
600 : if (!is_existing || !totalmatrix.load(path_mat.string())) {
601 :
602 : // --- Build matrix ---
603 : totalmatrix = hamiltonian_energy - hamiltonian_electricMomentum_0 * E_0 +
604 : hamiltonian_electricMomentum_p * E_m + hamiltonian_electricMomentum_m * E_p +
605 : hamiltonian_magneticMomentum_0 * B_0 - hamiltonian_magneticMomentum_p * B_m -
606 : hamiltonian_magneticMomentum_m * B_p +
607 : hamiltonian_diamagnetism_00 * (B_0 * B_0 - 2. * B_p * B_m) -
608 : hamiltonian_diamagnetism_20 * (B_0 * B_0 + B_p * B_m) +
609 : std::sqrt(3) * hamiltonian_diamagnetism_2p * B_0 * B_m +
610 : std::sqrt(3) * hamiltonian_diamagnetism_2m * B_0 * B_p -
611 : std::sqrt(1.5) * hamiltonian_diamagnetism_2pp * B_m * B_m -
612 : std::sqrt(1.5) * hamiltonian_diamagnetism_2mm * B_p * B_p;
613 :
614 : // Stdout: Hamiltonian assembled
615 : #pragma omp critical(textoutput)
616 : {
617 : std::cout << fmt::format(">>DIM{:7d}", totalmatrix.num_basisvectors()) << std::endl;
618 : std::cout << "One-atom Hamiltonian, " << step + 1 << ". Hamiltonian assembled"
619 : << std::endl;
620 : }
621 :
622 : // --- Diagonalize matrix and save diagonalized matrix ---
623 : totalmatrix.diagonalize();
624 : totalmatrix.save(path_mat.string());
625 :
626 : // Stdout: Hamiltonian diagonalized
627 : #pragma omp critical(textoutput)
628 : {
629 : std::cout << fmt::format(">>OUT{:7d}{:7d}{:7d}{:7d} {:s}", (step + 1), step, 1, 0,
630 : path.string())
631 : << std::endl;
632 : std::cout << "One-atom Hamiltonian, " << step + 1 << ". Hamiltonian diagonalized"
633 : << std::endl;
634 : }
635 : } else {
636 : // Stdout: Hamiltonian loaded
637 : #pragma omp critical(textoutput)
638 : {
639 : std::cout << fmt::format(">>DIM{:7d}", totalmatrix.num_basisvectors()) << std::endl;
640 : std::cout << fmt::format(">>OUT{:7d}{:7d}{:7d}{:7d} {:s}", (step + 1), step, 1, 0,
641 : path.string())
642 : << std::endl;
643 : std::cout << "One-atom Hamiltonian, " << step + 1 << ". Hamiltonian loaded"
644 : << std::endl;
645 : }
646 : }
647 :
648 : // === Store path to configuration and diagonalized matrix ===
649 : this->matrix_path[step] = path.string();
650 : this->matrix_diag[step] =
651 : std::make_shared<Hamiltonianmatrix<Scalar>>(totalmatrix); // TODO maybe remove
652 : this->params[step] = std::make_shared<Configuration>(conf); // TODO maybe remove
653 : }
654 :
655 3 : std::cout << "One-atom Hamiltonian, all Hamiltonians processed" << std::endl;
656 3 : }
657 :
658 : template class HamiltonianOne<std::complex<double>>;
659 : template class HamiltonianOne<double>;
|
