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