Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 PairInteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/basis/BasisPairCreator.hpp"
       5             : 
       6             : #include "pairinteraction/basis/BasisAtom.hpp"
       7             : #include "pairinteraction/basis/BasisPair.hpp"
       8             : #include "pairinteraction/enums/Parity.hpp"
       9             : #include "pairinteraction/enums/TransformationType.hpp"
      10             : #include "pairinteraction/ket/KetAtom.hpp"
      11             : #include "pairinteraction/ket/KetPair.hpp"
      12             : #include "pairinteraction/system/SystemAtom.hpp"
      13             : 
      14             : #include <algorithm>
      15             : #include <cassert>
      16             : #include <limits>
      17             : #include <memory>
      18             : 
      19             : namespace pairinteraction {
      20             : template <typename Scalar>
      21         103 : BasisPairCreator<Scalar>::BasisPairCreator() : product_of_parities(Parity::UNKNOWN) {}
      22             : 
      23             : template <typename Scalar>
      24         206 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::add(const SystemAtom<Scalar> &system_atom) {
      25         206 :     if (!system_atom.is_diagonal()) {
      26           0 :         throw std::invalid_argument("The system must be diagonalized before it can be added.");
      27             :     }
      28         206 :     systems_atom.push_back(system_atom);
      29         206 :     return *this;
      30             : }
      31             : 
      32             : template <typename Scalar>
      33          94 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_energy(real_t min, real_t max) {
      34          94 :     range_energy = {min, max};
      35          94 :     return *this;
      36             : }
      37             : 
      38             : template <typename Scalar>
      39          16 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_quantum_number_m(real_t min,
      40             :                                                                               real_t max) {
      41          16 :     range_quantum_number_m = {min, max};
      42          16 :     return *this;
      43             : }
      44             : 
      45             : template <typename Scalar>
      46           0 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_product_of_parities(Parity value) {
      47           0 :     product_of_parities = value;
      48           0 :     return *this;
      49             : }
      50             : 
      51             : template <typename Scalar>
      52         103 : std::shared_ptr<const BasisPair<Scalar>> BasisPairCreator<Scalar>::create() const {
      53         103 :     if (systems_atom.size() != 2) {
      54           0 :         throw std::invalid_argument("Two SystemAtom must be added before creating the BasisPair.");
      55             :     }
      56             : 
      57         103 :     constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
      58             : 
      59             :     // Sort the states, which are eigenstates, by their energy
      60         103 :     auto system1 = systems_atom[0].get();
      61         103 :     auto system2 = systems_atom[1].get();
      62         103 :     system1.transform(system1.get_sorter({TransformationType::SORT_BY_ENERGY}));
      63         103 :     system2.transform(system2.get_sorter({TransformationType::SORT_BY_ENERGY}));
      64             : 
      65             :     // Construct the canonical basis that contains all KetPair objects with allowed energies and
      66             :     // quantum numbers
      67         103 :     auto basis1 = system1.get_basis();
      68         103 :     auto basis2 = system2.get_basis();
      69         103 :     auto eigenenergies1 = system1.get_eigenenergies();
      70         103 :     auto eigenenergies2 = system2.get_eigenenergies();
      71         103 :     real_t *eigenenergies2_begin = eigenenergies2.data();
      72         103 :     real_t *eigenenergies2_end = eigenenergies2_begin + eigenenergies2.size();
      73             : 
      74         103 :     ketvec_t kets;
      75         103 :     kets.reserve(eigenenergies1.size() * eigenenergies2.size());
      76             : 
      77         103 :     typename basis_t::map_range_t map_range_of_state_index2;
      78         103 :     map_range_of_state_index2.reserve(eigenenergies1.size());
      79             : 
      80         103 :     typename basis_t::map_indices_t state_indices_to_ket_index;
      81             : 
      82             :     // Loop only over states with an allowed energy
      83         103 :     size_t ket_index = 0;
      84       23760 :     for (size_t idx1 = 0; idx1 < static_cast<size_t>(eigenenergies1.size()); ++idx1) {
      85             :         // Get the energetically allowed range of the second index
      86       23657 :         size_t min = 0;
      87       23657 :         size_t max = eigenenergies2.size();
      88       23657 :         if (range_energy.is_finite()) {
      89       23351 :             real_t min_val2 = range_energy.min() - eigenenergies1[idx1];
      90       23351 :             real_t max_val2 = range_energy.max() - eigenenergies1[idx1];
      91       23351 :             min =
      92       23351 :                 std::distance(eigenenergies2_begin,
      93             :                               std::lower_bound(eigenenergies2_begin, eigenenergies2_end, min_val2));
      94       23351 :             max =
      95       23351 :                 std::distance(eigenenergies2_begin,
      96             :                               std::upper_bound(eigenenergies2_begin, eigenenergies2_end, max_val2));
      97             :         }
      98       23657 :         map_range_of_state_index2.emplace(idx1, typename basis_t::range_t(min, max));
      99             : 
     100             :         // Loop over the energetically allowed range of the second index
     101      172591 :         for (size_t idx2 = min; idx2 < max; ++idx2) {
     102             :             // Get energy
     103       74467 :             const real_t energy = eigenenergies1[idx1] + eigenenergies2[idx2];
     104       74467 :             assert(!range_energy.is_finite() ||
     105             :                    (energy >= range_energy.min() && energy <= range_energy.max()));
     106             : 
     107             :             // Check the parity of the sum of the parities
     108       74467 :             if (product_of_parities != Parity::UNKNOWN) {
     109           0 :                 if (static_cast<int>(basis1->get_parity(idx1)) *
     110           0 :                         static_cast<int>(basis2->get_parity(idx2)) !=
     111           0 :                     static_cast<int>(product_of_parities)) {
     112        4820 :                     continue;
     113             :                 }
     114             :             }
     115             : 
     116             :             // Create a KetPair object
     117      223401 :             auto ket = std::make_shared<ket_t>(
     118       74467 :                 typename ket_t::Private(), std::initializer_list<size_t>{idx1, idx2},
     119             :                 std::initializer_list<std::shared_ptr<const BasisAtom<Scalar>>>{basis1, basis2},
     120             :                 energy);
     121             : 
     122             :             // Check the quantum number m
     123       74467 :             if (ket->has_quantum_number_m()) {
     124       80323 :                 if (range_quantum_number_m.is_finite() &&
     125        5856 :                     (ket->get_quantum_number_m() <
     126        5856 :                          range_quantum_number_m.min() - numerical_precision ||
     127        2332 :                      ket->get_quantum_number_m() >
     128        2332 :                          range_quantum_number_m.max() + numerical_precision)) {
     129        4820 :                     continue;
     130             :                 }
     131           0 :             } else if (range_quantum_number_m.is_finite()) {
     132           0 :                 throw std::invalid_argument(
     133             :                     "The quantum number m must not be restricted because it is not well-defined.");
     134             :             }
     135             : 
     136             :             // Store the KetPair object as a ket
     137       69647 :             kets.emplace_back(std::move(ket));
     138             : 
     139             :             // Store the ket index
     140       69647 :             state_indices_to_ket_index.emplace(std::vector<size_t>{idx1, idx2}, ket_index++);
     141             :         }
     142             :     }
     143             : 
     144         103 :     kets.shrink_to_fit();
     145             : 
     146         103 :     return std::make_shared<basis_t>(typename basis_t::Private(), std::move(kets),
     147         103 :                                      std::move(map_range_of_state_index2),
     148         309 :                                      std::move(state_indices_to_ket_index), basis1, basis2);
     149       74570 : }
     150             : 
     151             : // Explicit instantiations
     152             : template class BasisPairCreator<double>;
     153             : template class BasisPairCreator<std::complex<double>>;
     154             : } // namespace pairinteraction