From 2f0260cd0dfdd280bc6ea9d956b5d1946193e566 Mon Sep 17 00:00:00 2001
From: Conrad Huebler <Conrad.Huebler@gmx.net>
Date: Tue, 21 Nov 2023 17:28:11 +0100
Subject: [PATCH] make random seed changeable for md runs

Signed-off-by: Conrad Huebler <Conrad.Huebler@gmx.net>
---
 src/capabilities/simplemd.cpp | 11 ++++++++++-
 src/capabilities/simplemd.h   |  4 +++-
 2 files changed, 13 insertions(+), 2 deletions(-)

diff --git a/src/capabilities/simplemd.cpp b/src/capabilities/simplemd.cpp
index 3cdf402..93069a8 100644
--- a/src/capabilities/simplemd.cpp
+++ b/src/capabilities/simplemd.cpp
@@ -71,6 +71,7 @@ void SimpleMD::LoadControlJson()
     m_dumb = Json2KeyWord<int>(m_defaults, "dump");
     m_print = Json2KeyWord<int>(m_defaults, "print");
     m_max_top_diff = Json2KeyWord<int>(m_defaults, "MaxTopoDiff");
+    m_seed = Json2KeyWord<int>(m_defaults, "seed");
 
     m_rmsd = Json2KeyWord<double>(m_defaults, "rmsd");
     m_impuls = Json2KeyWord<double>(m_defaults, "impuls");
@@ -247,7 +248,14 @@ void SimpleMD::InitVelocities(double scaling)
 {
     std::random_device rd{};
     std::mt19937 gen{ rd() };
-    std::normal_distribution<> d{ 0, 1 };
+    if (m_seed == -1) {
+        const auto start = std::chrono::high_resolution_clock::now();
+        m_seed = std::chrono::duration_cast<std::chrono::seconds>(start.time_since_epoch()).count();
+    } else if (m_seed == 0)
+        m_seed = m_T0 * m_mass.size();
+    std::cout << "Random seed is " << m_seed << std::endl;
+    gen.seed(m_seed);
+    std::normal_distribution<> d{ -1, 1 };
     double Px = 0.0, Py = 0.0, Pz = 0.0;
     for (int i = 0; i < m_natoms; ++i) {
         double v0 = sqrt(kb * m_T0 * amu2au / (m_mass[i])) * scaling / fs2amu;
@@ -897,6 +905,7 @@ void SimpleMD::CSVR()
     double c = exp(-(m_timestep * m_respa) / m_coupling);
     std::random_device rd{};
     std::mt19937 gen{ rd() };
+    gen.seed(m_seed);
     std::normal_distribution<> d{ 0, 1 };
     std::chi_squared_distribution<float> dchi{ dof };
 
diff --git a/src/capabilities/simplemd.h b/src/capabilities/simplemd.h
index 85f7ade..ed0dbea 100644
--- a/src/capabilities/simplemd.h
+++ b/src/capabilities/simplemd.h
@@ -63,7 +63,8 @@ static json CurcumaMDJson{
     { "rattle_tolerance", 1e-5 },
     { "thermostat", "csvr" },
     { "respa", 1 },
-    { "dipole", false }
+    { "dipole", false },
+    { "seed", -1 }
 };
 
 class SimpleMD : public CurcumaMethod {
@@ -164,6 +165,7 @@ class SimpleMD : public CurcumaMethod {
     std::string m_method = "UFF", m_initfile = "none", m_thermostat = "csvr";
     bool m_unstable = false;
     bool m_dipole = false;
+    int m_seed = -1;
 };
 
 class MDThread : public CxxThread {