fix some bugs in MD simulation function and edit ORCA interface

example/input/orca_dyn.json

@@ -5,10 +5,9 @@
     },
     "orca" : {
         "interface" : "../../interface/orca.sh",
-        "method" : "hf",
         "dynamics" : {
             "iters" : 100,
-            "step" : 1
+            "step" : 20
         }
     }
 }

example/input/rfci_dyn.json

@@ -7,7 +7,7 @@
     "rci" : {
         "dynamics" : {
             "iters" : 100,
-            "step" : 1
+            "step" : 20
         },
         "gradient" : {
             "step" : 1e-5

example/input/rhf_dyn.json

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

example/input/rmp_dyn.json

@@ -7,7 +7,7 @@
     "rmp" : {
         "dynamics" : {
             "iters" : 100,
-            "step" : 1
+            "step" : 20
         },
         "gradient" : {
             "step" : 1e-5

include/constant.h

@@ -2,6 +2,8 @@
 
 #include <iostream>
 
+#define AMU 1822.888486192
+#define AU2K 315774.67
 #define A2BOHR 1.889726124626
 #define AU2FS 0.02418884254
 #define BOHR2A 0.529177210903

include/default.h

@@ -45,5 +45,5 @@ inline std::string rhfoptstr = R"({
 })";
 
 inline std::string orcoptstr = R"({
-    "interface" : "orca.sh", "method" : "hf", "dynamics" : {"iters" : 100, "step" : 1}
+    "interface" : "orca.sh", "dynamics" : {"iters" : 100, "step" : 1}
 })";

interface/orca.sh

@@ -1,17 +1,14 @@
 #!/bin/bash
 
-# create the input file
+# ARGUMENTS: $1 = charge, $2 = multiplicity, $3 = basis
+
+# specify the input file
 cat << EOT > orca.inp
-! ${4^^} ENGRAD ${3^^}
+! HF ENGRAD ${3^^}
 
 *xyzfile $1 $2 molecule.xyz
 EOT
 
-# numerical gradients for some methods
-if [[ ${4^^} == "CISD" ]] || [[ ${4^^} == "CISD(T)" ]] || [[ ${4^^} == "CCSD" ]] || [[ ${4^^} == "CCSD(T)" ]]; then
-    sed -i 's/ENGRAD/ENGRAD NUMGRAD/' orca.inp
-fi
-
 # run the calculation
 orca orca.inp > orca.out
 

src/main.cpp

