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/system/GreenTensorInterpolator.hpp"
12 : #include "pairinteraction/utils/eigen_assertion.hpp"
13 : #include "pairinteraction/utils/eigen_compat.hpp"
14 : #include "pairinteraction/utils/operator.hpp"
15 : #include "pairinteraction/utils/spherical.hpp"
16 :
17 : #include <Eigen/Dense>
18 : #include <Eigen/SparseCore>
19 : #include <algorithm>
20 : #include <cmath>
21 : #include <initializer_list>
22 : #include <iterator>
23 : #include <limits>
24 : #include <memory>
25 : #include <set>
26 : #include <spdlog/spdlog.h>
27 : #include <variant>
28 : #include <vector>
29 :
30 : namespace pairinteraction {
31 : template <typename Scalar>
32 663 : SystemAtom<Scalar>::SystemAtom(std::shared_ptr<const basis_t> basis)
33 663 : : System<SystemAtom<Scalar>>(std::move(basis)) {}
34 :
35 : template <typename Scalar>
36 160 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_electric_field(const std::array<real_t, 3> &field) {
37 160 : this->hamiltonian_requires_construction = true;
38 :
39 104 : if (!traits::NumTraits<Scalar>::is_complex_v && field[1] != 0) {
40 0 : throw std::invalid_argument(
41 : "The field must not have a y-component if the scalar type is real.");
42 : }
43 :
44 160 : electric_field = field;
45 :
46 160 : return *this;
47 : }
48 :
49 : template <typename Scalar>
50 110 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_magnetic_field(const std::array<real_t, 3> &field) {
51 110 : this->hamiltonian_requires_construction = true;
52 :
53 76 : if (!traits::NumTraits<Scalar>::is_complex_v && field[1] != 0) {
54 0 : throw std::invalid_argument(
55 : "The field must not have a y-component if the scalar type is real.");
56 : }
57 :
58 110 : magnetic_field = field;
59 :
60 110 : return *this;
61 : }
62 :
63 : template <typename Scalar>
64 97 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_diamagnetism_enabled(bool enable) {
65 97 : this->hamiltonian_requires_construction = true;
66 97 : diamagnetism_enabled = enable;
67 97 : return *this;
68 : }
69 :
70 : template <typename Scalar>
71 : SystemAtom<Scalar> &
72 21 : SystemAtom<Scalar>::set_ion_distance_vector(const std::array<real_t, 3> &vector) {
73 21 : this->hamiltonian_requires_construction = true;
74 :
75 4 : if (!traits::NumTraits<Scalar>::is_complex_v && vector[1] != 0) {
76 0 : throw std::invalid_argument(
77 : "The distance vector must not have a y-component if the scalar type is real.");
78 : }
79 :
80 21 : ion_distance_vector = vector;
81 :
82 21 : return *this;
83 : }
84 :
85 : template <typename Scalar>
86 19 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_ion_charge(real_t charge) {
87 19 : this->hamiltonian_requires_construction = true;
88 19 : ion_charge = charge;
89 19 : return *this;
90 : }
91 :
92 : template <typename Scalar>
93 21 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_ion_interaction_order(int value) {
94 21 : this->hamiltonian_requires_construction = true;
95 21 : if (value < 2 || value > 3) {
96 0 : throw std::invalid_argument("The order of the Rydberg-ion interaction must be 2 or 3");
97 : }
98 21 : ion_interaction_order = value;
99 21 : return *this;
100 : }
101 :
102 : template <typename Scalar>
103 2 : SystemAtom<Scalar> &SystemAtom<Scalar>::set_green_tensor_interpolator(
104 : const std::shared_ptr<const GreenTensorInterpolator<Scalar>> &green_tensor_interpolator) {
105 2 : this->hamiltonian_requires_construction = true;
106 2 : this->green_tensor_interpolator = green_tensor_interpolator;
107 2 : return *this;
108 : }
109 :
110 : template <typename Scalar>
111 663 : void SystemAtom<Scalar>::construct_hamiltonian() const {
112 1686 : auto get_operator_matrix = [this](OperatorType type, int q = 0) {
113 341 : return this->basis->get_database().get_matrix_elements_in_canonical_basis(
114 341 : this->basis, this->basis, type, q);
115 : };
116 : // Helper function for constructing matrices of spherical harmonics operators in the
117 : // canonical basis. The full Hamiltonian is transformed into the actual basis at the end.
118 4 : auto get_matrices = [](auto basis, OperatorType type, std::initializer_list<int> m,
119 : bool conjugate) {
120 4 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> matrices;
121 4 : matrices.reserve(m.size());
122 4 : int factor = conjugate ? -1 : 1;
123 52 : std::transform(m.begin(), m.end(), std::back_inserter(matrices), [&](int q) {
124 24 : return (std::pow(factor, q) *
125 : basis->get_database().get_matrix_elements_in_canonical_basis(basis, basis, type,
126 : factor * q))
127 36 : .eval();
128 : });
129 8 : return matrices;
130 0 : };
131 :
132 : // Construct the unperturbed Hamiltonian in the canonical atomic basis
133 663 : this->matrix = utils::get_energies_in_canonical_basis(this->basis);
134 :
135 663 : this->hamiltonian_is_diagonal = false;
136 663 : bool sort_by_quantum_number_f = this->basis->has_quantum_number_f();
137 663 : bool sort_by_quantum_number_m = this->basis->has_quantum_number_m();
138 663 : bool sort_by_parity = this->basis->has_parity();
139 :
140 : // Estimate the numerical precision so that we can decide which terms to keep
141 663 : Eigen::VectorX<real_t> diag = this->matrix.diagonal().real();
142 663 : real_t scale = (diag - diag.mean() * Eigen::VectorX<real_t>::Ones(diag.size())).norm();
143 663 : real_t numerical_precision = 100 * scale * std::numeric_limits<real_t>::epsilon();
144 :
145 663 : real_t typical_magnetic_dipole = 1e2; // ~n^1
146 663 : real_t typical_electric_dipole = 1e4; // ~n^2
147 663 : real_t typical_electric_quadrupole = 1e8; // ~n^4
148 :
149 : // Add the interaction with the field of an ion (see
150 : // https://en.wikipedia.org/wiki/Multipole_expansion#Spherical_form for details)
151 :
152 1326 : Eigen::Map<const Eigen::Vector3<real_t>> vector_map(ion_distance_vector.data(),
153 663 : ion_distance_vector.size());
154 663 : real_t distance = vector_map.norm();
155 663 : SPDLOG_DEBUG("Distance to the ion: {}", distance);
156 :
157 : // Dipole order
158 663 : if (std::isfinite(distance) && ion_interaction_order >= 2) {
159 : // Calculate sqrt(4pi/3) * r * Y_1,q with q = -1, 0, 1 for the first three elements. Take
160 : // the conjugate of the result and scale it by 1/distance^2.
161 21 : Eigen::Vector3<Scalar> vector_dipole_order =
162 21 : spherical::get_transformator<Scalar>(1) * vector_map / distance;
163 21 : vector_dipole_order = vector_dipole_order.conjugate() / std::pow(distance, 2);
164 :
165 84 : for (int q = -1; q <= 1; ++q) {
166 63 : if (std::abs(vector_dipole_order[q + 1]) * typical_electric_dipole >
167 : numerical_precision) {
168 33 : this->matrix -= ion_charge * vector_dipole_order[q + 1] *
169 : get_operator_matrix(OperatorType::ELECTRIC_DIPOLE, q);
170 33 : sort_by_quantum_number_f = false;
171 33 : sort_by_parity = false;
172 33 : sort_by_quantum_number_m &= (q == 0);
173 : }
174 : }
175 : }
176 :
177 : // Quadrupole order (the last entry of vector_quadrupole_order would correspond to
178 : // ELECTRIC_QUADRUPOLE_ZERO and does not contribute to a *traceless* quadrupole)
179 663 : if (std::isfinite(distance) && ion_interaction_order >= 3) {
180 : // Calculate sqrt(4pi/5) * r^2 * Y_2,q / 3 with q = -2, -1, 0, 1, 2 for the first five
181 : // elements and (x^2 + y^2 + z^2) / 6 as the last element. Take the conjugate of the
182 : // result and scale it by 1/distance^3.
183 20 : Eigen::Vector<Scalar, 6> vector_quadrupole_order = spherical::get_transformator<Scalar>(2) *
184 20 : Eigen::KroneckerProduct(vector_map / distance, vector_map / distance);
185 20 : vector_quadrupole_order = 3 * vector_quadrupole_order.conjugate() / std::pow(distance, 3);
186 :
187 120 : for (int q = -2; q <= 2; ++q) {
188 100 : if (std::abs(vector_quadrupole_order[q + 2]) * typical_electric_quadrupole >
189 : numerical_precision) {
190 44 : this->matrix -= ion_charge * vector_quadrupole_order[q + 2] *
191 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, q);
192 44 : sort_by_quantum_number_f = false;
193 44 : sort_by_quantum_number_m &= (q == 0);
194 : }
195 : }
196 : }
197 :
198 : // Add external fields (see https://arxiv.org/abs/1612.08053 for details)
199 :
200 : // Transform the electric field to spherical coordinates and take the conjugate
201 663 : Eigen::Vector3<Scalar> electric_field_spherical = spherical::get_transformator<Scalar>(1) *
202 663 : Eigen::Map<const Eigen::Vector3<real_t>>(electric_field.data(), electric_field.size());
203 663 : electric_field_spherical = electric_field_spherical.conjugate();
204 :
205 : // Transform the magnetic field to spherical coordinates and take the conjugate
206 663 : Eigen::Vector3<Scalar> magnetic_field_spherical = spherical::get_transformator<Scalar>(1) *
207 663 : Eigen::Map<const Eigen::Vector3<real_t>>(magnetic_field.data(), magnetic_field.size());
208 663 : magnetic_field_spherical = magnetic_field_spherical.conjugate();
209 :
210 : // Stark effect: - \vec{d} \vec{E} = - d_{1,0} E_{0} + d_{1,1} E_{-} + d_{1,-1} E_{+}
211 : // = - d_{1,0} E_{0} - d_{1,1} E_{+}^* - d_{1,-1} E_{-}^*
212 : // with the electric dipole operator: d_{1,q} = - e r sqrt{4 pi / 3} Y_{1,q}(\theta, \phi)
213 : // where electric_field_spherical=[E_{-}^*, E_{0}, E_{+}^*]
214 2652 : for (int q = -1; q <= 1; ++q) {
215 1989 : if (std::abs(electric_field_spherical[q + 1]) * typical_electric_dipole >
216 : numerical_precision) {
217 186 : this->matrix -= electric_field_spherical[q + 1] *
218 : get_operator_matrix(OperatorType::ELECTRIC_DIPOLE, q);
219 186 : sort_by_quantum_number_f = false;
220 186 : sort_by_parity = false;
221 186 : sort_by_quantum_number_m &= (q == 0);
222 : }
223 : }
224 :
225 : // Zeeman effect: - \vec{\mu} \vec{B} = - \mu_{1,0} B_{0} + \mu_{1,1} B_{-} + \mu_{1,-1} B_{+}
226 : // = - \mu_{1,0} B_{0} - \mu_{1,1} B_{+}^* - \mu_{1,-1} B_{-}^*
227 : // with the magnetic dipole operator: \vec{\mu} = - \mu_B / \hbar (g_l \vec{l} + g_s \vec{s})
228 : // where magnetic_field_spherical=[B_{-}^*, B_{0}, B_{+}^*]
229 2652 : for (int q = -1; q <= 1; ++q) {
230 1989 : if (std::abs(magnetic_field_spherical[q + 1]) * typical_magnetic_dipole >
231 : numerical_precision) {
232 56 : this->matrix -= magnetic_field_spherical[q + 1] *
233 : get_operator_matrix(OperatorType::MAGNETIC_DIPOLE, q);
234 56 : sort_by_quantum_number_f = false;
235 56 : sort_by_quantum_number_m &= (q == 0);
236 : }
237 : }
238 :
239 : // Diamagnetism: 1 / (8 m_e) abs(\vec{d} \times \vec{B})^2
240 : // = (1/12) \left[ B_0^2 (q0 - d_{2,0}) + B_+ B_- (-2 q0 - d_{2,0})
241 : // + \sqrt{3} B_0 B_- d_{2,1} + \sqrt{3} B_0 B_+ d_{2,-1}
242 : // - \sqrt{3/2} B_-^2 d_{2,2} - \sqrt{3/2} B_+^2 d_{2,-2} \right]
243 : // with the operator: q0 = e^2 r^2 sqrt{4 pi} Y_{0,0} = e^2 r^2
244 : // and the electric quadrupole operator: d_{2,q} = e^2 r^2 sqrt{4 pi / 5} Y_{2,q}(\theta, \phi)
245 : // where magnetic_field_spherical=[B_{-}^*, B_{0}, B_{+}^*]
246 663 : if (diamagnetism_enabled) {
247 59 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
248 : numerical_precision) {
249 14 : this->matrix += static_cast<real_t>(1 / 12.) *
250 11 : static_cast<Scalar>(std::pow(magnetic_field_spherical[1], 2)) *
251 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0);
252 14 : this->matrix -= static_cast<real_t>(1 / 12.) *
253 11 : static_cast<Scalar>(std::pow(magnetic_field_spherical[1], 2)) *
254 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, 0);
255 7 : sort_by_quantum_number_f = false;
256 : }
257 59 : if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
258 61 : numerical_precision &&
259 2 : std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
260 : numerical_precision) {
261 4 : this->matrix -= static_cast<real_t>(2 / 12.) * magnetic_field_spherical[0] *
262 2 : magnetic_field_spherical[2] *
263 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0);
264 4 : this->matrix -= static_cast<real_t>(1 / 12.) * magnetic_field_spherical[0] *
265 2 : magnetic_field_spherical[2] *
266 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, 0);
267 2 : sort_by_quantum_number_f = false;
268 : }
269 59 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
270 66 : numerical_precision &&
271 7 : std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
272 : numerical_precision) {
273 0 : this->matrix -= static_cast<real_t>(std::sqrt(3.0) / 12.) *
274 0 : magnetic_field_spherical[1] * magnetic_field_spherical[2] *
275 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, 1);
276 0 : sort_by_quantum_number_f = false;
277 0 : sort_by_quantum_number_m = false;
278 : }
279 59 : if (std::abs(magnetic_field_spherical[1]) * typical_electric_quadrupole >
280 66 : numerical_precision &&
281 7 : std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
282 : numerical_precision) {
283 0 : this->matrix -= static_cast<real_t>(std::sqrt(3.0) / 12.) *
284 0 : magnetic_field_spherical[1] * magnetic_field_spherical[0] *
285 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, -1);
286 0 : sort_by_quantum_number_f = false;
287 0 : sort_by_quantum_number_m = false;
288 : }
289 59 : if (std::abs(magnetic_field_spherical[2]) * typical_electric_quadrupole >
290 : numerical_precision) {
291 4 : this->matrix -= static_cast<real_t>(std::sqrt(1.5) / 12.) *
292 4 : static_cast<Scalar>(std::pow(magnetic_field_spherical[2], 2)) *
293 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, 2);
294 2 : sort_by_quantum_number_f = false;
295 2 : sort_by_quantum_number_m = false;
296 : }
297 59 : if (std::abs(magnetic_field_spherical[0]) * typical_electric_quadrupole >
298 : numerical_precision) {
299 4 : this->matrix -= static_cast<real_t>(std::sqrt(1.5) / 12.) *
300 4 : static_cast<Scalar>(std::pow(magnetic_field_spherical[0], 2)) *
301 : get_operator_matrix(OperatorType::ELECTRIC_QUADRUPOLE, -2);
302 2 : sort_by_quantum_number_f = false;
303 2 : sort_by_quantum_number_m = false;
304 : }
305 : }
306 :
307 : // Add self-interaction via Green tensor
308 : // H_self = (1/2) * Σ_{ij} D_left[i] * G_{ij} * D_right[j]
309 : // where D_left uses conjugated convention and
310 : // D_right uses normal convention.
311 :
312 : // Helper function for adding self interaction
313 665 : auto add_interaction = [this, &sort_by_quantum_number_f,
314 : &sort_by_quantum_number_m](const auto &entries, const auto &op_left,
315 : const auto &op_right, int delta) {
316 4 : for (const auto &entry : entries) {
317 2 : if (std::holds_alternative<
318 2 : typename GreenTensorInterpolator<Scalar>::OmegaDependentEntry>(entry)) {
319 0 : throw std::logic_error(
320 : "Green tensor with omega dependent entries is currently not supported.");
321 : }
322 :
323 : const auto &constant_entry =
324 2 : std::get<typename GreenTensorInterpolator<Scalar>::ConstantEntry>(entry);
325 2 : Eigen::SparseMatrix<Scalar, Eigen::RowMajor> self_interaction =
326 2 : (constant_entry.val() / 2. * op_left[constant_entry.row()] *
327 2 : op_right[constant_entry.col()])
328 : .eval();
329 2 : this->matrix += self_interaction;
330 :
331 2 : sort_by_quantum_number_f = false;
332 2 : if (constant_entry.row() != constant_entry.col() + delta) {
333 0 : sort_by_quantum_number_m = false;
334 : }
335 : }
336 : };
337 :
338 663 : if (green_tensor_interpolator) {
339 2 : if (!green_tensor_interpolator->get_spherical_entries(1, 1).empty()) {
340 2 : auto dipole_left =
341 2 : get_matrices(this->basis, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, true);
342 2 : auto dipole_right =
343 2 : get_matrices(this->basis, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, false);
344 :
345 2 : add_interaction(green_tensor_interpolator->get_spherical_entries(1, 1), dipole_left,
346 : dipole_right, 0);
347 2 : }
348 : }
349 :
350 : // Transform from the canonical basis into the actual basis
351 663 : this->matrix =
352 663 : this->basis->get_coefficients().adjoint() * this->matrix * this->basis->get_coefficients();
353 :
354 : // Store which labels can be used to block-diagonalize the Hamiltonian
355 663 : this->blockdiagonalizing_labels.clear();
356 663 : if (sort_by_quantum_number_f) {
357 470 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_F);
358 : }
359 663 : if (sort_by_quantum_number_m) {
360 617 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_M);
361 : }
362 663 : if (sort_by_parity) {
363 518 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_PARITY);
364 : }
365 663 : }
366 :
367 : // Explicit instantiations
368 : template class SystemAtom<double>;
369 : template class SystemAtom<std::complex<double>>;
370 : } // namespace pairinteraction
|
