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