From 587e985410b2072049cafcbf3f60487b117880ff Mon Sep 17 00:00:00 2001 From: Tomas Jira Date: Mon, 12 Feb 2024 10:56:36 +0100 Subject: [PATCH] implement Berendsen thermostat --- TODO.md | 2 ++ example/input/orca_dyn_thermo.json | 16 +++++++++++++++ example/input/rfci_dyn_thermo.json | 19 ++++++++++++++++++ example/input/rhf_dyn.json | 4 +--- example/input/rhf_dyn_thermo.json | 16 +++++++++++++++ example/input/rmp_dyn_thermo.json | 19 ++++++++++++++++++ example/input/uhf_dyn_thermo.json | 20 +++++++++++++++++++ include/default.h | 16 ++++++++------- include/modelsolver.h | 14 +++++++++---- include/orca.h | 7 +++++-- include/restrictedconfigurationinteraction.h | 9 +++++++-- include/restrictedhartreefock.h | 9 +++++++-- include/restrictedmollerplesset.h | 9 +++++++-- include/unrestrictedhartreefock.h | 9 +++++++-- src/main.cpp | 21 +++++++++++++------- src/method.cpp | 10 +++++++++- 16 files changed, 168 insertions(+), 32 deletions(-) create mode 100644 example/input/orca_dyn_thermo.json create mode 100644 example/input/rfci_dyn_thermo.json create mode 100644 example/input/rhf_dyn_thermo.json create mode 100644 example/input/rmp_dyn_thermo.json create mode 100644 example/input/uhf_dyn_thermo.json diff --git a/TODO.md b/TODO.md index 43bb4ed..a591ddb 100644 --- a/TODO.md +++ b/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 diff --git a/example/input/orca_dyn_thermo.json b/example/input/orca_dyn_thermo.json new file mode 100644 index 0000000..730dcf7 --- /dev/null +++ b/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 + } + } + } +} diff --git a/example/input/rfci_dyn_thermo.json b/example/input/rfci_dyn_thermo.json new file mode 100644 index 0000000..c12c2c7 --- /dev/null +++ b/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 + } + } +} diff --git a/example/input/rhf_dyn.json b/example/input/rhf_dyn.json index f4bdc63..8f760d5 100644 --- a/example/input/rhf_dyn.json +++ b/example/input/rhf_dyn.json @@ -8,8 +8,6 @@ "iters" : 100, "step" : 20 }, - "gradient" : { - "step" : 1e-5 - } + "gradient" : {} } } diff --git a/example/input/rhf_dyn_thermo.json b/example/input/rhf_dyn_thermo.json new file mode 100644 index 0000000..9c27f21 --- /dev/null +++ b/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" : {} + } +} diff --git a/example/input/rmp_dyn_thermo.json b/example/input/rmp_dyn_thermo.json new file mode 100644 index 0000000..d92ba7f --- /dev/null +++ b/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 + } + } +} diff --git a/example/input/uhf_dyn_thermo.json b/example/input/uhf_dyn_thermo.json new file mode 100644 index 0000000..be3cd22 --- /dev/null +++ b/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 + } + } +} diff --git a/include/default.h b/include/default.h index 3c9af30..7299167 100644 --- a/include/default.h +++ b/include/default.h @@ -3,7 +3,7 @@ #include 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} })"; diff --git a/include/modelsolver.h b/include/modelsolver.h index 358ee2f..39b086d 100644 --- a/include/modelsolver.h +++ b/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 guess; std::string folder; double step=0.1; int iters=1000; diff --git a/include/orca.h b/include/orca.h index f375b73..6b80ae0 100644 --- a/include/orca.h +++ b/include/orca.h @@ -6,8 +6,11 @@ class Orca : public Method { friend class Method; 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; diff --git a/include/restrictedconfigurationinteraction.h b/include/restrictedconfigurationinteraction.h index b0c9f8b..cbf4992 100644 --- a/include/restrictedconfigurationinteraction.h +++ b/include/restrictedconfigurationinteraction.h @@ -7,9 +7,14 @@ class RestrictedConfigurationInteraction : public RestrictedMethod; 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 diff --git a/include/restrictedhartreefock.h b/include/restrictedhartreefock.h index b07dce4..05a6079 100644 --- a/include/restrictedhartreefock.h +++ b/include/restrictedhartreefock.h @@ -6,9 +6,14 @@ class RestrictedHartreeFock : public RestrictedMethod { friend class Method; 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; diff --git a/include/restrictedmollerplesset.h b/include/restrictedmollerplesset.h index cdbee25..93e2fa0 100644 --- a/include/restrictedmollerplesset.h +++ b/include/restrictedmollerplesset.h @@ -7,9 +7,14 @@ class RestrictedMollerPlesset : public RestrictedMethod friend class Method; 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; diff --git a/include/unrestrictedhartreefock.h b/include/unrestrictedhartreefock.h index 25e6f94..ac4f634 100644 --- a/include/unrestrictedhartreefock.h +++ b/include/unrestrictedhartreefock.h @@ -6,9 +6,14 @@ class UnrestrictedHartreeFock : public UnrestrictedMethod; 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; diff --git a/src/main.cpp b/src/main.cpp index 59a15e5..8113cf6 100644 --- a/src/main.cpp +++ b/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); diff --git a/src/method.cpp b/src/method.cpp index bf87f48..8411aa6 100644 --- a/src/method.cpp +++ b/src/method.cpp @@ -48,7 +48,15 @@ void Method::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