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 <cmath>
      21             : #include <doctest/doctest.h>
      22             : #include <fmt/ranges.h>
      23             : 
      24             : namespace pairinteraction {
      25             : 
      26             : constexpr double VOLT_PER_CM_IN_ATOMIC_UNITS = 1 / 5.14220675112e9;
      27             : constexpr double UM_IN_ATOMIC_UNITS = 1 / 5.29177210544e-5;
      28             : constexpr double HARTREE_IN_GHZ = 6579683.920501762;
      29             : 
      30           1 : DOCTEST_TEST_CASE("construct a pair Hamiltonian") {
      31           1 :     auto &database = Database::get_global_instance();
      32           1 :     auto diagonalizer = DiagonalizerEigen<double>();
      33             : 
      34             :     // Construct and diagonalize the constituent systems
      35           1 :     auto basis = BasisAtomCreator<double>()
      36           2 :                      .set_species("Rb")
      37           2 :                      .restrict_quantum_number("n", 60, 61)
      38           2 :                      .restrict_quantum_number("l", 0, 1)
      39           2 :                      .restrict_quantum_number("m", -0.5, 0.5)
      40           1 :                      .create(database);
      41             : 
      42           1 :     SystemAtom<double> system1(basis);
      43           1 :     SystemAtom<double> system2(basis);
      44           1 :     system1.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
      45           1 :     system2.set_electric_field({0, 0, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS});
      46           1 :     diagonalize<SystemAtom<double>>({system1, system2}, diagonalizer);
      47             : 
      48             :     // Construct and diagonalize the system_pair
      49           1 :     auto basis_pair = BasisPairCreator<double>().add(system1).add(system2).create();
      50           1 :     DOCTEST_MESSAGE("Number of states in pair basis: ", basis_pair->get_number_of_states());
      51             : 
      52           1 :     auto system_pair = SystemPair<double>(basis_pair);
      53           1 :     system_pair.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
      54           1 :     system_pair.diagonalize(diagonalizer);
      55             : 
      56             :     // Print the largest and smallest eigenenergies
      57           1 :     auto eigenenergies = system_pair.get_eigenenergies();
      58           1 :     DOCTEST_MESSAGE("Lowest energy: ", eigenenergies.minCoeff());
      59           1 :     DOCTEST_MESSAGE("Highest energy: ", eigenenergies.maxCoeff());
      60           1 : }
      61             : 
      62           1 : DOCTEST_TEST_CASE("construct a pair Hamiltonian in a non-canonical pair basis") {
      63           1 :     auto &database = Database::get_global_instance();
      64             : 
      65           1 :     auto basis = BasisAtomCreator<double>()
      66           2 :                      .set_species("Rb")
      67           2 :                      .restrict_quantum_number("n", 60, 61)
      68           2 :                      .restrict_quantum_number("l", 0, 1)
      69           2 :                      .restrict_quantum_number("m", -0.5, 0.5)
      70           1 :                      .create(database);
      71             : 
      72           1 :     auto diagonalizer = DiagonalizerEigen<double>();
      73           1 :     SystemAtom<double> system1(basis);
      74           1 :     SystemAtom<double> system2(basis);
      75           1 :     system1.set_electric_field({0, 0, 1 * VOLT_PER_CM_IN_ATOMIC_UNITS});
      76           1 :     system2.set_electric_field({0, 0, 2 * VOLT_PER_CM_IN_ATOMIC_UNITS});
      77           1 :     diagonalize<SystemAtom<double>>({system1, system2}, diagonalizer);
      78             : 
      79           1 :     auto pair_basis = BasisPairCreator<double>().add(system1).add(system2).create();
      80           1 :     DOCTEST_REQUIRE(pair_basis->get_number_of_states() >= 2);
      81             : 
      82           1 :     SystemPair<double> reference_system(pair_basis);
      83           1 :     reference_system.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
      84           1 :     const auto &reference_matrix = reference_system.get_matrix();
      85             : 
      86             :     Eigen::SparseMatrix<double, Eigen::RowMajor> transformation(
      87           1 :         static_cast<Eigen::Index>(pair_basis->get_number_of_states()),
      88           2 :         static_cast<Eigen::Index>(pair_basis->get_number_of_states()));
      89           1 :     transformation.setIdentity();
      90             : 
      91           1 :     double inverse_sqrt_two = 1 / std::sqrt(2.0);
      92           1 :     transformation.coeffRef(0, 0) = inverse_sqrt_two;
      93           1 :     transformation.coeffRef(1, 0) = inverse_sqrt_two;
      94           1 :     transformation.coeffRef(0, 1) = inverse_sqrt_two;
      95           1 :     transformation.coeffRef(1, 1) = -inverse_sqrt_two;
      96           1 :     transformation.makeCompressed();
      97             : 
      98           1 :     auto transformed_pair_basis = pair_basis->transformed(transformation);
      99           1 :     SystemPair<double> transformed_system(transformed_pair_basis);
     100           1 :     transformed_system.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
     101             : 
     102             :     Eigen::SparseMatrix<double, Eigen::RowMajor> expected_matrix =
     103           1 :         transformation.adjoint() * reference_matrix * transformation;
     104             : 
     105           1 :     DOCTEST_CHECK(transformed_system.get_matrix().isApprox(expected_matrix, 1e-11));
     106           1 : }
     107             : 
     108             : #ifdef WITH_LAPACKE
     109           1 : DOCTEST_TEST_CASE("diagonalize with lapacke_evr") {
     110           1 :     auto &database = Database::get_global_instance();
     111           1 :     auto diagonalizer = DiagonalizerLapackeEvr<double>();
     112             : 
     113             :     // Construct the state of interest
     114           1 :     auto ket = KetAtomCreator()
     115           2 :                    .set_species("Rb")
     116           2 :                    .set_quantum_number("n", 60)
     117           2 :                    .set_quantum_number("l", 0)
     118           2 :                    .set_quantum_number("j", 0.5)
     119           2 :                    .set_quantum_number("m", 0.5)
     120           1 :                    .create(database);
     121             : 
     122             :     // Construct and diagonalize the constituent systems
     123           1 :     auto basis = BasisAtomCreator<double>()
     124           2 :                      .set_species("Rb")
     125           2 :                      .restrict_quantum_number("n", ket->get_quantum_number("n") - 3,
     126           2 :                                               ket->get_quantum_number("n") + 3)
     127           2 :                      .restrict_quantum_number("l", ket->get_quantum_number("l") - 1,
     128           2 :                                               ket->get_quantum_number("l") + 1)
     129           1 :                      .create(database);
     130           1 :     SystemAtom<double> system(basis);
     131             : 
     132             :     // Construct and diagonalize the system_pair
     133           1 :     auto basis_pair = BasisPairCreator<double>()
     134           1 :                           .add(system)
     135           1 :                           .add(system)
     136           1 :                           .restrict_energy(2 * ket->get_energy() - 2 / HARTREE_IN_GHZ,
     137           1 :                                            2 * ket->get_energy() + 2 / HARTREE_IN_GHZ)
     138           1 :                           .create();
     139           1 :     DOCTEST_MESSAGE("Number of states in pair basis: ", basis_pair->get_number_of_states());
     140             : 
     141           1 :     auto system_pair = SystemPair<double>(basis_pair);
     142           1 :     system_pair.set_distance_vector({0, 0, 3 * UM_IN_ATOMIC_UNITS});
     143           1 :     system_pair.diagonalize(diagonalizer, 2 * ket->get_energy() - 0.5 / HARTREE_IN_GHZ,
     144           1 :                             2 * ket->get_energy() + 0.5 / HARTREE_IN_GHZ);
     145             : 
     146             :     // Print the largest and smallest eigenenergies
     147           1 :     auto eigenenergies = system_pair.get_eigenenergies();
     148           1 :     DOCTEST_MESSAGE("Lowest energy: ", eigenenergies.minCoeff());
     149           1 :     DOCTEST_MESSAGE("Highest energy: ", eigenenergies.maxCoeff());
     150           1 : }
     151             : #endif
     152             : } // namespace pairinteraction