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 <string>
      19             : #include <vector>
      20             : 
      21             : namespace pairinteraction {
      22             : template <typename Scalar>
      23             : // NOLINTNEXTLINE(cert-err58-cpp)
      24             : const std::map<SorterType, std::string> BasisPair<Scalar>::sorter_type_to_quantum_number_name{
      25             :     {SorterType::QUANTUM_NUMBER_F, "f"},
      26             :     {SorterType::QUANTUM_NUMBER_M, "m"},
      27             :     {SorterType::PARITY, "product_of_parities"},
      28             : };
      29             : 
      30             : template <typename Scalar>
      31         552 : BasisPair<Scalar>::BasisPair(Private /*unused*/, ketvec_t &&kets,
      32             :                              map_range_t &&state_index1_to_state_index_range2,
      33             :                              map_indices_t &&state_indices_to_ket_index,
      34             :                              std::shared_ptr<const BasisAtom<Scalar>> basis1,
      35             :                              std::shared_ptr<const BasisAtom<Scalar>> basis2)
      36         552 :     : Basis<BasisPair<Scalar>>(std::move(kets)),
      37         552 :       state_index1_to_state_index_range2(std::move(state_index1_to_state_index_range2)),
      38         552 :       state_indices_to_ket_index(std::move(state_indices_to_ket_index)), basis1(std::move(basis1)),
      39        1104 :       basis2(std::move(basis2)) {}
      40             : 
      41             : template <typename Scalar>
      42             : const typename BasisPair<Scalar>::range_t &
      43     7081626 : BasisPair<Scalar>::get_index_range(size_t state_index1) const {
      44     7081626 :     return state_index1_to_state_index_range2.at(state_index1);
      45             : }
      46             : 
      47             : template <typename Scalar>
      48        1095 : std::shared_ptr<const BasisAtom<Scalar>> BasisPair<Scalar>::get_basis1() const {
      49        1095 :     return basis1;
      50             : }
      51             : 
      52             : template <typename Scalar>
      53        1095 : std::shared_ptr<const BasisAtom<Scalar>> BasisPair<Scalar>::get_basis2() const {
      54        1095 :     return basis2;
      55             : }
      56             : 
      57             : template <typename Scalar>
      58     8634098 : int BasisPair<Scalar>::get_ket_index_from_tuple(size_t state_index1, size_t state_index2) const {
      59     8634098 :     if (!state_indices_to_ket_index.contains({state_index1, state_index2})) {
      60     1335392 :         return -1;
      61             :     }
      62     7218868 :     return state_indices_to_ket_index.at({state_index1, state_index2});
      63             : }
      64             : 
      65             : template <typename Scalar>
      66             : std::shared_ptr<const typename BasisPair<Scalar>::Type>
      67           7 : BasisPair<Scalar>::merge(std::shared_ptr<const Type> other) const {
      68           7 :     if (basis1 != other->basis1 || basis2 != other->basis2) {
      69           1 :         throw std::invalid_argument(
      70             :             "Cannot merge two pair bases whose underlying atomic bases are not identical objects. "
      71             :             "Both pair bases must be constructed from the same pair of SystemAtom instances, "
      72             :             "and these instances must not have been changed in between.");
      73             :     }
      74           6 :     if (!this->is_canonical() || !other->is_canonical()) {
      75           0 :         throw std::invalid_argument(
      76             :             "Cannot merge non-canonical bases (i.e., bases with non-identity coefficients). "
      77             :             "Canonicalize the bases first.");
      78             :     }
      79             : 
      80           6 :     std::map<std::vector<size_t>, std::shared_ptr<const ket_t>> ket_by_state_indices;
      81          18 :     for (const auto &basis : {this, other.get()}) {
      82         364 :         for (const auto &[state_indices, ket_index] : basis->state_indices_to_ket_index) {
      83         352 :             if (!ket_by_state_indices.contains(state_indices)) {
      84         288 :                 ket_by_state_indices.emplace(state_indices, basis->kets.at(ket_index));
      85             :             }
      86             :         }
      87             :     }
      88             : 
      89           6 :     ketvec_t merged_kets;
      90           6 :     merged_kets.reserve(ket_by_state_indices.size());
      91           6 :     map_indices_t merged_state_indices;
      92           6 :     merged_state_indices.reserve(ket_by_state_indices.size());
      93         294 :     for (const auto &[state_indices, ket] : ket_by_state_indices) {
      94         288 :         merged_state_indices.try_emplace(state_indices, merged_kets.size());
      95         288 :         merged_kets.push_back(ket);
      96             :     }
      97             : 
      98           6 :     const size_t number_of_states1 = basis1->get_number_of_states();
      99           6 :     const size_t number_of_states2 = basis2->get_number_of_states();
     100           6 :     std::vector<size_t> minimum_indices2(number_of_states1, number_of_states2);
     101           6 :     std::vector<size_t> maximum_indices2(number_of_states1, 0);
     102           6 :     std::vector<bool> has_index(number_of_states1, false);
     103         294 :     for (const auto &[state_indices, ket_index] : merged_state_indices) {
     104             :         static_cast<void>(ket_index);
     105         288 :         const size_t state_index1 = state_indices[0];
     106         288 :         const size_t state_index2 = state_indices[1];
     107         288 :         minimum_indices2[state_index1] = std::min(minimum_indices2[state_index1], state_index2);
     108         288 :         maximum_indices2[state_index1] = std::max(maximum_indices2[state_index1], state_index2);
     109         288 :         has_index[state_index1] = true;
     110             :     }
     111             : 
     112           6 :     map_range_t merged_ranges;
     113           6 :     merged_ranges.reserve(number_of_states1);
     114         546 :     for (size_t state_index1 = 0; state_index1 < number_of_states1; ++state_index1) {
     115         540 :         if (has_index[state_index1]) {
     116          40 :             merged_ranges.try_emplace(
     117             :                 state_index1,
     118          80 :                 range_t(minimum_indices2[state_index1], maximum_indices2[state_index1] + 1));
     119             :         } else {
     120         500 :             merged_ranges.try_emplace(state_index1, range_t(0, 0));
     121             :         }
     122             :     }
     123             : 
     124           6 :     return std::make_shared<const Type>(Private(), std::move(merged_kets), std::move(merged_ranges),
     125          18 :                                         std::move(merged_state_indices), basis1, basis2);
     126           6 : }
     127             : 
     128             : template <typename Scalar>
     129             : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
     130         210 : BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const Type> final_state, OperatorType type1,
     131             :                                        OperatorType type2, int q1, int q2) const {
     132         210 :     auto initial1 = this->get_basis1();
     133         210 :     auto initial2 = this->get_basis2();
     134         210 :     auto final1 = final_state->get_basis1();
     135         210 :     auto final2 = final_state->get_basis2();
     136             : 
     137         210 :     auto matrix_elements1 = initial1->get_database().get_matrix_elements_in_canonical_basis(
     138             :         initial1, final1, type1, q1);
     139         210 :     matrix_elements1 =
     140         210 :         final1->get_coefficients().adjoint() * matrix_elements1 * initial1->get_coefficients();
     141         210 :     auto matrix_elements2 = initial2->get_database().get_matrix_elements_in_canonical_basis(
     142             :         initial2, final2, type2, q2);
     143         210 :     matrix_elements2 =
     144         210 :         final2->get_coefficients().adjoint() * matrix_elements2 * initial2->get_coefficients();
     145         210 :     auto matrix_elements = utils::calculate_tensor_product_in_canonical_basis(
     146             :         this->shared_from_this(), final_state, matrix_elements1, matrix_elements2);
     147         210 :     assert(static_cast<size_t>(matrix_elements.rows()) == final_state->get_number_of_kets());
     148         210 :     assert(static_cast<size_t>(matrix_elements.cols()) == this->get_number_of_kets());
     149             : 
     150         420 :     return final_state->get_coefficients().adjoint() * matrix_elements * this->get_coefficients();
     151         210 : }
     152             : 
     153             : // Explicit instantiations
     154             : template class BasisPair<double>;
     155             : template class BasisPair<std::complex<double>>;
     156             : } // namespace pairinteraction