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