Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/basis/BasisAtom.hpp"
       5             : 
       6             : #include "pairinteraction/basis/BasisAtomCreator.hpp"
       7             : #include "pairinteraction/database/Database.hpp"
       8             : #include "pairinteraction/ket/KetAtom.hpp"
       9             : 
      10             : #include <cassert>
      11             : #include <map>
      12             : #include <stdexcept>
      13             : #include <string>
      14             : #include <unordered_set>
      15             : 
      16             : namespace pairinteraction {
      17             : template <typename Scalar>
      18             : // NOLINTNEXTLINE(cert-err58-cpp)
      19             : const std::map<SorterType, std::string> BasisAtom<Scalar>::sorter_type_to_quantum_number_name{
      20             :     {SorterType::QUANTUM_NUMBER_F, "f"},
      21             :     {SorterType::QUANTUM_NUMBER_M, "m"},
      22             :     {SorterType::PARITY, "parity"},
      23             : };
      24             : 
      25             : template <typename Scalar>
      26        1666 : BasisAtom<Scalar>::BasisAtom(Private /*unused*/, ketvec_t &&kets, std::string &&canonical_basis_id,
      27             :                              Database &database)
      28        1666 :     : Basis<BasisAtom<Scalar>>(std::move(kets)), canonical_basis_id(std::move(canonical_basis_id)),
      29        3332 :       database(database) {
      30       48993 :     for (size_t i = 0; i < this->kets.size(); ++i) {
      31       47327 :         ket_id_to_ket_index[this->kets[i]->get_id_in_database()] = i;
      32             :     }
      33        1666 : }
      34             : 
      35             : template <typename Scalar>
      36       18926 : Database &BasisAtom<Scalar>::get_database() const {
      37       18926 :     return database;
      38             : }
      39             : 
      40             : template <typename Scalar>
      41       14948 : const std::string &BasisAtom<Scalar>::get_species() const {
      42       14948 :     return this->kets[0]->get_species();
      43             : }
      44             : 
      45             : template <typename Scalar>
      46     1752434 : int BasisAtom<Scalar>::get_ket_index_from_id(size_t ket_id) const {
      47     1752434 :     if (!ket_id_to_ket_index.contains(ket_id)) {
      48           0 :         return -1;
      49             :     }
      50     1752434 :     return ket_id_to_ket_index.at(ket_id);
      51             : }
      52             : 
      53             : template <typename Scalar>
      54       11805 : const std::string &BasisAtom<Scalar>::get_canonical_basis_id() const {
      55       11805 :     return canonical_basis_id;
      56             : }
      57             : 
      58             : template <typename Scalar>
      59             : std::shared_ptr<const typename BasisAtom<Scalar>::Type>
      60          48 : BasisAtom<Scalar>::merge(std::shared_ptr<const Type> other) const {
      61          48 :     if (&database != &other->database) {
      62           0 :         throw std::invalid_argument("Cannot merge atomic bases from different Database instances.");
      63             :     }
      64          48 :     if (get_species() != other->get_species()) {
      65           0 :         throw std::invalid_argument("Cannot merge atomic bases with different species.");
      66             :     }
      67          48 :     if (!this->is_canonical() || !other->is_canonical()) {
      68           0 :         throw std::invalid_argument(
      69             :             "Cannot merge non-canonical bases (i.e., bases with non-identity coefficients). "
      70             :             "Canonicalize the bases first.");
      71             :     }
      72             : 
      73          48 :     BasisAtomCreator<Scalar> creator(database);
      74          48 :     std::unordered_set<size_t> ket_ids;
      75          48 :     ket_ids.reserve(this->kets.size() + other->kets.size());
      76         144 :     for (const auto &basis : {this, other.get()}) {
      77        1680 :         for (const auto &ket : basis->kets) {
      78        1584 :             if (ket_ids.insert(ket->get_id_in_database()).second) {
      79         910 :                 creator.add_ket(ket);
      80             :             }
      81             :         }
      82             :     }
      83          96 :     return creator.create();
      84          48 : }
      85             : 
      86             : template <typename Scalar>
      87             : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
      88         999 : BasisAtom<Scalar>::get_matrix_elements(std::shared_ptr<const Type> other, OperatorType type,
      89             :                                        int q) const {
      90         999 :     auto matrix_elements = this->get_database().get_matrix_elements_in_canonical_basis(
      91             :         this->shared_from_this(), other, type, q);
      92         999 :     matrix_elements =
      93         999 :         other->get_coefficients().adjoint() * matrix_elements * this->get_coefficients();
      94             : 
      95         999 :     assert(static_cast<size_t>(matrix_elements.rows()) == other->get_number_of_states());
      96         999 :     assert(static_cast<size_t>(matrix_elements.cols()) == this->get_number_of_states());
      97             : 
      98         999 :     return matrix_elements;
      99           0 : }
     100             : 
     101             : // Explicit instantiations
     102             : template class BasisAtom<double>;
     103             : template class BasisAtom<std::complex<double>>;
     104             : } // namespace pairinteraction