Line data Source code
1 : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/basis/BasisAtomCreator.hpp"
5 :
6 : #include "pairinteraction/basis/BasisAtom.hpp"
7 : #include "pairinteraction/database/Database.hpp"
8 : #include "pairinteraction/diagonalize/DiagonalizerEigen.hpp"
9 : #include "pairinteraction/enums/OperatorType.hpp"
10 : #include "pairinteraction/enums/Parity.hpp"
11 : #include "pairinteraction/enums/TransformationType.hpp"
12 : #include "pairinteraction/ket/KetAtom.hpp"
13 : #include "pairinteraction/ket/KetAtomCreator.hpp"
14 : #include "pairinteraction/system/SystemAtom.hpp"
15 :
16 : #include <doctest/doctest.h>
17 :
18 : namespace pairinteraction {
19 :
20 : constexpr double VOLT_PER_CM_IN_ATOMIC_UNITS = 1 / 5.14220675112e9;
21 :
22 1 : DOCTEST_TEST_CASE("create a basis for strontium 88") {
23 1 : Database &database = Database::get_global_instance();
24 1 : auto basis = BasisAtomCreator<double>()
25 2 : .set_species("Sr88_singlet")
26 1 : .restrict_quantum_number_n(60, 60)
27 1 : .restrict_quantum_number_l(0, 2)
28 1 : .create(database);
29 19 : for (const auto &ket : *basis) {
30 9 : DOCTEST_CHECK(ket->get_species() == "Sr88_singlet");
31 9 : DOCTEST_MESSAGE("Ket: ", *ket);
32 9 : }
33 1 : }
34 :
35 1 : DOCTEST_TEST_CASE("create a basis for strontium 87") {
36 1 : Database &database = Database::get_global_instance();
37 1 : auto basis = BasisAtomCreator<double>()
38 2 : .set_species("Sr87_mqdt")
39 1 : .restrict_quantum_number_nu(59, 61)
40 1 : .restrict_quantum_number_l(0, 0)
41 1 : .create(database);
42 113 : for (const auto &ket : *basis) {
43 56 : DOCTEST_CHECK(ket->get_species() == "Sr87_mqdt");
44 56 : DOCTEST_MESSAGE("Ket: ", *ket);
45 56 : }
46 1 : }
47 :
48 1 : DOCTEST_TEST_CASE("create a basis from kets") {
49 1 : Database &database = Database::get_global_instance();
50 1 : auto ket1 = KetAtomCreator("Sr88_singlet", 59, 0, 0, 0).create(database);
51 1 : auto ket2 = KetAtomCreator("Sr88_singlet", 60, 0, 0, 0).create(database);
52 1 : auto ket3 = KetAtomCreator("Sr88_singlet", 61, 0, 0, 0).create(database);
53 : auto basis =
54 2 : BasisAtomCreator<double>().append_ket(ket1).append_ket(ket2).append_ket(ket3).create(
55 1 : database);
56 7 : for (const auto &ket : *basis) {
57 3 : DOCTEST_CHECK(ket->get_species() == "Sr88_singlet");
58 3 : DOCTEST_MESSAGE("Ket: ", *ket);
59 3 : }
60 1 : }
61 :
62 1 : DOCTEST_TEST_CASE("create a basis and sort it according to parity and m") {
63 1 : Database &database = Database::get_global_instance();
64 1 : auto basis_unsorted = BasisAtomCreator<double>()
65 2 : .set_species("Rb")
66 1 : .restrict_quantum_number_n(60, 60)
67 1 : .restrict_quantum_number_l(0, 3)
68 1 : .restrict_quantum_number_m(-0.5, 0.5)
69 1 : .create(database);
70 :
71 : // Sort the basis by parity and the m quantum number
72 1 : auto sorter = basis_unsorted->get_sorter(
73 1 : {TransformationType::SORT_BY_PARITY, TransformationType::SORT_BY_QUANTUM_NUMBER_M});
74 1 : auto basis = basis_unsorted->transformed(sorter);
75 :
76 : // Check if the basis is properly sorted
77 1 : auto parity = Parity::ODD;
78 1 : auto quantum_number_m = std::numeric_limits<double>::lowest();
79 15 : for (size_t i = 0; i < basis->get_number_of_states(); ++i) {
80 14 : DOCTEST_MESSAGE("State ", i, ": Parity = ", basis->get_parity(i),
81 : ", M = ", basis->get_quantum_number_m(i));
82 14 : DOCTEST_CHECK(basis->get_parity(i) >= parity);
83 14 : if (basis->get_parity(i) != parity) {
84 1 : parity = basis->get_parity(i);
85 1 : quantum_number_m = std::numeric_limits<double>::lowest();
86 : }
87 14 : DOCTEST_CHECK(basis->get_quantum_number_m(i) >= quantum_number_m);
88 14 : quantum_number_m = basis->get_quantum_number_m(i);
89 : }
90 :
91 : // Check that the blocks are correctly determined
92 1 : auto blocks = basis->get_indices_of_blocks(
93 1 : {TransformationType::SORT_BY_PARITY, TransformationType::SORT_BY_QUANTUM_NUMBER_M});
94 1 : std::vector<size_t> expected_start = {0, 4, 8, 11};
95 :
96 1 : DOCTEST_CHECK(blocks.size() == expected_start.size());
97 :
98 1 : size_t idx = 0;
99 5 : for (const auto &block : blocks) {
100 4 : DOCTEST_MESSAGE("Block ", idx, " starts at ", block.start);
101 4 : DOCTEST_CHECK(block.start == expected_start[idx]);
102 4 : idx++;
103 : }
104 :
105 : // Test implicit conversion of an eigen matrix to a transformator
106 1 : size_t dim = basis->get_number_of_states();
107 : Eigen::SparseMatrix<double, Eigen::RowMajor> matrix(static_cast<long>(dim),
108 1 : static_cast<long>(dim));
109 1 : matrix.setIdentity();
110 1 : auto transformed = basis->transformed(matrix);
111 1 : auto transformation = transformed->get_transformation();
112 1 : DOCTEST_CHECK(transformation.transformation_type.back() == TransformationType::ARBITRARY);
113 1 : }
114 :
115 3 : DOCTEST_TEST_CASE("calculation of matrix elements") {
116 3 : auto &database = Database::get_global_instance();
117 :
118 3 : auto ket_s = KetAtomCreator()
119 6 : .set_species("Rb")
120 3 : .set_quantum_number_n(60)
121 3 : .set_quantum_number_l(0)
122 3 : .set_quantum_number_j(0.5)
123 3 : .set_quantum_number_m(0.5)
124 3 : .create(database);
125 :
126 3 : auto ket_p = KetAtomCreator()
127 6 : .set_species("Rb")
128 3 : .set_quantum_number_n(60)
129 3 : .set_quantum_number_l(1)
130 3 : .set_quantum_number_j(0.5)
131 3 : .set_quantum_number_m(0.5)
132 3 : .create(database);
133 :
134 3 : auto basis = BasisAtomCreator<double>()
135 6 : .set_species("Rb")
136 3 : .restrict_quantum_number_n(59, 61)
137 3 : .restrict_quantum_number_l(0, 1)
138 3 : .restrict_quantum_number_m(0.5, 0.5)
139 3 : .create(database);
140 :
141 3 : SystemAtom<double> system(basis);
142 :
143 3 : DOCTEST_SUBCASE("calculate energy") {
144 2 : auto m1 = basis->get_canonical_state_from_ket(ket_s)->get_matrix_elements(
145 1 : ket_s, OperatorType::ENERGY, 0);
146 1 : DOCTEST_CHECK(m1.size() == 1);
147 1 : double energy1 = m1[0];
148 :
149 1 : auto m2 = basis->get_matrix_elements(ket_s, OperatorType::ENERGY, 0);
150 1 : DOCTEST_CHECK(m2.size() == basis->get_number_of_states());
151 1 : double energy2 = m2[static_cast<int>(basis->get_corresponding_state_index(ket_s))];
152 :
153 1 : double reference = ket_s->get_energy();
154 1 : DOCTEST_CHECK(std::abs(energy1 - reference) < 1e-11);
155 1 : DOCTEST_CHECK(std::abs(energy2 - reference) < 1e-11);
156 4 : }
157 :
158 3 : DOCTEST_SUBCASE("calculate electric dipole matrix element") {
159 1 : auto m = basis->get_matrix_elements(ket_p, OperatorType::ELECTRIC_DIPOLE, 0);
160 1 : DOCTEST_CHECK(m.size() == basis->get_number_of_states());
161 1 : double dipole = m[static_cast<int>(basis->get_corresponding_state_index(ket_s))];
162 :
163 1 : DOCTEST_CHECK(std::abs(dipole - 1247.5985544327215848) < 1e-6);
164 4 : }
165 :
166 3 : DOCTEST_SUBCASE("calculate electric dipole matrix element with and without an induced dipole") {
167 : {
168 1 : auto state = basis->get_corresponding_state(ket_s);
169 :
170 1 : auto m = state->get_matrix_elements(state, OperatorType::ELECTRIC_DIPOLE, 0);
171 1 : DOCTEST_CHECK(m.rows() == 1);
172 1 : DOCTEST_CHECK(m.cols() == 1);
173 1 : double dipole = m.coeff(0, 0);
174 :
175 1 : DOCTEST_CHECK(std::abs(dipole - 0) < 1e-6);
176 1 : }
177 :
178 : {
179 1 : system.set_electric_field({0, 0, VOLT_PER_CM_IN_ATOMIC_UNITS});
180 1 : system.diagonalize(DiagonalizerEigen<double>());
181 1 : auto state = system.get_eigenbasis()->get_corresponding_state(ket_s);
182 :
183 1 : auto m = state->get_matrix_elements(state, OperatorType::ELECTRIC_DIPOLE, 0);
184 1 : DOCTEST_CHECK(m.rows() == 1);
185 1 : DOCTEST_CHECK(m.cols() == 1);
186 1 : double dipole = m.coeff(0, 0);
187 :
188 1 : DOCTEST_CHECK(std::abs(dipole - 135.01903532588247003) < 1e-6);
189 1 : }
190 3 : }
191 3 : }
192 :
193 : } // namespace pairinteraction
|
