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         546 : 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         546 :     GreenTensorInterpolator<Scalar> green_tensor_interpolator;
      54             : 
      55             :     // Normalize the distance vector, return zero green tensor if the distance is infinity
      56        1092 :     Eigen::Map<const Eigen::Vector3<real_t>> vector_map(distance_vector.data(),
      57         546 :                                                         distance_vector.size());
      58         546 :     real_t distance = vector_map.norm();
      59         546 :     SPDLOG_DEBUG("Interatomic distance: {}", distance);
      60         546 :     if (!std::isfinite(distance)) {
      61          22 :         return green_tensor_interpolator;
      62             :     }
      63         524 :     Eigen::Vector3<real_t> unitvec = vector_map / distance;
      64             : 
      65             :     // Dyadic green function of dipole-dipole interaction
      66         524 :     if (interaction_order >= 3) {
      67         524 :         Eigen::Matrix3<Scalar> entries =
      68         524 :             Eigen::Matrix3<real_t>::Identity() - 3 * unitvec * unitvec.transpose();
      69             : 
      70         524 :         green_tensor_interpolator.create_entries_from_cartesian(
      71        1048 :             1, 1, (entries / std::pow(distance, 3)).template cast<Scalar>());
      72             :     }
      73             : 
      74             :     // Dyadic green function of dipole-quadrupole interaction
      75         524 :     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         524 :     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         524 :     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         524 :     return green_tensor_interpolator;
     151           0 : }
     152             : 
     153             : template <typename Scalar>
     154             : OperatorMatrices<Scalar>
     155         566 : 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        1168 :     auto get_matrices = [](auto basis, OperatorType type, std::initializer_list<int> m,
     160             :                            bool conjugate) {
     161        1168 :         std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> matrices;
     162        1168 :         matrices.reserve(m.size());
     163        1168 :         int factor = conjugate ? -1 : 1;
     164       22192 :         std::transform(m.begin(), m.end(), std::back_inserter(matrices), [&](int q) {
     165        3504 :             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        5184 :             return (basis->get_coefficients().adjoint() * matrix_elements *
     170        3504 :                     basis->get_coefficients())
     171        8688 :                 .eval();
     172        3504 :         });
     173        2336 :         return matrices;
     174           0 :     };
     175             : 
     176         566 :     OperatorMatrices<Scalar> op;
     177             : 
     178             :     // Operator matrices for Rydberg-Rydberg interaction
     179         588 :     if (!green_tensor_interpolator.get_spherical_entries(1, 1).empty() ||
     180          22 :         !green_tensor_interpolator.get_spherical_entries(1, 2).empty()) {
     181         544 :         op.d1 = get_matrices(basis1, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, true);
     182             :     }
     183         588 :     if (!green_tensor_interpolator.get_spherical_entries(1, 1).empty() ||
     184          22 :         !green_tensor_interpolator.get_spherical_entries(2, 1).empty()) {
     185         544 :         op.d2 = get_matrices(basis2, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, false);
     186             :     }
     187        1132 :     if (!green_tensor_interpolator.get_spherical_entries(2, 2).empty() ||
     188         566 :         !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        1132 :     if (!green_tensor_interpolator.get_spherical_entries(2, 2).empty() ||
     193         566 :         !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        1132 :     return op;
     200           0 : }
     201             : 
     202             : template <typename Scalar>
     203         567 : SystemPair<Scalar>::SystemPair(std::shared_ptr<const basis_t> basis)
     204         567 :     : System<SystemPair<Scalar>>(std::move(basis)) {}
     205             : 
     206             : template <typename Scalar>
     207          94 : SystemPair<Scalar> &SystemPair<Scalar>::set_interaction_order(int value) {
     208          94 :     this->hamiltonian_requires_construction = true;
     209             : 
     210          94 :     if (value < 3 || value > 5) {
     211           0 :         throw std::invalid_argument("The order must be 3, 4, or 5.");
     212             :     }
     213             : 
     214          94 :     if (green_tensor_interpolator) {
     215           0 :         throw std::invalid_argument(
     216             :             "Cannot set interaction order if a user-defined green tensor interpolator is set.");
     217             :     }
     218             : 
     219          94 :     interaction_order = value;
     220             : 
     221          94 :     return *this;
     222             : }
     223             : 
     224             : template <typename Scalar>
     225         525 : SystemPair<Scalar> &SystemPair<Scalar>::set_distance_vector(const std::array<real_t, 3> &vector) {
     226         525 :     this->hamiltonian_requires_construction = true;
     227             : 
     228         276 :     if (!traits::NumTraits<Scalar>::is_complex_v && vector[1] != 0) {
     229           1 :         throw std::invalid_argument(
     230             :             "The distance vector must not have a y-component if the scalar type is real.");
     231             :     }
     232             : 
     233         524 :     if (green_tensor_interpolator) {
     234           0 :         throw std::invalid_argument(
     235             :             "Cannot set distance vector if a user-defined green tensor interpolator is set.");
     236             :     }
     237             : 
     238         524 :     distance_vector = vector;
     239             : 
     240         524 :     return *this;
     241             : }
     242             : 
     243             : template <typename Scalar>
     244          20 : SystemPair<Scalar> &SystemPair<Scalar>::set_green_tensor_interpolator(
     245             :     const std::shared_ptr<const GreenTensorInterpolator<Scalar>> &green_tensor_interpolator) {
     246          20 :     this->hamiltonian_requires_construction = true;
     247             : 
     248          40 :     if (std::isfinite(distance_vector[0]) && std::isfinite(distance_vector[1]) &&
     249          20 :         std::isfinite(distance_vector[2])) {
     250           0 :         throw std::invalid_argument(
     251             :             "Cannot set green tensor interpolator if a finite distance vector is set.");
     252             :     }
     253             : 
     254          20 :     this->green_tensor_interpolator = green_tensor_interpolator;
     255             : 
     256          20 :     return *this;
     257             : }
     258             : 
     259             : template <typename Scalar>
     260         566 : void SystemPair<Scalar>::construct_hamiltonian() const {
     261         566 :     auto basis1 = this->basis->get_basis1();
     262         566 :     auto basis2 = this->basis->get_basis2();
     263             : 
     264         566 :     std::shared_ptr<const GreenTensorInterpolator<Scalar>> green_tensor_interpolator_ptr;
     265         566 :     if (green_tensor_interpolator) {
     266          20 :         green_tensor_interpolator_ptr = green_tensor_interpolator;
     267             :     } else {
     268         546 :         green_tensor_interpolator_ptr = std::make_shared<const GreenTensorInterpolator<Scalar>>(
     269         546 :             construct_green_tensor_interpolator<Scalar>(distance_vector, interaction_order));
     270             :     }
     271             : 
     272         566 :     auto op = construct_operator_matrices(*green_tensor_interpolator_ptr, basis1, basis2);
     273             : 
     274             :     // Construct the unperturbed Hamiltonian in the canonical pair basis
     275         566 :     this->matrix = utils::get_energies_in_canonical_basis(this->basis);
     276             : 
     277         566 :     this->hamiltonian_is_diagonal = false;
     278         566 :     bool sort_by_quantum_number_f = this->basis->has_quantum_number_f();
     279         566 :     bool sort_by_quantum_number_m = this->basis->has_quantum_number_m();
     280         566 :     bool sort_by_parity = this->basis->has_parity();
     281             : 
     282             :     // Add Rydberg-Rydberg interaction via Green tensor
     283             :     // H_RR = Σ_{ij} D_1,left[i] * G_{ij} * D_2,right[j]
     284             :     // where D_1,left uses conjugated convention and
     285             :     // D_2,right uses normal convention.
     286             : 
     287             :     // Helper function for adding Rydberg-Rydberg interaction
     288        2830 :     auto add_interaction = [this, &sort_by_quantum_number_f, &sort_by_quantum_number_m](
     289             :                                const auto &entries, const auto &op1, const auto &op2, int delta) {
     290        4106 :         for (const auto &entry : entries) {
     291        1842 :             if (std::holds_alternative<
     292        1842 :                     typename GreenTensorInterpolator<Scalar>::OmegaDependentEntry>(entry)) {
     293           0 :                 throw std::logic_error(
     294             :                     "Green tensor with omega dependent entries is currently not supported.");
     295             :             }
     296             : 
     297             :             const auto &constant_entry =
     298        1842 :                 std::get<typename GreenTensorInterpolator<Scalar>::ConstantEntry>(entry);
     299        1842 :             this->matrix += constant_entry.val() *
     300        1842 :                 utils::calculate_tensor_product_in_canonical_basis(this->basis, this->basis,
     301        1842 :                                                                    op1[constant_entry.row()],
     302        1842 :                                                                    op2[constant_entry.col()]);
     303             : 
     304        1842 :             sort_by_quantum_number_f = false;
     305        1842 :             if (constant_entry.row() != constant_entry.col() + delta) {
     306          90 :                 sort_by_quantum_number_m = false;
     307             :             }
     308             :         }
     309             :     };
     310             : 
     311             :     // Dipole-dipole interaction
     312         566 :     add_interaction(green_tensor_interpolator_ptr->get_spherical_entries(1, 1), op.d1, op.d2, 0);
     313             : 
     314             :     // Dipole-quadrupole interaction
     315         566 :     add_interaction(green_tensor_interpolator_ptr->get_spherical_entries(1, 2), op.d1, op.q2, -1);
     316             : 
     317             :     // Quadrupole-dipole interaction
     318         566 :     add_interaction(green_tensor_interpolator_ptr->get_spherical_entries(2, 1), op.q1, op.d2, +1);
     319             : 
     320             :     // Quadrupole-quadrupole interaction
     321         566 :     add_interaction(green_tensor_interpolator_ptr->get_spherical_entries(2, 2), op.q1, op.q2, 0);
     322             : 
     323             :     // Transform from the canonical basis into the actual basis
     324         566 :     this->matrix =
     325         566 :         this->basis->get_coefficients().adjoint() * this->matrix * this->basis->get_coefficients();
     326             : 
     327             :     // Store which labels can be used to block-diagonalize the Hamiltonian
     328         566 :     this->blockdiagonalizing_labels.clear();
     329         566 :     if (sort_by_quantum_number_f) {
     330           0 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_F);
     331             :     }
     332         566 :     if (sort_by_quantum_number_m) {
     333         542 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_M);
     334             :     }
     335         566 :     if (sort_by_parity) {
     336           0 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_PARITY);
     337             :     }
     338         566 : }
     339             : 
     340             : // Explicit instantiations
     341             : template class SystemPair<double>;
     342             : template class SystemPair<std::complex<double>>;
     343             : } // namespace pairinteraction