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/KetPair.hpp"
       9             : #include "pairinteraction/utils/Range.hpp"
      10             : #include "pairinteraction/utils/tensor.hpp"
      11             : 
      12             : #include <cassert>
      13             : #include <memory>
      14             : #include <vector>
      15             : 
      16             : namespace pairinteraction {
      17             : template <typename Scalar>
      18         483 : BasisPair<Scalar>::BasisPair(Private /*unused*/, ketvec_t &&kets,
      19             :                              map_range_t &&map_range_of_state_index2,
      20             :                              map_indices_t &&state_indices_to_ket_index,
      21             :                              std::shared_ptr<const BasisAtom<Scalar>> basis1,
      22             :                              std::shared_ptr<const BasisAtom<Scalar>> basis2)
      23         483 :     : Basis<BasisPair<Scalar>>(std::move(kets)),
      24         483 :       map_range_of_state_index2(std::move(map_range_of_state_index2)),
      25         483 :       state_indices_to_ket_index(std::move(state_indices_to_ket_index)), basis1(std::move(basis1)),
      26         966 :       basis2(std::move(basis2)) {}
      27             : 
      28             : template <typename Scalar>
      29             : const typename BasisPair<Scalar>::range_t &
      30     7365236 : BasisPair<Scalar>::get_index_range(size_t state_index1) const {
      31     7365236 :     return map_range_of_state_index2.at(state_index1);
      32             : }
      33             : 
      34             : template <typename Scalar>
      35         976 : std::shared_ptr<const BasisAtom<Scalar>> BasisPair<Scalar>::get_basis1() const {
      36         976 :     return basis1;
      37             : }
      38             : 
      39             : template <typename Scalar>
      40         976 : std::shared_ptr<const BasisAtom<Scalar>> BasisPair<Scalar>::get_basis2() const {
      41         976 :     return basis2;
      42             : }
      43             : 
      44             : template <typename Scalar>
      45     8573111 : int BasisPair<Scalar>::get_ket_index_from_tuple(size_t state_index1, size_t state_index2) const {
      46     8573111 :     if (!state_indices_to_ket_index.contains({state_index1, state_index2})) {
      47     1177519 :         return -1;
      48             :     }
      49     7395592 :     return state_indices_to_ket_index.at({state_index1, state_index2});
      50             : }
      51             : 
      52             : template <typename Scalar>
      53             : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
      54         188 : BasisPair<Scalar>::get_matrix_elements(std::shared_ptr<const Type> final_state, OperatorType type1,
      55             :                                        OperatorType type2, int q1, int q2) const {
      56         188 :     auto initial1 = this->get_basis1();
      57         188 :     auto initial2 = this->get_basis2();
      58         188 :     auto final1 = final_state->get_basis1();
      59         188 :     auto final2 = final_state->get_basis2();
      60             : 
      61         188 :     auto matrix_elements1 = initial1->get_database().get_matrix_elements_in_canonical_basis(
      62             :         initial1, final1, type1, q1);
      63         188 :     matrix_elements1 =
      64         188 :         final1->get_coefficients().adjoint() * matrix_elements1 * initial1->get_coefficients();
      65         188 :     auto matrix_elements2 = initial2->get_database().get_matrix_elements_in_canonical_basis(
      66             :         initial2, final2, type2, q2);
      67         188 :     matrix_elements2 =
      68         188 :         final2->get_coefficients().adjoint() * matrix_elements2 * initial2->get_coefficients();
      69         188 :     auto matrix_elements = utils::calculate_tensor_product_in_canonical_basis(
      70             :         this->shared_from_this(), final_state, matrix_elements1, matrix_elements2);
      71         188 :     assert(static_cast<size_t>(matrix_elements.rows()) == final_state->get_number_of_kets());
      72         188 :     assert(static_cast<size_t>(matrix_elements.cols()) == this->get_number_of_kets());
      73             : 
      74         376 :     return final_state->get_coefficients().adjoint() * matrix_elements * this->get_coefficients();
      75         188 : }
      76             : 
      77             : // Explicit instantiations
      78             : template class BasisPair<double>;
      79             : template class BasisPair<std::complex<double>>;
      80             : } // namespace pairinteraction