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