From f8c9a0532374c182b47b34de6c77237e1156a5d7 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Sun, 21 Oct 2018 21:30:04 +0200 Subject: [PATCH 01/13] added the selected fock space --- cmake/SetCMakeEnvironment.cmake | 4 + include/Configuration.hpp | 21 +++++ include/FockSpace/BaseFockSpace.hpp | 2 +- include/FockSpace/FockSpaceType.hpp | 3 +- include/FockSpace/SelectedFockSpace.hpp | 82 +++++++++++++++++ include/ONV.hpp | 18 ++++ src/FockSpace/SelectedFockSpace.cpp | 102 +++++++++++++++++++++ src/ONV.cpp | 47 +++++++++- tests/FockSpace/SelectedFockSpace_test.cpp | 55 +++++++++++ tests/ONV_test.cpp | 33 +++++++ 10 files changed, 363 insertions(+), 4 deletions(-) create mode 100644 include/Configuration.hpp create mode 100644 include/FockSpace/SelectedFockSpace.hpp create mode 100644 src/FockSpace/SelectedFockSpace.cpp create mode 100644 tests/FockSpace/SelectedFockSpace_test.cpp diff --git a/cmake/SetCMakeEnvironment.cmake b/cmake/SetCMakeEnvironment.cmake index 546bb74d6e..9ceaf1faaf 100755 --- a/cmake/SetCMakeEnvironment.cmake +++ b/cmake/SetCMakeEnvironment.cmake @@ -31,6 +31,7 @@ set(PROJECT_SOURCE_FILES ${PROJECT_SOURCE_FOLDER}/FockSpace/BaseFockSpace.cpp ${PROJECT_SOURCE_FOLDER}/FockSpace/FockSpace.cpp ${PROJECT_SOURCE_FOLDER}/FockSpace/FockSpaceProduct.cpp + ${PROJECT_SOURCE_FOLDER}/FockSpace/SelectedFockSpace.cpp ${PROJECT_SOURCE_FOLDER}/HamiltonianBuilder/DOCI.cpp ${PROJECT_SOURCE_FOLDER}/HamiltonianBuilder/FCI.cpp ${PROJECT_SOURCE_FOLDER}/HamiltonianBuilder/HamiltonianBuilder.cpp @@ -78,6 +79,7 @@ set(PROJECT_INCLUDE_FILES ${PROJECT_INCLUDE_FOLDER}/FockSpace/FockSpace.hpp ${PROJECT_INCLUDE_FOLDER}/FockSpace/FockSpaceProduct.hpp ${PROJECT_INCLUDE_FOLDER}/FockSpace/FockSpaceType.hpp + ${PROJECT_INCLUDE_FOLDER}/FockSpace/SelectedFockSpace.hpp ${PROJECT_INCLUDE_FOLDER}/HamiltonianBuilder/DOCI.hpp ${PROJECT_INCLUDE_FOLDER}/HamiltonianBuilder/FCI.hpp ${PROJECT_INCLUDE_FOLDER}/HamiltonianBuilder/HamiltonianBuilder.hpp @@ -102,6 +104,7 @@ set(PROJECT_INCLUDE_FILES ${PROJECT_INCLUDE_FOLDER}/AOBasis.hpp ${PROJECT_INCLUDE_FOLDER}/Atom.hpp ${PROJECT_INCLUDE_FOLDER}/common.hpp + ${PROJECT_INCLUDE_FOLDER}/Configuration.hpp ${PROJECT_INCLUDE_FOLDER}/DOCINewtonOrbitalOptimizer.hpp ${PROJECT_INCLUDE_FOLDER}/elements.hpp ${PROJECT_INCLUDE_FOLDER}/JacobiRotationParameters.hpp @@ -129,6 +132,7 @@ set(PROJECT_TEST_SOURCE_FILES ${PROJECT_TESTS_FOLDER}/AP1roG/AP1roGPSESolver_test.cpp ${PROJECT_TESTS_FOLDER}/FockSpace/FockSpace_test.cpp ${PROJECT_TESTS_FOLDER}/FockSpace/FockSpaceProduct_test.cpp + ${PROJECT_TESTS_FOLDER}/FockSpace/SelectedFockSpace_test.cpp ${PROJECT_TESTS_FOLDER}/HamiltonianBuilder/DOCI_test.cpp ${PROJECT_TESTS_FOLDER}/HamiltonianBuilder/FCI_test.cpp ${PROJECT_TESTS_FOLDER}/HamiltonianParameters/HamiltonianParameters_test.cpp diff --git a/include/Configuration.hpp b/include/Configuration.hpp new file mode 100644 index 0000000000..ef420c3150 --- /dev/null +++ b/include/Configuration.hpp @@ -0,0 +1,21 @@ +#ifndef GQCP_CONFIGURATION_HPP +#define GQCP_CONFIGURATION_HPP + + +#include "ONV.hpp" +#include "common.hpp" + + +namespace GQCP { + + +struct Configuration { + ONV onv_alpha; + ONV onv_beta; +}; + + +} // namespace GQCP + + +#endif // GQCP_CONFIGURATION_HPP diff --git a/include/FockSpace/BaseFockSpace.hpp b/include/FockSpace/BaseFockSpace.hpp index f1bdf14374..5350d3f4e2 100644 --- a/include/FockSpace/BaseFockSpace.hpp +++ b/include/FockSpace/BaseFockSpace.hpp @@ -37,7 +37,7 @@ namespace GQCP { class BaseFockSpace { protected: const size_t K; // number of spatial orbitals - const size_t dim; // dimension of the Fock space + size_t dim; // dimension of the Fock space // PROTECTED CONSTRUCTORS diff --git a/include/FockSpace/FockSpaceType.hpp b/include/FockSpace/FockSpaceType.hpp index 4cbfb73717..6d24c5bd3c 100644 --- a/include/FockSpace/FockSpaceType.hpp +++ b/include/FockSpace/FockSpaceType.hpp @@ -9,7 +9,8 @@ namespace GQCP { */ enum class FockSpaceType { FockSpace, - FockSpaceProduct + FockSpaceProduct, + SelectedFockSpace }; diff --git a/include/FockSpace/SelectedFockSpace.hpp b/include/FockSpace/SelectedFockSpace.hpp new file mode 100644 index 0000000000..9b4016134b --- /dev/null +++ b/include/FockSpace/SelectedFockSpace.hpp @@ -0,0 +1,82 @@ +// This file is part of GQCG-gqcp. +// +// Copyright (C) 2017-2018 the GQCG developers +// +// GQCG-gqcp is free software: you can redistribute it and/or modify +// it under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// GQCG-gqcp is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with GQCG-gqcp. If not, see . +// +#ifndef GQCP_SELECTEDSPACEPRODUCT_HPP +#define GQCP_SELECTEDSPACEPRODUCT_HPP + + +#include "FockSpace/BaseFockSpace.hpp" +#include "Configuration.hpp" + +#include +#include + + +namespace GQCP { + + +/** + * A Fock space for a given set of orbitals and number of alpha and beta electrons. + * Where the considered configurations are manually selected and represented as an ONV pair : + * a combination of two ONVs, one holding the alpha configuration and holding the bea configuration. + */ +class SelectedFockSpace : public GQCG::BaseFockSpace { +private: + const size_t N_alpha; // number of alpha electrons + const size_t N_beta; // number of beta electrons + + std::vector selection; // the selected configurations + + /** + * Given @param onv1 and onv2, make a configuration + */ + Configuration makeConfiguration(std::string onv1, std::string onv2); + +public: + // CONSTRUCTORS + /** + * Constructor given a @param K (spatial orbitals), N_alpha and N_beta (electrons); + * the initial dimension of the space is 0 as no selections are made. + */ + SelectedFockSpace(size_t K, size_t N_alpha, size_t N_beta); + + + // DESTRUCTORS + ~SelectedFockSpace() override = default; + + + // GETTERS + size_t get_N_alpha() const { return this->N_alpha; } + size_t get_N_beta() const { return this->N_beta; } + Configuration get_Configuration(size_t index) const { return this->selection[index]; } + FockSpaceType get_type() const override { return FockSpaceType::SelectedFockSpace; } + + + // PUBLIC METHODS + /** + * Add configuration(s) (@param onv1(s) and onv2(s)) to the current selected Fock space + * increasing its selected dimension. + */ + void addConfiguration(std::string onv1, std::string onv2); + void addConfiguration(std::vector onv1s, std::vector onv2s); +}; + + +} // namespace GQCP + + +#endif // GQCP_SELECTEDSPACEPRODUCT_HPP diff --git a/include/ONV.hpp b/include/ONV.hpp index 05561e39b2..5dd040abfc 100644 --- a/include/ONV.hpp +++ b/include/ONV.hpp @@ -55,6 +55,11 @@ class ONV { */ ONV(size_t K, size_t N, size_t unsigned_representation); + /** + * Constructor from a @param K orbitals and an @param unsigned_representation + */ + ONV(size_t K, size_t unsigned_representation); + // OPERATORS /** @@ -159,6 +164,19 @@ class ONV { size_t slice(size_t index_start, size_t index_end) const; + /** + * @return the number of different bits between this and @param other, i.e. two times the number of electron excitations + */ + size_t countNumberOfDifferences(const ONV& other) const; + + + /** + * @return the positions of the bits (lexical indices) that are occupied in this, but unoccupied in @param other + */ + std::vector findOccupiedDifferences(const ONV& other) const; + + + // FRIEND CLASSES friend class FockSpace; }; diff --git a/src/FockSpace/SelectedFockSpace.cpp b/src/FockSpace/SelectedFockSpace.cpp new file mode 100644 index 0000000000..b012856653 --- /dev/null +++ b/src/FockSpace/SelectedFockSpace.cpp @@ -0,0 +1,102 @@ +// This file is part of GQCG-gqcp. +// +// Copyright (C) 2017-2018 the GQCG developers +// +// GQCG-gqcp is free software: you can redistribute it and/or modify +// it under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// GQCG-gqcp is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with GQCG-gqcp. If not, see . +// +#include "FockSpace/SelectedFockSpace.hpp" + + +#include "boost/dynamic_bitset.hpp" + + +namespace GQCP { + + +/* + * PRIVATE METHODS + */ + +/** + * Given @param onv1 and onv2, make a configuration + */ +Configuration SelectedFockSpace::makeConfiguration(std::string onv1, std::string onv2){ + + boost::dynamic_bitset<> alpha_transfer (onv1); + boost::dynamic_bitset<> beta_transfer (onv2); + + if (alpha_transfer.size() != this->K | beta_transfer.size() != this->K) { + throw std::invalid_argument("Added configuration is not compatible with the orbitals number of the Fock space"); + } + + if (alpha_transfer.count() != this->N_alpha | beta_transfer.count() != this->N_beta) { + throw std::invalid_argument("Added configuration is not compatible with the electron numbers of the Fock space"); + } + + size_t alpha_s = alpha_transfer.to_ulong(); + size_t beta_s = beta_transfer.to_ulong(); + + ONV alpha (alpha_transfer.size(), alpha_s); + ONV beta (beta_transfer.size(), beta_s); + + return Configuration { alpha, beta }; +} + + + +/* + * CONSTRUCTORS + */ + +/** + * Constructor given a @param K (spatial orbitals), N_alpha and N_beta (electrons); + * the initial dimension of the space is 0 as no selections are made. + */ +SelectedFockSpace::SelectedFockSpace(size_t K, size_t N_alpha, size_t N_beta) : + BaseFockSpace(K, 0), + N_alpha (N_alpha), + N_beta (N_beta) +{} + + + +/* + * PUBLIC METHODS + */ + +/** + * Add configuration(s) (@param onv1(s) and onv2(s)) to the current selected Fock space + * increasing its selected dimension. + */ +void SelectedFockSpace::addConfiguration(std::string onv1, std::string onv2){ + + this->dim++; + + Configuration configuration = makeConfiguration(onv1, onv2); + selection.push_back(configuration); +} + +void SelectedFockSpace::addConfiguration(std::vector onv1s, std::vector onv2s){ + + if (onv1s.size() != onv2s.size()) { + throw std::invalid_argument("Size of both ONV entry vectors do not match"); + } + + for (int i = 0; iaddConfiguration(onv1s[i], onv2s[i]); + } +} + + +} // namespace GQCP diff --git a/src/ONV.cpp b/src/ONV.cpp index f27ff7fe23..0daa34a38a 100644 --- a/src/ONV.cpp +++ b/src/ONV.cpp @@ -15,7 +15,7 @@ // You should have received a copy of the GNU Lesser General Public License // along with GQCG-gqcp. If not, see . // -#include "ONV.hpp" +#include "../../gqcg/include/ONV.hpp" #include @@ -28,7 +28,7 @@ namespace GQCP { */ /** - * Constructor from a @param K orbitals, N electrons and a representation for the ONV + * Constructor from a @param K orbitals, N electrons and an @param unsigned_representation */ ONV::ONV(size_t K, size_t N, size_t representation): K(K), @@ -39,6 +39,17 @@ ONV::ONV(size_t K, size_t N, size_t representation): this->updateOccupationIndices(); // throws error if the representation and N are not compatible } +/** + * Constructor from a @param K orbitals and an @param unsigned_representation + */ +ONV::ONV(size_t K, size_t representation): + K(K), + N(__builtin_popcountl(representation)), + unsigned_representation(representation) +{ + occupation_indices = VectorXs::Zero(N); + this->updateOccupationIndices(); // throws error if the representation and N are not compatible +} /* @@ -246,4 +257,36 @@ size_t ONV::slice(size_t index_start, size_t index_end) const { } +/** + * @return the number of different bits between this and @param other, i.e. two times the number of electron excitations + */ +size_t ONV::countNumberOfDifferences(const ONV& other) const { + return __builtin_popcountl(this->unsigned_representation ^ other.unsigned_representation); +} + + +/** + * @return the positions of the bits (lexical indices) that are occupied in this, but unoccupied in @param other + */ +std::vector ONV::findOccupiedDifferences(const ONV& other) const { + + size_t differences = this->unsigned_representation ^ other.unsigned_representation; + size_t occupied_differences = differences & this->unsigned_representation; // this holds all indices occupied in this, but unoccupied in other + + size_t number_of_occupied_differences = __builtin_popcountl(occupied_differences); + std::vector positions (number_of_occupied_differences); + + + // Find the positions of the set bits in occupied_differences + for (size_t counter = 0; counter < number_of_occupied_differences; counter++) { // counts the number of occupied differences we have already encountered + size_t position = __builtin_ctzl(occupied_differences); // count trailing zeros + positions[counter] = position; + + occupied_differences ^= occupied_differences & -occupied_differences; // annihilate the least significant set bit + } + + return positions; +} + + } // namespace GQCP diff --git a/tests/FockSpace/SelectedFockSpace_test.cpp b/tests/FockSpace/SelectedFockSpace_test.cpp new file mode 100644 index 0000000000..6783b5b96b --- /dev/null +++ b/tests/FockSpace/SelectedFockSpace_test.cpp @@ -0,0 +1,55 @@ +// This file is part of GQCG-gqcp. +// +// Copyright (C) 2017-2018 the GQCG developers +// +// GQCG-gqcp is free software: you can redistribute it and/or modify +// it under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// GQCG-gqcp is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with GQCG-gqcp. If not, see . +// +#define BOOST_TEST_MODULE "SelectedFockSpace" + + +#include "FockSpace/SelectedFockSpace.hpp" + + +#include +#include // include this to get main(), otherwise the compiler will complain + + +BOOST_AUTO_TEST_CASE ( constructor ) { + + BOOST_CHECK_NO_THROW(GQCP::SelectedFockSpace (10, 5, 5)); +} + + +BOOST_AUTO_TEST_CASE ( addConfiguration ) { + + // Create a faulty expansion: one of the orbitals is different + GQCP::SelectedFockSpace fock_space (3, 1, 1); + + std::vector alpha_set = {"001", "010"}; + std::vector beta_set = {"001", "010"}; + + BOOST_CHECK_NO_THROW(fock_space.addConfiguration(alpha_set, beta_set)); + + // Test throw with one of the sets is not the same size + std::vector beta_set_long = {"001", "010", "100"}; + BOOST_CHECK_THROW(fock_space.addConfiguration(alpha_set, beta_set_long), std::invalid_argument); + + // Test throw with incompatible orbital numbers + BOOST_CHECK_THROW(fock_space.addConfiguration("0001", "0100"), std::invalid_argument); + + // Test throw with incompatible electron numbers + BOOST_CHECK_THROW(fock_space.addConfiguration("011", "011"), std::invalid_argument); +} + + diff --git a/tests/ONV_test.cpp b/tests/ONV_test.cpp index 2f48666ac3..29d543baf0 100644 --- a/tests/ONV_test.cpp +++ b/tests/ONV_test.cpp @@ -274,3 +274,36 @@ BOOST_AUTO_TEST_CASE ( ONV_address_setNext_fullspace ) { // Checks if no unexpected results occured in a full iteration BOOST_CHECK(is_correct); } + + +BOOST_AUTO_TEST_CASE ( countNumberOfDifferences_unsigned ) { + + GQCP::ONV spin_string1 (5, 3, 21); // "10101" (21) + GQCP::ONV spin_string2 (5, 3, 22); // "10110" (22) + GQCP::ONV spin_string3 (5, 3, 26); // "11010" (26) + + BOOST_CHECK_EQUAL(spin_string1.countNumberOfDifferences(spin_string1), 0); + BOOST_CHECK_EQUAL(spin_string2.countNumberOfDifferences(spin_string2), 0); + BOOST_CHECK_EQUAL(spin_string3.countNumberOfDifferences(spin_string3), 0); + + BOOST_CHECK_EQUAL(spin_string1.countNumberOfDifferences(spin_string2), 2); + BOOST_CHECK_EQUAL(spin_string1.countNumberOfDifferences(spin_string3), 4); + BOOST_CHECK_EQUAL(spin_string2.countNumberOfDifferences(spin_string3), 2); +} + + +BOOST_AUTO_TEST_CASE ( findOccupiedDifferences_unsigned ) { + + GQCP::ONV spin_string1 (5, 3, 21); // "10101" (21) + GQCP::ONV spin_string2 (5, 3, 22); // "10110" (22) + GQCP::ONV spin_string3 (5, 3, 26); // "11010" (26) + + BOOST_TEST(spin_string1.findOccupiedDifferences(spin_string2) == (std::vector {0}), boost::test_tools::per_element()); + BOOST_TEST(spin_string2.findOccupiedDifferences(spin_string1) == (std::vector {1}), boost::test_tools::per_element()); + + BOOST_TEST(spin_string1.findOccupiedDifferences(spin_string3) == (std::vector {0, 2}), boost::test_tools::per_element()); + BOOST_TEST(spin_string3.findOccupiedDifferences(spin_string1) == (std::vector {1, 3}), boost::test_tools::per_element()); + + BOOST_TEST(spin_string2.findOccupiedDifferences(spin_string3) == (std::vector {2}), boost::test_tools::per_element()); + BOOST_TEST(spin_string3.findOccupiedDifferences(spin_string2) == (std::vector {3}), boost::test_tools::per_element()); +} \ No newline at end of file From 462e1742168dc67951df276c0a07a8334299ccc3 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Sun, 21 Oct 2018 21:49:05 +0200 Subject: [PATCH 02/13] fixed namespace, added new functionality --- cmake/SetCMakeEnvironment.cmake | 8 ++++---- include/CISolver/CISolver.hpp | 2 +- include/DOCINewtonOrbitalOptimizer.hpp | 2 +- include/FockSpace/SelectedFockSpace.hpp | 2 +- include/{ => WaveFunction}/WaveFunction.hpp | 0 include/WaveFunction/WaveFunctionReader.hpp | 22 +++++++++++++++++++++ src/{ => WaveFunction}/WaveFunction.cpp | 2 +- 7 files changed, 30 insertions(+), 8 deletions(-) rename include/{ => WaveFunction}/WaveFunction.hpp (100%) create mode 100644 include/WaveFunction/WaveFunctionReader.hpp rename src/{ => WaveFunction}/WaveFunction.cpp (96%) diff --git a/cmake/SetCMakeEnvironment.cmake b/cmake/SetCMakeEnvironment.cmake index 9ceaf1faaf..d73953a65e 100755 --- a/cmake/SetCMakeEnvironment.cmake +++ b/cmake/SetCMakeEnvironment.cmake @@ -53,6 +53,7 @@ set(PROJECT_SOURCE_FILES ${PROJECT_SOURCE_FOLDER}/RHF/PlainRHFSCFSolver.cpp ${PROJECT_SOURCE_FOLDER}/RHF/RHF.cpp ${PROJECT_SOURCE_FOLDER}/RHF/RHFSCFSolver.cpp + ${PROJECT_SOURCE_FOLDER}/WaveFunction/WaveFunction.cpp ${PROJECT_SOURCE_FOLDER}/AOBasis.cpp ${PROJECT_SOURCE_FOLDER}/Atom.cpp ${PROJECT_SOURCE_FOLDER}/DOCINewtonOrbitalOptimizer.cpp @@ -62,8 +63,7 @@ set(PROJECT_SOURCE_FILES ${PROJECT_SOURCE_FOLDER}/miscellaneous.cpp ${PROJECT_SOURCE_FOLDER}/Molecule.cpp ${PROJECT_SOURCE_FOLDER}/ONV.cpp - ${PROJECT_SOURCE_FOLDER}/RMP2.cpp - ${PROJECT_SOURCE_FOLDER}/WaveFunction.cpp) + ${PROJECT_SOURCE_FOLDER}/RMP2.cpp) # Find the header folder set(PROJECT_INCLUDE_FOLDER ${CMAKE_SOURCE_DIR}/include) @@ -101,6 +101,7 @@ set(PROJECT_INCLUDE_FILES ${PROJECT_INCLUDE_FOLDER}/RHF/PlainRHFSCFSolver.hpp ${PROJECT_INCLUDE_FOLDER}/RHF/RHF.hpp ${PROJECT_INCLUDE_FOLDER}/RHF/RHFSCFSolver.hpp + ${PROJECT_INCLUDE_FOLDER}/WaveFunction/WaveFunction.hpp ${PROJECT_INCLUDE_FOLDER}/AOBasis.hpp ${PROJECT_INCLUDE_FOLDER}/Atom.hpp ${PROJECT_INCLUDE_FOLDER}/common.hpp @@ -113,8 +114,7 @@ set(PROJECT_INCLUDE_FILES ${PROJECT_INCLUDE_FOLDER}/Molecule.hpp ${PROJECT_INCLUDE_FOLDER}/ONV.hpp ${PROJECT_INCLUDE_FOLDER}/RMP2.hpp - ${PROJECT_INCLUDE_FOLDER}/units.hpp - ${PROJECT_INCLUDE_FOLDER}/WaveFunction.hpp) + ${PROJECT_INCLUDE_FOLDER}/units.hpp) # Find the tests folder set(PROJECT_TESTS_FOLDER ${CMAKE_SOURCE_DIR}/tests) diff --git a/include/CISolver/CISolver.hpp b/include/CISolver/CISolver.hpp index ac9fcf9f8e..3a5fa87e41 100644 --- a/include/CISolver/CISolver.hpp +++ b/include/CISolver/CISolver.hpp @@ -21,7 +21,7 @@ #include "HamiltonianBuilder/HamiltonianBuilder.hpp" #include "HamiltonianParameters/HamiltonianParameters.hpp" -#include "WaveFunction.hpp" +#include "WaveFunction/WaveFunction.hpp" #include diff --git a/include/DOCINewtonOrbitalOptimizer.hpp b/include/DOCINewtonOrbitalOptimizer.hpp index cd05414ec7..ed9c85a635 100644 --- a/include/DOCINewtonOrbitalOptimizer.hpp +++ b/include/DOCINewtonOrbitalOptimizer.hpp @@ -4,7 +4,7 @@ #include "HamiltonianParameters/HamiltonianParameters.hpp" #include "HamiltonianBuilder/DOCI.hpp" #include "OrbitalOptimizationOptions.hpp" -#include "WaveFunction.hpp" +#include "WaveFunction/WaveFunction.hpp" #include diff --git a/include/FockSpace/SelectedFockSpace.hpp b/include/FockSpace/SelectedFockSpace.hpp index 9b4016134b..1da43cff43 100644 --- a/include/FockSpace/SelectedFockSpace.hpp +++ b/include/FockSpace/SelectedFockSpace.hpp @@ -34,7 +34,7 @@ namespace GQCP { * Where the considered configurations are manually selected and represented as an ONV pair : * a combination of two ONVs, one holding the alpha configuration and holding the bea configuration. */ -class SelectedFockSpace : public GQCG::BaseFockSpace { +class SelectedFockSpace : public GQCP::BaseFockSpace { private: const size_t N_alpha; // number of alpha electrons const size_t N_beta; // number of beta electrons diff --git a/include/WaveFunction.hpp b/include/WaveFunction/WaveFunction.hpp similarity index 100% rename from include/WaveFunction.hpp rename to include/WaveFunction/WaveFunction.hpp diff --git a/include/WaveFunction/WaveFunctionReader.hpp b/include/WaveFunction/WaveFunctionReader.hpp new file mode 100644 index 0000000000..18c79a7f4c --- /dev/null +++ b/include/WaveFunction/WaveFunctionReader.hpp @@ -0,0 +1,22 @@ +#ifndef GQCP_WAVEFUNCTIONREADER_HPP +#define GQCP_WAVEFUNCTIONREADER_HPP + + + + + + + + + + + + + + + + + + + +#endif //GQCP_WAVEFUNCTIONREADER_HPP diff --git a/src/WaveFunction.cpp b/src/WaveFunction/WaveFunction.cpp similarity index 96% rename from src/WaveFunction.cpp rename to src/WaveFunction/WaveFunction.cpp index f8d1ceb0e1..1079c562c9 100644 --- a/src/WaveFunction.cpp +++ b/src/WaveFunction/WaveFunction.cpp @@ -15,7 +15,7 @@ // You should have received a copy of the GNU Lesser General Public License // along with GQCG-gqcp. If not, see . // -#include "WaveFunction.hpp" +#include "WaveFunction/WaveFunction.hpp" namespace GQCP { From 0227d659264fbb45fc081f3bfdef8d78676457d5 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Sun, 21 Oct 2018 22:49:53 +0200 Subject: [PATCH 03/13] fixed pathnames, added selectedReader --- cmake/SetCMakeEnvironment.cmake | 1 + include/FockSpace/BaseFockSpace.hpp | 4 +- include/FockSpace/SelectedFockSpace.hpp | 12 +- include/WaveFunction/WaveFunction.hpp | 2 +- include/WaveFunction/WaveFunctionReader.hpp | 25 +++- src/ONV.cpp | 2 +- src/WaveFunction/WaveFunctionReader.cpp | 135 ++++++++++++++++++++ 7 files changed, 169 insertions(+), 12 deletions(-) create mode 100644 src/WaveFunction/WaveFunctionReader.cpp diff --git a/cmake/SetCMakeEnvironment.cmake b/cmake/SetCMakeEnvironment.cmake index d73953a65e..2cb40c2718 100755 --- a/cmake/SetCMakeEnvironment.cmake +++ b/cmake/SetCMakeEnvironment.cmake @@ -54,6 +54,7 @@ set(PROJECT_SOURCE_FILES ${PROJECT_SOURCE_FOLDER}/RHF/RHF.cpp ${PROJECT_SOURCE_FOLDER}/RHF/RHFSCFSolver.cpp ${PROJECT_SOURCE_FOLDER}/WaveFunction/WaveFunction.cpp + ${PROJECT_SOURCE_FOLDER}/WaveFunction/WaveFunctionReader.cpp ${PROJECT_SOURCE_FOLDER}/AOBasis.cpp ${PROJECT_SOURCE_FOLDER}/Atom.cpp ${PROJECT_SOURCE_FOLDER}/DOCINewtonOrbitalOptimizer.cpp diff --git a/include/FockSpace/BaseFockSpace.hpp b/include/FockSpace/BaseFockSpace.hpp index 5350d3f4e2..7fe5ee586a 100644 --- a/include/FockSpace/BaseFockSpace.hpp +++ b/include/FockSpace/BaseFockSpace.hpp @@ -36,7 +36,7 @@ namespace GQCP { */ class BaseFockSpace { protected: - const size_t K; // number of spatial orbitals + size_t K; // number of spatial orbitals size_t dim; // dimension of the Fock space @@ -45,7 +45,7 @@ class BaseFockSpace { * Protected constructor given a @param K and @param dim */ explicit BaseFockSpace(size_t K, size_t dim); - + BaseFockSpace() = default; public: // DESTRUCTOR diff --git a/include/FockSpace/SelectedFockSpace.hpp b/include/FockSpace/SelectedFockSpace.hpp index 1da43cff43..9ab71c4510 100644 --- a/include/FockSpace/SelectedFockSpace.hpp +++ b/include/FockSpace/SelectedFockSpace.hpp @@ -31,13 +31,13 @@ namespace GQCP { /** * A Fock space for a given set of orbitals and number of alpha and beta electrons. - * Where the considered configurations are manually selected and represented as an ONV pair : - * a combination of two ONVs, one holding the alpha configuration and holding the bea configuration. + * Where the considered configurations are manually selected and represented as a Configuration : + * a combination of two ONVs, one holding the alpha configuration and holding the beta configuration. */ class SelectedFockSpace : public GQCP::BaseFockSpace { private: - const size_t N_alpha; // number of alpha electrons - const size_t N_beta; // number of beta electrons + size_t N_alpha; // number of alpha electrons + size_t N_beta; // number of beta electrons std::vector selection; // the selected configurations @@ -50,10 +50,10 @@ class SelectedFockSpace : public GQCP::BaseFockSpace { // CONSTRUCTORS /** * Constructor given a @param K (spatial orbitals), N_alpha and N_beta (electrons); - * the initial dimension of the space is 0 as no selections are made. + * the initial dimension of the space is 0, as no selections are made. */ SelectedFockSpace(size_t K, size_t N_alpha, size_t N_beta); - + SelectedFockSpace() = default; // DESTRUCTORS ~SelectedFockSpace() override = default; diff --git a/include/WaveFunction/WaveFunction.hpp b/include/WaveFunction/WaveFunction.hpp index c32eb099ad..c318c4a483 100644 --- a/include/WaveFunction/WaveFunction.hpp +++ b/include/WaveFunction/WaveFunction.hpp @@ -39,7 +39,7 @@ class WaveFunction { public: // CONSTRUCTORS WaveFunction(BaseFockSpace& base_fock_space, const Eigen::VectorXd& coefficients); - + WaveFunction() = default; // GETTERS Eigen::VectorXd get_coefficients() const { return coefficients; } diff --git a/include/WaveFunction/WaveFunctionReader.hpp b/include/WaveFunction/WaveFunctionReader.hpp index 18c79a7f4c..87aab2de2e 100644 --- a/include/WaveFunction/WaveFunctionReader.hpp +++ b/include/WaveFunction/WaveFunctionReader.hpp @@ -2,20 +2,41 @@ #define GQCP_WAVEFUNCTIONREADER_HPP +#include "FockSpace/SelectedFockSpace.hpp" +#include "WaveFunction/WaveFunction.hpp" +#include +namespace GQCP { +/** + * Class that reads and stores a selected wavefunction expansion + */ +class WaveFunctionReader { +private: + GQCP::SelectedFockSpace fock_space; + Eigen::VectorXd coefficients; + GQCP::WaveFunction wave_function; +public: + /** + * Constructor based on a given @param GAMESS_filename + */ + explicit WaveFunctionReader(const std::string& GAMESS_filename); + // GETTERS + SelectedFockSpace get_fock_space() const { return this->fock_space; } + Eigen::VectorXd get_coefficients() const { return this->coefficients; } + WaveFunction get_wave_function() const { return this->wave_function; } +}; - - +} // namespace GQCP diff --git a/src/ONV.cpp b/src/ONV.cpp index 0daa34a38a..6c3d5a9fc6 100644 --- a/src/ONV.cpp +++ b/src/ONV.cpp @@ -15,7 +15,7 @@ // You should have received a copy of the GNU Lesser General Public License // along with GQCG-gqcp. If not, see . // -#include "../../gqcg/include/ONV.hpp" +#include "ONV.hpp" #include diff --git a/src/WaveFunction/WaveFunctionReader.cpp b/src/WaveFunction/WaveFunctionReader.cpp new file mode 100644 index 0000000000..7661f482fe --- /dev/null +++ b/src/WaveFunction/WaveFunctionReader.cpp @@ -0,0 +1,135 @@ +// This file is part of GQCG-gqcp. +// +// Copyright (C) 2017-2018 the GQCG developers +// +// GQCG-gqcp is free software: you can redistribute it and/or modify +// it under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// GQCG-gqcp is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with GQCG-gqcp. If not, see . +// + +#include "WaveFunction/WaveFunctionReader.hpp" +#include + +namespace GQCP { + + +/* + * CONSTRUCTORS + */ + +/** + * Constructor based on a given @param GAMESS_filename + */ +WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) +{ + // If the filename isn't properly converted into an input file stream, we assume the user supplied a wrong file + std::ifstream input_file_stream(GAMESS_filename); + std::ifstream input_file_stream_count(GAMESS_filename); // made to count the expansion size + + if (!input_file_stream.good()) { + throw std::runtime_error("WaveFunctionReader(): The provided GAMESS file is illegible. Maybe you specified a wrong path?"); + } + + // Do the actual parsing + // Read in dummy lines up until we actually get to the ONVs and coefficients + std::string line; + std::string buffer; // dummy for the getline() TODO: find "correcter" way if possible + while (std::getline(input_file_stream, line)) { + std::getline(input_file_stream_count, buffer); + if (line.find("ALPHA") != std::string::npos + && line.find("BETA") != std::string::npos + && line.find("COEFFICIENT") != std::string::npos) { // if find returns an index that's different from the 'not-found' index + + // this line should have dashes and we skip it + std::getline(input_file_stream, line); + std::getline(input_file_stream_count, buffer); + break; + } + } + + size_t space_size = 0; + // Count the amount of configurations + while (std::getline(input_file_stream, buffer)) { + if (buffer.length() != 0 | buffer != '\n') { + space_size++; + } + } + + this->coefficients = Eigen::VectorXd::Zero(space_size); + + + std::getline(input_file_stream, line); + // Read the first line containing the configurations + std::vector splitted_line; + boost::split(splitted_line, line, boost::is_any_of("|")); // split on '|' + + + // Create an alpha SpinString for the first field + std::string trimmed_alpha = boost::algorithm::trim_copy(splitted_line[0]); + std::string reversed_alpha (trimmed_alpha.rbegin(), trimmed_alpha.rend()); + + // Create a beta SpinString for the second field + std::string trimmed_beta = boost::algorithm::trim_copy(splitted_line[1]); + std::string reversed_beta(trimmed_beta.rbegin(), trimmed_beta.rend()); + + size_t counter = 0; + // add the double + coefficient(counter) = std::stod(splitted_line[2]); + + // dynamic bitset provides us with functionallity + boost::dynamic_bitset<> alpha_transfer (trimmed_alpha); + boost::dynamic_bitset<> beta_transfer (trimmed_beta); + + size_t K = alpha_transfer.size(); + size_t N_alpha = alpha_transfer.count(); + size_t N_beta = beta_transfer.count(); + + fock_space = SelectedFockSpace(K, N_alpha, N_beta); + + fock_space.addConfiguration(trimmed_alpha, trimmed_beta); + + + // Read in the ONVs and the coefficients by splitting the line on '|', and then trimming whitespace + // In the GAMESS format, the bit strings are ordered in reverse + while (std::getline(input_file_stream, line)) { + + counter++; + std::vector splitted_line; + boost::split(splitted_line, line, boost::is_any_of("|")); // split on '|' + + + // Create an alpha SpinString for the first field + std::string trimmed_alpha = boost::algorithm::trim_copy(splitted_line[0]); + if (trimmed_alpha.length() != K) { + throw std::invalid_argument("WaveFunctionReader(): One of the provided alpha ONVs does not have the correct number of orbitals."); + } + std::string reversed_alpha (trimmed_alpha.rbegin(), trimmed_alpha.rend()); + + // Create a beta SpinString for the second field + std::string trimmed_beta = boost::algorithm::trim_copy(splitted_line[1]); + if (trimmed_beta.length() != K) { + throw std::invalid_argument("WaveFunctionReader(): One of the provided beta ONVs does not have the correct number of orbitals."); + } + std::string reversed_beta(trimmed_beta.rbegin(), trimmed_beta.rend()); + + + // Create a double for the third field + coefficient(counter) = std::stod(splitted_line[2]); + fock_space.addConfiguration(trimmed_alpha, trimmed_beta); + + } // while getline + + this->wave_function = WaveFunction(fock_space, coefficients); +} + + +} // namespace GQCP From bd7754a28cb2b7e3f0fb6716e951ccfc2003a544 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Sun, 21 Oct 2018 22:51:20 +0200 Subject: [PATCH 04/13] fixed header, and typos --- src/WaveFunction/WaveFunctionReader.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/WaveFunction/WaveFunctionReader.cpp b/src/WaveFunction/WaveFunctionReader.cpp index 7661f482fe..39267d06f7 100644 --- a/src/WaveFunction/WaveFunctionReader.cpp +++ b/src/WaveFunction/WaveFunctionReader.cpp @@ -17,6 +17,10 @@ // #include "WaveFunction/WaveFunctionReader.hpp" + +#include +#include +#include "boost/dynamic_bitset.hpp" #include namespace GQCP { @@ -59,7 +63,7 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) size_t space_size = 0; // Count the amount of configurations while (std::getline(input_file_stream, buffer)) { - if (buffer.length() != 0 | buffer != '\n') { + if (buffer.length() != 0 | buffer[0] != '\n') { space_size++; } } @@ -83,7 +87,7 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) size_t counter = 0; // add the double - coefficient(counter) = std::stod(splitted_line[2]); + this->coefficients(counter) = std::stod(splitted_line[2]); // dynamic bitset provides us with functionallity boost::dynamic_bitset<> alpha_transfer (trimmed_alpha); @@ -93,9 +97,9 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) size_t N_alpha = alpha_transfer.count(); size_t N_beta = beta_transfer.count(); - fock_space = SelectedFockSpace(K, N_alpha, N_beta); + this->fock_space = SelectedFockSpace(K, N_alpha, N_beta); - fock_space.addConfiguration(trimmed_alpha, trimmed_beta); + this->fock_space.addConfiguration(trimmed_alpha, trimmed_beta); // Read in the ONVs and the coefficients by splitting the line on '|', and then trimming whitespace @@ -123,12 +127,12 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) // Create a double for the third field - coefficient(counter) = std::stod(splitted_line[2]); - fock_space.addConfiguration(trimmed_alpha, trimmed_beta); + this->coefficients(counter) = std::stod(splitted_line[2]); + this->fock_space.addConfiguration(trimmed_alpha, trimmed_beta); } // while getline - this->wave_function = WaveFunction(fock_space, coefficients); + this->wave_function = WaveFunction(this->fock_space, this->coefficients); } From 3392bbde2cf14812064fd09b618df3cfad937daa Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Sun, 21 Oct 2018 23:10:44 +0200 Subject: [PATCH 05/13] boost/split, libint error --- src/WaveFunction/WaveFunctionReader.cpp | 10 +++++++--- tests/FockSpace/SelectedFockSpace_test.cpp | 10 ++++++++++ tests/data/test_GAMESS_expansion | 6 ++++++ 3 files changed, 23 insertions(+), 3 deletions(-) create mode 100644 tests/data/test_GAMESS_expansion diff --git a/src/WaveFunction/WaveFunctionReader.cpp b/src/WaveFunction/WaveFunctionReader.cpp index 39267d06f7..f94ec227a5 100644 --- a/src/WaveFunction/WaveFunctionReader.cpp +++ b/src/WaveFunction/WaveFunctionReader.cpp @@ -18,11 +18,15 @@ #include "WaveFunction/WaveFunctionReader.hpp" -#include -#include -#include "boost/dynamic_bitset.hpp" #include +#include +#include + + +#include + + namespace GQCP { diff --git a/tests/FockSpace/SelectedFockSpace_test.cpp b/tests/FockSpace/SelectedFockSpace_test.cpp index 6783b5b96b..2add9222e4 100644 --- a/tests/FockSpace/SelectedFockSpace_test.cpp +++ b/tests/FockSpace/SelectedFockSpace_test.cpp @@ -20,6 +20,8 @@ #include "FockSpace/SelectedFockSpace.hpp" +#include "WaveFunction/WaveFunctionReader.hpp" + #include #include // include this to get main(), otherwise the compiler will complain @@ -53,3 +55,11 @@ BOOST_AUTO_TEST_CASE ( addConfiguration ) { } +BOOST_AUTO_TEST_CASE ( reader_test ) { + + GQCP::WaveFunctionReader test_reader ("../tests/data/test_GAMESS_expansion"); + Eigen::VectorXd test = {1.0000, 0.0000}; + + // Check if both expansions are considered equal + BOOST_CHECK(test.isApprox(test_reader.get_coefficients())); +} \ No newline at end of file diff --git a/tests/data/test_GAMESS_expansion b/tests/data/test_GAMESS_expansion new file mode 100644 index 0000000000..d133e7a45f --- /dev/null +++ b/tests/data/test_GAMESS_expansion @@ -0,0 +1,6 @@ + STATE 1 ENERGY= -1.0502 + + ALPHA | BETA | COEFFICIENT +----------------------------------------------|----------------------------------------------|--------------- + 1000000000000000000000000000000000000000000000 | 1000000000000000000000000000000000000000000000 | 1.0000000 + 1000000000000000000000000000000000000000000000 | 0100000000000000000000000000000000000000000000 | 0.0000000 From fffab31eb327fbbc5784cb6bc5a5ad6b668e00a3 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Mon, 22 Oct 2018 10:50:23 +0200 Subject: [PATCH 06/13] added onv default --- include/ONV.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/include/ONV.hpp b/include/ONV.hpp index 5dd040abfc..0f13d5c0b8 100644 --- a/include/ONV.hpp +++ b/include/ONV.hpp @@ -40,8 +40,8 @@ namespace GQCP { */ class ONV { private: - const size_t K; // number of spatial orbitals - const size_t N; // number of electrons + size_t K; // number of spatial orbitals + size_t N; // number of electrons size_t unsigned_representation; // unsigned representation VectorXs occupation_indices; // the occupied orbital electron indexes // it is a vector of N elements in which occupation_indices[j] @@ -50,6 +50,7 @@ class ONV { public: // CONSTRUCTORS + ONV() = default; /** * Constructor from a @param K orbitals, N electrons and an @param unsigned_representation */ From 4f14447905187b5e1772f9249cd9f7e3efe94f73 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Mon, 22 Oct 2018 12:14:15 +0200 Subject: [PATCH 07/13] disabled WaveFunctionReader temp --- cmake/FindPackages.cmake | 2 +- include/CISolver/CISolver.hpp | 1 + include/WaveFunction/WaveFunctionReader.hpp | 4 +++- src/WaveFunction/WaveFunctionReader.cpp | 10 +++++----- tests/FockSpace/SelectedFockSpace_test.cpp | 7 ++++--- 5 files changed, 14 insertions(+), 10 deletions(-) diff --git a/cmake/FindPackages.cmake b/cmake/FindPackages.cmake index d69940454c..c4bf25fd0a 100755 --- a/cmake/FindPackages.cmake +++ b/cmake/FindPackages.cmake @@ -2,7 +2,7 @@ # Find the Boost package -find_package(Boost REQUIRED) +find_package(Boost REQUIRED COMPONENTS system thread) # Find Eigen3 find_package(Eigen3 3.3.4 REQUIRED) diff --git a/include/CISolver/CISolver.hpp b/include/CISolver/CISolver.hpp index 3a5fa87e41..b56776e242 100644 --- a/include/CISolver/CISolver.hpp +++ b/include/CISolver/CISolver.hpp @@ -19,6 +19,7 @@ #define GQCP_CISOLVER_HPP + #include "HamiltonianBuilder/HamiltonianBuilder.hpp" #include "HamiltonianParameters/HamiltonianParameters.hpp" #include "WaveFunction/WaveFunction.hpp" diff --git a/include/WaveFunction/WaveFunctionReader.hpp b/include/WaveFunction/WaveFunctionReader.hpp index 87aab2de2e..d5047b5b37 100644 --- a/include/WaveFunction/WaveFunctionReader.hpp +++ b/include/WaveFunction/WaveFunctionReader.hpp @@ -2,10 +2,12 @@ #define GQCP_WAVEFUNCTIONREADER_HPP + +#include + #include "FockSpace/SelectedFockSpace.hpp" #include "WaveFunction/WaveFunction.hpp" - #include diff --git a/src/WaveFunction/WaveFunctionReader.cpp b/src/WaveFunction/WaveFunctionReader.cpp index f94ec227a5..ffee161ee5 100644 --- a/src/WaveFunction/WaveFunctionReader.cpp +++ b/src/WaveFunction/WaveFunctionReader.cpp @@ -16,15 +16,13 @@ // along with GQCG-gqcp. If not, see . // -#include "WaveFunction/WaveFunctionReader.hpp" - #include -#include -#include -#include +#include "WaveFunction/WaveFunctionReader.hpp" + + namespace GQCP { @@ -39,6 +37,7 @@ namespace GQCP { */ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) { + /* // If the filename isn't properly converted into an input file stream, we assume the user supplied a wrong file std::ifstream input_file_stream(GAMESS_filename); std::ifstream input_file_stream_count(GAMESS_filename); // made to count the expansion size @@ -137,6 +136,7 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) } // while getline this->wave_function = WaveFunction(this->fock_space, this->coefficients); + */ } diff --git a/tests/FockSpace/SelectedFockSpace_test.cpp b/tests/FockSpace/SelectedFockSpace_test.cpp index 2add9222e4..fb39c9d0d3 100644 --- a/tests/FockSpace/SelectedFockSpace_test.cpp +++ b/tests/FockSpace/SelectedFockSpace_test.cpp @@ -20,7 +20,7 @@ #include "FockSpace/SelectedFockSpace.hpp" -#include "WaveFunction/WaveFunctionReader.hpp" +//#include "WaveFunction/WaveFunctionReader.hpp" #include @@ -54,7 +54,7 @@ BOOST_AUTO_TEST_CASE ( addConfiguration ) { BOOST_CHECK_THROW(fock_space.addConfiguration("011", "011"), std::invalid_argument); } - +/* BOOST_AUTO_TEST_CASE ( reader_test ) { GQCP::WaveFunctionReader test_reader ("../tests/data/test_GAMESS_expansion"); @@ -62,4 +62,5 @@ BOOST_AUTO_TEST_CASE ( reader_test ) { // Check if both expansions are considered equal BOOST_CHECK(test.isApprox(test_reader.get_coefficients())); -} \ No newline at end of file +} + */ \ No newline at end of file From b50c18df937ff2ad25e7ea54b489ec7d2978e725 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Mon, 22 Oct 2018 16:08:25 +0200 Subject: [PATCH 08/13] added tests for the selected RDMs --- cmake/SetCMakeEnvironment.cmake | 4 + include/FockSpace/SelectedFockSpace.hpp | 6 + include/RDM/OneRDM.hpp | 3 + include/RDM/RDMCalculator.hpp | 4 +- include/RDM/SelectedRDMBuilder.hpp | 68 +++++ src/FockSpace/SelectedFockSpace.cpp | 38 +++ src/RDM/RDMCalculator.cpp | 15 + src/RDM/SelectedRDMBuilder.cpp | 360 ++++++++++++++++++++++++ tests/RDM/SelectedRDMBuilder_test.cpp | 129 +++++++++ 9 files changed, 626 insertions(+), 1 deletion(-) create mode 100644 include/RDM/SelectedRDMBuilder.hpp create mode 100644 src/RDM/SelectedRDMBuilder.cpp create mode 100644 tests/RDM/SelectedRDMBuilder_test.cpp diff --git a/cmake/SetCMakeEnvironment.cmake b/cmake/SetCMakeEnvironment.cmake index 2cb40c2718..0566aae695 100755 --- a/cmake/SetCMakeEnvironment.cmake +++ b/cmake/SetCMakeEnvironment.cmake @@ -46,8 +46,10 @@ set(PROJECT_SOURCE_FILES ${PROJECT_SOURCE_FOLDER}/RDM/DOCIRDMBuilder.cpp ${PROJECT_SOURCE_FOLDER}/RDM/FCIRDMBuilder.cpp ${PROJECT_SOURCE_FOLDER}/RDM/OneRDM.cpp + ${PROJECT_SOURCE_FOLDER}/RDM/SelectedRDMBuilder.cpp ${PROJECT_SOURCE_FOLDER}/RDM/RDMCalculator.cpp ${PROJECT_SOURCE_FOLDER}/RDM/RDMs.cpp + ${PROJECT_SOURCE_FOLDER}/RDM/SelectedRDMBuilder.cpp ${PROJECT_SOURCE_FOLDER}/RDM/TwoRDM.cpp ${PROJECT_SOURCE_FOLDER}/RHF/DIISRHFSCFSolver.cpp ${PROJECT_SOURCE_FOLDER}/RHF/PlainRHFSCFSolver.cpp @@ -95,6 +97,7 @@ set(PROJECT_INCLUDE_FILES ${PROJECT_INCLUDE_FOLDER}/RDM/DOCIRDMBuilder.hpp ${PROJECT_INCLUDE_FOLDER}/RDM/FCIRDMBuilder.hpp ${PROJECT_INCLUDE_FOLDER}/RDM/OneRDM.hpp + ${PROJECT_INCLUDE_FOLDER}/RDM/SelectedRDMBuilder.hpp ${PROJECT_INCLUDE_FOLDER}/RDM/RDMCalculator.hpp ${PROJECT_INCLUDE_FOLDER}/RDM/RDMs.hpp ${PROJECT_INCLUDE_FOLDER}/RDM/TwoRDM.hpp @@ -142,6 +145,7 @@ set(PROJECT_TEST_SOURCE_FILES ${PROJECT_TESTS_FOLDER}/Operator/TwoElectronOperator_test.cpp ${PROJECT_TESTS_FOLDER}/RDM/DOCIRDMBuilder_test.cpp ${PROJECT_TESTS_FOLDER}/RDM/FCIRDMBuilder_test.cpp + ${PROJECT_TESTS_FOLDER}/RDM/SelectedRDMBuilder_test.cpp ${PROJECT_TESTS_FOLDER}/RDM/RDMCalculator_test.cpp ${PROJECT_TESTS_FOLDER}/RHF/DIISRHFSCFSolver_test.cpp ${PROJECT_TESTS_FOLDER}/RHF/PlainRHFSCFSolver_test.cpp diff --git a/include/FockSpace/SelectedFockSpace.hpp b/include/FockSpace/SelectedFockSpace.hpp index 9ab71c4510..e78fd670ea 100644 --- a/include/FockSpace/SelectedFockSpace.hpp +++ b/include/FockSpace/SelectedFockSpace.hpp @@ -20,6 +20,7 @@ #include "FockSpace/BaseFockSpace.hpp" +#include "FockSpace/FockSpaceProduct.hpp" #include "Configuration.hpp" #include @@ -73,6 +74,11 @@ class SelectedFockSpace : public GQCP::BaseFockSpace { */ void addConfiguration(std::string onv1, std::string onv2); void addConfiguration(std::vector onv1s, std::vector onv2s); + + /** + * Sets the current expansion to that of @param fock_space + */ + void setExpansion(const FockSpaceProduct& fock_space); }; diff --git a/include/RDM/OneRDM.hpp b/include/RDM/OneRDM.hpp index bf001efcfb..dc036e9376 100644 --- a/include/RDM/OneRDM.hpp +++ b/include/RDM/OneRDM.hpp @@ -43,6 +43,9 @@ class OneRDM : public BaseRDM { Eigen::MatrixXd get_matrix_representation() const { return this->D; } double get(size_t p, size_t q) const { return this->D(p, q); } + // OPERATOR + + // PUBLIC METHODS /** diff --git a/include/RDM/RDMCalculator.hpp b/include/RDM/RDMCalculator.hpp index 9b5b197670..b936b2a621 100644 --- a/include/RDM/RDMCalculator.hpp +++ b/include/RDM/RDMCalculator.hpp @@ -4,6 +4,7 @@ #include "RDM/BaseRDMBuilder.hpp" #include "FockSpace/FockSpace.hpp" #include "FockSpace/FockSpaceProduct.hpp" +#include "FockSpace/SelectedFockSpace.hpp" #include @@ -26,7 +27,8 @@ class RDMCalculator { */ RDMCalculator(const FockSpace& fock_space); // DOCIRDMBuilder RDMCalculator(const FockSpaceProduct& fock_space); // FCIRDMBuilder - RDMCalculator(const BaseFockSpace& fock_space); + RDMCalculator(const SelectedFockSpace& fock_space); // SelectedRDMBuilder + RDMCalculator(const BaseFockSpace& fock_space); // runtime // PUBLIC METHODS diff --git a/include/RDM/SelectedRDMBuilder.hpp b/include/RDM/SelectedRDMBuilder.hpp new file mode 100644 index 0000000000..3c49e49f43 --- /dev/null +++ b/include/RDM/SelectedRDMBuilder.hpp @@ -0,0 +1,68 @@ +// This file is part of GQCG-gqcp. +// +// Copyright (C) 2017-2018 the GQCG developers +// +// GQCG-gqcp is free software: you can redistribute it and/or modify +// it under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// GQCG-gqcp is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with GQCG-gqcp. If not, see . +// +#ifndef GQCP_SELECTEDRDMBUILDER_HPP +#define GQCP_SELECTEDRDMBUILDER_HPP + + +#include "FockSpace/SelectedFockSpace.hpp" +#include "RDM/BaseRDMBuilder.hpp" +#include "RDM/RDMs.hpp" + + +namespace GQCP { + + +/** + * SelectedRDMBuilder is a class for the calculation of a density matrix from a given wave function + * or coefficient expansion in the full CI Fock space + */ +class SelectedRDMBuilder : public BaseRDMBuilder { + SelectedFockSpace fock_space; // Fock space containing the selected configuration + + +public: + // CONSTRUCTOR + explicit SelectedRDMBuilder (const SelectedFockSpace& fock_space); + + + // DESTRUCTOR + ~SelectedRDMBuilder() = default; + + + // OVERRIDDEN PUBLIC METHODS + /** + * @return all 1-RDMs from a coefficient vector @param x + */ + OneRDMs calculate1RDMs(const Eigen::VectorXd& x) override; + + /** + * @return all 2-RDMs from a coefficient vector @param x + */ + TwoRDMs calculate2RDMs(const Eigen::VectorXd& x) override; + + /** + * @return the Fock space of the RDMBuilder + */ + BaseFockSpace* get_fock_space() override { return &fock_space; } +}; + + +} // namespace GQCP + + +#endif // GQCP_SELECTEDRDMBUILDER_HPP diff --git a/src/FockSpace/SelectedFockSpace.cpp b/src/FockSpace/SelectedFockSpace.cpp index b012856653..8b14ce7c1e 100644 --- a/src/FockSpace/SelectedFockSpace.cpp +++ b/src/FockSpace/SelectedFockSpace.cpp @@ -98,5 +98,43 @@ void SelectedFockSpace::addConfiguration(std::vector onv1s, std::ve } } +/** + * Sets the current expansion to that of @param fock_space + */ +void SelectedFockSpace::setExpansion(const FockSpaceProduct& fock_space){ + + std::vector configurations; + + FockSpace fock_space_alpha = fock_space.get_fock_space_alpha(); + FockSpace fock_space_beta = fock_space.get_fock_space_beta(); + + auto dim_alpha = fock_space_alpha.get_dimension(); + auto dim_beta = fock_space_beta.get_dimension(); + + ONV alpha = fock_space_alpha.get_ONV(0); + for (size_t I_alpha = 0; I_alpha < dim_alpha; I_alpha++) { + + ONV beta = fock_space_beta.get_ONV(0); + for (size_t I_beta = 0; I_beta < dim_beta; I_beta++) { + + configurations.push_back(Configuration {alpha, beta}); + + + if (I_beta < dim_beta - 1) { // prevent the last permutation to occur + fock_space_beta.setNext(beta); + } + + } + + if (I_alpha < dim_alpha - 1) { // prevent the last permutation to occur + fock_space_alpha.setNext(alpha); + } + + } + this->dim = fock_space.get_dimension(); + this->selection = configurations; + + +} } // namespace GQCP diff --git a/src/RDM/RDMCalculator.cpp b/src/RDM/RDMCalculator.cpp index c505858fa1..deafa89082 100644 --- a/src/RDM/RDMCalculator.cpp +++ b/src/RDM/RDMCalculator.cpp @@ -19,6 +19,7 @@ #include "RDM/DOCIRDMBuilder.hpp" #include "RDM/FCIRDMBuilder.hpp" +#include "RDM/SelectedRDMBuilder.hpp" @@ -45,6 +46,14 @@ RDMCalculator::RDMCalculator(const FockSpaceProduct& fock_space) { } +/** + * Allocates a SelectedRDMBuilder based on @param fock_space + */ +RDMCalculator::RDMCalculator(const SelectedFockSpace& fock_space) { + rdm_builder = std::make_shared(fock_space); +} + + /** * Allocates the correct derived BaseRDMBuilder based on @param fock_space */ @@ -63,6 +72,12 @@ RDMCalculator::RDMCalculator(const BaseFockSpace& fock_space) { break; } + + case FockSpaceType::SelectedFockSpace: { + rdm_builder = std::make_shared(dynamic_cast(fock_space)); + + break; + } } } diff --git a/src/RDM/SelectedRDMBuilder.cpp b/src/RDM/SelectedRDMBuilder.cpp new file mode 100644 index 0000000000..75cc5f5e38 --- /dev/null +++ b/src/RDM/SelectedRDMBuilder.cpp @@ -0,0 +1,360 @@ +// This file is part of GQCG-gqcp. +// +// Copyright (C) 2017-2018 the GQCG developers +// +// GQCG-gqcp is free software: you can redistribute it and/or modify +// it under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// GQCG-gqcp is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with GQCG-gqcp. If not, see . +// +#include "RDM/SelectedRDMBuilder.hpp" + +#include "Configuration.hpp" + + +namespace GQCP { + + +/* + * CONSTRUCTOR + */ +SelectedRDMBuilder::SelectedRDMBuilder(const SelectedFockSpace& fock_space) : + fock_space (fock_space) +{} + + +/* + * OVERRIDDEN PUBLIC METHODS + */ + +/** + * @return 1RDM from a coefficient vector @param x + */ +OneRDMs SelectedRDMBuilder::calculate1RDMs(const Eigen::VectorXd& x) { + + // Initialize as zero matrices + size_t K = this->fock_space.get_K(); + + Eigen::MatrixXd D_aa = Eigen::MatrixXd::Zero(K, K); + Eigen::MatrixXd D_bb = Eigen::MatrixXd::Zero(K, K); + + + size_t dim = fock_space.get_dimension(); + + for (size_t I = 0; I < dim; I++) { // loop over all addresses (1) + Configuration configuration_I = this->fock_space.get_Configuration(I); + ONV alpha_I = configuration_I.onv_alpha; + ONV beta_I = configuration_I.onv_beta; + + double c_I = x(I); + + + // Calculate the diagonal of the 1-RDMs + for (size_t p = 0; p < K; p++) { + + if (alpha_I.isOccupied(p)) { + D_aa(p,p) += std::pow(c_I, 2); + } + + if (beta_I.isOccupied(p)) { + D_bb(p,p) += std::pow(c_I, 2); + } + } + + + // Calculate the off-diagonal elements, by going over all other ONVs + for (size_t J = I+1; J < dim; J++) { + + Configuration configuration_J = this->fock_space.get_Configuration(J); + ONV alpha_J = configuration_J.onv_alpha; + ONV beta_J = configuration_J.onv_beta; + + double c_J = x(J); + + + // 1 electron excitation in alpha (i.e. 2 differences), 0 in beta + if ((alpha_I.countNumberOfDifferences(alpha_J) == 2) && (beta_I.countNumberOfDifferences(beta_J) == 0)) { + + // Find the orbitals that are occupied in one string, and aren't in the other + size_t p = alpha_I.findOccupiedDifferences(alpha_J)[0]; // we're sure that there is only 1 element in the std::vector + size_t q = alpha_J.findOccupiedDifferences(alpha_I)[0]; // we're sure that there is only 1 element in the std::vector + + // Calculate the total sign, and include it in the RDM contribution + int sign = alpha_I.operatorPhaseFactor(p) * alpha_J.operatorPhaseFactor(q); + D_aa(p,q) += sign * c_I * c_J; + D_aa(q,p) += sign * c_I * c_J; + } + + + // 1 electron excitation in beta, 0 in alpha + if ((alpha_I.countNumberOfDifferences(alpha_J) == 0) && (beta_I.countNumberOfDifferences(beta_J) == 2)) { + + // Find the orbitals that are occupied in one string, and aren't in the other + size_t p = beta_I.findOccupiedDifferences(beta_J)[0]; // we're sure that there is only 1 element in the std::vector + size_t q = beta_J.findOccupiedDifferences(beta_I)[0]; // we're sure that there is only 1 element in the std::vector + + // Calculate the total sign, and include it in the RDM contribution + int sign = beta_I.operatorPhaseFactor(p) * beta_J.operatorPhaseFactor(q); + D_bb(p,q) += sign * c_I * c_J; + D_bb(q,p) += sign * c_I * c_J; + } + + } // loop over addresses J > I + } // loop over addresses I + + OneRDM one_rdm_aa (D_aa); + OneRDM one_rdm_bb (D_bb); + // The total 1-RDM is the sum of the spin components + return OneRDMs (one_rdm_aa, one_rdm_bb); +} + + +/** + * @return 2RDM from a coefficient vector @param x + */ +TwoRDMs SelectedRDMBuilder::calculate2RDMs(const Eigen::VectorXd& x) { + + + // Initialize as zero matrices + size_t K = this->fock_space.get_K(); + + size_t dim = fock_space.get_dimension(); + + Eigen::Tensor d_aaaa (K,K,K,K); + d_aaaa.setZero(); + Eigen::Tensor d_aabb (K,K,K,K); + d_aabb.setZero(); + Eigen::Tensor d_bbaa (K,K,K,K); + d_bbaa.setZero(); + Eigen::Tensor d_bbbb (K,K,K,K); + d_bbbb.setZero(); + + for (size_t I = 0; I < dim; I++) { // loop over all addresses I + + Configuration configuration_I = this->fock_space.get_Configuration(I); + ONV alpha_I = configuration_I.onv_alpha; + ONV beta_I = configuration_I.onv_beta; + + double c_I = x(I); + + for (size_t p = 0; p < K; p++) { + + // 'Diagonal' elements of the 2-RDM: aaaa and aabb + if (alpha_I.isOccupied(p)) { + for (size_t q = 0; q < K; q++) { + if (beta_I.isOccupied(q)) { + d_aabb(p,p,q,q) += std::pow(c_I, 2); + } + + if (p != q) { // can't create/annihilate the same orbital twice + if (alpha_I.isOccupied(q)) { + d_aaaa(p,p,q,q) += std::pow(c_I, 2); + d_aaaa(p,q,q,p) -= std::pow(c_I, 2); + } + } + + } // loop over q + } + + // 'Diagonal' elements of the 2-RDM: bbbb and bbaa + if (beta_I.isOccupied(p)) { + for (size_t q = 0; q < K; q++) { + if (alpha_I.isOccupied(q)) { + d_bbaa(p,p,q,q) += std::pow(c_I, 2); + } + + if (p != q) { // can't create/annihilate the same orbital twice + if (beta_I.isOccupied(q)) { + d_bbbb(p,p,q,q) += std::pow(c_I, 2); + d_bbbb(p,q,q,p) -= std::pow(c_I, 2); + } + } + } // loop over q + } + } // loop over q + + + for (size_t J = I+1; J < dim; J++) { + + Configuration configuration_J = this->fock_space.get_Configuration(J); + ONV alpha_J = configuration_J.onv_alpha; + ONV beta_J = configuration_J.onv_beta; + + double c_J = x(J); + + // 1 electron excitation in alpha, 0 in beta + if ((alpha_I.countNumberOfDifferences(alpha_J) == 2) && (beta_I.countNumberOfDifferences(beta_J) == 0)) { + + // Find the orbitals that are occupied in one string, and aren't in the other + size_t p = alpha_I.findOccupiedDifferences(alpha_J)[0]; // we're sure that there is only 1 element in the std::vector + size_t q = alpha_J.findOccupiedDifferences(alpha_I)[0]; // we're sure that there is only 1 element in the std::vector + + // Calculate the total sign + int sign = alpha_I.operatorPhaseFactor(p) * alpha_J.operatorPhaseFactor(q); + + + for (size_t r = 0; r < K; r++) { // r loops over spatial orbitals + + if (alpha_I.isOccupied(r) && alpha_J.isOccupied(r)) { // r must be occupied on the left and on the right + if ((p != r) && (q != r)) { // can't create or annihilate the same orbital + // Fill in the 2-RDM contributions + d_aaaa(p,q,r,r) += sign * c_I * c_J; + d_aaaa(r,q,p,r) -= sign * c_I * c_J; + d_aaaa(p,r,r,q) -= sign * c_I * c_J; + d_aaaa(r,r,p,q) += sign * c_I * c_J; + + d_aaaa(q,p,r,r) += sign * c_I * c_J; + d_aaaa(q,r,r,p) -= sign * c_I * c_J; + d_aaaa(r,p,q,r) -= sign * c_I * c_J; + d_aaaa(r,r,q,p) += sign * c_I * c_J; + } + } + + if (beta_I.isOccupied(r)) { // beta_I == beta_J from the previous if-branch + + // Fill in the 2-RDM contributions + d_aabb(p,q,r,r) += sign * c_I * c_J; + d_aabb(q,p,r,r) += sign * c_I * c_J; + + d_bbaa(r,r,p,q) += sign * c_I * c_J; + d_bbaa(r,r,q,p) += sign * c_I * c_J; + } + } + } + + + // 0 electron excitations in alpha, 1 in beta + if ((alpha_I.countNumberOfDifferences(alpha_J) == 0) && (beta_I.countNumberOfDifferences(beta_J) == 2)) { + + // Find the orbitals that are occupied in one string, and aren't in the other + size_t p = beta_I.findOccupiedDifferences(beta_J)[0]; // we're sure that there is only 1 element in the std::vector + size_t q = beta_J.findOccupiedDifferences(beta_I)[0]; // we're sure that there is only 1 element in the std::vector + + // Calculate the total sign + int sign = beta_I.operatorPhaseFactor(p) * beta_J.operatorPhaseFactor(q); + + + for (size_t r = 0; r < K; r++) { // r loops over spatial orbitals + + if (beta_I.isOccupied(r) && beta_J.isOccupied(r)) { // r must be occupied on the left and on the right + if ((p != r) && (q != r)) { // can't create or annihilate the same orbital + // Fill in the 2-RDM contributions + d_bbbb(p,q,r,r) += sign * c_I * c_J; + d_bbbb(r,q,p,r) -= sign * c_I * c_J; + d_bbbb(p,r,r,q) -= sign * c_I * c_J; + d_bbbb(r,r,p,q) += sign * c_I * c_J; + + d_bbbb(q,p,r,r) += sign * c_I * c_J; + d_bbbb(q,r,r,p) -= sign * c_I * c_J; + d_bbbb(r,p,q,r) -= sign * c_I * c_J; + d_bbbb(r,r,q,p) += sign * c_I * c_J; + } + } + + if (alpha_I.isOccupied(r)) { // alpha_I == alpha_J from the previous if-branch + + // Fill in the 2-RDM contributions + d_bbaa(p,q,r,r) += sign * c_I * c_J; + d_bbaa(q,p,r,r) += sign * c_I * c_J; + + d_aabb(r,r,p,q) += sign * c_I * c_J; + d_aabb(r,r,q,p) += sign * c_I * c_J; + } + } + } + + + // 1 electron excitation in alpha, 1 in beta + if ((alpha_I.countNumberOfDifferences(alpha_J) == 2) && (beta_I.countNumberOfDifferences(beta_J) == 2)) { + + // Find the orbitals that are occupied in one string, and aren't in the other + size_t p = alpha_I.findOccupiedDifferences(alpha_J)[0]; // we're sure that there is only 1 element in the std::vector + size_t q = alpha_J.findOccupiedDifferences(alpha_I)[0]; // we're sure that there is only 1 element in the std::vector + + size_t r = beta_I.findOccupiedDifferences(beta_J)[0]; // we're sure that there is only 1 element in the std::vector + size_t s = beta_J.findOccupiedDifferences(beta_I)[0]; // we're sure that there is only 1 element in the std::vector + + // Calculate the total sign, and include it in the 2-RDM contribution + int sign = alpha_I.operatorPhaseFactor(p) * alpha_J.operatorPhaseFactor(q) * beta_I.operatorPhaseFactor(r) * beta_J.operatorPhaseFactor(s); + d_aabb(p,q,r,s) += sign * c_I * c_J; + d_aabb(q,p,s,r) += sign * c_I * c_J; + + d_bbaa(r,s,p,q) += sign * c_I * c_J; + d_bbaa(s,r,q,p) += sign * c_I * c_J; + } + + + // 2 electron excitations in alpha, 0 in beta + if ((alpha_I.countNumberOfDifferences(alpha_J) == 4) && (beta_I.countNumberOfDifferences(beta_J) == 0)) { + + // Find the orbitals that are occupied in one string, and aren't in the other + std::vector occupied_indices_I = alpha_I.findOccupiedDifferences(alpha_J); // we're sure this has two elements + size_t p = occupied_indices_I[0]; + size_t r = occupied_indices_I[1]; + + std::vector occupied_indices_J = alpha_J.findOccupiedDifferences(alpha_I); // we're sure this has two elements + size_t q = occupied_indices_J[0]; + size_t s = occupied_indices_J[1]; + + + // Calculate the total sign, and include it in the 2-RDM contribution + int sign = alpha_I.operatorPhaseFactor(p) * alpha_I.operatorPhaseFactor(r) * alpha_J.operatorPhaseFactor(q) * alpha_J.operatorPhaseFactor(s); + d_aaaa(p,q,r,s) += sign * c_I * c_J; + d_aaaa(p,s,r,q) -= sign * c_I * c_J; + d_aaaa(r,q,p,s) -= sign * c_I * c_J; + d_aaaa(r,s,p,q) += sign * c_I * c_J; + + d_aaaa(q,p,s,r) += sign * c_I * c_J; + d_aaaa(s,p,q,r) -= sign * c_I * c_J; + d_aaaa(q,r,s,p) -= sign * c_I * c_J; + d_aaaa(s,r,q,p) += sign * c_I * c_J; + } + + + // 0 electron excitations in alpha, 2 in beta + if ((alpha_I.countNumberOfDifferences(alpha_J) == 0) && (beta_I.countNumberOfDifferences(beta_J) == 4)) { + + // Find the orbitals that are occupied in one string, and aren't in the other + std::vector occupied_indices_I = beta_I.findOccupiedDifferences(beta_J); // we're sure this has two elements + size_t p = occupied_indices_I[0]; + size_t r = occupied_indices_I[1]; + + std::vector occupied_indices_J = beta_J.findOccupiedDifferences(beta_I); // we're sure this has two elements + size_t q = occupied_indices_J[0]; + size_t s = occupied_indices_J[1]; + + + // Calculate the total sign, and include it in the 2-RDM contribution + int sign = beta_I.operatorPhaseFactor(p) * beta_I.operatorPhaseFactor(r) * beta_J.operatorPhaseFactor(q) * beta_J.operatorPhaseFactor(s); + d_bbbb(p,q,r,s) += sign * c_I * c_J; + d_bbbb(p,s,r,q) -= sign * c_I * c_J; + d_bbbb(r,q,p,s) -= sign * c_I * c_J; + d_bbbb(r,s,p,q) += sign * c_I * c_J; + + d_bbbb(q,p,s,r) += sign * c_I * c_J; + d_bbbb(s,p,q,r) -= sign * c_I * c_J; + d_bbbb(q,r,s,p) -= sign * c_I * c_J; + d_bbbb(s,r,q,p) += sign * c_I * c_J; + } + + } // loop over all addresses J > I + + } // loop over all addresses I + + TwoRDM two_rdm_aaaa (d_aaaa); + TwoRDM two_rdm_aabb (d_aabb); + TwoRDM two_rdm_bbaa (d_bbaa); + TwoRDM two_rdm_bbbb (d_bbbb); + return TwoRDMs (two_rdm_aaaa, two_rdm_aabb, two_rdm_bbaa, two_rdm_bbbb); +} + + +} // namespace GQCP diff --git a/tests/RDM/SelectedRDMBuilder_test.cpp b/tests/RDM/SelectedRDMBuilder_test.cpp new file mode 100644 index 0000000000..7a8c6bb1e6 --- /dev/null +++ b/tests/RDM/SelectedRDMBuilder_test.cpp @@ -0,0 +1,129 @@ +// This file is part of GQCG-gqcp. +// +// Copyright (C) 2017-2018 the GQCG developers +// +// GQCG-gqcp is free software: you can redistribute it and/or modify +// it under the terms of the GNU Lesser General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// GQCG-gqcp is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with GQCG-gqcp. If not, see . +// +#define BOOST_TEST_MODULE "Selected_RDM_test" + + + +#include "RDM/RDMCalculator.hpp" + +#include "CISolver/CISolver.hpp" +#include "HamiltonianBuilder/FCI.hpp" +#include "HamiltonianParameters/HamiltonianParameters_constructors.hpp" + +#include + +#include +#include + + +BOOST_AUTO_TEST_CASE ( one_rdms_fci_H2_6_31G ) { + + // Do an H2@FCI//6-31G calculation + // test if 1-RDM SelectedRDM and FCIRDM are equal + size_t N_a = 1; + size_t N_b = 1; + + // Create a Molecule and an AOBasis + GQCP::Molecule h2("../tests/data/h2.xyz"); + auto ao_basis = std::make_shared(h2, "6-31G"); + + // Create the molecular Hamiltonian parameters for this molecule and basis + auto ham_par = GQCP::constructMolecularHamiltonianParameters(ao_basis); + size_t K = ham_par.get_K(); // 4 + + GQCP::FockSpaceProduct fock_space(K, N_a, N_b); // dim = 16 + GQCP::FCI fci(fock_space); + + // Specify solver options and solve the eigenvalue problem + // Solve the dense FCI eigenvalue problem + GQCP::CISolver ci_solver(fci, ham_par); + numopt::eigenproblem::DenseSolverOptions solver_options; + ci_solver.solve(solver_options); + + Eigen::VectorXd coef = ci_solver.get_eigenpair().get_eigenvector(); + + // Get the 1-RDM from FCI + GQCP::RDMCalculator fci_rdm(fock_space); + GQCP::OneRDMs one_rdms = fci_rdm.calculate1RDMs(coef); + + + GQCP::SelectedFockSpace selected_fock_space(K, N_a, N_b); + selected_fock_space.setExpansion(fock_space); + + // Get the 1-RDM from SelectedCI + GQCP::RDMCalculator selected_rdm(selected_fock_space); + GQCP::OneRDMs one_rdms_s = selected_rdm.calculate1RDMs(coef); + + + BOOST_CHECK(one_rdms_s.one_rdm.get_matrix_representation().isApprox( + one_rdms.one_rdm.get_matrix_representation())); + + BOOST_CHECK(one_rdms_s.one_rdm_aa.get_matrix_representation().isApprox( + one_rdms.one_rdm_aa.get_matrix_representation())); + + BOOST_CHECK(one_rdms_s.one_rdm_bb.get_matrix_representation().isApprox( + one_rdms.one_rdm_bb.get_matrix_representation())); +} + + + +BOOST_AUTO_TEST_CASE ( two_rdms_fci_H2_6_31G ) { + + // Do an H2@FCI//6-31G calculation + // test if 2-RDM SelectedRDM and FCIRDM are equal + size_t N_a = 1; + size_t N_b = 1; + + // Create a Molecule and an AOBasis + GQCP::Molecule h2("../tests/data/h2.xyz"); + auto ao_basis = std::make_shared(h2, "6-31G"); + + // Create the molecular Hamiltonian parameters for this molecule and basis + auto ham_par = GQCP::constructMolecularHamiltonianParameters(ao_basis); + size_t K = ham_par.get_K(); // 4 + + GQCP::FockSpaceProduct fock_space(K, N_a, N_b); // dim = 16 + GQCP::FCI fci(fock_space); + + // Specify solver options and solve the eigenvalue problem + // Solve the dense FCI eigenvalue problem + GQCP::CISolver ci_solver(fci, ham_par); + numopt::eigenproblem::DenseSolverOptions solver_options; + ci_solver.solve(solver_options); + + Eigen::VectorXd coef = ci_solver.get_eigenpair().get_eigenvector(); + + // Get the 1-RDM from FCI + GQCP::RDMCalculator fci_rdm(fock_space); + GQCP::TwoRDMs two_rdms = fci_rdm.calculate2RDMs(coef); + + + GQCP::SelectedFockSpace selected_fock_space(K, N_a, N_b); + selected_fock_space.setExpansion(fock_space); + + // Get the 1-RDM from SelectedCI + GQCP::RDMCalculator selected_rdm(selected_fock_space); + GQCP::TwoRDMs two_rdms_s = selected_rdm.calculate2RDMs(coef); + + + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_aaaa.get_matrix_representation(), two_rdms.two_rdm_aaaa.get_matrix_representation(), 1.0e-06)); + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_aabb.get_matrix_representation(), two_rdms.two_rdm_aabb.get_matrix_representation(), 1.0e-06)); + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_bbaa.get_matrix_representation(), two_rdms.two_rdm_bbaa.get_matrix_representation(), 1.0e-06)); + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_bbbb.get_matrix_representation(), two_rdms.two_rdm_bbbb.get_matrix_representation(), 1.0e-06)); + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm.get_matrix_representation(), two_rdms.two_rdm.get_matrix_representation(), 1.0e-06)); +} \ No newline at end of file From e15f884c0da7f796185e9aeeea551227845b764a Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Mon, 22 Oct 2018 17:50:22 +0200 Subject: [PATCH 09/13] added tests for selectedfockspace and rdms --- include/WaveFunction/WaveFunctionReader.hpp | 3 -- src/WaveFunction/WaveFunctionReader.cpp | 33 ++++++++------- tests/FockSpace/SelectedFockSpace_test.cpp | 47 ++++++++++++++++++--- tests/RDM/SelectedRDMBuilder_test.cpp | 7 ++- 4 files changed, 61 insertions(+), 29 deletions(-) diff --git a/include/WaveFunction/WaveFunctionReader.hpp b/include/WaveFunction/WaveFunctionReader.hpp index d5047b5b37..a6a1f541af 100644 --- a/include/WaveFunction/WaveFunctionReader.hpp +++ b/include/WaveFunction/WaveFunctionReader.hpp @@ -2,9 +2,6 @@ #define GQCP_WAVEFUNCTIONREADER_HPP - -#include - #include "FockSpace/SelectedFockSpace.hpp" #include "WaveFunction/WaveFunction.hpp" diff --git a/src/WaveFunction/WaveFunctionReader.cpp b/src/WaveFunction/WaveFunctionReader.cpp index ffee161ee5..70614d588a 100644 --- a/src/WaveFunction/WaveFunctionReader.cpp +++ b/src/WaveFunction/WaveFunctionReader.cpp @@ -18,13 +18,14 @@ #include - +#include +#include +#include #include "WaveFunction/WaveFunctionReader.hpp" - namespace GQCP { @@ -37,7 +38,7 @@ namespace GQCP { */ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) { - /* + // If the filename isn't properly converted into an input file stream, we assume the user supplied a wrong file std::ifstream input_file_stream(GAMESS_filename); std::ifstream input_file_stream_count(GAMESS_filename); // made to count the expansion size @@ -65,7 +66,7 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) size_t space_size = 0; // Count the amount of configurations - while (std::getline(input_file_stream, buffer)) { + while (std::getline(input_file_stream_count, buffer)) { if (buffer.length() != 0 | buffer[0] != '\n') { space_size++; } @@ -76,25 +77,25 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) std::getline(input_file_stream, line); // Read the first line containing the configurations - std::vector splitted_line; - boost::split(splitted_line, line, boost::is_any_of("|")); // split on '|' + std::vector splitted_linef; + boost::split(splitted_linef, line, boost::is_any_of("|")); // split on '|' // Create an alpha SpinString for the first field - std::string trimmed_alpha = boost::algorithm::trim_copy(splitted_line[0]); - std::string reversed_alpha (trimmed_alpha.rbegin(), trimmed_alpha.rend()); + std::string trimmed_alphaf = boost::algorithm::trim_copy(splitted_linef[0]); + std::string reversed_alphaf (trimmed_alphaf.rbegin(), trimmed_alphaf.rend()); // Create a beta SpinString for the second field - std::string trimmed_beta = boost::algorithm::trim_copy(splitted_line[1]); - std::string reversed_beta(trimmed_beta.rbegin(), trimmed_beta.rend()); + std::string trimmed_betaf = boost::algorithm::trim_copy(splitted_linef[1]); + std::string reversed_betaf(trimmed_betaf.rbegin(), trimmed_betaf.rend()); size_t counter = 0; // add the double - this->coefficients(counter) = std::stod(splitted_line[2]); + this->coefficients(counter) = std::stod(splitted_linef[2]); // dynamic bitset provides us with functionallity - boost::dynamic_bitset<> alpha_transfer (trimmed_alpha); - boost::dynamic_bitset<> beta_transfer (trimmed_beta); + boost::dynamic_bitset<> alpha_transfer (reversed_alphaf); + boost::dynamic_bitset<> beta_transfer (reversed_betaf); size_t K = alpha_transfer.size(); size_t N_alpha = alpha_transfer.count(); @@ -102,7 +103,7 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) this->fock_space = SelectedFockSpace(K, N_alpha, N_beta); - this->fock_space.addConfiguration(trimmed_alpha, trimmed_beta); + this->fock_space.addConfiguration(reversed_alphaf, reversed_betaf); // Read in the ONVs and the coefficients by splitting the line on '|', and then trimming whitespace @@ -131,12 +132,12 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) // Create a double for the third field this->coefficients(counter) = std::stod(splitted_line[2]); - this->fock_space.addConfiguration(trimmed_alpha, trimmed_beta); + this->fock_space.addConfiguration(reversed_alpha, reversed_beta); } // while getline this->wave_function = WaveFunction(this->fock_space, this->coefficients); - */ + } diff --git a/tests/FockSpace/SelectedFockSpace_test.cpp b/tests/FockSpace/SelectedFockSpace_test.cpp index fb39c9d0d3..f403886dc2 100644 --- a/tests/FockSpace/SelectedFockSpace_test.cpp +++ b/tests/FockSpace/SelectedFockSpace_test.cpp @@ -20,7 +20,7 @@ #include "FockSpace/SelectedFockSpace.hpp" -//#include "WaveFunction/WaveFunctionReader.hpp" +#include "WaveFunction/WaveFunctionReader.hpp" #include @@ -54,13 +54,48 @@ BOOST_AUTO_TEST_CASE ( addConfiguration ) { BOOST_CHECK_THROW(fock_space.addConfiguration("011", "011"), std::invalid_argument); } -/* + BOOST_AUTO_TEST_CASE ( reader_test ) { GQCP::WaveFunctionReader test_reader ("../tests/data/test_GAMESS_expansion"); - Eigen::VectorXd test = {1.0000, 0.0000}; + Eigen::Vector2d test_vector; + test_vector << 1.0000, 0.0000; + + // Check read vector is correct + BOOST_CHECK(test_vector.isApprox(test_reader.get_coefficients())); + + // Check if the expansions are equal + std::string alpha1 = "0000000000000000000000000000000000000000000001"; + std::string alpha2 = "0000000000000000000000000000000000000000000001"; + std::string beta1 = "0000000000000000000000000000000000000000000001"; + std::string beta2 = "0000000000000000000000000000000000000000000010"; + + GQCP::Configuration configuration1 = test_reader.get_fock_space().get_Configuration(0); + GQCP::Configuration configuration2 = test_reader.get_fock_space().get_Configuration(1); + + std::string alpha1_test; + std::string alpha2_test; + std::string beta1_test; + std::string beta2_test; + + std::ostringstream ss; + + ss << configuration1.onv_alpha; + alpha1_test = ss.str(); + ss.str(""); + ss << configuration1.onv_beta; + beta1_test = ss.str(); + ss.str(""); + ss << configuration2.onv_alpha; + alpha2_test = ss.str(); + ss.str(""); + ss << configuration2.onv_beta; + beta2_test = ss.str(); + ss.str(""); + + BOOST_CHECK(alpha1_test == alpha1); + BOOST_CHECK(alpha2_test == alpha2); + BOOST_CHECK(beta1_test == beta1); + BOOST_CHECK(beta2_test == beta2); - // Check if both expansions are considered equal - BOOST_CHECK(test.isApprox(test_reader.get_coefficients())); } - */ \ No newline at end of file diff --git a/tests/RDM/SelectedRDMBuilder_test.cpp b/tests/RDM/SelectedRDMBuilder_test.cpp index 7a8c6bb1e6..8d97e8c9b4 100644 --- a/tests/RDM/SelectedRDMBuilder_test.cpp +++ b/tests/RDM/SelectedRDMBuilder_test.cpp @@ -72,10 +72,8 @@ BOOST_AUTO_TEST_CASE ( one_rdms_fci_H2_6_31G ) { BOOST_CHECK(one_rdms_s.one_rdm.get_matrix_representation().isApprox( one_rdms.one_rdm.get_matrix_representation())); - BOOST_CHECK(one_rdms_s.one_rdm_aa.get_matrix_representation().isApprox( one_rdms.one_rdm_aa.get_matrix_representation())); - BOOST_CHECK(one_rdms_s.one_rdm_bb.get_matrix_representation().isApprox( one_rdms.one_rdm_bb.get_matrix_representation())); } @@ -119,11 +117,12 @@ BOOST_AUTO_TEST_CASE ( two_rdms_fci_H2_6_31G ) { // Get the 1-RDM from SelectedCI GQCP::RDMCalculator selected_rdm(selected_fock_space); GQCP::TwoRDMs two_rdms_s = selected_rdm.calculate2RDMs(coef); - + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_aaaa.get_matrix_representation(), two_rdms.two_rdm_aaaa.get_matrix_representation(), 1.0e-06)); BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_aabb.get_matrix_representation(), two_rdms.two_rdm_aabb.get_matrix_representation(), 1.0e-06)); BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_bbaa.get_matrix_representation(), two_rdms.two_rdm_bbaa.get_matrix_representation(), 1.0e-06)); BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_bbbb.get_matrix_representation(), two_rdms.two_rdm_bbbb.get_matrix_representation(), 1.0e-06)); BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm.get_matrix_representation(), two_rdms.two_rdm.get_matrix_representation(), 1.0e-06)); -} \ No newline at end of file +} + From a186ce08a8f559930b235c9a6b1aa6e40fe0b62d Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Tue, 23 Oct 2018 13:58:04 +0200 Subject: [PATCH 10/13] performed requested changes, added additional tests and functionality --- cmake/FindPackages.cmake | 2 +- cmake/SetCMakeEnvironment.cmake | 9 +- include/FockSpace/BaseFockSpace.hpp | 2 +- include/{ => FockSpace}/Configuration.hpp | 5 +- include/{ => FockSpace}/ONV.hpp | 25 ++-- include/FockSpace/SelectedFockSpace.hpp | 53 ++++++--- include/RDM/OneRDM.hpp | 3 - include/RDM/SelectedRDMBuilder.hpp | 8 +- include/WaveFunction/WaveFunction.hpp | 3 +- include/WaveFunction/WaveFunctionReader.hpp | 4 +- src/{ => FockSpace}/ONV.cpp | 35 +++--- src/FockSpace/SelectedFockSpace.cpp | 123 +++++++++++++------- src/RDM/SelectedRDMBuilder.cpp | 2 +- src/WaveFunction/WaveFunctionReader.cpp | 54 ++++----- tests/{ => FockSpace}/ONV_test.cpp | 4 +- tests/FockSpace/SelectedFockSpace_test.cpp | 90 ++++++++------ tests/RDM/SelectedRDMBuilder_test.cpp | 118 ++++++++++++++++--- 17 files changed, 361 insertions(+), 179 deletions(-) rename include/{ => FockSpace}/Configuration.hpp (69%) rename include/{ => FockSpace}/ONV.hpp (87%) rename src/{ => FockSpace}/ONV.cpp (91%) rename tests/{ => FockSpace}/ONV_test.cpp (99%) diff --git a/cmake/FindPackages.cmake b/cmake/FindPackages.cmake index c4bf25fd0a..d69940454c 100755 --- a/cmake/FindPackages.cmake +++ b/cmake/FindPackages.cmake @@ -2,7 +2,7 @@ # Find the Boost package -find_package(Boost REQUIRED COMPONENTS system thread) +find_package(Boost REQUIRED) # Find Eigen3 find_package(Eigen3 3.3.4 REQUIRED) diff --git a/cmake/SetCMakeEnvironment.cmake b/cmake/SetCMakeEnvironment.cmake index 0566aae695..23d0dffb89 100755 --- a/cmake/SetCMakeEnvironment.cmake +++ b/cmake/SetCMakeEnvironment.cmake @@ -31,6 +31,7 @@ set(PROJECT_SOURCE_FILES ${PROJECT_SOURCE_FOLDER}/FockSpace/BaseFockSpace.cpp ${PROJECT_SOURCE_FOLDER}/FockSpace/FockSpace.cpp ${PROJECT_SOURCE_FOLDER}/FockSpace/FockSpaceProduct.cpp + ${PROJECT_SOURCE_FOLDER}/FockSpace/ONV.cpp ${PROJECT_SOURCE_FOLDER}/FockSpace/SelectedFockSpace.cpp ${PROJECT_SOURCE_FOLDER}/HamiltonianBuilder/DOCI.cpp ${PROJECT_SOURCE_FOLDER}/HamiltonianBuilder/FCI.cpp @@ -46,7 +47,6 @@ set(PROJECT_SOURCE_FILES ${PROJECT_SOURCE_FOLDER}/RDM/DOCIRDMBuilder.cpp ${PROJECT_SOURCE_FOLDER}/RDM/FCIRDMBuilder.cpp ${PROJECT_SOURCE_FOLDER}/RDM/OneRDM.cpp - ${PROJECT_SOURCE_FOLDER}/RDM/SelectedRDMBuilder.cpp ${PROJECT_SOURCE_FOLDER}/RDM/RDMCalculator.cpp ${PROJECT_SOURCE_FOLDER}/RDM/RDMs.cpp ${PROJECT_SOURCE_FOLDER}/RDM/SelectedRDMBuilder.cpp @@ -65,7 +65,6 @@ set(PROJECT_SOURCE_FILES ${PROJECT_SOURCE_FOLDER}/LibintCommunicator.cpp ${PROJECT_SOURCE_FOLDER}/miscellaneous.cpp ${PROJECT_SOURCE_FOLDER}/Molecule.cpp - ${PROJECT_SOURCE_FOLDER}/ONV.cpp ${PROJECT_SOURCE_FOLDER}/RMP2.cpp) # Find the header folder @@ -79,9 +78,11 @@ set(PROJECT_INCLUDE_FILES ${PROJECT_INCLUDE_FOLDER}/AP1roG/AP1roGJacobiOrbitalOptimizer.hpp ${PROJECT_INCLUDE_FOLDER}/AP1roG/AP1roGPSESolver.hpp ${PROJECT_INCLUDE_FOLDER}/FockSpace/BaseFockSpace.hpp + ${PROJECT_INCLUDE_FOLDER}/FockSpace/Configuration.hpp ${PROJECT_INCLUDE_FOLDER}/FockSpace/FockSpace.hpp ${PROJECT_INCLUDE_FOLDER}/FockSpace/FockSpaceProduct.hpp ${PROJECT_INCLUDE_FOLDER}/FockSpace/FockSpaceType.hpp + ${PROJECT_INCLUDE_FOLDER}/FockSpace/ONV.hpp ${PROJECT_INCLUDE_FOLDER}/FockSpace/SelectedFockSpace.hpp ${PROJECT_INCLUDE_FOLDER}/HamiltonianBuilder/DOCI.hpp ${PROJECT_INCLUDE_FOLDER}/HamiltonianBuilder/FCI.hpp @@ -109,14 +110,12 @@ set(PROJECT_INCLUDE_FILES ${PROJECT_INCLUDE_FOLDER}/AOBasis.hpp ${PROJECT_INCLUDE_FOLDER}/Atom.hpp ${PROJECT_INCLUDE_FOLDER}/common.hpp - ${PROJECT_INCLUDE_FOLDER}/Configuration.hpp ${PROJECT_INCLUDE_FOLDER}/DOCINewtonOrbitalOptimizer.hpp ${PROJECT_INCLUDE_FOLDER}/elements.hpp ${PROJECT_INCLUDE_FOLDER}/JacobiRotationParameters.hpp ${PROJECT_INCLUDE_FOLDER}/LibintCommunicator.hpp ${PROJECT_INCLUDE_FOLDER}/miscellaneous.hpp ${PROJECT_INCLUDE_FOLDER}/Molecule.hpp - ${PROJECT_INCLUDE_FOLDER}/ONV.hpp ${PROJECT_INCLUDE_FOLDER}/RMP2.hpp ${PROJECT_INCLUDE_FOLDER}/units.hpp) @@ -136,6 +135,7 @@ set(PROJECT_TEST_SOURCE_FILES ${PROJECT_TESTS_FOLDER}/AP1roG/AP1roGPSESolver_test.cpp ${PROJECT_TESTS_FOLDER}/FockSpace/FockSpace_test.cpp ${PROJECT_TESTS_FOLDER}/FockSpace/FockSpaceProduct_test.cpp + ${PROJECT_TESTS_FOLDER}/FockSpace/ONV_test.cpp ${PROJECT_TESTS_FOLDER}/FockSpace/SelectedFockSpace_test.cpp ${PROJECT_TESTS_FOLDER}/HamiltonianBuilder/DOCI_test.cpp ${PROJECT_TESTS_FOLDER}/HamiltonianBuilder/FCI_test.cpp @@ -157,7 +157,6 @@ set(PROJECT_TEST_SOURCE_FILES ${PROJECT_TESTS_FOLDER}/LibintCommunicator_test.cpp ${PROJECT_TESTS_FOLDER}/miscellaneous_test.cpp ${PROJECT_TESTS_FOLDER}/Molecule_test.cpp - ${PROJECT_TESTS_FOLDER}/ONV_test.cpp ${PROJECT_TESTS_FOLDER}/OO_DOCI_test.cpp ${PROJECT_TESTS_FOLDER}/RMP2_test.cpp ${PROJECT_TESTS_FOLDER}/units_test.cpp) diff --git a/include/FockSpace/BaseFockSpace.hpp b/include/FockSpace/BaseFockSpace.hpp index 7fe5ee586a..3d2578bc40 100644 --- a/include/FockSpace/BaseFockSpace.hpp +++ b/include/FockSpace/BaseFockSpace.hpp @@ -41,11 +41,11 @@ class BaseFockSpace { // PROTECTED CONSTRUCTORS + BaseFockSpace() = default; /** * Protected constructor given a @param K and @param dim */ explicit BaseFockSpace(size_t K, size_t dim); - BaseFockSpace() = default; public: // DESTRUCTOR diff --git a/include/Configuration.hpp b/include/FockSpace/Configuration.hpp similarity index 69% rename from include/Configuration.hpp rename to include/FockSpace/Configuration.hpp index ef420c3150..04751a138c 100644 --- a/include/Configuration.hpp +++ b/include/FockSpace/Configuration.hpp @@ -3,12 +3,15 @@ #include "ONV.hpp" -#include "common.hpp" namespace GQCP { +/** + * Struct containing an electron distribution + * represented by the alpha and beta ONV. + */ struct Configuration { ONV onv_alpha; ONV onv_beta; diff --git a/include/ONV.hpp b/include/FockSpace/ONV.hpp similarity index 87% rename from include/ONV.hpp rename to include/FockSpace/ONV.hpp index 0f13d5c0b8..f06c092d60 100644 --- a/include/ONV.hpp +++ b/include/FockSpace/ONV.hpp @@ -33,7 +33,7 @@ namespace GQCP { * An example for 3 alpha electrons in a Fock space spanned by 4 spatial orbitals is * a_1^\dagger a_2^\dagger a_3^\dagger |vac> = |1,1,1,0> * - * In this code, we are using REVERSE LEXICAL notation, i.e. bitstrings are read from right to left. This means that the + * In this code bitstrings are read from right to left. This means that the * least significant bit relates to the first orbital. Using this notation is how normally bits are read, leading * to more efficient code. As is also usual, the least significant bit has index 0. * The previous example is then represented by the bit string "0111" (7). @@ -52,12 +52,17 @@ class ONV { // CONSTRUCTORS ONV() = default; /** - * Constructor from a @param K orbitals, N electrons and an @param unsigned_representation + * Constructor + * @param K a given number of orbitals + * @param N a given number of electrons + * @param unsigned_representation a representation for the ONV */ ONV(size_t K, size_t N, size_t unsigned_representation); /** - * Constructor from a @param K orbitals and an @param unsigned_representation + * Constructor + * @param K a given number of orbitals + * @param unsigned_representation a representation for the ONV */ ONV(size_t K, size_t unsigned_representation); @@ -103,7 +108,7 @@ class ONV { /** * @return if the @param p-th spatial orbital is occupied, starting from 0 - * @param p is the lexical index (i.e. read from right to left) + * @param p is counted from right to left */ bool isOccupied(size_t p) const; @@ -144,7 +149,7 @@ class ONV { bool create(size_t p, int& sign); /** - * @return the phase factor (+1 or -1) that arises by applying an annihilation or creation operator on orbital @param p, starting from 0 in reverse lexical ordering. + * @return the phase factor (+1 or -1) that arises by applying an annihilation or creation operator on orbital @param p, starting from 0, read from right to left. * * Let's say that there are m electrons in the orbitals up to p (not included). If m is even, the phase factor is (+1) and if m is odd, the phase factor is (-1), since electrons are fermions. */ @@ -155,8 +160,8 @@ class ONV { * @return the representation of a slice (i.e. a subset) of the spin string between @param index_start (included) * and @param index_end (not included). * - * Both @param index_start and @param index_end are 'lexical' (i.e. from right to left), which means that the slice - * occurs 'lexically' as well (i.e. from right to left). + * Both @param index_start and @param index_end are read from right to left, which means that the slice + * is from right to left as well. * * Example: * "010011".slice(1, 4) => "01[001]1" -> "001" @@ -172,11 +177,15 @@ class ONV { /** - * @return the positions of the bits (lexical indices) that are occupied in this, but unoccupied in @param other + * @return the positions of the bits (from right to left) that are occupied in this, but unoccupied in @param other */ std::vector findOccupiedDifferences(const ONV& other) const; + /** + * @return std::string containing the ONV representation + */ + std::string string_representation() const; // FRIEND CLASSES friend class FockSpace; diff --git a/include/FockSpace/SelectedFockSpace.hpp b/include/FockSpace/SelectedFockSpace.hpp index e78fd670ea..669fe98f9d 100644 --- a/include/FockSpace/SelectedFockSpace.hpp +++ b/include/FockSpace/SelectedFockSpace.hpp @@ -15,8 +15,8 @@ // You should have received a copy of the GNU Lesser General Public License // along with GQCG-gqcp. If not, see . // -#ifndef GQCP_SELECTEDSPACEPRODUCT_HPP -#define GQCP_SELECTEDSPACEPRODUCT_HPP +#ifndef GQCP_SELECTEDFOCKSPACE_HPP +#define GQCP_SELECTEDFOCKSPACE_HPP #include "FockSpace/BaseFockSpace.hpp" @@ -31,58 +31,73 @@ namespace GQCP { /** - * A Fock space for a given set of orbitals and number of alpha and beta electrons. - * Where the considered configurations are manually selected and represented as a Configuration : - * a combination of two ONVs, one holding the alpha configuration and holding the beta configuration. + * A Fock space that is flexible in the number of states that span it. + * The configurations are represented as a Configuration: + * a combination of two ONVs, holding the alpha and beta ONVs. */ class SelectedFockSpace : public GQCP::BaseFockSpace { private: size_t N_alpha; // number of alpha electrons size_t N_beta; // number of beta electrons - std::vector selection; // the selected configurations + std::vector configurations; /** - * Given @param onv1 and onv2, make a configuration + * Member taking two string arguments and creating a Configuration + * @param onv1 a string representation read from right to left + * @param onv2 a string representation read from right to left + * @return a Configuration + * !!! only works for up to 64 bits !!! */ - Configuration makeConfiguration(std::string onv1, std::string onv2); + Configuration makeConfiguration(const std::string& onv1, const std::string& onv2); public: // CONSTRUCTORS + SelectedFockSpace() = default; /** * Constructor given a @param K (spatial orbitals), N_alpha and N_beta (electrons); * the initial dimension of the space is 0, as no selections are made. */ SelectedFockSpace(size_t K, size_t N_alpha, size_t N_beta); - SelectedFockSpace() = default; - // DESTRUCTORS - ~SelectedFockSpace() override = default; + /** + * Constructor that generates expansion of a given FockSpaceProduct + * @param fock_space generated Fock space + */ + explicit SelectedFockSpace(const FockSpaceProduct& fock_space); + + /** + * Constructor that generates expansion of a given FockSpace + * @param fock_space generated Fock space + */ + explicit SelectedFockSpace(const FockSpace& fock_space); // GETTERS size_t get_N_alpha() const { return this->N_alpha; } size_t get_N_beta() const { return this->N_beta; } - Configuration get_Configuration(size_t index) const { return this->selection[index]; } + Configuration get_configuration(size_t index) const { return this->configurations[index]; } FockSpaceType get_type() const override { return FockSpaceType::SelectedFockSpace; } // PUBLIC METHODS /** - * Add configuration(s) (@param onv1(s) and onv2(s)) to the current selected Fock space - * increasing its selected dimension. + * Member taking two string arguments to add a Configuration + * @see makeConfiguration() + * add a Configuration to @var configurations */ - void addConfiguration(std::string onv1, std::string onv2); - void addConfiguration(std::vector onv1s, std::vector onv2s); + void addConfiguration(const std::string& onv1, const std::string& onv2); /** - * Sets the current expansion to that of @param fock_space + * Member taking two vector arguments to add Configurations + * @see makeConfiguration() + * add multiple Configurations to @var configurations */ - void setExpansion(const FockSpaceProduct& fock_space); + void addConfiguration(const std::vector& onv1s, const std::vector& onv2s); }; } // namespace GQCP -#endif // GQCP_SELECTEDSPACEPRODUCT_HPP +#endif // GQCP_SELECTEDFOCKSPACE_HPP diff --git a/include/RDM/OneRDM.hpp b/include/RDM/OneRDM.hpp index dc036e9376..bf001efcfb 100644 --- a/include/RDM/OneRDM.hpp +++ b/include/RDM/OneRDM.hpp @@ -43,9 +43,6 @@ class OneRDM : public BaseRDM { Eigen::MatrixXd get_matrix_representation() const { return this->D; } double get(size_t p, size_t q) const { return this->D(p, q); } - // OPERATOR - - // PUBLIC METHODS /** diff --git a/include/RDM/SelectedRDMBuilder.hpp b/include/RDM/SelectedRDMBuilder.hpp index 3c49e49f43..ff90c23136 100644 --- a/include/RDM/SelectedRDMBuilder.hpp +++ b/include/RDM/SelectedRDMBuilder.hpp @@ -29,10 +29,10 @@ namespace GQCP { /** * SelectedRDMBuilder is a class for the calculation of a density matrix from a given wave function - * or coefficient expansion in the full CI Fock space + * or coefficient expansion in a selected Fock space */ class SelectedRDMBuilder : public BaseRDMBuilder { - SelectedFockSpace fock_space; // Fock space containing the selected configuration + SelectedFockSpace fock_space; // Fock space containing the selected configurations public: @@ -40,10 +40,6 @@ class SelectedRDMBuilder : public BaseRDMBuilder { explicit SelectedRDMBuilder (const SelectedFockSpace& fock_space); - // DESTRUCTOR - ~SelectedRDMBuilder() = default; - - // OVERRIDDEN PUBLIC METHODS /** * @return all 1-RDMs from a coefficient vector @param x diff --git a/include/WaveFunction/WaveFunction.hpp b/include/WaveFunction/WaveFunction.hpp index c318c4a483..948956dfbb 100644 --- a/include/WaveFunction/WaveFunction.hpp +++ b/include/WaveFunction/WaveFunction.hpp @@ -38,8 +38,9 @@ class WaveFunction { public: // CONSTRUCTORS + WaveFunction() = default; WaveFunction(BaseFockSpace& base_fock_space, const Eigen::VectorXd& coefficients); - WaveFunction() = default; + // GETTERS Eigen::VectorXd get_coefficients() const { return coefficients; } diff --git a/include/WaveFunction/WaveFunctionReader.hpp b/include/WaveFunction/WaveFunctionReader.hpp index a6a1f541af..b16949a729 100644 --- a/include/WaveFunction/WaveFunctionReader.hpp +++ b/include/WaveFunction/WaveFunctionReader.hpp @@ -5,8 +5,6 @@ #include "FockSpace/SelectedFockSpace.hpp" #include "WaveFunction/WaveFunction.hpp" -#include - namespace GQCP { @@ -39,4 +37,4 @@ class WaveFunctionReader { -#endif //GQCP_WAVEFUNCTIONREADER_HPP +#endif // GQCP_WAVEFUNCTIONREADER_HPP diff --git a/src/ONV.cpp b/src/FockSpace/ONV.cpp similarity index 91% rename from src/ONV.cpp rename to src/FockSpace/ONV.cpp index 6c3d5a9fc6..820747c6be 100644 --- a/src/ONV.cpp +++ b/src/FockSpace/ONV.cpp @@ -15,7 +15,7 @@ // You should have received a copy of the GNU Lesser General Public License // along with GQCG-gqcp. If not, see . // -#include "ONV.hpp" +#include "FockSpace/ONV.hpp" #include @@ -31,9 +31,9 @@ namespace GQCP { * Constructor from a @param K orbitals, N electrons and an @param unsigned_representation */ ONV::ONV(size_t K, size_t N, size_t representation): - K(K), - N(N), - unsigned_representation(representation) + K (K), + N (N), + unsigned_representation (representation) { occupation_indices = VectorXs::Zero(N); this->updateOccupationIndices(); // throws error if the representation and N are not compatible @@ -43,9 +43,7 @@ ONV::ONV(size_t K, size_t N, size_t representation): * Constructor from a @param K orbitals and an @param unsigned_representation */ ONV::ONV(size_t K, size_t representation): - K(K), - N(__builtin_popcountl(representation)), - unsigned_representation(representation) + ONV(K, __builtin_popcountl(representation), representation) { occupation_indices = VectorXs::Zero(N); this->updateOccupationIndices(); // throws error if the representation and N are not compatible @@ -60,7 +58,7 @@ ONV::ONV(size_t K, size_t representation): * Overloading of operator<< for a GQCP::ONV to be used with streams */ std::ostream& operator<<(std::ostream& os, const GQCP::ONV& onv) { - return os< (onv.K, onv.unsigned_representation); + return os< "01[001]1" -> "001" @@ -266,7 +264,7 @@ size_t ONV::countNumberOfDifferences(const ONV& other) const { /** - * @return the positions of the bits (lexical indices) that are occupied in this, but unoccupied in @param other + * @return the positions of the bits (from right to left) that are occupied in this, but unoccupied in @param other */ std::vector ONV::findOccupiedDifferences(const ONV& other) const { @@ -289,4 +287,15 @@ std::vector ONV::findOccupiedDifferences(const ONV& other) const { } +/** + * @return std::string containing the ONV representation + */ +std::string ONV::string_representation() const { + boost::dynamic_bitset<> transfer_set (this->K, this->unsigned_representation); + std::string buffer; + boost::to_string(transfer_set, buffer); + return buffer; +} + + } // namespace GQCP diff --git a/src/FockSpace/SelectedFockSpace.cpp b/src/FockSpace/SelectedFockSpace.cpp index 8b14ce7c1e..bdf1bf7c11 100644 --- a/src/FockSpace/SelectedFockSpace.cpp +++ b/src/FockSpace/SelectedFockSpace.cpp @@ -29,19 +29,23 @@ namespace GQCP { */ /** - * Given @param onv1 and onv2, make a configuration + * Member taking two string arguments and creating a Configuration + * @param onv1 a string representation read from right to left + * @param onv2 a string representation read from right to left + * @return a Configuration + * !!! only works for up to 64 bits !!! */ -Configuration SelectedFockSpace::makeConfiguration(std::string onv1, std::string onv2){ +Configuration SelectedFockSpace::makeConfiguration(const std::string& onv1, const std::string& onv2){ boost::dynamic_bitset<> alpha_transfer (onv1); boost::dynamic_bitset<> beta_transfer (onv2); if (alpha_transfer.size() != this->K | beta_transfer.size() != this->K) { - throw std::invalid_argument("Added configuration is not compatible with the orbitals number of the Fock space"); + throw std::invalid_argument("Given string representations for ONVs are not compatible with the number of orbitals of the Fock space"); } if (alpha_transfer.count() != this->N_alpha | beta_transfer.count() != this->N_beta) { - throw std::invalid_argument("Added configuration is not compatible with the electron numbers of the Fock space"); + throw std::invalid_argument("Given string representations for ONVs are not compatible with the number of orbitals of the Fock space"); } size_t alpha_s = alpha_transfer.to_ulong(); @@ -50,7 +54,7 @@ Configuration SelectedFockSpace::makeConfiguration(std::string onv1, std::string ONV alpha (alpha_transfer.size(), alpha_s); ONV beta (beta_transfer.size(), beta_s); - return Configuration { alpha, beta }; + return Configuration {alpha, beta}; } @@ -69,39 +73,13 @@ SelectedFockSpace::SelectedFockSpace(size_t K, size_t N_alpha, size_t N_beta) : N_beta (N_beta) {} - - -/* - * PUBLIC METHODS - */ - /** - * Add configuration(s) (@param onv1(s) and onv2(s)) to the current selected Fock space - * increasing its selected dimension. + * Constructor that generates expansion of a given FockSpaceProduct + * @param fock_space generated Fock space */ -void SelectedFockSpace::addConfiguration(std::string onv1, std::string onv2){ - - this->dim++; - - Configuration configuration = makeConfiguration(onv1, onv2); - selection.push_back(configuration); -} - -void SelectedFockSpace::addConfiguration(std::vector onv1s, std::vector onv2s){ - - if (onv1s.size() != onv2s.size()) { - throw std::invalid_argument("Size of both ONV entry vectors do not match"); - } - - for (int i = 0; iaddConfiguration(onv1s[i], onv2s[i]); - } -} - -/** - * Sets the current expansion to that of @param fock_space - */ -void SelectedFockSpace::setExpansion(const FockSpaceProduct& fock_space){ +SelectedFockSpace::SelectedFockSpace(const FockSpaceProduct& fock_space) : + SelectedFockSpace (fock_space.get_K(), fock_space.get_N_alpha(), fock_space.get_N_beta()) +{ std::vector configurations; @@ -119,22 +97,85 @@ void SelectedFockSpace::setExpansion(const FockSpaceProduct& fock_space){ configurations.push_back(Configuration {alpha, beta}); - if (I_beta < dim_beta - 1) { // prevent the last permutation to occur fock_space_beta.setNext(beta); } - } - if (I_alpha < dim_alpha - 1) { // prevent the last permutation to occur fock_space_alpha.setNext(alpha); } - } this->dim = fock_space.get_dimension(); - this->selection = configurations; + this->configurations = configurations; + +} + +/** + * Constructor that generates expansion of a given FockSpace + * @param fock_space generated Fock space + */ +SelectedFockSpace::SelectedFockSpace(const FockSpace& fock_space) : + SelectedFockSpace (fock_space.get_K(), fock_space.get_N(), fock_space.get_N()) +{ + + std::vector configurations; + // Current workaround to call non-const functions + FockSpace fock_space_single = fock_space; + + + auto dim = fock_space_single.get_dimension(); + + // Iterate over the Fock space and add all onvs as doubly occupied configurations + ONV onv = fock_space_single.get_ONV(0); + for (size_t I = 0; I < dim; I++) { + + configurations.push_back(Configuration {onv, onv}); + + if (I < dim - 1) { // prevent the last permutation to occur + fock_space_single.setNext(onv); + } + + } + + this->dim = dim; + this->configurations = configurations; } + +/* + * PUBLIC METHODS + */ + +/** + * Member taking two string arguments to add a Configuration + * @see makeConfiguration() + * add a Configuration to @var configurations + */ +void SelectedFockSpace::addConfiguration(const std::string& onv1, const std::string& onv2){ + + this->dim++; + + Configuration configuration = makeConfiguration(onv1, onv2); + configurations.push_back(configuration); +} + +/** + * Member taking two vector arguments to add Configurations + * @see makeConfiguration() + * add multiple Configurations to @var configurations + */ +void SelectedFockSpace::addConfiguration(const std::vector& onv1s, const std::vector& onv2s){ + + if (onv1s.size() != onv2s.size()) { + throw std::invalid_argument("Size of both ONV entry vectors do not match"); + } + + for (int i = 0; iaddConfiguration(onv1s[i], onv2s[i]); + } +} + + } // namespace GQCP diff --git a/src/RDM/SelectedRDMBuilder.cpp b/src/RDM/SelectedRDMBuilder.cpp index 75cc5f5e38..0b52455695 100644 --- a/src/RDM/SelectedRDMBuilder.cpp +++ b/src/RDM/SelectedRDMBuilder.cpp @@ -17,7 +17,7 @@ // #include "RDM/SelectedRDMBuilder.hpp" -#include "Configuration.hpp" +#include "FockSpace/Configuration.hpp" namespace GQCP { diff --git a/src/WaveFunction/WaveFunctionReader.cpp b/src/WaveFunction/WaveFunctionReader.cpp index 70614d588a..0a96f73833 100644 --- a/src/WaveFunction/WaveFunctionReader.cpp +++ b/src/WaveFunction/WaveFunctionReader.cpp @@ -40,8 +40,8 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) { // If the filename isn't properly converted into an input file stream, we assume the user supplied a wrong file - std::ifstream input_file_stream(GAMESS_filename); - std::ifstream input_file_stream_count(GAMESS_filename); // made to count the expansion size + std::ifstream input_file_stream (GAMESS_filename); + std::ifstream input_file_stream_count (GAMESS_filename); // made to count the expansion size if (!input_file_stream.good()) { throw std::runtime_error("WaveFunctionReader(): The provided GAMESS file is illegible. Maybe you specified a wrong path?"); @@ -50,7 +50,7 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) // Do the actual parsing // Read in dummy lines up until we actually get to the ONVs and coefficients std::string line; - std::string buffer; // dummy for the getline() TODO: find "correcter" way if possible + std::string buffer; // dummy for the counting stream TODO: find "correcter" way if possible while (std::getline(input_file_stream, line)) { std::getline(input_file_stream_count, buffer); if (line.find("ALPHA") != std::string::npos @@ -77,25 +77,27 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) std::getline(input_file_stream, line); // Read the first line containing the configurations - std::vector splitted_linef; - boost::split(splitted_linef, line, boost::is_any_of("|")); // split on '|' + std::vector splitted_line; + boost::split(splitted_line, line, boost::is_any_of("|")); // split on '|' - // Create an alpha SpinString for the first field - std::string trimmed_alphaf = boost::algorithm::trim_copy(splitted_linef[0]); - std::string reversed_alphaf (trimmed_alphaf.rbegin(), trimmed_alphaf.rend()); + // Create an alpha ONV for the first field + std::string trimmed_alpha = boost::algorithm::trim_copy(splitted_line[0]); + std::string reversed_alpha (trimmed_alpha.rbegin(), trimmed_alpha.rend()); - // Create a beta SpinString for the second field - std::string trimmed_betaf = boost::algorithm::trim_copy(splitted_linef[1]); - std::string reversed_betaf(trimmed_betaf.rbegin(), trimmed_betaf.rend()); + // Create a beta ONV for the second field + std::string trimmed_beta = boost::algorithm::trim_copy(splitted_line[1]); + std::string reversed_beta (trimmed_beta.rbegin(), trimmed_beta.rend()); + + + size_t index_count = 0; // counts the configurations added - size_t counter = 0; // add the double - this->coefficients(counter) = std::stod(splitted_linef[2]); + this->coefficients(index_count) = std::stod(splitted_line[2]); - // dynamic bitset provides us with functionallity - boost::dynamic_bitset<> alpha_transfer (reversed_alphaf); - boost::dynamic_bitset<> beta_transfer (reversed_betaf); + // Construct dynamic bitsets to use its functionality + boost::dynamic_bitset<> alpha_transfer (reversed_alpha); + boost::dynamic_bitset<> beta_transfer (reversed_beta); size_t K = alpha_transfer.size(); size_t N_alpha = alpha_transfer.count(); @@ -103,41 +105,39 @@ WaveFunctionReader::WaveFunctionReader(const std::string& GAMESS_filename) this->fock_space = SelectedFockSpace(K, N_alpha, N_beta); - this->fock_space.addConfiguration(reversed_alphaf, reversed_betaf); + this->fock_space.addConfiguration(reversed_alpha, reversed_beta); // Read in the ONVs and the coefficients by splitting the line on '|', and then trimming whitespace // In the GAMESS format, the bit strings are ordered in reverse while (std::getline(input_file_stream, line)) { - counter++; - std::vector splitted_line; + index_count++; boost::split(splitted_line, line, boost::is_any_of("|")); // split on '|' - // Create an alpha SpinString for the first field - std::string trimmed_alpha = boost::algorithm::trim_copy(splitted_line[0]); + // Create an alpha ONV for the first field + trimmed_alpha = boost::algorithm::trim_copy(splitted_line[0]); if (trimmed_alpha.length() != K) { throw std::invalid_argument("WaveFunctionReader(): One of the provided alpha ONVs does not have the correct number of orbitals."); } - std::string reversed_alpha (trimmed_alpha.rbegin(), trimmed_alpha.rend()); + reversed_alpha = std::string(trimmed_alpha.rbegin(), trimmed_alpha.rend()); - // Create a beta SpinString for the second field - std::string trimmed_beta = boost::algorithm::trim_copy(splitted_line[1]); + // Create a beta ONV for the second field + trimmed_beta = boost::algorithm::trim_copy(splitted_line[1]); if (trimmed_beta.length() != K) { throw std::invalid_argument("WaveFunctionReader(): One of the provided beta ONVs does not have the correct number of orbitals."); } - std::string reversed_beta(trimmed_beta.rbegin(), trimmed_beta.rend()); + reversed_beta = std::string(trimmed_beta.rbegin(), trimmed_beta.rend()); // Create a double for the third field - this->coefficients(counter) = std::stod(splitted_line[2]); + this->coefficients(index_count) = std::stod(splitted_line[2]); this->fock_space.addConfiguration(reversed_alpha, reversed_beta); } // while getline this->wave_function = WaveFunction(this->fock_space, this->coefficients); - } diff --git a/tests/ONV_test.cpp b/tests/FockSpace/ONV_test.cpp similarity index 99% rename from tests/ONV_test.cpp rename to tests/FockSpace/ONV_test.cpp index 29d543baf0..669b64e2ec 100644 --- a/tests/ONV_test.cpp +++ b/tests/FockSpace/ONV_test.cpp @@ -18,7 +18,7 @@ #define BOOST_TEST_MODULE "ONV" -#include "ONV.hpp" +#include "FockSpace/ONV.hpp" #include "FockSpace/FockSpace.hpp" #include @@ -306,4 +306,4 @@ BOOST_AUTO_TEST_CASE ( findOccupiedDifferences_unsigned ) { BOOST_TEST(spin_string2.findOccupiedDifferences(spin_string3) == (std::vector {2}), boost::test_tools::per_element()); BOOST_TEST(spin_string3.findOccupiedDifferences(spin_string2) == (std::vector {3}), boost::test_tools::per_element()); -} \ No newline at end of file +} diff --git a/tests/FockSpace/SelectedFockSpace_test.cpp b/tests/FockSpace/SelectedFockSpace_test.cpp index f403886dc2..6996ec7cc5 100644 --- a/tests/FockSpace/SelectedFockSpace_test.cpp +++ b/tests/FockSpace/SelectedFockSpace_test.cpp @@ -30,6 +30,12 @@ BOOST_AUTO_TEST_CASE ( constructor ) { BOOST_CHECK_NO_THROW(GQCP::SelectedFockSpace (10, 5, 5)); + + GQCP::FockSpaceProduct fock_space_product (10, 5, 5); + GQCP::FockSpace fock_space (10, 5); + + BOOST_CHECK_NO_THROW(GQCP::SelectedFockSpace fock (fock_space_product)); + BOOST_CHECK_NO_THROW(GQCP::SelectedFockSpace fock (fock_space)); } @@ -52,50 +58,68 @@ BOOST_AUTO_TEST_CASE ( addConfiguration ) { // Test throw with incompatible electron numbers BOOST_CHECK_THROW(fock_space.addConfiguration("011", "011"), std::invalid_argument); + + fock_space.addConfiguration(alpha_set, beta_set); + + + // Check if the expansions are equal + // Generate the expected results + std::string alpha1_ref = "001"; + std::string alpha2_ref = "010"; + std::string beta1_ref = "001"; + std::string beta2_ref = "010"; + + // Retrieve the added results + GQCP::Configuration configuration1 = fock_space.get_configuration(0); + GQCP::Configuration configuration2 = fock_space.get_configuration(1); + + // Retrieve the string representation of the ONVs + std::string alpha1_test = configuration1.onv_alpha.string_representation(); + std::string alpha2_test = configuration2.onv_alpha.string_representation(); + std::string beta1_test = configuration1.onv_beta.string_representation(); + std::string beta2_test = configuration2.onv_beta.string_representation(); + + BOOST_CHECK(alpha1_test == alpha1_ref); + BOOST_CHECK(alpha2_test == alpha2_ref); + BOOST_CHECK(beta1_test == beta1_ref); + BOOST_CHECK(beta2_test == beta2_ref); } BOOST_AUTO_TEST_CASE ( reader_test ) { + // We will test if we can construct a selected fock space and a corresponding coefficients + // through an input file GQCP::WaveFunctionReader test_reader ("../tests/data/test_GAMESS_expansion"); + + + // Check read vector is correct + // Gerenate the expected result Eigen::Vector2d test_vector; test_vector << 1.0000, 0.0000; - // Check read vector is correct BOOST_CHECK(test_vector.isApprox(test_reader.get_coefficients())); // Check if the expansions are equal - std::string alpha1 = "0000000000000000000000000000000000000000000001"; - std::string alpha2 = "0000000000000000000000000000000000000000000001"; - std::string beta1 = "0000000000000000000000000000000000000000000001"; - std::string beta2 = "0000000000000000000000000000000000000000000010"; - - GQCP::Configuration configuration1 = test_reader.get_fock_space().get_Configuration(0); - GQCP::Configuration configuration2 = test_reader.get_fock_space().get_Configuration(1); - - std::string alpha1_test; - std::string alpha2_test; - std::string beta1_test; - std::string beta2_test; - - std::ostringstream ss; - - ss << configuration1.onv_alpha; - alpha1_test = ss.str(); - ss.str(""); - ss << configuration1.onv_beta; - beta1_test = ss.str(); - ss.str(""); - ss << configuration2.onv_alpha; - alpha2_test = ss.str(); - ss.str(""); - ss << configuration2.onv_beta; - beta2_test = ss.str(); - ss.str(""); - - BOOST_CHECK(alpha1_test == alpha1); - BOOST_CHECK(alpha2_test == alpha2); - BOOST_CHECK(beta1_test == beta1); - BOOST_CHECK(beta2_test == beta2); + // Generate the expected results + std::string alpha1_ref = "0000000000000000000000000000000000000000000001"; + std::string alpha2_ref = "0000000000000000000000000000000000000000000001"; + std::string beta1_ref = "0000000000000000000000000000000000000000000001"; + std::string beta2_ref = "0000000000000000000000000000000000000000000010"; + + // Retrieve the read results + GQCP::Configuration configuration1 = test_reader.get_fock_space().get_configuration(0); + GQCP::Configuration configuration2 = test_reader.get_fock_space().get_configuration(1); + + // Retrieve the string representation of the ONVs + std::string alpha1_test = configuration1.onv_alpha.string_representation(); + std::string alpha2_test = configuration2.onv_alpha.string_representation(); + std::string beta1_test = configuration1.onv_beta.string_representation(); + std::string beta2_test = configuration2.onv_beta.string_representation(); + + BOOST_CHECK(alpha1_test == alpha1_ref); + BOOST_CHECK(alpha2_test == alpha2_ref); + BOOST_CHECK(beta1_test == beta1_ref); + BOOST_CHECK(beta2_test == beta2_ref); } diff --git a/tests/RDM/SelectedRDMBuilder_test.cpp b/tests/RDM/SelectedRDMBuilder_test.cpp index 8d97e8c9b4..5f5e8adfeb 100644 --- a/tests/RDM/SelectedRDMBuilder_test.cpp +++ b/tests/RDM/SelectedRDMBuilder_test.cpp @@ -23,6 +23,7 @@ #include "CISolver/CISolver.hpp" #include "HamiltonianBuilder/FCI.hpp" +#include "HamiltonianBuilder/DOCI.hpp" #include "HamiltonianParameters/HamiltonianParameters_constructors.hpp" #include @@ -39,19 +40,19 @@ BOOST_AUTO_TEST_CASE ( one_rdms_fci_H2_6_31G ) { size_t N_b = 1; // Create a Molecule and an AOBasis - GQCP::Molecule h2("../tests/data/h2.xyz"); + GQCP::Molecule h2 ("../tests/data/h2.xyz"); auto ao_basis = std::make_shared(h2, "6-31G"); // Create the molecular Hamiltonian parameters for this molecule and basis auto ham_par = GQCP::constructMolecularHamiltonianParameters(ao_basis); size_t K = ham_par.get_K(); // 4 - GQCP::FockSpaceProduct fock_space(K, N_a, N_b); // dim = 16 - GQCP::FCI fci(fock_space); + GQCP::FockSpaceProduct fock_space (K, N_a, N_b); // dim = 16 + GQCP::FCI fci (fock_space); // Specify solver options and solve the eigenvalue problem // Solve the dense FCI eigenvalue problem - GQCP::CISolver ci_solver(fci, ham_par); + GQCP::CISolver ci_solver (fci, ham_par); numopt::eigenproblem::DenseSolverOptions solver_options; ci_solver.solve(solver_options); @@ -62,11 +63,10 @@ BOOST_AUTO_TEST_CASE ( one_rdms_fci_H2_6_31G ) { GQCP::OneRDMs one_rdms = fci_rdm.calculate1RDMs(coef); - GQCP::SelectedFockSpace selected_fock_space(K, N_a, N_b); - selected_fock_space.setExpansion(fock_space); + GQCP::SelectedFockSpace selected_fock_space (fock_space); // Get the 1-RDM from SelectedCI - GQCP::RDMCalculator selected_rdm(selected_fock_space); + GQCP::RDMCalculator selected_rdm (selected_fock_space); GQCP::OneRDMs one_rdms_s = selected_rdm.calculate1RDMs(coef); @@ -88,19 +88,19 @@ BOOST_AUTO_TEST_CASE ( two_rdms_fci_H2_6_31G ) { size_t N_b = 1; // Create a Molecule and an AOBasis - GQCP::Molecule h2("../tests/data/h2.xyz"); + GQCP::Molecule h2 ("../tests/data/h2.xyz"); auto ao_basis = std::make_shared(h2, "6-31G"); // Create the molecular Hamiltonian parameters for this molecule and basis auto ham_par = GQCP::constructMolecularHamiltonianParameters(ao_basis); size_t K = ham_par.get_K(); // 4 - GQCP::FockSpaceProduct fock_space(K, N_a, N_b); // dim = 16 - GQCP::FCI fci(fock_space); + GQCP::FockSpaceProduct fock_space (K, N_a, N_b); // dim = 16 + GQCP::FCI fci (fock_space); // Specify solver options and solve the eigenvalue problem // Solve the dense FCI eigenvalue problem - GQCP::CISolver ci_solver(fci, ham_par); + GQCP::CISolver ci_solver (fci, ham_par); numopt::eigenproblem::DenseSolverOptions solver_options; ci_solver.solve(solver_options); @@ -111,11 +111,10 @@ BOOST_AUTO_TEST_CASE ( two_rdms_fci_H2_6_31G ) { GQCP::TwoRDMs two_rdms = fci_rdm.calculate2RDMs(coef); - GQCP::SelectedFockSpace selected_fock_space(K, N_a, N_b); - selected_fock_space.setExpansion(fock_space); + GQCP::SelectedFockSpace selected_fock_space (fock_space); // Get the 1-RDM from SelectedCI - GQCP::RDMCalculator selected_rdm(selected_fock_space); + GQCP::RDMCalculator selected_rdm (selected_fock_space); GQCP::TwoRDMs two_rdms_s = selected_rdm.calculate2RDMs(coef); @@ -126,3 +125,94 @@ BOOST_AUTO_TEST_CASE ( two_rdms_fci_H2_6_31G ) { BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm.get_matrix_representation(), two_rdms.two_rdm.get_matrix_representation(), 1.0e-06)); } + +BOOST_AUTO_TEST_CASE ( one_rdms_doci_H2_6_31G ) { + + // Do an H2@doci//6-31G calculation + // test if 1-RDM SelectedRDM and dociRDM are equal + size_t N = 1; + + // Create a Molecule and an AOBasis + GQCP::Molecule h2 ("../tests/data/h2.xyz"); + auto ao_basis = std::make_shared(h2, "6-31G"); + + // Create the molecular Hamiltonian parameters for this molecule and basis + auto ham_par = GQCP::constructMolecularHamiltonianParameters(ao_basis); + size_t K = ham_par.get_K(); // 4 + + GQCP::FockSpace fock_space (K, N); // dim = 4 + GQCP::DOCI doci (fock_space); + + // Specify solver options and solve the eigenvalue problem + // Solve the dense doci eigenvalue problem + GQCP::CISolver ci_solver (doci, ham_par); + numopt::eigenproblem::DenseSolverOptions solver_options; + ci_solver.solve(solver_options); + + Eigen::VectorXd coef = ci_solver.get_eigenpair().get_eigenvector(); + + // Get the 1-RDM from doci + GQCP::RDMCalculator doci_rdm(fock_space); + GQCP::OneRDMs one_rdms = doci_rdm.calculate1RDMs(coef); + + + GQCP::SelectedFockSpace selected_fock_space (fock_space); + + // Get the 1-RDM from SelectedCI + GQCP::RDMCalculator selected_rdm (selected_fock_space); + GQCP::OneRDMs one_rdms_s = selected_rdm.calculate1RDMs(coef); + + + BOOST_CHECK(one_rdms_s.one_rdm.get_matrix_representation().isApprox( + one_rdms.one_rdm.get_matrix_representation())); + BOOST_CHECK(one_rdms_s.one_rdm_aa.get_matrix_representation().isApprox( + one_rdms.one_rdm_aa.get_matrix_representation())); + BOOST_CHECK(one_rdms_s.one_rdm_bb.get_matrix_representation().isApprox( + one_rdms.one_rdm_bb.get_matrix_representation())); +} + + + +BOOST_AUTO_TEST_CASE ( two_rdms_doci_H2_6_31G ) { + + // Do an H2@doci//6-31G calculation + // test if 2-RDM SelectedRDM and dociRDM are equal + size_t N = 1; + + // Create a Molecule and an AOBasis + GQCP::Molecule h2 ("../tests/data/h2.xyz"); + auto ao_basis = std::make_shared(h2, "6-31G"); + + // Create the molecular Hamiltonian parameters for this molecule and basis + auto ham_par = GQCP::constructMolecularHamiltonianParameters(ao_basis); + size_t K = ham_par.get_K(); // 4 + + GQCP::FockSpace fock_space (K, N); // dim = 4 + GQCP::DOCI doci (fock_space); + + // Specify solver options and solve the eigenvalue problem + // Solve the dense doci eigenvalue problem + GQCP::CISolver ci_solver (doci, ham_par); + numopt::eigenproblem::DenseSolverOptions solver_options; + ci_solver.solve(solver_options); + + Eigen::VectorXd coef = ci_solver.get_eigenpair().get_eigenvector(); + + // Get the 1-RDM from doci + GQCP::RDMCalculator doci_rdm(fock_space); + GQCP::TwoRDMs two_rdms = doci_rdm.calculate2RDMs(coef); + + + GQCP::SelectedFockSpace selected_fock_space (fock_space); + + // Get the 1-RDM from SelectedCI + GQCP::RDMCalculator selected_rdm (selected_fock_space); + GQCP::TwoRDMs two_rdms_s = selected_rdm.calculate2RDMs(coef); + + + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_aaaa.get_matrix_representation(), two_rdms.two_rdm_aaaa.get_matrix_representation(), 1.0e-06)); + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_aabb.get_matrix_representation(), two_rdms.two_rdm_aabb.get_matrix_representation(), 1.0e-06)); + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_bbaa.get_matrix_representation(), two_rdms.two_rdm_bbaa.get_matrix_representation(), 1.0e-06)); + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm_bbbb.get_matrix_representation(), two_rdms.two_rdm_bbbb.get_matrix_representation(), 1.0e-06)); + BOOST_CHECK(cpputil::linalg::areEqual(two_rdms_s.two_rdm.get_matrix_representation(), two_rdms.two_rdm.get_matrix_representation(), 1.0e-06)); +} From b16c4b78d8f241f8f629c408b773fa2b985c74d0 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Tue, 23 Oct 2018 14:13:39 +0200 Subject: [PATCH 11/13] fixed typo in RDMBuilder --- src/RDM/SelectedRDMBuilder.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/RDM/SelectedRDMBuilder.cpp b/src/RDM/SelectedRDMBuilder.cpp index 0b52455695..a225fa3de0 100644 --- a/src/RDM/SelectedRDMBuilder.cpp +++ b/src/RDM/SelectedRDMBuilder.cpp @@ -50,7 +50,7 @@ OneRDMs SelectedRDMBuilder::calculate1RDMs(const Eigen::VectorXd& x) { size_t dim = fock_space.get_dimension(); for (size_t I = 0; I < dim; I++) { // loop over all addresses (1) - Configuration configuration_I = this->fock_space.get_Configuration(I); + Configuration configuration_I = this->fock_space.get_configuration(I); ONV alpha_I = configuration_I.onv_alpha; ONV beta_I = configuration_I.onv_beta; @@ -73,7 +73,7 @@ OneRDMs SelectedRDMBuilder::calculate1RDMs(const Eigen::VectorXd& x) { // Calculate the off-diagonal elements, by going over all other ONVs for (size_t J = I+1; J < dim; J++) { - Configuration configuration_J = this->fock_space.get_Configuration(J); + Configuration configuration_J = this->fock_space.get_configuration(J); ONV alpha_J = configuration_J.onv_alpha; ONV beta_J = configuration_J.onv_beta; @@ -139,7 +139,7 @@ TwoRDMs SelectedRDMBuilder::calculate2RDMs(const Eigen::VectorXd& x) { for (size_t I = 0; I < dim; I++) { // loop over all addresses I - Configuration configuration_I = this->fock_space.get_Configuration(I); + Configuration configuration_I = this->fock_space.get_configuration(I); ONV alpha_I = configuration_I.onv_alpha; ONV beta_I = configuration_I.onv_beta; @@ -184,7 +184,7 @@ TwoRDMs SelectedRDMBuilder::calculate2RDMs(const Eigen::VectorXd& x) { for (size_t J = I+1; J < dim; J++) { - Configuration configuration_J = this->fock_space.get_Configuration(J); + Configuration configuration_J = this->fock_space.get_configuration(J); ONV alpha_J = configuration_J.onv_alpha; ONV beta_J = configuration_J.onv_beta; From 6586f50ae345f238bcdb7baf428596a031c8e661 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Tue, 23 Oct 2018 14:28:23 +0200 Subject: [PATCH 12/13] switched cmake order --- cmake/SetCMakeEnvironment.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/SetCMakeEnvironment.cmake b/cmake/SetCMakeEnvironment.cmake index 23d0dffb89..c12dc68392 100755 --- a/cmake/SetCMakeEnvironment.cmake +++ b/cmake/SetCMakeEnvironment.cmake @@ -145,8 +145,8 @@ set(PROJECT_TEST_SOURCE_FILES ${PROJECT_TESTS_FOLDER}/Operator/TwoElectronOperator_test.cpp ${PROJECT_TESTS_FOLDER}/RDM/DOCIRDMBuilder_test.cpp ${PROJECT_TESTS_FOLDER}/RDM/FCIRDMBuilder_test.cpp - ${PROJECT_TESTS_FOLDER}/RDM/SelectedRDMBuilder_test.cpp ${PROJECT_TESTS_FOLDER}/RDM/RDMCalculator_test.cpp + ${PROJECT_TESTS_FOLDER}/RDM/SelectedRDMBuilder_test.cpp ${PROJECT_TESTS_FOLDER}/RHF/DIISRHFSCFSolver_test.cpp ${PROJECT_TESTS_FOLDER}/RHF/PlainRHFSCFSolver_test.cpp ${PROJECT_TESTS_FOLDER}/RHF/RHF_test.cpp From c42d0cab32c291bc182c259bc42f6d9d4d7a0db1 Mon Sep 17 00:00:00 2001 From: tmhuysen Date: Tue, 23 Oct 2018 14:41:53 +0200 Subject: [PATCH 13/13] asString --- include/FockSpace/ONV.hpp | 2 +- src/FockSpace/ONV.cpp | 4 ++-- tests/FockSpace/SelectedFockSpace_test.cpp | 16 ++++++++-------- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/include/FockSpace/ONV.hpp b/include/FockSpace/ONV.hpp index f06c092d60..00a7b12dd9 100644 --- a/include/FockSpace/ONV.hpp +++ b/include/FockSpace/ONV.hpp @@ -185,7 +185,7 @@ class ONV { /** * @return std::string containing the ONV representation */ - std::string string_representation() const; + std::string asString() const; // FRIEND CLASSES friend class FockSpace; diff --git a/src/FockSpace/ONV.cpp b/src/FockSpace/ONV.cpp index 820747c6be..894867159e 100644 --- a/src/FockSpace/ONV.cpp +++ b/src/FockSpace/ONV.cpp @@ -58,7 +58,7 @@ ONV::ONV(size_t K, size_t representation): * Overloading of operator<< for a GQCP::ONV to be used with streams */ std::ostream& operator<<(std::ostream& os, const GQCP::ONV& onv) { - return os< ONV::findOccupiedDifferences(const ONV& other) const { /** * @return std::string containing the ONV representation */ -std::string ONV::string_representation() const { +std::string ONV::asString() const { boost::dynamic_bitset<> transfer_set (this->K, this->unsigned_representation); std::string buffer; boost::to_string(transfer_set, buffer); diff --git a/tests/FockSpace/SelectedFockSpace_test.cpp b/tests/FockSpace/SelectedFockSpace_test.cpp index 6996ec7cc5..dc1158f8e5 100644 --- a/tests/FockSpace/SelectedFockSpace_test.cpp +++ b/tests/FockSpace/SelectedFockSpace_test.cpp @@ -74,10 +74,10 @@ BOOST_AUTO_TEST_CASE ( addConfiguration ) { GQCP::Configuration configuration2 = fock_space.get_configuration(1); // Retrieve the string representation of the ONVs - std::string alpha1_test = configuration1.onv_alpha.string_representation(); - std::string alpha2_test = configuration2.onv_alpha.string_representation(); - std::string beta1_test = configuration1.onv_beta.string_representation(); - std::string beta2_test = configuration2.onv_beta.string_representation(); + std::string alpha1_test = configuration1.onv_alpha.asString(); + std::string alpha2_test = configuration2.onv_alpha.asString(); + std::string beta1_test = configuration1.onv_beta.asString(); + std::string beta2_test = configuration2.onv_beta.asString(); BOOST_CHECK(alpha1_test == alpha1_ref); BOOST_CHECK(alpha2_test == alpha2_ref); @@ -112,10 +112,10 @@ BOOST_AUTO_TEST_CASE ( reader_test ) { GQCP::Configuration configuration2 = test_reader.get_fock_space().get_configuration(1); // Retrieve the string representation of the ONVs - std::string alpha1_test = configuration1.onv_alpha.string_representation(); - std::string alpha2_test = configuration2.onv_alpha.string_representation(); - std::string beta1_test = configuration1.onv_beta.string_representation(); - std::string beta2_test = configuration2.onv_beta.string_representation(); + std::string alpha1_test = configuration1.onv_alpha.asString(); + std::string alpha2_test = configuration2.onv_alpha.asString(); + std::string beta1_test = configuration1.onv_beta.asString(); + std::string beta2_test = configuration2.onv_beta.asString(); BOOST_CHECK(alpha1_test == alpha1_ref); BOOST_CHECK(alpha2_test == alpha2_ref);