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 "SystemBase.hpp"
21 : #include "Symmetry.hpp"
22 :
23 : #include "EigenCompat.hpp"
24 : #include <Eigen/Geometry>
25 : #include <unsupported/Eigen/MatrixFunctions>
26 :
27 : #include <numeric>
28 :
29 : #ifdef WITH_INTEL_MKL
30 : #define MKL_Complex8 std::complex<float>
31 : #define MKL_Complex16 std::complex<double>
32 : #include <mkl_lapacke.h>
33 : #include <mkl_solvers_ee.h>
34 : #elif defined EIGEN_USE_LAPACKE
35 : #include <lapacke.h>
36 : #endif // WITH_INTEL_MKL
37 :
38 : #ifdef WITH_INTEL_MKL
39 0 : static void feast_csrev(const char *uplo, const MKL_INT *n, const double *a, const MKL_INT *ia,
40 : const MKL_INT *ja, MKL_INT *fpm, double *epsout, MKL_INT *loop,
41 : const double *emin, const double *emax, MKL_INT *m0, double *e, double *x,
42 : MKL_INT *m, double *res, MKL_INT *info) {
43 0 : dfeast_scsrev(uplo, n, a, ia, ja, fpm, epsout, loop, emin, emax, m0, e, x, m, res, info);
44 0 : }
45 :
46 2 : static void feast_csrev(const char *uplo, const MKL_INT *n, const MKL_Complex16 *a,
47 : const MKL_INT *ia, const MKL_INT *ja, MKL_INT *fpm, double *epsout,
48 : MKL_INT *loop, const double *emin, const double *emax, MKL_INT *m0,
49 : double *e, MKL_Complex16 *x, MKL_INT *m, double *res, MKL_INT *info) {
50 2 : zfeast_hcsrev(uplo, n, a, ia, ja, fpm, epsout, loop, emin, emax, m0, e, x, m, res, info);
51 2 : }
52 : #endif // WITH_INTEL_MKL
53 :
54 : #if defined EIGEN_USE_LAPACKE || WITH_INTEL_MKL
55 259 : static int LAPACKE_evd(const int matrix_layout, const char jobz, const char uplo,
56 : const lapack_int n, double *a, const lapack_int lda, double *w) {
57 259 : return LAPACKE_dsyevd(matrix_layout, jobz, uplo, n, a, lda, w);
58 : }
59 :
60 74 : static int LAPACKE_evd(const int matrix_layout, const char jobz, const char uplo,
61 : const lapack_int n, lapack_complex_double *a, const lapack_int lda,
62 : double *w) {
63 74 : return LAPACKE_zheevd(matrix_layout, jobz, uplo, n, a, lda, w);
64 : }
65 : #endif // EIGEN_USE_LAPACKE || WITH_INTEL_MKL
66 :
67 : template <typename State>
68 : State createStateFromLabel(const std::string &label);
69 : template <>
70 0 : StateOne createStateFromLabel<StateOne>(const std::string &label) {
71 0 : return StateOne(label);
72 : }
73 : template <>
74 28 : StateTwo createStateFromLabel<StateTwo>(const std::string &label) {
75 56 : return StateTwo(std::array<std::string, 2>({{"0_" + label, "1_" + label}}));
76 : }
77 :
78 : // Implementation of SystemBase
79 :
80 : template <class Scalar, class State>
81 8 : void SystemBase<Scalar, State>::setMinimalNorm(const double &threshold) {
82 8 : threshold_for_sqnorm = threshold;
83 8 : }
84 :
85 : ////////////////////////////////////////////////////////////////////
86 : /// Methods to restrict the number of states inside the basis //////
87 : ////////////////////////////////////////////////////////////////////
88 :
89 : template <class Scalar, class State>
90 88 : void SystemBase<Scalar, State>::restrictEnergy(double e_min,
91 : double e_max) { // restricts diagonal entries
92 88 : energy_min = e_min;
93 88 : energy_max = e_max;
94 88 : }
95 :
96 : template <class Scalar, class State>
97 43 : void SystemBase<Scalar, State>::restrictN(int n_min, int n_max) {
98 43 : this->range(range_n, n_min, n_max);
99 43 : }
100 :
101 : template <class Scalar, class State>
102 150 : void SystemBase<Scalar, State>::restrictN(std::set<int> n) {
103 150 : range_n = n;
104 150 : }
105 :
106 : template <class Scalar, class State>
107 42 : void SystemBase<Scalar, State>::restrictL(int l_min, int l_max) {
108 42 : this->range(range_l, l_min, l_max);
109 42 : }
110 :
111 : template <class Scalar, class State>
112 150 : void SystemBase<Scalar, State>::restrictL(std::set<int> l) {
113 150 : range_l = l;
114 150 : }
115 :
116 : template <class Scalar, class State>
117 5 : void SystemBase<Scalar, State>::restrictJ(float j_min, float j_max) {
118 5 : this->range(range_j, j_min, j_max);
119 5 : }
120 :
121 : template <class Scalar, class State>
122 150 : void SystemBase<Scalar, State>::restrictJ(std::set<float> j) {
123 150 : range_j = j;
124 150 : }
125 :
126 : template <class Scalar, class State>
127 2 : void SystemBase<Scalar, State>::restrictM(float m_min, float m_max) {
128 2 : this->range(range_m, m_min, m_max);
129 2 : }
130 :
131 : template <class Scalar, class State>
132 150 : void SystemBase<Scalar, State>::restrictM(std::set<float> m) {
133 150 : range_m = m;
134 150 : }
135 :
136 : ////////////////////////////////////////////////////////////////////
137 : /// Method for adding user-defined states //////////////////////////
138 : ////////////////////////////////////////////////////////////////////
139 :
140 : template <class Scalar, class State>
141 0 : void SystemBase<Scalar, State>::addStates(const State &s) {
142 0 : states_to_add.insert(s);
143 0 : }
144 :
145 : template <class Scalar, class State>
146 0 : void SystemBase<Scalar, State>::addStates(const std::set<State> &s) {
147 0 : states_to_add.insert(s.begin(), s.end());
148 0 : }
149 :
150 : // TODO make it possible to just use added states, i.e. use no restrictions on quantum numbers
151 : // and energies
152 :
153 : ////////////////////////////////////////////////////////////////////
154 : /// Methods to get overlaps ////////////////////////////////////////
155 : ////////////////////////////////////////////////////////////////////
156 :
157 : template <class Scalar, class State>
158 209 : Eigen::VectorX<double> SystemBase<Scalar, State>::getOverlap(const State &generalizedstate) {
159 209 : return this->getOverlap(generalizedstate, 0, 0, 0);
160 : }
161 :
162 : template <class Scalar, class State>
163 : Eigen::VectorX<double>
164 3 : SystemBase<Scalar, State>::getOverlap(const std::vector<State> &generalizedstates) {
165 3 : return this->getOverlap(generalizedstates, 0, 0, 0);
166 : }
167 :
168 : template <class Scalar, class State>
169 0 : Eigen::VectorX<double> SystemBase<Scalar, State>::getOverlap(const size_t &state_index) {
170 0 : return this->getOverlap(state_index, 0, 0, 0);
171 : }
172 :
173 : template <class Scalar, class State>
174 : Eigen::VectorX<double>
175 0 : SystemBase<Scalar, State>::getOverlap(const std::vector<size_t> &states_indices) {
176 0 : return this->getOverlap(states_indices, 0, 0, 0);
177 : }
178 :
179 : template <class Scalar, class State>
180 0 : Eigen::VectorX<double> SystemBase<Scalar, State>::getOverlap(const State &generalizedstate,
181 : std::array<double, 3> to_z_axis,
182 : std::array<double, 3> to_y_axis) {
183 0 : auto euler_zyz = this->getEulerAngles(to_z_axis, to_y_axis);
184 0 : return this->getOverlap(generalizedstate, euler_zyz[0], euler_zyz[1], euler_zyz[2]);
185 : }
186 :
187 : template <class Scalar, class State>
188 : Eigen::VectorX<double>
189 0 : SystemBase<Scalar, State>::getOverlap(const std::vector<State> &generalizedstates,
190 : std::array<double, 3> to_z_axis,
191 : std::array<double, 3> to_y_axis) {
192 0 : auto euler_zyz = this->getEulerAngles(to_z_axis, to_y_axis);
193 0 : return this->getOverlap(generalizedstates, euler_zyz[0], euler_zyz[1], euler_zyz[2]);
194 : }
195 :
196 : template <class Scalar, class State>
197 0 : Eigen::VectorX<double> SystemBase<Scalar, State>::getOverlap(const size_t &state_index,
198 : std::array<double, 3> to_z_axis,
199 : std::array<double, 3> to_y_axis) {
200 0 : auto euler_zyz = this->getEulerAngles(to_z_axis, to_y_axis);
201 0 : return this->getOverlap(state_index, euler_zyz[0], euler_zyz[1], euler_zyz[2]);
202 : }
203 :
204 : template <class Scalar, class State>
205 : Eigen::VectorX<double>
206 0 : SystemBase<Scalar, State>::getOverlap(const std::vector<size_t> &states_indices,
207 : std::array<double, 3> to_z_axis,
208 : std::array<double, 3> to_y_axis) {
209 0 : auto euler_zyz = this->getEulerAngles(to_z_axis, to_y_axis);
210 0 : return this->getOverlap(states_indices, euler_zyz[0], euler_zyz[1], euler_zyz[2]);
211 : }
212 :
213 : template <class Scalar, class State>
214 249 : Eigen::VectorX<double> SystemBase<Scalar, State>::getOverlap(const State &generalizedstate,
215 : double alpha, double beta,
216 : double gamma) {
217 996 : std::vector<State> generalizedstates({generalizedstate});
218 498 : return this->getOverlap(generalizedstates, alpha, beta, gamma);
219 : }
220 :
221 : template <class Scalar, class State>
222 0 : Eigen::VectorX<double> SystemBase<Scalar, State>::getOverlap(const size_t &state_index,
223 : double alpha, double beta,
224 : double gamma) {
225 0 : std::vector<size_t> states_indices({state_index});
226 0 : return this->getOverlap(states_indices, alpha, beta, gamma);
227 : }
228 :
229 : template <class Scalar, class State>
230 : Eigen::VectorX<double>
231 254 : SystemBase<Scalar, State>::getOverlap(const std::vector<State> &generalizedstates, double alpha,
232 : double beta, double gamma) {
233 : // Build basis
234 254 : this->buildBasis();
235 :
236 : // Determine indices of the specified states
237 508 : std::vector<size_t> states_indices;
238 254 : states_indices.reserve(generalizedstates.size());
239 :
240 513 : for (const auto &searched_state : generalizedstates) {
241 259 : if (utils::is_true(searched_state.isGeneralized())) {
242 305 : for (const auto &state : states) {
243 150 : if (state.state ^ searched_state) { // Check whether state.state is contained in
244 : // searched_state
245 10 : states_indices.push_back(state.idx);
246 : }
247 : }
248 : } else {
249 254 : auto state_iter = states.template get<1>().find(searched_state);
250 254 : if (state_iter != states.template get<1>().end()) {
251 249 : states_indices.push_back(state_iter->idx);
252 : }
253 : }
254 : }
255 :
256 : // Get the overlap
257 508 : return this->getOverlap(states_indices, alpha, beta, gamma);
258 : }
259 :
260 : template <class Scalar, class State>
261 : Eigen::VectorX<double>
262 254 : SystemBase<Scalar, State>::getOverlap(const std::vector<size_t> &states_indices, double alpha,
263 : double beta, double gamma) {
264 : // Build basis
265 254 : this->buildBasis();
266 :
267 : // Build state vectors
268 508 : Eigen::SparseMatrix<Scalar> overlap_states;
269 254 : if (alpha == 0 && beta == 0 && gamma == 0) {
270 212 : std::vector<Eigen::Triplet<Scalar>> overlap_states_triplets;
271 212 : overlap_states_triplets.reserve(states_indices.size());
272 :
273 212 : size_t current = 0;
274 427 : for (auto const &idx : states_indices) {
275 215 : overlap_states_triplets.emplace_back(idx, current++, 1);
276 : }
277 :
278 212 : overlap_states.resize(states.size(), states_indices.size());
279 212 : overlap_states.setFromTriplets(overlap_states_triplets.begin(),
280 424 : overlap_states_triplets.end());
281 : } else {
282 42 : overlap_states = this->rotateStates(states_indices, alpha, beta, gamma);
283 : }
284 :
285 : // Calculate overlap
286 508 : Eigen::SparseMatrix<Scalar> product = basisvectors.adjoint() * overlap_states;
287 254 : Eigen::VectorX<double> overlap = Eigen::VectorX<double>::Zero(product.rows());
288 513 : for (int k = 0; k < product.outerSize(); ++k) {
289 11689 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(product, k); triple;
290 11430 : ++triple) {
291 11430 : overlap[triple.row()] += std::pow(std::abs(triple.value()), 2);
292 : }
293 : }
294 :
295 508 : return overlap;
296 : }
297 :
298 : ////////////////////////////////////////////////////////////////////
299 : /// Methods to get properties of the system ////////////////////////
300 : ////////////////////////////////////////////////////////////////////
301 :
302 : template <class Scalar, class State>
303 76 : std::vector<State> SystemBase<Scalar, State>::getStates() {
304 76 : this->buildBasis();
305 76 : std::vector<State> states_converted;
306 76 : states_converted.reserve(states.size());
307 11284 : for (const auto &s : states) {
308 5604 : states_converted.push_back(std::move(s.state));
309 : }
310 76 : return states_converted;
311 : }
312 :
313 : template <class Scalar, class State>
314 : const typename states_set<State>::type &
315 1307712 : SystemBase<Scalar, State>::getStatesMultiIndex() { // TODO @hmenke typemap for "const typename
316 : // states_set<T>::type &"
317 1307712 : return states;
318 : }
319 :
320 : template <class Scalar, class State>
321 74543 : Eigen::SparseMatrix<Scalar> &SystemBase<Scalar, State>::getBasisvectors() {
322 74543 : this->buildBasis();
323 74543 : return basisvectors;
324 : }
325 :
326 : template <class Scalar, class State>
327 306548 : Eigen::SparseMatrix<Scalar> &SystemBase<Scalar, State>::getHamiltonian() {
328 306548 : this->buildHamiltonian();
329 306548 : return hamiltonian;
330 : }
331 :
332 : template <class Scalar, class State>
333 295088 : size_t SystemBase<Scalar, State>::getNumBasisvectors() {
334 : // Build basis
335 295088 : this->buildBasis();
336 :
337 : // Check variables for consistency
338 590176 : if ((basisvectors.outerSize() != basisvectors.cols()) ||
339 590176 : (basisvectors.outerSize() != hamiltonian.rows()) ||
340 295088 : (basisvectors.outerSize() != hamiltonian.cols())) {
341 0 : throw std::runtime_error("Inconsistent variables at " + std::string(__FILE__) + ":" +
342 0 : std::to_string(__LINE__) + ".");
343 : }
344 :
345 295088 : return basisvectors.cols();
346 : }
347 :
348 : template <class Scalar, class State>
349 325 : size_t SystemBase<Scalar, State>::getNumStates() {
350 : // Build basis
351 325 : this->buildBasis();
352 :
353 : // Check variables for consistency
354 650 : if ((basisvectors.innerSize() != basisvectors.rows()) ||
355 325 : (static_cast<size_t>(basisvectors.innerSize()) != states.size())) {
356 0 : throw std::runtime_error("Inconsistent variables at " + std::string(__FILE__) + ":" +
357 0 : std::to_string(__LINE__) + ".");
358 : }
359 :
360 325 : return basisvectors.rows();
361 : }
362 :
363 : template <class Scalar, class State>
364 25 : std::vector<State> SystemBase<Scalar, State>::getMainStates() {
365 : // Build basis
366 25 : this->buildBasis();
367 :
368 : // Get states with the main contribution
369 25 : std::vector<State> states_with_maxval;
370 25 : states_with_maxval.reserve(basisvectors.cols());
371 :
372 3383 : for (int k = 0; k < basisvectors.outerSize(); ++k) { // col == idx_vector
373 3358 : double maxval = -1;
374 : size_t row_with_maxval;
375 :
376 27612 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basisvectors, k); triple;
377 24254 : ++triple) {
378 24254 : if (std::abs(triple.value()) > maxval) {
379 4791 : row_with_maxval = triple.row();
380 4791 : maxval = std::abs(triple.value());
381 : }
382 : }
383 :
384 3358 : states_with_maxval.push_back(states[row_with_maxval].state);
385 : }
386 :
387 25 : return states_with_maxval;
388 : }
389 :
390 : template <class Scalar, class State>
391 : std::array<std::vector<size_t>, 2>
392 8 : SystemBase<Scalar, State>::getConnections(SystemBase<Scalar, State> &system_to, double threshold) {
393 : // Build basis
394 8 : this->buildBasis();
395 8 : system_to.buildBasis();
396 :
397 8 : double threshold_sqrt = std::sqrt(threshold);
398 :
399 : // Calculate transformator between the set of states
400 16 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
401 8 : triplets_transformator.reserve(std::min(states.size(), system_to.states.size()));
402 :
403 768 : for (const auto &s : system_to.states) {
404 380 : auto state_iter = states.template get<1>().find(s.state);
405 380 : if (state_iter != states.template get<1>().end()) {
406 346 : size_t idx_from = state_iter->idx;
407 346 : triplets_transformator.emplace_back(idx_from, s.idx, 1);
408 : }
409 : }
410 :
411 16 : Eigen::SparseMatrix<Scalar> transformator(states.size(), system_to.states.size());
412 8 : transformator.setFromTriplets(triplets_transformator.begin(), triplets_transformator.end());
413 :
414 : // Calculate overlap
415 16 : Eigen::SparseMatrix<double> overlap_sqrt =
416 8 : (basisvectors.adjoint() * transformator * system_to.basisvectors).cwiseAbs();
417 :
418 : // Determine the indices of the maximal values within the rows
419 16 : std::vector<size_t> rows_with_maxval;
420 8 : rows_with_maxval.reserve(overlap_sqrt.cols());
421 :
422 168 : for (int k = 0; k < overlap_sqrt.outerSize(); ++k) {
423 160 : double maxval = -1;
424 : size_t row_with_maxval;
425 :
426 1311 : for (typename Eigen::SparseMatrix<double>::InnerIterator triple(overlap_sqrt, k); triple;
427 1151 : ++triple) {
428 1151 : if (triple.value() > maxval) {
429 445 : row_with_maxval = triple.row();
430 445 : maxval = triple.value();
431 : }
432 : }
433 :
434 160 : rows_with_maxval.push_back(row_with_maxval);
435 : }
436 :
437 : // Determine the maximal values within the columns and construct connections
438 16 : Eigen::SparseMatrix<double> overlap_sqrt_transposed = overlap_sqrt.transpose();
439 8 : std::array<std::vector<size_t>, 2> connections;
440 8 : connections[0].reserve(std::max(overlap_sqrt.rows(), overlap_sqrt.cols()));
441 8 : connections[1].reserve(std::max(overlap_sqrt.rows(), overlap_sqrt.cols()));
442 :
443 168 : for (int k = 0; k < overlap_sqrt_transposed.outerSize(); ++k) { // cols
444 160 : double maxval = threshold_sqrt;
445 : size_t idx_from;
446 : size_t idx_to;
447 :
448 1311 : for (typename Eigen::SparseMatrix<double>::InnerIterator triple(overlap_sqrt_transposed, k);
449 1151 : triple; ++triple) {
450 1151 : if (triple.value() > maxval) {
451 335 : idx_from = triple.col();
452 335 : idx_to = triple.row();
453 335 : maxval = triple.value();
454 : }
455 : }
456 :
457 160 : if (maxval > threshold_sqrt && rows_with_maxval[idx_to] == idx_from) {
458 158 : connections[0].push_back(idx_from);
459 158 : connections[1].push_back(idx_to);
460 : }
461 : }
462 :
463 8 : connections[0].shrink_to_fit();
464 8 : connections[1].shrink_to_fit();
465 :
466 16 : return connections;
467 : }
468 :
469 : ////////////////////////////////////////////////////////////////////
470 : /// Methods to build, transform, and destroy the system ////////////
471 : ////////////////////////////////////////////////////////////////////
472 :
473 : template <class Scalar, class State>
474 307047 : void SystemBase<Scalar, State>::buildHamiltonian() {
475 : // Build basis, also constructs the Hamiltonian matrix without interaction
476 307047 : this->buildBasis();
477 :
478 : // Initialize Hamiltonian matrix with interaction if required
479 307047 : if (is_new_hamiltonian_required) {
480 346 : if (is_interaction_already_contained) {
481 :
482 : // Check variables for consistency
483 474 : if (memory_saving || basisvectors_unperturbed_cache.size() == 0 ||
484 237 : hamiltonian_unperturbed_cache.size() == 0) {
485 0 : throw std::runtime_error("Inconsistent variables at " + std::string(__FILE__) +
486 0 : ":" + std::to_string(__LINE__) + ".");
487 : }
488 :
489 : // Reset the Hamiltonian if it already contains the interaction
490 237 : basisvectors = basisvectors_unperturbed_cache;
491 237 : hamiltonian = hamiltonian_unperturbed_cache;
492 109 : } else if (!memory_saving) {
493 :
494 : // Store the Hamiltonian without interaction
495 109 : basisvectors_unperturbed_cache = basisvectors;
496 109 : hamiltonian_unperturbed_cache = hamiltonian;
497 : }
498 :
499 : // Build interaction
500 346 : this->initializeInteraction();
501 :
502 : // Add interaction to the Hamiltonian
503 346 : this->addInteraction();
504 :
505 346 : if (memory_saving) {
506 : // Delete the variables used to add the interaction to the Hamiltonian
507 0 : this->deleteInteraction();
508 : }
509 :
510 346 : is_interaction_already_contained = true;
511 346 : is_new_hamiltonian_required = false;
512 : }
513 307047 : }
514 :
515 : template <class Scalar, class State>
516 9 : void SystemBase<Scalar, State>::buildInteraction() {
517 : // Build basis
518 9 : this->buildBasis();
519 :
520 : // Initialize interaction
521 9 : this->initializeInteraction(); // this method checks by itself whether a new initialization
522 : // is required
523 9 : }
524 :
525 : template <class Scalar, class State>
526 677929 : void SystemBase<Scalar, State>::buildBasis() {
527 : // Check variables for consistency
528 1355858 : if (((hamiltonian.size() == 0) != states.empty()) ||
529 677929 : ((hamiltonian.size() == 0) != (basisvectors.size() == 0))) {
530 0 : throw std::runtime_error("Inconsistent variables at " + std::string(__FILE__) + ":" +
531 0 : std::to_string(__LINE__) + ".");
532 : }
533 :
534 : // In case of no new basis restrictions and already initialized basis, there is nothing to
535 : // do
536 2033459 : if (!states.empty() && states_to_add.empty() && range_n.empty() && range_l.empty() &&
537 2711224 : range_j.empty() && range_m.empty() && energy_min == std::numeric_limits<double>::lowest() &&
538 677765 : energy_max == std::numeric_limits<double>::max()) { // TODO check for new threshold, too
539 677765 : return;
540 : }
541 :
542 : // Check whether the basis does not exist
543 164 : if (states.empty()) {
544 :
545 : // Initialize the basis
546 164 : this->initializeBasis();
547 :
548 : // Forget basis restrictions as they were applied now
549 164 : this->forgetRestrictions();
550 :
551 : } else {
552 0 : this->updateEverything();
553 : }
554 :
555 : // Check whether the basis is empty
556 164 : if (basisvectors.rows() == 0) {
557 0 : throw std::runtime_error("The basis contains no states.");
558 : }
559 164 : if (basisvectors.cols() == 0) {
560 0 : throw std::runtime_error("The basis contains no vectors.");
561 : }
562 : }
563 :
564 : template <class Scalar, class State>
565 1 : void SystemBase<Scalar, State>::diagonalize(double energy_lower_bound, double energy_upper_bound) {
566 1 : this->diagonalize(energy_lower_bound, energy_upper_bound, 0);
567 1 : }
568 :
569 : template <class Scalar, class State>
570 2 : void SystemBase<Scalar, State>::diagonalize(double energy_lower_bound, double energy_upper_bound,
571 : double threshold) {
572 : #ifdef WITH_INTEL_MKL
573 2 : this->buildHamiltonian();
574 :
575 : // Check if already diagonal
576 2 : if (checkIsDiagonal(hamiltonian)) {
577 0 : return;
578 : }
579 :
580 : // Estimate number of found eigenvalues
581 4 : std::vector<double> diagonal_max(hamiltonian.outerSize(), 0);
582 4 : std::vector<double> diagonal_val;
583 4 : std::vector<int> diagonal_idx;
584 2 : diagonal_val.reserve(hamiltonian.outerSize());
585 2 : diagonal_idx.reserve(hamiltonian.outerSize());
586 25 : for (int k = 0; k < hamiltonian.outerSize(); ++k) {
587 23 : double val = 0;
588 268 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(hamiltonian, k); triple;
589 245 : ++triple) {
590 245 : if (triple.row() == triple.col()) {
591 23 : val = std::real(triple.value());
592 222 : } else if (triple.row() != triple.col()) {
593 222 : diagonal_max[k] = std::max(diagonal_max[k], std::abs(triple.value()));
594 : }
595 : }
596 23 : diagonal_idx.push_back(k);
597 23 : diagonal_val.push_back(val);
598 : }
599 :
600 2 : MKL_INT m0 = std::count_if(diagonal_val.begin(), diagonal_val.end(),
601 43 : [&energy_lower_bound, &energy_upper_bound](const double &val) {
602 23 : return (val < energy_upper_bound) && (val > energy_lower_bound);
603 4 : }) +
604 2 : std::count_if(diagonal_idx.begin(), diagonal_idx.end(),
605 100 : [&energy_lower_bound, &energy_upper_bound, &diagonal_val,
606 : &diagonal_max](const int &i) {
607 23 : return ((diagonal_val[i] >= energy_upper_bound) ||
608 20 : (diagonal_val[i] <= energy_lower_bound)) &&
609 45 : (diagonal_val[i] < energy_upper_bound + diagonal_max[i]) &&
610 25 : (diagonal_val[i] > energy_lower_bound - diagonal_max[i]);
611 : });
612 :
613 : // Inplace conversion of the Hamiltonian to CSR with one-based indexing
614 2 : hamiltonian = hamiltonian.transpose();
615 2 : hamiltonian.makeCompressed();
616 27 : for (int i = 0; i < hamiltonian.rows() + 1; ++i) {
617 25 : hamiltonian.outerIndexPtr()[i] += 1;
618 : }
619 247 : for (int i = 0; i < hamiltonian.nonZeros(); ++i) {
620 245 : hamiltonian.innerIndexPtr()[i] += 1;
621 : }
622 :
623 : // Set default parameters for the diagonalization as described at
624 : // https://software.intel.com/en-us/mkl-developer-reference-c-extended-eigensolver-input-parameters.
625 :
626 4 : std::vector<MKL_INT> fpm(128);
627 2 : feastinit(&fpm[0]);
628 2 : fpm[0] = 1; // enables terminal output
629 2 : fpm[1] = 6; // number of contour points
630 2 : fpm[26] = 0; // disables matrix checker
631 2 : fpm[3] = 5; // maximum number of refinement loops allowed
632 2 : if (threshold != 0) {
633 : // Adapt the error trace stopping criteria (10-fpm[2])
634 1 : fpm[2] = std::min(std::round(-std::log10(threshold)), 12.);
635 : }
636 :
637 : // Do the diagonalization
638 : {
639 2 : MKL_INT n = hamiltonian.rows(); // size of the matrix
640 : MKL_INT m; // will contain the number of eigenvalues
641 4 : std::vector<Scalar> x(m0 * n); // the first m columns will contain the eigenvectors
642 : {
643 4 : std::vector<double> e(m0); // will contain the first m eigenvalues
644 : {
645 2 : char uplo = 'F'; // full matrix is stored
646 : MKL_INT info; // will contain return codes
647 : double epsout; // will contain relative error
648 : MKL_INT loop; // will contain number of used refinement
649 4 : std::vector<double> res(m0); // will contain the residual errors
650 :
651 2 : feast_csrev(&uplo, &n, hamiltonian.valuePtr(), hamiltonian.outerIndexPtr(),
652 2 : hamiltonian.innerIndexPtr(), &fpm[0], &epsout, &loop,
653 2 : &energy_lower_bound, &energy_upper_bound, &m0, &e[0], &x[0], &m,
654 2 : &res[0], &info);
655 2 : if (info != 0) {
656 0 : throw std::runtime_error("Diagonalization with FEAST failed.");
657 : }
658 : }
659 :
660 : // Build the new hamiltonian
661 2 : hamiltonian.resize(m, m);
662 2 : hamiltonian.setZero();
663 2 : hamiltonian.reserve(m);
664 20 : for (int idx = 0; idx < m; ++idx) {
665 18 : hamiltonian.insert(idx, idx) = e[idx];
666 : }
667 2 : hamiltonian.makeCompressed();
668 : }
669 :
670 : // Transform the basis vectors
671 4 : Eigen::SparseMatrix<Scalar> evecs =
672 4 : Eigen::Map<Eigen::MatrixX<Scalar>>(&x[0], n, m).sparseView();
673 2 : if (threshold == 0) {
674 1 : basisvectors = basisvectors * evecs;
675 : } else {
676 1 : basisvectors = (basisvectors * evecs).pruned(threshold, 1);
677 : }
678 : }
679 : #else // WITH_INTEL_MKL
680 : (void)energy_lower_bound;
681 : (void)energy_upper_bound;
682 : (void)threshold;
683 : throw std::runtime_error(
684 : "The method does not work because the program was compiled without MKL support.");
685 : #endif // WITH_INTEL_MKL
686 : }
687 :
688 : template <class Scalar, class State>
689 216 : void SystemBase<Scalar, State>::diagonalize() {
690 216 : this->diagonalize(0);
691 216 : }
692 :
693 : template <class Scalar, class State>
694 477 : void SystemBase<Scalar, State>::diagonalize(double threshold) {
695 477 : this->buildHamiltonian();
696 :
697 : // Check if already diagonal
698 477 : if (checkIsDiagonal(hamiltonian)) {
699 144 : return;
700 : }
701 :
702 : #if defined EIGEN_USE_LAPACKE || WITH_INTEL_MKL
703 :
704 : // Diagonalize hamiltonian
705 333 : char jobz = 'V'; // eigenvalues and eigenvectors are computed
706 333 : char uplo = 'U'; // full matrix is stored, upper is used
707 333 : int n = hamiltonian.cols(); // size of the matrix
708 666 : Eigen::MatrixX<Scalar> mat = hamiltonian; // matrix
709 333 : int lda = mat.outerStride(); // leading dimension
710 666 : Eigen::VectorX<double> evals(n); // eigenvalues
711 333 : int info = LAPACKE_evd(LAPACK_COL_MAJOR, jobz, uplo, n, mat.data(), lda, evals.data());
712 333 : if (info != 0) {
713 0 : throw std::runtime_error("Diagonalization with LAPACKE failed.");
714 : }
715 :
716 : // Get eigenvectors
717 666 : Eigen::SparseMatrix<Scalar> evecs = mat.sparseView();
718 :
719 : #else // EIGEN_USE_LAPACKE || WITH_INTEL_MKL
720 :
721 : // Diagonalize hamiltonian
722 : Eigen::SelfAdjointEigenSolver<Eigen::MatrixX<Scalar>> eigensolver(hamiltonian);
723 :
724 : // Get eigenvalues and eigenvectors
725 : Eigen::VectorX<double> evals = eigensolver.eigenvalues();
726 : Eigen::SparseMatrix<Scalar> evecs = eigensolver.eigenvectors().sparseView();
727 :
728 : #endif // EIGEN_USE_LAPACKE || WITH_INTEL_MKL
729 :
730 : // Build the new hamiltonian
731 333 : hamiltonian.setZero();
732 333 : hamiltonian.reserve(evals.size());
733 17284 : for (int idx = 0; idx < evals.size(); ++idx) {
734 16951 : hamiltonian.insert(idx, idx) = evals.coeffRef(idx);
735 : }
736 333 : hamiltonian.makeCompressed();
737 :
738 : // Transform the basis vectors
739 333 : if (threshold == 0) {
740 79 : basisvectors = basisvectors * evecs;
741 : } else {
742 254 : basisvectors = (basisvectors * evecs).pruned(threshold, 1);
743 : }
744 :
745 : // TODO call transformInteraction (see applyRightsideTransformator), perhaps not?
746 : }
747 :
748 : template <class Scalar, class State>
749 0 : void SystemBase<Scalar, State>::canonicalize() {
750 0 : this->buildHamiltonian();
751 :
752 : // Transform the hamiltonian
753 0 : hamiltonian = basisvectors * hamiltonian * basisvectors.adjoint();
754 :
755 : // Transform the basis vectors
756 0 : basisvectors = basisvectors * basisvectors.adjoint();
757 :
758 : // TODO call transformInteraction (see applyRightsideTransformator), perhaps yes?
759 0 : }
760 :
761 : template <class Scalar, class State>
762 2 : void SystemBase<Scalar, State>::unitarize() {
763 2 : this->buildHamiltonian();
764 :
765 : // Determine hash of the list of states
766 : size_t hashvalue_states =
767 2 : utils::hash_range(states.template get<0>().begin(), states.template get<0>().end());
768 :
769 : // Loop over basisvectors
770 2 : size_t num_basisvectors = basisvectors.outerSize();
771 :
772 2 : states.clear();
773 2 : size_t idx_new = 0;
774 30 : for (int k = 0; k < basisvectors.outerSize(); ++k) { // col == idx_vector
775 28 : size_t hashvalue = hashvalue_states;
776 160 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basisvectors, k); triple;
777 132 : ++triple) {
778 132 : utils::hash_combine(hashvalue, triple.row());
779 132 : utils::hash_combine(hashvalue, triple.value());
780 : }
781 :
782 : // Rewrite basisvectors as states
783 28 : std::stringstream ss;
784 28 : ss << std::hex << hashvalue;
785 28 : states.push_back(enumerated_state<State>(idx_new++, createStateFromLabel<State>(ss.str())));
786 : }
787 2 : states.shrink_to_fit();
788 :
789 2 : if (num_basisvectors != states.size()) {
790 0 : throw std::runtime_error("A hash collision occurred.");
791 : }
792 :
793 : // Adapt basisvectors
794 2 : basisvectors.resize(states.size(), states.size());
795 2 : basisvectors.setZero();
796 2 : basisvectors.reserve(states.size());
797 30 : for (size_t idx = 0; idx < states.size(); ++idx) {
798 28 : basisvectors.insert(idx, idx) = 1;
799 : }
800 2 : basisvectors.makeCompressed();
801 :
802 : // Delete caches as they are no longer meaningful
803 2 : basisvectors_unperturbed_cache.resize(0, 0);
804 2 : hamiltonian_unperturbed_cache.resize(0, 0);
805 2 : }
806 :
807 : template <class Scalar, class State>
808 1 : void SystemBase<Scalar, State>::rotate(std::array<double, 3> to_z_axis,
809 : std::array<double, 3> to_y_axis) {
810 1 : auto euler_zyz = this->getEulerAngles(to_z_axis, to_y_axis);
811 1 : this->rotate(euler_zyz[0], euler_zyz[1], euler_zyz[2]);
812 1 : }
813 :
814 : template <class Scalar, class State>
815 10 : void SystemBase<Scalar, State>::rotate(double alpha, double beta,
816 : double gamma) { // TODO applyRightsideTransformator ?
817 : // Build Hamiltonian and basis
818 10 : this->buildHamiltonian();
819 :
820 : // Get the rotator for the basis states
821 20 : Eigen::SparseMatrix<Scalar> transformator = this->buildStaterotator(alpha, beta, gamma);
822 :
823 : // Rotate basis
824 10 : this->transformInteraction(basisvectors.adjoint());
825 :
826 10 : basisvectors = transformator * basisvectors;
827 10 : if (basisvectors_unperturbed_cache.size() != 0) {
828 7 : basisvectors_unperturbed_cache = transformator * basisvectors_unperturbed_cache;
829 : }
830 :
831 10 : this->transformInteraction(basisvectors);
832 10 : }
833 :
834 : template <class Scalar, class State>
835 38 : void SystemBase<Scalar, State>::add(SystemBase<Scalar, State> &system) {
836 : // --- Build bases ---
837 38 : this->buildBasis();
838 38 : system.buildBasis();
839 :
840 38 : size_t size_this = basisvectors.cols();
841 38 : size_t size_other = system.basisvectors.cols();
842 :
843 : // --- Combine system specific variables ---
844 38 : this->incorporate(system);
845 :
846 : // --- Combine universal variables ---
847 38 : if (memory_saving != system.memory_saving) {
848 0 : throw std::runtime_error(
849 : "The value of the variable 'memory_saving' must be the same for both systems.");
850 : }
851 38 : if (is_interaction_already_contained != system.is_interaction_already_contained) {
852 0 : throw std::runtime_error("The value of the variable 'is_interaction_already_contained' "
853 : "must be the same for both systems.");
854 : }
855 38 : if (is_new_hamiltonian_required != system.is_new_hamiltonian_required) {
856 0 : throw std::runtime_error(
857 : "The value of the variable 'is_new_hamiltonian_required' must be the same "
858 : "for both systems.");
859 : }
860 :
861 : // --- Combine states and build transformators ---
862 76 : Eigen::SparseMatrix<Scalar> transformator;
863 : {
864 38 : std::vector<Eigen::Triplet<Scalar>> transformator_triplets;
865 38 : transformator_triplets.reserve(system.states.size());
866 :
867 2115 : for (const auto &entry : system.states) {
868 2077 : auto state_iter = states.template get<1>().find(entry.state);
869 :
870 2077 : size_t newidx = states.size();
871 2077 : if (state_iter == states.template get<1>().end()) {
872 495 : states.push_back(enumerated_state<State>(newidx, entry.state));
873 : } else {
874 1582 : newidx = state_iter->idx;
875 : }
876 :
877 2077 : transformator_triplets.emplace_back(newidx, entry.idx, 1);
878 : }
879 :
880 38 : transformator.resize(states.size(), system.basisvectors.rows());
881 38 : transformator.setFromTriplets(transformator_triplets.begin(), transformator_triplets.end());
882 : }
883 :
884 : // --- Combine basisvectors and basisvectors_unperturbed_cache ---
885 38 : basisvectors.conservativeResize(states.size(), size_this + size_other);
886 38 : basisvectors.rightCols(size_other) = transformator * system.basisvectors;
887 :
888 38 : if ((basisvectors_unperturbed_cache.size() != 0) !=
889 38 : (system.basisvectors_unperturbed_cache.size() != 0)) {
890 0 : throw std::runtime_error("Inconsistent variables at " + std::string(__FILE__) + ":" +
891 0 : std::to_string(__LINE__) + ".");
892 : }
893 :
894 38 : if (basisvectors_unperturbed_cache.size() != 0) {
895 27 : basisvectors_unperturbed_cache.conservativeResize(
896 27 : states.size(),
897 27 : basisvectors_unperturbed_cache.cols() + system.basisvectors_unperturbed_cache.cols());
898 27 : basisvectors_unperturbed_cache.rightCols(system.basisvectors_unperturbed_cache.cols()) =
899 54 : transformator * system.basisvectors_unperturbed_cache;
900 : }
901 :
902 : // --- Check if basis vectors are orthogonal ---
903 76 : Eigen::SparseMatrix<Scalar> tmp =
904 38 : (basisvectors.leftCols(size_this).adjoint() * basisvectors.rightCols(size_other))
905 38 : .pruned(1e-12, 1);
906 38 : if (tmp.nonZeros() != 0) {
907 0 : throw std::runtime_error(
908 : "Two systems cannot be combined if their basis vectors are not orthogonal.");
909 : }
910 :
911 : // --- Combine hamiltonian and hamiltonian_unperturbed_cache ---
912 76 : Eigen::SparseMatrix<Scalar> shifter(hamiltonian.rows() + system.hamiltonian.rows(),
913 : system.hamiltonian.rows());
914 : {
915 38 : std::vector<Eigen::Triplet<Scalar>> shifter_triplets;
916 38 : shifter_triplets.reserve(system.hamiltonian.rows());
917 :
918 960 : for (Eigen::Index idx = 0; idx < system.hamiltonian.rows(); ++idx) {
919 922 : shifter_triplets.emplace_back(hamiltonian.rows() + idx, idx, 1);
920 : }
921 :
922 38 : shifter.setFromTriplets(shifter_triplets.begin(), shifter_triplets.end());
923 : }
924 :
925 38 : hamiltonian.conservativeResize(hamiltonian.rows() + system.hamiltonian.rows(),
926 38 : hamiltonian.cols() + system.hamiltonian.cols());
927 38 : hamiltonian.rightCols(system.hamiltonian.cols()) = shifter * system.hamiltonian;
928 :
929 38 : if ((hamiltonian_unperturbed_cache.size() != 0) !=
930 38 : (system.hamiltonian_unperturbed_cache.size() != 0)) {
931 0 : throw std::runtime_error("Inconsistent variables at " + std::string(__FILE__) + ":" +
932 0 : std::to_string(__LINE__) + ".");
933 : }
934 :
935 38 : if (hamiltonian_unperturbed_cache.size() != 0) {
936 27 : hamiltonian_unperturbed_cache.conservativeResize(
937 27 : hamiltonian_unperturbed_cache.rows() + system.hamiltonian_unperturbed_cache.rows(),
938 27 : hamiltonian_unperturbed_cache.cols() + system.hamiltonian_unperturbed_cache.cols());
939 27 : hamiltonian_unperturbed_cache.rightCols(system.hamiltonian_unperturbed_cache.cols()) =
940 54 : shifter * system.hamiltonian_unperturbed_cache;
941 : }
942 38 : }
943 :
944 : template <class Scalar, class State>
945 4 : void SystemBase<Scalar, State>::constrainBasisvectors(
946 : std::vector<size_t> indices_of_wanted_basisvectors) {
947 4 : this->buildHamiltonian();
948 :
949 : // Check if indices are unique
950 : {
951 8 : std::set<size_t> tmp(indices_of_wanted_basisvectors.begin(),
952 : indices_of_wanted_basisvectors.end());
953 4 : if (tmp.size() < indices_of_wanted_basisvectors.size()) {
954 0 : throw std::runtime_error("Indices are occuring multiple times.");
955 : }
956 : }
957 :
958 : // Build transformator and remove vectors
959 8 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
960 4 : triplets_transformator.reserve(indices_of_wanted_basisvectors.size());
961 :
962 4 : size_t idx_new = 0;
963 20 : for (const auto &idx : indices_of_wanted_basisvectors) {
964 16 : if (Eigen::Index(idx) >= basisvectors.cols()) {
965 0 : throw std::runtime_error("A basis vector with index " + std::to_string(idx) +
966 : " could not be found.");
967 : }
968 16 : triplets_transformator.emplace_back(idx, idx_new++, 1);
969 : }
970 :
971 4 : this->applyRightsideTransformator(triplets_transformator);
972 4 : }
973 :
974 : template <class Scalar, class State>
975 4 : void SystemBase<Scalar, State>::applySchriefferWolffTransformation(
976 : SystemBase<Scalar, State> &system0) {
977 4 : this->diagonalize();
978 4 : system0.buildHamiltonian();
979 :
980 : // Check that system, on which applySchriefferWolffTransformation() is called, is unitary
981 4 : if (!this->checkIsUnitary(basisvectors)) {
982 0 : throw std::runtime_error("The system, on which applySchriefferWolffTransformation() is "
983 : "called, is not unitary. Call unitarize() on the system.");
984 : }
985 :
986 : // Check that system0, i.e. the unperturbed system, is unitary
987 4 : if (!this->checkIsUnitary(system0.basisvectors)) {
988 0 : throw std::runtime_error(
989 : "The unperturbed system passed to applySchriefferWolffTransformation() is not "
990 : "unitary. Call unitarize() on the system.");
991 : }
992 :
993 : // Check that the unperturbed system is diagonal
994 4 : if (!this->checkIsDiagonal(system0.hamiltonian)) {
995 0 : throw std::runtime_error("The unperturbed system passed to "
996 : "applySchriefferWolffTransformation() is not diagonal.");
997 : }
998 :
999 : // --- Express the basis vectors of system0 as a linearcombination of all states ---
1000 :
1001 : // Combine states and build transformators
1002 8 : Eigen::SparseMatrix<Scalar> transformator;
1003 : {
1004 4 : std::vector<Eigen::Triplet<Scalar>> transformator_triplets;
1005 4 : transformator_triplets.reserve(system0.states.size());
1006 :
1007 200 : for (const auto &entry : system0.states) {
1008 196 : auto state_iter = states.template get<1>().find(entry.state);
1009 :
1010 : // Check that all the states of system0 occur in the perturbed system
1011 196 : if (state_iter == states.template get<1>().end()) {
1012 0 : throw std::runtime_error(
1013 : "The unperturbed system contains states that are not occuring. "
1014 : "This might happen if you did not use the same system as basis for the "
1015 : "unperturbed and perturbed system.");
1016 : }
1017 196 : size_t newidx = state_iter->idx;
1018 :
1019 196 : transformator_triplets.emplace_back(newidx, entry.idx, 1);
1020 : }
1021 :
1022 4 : transformator.resize(states.size(), system0.states.size());
1023 4 : transformator.setFromTriplets(transformator_triplets.begin(), transformator_triplets.end());
1024 : }
1025 :
1026 : // Get basisvectors of system0
1027 8 : Eigen::SparseMatrix<Scalar> low_energy_basis0 = transformator * system0.basisvectors;
1028 :
1029 : // --- Find basisvectors which corresponds to the basis vectors of system0 ---
1030 :
1031 : // Get overlaps between basis vectors
1032 8 : Eigen::VectorX<double> overlap = (basisvectors.adjoint() * low_energy_basis0).cwiseAbs2() *
1033 4 : Eigen::VectorX<double>::Ones(low_energy_basis0.cols());
1034 :
1035 : // Get indices of the low_energy_basis0.cols() largest entries and build transformator
1036 : {
1037 8 : std::vector<int> indices(basisvectors.cols());
1038 4 : std::iota(indices.begin(), indices.end(), 0);
1039 4 : std::nth_element(indices.begin(), indices.begin() + low_energy_basis0.cols(), indices.end(),
1040 368 : [&overlap](int a, int b) { return overlap[a] > overlap[b]; });
1041 :
1042 4 : std::vector<Eigen::Triplet<Scalar>> transformator_triplets;
1043 4 : transformator_triplets.reserve(low_energy_basis0.cols());
1044 :
1045 20 : for (int i = 0; i < low_energy_basis0.cols(); ++i) {
1046 16 : transformator_triplets.emplace_back(indices[i], i, 1);
1047 : }
1048 :
1049 4 : transformator.resize(basisvectors.cols(), low_energy_basis0.cols());
1050 4 : transformator.setFromTriplets(transformator_triplets.begin(), transformator_triplets.end());
1051 : }
1052 :
1053 : // Get corresponding basis vectors
1054 8 : Eigen::SparseMatrix<Scalar> low_energy_basis = basisvectors * transformator;
1055 :
1056 : // --- Perform the Schrieffer Wolff transformation ---
1057 :
1058 : // Projector on the selected basis vectors (typically low-energy subspace)
1059 8 : Eigen::SparseMatrix<Scalar> projector0 = low_energy_basis0 * low_energy_basis0.adjoint();
1060 8 : Eigen::SparseMatrix<Scalar> projector = low_energy_basis * low_energy_basis.adjoint();
1061 :
1062 : // Reflection operator which change signes of the selectes basis vectors but act
1063 : // trivially on others
1064 8 : Eigen::MatrixX<Scalar> reflection0 =
1065 4 : Eigen::MatrixX<Scalar>::Identity(states.size(), states.size()) - 2 * projector0;
1066 8 : Eigen::MatrixX<Scalar> reflection =
1067 4 : Eigen::MatrixX<Scalar>::Identity(states.size(), states.size()) - 2 * projector;
1068 :
1069 : // Direct rotator from low_energy_basis to low_energy_basis0
1070 8 : Eigen::SparseMatrix<Scalar> rotator = (reflection0 * reflection).sqrt().sparseView();
1071 :
1072 4 : if (std::isnan(std::abs(rotator.coeffRef(0, 0)))) {
1073 0 : throw std::runtime_error(
1074 : "Error in calculating the matrix square root. Try to use the picomplex module or "
1075 : "increase accuracy by choosing smaller thresholds.");
1076 : }
1077 :
1078 : // Calculate the effective Hamiltonian
1079 4 : transformator = basisvectors.adjoint() * rotator.adjoint() * projector0 * low_energy_basis0;
1080 :
1081 4 : this->applyRightsideTransformator(transformator);
1082 4 : }
1083 :
1084 : ////////////////////////////////////////////////////////////////////
1085 : /// Methods to manipulate entries of the Hamiltonian ///////////////
1086 : ////////////////////////////////////////////////////////////////////
1087 :
1088 : template <class Scalar, class State>
1089 7 : size_t SystemBase<Scalar, State>::getStateIndex(const State &searched_state) {
1090 7 : this->buildBasis();
1091 :
1092 7 : if (utils::is_true(searched_state.isGeneralized())) {
1093 0 : throw std::runtime_error("The method must not be called on a generalized state.");
1094 : }
1095 :
1096 7 : auto state_iter = states.template get<1>().find(searched_state);
1097 7 : if (state_iter == states.template get<1>().end()) {
1098 0 : throw std::runtime_error("The method is called on a non-existing state.");
1099 : }
1100 :
1101 7 : return state_iter->idx;
1102 : }
1103 :
1104 : template <class Scalar, class State>
1105 : std::vector<size_t>
1106 7 : SystemBase<Scalar, State>::getStateIndex(const std::vector<State> &searched_states) {
1107 7 : this->buildBasis();
1108 :
1109 7 : std::vector<size_t> indices;
1110 7 : indices.reserve(searched_states.size());
1111 :
1112 35 : for (const auto &state : searched_states) {
1113 28 : if (utils::is_true(state.isGeneralized())) {
1114 0 : throw std::runtime_error("The method must not be called on a generalized state.");
1115 : }
1116 :
1117 28 : auto state_iter = states.template get<1>().find(state);
1118 28 : if (state_iter == states.template get<1>().end()) {
1119 0 : throw std::runtime_error("The method is called on a non-existing state.");
1120 : }
1121 28 : indices.push_back(state_iter->idx);
1122 : }
1123 :
1124 7 : return indices;
1125 : }
1126 :
1127 : template <class Scalar, class State>
1128 7 : size_t SystemBase<Scalar, State>::getBasisvectorIndex(const State &searched_state) {
1129 7 : this->buildBasis();
1130 :
1131 7 : Eigen::Index stateidx = this->getStateIndex(searched_state);
1132 :
1133 7 : double maxval = -1;
1134 7 : size_t col_with_maxval = -1;
1135 365 : for (int k = 0; k < basisvectors.outerSize(); ++k) { // col == idx_vector
1136 30414 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basisvectors, k); triple;
1137 30056 : ++triple) {
1138 30298 : if (triple.row() == stateidx) {
1139 242 : if (std::abs(triple.value()) > maxval) {
1140 23 : col_with_maxval = triple.col();
1141 23 : maxval = std::abs(triple.value());
1142 : }
1143 242 : break;
1144 : }
1145 : }
1146 : }
1147 :
1148 7 : return col_with_maxval;
1149 : }
1150 :
1151 : template <class Scalar, class State>
1152 : std::vector<size_t>
1153 5 : SystemBase<Scalar, State>::getBasisvectorIndex(const std::vector<State> &searched_states) {
1154 5 : this->buildBasis();
1155 :
1156 : // Ensure that the states within searched_states are unique
1157 : {
1158 10 : std::set<State> tmp(searched_states.begin(), searched_states.end());
1159 5 : if (tmp.size() < searched_states.size()) {
1160 0 : throw std::runtime_error("States are occuring multiple times.");
1161 : }
1162 : }
1163 :
1164 : // Construct canonical basis vectors
1165 10 : Eigen::SparseMatrix<Scalar> canonicalbasis;
1166 : {
1167 10 : std::vector<size_t> state_indices = this->getStateIndex(searched_states);
1168 5 : std::vector<Eigen::Triplet<Scalar>> canonicalbasis_triplets;
1169 5 : canonicalbasis_triplets.reserve(searched_states.size());
1170 25 : for (size_t i = 0; i < state_indices.size(); ++i) {
1171 20 : canonicalbasis_triplets.emplace_back(state_indices[i], i, 1);
1172 : }
1173 5 : canonicalbasis.resize(states.size(), searched_states.size());
1174 5 : canonicalbasis.setFromTriplets(canonicalbasis_triplets.begin(),
1175 10 : canonicalbasis_triplets.end());
1176 : }
1177 :
1178 : // Get overlaps between basis vectors
1179 10 : Eigen::SparseMatrix<double> overlap = (canonicalbasis.adjoint() * basisvectors).cwiseAbs2();
1180 10 : Eigen::VectorX<double> overlap_total =
1181 5 : Eigen::VectorX<double>::Ones(canonicalbasis.cols()).transpose() * overlap;
1182 :
1183 : // Get indices of the canonicalbasis.cols() largest entries
1184 10 : std::vector<size_t> indices_available(basisvectors.cols());
1185 5 : std::iota(indices_available.begin(), indices_available.end(), 0);
1186 5 : std::nth_element(indices_available.begin(), indices_available.begin() + canonicalbasis.cols(),
1187 1260 : indices_available.end(), [&overlap_total](size_t a, size_t b) {
1188 630 : return overlap_total[a] > overlap_total[b];
1189 : });
1190 5 : indices_available.resize(canonicalbasis.cols());
1191 :
1192 : // Resort the indices to match the order of searched_states, TODO use Hungarian algorithm
1193 5 : std::vector<size_t> indices(canonicalbasis.cols(), std::numeric_limits<size_t>::max());
1194 25 : for (const auto &k : indices_available) {
1195 : size_t row_with_maxval;
1196 20 : double maxval = -1;
1197 56 : for (typename Eigen::SparseMatrix<double>::InnerIterator triple(overlap, k); triple;
1198 36 : ++triple) {
1199 64 : if (indices[triple.row()] == std::numeric_limits<size_t>::max() &&
1200 28 : triple.value() > maxval) {
1201 22 : row_with_maxval = triple.row();
1202 22 : maxval = triple.value();
1203 : }
1204 : }
1205 :
1206 20 : if (maxval == -1) {
1207 0 : throw std::runtime_error("There is a state for which no basis vector could be found.");
1208 : }
1209 :
1210 20 : indices[row_with_maxval] = k;
1211 : }
1212 :
1213 10 : return indices;
1214 : }
1215 :
1216 : template <class Scalar, class State>
1217 0 : void SystemBase<Scalar, State>::forgetStatemixing() {
1218 0 : this->diagonalize();
1219 :
1220 0 : std::vector<Eigen::Triplet<Scalar>> basisvectors_triplets;
1221 0 : basisvectors_triplets.reserve(basisvectors.cols());
1222 0 : double threshold = std::sqrt(0.5);
1223 :
1224 0 : for (int k = 0; k < basisvectors.outerSize(); ++k) { // col == idx_vector
1225 0 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basisvectors, k); triple;
1226 0 : ++triple) {
1227 0 : if (std::abs(triple.value()) > threshold) {
1228 0 : basisvectors_triplets.emplace_back(triple.row(), triple.col(), 1);
1229 0 : break;
1230 : }
1231 : }
1232 : }
1233 :
1234 0 : if (Eigen::Index(basisvectors_triplets.size()) < basisvectors.cols()) {
1235 0 : throw std::runtime_error(
1236 : "The states are mixed too strongly for calling forgetStatemixing().");
1237 : }
1238 :
1239 0 : basisvectors.setFromTriplets(basisvectors_triplets.begin(), basisvectors_triplets.end());
1240 0 : }
1241 :
1242 : template <class Scalar, class State>
1243 0 : Scalar SystemBase<Scalar, State>::getHamiltonianEntry(const State &state_row,
1244 : const State &state_col) {
1245 0 : this->buildHamiltonian();
1246 :
1247 0 : size_t idx_row = this->getStateIndex(state_row);
1248 0 : size_t idx_col = this->getStateIndex(state_col);
1249 :
1250 0 : return (basisvectors.row(idx_row) * hamiltonian).dot(basisvectors.row(idx_col));
1251 : }
1252 :
1253 : template <class Scalar, class State>
1254 0 : void SystemBase<Scalar, State>::setHamiltonianEntry(const State &state_row, const State &state_col,
1255 : Scalar value) {
1256 0 : this->buildHamiltonian();
1257 :
1258 0 : size_t idx_row = this->getStateIndex(state_row);
1259 0 : size_t idx_col = this->getStateIndex(state_col);
1260 :
1261 0 : value -= (basisvectors.row(idx_row) * hamiltonian).dot(basisvectors.row(idx_col));
1262 :
1263 0 : Eigen::SparseMatrix<Scalar> tmp(states.size(), states.size());
1264 0 : tmp.reserve(2);
1265 0 : tmp.insert(idx_row, idx_col) = value;
1266 0 : if (idx_row != idx_col) {
1267 0 : tmp.insert(idx_col, idx_row) = utils::conjugate(value);
1268 : }
1269 0 : tmp.makeCompressed();
1270 :
1271 0 : hamiltonian += basisvectors.adjoint() * tmp * basisvectors;
1272 0 : }
1273 :
1274 : template <class Scalar, class State>
1275 0 : void SystemBase<Scalar, State>::addHamiltonianEntry(const State &state_row, const State &state_col,
1276 : Scalar value) {
1277 0 : this->buildHamiltonian();
1278 :
1279 0 : size_t idx_row = this->getStateIndex(state_row);
1280 0 : size_t idx_col = this->getStateIndex(state_col);
1281 :
1282 0 : Eigen::SparseMatrix<Scalar> tmp(states.size(), states.size());
1283 0 : tmp.reserve(2);
1284 0 : tmp.insert(idx_row, idx_col) = value;
1285 0 : if (idx_row != idx_col) {
1286 0 : tmp.insert(idx_col, idx_row) = utils::conjugate(value);
1287 : }
1288 0 : tmp.makeCompressed();
1289 :
1290 0 : hamiltonian += basisvectors.adjoint() * tmp * basisvectors;
1291 0 : }
1292 :
1293 : template <class Scalar, class State>
1294 0 : SystemBase<Scalar, State>::SystemBase()
1295 0 : : energy_min(std::numeric_limits<double>::lowest()),
1296 0 : energy_max(std::numeric_limits<double>::max()) {}
1297 :
1298 : template <class Scalar, class State>
1299 121 : SystemBase<Scalar, State>::SystemBase(MatrixElementCache &cache)
1300 242 : : m_cache(std::addressof(cache)), energy_min(std::numeric_limits<double>::lowest()),
1301 121 : energy_max(std::numeric_limits<double>::max()) {}
1302 :
1303 : template <class Scalar, class State>
1304 0 : SystemBase<Scalar, State>::SystemBase(MatrixElementCache &cache, bool memory_saving)
1305 0 : : m_cache(std::addressof(cache)), energy_min(std::numeric_limits<double>::lowest()),
1306 0 : energy_max(std::numeric_limits<double>::max()), memory_saving(memory_saving) {}
1307 :
1308 : ////////////////////////////////////////////////////////////////////
1309 : /// Helper method that shoul be called by the derived classes //////
1310 : ////////////////////////////////////////////////////////////////////
1311 :
1312 : template <class Scalar, class State>
1313 435 : void SystemBase<Scalar, State>::onParameterChange() {
1314 : // Check variables for consistency
1315 435 : if ((basisvectors_unperturbed_cache.size() == 0) !=
1316 435 : (hamiltonian_unperturbed_cache.size() == 0)) {
1317 0 : throw std::runtime_error("Inconsistent variables at " + std::string(__FILE__) + ":" +
1318 0 : std::to_string(__LINE__) + ".");
1319 : }
1320 :
1321 : // Throw error if the Hamiltonian cannot be changed anymore
1322 435 : if (is_interaction_already_contained && basisvectors_unperturbed_cache.size() == 0) {
1323 0 : throw std::runtime_error(
1324 : "If memory saving is activated or unitarize() has been called, one cannot change "
1325 : "parameters after interaction was added to the Hamiltonian.");
1326 : }
1327 :
1328 435 : is_new_hamiltonian_required = true;
1329 435 : }
1330 :
1331 : template <class Scalar, class State>
1332 83 : void SystemBase<Scalar, State>::onSymmetryChange() {
1333 : // Throw error if the Symmetry cannot be changed anymore
1334 83 : if (!states.empty()) {
1335 0 : throw std::runtime_error("One cannot change symmetries after the basis was built.");
1336 : }
1337 83 : }
1338 :
1339 : ////////////////////////////////////////////////////////////////////
1340 : /// Helper methods to check diagonality and unitarity of a matrix //
1341 : ////////////////////////////////////////////////////////////////////
1342 :
1343 : template <class Scalar, class State>
1344 483 : bool SystemBase<Scalar, State>::checkIsDiagonal(const Eigen::SparseMatrix<Scalar> &mat) {
1345 966 : Eigen::SparseMatrix<Scalar> tmp = mat;
1346 483 : tmp.prune(1e-12, 1);
1347 :
1348 6832 : for (int k = 0; k < tmp.outerSize(); ++k) {
1349 13368 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(tmp, k); triple; ++triple) {
1350 7019 : if (triple.row() != triple.col()) {
1351 335 : return false;
1352 : }
1353 : }
1354 : }
1355 148 : return true;
1356 : }
1357 :
1358 : template <class Scalar, class State>
1359 8 : bool SystemBase<Scalar, State>::checkIsUnitary(const Eigen::SparseMatrix<Scalar> &mat) {
1360 16 : Eigen::SparseMatrix<Scalar> tmp = (mat.adjoint() * mat).pruned(1e-12, 1);
1361 :
1362 8 : if (tmp.nonZeros() != tmp.outerSize()) {
1363 0 : return false;
1364 : }
1365 :
1366 240 : for (int k = 0; k < tmp.outerSize(); ++k) {
1367 464 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(tmp, k); triple; ++triple) {
1368 232 : if (triple.row() != triple.col()) {
1369 0 : return false;
1370 : }
1371 464 : if (std::abs(std::real(triple.value()) - 1) > 1e-12 ||
1372 232 : std::abs(std::imag(triple.value())) > 1e-12) {
1373 0 : return false;
1374 : }
1375 : }
1376 : }
1377 8 : return true;
1378 : }
1379 :
1380 : ////////////////////////////////////////////////////////////////////
1381 : /// Helper methods to check the validity of states /////////////////
1382 : ////////////////////////////////////////////////////////////////////
1383 :
1384 : template <class Scalar, class State>
1385 0 : bool SystemBase<Scalar, State>::checkIsQuantumstateValid(const State &state) {
1386 0 : return checkIsQuantumstateValid(state, state.isArtificial());
1387 : }
1388 :
1389 : template <class Scalar, class State>
1390 0 : bool SystemBase<Scalar, State>::checkIsQuantumstateValid(const StateOne &state, bool a) {
1391 0 : return a ||
1392 0 : (checkIsQuantumnumberValid(state.getN(), range_n) &&
1393 0 : checkIsQuantumnumberValid(state.getL(), range_l) &&
1394 0 : checkIsQuantumnumberValid(state.getJ(), range_j) &&
1395 0 : checkIsQuantumnumberValid(state.getM(), range_m));
1396 : }
1397 :
1398 : template <class Scalar, class State>
1399 0 : bool SystemBase<Scalar, State>::checkIsQuantumstateValid(const StateTwo &state,
1400 : std::array<bool, 2> a) {
1401 0 : bool valid = true;
1402 0 : for (int idx = 0; idx < 2; ++idx) {
1403 0 : valid = valid &&
1404 0 : (a[idx] ||
1405 0 : (checkIsQuantumnumberValid(state.getN(idx), range_n) &&
1406 0 : checkIsQuantumnumberValid(state.getL(idx), range_l) &&
1407 0 : checkIsQuantumnumberValid(state.getJ(idx), range_j) &&
1408 0 : checkIsQuantumnumberValid(state.getM(idx), range_m)));
1409 : }
1410 0 : return valid;
1411 : }
1412 :
1413 : template <class Scalar, class State>
1414 155696 : bool SystemBase<Scalar, State>::checkIsEnergyValid(double e) {
1415 265643 : return (e > energy_min || energy_min == std::numeric_limits<double_t>::lowest()) &&
1416 207247 : (e < energy_max ||
1417 97300 : energy_max ==
1418 252996 : std::numeric_limits<double_t>::max()); // TODO take into account numerical errors
1419 : }
1420 :
1421 : ////////////////////////////////////////////////////////////////////
1422 : /// Helper methods to change the set of basis vectors //////////////
1423 : ////////////////////////////////////////////////////////////////////
1424 :
1425 : template <class Scalar, class State>
1426 75 : void SystemBase<Scalar, State>::applyLeftsideTransformator(
1427 : std::vector<Eigen::Triplet<Scalar>> &triplets_transformator) {
1428 150 : Eigen::SparseMatrix<Scalar> transformator(triplets_transformator.size(), basisvectors.rows());
1429 75 : transformator.setFromTriplets(triplets_transformator.begin(), triplets_transformator.end());
1430 :
1431 75 : this->applyLeftsideTransformator(transformator);
1432 75 : }
1433 :
1434 : template <class Scalar, class State>
1435 75 : void SystemBase<Scalar, State>::applyLeftsideTransformator(
1436 : Eigen::SparseMatrix<Scalar> &transformator) {
1437 : // Apply transformator in order to remove rows from the basisvector matrix (i.e. states)
1438 75 : basisvectors = transformator * basisvectors;
1439 75 : if (basisvectors_unperturbed_cache.size() != 0) {
1440 0 : basisvectors_unperturbed_cache = transformator * basisvectors_unperturbed_cache;
1441 : }
1442 75 : }
1443 :
1444 : template <class Scalar, class State>
1445 79 : void SystemBase<Scalar, State>::applyRightsideTransformator(
1446 : std::vector<Eigen::Triplet<Scalar>> &triplets_transformator) {
1447 158 : Eigen::SparseMatrix<Scalar> transformator(basisvectors.cols(), triplets_transformator.size());
1448 79 : transformator.setFromTriplets(triplets_transformator.begin(), triplets_transformator.end());
1449 :
1450 79 : this->applyRightsideTransformator(transformator);
1451 79 : }
1452 :
1453 : template <class Scalar, class State>
1454 83 : void SystemBase<Scalar, State>::applyRightsideTransformator(
1455 : Eigen::SparseMatrix<Scalar> &transformator) {
1456 : // Apply transformator in order to remove columns from the basisvector matrix (i.e. basis
1457 : // vectors)
1458 83 : basisvectors = basisvectors * transformator;
1459 83 : if (basisvectors_unperturbed_cache.size() != 0) {
1460 3 : basisvectors_unperturbed_cache = basisvectors_unperturbed_cache * transformator;
1461 : }
1462 :
1463 : // Apply transformator in order to remove rows and columns from the matrices that help
1464 : // constructing the total Hamiltonian
1465 83 : this->transformInteraction(transformator);
1466 :
1467 : // Apply transformator in order to remove rows and columns from the total Hamiltonian
1468 83 : hamiltonian = transformator.adjoint() * hamiltonian * transformator;
1469 83 : if (hamiltonian_unperturbed_cache.size() != 0) {
1470 3 : hamiltonian_unperturbed_cache =
1471 6 : transformator.adjoint() * hamiltonian_unperturbed_cache * transformator;
1472 : }
1473 83 : }
1474 :
1475 : template <class Scalar, class State>
1476 75 : void SystemBase<Scalar, State>::removeRestrictedStates(
1477 : std::function<bool(const enumerated_state<State> &)> checkIsValidEntry) {
1478 : // Build transformator and remove states
1479 150 : typename states_set<State>::type states_new;
1480 75 : states_new.reserve(states.size());
1481 150 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
1482 75 : triplets_transformator.reserve(states.size());
1483 :
1484 75 : size_t idx_new = 0;
1485 :
1486 29999 : for (const auto &entry : states) {
1487 14962 : if (checkIsValidEntry(entry)) {
1488 8284 : states_new.push_back(enumerated_state<State>(idx_new, entry.state));
1489 8284 : triplets_transformator.emplace_back(idx_new, entry.idx, 1);
1490 8284 : ++idx_new;
1491 : }
1492 : }
1493 :
1494 75 : states_new.shrink_to_fit();
1495 75 : states = states_new;
1496 :
1497 75 : this->applyLeftsideTransformator(triplets_transformator);
1498 75 : }
1499 :
1500 : ////////////////////////////////////////////////////////////////////
1501 : /// Utility methods ////////////////////////////////////////////////
1502 : ////////////////////////////////////////////////////////////////////
1503 :
1504 : template <class Scalar, class State>
1505 : Eigen::Matrix<double, 3, 3>
1506 1 : SystemBase<Scalar, State>::buildRotator(std::array<double, 3> to_z_axis,
1507 : std::array<double, 3> to_y_axis) {
1508 1 : auto to_z_axis_mapped = Eigen::Map<Eigen::Matrix<double, 3, 1>>(&to_z_axis[0]).normalized();
1509 1 : auto to_y_axis_mapped = Eigen::Map<Eigen::Matrix<double, 3, 1>>(&to_y_axis[0]).normalized();
1510 :
1511 1 : double tolerance = 1e-16;
1512 1 : if (std::abs(to_z_axis_mapped.dot(to_y_axis_mapped)) > tolerance) {
1513 0 : throw std::runtime_error("The z-axis and the y-axis are not orhogonal.");
1514 : }
1515 :
1516 1 : Eigen::Matrix<double, 3, 3> transformator;
1517 1 : transformator << to_y_axis_mapped.cross(to_z_axis_mapped), to_y_axis_mapped, to_z_axis_mapped;
1518 :
1519 2 : return transformator;
1520 : }
1521 :
1522 : template <class Scalar, class State>
1523 3 : Eigen::Matrix<double, 3, 3> SystemBase<Scalar, State>::buildRotator(const double &alpha,
1524 : const double &beta,
1525 : const double &gamma) {
1526 3 : Eigen::Matrix<double, 3, 3> transformator;
1527 :
1528 6 : transformator = Eigen::AngleAxisd(alpha, Eigen::Matrix<double, 3, 1>::UnitZ()) *
1529 6 : Eigen::AngleAxisd(beta, Eigen::Matrix<double, 3, 1>::UnitY()) *
1530 3 : Eigen::AngleAxisd(gamma, Eigen::Matrix<double, 3, 1>::UnitZ());
1531 :
1532 3 : return transformator;
1533 : }
1534 :
1535 : template <class Scalar, class State>
1536 : Eigen::Matrix<double, 3, 1>
1537 1 : SystemBase<Scalar, State>::getEulerAngles(const std::array<double, 3> &to_z_axis,
1538 : const std::array<double, 3> &to_y_axis) {
1539 1 : Eigen::Matrix<double, 3, 3> rotator = this->buildRotator(to_z_axis, to_y_axis);
1540 1 : Eigen::Matrix<double, 3, 1> euler_zyz = rotator.eulerAngles(2, 1, 2);
1541 2 : return euler_zyz; // alpha, beta, gamma
1542 : }
1543 :
1544 : template <class Scalar, class State>
1545 164 : void SystemBase<Scalar, State>::forgetRestrictions() {
1546 164 : energy_min = std::numeric_limits<double>::lowest();
1547 164 : energy_max = std::numeric_limits<double>::max();
1548 164 : range_n.clear();
1549 164 : range_l.clear();
1550 164 : range_j.clear();
1551 164 : range_m.clear();
1552 164 : states_to_add.clear();
1553 164 : }
1554 :
1555 : ////////////////////////////////////////////////////////////////////
1556 : /// Method to update the system ////////////////////////////////////
1557 : ////////////////////////////////////////////////////////////////////
1558 :
1559 : template <class Scalar, class State>
1560 0 : void SystemBase<Scalar, State>::updateEverything() {
1561 :
1562 0 : if (!range_n.empty() || !range_l.empty() || !range_j.empty() || !range_m.empty()) {
1563 :
1564 : ////////////////////////////////////////////////////////////////////
1565 : /// Remove restricted states ///////////////////////////////////////
1566 : ////////////////////////////////////////////////////////////////////
1567 :
1568 : // Remove states if the quantum numbers are not allowed
1569 0 : removeRestrictedStates([=](const enumerated_state<State> &entry) -> bool {
1570 0 : return this->checkIsQuantumstateValid(entry.state);
1571 : });
1572 :
1573 : // Comunicate that the list of states has changed
1574 0 : this->onStatesChange();
1575 : }
1576 :
1577 0 : if (!range_n.empty() || !range_l.empty() || !range_j.empty() || !range_m.empty() ||
1578 0 : energy_min != std::numeric_limits<double>::lowest() ||
1579 0 : energy_max != std::numeric_limits<double>::max()) {
1580 :
1581 : ////////////////////////////////////////////////////////////////////
1582 : /// Remove vectors with restricted energies or too small norm //////
1583 : ////////////////////////////////////////////////////////////////////
1584 :
1585 : // Build transformator and remove vectors (if their energy is not allowed or the squared
1586 : // norm too small)
1587 0 : std::vector<Eigen::Triplet<Scalar>> triplets_transformator;
1588 0 : triplets_transformator.reserve(basisvectors.cols());
1589 :
1590 0 : size_t idx_new = 0;
1591 0 : for (int idx = 0; idx < basisvectors.cols(); ++idx) { // idx = col = num basis vector
1592 :
1593 0 : if (checkIsEnergyValid(std::real(hamiltonian.coeff(idx, idx)))) {
1594 0 : double_t sqnorm = 0;
1595 :
1596 : // Calculate the square norm of the columns of the basisvector matrix
1597 0 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basisvectors, idx);
1598 0 : triple; ++triple) {
1599 0 : sqnorm += std::pow(std::abs(triple.value()), 2);
1600 : }
1601 0 : if (sqnorm > threshold_for_sqnorm) {
1602 0 : triplets_transformator.emplace_back(idx, idx_new++, 1);
1603 : }
1604 : }
1605 : }
1606 :
1607 0 : this->applyRightsideTransformator(triplets_transformator);
1608 :
1609 : ////////////////////////////////////////////////////////////////////
1610 : /// Remove states that barely occur within the vectors /////////////
1611 : ////////////////////////////////////////////////////////////////////
1612 :
1613 : // Calculate the square norm of the rows of the basisvector matrix
1614 0 : std::vector<double> sqnorm_list(basisvectors.rows(), 0);
1615 0 : for (int k = 0; k < this->basisvectors.cols(); ++k) {
1616 0 : for (typename Eigen::SparseMatrix<Scalar>::InnerIterator triple(basisvectors, k);
1617 0 : triple; ++triple) {
1618 0 : sqnorm_list[triple.row()] += std::pow(std::abs(triple.value()), 2);
1619 : }
1620 : }
1621 :
1622 : // Remove states if the squared norm is to small
1623 0 : removeRestrictedStates([=](const enumerated_state<State> &entry) -> bool {
1624 0 : return sqnorm_list[entry.idx] > threshold_for_sqnorm;
1625 : });
1626 :
1627 : // Comunicate that the list of states has changed
1628 0 : this->onStatesChange();
1629 : }
1630 :
1631 0 : if (!states_to_add.empty()) {
1632 0 : throw std::runtime_error(
1633 : "States cannot be added anymore once the basis vectors have been created.");
1634 : }
1635 :
1636 : // Forget basis restrictions as they were applied now
1637 0 : this->forgetRestrictions();
1638 0 : }
1639 :
1640 : template class SystemBase<std::complex<double>, StateOne>;
1641 : template class SystemBase<std::complex<double>, StateTwo>;
1642 : template class SystemBase<double, StateOne>;
1643 : template class SystemBase<double, StateTwo>;
|
