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/enums/SorterType.hpp"
16 : #include "pairinteraction/ket/KetAtom.hpp"
17 : #include "pairinteraction/ket/KetAtomCreator.hpp"
18 :
19 : #include <Eigen/Eigenvalues>
20 : #include <cmath>
21 : #include <doctest/doctest.h>
22 : #include <fmt/ranges.h>
23 : #include <stdexcept>
24 :
25 : namespace pairinteraction {
26 :
27 : constexpr double VOLT_PER_CM_IN_ATOMIC_UNITS = 1 / 5.14220675112e9;
28 : constexpr double UM_IN_ATOMIC_UNITS = 1 / 5.29177210544e-5;
29 : constexpr double HARTREE_IN_GHZ = 6579683.920501762;
30 :
31 1 : DOCTEST_TEST_CASE("construct and diagonalize a small Hamiltonian") {
32 1 : auto &database = Database::get_global_instance();
33 1 : auto diagonalizer = DiagonalizerEigen<double>();
34 :
35 1 : auto ket1 = KetAtomCreator(database)
36 2 : .set_species("Rb")
37 2 : .set_quantum_number("n", 60)
38 2 : .set_quantum_number("l", 0)
39 2 : .set_quantum_number("j", 0.5)
40 2 : .set_quantum_number("m", 0.5)
41 1 : .create();
42 1 : auto ket2 = KetAtomCreator(database)
43 2 : .set_species("Rb")
44 2 : .set_quantum_number("n", 60)
45 2 : .set_quantum_number("l", 1)
46 2 : .set_quantum_number("j", 0.5)
47 2 : .set_quantum_number("m", 0.5)
48 1 : .create();
49 1 : auto basis = BasisAtomCreator<double>(database).add_ket(ket1).add_ket(ket2).create();
50 :
51 1 : auto system = SystemAtom<double>(basis);
52 1 : system.set_electric_field({0, 0, 0.0001});
53 :
54 1 : Eigen::MatrixXd tmp = Eigen::MatrixXd(1e5 * system.get_matrix()).array().round() / 1e5;
55 1 : std::vector<double> matrix_vector(tmp.data(), tmp.data() + tmp.size());
56 3 : DOCTEST_MESSAGE(fmt::format("Constructed: {}", fmt::join(matrix_vector, ", ")));
57 :
58 1 : system.diagonalize(diagonalizer);
59 1 : tmp = Eigen::MatrixXd(1e5 * system.get_matrix()).array().round() / 1e5;
60 1 : matrix_vector = std::vector<double>(tmp.data(), tmp.data() + tmp.size());
61 3 : DOCTEST_MESSAGE(fmt::format("Diagonalized: {}", fmt::join(matrix_vector, ", ")));
62 :
63 1 : Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver;
64 1 : eigensolver.compute(system.get_matrix());
65 1 : auto eigenenergies_eigen = eigensolver.eigenvalues();
66 1 : auto eigenenergies_pairinteraction = system.get_eigenenergies();
67 3 : for (int i = 0; i < eigenenergies_eigen.size(); ++i) {
68 2 : DOCTEST_CHECK(std::abs(eigenenergies_eigen(i) - eigenenergies_pairinteraction(i)) < 1e-10);
69 : }
70 1 : }
71 :
72 1 : DOCTEST_TEST_CASE("construct and diagonalize two Hamiltonians in parallel") {
73 1 : auto &database = Database::get_global_instance();
74 1 : auto diagonalizer = DiagonalizerEigen<std::complex<double>>();
75 :
76 1 : auto basis = BasisAtomCreator<std::complex<double>>(database)
77 2 : .set_species("Rb")
78 2 : .restrict_quantum_number("n", 59, 61)
79 2 : .restrict_quantum_number("l", 0, 1)
80 1 : .create();
81 :
82 1 : auto system1 = SystemAtom<std::complex<double>>(basis);
83 1 : system1.set_electric_field({0, 0, 0.0001});
84 :
85 1 : auto system2 = SystemAtom<std::complex<double>>(basis);
86 1 : system2.set_electric_field({0, 0, 0.0002});
87 :
88 1 : diagonalize<SystemAtom<std::complex<double>>>({system1, system2}, diagonalizer);
89 :
90 1 : auto matrix1 = system1.get_matrix();
91 1 : auto matrix2 = system2.get_matrix();
92 25 : for (int i = 0; i < matrix1.rows(); ++i) {
93 600 : for (int j = 0; j < matrix1.cols(); ++j) {
94 576 : if (i != j) {
95 552 : DOCTEST_CHECK(std::abs(matrix1.coeff(i, j)) < 1e-10);
96 552 : DOCTEST_CHECK(std::abs(matrix2.coeff(i, j)) < 1e-10);
97 : }
98 : }
99 : }
100 1 : }
101 :
102 1 : DOCTEST_TEST_CASE("construct an atomic Hamiltonian in a non-canonical atomic basis") {
103 1 : auto &database = Database::get_global_instance();
104 :
105 1 : auto basis = BasisAtomCreator<double>(database)
106 2 : .set_species("Rb")
107 2 : .restrict_quantum_number("n", 60, 61)
108 2 : .restrict_quantum_number("l", 0, 1)
109 2 : .restrict_quantum_number("m", 0.5, 0.5)
110 1 : .create();
111 1 : DOCTEST_REQUIRE(basis->get_number_of_states() >= 2);
112 :
113 1 : SystemAtom<double> reference_system(basis);
114 1 : reference_system.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
115 1 : const auto &reference_matrix = reference_system.get_matrix();
116 :
117 : Eigen::SparseMatrix<double, Eigen::RowMajor> transformation(
118 1 : static_cast<Eigen::Index>(basis->get_number_of_states()),
119 2 : static_cast<Eigen::Index>(basis->get_number_of_states()));
120 1 : transformation.setIdentity();
121 :
122 1 : double inverse_sqrt_two = 1 / std::sqrt(2.0);
123 1 : transformation.coeffRef(0, 0) = inverse_sqrt_two;
124 1 : transformation.coeffRef(1, 0) = inverse_sqrt_two;
125 1 : transformation.coeffRef(0, 1) = inverse_sqrt_two;
126 1 : transformation.coeffRef(1, 1) = -inverse_sqrt_two;
127 1 : transformation.makeCompressed();
128 :
129 1 : auto transformed_basis = basis->transformed(transformation);
130 1 : SystemAtom<double> transformed_system(transformed_basis);
131 1 : transformed_system.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
132 :
133 : Eigen::SparseMatrix<double, Eigen::RowMajor> expected_matrix =
134 1 : transformation.adjoint() * reference_matrix * transformation;
135 :
136 1 : DOCTEST_CHECK(transformed_system.get_matrix().isApprox(expected_matrix, 1e-11));
137 1 : }
138 :
139 0 : DOCTEST_TEST_CASE("construct and diagonalize multiple Hamiltonians in parallel" *
140 : doctest::skip(true)) {
141 : // TODO For a slow database, the fast parallelized construction of the tiny Hamiltonians seems
142 : // to lead to the wrong Hamiltonians (visible in the warning "The floating point error (5e-324
143 : // Hartree) is similar or larger than error estimated from the specified tolerance (0
144 : // Hartree)."). This could be caused by an issue with thread safety or memory access.
145 0 : int n = 10;
146 :
147 0 : auto &database = Database::get_global_instance();
148 0 : auto diagonalizer = DiagonalizerEigen<std::complex<double>>();
149 :
150 0 : auto basis = BasisAtomCreator<std::complex<double>>(database)
151 0 : .set_species("Sr87_mqdt")
152 0 : .restrict_quantum_number("nu", 60, 61)
153 0 : .restrict_quantum_number("l", 0, 1)
154 0 : .create();
155 :
156 0 : std::vector<SystemAtom<std::complex<double>>> systems;
157 0 : systems.reserve(n);
158 0 : for (int i = 0; i < n; ++i) {
159 0 : auto system = SystemAtom<std::complex<double>>(basis);
160 0 : system.set_electric_field({0, 0, 0.0001 * i});
161 0 : systems.push_back(std::move(system));
162 0 : }
163 :
164 0 : DOCTEST_MESSAGE("Basis size: ", basis->get_number_of_states());
165 :
166 0 : diagonalize<SystemAtom<std::complex<double>>>(systems, diagonalizer);
167 0 : }
168 :
169 2 : DOCTEST_TEST_CASE("construct and diagonalize a Hamiltonian using different methods") {
170 2 : auto &database = Database::get_global_instance();
171 :
172 2 : auto basis = BasisAtomCreator<std::complex<double>>(database)
173 4 : .set_species("Rb")
174 4 : .restrict_quantum_number("n", 60, 61)
175 4 : .restrict_quantum_number("l", 0, 1)
176 2 : .create();
177 :
178 : // Diagonalize using the Eigen library
179 2 : auto system = SystemAtom<std::complex<double>>(basis);
180 2 : system.set_electric_field({1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS,
181 : 3 * VOLT_PER_CM_IN_ATOMIC_UNITS});
182 :
183 2 : Eigen::MatrixXcd matrix = system.get_matrix();
184 2 : Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> eigensolver;
185 2 : eigensolver.compute(matrix);
186 2 : auto eigenenergies_eigen = eigensolver.eigenvalues();
187 :
188 : // Create diagonalizers
189 2 : std::vector<std::unique_ptr<DiagonalizerInterface<std::complex<double>>>> diagonalizers;
190 2 : std::vector<double> rtols;
191 2 : double eps{};
192 2 : DOCTEST_SUBCASE("Double precision") {
193 1 : diagonalizers.push_back(std::make_unique<DiagonalizerEigen<std::complex<double>>>());
194 : #ifdef WITH_LAPACKE
195 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvd<std::complex<double>>>());
196 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvr<std::complex<double>>>());
197 : #endif
198 : #ifdef WITH_MKL
199 1 : diagonalizers.push_back(std::make_unique<DiagonalizerFeast<std::complex<double>>>(300));
200 : #endif
201 1 : rtols = {1e-1, 1e-6, 1e-14};
202 1 : eps = std::numeric_limits<double>::epsilon();
203 2 : }
204 :
205 2 : DOCTEST_SUBCASE("Single precision") {
206 1 : diagonalizers.push_back(
207 2 : std::make_unique<DiagonalizerEigen<std::complex<double>>>(FloatType::FLOAT32));
208 : #ifdef WITH_LAPACKE
209 1 : diagonalizers.push_back(
210 2 : std::make_unique<DiagonalizerLapackeEvd<std::complex<double>>>(FloatType::FLOAT32));
211 1 : diagonalizers.push_back(
212 2 : std::make_unique<DiagonalizerLapackeEvr<std::complex<double>>>(FloatType::FLOAT32));
213 : #endif
214 : #ifdef WITH_MKL
215 1 : diagonalizers.push_back(
216 2 : std::make_unique<DiagonalizerFeast<std::complex<double>>>(300, FloatType::FLOAT32));
217 : #endif
218 1 : rtols = {1e-1, 1e-6};
219 1 : eps = std::numeric_limits<float>::epsilon();
220 2 : }
221 :
222 : // Diagonalize using pairinteraction
223 7 : for (double rtol_eigenenergies : rtols) {
224 : double atol_eigenvectors =
225 5 : std::max(0.5 * rtol_eigenenergies / std::sqrt(basis->get_number_of_states()), 5 * eps);
226 5 : DOCTEST_MESSAGE("Precision: " << rtol_eigenenergies << " (rtol eigenenergies), "
227 : << atol_eigenvectors << " (atol eigenvectors)");
228 :
229 25 : for (const auto &diagonalizer : diagonalizers) {
230 20 : auto system = SystemAtom<std::complex<double>>(basis);
231 20 : system.set_electric_field({1 * VOLT_PER_CM_IN_ATOMIC_UNITS,
232 : 2 * VOLT_PER_CM_IN_ATOMIC_UNITS,
233 : 3 * VOLT_PER_CM_IN_ATOMIC_UNITS});
234 :
235 : // We specify a search interval because this is required if the FEAST routine is
236 : // used. To avoid overflows, the interval ranges from half the smallest possible
237 : // value to half the largest possible value.
238 20 : system.diagonalize(*diagonalizer, std::numeric_limits<float>::lowest() / 2,
239 20 : std::numeric_limits<float>::max() / 2, rtol_eigenenergies);
240 20 : auto eigenenergies_pairinteraction = system.get_eigenenergies();
241 20 : auto eigenvectors_pairinteraction = system.get_eigenbasis()->get_coefficients();
242 :
243 20 : DOCTEST_CHECK(
244 : (eigenenergies_eigen - eigenenergies_pairinteraction).array().abs().maxCoeff() <
245 : rtol_eigenenergies * matrix.norm());
246 :
247 340 : for (int i = 0; i < eigenvectors_pairinteraction.cols(); ++i) {
248 320 : DOCTEST_CHECK(abs(1 - eigenvectors_pairinteraction.col(i).norm()) <
249 : atol_eigenvectors * eigenvectors_pairinteraction.rows());
250 : }
251 20 : }
252 : }
253 2 : }
254 :
255 1 : DOCTEST_TEST_CASE("construct and diagonalize a Hamiltonian with energy restrictions") {
256 1 : double min_energy = 0.153355;
257 1 : double max_energy = 0.153360;
258 :
259 1 : auto &database = Database::get_global_instance();
260 :
261 1 : auto basis = BasisAtomCreator<double>(database)
262 2 : .set_species("Rb")
263 2 : .restrict_quantum_number("n", 58, 62)
264 2 : .restrict_quantum_number("l", 0, 1)
265 1 : .create();
266 :
267 : // Diagonalize using the Eigen library
268 1 : auto system = SystemAtom<double>(basis);
269 1 : system.set_electric_field(
270 : {1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
271 :
272 1 : Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver;
273 1 : eigensolver.compute(system.get_matrix());
274 1 : auto eigenenergies_all = eigensolver.eigenvalues();
275 1 : std::vector<double> eigenenergies_eigen;
276 41 : for (int i = 0; i < eigenenergies_all.size(); ++i) {
277 40 : if (eigenenergies_all[i] > min_energy && eigenenergies_all[i] < max_energy) {
278 8 : eigenenergies_eigen.push_back(eigenenergies_all[i]);
279 : }
280 : }
281 :
282 : // Create diagonalizer
283 1 : std::vector<std::unique_ptr<DiagonalizerInterface<double>>> diagonalizers;
284 1 : diagonalizers.push_back(std::make_unique<DiagonalizerEigen<double>>(FloatType::FLOAT64));
285 : #ifdef WITH_LAPACKE
286 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvd<double>>(FloatType::FLOAT64));
287 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvr<double>>(FloatType::FLOAT64));
288 : #endif
289 : #ifdef WITH_MKL
290 1 : diagonalizers.push_back(std::make_unique<DiagonalizerFeast<double>>(10, FloatType::FLOAT64));
291 : #endif
292 :
293 : // Diagonalize using pairinteraction
294 5 : for (const auto &diagonalizer : diagonalizers) {
295 4 : auto system = SystemAtom<double>(basis);
296 4 : system.set_electric_field(
297 : {1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
298 :
299 4 : system.diagonalize(*diagonalizer, min_energy, max_energy, 1e-6);
300 4 : auto eigenenergies_pairinteraction = system.get_eigenenergies();
301 :
302 4 : Eigen::MatrixXd tmp = (1e5 * eigenenergies_pairinteraction).array().round() / 1e5;
303 4 : std::vector<double> eigenenergies_vector(tmp.data(), tmp.data() + tmp.size());
304 12 : DOCTEST_MESSAGE(fmt::format("Eigenenergies: {}", fmt::join(eigenenergies_vector, ", ")));
305 :
306 4 : DOCTEST_CHECK(eigenenergies_eigen.size() == 8);
307 4 : DOCTEST_CHECK(eigenenergies_pairinteraction.size() == 8);
308 36 : for (size_t i = 0; i < eigenenergies_eigen.size(); ++i) {
309 32 : DOCTEST_CHECK(std::abs(eigenenergies_eigen[i] - eigenenergies_pairinteraction[i]) <
310 : 1e-10);
311 : }
312 4 : }
313 1 : }
314 :
315 2 : DOCTEST_TEST_CASE("the eigenbasis knows the eigenenergies") {
316 2 : auto &database = Database::get_global_instance();
317 2 : auto diagonalizer = DiagonalizerEigen<double>();
318 :
319 2 : DOCTEST_SUBCASE("diagonalized Hamiltonian") {
320 1 : auto basis = BasisAtomCreator<double>(database)
321 2 : .set_species("Rb")
322 2 : .restrict_quantum_number("n", 58, 62)
323 2 : .restrict_quantum_number("l", 0, 1)
324 1 : .create();
325 :
326 3 : for (bool restrict_energy : {false, true}) {
327 2 : auto system = SystemAtom<double>(basis);
328 2 : system.set_electric_field(
329 : {1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
330 2 : if (restrict_energy) {
331 1 : system.diagonalize(diagonalizer, 0.153355, 0.153360);
332 : } else {
333 1 : system.diagonalize(diagonalizer);
334 : }
335 :
336 2 : auto eigenbasis = system.get_eigenbasis();
337 2 : auto eigenenergies = system.get_eigenenergies();
338 2 : DOCTEST_REQUIRE(eigenbasis->get_number_of_states() ==
339 : static_cast<size_t>(eigenenergies.size()));
340 50 : for (Eigen::Index i = 0; i < eigenenergies.size(); ++i) {
341 48 : DOCTEST_CHECK(eigenbasis->get_energy(static_cast<size_t>(i)) == eigenenergies[i]);
342 : }
343 2 : }
344 3 : }
345 :
346 2 : DOCTEST_SUBCASE("Hamiltonian that is diagonal without diagonalization") {
347 : // For a single fine structure level, the Zeeman Hamiltonian is diagonal in the canonical
348 : // basis
349 1 : auto basis = BasisAtomCreator<double>(database)
350 2 : .set_species("Rb")
351 2 : .restrict_quantum_number("n", 60, 60)
352 2 : .restrict_quantum_number("l", 0, 0)
353 1 : .create();
354 1 : auto system = SystemAtom<double>(basis);
355 1 : system.set_magnetic_field({0, 0, 1e-6});
356 :
357 1 : auto eigenbasis = system.get_eigenbasis();
358 1 : auto eigenenergies = system.get_eigenenergies();
359 1 : DOCTEST_REQUIRE(eigenenergies.size() == 2);
360 3 : for (Eigen::Index i = 0; i < eigenenergies.size(); ++i) {
361 2 : auto state_index = static_cast<size_t>(i);
362 2 : DOCTEST_CHECK(eigenbasis->get_energy(state_index) == eigenenergies[i]);
363 2 : DOCTEST_CHECK(eigenbasis->get_energy(state_index) !=
364 : eigenbasis->get_ket(state_index)->get_energy());
365 : }
366 3 : }
367 2 : }
368 :
369 : #ifdef WITH_MKL
370 : #include <Eigen/Dense>
371 : #include <mkl.h>
372 1 : DOCTEST_TEST_CASE("diagonalization with mkl") {
373 : // We loop several times to check for memory errors
374 11 : for (size_t i = 0; i < 10; ++i) {
375 : // Create a symmetric matrix
376 10 : int n = 100;
377 10 : Eigen::MatrixXd matrix = Eigen::MatrixXd::Random(n, n);
378 10 : matrix = (matrix + matrix.transpose()).eval();
379 10 : Eigen::VectorXd eigenenergies(n);
380 :
381 : // Diagonalize the matrix
382 : int info =
383 10 : LAPACKE_dsyev(LAPACK_COL_MAJOR, 'V', 'U', n, matrix.data(), n, eigenenergies.data());
384 10 : DOCTEST_CHECK(info == 0);
385 10 : }
386 1 : }
387 : #endif
388 :
389 1 : DOCTEST_TEST_CASE("handle it gracefully if no eigenenergies are within energy restrictions") {
390 1 : double min_energy = -1;
391 1 : double max_energy = -1;
392 :
393 1 : auto &database = Database::get_global_instance();
394 :
395 1 : auto basis = BasisAtomCreator<double>(database)
396 2 : .set_species("Rb")
397 2 : .restrict_quantum_number("n", 58, 62)
398 2 : .restrict_quantum_number("l", 0, 1)
399 1 : .create();
400 :
401 1 : std::vector<std::unique_ptr<DiagonalizerInterface<double>>> diagonalizers;
402 1 : diagonalizers.push_back(std::make_unique<DiagonalizerEigen<double>>());
403 : #ifdef WITH_LAPACKE
404 1 : diagonalizers.push_back(std::make_unique<DiagonalizerLapackeEvd<double>>());
405 : #endif
406 :
407 3 : for (const auto &diagonalizer : diagonalizers) {
408 2 : auto system = SystemAtom<double>(basis);
409 2 : system.set_electric_field(
410 : {1 * VOLT_PER_CM_IN_ATOMIC_UNITS, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
411 :
412 2 : system.diagonalize(*diagonalizer, min_energy, max_energy, 1e-6);
413 2 : auto eigenenergies_pairinteraction = system.get_eigenenergies();
414 :
415 2 : DOCTEST_CHECK(eigenenergies_pairinteraction.size() == 0);
416 2 : }
417 1 : }
418 :
419 1 : DOCTEST_TEST_CASE("atom ion interaction") {
420 1 : auto &database = Database::get_global_instance();
421 1 : DiagonalizerEigen<double> diagonalizer;
422 :
423 1 : auto ket = KetAtomCreator(database)
424 2 : .set_species("Rb")
425 2 : .set_quantum_number("n", 60)
426 2 : .set_quantum_number("l", 1)
427 2 : .set_quantum_number("j", 0.5)
428 2 : .set_quantum_number("m", 0.5)
429 1 : .create();
430 1 : double energy = ket->get_energy();
431 1 : double min_energy = energy - 50 / HARTREE_IN_GHZ;
432 1 : double max_energy = energy + 50 / HARTREE_IN_GHZ;
433 :
434 1 : auto basis = BasisAtomCreator<double>(database)
435 2 : .set_species("Rb")
436 2 : .restrict_quantum_number("n", 58, 62)
437 2 : .restrict_quantum_number("l", 0, 3)
438 2 : .restrict_quantum_number("m", 0.5, 0.5)
439 1 : .create();
440 :
441 1 : auto system3 = SystemAtom<double>(basis);
442 1 : system3.set_ion_interaction_order(3);
443 1 : system3.set_ion_distance_vector({0, 0, 10 * UM_IN_ATOMIC_UNITS});
444 1 : system3.diagonalize(diagonalizer, min_energy, max_energy, 1e-6);
445 1 : auto energies3 = system3.get_eigenenergies();
446 :
447 1 : auto system2 = SystemAtom<double>(basis);
448 1 : system2.set_ion_interaction_order(2);
449 1 : system2.set_ion_distance_vector({0, 0, 10 * UM_IN_ATOMIC_UNITS});
450 1 : system2.diagonalize(diagonalizer, min_energy, max_energy, 1e-6);
451 1 : auto energies2 = system2.get_eigenenergies();
452 :
453 : // Ensure that the quadrupole order has a significant effect
454 1 : size_t num_energies = std::min(energies2.size(), energies3.size());
455 17 : for (size_t i = 0; i < num_energies; ++i) {
456 16 : DOCTEST_CHECK(std::abs(energies3[i] - energies2[i]) * HARTREE_IN_GHZ > 1e-6);
457 : }
458 1 : }
459 :
460 1 : DOCTEST_TEST_CASE("obtain the blocks of a Hamiltonian") {
461 1 : auto &database = Database::get_global_instance();
462 :
463 1 : auto basis = BasisAtomCreator<double>(database)
464 2 : .set_species("Rb")
465 2 : .restrict_quantum_number("n", 60, 60)
466 2 : .restrict_quantum_number("l", 0, 1)
467 1 : .create();
468 :
469 1 : auto system = SystemAtom<double>(basis);
470 1 : system.set_electric_field({0, 0, 0.0001});
471 :
472 : // The electric field conserves the quantum number m, thus the Hamiltonian can be
473 : // block-diagonalized by m
474 1 : system.transform(system.get_sorter({SorterType::QUANTUM_NUMBER_M}));
475 1 : auto blocks = system.get_indices_of_blocks({SorterType::QUANTUM_NUMBER_M});
476 1 : DOCTEST_CHECK(blocks.size() > 1);
477 :
478 1 : size_t expected_start = 0;
479 5 : for (const auto &block : blocks) {
480 4 : DOCTEST_MESSAGE("Block from ", block.start, " to ", block.end);
481 4 : DOCTEST_CHECK(block.start == expected_start);
482 4 : expected_start = block.end;
483 : }
484 1 : DOCTEST_CHECK(expected_start == basis->get_number_of_states());
485 :
486 : // Blocks are obtained to diagonalize the Hamiltonian block by block. Because energy blocks
487 : // would require an already diagonal Hamiltonian, they cannot be obtained ...
488 3 : DOCTEST_CHECK_THROWS_AS(system.get_indices_of_blocks({SorterType::ENERGY}),
489 : std::invalid_argument);
490 :
491 : // ... whereas sorting by the energy is supported
492 1 : DOCTEST_CHECK_NOTHROW(system.get_sorter({SorterType::ENERGY}));
493 :
494 : // A repeated label cannot influence the order anymore and must be ignored instead of
495 : // terminating the program
496 1 : DOCTEST_CHECK_NOTHROW(system.get_sorter({SorterType::ENERGY, SorterType::ENERGY}));
497 1 : DOCTEST_CHECK_NOTHROW(
498 : system.get_sorter({SorterType::QUANTUM_NUMBER_M, SorterType::ENERGY, SorterType::ENERGY}));
499 1 : }
500 :
501 1 : DOCTEST_TEST_CASE("sort a Hamiltonian by several labels") {
502 1 : auto &database = Database::get_global_instance();
503 :
504 1 : auto basis = BasisAtomCreator<double>(database)
505 2 : .set_species("Rb")
506 2 : .restrict_quantum_number("n", 60, 60)
507 2 : .restrict_quantum_number("l", 0, 1)
508 1 : .create();
509 :
510 : // The first label is the primary sorting criterion, the following labels only break ties
511 : {
512 1 : auto system = SystemAtom<double>(basis);
513 1 : system.set_electric_field({0, 0, 0.0001});
514 1 : system.transform(system.get_sorter({SorterType::QUANTUM_NUMBER_M, SorterType::ENERGY}));
515 :
516 : // Sorting by m first means that states of equal m are contiguous ...
517 1 : auto blocks = system.get_indices_of_blocks({SorterType::QUANTUM_NUMBER_M});
518 1 : DOCTEST_CHECK(blocks.size() > 1);
519 :
520 : // ... and that within each block, the energies are ascending
521 1 : const auto &matrix = system.get_matrix();
522 5 : for (const auto &block : blocks) {
523 8 : for (auto i = static_cast<long>(block.start) + 1; i < static_cast<long>(block.end);
524 : ++i) {
525 4 : DOCTEST_CHECK(matrix.coeff(i - 1, i - 1) <= matrix.coeff(i, i));
526 : }
527 : }
528 1 : }
529 :
530 : // If the energy is the first label, the states are sorted by the energy globally
531 : {
532 1 : auto system = SystemAtom<double>(basis);
533 1 : system.set_electric_field({0, 0, 0.0001});
534 1 : system.transform(system.get_sorter({SorterType::ENERGY, SorterType::QUANTUM_NUMBER_M}));
535 :
536 1 : const auto &matrix = system.get_matrix();
537 8 : for (long i = 1; i < matrix.rows(); ++i) {
538 7 : DOCTEST_CHECK(matrix.coeff(i - 1, i - 1) <= matrix.coeff(i, i));
539 : }
540 1 : }
541 :
542 : // If the energy is enclosed by other labels, the labels before the energy take precedence over
543 : // the energy and the energy takes precedence over the labels after it
544 : {
545 1 : auto system = SystemAtom<double>(basis);
546 1 : system.set_electric_field({0, 0, 0.0001});
547 1 : system.transform(system.get_sorter(
548 : {SorterType::PARITY, SorterType::ENERGY, SorterType::QUANTUM_NUMBER_M}));
549 :
550 1 : auto sorted_basis = system.get_basis();
551 1 : size_t num_states = sorted_basis->get_number_of_states();
552 1 : const auto &matrix = system.get_matrix();
553 :
554 42 : auto parity = [&](size_t i) { return sorted_basis->get_quantum_number("parity", i); };
555 24 : auto energy = [&](size_t i) {
556 24 : return matrix.coeff(static_cast<long>(i), static_cast<long>(i));
557 1 : };
558 24 : auto quantum_number_m = [&](size_t i) { return sorted_basis->get_quantum_number("m", i); };
559 :
560 : // The parity is the primary criterion, thus it is ascending globally ...
561 8 : for (size_t i = 1; i < num_states; ++i) {
562 7 : DOCTEST_CHECK(parity(i - 1) <= parity(i));
563 : }
564 :
565 : // ... the energy is only ascending within a block of equal parity ...
566 8 : for (size_t i = 1; i < num_states; ++i) {
567 7 : if (parity(i - 1) == parity(i)) {
568 6 : DOCTEST_CHECK(energy(i - 1) <= energy(i));
569 : }
570 : }
571 :
572 : // ... and the quantum number m only breaks ties between states of equal parity and energy
573 1 : size_t num_tie_breaks = 0;
574 8 : for (size_t i = 1; i < num_states; ++i) {
575 7 : if (parity(i - 1) == parity(i) && energy(i - 1) == energy(i)) {
576 5 : DOCTEST_CHECK(quantum_number_m(i - 1) <= quantum_number_m(i));
577 5 : ++num_tie_breaks;
578 : }
579 : }
580 :
581 : // Check that the basis contains degenerate states so that the tie breaking by m above is
582 : // not vacuous
583 1 : DOCTEST_CHECK(num_tie_breaks > 0);
584 :
585 : // If m was the primary criterion, m would be ascending globally instead
586 1 : bool m_ascending_globally = true;
587 8 : for (size_t i = 1; i < num_states; ++i) {
588 7 : if (quantum_number_m(i - 1) > quantum_number_m(i)) {
589 2 : m_ascending_globally = false;
590 : }
591 : }
592 1 : DOCTEST_CHECK_FALSE(m_ascending_globally);
593 1 : }
594 1 : }
595 :
596 : } // namespace pairinteraction
|