@@ -46,7 +46,7 @@ NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedHartreeFock::Options, dynamics, gra
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedConfigurationInteraction::Options, dynamics, gradient, hessian);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(ModelSolver::OptionsNonadiabatic, dynamics, step, iters, guess, folder);
 NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(RestrictedMollerPlesset::Options, dynamics, gradient, hessian, order);
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Orca::Options, dynamics, interface, method, folder);
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Orca::Options, dynamics, interface, folder);
 
 int main(int argc, char** argv) {
     // create the argument parser and the program timer

src/method.cpp

@@ -6,21 +6,21 @@
 template <class M>
 void Method<M>::dynamics(System system, Integrals ints, Result res, bool print) const {
     // print the header
-    if (print) std::printf(" ITER  TIME [fs]        E [Eh]              KIN [Eh]              TK [K]         |GRAD|      TIME\n");
+    if (print) std::printf(" ITER  TIME [fs]    POT [Eh]       KIN [Eh]        E [Eh]         TK [K]      |GRAD|      TIME\n");
 
     // create position, velocity and mass matrices
     Matrix<> v(system.getAtoms().size(), 3), a(system.getAtoms().size(), 3), m(system.getAtoms().size(), 3);
 
     // fill the mass matrix
     for(size_t i = 0; i < system.getAtoms().size(); i++) {
-        m.row(i) = [](double m) {return Vector<>::Constant(3, m);}(an2m.at(system.getAtoms().at(i).atomic_number));
+        m.row(i) = [](double m) {return Vector<>::Constant(3, m);}(an2m.at(system.getAtoms().at(i).atomic_number) * AMU);
     }
 
     // get the degrees of freedom and write the initial geometry
     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();
+    // get the dynamics step and start the timer
+    double step = get()->opt.dynamics.step; auto start = Timer::Now();
 
     // calculate the initial energy with gradient
     if constexpr (std::is_same_v<M, Orca>) {
@@ -38,20 +38,20 @@ void Method<M>::dynamics(System system, Integrals ints, Result res, bool print)
     double Ekin = 0.5 * (m.array() * v.array() * v.array()).sum(); double T = 2 * Ekin / Nf;
 
     // print the zeroth iteration
-    if (print) std::printf("%6d %9.4f %20.14f %20.14f %20.14f %.2e %s\n", 0, 0.0, res.Etot, Ekin, T / BOLTZMANN, res.G.norm(), Timer::Format(Timer::Elapsed(start)).c_str());
+    if (print) std::printf("%6d %9.4f %14.8f %14.8f %14.8f %14.8f %.2e %s\n", 0, 0.0, res.Etot, Ekin, res.Etot + Ekin, T * AU2K, res.G.norm(), Timer::Format(Timer::Elapsed(start)).c_str());
 
     for (int i = 0; i < get()->opt.dynamics.iters; i++) {
         // start the timer and store the previous v and a
         start = Timer::Now(); Matrix<> vp = v, ap = a;
 
         // calculate the velocity and accceleration
-        a = -res.G.array() / m.array(); v = vp + (ap + a) * get()->opt.dynamics.step / 2;
+        a = -res.G.array() / m.array(); v = vp + 0.5 * (ap + a) * step;
 
         // calculate the kinetic energy and temperature
         Ekin = 0.5 * (m.array() * v.array() * v.array()).sum(); T = 2 * Ekin / Nf;
 
         // move the system
-        system.move(get()->opt.dynamics.step * (v + 0.5 * a * get()->opt.dynamics.step));
+        system.move(step * (v + 0.5 * a * step));
 
         // write the current geometry
         system.save(get()->opt.dynamics.folder / std::filesystem::path("trajectory.xyz"), std::ios::app);
@@ -69,9 +69,8 @@ void Method<M>::dynamics(System system, Integrals ints, Result res, bool print)
         }
 
         // print the iteration info
-        if (print) std::printf("%6d %9.4f %20.14f %20.14f %20.14f %.2e %s\n", i + 1, AU2FS * get()->opt.dynamics.step * (i + 1), res.Etot, Ekin, T / BOLTZMANN, res.G.norm(), Timer::Format(Timer::Elapsed(start)).c_str());
+        if (print) std::printf("%6d %9.4f %14.8f %14.8f %14.8f %14.8f %.2e %s\n", i + 1, AU2FS * step * (i + 1), res.Etot, Ekin, res.Etot + Ekin, T * AU2K, res.G.norm(), Timer::Format(Timer::Elapsed(start)).c_str());
     }
-
 }
 
 template <class M>

src/orca.cpp

@@ -5,7 +5,7 @@ Result Orca::gradient(const System& system, const Integrals&, Result res, bool)
     std::filesystem::path orcadir = opt.folder / (".orca." + std::to_string(Timer::Now().time_since_epoch().count())); res.G = Matrix<>(system.getAtoms().size(), 3);
 
     // define the execution command and create the execution directory
-    std::stringstream cmd; cmd << "./orca.sh " << system.getCharge() << " " << system.getMulti() << " " << system.getBasis() << " \"" << opt.method << "\"", std::filesystem::create_directory(orcadir);
+    std::stringstream cmd; cmd << "./orca.sh " << system.getCharge() << " " << system.getMulti() << " " << system.getBasis(), std::filesystem::create_directory(orcadir);
 
     // save the system and copy the interface
     system.save(orcadir / "molecule.xyz"), std::filesystem::copy_file(opt.interface, orcadir / "orca.sh", std::filesystem::copy_options::overwrite_existing);