add a lot of tests for RHF, MP2 and FCI

CMakeLists.txt

@@ -26,7 +26,7 @@ FetchContent_Declare(json GIT_REPOSITORY https://github.com/nlohmann/json.git GI
 FetchContent_MakeAvailable(argparse exprtk json)
 
 # include directories
-include_directories(include ${exprtk_SOURCE_DIR} ${argparse_SOURCE_DIR}/include ${json_SOURCE_DIR}/include)
+include_directories(include ${PROJECT_SOURCE_DIR}/bin ${exprtk_SOURCE_DIR} ${argparse_SOURCE_DIR}/include ${json_SOURCE_DIR}/include)
 
 # find OpenMP
 find_package(Eigen3 REQUIRED)
@@ -73,19 +73,19 @@ enable_testing()
 
 if (BUILD_TESTING)
     # add test executables
+    add_executable(test_rmp2 test/rmp2.cpp)
+    add_executable(test_rfci test/rfci.cpp)
     add_executable(test_rhf test/rhf.cpp)
-    add_executable(test_mp2 test/mp2.cpp)
-    add_executable(test_fci test/fci.cpp)
 
     # link the test executables to the library
+    target_link_libraries(test_rmp2 acorn)
+    target_link_libraries(test_rfci acorn)
     target_link_libraries(test_rhf acorn)
-    target_link_libraries(test_mp2 acorn)
-    target_link_libraries(test_fci acorn)
 
     # change the output folder for the test executables
+    set_target_properties(test_rmp2 PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/bin/test)
+    set_target_properties(test_rfci PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/bin/test)
     set_target_properties(test_rhf PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/bin/test)
-    set_target_properties(test_mp2 PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/bin/test)
-    set_target_properties(test_fci PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/bin/test)
 endif()
 
 # include test file

include/default.h

@@ -3,7 +3,7 @@
 #include <string>
 
 inline std::string rcioptstr = R"({
-    "dynamics" : {"iters" : 100, "output" : "trajectory.xyz", "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
+    "dynamics" : {"iters" : 100, "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
     "export" : {"coulombms" : false, "kineticms" : false, "nuclearms" : false, "hcorems" : false, "hamiltonian" : false, "energies" : false},
     "print" : {"coulombms" : false, "kineticms" : false, "nuclearms" : false, "hcorems" : false, "hamiltonian" : false, "energies" : false}
 })";
@@ -17,29 +17,33 @@ inline std::string mdloptstr = R"({
     "mass" : 1, "ngrid" : 512, "limits" : [[-16, 16]]
 })";
 
+inline std::string moloptstr = R"({
+    "charge" : 0, "multiplicity" : 1
+})";
+
 inline std::string msaoptstr = R"({
     "real" : false, "step" : 0.1, "optimize" : true, "guess" : "-x^2", "nstate" : 3, "iters" : 1000, "thresh" : 1e-12,
-    "dynamics" : {"iters" : 100, "output" : "trajectory.xyz", "step" : 1}
+    "dynamics" : {"iters" : 100, "step" : 1}
 })";
 
 inline std::string msnoptstr = R"({
-    "step" : 0.1, "guess" : ["-x^2"], "iters" : 1000, "dynamics" : {"iters" : 100, "output" : "trajectory.xyz", "step" : 1}
+    "step" : 0.1, "guess" : ["-x^2"], "iters" : 1000, "dynamics" : {"iters" : 100, "step" : 1}
 })";
 
 inline std::string rmpoptstr = R"({
     "order" : 2,
-    "dynamics" : {"iters" : 100, "output" : "trajectory.xyz", "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
+    "dynamics" : {"iters" : 100, "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
     "export" : {"coulombmo" : false},
     "print" : {"coulombmo" : false}
 })";
 
 inline std::string rhfoptstr = R"({
     "maxiter" : 1000, "thresh" : 1e-8,
-    "dynamics" : {"iters" : 100, "output" : "trajectory.xyz", "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
+    "dynamics" : {"iters" : 100, "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
     "export" : {"coef" : false, "density" : false, "orben" : false, "hcore" : false},
     "print" : {"coef" : false, "density" : false, "orben" : false, "hcore" : false}
 })";
 
 inline std::string orcoptstr = R"({
