Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/system/SystemPair.hpp"
       5             : 
       6             : #include "pairinteraction/basis/BasisAtom.hpp"
       7             : #include "pairinteraction/basis/BasisAtomCreator.hpp"
       8             : #include "pairinteraction/basis/BasisPair.hpp"
       9             : #include "pairinteraction/basis/BasisPairCreator.hpp"
      10             : #include "pairinteraction/database/Database.hpp"
      11             : #include "pairinteraction/diagonalize/DiagonalizerEigen.hpp"
      12             : #include "pairinteraction/diagonalize/DiagonalizerFeast.hpp"
      13             : #include "pairinteraction/diagonalize/DiagonalizerLapackeEvr.hpp"
      14             : #include "pairinteraction/diagonalize/diagonalize.hpp"
      15             : #include "pairinteraction/ket/KetAtom.hpp"
      16             : #include "pairinteraction/ket/KetAtomCreator.hpp"
      17             : #include "pairinteraction/system/SystemAtom.hpp"
      18             : #include "pairinteraction/utils/Range.hpp"
      19             : 
      20             : #include <doctest/doctest.h>
      21             : #include <fmt/ranges.h>
      22             : 
      23             : namespace pairinteraction {
      24             : 
      25             : constexpr double VOLT_PER_CM_IN_ATOMIC_UNITS = 1 / 5.14220675112e9;
      26             : constexpr double UM_IN_ATOMIC_UNITS = 1 / 5.29177210544e-5;
      27             : constexpr double HARTREE_IN_GHZ = 6579683.920501762;
      28             : 
      29           1 : DOCTEST_TEST_CASE("construct a pair Hamiltonian") {
      30           1 :     auto &database = Database::get_global_instance();
      31           1 :     auto diagonalizer = DiagonalizerEigen<double>();
      32             : 
      33             :     // Construct and diagonalize the constituent systems
      34           1 :     auto basis = BasisAtomCreator<double>()
      35           2 :                      .set_species("Rb")
      36           1 :                      .restrict_quantum_number_n(60, 61)
      37           1 :                      .restrict_quantum_number_l(0, 1)
      38           1 :                      .restrict_quantum_number_m(-0.5, 0.5)
      39           1 :                      .create(database);
      40             : 
      41           1 :     SystemAtom<double> system1(basis);
      42           1 :     SystemAtom<double> system2(basis);
      43           1 :     system1.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
      44           1 :     system2.set_electric_field({0, 0, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS});
      45           1 :     diagonalize<SystemAtom<double>>({system1, system2}, diagonalizer);
      46             : 
      47             :     // Construct and diagonalize the system_pair
      48           1 :     auto basis_pair = BasisPairCreator<double>().add(system1).add(system2).create();
      49           1 :     DOCTEST_MESSAGE("Number of states in pair basis: ", basis_pair->get_number_of_states());
      50             : 
      51           1 :     auto system_pair = SystemPair<double>(basis_pair);
      52           1 :     system_pair.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
      53           1 :     system_pair.diagonalize(diagonalizer);
      54             : 
      55             :     // Print the largest and smallest eigenenergies
      56           1 :     auto eigenenergies = system_pair.get_eigenenergies();
      57           1 :     DOCTEST_MESSAGE("Lowest energy: ", eigenenergies.minCoeff());
      58           1 :     DOCTEST_MESSAGE("Highest energy: ", eigenenergies.maxCoeff());
      59           1 : }
      60             : 
      61             : #ifdef WITH_LAPACKE
      62           1 : DOCTEST_TEST_CASE("diagonalize with lapacke_evr") {
      63           1 :     auto &database = Database::get_global_instance();
      64           1 :     auto diagonalizer = DiagonalizerLapackeEvr<double>();
      65             : 
      66             :     // Construct the state of interest
      67           1 :     auto ket = KetAtomCreator()
      68           2 :                    .set_species("Rb")
      69           1 :                    .set_quantum_number_n(60)
      70           1 :                    .set_quantum_number_l(0)
      71           1 :                    .set_quantum_number_j(0.5)
      72           1 :                    .set_quantum_number_m(0.5)
      73           1 :                    .create(database);
      74             : 
      75             :     // Construct and diagonalize the constituent systems
      76           1 :     auto basis = BasisAtomCreator<double>()
      77           2 :                      .set_species("Rb")
      78           1 :                      .restrict_quantum_number_n(ket->get_quantum_number_n() - 3,
      79           1 :                                                 ket->get_quantum_number_n() + 3)
      80           1 :                      .restrict_quantum_number_l(ket->get_quantum_number_l() - 1,
      81           1 :                                                 ket->get_quantum_number_l() + 1)
      82           1 :                      .create(database);
      83           1 :     SystemAtom<double> system(basis);
      84             : 
      85             :     // Construct and diagonalize the system_pair
      86           2 :     auto basis_pair = BasisPairCreator<double>()
      87           1 :                           .add(system)
      88           1 :                           .add(system)
      89           1 :                           .restrict_energy(2 * ket->get_energy() - 2 / HARTREE_IN_GHZ,
      90           1 :                                            2 * ket->get_energy() + 2 / HARTREE_IN_GHZ)
      91           1 :                           .create();
      92           1 :     DOCTEST_MESSAGE("Number of states in pair basis: ", basis_pair->get_number_of_states());
      93             : 
      94           1 :     auto system_pair = SystemPair<double>(basis_pair);
      95           1 :     system_pair.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
      96           1 :     system_pair.diagonalize(diagonalizer, 2 * ket->get_energy() - 0.5 / HARTREE_IN_GHZ,
      97           1 :                             2 * ket->get_energy() + 0.5 / HARTREE_IN_GHZ);
      98             : 
      99             :     // Print the largest and smallest eigenenergies
     100           1 :     auto eigenenergies = system_pair.get_eigenenergies();
     101           1 :     DOCTEST_MESSAGE("Lowest energy: ", eigenenergies.minCoeff());
     102           1 :     DOCTEST_MESSAGE("Highest energy: ", eigenenergies.maxCoeff());
     103           1 : }
     104             : #endif
     105             : } // namespace pairinteraction