Line data Source code
1 : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/system/SystemAtom.hpp"
5 :
6 : #include "pairinteraction/basis/BasisAtom.hpp"
7 : #include "pairinteraction/basis/BasisAtomCreator.hpp"
8 : #include "pairinteraction/database/Database.hpp"
9 : #include "pairinteraction/diagonalize/DiagonalizerEigen.hpp"
10 : #include "pairinteraction/diagonalize/DiagonalizerFeast.hpp"
11 : #include "pairinteraction/diagonalize/DiagonalizerLapackeEvd.hpp"
12 : #include "pairinteraction/diagonalize/DiagonalizerLapackeEvr.hpp"
13 : #include "pairinteraction/diagonalize/diagonalize.hpp"
14 : #include "pairinteraction/enums/FloatType.hpp"
15 : #include "pairinteraction/ket/KetAtom.hpp"
16 : #include "pairinteraction/ket/KetAtomCreator.hpp"
17 :
18 : #include <Eigen/Eigenvalues>
19 : #include <doctest/doctest.h>
20 : #include <fmt/ranges.h>
21 :
22 : namespace pairinteraction {
23 :
24 : constexpr double VOLT_PER_CM_IN_ATOMIC_UNITS = 1 / 5.14220675112e9;
25 : constexpr double UM_IN_ATOMIC_UNITS = 1 / 5.29177210544e-5;
26 : constexpr double HARTREE_IN_GHZ = 6579683.920501762;
27 :
28 1 : DOCTEST_TEST_CASE("construct and diagonalize a small Hamiltonian") {
29 1 : auto &database = Database::get_global_instance();
30 1 : auto diagonalizer = DiagonalizerEigen<double>();
31 :
32 1 : auto ket1 = KetAtomCreator()
33 2 : .set_species("Rb")
34 1 : .set_quantum_number_n(60)
35 1 : .set_quantum_number_l(0)
36 1 : .set_quantum_number_j(0.5)
37 1 : .set_quantum_number_m(0.5)
38 1 : .create(database);
39 1 : auto ket2 = KetAtomCreator()
40 2 : .set_species("Rb")
41 1 : .set_quantum_number_n(60)
42 1 : .set_quantum_number_l(1)
43 1 : .set_quantum_number_j(0.5)
44 1 : .set_quantum_number_m(0.5)
45 1 : .create(database);
46 1 : auto basis = BasisAtomCreator<double>().append_ket(ket1).append_ket(ket2).create(database);
47 :
48 1 : auto system = SystemAtom<double>(basis);
49 1 : system.set_electric_field({0, 0, 0.0001});
50 :
51 1 : Eigen::MatrixXd tmp = Eigen::MatrixXd(1e5 * system.get_matrix()).array().round() / 1e5;
52 1 : std::vector<double> matrix_vector(tmp.data(), tmp.data() + tmp.size());
53 3 : DOCTEST_MESSAGE(fmt::format("Constructed: {}", fmt::join(matrix_vector, ", ")));
54 :
55 1 : system.diagonalize(diagonalizer);
56 1 : tmp = Eigen::MatrixXd(1e5 * system.get_matrix()).array().round() / 1e5;
57 1 : matrix_vector = std::vector<double>(tmp.data(), tmp.data() + tmp.size());
58 3 : DOCTEST_MESSAGE(fmt::format("Diagonalized: {}", fmt::join(matrix_vector, ", ")));
59 :
60 1 : Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver;
61 1 : eigensolver.compute(system.get_matrix());
62 1 : auto eigenenergies_eigen = eigensolver.eigenvalues();
63 1 : auto eigenenergies_pairinteraction = system.get_eigenenergies();
64 3 : for (int i = 0; i < eigenenergies_eigen.size(); ++i) {
65 2 : DOCTEST_CHECK(std::abs(eigenenergies_eigen(i) - eigenenergies_pairinteraction(i)) < 1e-10);
66 : }
67 1 : }
68 :
69 1 : DOCTEST_TEST_CASE("construct and diagonalize two Hamiltonians in parallel") {
70 1 : auto &database = Database::get_global_instance();
71 1 : auto diagonalizer = DiagonalizerEigen<std::complex<double>>();
72 :
73 1 : auto basis = BasisAtomCreator<std::complex<double>>()
74 2 : .set_species("Rb")
75 1 : .restrict_quantum_number_n(59, 61)
76 1 : .restrict_quantum_number_l(0, 1)
77 1 : .create(database);
78 :
79 1 : auto system1 = SystemAtom<std::complex<double>>(basis);
80 1 : system1.set_electric_field({0, 0, 0.0001});
81 :
82 1 : auto system2 = SystemAtom<std::complex<double>>(basis);
83 1 : system2.set_electric_field({0, 0, 0.0002});
84 :
85 1 : diagonalize<SystemAtom<std::complex<double>>>({system1, system2}, diagonalizer);
86 :
87 1 : auto matrix1 = system1.get_matrix();
88 1 : auto matrix2 = system2.get_matrix();
89 25 : for (int i = 0; i < matrix1.rows(); ++i) {
90 600 : for (int j = 0; j < matrix1.cols(); ++j) {
91 576 : if (i != j) {
92 552 : DOCTEST_CHECK(std::abs(matrix1.coeff(i, j)) < 1e-10);
93 552 : DOCTEST_CHECK(std::abs(matrix2.coeff(i, j)) < 1e-10);
94 : }
95 : }
96 : }
97 1 : }
98 :
99 0 : DOCTEST_TEST_CASE("construct and diagonalize multiple Hamiltonians in parallel" *
100 : doctest::skip(true)) {
101 : // TODO For a slow database, the fast parallelized construction of the tiny Hamiltonians seems
102 : // to lead to the wrong Hamiltonians (visible in the warning "The floating point error (5e-324
103 : // Hartree) is similar or larger than error estimated from the specified tolerance (0
104 : // Hartree)."). This could be caused by an issue with thread safety or memory access.
105 0 : int n = 10;
106 :
107 0 : auto &database = Database::get_global_instance();
108 0 : auto diagonalizer = DiagonalizerEigen<std::complex<double>>();
109 :
110 0 : auto basis = BasisAtomCreator<std::complex<double>>()
111 0 : .set_species("Sr87_mqdt")
112 0 : .restrict_quantum_number_nu(60, 61)
113 0 : .restrict_quantum_number_l(0, 1)
114 0 : .create(database);
115 :
116 0 : std::vector<SystemAtom<std::complex<double>>> systems;
117 0 : systems.reserve(n);
118 0 : for (int i = 0; i < n; ++i) {
119 0 : auto system = SystemAtom<std::complex<double>>(basis);
120 0 : system.set_electric_field({0, 0, 0.0001 * i});
121 0 : systems.push_back(std::move(system));
122 0 : }
123 :
124 0 : DOCTEST_MESSAGE("Basis size: ", basis->get_number_of_states());
125 :
126 0 : diagonalize<SystemAtom<std::complex<double>>>(systems, diagonalizer);
127 0 : }
128 :
129 2 : DOCTEST_TEST_CASE("construct and diagonalize a Hamiltonian using different methods") {
130 2 : auto &database = Database::get_global_instance();
131 :
132 2 : auto basis = BasisAtomCreator<std::complex<double>>()
133 4 : .set_species("Rb")
134 2 : .restrict_quantum_number_n(60, 61)
135 2 : .restrict_quantum_number_l(0, 1)
136 2 : .create(database);
137 :
138 : // Diagonalize using the Eigen library
139 2 : auto system = SystemAtom<std::complex<double>>(basis);
140 2 : system.set_electric_field({1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS,
141 : 3 * VOLT_PER_CM_IN_ATOMIC_UNITS});
142 :
143 2 : Eigen::MatrixXcd matrix = system.get_matrix();
144 2 : Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver;
145 2 : eigensolver.compute(matrix);
146 2 : auto eigenenergies_eigen = eigensolver.eigenvalues();
147 :
148 : // Create diagonalizers
149 2 : std::vector<std::unique_ptr<DiagonalizerInterface<std::complex<double>>>> diagonalizers;
150 2 : std::vector<double> rtols;
151 2 : double eps{};
152 2 : DOCTEST_SUBCASE("Double precision") {
153 1 : diagonalizers.push_back(std::make_unique<DiagonalizerEigen<std::complex<double>>>());
154 : #ifdef WITH_LAPACKE
155 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvd<std::complex<double>>>());
156 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvr<std::complex<double>>>());
157 : #endif
158 : #ifdef WITH_MKL
159 1 : diagonalizers.push_back(std::make_unique<DiagonalizerFeast<std::complex<double>>>(300));
160 : #endif
161 1 : rtols = {1e-1, 1e-6, 1e-14};
162 1 : eps = std::numeric_limits<double>::epsilon();
163 2 : }
164 :
165 2 : DOCTEST_SUBCASE("Single precision") {
166 1 : diagonalizers.push_back(
167 2 : std::make_unique<DiagonalizerEigen<std::complex<double>>>(FloatType::FLOAT32));
168 : #ifdef WITH_LAPACKE
169 1 : diagonalizers.push_back(
170 2 : std::make_unique<DiagonalizerLapackeEvd<std::complex<double>>>(FloatType::FLOAT32));
171 1 : diagonalizers.push_back(
172 2 : std::make_unique<DiagonalizerLapackeEvr<std::complex<double>>>(FloatType::FLOAT32));
173 : #endif
174 : #ifdef WITH_MKL
175 1 : diagonalizers.push_back(
176 2 : std::make_unique<DiagonalizerFeast<std::complex<double>>>(300, FloatType::FLOAT32));
177 : #endif
178 1 : rtols = {1e-1, 1e-6};
179 1 : eps = std::numeric_limits<float>::epsilon();
180 2 : }
181 :
182 : // Diagonalize using pairinteraction
183 7 : for (double rtol_eigenenergies : rtols) {
184 : double atol_eigenvectors =
185 5 : std::max(0.5 * rtol_eigenenergies / std::sqrt(basis->get_number_of_states()), 4 * eps);
186 5 : DOCTEST_MESSAGE("Precision: " << rtol_eigenenergies << " (rtol eigenenergies), "
187 : << atol_eigenvectors << " (atol eigenvectors)");
188 :
189 25 : for (const auto &diagonalizer : diagonalizers) {
190 20 : auto system = SystemAtom<std::complex<double>>(basis);
191 20 : system.set_electric_field({1 * VOLT_PER_CM_IN_ATOMIC_UNITS,
192 : 2 * VOLT_PER_CM_IN_ATOMIC_UNITS,
193 : 3 * VOLT_PER_CM_IN_ATOMIC_UNITS});
194 :
195 : // We specify a search interval because this is required if the FEAST routine is
196 : // used. To avoid overflows, the interval ranges from half the smallest possible
197 : // value to half the largest possible value.
198 20 : system.diagonalize(*diagonalizer, std::numeric_limits<float>::lowest() / 2,
199 20 : std::numeric_limits<float>::max() / 2, rtol_eigenenergies);
200 20 : auto eigenenergies_pairinteraction = system.get_eigenenergies();
201 20 : auto eigenvectors_pairinteraction = system.get_eigenbasis()->get_coefficients();
202 :
203 20 : DOCTEST_CHECK(
204 : (eigenenergies_eigen - eigenenergies_pairinteraction).array().abs().maxCoeff() <
205 : rtol_eigenenergies * matrix.norm());
206 :
207 340 : for (int i = 0; i < eigenvectors_pairinteraction.cols(); ++i) {
208 320 : DOCTEST_CHECK(abs(1 - eigenvectors_pairinteraction.col(i).norm()) <
209 : atol_eigenvectors * eigenvectors_pairinteraction.rows());
210 : }
211 20 : }
212 : }
213 2 : }
214 :
215 1 : DOCTEST_TEST_CASE("construct and diagonalize a Hamiltonian with energy restrictions") {
216 1 : double min_energy = 0.153355;
217 1 : double max_energy = 0.153360;
218 :
219 1 : auto &database = Database::get_global_instance();
220 :
221 1 : auto basis = BasisAtomCreator<double>()
222 2 : .set_species("Rb")
223 1 : .restrict_quantum_number_n(58, 62)
224 1 : .restrict_quantum_number_l(0, 1)
225 1 : .create(database);
226 :
227 : // Diagonalize using the Eigen library
228 1 : auto system = SystemAtom<double>(basis);
229 1 : system.set_electric_field(
230 : {1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
231 :
232 1 : Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver;
233 1 : eigensolver.compute(system.get_matrix());
234 1 : auto eigenenergies_all = eigensolver.eigenvalues();
235 1 : std::vector<double> eigenenergies_eigen;
236 41 : for (int i = 0; i < eigenenergies_all.size(); ++i) {
237 40 : if (eigenenergies_all[i] > min_energy && eigenenergies_all[i] < max_energy) {
238 8 : eigenenergies_eigen.push_back(eigenenergies_all[i]);
239 : }
240 : }
241 :
242 : // Create diagonalizer
243 1 : std::vector<std::unique_ptr<DiagonalizerInterface<double>>> diagonalizers;
244 1 : diagonalizers.push_back(std::make_unique<DiagonalizerEigen<double>>(FloatType::FLOAT64));
245 : #ifdef WITH_LAPACKE
246 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvd<double>>(FloatType::FLOAT64));
247 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvr<double>>(FloatType::FLOAT64));
248 : #endif
249 : #ifdef WITH_MKL
250 1 : diagonalizers.push_back(std::make_unique<DiagonalizerFeast<double>>(10, FloatType::FLOAT64));
251 : #endif
252 :
253 : // Diagonalize using pairinteraction
254 5 : for (const auto &diagonalizer : diagonalizers) {
255 4 : auto system = SystemAtom<double>(basis);
256 4 : system.set_electric_field(
257 : {1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
258 :
259 4 : system.diagonalize(*diagonalizer, min_energy, max_energy, 1e-6);
260 4 : auto eigenenergies_pairinteraction = system.get_eigenenergies();
261 :
262 4 : Eigen::MatrixXd tmp = (1e5 * eigenenergies_pairinteraction).array().round() / 1e5;
263 4 : std::vector<double> eigenenergies_vector(tmp.data(), tmp.data() + tmp.size());
264 12 : DOCTEST_MESSAGE(fmt::format("Eigenenergies: {}", fmt::join(eigenenergies_vector, ", ")));
265 :
266 4 : DOCTEST_CHECK(eigenenergies_eigen.size() == 8);
267 4 : DOCTEST_CHECK(eigenenergies_pairinteraction.size() == 8);
268 36 : for (size_t i = 0; i < eigenenergies_eigen.size(); ++i) {
269 32 : DOCTEST_CHECK(std::abs(eigenenergies_eigen[i] - eigenenergies_pairinteraction[i]) <
270 : 1e-10);
271 : }
272 4 : }
273 1 : }
274 :
275 : #ifdef WITH_MKL
276 : #include <Eigen/Dense>
277 : #include <mkl.h>
278 1 : DOCTEST_TEST_CASE("diagonalization with mkl") {
279 : // We loop several times to check for memory errors
280 11 : for (size_t i = 0; i < 10; ++i) {
281 : // Create a symmetric matrix
282 10 : int n = 100;
283 10 : Eigen::MatrixXd matrix = Eigen::MatrixXd::Random(n, n);
284 10 : matrix = (matrix + matrix.transpose()).eval();
285 10 : Eigen::VectorXd eigenenergies(n);
286 :
287 : // Diagonalize the matrix
288 : int info =
289 10 : LAPACKE_dsyev(LAPACK_COL_MAJOR, 'V', 'U', n, matrix.data(), n, eigenenergies.data());
290 10 : DOCTEST_CHECK(info == 0);
291 10 : }
292 1 : }
293 : #endif
294 :
295 1 : DOCTEST_TEST_CASE("handle it gracefully if no eigenenergies are within energy restrictions") {
296 1 : double min_energy = -1;
297 1 : double max_energy = -1;
298 :
299 1 : auto &database = Database::get_global_instance();
300 :
301 1 : auto basis = BasisAtomCreator<double>()
302 2 : .set_species("Rb")
303 1 : .restrict_quantum_number_n(58, 62)
304 1 : .restrict_quantum_number_l(0, 1)
305 1 : .create(database);
306 :
307 1 : std::vector<std::unique_ptr<DiagonalizerInterface<double>>> diagonalizers;
308 1 : diagonalizers.push_back(std::make_unique<DiagonalizerEigen<double>>());
309 : #ifdef WITH_LAPACKE
310 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvd<double>>());
311 : #endif
312 :
313 3 : for (const auto &diagonalizer : diagonalizers) {
314 2 : auto system = SystemAtom<double>(basis);
315 2 : system.set_electric_field(
316 : {1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
317 :
318 2 : system.diagonalize(*diagonalizer, min_energy, max_energy, 1e-6);
319 2 : auto eigenenergies_pairinteraction = system.get_eigenenergies();
320 :
321 2 : DOCTEST_CHECK(eigenenergies_pairinteraction.size() == 0);
322 2 : }
323 1 : }
324 :
325 1 : DOCTEST_TEST_CASE("atom ion interaction") {
326 1 : auto &database = Database::get_global_instance();
327 1 : DiagonalizerEigen<double> diagonalizer;
328 :
329 1 : auto ket = KetAtomCreator()
330 2 : .set_species("Rb")
331 1 : .set_quantum_number_n(60)
332 1 : .set_quantum_number_l(1)
333 1 : .set_quantum_number_j(0.5)
334 1 : .set_quantum_number_m(0.5)
335 1 : .create(database);
336 1 : double energy = ket->get_energy();
337 1 : double min_energy = energy - 50 / HARTREE_IN_GHZ;
338 1 : double max_energy = energy + 50 / HARTREE_IN_GHZ;
339 :
340 1 : auto basis = BasisAtomCreator<double>()
341 2 : .set_species("Rb")
342 1 : .restrict_quantum_number_n(58, 62)
343 1 : .restrict_quantum_number_l(0, 3)
344 1 : .restrict_quantum_number_m(0.5, 0.5)
345 1 : .create(database);
346 :
347 1 : auto system3 = SystemAtom<double>(basis);
348 1 : system3.set_ion_interaction_order(3);
349 1 : system3.set_ion_distance_vector({0, 0, 10 * UM_IN_ATOMIC_UNITS});
350 1 : system3.diagonalize(diagonalizer, min_energy, max_energy, 1e-6);
351 1 : auto energies3 = system3.get_eigenenergies();
352 :
353 1 : auto system2 = SystemAtom<double>(basis);
354 1 : system2.set_ion_interaction_order(2);
355 1 : system2.set_ion_distance_vector({0, 0, 10 * UM_IN_ATOMIC_UNITS});
356 1 : system2.diagonalize(diagonalizer, min_energy, max_energy, 1e-6);
357 1 : auto energies2 = system2.get_eigenenergies();
358 :
359 : // Ensure that the quadrupole order has a significant effect
360 1 : size_t num_energies = std::min(energies2.size(), energies3.size());
361 17 : for (size_t i = 0; i < num_energies; ++i) {
362 16 : DOCTEST_CHECK(std::abs(energies3[i] - energies2[i]) * HARTREE_IN_GHZ > 1e-6);
363 : }
364 1 : }
365 :
366 : } // namespace pairinteraction
|
