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 "SystemOne.hpp"
21 : #include "Constants.hpp"
22 : #include "MatrixElementCache.hpp"
23 : #include "Symmetry.hpp"
24 :
25 : #include <cctype>
26 : #include <cmath>
27 : #include <limits>
28 : #include <numeric>
29 : #include <string>
30 : #include <type_traits>
31 : #include <unordered_set>
32 : #include <utility>
33 : #include <vector>
34 :
35 : template <typename Scalar>
36 0 : SystemOne<Scalar>::SystemOne()
37 0 : : efield({{0, 0, 0}}), bfield({{0, 0, 0}}), sym_rotation({static_cast<float>(ARB)}) {}
38 :
39 : template <typename Scalar>
40 59 : SystemOne<Scalar>::SystemOne(std::string species, MatrixElementCache &cache)
41 : : SystemBase<Scalar, StateOne>(cache), efield({{0, 0, 0}}), bfield({{0, 0, 0}}),
42 59 : diamagnetism(true), charge(0), ordermax(0), distance(std::numeric_limits<double>::max()),
43 59 : species(std::move(species)), sym_reflection(NA), sym_rotation({static_cast<float>(ARB)}) {}
44 :
45 : template <typename Scalar>
46 0 : SystemOne<Scalar>::SystemOne(std::string species, MatrixElementCache &cache, bool memory_saving)
47 : : SystemBase<Scalar, StateOne>(cache, memory_saving), efield({{0, 0, 0}}), bfield({{0, 0, 0}}),
48 0 : diamagnetism(true), charge(0), ordermax(0), distance(std::numeric_limits<double>::max()),
49 0 : species(std::move(species)), sym_reflection(NA), sym_rotation({static_cast<float>(ARB)}) {}
50 :
51 : template <typename Scalar>
52 124 : const std::string &SystemOne<Scalar>::getSpecies() const {
53 124 : return species;
54 : }
55 :
56 : template <typename Scalar>
57 122 : void SystemOne<Scalar>::setEfield(std::array<double, 3> field) {
58 122 : this->onParameterChange();
59 122 : efield = field;
60 :
61 : // Transform the electric field into spherical coordinates
62 122 : this->changeToSphericalbasis(efield, efield_spherical);
63 122 : }
64 :
65 : template <typename Scalar>
66 24 : void SystemOne<Scalar>::setBfield(std::array<double, 3> field) {
67 24 : this->onParameterChange();
68 24 : bfield = field;
69 :
70 : // Transform the magnetic field into spherical coordinates
71 24 : this->changeToSphericalbasis(bfield, bfield_spherical);
72 :
73 24 : diamagnetism_terms[{{0, +0}}] = bfield_spherical[+0] * bfield_spherical[+0] -
74 24 : bfield_spherical[+1] * bfield_spherical[-1] * 2.;
75 20 : diamagnetism_terms[{{2, +0}}] =
76 24 : bfield_spherical[+0] * bfield_spherical[+0] + bfield_spherical[+1] * bfield_spherical[-1];
77 24 : diamagnetism_terms[{{2, +1}}] = bfield_spherical[+0] * bfield_spherical[-1];
78 24 : diamagnetism_terms[{{2, -1}}] = bfield_spherical[+0] * bfield_spherical[+1];
79 24 : diamagnetism_terms[{{2, +2}}] = bfield_spherical[-1] * bfield_spherical[-1];
80 24 : diamagnetism_terms[{{2, -2}}] = bfield_spherical[+1] * bfield_spherical[+1];
81 24 : }
82 :
83 : template <typename Scalar>
84 0 : void SystemOne<Scalar>::setEfield(std::array<double, 3> field, std::array<double, 3> to_z_axis,
85 : std::array<double, 3> to_y_axis) {
86 0 : this->rotateVector(field, to_z_axis, to_y_axis);
87 0 : this->setEfield(field);
88 0 : }
89 :
90 : template <typename Scalar>
91 0 : void SystemOne<Scalar>::setBfield(std::array<double, 3> field, std::array<double, 3> to_z_axis,
92 : std::array<double, 3> to_y_axis) {
93 0 : this->rotateVector(field, to_z_axis, to_y_axis);
94 0 : this->setBfield(field);
95 0 : }
96 :
97 : template <typename Scalar>
98 1 : void SystemOne<Scalar>::setEfield(std::array<double, 3> field, double alpha, double beta,
99 : double gamma) {
100 1 : this->rotateVector(field, alpha, beta, gamma);
101 1 : this->setEfield(field);
102 1 : }
103 :
104 : template <typename Scalar>
105 2 : void SystemOne<Scalar>::setBfield(std::array<double, 3> field, double alpha, double beta,
106 : double gamma) {
107 2 : this->rotateVector(field, alpha, beta, gamma);
108 2 : this->setBfield(field);
109 2 : }
110 :
111 : template <typename Scalar>
112 11 : void SystemOne<Scalar>::enableDiamagnetism(bool enable) {
113 11 : this->onParameterChange();
114 11 : diamagnetism = enable;
115 11 : }
116 : template <typename Scalar>
117 1 : void SystemOne<Scalar>::setIonCharge(int c) {
118 1 : this->onParameterChange();
119 1 : charge = c;
120 1 : }
121 :
122 : template <typename Scalar>
123 1 : void SystemOne<Scalar>::setRydIonOrder(unsigned int o) {
124 1 : this->onParameterChange();
125 1 : ordermax = o;
126 1 : }
127 :
128 : template <typename Scalar>
129 1 : void SystemOne<Scalar>::setRydIonDistance(double d) {
130 1 : this->onParameterChange();
131 1 : distance = d;
132 1 : }
133 :
134 : template <typename Scalar>
135 6 : void SystemOne<Scalar>::setConservedParityUnderReflection(parity_t parity) {
136 6 : this->onSymmetryChange();
137 6 : sym_reflection = parity;
138 6 : if (!this->isRefelectionAndRotationCompatible()) {
139 0 : throw std::runtime_error("The conserved parity under reflection is not compatible to the "
140 : "previously specified conserved momenta.");
141 : }
142 6 : }
143 :
144 : template <typename Scalar>
145 9 : void SystemOne<Scalar>::setConservedMomentaUnderRotation(const std::set<float> &momenta) {
146 9 : if (momenta.count(static_cast<float>(ARB)) != 0 && momenta.size() > 1) {
147 0 : throw std::runtime_error(
148 : "If ARB (=arbitrary momentum) is specified, momenta must not be passed explicitely.");
149 : }
150 9 : this->onSymmetryChange();
151 9 : sym_rotation = momenta;
152 9 : if (!this->isRefelectionAndRotationCompatible()) {
153 0 : throw std::runtime_error("The conserved momenta are not compatible to the previously "
154 : "specified conserved parity under reflection.");
155 : }
156 9 : }
157 :
158 : ////////////////////////////////////////////////////////////////////
159 : /// Method that allows base class to initialize Basis //////////////
160 : ////////////////////////////////////////////////////////////////////
161 :
162 : template <typename Scalar>
163 89 : void SystemOne<Scalar>::initializeBasis() {
164 : // If the basis is infinite, throw an error
165 89 : if (this->range_n.empty() &&
166 0 : (this->energy_min == std::numeric_limits<double>::lowest() ||
167 0 : this->energy_max == std::numeric_limits<double>::max())) {
168 0 : throw std::runtime_error(
169 : "The number of basis elements is infinite. The basis has to be restricted.");
170 : }
171 :
172 89 : auto &cache = *this->m_cache;
173 :
174 : ////////////////////////////////////////////////////////////////////
175 : /// Build one atom states //////////////////////////////////////////
176 : ////////////////////////////////////////////////////////////////////
177 :
178 : // TODO check whether symmetries are applicable
179 : // TODO check whether range_j, range_m is half-integer or integer valued
180 :
181 89 : float s = 0.5;
182 89 : if (std::isdigit(species.back()) != 0) {
183 0 : s = ((species.back() - '0') - 1) / 2.; // TODO think of a better solution
184 : }
185 :
186 89 : size_t idx = 0;
187 : std::vector<Eigen::Triplet<Scalar>>
188 178 : basisvectors_triplets; // TODO reserve states, basisvectors_triplets,
189 : // hamiltonian_triplets
190 :
191 178 : std::vector<Eigen::Triplet<Scalar>> hamiltonian_triplets;
192 :
193 : /// Loop over specified quantum numbers ////////////////////////////
194 :
195 178 : std::set<int> range_adapted_n, range_adapted_l;
196 178 : std::set<float> range_adapted_j, range_adapted_m;
197 :
198 89 : if (this->range_n.empty()) {
199 0 : throw std::runtime_error(
200 : "The calculation of range_n via energy restrictions is not yet implemented."); // TODO
201 : }
202 89 : range_adapted_n = this->range_n;
203 :
204 482 : for (auto n : range_adapted_n) {
205 :
206 393 : if (this->range_l.empty()) {
207 7 : this->range(range_adapted_l, 0, n - 1);
208 : } else {
209 386 : range_adapted_l = this->range_l;
210 : }
211 2238 : for (auto l : range_adapted_l) {
212 1845 : if (l > n - 1 || l < 0) {
213 348 : continue;
214 : }
215 :
216 1497 : if (this->range_j.empty()) {
217 1281 : this->range(range_adapted_j, std::fabs(l - s), l + s);
218 : } else {
219 216 : range_adapted_j = this->range_j;
220 : }
221 4434 : for (auto j : range_adapted_j) {
222 2937 : if (std::fabs(j - l) > s || j < 0) {
223 1710 : continue;
224 : }
225 :
226 2577 : double energy = StateOne(species, n, l, j, s).getEnergy(cache);
227 2577 : if (!this->checkIsEnergyValid(energy)) {
228 1350 : continue;
229 : }
230 :
231 1227 : if (this->range_m.empty()) {
232 1048 : this->range(range_adapted_m, -j, j);
233 : } else {
234 179 : range_adapted_m = this->range_m;
235 : }
236 :
237 : // Consider rotation symmetry
238 2454 : std::set<float> range_allowed_m;
239 1227 : if (sym_rotation.count(static_cast<float>(ARB)) == 0) {
240 126 : std::set_intersection(sym_rotation.begin(), sym_rotation.end(),
241 : range_adapted_m.begin(), range_adapted_m.end(),
242 : std::inserter(range_allowed_m, range_allowed_m.begin()));
243 : } else {
244 1101 : range_allowed_m = range_adapted_m;
245 : }
246 :
247 5157 : for (auto m : range_allowed_m) {
248 3930 : if (std::fabs(m) > j) {
249 0 : continue;
250 : }
251 :
252 : // Create state
253 7860 : StateOne state(species, n, l, j, m);
254 :
255 : // Check whether reflection symmetry can be realized with the states available
256 3986 : if (sym_reflection != NA && state.getM() != 0 &&
257 3986 : range_allowed_m.count(-state.getM()) == 0) {
258 0 : throw std::runtime_error("The momentum " + std::to_string(-state.getM()) +
259 : " required by symmetries cannot be found.");
260 : }
261 :
262 : // Add symmetrized basis vectors
263 3930 : this->addSymmetrizedBasisvectors(state, idx, energy, basisvectors_triplets,
264 3930 : hamiltonian_triplets, sym_reflection);
265 : }
266 : }
267 : }
268 : }
269 :
270 : /// Loop over user-defined states //////////////////////////////////
271 :
272 : // Check that the user-defined states are not already contained in the list of states
273 89 : for (const auto &state : this->states_to_add) {
274 0 : if (this->states.template get<1>().find(state) != this->states.template get<1>().end()) {
275 0 : throw std::runtime_error("The state " + state.str() +
276 : " is already contained in the list of states.");
277 : }
278 0 : if (!state.isArtificial() && state.getSpecies() != species) {
279 0 : throw std::runtime_error("The state " + state.str() + " is of the wrong species.");
280 : }
281 : }
282 :
283 : // Add user-defined states
284 89 : for (const auto &state : this->states_to_add) {
285 : // Get energy of the state
286 0 : double energy = state.isArtificial() ? 0 : state.getEnergy(cache);
287 :
288 : // In case of artificial states, symmetries won't work
289 0 : auto sym_reflection_local = sym_reflection;
290 0 : auto sym_rotation_local = sym_rotation;
291 0 : if (state.isArtificial()) {
292 0 : if (sym_reflection_local != NA ||
293 0 : sym_rotation_local.count(static_cast<float>(ARB)) == 0) {
294 : std::cerr
295 0 : << "WARNING: Only permutation symmetry can be applied to artificial states."
296 0 : << std::endl;
297 : }
298 0 : sym_reflection_local = NA;
299 0 : sym_rotation_local = std::set<float>({static_cast<float>(ARB)});
300 : }
301 :
302 : // Consider rotation symmetry
303 0 : if (sym_rotation_local.count(static_cast<float>(ARB)) == 0 &&
304 0 : sym_rotation_local.count(state.getM()) == 0) {
305 0 : continue;
306 : }
307 :
308 : // Check whether reflection symmetry can be realized with the states available
309 0 : if (sym_reflection_local != NA && state.getM() != 0) {
310 0 : auto state_reflected = state.getReflected();
311 0 : if (this->states_to_add.find(state_reflected) == this->states_to_add.end()) {
312 0 : throw std::runtime_error("The state " + state_reflected.str() +
313 : " required by symmetries cannot be found.");
314 : }
315 : }
316 :
317 : // Add symmetrized basis vectors
318 0 : this->addSymmetrizedBasisvectors(state, idx, energy, basisvectors_triplets,
319 : hamiltonian_triplets, sym_reflection_local);
320 : }
321 :
322 : /// Build data /////////////////////////////////////////////////////
323 :
324 89 : this->basisvectors.resize(this->states.size(), idx);
325 89 : this->basisvectors.setFromTriplets(basisvectors_triplets.begin(), basisvectors_triplets.end());
326 89 : basisvectors_triplets.clear();
327 :
328 89 : this->hamiltonian.resize(idx, idx);
329 89 : this->hamiltonian.setFromTriplets(hamiltonian_triplets.begin(), hamiltonian_triplets.end());
330 89 : hamiltonian_triplets.clear();
331 89 : }
332 :
333 : ////////////////////////////////////////////////////////////////////
334 : /// Method that allows base class to calculate the interaction /////
335 : ////////////////////////////////////////////////////////////////////
336 :
337 : template <typename Scalar>
338 157 : void SystemOne<Scalar>::initializeInteraction() {
339 : ////////////////////////////////////////////////////////////////////
340 : /// Prepare the calculation of the interaction /////////////////////
341 : ////////////////////////////////////////////////////////////////////
342 :
343 157 : auto &cache = *this->m_cache;
344 :
345 : // Check if something to do
346 157 : double tolerance = 1e-24;
347 :
348 157 : std::vector<int> erange, brange;
349 157 : std::vector<std::array<int, 2>> drange;
350 157 : std::vector<int> orange;
351 577 : for (const auto &entry : efield_spherical) {
352 420 : if (entry.first < 0) {
353 140 : continue;
354 : }
355 433 : if (std::abs(entry.second) > tolerance &&
356 433 : interaction_efield.find(-entry.first) == interaction_efield.end()) {
357 51 : erange.push_back(entry.first);
358 : }
359 : }
360 589 : for (const auto &entry : bfield_spherical) {
361 432 : if (entry.first < 0) {
362 144 : continue;
363 : }
364 351 : if (std::abs(entry.second) > tolerance &&
365 351 : interaction_bfield.find(-entry.first) == interaction_bfield.end()) {
366 60 : brange.push_back(entry.first);
367 : }
368 : }
369 1021 : for (const auto &entry : diamagnetism_terms) {
370 864 : if (entry.first[1] < 0) {
371 288 : continue;
372 : }
373 708 : if (diamagnetism && std::abs(entry.second) > tolerance &&
374 708 : interaction_diamagnetism.find(entry.first) == interaction_diamagnetism.end()) {
375 123 : drange.push_back(entry.first);
376 : }
377 : }
378 :
379 157 : if (charge != 0) {
380 2 : for (unsigned int order = 1; order <= ordermax; ++order) {
381 1 : if (interaction_multipole.find(order) == interaction_multipole.end()) {
382 1 : orange.push_back(order);
383 : }
384 : }
385 : }
386 : // Return if there is nothing to do
387 157 : if (erange.empty() && brange.empty() && drange.empty() && orange.empty()) {
388 107 : return;
389 : }
390 :
391 : // Precalculate matrix elements
392 100 : auto states_converted = this->getStates();
393 101 : for (const auto &i : erange) {
394 51 : cache.precalculateElectricMomentum(states_converted, i);
395 51 : if (i != 0) {
396 19 : cache.precalculateElectricMomentum(states_converted, -i);
397 : }
398 : }
399 110 : for (const auto &i : brange) {
400 60 : cache.precalculateMagneticMomentum(states_converted, i);
401 60 : if (i != 0) {
402 30 : cache.precalculateMagneticMomentum(states_converted, -i);
403 : }
404 : }
405 173 : for (const auto &i : drange) {
406 123 : cache.precalculateDiamagnetism(states_converted, i[0], i[1]);
407 123 : if (i[1] != 0) {
408 47 : cache.precalculateDiamagnetism(states_converted, i[0], -i[1]);
409 : }
410 : }
411 50 : if (charge != 0) {
412 2 : for (unsigned int order = 1; order <= ordermax; ++order) {
413 1 : cache.precalculateMultipole(states_converted, order);
414 : }
415 : }
416 :
417 : ////////////////////////////////////////////////////////////////////
418 : /// Calculate the interaction in the canonical basis ///////////////
419 : ////////////////////////////////////////////////////////////////////
420 :
421 : std::unordered_map<int, std::vector<Eigen::Triplet<Scalar>>>
422 100 : interaction_efield_triplets; // TODO reserve
423 : std::unordered_map<int, std::vector<Eigen::Triplet<Scalar>>>
424 100 : interaction_bfield_triplets; // TODO reserve
425 : std::unordered_map<std::array<int, 2>, std::vector<Eigen::Triplet<Scalar>>,
426 : utils::hash<std::array<int, 2>>>
427 100 : interaction_diamagnetism_triplets; // TODO reserve
428 : std::unordered_map<int, std::vector<Eigen::Triplet<Scalar>>>
429 100 : interaction_multipole_triplets; // TODO reserve
430 : // Loop over column entries
431 3042 : for (const auto &c : this->states) { // TODO parallelization
432 1496 : if (c.state.isArtificial()) {
433 0 : continue;
434 : }
435 :
436 : // Loop over row entries
437 211800 : for (const auto &r : this->states) {
438 105152 : if (r.state.isArtificial()) {
439 0 : continue;
440 : }
441 :
442 : // E-field interaction
443 220416 : for (const auto &i : erange) {
444 119574 : if (i == 0 && r.idx < c.idx) {
445 34288 : continue;
446 : }
447 :
448 85286 : if (selectionRulesMultipoleNew(r.state, c.state, 1, i)) {
449 4310 : Scalar value = cache.getElectricDipole(r.state, c.state);
450 4310 : this->addTriplet(interaction_efield_triplets[i], r.idx, c.idx, value);
451 4310 : break; // because for the other operators, the selection rule for the magnetic
452 : // quantum numbers will not be fulfilled
453 : }
454 : }
455 :
456 : // B-field interaction
457 219095 : for (const auto &i : brange) {
458 118353 : if (i == 0 && r.idx < c.idx) {
459 31195 : continue;
460 : }
461 :
462 87158 : if (selectionRulesMomentumNew(r.state, c.state, i)) {
463 4410 : Scalar value = cache.getMagneticDipole(r.state, c.state);
464 4410 : this->addTriplet(interaction_bfield_triplets[i], r.idx, c.idx, value);
465 4410 : break; // because for the other operators, the selection rule for the magnetic
466 : // quantum numbers will not be fulfilled
467 : }
468 : }
469 :
470 : // Diamagnetic interaction
471 329386 : for (const auto &i : drange) {
472 224234 : if (i[1] == 0 && r.idx < c.idx) {
473 57962 : continue;
474 : }
475 :
476 166272 : if (selectionRulesMultipoleNew(r.state, c.state, i[0], i[1])) {
477 9400 : Scalar value = 1. / (8 * electron_rest_mass) *
478 9400 : cache.getDiamagnetism(r.state, c.state, i[0]);
479 9400 : this->addTriplet(interaction_diamagnetism_triplets[i], r.idx, c.idx, value);
480 : }
481 : }
482 :
483 : // Multipole interaction with an ion
484 105152 : if (charge != 0) {
485 28900 : int q = r.state.getM() - c.state.getM();
486 28900 : if (q == 0) { // total momentum consreved
487 57800 : for (const auto &order : orange) {
488 28900 : if (selectionRulesMultipoleNew(r.state, c.state, order)) {
489 1008 : double val = -coulombs_constant * elementary_charge *
490 1008 : cache.getElectricMultipole(r.state, c.state, order);
491 1008 : this->addTriplet(interaction_multipole_triplets[order], r.idx, c.idx,
492 : val);
493 : }
494 : }
495 : }
496 : }
497 : }
498 : }
499 : ////////////////////////////////////////////////////////////////////
500 : /// Build and transform the interaction to the used basis //////////
501 : ////////////////////////////////////////////////////////////////////
502 :
503 101 : for (const auto &i : erange) {
504 51 : interaction_efield[i].resize(this->states.size(), this->states.size());
505 102 : interaction_efield[i].setFromTriplets(interaction_efield_triplets[i].begin(),
506 51 : interaction_efield_triplets[i].end());
507 51 : interaction_efield_triplets[i].clear();
508 :
509 51 : if (i == 0) {
510 64 : interaction_efield[i] = this->basisvectors.adjoint() *
511 96 : interaction_efield[i].template selfadjointView<Eigen::Lower>() * this->basisvectors;
512 : } else {
513 19 : interaction_efield[i] =
514 19 : this->basisvectors.adjoint() * interaction_efield[i] * this->basisvectors;
515 19 : interaction_efield[-i] = std::pow(-1, i) * interaction_efield[i].adjoint();
516 : }
517 : }
518 :
519 110 : for (const auto &i : brange) {
520 60 : interaction_bfield[i].resize(this->states.size(), this->states.size());
521 120 : interaction_bfield[i].setFromTriplets(interaction_bfield_triplets[i].begin(),
522 60 : interaction_bfield_triplets[i].end());
523 60 : interaction_bfield_triplets[i].clear();
524 :
525 60 : if (i == 0) {
526 60 : interaction_bfield[i] = this->basisvectors.adjoint() *
527 90 : interaction_bfield[i].template selfadjointView<Eigen::Lower>() * this->basisvectors;
528 : } else {
529 30 : interaction_bfield[i] =
530 30 : this->basisvectors.adjoint() * interaction_bfield[i] * this->basisvectors;
531 30 : interaction_bfield[-i] = std::pow(-1, i) * interaction_bfield[i].adjoint();
532 : }
533 : }
534 :
535 173 : for (const auto &i : drange) {
536 123 : interaction_diamagnetism[i].resize(this->states.size(), this->states.size());
537 246 : interaction_diamagnetism[i].setFromTriplets(interaction_diamagnetism_triplets[i].begin(),
538 123 : interaction_diamagnetism_triplets[i].end());
539 123 : interaction_diamagnetism_triplets[i].clear();
540 :
541 123 : if (i[1] == 0) {
542 152 : interaction_diamagnetism[i] = this->basisvectors.adjoint() *
543 228 : interaction_diamagnetism[i].template selfadjointView<Eigen::Lower>() *
544 : this->basisvectors;
545 : } else {
546 47 : interaction_diamagnetism[i] =
547 47 : this->basisvectors.adjoint() * interaction_diamagnetism[i] * this->basisvectors;
548 47 : interaction_diamagnetism[{{i[0], -i[1]}}] =
549 94 : std::pow(-1, i[1]) * interaction_diamagnetism[i].adjoint();
550 : }
551 : }
552 50 : if (charge != 0) {
553 2 : for (const auto &i : orange) {
554 1 : interaction_multipole[i].resize(this->states.size(), this->states.size());
555 2 : interaction_multipole[i].setFromTriplets(interaction_multipole_triplets[i].begin(),
556 1 : interaction_multipole_triplets[i].end());
557 1 : interaction_multipole_triplets[i].clear();
558 1 : if (i == 0) {
559 0 : interaction_multipole[i] = this->basisvectors.adjoint() *
560 0 : interaction_multipole[i].template selfadjointView<Eigen::Lower>() *
561 : this->basisvectors;
562 : } else {
563 1 : interaction_multipole[i] =
564 1 : this->basisvectors.adjoint() * interaction_multipole[i] * this->basisvectors;
565 1 : interaction_multipole[-i] = std::pow(-1, i) * interaction_multipole[i].adjoint();
566 : }
567 : }
568 : }
569 : }
570 :
571 : ////////////////////////////////////////////////////////////////////
572 : /// Method that allows base class to construct Hamiltonian /////////
573 : ////////////////////////////////////////////////////////////////////
574 :
575 : template <typename Scalar>
576 154 : void SystemOne<Scalar>::addInteraction() {
577 : // Build the total Hamiltonian
578 154 : double tolerance = 1e-24;
579 :
580 154 : if (std::abs(efield_spherical[+0]) > tolerance) {
581 130 : this->hamiltonian -= interaction_efield[+0] * efield_spherical[+0];
582 : }
583 154 : if (std::abs(efield_spherical[-1]) > tolerance) {
584 18 : this->hamiltonian += interaction_efield[+1] * efield_spherical[-1];
585 : }
586 154 : if (std::abs(efield_spherical[+1]) > tolerance) {
587 18 : this->hamiltonian += interaction_efield[-1] * efield_spherical[+1];
588 : }
589 154 : if (std::abs(bfield_spherical[+0]) > tolerance) {
590 29 : this->hamiltonian -= interaction_bfield[+0] * bfield_spherical[+0];
591 : }
592 154 : if (std::abs(bfield_spherical[-1]) > tolerance) {
593 30 : this->hamiltonian += interaction_bfield[+1] * bfield_spherical[-1];
594 : }
595 154 : if (std::abs(bfield_spherical[+1]) > tolerance) {
596 30 : this->hamiltonian += interaction_bfield[-1] * bfield_spherical[+1];
597 : }
598 :
599 154 : if (diamagnetism && std::abs(diamagnetism_terms[{{0, +0}}]) > tolerance) {
600 39 : this->hamiltonian += interaction_diamagnetism[{{0, +0}}] * diamagnetism_terms[{{0, +0}}];
601 : }
602 154 : if (diamagnetism && std::abs(diamagnetism_terms[{{2, +0}}]) > tolerance) {
603 39 : this->hamiltonian -= interaction_diamagnetism[{{2, +0}}] * diamagnetism_terms[{{2, +0}}];
604 : }
605 154 : if (diamagnetism && std::abs(diamagnetism_terms[{{2, +1}}]) > tolerance) {
606 16 : this->hamiltonian +=
607 32 : interaction_diamagnetism[{{2, +1}}] * diamagnetism_terms[{{2, +1}}] * std::sqrt(3);
608 : }
609 154 : if (diamagnetism && std::abs(diamagnetism_terms[{{2, -1}}]) > tolerance) {
610 16 : this->hamiltonian +=
611 32 : interaction_diamagnetism[{{2, -1}}] * diamagnetism_terms[{{2, -1}}] * std::sqrt(3);
612 : }
613 154 : if (diamagnetism && std::abs(diamagnetism_terms[{{2, +2}}]) > tolerance) {
614 30 : this->hamiltonian -=
615 60 : interaction_diamagnetism[{{2, +2}}] * diamagnetism_terms[{{2, +2}}] * std::sqrt(1.5);
616 : }
617 154 : if (diamagnetism && std::abs(diamagnetism_terms[{{2, -2}}]) > tolerance) {
618 30 : this->hamiltonian -=
619 60 : interaction_diamagnetism[{{2, -2}}] * diamagnetism_terms[{{2, -2}}] * std::sqrt(1.5);
620 : }
621 154 : if (charge != 0 && distance != std::numeric_limits<double>::max()) {
622 2 : for (unsigned int order = 1; order <= ordermax; ++order) {
623 1 : double powerlaw = 1. / std::pow(distance, order + 1);
624 1 : this->hamiltonian += interaction_multipole[order] * charge * powerlaw;
625 : }
626 : }
627 154 : }
628 :
629 : ////////////////////////////////////////////////////////////////////
630 : /// Method that allows base class to transform the interaction /////
631 : ////////////////////////////////////////////////////////////////////
632 :
633 : template <typename Scalar>
634 18 : void SystemOne<Scalar>::transformInteraction(const Eigen::SparseMatrix<Scalar> &transformator) {
635 34 : for (auto &entry : interaction_efield) {
636 16 : entry.second = transformator.adjoint() * entry.second * transformator;
637 : }
638 42 : for (auto &entry : interaction_bfield) {
639 24 : entry.second = transformator.adjoint() * entry.second * transformator;
640 : }
641 66 : for (auto &entry : interaction_diamagnetism) {
642 48 : entry.second = transformator.adjoint() * entry.second * transformator; // NOLINT
643 : }
644 18 : for (auto &entry : interaction_multipole) {
645 0 : entry.second = transformator.adjoint() * entry.second * transformator; // NOLINT
646 : }
647 18 : }
648 :
649 : ////////////////////////////////////////////////////////////////////
650 : /// Method that allows base class to delete the interaction ////////
651 : ////////////////////////////////////////////////////////////////////
652 :
653 : template <typename Scalar>
654 5 : void SystemOne<Scalar>::deleteInteraction() {
655 5 : interaction_efield.clear();
656 5 : interaction_bfield.clear();
657 5 : interaction_diamagnetism.clear();
658 5 : interaction_multipole.clear();
659 5 : }
660 :
661 : ////////////////////////////////////////////////////////////////////
662 : /// Methods that allows base class to rotate states ////////////////
663 : ////////////////////////////////////////////////////////////////////
664 :
665 : template <typename Scalar>
666 : Eigen::SparseMatrix<Scalar>
667 3 : SystemOne<Scalar>::rotateStates(const std::vector<size_t> &states_indices, double alpha,
668 : double beta, double gamma) {
669 : // Initialize Wigner D matrix
670 3 : WignerD wigner;
671 :
672 : // Rotate state
673 6 : std::vector<Eigen::Triplet<Scalar>> states_rotated_triplets;
674 3 : states_rotated_triplets.reserve(
675 3 : std::min(static_cast<size_t>(10), this->states.size()) *
676 3 : states_indices.size()); // TODO std::min( 2*jmax+1, states.size() ) * states_indices.size()
677 :
678 3 : size_t current = 0;
679 8 : for (auto const &idx : states_indices) {
680 5 : this->addRotated(this->states[idx].state, current++, states_rotated_triplets, wigner, alpha,
681 : beta, gamma);
682 : }
683 :
684 3 : Eigen::SparseMatrix<Scalar> states_rotated(this->states.size(), states_indices.size());
685 3 : states_rotated.setFromTriplets(states_rotated_triplets.begin(), states_rotated_triplets.end());
686 3 : states_rotated_triplets.clear();
687 :
688 6 : return states_rotated;
689 : }
690 :
691 : template <typename Scalar>
692 9 : Eigen::SparseMatrix<Scalar> SystemOne<Scalar>::buildStaterotator(double alpha, double beta,
693 : double gamma) {
694 : // Initialize Wigner D matrix
695 9 : WignerD wigner;
696 :
697 : // Build rotator
698 18 : std::vector<Eigen::Triplet<Scalar>> rotator_triplets;
699 9 : rotator_triplets.reserve(
700 9 : std::min(static_cast<size_t>(10), this->states.size()) *
701 9 : this->states.size()); // TODO std::min( 2*jmax+1, states.size() ) * states.size()
702 :
703 549 : for (auto const &entry : this->states) {
704 270 : this->addRotated(entry.state, entry.idx, rotator_triplets, wigner, alpha, beta, gamma);
705 : }
706 :
707 9 : Eigen::SparseMatrix<Scalar> rotator(this->states.size(), this->states.size());
708 9 : rotator.setFromTriplets(rotator_triplets.begin(), rotator_triplets.end()); // NOLINT
709 9 : rotator_triplets.clear();
710 :
711 18 : return rotator;
712 : }
713 :
714 : ////////////////////////////////////////////////////////////////////
715 : /// Method that allows base class to combine systems ///////////////
716 : ////////////////////////////////////////////////////////////////////
717 :
718 : template <typename Scalar>
719 5 : void SystemOne<Scalar>::incorporate(SystemBase<Scalar, StateOne> &system) {
720 : // Combine parameters
721 5 : if (species != dynamic_cast<SystemOne<Scalar> &>(system).species) {
722 0 : throw std::runtime_error(
723 : "The value of the variable 'element' must be the same for both systems.");
724 : }
725 5 : if (efield != dynamic_cast<SystemOne<Scalar> &>(system).efield) {
726 0 : throw std::runtime_error("The value of the variable 'distance' must be the same for both "
727 : "systems."); // implies that efield_spherical is the same, too
728 : }
729 5 : if (bfield != dynamic_cast<SystemOne<Scalar> &>(system).bfield) {
730 0 : throw std::runtime_error("The value of the variable 'angle' must be the same for both "
731 : "systems."); // implies that
732 : // bfield_spherical
733 : // is the same,
734 : // too
735 : }
736 5 : if (diamagnetism != dynamic_cast<SystemOne<Scalar> &>(system).diamagnetism) {
737 0 : throw std::runtime_error(
738 : "The value of the variable 'ordermax' must be the same for both systems.");
739 : }
740 :
741 : // Combine symmetries
742 5 : unsigned int num_different_symmetries = 0;
743 5 : if (sym_reflection != dynamic_cast<SystemOne<Scalar> &>(system).sym_reflection) {
744 2 : sym_reflection = NA;
745 2 : ++num_different_symmetries;
746 : }
747 8 : if (!(sym_rotation.size() == dynamic_cast<SystemOne<Scalar> &>(system).sym_rotation.size() &&
748 3 : std::equal(sym_rotation.begin(), sym_rotation.end(),
749 3 : dynamic_cast<SystemOne<Scalar> &>(system).sym_rotation.begin()))) {
750 6 : if (sym_rotation.count(static_cast<float>(ARB)) != 0 ||
751 6 : dynamic_cast<SystemOne<Scalar> &>(system).sym_rotation.count(static_cast<float>(ARB)) !=
752 : 0) {
753 0 : sym_rotation = {static_cast<float>(ARB)};
754 : } else {
755 3 : sym_rotation.insert(dynamic_cast<SystemOne<Scalar> &>(system).sym_rotation.begin(),
756 3 : dynamic_cast<SystemOne<Scalar> &>(system).sym_rotation.end());
757 : }
758 3 : ++num_different_symmetries;
759 : }
760 5 : if (num_different_symmetries > 1) {
761 0 : std::cerr << "Warning: The systems differ in more than one symmetry. For the combined "
762 : "system, the notion of symmetries might be meaningless."
763 0 : << std::endl;
764 : }
765 :
766 : // Clear cached interaction
767 5 : this->deleteInteraction();
768 5 : }
769 :
770 : ////////////////////////////////////////////////////////////////////
771 : /// Utility methods ////////////////////////////////////////////////
772 : ////////////////////////////////////////////////////////////////////
773 :
774 : template <typename Scalar>
775 3930 : void SystemOne<Scalar>::addSymmetrizedBasisvectors(
776 : const StateOne &state, size_t &idx, const double &energy,
777 : std::vector<Eigen::Triplet<Scalar>> &basisvectors_triplets,
778 : std::vector<Eigen::Triplet<Scalar>> &hamiltonian_triplets, parity_t &sym_reflection_local) {
779 : // In case of reflection symmetry, skip half of the basis vectors
780 3930 : if (sym_reflection_local != NA && state.getM() != 0) {
781 56 : if (state.getM() < 0) {
782 28 : return;
783 : }
784 : }
785 :
786 : // Store the energy of the unperturbed one atom state
787 3902 : hamiltonian_triplets.emplace_back(idx, idx, energy);
788 :
789 : // Adapt the normalization if required by symmetries
790 3902 : Scalar value = 1;
791 3902 : if (sym_reflection_local != NA && state.getM() != 0) {
792 28 : value /= std::sqrt(2);
793 : }
794 :
795 : // Add an entry to the current basis vector
796 3902 : this->addBasisvectors(state, idx, value, basisvectors_triplets);
797 :
798 : // Add further entries to the current basis vector if required by symmetries
799 3902 : if (sym_reflection_local != NA && state.getM() != 0) {
800 28 : value *= std::pow(-1, state.getL() + state.getM() - state.getJ()) *
801 28 : utils::imaginary_unit<Scalar>();
802 28 : value *= (sym_reflection_local == EVEN) ? 1 : -1;
803 : // S_y is invariant under reflection through xz-plane
804 : // TODO is the s quantum number of importance here?
805 28 : this->addBasisvectors(state.getReflected(), idx, value, basisvectors_triplets);
806 : }
807 :
808 3902 : ++idx;
809 : }
810 :
811 : template <typename Scalar>
812 3930 : void SystemOne<Scalar>::addBasisvectors(
813 : const StateOne &state, const size_t &idx, const Scalar &value,
814 : std::vector<Eigen::Triplet<Scalar>> &basisvectors_triplets) {
815 3930 : auto state_iter = this->states.template get<1>().find(state);
816 :
817 : size_t row;
818 3930 : if (state_iter != this->states.template get<1>().end()) {
819 0 : row = state_iter->idx;
820 : } else {
821 3930 : row = this->states.size();
822 3930 : this->states.push_back(enumerated_state<StateOne>(row, state));
823 : }
824 :
825 3930 : basisvectors_triplets.emplace_back(row, idx, value);
826 3930 : }
827 :
828 : template <typename Scalar>
829 113 : void SystemOne<Scalar>::changeToSphericalbasis(std::array<double, 3> field,
830 : std::unordered_map<int, double> &field_spherical) {
831 113 : if (field[1] != 0) {
832 0 : throw std::runtime_error(
833 : "For fields with non-zero y-coordinates, a complex data type is needed.");
834 : }
835 113 : field_spherical[1] = -field[0] / std::sqrt(2);
836 113 : field_spherical[-1] = field[0] / std::sqrt(2);
837 113 : field_spherical[0] = field[2];
838 113 : }
839 :
840 : template <typename Scalar>
841 33 : void SystemOne<Scalar>::changeToSphericalbasis(
842 : std::array<double, 3> field, std::unordered_map<int, std::complex<double>> &field_spherical) {
843 33 : field_spherical[1] = std::complex<double>(-field[0] / std::sqrt(2), -field[1] / std::sqrt(2));
844 33 : field_spherical[-1] = std::complex<double>(field[0] / std::sqrt(2), -field[1] / std::sqrt(2));
845 33 : field_spherical[0] = std::complex<double>(field[2], 0);
846 33 : }
847 :
848 : template <typename Scalar>
849 19128 : void SystemOne<Scalar>::addTriplet(std::vector<Eigen::Triplet<Scalar>> &triplets,
850 : const size_t r_idx, const size_t c_idx, const Scalar val) {
851 19128 : triplets.emplace_back(r_idx, c_idx, val);
852 19128 : }
853 :
854 : template <typename Scalar>
855 0 : void SystemOne<Scalar>::rotateVector(std::array<double, 3> &field, std::array<double, 3> &to_z_axis,
856 : std::array<double, 3> &to_y_axis) {
857 0 : auto field_mapped = Eigen::Map<Eigen::Matrix<double, 3, 1>>(&field[0]);
858 :
859 0 : if (field_mapped.norm() != 0) {
860 0 : Eigen::Matrix<double, 3, 3> rotator = this->buildRotator(to_z_axis, to_y_axis);
861 0 : field_mapped = rotator.transpose() * field_mapped;
862 : }
863 0 : }
864 :
865 : template <typename Scalar>
866 3 : void SystemOne<Scalar>::rotateVector(std::array<double, 3> &field, double alpha, double beta,
867 : double gamma) {
868 3 : auto field_mapped = Eigen::Map<Eigen::Matrix<double, 3, 1>>(&field[0]);
869 :
870 3 : if (field_mapped.norm() != 0) {
871 3 : Eigen::Matrix<double, 3, 3> rotator = this->buildRotator(alpha, beta, gamma);
872 3 : field_mapped = rotator.transpose() * field_mapped;
873 : }
874 3 : }
875 :
876 : template <typename Scalar>
877 15 : bool SystemOne<Scalar>::isRefelectionAndRotationCompatible() {
878 15 : if (sym_rotation.count(static_cast<float>(ARB)) != 0 || sym_reflection == NA) {
879 15 : return true;
880 : }
881 :
882 0 : for (const auto &s : sym_rotation) {
883 0 : if (sym_rotation.count(-s) == 0) {
884 0 : return false;
885 : }
886 : }
887 :
888 0 : return true;
889 : }
890 :
891 : template class SystemOne<std::complex<double>>;
892 : template class SystemOne<double>;
|
