Line data Source code
1 : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/basis/BasisPairCreator.hpp"
5 :
6 : #include "pairinteraction/basis/BasisAtom.hpp"
7 : #include "pairinteraction/basis/BasisPair.hpp"
8 : #include "pairinteraction/ket/KetPair.hpp"
9 : #include "pairinteraction/utils/TaskControl.hpp"
10 : #include "pairinteraction/utils/hash.hpp"
11 :
12 : #include <algorithm>
13 : #include <array>
14 : #include <cassert>
15 : #include <cmath>
16 : #include <limits>
17 : #include <memory>
18 : #include <stdexcept>
19 : #include <string>
20 : #include <unordered_map>
21 : #include <vector>
22 :
23 : namespace pairinteraction {
24 : template <typename Scalar>
25 553 : BasisPairCreator<Scalar>::BasisPairCreator(std::shared_ptr<const BasisAtom<Scalar>> basis1,
26 : std::shared_ptr<const BasisAtom<Scalar>> basis2)
27 553 : : basis1(std::move(basis1)), basis2(std::move(basis2)) {}
28 :
29 : template <typename Scalar>
30 232 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_energy(real_t min, real_t max) {
31 232 : range_energy = {min, max};
32 232 : return *this;
33 : }
34 :
35 : template <typename Scalar>
36 57 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_quantum_number_m(real_t min,
37 : real_t max) {
38 57 : range_quantum_number_m = {min, max};
39 57 : return *this;
40 : }
41 :
42 : template <typename Scalar>
43 8 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_parity_under_inversion(int value) {
44 8 : if (value != 1 && value != -1) {
45 1 : throw std::invalid_argument("The parity must be +1 or -1.");
46 : }
47 7 : parity_under_inversion = value;
48 7 : return *this;
49 : }
50 :
51 : template <typename Scalar>
52 10 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_parity_under_permutation(int value) {
53 10 : if (value != 1 && value != -1) {
54 1 : throw std::invalid_argument("The parity must be +1 or -1.");
55 : }
56 9 : parity_under_permutation = value;
57 9 : return *this;
58 : }
59 :
60 : template <typename Scalar>
61 551 : std::shared_ptr<const BasisPair<Scalar>> BasisPairCreator<Scalar>::create() const {
62 551 : set_task_status("Constructing pair basis...");
63 :
64 551 : constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
65 551 : const bool has_symmetry_restriction =
66 551 : parity_under_inversion.has_value() || parity_under_permutation.has_value();
67 :
68 : // This ensures that a one-atom state can be identified across both atoms by its state index
69 551 : if (has_symmetry_restriction && basis1 != basis2) {
70 2 : throw std::invalid_argument(
71 : "Parity restrictions require the same BasisAtom for both atoms, because "
72 : "symmetrization is only defined for two identical atoms.");
73 : }
74 :
75 549 : std::optional<int> inferred_product_of_parities;
76 549 : if (parity_under_inversion.has_value() && parity_under_permutation.has_value()) {
77 1 : inferred_product_of_parities = *parity_under_inversion * *parity_under_permutation;
78 : }
79 :
80 : // Sorting is required for the binary search of the energetically allowed range
81 781 : if (range_energy.is_finite() &&
82 232 : (!basis1->is_sorted_by_energy() || !basis2->is_sorted_by_energy())) {
83 3 : throw std::invalid_argument(
84 : "The energy can only be restricted if the states of the BasisAtom have well-defined "
85 : "energies and are sorted by energy. Consider using the eigenbasis of a diagonalized "
86 : "SystemAtom.");
87 : }
88 :
89 : // Construct the canonical basis that contains all KetPair objects with allowed energies and
90 : // quantum numbers. States without a well-defined energy are marked by the maximum value.
91 : static constexpr real_t unknown_energy = std::numeric_limits<real_t>::max();
92 1092 : auto get_energies = [](const auto &basis) {
93 1092 : std::vector<real_t> energies(basis->get_number_of_states(), unknown_energy);
94 1092 : if (basis->has_energy()) {
95 112903 : for (size_t state_index = 0; state_index < energies.size(); ++state_index) {
96 111824 : energies[state_index] = basis->get_energy(state_index);
97 : }
98 : }
99 1092 : return energies;
100 0 : };
101 546 : const std::vector<real_t> eigenenergies1 = get_energies(basis1);
102 546 : const std::vector<real_t> eigenenergies2 = get_energies(basis2);
103 546 : const real_t *eigenenergies2_begin = eigenenergies2.data();
104 546 : const real_t *eigenenergies2_end = eigenenergies2_begin + eigenenergies2.size();
105 :
106 : // The quantum number m of a pair state is only well-defined if it is well-defined for both
107 : // atoms
108 546 : const bool has_quantum_number_m =
109 546 : basis1->has_quantum_number("m") && basis2->has_quantum_number("m");
110 546 : if (!has_quantum_number_m && range_quantum_number_m.is_finite()) {
111 0 : throw std::invalid_argument(
112 : "The quantum number m must not be restricted because it is not well-defined.");
113 : }
114 :
115 546 : ketvec_t kets;
116 546 : kets.reserve(eigenenergies1.size() * eigenenergies2.size());
117 :
118 546 : typename basis_t::map_range_t state_index1_to_state_index_range2;
119 546 : state_index1_to_state_index_range2.reserve(eigenenergies1.size());
120 :
121 546 : typename basis_t::map_indices_t state_indices_to_ket_index;
122 :
123 546 : Eigen::Index state_index = 0;
124 : std::unordered_map<std::array<size_t, 2>, Eigen::Index, utils::hash<std::array<size_t, 2>>>
125 546 : ket_indices2state_index;
126 546 : std::vector<Eigen::Triplet<Scalar>> transformation_triplets;
127 546 : if (has_symmetry_restriction) {
128 13 : transformation_triplets.reserve(eigenenergies1.size() * eigenenergies2.size());
129 : }
130 546 : const double inverse_sqrt_two = 1 / std::sqrt(2.0);
131 :
132 : // Construct the symmetry transformation by recording the contribution of the pair state
133 : // |idx1, idx2> to the symmetrized basis. Because the two atoms are identical (enforced above),
134 : // a one-atom state is uniquely identified across both atoms by its state index alone.
135 327144 : auto construct_symmetry_transformation = [&](Eigen::Index row_index, size_t idx1, size_t idx2) {
136 : // Following https://doi.org/10.1088/1361-6455/aa743a, pair states |a, a> cannot be of even
137 : // parity.
138 65456 : if (idx1 == idx2 && (parity_under_inversion == 1 || parity_under_permutation == 1)) {
139 33118 : return;
140 : }
141 :
142 : // Map the (unordered) pair of one-atom state indices to the column index of the symmetrized
143 : // state it contributes to, creating a new column the first time the pair is encountered.
144 65084 : std::array<size_t, 2> ordered_indices{std::max(idx1, idx2), std::min(idx1, idx2)};
145 65084 : auto [iterator, inserted] =
146 65084 : ket_indices2state_index.try_emplace(ordered_indices, state_index);
147 65084 : if (inserted) {
148 32746 : ++state_index;
149 : }
150 65084 : Eigen::Index column_index = iterator->second;
151 :
152 : // A pair state |a, a> contributes with coefficient one.
153 65084 : if (idx1 == idx2) {
154 408 : transformation_triplets.emplace_back(row_index, column_index, 1);
155 408 : return;
156 : }
157 :
158 : // We let pair states with idx1 > idx2 contribute with coefficient 1/sqrt(2) and put the
159 : // phase into the contribution of the partner state with idx1 < idx2.
160 64676 : if (idx1 > idx2) {
161 32338 : transformation_triplets.emplace_back(row_index, column_index, inverse_sqrt_two);
162 32338 : return;
163 : }
164 :
165 : // Determine the phase of the contribution of the partner state with idx1 < idx2.
166 32338 : int phase = 0;
167 32338 : if (parity_under_permutation.has_value()) {
168 : // If both inversion and permutation are restricted, the earlier filter on the product
169 : // of parities already guarantees that the phases in the inversion- and
170 : // permutation-symmetric states are the same.
171 16252 : phase = -*parity_under_permutation;
172 : } else {
173 : // If only inversion is restricted, the phase is determined by the product of the
174 : // parities of the one-atom states and the specified inversion parity.
175 16086 : phase = -*parity_under_inversion *
176 32172 : static_cast<int>(basis1->get_quantum_number("parity", idx1)) *
177 16086 : static_cast<int>(basis2->get_quantum_number("parity", idx2));
178 : }
179 32338 : transformation_triplets.emplace_back(row_index, column_index, phase * inverse_sqrt_two);
180 : };
181 :
182 : // Loop only over states with an allowed energy
183 546 : size_t ket_index = 0;
184 55791 : for (size_t idx1 = 0; idx1 < eigenenergies1.size(); ++idx1) {
185 55245 : set_task_status("Constructing pair basis...");
186 :
187 : // Get the energetically allowed range of the second index
188 55245 : size_t min = 0;
189 55245 : size_t max = eigenenergies2.size();
190 55245 : if (range_energy.is_finite()) {
191 37003 : real_t min_val2 = range_energy.min() - eigenenergies1[idx1];
192 37003 : real_t max_val2 = range_energy.max() - eigenenergies1[idx1];
193 37003 : min =
194 37003 : std::distance(eigenenergies2_begin,
195 : std::lower_bound(eigenenergies2_begin, eigenenergies2_end, min_val2));
196 37003 : max =
197 37003 : std::distance(eigenenergies2_begin,
198 : std::upper_bound(eigenenergies2_begin, eigenenergies2_end, max_val2));
199 : }
200 55245 : state_index1_to_state_index_range2.try_emplace(idx1, typename basis_t::range_t(min, max));
201 :
202 : // Loop over the energetically allowed range of the second index
203 5412543 : for (size_t idx2 = min; idx2 < max; ++idx2) {
204 : // Get energy
205 2678681 : const real_t energy =
206 5357356 : (eigenenergies1[idx1] == unknown_energy || eigenenergies2[idx2] == unknown_energy)
207 2678937 : ? unknown_energy
208 2678419 : : eigenenergies1[idx1] + eigenenergies2[idx2];
209 2678681 : assert(!range_energy.is_finite() ||
210 : (energy >= range_energy.min() && energy <= range_energy.max()));
211 :
212 : // Check the parity of the product of the parities
213 2678681 : if (inferred_product_of_parities.has_value()) {
214 288 : if (static_cast<int>(basis1->get_quantum_number("parity", idx1)) *
215 288 : static_cast<int>(basis2->get_quantum_number("parity", idx2)) !=
216 144 : *inferred_product_of_parities) {
217 70303 : continue;
218 : }
219 : }
220 :
221 : // Get the quantum numbers and check the quantum number m
222 2678617 : std::unordered_map<std::string, double> quantum_numbers;
223 2678617 : if (has_quantum_number_m) {
224 2678355 : const real_t m =
225 2678355 : basis1->get_quantum_number("m", idx1) + basis2->get_quantum_number("m", idx2);
226 2808727 : if (range_quantum_number_m.is_finite() &&
227 130372 : (m < range_quantum_number_m.min() - numerical_precision ||
228 74569 : m > range_quantum_number_m.max() + numerical_precision)) {
229 70239 : continue;
230 : }
231 2608116 : quantum_numbers["m"] = m;
232 : }
233 :
234 : // Create a KetPair object
235 7825134 : auto ket = std::make_shared<ket_t>(
236 2608378 : typename ket_t::Private(), std::initializer_list<size_t>{idx1, idx2},
237 5216756 : std::initializer_list<std::shared_ptr<const BasisAtom<Scalar>>>{basis1, basis2},
238 2608378 : energy, std::move(quantum_numbers));
239 :
240 : // Store the KetPair object as a ket
241 2608378 : kets.emplace_back(std::move(ket));
242 2608378 : state_indices_to_ket_index.try_emplace(std::vector<size_t>{idx1, idx2}, ket_index);
243 :
244 2608378 : auto row_index = static_cast<Eigen::Index>(ket_index++);
245 2608378 : if (has_symmetry_restriction) {
246 65456 : construct_symmetry_transformation(row_index, idx1, idx2);
247 : }
248 : }
249 : }
250 :
251 546 : kets.shrink_to_fit();
252 :
253 546 : std::shared_ptr<const basis_t> basis = std::make_shared<basis_t>(
254 546 : typename basis_t::Private(), std::move(kets), std::move(state_index1_to_state_index_range2),
255 546 : std::move(state_indices_to_ket_index), basis1, basis2);
256 :
257 546 : if (!has_symmetry_restriction) {
258 533 : return basis;
259 : }
260 :
261 13 : transformation_triplets.shrink_to_fit();
262 :
263 26 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor> transformation_matrix(
264 13 : basis->get_number_of_states(), state_index);
265 13 : transformation_matrix.setFromTriplets(transformation_triplets.begin(),
266 13 : transformation_triplets.end());
267 :
268 13 : Eigen::Matrix<real_t, Eigen::Dynamic, 1> sum_of_squared_coefficients =
269 : transformation_matrix.cwiseAbs2().transpose() *
270 : Eigen::Matrix<real_t, Eigen::Dynamic, 1>::Ones(transformation_matrix.rows());
271 32759 : for (Eigen::Index column_index = 0; column_index < sum_of_squared_coefficients.size();
272 : ++column_index) {
273 32746 : if (std::abs(sum_of_squared_coefficients[column_index] - 1) > numerical_precision) {
274 0 : throw std::invalid_argument(
275 : "The basis could not be symmetrized. This likely means that the specified parity "
276 : "restrictions are invalid for the given one-atom systems.");
277 : }
278 : }
279 :
280 : // TODO: on the long run, construct the coefficient matrix directly
281 13 : return basis->transformed(transformation_matrix);
282 2608924 : }
283 :
284 : // Explicit instantiations
285 : template class BasisPairCreator<double>;
286 : template class BasisPairCreator<std::complex<double>>;
287 : } // namespace pairinteraction
|
