update md and xtb

Signed-off-by: Conrad Huebler <[email protected]>

Author: conradhuebler

Diff for: external/tblite

@@ -226,6 +226,9 @@ void SimpleMD::LoadControlJson()
     m_rm_COM = Json2KeyWord<double>(m_defaults, "rm_COM");
     int rattle = Json2KeyWord<int>(m_defaults, "rattle");
     m_rattle_maxiter = Json2KeyWord<int>(m_defaults, "rattle_maxiter");
+    m_rattle_dynamic_tol_iter = Json2KeyWord<int>(m_defaults, "rattle_dynamic_tol_iter");
+    m_rattle_dynamic_tol = Json2KeyWord<bool>(m_defaults, "rattle_dynamic_tol");
+
     if (rattle == 1) {
         Integrator = [=](double* grad) {
             this->Rattle(grad);
@@ -718,6 +721,11 @@ nlohmann::json SimpleMD::WriteRestartInformation()
     restart["average_Virial"] = m_average_virial_correction;
     restart["average_Wall"] = m_average_wall_potential;
 
+    restart["rattle"] = m_rattle;
+    restart["rattle_maxiter"] = m_rattle_maxiter;
+    restart["rattle_dynamic_tol"] = m_rattle_tolerance;
+    restart["rattle_dynamic_tol_iter"] = m_rattle_dynamic_tol_iter;
+
     restart["coupling"] = m_coupling;
     restart["MaxTopoDiff"] = m_max_top_diff;
     restart["impuls"] = m_impuls;
@@ -789,7 +797,7 @@ bool SimpleMD::LoadRestartInformation()
 
 bool SimpleMD::LoadRestartInformation(const json& state)
 {
-    std::string geometry, velocities, constrains, zeta, xi, Q;
+    std::string geometry, velocities, constrains, xi, Q;
 
     try {
         m_method = state["method"];
@@ -879,10 +887,6 @@ bool SimpleMD::LoadRestartInformation(const json& state)
     } catch (json::type_error& e) {
     }
 
-    try {
-        zeta = state["zeta"];
-    } catch (json::type_error& e) {
-    }
     try {
         xi = state["xi"];
     } catch (json::type_error& e) {
@@ -897,6 +901,30 @@ bool SimpleMD::LoadRestartInformation(const json& state)
     } catch (json::type_error& e) {
     }
 
+    try {
+        m_rattle = state["rattle"];
+    } catch (json::type_error& e) {
+    }
+
+    try {
+        m_rattle_tolerance = state["rattle_tolerance"];
+    } catch (json::type_error& e) {
+    }
+
+    try {
+        m_rattle_maxiter = state["rattle_maxiter"];
+    } catch (json::type_error& e) {
+    }
+
+    try {
+        m_rattle_dynamic_tol = state["rattle_dynamic_tol"];
+    } catch (json::type_error& e) {
+    }
+
+    try {
+        m_rattle_dynamic_tol_iter = state["rattle_dynamic_tol_iter"];
+    } catch (json::type_error& e) {
+    }
     try {
         m_rmsd_mtd = state["rmsd_mtd"];
         if (m_rmsd_mtd) {
@@ -1112,7 +1140,10 @@ void SimpleMD::start()
             fmt::print(fg(fmt::color::salmon) | fmt::emphasis::bold, "Simulation got unstable, exiting!\n");
 
             std::ofstream restart_file("unstable_curcuma.json");
-            restart_file << WriteRestartInformation() << std::endl;
+            nlohmann::json restart;
+            restart[MethodName()[0]] = WriteRestartInformation();
+            restart_file << restart << std::endl;
+
             m_time_step = 0;
             aborted = true;
 
@@ -1127,14 +1158,20 @@ void SimpleMD::start()
         if (m_writerestart > -1 && m_step % m_writerestart == 0) {
             std::ofstream restart_file("curcuma_step_" + std::to_string(int(m_step * m_dT)) + ".json");
             nlohmann::json restart;
-            restart_file << WriteRestartInformation() << std::endl;
+            restart[MethodName()[0]] = WriteRestartInformation();
+            restart_file << restart << std::endl;
         }
         if ((m_step && int(m_step * m_dT) % m_print == 0)) {
             m_Etot = m_Epot + m_Ekin;
             PrintStatus();
             m_time_step = 0;
         }
-
+        if (m_rattle && m_rattle_dynamic_tol) {
+            m_aver_rattle_Temp += m_T;
+            m_rattle_counter++;
+            if (m_rattle_counter == m_rattle_dynamic_tol_iter)
+                AdjustRattleTolerance();
+        }
         if (m_impuls > m_T) {
             InitVelocities(m_scale_velo * m_impuls_scaling);
             EKin();
@@ -1187,12 +1224,30 @@ void SimpleMD::start()
         }
     }
     std::ofstream restart_file("curcuma_final.json");
-    restart_file << WriteRestartInformation() << std::endl;
+    nlohmann::json restart;
+    restart[MethodName()[0]] = WriteRestartInformation();
+    restart_file << restart << std::endl;
     if (aborted == false)
         std::remove("curcuma_restart.json");
     delete[] gradient;
 }
 
+void SimpleMD::AdjustRattleTolerance()
+{
+    m_aver_rattle_Temp /= double(m_rattle_counter);
+
+    // std::pair<double, double> pair(m_rattle_tolerance, m_aver_Temp);
+
+    if (m_aver_rattle_Temp > m_T0)
+        m_rattle_tolerance -= 0.01;
+    else if (m_aver_rattle_Temp < m_T0)
+        m_rattle_tolerance += 0.01;
+    std::cout << m_rattle_counter << " " << m_aver_rattle_Temp << " " << m_rattle_tolerance << std::endl;
+    m_rattle_tolerance = std::abs(m_rattle_tolerance);
+    m_rattle_counter = 0;
+    m_aver_rattle_Temp = 0;
+}
+
 void SimpleMD::Verlet(double* grad)
 {
     double ekin = 0;
@@ -1254,7 +1309,7 @@ void SimpleMD::Verlet(double* grad)
     }
     ekin *= 0.5;
     double T = 2.0 * ekin / (kb_Eh * m_dof);
-    m_unstable = T > 10000 * m_T;
+    m_unstable = T > 10000 * m_T || std::isnan(T);
     m_T = T;
     m_Ekin = ekin;
     ThermostatFunction();
@@ -1275,6 +1330,7 @@ void SimpleMD::Rattle(double* grad)
      * like dT^3 -> dT^2 and
      * updated velocities of the second atom (minus instead of plus)
      */
+    TriggerWriteRestart();
     double* coord = new double[3 * m_natoms];
     double m_dT_inverse = 1 / m_dT;
     std::vector<int> moved(m_natoms, 0);
@@ -1301,8 +1357,7 @@ void SimpleMD::Rattle(double* grad)
                 + (coord[3 * i + 1] - coord[3 * j + 1]) * (coord[3 * i + 1] - coord[3 * j + 1])
                 + (coord[3 * i + 2] - coord[3 * j + 2]) * (coord[3 * i + 2] - coord[3 * j + 2]));
 
-            if (std::abs(distance - distance_current) > 2 * m_rattle_tolerance * distance) {
-                active++;
+            if (std::abs(distance - distance_current) > m_rattle_tolerance) {
                 move = true;
                 double r = distance - distance_current;
                 double dx = m_current_geometry[3 * i + 0] - m_current_geometry[3 * j + 0];
@@ -1315,6 +1370,7 @@ void SimpleMD::Rattle(double* grad)
                 if (scalarproduct >= m_rattle_tolerance * distance) {
                     moved[i] = 1;
                     moved[j] = 1;
+                    active++;
 
                     double lambda = r / (1 * (m_rmass[i] + m_rmass[j]) * scalarproduct);
                     while (std::abs(lambda) > max_mu)
@@ -1340,6 +1396,13 @@ void SimpleMD::Rattle(double* grad)
         if (active == 0)
             break;
     }
+    if (iter >= m_rattle_maxiter) {
+        std::cout << "numeric difficulties - 1st step in rattle velocity verlet" << std::endl;
+        std::ofstream restart_file("unstable_curcuma_" + std::to_string(m_currentStep) + ".json");
+        nlohmann::json restart;
+        restart[MethodName()[0]] = WriteRestartInformation();
+        restart_file << restart << std::endl;
+    }
     double ekin = 0;
 
     for (int i = 0; i < m_natoms; ++i) {
@@ -1397,6 +1460,7 @@ void SimpleMD::Rattle(double* grad)
     ekin = 0.0;
     while (iter < m_rattle_maxiter) {
         iter++;
+        int active = 0;
         for (auto bond : m_bond_constrained) {
             int i = bond.first.first, j = bond.first.second;
             if (moved[i] == 1 && moved[j] == 1) {
@@ -1414,7 +1478,8 @@ void SimpleMD::Rattle(double* grad)
                 double mu = -1 * r / ((m_rmass[i] + m_rmass[j]) * distance);
                 while (std::abs(mu) > max_mu)
                     mu /= 2;
-                if (std::abs(mu) > m_rattle_tolerance && std::abs(mu) < max_mu) {
+                if (std::abs(mu) > m_rattle_tolerance) {
+                    active = 1;
                     m_virial_correction += mu * distance;
                     m_velocities[3 * i + 0] += dx * mu * m_rmass[i];
                     m_velocities[3 * i + 1] += dy * mu * m_rmass[i];
@@ -1426,7 +1491,18 @@ void SimpleMD::Rattle(double* grad)
                 }
             }
         }
+        if (active == 0)
+            break;
     }
+
+    if (iter >= m_rattle_maxiter) {
+        std::cout << "numeric difficulties - 2nd in rattle velocity verlet" << iter << std::endl;
+        std::ofstream restart_file("unstable_curcuma_" + std::to_string(m_currentStep) + ".json");
+        nlohmann::json restart;
+        restart[MethodName()[0]] = WriteRestartInformation();
+        restart_file << restart << std::endl;
+    }
+
     if (move)
         RemoveRotations(m_velocities);
 
@@ -1436,7 +1512,7 @@ void SimpleMD::Rattle(double* grad)
     }
     ekin *= 0.5;
     double T = 2.0 * ekin / (kb_Eh * m_dof);
-    m_unstable = T > 10000 * m_T;
+    m_unstable = T > 10000 * m_T || std::isnan(T);
     m_T = T;
     ThermostatFunction();
     EKin();

Diff for: external/xtb

@@ -151,7 +151,9 @@ static json CurcumaMDJson{
     { "rattle", false },
     { "rattle_tolerance", 1e-1 },
     { "rattle_maxiter", 100 },
-    { "thermostat", "csvr" },
+    { "rattle_dynamic_tol", false },
+    { "rattle_dynamic_tol_iter", 100 },
+    { "thermostat", "csvr" }, // can be csvr (default), berendson, none, anderson or nosehover
     { "respa", 1 },
     { "threads", 1 },
     { "dipole", false },
@@ -254,6 +256,8 @@ class SimpleMD : public CurcumaMethod {
     void Rattle_Verlet_Second(double* coord, double* grad);
     void Rattle_Constrain_Second(double* coord, double* grad);
 
+    void AdjustRattleTolerance();
+
     void RemoveRotation(std::vector<double>& velo);
     void RemoveRotations(std::vector<double>& velo);
 
@@ -294,7 +298,7 @@ class SimpleMD : public CurcumaMethod {
     double m_single_step = 1;
     double m_dT = 0.5, m_currentStep = 0, m_maxtime = 1000;
     int m_spin = 0, m_charge = 0, m_print = 100;
-    double m_T0 = 298.13, m_aver_Temp = 0, m_rmsd = 1.5;
+    double m_T0 = 298.13, m_aver_Temp = 0, m_aver_rattle_Temp = 0, m_rmsd = 1.5;
     double m_x0 = 0, m_y0 = 0, m_z0 = 0;
     double m_Ekin_exchange = 0.0;
     std::vector<double> m_current_geometry, m_mass, m_velocities, m_gradient, m_rmass, m_virial, m_gradient_bias;
@@ -318,6 +322,7 @@ class SimpleMD : public CurcumaMethod {
     int m_unix_started = 0, m_prev_index = 0, m_max_rescue = 10, m_current_rescue = 0, m_currentTime = 0, m_max_top_diff = 15, m_step = 0;
     int m_writerestart = -1;
     int m_respa = 1;
+    int m_rattle_dynamic_tol_iter = 100;
     double m_pos_conv = 0, m_scale_velo = 1.0, m_coupling = 10;
     double m_impuls = 0, m_impuls_scaling = 0.75, m_dt2 = 0;
     double m_rattle_tolerance = 1;
@@ -340,6 +345,7 @@ class SimpleMD : public CurcumaMethod {
     int m_bias_structure_count = 0;
     int m_rmsd_fragment_count = 0;
     int m_wall_type = 0;
+    int m_rattle_counter = 0;
     std::vector<double> m_collected_dipole;
     Matrix m_topo_initial;
     std::vector<Molecule*> m_unique_structures;
@@ -352,7 +358,7 @@ class SimpleMD : public CurcumaMethod {
     bool m_rmsd_mtd = false;
     bool m_wtmtd = false;
     bool m_rmsd_fix_structure = false;
-
+    bool m_rattle_dynamic_tol = false;
     int m_mtd_dT = -1;
     int m_seed = -1;
     int m_time_step = 0;
@@ -367,6 +373,8 @@ class SimpleMD : public CurcumaMethod {
     int m_chain_length = 3; // Länge der Thermostatkette
 
     double m_anderson = 0.01;
+
+    std::vector<std::pair<double, double>> m_rattle_tol_temp;
 };
 
 class MDThread : public CxxThread {

Diff for: src/capabilities/simplemd.cpp

@@ -902,6 +902,8 @@ int main(int argc, char **argv) {
             while (!file.AtEnd()) {
                 Molecule mol = file.Next();
                 std::pair<double, double> gyr = mol.GyrationRadius(hmass);
+                if (std::isnan(gyr.first) || std::isnan(gyr.second))
+                    continue;
                 sum += gyr.first;
                 sum_mass += gyr.second;
                 sqrt_sum += sqrt(gyr.first);