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