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