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/database/Database.hpp"
8 : #include "pairinteraction/enums/OperatorType.hpp"
9 : #include "pairinteraction/enums/TransformationType.hpp"
10 : #include "pairinteraction/ket/KetAtom.hpp"
11 : #include "pairinteraction/utils/eigen_assertion.hpp"
12 : #include "pairinteraction/utils/eigen_compat.hpp"
13 : #include "pairinteraction/utils/operator.hpp"
14 : #include "pairinteraction/utils/spherical.hpp"
15 :
16 : #include <Eigen/Dense>
17 : #include <algorithm>
18 : #include <cmath>
19 : #include <limits>
20 : #include <memory>
21 : #include <set>
22 : #include <spdlog/spdlog.h>
23 :
24 : namespace pairinteraction {
25 : template <typename Scalar>
26 284 : SystemAtom<Scalar>::SystemAtom(std::shared_ptr<const basis_t> basis)
27 284 : : System<SystemAtom<Scalar>>(std::move(basis)) {}
28 :
29 : template <typename Scalar>
30 157 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_electric_field(const std::array<real_t, 3> &field) {
31 157 : this->hamiltonian_requires_construction = true;
32 :
33 102 : if (!traits::NumTraits<Scalar>::is_complex_v && field[1] != 0) {
34 0 : throw std::invalid_argument(
35 : "The field must not have a y-component if the scalar type is real.");
36 : }
37 :
38 157 : electric_field = field;
39 :
40 157 : return *this;
41 : }
42 :
43 : template <typename Scalar>
44 111 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_magnetic_field(const std::array<real_t, 3> &field) {
45 111 : this->hamiltonian_requires_construction = true;
46 :
47 77 : if (!traits::NumTraits<Scalar>::is_complex_v && field[1] != 0) {
48 0 : throw std::invalid_argument(
49 : "The field must not have a y-component if the scalar type is real.");
50 : }
51 :
52 111 : magnetic_field = field;
53 :
54 111 : return *this;
55 : }
56 :
57 : template <typename Scalar>
58 97 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_diamagnetism_enabled(bool enable) {
59 97 : this->hamiltonian_requires_construction = true;
60 97 : diamagnetism_enabled = enable;
61 97 : return *this;
62 : }
63 :
64 : template <typename Scalar>
65 : SystemAtom<Scalar> &
66 21 : SystemAtom<Scalar>::set_ion_distance_vector(const std::array<real_t, 3> &vector) {
67 21 : this->hamiltonian_requires_construction = true;
68 :
69 4 : if (!traits::NumTraits<Scalar>::is_complex_v && vector[1] != 0) {
70 0 : throw std::invalid_argument(
71 : "The distance vector must not have a y-component if the scalar type is real.");
72 : }
73 :
74 21 : ion_distance_vector = vector;
75 :
76 21 : return *this;
77 : }
78 :
79 : template <typename Scalar>
80 19 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_ion_charge(real_t charge) {
81 19 : this->hamiltonian_requires_construction = true;
82 19 : ion_charge = charge;
83 19 : return *this;
84 : }
85 :
86 : template <typename Scalar>
87 21 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_ion_interaction_order(int value) {
88 21 : this->hamiltonian_requires_construction = true;
89 21 : if (value < 2 || value > 3) {
90 0 : throw std::invalid_argument("The order of the Rydberg-ion interaction must be 2 or 3");
91 : }
92 21 : ion_interaction_order = value;
93 21 : return *this;
94 : }
95 :
96 : template <typename Scalar>
97 282 : void SystemAtom<Scalar>::construct_hamiltonian() const {
98 1299 : auto get_operator_matrix = [this](OperatorType type, int q = 0) {
99 339 : return this->basis->get_database().get_matrix_elements_in_canonical_basis(
100 339 : this->basis, this->basis, type, q);
101 : };
102 :
103 : // Construct the unperturbed Hamiltonian in the canonical atomic basis
104 282 : this->matrix = utils::get_energies_in_canonical_basis(this->basis);
105 :
106 282 : this->hamiltonian_is_diagonal = false;
107 282 : bool sort_by_quantum_number_f = this->basis->has_quantum_number_f();
108 282 : bool sort_by_quantum_number_m = this->basis->has_quantum_number_m();
109 282 : bool sort_by_parity = this->basis->has_parity();
110 :
111 : // Estimate the numerical precision so that we can decide which terms to keep
112 282 : Eigen::VectorX<real_t> diag = this->matrix.diagonal().real();
113 282 : real_t scale = (diag - diag.mean() * Eigen::VectorX<real_t>::Ones(diag.size())).norm();
114 282 : real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
115 :
116 282 : real_t typical_magnetic_dipole = 1e2; // ~n^1
117 282 : real_t typical_electric_dipole = 1e4; // ~n^2
118 282 : real_t typical_electric_quadrupole = 1e8; // ~n^4
119 :
120 : // Add the interaction with the field of an ion (see
121 : // https://en.wikipedia.org/wiki/Multipole_expansion#Spherical_form for details)
122 :
123 564 : Eigen::Map<const Eigen::Vector3<real_t>> vector_map(ion_distance_vector.data(),
124 282 : ion_distance_vector.size());
125 282 : real_t distance = vector_map.norm();
126 282 : SPDLOG_DEBUG("Distance to the ion: {}", distance);
127 :
128 : // Dipole order
129 282 : if (std::isfinite(distance) && ion_interaction_order >= 2) {
130 : // Calculate sqrt(4pi/3) * r * Y_1,q with q = -1, 0, 1 for the first three elements. Take
131 : // the conjugate of the result and scale it by 1/distance^2.
132 21 : Eigen::Vector3<Scalar> vector_dipole_order =
133 21 : spherical::get_transformator<Scalar>(1) * vector_map / distance;
134 21 : vector_dipole_order = vector_dipole_order.conjugate() / std::pow(distance, 2);
135 :
136 84 : for (int q = -1; q <= 1; ++q) {
137 63 : if (std::abs(vector_dipole_order[q + 1]) * typical_electric_dipole >
138 : numerical_precision) {
139 33 : this->matrix -= ion_charge * vector_dipole_order[q + 1] *
140 : get_operator_matrix(OperatorType::ELECTRIC_DIPOLE, q);
141 33 : sort_by_quantum_number_f = false;
142 33 : sort_by_parity = false;
143 33 : sort_by_quantum_number_m &= (q == 0);
144 : }
145 : }
146 : }
147 :
148 : // Quadrupole order (the last entry of vector_quadrupole_order would correspond to
149 : // ELECTRIC_QUADRUPOLE_ZERO and does not contribute to a *traceless* quadrupole)
150 282 : if (std::isfinite(distance) && ion_interaction_order >= 3) {
151 : // Calculate sqrt(4pi/5) * r^2 * Y_2,q / 3 with q = -2, -1, 0, 1, 2 for the first five
152 : // elements and (x^2 + y^2 + z^2) / 6 as the last element. Take the conjugate of the
153 : // result and scale it by 1/distance^3.
154 20 : Eigen::Vector<Scalar, 6> vector_quadrupole_order = spherical::get_transformator<Scalar>(2) *
155 20 : Eigen::KroneckerProduct(vector_map / distance, vector_map / distance);
156 20 : vector_quadrupole_order = 3 * vector_quadrupole_order.conjugate() / std::pow(distance, 3);
157 :
158 120 : for (int q = -2; q <= 2; ++q) {
159 100 : if (std::abs(vector_quadrupole_order[q + 2]) * typical_electric_quadrupole >
160 : numerical_precision) {
161 44 : this->matrix -= ion_charge * vector_quadrupole_order[q + 2] *
162 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, q);
163 44 : sort_by_quantum_number_f = false;
164 44 : sort_by_quantum_number_m &= (q == 0);
165 : }
166 : }
167 : }
168 :
169 : // Add external fields (see https://arxiv.org/abs/1612.08053 for details)
170 :
171 : // Transform the electric field to spherical coordinates and take the conjugate
172 282 : Eigen::Vector3<Scalar> electric_field_spherical = spherical::get_transformator<Scalar>(1) *
173 282 : Eigen::Map<const Eigen::Vector3<real_t>>(electric_field.data(), electric_field.size());
174 282 : electric_field_spherical = electric_field_spherical.conjugate();
175 :
176 : // Transform the magnetic field to spherical coordinates and take the conjugate
177 282 : Eigen::Vector3<Scalar> magnetic_field_spherical = spherical::get_transformator<Scalar>(1) *
178 282 : Eigen::Map<const Eigen::Vector3<real_t>>(magnetic_field.data(), magnetic_field.size());
179 282 : magnetic_field_spherical = magnetic_field_spherical.conjugate();
180 :
181 : // Stark effect: - \vec{d} \vec{E} = - d_{1,0} E_{0} + d_{1,1} E_{-} + d_{1,-1} E_{+}
182 : // = - d_{1,0} E_{0} - d_{1,1} E_{+}^* - d_{1,-1} E_{-}^*
183 : // with the electric dipole operator: d_{1,q} = - e r sqrt{4 pi / 3} Y_{1,q}(\theta, \phi)
184 : // where electric_field_spherical=[E_{-}^*, E_{0}, E_{+}^*]
185 1128 : for (int q = -1; q <= 1; ++q) {
186 846 : if (std::abs(electric_field_spherical[q + 1]) * typical_electric_dipole >
187 : numerical_precision) {
188 183 : this->matrix -= electric_field_spherical[q + 1] *
189 : get_operator_matrix(OperatorType::ELECTRIC_DIPOLE, q);
190 183 : sort_by_quantum_number_f = false;
191 183 : sort_by_parity = false;
192 183 : sort_by_quantum_number_m &= (q == 0);
193 : }
194 : }
195 :
196 : // Zeeman effect: - \vec{\mu} \vec{B} = - \mu_{1,0} B_{0} + \mu_{1,1} B_{-} + \mu_{1,-1} B_{+}
197 : // = - \mu_{1,0} B_{0} - \mu_{1,1} B_{+}^* - \mu_{1,-1} B_{-}^*
198 : // with the magnetic dipole operator: \vec{\mu} = - \mu_B / \hbar (g_l \vec{l} + g_s \vec{s})
199 : // where magnetic_field_spherical=[B_{-}^*, B_{0}, B_{+}^*]
200 1128 : for (int q = -1; q <= 1; ++q) {
201 846 : if (std::abs(magnetic_field_spherical[q + 1]) * typical_magnetic_dipole >
202 : numerical_precision) {
203 57 : this->matrix -= magnetic_field_spherical[q + 1] *
204 : get_operator_matrix(OperatorType::MAGNETIC_DIPOLE, q);
205 57 : sort_by_quantum_number_f = false;
206 57 : sort_by_quantum_number_m &= (q == 0);
207 : }
208 : }
209 :
210 : // Diamagnetism: 1 / (8 m_e) abs(\vec{d} \times \vec{B})^2
211 : // = (1/12) \left[ B_0^2 (q0 - d_{2,0}) + B_+ B_- (-2 q0 - d_{2,0})
212 : // + \sqrt{3} B_0 B_- d_{2,1} + \sqrt{3} B_0 B_+ d_{2,-1}
213 : // - \sqrt{3/2} B_-^2 d_{2,2} - \sqrt{3/2} B_+^2 d_{2,-2} \right]
214 : // with the operator: q0 = e^2 r^2 sqrt{4 pi} Y_{0,0} = e^2 r^2
215 : // and the electric quadrupole operator: d_{2,q} = e^2 r^2 sqrt{4 pi / 5} Y_{2,q}(\theta, \phi)
216 : // where magnetic_field_spherical=[B_{-}^*, B_{0}, B_{+}^*]
217 282 : if (diamagnetism_enabled) {
218 59 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
219 : numerical_precision) {
220 14 : this->matrix += static_cast<real_t>(1 / 12.) *
221 11 : static_cast<Scalar>(std::pow(magnetic_field_spherical[1], 2)) *
222 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0);
223 14 : this->matrix -= static_cast<real_t>(1 / 12.) *
224 11 : static_cast<Scalar>(std::pow(magnetic_field_spherical[1], 2)) *
225 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, 0);
226 7 : sort_by_quantum_number_f = false;
227 : }
228 59 : if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
229 61 : numerical_precision &&
230 2 : std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
231 : numerical_precision) {
232 4 : this->matrix -= static_cast<real_t>(2 / 12.) * magnetic_field_spherical[0] *
233 2 : magnetic_field_spherical[2] *
234 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0);
235 4 : this->matrix -= static_cast<real_t>(1 / 12.) * magnetic_field_spherical[0] *
236 2 : magnetic_field_spherical[2] *
237 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, 0);
238 2 : sort_by_quantum_number_f = false;
239 : }
240 59 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
241 66 : numerical_precision &&
242 7 : std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
243 : numerical_precision) {
244 0 : this->matrix -= static_cast<real_t>(std::sqrt(3.0) / 12.) *
245 0 : magnetic_field_spherical[1] * magnetic_field_spherical[2] *
246 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, 1);
247 0 : sort_by_quantum_number_f = false;
248 0 : sort_by_quantum_number_m = false;
249 : }
250 59 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
251 66 : numerical_precision &&
252 7 : std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
253 : numerical_precision) {
254 0 : this->matrix -= static_cast<real_t>(std::sqrt(3.0) / 12.) *
255 0 : magnetic_field_spherical[1] * magnetic_field_spherical[0] *
256 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, -1);
257 0 : sort_by_quantum_number_f = false;
258 0 : sort_by_quantum_number_m = false;
259 : }
260 59 : if (std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
261 : numerical_precision) {
262 4 : this->matrix -= static_cast<real_t>(std::sqrt(1.5) / 12.) *
263 4 : static_cast<Scalar>(std::pow(magnetic_field_spherical[2], 2)) *
264 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, 2);
265 2 : sort_by_quantum_number_f = false;
266 2 : sort_by_quantum_number_m = false;
267 : }
268 59 : if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
269 : numerical_precision) {
270 4 : this->matrix -= static_cast<real_t>(std::sqrt(1.5) / 12.) *
271 4 : static_cast<Scalar>(std::pow(magnetic_field_spherical[0], 2)) *
272 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, -2);
273 2 : sort_by_quantum_number_f = false;
274 2 : sort_by_quantum_number_m = false;
275 : }
276 : }
277 :
278 : // Transform from the canonical basis into the actual basis
279 282 : this->matrix =
280 282 : this->basis->get_coefficients().adjoint() * this->matrix * this->basis->get_coefficients();
281 :
282 : // Store which labels can be used to block-diagonalize the Hamiltonian
283 282 : this->blockdiagonalizing_labels.clear();
284 282 : if (sort_by_quantum_number_f) {
285 96 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_F);
286 : }
287 282 : if (sort_by_quantum_number_m) {
288 236 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_M);
289 : }
290 282 : if (sort_by_parity) {
291 142 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_PARITY);
292 : }
293 282 : }
294 :
295 : // Explicit instantiations
296 : template class SystemAtom<double>;
297 : template class SystemAtom<std::complex<double>>;
298 : } // namespace pairinteraction
|
