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