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 UM_IN_ATOMIC_UNITS = 1 / 5.29177210544e-5;
      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 and diagonalize systems for two atoms
      38           1 :     auto basis = pairinteraction::BasisAtomCreator<double>()
      39           2 :                      .set_species("Rb")
      40           1 :                      .restrict_quantum_number_n(58, 62)
      41           1 :                      .restrict_quantum_number_l(0, 2)
      42           1 :                      .create(database);
      43           2 :     SPDLOG_INFO("Number of single-atom basis states: {}", basis->get_number_of_states());
      44             : 
      45           1 :     pairinteraction::SystemAtom<double> system(basis);
      46             : 
      47             :     // Create two-atom systems for different interatomic distances
      48           1 :     auto ket = pairinteraction::KetAtomCreator()
      49           2 :                    .set_species("Rb")
      50           1 :                    .set_quantum_number_n(60)
      51           1 :                    .set_quantum_number_l(0)
      52           1 :                    .set_quantum_number_m(0.5)
      53           1 :                    .create(database);
      54           1 :     double min_energy = 2 * ket->get_energy() - 3 / HARTREE_IN_GHZ;
      55           1 :     double max_energy = 2 * ket->get_energy() + 3 / HARTREE_IN_GHZ;
      56             : 
      57           2 :     auto basis_pair = pairinteraction::BasisPairCreator<double>()
      58           1 :                           .add(system)
      59           1 :                           .add(system)
      60           1 :                           .restrict_energy(min_energy, max_energy)
      61           1 :                           .restrict_quantum_number_m(1, 1)
      62           1 :                           .create();
      63           2 :     SPDLOG_INFO("Number of two-atom basis states: {}", basis_pair->get_number_of_states());
      64             : 
      65           1 :     std::vector<pairinteraction::SystemPair<double>> system_pairs;
      66           1 :     system_pairs.reserve(5);
      67           6 :     for (int i = 1; i < 6; ++i) {
      68           5 :         pairinteraction::SystemPair<double> system(basis_pair);
      69           5 :         system.set_distance_vector({0, 0, i * UM_IN_ATOMIC_UNITS});
      70           5 :         system_pairs.push_back(std::move(system));
      71           5 :     }
      72             : 
      73             :     // Diagonalize the systems in parallel
      74           1 :     pairinteraction::DiagonalizerEigen<double> diagonalizer(pairinteraction::FloatType::FLOAT32);
      75           1 :     pairinteraction::diagonalize(system_pairs, diagonalizer);
      76             : 
      77             :     // Sort by the eigenenergies
      78           6 :     for (auto &system : system_pairs) {
      79           5 :         auto sorter = system.get_sorter({pairinteraction::TransformationType::SORT_BY_ENERGY});
      80           5 :         system.transform(sorter);
      81           5 :     }
      82             : 
      83             :     // Extract results
      84           1 :     std::vector<std::string> kets;
      85           1 :     Eigen::MatrixX<double> eigenenergies(system_pairs.size(), basis_pair->get_number_of_states());
      86           0 :     Eigen::MatrixX<double> eigenstates(system_pairs.size(),
      87           1 :                                        basis_pair->get_number_of_states() *
      88           1 :                                            basis_pair->get_number_of_states());
      89           1 :     Eigen::MatrixX<double> overlaps(system_pairs.size(), basis_pair->get_number_of_states());
      90             : 
      91           1 :     kets.reserve(basis->get_number_of_states());
      92          20 :     for (const auto &ket : *basis_pair) {
      93          19 :         std::stringstream ss;
      94          19 :         ss << *ket;
      95          19 :         kets.push_back(ss.str());
      96          19 :     }
      97             : 
      98           6 :     for (Eigen::Index i = 0; i < static_cast<Eigen::Index>(system_pairs.size()); ++i) {
      99           5 :         eigenenergies.row(i) = system_pairs[i].get_eigenenergies() * HARTREE_IN_GHZ;
     100             : 
     101             :         Eigen::MatrixX<double> tmp =
     102           5 :             system_pairs[i].get_eigenbasis()->get_coefficients().toDense().transpose();
     103           5 :         eigenstates.row(i) = Eigen::Map<Eigen::VectorXd>(tmp.data(), tmp.size());
     104             : 
     105           5 :         overlaps.row(i) = system_pairs[i].get_eigenbasis()->get_overlaps(ket, ket);
     106           5 :     }
     107             : 
     108             :     // Compare with reference data
     109           1 :     bool success = true;
     110             : 
     111             :     // Check kets
     112             :     const std::filesystem::path reference_kets_file =
     113           1 :         data_dir / "reference_pair_potential/kets.txt";
     114           2 :     SPDLOG_INFO("Reference kets: {}", reference_kets_file.string());
     115             : 
     116           1 :     std::vector<std::string> reference_kets;
     117           1 :     std::string line;
     118           1 :     std::ifstream stream(reference_kets_file);
     119           1 :     if (!stream) {
     120           0 :         SPDLOG_ERROR("Could not open reference kets file");
     121           0 :         success = false;
     122             :     } else {
     123           1 :         reference_kets.reserve(basis_pair->get_number_of_states());
     124          20 :         while (std::getline(stream, line)) {
     125          19 :             reference_kets.push_back(line);
     126             :         }
     127           1 :         stream.close();
     128           1 :         if (kets.size() != reference_kets.size()) {
     129           0 :             SPDLOG_ERROR("Number of kets does not match reference data");
     130           0 :             success = false;
     131           1 :         } else if (kets != reference_kets) {
     132           0 :             for (size_t i = 0; i < kets.size(); ++i) {
     133           0 :                 SPDLOG_DEBUG("Ket: {} vs {}, match: {}", kets[i], reference_kets[i],
     134             :                              kets[i] == reference_kets[i]);
     135             :             }
     136           0 :             SPDLOG_ERROR("Kets do not match reference data");
     137           0 :             success = false;
     138             :         }
     139             :     }
     140             : 
     141             :     // Check eigenenergies
     142             :     const std::filesystem::path reference_eigenenergies_file =
     143           1 :         data_dir / "reference_pair_potential/eigenenergies.txt";
     144           2 :     SPDLOG_INFO("Reference eigenenergies: {}", reference_eigenenergies_file.string());
     145             : 
     146           0 :     Eigen::MatrixXd reference_eigenenergies(system_pairs.size(),
     147           1 :                                             basis_pair->get_number_of_states());
     148           1 :     stream = std::ifstream(reference_eigenenergies_file);
     149           1 :     if (!stream) {
     150           0 :         SPDLOG_ERROR("Could not open reference eigenenergies file");
     151           0 :         success = false;
     152             :     } else {
     153           6 :         for (Eigen::Index i = 0; i < static_cast<Eigen::Index>(system_pairs.size()); ++i) {
     154         100 :             for (Eigen::Index j = 0;
     155         100 :                  j < static_cast<Eigen::Index>(basis_pair->get_number_of_states()); ++j) {
     156          95 :                 stream >> reference_eigenenergies(i, j);
     157             :             }
     158             :         }
     159           1 :         stream.close();
     160           1 :         if ((eigenenergies - reference_eigenenergies).array().abs().maxCoeff() >
     161             :             100e-6) { // 100 kHz precision
     162           0 :             for (Eigen::Index i = 0; i < eigenenergies.size(); ++i) {
     163           0 :                 SPDLOG_DEBUG("Eigenenergy: {} vs {}, delta: {}", eigenenergies(i),
     164             :                              reference_eigenenergies(i),
     165             :                              std::abs(eigenenergies(i) - reference_eigenenergies(i)));
     166             :             }
     167           0 :             SPDLOG_ERROR("Eigenenergies do not match reference data");
     168           0 :             success = false;
     169             :         }
     170             :     }
     171             : 
     172             :     // Check overlaps
     173             :     const std::filesystem::path reference_overlaps_file =
     174           1 :         data_dir / "reference_pair_potential/overlaps.txt";
     175           2 :     SPDLOG_INFO("Reference overlaps: {}", reference_overlaps_file.string());
     176             : 
     177           1 :     Eigen::MatrixXd reference_overlaps(system_pairs.size(), basis_pair->get_number_of_states());
     178           1 :     stream = std::ifstream(reference_overlaps_file);
     179           1 :     if (!stream) {
     180           0 :         SPDLOG_ERROR("Could not open reference overlap file");
     181           0 :         success = false;
     182             :     } else {
     183           6 :         for (Eigen::Index i = 0; i < static_cast<Eigen::Index>(system_pairs.size()); ++i) {
     184         100 :             for (Eigen::Index j = 0;
     185         100 :                  j < static_cast<Eigen::Index>(basis_pair->get_number_of_states()); ++j) {
     186          95 :                 stream >> reference_overlaps(i, j);
     187             :             }
     188             :         }
     189           1 :         stream.close();
     190           1 :         if ((overlaps - reference_overlaps).array().abs().maxCoeff() > 1e-5) {
     191           0 :             for (Eigen::Index i = 0; i < overlaps.size(); ++i) {
     192           0 :                 SPDLOG_DEBUG("Overlap: {} vs {}, delta: {}", overlaps(i), reference_overlaps(i),
     193             :                              std::abs(overlaps(i) - reference_overlaps(i)));
     194             :             }
     195           0 :             SPDLOG_ERROR("Overlaps do not match reference data");
     196           0 :             success = false;
     197             :         }
     198             :     }
     199             : 
     200             :     // Check eigenstates
     201             :     // Because of degeneracies, checking the eigenstates against reference data is complicated.
     202             :     // Thus, we only check their normalization and orthogonality.
     203             :     Eigen::VectorXd cumulative_norm =
     204           1 :         (eigenstates.adjoint().array() * eigenstates.transpose().array()).colwise().sum();
     205           1 :     if (!cumulative_norm.isApprox(Eigen::VectorXd::Constant(5, 19), 1e-5)) {
     206           0 :         SPDLOG_ERROR("Eigenvectors are not orthonormal.");
     207           0 :         success = false;
     208             :     }
     209             : 
     210           1 :     return success ? 0 : 1;
     211           1 : }