Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/basis/BasisPair.hpp"
       5             : 
       6             : #include "pairinteraction/basis/BasisAtom.hpp"
       7             : #include "pairinteraction/database/Database.hpp"
       8             : #include "pairinteraction/ket/KetAtom.hpp"
       9             : #include "pairinteraction/ket/KetPair.hpp"
      10             : #include "pairinteraction/utils/Range.hpp"
      11             : #include "pairinteraction/utils/tensor.hpp"
      12             : 
      13             : #include <algorithm>
      14             : #include <cassert>
      15             : #include <map>
      16             : #include <memory>
      17             : #include <stdexcept>
      18             : #include <vector>
      19             : 
      20             : namespace pairinteraction {
      21             : template <typename Scalar>
      22         513 : BasisPair<Scalar>::BasisPair(Private /*unused*/, ketvec_t &&kets,
      23             :                              map_range_t &&map_range_of_state_index2,
      24             :                              map_indices_t &&state_indices_to_ket_index,
      25             :                              std::shared_ptr<const BasisAtom<Scalar>> basis1,
      26             :                              std::shared_ptr<const BasisAtom<Scalar>> basis2)
      27         513 :     : Basis<BasisPair<Scalar>>(std::move(kets)),
      28         513 :       map_range_of_state_index2(std::move(map_range_of_state_index2)),
      29         513 :       state_indices_to_ket_index(std::move(state_indices_to_ket_index)), basis1(std::move(basis1)),
      30        1026 :       basis2(std::move(basis2)) {}
      31             : 
      32             : template <typename Scalar>
      33             : const typename BasisPair<Scalar>::range_t &
      34     6238765 : BasisPair<Scalar>::get_index_range(size_t state_index1) const {
      35     6238765 :     return map_range_of_state_index2.at(state_index1);
      36             : }
      37             : 
      38             : template <typename Scalar>
      39        1008 : std::shared_ptr<const BasisAtom<Scalar>> BasisPair<Scalar>::get_basis1() const {
      40        1008 :     return basis1;
      41             : }
      42             : 
      43             : template <typename Scalar>
      44        1008 : std::shared_ptr<const BasisAtom<Scalar>> BasisPair<Scalar>::get_basis2() const {
      45        1008 :     return basis2;
      46             : }
      47             : 
      48             : template <typename Scalar>
      49     7105169 : int BasisPair<Scalar>::get_ket_index_from_tuple(size_t state_index1, size_t state_index2) const {
      50     7105169 :     if (!state_indices_to_ket_index.contains({state_index1, state_index2})) {
      51      993846 :         return -1;
      52             :     }
      53     6072958 :     return state_indices_to_ket_index.at({state_index1, state_index2});
      54             : }
      55             : 
      56             : template <typename Scalar>
      57             : std::shared_ptr<const typename BasisPair<Scalar>::Type>
      58           7 : BasisPair<Scalar>::merge(std::shared_ptr<const Type> other) const {
      59           7 :     if (basis1 != other->basis1 || basis2 != other->basis2) {
      60           1 :         throw std::invalid_argument(
      61             :             "Cannot merge two pair bases whose underlying atomic bases are not identical objects. "
      62             :             "Both pair bases must be constructed from the same pair of SystemAtom instances, "
      63             :             "and these instances must not have been changed in between.");
      64             :     }
      65           6 :     if (!this->is_canonical() || !other->is_canonical()) {
      66           0 :         throw std::invalid_argument(
      67             :             "Cannot merge non-canonical bases (i.e., bases with non-identity coefficients). "
      68             :             "Canonicalize the bases first.");
      69             :     }
      70             : 
      71           6 :     std::map<std::vector<size_t>, std::shared_ptr<const ket_t>> ket_by_state_indices;
      72          18 :     for (const auto &basis : {this, other.get()}) {
      73         364 :         for (const auto &[state_indices, ket_index] : basis->state_indices_to_ket_index) {
      74         352 :             if (!ket_by_state_indices.contains(state_indices)) {
      75         288 :                 ket_by_state_indices.emplace(state_indices, basis->kets.at(ket_index));
      76             :             }
      77             :         }
      78             :     }
      79             : 
      80           6 :     ketvec_t merged_kets;
      81           6 :     merged_kets.reserve(ket_by_state_indices.size());
      82           6 :     map_indices_t merged_state_indices;
      83           6 :     merged_state_indices.reserve(ket_by_state_indices.size());
      84         294 :     for (const auto &[state_indices, ket] : ket_by_state_indices) {
      85         288 :         merged_state_indices.try_emplace(state_indices, merged_kets.size());
      86         288 :         merged_kets.push_back(ket);
      87             :     }
      88             : 
      89           6 :     const size_t number_of_states1 = basis1->get_number_of_states();
      90           6 :     const size_t number_of_states2 = basis2->get_number_of_states();
      91           6 :     std::vector<size_t> minimum_indices2(number_of_states1, number_of_states2);
      92           6 :     std::vector<size_t> maximum_indices2(number_of_states1, 0);
      93           6 :     std::vector<bool> has_index(number_of_states1, false);
      94         294 :     for (const auto &[state_indices, ket_index] : merged_state_indices) {
      95             :         static_cast<void>(ket_index);
      96         288 :         const size_t state_index1 = state_indices[0];
      97         288 :         const size_t state_index2 = state_indices[1];
      98         288 :         minimum_indices2[state_index1] = std::min(minimum_indices2[state_index1], state_index2);
      99         288 :         maximum_indices2[state_index1] = std::max(maximum_indices2[state_index1], state_index2);
     100         288 :         has_index[state_index1] = true;
     101             :     }
     102             : 
     103           6 :     map_range_t merged_ranges;
     104           6 :     merged_ranges.reserve(number_of_states1);
     105         546 :     for (size_t state_index1 = 0; state_index1 < number_of_states1; ++state_index1) {
     106         540 :         if (has_index[state_index1]) {
     107          40 :             merged_ranges.try_emplace(
     108             :                 state_index1,
     109          80 :                 range_t(minimum_indices2[state_index1], maximum_indices2[state_index1] + 1));
     110             :         } else {
     111         500 :             merged_ranges.try_emplace(state_index1, range_t(0, 0));
     112             :         }
     113             :     }
     114             : 
     115           6 :     return std::make_shared<const Type>(Private(), std::move(merged_kets), std::move(merged_ranges),
     116          18 :                                         std::move(merged_state_indices), basis1, basis2);
     117           6 : }
     118             : 
     119             : template <typename Scalar>
     120             : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
     121         204 : BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const Type> final_state, OperatorType type1,
     122             :                                        OperatorType type2, int q1, int q2) const {
     123         204 :     auto initial1 = this->get_basis1();
     124         204 :     auto initial2 = this->get_basis2();
     125         204 :     auto final1 = final_state->get_basis1();
     126         204 :     auto final2 = final_state->get_basis2();
     127             : 
     128         204 :     auto matrix_elements1 = initial1->get_database().get_matrix_elements_in_canonical_basis(
     129             :         initial1, final1, type1, q1);
     130         204 :     matrix_elements1 =
     131         204 :         final1->get_coefficients().adjoint() * matrix_elements1 * initial1->get_coefficients();
     132         204 :     auto matrix_elements2 = initial2->get_database().get_matrix_elements_in_canonical_basis(
     133             :         initial2, final2, type2, q2);
     134         204 :     matrix_elements2 =
     135         204 :         final2->get_coefficients().adjoint() * matrix_elements2 * initial2->get_coefficients();
     136         204 :     auto matrix_elements = utils::calculate_tensor_product_in_canonical_basis(
     137             :         this->shared_from_this(), final_state, matrix_elements1, matrix_elements2);
     138         204 :     assert(static_cast<size_t>(matrix_elements.rows()) == final_state->get_number_of_kets());
     139         204 :     assert(static_cast<size_t>(matrix_elements.cols()) == this->get_number_of_kets());
     140             : 
     141         408 :     return final_state->get_coefficients().adjoint() * matrix_elements * this->get_coefficients();
     142         204 : }
     143             : 
     144             : // Explicit instantiations
     145             : template class BasisPair<double>;
     146             : template class BasisPair<std::complex<double>>;
     147             : } // namespace pairinteraction