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/BasisPair.hpp"
8 : #include "pairinteraction/database/Database.hpp"
9 : #include "pairinteraction/enums/OperatorType.hpp"
10 : #include "pairinteraction/enums/Parity.hpp"
11 : #include "pairinteraction/enums/TransformationType.hpp"
12 : #include "pairinteraction/ket/KetAtom.hpp"
13 : #include "pairinteraction/ket/KetPair.hpp"
14 : #include "pairinteraction/system/GreenTensorInterpolator.hpp"
15 : #include "pairinteraction/system/SystemAtom.hpp"
16 : #include "pairinteraction/utils/Range.hpp"
17 : #include "pairinteraction/utils/eigen_assertion.hpp"
18 : #include "pairinteraction/utils/eigen_compat.hpp"
19 : #include "pairinteraction/utils/operator.hpp"
20 : #include "pairinteraction/utils/spherical.hpp"
21 : #include "pairinteraction/utils/streamed.hpp"
22 : #include "pairinteraction/utils/tensor.hpp"
23 : #include "pairinteraction/utils/traits.hpp"
24 :
25 : #include <Eigen/SparseCore>
26 : #include <array>
27 : #include <cmath>
28 : #include <complex>
29 : #include <limits>
30 : #include <memory>
31 : #include <spdlog/spdlog.h>
32 : #include <vector>
33 :
34 : namespace pairinteraction {
35 : template <typename Scalar>
36 : struct OperatorMatrices {
37 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> d1;
38 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> d2;
39 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> q1;
40 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> q2;
41 : };
42 :
43 : template <typename Scalar>
44 539 : GreenTensorInterpolator<Scalar> construct_green_tensor_interpolator(
45 : const std::array<typename traits::NumTraits<Scalar>::real_t, 3> &distance_vector,
46 : int interaction_order) {
47 : // https://doi.org/10.1103/PhysRevA.96.062509
48 : // https://doi.org/10.1103/PhysRevA.82.010901
49 : // https://en.wikipedia.org/wiki/Table_of_spherical_harmonics
50 :
51 : using real_t = typename traits::NumTraits<Scalar>::real_t;
52 :
53 539 : GreenTensorInterpolator<Scalar> green_tensor_interpolator;
54 :
55 : // Normalize the distance vector, return zero green tensor if the distance is infinity
56 1078 : Eigen::Map<const Eigen::Vector3<real_t>> vector_map(distance_vector.data(),
57 539 : distance_vector.size());
58 539 : real_t distance = vector_map.norm();
59 539 : SPDLOG_DEBUG("Interatomic distance: {}", distance);
60 539 : if (!std::isfinite(distance)) {
61 20 : return green_tensor_interpolator;
62 : }
63 519 : Eigen::Vector3<real_t> unitvec = vector_map / distance;
64 :
65 : // Dyadic green function of dipole-dipole interaction
66 519 : if (interaction_order >= 3) {
67 519 : Eigen::Matrix3<Scalar> entries =
68 519 : Eigen::Matrix3<real_t>::Identity() - 3 * unitvec * unitvec.transpose();
69 :
70 519 : green_tensor_interpolator.create_entries_from_cartesian(
71 1038 : 1, 1, (entries / std::pow(distance, 3)).template cast<Scalar>());
72 : }
73 :
74 : // Dyadic green function of dipole-quadrupole interaction
75 519 : if (interaction_order >= 4) {
76 20 : Eigen::Matrix<real_t, 3, 9> entries = Eigen::Matrix<real_t, 3, 9>::Zero();
77 80 : for (Eigen::Index q = 0; q < 3; ++q) {
78 60 : Eigen::Index row = q;
79 240 : for (Eigen::Index j = 0; j < 3; ++j) {
80 720 : for (Eigen::Index i = 0; i < 3; ++i) {
81 540 : Eigen::Index col = 3 * j + i;
82 540 : real_t v = 15 * unitvec[q] * unitvec[j] * unitvec[i];
83 540 : if (i == j) v += -3 * unitvec[q];
84 540 : if (i == q) v += -3 * unitvec[j];
85 540 : if (j == q) v += -3 * unitvec[i];
86 540 : entries(row, col) += v;
87 : }
88 : }
89 : }
90 :
91 20 : green_tensor_interpolator.create_entries_from_cartesian(
92 40 : 1, 2, (entries / std::pow(distance, 4)).template cast<Scalar>());
93 : }
94 :
95 : // Dyadic green function of quadrupole-dipole interaction
96 519 : if (interaction_order >= 4) {
97 20 : Eigen::Matrix<real_t, 9, 3> entries = Eigen::Matrix<real_t, 9, 3>::Zero();
98 80 : for (Eigen::Index q = 0; q < 3; ++q) {
99 240 : for (Eigen::Index j = 0; j < 3; ++j) {
100 180 : Eigen::Index row = 3 * q + j;
101 720 : for (Eigen::Index i = 0; i < 3; ++i) {
102 540 : Eigen::Index col = i;
103 540 : real_t v = -15 * unitvec[q] * unitvec[j] * unitvec[i];
104 540 : if (i == j) v += 3 * unitvec[q];
105 540 : if (i == q) v += 3 * unitvec[j];
106 540 : if (j == q) v += 3 * unitvec[i];
107 540 : entries(row, col) += v;
108 : }
109 : }
110 : }
111 :
112 20 : green_tensor_interpolator.create_entries_from_cartesian(
113 40 : 2, 1, (entries / std::pow(distance, 4)).template cast<Scalar>());
114 : }
115 :
116 : // Dyadic green function of quadrupole-quadrupole interaction
117 519 : if (interaction_order >= 5) {
118 0 : SPDLOG_WARN("Quadrupole-quadrupole interaction is considered but "
119 : "not dipole-octupole interaction although this interaction would be "
120 : "of the same order. We plan to implement dipole-octupole interaction "
121 : "in the future.");
122 :
123 0 : Eigen::Matrix<real_t, 9, 9> entries = Eigen::Matrix<real_t, 9, 9>::Zero();
124 0 : for (Eigen::Index q = 0; q < 3; ++q) {
125 0 : for (Eigen::Index j = 0; j < 3; ++j) {
126 0 : Eigen::Index row = 3 * q + j;
127 0 : for (Eigen::Index i = 0; i < 3; ++i) {
128 0 : for (Eigen::Index k = 0; k < 3; ++k) {
129 0 : Eigen::Index col = 3 * i + k;
130 0 : real_t v = 105 * unitvec[q] * unitvec[j] * unitvec[i] * unitvec[k];
131 0 : if (i == j) v += -15 * unitvec[q] * unitvec[k];
132 0 : if (i == q) v += -15 * unitvec[j] * unitvec[k];
133 0 : if (j == q) v += -15 * unitvec[i] * unitvec[k];
134 0 : if (k == q) v += -15 * unitvec[j] * unitvec[i];
135 0 : if (k == j) v += -15 * unitvec[q] * unitvec[i];
136 0 : if (k == i) v += -15 * unitvec[q] * unitvec[j];
137 0 : if (q == k && i == j) v += 3;
138 0 : if (i == k && j == q) v += 3;
139 0 : if (j == k && i == q) v += 3;
140 0 : entries(row, col) += v;
141 : }
142 : }
143 : }
144 : }
145 :
146 0 : green_tensor_interpolator.create_entries_from_cartesian(
147 0 : 2, 2, (entries / std::pow(distance, 5)).template cast<Scalar>());
148 : }
149 :
150 519 : return green_tensor_interpolator;
151 0 : }
152 :
153 : template <typename Scalar>
154 : OperatorMatrices<Scalar>
155 561 : construct_operator_matrices(const GreenTensorInterpolator<Scalar> &green_tensor_interpolator,
156 : const std::shared_ptr<const BasisAtom<Scalar>> &basis1,
157 : const std::shared_ptr<const BasisAtom<Scalar>> &basis2) {
158 : // Helper function for constructing matrices of spherical harmonics operators
159 1162 : auto get_matrices = [](auto basis, OperatorType type, std::initializer_list<int> m,
160 : bool conjugate) {
161 1162 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> matrices;
162 1162 : matrices.reserve(m.size());
163 1162 : int factor = conjugate ? -1 : 1;
164 22078 : std::transform(m.begin(), m.end(), std::back_inserter(matrices), [&](int q) {
165 3486 : auto matrix_elements = (std::pow(factor, q) *
166 : basis->get_database().get_matrix_elements_in_canonical_basis(
167 : basis, basis, type, factor * q))
168 : .eval();
169 5166 : return (basis->get_coefficients().adjoint() * matrix_elements *
170 3486 : basis->get_coefficients())
171 8652 : .eval();
172 3486 : });
173 2324 : return matrices;
174 0 : };
175 :
176 561 : OperatorMatrices<Scalar> op;
177 :
178 : // Operator matrices for Rydberg-Rydberg interaction
179 581 : if (!green_tensor_interpolator.get_spherical_entries(1, 1).empty() ||
180 20 : !green_tensor_interpolator.get_spherical_entries(1, 2).empty()) {
181 541 : op.d1 = get_matrices(basis1, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, true);
182 : }
183 581 : if (!green_tensor_interpolator.get_spherical_entries(1, 1).empty() ||
184 20 : !green_tensor_interpolator.get_spherical_entries(2, 1).empty()) {
185 541 : op.d2 = get_matrices(basis2, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, false);
186 : }
187 1122 : if (!green_tensor_interpolator.get_spherical_entries(2, 2).empty() ||
188 561 : !green_tensor_interpolator.get_spherical_entries(2, 1).empty()) {
189 20 : op.q1 = get_matrices(basis1, OperatorType::ELECTRIC_QUADRUPOLE, {-2, -1, 0, +1, +2}, true);
190 20 : op.q1.push_back(get_matrices(basis1, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, {0}, true)[0]);
191 : }
192 1122 : if (!green_tensor_interpolator.get_spherical_entries(2, 2).empty() ||
193 561 : !green_tensor_interpolator.get_spherical_entries(1, 2).empty()) {
194 20 : op.q2 = get_matrices(basis2, OperatorType::ELECTRIC_QUADRUPOLE, {-2, -1, 0, +1, +2}, false);
195 20 : op.q2.push_back(
196 40 : get_matrices(basis2, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, {0}, false)[0]);
197 : }
198 :
199 1122 : return op;
200 0 : }
201 :
202 : // "overloaded" pattern for std::visit
203 : template <class... Ts>
204 : struct overloaded : Ts... {
205 : using Ts::operator()...;
206 : };
207 : template <class... Ts>
208 : overloaded(Ts...) -> overloaded<Ts...>;
209 :
210 : template <typename Scalar>
211 561 : SystemPair<Scalar>::SystemPair(std::shared_ptr<const basis_t> basis)
212 561 : : System<SystemPair<Scalar>>(std::move(basis)) {}
213 :
214 : template <typename Scalar>
215 95 : SystemPair<Scalar> &SystemPair<Scalar>::set_interaction_order(int value) {
216 95 : this->hamiltonian_requires_construction = true;
217 :
218 95 : if (value < 3 || value > 5) {
219 0 : throw std::invalid_argument("The order must be 3, 4, or 5.");
220 : }
221 :
222 95 : if (green_tensor_interpolator) {
223 0 : throw std::invalid_argument(
224 : "Cannot set interaction order if a user-defined green tensor interpolator is set.");
225 : }
226 :
227 95 : interaction_order = value;
228 :
229 95 : return *this;
230 : }
231 :
232 : template <typename Scalar>
233 519 : SystemPair<Scalar> &SystemPair<Scalar>::set_distance_vector(const std::array<real_t, 3> &vector) {
234 519 : this->hamiltonian_requires_construction = true;
235 :
236 270 : if (!traits::NumTraits<Scalar>::is_complex_v && distance_vector[1] != 0) {
237 0 : throw std::invalid_argument(
238 : "The distance vector must not have a y-component if the scalar type is real.");
239 : }
240 :
241 519 : if (green_tensor_interpolator) {
242 0 : throw std::invalid_argument(
243 : "Cannot set distance vector if a user-defined green tensor interpolator is set.");
244 : }
245 :
246 519 : distance_vector = vector;
247 :
248 519 : return *this;
249 : }
250 :
251 : template <typename Scalar>
252 22 : SystemPair<Scalar> &SystemPair<Scalar>::set_green_tensor_interpolator(
253 : std::shared_ptr<const GreenTensorInterpolator<Scalar>> &green_tensor_interpolator) {
254 22 : this->hamiltonian_requires_construction = true;
255 :
256 44 : if (std::isfinite(distance_vector[0]) && std::isfinite(distance_vector[1]) &&
257 22 : std::isfinite(distance_vector[2])) {
258 0 : throw std::invalid_argument(
259 : "Cannot set green tensor interpolator if a finite distance vector is set.");
260 : }
261 :
262 22 : this->green_tensor_interpolator = green_tensor_interpolator;
263 :
264 22 : return *this;
265 : }
266 :
267 : template <typename Scalar>
268 561 : void SystemPair<Scalar>::construct_hamiltonian() const {
269 561 : auto basis1 = this->basis->get_basis1();
270 561 : auto basis2 = this->basis->get_basis2();
271 :
272 561 : std::shared_ptr<const GreenTensorInterpolator<Scalar>> green_tensor_interpolator_ptr;
273 561 : if (green_tensor_interpolator) {
274 22 : green_tensor_interpolator_ptr = green_tensor_interpolator;
275 : } else {
276 539 : green_tensor_interpolator_ptr = std::make_shared<const GreenTensorInterpolator<Scalar>>(
277 539 : construct_green_tensor_interpolator<Scalar>(distance_vector, interaction_order));
278 : }
279 :
280 561 : auto op = construct_operator_matrices(*green_tensor_interpolator_ptr, basis1, basis2);
281 :
282 : // Construct the unperturbed Hamiltonian in the canonical pair basis
283 561 : this->matrix = utils::get_energies_in_canonical_basis(this->basis);
284 :
285 561 : this->hamiltonian_is_diagonal = false;
286 561 : bool sort_by_quantum_number_f = this->basis->has_quantum_number_f();
287 561 : bool sort_by_quantum_number_m = this->basis->has_quantum_number_m();
288 561 : bool sort_by_parity = this->basis->has_parity();
289 :
290 : // Store the energies (they are needed in case of Rydberg-Rydberg interaction with an
291 : // OmegaDependentEntry)
292 561 : auto energies = this->matrix.diagonal().real();
293 :
294 : // Helper function for adding Rydberg-Rydberg interaction
295 2805 : auto add_interaction = [this, &energies, &sort_by_quantum_number_f, &sort_by_quantum_number_m](
296 : const auto &entries, const auto &op1, const auto &op2, int delta) {
297 4085 : for (const auto &entry : entries) {
298 1841 : std::visit(
299 : overloaded{
300 1841 : [this, &sort_by_quantum_number_m, &op1, &op2,
301 : &delta](const typename GreenTensorInterpolator<Scalar>::ConstantEntry &ce) {
302 1841 : this->matrix += ce.val() *
303 : utils::calculate_tensor_product_in_canonical_basis(
304 1841 : this->basis, this->basis, op1[ce.row()], op2[ce.col()]);
305 1841 : if (ce.row() != ce.col() + delta) {
306 98 : sort_by_quantum_number_m = false;
307 : }
308 : },
309 :
310 0 : [this, &energies, &sort_by_quantum_number_m, &op1, &op2, &delta](
311 : const typename GreenTensorInterpolator<Scalar>::OmegaDependentEntry &oe) {
312 0 : auto tensor_product = utils::calculate_tensor_product_in_canonical_basis(
313 0 : this->basis, this->basis, op1[oe.row()], op2[oe.col()]);
314 0 : for (int k = 0; k < tensor_product.outerSize(); ++k) {
315 0 : for (typename Eigen::SparseMatrix<
316 0 : Scalar, Eigen::RowMajor>::InnerIterator it(tensor_product, k);
317 0 : it; ++it) {
318 0 : it.valueRef() *= oe.val(energies(it.row()) - energies(it.col()));
319 : }
320 : }
321 0 : this->matrix += tensor_product;
322 0 : if (oe.row() != oe.col() + delta) {
323 0 : sort_by_quantum_number_m = false;
324 : }
325 0 : }},
326 : entry);
327 :
328 1841 : sort_by_quantum_number_f = false;
329 : }
330 : };
331 :
332 : // Dipole-dipole interaction
333 561 : add_interaction(green_tensor_interpolator_ptr->get_spherical_entries(1, 1), op.d1, op.d2, 0);
334 :
335 : // Dipole-quadrupole interaction
336 561 : add_interaction(green_tensor_interpolator_ptr->get_spherical_entries(1, 2), op.d1, op.q2, -1);
337 :
338 : // Quadrupole-dipole interaction
339 561 : add_interaction(green_tensor_interpolator_ptr->get_spherical_entries(2, 1), op.q1, op.d2, +1);
340 :
341 : // Quadrupole-quadrupole interaction
342 561 : add_interaction(green_tensor_interpolator_ptr->get_spherical_entries(2, 2), op.q1, op.q2, 0);
343 :
344 : // Transform from the canonical basis into the actual basis
345 561 : this->matrix =
346 561 : this->basis->get_coefficients().adjoint() * this->matrix * this->basis->get_coefficients();
347 :
348 : // Store which labels can be used to block-diagonalize the Hamiltonian
349 561 : this->blockdiagonalizing_labels.clear();
350 561 : if (sort_by_quantum_number_f) {
351 0 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_F);
352 : }
353 561 : if (sort_by_quantum_number_m) {
354 533 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_M);
355 : }
356 561 : if (sort_by_parity) {
357 0 : this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_PARITY);
358 : }
359 561 : }
360 :
361 : // Explicit instantiations
362 : template class SystemPair<double>;
363 : template class SystemPair<std::complex<double>>;
364 : } // namespace pairinteraction
|
