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