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           9 : 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           9 :     GreenTensor<Scalar> green_tensor;
      53             : 
      54             :     // Normalize the distance vector, return zero green tensor if the distance is infinity
      55          18 :     Eigen::Map<const Eigen::Vector3<real_t>> vector_map(distance_vector.data(),
      56           9 :                                                         distance_vector.size());
      57           9 :     real_t distance = vector_map.norm();
      58          18 :     SPDLOG_DEBUG("Interatomic distance: {}", distance);
      59           9 :     if (distance == std::numeric_limits<real_t>::infinity()) {
      60           0 :         return green_tensor;
      61             :     }
      62           9 :     Eigen::Vector3<real_t> vector_normalized = vector_map / distance;
      63             : 
      64             :     // Dyadic green function of dipole-dipole interaction
      65           9 :     if (interaction_order >= 3) {
      66           9 :         Eigen::Matrix3<Scalar> entries = Eigen::Matrix3<real_t>::Identity() -
      67           9 :             3 * vector_normalized * vector_normalized.transpose();
      68           9 :         entries /= std::pow(distance, 3);
      69             : 
      70           9 :         green_tensor.set_entries(1, 1, entries.template cast<Scalar>());
      71             :     }
      72             : 
      73             :     // Dyadic green function of dipole-quadrupole interaction
      74           9 :     if (interaction_order >= 4) {
      75           0 :         Eigen::Matrix<real_t, 3, 9> entries = Eigen::Matrix<real_t, 3, 9>::Zero();
      76           0 :         for (Eigen::Index q = 0; q < 3; ++q) {
      77           0 :             Eigen::Index row = q;
      78           0 :             for (Eigen::Index j = 0; j < 3; ++j) {
      79           0 :                 for (Eigen::Index i = 0; i < 3; ++i) {
      80           0 :                     Eigen::Index col = 3 * j + i;
      81           0 :                     entries(row, col) +=
      82           0 :                         15 * vector_normalized[q] * vector_normalized[j] * vector_normalized[i];
      83           0 :                     if (i == j) {
      84           0 :                         entries(row, col) += -3 * vector_normalized[q];
      85             :                     }
      86           0 :                     if (i == q) {
      87           0 :                         entries(row, col) += -3 * vector_normalized[j];
      88             :                     }
      89           0 :                     if (j == q) {
      90           0 :                         entries(row, col) += -3 * vector_normalized[i];
      91             :                     }
      92             :                 }
      93             :             }
      94             :         }
      95           0 :         entries /= std::pow(distance, 4);
      96             : 
      97           0 :         green_tensor.set_entries(1, 2, entries.template cast<Scalar>());
      98             :     }
      99             : 
     100             :     // Dyadic green function of quadrupole-dipole interaction
     101           9 :     if (interaction_order >= 4) {
     102           0 :         Eigen::Matrix<real_t, 9, 3> entries = Eigen::Matrix<real_t, 9, 3>::Zero();
     103           0 :         for (Eigen::Index q = 0; q < 3; ++q) {
     104           0 :             for (Eigen::Index j = 0; j < 3; ++j) {
     105           0 :                 Eigen::Index row = 3 * q + j;
     106           0 :                 for (Eigen::Index i = 0; i < 3; ++i) {
     107           0 :                     Eigen::Index col = i;
     108           0 :                     entries(row, col) +=
     109           0 :                         -15 * vector_normalized[q] * vector_normalized[j] * vector_normalized[i];
     110           0 :                     if (i == j) {
     111           0 :                         entries(row, col) += 3 * vector_normalized[q];
     112             :                     }
     113           0 :                     if (i == q) {
     114           0 :                         entries(row, col) += 3 * vector_normalized[j];
     115             :                     }
     116           0 :                     if (j == q) {
     117           0 :                         entries(row, col) += 3 * vector_normalized[i];
     118             :                     }
     119             :                 }
     120             :             }
     121             :         }
     122           0 :         entries /= std::pow(distance, 4);
     123             : 
     124           0 :         green_tensor.set_entries(2, 1, entries.template cast<Scalar>());
     125             :     }
     126             : 
     127             :     // Dyadic green function of quadrupole-quadrupole interaction
     128           9 :     if (interaction_order >= 5) {
     129           0 :         SPDLOG_WARN("Quadrupole-quadrupole interaction is considered but "
     130             :                     "not dipole-octupole interaction although this interaction would be "
     131             :                     "of the same order. We plan to implement dipole-octupole interaction "
     132             :                     "in the future.");
     133             : 
     134           0 :         Eigen::Matrix<real_t, 9, 9> entries = Eigen::Matrix<real_t, 9, 9>::Zero();
     135           0 :         for (Eigen::Index q = 0; q < 3; ++q) {
     136           0 :             for (Eigen::Index j = 0; j < 3; ++j) {
     137           0 :                 Eigen::Index row = 3 * q + j;
     138           0 :                 for (Eigen::Index i = 0; i < 3; ++i) {
     139           0 :                     for (Eigen::Index k = 0; k < 3; ++k) {
     140           0 :                         Eigen::Index col = 3 * i + k;
     141           0 :                         entries(row, col) += 105 * vector_normalized[q] * vector_normalized[j] *
     142           0 :                             vector_normalized[i] * vector_normalized[k];
     143           0 :                         if (i == j) {
     144           0 :                             entries(row, col) += -15 * vector_normalized[q] * vector_normalized[k];
     145             :                         }
     146           0 :                         if (i == q) {
     147           0 :                             entries(row, col) += -15 * vector_normalized[j] * vector_normalized[k];
     148             :                         }
     149           0 :                         if (j == q) {
     150           0 :                             entries(row, col) += -15 * vector_normalized[i] * vector_normalized[k];
     151             :                         }
     152           0 :                         if (k == q) {
     153           0 :                             entries(row, col) += -15 * vector_normalized[j] * vector_normalized[i];
     154             :                         }
     155           0 :                         if (k == j) {
     156           0 :                             entries(row, col) += -15 * vector_normalized[q] * vector_normalized[i];
     157             :                         }
     158           0 :                         if (k == i) {
     159           0 :                             entries(row, col) += -15 * vector_normalized[q] * vector_normalized[j];
     160             :                         }
     161           0 :                         if (q == k && i == j) {
     162           0 :                             entries(row, col) += 3;
     163             :                         }
     164           0 :                         if (i == k && j == q) {
     165           0 :                             entries(row, col) += 3;
     166             :                         }
     167           0 :                         if (j == k && i == q) {
     168           0 :                             entries(row, col) += 3;
     169             :                         }
     170             :                     }
     171             :                 }
     172             :             }
     173             :         }
     174           0 :         entries /= std::pow(distance, 5);
     175             : 
     176           0 :         green_tensor.set_entries(2, 2, entries.template cast<Scalar>());
     177             :     }
     178             : 
     179           9 :     return green_tensor;
     180           0 : }
     181             : 
     182             : template <typename Scalar>
     183             : OperatorMatrices<Scalar>
     184           9 : construct_operator_matrices(const GreenTensor<Scalar> &green_tensor,
     185             :                             const std::shared_ptr<const BasisAtom<Scalar>> &basis1,
     186             :                             const std::shared_ptr<const BasisAtom<Scalar>> &basis2) {
     187           9 :     OperatorMatrices<Scalar> op;
     188             : 
     189           9 :     if (!green_tensor.get_entries(1, 1).empty() || !green_tensor.get_entries(1, 2).empty()) {
     190           9 :         op.d1.push_back(
     191          17 :             -OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_DIPOLE, 1).get_matrix());
     192          18 :         op.d1.push_back(
     193           9 :             OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_DIPOLE, 0).get_matrix());
     194           9 :         op.d1.push_back(
     195          18 :             -OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_DIPOLE, -1).get_matrix());
     196             :     }
     197             : 
     198           9 :     if (!green_tensor.get_entries(1, 1).empty() || !green_tensor.get_entries(2, 1).empty()) {
     199          18 :         op.d2.push_back(
     200           9 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_DIPOLE, -1).get_matrix());
     201          18 :         op.d2.push_back(
     202           9 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_DIPOLE, 0).get_matrix());
     203          18 :         op.d2.push_back(
     204           9 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_DIPOLE, 1).get_matrix());
     205             :     }
     206             : 
     207           9 :     if (!green_tensor.get_entries(2, 2).empty() || !green_tensor.get_entries(2, 1).empty()) {
     208           0 :         op.q1.push_back(
     209           0 :             OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_QUADRUPOLE, 2).get_matrix());
     210           0 :         op.q1.push_back(
     211           0 :             -OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_QUADRUPOLE, 1).get_matrix());
     212           0 :         op.q1.push_back(
     213           0 :             OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_QUADRUPOLE, 0).get_matrix());
     214           0 :         op.q1.push_back(
     215           0 :             -OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_QUADRUPOLE, -1).get_matrix());
     216           0 :         op.q1.push_back(
     217           0 :             OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_QUADRUPOLE, -2).get_matrix());
     218           0 :         op.q1.push_back(
     219           0 :             OperatorAtom<Scalar>(basis1, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0).get_matrix());
     220             :     }
     221             : 
     222           9 :     if (!green_tensor.get_entries(2, 2).empty() || !green_tensor.get_entries(1, 2).empty()) {
     223           0 :         op.q2.push_back(
     224           0 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_QUADRUPOLE, -2).get_matrix());
     225           0 :         op.q2.push_back(
     226           0 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_QUADRUPOLE, -1).get_matrix());
     227           0 :         op.q2.push_back(
     228           0 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_QUADRUPOLE, 0).get_matrix());
     229           0 :         op.q2.push_back(
     230           0 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_QUADRUPOLE, 1).get_matrix());
     231           0 :         op.q2.push_back(
     232           0 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_QUADRUPOLE, 2).get_matrix());
     233           0 :         op.q2.push_back(
     234           0 :             OperatorAtom<Scalar>(basis2, OperatorType::ELECTRIC_QUADRUPOLE_ZERO, 0).get_matrix());
     235             :     }
     236             : 
     237           9 :     return op;
     238           0 : }
     239             : 
     240             : template <typename Scalar>
     241           9 : SystemPair<Scalar>::SystemPair(std::shared_ptr<const basis_t> basis)
     242           9 :     : System<SystemPair<Scalar>>(std::move(basis)) {}
     243             : 
     244             : template <typename Scalar>
     245           0 : SystemPair<Scalar> &SystemPair<Scalar>::set_interaction_order(int value) {
     246           0 :     this->hamiltonian_requires_construction = true;
     247           0 :     if (value < 3 || value > 5) {
     248           0 :         throw std::invalid_argument("The order must be 3, 4, or 5.");
     249             :     }
     250           0 :     interaction_order = value;
     251           0 :     return *this;
     252             : }
     253             : 
     254             : template <typename Scalar>
     255           9 : SystemPair<Scalar> &SystemPair<Scalar>::set_distance_vector(const std::array<real_t, 3> &vector) {
     256           9 :     this->hamiltonian_requires_construction = true;
     257             : 
     258           9 :     constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
     259           9 :     if (!traits::NumTraits<Scalar>::is_complex_v &&
     260           9 :         std::abs(distance_vector[1]) > numerical_precision) {
     261           0 :         throw std::invalid_argument(
     262             :             "The distance vector must not have a y-component if the scalar type is real.");
     263             :     }
     264             : 
     265           9 :     distance_vector = vector;
     266             : 
     267           9 :     return *this;
     268             : }
     269             : 
     270             : template <typename Scalar>
     271           9 : void SystemPair<Scalar>::construct_hamiltonian() const {
     272           9 :     auto basis = this->hamiltonian->get_basis();
     273           9 :     auto basis1 = basis->get_basis1();
     274           9 :     auto basis2 = basis->get_basis2();
     275             : 
     276           9 :     auto green_tensor = construct_green_tensor<Scalar>(distance_vector, interaction_order);
     277           9 :     auto op = construct_operator_matrices(green_tensor, basis1, basis2);
     278             : 
     279             :     // Construct the unperturbed Hamiltonian
     280           9 :     this->hamiltonian = std::make_unique<OperatorPair<Scalar>>(basis, OperatorType::ENERGY);
     281           9 :     this->hamiltonian_is_diagonal = true;
     282           9 :     bool sort_by_quantum_number_f = basis->has_quantum_number_f();
     283           9 :     bool sort_by_quantum_number_m = basis->has_quantum_number_m();
     284           9 :     bool sort_by_parity = basis->has_parity();
     285             : 
     286             :     // Dipole-dipole interaction
     287          36 :     for (const auto &entry : green_tensor.get_entries(1, 1)) {
     288          27 :         if (auto ce = std::get_if<typename GreenTensor<Scalar>::ConstantEntry>(&entry)) {
     289          27 :             this->hamiltonian->get_matrix() += ce->val() *
     290          27 :                 utils::calculate_tensor_product(basis, basis, op.d1[ce->row()], op.d2[ce->col()]);
     291          27 :             if (ce->row() != ce->col()) {
     292           0 :                 sort_by_quantum_number_m = false;
     293             :             }
     294          27 :             this->hamiltonian_is_diagonal = false;
     295          27 :             sort_by_quantum_number_f = false;
     296             :         } else {
     297           0 :             throw std::invalid_argument("Omega-dependent entries are not yet supported.");
     298             :         }
     299             :     }
     300             : 
     301             :     // Dipole-quadrupole interaction
     302           9 :     for (const auto &entry : green_tensor.get_entries(1, 2)) {
     303           0 :         if (auto ce = std::get_if<typename GreenTensor<Scalar>::ConstantEntry>(&entry)) {
     304           0 :             this->hamiltonian->get_matrix() += ce->val() *
     305           0 :                 utils::calculate_tensor_product(basis, basis, op.d1[ce->row()], op.q2[ce->col()]);
     306           0 :             if (ce->row() != ce->col() - 1) {
     307           0 :                 sort_by_quantum_number_m = false;
     308             :             }
     309           0 :             this->hamiltonian_is_diagonal = false;
     310           0 :             sort_by_quantum_number_f = false;
     311             :         } else {
     312           0 :             throw std::invalid_argument("Omega-dependent entries are not yet supported.");
     313             :         }
     314             :     }
     315             : 
     316             :     // Quadrupole-dipole interaction
     317           9 :     for (const auto &entry : green_tensor.get_entries(2, 1)) {
     318           0 :         if (auto ce = std::get_if<typename GreenTensor<Scalar>::ConstantEntry>(&entry)) {
     319           0 :             this->hamiltonian->get_matrix() += ce->val() *
     320           0 :                 utils::calculate_tensor_product(basis, basis, op.q1[ce->row()], op.d2[ce->col()]);
     321           0 :             if (ce->row() - 1 != ce->col()) {
     322           0 :                 sort_by_quantum_number_m = false;
     323             :             }
     324           0 :             this->hamiltonian_is_diagonal = false;
     325           0 :             sort_by_quantum_number_f = false;
     326             :         } else {
     327           0 :             throw std::invalid_argument("Omega-dependent entries are not yet supported.");
     328             :         }
     329             :     }
     330             : 
     331             :     // Quadrupole-quadrupole interaction
     332           9 :     for (const auto &entry : green_tensor.get_entries(2, 2)) {
     333           0 :         if (auto ce = std::get_if<typename GreenTensor<Scalar>::ConstantEntry>(&entry)) {
     334           0 :             this->hamiltonian->get_matrix() += ce->val() *
     335           0 :                 utils::calculate_tensor_product(basis, basis, op.q1[ce->row()], op.q2[ce->col()]);
     336           0 :             if (ce->row() != ce->col()) {
     337           0 :                 sort_by_quantum_number_m = false;
     338             :             }
     339           0 :             this->hamiltonian_is_diagonal = false;
     340           0 :             sort_by_quantum_number_f = false;
     341             :         } else {
     342           0 :             throw std::invalid_argument("Omega-dependent entries are not yet supported.");
     343             :         }
     344             :     }
     345             : 
     346             :     // Store which labels can be used to block-diagonalize the Hamiltonian
     347           9 :     this->blockdiagonalizing_labels.clear();
     348           9 :     if (sort_by_quantum_number_f) {
     349           0 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_F);
     350             :     }
     351           9 :     if (sort_by_quantum_number_m) {
     352           9 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_QUANTUM_NUMBER_M);
     353             :     }
     354           9 :     if (sort_by_parity) {
     355           0 :         this->blockdiagonalizing_labels.push_back(TransformationType::SORT_BY_PARITY);
     356             :     }
     357           9 : }
     358             : 
     359             : // Explicit instantiations
     360             : template class SystemPair<double>;
     361             : template class SystemPair<std::complex<double>>;
     362             : } // namespace pairinteraction