implement Berendsen thermostat

TODO.md

@@ -2,12 +2,14 @@
 
 - [ ] Add BAGEL interface with excited state dynamics
 - [ ] Add Mulliken analysis
+- [ ] Add excitation specifications to FCI
 - [x] Add some tests
 - [x] Add the ability to export the atomic integrals
 - [x] Add the ability to use orca to perform ground state dynamics
 - [x] Check the correct working of the library before merging the fci branch
 - [x] Fix gradient calculation
 - [x] Fix the shells variable in system class that reads the basis each time the system is created
+- [x] Implement Berendsen thermostat
 - [x] Implement MP2 method
 - [x] Implement passing of the density matrix to HF method
 - [x] Implement the analytical gradient for RHF

example/input/orca_dyn_thermo.json

@@ -0,0 +1,16 @@
+{
+    "molecule" : {
+        "file" : "../molecule/water.xyz",
+        "basis" : "sto-3g"
+    },
+    "orca" : {
+        "interface" : "../../interface/orca.sh",
+        "dynamics" : {
+            "iters" : 100,
+            "step" : 20,
+            "berendsen" : {
+                "temp" : 298.15
+            }
+        }
+    }
+}

example/input/rfci_dyn_thermo.json

@@ -0,0 +1,19 @@
+{
+    "molecule" : {
+        "file" : "../molecule/water.xyz",
+        "basis" : "sto-3g"
+    },
+    "rhf" : {},
+    "rci" : {
+        "dynamics" : {
+            "iters" : 100,
+            "step" : 20,
+            "berendsen" : {
+                "temp" : 298.15
+            }
+        },
+        "gradient" : {
+            "step" : 1e-5
+        }
+    }
+}

example/input/rhf_dyn.json

@@ -8,8 +8,6 @@
             "iters" : 100,
             "step" : 20
         },
-        "gradient" : {
-            "step" : 1e-5
-        }
+        "gradient" : {}
     }
 }

example/input/rhf_dyn_thermo.json

@@ -0,0 +1,16 @@
+{
+    "molecule" : {
+        "file" : "../molecule/water.xyz",
+        "basis" : "sto-3g"
+    },
+    "rhf" : {
+        "dynamics" : {
+            "iters" : 100,
+            "step" : 20,
+            "berendsen" : {
+                "temp" : 298.15
+            }
+        },
+        "gradient" : {}
+    }
+}

example/input/rmp_dyn_thermo.json

@@ -0,0 +1,19 @@
+{
+    "molecule" : {
+        "file" : "../molecule/water.xyz",
+        "basis" : "sto-3g"
+    },
+    "rhf" : {},
+    "rmp" : {
+        "dynamics" : {
+            "iters" : 100,
+            "step" : 20,
+            "berendsen" : {
+                "temp" : 298.15
+            }
+        },
+        "gradient" : {
+            "step" : 1e-5
+        }
+    }
+}

example/input/uhf_dyn_thermo.json

@@ -0,0 +1,20 @@
+{
+    "molecule" : {
+        "file" : "../molecule/water.xyz",
+        "basis" : "sto-3g",
+        "charge" : 1,
+        "multiplicity" : 2
+    },
+    "uhf" : {
+        "dynamics" : {
+            "iters" : 100,
+            "step" : 20,
+            "berendsen" : {
+                "temp" : 298.15
+            }
+        },
+        "gradient" : {
+            "step" : 1e-5
+        }
+    }
+}

include/default.h

@@ -3,7 +3,7 @@
 #include <string>
 
 inline std::string rcioptstr = R"({
-    "dynamics" : {"iters" : 100, "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
+    "dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}, "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}
 })";
@@ -23,34 +23,36 @@ inline std::string moloptstr = R"({
 
 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, "step" : 1}
+    "dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}
 })";
 
 inline std::string msnoptstr = R"({
