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