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 <stdexcept>
      12             : #include <unordered_set>
      13             : 
      14             : namespace pairinteraction {
      15             : template <typename Scalar>
      16        1541 : BasisAtom<Scalar>::BasisAtom(Private /*unused*/, ketvec_t &&kets, std::string &&canonical_basis_id,
      17             :                              Database &database)
      18        1541 :     : Basis<BasisAtom<Scalar>>(std::move(kets)), canonical_basis_id(std::move(canonical_basis_id)),
      19        3082 :       database(database) {
      20       44334 :     for (size_t i = 0; i < this->kets.size(); ++i) {
      21       42793 :         ket_id_to_ket_index[this->kets[i]->get_id_in_database()] = i;
      22             :     }
      23        1541 : }
      24             : 
      25             : template <typename Scalar>
      26       15595 : Database &BasisAtom<Scalar>::get_database() const {
      27       15595 :     return database;
      28             : }
      29             : 
      30             : template <typename Scalar>
      31       12118 : const std::string &BasisAtom<Scalar>::get_species() const {
      32       12118 :     return this->kets[0]->get_species();
      33             : }
      34             : 
      35             : template <typename Scalar>
      36     1561128 : int BasisAtom<Scalar>::get_ket_index_from_id(size_t ket_id) const {
      37     1561128 :     if (!ket_id_to_ket_index.contains(ket_id)) {
      38           0 :         return -1;
      39             :     }
      40     1561128 :     return ket_id_to_ket_index.at(ket_id);
      41             : }
      42             : 
      43             : template <typename Scalar>
      44       10424 : const std::string &BasisAtom<Scalar>::get_canonical_basis_id() const {
      45       10424 :     return canonical_basis_id;
      46             : }
      47             : 
      48             : template <typename Scalar>
      49             : std::shared_ptr<const typename BasisAtom<Scalar>::Type>
      50          48 : BasisAtom<Scalar>::merge(std::shared_ptr<const Type> other) const {
      51          48 :     if (&database != &other->database) {
      52           0 :         throw std::invalid_argument("Cannot merge atomic bases from different Database instances.");
      53             :     }
      54          48 :     if (get_species() != other->get_species()) {
      55           0 :         throw std::invalid_argument("Cannot merge atomic bases with different species.");
      56             :     }
      57          48 :     if (!this->is_canonical() || !other->is_canonical()) {
      58           0 :         throw std::invalid_argument(
      59             :             "Cannot merge non-canonical bases (i.e., bases with non-identity coefficients). "
      60             :             "Canonicalize the bases first.");
      61             :     }
      62             : 
      63          48 :     BasisAtomCreator<Scalar> creator;
      64          48 :     std::unordered_set<size_t> ket_ids;
      65          48 :     ket_ids.reserve(this->kets.size() + other->kets.size());
      66         144 :     for (const auto &basis : {this, other.get()}) {
      67        1680 :         for (const auto &ket : basis->kets) {
      68        1584 :             if (ket_ids.insert(ket->get_id_in_database()).second) {
      69         910 :                 creator.add_ket(ket);
      70             :             }
      71             :         }
      72             :     }
      73          96 :     return creator.create(database);
      74          48 : }
      75             : 
      76             : template <typename Scalar>
      77             : Eigen::SparseMatrix<Scalar, Eigen::RowMajor>
      78         987 : BasisAtom<Scalar>::get_matrix_elements(std::shared_ptr<const Type> other, OperatorType type,
      79             :                                        int q) const {
      80         987 :     auto matrix_elements = this->get_database().get_matrix_elements_in_canonical_basis(
      81             :         this->shared_from_this(), other, type, q);
      82         987 :     matrix_elements =
      83         987 :         other->get_coefficients().adjoint() * matrix_elements * this->get_coefficients();
      84             : 
      85         987 :     assert(static_cast<size_t>(matrix_elements.rows()) == other->get_number_of_states());
      86         987 :     assert(static_cast<size_t>(matrix_elements.cols()) == this->get_number_of_states());
      87             : 
      88         987 :     return matrix_elements;
      89           0 : }
      90             : 
      91             : // Explicit instantiations
      92             : template class BasisAtom<double>;
      93             : template class BasisAtom<std::complex<double>>;
      94             : } // namespace pairinteraction