-    "step" : 0.1, "guess" : ["-x^2"], "iters" : 1000, "dynamics" : {"iters" : 100, "step" : 1}
+    "step" : 0.1, "guess" : ["-x^2"], "iters" : 1000,
+    "dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}
 })";
 
 inline std::string rmpoptstr = R"({
     "order" : 2,
-    "dynamics" : {"iters" : 100, "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
+    "dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
     "export" : {"coulombmo" : false},
     "print" : {"coulombmo" : false}
 })";
 
 inline std::string rhfoptstr = R"({
     "maxiter" : 1000, "thresh" : 1e-12,
-    "dynamics" : {"iters" : 100, "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
+    "dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}, "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", "dynamics" : {"iters" : 100, "step" : 1}
+    "interface" : "orca.sh",
+    "dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}
 })";
 
 inline std::string uhfoptstr = R"({
     "maxiter" : 1000, "thresh" : 1e-12,
-    "dynamics" : {"iters" : 100, "step" : 1}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
+    "dynamics" : {"iters" : 100, "step" : 20, "berendsen" : {"tau" : 15, "temp" : 0, "timeout" : 20}}, "gradient" : {"step" : 1e-5}, "hessian" : {"step" : 1e-5},
     "export" : {"coefa" : false, "coefb" : false, "densitya" : false, "densityb" : false, "orbena" : false, "orbenb" : false, "hcore" : false},
     "print" : {"coefa" : false, "coefb" : false, "densitya" : false, "densityb" : false, "orbena" : false, "orbenb" : false, "hcore" : false}
 })";

include/modelsolver.h

