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/enums/OperatorType.hpp"
       9             : #include "pairinteraction/enums/Parity.hpp"
      10             : #include "pairinteraction/enums/TransformationType.hpp"
      11             : #include "pairinteraction/ket/KetAtom.hpp"
      12             : #include "pairinteraction/ket/KetPair.hpp"
      13             : #include "pairinteraction/operator/OperatorAtom.hpp"
      14             : #include "pairinteraction/operator/OperatorPair.hpp"
      15             : #include "pairinteraction/system/GreenTensor.hpp"
      16             : #include "pairinteraction/system/SystemAtom.hpp"
      17             : #include "pairinteraction/utils/Range.hpp"
      18             : #include "pairinteraction/utils/eigen_assertion.hpp"
      19             : #include "pairinteraction/utils/eigen_compat.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 <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>> 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             : 
      42             : template <typename Scalar>
      43         147 : GreenTensor<Scalar> construct_green_tensor(
      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             : 
      50             :     using real_t = typename traits::NumTraits<Scalar>::real_t;
      51             : 
      52         147 :     GreenTensor<Scalar> green_tensor;
      53             : 
      54             :     // Normalize the distance vector, return zero green tensor if the distance is infinity
      55         294 :     Eigen::Map<const Eigen::Vector3<real_t>> vector_map(distance_vector.data(),
      56         147 :                                                         distance_vector.size());
      57         147 :     real_t distance = vector_map.norm();
      58         294 :     SPDLOG_DEBUG("Interatomic distance: {}", distance);
      59         147 :     if (!std::isfinite(distance)) {
      60          78 :         return green_tensor;
      61             :     }
      62          69 :     Eigen::Vector3<real_t> unitvec = vector_map / distance;
      63             : 
      64             :     // Dyadic green function of dipole-dipole interaction
      65          69 :     if (interaction_order >= 3) {
      66          69 :         Eigen::Matrix3<Scalar> entries =
      67          69 :             Eigen::Matrix3<real_t>::Identity() - 3 * unitvec * unitvec.transpose();
      68             : 
      69          69 :         green_tensor.create_entries_from_cartesian(
      70         138 :             1, 1, (entries / std::pow(distance, 3)).template cast<Scalar>());
      71             :     }
      72             : 
      73             :     // Dyadic green function of dipole-quadrupole interaction
      74          69 :     if (interaction_order >= 4) {
      75          10 :         Eigen::Matrix<real_t, 3, 9> entries = Eigen::Matrix<real_t, 3, 9>::Zero();
      76          40 :         for (Eigen::Index q = 0; q < 3; ++q) {
      77          30 :             Eigen::Index row = q;
      78         120 :             for (Eigen::Index j = 0; j < 3; ++j) {
      79         360 :                 for (Eigen::Index i = 0; i < 3; ++i) {
      80         270 :                     Eigen::Index col = 3 * j + i;
      81         270 :                     real_t v = 15 * unitvec[q] * unitvec[j] * unitvec[i];
      82         270 :                     if (i == j) v += -3 * unitvec[q];
      83         270 :                     if (i == q) v += -3 * unitvec[j];
      84         270 :                     if (j == q) v += -3 * unitvec[i];
      85         270 :                     entries(row, col) += v;
      86             :                 }
      87             :             }
      88             :         }
      89             : 
      90          10 :         green_tensor.create_entries_from_cartesian(
      91          20 :             1, 2, (entries / std::pow(distance, 4)).template cast<Scalar>());
      92             :     }
      93             : 
      94             :     // Dyadic green function of quadrupole-dipole interaction
      95          69 :     if (interaction_order >= 4) {
      96          10 :         Eigen::Matrix<real_t, 9, 3> entries = Eigen::Matrix<real_t, 9, 3>::Zero();
      97          40 :         for (Eigen::Index q = 0; q < 3; ++q) {
      98         120 :             for (Eigen::Index j = 0; j < 3; ++j) {
      99          90 :                 Eigen::Index row = 3 * q + j;
     100         360 :                 for (Eigen::Index i = 0; i < 3; ++i) {
     101         270 :                     Eigen::Index col = i;
     102         270 :                     real_t v = -15 * unitvec[q] * unitvec[j] * unitvec[i];
     103         270 :                     if (i == j) v += 3 * unitvec[q];
     104         270 :                     if (i == q) v += 3 * unitvec[j];
     105         270 :                     if (j == q) v += 3 * unitvec[i];
     106         270 :                     entries(row, col) += v;
     107             :                 }
     108             :             }
     109             :         }
     110             : 
     111          10 :         green_tensor.create_entries_from_cartesian(
     112          20 :             2, 1, (entries / std::pow(distance, 4)).template cast<Scalar>());
     113             :     }
     114             : 
     115             :     // Dyadic green function of quadrupole-quadrupole interaction
     116          69 :     if (interaction_order >= 5) {
     117           0 :         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           0 :         Eigen::Matrix<real_t, 9, 9> entries = Eigen::Matrix<real_t, 9, 9>::Zero();
     123           0 :         for (Eigen::Index q = 0; q < 3; ++q) {
     124           0 :             for (Eigen::Index j = 0; j < 3; ++j) {
     125           0 :                 Eigen::Index row = 3 * q + j;
     126           0 :                 for (Eigen::Index i = 0; i < 3; ++i) {
     127           0 :                     for (Eigen::Index k = 0; k < 3; ++k) {
     128           0 :                         Eigen::Index col = 3 * i + k;
     129           0 :                         real_t v = 105 * unitvec[q] * unitvec[j] * unitvec[i] * unitvec[k];
     130           0 :                         if (i == j) v += -15 * unitvec[q] * unitvec[k];
     131           0 :                         if (i == q) v += -15 * unitvec[j] * unitvec[k];
     132           0 :                         if (j == q) v += -15 * unitvec[i] * unitvec[k];
     133           0 :                         if (k == q) v += -15 * unitvec[j] * unitvec[i];
     134           0 :                         if (k == j) v += -15 * unitvec[q] * unitvec[i];
     135           0 :                         if (k == i) v += -15 * unitvec[q] * unitvec[j];
     136           0 :                         if (q == k && i == j) v += 3;
     137           0 :                         if (i == k && j == q) v += 3;
     138           0 :                         if (j == k && i == q) v += 3;
     139           0 :                         entries(row, col) += v;
     140             :                     }
     141             :                 }
     142             :             }
     143             :         }
     144             : 
     145           0 :         green_tensor.create_entries_from_cartesian(
     146           0 :             2, 2, (entries / std::pow(distance, 5)).template cast<Scalar>());
     147             :     }
     148             : 
     149          69 :     return green_tensor;
     150           0 : }
     151             : 
     152             : template <typename Scalar>
     153             : OperatorMatrices<Scalar>
     154         150 : construct_operator_matrices(const GreenTensor<Scalar> &green_tensor,
     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         184 :     auto get_matrices = [](auto basis, OperatorType type, std::initializer_list<int> m,
     159             :                            bool conjugate) {
     160         184 :         std::vector<Eigen::SparseMatrix<Scalar, Eigen::RowMajor>> matrices;
     161         184 :         matrices.reserve(m.size());
     162         184 :         int factor = conjugate ? -1 : 1;
     163        1288 :         std::transform(m.begin(), m.end(), std::back_inserter(matrices), [&](int q) {
     164           0 :             return (std::pow(factor, q) *
     165        1104 :                     OperatorAtom<Scalar>(basis, type, factor * q).get_matrix())
     166        1104 :                 .eval();
     167             :         });
     168         368 :         return matrices;
     169           0 :     };
     170             : 
     171         150 :     OperatorMatrices<Scalar> op;
     172             : 
     173             :     // Operator matrices for Rydberg-Rydberg interaction
     174         228 :     if (!green_tensor.get_spherical_entries(1, 1).empty() ||
     175          78 :         !green_tensor.get_spherical_entries(1, 2).empty()) {
     176          72 :         op.d1 = get_matrices(basis1, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, true);
     177             :     }
     178         228 :     if (!green_tensor.get_spherical_entries(1, 1).empty() ||
     179          78 :         !green_tensor.get_spherical_entries(2, 1).empty()) {
     180          72 :         op.d2 = get_matrices(basis2, OperatorType::ELECTRIC_DIPOLE, {-1, 0, +1}, false);
     181             :     }
     182         300 :     if (!green_tensor.get_spherical_entries(2, 2).empty() ||
     183         150 :         !green_tensor.get_spherical_entries(2, 1).empty()) {
     184          10 :         op.q1 = get_matrices(basis1, OperatorType::ELECTRIC_QUADRUPOLE, {-2, -1, 0, +1, +2}, true);
     185          10 :         op.q1.push_back(get_matrices(basis1, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, {0}, true)[0]);
     186             :     }
     187         300 :     if (!green_tensor.get_spherical_entries(2, 2).empty() ||
     188         150 :         !green_tensor.get_spherical_entries(1, 2).empty()) {
     189          10 :         op.q2 = get_matrices(basis2, OperatorType::ELECTRIC_QUADRUPOLE, {-2, -1, 0, +1, +2}, false);
     190          10 :         op.q2.push_back(
     191          20 :             get_matrices(basis2, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, {0}, false)[0]);
     192             :     }
     193             : 
     194         300 :     return op;
     195           0 : }
     196             : 
     197             : // "overloaded" pattern for std::visit
     198             : template <class... Ts>
     199             : struct overloaded : Ts... {
     200             :     using Ts::operator()...;
     201             : };
     202             : template <class... Ts>
     203             : overloaded(Ts...) -> overloaded<Ts...>;
     204             : 
     205             : template <typename Scalar>
     206          82 : SystemPair<Scalar>::SystemPair(std::shared_ptr<const basis_t> basis)
     207          82 :     : System<SystemPair<Scalar>>(std::move(basis)) {}
     208             : 
     209             : template <typename Scalar>
     210          40 : SystemPair<Scalar> &SystemPair<Scalar>::set_interaction_order(int value) {
     211          40 :     this->hamiltonian_requires_construction = true;
     212             : 
     213          40 :     if (value < 3 || value > 5) {
     214           0 :         throw std::invalid_argument("The order must be 3, 4, or 5.");
     215             :     }
     216             : 
     217          40 :     if (user_defined_green_tensor) {
     218           0 :         throw std::invalid_argument(
     219             :             "Cannot set interaction order if a user-defined green tensor is set.");
     220             :     }
     221             : 
     222          40 :     interaction_order = value;
     223             : 
     224          40 :     return *this;
     225             : }
     226             : 
     227             : template <typename Scalar>
     228          69 : SystemPair<Scalar> &SystemPair<Scalar>::set_distance_vector(const std::array<real_t, 3> &vector) {
     229          69 :     this->hamiltonian_requires_construction = true;
     230             : 
     231          69 :     if (!traits::NumTraits<Scalar>::is_complex_v && distance_vector[1] != 0) {
     232           0 :         throw std::invalid_argument(
     233             :             "The distance vector must not have a y-component if the scalar type is real.");
     234             :     }
     235             : 
     236          69 :     if (user_defined_green_tensor) {
     237           0 :         throw std::invalid_argument(
     238             :             "Cannot set distance vector if a user-defined green tensor is set.");
     239             :     }
     240             : 
     241          69 :     distance_vector = vector;
     242             : 
     243          69 :     return *this;
     244             : }
     245             : 
     246             : template <typename Scalar>
     247             : SystemPair<Scalar> &
     248           3 : SystemPair<Scalar>::set_green_tensor(std::shared_ptr<const GreenTensor<Scalar>> &green_tensor) {
     249           3 :     this->hamiltonian_requires_construction = true;
     250             : 
     251           6 :     if (std::isfinite(distance_vector[0]) && std::isfinite(distance_vector[1]) &&
     252           3 :         std::isfinite(distance_vector[2])) {
     253           0 :         throw std::invalid_argument("Cannot set green tensor if a finite distance vector is set.");
     254             :     }
     255             : 
     256           3 :     user_defined_green_tensor = green_tensor;
     257             : 
     258           3 :     return *this;
     259             : }
     260             : 
     261             : template <typename Scalar>
     262         150 : void SystemPair<Scalar>::construct_hamiltonian() const {
     263         150 :     auto basis = this->hamiltonian->get_basis();
     264         150 :     auto basis1 = basis->get_basis1();
     265         150 :     auto basis2 = basis->get_basis2();
     266             : 
     267         150 :     std::shared_ptr<const GreenTensor<Scalar>> green_tensor_ptr;
     268         150 :     if (user_defined_green_tensor) {
     269           3 :         green_tensor_ptr = user_defined_green_tensor;
     270             :     } else {
     271         147 :         green_tensor_ptr = std::make_shared<const GreenTensor<Scalar>>(
     272         147 :             construct_green_tensor<Scalar>(distance_vector, interaction_order));
     273             :     }
     274             : 
     275         150 :     auto op = construct_operator_matrices(*green_tensor_ptr, basis1, basis2);
     276             : 
     277             :     // Construct the unperturbed Hamiltonian
     278         150 :     this->hamiltonian = std::make_unique<OperatorPair<Scalar>>(basis, OperatorType::ENERGY);
     279         150 :     this->hamiltonian_is_diagonal = true;
     280         150 :     bool sort_by_quantum_number_f = basis->has_quantum_number_f();
     281         150 :     bool sort_by_quantum_number_m = basis->has_quantum_number_m();
     282         150 :     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         150 :     auto energies = this->hamiltonian->get_matrix().diagonal().real();
     287             : 
     288             :     // Helper function for adding Rydberg-Rydberg interaction
     289         750 :     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         876 :         for (const auto &entry : entries) {
     293         276 :             std::visit(
     294             :                 overloaded{
     295         276 :                     [this, &basis, &sort_by_quantum_number_m, &op1, &op2,
     296             :                      &delta](const typename GreenTensor<Scalar>::ConstantEntry &ce) {
     297         276 :                         this->hamiltonian->get_matrix() += ce.val() *
     298         276 :                             utils::calculate_tensor_product(basis, basis, op1[ce.row()],
     299         276 :                                                             op2[ce.col()]);
     300         276 :                         if (ce.row() != ce.col() + delta) {
     301           0 :                             sort_by_quantum_number_m = false;
     302             :                         }
     303             :                     },
     304             : 
     305           0 :                     [this, &basis, &energies, &sort_by_quantum_number_m, &op1, &op2,
     306             :                      &delta](const typename GreenTensor<Scalar>::OmegaDependentEntry &oe) {
     307           0 :                         auto tensor_product = utils::calculate_tensor_product(
     308           0 :                             basis, basis, op1[oe.row()], op2[oe.col()]);
     309           0 :                         for (int k = 0; k < tensor_product.outerSize(); ++k) {
     310           0 :                             for (typename Eigen::SparseMatrix<
     311           0 :                                      Scalar, Eigen::RowMajor>::InnerIterator it(tensor_product, k);
     312           0 :                                  it; ++it) {
     313           0 :                                 it.valueRef() *= oe.val(energies(it.row()) - energies(it.col()));
     314             :                             }
     315             :                         }
     316           0 :                         this->hamiltonian->get_matrix() += tensor_product;
     317           0 :                         if (oe.row() != oe.col() + delta) {
     318           0 :                             sort_by_quantum_number_m = false;
     319             :                         }
     320           0 :                     }},
     321             :                 entry);
     322             : 
     323         276 :             this->hamiltonian_is_diagonal = false;
     324         276 :             sort_by_quantum_number_f = false;
     325             :         }
     326             :     };
     327             : 
     328             :     // Dipole-dipole interaction
     329         150 :     add_interaction(green_tensor_ptr->get_spherical_entries(1, 1), op.d1, op.d2, 0);
     330             : 
     331             :     // Dipole-quadrupole interaction
     332         150 :     add_interaction(green_tensor_ptr->get_spherical_entries(1, 2), op.d1, op.q2, -1);
     333             : 
     334             :     // Quadrupole-dipole interaction
     335         150 :     add_interaction(green_tensor_ptr->get_spherical_entries(2, 1), op.q1, op.d2, +1);
     336             : 
     337             :     // Quadrupole-quadrupole interaction
     338         150 :     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         150 :     this->blockdiagonalizing_labels.clear();
     342         150 :     if (sort_by_quantum_number_f) {
     343           0 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_F);
     344             :     }
     345         150 :     if (sort_by_quantum_number_m) {
     346         150 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_M);
     347             :     }
     348         150 :     if (sort_by_parity) {
     349           0 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_PARITY);
     350             :     }
     351         150 : }
     352             : 
     353             : // Explicit instantiations
     354             : template class SystemPair<double>;
     355             : template class SystemPair<std::complex<double>>;
     356             : } // namespace pairinteraction