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 <Eigen/Dense>
       8             : #include <filesystem>
       9             : #include <fstream>
      10             : #include <spdlog/spdlog.h>
      11             : #include <vector>
      12             : 
      13             : constexpr double HARTREE_IN_GHZ = 6579683.920501762;
      14             : constexpr double VOLT_PER_CM_IN_ATOMIC_UNITS = 1 / 5.14220675112e9;
      15             : 
      16           1 : int main(int argc, char **argv) {
      17             :     // Call the setup function to configure logging
      18           1 :     pairinteraction::setup();
      19             : 
      20             :     // Create a database instance
      21           1 :     std::filesystem::path database_dir;
      22           1 :     std::filesystem::path data_dir;
      23           1 :     bool download_missing = false;
      24             : 
      25           3 :     for (int i = 1; i < argc; ++i) {
      26           2 :         bool found = pairinteraction::args::parse_download_missing(i, argc, argv, download_missing);
      27           2 :         if (!found) {
      28           2 :             found = pairinteraction::args::parse_database_dir(i, argc, argv, database_dir);
      29             :         }
      30           2 :         if (!found) {
      31           1 :             pairinteraction::args::parse_data_dir(i, argc, argv, data_dir);
      32             :         }
      33             :     }
      34             : 
      35           1 :     pairinteraction::Database database(download_missing, true, database_dir);
      36             : 
      37             :     // Create a basis
      38           1 :     auto ket = pairinteraction::KetAtomCreator()
      39           2 :                    .set_species("Rb")
      40           1 :                    .set_quantum_number_n(60)
      41           1 :                    .set_quantum_number_l(0)
      42           1 :                    .set_quantum_number_m(0.5)
      43           1 :                    .create(database);
      44             : 
      45           1 :     auto basis = pairinteraction::BasisAtomCreator<double>()
      46           2 :                      .set_species("Rb")
      47           1 :                      .restrict_quantum_number_n(58, 62)
      48           1 :                      .restrict_quantum_number_l(0, 2)
      49           1 :                      .create(database);
      50             : 
      51           2 :     SPDLOG_INFO("Number of basis states: {}", basis->get_number_of_states());
      52             : 
      53             :     // Create systems for different values of the electric field
      54           1 :     std::vector<pairinteraction::SystemAtom<double>> systems;
      55           1 :     systems.reserve(11);
      56          12 :     for (int i = 0; i < 11; ++i) {
      57          11 :         pairinteraction::SystemAtom<double> system(basis);
      58          11 :         system.set_electric_field({0, 0, i * VOLT_PER_CM_IN_ATOMIC_UNITS});
      59          11 :         systems.push_back(std::move(system));
      60          11 :     }
      61             : 
      62             :     // Diagonalize the systems in parallel
      63           1 :     pairinteraction::DiagonalizerEigen<double> diagonalizer(pairinteraction::FloatType::FLOAT32);
      64           1 :     pairinteraction::diagonalize(systems, diagonalizer);
      65             : 
      66             :     // Sort by the eigenenergies
      67          12 :     for (auto &system : systems) {
      68          11 :         auto sorter = system.get_sorter({pairinteraction::TransformationType::SORT_BY_ENERGY});
      69          11 :         system.transform(sorter);
      70          11 :     }
      71             : 
      72             :     // Extract results
      73           1 :     std::vector<std::string> kets;
      74           1 :     Eigen::MatrixX<double> eigenenergies(systems.size(), basis->get_number_of_states());
      75             :     Eigen::MatrixX<double> eigenstates(
      76           1 :         systems.size(), basis->get_number_of_states() * basis->get_number_of_states());
      77           1 :     Eigen::MatrixX<double> overlaps(systems.size(), basis->get_number_of_states());
      78             : 
      79           1 :     kets.reserve(basis->get_number_of_states());
      80          91 :     for (const auto &ket : *basis) {
      81          90 :         std::stringstream ss;
      82          90 :         ss << *ket;
      83          90 :         kets.push_back(ss.str());
      84          90 :     }
      85             : 
      86          12 :     for (Eigen::Index i = 0; i < static_cast<Eigen::Index>(systems.size()); ++i) {
      87          11 :         eigenenergies.row(i) = systems[i].get_eigenenergies() * HARTREE_IN_GHZ;
      88             : 
      89             :         Eigen::MatrixX<double> tmp =
      90          11 :             systems[i].get_eigenbasis()->get_coefficients().toDense().transpose();
      91          11 :         eigenstates.row(i) = Eigen::Map<Eigen::VectorXd>(tmp.data(), tmp.size());
      92             : 
      93          11 :         overlaps.row(i) = systems[i].get_eigenbasis()->get_overlaps(ket);
      94          11 :     }
      95             : 
      96             :     // Compare with reference data
      97           1 :     bool success = true;
      98             : 
      99             :     // Check kets
     100           1 :     const std::filesystem::path reference_kets_file = data_dir / "reference_stark_map/kets.txt";
     101           2 :     SPDLOG_INFO("Reference kets: {}", reference_kets_file.string());
     102             : 
     103           1 :     std::vector<std::string> reference_kets;
     104           1 :     std::string line;
     105           1 :     std::ifstream stream(reference_kets_file);
     106           1 :     if (!stream) {
     107           0 :         SPDLOG_ERROR("Could not open reference kets file");
     108           0 :         success = false;
     109             :     } else {
     110           1 :         reference_kets.reserve(basis->get_number_of_states());
     111          91 :         while (std::getline(stream, line)) {
     112          90 :             reference_kets.push_back(line);
     113             :         }
     114           1 :         stream.close();
     115           1 :         if (kets.size() != reference_kets.size()) {
     116           0 :             SPDLOG_ERROR("Number of kets does not match reference data");
     117           0 :             success = false;
     118           1 :         } else if (kets != reference_kets) {
     119           0 :             for (size_t i = 0; i < kets.size(); ++i) {
     120           0 :                 SPDLOG_DEBUG("Ket: {} vs {}, match: {}", kets[i], reference_kets[i],
     121             :                              kets[i] == reference_kets[i]);
     122             :             }
     123           0 :             SPDLOG_ERROR("Kets do not match reference data");
     124           0 :             success = false;
     125             :         }
     126             :     }
     127             : 
     128             :     // Check eigenenergies
     129             :     const std::filesystem::path reference_eigenenergies_file =
     130           1 :         data_dir / "reference_stark_map/eigenenergies.txt";
     131           2 :     SPDLOG_INFO("Reference eigenenergies: {}", reference_eigenenergies_file.string());
     132             : 
     133           1 :     Eigen::MatrixXd reference_eigenenergies(systems.size(), basis->get_number_of_states());
     134           1 :     stream = std::ifstream(reference_eigenenergies_file);
     135           1 :     if (!stream) {
     136           0 :         SPDLOG_ERROR("Could not open reference eigenenergies file");
     137           0 :         success = false;
     138             :     } else {
     139          12 :         for (Eigen::Index i = 0; i < static_cast<Eigen::Index>(systems.size()); ++i) {
     140        1001 :             for (Eigen::Index j = 0; j < static_cast<Eigen::Index>(basis->get_number_of_states());
     141             :                  ++j) {
     142         990 :                 stream >> reference_eigenenergies(i, j);
     143             :             }
     144             :         }
     145           1 :         stream.close();
     146           1 :         if ((eigenenergies - reference_eigenenergies).array().abs().maxCoeff() >
     147             :             500e-6) { // 500 kHz precision
     148           0 :             for (Eigen::Index i = 0; i < eigenenergies.size(); ++i) {
     149           0 :                 SPDLOG_DEBUG("Eigenenergy: {} vs {}, delta: {}", eigenenergies(i),
     150             :                              reference_eigenenergies(i),
     151             :                              std::abs(eigenenergies(i) - reference_eigenenergies(i)));
     152             :             }
     153           0 :             SPDLOG_ERROR("Eigenenergies do not match reference data");
     154           0 :             success = false;
     155             :         }
     156             :     }
     157             : 
     158             :     // Check overlaps
     159             :     const std::filesystem::path reference_overlaps_file =
     160           1 :         data_dir / "reference_stark_map/overlaps.txt";
     161           2 :     SPDLOG_INFO("Reference overlaps: {}", reference_overlaps_file.string());
     162             : 
     163           1 :     Eigen::MatrixXd reference_overlaps(systems.size(), basis->get_number_of_states());
     164           1 :     stream = std::ifstream(reference_overlaps_file);
     165           1 :     if (!stream) {
     166           0 :         SPDLOG_ERROR("Could not open reference overlap file");
     167           0 :         success = false;
     168             :     } else {
     169          12 :         for (Eigen::Index i = 0; i < static_cast<Eigen::Index>(systems.size()); ++i) {
     170        1001 :             for (Eigen::Index j = 0; j < static_cast<Eigen::Index>(basis->get_number_of_states());
     171             :                  ++j) {
     172         990 :                 stream >> reference_overlaps(i, j);
     173             :             }
     174             :         }
     175           1 :         stream.close();
     176           1 :         if ((overlaps - reference_overlaps).array().abs().maxCoeff() > 5e-5) {
     177           0 :             for (Eigen::Index i = 0; i < overlaps.size(); ++i) {
     178           0 :                 SPDLOG_DEBUG("Overlap: {} vs {}, delta: {}", overlaps(i), reference_overlaps(i),
     179             :                              std::abs(overlaps(i) - reference_overlaps(i)));
     180             :             }
     181           0 :             SPDLOG_ERROR("Overlaps do not match reference data");
     182           0 :             success = false;
     183             :         }
     184             :     }
     185             : 
     186             :     // Check eigenstates
     187             :     // Because of degeneracies, checking the eigenstates against reference data is complicated.
     188             :     // Thus, we only check their normalization and orthogonality.
     189             :     Eigen::VectorXd cumulative_norm =
     190           1 :         (eigenstates.adjoint().array() * eigenstates.transpose().array()).colwise().sum();
     191           1 :     if (!cumulative_norm.isApprox(Eigen::VectorXd::Constant(11, 90), 5e-5)) {
     192           0 :         SPDLOG_ERROR("Eigenvectors are not orthonormal.");
     193           0 :         success = false;
     194             :     }
     195             : 
     196           1 :     return success ? 0 : 1;
     197           1 : }