@@ -6,16 +6,22 @@
 class ModelSolver {
 public:
     struct OptionsAdiabatic {OptionsAdiabatic(){};
-        // structures
-        struct Dynamics {int iters=100; double step=1.0; std::string folder;} dynamics={};
+        // dynamics structure
+        struct Dynamics {
+            struct Berendsen {double tau=1, temp=0, timeout=10;} berendsen={};
+            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 folder;} dynamics={};
+        // dynamics structure
+        struct Dynamics {
+            struct Berendsen {double tau=1, temp=0, timeout=10;} berendsen={};
+            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

@@ -6,8 +6,11 @@ class Orca : public Method<Orca> {
     friend class Method<Orca>;
 public:
     struct Options {Options(){};
-        // structures
-        struct Dynamics {int iters=100; double step=1.0; std::string folder;} dynamics={};
+        // dynamics structure
+        struct Dynamics {
+            struct Berendsen {bool enabled=false; double tau = 1, temp=298.15, timeout=10;} berendsen={};
+            int iters=100; double step=1.0; std::string folder;
+        } dynamics={};
 
         // variables
         std::string method; std::filesystem::path interface, folder;

include/restrictedconfigurationinteraction.h

@@ -7,9 +7,14 @@ class RestrictedConfigurationInteraction : public RestrictedMethod<RestrictedCon
     friend class Method<RestrictedConfigurationInteraction>;
 public:
     struct Options {Options(){};
-        // structure of the options
+        // gradient and hessian structures
         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;
+
+        // dynamics structure
+        struct Dynamics {
+            struct Berendsen {double tau=1, temp=0, timeout=10;} berendsen={};
+            int iters=100; double step=1.0; std::string folder;
+        } dynamics={};
     };
 public:
     // constructors and destructors

include/restrictedhartreefock.h

@@ -6,9 +6,14 @@ class RestrictedHartreeFock : public RestrictedMethod<RestrictedHartreeFock> {
     friend class Method<RestrictedHartreeFock>;
 public:
     struct Options {Options(){};
-        // structures
+        // gradient and hessian structrures
         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={};
+
+        // dynamics structure
+        struct Dynamics {
+            struct Berendsen {double tau=1, temp=0, timeout=10;} berendsen={};
+            int iters=100; double step=1.0; std::string folder;
+        } dynamics={};
 
         // variables
         int maxiter=100; double thresh=1e-8;

include/restrictedmollerplesset.h

@@ -7,9 +7,14 @@ class RestrictedMollerPlesset : public RestrictedMethod<RestrictedMollerPlesset>
     friend class Method<RestrictedMollerPlesset>;
 public:
     struct Options {Options(){};
-        // structure of the options
+        // gradient and hessian structures
         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;
+
+        // dynamics structure
+        struct Dynamics {
+            struct Berendsen {double tau=1, temp=0, timeout=10;} berendsen={};
+            int iters=100; double step=1.0; std::string folder;
+        } dynamics={};
 
         // values of simple options
         int order=2;

include/unrestrictedhartreefock.h

@@ -6,9 +6,14 @@ class UnrestrictedHartreeFock : public UnrestrictedMethod<UnrestrictedHartreeFoc
     friend class Method<UnrestrictedHartreeFock>;
 public:
     struct Options {Options(){};
-        // structures
+        // gradient and hessian structrures
         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={};
+
+        // dynamics structure
+        struct Dynamics {
+            struct Berendsen {double tau=1, temp=0, timeout=10;} berendsen={};
+            int iters=100; double step=1.0; std::string folder;
+        } dynamics={};
 
         // variables
         int maxiter=100; double thresh=1e-8;

src/main.cpp

@@ -20,31 +20,38 @@
 #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, folder);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedHartreeFock::Options::Dynamics::Berendsen, tau, temp, timeout);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedHartreeFock::Options::Dynamics, iters, step, folder, berendsen);
 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, folder);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedConfigurationInteraction::Options::Dynamics::Berendsen, tau, temp, timeout);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedConfigurationInteraction::Options::Dynamics, iters, step, folder, berendsen);
 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, folder);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options::Dynamics::Berendsen, tau, temp, timeout);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options::Dynamics, iters, step, folder, berendsen);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options::Gradient, step);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options::Hessian, step);
 
 // option structures loaders for UHF
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(UnrestrictedHartreeFock::Options::Dynamics, iters, step, folder);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(UnrestrictedHartreeFock::Options::Dynamics::Berendsen, tau, temp, timeout);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(UnrestrictedHartreeFock::Options::Dynamics, iters, step, folder, berendsen);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(UnrestrictedHartreeFock::Options::Gradient, step);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(UnrestrictedHartreeFock::Options::Hessian, step);
 
 // option structures loaders for ORCA
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Orca::Options::Dynamics, iters, step, folder);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Orca::Options::Dynamics::Berendsen, tau, temp, timeout);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Orca::Options::Dynamics, iters, step, folder, berendsen);
 
 // option structures loaders for model method
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic::Dynamics, iters, step, folder);
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic::Dynamics, iters, step, folder);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic::Dynamics::Berendsen, tau, temp, timeout);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic::Dynamics, iters, step, folder, berendsen);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic::Dynamics::Berendsen, tau, temp, timeout);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic::Dynamics, iters, step, folder, berendsen);
 
 // option loaders
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsAdiabatic, dynamics, real, step, iters, nstate, thresh, optimize, guess, folder);

src/method.cpp

@@ -48,7 +48,15 @@ void Method<M>::dynamics(System system, Integrals ints, Result res, bool print)
         // calculate the velocity and accceleration
         a = -res.G.array() / m.array(); v = vp + 0.5 * (ap + a) * step;
 
-        // calculate the kinetic energy and temperature
+        // calculate the kinetic energy and temperature before thermostatting
+        Ekin = 0.5 * (m.array() * v.array() * v.array()).sum(); T = 2 * Ekin / Nf;
+
+        // apply the berendsen thermostat
+        if (get()->opt.dynamics.berendsen.temp > 0 && AU2FS * step * i < get()->opt.dynamics.berendsen.timeout) {
+            v.array() *= std::sqrt(1 + (get()->opt.dynamics.berendsen.temp / T / AU2K - 1) * step / get()->opt.dynamics.berendsen.tau);
+        }
+
+        // calculate the kinetic energy and temperature after thermostatting
         Ekin = 0.5 * (m.array() * v.array() * v.array()).sum(); T = 2 * Ekin / Nf;
 
         // move the system