Line data    Source code

       1             : /*
       2             :  * Copyright (c) 2017 Sebastian Weber, Henri Menke. All rights reserved.
       3             :  *
       4             :  * This file is part of the pairinteraction library.
       5             :  *
       6             :  * The pairinteraction library is free software: you can redistribute it and/or modify
       7             :  * it under the terms of the GNU Lesser General Public License as published by
       8             :  * the Free Software Foundation, either version 3 of the License, or
       9             :  * (at your option) any later version.
      10             :  *
      11             :  * The pairinteraction library is distributed in the hope that it will be useful,
      12             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      13             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      14             :  * GNU Lesser General Public License for more details.
      15             :  *
      16             :  * You should have received a copy of the GNU Lesser General Public License
      17             :  * along with the pairinteraction library. If not, see <http://www.gnu.org/licenses/>.
      18             :  */
      19             : 
      20             : #include "MatrixElementCache.hpp"
      21             : #include "State.hpp"
      22             : #include "SystemOne.hpp"
      23             : #include "SystemTwo.hpp"
      24             : #include "filesystem.hpp"
      25             : 
      26             : #include <cereal/archives/json.hpp>
      27             : #define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
      28             : #include <doctest/doctest.h>
      29             : 
      30             : #include <algorithm>
      31             : #include <fstream>
      32             : #include <iostream>
      33             : #include <memory>
      34             : #include <stdexcept>
      35             : 
      36             : struct F {
      37           1 :     F() : path_cache(fs::create_temp_directory()) {
      38           2 :         std::cout << "Cache directory " << fs::absolute(path_cache).string() << " created."
      39           1 :                   << std::endl;
      40           1 :     }
      41           1 :     ~F() {
      42             :         // Delete cache directory
      43           1 :         fs::remove_all(path_cache);
      44           1 :     }
      45             :     fs::path path_cache;
      46             : };
      47             : 
      48           2 : TEST_CASE_FIXTURE(F, "integration_test") // NOLINT
      49             : {
      50             :     using Scalar = double;
      51           1 :     constexpr bool dump_new_reference_data = false;
      52           1 :     constexpr double tolerance = 1e-6;
      53             : 
      54             :     ////////////////////////////////////////////////////////////////////
      55             :     /// Preparations ///////////////////////////////////////////////////
      56             :     ////////////////////////////////////////////////////////////////////
      57             : 
      58             :     // Set up cache
      59           2 :     MatrixElementCache cache(path_cache.string());
      60             : 
      61             :     // Load reference data
      62           2 :     Eigen::SparseMatrix<Scalar> hamiltonian_one_reference, hamiltonian_two_reference;
      63           2 :     Eigen::SparseMatrix<Scalar> basis_one_reference, basis_two_reference;
      64             : 
      65             :     if (!dump_new_reference_data) {
      66           2 :         std::ifstream ifs("./pairinteraction/unit_test/integration_test_referencedata.json");
      67           1 :         cereal::JSONInputArchive ia(ifs);
      68           1 :         ia >> cereal::make_nvp("hamiltonian_one", hamiltonian_one_reference) >>
      69           1 :             cereal::make_nvp("basis_one", basis_one_reference) >>
      70           1 :             cereal::make_nvp("hamiltonian_two", hamiltonian_two_reference) >>
      71           1 :             cereal::make_nvp("basis_two", basis_two_reference);
      72             :     }
      73             : 
      74             :     // Setup states
      75           3 :     StateOne state_one("Rb", 61, 2, 1.5, 1.5);
      76           3 :     StateTwo state_two(state_one, state_one);
      77             : 
      78             :     ////////////////////////////////////////////////////////////////////
      79             :     /// Test system consisting of one atom /////////////////////////////
      80             :     ////////////////////////////////////////////////////////////////////
      81             : 
      82             :     // Build one-atom system
      83           2 :     SystemOne<Scalar> system_one(state_one.getSpecies(), cache);
      84           1 :     system_one.restrictEnergy(state_one.getEnergy() - 40, state_one.getEnergy() + 40);
      85           1 :     system_one.restrictN(state_one.getN() - 1, state_one.getN() + 1);
      86           1 :     system_one.restrictL(state_one.getL() - 1, state_one.getL() + 1);
      87           1 :     system_one.setEfield({{0, 0, 0.1}});
      88           1 :     system_one.setBfield({{0, 0, 1}});
      89           1 :     system_one.enableDiamagnetism(false);
      90             : 
      91             :     // Check for correct dimensions
      92           1 :     CHECK(system_one.getNumBasisvectors() == 64);
      93           1 :     CHECK(system_one.getNumStates() == 64);
      94             : 
      95             :     // Compare current results to the reference data (the results have to be
      96             :     // compared before diagonalization as the order of the eigenvectors is not
      97             :     // fixed)
      98           2 :     Eigen::SparseMatrix<Scalar> hamiltonian_one = system_one.getHamiltonian();
      99           2 :     Eigen::SparseMatrix<Scalar> basis_one = system_one.getBasisvectors();
     100             : 
     101           2 :     Eigen::SparseMatrix<double> diff;
     102             :     double max_diff_hamiltonian, max_diff_basis;
     103             : 
     104             :     if (!dump_new_reference_data) {
     105           1 :         diff = (hamiltonian_one - hamiltonian_one_reference)
     106           1 :                    .pruned(tolerance, 1) // without pruning, max_diff_hamiltonian
     107             :                                          // might be infinity due to division by
     108             :                                          // zero
     109           1 :                    .cwiseQuotient(hamiltonian_one.cwiseMin(hamiltonian_one_reference))
     110           1 :                    .cwiseAbs();
     111           1 :         max_diff_hamiltonian =
     112           1 :             *std::max_element(diff.valuePtr(), diff.valuePtr() + diff.nonZeros());
     113             :         std::cout << "One-atom system, relative maximum deviation from "
     114           1 :                      "reference Hamiltonian: "
     115           1 :                   << max_diff_hamiltonian << std::endl;
     116           1 :         CHECK(max_diff_hamiltonian == doctest::Approx(0.0).epsilon(tolerance));
     117             : 
     118           1 :         diff = (basis_one - basis_one_reference)
     119           1 :                    .pruned(tolerance, 1) // without pruning, max_diff_hamiltonian
     120             :                                          // might be infinity due to division by
     121             :                                          // zero
     122           1 :                    .cwiseQuotient(basis_one.cwiseMin(basis_one_reference))
     123           1 :                    .cwiseAbs();
     124           1 :         max_diff_basis = *std::max_element(diff.valuePtr(), diff.valuePtr() + diff.nonZeros());
     125             :         std::cout << "One-atom system, relative maximum deviation from "
     126           1 :                      "reference basis: "
     127           1 :                   << max_diff_basis << std::endl;
     128           1 :         CHECK(max_diff_basis == doctest::Approx(0.0).epsilon(tolerance));
     129             :     }
     130             : 
     131             :     // Diagonalize one-atom system
     132           1 :     system_one.diagonalize();
     133             : 
     134             :     ////////////////////////////////////////////////////////////////////
     135             :     /// Test system consisting of two atoms ////////////////////////////
     136             :     ////////////////////////////////////////////////////////////////////
     137             : 
     138             :     // Build one-atom system (for this test, system_one has to be diagonal by
     139             :     // itself because diagonalization can lead to different order of
     140             :     // eigenvectors)
     141             :     // system_one = SystemOne(state_one.species, cache); // TODO  object of type 'SystemOne' cannot
     142             :     // be assigned because its copy assignment operator is implicitly deleted
     143           2 :     SystemOne<Scalar> system_one_new(state_one.getSpecies(), cache);
     144           1 :     system_one_new.restrictEnergy(state_one.getEnergy() - 40, state_one.getEnergy() + 40);
     145           1 :     system_one_new.restrictN(state_one.getN() - 1, state_one.getN() + 1);
     146           1 :     system_one_new.restrictL(state_one.getL() - 1, state_one.getL() + 1);
     147             : 
     148             :     // Build two-atom system
     149           2 :     SystemTwo<Scalar> system_two(system_one_new, system_one_new, cache);
     150           1 :     system_two.restrictEnergy(state_two.getEnergy() - 2, state_two.getEnergy() + 2);
     151           1 :     system_two.setConservedParityUnderPermutation(ODD);
     152           1 :     system_two.setDistance(6);
     153           1 :     system_two.setAngle(0.9);
     154             : 
     155             :     // Check for correct dimensions
     156           1 :     CHECK(system_two.getNumBasisvectors() == 239);
     157           1 :     CHECK(system_two.getNumStates() == 468);
     158             : 
     159             :     // Compare current results to the reference data (the results have to be
     160             :     // compared before diagonalization as the order of the eigenvectors is not
     161             :     // fixed)
     162           2 :     Eigen::SparseMatrix<Scalar> hamiltonian_two = system_two.getHamiltonian();
     163           2 :     Eigen::SparseMatrix<Scalar> basis_two = system_two.getBasisvectors();
     164             : 
     165             :     if (!dump_new_reference_data) {
     166           1 :         diff = (hamiltonian_two - hamiltonian_two_reference)
     167           1 :                    .pruned(tolerance, 1) // without pruning, max_diff_hamiltonian
     168             :                                          // might be infinity due to division by
     169             :                                          // zero
     170           1 :                    .cwiseQuotient(hamiltonian_two.cwiseMin(hamiltonian_two_reference))
     171           1 :                    .cwiseAbs();
     172           1 :         max_diff_hamiltonian =
     173           1 :             *std::max_element(diff.valuePtr(), diff.valuePtr() + diff.nonZeros());
     174             :         std::cout << "Two-atom system, relative maximum deviation from "
     175           1 :                      "reference Hamiltonian: "
     176           1 :                   << max_diff_hamiltonian << std::endl;
     177           1 :         CHECK(max_diff_hamiltonian == doctest::Approx(0.0).epsilon(tolerance));
     178             : 
     179           1 :         diff = (basis_two - basis_two_reference)
     180           1 :                    .pruned(tolerance, 1) // without pruning, max_diff_hamiltonian
     181             :                                          // might be infinity due to division by
     182             :                                          // zero
     183           1 :                    .cwiseQuotient(basis_two.cwiseMin(basis_two_reference))
     184           1 :                    .cwiseAbs();
     185           1 :         max_diff_basis = *std::max_element(diff.valuePtr(), diff.valuePtr() + diff.nonZeros());
     186             :         std::cout << "Two-atom system, relative maximum deviation from "
     187           1 :                      "reference basis: "
     188           1 :                   << max_diff_basis << std::endl;
     189           1 :         CHECK(max_diff_basis == doctest::Approx(0.0).epsilon(tolerance));
     190             :     }
     191             : 
     192             :     // Diagonalize two-atom system
     193           1 :     system_two.diagonalize();
     194             : 
     195             :     ////////////////////////////////////////////////////////////////////
     196             :     /// Clean up ///////////////////////////////////////////////////////
     197             :     ////////////////////////////////////////////////////////////////////
     198             : 
     199             :     if (dump_new_reference_data) {
     200             :         {
     201             :             std::ofstream ofs("../pairinteraction/unit_test/integration_test_referencedata.json");
     202             :             cereal::JSONOutputArchive oa(ofs);
     203             :             oa << CEREAL_NVP(hamiltonian_one) << CEREAL_NVP(basis_one)
     204             :                << CEREAL_NVP(hamiltonian_two) << CEREAL_NVP(basis_two);
     205             :         }
     206             : 
     207             :         // ATTENTION
     208             :         // After generating integration_test_referencedata.txt, we possibly have to manually modify
     209             :         // the Boost serialization library version - the beginning of the file should read "22
     210             :         // serialization::archive 12". Otherwise, unsupported version exceptions are thrown with
     211             :         // older versions of Boost even though they are compatible.
     212             : 
     213             :         // We deliberately crash the test in case we are dumping new
     214             :         // data because we want to prevent accidentally dumping new
     215             :         // data when running in Travis CI.
     216             :         FAIL("No tests were executed. Only dumping data!");
     217             :     }
     218             : 
     219             :     // TODO call more methods to increase code covering
     220             :     // TODO cause exceptions and check whether they are handled correctly
     221           1 : }