Feature/selected fock space

cmake/SetCMakeEnvironment.cmake

@@ -31,6 +31,8 @@ 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
         ${PROJECT_SOURCE_FOLDER}/HamiltonianBuilder/HamiltonianBuilder.cpp
@@ -47,11 +49,14 @@ set(PROJECT_SOURCE_FILES
         ${PROJECT_SOURCE_FOLDER}/RDM/OneRDM.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
         ${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
@@ -60,9 +65,7 @@ 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
-        ${PROJECT_SOURCE_FOLDER}/WaveFunction.cpp)
+        ${PROJECT_SOURCE_FOLDER}/RMP2.cpp)
 
 # Find the header folder
 set(PROJECT_INCLUDE_FOLDER ${CMAKE_SOURCE_DIR}/include)
@@ -75,9 +78,12 @@ 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
         ${PROJECT_INCLUDE_FOLDER}/HamiltonianBuilder/HamiltonianBuilder.hpp
@@ -92,13 +98,15 @@ 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
         ${PROJECT_INCLUDE_FOLDER}/RHF/DIISRHFSCFSolver.hpp
         ${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
@@ -108,10 +116,8 @@ set(PROJECT_INCLUDE_FILES
         ${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
-        ${PROJECT_INCLUDE_FOLDER}/WaveFunction.hpp)
+        ${PROJECT_INCLUDE_FOLDER}/units.hpp)
 
 # Find the tests folder
 set(PROJECT_TESTS_FOLDER ${CMAKE_SOURCE_DIR}/tests)
@@ -129,6 +135,8 @@ 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
         ${PROJECT_TESTS_FOLDER}/HamiltonianParameters/HamiltonianParameters_test.cpp
@@ -138,6 +146,7 @@ set(PROJECT_TEST_SOURCE_FILES
         ${PROJECT_TESTS_FOLDER}/RDM/DOCIRDMBuilder_test.cpp
         ${PROJECT_TESTS_FOLDER}/RDM/FCIRDMBuilder_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
@@ -148,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)

include/CISolver/CISolver.hpp

@@ -19,9 +19,10 @@
 #define GQCP_CISOLVER_HPP
 
 
+
 #include "HamiltonianBuilder/HamiltonianBuilder.hpp"
 #include "HamiltonianParameters/HamiltonianParameters.hpp"
-#include "WaveFunction.hpp"
+#include "WaveFunction/WaveFunction.hpp"
 
 #include <numopt.hpp>
 

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 <numopt.hpp>
 

include/FockSpace/BaseFockSpace.hpp

@@ -36,17 +36,17 @@ namespace GQCP {
  */
 class BaseFockSpace {
 protected:
-    const size_t K;  // number of spatial orbitals
-    const size_t dim;  // dimension of the Fock space
+    size_t K;  // number of spatial orbitals
+    size_t dim;  // dimension of the Fock space
 
 
     // PROTECTED CONSTRUCTORS
+    BaseFockSpace() = default;
     /**
      *  Protected constructor given a @param K and @param dim
      */
     explicit BaseFockSpace(size_t K, size_t dim);
 
-
 public:
     // DESTRUCTOR
     /**

include/FockSpace/Configuration.hpp

@@ -0,0 +1,24 @@
+#ifndef GQCP_CONFIGURATION_HPP
+#define GQCP_CONFIGURATION_HPP
+
+
+#include "ONV.hpp"
+
+
+namespace GQCP {
+
+
+/**
+ *  Struct containing an electron distribution
+ *  represented by the alpha and beta ONV.
+ */
+struct Configuration {
+    ONV onv_alpha;
+    ONV onv_beta;
+};
+
+
+}  // namespace GQCP
+
+
+#endif  // GQCP_CONFIGURATION_HPP

include/FockSpace/FockSpaceType.hpp

@@ -9,7 +9,8 @@ namespace GQCP {
  */
 enum class FockSpaceType {
     FockSpace,
-    FockSpaceProduct
+    FockSpaceProduct,
+    SelectedFockSpace
 };
 
 

include/ONV.hpp → include/FockSpace/ONV.hpp

@@ -33,15 +33,15 @@ 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).
  */
 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,11 +50,22 @@ class ONV {
 
 public:
     // 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
+     *  @param K a given number of orbitals
+     *  @param unsigned_representation a representation for the ONV
+     */
+    ONV(size_t K, size_t unsigned_representation);
+
 
     // OPERATORS
     /**
@@ -97,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;
 
@@ -138,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.
      */
@@ -149,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"
@@ -159,6 +170,23 @@ 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 (from right to left) that are occupied in this, but unoccupied in @param other
+     */
+    std::vector<size_t> findOccupiedDifferences(const ONV& other) const;
+
+
+    /**
+     * @return std::string containing the ONV representation
+     */
+    std::string asString() const;
+
     // FRIEND CLASSES
     friend class FockSpace;
 };

include/FockSpace/SelectedFockSpace.hpp

@@ -0,0 +1,103 @@
+// 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 <http://www.gnu.org/licenses/>.
+// 
+#ifndef GQCP_SELECTEDFOCKSPACE_HPP
+#define GQCP_SELECTEDFOCKSPACE_HPP
+
+
+#include "FockSpace/BaseFockSpace.hpp"
+#include "FockSpace/FockSpaceProduct.hpp"
+#include "Configuration.hpp"
+
+#include <boost/numeric/conversion/converter.hpp>
+#include <boost/math/special_functions.hpp>
+
+
+namespace GQCP {
+
+
+/**
+ *  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<GQCP::Configuration> configurations;
+
+    /**
+     *  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(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);
+
+    /**
+     * 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->configurations[index]; }
+    FockSpaceType get_type() const override { return FockSpaceType::SelectedFockSpace; }
+
+
+    // PUBLIC METHODS
+    /**
+     *  Member taking two string arguments to add a Configuration
+     *  @see makeConfiguration()
+     *  add a Configuration to @var configurations
+     */
+    void addConfiguration(const std::string& onv1, const std::string& onv2);
+
+    /**
+     *  Member taking two vector<string> arguments to add Configurations
+     *  @see makeConfiguration()
+     *  add multiple Configurations to @var configurations
+     */
+    void addConfiguration(const std::vector<std::string>& onv1s, const std::vector<std::string>& onv2s);
+};
+
+
+}  // namespace GQCP
+
+
+#endif  // GQCP_SELECTEDFOCKSPACE_HPP

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 <memory>
 
@@ -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