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 2 : .restrict_quantum_number("n", 60, 60)
27 2 : .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 : }
32 1 : }
33 :
34 1 : DOCTEST_TEST_CASE("create a basis for strontium 87") {
35 1 : Database &database = Database::get_global_instance();
36 1 : auto basis = BasisAtomCreator<double>()
37 2 : .set_species("Sr87_mqdt")
38 2 : .restrict_quantum_number("nu", 59, 61)
39 2 : .restrict_quantum_number("l", 0, 0)
40 1 : .create(database);
41 113 : for (const auto &ket : *basis) {
42 56 : DOCTEST_CHECK(ket->get_species() == "Sr87_mqdt");
43 56 : }
44 1 : }
45 :
46 1 : DOCTEST_TEST_CASE("create a basis from kets") {
47 1 : Database &database = Database::get_global_instance();
48 1 : auto ket1 = KetAtomCreator("Sr88_singlet", 59, 0, 0, 0).create(database);
49 1 : auto ket2 = KetAtomCreator("Sr88_singlet", 60, 0, 0, 0).create(database);
50 1 : auto ket3 = KetAtomCreator("Sr88_singlet", 61, 0, 0, 0).create(database);
51 : auto basis =
52 1 : BasisAtomCreator<double>().add_ket(ket1).add_ket(ket2).add_ket(ket3).create(database);
53 7 : for (const auto &ket : *basis) {
54 3 : DOCTEST_CHECK(ket->get_species() == "Sr88_singlet");
55 3 : }
56 1 : }
57 :
58 1 : DOCTEST_TEST_CASE("create a basis and sort it according to parity and m") {
59 1 : Database &database = Database::get_global_instance();
60 1 : auto basis_unsorted = BasisAtomCreator<double>()
61 2 : .set_species("Rb")
62 2 : .restrict_quantum_number("n", 60, 60)
63 2 : .restrict_quantum_number("l", 0, 3)
64 2 : .restrict_quantum_number("m", -0.5, 0.5)
65 1 : .create(database);
66 :
67 : // Sort the basis by parity and the m quantum number
68 1 : auto sorter = basis_unsorted->get_sorter(
69 1 : {TransformationType::SORT_BY_PARITY, TransformationType::SORT_BY_QUANTUM_NUMBER_M});
70 1 : auto basis = basis_unsorted->transformed(sorter);
71 :
72 : // Check if the basis is properly sorted
73 1 : auto parity = Parity::ODD;
74 1 : auto quantum_number_m = std::numeric_limits<double>::lowest();
75 15 : for (size_t i = 0; i < basis->get_number_of_states(); ++i) {
76 14 : DOCTEST_MESSAGE("State ", i, ": Parity = ", basis->get_parity(i),
77 : ", M = ", basis->get_quantum_number_m(i));
78 14 : DOCTEST_CHECK(basis->get_parity(i) >= parity);
79 14 : if (basis->get_parity(i) != parity) {
80 1 : parity = basis->get_parity(i);
81 1 : quantum_number_m = std::numeric_limits<double>::lowest();
82 : }
83 14 : DOCTEST_CHECK(basis->get_quantum_number_m(i) >= quantum_number_m);
84 14 : quantum_number_m = basis->get_quantum_number_m(i);
85 : }
86 :
87 : // Check that the blocks are correctly determined
88 1 : auto blocks = basis->get_indices_of_blocks(
89 1 : {TransformationType::SORT_BY_PARITY, TransformationType::SORT_BY_QUANTUM_NUMBER_M});
90 1 : std::vector<size_t> expected_start = {0, 4, 8, 11};
91 :
92 1 : DOCTEST_CHECK(blocks.size() == expected_start.size());
93 :
94 1 : size_t idx = 0;
95 5 : for (const auto &block : blocks) {
96 4 : DOCTEST_MESSAGE("Block ", idx, " starts at ", block.start);
97 4 : DOCTEST_CHECK(block.start == expected_start[idx]);
98 4 : idx++;
99 : }
100 :
101 : // Test implicit conversion of an eigen matrix to a transformator
102 1 : size_t dim = basis->get_number_of_states();
103 : Eigen::SparseMatrix<double, Eigen::RowMajor> matrix(static_cast<long>(dim),
104 1 : static_cast<long>(dim));
105 1 : matrix.setIdentity();
106 1 : auto transformed = basis->transformed(matrix);
107 1 : auto transformation = transformed->get_transformation();
108 1 : DOCTEST_CHECK(transformation.transformation_type.back() == TransformationType::ARBITRARY);
109 1 : }
110 :
111 3 : DOCTEST_TEST_CASE("calculation of matrix elements") {
112 3 : auto &database = Database::get_global_instance();
113 :
114 3 : auto ket_s = KetAtomCreator()
115 6 : .set_species("Rb")
116 6 : .set_quantum_number("n", 60)
117 6 : .set_quantum_number("l", 0)
118 6 : .set_quantum_number("j", 0.5)
119 6 : .set_quantum_number("m", 0.5)
120 3 : .create(database);
121 :
122 3 : auto ket_p = KetAtomCreator()
123 6 : .set_species("Rb")
124 6 : .set_quantum_number("n", 60)
125 6 : .set_quantum_number("l", 1)
126 6 : .set_quantum_number("j", 0.5)
127 6 : .set_quantum_number("m", 0.5)
128 3 : .create(database);
129 :
130 3 : auto basis = BasisAtomCreator<double>()
131 6 : .set_species("Rb")
132 6 : .restrict_quantum_number("n", 59, 61)
133 6 : .restrict_quantum_number("l", 0, 1)
134 6 : .restrict_quantum_number("m", 0.5, 0.5)
135 3 : .create(database);
136 :
137 3 : SystemAtom<double> system(basis);
138 :
139 4 : auto get_corresponding_state_index = [&database](const auto &b,
140 4 : const std::shared_ptr<const KetAtom> &ket) {
141 4 : auto basis_ket = BasisAtomCreator<double>().add_ket(ket).create(database);
142 4 : Eigen::MatrixXd overlaps =
143 8 : Eigen::MatrixXd(b->get_matrix_elements(basis_ket, OperatorType::IDENTITY, 0))
144 : .cwiseAbs();
145 4 : Eigen::Index idx = 0;
146 4 : overlaps.row(0).maxCoeff(&idx);
147 4 : return idx;
148 4 : };
149 :
150 3 : DOCTEST_SUBCASE("calculate energy") {
151 1 : auto basis_ket_s = BasisAtomCreator<double>().add_ket(ket_s).create(database);
152 :
153 1 : auto m1 = basis_ket_s->get_matrix_elements(basis_ket_s, OperatorType::ENERGY, 0);
154 1 : DOCTEST_CHECK(m1.rows() == 1);
155 1 : DOCTEST_CHECK(m1.cols() == 1);
156 1 : double energy1 = m1.coeff(0, 0);
157 :
158 1 : auto m2 = basis->get_matrix_elements(basis_ket_s, OperatorType::ENERGY, 0);
159 1 : DOCTEST_CHECK(m2.rows() == 1);
160 1 : DOCTEST_CHECK(m2.cols() == basis->get_number_of_states());
161 1 : double energy2 = m2.coeff(0, static_cast<int>(get_corresponding_state_index(basis, ket_s)));
162 :
163 1 : double reference = ket_s->get_energy();
164 1 : DOCTEST_CHECK(std::abs(energy1 - reference) < 1e-11);
165 1 : DOCTEST_CHECK(std::abs(energy2 - reference) < 1e-11);
166 4 : }
167 :
168 3 : DOCTEST_SUBCASE("calculate electric dipole matrix element") {
169 1 : auto basis_ket_p = BasisAtomCreator<double>().add_ket(ket_p).create(database);
170 :
171 1 : auto m = basis->get_matrix_elements(basis_ket_p, OperatorType::ELECTRIC_DIPOLE, 0);
172 1 : DOCTEST_CHECK(m.rows() == 1);
173 1 : DOCTEST_CHECK(m.cols() == basis->get_number_of_states());
174 1 : double dipole = m.coeff(0, static_cast<int>(get_corresponding_state_index(basis, ket_s)));
175 :
176 1 : DOCTEST_CHECK(std::abs(dipole - 1247.6043831131365) < 1e-6);
177 4 : }
178 :
179 3 : DOCTEST_SUBCASE("calculate electric dipole matrix element with and without an induced dipole") {
180 : {
181 1 : auto state = basis->get_state(get_corresponding_state_index(basis, 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 - 0) < 1e-6);
189 1 : }
190 :
191 : {
192 1 : system.set_electric_field({0, 0, VOLT_PER_CM_IN_ATOMIC_UNITS});
193 1 : system.diagonalize(DiagonalizerEigen<double>());
194 1 : auto eigenbasis = system.get_eigenbasis();
195 1 : auto state = eigenbasis->get_state(get_corresponding_state_index(eigenbasis, ket_s));
196 :
197 1 : auto m = state->get_matrix_elements(state, OperatorType::ELECTRIC_DIPOLE, 0);
198 1 : DOCTEST_CHECK(m.rows() == 1);
199 1 : DOCTEST_CHECK(m.cols() == 1);
200 1 : double dipole = m.coeff(0, 0);
201 :
202 1 : DOCTEST_CHECK(std::abs(dipole - 135.04130863117354) < 1e-6);
203 1 : }
204 3 : }
205 3 : }
206 :
207 : } // namespace pairinteraction
|
