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/BasisAtomCreator.hpp"
8 : #include "pairinteraction/basis/BasisPair.hpp"
9 : #include "pairinteraction/database/Database.hpp"
10 : #include "pairinteraction/diagonalize/DiagonalizerEigen.hpp"
11 : #include "pairinteraction/enums/OperatorType.hpp"
12 : #include "pairinteraction/enums/SorterType.hpp"
13 : #include "pairinteraction/ket/KetAtom.hpp"
14 : #include "pairinteraction/ket/KetAtomCreator.hpp"
15 : #include "pairinteraction/ket/KetPair.hpp"
16 : #include "pairinteraction/system/SystemAtom.hpp"
17 : #include "pairinteraction/system/SystemPair.hpp"
18 : #include "pairinteraction/utils/hash.hpp"
19 : #include "pairinteraction/utils/streamed.hpp"
20 :
21 : #include <algorithm>
22 : #include <array>
23 : #include <cmath>
24 : #include <doctest/doctest.h>
25 : #include <optional>
26 : #include <unordered_map>
27 : #include <utility>
28 : #include <vector>
29 :
30 : namespace pairinteraction {
31 :
32 : constexpr double HARTREE_IN_GHZ = 6579683.920501762;
33 : constexpr double VOLT_PER_CM_IN_ATOMIC_UNITS = 1 / 5.14220675112e9;
34 : constexpr double UM_IN_ATOMIC_UNITS = 1 / 5.29177210544e-5;
35 :
36 : namespace {
37 : template <typename Scalar>
38 : std::shared_ptr<const typename BasisAtom<Scalar>::ket_t>
39 864 : get_dominant_ket(const std::shared_ptr<const BasisAtom<Scalar>> &basis, size_t state_index) {
40 864 : Eigen::MatrixX<Scalar> column(
41 864 : basis->get_coefficients().col(static_cast<Eigen::Index>(state_index)));
42 864 : Eigen::Index ket_index = 0;
43 864 : column.cwiseAbs().col(0).maxCoeff(&ket_index);
44 1728 : return basis->get_ket(static_cast<size_t>(ket_index));
45 864 : }
46 :
47 : template <typename Scalar>
48 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
49 2 : build_manual_symmetrizer(const std::shared_ptr<const BasisPair<Scalar>> &basis,
50 : std::optional<int> parity_under_inversion,
51 : std::optional<int> parity_under_permutation) {
52 : using real_t = typename BasisPair<Scalar>::real_t;
53 :
54 2 : const auto basis1 = basis->get_basis1();
55 2 : const auto basis2 = basis->get_basis2();
56 2 : const auto inv_sqrt_two = static_cast<real_t>(1 / std::sqrt(2.0));
57 :
58 2 : std::vector<Eigen::Triplet<Scalar>> triplets;
59 2 : triplets.reserve(2 * basis->get_number_of_states());
60 :
61 2 : Eigen::Index state_index = 0;
62 : std::unordered_map<std::array<size_t, 2>, Eigen::Index, utils::hash<std::array<size_t, 2>>>
63 2 : ket_ids_to_state_index;
64 26 : for (size_t idx1 = 0; idx1 < basis1->get_number_of_states(); ++idx1) {
65 312 : for (size_t idx2 = 0; idx2 < basis2->get_number_of_states(); ++idx2) {
66 288 : int ket_index = basis->get_ket_index_from_tuple(idx1, idx2);
67 288 : if (ket_index < 0) {
68 24 : continue;
69 : }
70 :
71 288 : size_t id1 = get_dominant_ket(basis1, idx1)->get_id_in_database();
72 288 : size_t id2 = get_dominant_ket(basis2, idx2)->get_id_in_database();
73 :
74 288 : if (id1 == id2) {
75 24 : if (parity_under_inversion == 1 || parity_under_permutation == 1) {
76 0 : continue;
77 : }
78 24 : triplets.emplace_back(ket_index, state_index++, Scalar{1});
79 24 : continue;
80 : }
81 :
82 264 : std::array<size_t, 2> ordered_ids{std::max(id1, id2), std::min(id1, id2)};
83 264 : auto [it, inserted] = ket_ids_to_state_index.try_emplace(ordered_ids, state_index);
84 264 : if (inserted) {
85 132 : ++state_index;
86 : }
87 264 : Eigen::Index column_index = it->second;
88 :
89 264 : if (id1 > id2) {
90 132 : triplets.emplace_back(ket_index, column_index, static_cast<Scalar>(inv_sqrt_two));
91 : } else {
92 132 : int swapped_sign = 0;
93 132 : if (parity_under_inversion.has_value() && !parity_under_permutation.has_value()) {
94 66 : swapped_sign = -*parity_under_inversion *
95 132 : static_cast<int>(basis1->get_quantum_number("parity", idx1)) *
96 66 : static_cast<int>(basis2->get_quantum_number("parity", idx2));
97 : } else {
98 66 : swapped_sign = -*parity_under_permutation;
99 : }
100 132 : triplets.emplace_back(ket_index, column_index,
101 132 : static_cast<Scalar>(swapped_sign * inv_sqrt_two));
102 : }
103 : }
104 : }
105 :
106 4 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor> transformation(
107 2 : static_cast<Eigen::Index>(basis->get_number_of_states()), state_index);
108 2 : transformation.setFromTriplets(triplets.begin(), triplets.end());
109 4 : return transformation;
110 2 : }
111 :
112 : template <typename Scalar>
113 : std::shared_ptr<const BasisPair<Scalar>>
114 4 : build_pair_basis(std::shared_ptr<const BasisAtom<Scalar>> basis1,
115 : std::shared_ptr<const BasisAtom<Scalar>> basis2) {
116 8 : return BasisPairCreator<Scalar>(std::move(basis1), std::move(basis2)).create();
117 : }
118 :
119 : template <typename Scalar>
120 2 : void check_same_pair_eigenenergies(const std::shared_ptr<const BasisPair<Scalar>> &basis1,
121 : const std::shared_ptr<const BasisPair<Scalar>> &basis2,
122 : const DiagonalizerEigen<Scalar> &diagonalizer) {
123 2 : auto system_pair_1 = SystemPair<Scalar>(basis1).set_distance_vector(
124 : std::array<typename BasisPair<Scalar>::real_t, 3>{0, 0, 1 * UM_IN_ATOMIC_UNITS});
125 2 : auto system_pair_2 = SystemPair<Scalar>(basis2).set_distance_vector(
126 : std::array<typename BasisPair<Scalar>::real_t, 3>{0, 0, 1 * UM_IN_ATOMIC_UNITS});
127 :
128 2 : system_pair_1.diagonalize(diagonalizer);
129 2 : system_pair_2.diagonalize(diagonalizer);
130 :
131 2 : auto eigenenergies_1 = system_pair_1.get_eigenenergies();
132 2 : auto eigenenergies_2 = system_pair_2.get_eigenenergies();
133 :
134 2 : DOCTEST_REQUIRE(eigenenergies_1.size() == eigenenergies_2.size());
135 2 : DOCTEST_CHECK(eigenenergies_1.isApprox(eigenenergies_2, 1e-11));
136 2 : }
137 : } // namespace
138 :
139 3 : DOCTEST_TEST_CASE("create a BasisPair") {
140 : // Create single-atom system
141 3 : Database &database = Database::get_global_instance();
142 3 : auto basis = BasisAtomCreator<double>(database)
143 6 : .set_species("Rb")
144 6 : .restrict_quantum_number("n", 58, 62)
145 6 : .restrict_quantum_number("l", 0, 2)
146 3 : .create();
147 3 : SystemAtom<double> system(basis);
148 3 : system.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
149 :
150 3 : DiagonalizerEigen<double> diagonalizer;
151 3 : system.diagonalize(diagonalizer);
152 :
153 : // Get energy window for a two-atom basis
154 3 : auto ket = KetAtomCreator(database)
155 6 : .set_species("Rb")
156 6 : .set_quantum_number("n", 60)
157 6 : .set_quantum_number("l", 0)
158 6 : .set_quantum_number("m", 0.5)
159 3 : .create();
160 3 : double min_energy = 2 * ket->get_energy() - 3 / HARTREE_IN_GHZ;
161 3 : double max_energy = 2 * ket->get_energy() + 3 / HARTREE_IN_GHZ;
162 :
163 : // Create two-atom bases
164 : auto basis_pair_a =
165 6 : pairinteraction::BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
166 3 : .restrict_energy(min_energy, max_energy)
167 3 : .restrict_quantum_number_m(1, 1)
168 3 : .create();
169 : auto basis_pair_b =
170 6 : pairinteraction::BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
171 3 : .restrict_energy(min_energy, max_energy)
172 3 : .restrict_quantum_number_m(1, 1)
173 3 : .create();
174 :
175 3 : DOCTEST_SUBCASE("check equality of kets") {
176 : // Obtain kets from the two-atom bases and check for equality
177 1 : auto ket1a = basis_pair_a->get_kets()[0];
178 1 : auto ket1b = basis_pair_b->get_kets()[0];
179 1 : auto ket2a = basis_pair_a->get_kets()[1];
180 1 : auto ket2b = basis_pair_b->get_kets()[1];
181 1 : DOCTEST_CHECK(*ket1a == *ket1a);
182 1 : DOCTEST_CHECK(*ket2a == *ket2a);
183 1 : DOCTEST_CHECK(*ket1a != *ket2b);
184 1 : DOCTEST_CHECK(*ket2a != *ket1b);
185 1 : DOCTEST_CHECK(*ket1a == *ket1b);
186 1 : DOCTEST_CHECK(*ket2a == *ket2b);
187 4 : }
188 :
189 3 : DOCTEST_SUBCASE("check overlap") {
190 1 : auto basis_ket = BasisAtomCreator<double>(database).add_ket(ket).create();
191 1 : auto basis_pair_ket = build_pair_basis<double>(basis_ket, basis_ket);
192 : Eigen::RowVectorXd amplitudes =
193 : basis_pair_a
194 2 : ->get_matrix_elements(basis_pair_ket, OperatorType::IDENTITY,
195 : OperatorType::IDENTITY, 0, 0)
196 1 : .row(0);
197 1 : auto overlaps = amplitudes.cwiseAbs2().eval();
198 :
199 : // The total overlap is less than 1 because of the restricted energy window
200 1 : DOCTEST_CHECK(overlaps.sum() == doctest::Approx(0.9107819201));
201 4 : }
202 :
203 3 : DOCTEST_SUBCASE("get the atomic states constituting a ket of the basis_pair") {
204 1 : auto atomic_states = basis_pair_a->get_kets()[0]->get_atomic_states();
205 1 : DOCTEST_CHECK(atomic_states.size() == 2);
206 1 : DOCTEST_CHECK(atomic_states[0]->get_number_of_states() == 1);
207 1 : DOCTEST_CHECK(atomic_states[0]->get_number_of_kets() == basis->get_number_of_kets());
208 4 : }
209 3 : }
210 :
211 1 : DOCTEST_TEST_CASE("merge rejects differently transformed atomic bases") {
212 1 : auto &database = Database::get_global_instance();
213 1 : auto atomic_basis = BasisAtomCreator<double>(database)
214 2 : .set_species("Rb")
215 2 : .restrict_quantum_number("n", 60, 60)
216 2 : .restrict_quantum_number("l", 0, 1)
217 2 : .restrict_quantum_number("m", 0.5, 0.5)
218 1 : .create();
219 :
220 1 : auto system_a = SystemAtom<double>(atomic_basis);
221 1 : system_a.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
222 1 : system_a.diagonalize(DiagonalizerEigen<double>());
223 :
224 1 : auto system_b = SystemAtom<double>(atomic_basis);
225 1 : system_b.set_electric_field({0, 0, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS});
226 1 : system_b.diagonalize(DiagonalizerEigen<double>());
227 :
228 : auto pair_basis_a =
229 1 : BasisPairCreator<double>(system_a.get_eigenbasis(), system_a.get_eigenbasis()).create();
230 : auto pair_basis_b =
231 1 : BasisPairCreator<double>(system_b.get_eigenbasis(), system_b.get_eigenbasis()).create();
232 :
233 2 : DOCTEST_CHECK_THROWS_WITH_AS(pair_basis_a->merge(pair_basis_b),
234 : doctest::Contains("Cannot merge two pair bases"),
235 : std::invalid_argument);
236 1 : }
237 :
238 2 : DOCTEST_TEST_CASE("get matrix elements in the pair basis") {
239 2 : DiagonalizerEigen<double> diagonalizer;
240 :
241 : // Create single-atom system
242 2 : Database &database = Database::get_global_instance();
243 2 : auto basis = BasisAtomCreator<double>(database)
244 4 : .set_species("Rb")
245 4 : .restrict_quantum_number("n", 58, 62)
246 4 : .restrict_quantum_number("l", 0, 2)
247 2 : .create();
248 2 : SystemAtom<double> system(basis);
249 2 : system.set_electric_field({0, 0, 10 * VOLT_PER_CM_IN_ATOMIC_UNITS});
250 2 : system.diagonalize(diagonalizer);
251 :
252 : // Get energy window for a two-atom basis
253 2 : auto ket = KetAtomCreator(database)
254 4 : .set_species("Rb")
255 4 : .set_quantum_number("n", 60)
256 4 : .set_quantum_number("l", 0)
257 4 : .set_quantum_number("m", 0.5)
258 2 : .create();
259 2 : double min_energy = 2 * ket->get_energy() - 3 / HARTREE_IN_GHZ;
260 2 : double max_energy = 2 * ket->get_energy() + 3 / HARTREE_IN_GHZ;
261 :
262 : // Create two-atom system
263 : auto basis_pair_unperturbed =
264 4 : pairinteraction::BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
265 2 : .restrict_energy(min_energy, max_energy)
266 2 : .restrict_quantum_number_m(1, 1)
267 2 : .create();
268 4 : auto system_pair = SystemPair<double>(basis_pair_unperturbed)
269 2 : .set_distance_vector({0, 0, 1 * UM_IN_ATOMIC_UNITS});
270 2 : system_pair.diagonalize(diagonalizer);
271 :
272 2 : auto basis_pair = system_pair.get_eigenbasis();
273 :
274 2 : DOCTEST_SUBCASE("check dimensions") {
275 : // <basis_pair_unperturbed|d0d0|basis_pair_unperturbed>
276 : auto matrix_elements_all = basis_pair_unperturbed->get_matrix_elements(
277 : basis_pair_unperturbed, OperatorType::ELECTRIC_DIPOLE, OperatorType::ELECTRIC_DIPOLE, 0,
278 1 : 0);
279 1 : DOCTEST_CHECK(matrix_elements_all.rows() == basis_pair_unperturbed->get_number_of_states());
280 1 : DOCTEST_CHECK(matrix_elements_all.cols() == basis_pair_unperturbed->get_number_of_states());
281 :
282 : // <ket_pair|d0d0|basis_pair_unperturbed>
283 1 : auto atomic_states = basis_pair_unperturbed->get_kets()[0]->get_atomic_states();
284 1 : auto basis_pair_ket_pair = build_pair_basis<double>(atomic_states[0], atomic_states[1]);
285 : Eigen::RowVectorXd matrix_elements_ket_pair =
286 : basis_pair_unperturbed
287 2 : ->get_matrix_elements(basis_pair_ket_pair, OperatorType::ELECTRIC_DIPOLE,
288 : OperatorType::ELECTRIC_DIPOLE, 0, 0)
289 1 : .row(0);
290 1 : DOCTEST_CHECK(matrix_elements_ket_pair.size() ==
291 : basis_pair_unperturbed->get_number_of_states());
292 :
293 : {
294 1 : Eigen::RowVectorXd ref = matrix_elements_all.row(0);
295 1 : DOCTEST_CHECK(ref.isApprox(matrix_elements_ket_pair, 1e-11));
296 1 : }
297 :
298 : // <basis x basis|d0d0|basis_pair>
299 1 : auto basis_pair_product = build_pair_basis<double>(basis, basis);
300 : auto matrix_elements_product = basis_pair->get_matrix_elements(
301 1 : basis_pair_product, OperatorType::ELECTRIC_DIPOLE, OperatorType::ELECTRIC_DIPOLE, 0, 0);
302 1 : DOCTEST_CHECK(matrix_elements_product.rows() ==
303 : basis->get_number_of_states() * basis->get_number_of_states());
304 1 : DOCTEST_CHECK(matrix_elements_product.cols() == basis_pair->get_number_of_states());
305 :
306 : // <ket,ket|d0d0|basis_pair>
307 1 : auto basis_ket = BasisAtomCreator<double>(database).add_ket(ket).create();
308 1 : auto basis_pair_ket = build_pair_basis<double>(basis_ket, basis_ket);
309 : auto matrix_elements_ket = basis_pair->get_matrix_elements(
310 1 : basis_pair_ket, OperatorType::ELECTRIC_DIPOLE, OperatorType::ELECTRIC_DIPOLE, 0, 0);
311 1 : DOCTEST_CHECK(matrix_elements_ket.rows() == 1);
312 1 : DOCTEST_CHECK(matrix_elements_ket.cols() == basis_pair->get_number_of_states());
313 3 : }
314 :
315 2 : DOCTEST_SUBCASE("check matrix elements") {
316 : // energy
317 : auto hamiltonian = basis_pair->get_matrix_elements(basis_pair, OperatorType::ENERGY,
318 1 : OperatorType::IDENTITY, 0, 0);
319 2 : hamiltonian += basis_pair->get_matrix_elements(basis_pair, OperatorType::IDENTITY,
320 1 : OperatorType::ENERGY, 0, 0);
321 :
322 : // interaction with electric field
323 : {
324 2 : Eigen::SparseMatrix<double, Eigen::RowMajor> tmp = -basis_pair->get_matrix_elements(
325 1 : basis_pair, OperatorType::ELECTRIC_DIPOLE, OperatorType::IDENTITY, 0, 0);
326 2 : tmp += -basis_pair->get_matrix_elements(basis_pair, OperatorType::IDENTITY,
327 1 : OperatorType::ELECTRIC_DIPOLE, 0, 0);
328 1 : hamiltonian += 10 * VOLT_PER_CM_IN_ATOMIC_UNITS * tmp;
329 1 : }
330 :
331 : // dipole-dipole interaction
332 : {
333 1 : Eigen::SparseMatrix<double, Eigen::RowMajor> tmp = -2 *
334 2 : basis_pair->get_matrix_elements(basis_pair, OperatorType::ELECTRIC_DIPOLE,
335 1 : OperatorType::ELECTRIC_DIPOLE, 0, 0);
336 2 : tmp += -basis_pair->get_matrix_elements(basis_pair, OperatorType::ELECTRIC_DIPOLE,
337 1 : OperatorType::ELECTRIC_DIPOLE, 1, -1);
338 2 : tmp += -basis_pair->get_matrix_elements(basis_pair, OperatorType::ELECTRIC_DIPOLE,
339 1 : OperatorType::ELECTRIC_DIPOLE, -1, 1);
340 1 : hamiltonian += std::pow(UM_IN_ATOMIC_UNITS, -3) * tmp;
341 1 : }
342 :
343 : // compare to reference
344 1 : const auto &ref = system_pair.get_matrix();
345 1 : DOCTEST_CHECK(ref.isApprox(hamiltonian, 1e-11));
346 3 : }
347 2 : }
348 :
349 1 : DOCTEST_TEST_CASE("get amplitudes (via matrix elements) between different pair basis") {
350 1 : DiagonalizerEigen<double> diagonalizer;
351 :
352 1 : Database &database = Database::get_global_instance();
353 1 : auto atomic_basis = BasisAtomCreator<double>(database)
354 2 : .set_species("Rb")
355 2 : .restrict_quantum_number("n", 58, 62)
356 2 : .restrict_quantum_number("l", 0, 2)
357 2 : .restrict_quantum_number("m", 0.5, 0.5)
358 1 : .create();
359 :
360 1 : auto ket = KetAtomCreator(database)
361 2 : .set_species("Rb")
362 2 : .set_quantum_number("n", 60)
363 2 : .set_quantum_number("l", 0)
364 2 : .set_quantum_number("m", 0.5)
365 1 : .create();
366 :
367 1 : auto perturbed_system1 = SystemAtom<double>(atomic_basis);
368 1 : perturbed_system1.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
369 1 : perturbed_system1.diagonalize(diagonalizer);
370 :
371 1 : auto perturbed_system2 = SystemAtom<double>(atomic_basis);
372 1 : perturbed_system2.set_electric_field({0, 0, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS});
373 1 : perturbed_system2.diagonalize(diagonalizer);
374 :
375 1 : double min_energy = 2 * ket->get_energy() - 20 / HARTREE_IN_GHZ;
376 1 : double max_energy = 2 * ket->get_energy() + 20 / HARTREE_IN_GHZ;
377 2 : auto perturbed_basis = BasisPairCreator<double>(perturbed_system1.get_eigenbasis(),
378 2 : perturbed_system2.get_eigenbasis())
379 1 : .restrict_energy(min_energy, max_energy)
380 1 : .create();
381 1 : auto state_in_perturbed_basis = perturbed_basis->get_state(42);
382 :
383 1 : min_energy = 2 * ket->get_energy() - 10 / HARTREE_IN_GHZ;
384 1 : max_energy = 2 * ket->get_energy() + 10 / HARTREE_IN_GHZ;
385 2 : auto small_unperturbed_basis = BasisPairCreator<double>(atomic_basis, atomic_basis)
386 1 : .restrict_energy(min_energy, max_energy)
387 1 : .create();
388 1 : DOCTEST_CHECK(perturbed_basis->get_number_of_states() !=
389 : small_unperturbed_basis->get_number_of_states());
390 :
391 : auto amplitudes = state_in_perturbed_basis->get_matrix_elements(
392 1 : small_unperturbed_basis, OperatorType::IDENTITY, OperatorType::IDENTITY, 0, 0);
393 1 : DOCTEST_CHECK(amplitudes.rows() == small_unperturbed_basis->get_number_of_states());
394 1 : DOCTEST_CHECK(amplitudes.cols() == state_in_perturbed_basis->get_number_of_states());
395 1 : }
396 :
397 6 : DOCTEST_TEST_CASE("create a symmetrized BasisPair") {
398 6 : auto &database = Database::get_global_instance();
399 6 : auto diagonalizer = DiagonalizerEigen<double>();
400 :
401 6 : auto basis = BasisAtomCreator<double>(database)
402 12 : .set_species("Rb")
403 12 : .restrict_quantum_number("n", 60, 61)
404 12 : .restrict_quantum_number("l", 0, 1)
405 12 : .restrict_quantum_number("m", -0.5, 0.5)
406 6 : .create();
407 :
408 6 : SystemAtom<double> system(basis);
409 6 : system.diagonalize(diagonalizer);
410 :
411 : auto canonical_basis =
412 6 : BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis()).create();
413 :
414 6 : DOCTEST_SUBCASE("restrict permutation parity") {
415 : auto symmetrized_basis =
416 2 : BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
417 1 : .restrict_parity_under_permutation(-1)
418 1 : .create();
419 :
420 1 : auto expected_basis = canonical_basis->transformed(
421 1 : build_manual_symmetrizer(canonical_basis, std::nullopt, -1));
422 :
423 1 : check_same_pair_eigenenergies(symmetrized_basis, expected_basis, diagonalizer);
424 1 : DOCTEST_CHECK(symmetrized_basis->get_number_of_states() <
425 : canonical_basis->get_number_of_states());
426 7 : }
427 :
428 6 : DOCTEST_SUBCASE("restrict inversion parity") {
429 : auto symmetrized_basis =
430 2 : BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
431 1 : .restrict_parity_under_inversion(-1)
432 1 : .create();
433 :
434 1 : auto expected_basis = canonical_basis->transformed(
435 1 : build_manual_symmetrizer(canonical_basis, -1, std::nullopt));
436 :
437 1 : check_same_pair_eigenenergies(symmetrized_basis, expected_basis, diagonalizer);
438 1 : DOCTEST_CHECK(symmetrized_basis->get_number_of_states() <
439 : canonical_basis->get_number_of_states());
440 7 : }
441 :
442 6 : DOCTEST_SUBCASE("restrict permutation parity to EVEN includes identical-state kets") {
443 : auto symmetrized_basis_even =
444 2 : BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
445 1 : .restrict_parity_under_permutation(1)
446 1 : .create();
447 :
448 : auto symmetrized_basis_odd =
449 2 : BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
450 1 : .restrict_parity_under_permutation(-1)
451 1 : .create();
452 :
453 : // Count kets in the canonical basis where both atoms are in the same state (id1 == id2).
454 : // Such kets are always permutation-symmetric and must appear in EVEN but not ODD.
455 1 : size_t num_diagonal_kets = 0;
456 145 : for (const auto &ket : *canonical_basis) {
457 144 : auto atomic_states = ket->get_atomic_states();
458 144 : DOCTEST_REQUIRE(atomic_states.size() == 2);
459 144 : size_t id1 = get_dominant_ket(atomic_states[0], 0)->get_id_in_database();
460 144 : size_t id2 = get_dominant_ket(atomic_states[1], 0)->get_id_in_database();
461 144 : if (id1 == id2) {
462 12 : ++num_diagonal_kets;
463 : }
464 144 : }
465 1 : DOCTEST_REQUIRE(num_diagonal_kets > 0);
466 :
467 : // EVEN and ODD together partition all canonical kets: off-diagonal pairs contribute one
468 : // state to each sector, diagonal kets only contribute to EVEN.
469 1 : DOCTEST_CHECK(symmetrized_basis_even->get_number_of_states() +
470 : symmetrized_basis_odd->get_number_of_states() ==
471 : canonical_basis->get_number_of_states());
472 1 : DOCTEST_CHECK(symmetrized_basis_odd->get_number_of_states() -
473 : symmetrized_basis_even->get_number_of_states() ==
474 : num_diagonal_kets);
475 :
476 : // The eigenenergies of the EVEN and ODD bases together must reproduce the
477 : // eigenenergies of the canonical (non-symmetrized) basis.
478 : auto system_pair_canonical =
479 1 : SystemPair<double>(canonical_basis).set_distance_vector({0, 0, 1 * UM_IN_ATOMIC_UNITS});
480 2 : auto system_pair_even = SystemPair<double>(symmetrized_basis_even)
481 1 : .set_distance_vector({0, 0, 1 * UM_IN_ATOMIC_UNITS});
482 2 : auto system_pair_odd = SystemPair<double>(symmetrized_basis_odd)
483 1 : .set_distance_vector({0, 0, 1 * UM_IN_ATOMIC_UNITS});
484 :
485 1 : system_pair_canonical.diagonalize(diagonalizer);
486 1 : system_pair_even.diagonalize(diagonalizer);
487 1 : system_pair_odd.diagonalize(diagonalizer);
488 :
489 1 : auto canonical_eigenenergies = system_pair_canonical.get_eigenenergies();
490 1 : auto even_eigenenergies = system_pair_even.get_eigenenergies();
491 1 : auto odd_eigenenergies = system_pair_odd.get_eigenenergies();
492 :
493 1 : DOCTEST_REQUIRE(canonical_eigenenergies.size() ==
494 : even_eigenenergies.size() + odd_eigenenergies.size());
495 :
496 1 : Eigen::VectorXd combined_eigenenergies(canonical_eigenenergies.size());
497 1 : combined_eigenenergies << even_eigenenergies, odd_eigenenergies;
498 1 : std::sort(combined_eigenenergies.data(),
499 1 : combined_eigenenergies.data() + combined_eigenenergies.size());
500 1 : std::sort(canonical_eigenenergies.data(),
501 1 : canonical_eigenenergies.data() + canonical_eigenenergies.size());
502 :
503 1 : DOCTEST_CHECK(combined_eigenenergies.isApprox(canonical_eigenenergies, 1e-11));
504 7 : }
505 :
506 6 : DOCTEST_SUBCASE("combine inversion and permutation parity") {
507 : auto symmetrized_basis =
508 2 : BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
509 1 : .restrict_parity_under_inversion(-1)
510 1 : .restrict_parity_under_permutation(-1)
511 1 : .create();
512 :
513 1 : DOCTEST_CHECK(symmetrized_basis->get_number_of_states() <
514 : canonical_basis->get_number_of_states());
515 :
516 : Eigen::SparseMatrix<double, Eigen::ColMajor> coefficients =
517 1 : symmetrized_basis->get_coefficients();
518 1 : const double inv_sqrt_two = 1 / std::sqrt(2.0);
519 :
520 47 : for (int state_index = 0; state_index < coefficients.outerSize(); ++state_index) {
521 46 : std::vector<std::pair<int, double>> entries;
522 46 : for (Eigen::SparseMatrix<double, Eigen::ColMajor>::InnerIterator it(coefficients,
523 46 : state_index);
524 126 : it; ++it) {
525 80 : auto atomic_states = symmetrized_basis->get_kets()[it.row()]->get_atomic_states();
526 80 : DOCTEST_REQUIRE(atomic_states.size() == 2);
527 80 : DOCTEST_CHECK(
528 : static_cast<int>(atomic_states[0]->get_quantum_number("parity", 0)) *
529 : static_cast<int>(atomic_states[1]->get_quantum_number("parity", 0)) ==
530 : 1);
531 80 : entries.emplace_back(it.row(), it.value());
532 80 : }
533 :
534 46 : DOCTEST_CHECK(entries.size() >= 1);
535 46 : DOCTEST_CHECK(entries.size() <= 2);
536 46 : if (entries.size() == 1) {
537 12 : DOCTEST_CHECK(entries[0].second == doctest::Approx(1));
538 : } else {
539 34 : DOCTEST_CHECK(std::abs(entries[0].second) == doctest::Approx(inv_sqrt_two));
540 34 : DOCTEST_CHECK(std::abs(entries[1].second) == doctest::Approx(inv_sqrt_two));
541 34 : DOCTEST_CHECK(entries[0].second == doctest::Approx(entries[1].second));
542 : }
543 46 : }
544 7 : }
545 :
546 6 : DOCTEST_SUBCASE("parity restrictions require the same SystemAtom twice") {
547 : // A second, independently constructed system represents a different atom. Even though it
548 : // is built from the same basis, it is a distinct object, so symmetrization is rejected.
549 1 : SystemAtom<double> system_other(basis);
550 1 : system_other.diagonalize(diagonalizer);
551 :
552 4 : DOCTEST_CHECK_THROWS_AS(
553 : BasisPairCreator<double>(system.get_eigenbasis(), system_other.get_eigenbasis())
554 : .restrict_parity_under_permutation(-1)
555 : .create(),
556 : std::invalid_argument);
557 :
558 4 : DOCTEST_CHECK_THROWS_AS(
559 : BasisPairCreator<double>(system.get_eigenbasis(), system_other.get_eigenbasis())
560 : .restrict_parity_under_inversion(-1)
561 : .create(),
562 : std::invalid_argument);
563 :
564 : // Without a parity restriction, two different systems remain allowed.
565 1 : DOCTEST_CHECK_NOTHROW(
566 : BasisPairCreator<double>(system.get_eigenbasis(), system_other.get_eigenbasis())
567 : .create());
568 7 : }
569 :
570 6 : DOCTEST_SUBCASE("parity restrictions must be +1 or -1") {
571 1 : auto eigenbasis = system.get_eigenbasis();
572 4 : DOCTEST_CHECK_THROWS_AS(
573 : BasisPairCreator<double>(eigenbasis, eigenbasis).restrict_parity_under_inversion(0),
574 : std::invalid_argument);
575 4 : DOCTEST_CHECK_THROWS_AS(
576 : BasisPairCreator<double>(eigenbasis, eigenbasis).restrict_parity_under_permutation(2),
577 : std::invalid_argument);
578 7 : }
579 6 : }
580 :
581 4 : DOCTEST_TEST_CASE("BasisPairCreator checks the energies of the BasisAtom") {
582 4 : auto &database = Database::get_global_instance();
583 4 : auto basis = BasisAtomCreator<double>(database)
584 8 : .set_species("Rb")
585 8 : .restrict_quantum_number("n", 60, 61)
586 8 : .restrict_quantum_number("l", 0, 1)
587 4 : .create();
588 4 : double min_energy = 2 * basis->get_energy(0);
589 4 : double max_energy = 2 * basis->get_energy(basis->get_number_of_states() - 1);
590 :
591 4 : DOCTEST_SUBCASE("a canonical basis sorted by energy can be used") {
592 1 : DOCTEST_CHECK_NOTHROW(BasisPairCreator<double>(basis, basis).create());
593 1 : DOCTEST_CHECK_NOTHROW(BasisPairCreator<double>(basis, basis)
594 : .restrict_energy(min_energy, max_energy)
595 : .create());
596 4 : }
597 :
598 4 : DOCTEST_SUBCASE("the states must have well-defined energies for an energy restriction") {
599 1 : auto basis_without_energies = basis->copy_with_coefficients(basis->get_coefficients());
600 4 : DOCTEST_CHECK_THROWS_WITH_AS(BasisPairCreator<double>(basis, basis_without_energies)
601 : .restrict_energy(min_energy, max_energy)
602 : .create(),
603 : doctest::Contains("well-defined energies"),
604 : std::invalid_argument);
605 :
606 : // Without an energy restriction, the pair states do not have a well-defined energy
607 1 : auto pair_basis = BasisPairCreator<double>(basis, basis_without_energies).create();
608 1 : DOCTEST_CHECK(pair_basis->get_number_of_states() ==
609 : basis->get_number_of_states() * basis->get_number_of_states());
610 1 : DOCTEST_CHECK_FALSE(pair_basis->has_energy());
611 5 : }
612 :
613 4 : DOCTEST_SUBCASE("the states must be sorted by energy for an energy restriction") {
614 : auto basis_sorted_by_m =
615 1 : basis->transformed(basis->get_sorter({SorterType::QUANTUM_NUMBER_M}));
616 4 : DOCTEST_CHECK_THROWS_WITH_AS(BasisPairCreator<double>(basis_sorted_by_m, basis)
617 : .restrict_energy(min_energy, max_energy)
618 : .create(),
619 : doctest::Contains("sorted by energy"), std::invalid_argument);
620 1 : DOCTEST_CHECK_NOTHROW(BasisPairCreator<double>(basis_sorted_by_m, basis).create());
621 5 : }
622 :
623 4 : DOCTEST_SUBCASE("the eigenbasis of a diagonalized system can be used") {
624 1 : SystemAtom<double> system(basis);
625 1 : system.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
626 1 : system.diagonalize(DiagonalizerEigen<double>());
627 1 : auto eigenbasis = system.get_eigenbasis();
628 1 : auto pair_basis = BasisPairCreator<double>(eigenbasis, eigenbasis).create();
629 :
630 : // The energies of the pair states are the sums of the eigenenergies of the atoms
631 257 : for (const auto &ket : *pair_basis) {
632 256 : auto atomic_states = ket->get_atomic_states();
633 256 : DOCTEST_CHECK(
634 : ket->get_energy() ==
635 : doctest::Approx(atomic_states[0]->get_energy(0) + atomic_states[1]->get_energy(0)));
636 256 : }
637 :
638 : // The eigenbasis of a system that was not sorted by energy cannot be used with an energy
639 : // restriction
640 1 : SystemAtom<double> unsorted_system(basis);
641 1 : unsorted_system.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
642 1 : unsorted_system.diagonalize(DiagonalizerEigen<double>(), {}, {}, 1e-6, false);
643 4 : DOCTEST_CHECK_THROWS_WITH_AS(
644 : BasisPairCreator<double>(unsorted_system.get_eigenbasis(), eigenbasis)
645 : .restrict_energy(min_energy, max_energy)
646 : .create(),
647 : doctest::Contains("sorted by energy"), std::invalid_argument);
648 5 : }
649 4 : }
650 :
651 : } // namespace pairinteraction
|
