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/SorterType.hpp"
11 : #include "pairinteraction/ket/KetAtom.hpp"
12 : #include "pairinteraction/ket/KetPair.hpp"
13 : #include "pairinteraction/system/GreenTensorInterpolator.hpp"
14 : #include "pairinteraction/system/SystemAtom.hpp"
15 : #include "pairinteraction/utils/Range.hpp"
16 : #include "pairinteraction/utils/eigen_assertion.hpp"
17 : #include "pairinteraction/utils/eigen_compat.hpp"
18 : #include "pairinteraction/utils/operator.hpp"
19 : #include "pairinteraction/utils/spherical.hpp"
20 : #include "pairinteraction/utils/streamed.hpp"
21 : #include "pairinteraction/utils/tensor.hpp"
22 : #include "pairinteraction/utils/traits.hpp"
23 :
24 : #include <Eigen/SparseCore>
25 : #include <array>
26 : #include <cmath>
27 : #include <complex>
28 : #include <limits>
29 : #include <memory>
30 : #include <spdlog/spdlog.h>
31 : #include <vector>
32 :
33 : namespace pairinteraction {
34 : template <typename Scalar>
35 : struct OperatorMatrices {
36 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> m1;
37 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> m2;
38 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> d1;
39 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> d2;
40 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> q1;
41 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> q2;
42 : };
43 :
44 : template <typename Scalar>
45 609 : GreenTensorInterpolator<Scalar> construct_green_tensor_interpolator(
46 : const std::array<typename traits::NumTraits<Scalar>::real_t, 3> &distance_vector,
47 : int interaction_order, int minimal_kappa1, int minimal_kappa2) {
48 : // https://doi.org/10.1103/PhysRevA.96.062509
49 : // https://doi.org/10.1103/PhysRevA.82.010901
50 : // https://en.wikipedia.org/wiki/Table_of_spherical_harmonics
51 :
52 : using real_t = typename traits::NumTraits<Scalar>::real_t;
53 :
54 609 : GreenTensorInterpolator<Scalar> green_tensor_interpolator;
55 :
56 : // Normalize the distance vector, return zero green tensor if the distance is infinity
57 1191 : Eigen::Map<const Eigen::Vector3<real_t>> vector_map(distance_vector.data(),
58 601 : distance_vector.size());
59 584 : real_t distance = vector_map.norm();
60 592 : SPDLOG_DEBUG("Interatomic distance: {}", distance);
61 622 : if (!std::isfinite(distance)) {
62 26 : return green_tensor_interpolator;
63 : }
64 595 : Eigen::Vector3<real_t> unitvec = vector_map / distance;
65 :
66 10498 : auto is_required = [interaction_order, minimal_kappa1, minimal_kappa2](int kappa1, int kappa2) {
67 4680 : return kappa1 + kappa2 + 1 <= interaction_order && kappa1 >= minimal_kappa1 &&
68 4680 : kappa2 >= minimal_kappa2;
69 : };
70 :
71 : // The cartesian tensors follow from the Taylor expansion of the Coulomb interaction
72 : // 1/|R + r2 - r1| in the coordinates r1 and r2 of the two charge distributions. The tensor
73 : // for the orders (kappa1, kappa2) is (-1)^kappa1 * d^(kappa1+kappa2)/dR^(kappa1+kappa2) (1/R).
74 : // Note that the factor 1/(kappa1! * kappa2!) of the Taylor expansion is not included because
75 : // it is absorbed in the normalization of the cartesian-to-spherical transformation, see
76 : // spherical.hpp.
77 :
78 : // Green function of monopole-monopole interaction
79 595 : if (is_required(0, 0)) {
80 56 : Eigen::Matrix<real_t, 1, 1> entries = Eigen::Matrix<real_t, 1, 1>::Constant(1);
81 :
82 56 : green_tensor_interpolator.create_entries_from_cartesian(
83 58 : 0, 0, (entries / distance).template cast<Scalar>());
84 : }
85 :
86 : // Green function of monopole-dipole interaction
87 580 : if (is_required(0, 1)) {
88 38 : Eigen::Matrix<real_t, 1, 3> entries = -unitvec.transpose();
89 :
90 38 : green_tensor_interpolator.create_entries_from_cartesian(
91 76 : 0, 1, (entries / std::pow(distance, 2)).template cast<Scalar>());
92 : }
93 :
94 : // Green function of dipole-monopole interaction
95 591 : if (is_required(1, 0)) {
96 50 : const Eigen::Matrix<real_t, 3, 1> &entries = unitvec;
97 :
98 50 : green_tensor_interpolator.create_entries_from_cartesian(
99 100 : 1, 0, (entries / std::pow(distance, 2)).template cast<Scalar>());
100 : }
101 :
102 : // Green function of monopole-quadrupole interaction
103 586 : if (is_required(0, 2)) {
104 20 : Eigen::Matrix<real_t, 1, 9> entries = Eigen::Matrix<real_t, 1, 9>::Zero();
105 80 : for (Eigen::Index j = 0; j < 3; ++j) {
106 240 : for (Eigen::Index i = 0; i < 3; ++i) {
107 180 : Eigen::Index col = 3 * j + i;
108 180 : real_t v = 3 * unitvec[j] * unitvec[i];
109 180 : if (i == j) v += -1;
110 180 : entries(0, col) += v;
111 : }
112 : }
113 :
114 20 : green_tensor_interpolator.create_entries_from_cartesian(
115 40 : 0, 2, (entries / std::pow(distance, 3)).template cast<Scalar>());
116 : }
117 :
118 : // Green function of quadrupole-monopole interaction
119 587 : if (is_required(2, 0)) {
120 26 : Eigen::Matrix<real_t, 9, 1> entries = Eigen::Matrix<real_t, 9, 1>::Zero();
121 104 : for (Eigen::Index q = 0; q < 3; ++q) {
122 312 : for (Eigen::Index j = 0; j < 3; ++j) {
123 234 : Eigen::Index row = 3 * q + j;
124 234 : real_t v = 3 * unitvec[q] * unitvec[j];
125 234 : if (q == j) v += -1;
126 234 : entries(row, 0) += v;
127 : }
128 : }
129 :
130 26 : green_tensor_interpolator.create_entries_from_cartesian(
131 52 : 2, 0, (entries / std::pow(distance, 3)).template cast<Scalar>());
132 : }
133 :
134 : // Dyadic green function of dipole-dipole interaction
135 591 : if (is_required(1, 1)) {
136 544 : Eigen::Matrix3<Scalar> entries =
137 545 : Eigen::Matrix3<real_t>::Identity() - 3 * unitvec * unitvec.transpose();
138 :
139 528 : green_tensor_interpolator.create_entries_from_cartesian(
140 1061 : 1, 1, (entries / std::pow(distance, 3)).template cast<Scalar>());
141 : }
142 :
143 : // Dyadic green function of dipole-quadrupole interaction
144 562 : if (is_required(1, 2)) {
145 20 : Eigen::Matrix<real_t, 3, 9> entries = Eigen::Matrix<real_t, 3, 9>::Zero();
146 80 : for (Eigen::Index q = 0; q < 3; ++q) {
147 60 : Eigen::Index row = q;
148 240 : for (Eigen::Index j = 0; j < 3; ++j) {
149 720 : for (Eigen::Index i = 0; i < 3; ++i) {
150 540 : Eigen::Index col = 3 * j + i;
151 540 : real_t v = 15 * unitvec[q] * unitvec[j] * unitvec[i];
152 540 : if (i == j) v += -3 * unitvec[q];
153 540 : if (i == q) v += -3 * unitvec[j];
154 540 : if (j == q) v += -3 * unitvec[i];
155 540 : entries(row, col) += v;
156 : }
157 : }
158 : }
159 :
160 20 : green_tensor_interpolator.create_entries_from_cartesian(
161 40 : 1, 2, (entries / std::pow(distance, 4)).template cast<Scalar>());
162 : }
163 :
164 : // Dyadic green function of quadrupole-dipole interaction
165 574 : if (is_required(2, 1)) {
166 18 : Eigen::Matrix<real_t, 9, 3> entries = Eigen::Matrix<real_t, 9, 3>::Zero();
167 80 : for (Eigen::Index q = 0; q < 3; ++q) {
168 240 : for (Eigen::Index j = 0; j < 3; ++j) {
169 180 : Eigen::Index row = 3 * q + j;
170 720 : for (Eigen::Index i = 0; i < 3; ++i) {
171 540 : Eigen::Index col = i;
172 540 : real_t v = -15 * unitvec[q] * unitvec[j] * unitvec[i];
173 540 : if (i == j) v += 3 * unitvec[q];
174 540 : if (i == q) v += 3 * unitvec[j];
175 540 : if (j == q) v += 3 * unitvec[i];
176 540 : entries(row, col) += v;
177 : }
178 : }
179 : }
180 :
181 20 : green_tensor_interpolator.create_entries_from_cartesian(
182 40 : 2, 1, (entries / std::pow(distance, 4)).template cast<Scalar>());
183 : }
184 :
185 : // Dyadic green function of quadrupole-quadrupole interaction
186 585 : if (interaction_order >= 5) {
187 0 : SPDLOG_WARN("Quadrupole-quadrupole interaction is considered but "
188 : "not dipole-octupole interaction although this interaction would be "
189 : "of the same order. We plan to implement dipole-octupole interaction "
190 : "in the future.");
191 :
192 0 : Eigen::Matrix<real_t, 9, 9> entries = Eigen::Matrix<real_t, 9, 9>::Zero();
193 0 : for (Eigen::Index q = 0; q < 3; ++q) {
194 0 : for (Eigen::Index j = 0; j < 3; ++j) {
195 0 : Eigen::Index row = 3 * q + j;
196 0 : for (Eigen::Index i = 0; i < 3; ++i) {
197 0 : for (Eigen::Index k = 0; k < 3; ++k) {
198 0 : Eigen::Index col = 3 * i + k;
199 0 : real_t v = 105 * unitvec[q] * unitvec[j] * unitvec[i] * unitvec[k];
200 0 : if (i == j) v += -15 * unitvec[q] * unitvec[k];
201 0 : if (i == q) v += -15 * unitvec[j] * unitvec[k];
202 0 : if (j == q) v += -15 * unitvec[i] * unitvec[k];
203 0 : if (k == q) v += -15 * unitvec[j] * unitvec[i];
204 0 : if (k == j) v += -15 * unitvec[q] * unitvec[i];
205 0 : if (k == i) v += -15 * unitvec[q] * unitvec[j];
206 0 : if (q == k && i == j) v += 3;
207 0 : if (i == k && j == q) v += 3;
208 0 : if (j == k && i == q) v += 3;
209 0 : entries(row, col) += v;
210 : }
211 : }
212 : }
213 : }
214 :
215 0 : green_tensor_interpolator.create_entries_from_cartesian(
216 0 : 2, 2, (entries / std::pow(distance, 5)).template cast<Scalar>());
217 : }
218 :
219 585 : return green_tensor_interpolator;
220 0 : }
221 :
222 : template <typename Scalar>
223 : OperatorMatrices<Scalar>
224 619 : construct_operator_matrices(const GreenTensorInterpolator<Scalar> &green_tensor_interpolator,
225 : const std::shared_ptr<const BasisAtom<Scalar>> &basis1,
226 : const std::shared_ptr<const BasisAtom<Scalar>> &basis2) {
227 : // Helper function for constructing matrices of spherical harmonics operators
228 1461 : auto get_matrices = [](auto basis, OperatorType type, std::initializer_list<int> m,
229 : bool conjugate) {
230 1461 : std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> matrices;
231 1461 : matrices.reserve(m.size());
232 1436 : int factor = conjugate ? -1 : 1;
233 26313 : std::transform(m.begin(), m.end(), std::back_inserter(matrices), [&](int q) {
234 4157 : auto matrix_elements = (std::pow(factor, q) *
235 : basis->get_database().get_matrix_elements_in_canonical_basis(
236 : basis, basis, type, factor * q))
237 : .eval();
238 6134 : return (basis->get_coefficients().adjoint() * matrix_elements *
239 4163 : basis->get_coefficients())
240 10275 : .eval();
241 4188 : });
242 2958 : return matrices;
243 0 : };
244 :
245 619 : OperatorMatrices<Scalar> op;
246 :
247 10198 : auto has = [&green_tensor_interpolator](int kappa1, int kappa2) {
248 9579 : return !green_tensor_interpolator.get_spherical_entries(kappa1, kappa2).empty();
249 : };
250 :
251 : // Operator matrices for Rydberg-Rydberg interaction
252 619 : if (has(0, 0) || has(0, 1) || has(0, 2)) {
253 56 : op.m1 = get_matrices(basis1, OperatorType::ELECTRIC_MONOPOLE, {0}, true);
254 : }
255 600 : if (has(0, 0) || has(1, 0) || has(2, 0)) {
256 68 : op.m2 = get_matrices(basis2, OperatorType::ELECTRIC_MONOPOLE, {0}, false);
257 : }
258 609 : if (has(1, 0) || has(1, 1) || has(1, 2)) {
259 567 : op.d1 = get_matrices(basis1, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, true);
260 : }
261 625 : if (has(0, 1) || has(1, 1) || has(2, 1)) {
262 575 : op.d2 = get_matrices(basis2, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, false);
263 : }
264 623 : if (has(2, 0) || has(2, 1) || has(2, 2)) {
265 46 : op.q1 = get_matrices(basis1, OperatorType::ELECTRIC_QUADRUPOLE, {-2, -1, 0, +1, +2}, true);
266 46 : op.q1.push_back(get_matrices(basis1, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, {0}, true)[0]);
267 : }
268 635 : if (has(0, 2) || has(1, 2) || has(2, 2)) {
269 40 : op.q2 = get_matrices(basis2, OperatorType::ELECTRIC_QUADRUPOLE, {-2, -1, 0, +1, +2}, false);
270 40 : op.q2.push_back(
271 80 : get_matrices(basis2, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, {0}, false)[0]);
272 : }
273 :
274 1266 : return op;
275 0 : }
276 :
277 : template <typename Scalar>
278 643 : SystemPair<Scalar>::SystemPair(std::shared_ptr<const basis_t> basis)
279 643 : : System<SystemPair<Scalar>>(std::move(basis)) {}
280 :
281 : template <typename Scalar>
282 166 : SystemPair<Scalar> &SystemPair<Scalar>::set_interaction_order(int value) {
283 166 : this->hamiltonian_requires_construction = true;
284 :
285 166 : if (value > 5) {
286 0 : throw std::invalid_argument("Interaction orders larger than 5 are not yet supported.");
287 : }
288 :
289 166 : if (green_tensor_interpolator) {
290 0 : throw std::invalid_argument(
291 : "Cannot set interaction order if a user-defined green tensor interpolator is set.");
292 : }
293 :
294 166 : interaction_order = value;
295 :
296 166 : return *this;
297 : }
298 :
299 : template <typename Scalar>
300 597 : SystemPair<Scalar> &SystemPair<Scalar>::set_distance_vector(const std::array<real_t, 3> &vector) {
301 597 : this->hamiltonian_requires_construction = true;
302 :
303 311 : if (!traits::NumTraits<Scalar>::is_complex_v && vector[1] != 0) {
304 1 : throw std::invalid_argument(
305 : "The distance vector must not have a y-component if the scalar type is real.");
306 : }
307 :
308 596 : if (green_tensor_interpolator) {
309 0 : throw std::invalid_argument(
310 : "Cannot set distance vector if a user-defined green tensor interpolator is set.");
311 : }
312 :
313 596 : distance_vector = vector;
314 :
315 596 : return *this;
316 : }
317 :
318 : template <typename Scalar>
319 20 : SystemPair<Scalar> &SystemPair<Scalar>::set_green_tensor_interpolator(
320 : const std::shared_ptr<const GreenTensorInterpolator<Scalar>> &green_tensor_interpolator) {
321 20 : this->hamiltonian_requires_construction = true;
322 :
323 40 : if (std::isfinite(distance_vector[0]) && std::isfinite(distance_vector[1]) &&
324 20 : std::isfinite(distance_vector[2])) {
325 0 : throw std::invalid_argument(
326 : "Cannot set green tensor interpolator if a finite distance vector is set.");
327 : }
328 :
329 20 : this->green_tensor_interpolator = green_tensor_interpolator;
330 :
331 20 : return *this;
332 : }
333 :
334 : template <typename Scalar>
335 641 : void SystemPair<Scalar>::construct_hamiltonian() const {
336 641 : auto basis1 = this->basis->get_basis1();
337 642 : auto basis2 = this->basis->get_basis2();
338 :
339 641 : int minimal_kappa1 = basis1->get_database().get_minimal_kappa(basis1->get_species());
340 636 : int minimal_kappa2 = basis2->get_database().get_minimal_kappa(basis2->get_species());
341 :
342 640 : std::shared_ptr<const GreenTensorInterpolator<Scalar>> green_tensor_interpolator_ptr;
343 640 : if (green_tensor_interpolator) {
344 20 : green_tensor_interpolator_ptr = green_tensor_interpolator;
345 : } else {
346 601 : green_tensor_interpolator_ptr = std::make_shared<const GreenTensorInterpolator<Scalar>>(
347 612 : construct_green_tensor_interpolator<Scalar>(distance_vector, interaction_order,
348 : minimal_kappa1, minimal_kappa2));
349 : }
350 :
351 595 : auto op = construct_operator_matrices(*green_tensor_interpolator_ptr, basis1, basis2);
352 :
353 : // Construct the unperturbed Hamiltonian in the canonical pair basis
354 636 : this->matrix = utils::get_energies_in_canonical_basis(this->basis);
355 :
356 641 : this->hamiltonian_is_diagonal = false;
357 641 : bool sort_by_quantum_number_f = this->basis->has_quantum_number("f");
358 634 : bool sort_by_quantum_number_m = this->basis->has_quantum_number("m");
359 629 : bool sort_by_parity = this->basis->has_quantum_number("product_of_parities");
360 :
361 : // Add Rydberg-Rydberg interaction via Green tensor
362 : // H_RR = Σ_{ij} D_1,left[i] * G_{ij} * D_2,right[j]
363 : // where D_1,left uses conjugated convention and
364 : // D_2,right uses normal convention.
365 :
366 : // Helper function for adding Rydberg-Rydberg interaction.
367 6256 : auto add_interaction = [this, &green_tensor_interpolator_ptr, &sort_by_quantum_number_f,
368 : &sort_by_quantum_number_m, &sort_by_parity](
369 : const auto &op1, const auto &op2, int kappa1, int kappa2) {
370 5636 : const auto &entries = green_tensor_interpolator_ptr->get_spherical_entries(kappa1, kappa2);
371 :
372 7744 : for (const auto &entry : entries) {
373 2109 : if (std::holds_alternative<
374 2089 : typename GreenTensorInterpolator<Scalar>::OmegaDependentEntry>(entry)) {
375 0 : throw std::logic_error(
376 : "Green tensor with omega dependent entries is currently not supported.");
377 : }
378 :
379 : const auto &constant_entry =
380 2109 : std::get<typename GreenTensorInterpolator<Scalar>::ConstantEntry>(entry);
381 :
382 : // Skip vanishing operators (e.g., for small basis that comprises a single state)
383 : // so that they do not break the block-diagonal structure of the Hamiltonian
384 4096 : if (op1[constant_entry.row()].nonZeros() == 0 ||
385 2000 : op2[constant_entry.col()].nonZeros() == 0) {
386 178 : continue;
387 : }
388 :
389 1946 : this->matrix += constant_entry.val() *
390 1951 : utils::calculate_tensor_product_in_canonical_basis(this->basis, this->basis,
391 1952 : op1[constant_entry.row()],
392 1953 : op2[constant_entry.col()]);
393 :
394 1937 : if (kappa1 != 0 || kappa2 != 0) {
395 1870 : sort_by_quantum_number_f = false;
396 : }
397 1937 : if (constant_entry.row() != constant_entry.col() + kappa1 - kappa2) {
398 72 : sort_by_quantum_number_m = false;
399 : }
400 1937 : if ((kappa1 + kappa2) % 2 != 0) {
401 178 : sort_by_parity = false;
402 : }
403 : }
404 : };
405 :
406 : // Monopole-monopole interaction
407 630 : add_interaction(op.m1, op.m2, 0, 0);
408 :
409 : // Monopole-dipole interaction
410 617 : add_interaction(op.m1, op.d2, 0, 1);
411 :
412 : // Dipole-monopole interaction
413 620 : add_interaction(op.d1, op.m2, 1, 0);
414 :
415 : // Monopole-quadrupole interaction
416 623 : add_interaction(op.m1, op.q2, 0, 2);
417 :
418 : // Quadrupole-monopole interaction
419 619 : add_interaction(op.q1, op.m2, 2, 0);
420 :
421 : // Dipole-dipole interaction
422 617 : add_interaction(op.d1, op.d2, 1, 1);
423 :
424 : // Dipole-quadrupole interaction
425 642 : add_interaction(op.d1, op.q2, 1, 2);
426 :
427 : // Quadrupole-dipole interaction
428 626 : add_interaction(op.q1, op.d2, 2, 1);
429 :
430 : // Quadrupole-quadrupole interaction
431 638 : add_interaction(op.q1, op.q2, 2, 2);
432 :
433 : // Transform from the canonical basis into the actual basis
434 625 : this->matrix =
435 636 : this->basis->get_coefficients().adjoint() * this->matrix * this->basis->get_coefficients();
436 :
437 : // Store which labels can be used to block-diagonalize the Hamiltonian
438 642 : this->blockdiagonalizing_labels.clear();
439 633 : if (sort_by_quantum_number_f) {
440 0 : this->blockdiagonalizing_labels.push_back(SorterType::QUANTUM_NUMBER_F);
441 : }
442 633 : if (sort_by_quantum_number_m) {
443 614 : this->blockdiagonalizing_labels.push_back(SorterType::QUANTUM_NUMBER_M);
444 : }
445 621 : if (sort_by_parity) {
446 0 : this->blockdiagonalizing_labels.push_back(SorterType::PARITY);
447 : }
448 621 : }
449 :
450 : // Explicit instantiations
451 : template class SystemPair<double>;
452 : template class SystemPair<std::complex<double>>;
453 : } // namespace pairinteraction
|
