Line data Source code
1 : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/system/SystemPair.hpp"
5 :
6 : #include "pairinteraction/basis/BasisAtom.hpp"
7 : #include "pairinteraction/basis/BasisAtomCreator.hpp"
8 : #include "pairinteraction/basis/BasisPair.hpp"
9 : #include "pairinteraction/basis/BasisPairCreator.hpp"
10 : #include "pairinteraction/database/Database.hpp"
11 : #include "pairinteraction/diagonalize/DiagonalizerEigen.hpp"
12 : #include "pairinteraction/diagonalize/DiagonalizerFeast.hpp"
13 : #include "pairinteraction/diagonalize/DiagonalizerLapackeEvr.hpp"
14 : #include "pairinteraction/diagonalize/diagonalize.hpp"
15 : #include "pairinteraction/ket/KetAtom.hpp"
16 : #include "pairinteraction/ket/KetAtomCreator.hpp"
17 : #include "pairinteraction/system/SystemAtom.hpp"
18 : #include "pairinteraction/utils/Range.hpp"
19 :
20 : #include <array>
21 : #include <cmath>
22 : #include <doctest/doctest.h>
23 : #include <fmt/ranges.h>
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 a pair Hamiltonian") {
32 1 : auto &database = Database::get_global_instance();
33 1 : auto diagonalizer = DiagonalizerEigen<double>();
34 :
35 : // Construct and diagonalize the constituent systems
36 1 : auto basis = BasisAtomCreator<double>(database)
37 2 : .set_species("Rb")
38 2 : .restrict_quantum_number("n", 60, 61)
39 2 : .restrict_quantum_number("l", 0, 1)
40 2 : .restrict_quantum_number("m", -0.5, 0.5)
41 1 : .create();
42 :
43 1 : SystemAtom<double> system1(basis);
44 1 : SystemAtom<double> system2(basis);
45 1 : system1.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
46 1 : system2.set_electric_field({0, 0, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS});
47 1 : diagonalize<SystemAtom<double>>({system1, system2}, diagonalizer);
48 :
49 : // Construct and diagonalize the system_pair
50 : auto basis_pair =
51 1 : BasisPairCreator<double>(system1.get_eigenbasis(), system2.get_eigenbasis()).create();
52 1 : DOCTEST_MESSAGE("Number of states in pair basis: ", basis_pair->get_number_of_states());
53 :
54 1 : auto system_pair = SystemPair<double>(basis_pair);
55 1 : system_pair.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
56 1 : system_pair.diagonalize(diagonalizer);
57 :
58 : // Print the largest and smallest eigenenergies
59 1 : auto eigenenergies = system_pair.get_eigenenergies();
60 1 : DOCTEST_MESSAGE("Lowest energy: ", eigenenergies.minCoeff());
61 1 : DOCTEST_MESSAGE("Highest energy: ", eigenenergies.maxCoeff());
62 1 : }
63 :
64 1 : DOCTEST_TEST_CASE("construct a pair Hamiltonian in a non-canonical pair basis") {
65 1 : auto &database = Database::get_global_instance();
66 :
67 1 : auto basis = BasisAtomCreator<double>(database)
68 2 : .set_species("Rb")
69 2 : .restrict_quantum_number("n", 60, 61)
70 2 : .restrict_quantum_number("l", 0, 1)
71 2 : .restrict_quantum_number("m", -0.5, 0.5)
72 1 : .create();
73 :
74 1 : auto diagonalizer = DiagonalizerEigen<double>();
75 1 : SystemAtom<double> system1(basis);
76 1 : SystemAtom<double> system2(basis);
77 1 : system1.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
78 1 : system2.set_electric_field({0, 0, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS});
79 1 : diagonalize<SystemAtom<double>>({system1, system2}, diagonalizer);
80 :
81 : auto pair_basis =
82 1 : BasisPairCreator<double>(system1.get_eigenbasis(), system2.get_eigenbasis()).create();
83 1 : DOCTEST_REQUIRE(pair_basis->get_number_of_states() >= 2);
84 :
85 1 : SystemPair<double> reference_system(pair_basis);
86 1 : reference_system.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
87 1 : const auto &reference_matrix = reference_system.get_matrix();
88 :
89 : Eigen::SparseMatrix<double, Eigen::RowMajor> transformation(
90 1 : static_cast<Eigen::Index>(pair_basis->get_number_of_states()),
91 2 : static_cast<Eigen::Index>(pair_basis->get_number_of_states()));
92 1 : transformation.setIdentity();
93 :
94 1 : double inverse_sqrt_two = 1 / std::sqrt(2.0);
95 1 : transformation.coeffRef(0, 0) = inverse_sqrt_two;
96 1 : transformation.coeffRef(1, 0) = inverse_sqrt_two;
97 1 : transformation.coeffRef(0, 1) = inverse_sqrt_two;
98 1 : transformation.coeffRef(1, 1) = -inverse_sqrt_two;
99 1 : transformation.makeCompressed();
100 :
101 1 : auto transformed_pair_basis = pair_basis->transformed(transformation);
102 1 : SystemPair<double> transformed_system(transformed_pair_basis);
103 1 : transformed_system.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
104 :
105 : Eigen::SparseMatrix<double, Eigen::RowMajor> expected_matrix =
106 1 : transformation.adjoint() * reference_matrix * transformation;
107 :
108 1 : DOCTEST_CHECK(transformed_system.get_matrix().isApprox(expected_matrix, 1e-11));
109 1 : }
110 :
111 1 : DOCTEST_TEST_CASE("atom ion pair interaction") {
112 : // The interaction of a Rydberg atom with a Rydberg ion in a single state must reproduce the
113 : // interaction of the Rydberg atom with a classical point charge implemented in SystemAtom
114 1 : auto &database = Database::get_global_instance();
115 1 : auto diagonalizer = DiagonalizerEigen<double>();
116 :
117 1 : auto basis_atom = BasisAtomCreator<double>(database)
118 2 : .set_species("Rb")
119 2 : .restrict_quantum_number("n", 58, 62)
120 2 : .restrict_quantum_number("l", 0, 3)
121 2 : .restrict_quantum_number("m", 0.5, 0.5)
122 1 : .create();
123 :
124 1 : auto ket_ion = KetAtomCreator(database)
125 2 : .set_species("Sr88_ion")
126 2 : .set_quantum_number("n", 60)
127 2 : .set_quantum_number("l", 0)
128 2 : .set_quantum_number("j", 0.5)
129 2 : .set_quantum_number("m", 0.5)
130 1 : .create();
131 1 : auto basis_ion = BasisAtomCreator<double>(database).add_ket(ket_ion).create();
132 1 : DOCTEST_REQUIRE(basis_ion->get_number_of_states() == 1);
133 :
134 1 : std::array<double, 3> distance_vector{0, 0, 3 * UM_IN_ATOMIC_UNITS};
135 :
136 3 : for (int order : {2, 3}) {
137 : // Pair system of the atom and the ion
138 2 : auto basis_pair = BasisPairCreator<double>(basis_atom, basis_ion).create();
139 2 : DOCTEST_REQUIRE(basis_pair->get_number_of_states() == basis_atom->get_number_of_states());
140 :
141 2 : SystemPair<double> system_pair(basis_pair);
142 2 : system_pair.set_interaction_order(order);
143 2 : system_pair.set_distance_vector(distance_vector);
144 2 : system_pair.diagonalize(diagonalizer);
145 2 : Eigen::VectorXd energies_pair = system_pair.get_eigenenergies();
146 2 : energies_pair.array() -= ket_ion->get_energy();
147 :
148 : // Atom in the field of a classical point charge
149 2 : SystemAtom<double> system_reference(basis_atom);
150 2 : system_reference.set_ion_charge(1);
151 2 : system_reference.set_ion_interaction_order(order);
152 2 : system_reference.set_ion_distance_vector(distance_vector);
153 2 : system_reference.diagonalize(diagonalizer);
154 2 : Eigen::VectorXd energies_reference = system_reference.get_eigenenergies();
155 :
156 2 : DOCTEST_REQUIRE(energies_pair.size() == energies_reference.size());
157 2 : DOCTEST_CHECK((energies_pair - energies_reference).cwiseAbs().maxCoeff() < 1e-12);
158 :
159 : // The interaction must have a significant effect
160 2 : SystemAtom<double> system_unperturbed(basis_atom);
161 2 : system_unperturbed.diagonalize(diagonalizer);
162 2 : DOCTEST_CHECK((energies_pair - system_unperturbed.get_eigenenergies()).norm() *
163 : HARTREE_IN_GHZ >
164 : 1e-3);
165 2 : }
166 1 : }
167 :
168 1 : DOCTEST_TEST_CASE("ion ion pair interaction") {
169 : // Two ions in a single state each repel each other via the Coulomb interaction Z1*Z2/R
170 1 : auto &database = Database::get_global_instance();
171 :
172 1 : auto ket_ion = KetAtomCreator(database)
173 2 : .set_species("Sr88_ion")
174 2 : .set_quantum_number("n", 60)
175 2 : .set_quantum_number("l", 0)
176 2 : .set_quantum_number("j", 0.5)
177 2 : .set_quantum_number("m", 0.5)
178 1 : .create();
179 1 : auto basis_ion = BasisAtomCreator<double>(database).add_ket(ket_ion).create();
180 1 : auto basis_pair = BasisPairCreator<double>(basis_ion, basis_ion).create();
181 1 : DOCTEST_REQUIRE(basis_pair->get_number_of_states() == 1);
182 :
183 3 : for (double distance : {1 * UM_IN_ATOMIC_UNITS, 3 * UM_IN_ATOMIC_UNITS}) {
184 2 : SystemPair<double> system_pair(basis_pair);
185 2 : system_pair.set_interaction_order(3);
186 2 : system_pair.set_distance_vector({0, 0, distance});
187 2 : double energy = system_pair.get_matrix().coeff(0, 0) - 2 * ket_ion->get_energy();
188 2 : DOCTEST_CHECK(energy == doctest::Approx(1 / distance).epsilon(1e-6));
189 2 : }
190 1 : }
191 :
192 : #ifdef WITH_LAPACKE
193 1 : DOCTEST_TEST_CASE("diagonalize with lapacke_evr") {
194 1 : auto &database = Database::get_global_instance();
195 1 : auto diagonalizer = DiagonalizerLapackeEvr<double>();
196 :
197 : // Construct the state of interest
198 1 : auto ket = KetAtomCreator(database)
199 2 : .set_species("Rb")
200 2 : .set_quantum_number("n", 60)
201 2 : .set_quantum_number("l", 0)
202 2 : .set_quantum_number("j", 0.5)
203 2 : .set_quantum_number("m", 0.5)
204 1 : .create();
205 :
206 : // Construct and diagonalize the constituent systems
207 1 : auto basis = BasisAtomCreator<double>(database)
208 2 : .set_species("Rb")
209 2 : .restrict_quantum_number("n", ket->get_quantum_number("n") - 3,
210 2 : ket->get_quantum_number("n") + 3)
211 2 : .restrict_quantum_number("l", ket->get_quantum_number("l") - 1,
212 2 : ket->get_quantum_number("l") + 1)
213 1 : .create();
214 1 : SystemAtom<double> system(basis);
215 :
216 : // Construct and diagonalize the system_pair
217 2 : auto basis_pair = BasisPairCreator<double>(system.get_eigenbasis(), system.get_eigenbasis())
218 1 : .restrict_energy(2 * ket->get_energy() - 2 / HARTREE_IN_GHZ,
219 1 : 2 * ket->get_energy() + 2 / HARTREE_IN_GHZ)
220 1 : .create();
221 1 : DOCTEST_MESSAGE("Number of states in pair basis: ", basis_pair->get_number_of_states());
222 :
223 1 : auto system_pair = SystemPair<double>(basis_pair);
224 1 : system_pair.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
225 1 : system_pair.diagonalize(diagonalizer, 2 * ket->get_energy() - 0.5 / HARTREE_IN_GHZ,
226 1 : 2 * ket->get_energy() + 0.5 / HARTREE_IN_GHZ);
227 :
228 : // Print the largest and smallest eigenenergies
229 1 : auto eigenenergies = system_pair.get_eigenenergies();
230 1 : DOCTEST_MESSAGE("Lowest energy: ", eigenenergies.minCoeff());
231 1 : DOCTEST_MESSAGE("Highest energy: ", eigenenergies.maxCoeff());
232 1 : }
233 : #endif
234 : } // namespace pairinteraction
|