-    "interface" : "orca.sh", "method" : "hf", "dynamics" : {"iters" : 100, "output" : "trajectory.xyz", "step" : 1}
+    "interface" : "orca.sh", "method" : "hf", "dynamics" : {"iters" : 100, "step" : 1}
 })";

include/modelsolver.h

@@ -7,15 +7,15 @@ class ModelSolver {
 public:
     struct OptionsAdiabatic {OptionsAdiabatic(){};
         // structures
-        struct Dynamics {int iters=100; double step=1.0; std::string output="trajectory.xyz";} dynamics={};
+        struct Dynamics {int iters=100; double step=1.0; std::string folder;} dynamics={};
 
         // variables
         int nstate=3, iters=1000; bool real=false, optimize=false;
         std::string folder, guess; double step=0.1, thresh=1e-8;
     };
     struct OptionsNonadiabatic {OptionsNonadiabatic(){};
         // structures
-        struct Dynamics {int iters=100; double step=1.0; std::string output="trajectory.xyz";} dynamics={};
+        struct Dynamics {int iters=100; double step=1.0; std::string folder;} dynamics={};
 
         // variables
         std::vector<std::string> guess; std::string folder; double step=0.1; int iters=1000;

include/orca.h

@@ -7,7 +7,7 @@ class Orca : public Method<Orca> {
 public:
     struct Options {Options(){};
         // structures
-        struct Dynamics {int iters=100; double step=1.0; std::string output="trajectory.xyz";} dynamics={};
+        struct Dynamics {int iters=100; double step=1.0; std::string folder;} dynamics={};
 
         // variables
         std::string method; std::filesystem::path interface, folder;

include/restrictedconfigurationinteraction.h

@@ -8,8 +8,8 @@ class RestrictedConfigurationInteraction : public RestrictedMethod<RestrictedCon
 public:
     struct Options {Options(){};
         // structure of the options
-        struct Dynamics {int iters=100; double step=1.0; std::string output="trajectory.xyz";} dynamics;
         struct Gradient {double step=1e-5;} gradient; struct Hessian {double step=1e-5;} hessian;
+        struct Dynamics {int iters=100; double step=1.0; std::string folder;} dynamics;
     };
 public:
     // constructors and destructors

include/restrictedhartreefock.h

@@ -7,8 +7,8 @@ class RestrictedHartreeFock : public RestrictedMethod<RestrictedHartreeFock> {
 public:
     struct Options {Options(){};
         // structures
-        struct Dynamics {int iters=100; double step=1.0; std::string output="trajectory.xyz";} dynamics={};
         struct Gradient {double step=1e-5;} gradient={}; struct Hessian {double step=1e-5;} hessian={};
+        struct Dynamics {int iters=100; double step=1.0; std::string folder;} dynamics={};
 
         // variables
         int maxiter=100; double thresh=1e-8;

include/restrictedmollerplesset.h

@@ -8,8 +8,8 @@ class RestrictedMollerPlesset : public RestrictedMethod<RestrictedMollerPlesset>
 public:
     struct Options {Options(){};
         // structure of the options
-        struct Dynamics {int iters=100; double step=1.0; std::string output="trajectory.xyz";} dynamics;
         struct Gradient {double step=1e-5;} gradient; struct Hessian {double step=1e-5;} hessian;
+        struct Dynamics {int iters=100; double step=1.0; std::string folder;} dynamics;
 
         // values of simple options
         int order=2;

include/table.h

@@ -12,9 +12,15 @@ inline std::unordered_map<int, double> an2m = {
 
 inline std::unordered_map<int, std::string> an2sm = {
     { 1,  "H"},
+    { 2, "He"},
     { 6,  "C"},
     { 7,  "N"},
     { 8,  "O"},
     { 9,  "F"},
-    {17, "Cl"}
+    {10, "Ne"},
+    {14, "Si"},
+    {15,  "P"},
+    {16,  "S"},
+    {17, "Cl"},
+    {18, "Ar"}
 };

script/testing.sh

@@ -0,0 +1,68 @@
+#!/bin/bash
+
+read -d '' RHF << EOF
+add_test(NAME %s COMMAND \${PROJECT_SOURCE_DIR}/bin/test/test_rhf %s %s)
+set_tests_properties(%s PROPERTIES PASS_REGULAR_EXPRESSION "%s")
+EOF
+
+read -d '' RMP2 << EOF
+add_test(NAME %s COMMAND \${PROJECT_SOURCE_DIR}/bin/test/test_rmp2 %s %s)
+set_tests_properties(%s PROPERTIES PASS_REGULAR_EXPRESSION "%s")
+EOF
+
+read -d '' RFCI << EOF
+add_test(NAME %s COMMAND \${PROJECT_SOURCE_DIR}/bin/test/test_rfci %s %s)
+set_tests_properties(%s PROPERTIES PASS_REGULAR_EXPRESSION "%s")
+EOF
+
+create() {
+    # extract the variables from the function arguments
+    SYSTEM="$1"; CHARGE="$2"; MULT="$3"; BASIS="${4,,}"; METHOD="$5"; BASISFILE=$(echo "$BASIS" | sed -e 's/+/p/g ; s/*/s/g')
+
+    # create the molecule block
+    jq -n --arg file "example/molecule/$SYSTEM.xyz" --arg basis "$BASIS" --argjson charge "$CHARGE" --argjson multiplicity "$MULT" '{molecule: $ARGS.named}' > molecule.json
+
+    # create the input file
+    if [ "$METHOD" == "RHF" ]; then
+        jq '. + {rhf: {}}' molecule.json > input.json
+    elif [ "$METHOD" == "RMP2" ]; then
+        jq '. + {rhf: {}, rmp: {}}' molecule.json > input.json
+    elif [ "$METHOD" == "RFCI" ]; then
+        jq '. + {rhf: {}, rci: {}}' molecule.json > input.json
+    fi
+
+    # run the program with a timeout
+    $(timeout 1 acorn input.json > output.out 2> /dev/null)
+
+    # create the test
+    if [ "$METHOD" == "RHF" ]; then
+        EXPECT=$(grep "RESTRICTED HARTREE-FOCK ENERGY" output.out | awk '{print $4}')
+        [[ "$EXPECT" ]] && printf "\n\n# test-energy: %s %d %d %s RHF\n" "$SYSTEM" "$CHARGE" "$MULT" "${BASIS^^}"; TNAME="${SYSTEM}_${CHARGE}-${MULT}_${BASISFILE}_rhf_energy"
+        [[ "$EXPECT" ]] && printf "$RHF" "$TNAME" "$SYSTEM" "$BASIS" "$TNAME" "$(echo $EXPECT | awk '{printf "%.8f", $1}')"
+    elif [ "$METHOD" == "RMP2" ]; then
+        EXPECT=$(grep "RESTRICTED MOLLER-PLESSET ENERGY" output.out | awk '{print $4}')
+        [[ "$EXPECT" ]] && printf "\n\n# test-energy: %s %d %d %s RMP2\n" "$SYSTEM" "$CHARGE" "$MULT" "${BASIS^^}"; TNAME="${SYSTEM}_${CHARGE}-${MULT}_${BASISFILE}_rmp2_energy"
+        [[ "$EXPECT" ]] && printf "$RMP2" "$TNAME" "$SYSTEM" "$BASIS" "$TNAME" "$(echo $EXPECT | awk '{printf "%.8f", $1}')"
+    elif [ "$METHOD" == "RFCI" ]; then
+        EXPECT=$(grep "RESTRICTED CONFIGURATION INTERACTION ENERGY" output.out | awk '{print $5}')
+        [[ "$EXPECT" ]] && printf "\n\n# test-energy: %s %d %d %s RFCI\n" "$SYSTEM" "$CHARGE" "$MULT" "${BASIS^^}"; TNAME="${SYSTEM}_${CHARGE}-${MULT}_${BASISFILE}_rfci_energy"
+        [[ "$EXPECT" ]] && printf "$RFCI" "$TNAME" "$SYSTEM" "$BASIS" "$TNAME" "$(echo $EXPECT | awk '{printf "%.8f", $1}')"
+    fi
+
+    # remove the temporary files
+    rm molecule.json input.json output.out
+}
+
+printf "# this file is for definitions of tests and can be generated by the testing.sh script in the script folder"
+
+for SYSTEM in example/molecule/*.xyz; do
+for BASIS in basis/*.g94; do
+    # extract the system name and basis set name
+    SYSTEM=$(basename "$SYSTEM" .xyz); BASIS=$(basename "$BASIS" .g94 | sed 's/augmentation/aug/g')
+
+    # create the tests
+    create "$SYSTEM" 0 1 "$BASIS" "RHF"; create "$SYSTEM" 0 1 "$BASIS" "RMP2"; create "$SYSTEM" 0 1 "$BASIS" "RFCI"
+done
+done
+
+echo ""

src/main.cpp

@@ -19,26 +19,26 @@
 #define MEASURE(T, F) std::cout << T << std::flush; {auto t = Timer::Now(); F; std::printf("%s\n", Timer::Format(Timer::Elapsed(t)).c_str());}
 
 // option structures loaders for RHF
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedHartreeFock::Options::Dynamics, iters, step, output);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedHartreeFock::Options::Dynamics, iters, step, folder);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedHartreeFock::Options::Gradient, step);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedHartreeFock::Options::Hessian, step);
 
 // option structures loaders for RCI
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedConfigurationInteraction::Options::Dynamics, iters, step, output);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedConfigurationInteraction::Options::Dynamics, iters, step, folder);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedConfigurationInteraction::Options::Gradient, step);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedConfigurationInteraction::Options::Hessian, step);
 
 // option structures loaders for RMP
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options::Dynamics, iters, step, output);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options::Dynamics, iters, step, folder);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options::Gradient, step);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options::Hessian, step);
 
 // option structures loaders for ORCA
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Orca::Options::Dynamics, iters, step, output);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Orca::Options::Dynamics, iters, step, folder);
 
 // option structures loaders for model method
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic::Dynamics, iters, step, output);
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic::Dynamics, iters, step, output);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic::Dynamics, iters, step, folder);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic::Dynamics, iters, step, folder);
 
 // option loaders
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic, dynamics, real, step, iters, nstate, thresh, optimize, guess, folder);
@@ -88,11 +88,16 @@ int main(int argc, char** argv) {
     auto intopt = nlohmann::json::parse(intoptstr), rhfopt = nlohmann::json::parse(rhfoptstr);
     auto msaopt = nlohmann::json::parse(msaoptstr), mdlopt = nlohmann::json::parse(mdloptstr);
     auto msnopt = nlohmann::json::parse(msnoptstr), rciopt = nlohmann::json::parse(rcioptstr);
+    auto molopt = nlohmann::json::parse(moloptstr);
 
     // assign the output paths to some method options that output files
     orcopt["folder"] = inputpath, msaopt["folder"] = inputpath, msnopt["folder"] = inputpath;
+    rhfopt["dynamics"]["folder"] = inputpath, rciopt["dynamics"]["folder"] = inputpath;
+    msnopt["dynamics"]["folder"] = inputpath, msaopt["dynamics"]["folder"] = inputpath;
+    rmpopt["dynamics"]["folder"] = inputpath, orcopt["dynamics"]["folder"] = inputpath;
 
     // patch the input json file and apply defaults
+    if (input.contains("molecule")) molopt.merge_patch(input.at("molecule"));
     if (input.contains("integral")) intopt.merge_patch(input.at("integral"));
     if (input.contains("model")) mdlopt.merge_patch(input.at("model"));
     if (input.contains("orca")) orcopt.merge_patch(input.at("orca"));
@@ -111,15 +116,11 @@ int main(int argc, char** argv) {
         // define the integral container
         Integrals ints(true);
 
-        // make all the input and output paths absolute
-        orcopt.at("dynamics").at("output") = inputpath / orcopt.at("dynamics").at("output");
-        rhfopt.at("dynamics").at("output") = inputpath / rhfopt.at("dynamics").at("output");
-        rciopt.at("dynamics").at("output") = inputpath / rciopt.at("dynamics").at("output");
-        rmpopt.at("dynamics").at("output") = inputpath / rmpopt.at("dynamics").at("output");
+        // make all the interface paths absolute
         orcopt.at("interface") = inputpath / orcopt.at("interface");
 
         // create the system from the system file
-        std::ifstream mstream(syspath); System system(mstream, input.at("molecule").at("basis")); mstream.close();
+        std::ifstream mstream(syspath); System system(mstream, molopt.at("basis"), molopt.at("charge"), molopt.at("multiplicity")); mstream.close();
 
         // print the molecule specification
         Printer::Title("ATOM COORDINATE FILE AND BASIS"); std::cout << system << std::endl << std::endl;

src/method.cpp

@@ -17,7 +17,7 @@ void Method<M>::dynamics(System system, const Integrals& ints, Result res, bool
     }
 
     // get the degrees of freedom and write the initial geometry
-    double Nf = system.getAtoms().size() * 3 - 6; system.save(get()->opt.dynamics.output) ;
+    double Nf = system.getAtoms().size() * 3 - 6; system.save(get()->opt.dynamics.folder / std::filesystem::path("trajectory.xyz")) ;
 
     // start the timer
     auto start = Timer::Now();
@@ -49,7 +49,7 @@ void Method<M>::dynamics(System system, const Integrals& ints, Result res, bool
         system.move(get()->opt.dynamics.step * (v + 0.5 * a * get()->opt.dynamics.step));
 
         // write the current geometry
-        system.save(get()->opt.dynamics.output, std::ios::app);
+        system.save(get()->opt.dynamics.folder / std::filesystem::path("trajectory.xyz"), std::ios::app);
 
         // calculate the next energy with gradient
         if constexpr (std::is_same_v<M, Orca>) {
@@ -133,7 +133,7 @@ Result Method<M>::hessian(const System& system, const Integrals&, Result res, bo
     // print the header
     if (print) std::printf("  ELEM      dE [Eh/Bohr]        TIME\n");
 
-    #pragma omp parallel for num_threads(nthread) collapse(2)
+    #pragma omp parallel for num_threads(nthread)
     for (int i = 0; i < res.H.rows(); i++) {
         for (int j = i; j < res.H.cols(); j++) {
             // start the timer

src/system.cpp

@@ -19,6 +19,13 @@ System::System(std::ifstream& stream, std::string basis, int charge, int multi)
     // check for the input geometry file existence
     if (!stream.good()) throw std::runtime_error("SYSTEM FILE DOES NOT EXIST");
 
+    // throw an error if impossible combination of charge and multiplicity
+    if (std::abs(charge) % 2 == 0 && multi % 2 == 0) {
+        throw std::runtime_error("MOLECULE CAN'T HAVE AN EVEN CHARGE AND MULTIPLICITY AT THE SAME TIME.");
+    } else if (std::abs(charge) % 2 == 1 && multi % 2 == 1) {
+        throw std::runtime_error("MOLECULE CAN'T HAVE AN ODD CHARGE AND MULTIPLICITY AT THE SAME TIME.");
+    }
+
     // replace the basis placeholders
     std::replace(this->basis.begin(), this->basis.end(), '*', 's'), std::replace(this->basis.begin(), this->basis.end(), '+', 'p');
 

test/fci.cpp → test/rfci.cpp

@@ -1,4 +1,4 @@
-#include "../bin/acorn.h"
+#include "acorn.h"
 
 int main(int argc, char** argv) {
     // throw an error if the number of arguments is incorrect

test/rhf.cpp

@@ -1,4 +1,4 @@
-#include "../bin/acorn.h"
+#include "acorn.h"
 
 int main(int argc, char** argv) {
     // throw an error if the number of arguments is incorrect

test/mp2.cpp → test/rmp2.cpp

@@ -1,4 +1,4 @@
-#include "../bin/acorn.h"
+#include "acorn.h"
 
 int main(int argc, char** argv) {
     // throw an error if the number of arguments is incorrect