Line data Source code
1 : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/basis/Basis.hpp"
5 :
6 : #include "pairinteraction/basis/BasisAtom.hpp"
7 : #include "pairinteraction/basis/BasisPair.hpp"
8 : #include "pairinteraction/enums/Parity.hpp"
9 : #include "pairinteraction/enums/TransformationType.hpp"
10 : #include "pairinteraction/ket/KetAtom.hpp"
11 : #include "pairinteraction/ket/KetPair.hpp"
12 : #include "pairinteraction/utils/eigen_assertion.hpp"
13 : #include "pairinteraction/utils/eigen_compat.hpp"
14 : #include "pairinteraction/utils/wigner.hpp"
15 :
16 : #include <cassert>
17 : #include <numeric>
18 : #include <set>
19 :
20 : namespace pairinteraction {
21 :
22 : template <typename Scalar>
23 : class BasisAtom;
24 :
25 : template <typename Derived>
26 130 : void Basis<Derived>::perform_sorter_checks(const std::vector<TransformationType> &labels) const {
27 : // Check if the labels are valid sorting labels
28 236 : for (const auto &label : labels) {
29 106 : if (!utils::is_sorting(label)) {
30 0 : throw std::invalid_argument("One of the labels is not a valid sorting label.");
31 : }
32 : }
33 130 : }
34 :
35 : template <typename Derived>
36 57 : void Basis<Derived>::perform_blocks_checks(
37 : const std::set<TransformationType> &unique_labels) const {
38 : // Check if the states are sorted by the requested labels
39 57 : std::set<TransformationType> unique_labels_present;
40 89 : for (const auto &label : get_transformation().transformation_type) {
41 58 : if (!utils::is_sorting(label) || unique_labels_present.size() >= unique_labels.size()) {
42 26 : break;
43 : }
44 32 : unique_labels_present.insert(label);
45 : }
46 57 : if (unique_labels != unique_labels_present) {
47 0 : throw std::invalid_argument("The states are not sorted by the requested labels.");
48 : }
49 :
50 : // Throw a meaningful error if getting the blocks by energy is requested as this might be a
51 : // common mistake
52 57 : if (unique_labels.count(TransformationType::SORT_BY_ENERGY) > 0) {
53 0 : throw std::invalid_argument("States do not store the energy and thus no energy blocks can "
54 : "be obtained. Use an energy operator instead.");
55 : }
56 57 : }
57 :
58 : template <typename Derived>
59 40 : Basis<Derived>::Basis(ketvec_t &&kets)
60 80 : : kets(std::move(kets)), coefficients{{static_cast<Eigen::Index>(this->kets.size()),
61 40 : static_cast<Eigen::Index>(this->kets.size())},
62 80 : {TransformationType::SORT_BY_KET}} {
63 40 : if (this->kets.empty()) {
64 0 : throw std::invalid_argument("The basis must contain at least one element.");
65 : }
66 40 : state_index_to_quantum_number_f.reserve(this->kets.size());
67 40 : state_index_to_quantum_number_m.reserve(this->kets.size());
68 40 : state_index_to_parity.reserve(this->kets.size());
69 40 : ket_to_ket_index.reserve(this->kets.size());
70 40 : size_t index = 0;
71 9721 : for (const auto &ket : this->kets) {
72 9681 : state_index_to_quantum_number_f.push_back(ket->get_quantum_number_f());
73 9681 : state_index_to_quantum_number_m.push_back(ket->get_quantum_number_m());
74 9681 : state_index_to_parity.push_back(ket->get_parity());
75 9681 : ket_to_ket_index[ket] = index++;
76 9681 : if (ket->get_quantum_number_f() == std::numeric_limits<real_t>::max()) {
77 8494 : _has_quantum_number_f = false;
78 : }
79 9681 : if (ket->get_quantum_number_m() == std::numeric_limits<real_t>::max()) {
80 0 : _has_quantum_number_m = false;
81 : }
82 9681 : if (ket->get_parity() == Parity::UNKNOWN) {
83 8494 : _has_parity = false;
84 : }
85 : }
86 40 : state_index_to_ket_index.resize(this->kets.size());
87 40 : std::iota(state_index_to_ket_index.begin(), state_index_to_ket_index.end(), 0);
88 40 : ket_index_to_state_index.resize(this->kets.size());
89 40 : std::iota(ket_index_to_state_index.begin(), ket_index_to_state_index.end(), 0);
90 40 : coefficients.matrix.setIdentity();
91 40 : }
92 :
93 : template <typename Derived>
94 9 : bool Basis<Derived>::has_quantum_number_f() const {
95 9 : return _has_quantum_number_f;
96 : }
97 :
98 : template <typename Derived>
99 18211 : bool Basis<Derived>::has_quantum_number_m() const {
100 18211 : return _has_quantum_number_m;
101 : }
102 :
103 : template <typename Derived>
104 9 : bool Basis<Derived>::has_parity() const {
105 9 : return _has_parity;
106 : }
107 :
108 : template <typename Derived>
109 153 : const Derived &Basis<Derived>::derived() const {
110 153 : return static_cast<const Derived &>(*this);
111 : }
112 :
113 : template <typename Derived>
114 72 : const typename Basis<Derived>::ketvec_t &Basis<Derived>::get_kets() const {
115 72 : return kets;
116 : }
117 :
118 : template <typename Derived>
119 : const Eigen::SparseMatrix<typename Basis<Derived>::scalar_t, Eigen::RowMajor> &
120 608 : Basis<Derived>::get_coefficients() const {
121 608 : return coefficients.matrix;
122 : }
123 :
124 : template <typename Derived>
125 : Eigen::SparseMatrix<typename Basis<Derived>::scalar_t, Eigen::RowMajor> &
126 0 : Basis<Derived>::get_coefficients() {
127 0 : return coefficients.matrix;
128 : }
129 :
130 : template <typename Derived>
131 34 : int Basis<Derived>::get_ket_index_from_ket(std::shared_ptr<const ket_t> ket) const {
132 34 : if (ket_to_ket_index.count(ket) == 0) {
133 0 : return -1;
134 : }
135 34 : return ket_to_ket_index.at(ket);
136 : }
137 :
138 : template <typename Derived>
139 : Eigen::VectorX<typename Basis<Derived>::scalar_t>
140 11 : Basis<Derived>::get_amplitudes(std::shared_ptr<const ket_t> ket) const {
141 11 : int ket_index = get_ket_index_from_ket(ket);
142 11 : if (ket_index < 0) {
143 0 : throw std::invalid_argument("The ket does not belong to the basis.");
144 : }
145 : // The following line is a more efficient alternative to
146 : // "get_amplitudes(get_canonical_state_from_ket(ket)).transpose()"
147 22 : return coefficients.matrix.row(ket_index);
148 : }
149 :
150 : template <typename Derived>
151 : Eigen::SparseMatrix<typename Basis<Derived>::scalar_t, Eigen::RowMajor>
152 0 : Basis<Derived>::get_amplitudes(std::shared_ptr<const Derived> other) const {
153 0 : return other->coefficients.matrix.adjoint() * coefficients.matrix;
154 : }
155 :
156 : template <typename Derived>
157 : Eigen::VectorX<typename Basis<Derived>::real_t>
158 11 : Basis<Derived>::get_overlaps(std::shared_ptr<const ket_t> ket) const {
159 11 : return get_amplitudes(ket).cwiseAbs2();
160 : }
161 :
162 : template <typename Derived>
163 : Eigen::SparseMatrix<typename Basis<Derived>::real_t, Eigen::RowMajor>
164 0 : Basis<Derived>::get_overlaps(std::shared_ptr<const Derived> other) const {
165 0 : return get_amplitudes(other).cwiseAbs2();
166 : }
167 :
168 : template <typename Derived>
169 0 : typename Basis<Derived>::real_t Basis<Derived>::get_quantum_number_f(size_t state_index) const {
170 0 : real_t quantum_number_f = state_index_to_quantum_number_f.at(state_index);
171 0 : if (quantum_number_f == std::numeric_limits<real_t>::max()) {
172 0 : throw std::invalid_argument("The state does not have a well-defined quantum number f.");
173 : }
174 0 : return quantum_number_f;
175 : }
176 :
177 : template <typename Derived>
178 18244 : typename Basis<Derived>::real_t Basis<Derived>::get_quantum_number_m(size_t state_index) const {
179 18244 : real_t quantum_number_m = state_index_to_quantum_number_m.at(state_index);
180 18244 : if (quantum_number_m == std::numeric_limits<real_t>::max()) {
181 0 : throw std::invalid_argument("The state does not have a well-defined quantum number m.");
182 : }
183 18244 : return quantum_number_m;
184 : }
185 :
186 : template <typename Derived>
187 43 : Parity Basis<Derived>::get_parity(size_t state_index) const {
188 43 : Parity parity = state_index_to_parity.at(state_index);
189 43 : if (parity == Parity::UNKNOWN) {
190 0 : throw std::invalid_argument("The state does not have a well-defined parity.");
191 : }
192 43 : return parity;
193 : }
194 :
195 : template <typename Derived>
196 : std::shared_ptr<const typename Basis<Derived>::ket_t>
197 38 : Basis<Derived>::get_corresponding_ket(size_t state_index) const {
198 38 : size_t ket_index = state_index_to_ket_index.at(state_index);
199 38 : if (ket_index == std::numeric_limits<int>::max()) {
200 0 : throw std::invalid_argument("The state does not belong to a ket in a well-defined way.");
201 : }
202 38 : return kets[ket_index];
203 : }
204 :
205 : template <typename Derived>
206 : std::shared_ptr<const typename Basis<Derived>::ket_t>
207 0 : Basis<Derived>::get_corresponding_ket(std::shared_ptr<const Derived> /*state*/) const {
208 0 : throw std::runtime_error("Not implemented yet.");
209 : }
210 :
211 : template <typename Derived>
212 4 : std::shared_ptr<const Derived> Basis<Derived>::get_state(size_t state_index) const {
213 : // Create a copy of the current object
214 4 : auto restricted = std::make_shared<Derived>(derived());
215 :
216 : // Restrict the copy to the state with the largest overlap
217 4 : restricted->coefficients.matrix = restricted->coefficients.matrix.col(state_index);
218 :
219 4 : std::fill(restricted->ket_index_to_state_index.begin(),
220 4 : restricted->ket_index_to_state_index.end(), std::numeric_limits<int>::max());
221 4 : restricted->ket_index_to_state_index[state_index_to_ket_index[state_index]] = 0;
222 :
223 4 : restricted->state_index_to_quantum_number_f = {state_index_to_quantum_number_f[state_index]};
224 4 : restricted->state_index_to_quantum_number_m = {state_index_to_quantum_number_m[state_index]};
225 4 : restricted->state_index_to_parity = {state_index_to_parity[state_index]};
226 4 : restricted->state_index_to_ket_index = {state_index_to_ket_index[state_index]};
227 :
228 8 : restricted->_has_quantum_number_f =
229 4 : restricted->state_index_to_quantum_number_f[0] != std::numeric_limits<real_t>::max();
230 8 : restricted->_has_quantum_number_m =
231 4 : restricted->state_index_to_quantum_number_m[0] != std::numeric_limits<real_t>::max();
232 4 : restricted->_has_parity = restricted->state_index_to_parity[0] != Parity::UNKNOWN;
233 :
234 8 : return restricted;
235 4 : }
236 :
237 : template <typename Derived>
238 : std::shared_ptr<const typename Basis<Derived>::ket_t>
239 0 : Basis<Derived>::get_ket(size_t ket_index) const {
240 0 : return kets[ket_index];
241 : }
242 :
243 : template <typename Derived>
244 2 : std::shared_ptr<const Derived> Basis<Derived>::get_corresponding_state(size_t ket_index) const {
245 2 : size_t state_index = ket_index_to_state_index.at(ket_index);
246 2 : if (state_index == std::numeric_limits<int>::max()) {
247 0 : throw std::runtime_error("The ket does not belong to a state in a well-defined way.");
248 : }
249 2 : return get_state(state_index);
250 : }
251 :
252 : template <typename Derived>
253 : std::shared_ptr<const Derived>
254 2 : Basis<Derived>::get_corresponding_state(std::shared_ptr<const ket_t> ket) const {
255 2 : int ket_index = get_ket_index_from_ket(ket);
256 2 : if (ket_index < 0) {
257 0 : throw std::invalid_argument("The ket does not belong to the basis.");
258 : }
259 2 : return get_corresponding_state(ket_index);
260 : }
261 :
262 : template <typename Derived>
263 2 : size_t Basis<Derived>::get_corresponding_state_index(size_t ket_index) const {
264 2 : int state_index = ket_index_to_state_index.at(ket_index);
265 2 : if (state_index == std::numeric_limits<int>::max()) {
266 0 : throw std::runtime_error("The ket does not belong to a state in a well-defined way.");
267 : }
268 2 : return state_index;
269 : }
270 :
271 : template <typename Derived>
272 2 : size_t Basis<Derived>::get_corresponding_state_index(std::shared_ptr<const ket_t> ket) const {
273 2 : int ket_index = get_ket_index_from_ket(ket);
274 2 : if (ket_index < 0) {
275 0 : throw std::invalid_argument("The ket does not belong to the basis.");
276 : }
277 2 : return get_corresponding_state_index(ket_index);
278 : }
279 :
280 : template <typename Derived>
281 0 : size_t Basis<Derived>::get_corresponding_ket_index(size_t state_index) const {
282 0 : int ket_index = state_index_to_ket_index.at(state_index);
283 0 : if (ket_index == std::numeric_limits<int>::max()) {
284 0 : throw std::runtime_error("The state does not belong to a ket in a well-defined way.");
285 : }
286 0 : return ket_index;
287 : }
288 :
289 : template <typename Derived>
290 0 : size_t Basis<Derived>::get_corresponding_ket_index(std::shared_ptr<const Derived> /*state*/) const {
291 0 : throw std::runtime_error("Not implemented yet.");
292 : }
293 :
294 : template <typename Derived>
295 : std::shared_ptr<const Derived>
296 19 : Basis<Derived>::get_canonical_state_from_ket(size_t ket_index) const {
297 : // Create a copy of the current object
298 19 : auto created = std::make_shared<Derived>(derived());
299 :
300 : // Fill the copy with the state corresponding to the ket index
301 19 : created->coefficients.matrix =
302 38 : Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>(coefficients.matrix.rows(), 1);
303 19 : created->coefficients.matrix.coeffRef(ket_index, 0) = 1;
304 19 : created->coefficients.matrix.makeCompressed();
305 :
306 19 : std::fill(created->ket_index_to_state_index.begin(), created->ket_index_to_state_index.end(),
307 19 : std::numeric_limits<int>::max());
308 19 : created->ket_index_to_state_index[ket_index] = 0;
309 :
310 19 : created->state_index_to_quantum_number_f = {kets[ket_index]->get_quantum_number_f()};
311 19 : created->state_index_to_quantum_number_m = {kets[ket_index]->get_quantum_number_m()};
312 19 : created->state_index_to_parity = {kets[ket_index]->get_parity()};
313 19 : created->state_index_to_ket_index = {ket_index};
314 :
315 38 : created->_has_quantum_number_f =
316 19 : created->state_index_to_quantum_number_f[0] != std::numeric_limits<real_t>::max();
317 38 : created->_has_quantum_number_m =
318 19 : created->state_index_to_quantum_number_m[0] != std::numeric_limits<real_t>::max();
319 19 : created->_has_parity = created->state_index_to_parity[0] != Parity::UNKNOWN;
320 :
321 38 : return created;
322 19 : }
323 :
324 : template <typename Derived>
325 : std::shared_ptr<const Derived>
326 19 : Basis<Derived>::get_canonical_state_from_ket(std::shared_ptr<const ket_t> ket) const {
327 19 : int ket_index = get_ket_index_from_ket(ket);
328 19 : if (ket_index < 0) {
329 0 : throw std::invalid_argument("The ket does not belong to the basis.");
330 : }
331 19 : return get_canonical_state_from_ket(ket_index);
332 : }
333 :
334 : template <typename Derived>
335 6 : typename Basis<Derived>::Iterator Basis<Derived>::begin() const {
336 6 : return kets.begin();
337 : }
338 :
339 : template <typename Derived>
340 6 : typename Basis<Derived>::Iterator Basis<Derived>::end() const {
341 6 : return kets.end();
342 : }
343 :
344 : template <typename Derived>
345 12 : Basis<Derived>::Iterator::Iterator(typename ketvec_t::const_iterator it) : it{std::move(it)} {}
346 :
347 : template <typename Derived>
348 191 : bool Basis<Derived>::Iterator::operator!=(const Iterator &other) const {
349 191 : return other.it != it;
350 : }
351 :
352 : template <typename Derived>
353 185 : std::shared_ptr<const typename Basis<Derived>::ket_t> Basis<Derived>::Iterator::operator*() const {
354 185 : return *it;
355 : }
356 :
357 : template <typename Derived>
358 185 : typename Basis<Derived>::Iterator &Basis<Derived>::Iterator::operator++() {
359 185 : ++it;
360 185 : return *this;
361 : }
362 :
363 : template <typename Derived>
364 3272 : size_t Basis<Derived>::get_number_of_states() const {
365 3272 : return coefficients.matrix.cols();
366 : }
367 :
368 : template <typename Derived>
369 419 : size_t Basis<Derived>::get_number_of_kets() const {
370 419 : return coefficients.matrix.rows();
371 : }
372 :
373 : template <typename Derived>
374 : const Transformation<typename Basis<Derived>::scalar_t> &
375 58 : Basis<Derived>::get_transformation() const {
376 58 : return coefficients;
377 : }
378 :
379 : template <typename Derived>
380 : Transformation<typename Basis<Derived>::scalar_t>
381 0 : Basis<Derived>::get_rotator(real_t alpha, real_t beta, real_t gamma) const {
382 0 : Transformation<scalar_t> transformation{{static_cast<Eigen::Index>(coefficients.matrix.rows()),
383 : static_cast<Eigen::Index>(coefficients.matrix.rows())},
384 : {TransformationType::ROTATE}};
385 :
386 0 : std::vector<Eigen::Triplet<scalar_t>> entries;
387 :
388 0 : for (size_t idx_initial = 0; idx_initial < kets.size(); ++idx_initial) {
389 0 : real_t f = kets[idx_initial]->get_quantum_number_f();
390 0 : real_t m_initial = kets[idx_initial]->get_quantum_number_m();
391 :
392 0 : assert(2 * f == std::floor(2 * f) && f >= 0);
393 0 : assert(2 * m_initial == std::floor(2 * m_initial) && m_initial >= -f && m_initial <= f);
394 :
395 0 : for (real_t m_final = -f; m_final <= f; // NOSONAR m_final is precisely representable
396 : ++m_final) {
397 0 : auto val = wigner::wigner_uppercase_d_matrix<scalar_t>(f, m_initial, m_final, alpha,
398 : beta, gamma);
399 0 : size_t idx_final = get_ket_index_from_ket(
400 0 : kets[idx_initial]->get_ket_for_different_quantum_number_m(m_final));
401 0 : entries.emplace_back(idx_final, idx_initial, val);
402 : }
403 : }
404 :
405 0 : transformation.matrix.setFromTriplets(entries.begin(), entries.end());
406 0 : transformation.matrix.makeCompressed();
407 :
408 0 : return transformation;
409 0 : }
410 :
411 : template <typename Derived>
412 1 : Sorting Basis<Derived>::get_sorter(const std::vector<TransformationType> &labels) const {
413 1 : perform_sorter_checks(labels);
414 :
415 : // Throw a meaningful error if sorting by energy is requested as this might be a common mistake
416 1 : if (std::find(labels.begin(), labels.end(), TransformationType::SORT_BY_ENERGY) !=
417 2 : labels.end()) {
418 0 : throw std::invalid_argument("States do not store the energy and thus can not be sorted by "
419 : "the energy. Use an energy operator instead.");
420 : }
421 :
422 : // Initialize transformation
423 1 : Sorting transformation;
424 1 : transformation.matrix.resize(coefficients.matrix.cols());
425 1 : transformation.matrix.setIdentity();
426 :
427 : // Get the sorter
428 1 : get_sorter_without_checks(labels, transformation);
429 :
430 : // Check if all labels have been used for sorting
431 1 : if (labels != transformation.transformation_type) {
432 0 : throw std::invalid_argument("The states could not be sorted by all the requested labels.");
433 : }
434 :
435 1 : return transformation;
436 0 : }
437 :
438 : template <typename Derived>
439 : std::vector<IndicesOfBlock>
440 1 : Basis<Derived>::get_indices_of_blocks(const std::vector<TransformationType> &labels) const {
441 1 : perform_sorter_checks(labels);
442 :
443 1 : std::set<TransformationType> unique_labels(labels.begin(), labels.end());
444 1 : perform_blocks_checks(unique_labels);
445 :
446 : // Get the blocks
447 1 : IndicesOfBlocksCreator blocks_creator({0, static_cast<size_t>(coefficients.matrix.cols())});
448 1 : get_indices_of_blocks_without_checks(unique_labels, blocks_creator);
449 :
450 2 : return blocks_creator.create();
451 1 : }
452 :
453 : template <typename Derived>
454 31 : void Basis<Derived>::get_sorter_without_checks(const std::vector<TransformationType> &labels,
455 : Sorting &transformation) const {
456 31 : constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
457 :
458 31 : int *perm_begin = transformation.matrix.indices().data();
459 31 : int *perm_end = perm_begin + coefficients.matrix.cols();
460 31 : const int *perm_back = perm_end - 1;
461 :
462 : // Sort the vector based on the requested labels
463 40930 : std::stable_sort(perm_begin, perm_end, [&](int a, int b) {
464 14413 : for (const auto &label : labels) {
465 9974 : switch (label) {
466 51 : case TransformationType::SORT_BY_PARITY:
467 51 : if (state_index_to_parity[a] != state_index_to_parity[b]) {
468 5504 : return state_index_to_parity[a] < state_index_to_parity[b];
469 : }
470 31 : break;
471 9923 : case TransformationType::SORT_BY_QUANTUM_NUMBER_M:
472 19846 : if (std::abs(state_index_to_quantum_number_m[a] -
473 19846 : state_index_to_quantum_number_m[b]) > numerical_precision) {
474 5484 : return state_index_to_quantum_number_m[a] < state_index_to_quantum_number_m[b];
475 : }
476 4439 : break;
477 0 : case TransformationType::SORT_BY_QUANTUM_NUMBER_F:
478 0 : if (std::abs(state_index_to_quantum_number_f[a] -
479 0 : state_index_to_quantum_number_f[b]) > numerical_precision) {
480 0 : return state_index_to_quantum_number_f[a] < state_index_to_quantum_number_f[b];
481 : }
482 0 : break;
483 0 : case TransformationType::SORT_BY_KET:
484 0 : if (state_index_to_ket_index[a] != state_index_to_ket_index[b]) {
485 0 : return state_index_to_ket_index[a] < state_index_to_ket_index[b];
486 : }
487 0 : break;
488 0 : default:
489 0 : std::abort(); // Can't happen because of previous checks
490 : }
491 : }
492 4439 : return false; // Elements are equal
493 : });
494 :
495 : // Check for invalid values and add transformation types
496 63 : for (const auto &label : labels) {
497 32 : switch (label) {
498 1 : case TransformationType::SORT_BY_PARITY:
499 1 : if (state_index_to_parity[*perm_back] == Parity::UNKNOWN) {
500 0 : throw std::invalid_argument(
501 : "States cannot be labeled and thus not sorted by the parity.");
502 : }
503 1 : transformation.transformation_type.push_back(TransformationType::SORT_BY_PARITY);
504 1 : break;
505 31 : case TransformationType::SORT_BY_QUANTUM_NUMBER_M:
506 31 : if (state_index_to_quantum_number_m[*perm_back] == std::numeric_limits<real_t>::max()) {
507 0 : throw std::invalid_argument(
508 : "States cannot be labeled and thus not sorted by the quantum number m.");
509 : }
510 31 : transformation.transformation_type.push_back(
511 31 : TransformationType::SORT_BY_QUANTUM_NUMBER_M);
512 31 : break;
513 0 : case TransformationType::SORT_BY_QUANTUM_NUMBER_F:
514 0 : if (state_index_to_quantum_number_f[*perm_back] == std::numeric_limits<real_t>::max()) {
515 0 : throw std::invalid_argument(
516 : "States cannot be labeled and thus not sorted by the quantum number f.");
517 : }
518 0 : transformation.transformation_type.push_back(
519 0 : TransformationType::SORT_BY_QUANTUM_NUMBER_F);
520 0 : break;
521 0 : case TransformationType::SORT_BY_KET:
522 0 : if (state_index_to_ket_index[*perm_back] == std::numeric_limits<int>::max()) {
523 0 : throw std::invalid_argument(
524 : "States cannot be labeled and thus not sorted by kets.");
525 : }
526 0 : transformation.transformation_type.push_back(TransformationType::SORT_BY_KET);
527 0 : break;
528 0 : default:
529 0 : std::abort(); // Can't happen because of previous checks
530 : }
531 : }
532 31 : }
533 :
534 : template <typename Derived>
535 31 : void Basis<Derived>::get_indices_of_blocks_without_checks(
536 : const std::set<TransformationType> &unique_labels,
537 : IndicesOfBlocksCreator &blocks_creator) const {
538 31 : constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
539 :
540 31 : auto last_quantum_number_f = state_index_to_quantum_number_f[0];
541 31 : auto last_quantum_number_m = state_index_to_quantum_number_m[0];
542 31 : auto last_parity = state_index_to_parity[0];
543 31 : auto last_ket = state_index_to_ket_index[0];
544 :
545 1845 : for (int i = 0; i < coefficients.matrix.cols(); ++i) {
546 3545 : for (auto label : unique_labels) {
547 1825 : if (label == TransformationType::SORT_BY_QUANTUM_NUMBER_F) {
548 0 : if (std::abs(state_index_to_quantum_number_f[i] - last_quantum_number_f) >
549 : numerical_precision) {
550 0 : blocks_creator.add(i);
551 0 : break;
552 : }
553 1825 : } else if (label == TransformationType::SORT_BY_QUANTUM_NUMBER_M) {
554 1814 : if (std::abs(state_index_to_quantum_number_m[i] - last_quantum_number_m) >
555 : numerical_precision) {
556 94 : blocks_creator.add(i);
557 94 : break;
558 : }
559 11 : } else if (label == TransformationType::SORT_BY_PARITY) {
560 11 : if (state_index_to_parity[i] != last_parity) {
561 0 : blocks_creator.add(i);
562 0 : break;
563 : }
564 0 : } else if (label == TransformationType::SORT_BY_KET) {
565 0 : if (state_index_to_ket_index[i] != last_ket) {
566 0 : blocks_creator.add(i);
567 0 : break;
568 : }
569 : }
570 : }
571 1814 : last_quantum_number_f = state_index_to_quantum_number_f[i];
572 1814 : last_quantum_number_m = state_index_to_quantum_number_m[i];
573 1814 : last_parity = state_index_to_parity[i];
574 1814 : last_ket = state_index_to_ket_index[i];
575 : }
576 31 : }
577 :
578 : template <typename Derived>
579 73 : std::shared_ptr<const Derived> Basis<Derived>::transformed(const Sorting &transformation) const {
580 : // Create a copy of the current object
581 73 : auto transformed = std::make_shared<Derived>(derived());
582 :
583 : // Apply the transformation
584 73 : transformed->coefficients.matrix = transformed->coefficients.matrix * transformation.matrix;
585 73 : transformed->coefficients.transformation_type = transformation.transformation_type;
586 :
587 73 : transformed->state_index_to_quantum_number_f.resize(transformation.matrix.size());
588 73 : transformed->state_index_to_quantum_number_m.resize(transformation.matrix.size());
589 73 : transformed->state_index_to_parity.resize(transformation.matrix.size());
590 73 : transformed->state_index_to_ket_index.resize(transformation.matrix.size());
591 :
592 4910 : for (int i = 0; i < transformation.matrix.size(); ++i) {
593 4837 : transformed->state_index_to_quantum_number_f[i] =
594 4837 : state_index_to_quantum_number_f[transformation.matrix.indices()[i]];
595 4837 : transformed->state_index_to_quantum_number_m[i] =
596 4837 : state_index_to_quantum_number_m[transformation.matrix.indices()[i]];
597 4837 : transformed->state_index_to_parity[i] =
598 4837 : state_index_to_parity[transformation.matrix.indices()[i]];
599 4837 : transformed->state_index_to_ket_index[i] =
600 4837 : state_index_to_ket_index[transformation.matrix.indices()[i]];
601 4837 : transformed->ket_index_to_state_index
602 4837 : [state_index_to_ket_index[transformation.matrix.indices()[i]]] = i;
603 : }
604 :
605 146 : return transformed;
606 73 : }
607 :
608 : template <typename Derived>
609 : std::shared_ptr<const Derived>
610 57 : Basis<Derived>::transformed(const Transformation<scalar_t> &transformation) const {
611 : // TODO why is "numerical_precision = 100 * std::sqrt(coefficients.matrix.rows()) *
612 : // std::numeric_limits<real_t>::epsilon()" too small for figuring out whether m is conserved?
613 57 : real_t numerical_precision = 0.001;
614 :
615 : // If the transformation is a rotation, it should be a rotation and nothing else
616 57 : bool is_rotation = false;
617 114 : for (auto t : transformation.transformation_type) {
618 57 : if (t == TransformationType::ROTATE) {
619 0 : is_rotation = true;
620 0 : break;
621 : }
622 : }
623 57 : if (is_rotation && transformation.transformation_type.size() != 1) {
624 0 : throw std::invalid_argument("A rotation can not be combined with other transformations.");
625 : }
626 :
627 : // To apply a rotation, the object must only be sorted but other transformations are not allowed
628 57 : if (is_rotation) {
629 0 : for (auto t : coefficients.transformation_type) {
630 0 : if (!utils::is_sorting(t)) {
631 0 : throw std::runtime_error(
632 : "If the object was transformed by a different transformation "
633 : "than sorting, it can not be rotated.");
634 : }
635 : }
636 : }
637 :
638 : // Create a copy of the current object
639 57 : auto transformed = std::make_shared<Derived>(derived());
640 :
641 : // Apply the transformation
642 : // If a quantum number turns out to be conserved by the transformation, it will be
643 : // rounded to the nearest half integer to avoid loss of numerical_precision.
644 57 : transformed->coefficients.matrix = coefficients.matrix * transformation.matrix;
645 57 : transformed->coefficients.transformation_type = transformation.transformation_type;
646 :
647 57 : Eigen::SparseMatrix<real_t> probs = transformation.matrix.cwiseAbs2().transpose();
648 :
649 : {
650 114 : auto map = Eigen::Map<const Eigen::VectorX<real_t>>(state_index_to_quantum_number_f.data(),
651 57 : state_index_to_quantum_number_f.size());
652 57 : Eigen::VectorX<real_t> val = probs * map;
653 57 : Eigen::VectorX<real_t> sq = probs * map.cwiseAbs2();
654 57 : Eigen::VectorX<real_t> diff = (val.cwiseAbs2() - sq).cwiseAbs();
655 57 : transformed->state_index_to_quantum_number_f.resize(probs.rows());
656 :
657 2151 : for (size_t i = 0; i < transformed->state_index_to_quantum_number_f.size(); ++i) {
658 2094 : if (diff[i] < numerical_precision) {
659 376 : transformed->state_index_to_quantum_number_f[i] = std::round(val[i] * 2) / 2;
660 : } else {
661 1718 : transformed->state_index_to_quantum_number_f[i] =
662 1718 : std::numeric_limits<real_t>::max();
663 1718 : transformed->_has_quantum_number_f = false;
664 : }
665 : }
666 57 : }
667 :
668 : {
669 114 : auto map = Eigen::Map<const Eigen::VectorX<real_t>>(state_index_to_quantum_number_m.data(),
670 57 : state_index_to_quantum_number_m.size());
671 57 : Eigen::VectorX<real_t> val = probs * map;
672 57 : Eigen::VectorX<real_t> sq = probs * map.cwiseAbs2();
673 57 : Eigen::VectorX<real_t> diff = (val.cwiseAbs2() - sq).cwiseAbs();
674 57 : transformed->state_index_to_quantum_number_m.resize(probs.rows());
675 :
676 2151 : for (size_t i = 0; i < transformed->state_index_to_quantum_number_m.size(); ++i) {
677 2094 : if (diff[i] < numerical_precision) {
678 1742 : transformed->state_index_to_quantum_number_m[i] = std::round(val[i] * 2) / 2;
679 : } else {
680 352 : transformed->state_index_to_quantum_number_m[i] =
681 352 : std::numeric_limits<real_t>::max();
682 352 : transformed->_has_quantum_number_m = false;
683 : }
684 : }
685 57 : }
686 :
687 : {
688 : using utype = std::underlying_type<Parity>::type;
689 57 : Eigen::VectorX<real_t> map(state_index_to_parity.size());
690 2431 : for (size_t i = 0; i < state_index_to_parity.size(); ++i) {
691 2374 : map[i] = static_cast<utype>(state_index_to_parity[i]);
692 : }
693 57 : Eigen::VectorX<real_t> val = probs * map;
694 57 : Eigen::VectorX<real_t> sq = probs * map.cwiseAbs2();
695 57 : Eigen::VectorX<real_t> diff = (val.cwiseAbs2() - sq).cwiseAbs();
696 57 : transformed->state_index_to_parity.resize(probs.rows());
697 :
698 2151 : for (size_t i = 0; i < transformed->state_index_to_parity.size(); ++i) {
699 2094 : if (diff[i] < numerical_precision) {
700 800 : transformed->state_index_to_parity[i] = static_cast<Parity>(std::lround(val[i]));
701 : } else {
702 1293 : transformed->state_index_to_parity[i] = Parity::UNKNOWN;
703 1293 : transformed->_has_parity = false;
704 : }
705 : }
706 57 : }
707 :
708 : {
709 : // In the following, we obtain a bijective map between state index and ket index.
710 :
711 : // Find the maximum value in each row and column
712 57 : std::vector<real_t> max_in_row(transformed->coefficients.matrix.rows(), 0);
713 57 : std::vector<real_t> max_in_col(transformed->coefficients.matrix.cols(), 0);
714 2431 : for (int row = 0; row < transformed->coefficients.matrix.outerSize(); ++row) {
715 2372 : for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
716 2373 : transformed->coefficients.matrix, row);
717 43691 : it; ++it) {
718 41351 : real_t val = std::pow(std::abs(it.value()), 2);
719 41413 : max_in_row[row] = std::max(max_in_row[row], val);
720 41439 : max_in_col[it.col()] = std::max(max_in_col[it.col()], val);
721 : }
722 : }
723 :
724 : // Use the maximum values to define a cost for a sub-optimal mapping
725 57 : std::vector<real_t> costs;
726 57 : std::vector<std::pair<int, int>> mappings;
727 57 : costs.reserve(transformed->coefficients.matrix.nonZeros());
728 57 : mappings.reserve(transformed->coefficients.matrix.nonZeros());
729 2431 : for (int row = 0; row < transformed->coefficients.matrix.outerSize(); ++row) {
730 2374 : for (typename Eigen::SparseMatrix<scalar_t, Eigen::RowMajor>::InnerIterator it(
731 2374 : transformed->coefficients.matrix, row);
732 43832 : it; ++it) {
733 41442 : real_t val = std::pow(std::abs(it.value()), 2);
734 41513 : real_t cost = max_in_row[row] + max_in_col[it.col()] - 2 * val;
735 41430 : costs.push_back(cost);
736 41478 : mappings.push_back({row, it.col()});
737 : }
738 : }
739 :
740 : // Obtain from the costs in which order the mappings should be considered
741 57 : std::vector<size_t> order(costs.size());
742 57 : std::iota(order.begin(), order.end(), 0);
743 57 : std::sort(order.begin(), order.end(),
744 522153 : [&](size_t a, size_t b) { return costs[a] < costs[b]; });
745 :
746 : // Fill ket_index_to_state_index with invalid values as there can be more kets than states
747 57 : std::fill(transformed->ket_index_to_state_index.begin(),
748 57 : transformed->ket_index_to_state_index.end(), std::numeric_limits<int>::max());
749 :
750 : // Generate the bijective map
751 57 : std::vector<bool> row_used(transformed->coefficients.matrix.rows(), false);
752 57 : std::vector<bool> col_used(transformed->coefficients.matrix.cols(), false);
753 57 : int num_used = 0;
754 5321 : for (size_t idx : order) {
755 5318 : int row = mappings[idx].first; // corresponds to the ket index
756 5322 : int col = mappings[idx].second; // corresponds to the state index
757 5320 : if (!row_used[row] && !col_used[col]) {
758 2094 : row_used[row] = true;
759 2093 : col_used[col] = true;
760 2094 : num_used++;
761 2094 : transformed->state_index_to_ket_index[col] = row;
762 2094 : transformed->ket_index_to_state_index[row] = col;
763 : }
764 5317 : if (num_used == transformed->coefficients.matrix.cols()) {
765 55 : break;
766 : }
767 : }
768 57 : assert(num_used == transformed->coefficients.matrix.cols());
769 57 : }
770 :
771 114 : return transformed;
772 57 : }
773 :
774 : template <typename Derived>
775 9749 : size_t Basis<Derived>::hash::operator()(const std::shared_ptr<const ket_t> &k) const {
776 9749 : return typename ket_t::hash()(*k);
777 : }
778 :
779 : template <typename Derived>
780 68 : bool Basis<Derived>::equal_to::operator()(const std::shared_ptr<const ket_t> &lhs,
781 : const std::shared_ptr<const ket_t> &rhs) const {
782 68 : return *lhs == *rhs;
783 : }
784 :
785 : // Explicit instantiations
786 : template class Basis<BasisAtom<double>>;
787 : template class Basis<BasisAtom<std::complex<double>>>;
788 : template class Basis<BasisPair<double>>;
789 : template class Basis<BasisPair<std::complex<double>>>;
790 : } // namespace pairinteraction
|
