Line data    Source code

       1             : // SPDX-FileCopyrightText: 2024 Pairinteraction Developers
       2             : // SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : #include "pairinteraction/pairinteraction.hpp"
       5             : #include "pairinteraction/utils/args.hpp"
       6             : 
       7             : #include <filesystem>
       8             : #include <spdlog/spdlog.h>
       9             : 
      10           1 : int main(int argc, char **argv) {
      11             :     // Call the setup function to configure logging
      12           1 :     pairinteraction::setup();
      13             : 
      14             :     // Create a database instance
      15           1 :     std::filesystem::path database_dir;
      16           1 :     bool download_missing = false;
      17             : 
      18           4 :     for (int i = 1; i < argc; ++i) {
      19           3 :         bool found = pairinteraction::args::parse_download_missing(i, argc, argv, download_missing);
      20           3 :         if (!found) {
      21           3 :             pairinteraction::args::parse_database_dir(i, argc, argv, database_dir);
      22             :         }
      23             :     }
      24             : 
      25           1 :     pairinteraction::Database database(download_missing, true, database_dir);
      26             : 
      27             :     // Test the calculation of a dipole operator
      28           1 :     bool success = true;
      29             : 
      30             :     // Create a dipole operator coupling two specific states
      31           1 :     auto ket1 = pairinteraction::KetAtomCreator()
      32           2 :                     .set_species("Rb")
      33           1 :                     .set_quantum_number_n(60)
      34           1 :                     .set_quantum_number_l(0)
      35           1 :                     .set_quantum_number_j(0.5)
      36           1 :                     .set_quantum_number_m(0.5)
      37           1 :                     .create(database);
      38             : 
      39           1 :     auto ket2 = pairinteraction::KetAtomCreator()
      40           2 :                     .set_species("Rb")
      41           1 :                     .set_quantum_number_n(60)
      42           1 :                     .set_quantum_number_l(1)
      43           1 :                     .set_quantum_number_j(0.5)
      44           1 :                     .set_quantum_number_m(0.5)
      45           1 :                     .create(database);
      46             : 
      47             :     auto basis_ket1_ket2 =
      48           2 :         pairinteraction::BasisAtomCreator<double>().append_ket(ket1).append_ket(ket2).create(
      49           1 :             database);
      50             : 
      51             :     pairinteraction::OperatorAtom<double> dipole_ket1_ket2(
      52           1 :         basis_ket1_ket2, pairinteraction::OperatorType::ELECTRIC_DIPOLE, 0);
      53           1 :     double dipole_ket1_ket2_value = dipole_ket1_ket2.get_matrix().coeff(0, 1);
      54             : 
      55           1 :     double reference = 1247.5985544327225;
      56             : 
      57           1 :     if (std::abs(dipole_ket1_ket2_value - reference) >
      58           1 :         10 * std::numeric_limits<double>::epsilon()) {
      59           0 :         SPDLOG_ERROR("The dipole operator value is not correct. Value: {}", dipole_ket1_ket2_value);
      60           0 :         success = false;
      61             :     }
      62             : 
      63           1 :     dipole_ket1_ket2 = 2 * dipole_ket1_ket2;
      64           1 :     dipole_ket1_ket2_value = dipole_ket1_ket2.get_matrix().coeff(0, 1);
      65             : 
      66           1 :     if (std::abs(dipole_ket1_ket2_value - 2 * reference) >
      67           1 :         10 * std::numeric_limits<double>::epsilon()) {
      68           0 :         SPDLOG_ERROR("The dipole operator value is not correct after multiplication by a scalar.");
      69           0 :         success = false;
      70             :     }
      71             : 
      72           1 :     dipole_ket1_ket2 = dipole_ket1_ket2 + dipole_ket1_ket2;
      73           1 :     dipole_ket1_ket2_value = dipole_ket1_ket2.get_matrix().coeff(0, 1);
      74             : 
      75           1 :     if (std::abs(dipole_ket1_ket2_value - 4 * reference) >
      76           1 :         10 * std::numeric_limits<double>::epsilon()) {
      77           0 :         SPDLOG_ERROR("The dipole operator value is not correct after summation.");
      78           0 :         success = false;
      79             :     }
      80             : 
      81             :     // Create dipole operators in a typical basis
      82           1 :     auto basis = pairinteraction::BasisAtomCreator<std::complex<double>>()
      83           2 :                      .set_species("Sr88_singlet")
      84           1 :                      .restrict_quantum_number_n(60, 63)
      85           1 :                      .restrict_quantum_number_l(0, 3)
      86           1 :                      .create(database);
      87             : 
      88             :     pairinteraction::OperatorAtom<std::complex<double>> dipole_0(
      89           1 :         basis, pairinteraction::OperatorType::ELECTRIC_DIPOLE, 0);
      90             :     pairinteraction::OperatorAtom<std::complex<double>> dipole_p(
      91           1 :         basis, pairinteraction::OperatorType::ELECTRIC_DIPOLE, 1);
      92             :     pairinteraction::OperatorAtom<std::complex<double>> dipole_m(
      93           1 :         basis, pairinteraction::OperatorType::ELECTRIC_DIPOLE, -1);
      94             : 
      95           1 :     if (dipole_0.get_matrix().rows() != 64) {
      96           0 :         SPDLOG_ERROR("Wrong dimension.");
      97           0 :         success = false;
      98             :     }
      99             : 
     100           1 :     if (dipole_p.get_matrix().rows() != 64) {
     101           0 :         SPDLOG_ERROR("Wrong dimension.");
     102           0 :         success = false;
     103             :     }
     104             : 
     105           1 :     if (dipole_m.get_matrix().rows() != 64) {
     106           0 :         SPDLOG_ERROR("Wrong dimension.");
     107           0 :         success = false;
     108             :     }
     109             : 
     110           1 :     if (dipole_0.get_matrix().nonZeros() != 288) {
     111           0 :         SPDLOG_ERROR("Wrong number of non-zeros.");
     112           0 :         success = false;
     113             :     }
     114             : 
     115           1 :     if (dipole_p.get_matrix().nonZeros() != 288) {
     116           0 :         SPDLOG_ERROR("Wrong number of non-zeros.");
     117           0 :         success = false;
     118             :     }
     119             : 
     120           1 :     if (dipole_m.get_matrix().nonZeros() != 288) {
     121           0 :         SPDLOG_ERROR("Wrong number of non-zeros.");
     122           0 :         success = false;
     123             :     }
     124             : 
     125           1 :     return success ? 0 : 1;
     126           1 : }