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/enums/Parity.hpp"
9 : #include "pairinteraction/ket/KetPair.hpp"
10 : #include "pairinteraction/system/SystemAtom.hpp"
11 : #include "pairinteraction/utils/TaskControl.hpp"
12 : #include "pairinteraction/utils/hash.hpp"
13 :
14 : #include <algorithm>
15 : #include <array>
16 : #include <cassert>
17 : #include <cmath>
18 : #include <limits>
19 : #include <memory>
20 : #include <stdexcept>
21 : #include <unordered_map>
22 :
23 : namespace pairinteraction {
24 : template <typename Scalar>
25 1030 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::add(const SystemAtom<Scalar> &system_atom) {
26 : // The system must be diagonalized and its eigenstates sorted by energy.
27 : // Sorting is required for the binary search of the energetically allowed range in create().
28 : // By default, System::diagonalize ensures this.
29 1030 : if (!system_atom.is_diagonal_and_sorted_by_energy()) {
30 0 : throw std::invalid_argument(
31 : "The system must be diagonalized and sorted by energy before it can be added. "
32 : "Consider calling diagonalize() on the SystemAtom which also sorts the eigenstates.");
33 : }
34 1030 : systems_atom.push_back(system_atom);
35 1030 : return *this;
36 : }
37 :
38 : template <typename Scalar>
39 222 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_energy(real_t min, real_t max) {
40 222 : range_energy = {min, max};
41 222 : return *this;
42 : }
43 :
44 : template <typename Scalar>
45 55 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_quantum_number_m(real_t min,
46 : real_t max) {
47 55 : range_quantum_number_m = {min, max};
48 55 : return *this;
49 : }
50 :
51 : template <typename Scalar>
52 9 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_parity_under_inversion(Parity value) {
53 9 : parity_under_inversion = value;
54 9 : return *this;
55 : }
56 :
57 : template <typename Scalar>
58 : BasisPairCreator<Scalar> &
59 11 : BasisPairCreator<Scalar>::restrict_parity_under_permutation(Parity value) {
60 11 : parity_under_permutation = value;
61 11 : return *this;
62 : }
63 :
64 : template <typename Scalar>
65 515 : std::shared_ptr<const BasisPair<Scalar>> BasisPairCreator<Scalar>::create() const {
66 515 : set_task_status("Constructing pair basis...");
67 :
68 515 : if (systems_atom.size() != 2) {
69 0 : throw std::invalid_argument("Two SystemAtom must be added before creating the BasisPair.");
70 : }
71 :
72 : // Only references to the systems are stored, so a system might have been changed since add()
73 1545 : for (const auto &system_atom : systems_atom) {
74 1030 : if (!system_atom.get().is_diagonal_and_sorted_by_energy()) {
75 0 : throw std::invalid_argument(
76 : "The systems must still be diagonalized and sorted by energy when the BasisPair is "
77 : "created. Do not change a SystemAtom after it has been added.");
78 : }
79 : }
80 :
81 515 : constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
82 515 : const bool has_symmetry_restriction =
83 515 : parity_under_inversion != Parity::UNKNOWN || parity_under_permutation != Parity::UNKNOWN;
84 :
85 : // This ensures that a one-atom state can be identified across both atoms by its state index
86 515 : if (has_symmetry_restriction && &systems_atom[0].get() != &systems_atom[1].get()) {
87 6 : throw std::invalid_argument(
88 : "Parity restrictions require the same SystemAtom to be added twice, because "
89 : "symmetrization is only defined for two identical atoms.");
90 : }
91 :
92 509 : Parity inferred_product_of_parities = Parity::UNKNOWN;
93 509 : if (parity_under_inversion != Parity::UNKNOWN && parity_under_permutation != Parity::UNKNOWN) {
94 1 : inferred_product_of_parities = static_cast<Parity>(
95 1 : static_cast<int>(parity_under_inversion) * static_cast<int>(parity_under_permutation));
96 : }
97 :
98 509 : const auto &system1 = systems_atom[0].get();
99 509 : const auto &system2 = systems_atom[1].get();
100 :
101 : // Construct the canonical basis that contains all KetPair objects with allowed energies and
102 : // quantum numbers
103 509 : auto basis1 = system1.get_basis();
104 509 : auto basis2 = system2.get_basis();
105 509 : auto eigenenergies1 = system1.get_eigenenergies();
106 509 : auto eigenenergies2 = system2.get_eigenenergies();
107 509 : real_t *eigenenergies2_begin = eigenenergies2.data();
108 509 : real_t *eigenenergies2_end = eigenenergies2_begin + eigenenergies2.size();
109 :
110 509 : ketvec_t kets;
111 509 : kets.reserve(eigenenergies1.size() * eigenenergies2.size());
112 :
113 509 : typename basis_t::map_range_t map_range_of_state_index2;
114 509 : map_range_of_state_index2.reserve(eigenenergies1.size());
115 :
116 509 : typename basis_t::map_indices_t state_indices_to_ket_index;
117 :
118 509 : Eigen::Index state_index = 0;
119 : std::unordered_map<std::array<size_t, 2>, Eigen::Index, utils::hash<std::array<size_t, 2>>>
120 509 : ket_indices2state_index;
121 509 : std::vector<Eigen::Triplet<Scalar>> transformation_triplets;
122 509 : if (has_symmetry_restriction) {
123 13 : transformation_triplets.reserve(eigenenergies1.size() * eigenenergies2.size());
124 : }
125 509 : const double inverse_sqrt_two = 1 / std::sqrt(2.0);
126 :
127 : // Construct the symmetry transformation by recording the contribution of the pair state
128 : // |idx1, idx2> to the symmetrized basis. Because the two atoms are identical (enforced above),
129 : // a one-atom state is uniquely identified across both atoms by its state index alone.
130 311021 : auto construct_symmetry_transformation = [&](Eigen::Index row_index, size_t idx1, size_t idx2) {
131 : // Following https://doi.org/10.1088/1361-6455/aa743a, pair states |a, a> cannot be of even
132 : // parity.
133 65456 : if (idx1 == idx2 &&
134 780 : (parity_under_inversion == Parity::EVEN || parity_under_permutation == Parity::EVEN)) {
135 33118 : return;
136 : }
137 :
138 : // Map the (unordered) pair of one-atom state indices to the column index of the symmetrized
139 : // state it contributes to, creating a new column the first time the pair is encountered.
140 65084 : std::array<size_t, 2> ordered_indices{std::max(idx1, idx2), std::min(idx1, idx2)};
141 65084 : auto [iterator, inserted] =
142 65084 : ket_indices2state_index.try_emplace(ordered_indices, state_index);
143 65084 : if (inserted) {
144 32746 : ++state_index;
145 : }
146 65084 : Eigen::Index column_index = iterator->second;
147 :
148 : // A pair state |a, a> contributes with coefficient one.
149 65084 : if (idx1 == idx2) {
150 408 : transformation_triplets.emplace_back(row_index, column_index, 1);
151 408 : return;
152 : }
153 :
154 : // We let pair states with idx1 > idx2 contribute with coefficient 1/sqrt(2) and put the
155 : // phase into the contribution of the partner state with idx1 < idx2.
156 64676 : if (idx1 > idx2) {
157 32338 : transformation_triplets.emplace_back(row_index, column_index, inverse_sqrt_two);
158 32338 : return;
159 : }
160 :
161 : // Determine the phase of the contribution of the partner state with idx1 < idx2.
162 32338 : int phase = 0;
163 32338 : if (parity_under_permutation != Parity::UNKNOWN) {
164 : // If both inversion and permutation are restricted, the earlier filter on the product
165 : // of parities already guarantees that the phases in the inversion- and
166 : // permutation-symmetric states are the same.
167 16252 : phase = -static_cast<int>(parity_under_permutation);
168 : } else {
169 : // If only inversion is restricted, the phase is determined by the product of the
170 : // parities of the one-atom states and the specified inversion parity.
171 48258 : phase = -static_cast<int>(parity_under_inversion) *
172 16086 : static_cast<int>(basis1->get_parity(idx1)) *
173 16086 : static_cast<int>(basis2->get_parity(idx2));
174 : }
175 32338 : transformation_triplets.emplace_back(row_index, column_index, phase * inverse_sqrt_two);
176 : };
177 :
178 : // Loop only over states with an allowed energy
179 509 : size_t ket_index = 0;
180 57644 : for (size_t idx1 = 0; idx1 < static_cast<size_t>(eigenenergies1.size()); ++idx1) {
181 57135 : set_task_status("Constructing pair basis...");
182 :
183 : // Get the energetically allowed range of the second index
184 57135 : size_t min = 0;
185 57135 : size_t max = eigenenergies2.size();
186 57135 : if (range_energy.is_finite()) {
187 39884 : real_t min_val2 = range_energy.min() - eigenenergies1[idx1];
188 39884 : real_t max_val2 = range_energy.max() - eigenenergies1[idx1];
189 39884 : min =
190 39884 : std::distance(eigenenergies2_begin,
191 : std::lower_bound(eigenenergies2_begin, eigenenergies2_end, min_val2));
192 39884 : max =
193 39884 : std::distance(eigenenergies2_begin,
194 : std::upper_bound(eigenenergies2_begin, eigenenergies2_end, max_val2));
195 : }
196 57135 : map_range_of_state_index2.try_emplace(idx1, typename basis_t::range_t(min, max));
197 :
198 : // Loop over the energetically allowed range of the second index
199 5988259 : for (size_t idx2 = min; idx2 < max; ++idx2) {
200 : // Get energy
201 2965594 : const real_t energy = eigenenergies1[idx1] + eigenenergies2[idx2];
202 2965594 : assert(!range_energy.is_finite() ||
203 : (energy >= range_energy.min() && energy <= range_energy.max()));
204 :
205 : // Check the parity of the product of the parities
206 2965594 : if (inferred_product_of_parities != Parity::UNKNOWN) {
207 144 : if (static_cast<int>(basis1->get_parity(idx1)) *
208 144 : static_cast<int>(basis2->get_parity(idx2)) !=
209 : static_cast<int>(inferred_product_of_parities)) {
210 67385 : continue;
211 : }
212 : }
213 :
214 : // Create a KetPair object
215 8896590 : auto ket = std::make_shared<ket_t>(
216 2965530 : typename ket_t::Private(), std::initializer_list<size_t>{idx1, idx2},
217 : std::initializer_list<std::shared_ptr<const BasisAtom<Scalar>>>{basis1, basis2},
218 : energy);
219 :
220 : // Check the quantum number m
221 2965530 : if (ket->has_quantum_number_m()) {
222 3092270 : if (range_quantum_number_m.is_finite() &&
223 126740 : (ket->get_quantum_number_m() <
224 126740 : range_quantum_number_m.min() - numerical_precision ||
225 73169 : ket->get_quantum_number_m() >
226 73169 : range_quantum_number_m.max() + numerical_precision)) {
227 67321 : continue;
228 : }
229 0 : } else if (range_quantum_number_m.is_finite()) {
230 0 : throw std::invalid_argument(
231 : "The quantum number m must not be restricted because it is not well-defined.");
232 : }
233 :
234 : // Store the KetPair object as a ket
235 2898209 : kets.emplace_back(std::move(ket));
236 2898209 : state_indices_to_ket_index.try_emplace(std::vector<size_t>{idx1, idx2}, ket_index);
237 :
238 2898209 : auto row_index = static_cast<Eigen::Index>(ket_index++);
239 2898209 : if (has_symmetry_restriction) {
240 65456 : construct_symmetry_transformation(row_index, idx1, idx2);
241 : }
242 : }
243 : }
244 :
245 509 : kets.shrink_to_fit();
246 :
247 509 : std::shared_ptr<const basis_t> basis = std::make_shared<basis_t>(
248 509 : typename basis_t::Private(), std::move(kets), std::move(map_range_of_state_index2),
249 509 : std::move(state_indices_to_ket_index), basis1, basis2);
250 :
251 509 : if (!has_symmetry_restriction) {
252 496 : return basis;
253 : }
254 :
255 13 : transformation_triplets.shrink_to_fit();
256 :
257 26 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor> transformation_matrix(
258 13 : basis->get_number_of_states(), state_index);
259 13 : transformation_matrix.setFromTriplets(transformation_triplets.begin(),
260 13 : transformation_triplets.end());
261 :
262 13 : Eigen::Matrix<real_t, Eigen::Dynamic, 1> sum_of_squared_coefficients =
263 : transformation_matrix.cwiseAbs2().transpose() *
264 : Eigen::Matrix<real_t, Eigen::Dynamic, 1>::Ones(transformation_matrix.rows());
265 32759 : for (Eigen::Index column_index = 0; column_index < sum_of_squared_coefficients.size();
266 : ++column_index) {
267 32746 : if (std::abs(sum_of_squared_coefficients[column_index] - 1) > numerical_precision) {
268 0 : throw std::invalid_argument(
269 : "The basis could not be symmetrized. This likely means that the specified parity "
270 : "restrictions are invalid for the given one-atom systems.");
271 : }
272 : }
273 :
274 : // TODO: on the long run, construct the coefficient matrix directly
275 13 : auto transformation = Transformation<Scalar>(std::move(transformation_matrix));
276 13 : return basis->transformed(transformation);
277 2966039 : }
278 :
279 : // Explicit instantiations
280 : template class BasisPairCreator<double>;
281 : template class BasisPairCreator<std::complex<double>>;
282 : } // namespace pairinteraction
|
