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