Line data Source code
1 : /*
2 : * Copyright (c) 2016 Sebastian Weber, Henri Menke, Johannes Block. 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 "SystemTwo.hpp"
21 : #include "Constants.hpp"
22 : #include "GreenTensor.hpp"
23 : #include "MatrixElementCache.hpp"
24 :
25 : #include <algorithm>
26 : #include <cmath>
27 : #include <limits>
28 : #include <numeric>
29 : #include <string>
30 : #include <unordered_set>
31 : #include <vector>
32 :
33 : template <typename Scalar>
34 0 : SystemTwo<Scalar>::SystemTwo()
35 0 : : minimal_le_roy_radius(std::numeric_limits<double>::max()),
36 0 : distance(std::numeric_limits<double>::max()), distance_x(0), distance_y(0),
37 0 : distance_z(std::numeric_limits<double>::max()), GTbool(false),
38 0 : surface_distance(std::numeric_limits<double>::max()), ordermax(3), sym_permutation(NA),
39 0 : sym_inversion(NA), sym_reflection(NA), sym_rotation({ARB}) {}
40 :
41 : template <typename Scalar>
42 62 : SystemTwo<Scalar>::SystemTwo(const SystemOne<Scalar> &b1, const SystemOne<Scalar> &b2,
43 : MatrixElementCache &cache)
44 248 : : SystemBase<Scalar, StateTwo>(cache), species({{b1.getSpecies(), b2.getSpecies()}}),
45 62 : system1(b1), system2(b2), minimal_le_roy_radius(std::numeric_limits<double>::max()),
46 62 : distance(std::numeric_limits<double>::max()), distance_x(0), distance_y(0),
47 62 : distance_z(std::numeric_limits<double>::max()), GTbool(false),
48 62 : surface_distance(std::numeric_limits<double>::max()), ordermax(3), sym_permutation(NA),
49 124 : sym_inversion(NA), sym_reflection(NA), sym_rotation({ARB}) {}
50 :
51 : template <typename Scalar>
52 0 : SystemTwo<Scalar>::SystemTwo(const SystemOne<Scalar> &b1, const SystemOne<Scalar> &b2,
53 : MatrixElementCache &cache, bool memory_saving)
54 : : SystemBase<Scalar, StateTwo>(cache, memory_saving),
55 0 : species({{b1.getSpecies(), b2.getSpecies()}}), system1(b1), system2(b2),
56 0 : minimal_le_roy_radius(std::numeric_limits<double>::max()),
57 0 : distance(std::numeric_limits<double>::max()), distance_x(0), distance_y(0),
58 0 : distance_z(std::numeric_limits<double>::max()), GTbool(false),
59 0 : surface_distance(std::numeric_limits<double>::max()), ordermax(3), sym_permutation(NA),
60 0 : sym_inversion(NA), sym_reflection(NA), sym_rotation({ARB}) {}
61 :
62 : template <typename Scalar>
63 : std::vector<StateOne>
64 95 : SystemTwo<Scalar>::getStatesFirst() { // TODO @hmenke typemap for "state_set<StateOne>"
65 95 : this->buildBasis();
66 190 : std::unordered_set<StateOne> states_one;
67 14631 : for (const auto &entry : this->states) {
68 14536 : states_one.insert(StateOne(entry.state.getFirstState()));
69 : }
70 190 : return std::vector<StateOne>(states_one.begin(), states_one.end());
71 : }
72 :
73 : template <typename Scalar>
74 : std::vector<StateOne>
75 95 : SystemTwo<Scalar>::getStatesSecond() { // TODO @hmenke typemap for "state_set<StateOne>"
76 95 : this->buildBasis();
77 190 : std::unordered_set<StateOne> states_one;
78 14631 : for (const auto &entry : this->states) {
79 14536 : states_one.insert(StateOne(entry.state.getSecondState()));
80 : }
81 190 : return std::vector<StateOne>(states_one.begin(), states_one.end());
82 : }
83 :
84 : template <typename Scalar>
85 0 : const std::array<std::string, 2> &SystemTwo<Scalar>::getSpecies() {
86 0 : return species;
87 : }
88 :
89 : template <typename Scalar>
90 46 : void SystemTwo<Scalar>::enableGreenTensor(bool GTboolean) {
91 46 : this->onParameterChange();
92 46 : GTbool = GTboolean;
93 46 : if (!GTbool && surface_distance != std::numeric_limits<double>::max()) {
94 0 : throw std::runtime_error("If there is interaction with a surface, the Green tensor "
95 : "approach must not be disabled.");
96 : }
97 46 : }
98 :
99 : template <typename Scalar>
100 37 : void SystemTwo<Scalar>::setSurfaceDistance(double d) {
101 37 : this->onParameterChange();
102 37 : surface_distance = d;
103 37 : if (surface_distance != std::numeric_limits<double>::max()) {
104 37 : this->enableGreenTensor(true);
105 : }
106 37 : }
107 :
108 : template <typename Scalar>
109 149 : void SystemTwo<Scalar>::setDistance(double d) {
110 149 : this->onParameterChange();
111 149 : distance_x = distance_x / distance * d;
112 149 : distance_y = distance_y / distance * d;
113 149 : distance_z = distance_z / distance * d;
114 149 : distance = d;
115 149 : }
116 :
117 : template <typename Scalar>
118 0 : void SystemTwo<Scalar>::setDistanceVector(std::array<double, 3> d) {
119 0 : this->onParameterChange();
120 0 : distance_x = d[0];
121 0 : distance_y = d[1];
122 0 : distance_z = d[2];
123 0 : distance = std::sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
124 0 : if (distance_y != 0) {
125 0 : this->enableGreenTensor(true);
126 : }
127 0 : }
128 :
129 : template <typename Scalar>
130 15 : void SystemTwo<Scalar>::setAngle(double a) {
131 15 : this->onParameterChange();
132 15 : distance_x = distance * std::sin(a);
133 15 : distance_y = 0;
134 15 : distance_z = distance * std::cos(a);
135 15 : }
136 :
137 : template <typename Scalar>
138 28 : void SystemTwo<Scalar>::setOrder(double o) {
139 28 : this->onParameterChange();
140 28 : ordermax = o;
141 28 : }
142 :
143 : template <typename Scalar>
144 23 : void SystemTwo<Scalar>::setConservedParityUnderPermutation(parity_t parity) {
145 23 : this->onSymmetryChange();
146 23 : sym_permutation = parity;
147 23 : }
148 :
149 : template <typename Scalar>
150 19 : void SystemTwo<Scalar>::setConservedParityUnderInversion(parity_t parity) {
151 19 : this->onSymmetryChange();
152 19 : sym_inversion = parity;
153 19 : }
154 :
155 : template <typename Scalar>
156 10 : void SystemTwo<Scalar>::setConservedParityUnderReflection(parity_t parity) {
157 10 : this->onSymmetryChange();
158 10 : sym_reflection = parity;
159 10 : if (!this->isRefelectionAndRotationCompatible()) {
160 0 : throw std::runtime_error("The conserved parity under reflection is not compatible to the "
161 : "previously specified conserved momenta.");
162 : }
163 : // if (sym_reflection != NA) std::cerr << "Warning: The one-atom states must already be
164 : // reflection symmetric in order to build reflection symmetric two-atom states." << std::endl;
165 :
166 : // TODO make it work with one-atom states that are not pre-symmetrized
167 10 : }
168 :
169 : template <typename Scalar>
170 16 : void SystemTwo<Scalar>::setConservedMomentaUnderRotation(const std::set<int> &momenta) {
171 16 : if (momenta.count(ARB) != 0 && momenta.size() > 1) {
172 0 : throw std::runtime_error(
173 : "If ARB (=arbitrary momentum) is specified, momenta must not be passed explicitely.");
174 : }
175 16 : this->onSymmetryChange();
176 16 : sym_rotation = momenta;
177 16 : if (!this->isRefelectionAndRotationCompatible()) {
178 0 : throw std::runtime_error("The conserved momenta are not compatible to the previously "
179 : "specified conserved parity under reflection.");
180 : }
181 16 : }
182 :
183 : template <typename Scalar>
184 2 : void SystemTwo<Scalar>::setOneAtomBasisvectors(const std::vector<std::array<size_t, 2>> &indices) {
185 : // Check that all pairs of indices are unique
186 4 : std::vector<std::array<size_t, 2>> tmp(indices);
187 2 : std::sort(tmp.begin(), tmp.end());
188 2 : if (auto const it = std::adjacent_find(tmp.begin(), tmp.end()); it != tmp.end()) {
189 2 : std::string const which =
190 2 : "[" + std::to_string((*it)[0]) + "," + std::to_string((*it)[1]) + "]";
191 1 : throw std::runtime_error("The pairs of indices are not unique: " + which);
192 : }
193 :
194 1 : one_atom_basisvectors_indices = indices;
195 1 : }
196 :
197 : ////////////////////////////////////////////////////////////////////
198 : /// Method that allows base class to initialize Basis //////////////
199 : ////////////////////////////////////////////////////////////////////
200 :
201 : template <typename Scalar>
202 75 : void SystemTwo<Scalar>::initializeBasis() {
203 : ////////////////////////////////////////////////////////////////////
204 : /// Restrict one atom states to the allowed quantum numbers ////////
205 : ////////////////////////////////////////////////////////////////////
206 :
207 75 : system1.diagonalize(); // it is important to call this method here!
208 75 : system1.restrictN(this->range_n);
209 75 : system1.restrictL(this->range_l);
210 75 : system1.restrictJ(this->range_j);
211 75 : system1.restrictM(this->range_m);
212 75 : system2.diagonalize(); // it is important to call this method here!
213 75 : system2.restrictN(this->range_n);
214 75 : system2.restrictL(this->range_l);
215 75 : system2.restrictJ(this->range_j);
216 75 : system2.restrictM(this->range_m);
217 :
218 : ////////////////////////////////////////////////////////////////////
219 : /// Check whether the single atom states fit to the symmetries /////
220 : ////////////////////////////////////////////////////////////////////
221 :
222 75 : if (sym_permutation != NA) {
223 : // TODO check system1 == system2
224 : }
225 :
226 : // TODO consider further symmetries and check whether they are applicable
227 :
228 : ////////////////////////////////////////////////////////////////////
229 : /// Check which basis vectors contain artificial states ////////////
230 : ////////////////////////////////////////////////////////////////////
231 :
232 : // TODO check whether system1 == system2
233 :
234 150 : std::vector<bool> artificial1(system1.getNumBasisvectors(), false);
235 2641 : for (size_t col = 0; col < system1.getNumBasisvectors(); ++col) {
236 15158 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(system1.getBasisvectors(),
237 : col);
238 12592 : triple; ++triple) {
239 12592 : if (system1.getStatesMultiIndex()[triple.row()].state.isArtificial()) {
240 0 : artificial1[triple.col()] = true;
241 : }
242 : }
243 : }
244 :
245 150 : std::vector<bool> artificial2(system2.getNumBasisvectors(), false);
246 2641 : for (size_t col = 0; col < system2.getNumBasisvectors(); ++col) {
247 15158 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(system2.getBasisvectors(),
248 : col);
249 12592 : triple; ++triple) {
250 12592 : if (system2.getStatesMultiIndex()[triple.row()].state.isArtificial()) {
251 0 : artificial2[triple.col()] = true;
252 : }
253 : }
254 : }
255 :
256 : ////////////////////////////////////////////////////////////////////
257 : /// Build two atom states //////////////////////////////////////////
258 : ////////////////////////////////////////////////////////////////////
259 :
260 75 : auto &cache = *this->m_cache;
261 :
262 : /// Combine one atom states ////////////////////////////////////////
263 :
264 150 : std::vector<Eigen::Triplet<Scalar>> hamiltonian_triplets;
265 75 : hamiltonian_triplets.reserve(system1.getNumBasisvectors() * system2.getNumBasisvectors());
266 150 : this->states.reserve(system1.getNumStates() * system2.getNumStates() +
267 75 : this->states_to_add.size());
268 150 : std::vector<Eigen::Triplet<Scalar>> basisvectors_triplets; // TODO reserve
269 225 : std::vector<double> sqnorm_list(
270 75 : system1.getNumStates() * system2.getNumStates() + this->states_to_add.size(), 0);
271 :
272 75 : int M = 0;
273 75 : int parityL = 0;
274 75 : int parityJ = 0;
275 75 : int parityM = 0;
276 :
277 75 : size_t col_new = 0;
278 :
279 : // Generate vector containing pairs of indices of one atom basisvectors
280 150 : std::vector<std::array<size_t, 2>> pairs;
281 75 : if (!one_atom_basisvectors_indices.empty()) {
282 1 : pairs = one_atom_basisvectors_indices;
283 : } else {
284 74 : pairs.reserve(system1.getNumBasisvectors() * system2.getNumBasisvectors());
285 2626 : for (size_t col_1 = 0; col_1 < system1.getNumBasisvectors(); ++col_1) {
286 286712 : for (size_t col_2 = 0; col_2 < system2.getNumBasisvectors(); ++col_2) {
287 284160 : pairs.push_back({col_1, col_2});
288 : }
289 : }
290 : }
291 :
292 : // Loop over the pairs
293 284239 : for (const auto &[col_1, col_2] : pairs) {
294 :
295 : // In case of artificial states, some symmetries won't work
296 284164 : auto sym_inversion_local = sym_inversion;
297 284164 : auto sym_reflection_local = sym_reflection;
298 284164 : auto sym_rotation_local = sym_rotation;
299 284164 : if (artificial1[col_1] || artificial2[col_2]) {
300 0 : if (sym_inversion_local != NA || sym_reflection_local != NA ||
301 0 : sym_rotation_local.count(ARB) == 0) {
302 : std::cerr
303 0 : << "WARNING: Only permutation symmetry can be applied to artificial states."
304 0 : << std::endl;
305 : }
306 0 : sym_inversion_local = NA;
307 0 : sym_reflection_local = NA;
308 0 : sym_rotation_local = std::set<int>({ARB});
309 : }
310 :
311 : // In case of inversion or permutation symmetry: skip half of the basis vector pairs
312 284164 : if ((sym_inversion_local == EVEN && col_1 <= col_2) || // gerade
313 134840 : (sym_inversion_local == ODD && col_1 < col_2)) { // ungerade
314 123431 : continue;
315 : }
316 160733 : if ((sym_permutation == EVEN && col_1 <= col_2) || // sym
317 160128 : (sym_permutation == ODD && col_1 < col_2)) { // asym
318 7614 : continue;
319 : }
320 :
321 : // Continue if the pair statet energy is not valid
322 153119 : double energy = std::real(system1.getHamiltonian().coeff(col_1, col_1) +
323 153119 : system2.getHamiltonian().coeff(col_2, col_2));
324 153119 : if (!this->checkIsEnergyValid(energy)) {
325 141699 : continue;
326 : }
327 :
328 : // Store the pair state energy
329 11420 : hamiltonian_triplets.emplace_back(col_new, col_new, energy);
330 :
331 : // Build the basis vector that corresponds to the stored pair state energy
332 69387 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple_1(system1.getBasisvectors(),
333 : col_1);
334 57967 : triple_1; ++triple_1) {
335 57967 : size_t row_1 = triple_1.row();
336 115934 : StateOne state_1 = system1.getStatesMultiIndex()[row_1].state;
337 :
338 1282528 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple_2(
339 57967 : system2.getBasisvectors(), col_2);
340 1224561 : triple_2; ++triple_2) {
341 1224561 : size_t row_2 = triple_2.row();
342 1224561 : StateOne state_2 = system2.getStatesMultiIndex()[row_2].state;
343 :
344 1224561 : Scalar value_new = triple_1.value() * triple_2.value();
345 :
346 1224561 : if (!artificial1[col_1] && !artificial2[col_2]) {
347 1224561 : M = state_1.getM() + state_2.getM();
348 1224561 : parityL = std::pow(-1, state_1.getL() + state_2.getL());
349 1224561 : parityJ = std::pow(-1, state_1.getJ() + state_2.getJ());
350 1224561 : parityM = std::pow(-1, M);
351 : }
352 :
353 1224561 : bool different = col_1 != col_2;
354 :
355 : // Consider rotation symmetry
356 1224561 : if (sym_rotation_local.count(ARB) == 0 && sym_rotation_local.count(M) == 0) {
357 31378 : continue;
358 : }
359 :
360 : // Combine symmetries
361 1193183 : bool skip_reflection = false;
362 1193183 : if (different) {
363 : // In case of inversion and permutation symmetry: the inversion symmetric
364 : // state is already permutation symmetric
365 1107617 : if (sym_inversion_local != NA && sym_permutation != NA) {
366 24818 : if (((sym_inversion_local == EVEN) ? -parityL : parityL) !=
367 12409 : ((sym_permutation == EVEN) ? -1 : 1)) {
368 6162 : continue; // parity under inversion and permutation is different
369 : }
370 : }
371 :
372 : // In case of inversion or permutation and reflection symmetry: the
373 : // inversion or permutation symmetric state is already reflection symmetric
374 1107215 : if (sym_inversion_local != NA && sym_reflection_local != NA &&
375 1107215 : StateTwo(state_1.getReflected(), state_2.getReflected()) ==
376 1107215 : StateTwo(state_2, state_1)) {
377 800 : if (((sym_inversion_local == EVEN) ? -parityL : parityL) !=
378 400 : ((sym_reflection_local == EVEN) ? parityL * parityJ * parityM
379 200 : : -parityL * parityJ * parityM)) {
380 200 : continue; // parity under inversion and reflection is different
381 : }
382 200 : skip_reflection = true; // parity under inversion and reflection is the same
383 :
384 1106415 : } else if (sym_permutation != NA && sym_reflection_local != NA &&
385 1106415 : StateTwo(state_1.getReflected(), state_2.getReflected()) ==
386 1106415 : StateTwo(state_2, state_1)) {
387 0 : if (((sym_permutation == EVEN) ? -1 : 1) !=
388 0 : ((sym_reflection_local == EVEN) ? parityL * parityJ * parityM
389 0 : : -parityL * parityJ * parityM)) {
390 0 : continue; // parity under permutation and reflection is different
391 : }
392 0 : skip_reflection =
393 : true; // parity under permutation and reflection is the same
394 : }
395 : }
396 :
397 : // Adapt the normalization if required by symmetries
398 1186821 : if ((sym_inversion_local != NA || sym_permutation != NA) && different) {
399 873521 : value_new /= std::sqrt(2);
400 : }
401 1186821 : if (sym_reflection_local != NA && !skip_reflection) {
402 82704 : value_new /= std::sqrt(2) * std::sqrt(2);
403 : // the second factor std::sqrt(2) is because of double counting
404 : }
405 :
406 : // Add an entry to the current basis vector
407 1186821 : this->addBasisvectors(StateTwo(state_1, state_2), col_new, value_new,
408 : basisvectors_triplets, sqnorm_list);
409 :
410 : // Add further entries to the current basis vector if required by symmetries
411 1186821 : if (different) {
412 1101255 : if (sym_inversion_local != NA) {
413 41751 : Scalar v = value_new;
414 41751 : v *= (sym_inversion_local == EVEN) ? -parityL : parityL;
415 41751 : this->addBasisvectors(StateTwo(state_2, state_1), col_new, v,
416 : basisvectors_triplets, sqnorm_list);
417 1059504 : } else if (sym_permutation != NA) {
418 831770 : Scalar v = value_new;
419 831770 : v *= (sym_permutation == EVEN) ? -1 : 1;
420 831770 : this->addBasisvectors(StateTwo(state_2, state_1), col_new, v,
421 : basisvectors_triplets, sqnorm_list);
422 : }
423 : }
424 1186821 : if (sym_reflection_local != NA && !skip_reflection) {
425 82704 : Scalar v = value_new;
426 124056 : v *= (sym_reflection_local == EVEN) ? parityL * parityJ * parityM
427 41352 : : -parityL * parityJ * parityM;
428 82704 : this->addBasisvectors(StateTwo(state_1.getReflected(), state_2.getReflected()),
429 : col_new, v, basisvectors_triplets, sqnorm_list);
430 :
431 82704 : if (different) {
432 76704 : if (sym_inversion_local != NA) {
433 5360 : Scalar v = value_new;
434 8040 : v *= (sym_reflection_local == EVEN) ? parityL * parityJ * parityM
435 2680 : : -parityL * parityJ * parityM;
436 5360 : v *= (sym_inversion_local == EVEN) ? -parityL : parityL;
437 5360 : this->addBasisvectors(
438 10720 : StateTwo(state_2.getReflected(), state_1.getReflected()), col_new,
439 : v, basisvectors_triplets, sqnorm_list);
440 71344 : } else if (sym_permutation != NA) {
441 0 : Scalar v = value_new;
442 0 : v *= (sym_reflection_local == EVEN) ? parityL * parityJ * parityM
443 0 : : -parityL * parityJ * parityM;
444 0 : v *= (sym_permutation == EVEN) ? -1 : 1;
445 0 : this->addBasisvectors(
446 0 : StateTwo(state_2.getReflected(), state_1.getReflected()), col_new,
447 : v, basisvectors_triplets, sqnorm_list);
448 : }
449 : }
450 : }
451 : }
452 : }
453 :
454 11420 : ++col_new;
455 : }
456 :
457 : // Delete unecessary storage
458 : // system1 = SystemOne(species[0], cache); // TODO
459 : // system2 = SystemOne(species[1], cache); // TODO
460 :
461 : /// Loop over user-defined states //////////////////////////////////
462 :
463 : // Check that the user-defined states are not already contained in the list of states
464 75 : for (const auto &state : this->states_to_add) {
465 0 : if (this->states.template get<1>().find(state) != this->states.template get<1>().end()) {
466 0 : throw std::runtime_error("The state " + state.str() +
467 : " is already contained in the list of states.");
468 : }
469 0 : for (int idx = 0; idx < 2; ++idx) {
470 0 : if (!state.isArtificial(idx) && state.getSpecies(idx) != species[idx]) {
471 0 : throw std::runtime_error("The state " + state.str() + " is of the wrong species.");
472 : }
473 : }
474 : }
475 :
476 : // Warn if reflection symmetry is selected
477 75 : if (!this->states_to_add.empty() && sym_reflection != NA) {
478 0 : std::cerr << "WARNING: Reflection symmetry cannot be handled for user-defined states."
479 0 : << std::endl;
480 : }
481 :
482 : // Add user-defined states
483 75 : for (const auto &state : this->states_to_add) {
484 0 : bool different = state.getFirstState() != state.getSecondState();
485 :
486 : // Get energy of the state
487 0 : double energy = 0;
488 0 : for (int idx = 0; idx < 2; ++idx) {
489 0 : if (!state.isArtificial(idx)) {
490 0 : energy += state.getEnergy(idx, cache);
491 : }
492 : }
493 :
494 : // In case of artificial states, some symmetries won't work
495 0 : auto sym_inversion_local = sym_inversion;
496 0 : auto sym_rotation_local = sym_rotation;
497 0 : if (state.isArtificial(0) || state.isArtificial(1)) {
498 0 : if (sym_inversion_local != NA || sym_rotation_local.count(ARB) == 0) {
499 : std::cerr
500 0 : << "WARNING: Only permutation symmetry can be applied to artificial states."
501 0 : << std::endl;
502 : }
503 0 : sym_inversion_local = NA;
504 0 : sym_rotation_local = std::set<int>({ARB});
505 : } else {
506 0 : M = state.getM(0) + state.getM(1);
507 0 : parityL = std::pow(-1, state.getL(0) + state.getL(1));
508 : }
509 :
510 : // Consider rotation symmetry
511 0 : if (sym_rotation_local.count(ARB) == 0 && sym_rotation_local.count(M) == 0) {
512 0 : continue;
513 : }
514 :
515 : // Combine symmetries (in case of inversion and permutation symmetry: the inversion
516 : // symmetric state is already permutation symmetric)
517 0 : if (sym_inversion_local != NA && sym_permutation != NA && different) {
518 0 : if (((sym_inversion_local == EVEN) ? -parityL : parityL) !=
519 0 : ((sym_permutation == EVEN) ? -1 : 1)) {
520 0 : continue; // parity under inversion and permutation is different
521 : }
522 : }
523 :
524 : // Check whether the symmetries can be realized with the states available
525 0 : if ((sym_inversion_local != NA || sym_permutation != NA) && different) {
526 0 : auto state_changed = StateTwo(state.getSecondState(), state.getFirstState());
527 0 : if (this->states_to_add.find(state_changed) == this->states_to_add.end()) {
528 0 : throw std::runtime_error("The state " + state_changed.str() +
529 : " required by symmetries cannot be found.");
530 : }
531 : }
532 :
533 : // In case of inversion or permutation symmetry: skip half of the states
534 0 : if ((sym_inversion_local == EVEN &&
535 0 : state.getFirstState() <= state.getSecondState()) || // gerade
536 0 : (sym_inversion_local == ODD &&
537 0 : state.getFirstState() < state.getSecondState())) { // ungerade
538 0 : continue;
539 : }
540 0 : if ((sym_permutation == EVEN && state.getFirstState() <= state.getSecondState()) || // sym
541 0 : (sym_permutation == ODD && state.getFirstState() < state.getSecondState())) { // asym
542 0 : continue;
543 : }
544 :
545 : // Store the energy of the two atom state
546 0 : hamiltonian_triplets.emplace_back(col_new, col_new, energy);
547 :
548 : // Adapt the normalization if required by symmetries
549 0 : Scalar value_new = 1;
550 0 : if ((sym_inversion_local != NA || sym_permutation != NA) && different) {
551 0 : value_new /= std::sqrt(2);
552 : }
553 :
554 : // Add an entry to the current basis vector
555 0 : this->addBasisvectors(state, col_new, value_new, basisvectors_triplets, sqnorm_list);
556 :
557 : // Add further entries to the current basis vector if required by symmetries
558 0 : if (different) {
559 0 : if (sym_inversion_local != NA) {
560 0 : Scalar v = value_new;
561 0 : v *= (sym_inversion_local == EVEN) ? -parityL : parityL;
562 0 : this->addBasisvectors(StateTwo(state.getSecondState(), state.getFirstState()),
563 : col_new, v, basisvectors_triplets, sqnorm_list);
564 0 : } else if (sym_permutation != NA) {
565 0 : Scalar v = value_new;
566 0 : v *= (sym_permutation == EVEN) ? -1 : 1;
567 0 : this->addBasisvectors(StateTwo(state.getSecondState(), state.getFirstState()),
568 : col_new, v, basisvectors_triplets, sqnorm_list);
569 : }
570 : }
571 :
572 0 : ++col_new;
573 : }
574 :
575 : /// Build data /////////////////////////////////////////////////////
576 :
577 75 : this->states.shrink_to_fit();
578 :
579 75 : this->basisvectors.resize(this->states.size(), col_new);
580 75 : this->basisvectors.setFromTriplets(basisvectors_triplets.begin(), basisvectors_triplets.end());
581 75 : basisvectors_triplets.clear();
582 :
583 75 : this->hamiltonian.resize(col_new, col_new);
584 75 : this->hamiltonian.setFromTriplets(hamiltonian_triplets.begin(), hamiltonian_triplets.end());
585 75 : hamiltonian_triplets.clear();
586 :
587 : ////////////////////////////////////////////////////////////////////
588 : /// Remove vectors with too small norm /////////////////////////////
589 : ////////////////////////////////////////////////////////////////////
590 :
591 : // Build transformator and remove vectors (if the squared norm is too small)
592 75 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
593 75 : triplets_transformator.reserve(this->basisvectors.cols());
594 :
595 75 : size_t idx_new = 0;
596 11495 : for (int idx = 0; idx < this->basisvectors.cols(); ++idx) { // idx = col = num basis vector
597 11420 : double_t sqnorm = 0;
598 :
599 : // Calculate the square norm of the columns of the basisvector matrix
600 1223357 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(this->basisvectors, idx);
601 1211937 : triple; ++triple) {
602 1211937 : sqnorm += std::pow(std::abs(triple.value()), 2);
603 : }
604 :
605 11420 : if (sqnorm > this->threshold_for_sqnorm) {
606 4618 : triplets_transformator.emplace_back(idx, idx_new++, 1);
607 : }
608 : }
609 :
610 75 : this->applyRightsideTransformator(triplets_transformator);
611 :
612 : ////////////////////////////////////////////////////////////////////
613 : /// Remove states that barely occur ////////////////////////////////
614 : ////////////////////////////////////////////////////////////////////
615 :
616 : // Build transformator and remove states if the squared norm is to small
617 29999 : this->removeRestrictedStates([=](const enumerated_state<StateTwo> &entry) -> bool {
618 14962 : return sqnorm_list[entry.idx] > this->threshold_for_sqnorm;
619 : });
620 :
621 : /*////////////////////////////////////////////////////////////////////
622 : /// Sort states ////////////////////////////////////////////////////
623 : ////////////////////////////////////////////////////////////////////
624 :
625 : // TODO put this into an extra method
626 :
627 : // Initialize original index locations
628 : std::vector<size_t> idx_sorted(states.size());
629 : iota(idx_sorted.begin(), idx_sorted.end(), 0);
630 :
631 : // Sort indexes based on comparing values in v
632 : std::sort(idx_sorted.begin(), idx_sorted.end(), [=](size_t i1, size_t i2) {return
633 : this->states[i1].state < this->states[i2].state;});
634 :
635 : // Make use of the sorted indexes in order to sort the states and transform the basisvectors
636 : accordingly states_set<StateTwo> states_new; states_new.reserve(states.size());
637 : triplets_transformator.clear();
638 : triplets_transformator.reserve(states.size());
639 :
640 : idx_new = 0;
641 : for (size_t idx : idx_sorted) {
642 : states_new.push_back(enumerated_state<StateTwo>(idx_new, states[idx].state));
643 : triplets_transformator.push_back(Eigen::Triplet<Scalar>(idx_new++,idx,1));
644 : }
645 :
646 : states_new.shrink_to_fit();
647 : states = states_new;
648 : this->applyLeftsideTransformator(triplets_transformator);*/
649 75 : }
650 :
651 : ////////////////////////////////////////////////////////////////////
652 : /// Method that allows base class to calculate the interaction /////
653 : ////////////////////////////////////////////////////////////////////
654 :
655 : template <typename Scalar>
656 198 : void SystemTwo<Scalar>::initializeInteraction() {
657 198 : if (distance == std::numeric_limits<double>::max()) {
658 103 : return;
659 : }
660 :
661 : // Check whether distance is larger than the minimal Le Roy radius
662 198 : this->checkDistance(distance);
663 :
664 : ////////////////////////////////////////////////////////////////////
665 : /// Prepare the calculation of the interaction /////////////////////
666 : ////////////////////////////////////////////////////////////////////
667 :
668 198 : std::unordered_set<int> interaction_angulardipole_keys;
669 198 : std::vector<int> interaction_greentensor_dd_keys;
670 198 : std::vector<int> interaction_greentensor_dq_keys;
671 198 : std::vector<int> interaction_greentensor_qd_keys;
672 198 : std::vector<int> interaction_multipole_keys;
673 :
674 198 : greentensor_terms_dd.clear();
675 198 : greentensor_terms_dq.clear();
676 198 : greentensor_terms_qd.clear();
677 198 : angle_terms.clear();
678 :
679 198 : auto &cache = *this->m_cache;
680 :
681 198 : double tolerance = 1e-16;
682 :
683 198 : if (distance_y != 0 || surface_distance != std::numeric_limits<double>::max() ||
684 396 : (distance_x != 0 && ordermax > 3) || GTbool) {
685 :
686 : // Check compatibility of settings
687 42 : if (!GTbool) {
688 0 : throw std::runtime_error(
689 : "The Green tensor approach must be used if a) the distance vector has a non-zero "
690 : "y-coordinate, b) a surface is present, or c) the interaction angle is non-zero "
691 : "and interaction of higher order than dipole-dipole is considered.");
692 : }
693 42 : if (surface_distance != std::numeric_limits<double>::max() && distance_y != 0) {
694 0 : throw std::runtime_error("The atoms must be in the xz-plane if a surface is present");
695 : }
696 :
697 : // Calculate Green tensor
698 42 : GreenTensor GT(distance_x, distance_y, distance_z);
699 42 : if (surface_distance != std::numeric_limits<double>::max()) {
700 37 : GT.addSurface(surface_distance);
701 : }
702 :
703 : // Determine which interaction matrices have to be calculated
704 42 : if (ordermax >= 3) {
705 42 : const auto &dd_green_tensor = GT.getDDTensor();
706 168 : for (int i1 = 0; i1 < 3; ++i1) {
707 504 : for (int i2 = 0; i2 < 3; ++i2) {
708 378 : if (std::abs(dd_green_tensor(i1, i2)) > tolerance) {
709 196 : greentensor_terms_dd[3 * i2 + i1] = dd_green_tensor(i1, i2);
710 196 : interaction_greentensor_dd_keys.push_back(3 * i2 + i1);
711 : }
712 : }
713 : }
714 : }
715 42 : if (ordermax >= 4) {
716 1 : const auto &dq_green_tensor = GT.getDQTensor();
717 1 : const auto &qd_green_tensor = GT.getQDTensor();
718 4 : for (int i1 = 0; i1 < 3; ++i1) {
719 12 : for (int i2 = 0; i2 < 3; ++i2) {
720 36 : for (int i3 = 0; i3 < 3; i3++) {
721 27 : if (std::abs(dq_green_tensor(i1, i2, i3)) > tolerance) {
722 0 : greentensor_terms_dq[9 * i3 + 3 * i2 + i1] =
723 7 : dq_green_tensor(i1, i2, i3);
724 7 : interaction_greentensor_dq_keys.push_back(9 * i3 + 3 * i2 + i1);
725 : }
726 27 : if (std::abs(qd_green_tensor(i1, i2, i3)) > tolerance) {
727 0 : greentensor_terms_qd[9 * i3 + 3 * i2 + i1] =
728 7 : qd_green_tensor(i1, i2, i3);
729 7 : interaction_greentensor_qd_keys.push_back(9 * i3 + 3 * i2 + i1);
730 : }
731 : }
732 : }
733 : }
734 : }
735 42 : if (ordermax >= 5) {
736 0 : throw std::runtime_error(
737 : "If the Green tensor approach is used, multipole interaction of "
738 : "higher order than dipole-quadrupole cannot be considered.");
739 : }
740 :
741 156 : } else if (distance_x != 0) {
742 :
743 : // We know from before that distance_y is zero so that we can easily calculate the
744 : // interaction angle
745 7 : double angle = std::atan2(distance_x, distance_z);
746 :
747 : // Calculate angle terms
748 7 : double val = 1. - 3. * std::pow(std::cos(angle), 2); // 0,0
749 7 : if (std::abs(val) > tolerance) {
750 7 : angle_terms[3 * (0 + 1) + (0 + 1)] = val;
751 : }
752 7 : val = -1. + 1.5 * std::pow(std::sin(angle), 2); // 1,-1; -1,1
753 7 : if (std::abs(val) > tolerance) {
754 7 : angle_terms[3 * (-1 + 1) + (1 + 1)] = val;
755 : }
756 7 : val = 3. / std::sqrt(2) * std::sin(angle) * std::cos(angle); // 1,0; 0,1
757 7 : if (std::abs(val) > tolerance) {
758 7 : angle_terms[3 * (0 + 1) + (1 + 1)] = val;
759 : }
760 7 : val = -3. / std::sqrt(2) * std::sin(angle) * std::cos(angle); // -1,0; 0,-1
761 7 : if (std::abs(val) > tolerance) {
762 7 : angle_terms[3 * (-1 + 1) + (0 + 1)] = val;
763 : }
764 7 : val = -1.5 * std::pow(std::sin(angle), 2); // 1,1
765 7 : if (std::abs(val) > tolerance) {
766 7 : angle_terms[3 * (1 + 1) + (1 + 1)] = val;
767 : }
768 7 : val = -1.5 * std::pow(std::sin(angle), 2); // -1,-1
769 7 : if (std::abs(val) > tolerance) {
770 7 : angle_terms[3 * (-1 + 1) + (-1 + 1)] = val;
771 : }
772 :
773 : // Determine which interaction matrices have to be calculated
774 49 : for (const auto &a : angle_terms) {
775 42 : if (interaction_angulardipole.find(a.first) == interaction_angulardipole.end()) {
776 42 : interaction_angulardipole_keys.insert(a.first);
777 : }
778 : }
779 : } else {
780 : // Determine which interaction matrices have to be calculated
781 345 : for (unsigned int order = 3; order <= ordermax; ++order) {
782 196 : if (interaction_multipole.find(order) == interaction_multipole.end()) {
783 93 : interaction_multipole_keys.push_back(order);
784 : }
785 : }
786 : }
787 :
788 : // Is there something to do?
789 198 : bool interaction_multipole_uncalculated = !interaction_multipole_keys.empty();
790 198 : bool interaction_greentensor_dd_uncalculated = !interaction_greentensor_dd_keys.empty();
791 198 : bool interaction_greentensor_dq_uncalculated = !interaction_greentensor_dq_keys.empty();
792 198 : bool interaction_greentensor_qd_uncalculated = !interaction_greentensor_qd_keys.empty();
793 198 : bool interaction_angulardipole_uncalculated = !interaction_angulardipole_keys.empty();
794 :
795 : // Return if there is nothing to do
796 198 : if (!interaction_multipole_uncalculated && !interaction_greentensor_dd_uncalculated &&
797 110 : !interaction_greentensor_dq_uncalculated && !interaction_greentensor_qd_uncalculated &&
798 110 : !interaction_angulardipole_uncalculated) {
799 103 : return;
800 : }
801 :
802 : // Precalculate matrix elements
803 190 : auto states1 = this->getStatesFirst();
804 190 : auto states2 = this->getStatesSecond();
805 238 : for (unsigned int kappa = 1; kappa <= ordermax - 2; ++kappa) {
806 143 : cache.precalculateMultipole(states1, kappa);
807 143 : cache.precalculateMultipole(states2, kappa); // TODO check whether system1 == system2
808 : }
809 :
810 : ////////////////////////////////////////////////////////////////////
811 : /// Generate the interaction in the canonical basis ////////////////
812 : ////////////////////////////////////////////////////////////////////
813 :
814 : std::unordered_map<int, std::vector<Eigen::Triplet<double>>>
815 190 : interaction_angulardipole_triplets; // TODO reserve
816 : std::unordered_map<int, std::vector<Eigen::Triplet<double>>>
817 190 : interaction_multipole_triplets; // TODO reserve
818 : std::unordered_map<int, std::vector<Eigen::Triplet<Scalar>>>
819 190 : interaction_greentensor_dd_triplets; // TODO reserve
820 : std::unordered_map<int, std::vector<Eigen::Triplet<Scalar>>>
821 190 : interaction_greentensor_dq_triplets; // TODO reserve
822 : std::unordered_map<int, std::vector<Eigen::Triplet<Scalar>>>
823 190 : interaction_greentensor_qd_triplets; // TODO reserve
824 :
825 : /*// Categorize states // TODO
826 : std::unordered_map<ljm_t, std::vector<enumerated_state>> states_ordered;
827 : for (const auto &s : states) {
828 : states_categorized[ljm_t(s.getL(), s.getJ(), s.getM())].push_back(s)
829 : }
830 :
831 : // Reserve vectors // TODO
832 : for (const auto &state : states) {
833 : for (const auto &order : orange) {
834 : auto list_of_coupled_ljm = getCoupled(state, order); //, q, kappa
835 : for (const auto &coupled : list_of_coupled_ljm) {
836 : num_elements[order] += states_categorized[coupled].size();
837 : }
838 : }
839 : }
840 : for (const auto &order : orange) {
841 : interaction_multipole_triplets[order].reserve(num_elements[order]);
842 : }
843 : // TODO
844 : for (size_t i = 0; i < calculation_required.size(); ++i) {
845 : // TODO
846 : }*/
847 :
848 : // Loop over column entries
849 29167 : for (const auto &c : this->states) { // TODO parallelization
850 14536 : if (c.state.isArtificial(0) || c.state.isArtificial(1)) {
851 0 : continue;
852 : }
853 :
854 : // Loop over row entries
855 6220492 : for (const auto &r : this->states) {
856 3102978 : if (r.state.isArtificial(0) || r.state.isArtificial(1)) {
857 0 : continue;
858 : }
859 3102978 : if (r.idx < c.idx) {
860 1544221 : continue;
861 : }
862 :
863 1558757 : int q1 = r.state.getM(0) - c.state.getM(0);
864 1558757 : int q2 = r.state.getM(1) - c.state.getM(1);
865 :
866 : // Green tensor approach
867 :
868 : // Begin of dipole-dipole interaction
869 1558757 : if (interaction_greentensor_dd_uncalculated) {
870 767067 : if (selectionRulesMultipoleNew(r.state.getFirstState(), c.state.getFirstState(),
871 862226 : 1) &&
872 95159 : selectionRulesMultipoleNew(r.state.getSecondState(), c.state.getSecondState(),
873 : 1)) { // dipole-dipole interaction
874 :
875 451428 : for (const auto &key : interaction_greentensor_dd_keys) {
876 372171 : int i1 = key % 3;
877 372171 : int i2 = key / 3;
878 :
879 : // Dipole vector of first atom
880 372171 : double vec1_entry = 0;
881 372171 : if (i1 == 0 && std::abs(q1) == 1) { // Calculate x-entry
882 87876 : vec1_entry = 1 / std::sqrt(2.) *
883 87876 : cache.getElectricDipole(r.state.getFirstState(),
884 : c.state.getFirstState());
885 87876 : vec1_entry *= (q1 == 1) ? -1 : 1;
886 284295 : } else if (i1 == 1 && std::abs(q1) == 1) { // Calculate y-entry
887 45156 : vec1_entry = 1 / std::sqrt(2.) *
888 45156 : cache.getElectricDipole(r.state.getFirstState(),
889 : c.state.getFirstState());
890 239139 : } else if (i1 == 2 && q1 == 0) { // Calculate z-entry
891 60181 : vec1_entry = cache.getElectricDipole(r.state.getFirstState(),
892 : c.state.getFirstState());
893 : } else {
894 178958 : continue;
895 : }
896 :
897 : // Dipole vector of second atom
898 193213 : double vec2_entry = 0;
899 193213 : if (i2 == 0 && std::abs(q2) == 1) { // Calculate x-entry
900 44196 : vec2_entry = 1 / std::sqrt(2.) *
901 44196 : cache.getElectricDipole(r.state.getSecondState(),
902 : c.state.getSecondState());
903 44196 : vec2_entry *= (q2 == 1) ? -1 : 1;
904 149017 : } else if (i2 == 1 && std::abs(q2) == 1) { // Calculate y-entry
905 27108 : vec2_entry = 1 / std::sqrt(2.) *
906 27108 : cache.getElectricDipole(r.state.getSecondState(),
907 : c.state.getSecondState());
908 121909 : } else if (i2 == 2 && q2 == 0) { // Calculate z-entry
909 29781 : vec2_entry = cache.getElectricDipole(r.state.getSecondState(),
910 : c.state.getSecondState());
911 : } else {
912 92128 : continue;
913 : }
914 :
915 : // Combine everything
916 101085 : Scalar val = -coulombs_constant * vec1_entry * vec2_entry;
917 101085 : if (i1 == 1 && i2 == 1) {
918 27108 : val *= -1; //"-" from i^2
919 : }
920 101085 : if (key % 2 != 0) {
921 0 : val *= utils::imaginary_unit<Scalar>();
922 : }
923 101085 : this->addTriplet(interaction_greentensor_dd_triplets[3 * i2 + i1], r.idx,
924 101085 : c.idx, val);
925 : }
926 : }
927 : } // End of dipole-dipole interaction
928 :
929 : // Begin of dipole-quadrupole interaction
930 1558757 : if (interaction_greentensor_dq_uncalculated) {
931 561 : if (((selectionRulesMultipoleNew(r.state.getFirstState(), c.state.getFirstState(),
932 119 : 1) &&
933 119 : selectionRulesMultipoleNew(r.state.getSecondState(), c.state.getSecondState(),
934 1122 : 2)) ||
935 525 : (selectionRulesMultipoleNew(r.state.getFirstState(), c.state.getFirstState(),
936 83 : 1) &&
937 83 : selectionRulesMultipoleNew(r.state.getSecondState(), c.state.getSecondState(),
938 : 0)))) { // dipole-quadrupole
939 288 : for (const auto &key : interaction_greentensor_dq_keys) {
940 252 : int i1 = (key % 9) % 3;
941 252 : int i2 = (key % 9) / 3;
942 252 : int i3 = key / 9;
943 252 : double vec_entry = 0.;
944 252 : double matrix_entry = 0.;
945 : // dipole vector of atom1
946 252 : if (i1 == 0 && std::abs(q1) == 1) {
947 48 : vec_entry = 1. / std::sqrt(2.) *
948 48 : cache.getElectricDipole(r.state.getFirstState(),
949 : c.state.getFirstState());
950 48 : vec_entry *= (q1 == 1) ? -1. : 1.;
951 204 : } else if (i1 == 1 && std::abs(q1) == 1) {
952 48 : vec_entry = 1. / std::sqrt(2.) *
953 48 : cache.getElectricDipole(r.state.getFirstState(),
954 : c.state.getFirstState());
955 156 : } else if (i1 == 2 && q1 == 0) {
956 36 : vec_entry = cache.getElectricDipole(r.state.getFirstState(),
957 : c.state.getFirstState());
958 : } else {
959 120 : continue;
960 : }
961 : // quadrupole matrix of atom2
962 132 : if (i2 == 0 && i3 == 0) {
963 12 : if (std::abs(q2) == 2) {
964 0 : matrix_entry =
965 0 : cache.getElectricMultipole(r.state.getSecondState(),
966 : c.state.getSecondState(), 2) *
967 : std::sqrt(5. / 4.); // with std::sqrt(2kappa +1 /4pi)
968 12 : } else if (q2 == 0) {
969 24 : matrix_entry = -std::sqrt(2. / 3.) *
970 12 : cache.getElectricMultipole(r.state.getSecondState(),
971 12 : c.state.getSecondState(), 2) *
972 : std::sqrt(5. / 4.) +
973 12 : std::sqrt(10. / 3.) *
974 12 : cache.getElectricMultipole(r.state.getSecondState(),
975 12 : c.state.getSecondState(), 0) *
976 : std::sqrt(1. / 4.);
977 : }
978 120 : } else if (i2 == 0 && i3 == 1) {
979 0 : if (std::abs(q2) == 2) {
980 0 : matrix_entry =
981 0 : cache.getElectricMultipole(r.state.getSecondState(),
982 : c.state.getSecondState(), 2) *
983 : std::sqrt(5. / 4.);
984 0 : matrix_entry *= (q2 == 2) ? -1. : 1.;
985 : }
986 120 : } else if (i2 == 0 && i3 == 2) {
987 24 : if (std::abs(q2) == 1) {
988 24 : matrix_entry =
989 24 : cache.getElectricMultipole(r.state.getSecondState(),
990 : c.state.getSecondState(), 2) *
991 : std::sqrt(5. / 4.);
992 24 : matrix_entry *= (q2 == 1) ? -1. : 1.;
993 : }
994 96 : } else if (i2 == 1 && i3 == 0) {
995 0 : if (std::abs(q2) == 2) {
996 0 : matrix_entry =
997 0 : cache.getElectricMultipole(r.state.getSecondState(),
998 : c.state.getSecondState(), 2) *
999 : std::sqrt(5. / 4.);
1000 0 : matrix_entry *= (q2 == 2) ? -1. : 1.;
1001 : }
1002 96 : } else if (i2 == 1 && i3 == 1) {
1003 12 : if (std::abs(q2) == 2) {
1004 0 : matrix_entry =
1005 0 : -cache.getElectricMultipole(r.state.getSecondState(),
1006 0 : c.state.getSecondState(), 2) *
1007 : std::sqrt(5. / 4.); // with std::sqrt(2kappa +1 /4pi)
1008 12 : } else if (q2 == 0) {
1009 24 : matrix_entry = -std::sqrt(2. / 3.) *
1010 12 : cache.getElectricMultipole(r.state.getSecondState(),
1011 12 : c.state.getSecondState(), 2) *
1012 : std::sqrt(5. / 4.) +
1013 12 : std::sqrt(10. / 3.) *
1014 12 : cache.getElectricMultipole(r.state.getSecondState(),
1015 12 : c.state.getSecondState(), 0) *
1016 : std::sqrt(1. / 4.);
1017 : }
1018 84 : } else if (i2 == 1 && i3 == 2) {
1019 24 : if (std::abs(q2) == 1) {
1020 24 : matrix_entry =
1021 24 : cache.getElectricMultipole(r.state.getSecondState(),
1022 : c.state.getSecondState(), 2) *
1023 : std::sqrt(5. / 4.);
1024 : }
1025 60 : } else if (i2 == 2 && i3 == 0) {
1026 24 : if (std::abs(q2) == 1) {
1027 24 : matrix_entry =
1028 24 : cache.getElectricMultipole(r.state.getSecondState(),
1029 : c.state.getSecondState(), 2) *
1030 : std::sqrt(5. / 4.);
1031 24 : matrix_entry *= (q2 == 1) ? -1. : 1.;
1032 : }
1033 :
1034 36 : } else if (i2 == 2 && i3 == 1) {
1035 24 : if (std::abs(q2) == 1) {
1036 24 : matrix_entry =
1037 24 : cache.getElectricMultipole(r.state.getSecondState(),
1038 : c.state.getSecondState(), 2) *
1039 : std::sqrt(5. / 4.);
1040 : }
1041 12 : } else if (i2 == 2 && i3 == 2) {
1042 12 : if (q2 == 0) {
1043 24 : matrix_entry = std::sqrt(8. / 3.) *
1044 12 : cache.getElectricMultipole(r.state.getSecondState(),
1045 12 : c.state.getSecondState(), 2) *
1046 : std::sqrt(5. / 4.) +
1047 12 : std::sqrt(10. / 3.) *
1048 12 : cache.getElectricMultipole(r.state.getSecondState(),
1049 12 : c.state.getSecondState(), 0) *
1050 : std::sqrt(1. / 4.);
1051 : }
1052 : }
1053 132 : if (std::abs(vec_entry * matrix_entry) > tolerance) {
1054 132 : Scalar val =
1055 132 : coulombs_constant * vec_entry * matrix_entry / std::sqrt(30.);
1056 132 : if ((i1 == 1) && (i2 != i3 && (i2 == 1 || i3 == 1))) {
1057 48 : val *= -1.; // from i^2
1058 : }
1059 132 : if (key % 2 != 0) {
1060 0 : val *= utils::imaginary_unit<Scalar>();
1061 : }
1062 132 : this->addTriplet(
1063 132 : interaction_greentensor_dq_triplets[9 * i3 + 3 * i2 + i1], r.idx,
1064 132 : c.idx, val);
1065 : }
1066 : }
1067 : }
1068 : }
1069 :
1070 1558757 : if (interaction_greentensor_qd_uncalculated) {
1071 561 : if (((selectionRulesMultipoleNew(r.state.getFirstState(), c.state.getFirstState(),
1072 271 : 2) &&
1073 271 : selectionRulesMultipoleNew(r.state.getSecondState(), c.state.getSecondState(),
1074 1122 : 1)) ||
1075 525 : (selectionRulesMultipoleNew(r.state.getFirstState(), c.state.getFirstState(),
1076 53 : 0) &&
1077 53 : selectionRulesMultipoleNew(r.state.getSecondState(), c.state.getSecondState(),
1078 : 1)))) {
1079 :
1080 : // quadrupole dipole terms
1081 288 : for (const auto &key : interaction_greentensor_qd_keys) {
1082 252 : int i1 = (key % 9) % 3;
1083 252 : int i2 = (key % 9) / 3;
1084 252 : int i3 = key / 9;
1085 252 : double vec_entry = 0.;
1086 252 : double matrix_entry = 0.;
1087 :
1088 : // dipole vector of atom 2
1089 252 : if (i3 == 0 && std::abs(q2) == 1) {
1090 48 : vec_entry = 1. / std::sqrt(2.) *
1091 48 : cache.getElectricDipole(r.state.getSecondState(),
1092 : c.state.getSecondState());
1093 48 : vec_entry *= (q2 == 1) ? -1. : 1.;
1094 204 : } else if (i3 == 1 && std::abs(q2) == 1) {
1095 48 : vec_entry = 1. / std::sqrt(2.) *
1096 48 : cache.getElectricDipole(r.state.getSecondState(),
1097 : c.state.getSecondState());
1098 156 : } else if (i3 == 2 && q2 == 0) {
1099 36 : vec_entry = cache.getElectricDipole(r.state.getSecondState(),
1100 : c.state.getSecondState());
1101 : } else {
1102 120 : continue;
1103 : }
1104 : // quadrupole matrix of atom 1
1105 132 : if (i1 == 0 && i2 == 0) {
1106 12 : if (std::abs(q1) == 2) {
1107 0 : matrix_entry =
1108 0 : cache.getElectricMultipole(r.state.getFirstState(),
1109 : c.state.getFirstState(), 2) *
1110 : std::sqrt(5. / 4.);
1111 12 : } else if (q1 == 0) {
1112 24 : matrix_entry = -std::sqrt(2. / 3.) *
1113 12 : cache.getElectricMultipole(r.state.getFirstState(),
1114 12 : c.state.getFirstState(), 2) *
1115 : std::sqrt(5. / 4.) +
1116 12 : std::sqrt(10. / 3.) *
1117 12 : cache.getElectricMultipole(r.state.getFirstState(),
1118 12 : c.state.getFirstState(), 0) *
1119 : std::sqrt(1. / 4.);
1120 : }
1121 120 : } else if (i1 == 0 && i2 == 1) {
1122 0 : if (std::abs(q1) == 2) {
1123 0 : matrix_entry =
1124 0 : cache.getElectricMultipole(r.state.getFirstState(),
1125 : c.state.getFirstState(), 2) *
1126 : std::sqrt(5. / 4.);
1127 0 : matrix_entry *= (q1 == 2) ? -1. : 1.;
1128 : }
1129 120 : } else if (i1 == 0 && i2 == 2) {
1130 24 : if (std::abs(q1) == 1) {
1131 24 : matrix_entry =
1132 24 : cache.getElectricMultipole(r.state.getFirstState(),
1133 : c.state.getFirstState(), 2) *
1134 : std::sqrt(5. / 4.);
1135 24 : matrix_entry *= (q1 == 1) ? -1. : 1.;
1136 : }
1137 96 : } else if (i1 == 1 && i2 == 0) {
1138 0 : if (std::abs(q1) == 2) {
1139 0 : matrix_entry =
1140 0 : cache.getElectricMultipole(r.state.getFirstState(),
1141 : c.state.getFirstState(), 2) *
1142 : std::sqrt(5. / 4.);
1143 0 : matrix_entry *= (q1 == 2) ? -1. : 1.;
1144 : }
1145 96 : } else if (i1 == 1 && i2 == 1) {
1146 12 : if (std::abs(q1) == 2) {
1147 0 : matrix_entry =
1148 0 : -cache.getElectricMultipole(r.state.getFirstState(),
1149 0 : c.state.getFirstState(), 2) *
1150 : std::sqrt(5. / 4.);
1151 12 : } else if (q1 == 0) {
1152 24 : matrix_entry = -std::sqrt(2. / 3.) *
1153 12 : cache.getElectricMultipole(r.state.getFirstState(),
1154 12 : c.state.getFirstState(), 2) *
1155 : std::sqrt(5. / 4.) +
1156 12 : std::sqrt(10. / 3.) *
1157 12 : cache.getElectricMultipole(r.state.getFirstState(),
1158 12 : c.state.getFirstState(), 0) *
1159 : std::sqrt(1. / 4.);
1160 : }
1161 84 : } else if (i1 == 1 && i2 == 2) {
1162 24 : if (std::abs(q1) == 1) {
1163 24 : matrix_entry =
1164 24 : cache.getElectricMultipole(r.state.getFirstState(),
1165 : c.state.getFirstState(), 2) *
1166 : std::sqrt(5. / 4.);
1167 : }
1168 60 : } else if (i1 == 2 && i2 == 0) {
1169 24 : if (std::abs(q1) == 1) {
1170 24 : matrix_entry =
1171 24 : cache.getElectricMultipole(r.state.getFirstState(),
1172 : c.state.getFirstState(), 2) *
1173 : std::sqrt(5. / 4.);
1174 24 : matrix_entry *= (q1 == 1) ? -1. : 1.;
1175 : }
1176 :
1177 36 : } else if (i1 == 2 && i2 == 1) {
1178 24 : if (std::abs(q1) == 1) {
1179 24 : matrix_entry =
1180 24 : cache.getElectricMultipole(r.state.getFirstState(),
1181 : c.state.getFirstState(), 2) *
1182 : std::sqrt(5. / 4.);
1183 : }
1184 12 : } else if (i1 == 2 && i2 == 2) {
1185 12 : if (q1 == 0) {
1186 24 : matrix_entry = std::sqrt(8. / 3.) *
1187 12 : cache.getElectricMultipole(r.state.getFirstState(),
1188 12 : c.state.getFirstState(), 2) *
1189 : std::sqrt(5. / 4.) +
1190 12 : std::sqrt(10. / 3.) *
1191 12 : cache.getElectricMultipole(r.state.getFirstState(),
1192 12 : c.state.getFirstState(), 0) *
1193 : std::sqrt(1. / 4.);
1194 : }
1195 : }
1196 132 : if (std::abs(vec_entry * matrix_entry) > tolerance) {
1197 132 : Scalar val =
1198 132 : coulombs_constant * vec_entry * matrix_entry / std::sqrt(30.);
1199 132 : if ((i3 == 1) && (i1 != i2 && (i1 == 1 || i2 == 1))) {
1200 48 : val *= -1.; // from i^2
1201 : }
1202 132 : if (key % 2 != 0) {
1203 0 : val *= utils::imaginary_unit<Scalar>();
1204 : }
1205 132 : this->addTriplet(
1206 132 : interaction_greentensor_qd_triplets[9 * i3 + 3 * i2 + i1], r.idx,
1207 132 : c.idx, val);
1208 : }
1209 : }
1210 : }
1211 : } // End of dipole-quadrupole interaction
1212 :
1213 : // Angular dependent dipole-dipole interaction (setAngle and setOrder take care that
1214 : // a non-zero angle cannot occur for other interaction than dipole-dipole if the
1215 : // Green tensor approach is not used)
1216 1558757 : if (interaction_angulardipole_uncalculated) {
1217 276250 : if (selectionRulesMultipoleNew(r.state.getFirstState(), c.state.getFirstState(),
1218 325890 : 1) &&
1219 49640 : selectionRulesMultipoleNew(r.state.getSecondState(), c.state.getSecondState(),
1220 : 1)) {
1221 20432 : auto key = (q1 < q2) ? 3 * (q1 + 1) + (q2 + 1) : 3 * (q2 + 1) + (q1 + 1);
1222 :
1223 20432 : if (interaction_angulardipole_keys.count(key) != 0) {
1224 40864 : double val = coulombs_constant *
1225 20432 : cache.getElectricDipole(r.state.getFirstState(),
1226 : c.state.getFirstState()) *
1227 20432 : cache.getElectricDipole(r.state.getSecondState(),
1228 : c.state.getSecondState());
1229 :
1230 20432 : this->addTriplet(interaction_angulardipole_triplets[key], r.idx, c.idx,
1231 : val);
1232 : }
1233 : }
1234 : }
1235 :
1236 : // Multipole interaction
1237 1558757 : if (interaction_multipole_uncalculated) {
1238 515440 : if (q1 + q2 == 0) {
1239 328609 : for (const auto &order : interaction_multipole_keys) {
1240 :
1241 204896 : double val = 0;
1242 531286 : for (int kappa1 = 1; kappa1 <= order - 2; ++kappa1) {
1243 :
1244 326390 : int kappa2 = order - 1 - kappa1;
1245 326390 : if (selectionRulesMultipoleNew(r.state.getFirstState(),
1246 433245 : c.state.getFirstState(), kappa1) &&
1247 106855 : selectionRulesMultipoleNew(r.state.getSecondState(),
1248 : c.state.getSecondState(), kappa2)) {
1249 :
1250 125631 : double binomials = boost::math::binomial_coefficient<double>(
1251 41877 : kappa1 + kappa2, kappa1 + q1) *
1252 83754 : boost::math::binomial_coefficient<double>(kappa1 + kappa2,
1253 41877 : kappa2 - q2);
1254 :
1255 41877 : val += coulombs_constant * std::pow(-1, kappa2) *
1256 83754 : std::sqrt(binomials) *
1257 41877 : cache.getElectricMultipole(r.state.getFirstState(),
1258 41877 : c.state.getFirstState(), kappa1) *
1259 41877 : cache.getElectricMultipole(r.state.getSecondState(),
1260 : c.state.getSecondState(), kappa2);
1261 : }
1262 : }
1263 :
1264 204896 : this->addTriplet(interaction_multipole_triplets[order], r.idx, c.idx, val);
1265 : }
1266 : }
1267 : }
1268 : }
1269 : }
1270 :
1271 : ////////////////////////////////////////////////////////////////////
1272 : /// Build and transform the interaction to the used basis //////////
1273 : ////////////////////////////////////////////////////////////////////
1274 :
1275 : // Build the interaction and change it from the canonical to the symmetrized basis
1276 :
1277 291 : for (const auto &i : interaction_greentensor_dd_keys) {
1278 196 : interaction_greentensor_dd[i].resize(this->states.size(), this->states.size());
1279 392 : interaction_greentensor_dd[i].setFromTriplets(
1280 196 : interaction_greentensor_dd_triplets[i].begin(),
1281 196 : interaction_greentensor_dd_triplets[i].end());
1282 196 : interaction_greentensor_dd_triplets[i].clear();
1283 :
1284 392 : interaction_greentensor_dd[i] = this->basisvectors.adjoint() *
1285 588 : interaction_greentensor_dd[i].template selfadjointView<Eigen::Lower>() *
1286 : this->basisvectors;
1287 : }
1288 102 : for (const auto &i : interaction_greentensor_dq_keys) {
1289 7 : interaction_greentensor_dq[i].resize(this->states.size(), this->states.size());
1290 14 : interaction_greentensor_dq[i].setFromTriplets(
1291 7 : interaction_greentensor_dq_triplets[i].begin(),
1292 7 : interaction_greentensor_dq_triplets[i].end());
1293 7 : interaction_greentensor_dq_triplets[i].clear();
1294 :
1295 14 : interaction_greentensor_dq[i] = this->basisvectors.adjoint() *
1296 21 : interaction_greentensor_dq[i].template selfadjointView<Eigen::Lower>() *
1297 : this->basisvectors;
1298 : }
1299 102 : for (const auto &i : interaction_greentensor_qd_keys) {
1300 7 : interaction_greentensor_qd[i].resize(this->states.size(), this->states.size());
1301 14 : interaction_greentensor_qd[i].setFromTriplets(
1302 7 : interaction_greentensor_qd_triplets[i].begin(),
1303 7 : interaction_greentensor_qd_triplets[i].end());
1304 7 : interaction_greentensor_qd_triplets[i].clear();
1305 :
1306 14 : interaction_greentensor_qd[i] = this->basisvectors.adjoint() *
1307 21 : interaction_greentensor_qd[i].template selfadjointView<Eigen::Lower>() *
1308 : this->basisvectors;
1309 : }
1310 137 : for (const auto &i : interaction_angulardipole_keys) {
1311 42 : interaction_angulardipole[i].resize(this->states.size(), this->states.size());
1312 84 : interaction_angulardipole[i].setFromTriplets(interaction_angulardipole_triplets[i].begin(),
1313 42 : interaction_angulardipole_triplets[i].end());
1314 42 : interaction_angulardipole_triplets[i].clear();
1315 :
1316 84 : interaction_angulardipole[i] = this->basisvectors.adjoint() *
1317 126 : interaction_angulardipole[i].template selfadjointView<Eigen::Lower>() *
1318 : this->basisvectors;
1319 : }
1320 :
1321 188 : for (const auto &i : interaction_multipole_keys) {
1322 93 : interaction_multipole[i].resize(this->states.size(), this->states.size());
1323 186 : interaction_multipole[i].setFromTriplets(interaction_multipole_triplets[i].begin(),
1324 93 : interaction_multipole_triplets[i].end());
1325 93 : interaction_multipole_triplets[i].clear();
1326 :
1327 186 : interaction_multipole[i] = this->basisvectors.adjoint() *
1328 279 : interaction_multipole[i].template selfadjointView<Eigen::Lower>() * this->basisvectors;
1329 : }
1330 : }
1331 :
1332 : ////////////////////////////////////////////////////////////////////
1333 : /// Method that allows base class to construct Hamiltonian /////////
1334 : ////////////////////////////////////////////////////////////////////
1335 :
1336 : template <typename Scalar>
1337 192 : void SystemTwo<Scalar>::addInteraction() {
1338 192 : if (distance == std::numeric_limits<double>::max()) {
1339 0 : return;
1340 : }
1341 :
1342 : // Build the total Hamiltonian
1343 192 : if (GTbool) {
1344 236 : for (const auto &g : greentensor_terms_dd) {
1345 195 : this->hamiltonian += interaction_greentensor_dd[g.first] * g.second;
1346 : }
1347 48 : for (const auto &g : greentensor_terms_dq) {
1348 7 : this->hamiltonian += interaction_greentensor_dq[g.first] * g.second;
1349 : }
1350 48 : for (const auto &g : greentensor_terms_qd) {
1351 7 : this->hamiltonian += interaction_greentensor_qd[g.first] * g.second;
1352 : }
1353 151 : } else if (distance_x != 0) {
1354 7 : double powerlaw = 1. / std::pow(distance, 3);
1355 49 : for (const auto &a : angle_terms) {
1356 42 : this->hamiltonian += interaction_angulardipole[a.first] * a.second * powerlaw;
1357 : }
1358 : } else {
1359 335 : for (unsigned int order = 3; order <= ordermax; ++order) {
1360 191 : double powerlaw = 1. / std::pow(distance, order);
1361 191 : this->hamiltonian += interaction_multipole[order] * powerlaw;
1362 : }
1363 : }
1364 : }
1365 :
1366 : ////////////////////////////////////////////////////////////////////
1367 : /// Method that allows base class to transform the interaction /////
1368 : ////////////////////////////////////////////////////////////////////
1369 :
1370 : template <typename Scalar>
1371 85 : void SystemTwo<Scalar>::transformInteraction(const Eigen::SparseMatrix<Scalar> &transformator) {
1372 85 : for (auto &entry : interaction_greentensor_dd) {
1373 0 : entry.second = transformator.adjoint() * entry.second * transformator;
1374 : }
1375 85 : for (auto &entry : interaction_greentensor_dq) {
1376 0 : entry.second = transformator.adjoint() * entry.second * transformator;
1377 : }
1378 85 : for (auto &entry : interaction_greentensor_qd) {
1379 0 : entry.second = transformator.adjoint() * entry.second * transformator;
1380 : }
1381 97 : for (auto &entry : interaction_angulardipole) {
1382 12 : entry.second = transformator.adjoint() * entry.second * transformator;
1383 : }
1384 89 : for (auto &entry : interaction_multipole) {
1385 4 : entry.second = transformator.adjoint() * entry.second * transformator; // NOLINT
1386 : }
1387 85 : }
1388 :
1389 : ////////////////////////////////////////////////////////////////////
1390 : /// Method that allows base class to delete the interaction ////////
1391 : ////////////////////////////////////////////////////////////////////
1392 :
1393 : template <typename Scalar>
1394 33 : void SystemTwo<Scalar>::deleteInteraction() {
1395 33 : interaction_greentensor_dd.clear();
1396 33 : interaction_greentensor_dq.clear();
1397 33 : interaction_greentensor_qd.clear();
1398 33 : interaction_angulardipole.clear();
1399 33 : interaction_multipole.clear();
1400 33 : }
1401 :
1402 : ////////////////////////////////////////////////////////////////////
1403 : /// Methods that allows base class to rotate states ////////////////
1404 : ////////////////////////////////////////////////////////////////////
1405 :
1406 : template <typename Scalar>
1407 : Eigen::SparseMatrix<Scalar>
1408 39 : SystemTwo<Scalar>::rotateStates(const std::vector<size_t> &states_indices, double alpha,
1409 : double beta, double gamma) {
1410 : // Initialize Wigner D matrix
1411 39 : WignerD wigner;
1412 :
1413 : // Rotate state
1414 78 : std::vector<Eigen::Triplet<Scalar>> states_rotated_triplets;
1415 39 : states_rotated_triplets.reserve(
1416 39 : std::min(static_cast<size_t>(10 * 10), this->states.size()) *
1417 39 : states_indices.size()); // TODO std::min( (2*jmax+1)*(2*jmax+1), states.size() ) *
1418 : // states_indices.size()
1419 :
1420 39 : size_t current = 0;
1421 78 : for (auto const &idx : states_indices) {
1422 39 : this->addRotated(this->states[idx].state, current++, states_rotated_triplets, wigner, alpha,
1423 : beta, gamma);
1424 : }
1425 :
1426 39 : Eigen::SparseMatrix<Scalar> states_rotated(this->states.size(), states_indices.size());
1427 39 : states_rotated.setFromTriplets(states_rotated_triplets.begin(), states_rotated_triplets.end());
1428 39 : states_rotated_triplets.clear();
1429 :
1430 78 : return states_rotated;
1431 : }
1432 :
1433 : template <typename Scalar>
1434 1 : Eigen::SparseMatrix<Scalar> SystemTwo<Scalar>::buildStaterotator(double alpha, double beta,
1435 : double gamma) {
1436 : // Initialize Wigner D matrix
1437 1 : WignerD wigner;
1438 :
1439 : // Build rotator
1440 2 : std::vector<Eigen::Triplet<Scalar>> rotator_triplets;
1441 1 : rotator_triplets.reserve(
1442 1 : std::min(static_cast<size_t>(10 * 10), this->states.size()) *
1443 : this->states
1444 1 : .size()); // TODO std::min( (2*jmax+1)*(2*jmax+1), states.size() ) * states.size()
1445 :
1446 537 : for (auto const &entry : this->states) {
1447 268 : this->addRotated(entry.state, entry.idx, rotator_triplets, wigner, alpha, beta, gamma);
1448 : }
1449 :
1450 1 : Eigen::SparseMatrix<Scalar> rotator(this->states.size(), this->states.size());
1451 1 : rotator.setFromTriplets(rotator_triplets.begin(), rotator_triplets.end()); // NOLINT
1452 1 : rotator_triplets.clear();
1453 :
1454 2 : return rotator;
1455 : }
1456 :
1457 : ////////////////////////////////////////////////////////////////////
1458 : /// Method that allows base class to combine systems ///////////////
1459 : ////////////////////////////////////////////////////////////////////
1460 :
1461 : template <typename Scalar>
1462 33 : void SystemTwo<Scalar>::incorporate(SystemBase<Scalar, StateTwo> &system) {
1463 : // Combine parameters
1464 33 : if (species[0] != dynamic_cast<SystemTwo &>(system).species[0]) {
1465 0 : throw std::runtime_error(
1466 : "The value of the variable 'element' must be the same for both systems.");
1467 : }
1468 33 : if (species[1] != dynamic_cast<SystemTwo &>(system).species[1]) {
1469 0 : throw std::runtime_error(
1470 : "The value of the variable 'element' must be the same for both systems.");
1471 : }
1472 33 : if (distance_x != dynamic_cast<SystemTwo &>(system).distance_x) {
1473 0 : throw std::runtime_error(
1474 : "The value of the variable 'distance' must be the same for both systems.");
1475 : }
1476 33 : if (distance_y != dynamic_cast<SystemTwo &>(system).distance_y) {
1477 0 : throw std::runtime_error(
1478 : "The value of the variable 'distance' must be the same for both systems.");
1479 : }
1480 33 : if (distance_z != dynamic_cast<SystemTwo &>(system).distance_z) {
1481 0 : throw std::runtime_error(
1482 : "The value of the variable 'distance' must be the same for both systems.");
1483 : }
1484 33 : if (ordermax != dynamic_cast<SystemTwo &>(system).ordermax) {
1485 0 : throw std::runtime_error(
1486 : "The value of the variable 'ordermax' must be the same for both systems.");
1487 : }
1488 :
1489 : // Combine symmetries
1490 33 : unsigned int num_different_symmetries = 0;
1491 33 : if (sym_permutation != dynamic_cast<SystemTwo &>(system).sym_permutation) {
1492 8 : sym_permutation = NA;
1493 8 : ++num_different_symmetries;
1494 : }
1495 33 : if (sym_inversion != dynamic_cast<SystemTwo &>(system).sym_inversion) {
1496 7 : sym_inversion = NA;
1497 7 : ++num_different_symmetries;
1498 : }
1499 33 : if (sym_reflection != dynamic_cast<SystemTwo &>(system).sym_reflection) {
1500 4 : sym_reflection = NA;
1501 4 : ++num_different_symmetries;
1502 : }
1503 61 : if (!(sym_rotation.size() == dynamic_cast<SystemTwo &>(system).sym_rotation.size() &&
1504 28 : std::equal(sym_rotation.begin(), sym_rotation.end(),
1505 28 : dynamic_cast<SystemTwo &>(system).sym_rotation.begin()))) {
1506 12 : if (sym_rotation.count(ARB) != 0 ||
1507 6 : dynamic_cast<SystemTwo &>(system).sym_rotation.count(ARB) != 0) {
1508 0 : sym_rotation = {ARB};
1509 : } else {
1510 6 : sym_rotation.insert(dynamic_cast<SystemTwo &>(system).sym_rotation.begin(),
1511 6 : dynamic_cast<SystemTwo &>(system).sym_rotation.end());
1512 : }
1513 6 : ++num_different_symmetries;
1514 : }
1515 : // if (num_different_symmetries > 1) std::cerr << "Warning: The systems differ in more than one
1516 : // symmetry. For the combined system, the notion of symmetries might be meaningless." <<
1517 : // std::endl; // TODO let the user enable/disable this warning
1518 :
1519 : // Clear cached interaction
1520 33 : this->deleteInteraction();
1521 33 : }
1522 :
1523 : ////////////////////////////////////////////////////////////////////
1524 : /// Methods that allows base class to communicate with subclass ////
1525 : ////////////////////////////////////////////////////////////////////
1526 :
1527 : template <typename Scalar>
1528 0 : void SystemTwo<Scalar>::onStatesChange() {
1529 0 : minimal_le_roy_radius = std::numeric_limits<double>::max();
1530 0 : }
1531 :
1532 : ////////////////////////////////////////////////////////////////////
1533 : /// Utility methods ////////////////////////////////////////////////
1534 : ////////////////////////////////////////////////////////////////////
1535 :
1536 : template <typename Scalar>
1537 198 : void SystemTwo<Scalar>::checkDistance(const double &distance) {
1538 :
1539 198 : auto &cache = *this->m_cache;
1540 :
1541 : // Get the minimal Le Roy radius
1542 198 : if (minimal_le_roy_radius == std::numeric_limits<double>::max()) {
1543 :
1544 : // Estimate minimal Le Roy radius
1545 177 : StateTwo crucial_state{{{"None", "None"}}};
1546 15019 : for (const auto &e : this->states) {
1547 7480 : if (e.state.isArtificial(0) || e.state.isArtificial(1)) {
1548 0 : continue;
1549 : }
1550 :
1551 7480 : auto n = e.state.getNStar(cache);
1552 7480 : auto l = e.state.getL();
1553 :
1554 7480 : double le_roy_radius = 2 * au2um *
1555 7480 : (std::sqrt(0.5 * n[0] * n[0] * (5 * n[0] * n[0] + 1 - 3 * l[0] * (l[0] + 1))) +
1556 7480 : std::sqrt(0.5 * n[1] * n[1] * (5 * n[1] * n[1] + 1 - 3 * l[1] * (l[1] + 1))));
1557 :
1558 7480 : if (le_roy_radius < minimal_le_roy_radius) {
1559 146 : minimal_le_roy_radius = le_roy_radius;
1560 146 : crucial_state = e.state;
1561 : }
1562 : }
1563 :
1564 : // Calculate minimal Le Roy radius precisely
1565 59 : if (crucial_state.isArtificial(0) || crucial_state.isArtificial(1)) {
1566 0 : minimal_le_roy_radius = 0;
1567 : } else {
1568 59 : minimal_le_roy_radius = crucial_state.getLeRoyRadius(cache);
1569 : }
1570 : }
1571 :
1572 198 : if (distance < minimal_le_roy_radius) {
1573 1 : std::cerr << "WARNING: The distance " << distance
1574 1 : << " um is smaller than the Le Roy radius " << minimal_le_roy_radius << " um."
1575 1 : << std::endl;
1576 : }
1577 198 : }
1578 :
1579 : template <typename Scalar>
1580 2148406 : void SystemTwo<Scalar>::addBasisvectors(const StateTwo &state, const size_t &col_new,
1581 : const Scalar &value_new,
1582 : std::vector<Eigen::Triplet<Scalar>> &basisvectors_triplets,
1583 : std::vector<double> &sqnorm_list) {
1584 2148406 : auto state_iter = this->states.template get<1>().find(state);
1585 :
1586 : size_t row_new;
1587 2148406 : if (state_iter != this->states.template get<1>().end()) {
1588 2133444 : row_new = state_iter->idx;
1589 : } else {
1590 14962 : row_new = this->states.size();
1591 14962 : this->states.emplace_back(row_new, state);
1592 : }
1593 :
1594 2148406 : basisvectors_triplets.emplace_back(row_new, col_new, value_new);
1595 2148406 : sqnorm_list[row_new] += std::pow(std::abs(value_new), 2);
1596 2148406 : }
1597 :
1598 : template <typename Scalar>
1599 26 : bool SystemTwo<Scalar>::isRefelectionAndRotationCompatible() {
1600 26 : if (sym_rotation.count(ARB) != 0 || sym_reflection == NA) {
1601 26 : return true;
1602 : }
1603 :
1604 0 : for (const auto &s : sym_rotation) {
1605 0 : if (sym_rotation.count(-s) == 0) {
1606 0 : return false;
1607 : }
1608 : }
1609 :
1610 0 : return true;
1611 : }
1612 :
1613 : template class SystemTwo<std::complex<double>>;
1614 : template class SystemTwo<double>;
|
