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