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         177 : BasisPairCreator<Scalar>::BasisPairCreator() : product_of_parities(Parity::UNKNOWN) {}
      22             : 
      23             : template <typename Scalar>
      24         354 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::add(const SystemAtom<Scalar> &system_atom) {
      25         354 :     if (!system_atom.is_diagonal()) {
      26           0 :         throw std::invalid_argument("The system must be diagonalized before it can be added.");
      27             :     }
      28         354 :     systems_atom.push_back(system_atom);
      29         354 :     return *this;
      30             : }
      31             : 
      32             : template <typename Scalar>
      33         152 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_energy(real_t min, real_t max) {
      34         152 :     range_energy = {min, max};
      35         152 :     return *this;
      36             : }
      37             : 
      38             : template <typename Scalar>
      39          23 : BasisPairCreator<Scalar> &BasisPairCreator<Scalar>::restrict_quantum_number_m(real_t min,
      40             :                                                                               real_t max) {
      41          23 :     range_quantum_number_m = {min, max};
      42          23 :     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         177 : std::shared_ptr<const BasisPair<Scalar>> BasisPairCreator<Scalar>::create() const {
      53         177 :     if (systems_atom.size() != 2) {
      54           0 :         throw std::invalid_argument("Two SystemAtom must be added before creating the BasisPair.");
      55             :     }
      56             : 
      57         177 :     constexpr real_t numerical_precision = 100 * std::numeric_limits<real_t>::epsilon();
      58             : 
      59             :     // Sort the states, which are eigenstates, by their energy
      60         177 :     auto system1 = systems_atom[0].get();
      61         177 :     auto system2 = systems_atom[1].get();
      62         177 :     system1.transform(system1.get_sorter({TransformationType::SORT_BY_ENERGY}));
      63         177 :     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         177 :     auto basis1 = system1.get_basis();
      68         177 :     auto basis2 = system2.get_basis();
      69         177 :     auto eigenenergies1 = system1.get_eigenenergies();
      70         177 :     auto eigenenergies2 = system2.get_eigenenergies();
      71         177 :     real_t *eigenenergies2_begin = eigenenergies2.data();
      72         177 :     real_t *eigenenergies2_end = eigenenergies2_begin + eigenenergies2.size();
      73             : 
      74         177 :     ketvec_t kets;
      75         177 :     kets.reserve(eigenenergies1.size() * eigenenergies2.size());
      76             : 
      77         177 :     typename basis_t::map_range_t map_range_of_state_index2;
      78         177 :     map_range_of_state_index2.reserve(eigenenergies1.size());
      79             : 
      80         177 :     typename basis_t::map_indices_t state_indices_to_ket_index;
      81             : 
      82             :     // Loop only over states with an allowed energy
      83         177 :     size_t ket_index = 0;
      84       40113 :     for (size_t idx1 = 0; idx1 < static_cast<size_t>(eigenenergies1.size()); ++idx1) {
      85             :         // Get the energetically allowed range of the second index
      86       39936 :         size_t min = 0;
      87       39936 :         size_t max = eigenenergies2.size();
      88       39936 :         if (range_energy.is_finite()) {
      89       38627 :             real_t min_val2 = range_energy.min() - eigenenergies1[idx1];
      90       38627 :             real_t max_val2 = range_energy.max() - eigenenergies1[idx1];
      91       38627 :             min =
      92       38627 :                 std::distance(eigenenergies2_begin,
      93             :                               std::lower_bound(eigenenergies2_begin, eigenenergies2_end, min_val2));
      94       38627 :             max =
      95       38627 :                 std::distance(eigenenergies2_begin,
      96             :                               std::upper_bound(eigenenergies2_begin, eigenenergies2_end, max_val2));
      97             :         }
      98       39936 :         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      409928 :         for (size_t idx2 = min; idx2 < max; ++idx2) {
     102             :             // Get energy
     103      184996 :             const real_t energy = eigenenergies1[idx1] + eigenenergies2[idx2];
     104      184996 :             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      184996 :             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        9033 :                     continue;
     113             :                 }
     114             :             }
     115             : 
     116             :             // Create a KetPair object
     117      554988 :             auto ket = std::make_shared<ket_t>(
     118      184996 :                 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      184996 :             if (ket->has_quantum_number_m()) {
     124      131928 :                 if (range_quantum_number_m.is_finite() &&
     125       10932 :                     (ket->get_quantum_number_m() <
     126       10932 :                          range_quantum_number_m.min() - numerical_precision ||
     127        4377 :                      ket->get_quantum_number_m() >
     128        4377 :                          range_quantum_number_m.max() + numerical_precision)) {
     129        9033 :                     continue;
     130             :                 }
     131       64000 :             } 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      175963 :             kets.emplace_back(std::move(ket));
     138             : 
     139             :             // Store the ket index
     140      175963 :             state_indices_to_ket_index.emplace(std::vector<size_t>{idx1, idx2}, ket_index++);
     141             :         }
     142             :     }
     143             : 
     144         177 :     kets.shrink_to_fit();
     145             : 
     146         177 :     return std::make_shared<basis_t>(typename basis_t::Private(), std::move(kets),
     147         177 :                                      std::move(map_range_of_state_index2),
     148         531 :                                      std::move(state_indices_to_ket_index), basis1, basis2);
     149      185173 : }
     150             : 
     151             : // Explicit instantiations
     152             : template class BasisPairCreator<double>;
     153             : template class BasisPairCreator<std::complex<double>>;
     154             : } // namespace pairinteraction