Line data Source code
1 : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
2 : // SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : #include "pairinteraction/basis/BasisPairCreator.hpp"
5 :
6 : #include "pairinteraction/basis/BasisAtom.hpp"
7 : #include "pairinteraction/basis/BasisAtomCreator.hpp"
8 : #include "pairinteraction/basis/BasisPair.hpp"
9 : #include "pairinteraction/database/Database.hpp"
10 : #include "pairinteraction/diagonalize/DiagonalizerEigen.hpp"
11 : #include "pairinteraction/enums/OperatorType.hpp"
12 : #include "pairinteraction/enums/Parity.hpp"
13 : #include "pairinteraction/enums/TransformationType.hpp"
14 : #include "pairinteraction/ket/KetAtom.hpp"
15 : #include "pairinteraction/ket/KetAtomCreator.hpp"
16 : #include "pairinteraction/ket/KetPair.hpp"
17 : #include "pairinteraction/system/SystemAtom.hpp"
18 : #include "pairinteraction/system/SystemPair.hpp"
19 : #include "pairinteraction/utils/streamed.hpp"
20 :
21 : #include <doctest/doctest.h>
22 :
23 : namespace pairinteraction {
24 :
25 : constexpr double HARTREE_IN_GHZ = 6579683.920501762;
26 : constexpr double VOLT_PER_CM_IN_ATOMIC_UNITS = 1 / 5.14220675112e9;
27 : constexpr double UM_IN_ATOMIC_UNITS = 1 / 5.29177210544e-5;
28 :
29 3 : DOCTEST_TEST_CASE("create a BasisPair") {
30 : // Create single-atom system
31 3 : Database &database = Database::get_global_instance();
32 3 : auto basis = BasisAtomCreator<double>()
33 6 : .set_species("Rb")
34 3 : .restrict_quantum_number_n(58, 62)
35 3 : .restrict_quantum_number_l(0, 2)
36 3 : .create(database);
37 3 : SystemAtom<double> system(basis);
38 3 : system.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
39 :
40 3 : DiagonalizerEigen<double> diagonalizer;
41 3 : system.diagonalize(diagonalizer);
42 :
43 : // Get energy window for a two-atom basis
44 3 : auto ket = KetAtomCreator()
45 6 : .set_species("Rb")
46 3 : .set_quantum_number_n(60)
47 3 : .set_quantum_number_l(0)
48 3 : .set_quantum_number_m(0.5)
49 3 : .create(database);
50 3 : double min_energy = 2 * ket->get_energy() - 3 / HARTREE_IN_GHZ;
51 3 : double max_energy = 2 * ket->get_energy() + 3 / HARTREE_IN_GHZ;
52 :
53 : // Create two-atom bases
54 6 : auto basis_pair_a = pairinteraction::BasisPairCreator<double>()
55 3 : .add(system)
56 3 : .add(system)
57 3 : .restrict_energy(min_energy, max_energy)
58 3 : .restrict_quantum_number_m(1, 1)
59 3 : .create();
60 6 : auto basis_pair_b = pairinteraction::BasisPairCreator<double>()
61 3 : .add(system)
62 3 : .add(system)
63 3 : .restrict_energy(min_energy, max_energy)
64 3 : .restrict_quantum_number_m(1, 1)
65 3 : .create();
66 :
67 3 : DOCTEST_SUBCASE("check equality of kets") {
68 : // Obtain kets from the two-atom bases and check for equality
69 1 : auto ket1a = basis_pair_a->get_kets()[0];
70 1 : auto ket1b = basis_pair_b->get_kets()[0];
71 1 : auto ket2a = basis_pair_a->get_kets()[1];
72 1 : auto ket2b = basis_pair_b->get_kets()[1];
73 1 : DOCTEST_CHECK(*ket1a == *ket1a);
74 1 : DOCTEST_CHECK(*ket2a == *ket2a);
75 1 : DOCTEST_CHECK(*ket1a != *ket2b);
76 1 : DOCTEST_CHECK(*ket2a != *ket1b);
77 :
78 : // Currently, kets from different BasisPair are never equal
79 1 : DOCTEST_CHECK(*ket1a != *ket1b);
80 1 : DOCTEST_CHECK(*ket2a != *ket2b);
81 4 : }
82 :
83 3 : DOCTEST_SUBCASE("check overlap") {
84 1 : auto overlaps = basis_pair_a->get_overlaps(ket, ket);
85 :
86 : // The total overlap is less than 1 because of the restricted energy window
87 1 : DOCTEST_CHECK(overlaps.sum() == doctest::Approx(0.9107819201));
88 4 : }
89 :
90 3 : DOCTEST_SUBCASE("get the atomic states constituting a ket of the basis_pair") {
91 1 : auto atomic_states = basis_pair_a->get_kets()[0]->get_atomic_states();
92 1 : DOCTEST_CHECK(atomic_states.size() == 2);
93 1 : DOCTEST_CHECK(atomic_states[0]->get_number_of_states() == 1);
94 1 : DOCTEST_CHECK(atomic_states[0]->get_number_of_kets() == basis->get_number_of_kets());
95 4 : }
96 3 : }
97 :
98 2 : DOCTEST_TEST_CASE("get matrix elements in the pair basis") {
99 2 : DiagonalizerEigen<double> diagonalizer;
100 :
101 : // Create single-atom system
102 2 : Database &database = Database::get_global_instance();
103 2 : auto basis = BasisAtomCreator<double>()
104 4 : .set_species("Rb")
105 2 : .restrict_quantum_number_n(58, 62)
106 2 : .restrict_quantum_number_l(0, 2)
107 2 : .create(database);
108 2 : SystemAtom<double> system(basis);
109 2 : system.set_electric_field({0, 0, 10 * VOLT_PER_CM_IN_ATOMIC_UNITS});
110 2 : system.diagonalize(diagonalizer);
111 :
112 : // Get energy window for a two-atom basis
113 2 : auto ket = KetAtomCreator()
114 4 : .set_species("Rb")
115 2 : .set_quantum_number_n(60)
116 2 : .set_quantum_number_l(0)
117 2 : .set_quantum_number_m(0.5)
118 2 : .create(database);
119 2 : double min_energy = 2 * ket->get_energy() - 3 / HARTREE_IN_GHZ;
120 2 : double max_energy = 2 * ket->get_energy() + 3 / HARTREE_IN_GHZ;
121 :
122 : // Create two-atom system
123 4 : auto basis_pair_unperturbed = pairinteraction::BasisPairCreator<double>()
124 2 : .add(system)
125 2 : .add(system)
126 2 : .restrict_energy(min_energy, max_energy)
127 2 : .restrict_quantum_number_m(1, 1)
128 2 : .create();
129 4 : auto system_pair = SystemPair<double>(basis_pair_unperturbed)
130 2 : .set_distance_vector({0, 0, 1 * UM_IN_ATOMIC_UNITS});
131 2 : system_pair.diagonalize(diagonalizer);
132 :
133 2 : auto basis_pair = system_pair.get_eigenbasis();
134 :
135 2 : DOCTEST_SUBCASE("check dimensions") {
136 : // <basis_pair_unperturbed|d0d0|basis_pair_unperturbed>
137 : auto matrix_elements_all = basis_pair_unperturbed->get_matrix_elements(
138 : basis_pair_unperturbed, OperatorType::ELECTRIC_DIPOLE, OperatorType::ELECTRIC_DIPOLE, 0,
139 1 : 0);
140 1 : DOCTEST_CHECK(matrix_elements_all.rows() == basis_pair_unperturbed->get_number_of_states());
141 1 : DOCTEST_CHECK(matrix_elements_all.cols() == basis_pair_unperturbed->get_number_of_states());
142 :
143 : // <ket_pair|d0d0|basis_pair_unperturbed>
144 1 : auto ket_pair = basis_pair_unperturbed->get_kets()[0];
145 : auto matrix_elements_ket_pair = basis_pair_unperturbed->get_matrix_elements(
146 1 : ket_pair, OperatorType::ELECTRIC_DIPOLE, OperatorType::ELECTRIC_DIPOLE, 0, 0);
147 1 : DOCTEST_CHECK(matrix_elements_ket_pair.size() ==
148 : basis_pair_unperturbed->get_number_of_states());
149 :
150 : {
151 1 : Eigen::VectorX<double> ref = matrix_elements_all.row(0);
152 1 : DOCTEST_CHECK(ref.isApprox(matrix_elements_ket_pair, 1e-11));
153 1 : }
154 :
155 : // <basis x basis|d0d0|basis_pair>
156 : auto matrix_elements_product = basis_pair->get_matrix_elements(
157 1 : basis, basis, OperatorType::ELECTRIC_DIPOLE, OperatorType::ELECTRIC_DIPOLE, 0, 0);
158 1 : DOCTEST_CHECK(matrix_elements_product.rows() ==
159 : basis->get_number_of_states() * basis->get_number_of_states());
160 1 : DOCTEST_CHECK(matrix_elements_product.cols() == basis_pair->get_number_of_states());
161 :
162 : // <ket,ket|d0d0|basis_pair>
163 : auto matrix_elements_ket = basis_pair->get_matrix_elements(
164 1 : ket, ket, OperatorType::ELECTRIC_DIPOLE, OperatorType::ELECTRIC_DIPOLE, 0, 0);
165 1 : DOCTEST_CHECK(matrix_elements_ket.size() == basis_pair->get_number_of_states());
166 3 : }
167 :
168 2 : DOCTEST_SUBCASE("check matrix elements") {
169 : // energy
170 : auto hamiltonian = basis_pair->get_matrix_elements(basis_pair, OperatorType::ENERGY,
171 1 : OperatorType::IDENTITY);
172 2 : hamiltonian += basis_pair->get_matrix_elements(basis_pair, OperatorType::IDENTITY,
173 1 : OperatorType::ENERGY);
174 :
175 : // interaction with electric field
176 : {
177 2 : Eigen::SparseMatrix<double, Eigen::RowMajor> tmp = -basis_pair->get_matrix_elements(
178 1 : basis_pair, OperatorType::ELECTRIC_DIPOLE, OperatorType::IDENTITY, 0, 0);
179 2 : tmp += -basis_pair->get_matrix_elements(basis_pair, OperatorType::IDENTITY,
180 1 : OperatorType::ELECTRIC_DIPOLE, 0, 0);
181 1 : hamiltonian += 10 * VOLT_PER_CM_IN_ATOMIC_UNITS * tmp;
182 1 : }
183 :
184 : // dipole-dipole interaction
185 : {
186 1 : Eigen::SparseMatrix<double, Eigen::RowMajor> tmp = -2 *
187 2 : basis_pair->get_matrix_elements(basis_pair, OperatorType::ELECTRIC_DIPOLE,
188 1 : OperatorType::ELECTRIC_DIPOLE, 0, 0);
189 2 : tmp += -basis_pair->get_matrix_elements(basis_pair, OperatorType::ELECTRIC_DIPOLE,
190 1 : OperatorType::ELECTRIC_DIPOLE, 1, -1);
191 2 : tmp += -basis_pair->get_matrix_elements(basis_pair, OperatorType::ELECTRIC_DIPOLE,
192 1 : OperatorType::ELECTRIC_DIPOLE, -1, 1);
193 1 : hamiltonian += std::pow(UM_IN_ATOMIC_UNITS, -3) * tmp;
194 1 : }
195 :
196 : // compare to reference
197 1 : const auto &ref = system_pair.get_matrix();
198 1 : DOCTEST_CHECK(ref.isApprox(hamiltonian, 1e-11));
199 3 : }
200 2 : }
201 :
202 : } // namespace pairinteraction
|
