bug fix: use correct eksum

Author: zhi-wang

include/md/integrator.h

@@ -153,7 +153,7 @@ class BasicThermostat : virtual public IntegratorStaticData
    virtual ~BasicThermostat();
    virtual void printDetail(FILE*);
    virtual void control1(time_prec timeStep);
-   virtual void control2(time_prec timeStep, bool save);
+   virtual void control2(time_prec timeStep, bool calcEkin);
 
    static BasicThermostat* create(ThermostatEnum);
 };
@@ -187,15 +187,15 @@ class NhcDevice : public BasicThermostat
    int nnose, nhc_nc;
    double g0;
    double vnh[maxnose], qnh[maxnose];
+   double* m_kin_ptr;
    double (*f_kin)();
    void (*scale_vel)(double);
    std::string name;
 
-   void controlImpl(double timeStep);
+   void controlImpl(double timeStep, bool calcEkin);
 
 public:
-   NhcDevice(int nhclen, int nc, double dfree, //
-      double (*kin)(), void (*scale)(double), std::string str);
+   NhcDevice(int nhclen, int nc, double dfree, double* kin_ptr, double (*kin)(), void (*scale)(double), std::string str);
    void printDetail(FILE*) override;
    void control1(time_prec time_prec) override;
    void control2(time_prec time_prec, bool) override;
@@ -217,7 +217,7 @@ class Nhc06Thermostat : public BasicThermostat
 
    void printDetail(FILE*) override;
    void control1(time_prec dt) override;
-   void control2(time_prec dt, bool save) override;
+   void control2(time_prec dt, bool) override;
 
    static double kineticRattleGroup();
    static void scaleVelocityRattleGroup(double scale);
@@ -299,6 +299,7 @@ class IsoBaroDevice : public BasicBarostat
 protected:
    double* m_vir;
    double* m_eksum;
+   double (*f_kin)();
 
    double m_fric;
    double m_rnd;
@@ -322,6 +323,7 @@ class AnisoBaroDevice : public BasicBarostat
    double* m_vir;
    double* m_eksum;
    double (*m_ekin)[3];
+   void (*f_kin)();
 
    double m_fric;
    double m_rnd[3][3];
@@ -404,8 +406,7 @@ class BasicIntegrator : virtual public IntegratorStaticData
 
 public:
    void printDetail(FILE*);
-   BasicIntegrator(int nRespaLogV, PropagatorEnum pe, ThermostatEnum te,
-      BarostatEnum be);
+   BasicIntegrator(int nRespaLogV, PropagatorEnum pe, ThermostatEnum te, BarostatEnum be);
    BasicIntegrator();
    virtual ~BasicIntegrator();
    virtual void dynamic(int istep, time_prec dt);

src/md/barostat.cpp

@@ -173,7 +173,11 @@ void IsoBaroDevice::control_1_2(time_prec dt, int idx)
    double dim = 3.0;
    double al = 1.0 + dim / dofP;
    double tr_vir = m_vir[0] + m_vir[4] + m_vir[8];
-   double eksu1 = *m_eksum;
+   double eksu1;
+   if (idx == 1)
+      eksu1 = *m_eksum;
+   else
+      eksu1 = f_kin();
 
    if (idx == 2 and printPressure) {
       auto o = stdout;
@@ -203,16 +207,19 @@ IsoBaroDevice::IsoBaroDevice(double fric)
    : BasicBarostat()
    , m_vir(nullptr)
    , m_eksum(nullptr)
+   , f_kin(nullptr)
    , m_fric(fric)
    , m_rnd()
    , m_langevin(fric != 0.0)
 {
    if (atomic) {
       m_vir = vir;
       m_eksum = &eksum;
+      f_kin = NhcDevice::kineticAtomic;
    } else {
       m_vir = hc_vir;
       m_eksum = &hc_eksum;
+      f_kin = Nhc06Thermostat::kineticRattleGroup;
    }
 
    if (m_langevin)
@@ -273,6 +280,12 @@ void IsoBaroDevice::control3(time_prec dt)
 }
 
 namespace tinker {
+static void kineticForAnisoBaroDeviceAtomic()
+{
+   double temp;
+   kinetic(temp);
+}
+
 void AnisoBaroDevice::control_1_2(time_prec dt, int idx)
 {
    time_prec dt2 = dt * 0.5;
@@ -281,6 +294,8 @@ void AnisoBaroDevice::control_1_2(time_prec dt, int idx)
    if (m_langevin)
       b = std::sqrt(2 * m_fric * units::gasconst * bath::kelvin / qbar);
 
+   if (idx == 2)
+      f_kin();
    double eksu1 = *m_eksum;
    double gbar[3][3] = {0};
    for (int i = 0; i < 3; ++i)
@@ -354,6 +369,7 @@ AnisoBaroDevice::AnisoBaroDevice(double fric)
    , m_vir(nullptr)
    , m_eksum(nullptr)
    , m_ekin(nullptr)
+   , f_kin(nullptr)
    , m_fric(fric)
    , m_rnd()
    , m_langevin(fric != 0.0)
@@ -362,10 +378,12 @@ AnisoBaroDevice::AnisoBaroDevice(double fric)
       m_vir = vir;
       m_eksum = &eksum;
       m_ekin = ekin;
+      f_kin = kineticForAnisoBaroDeviceAtomic;
    } else {
       m_vir = hc_vir;
       m_eksum = &hc_eksum;
       m_ekin = hc_ekin;
+      f_kin = hcKinetic;
    }
 
    if (m_langevin) {
@@ -431,8 +449,7 @@ void AnisoBaroDevice::control3(time_prec dt)
 
    double scal[3][3];
    trimatExp(scal, vbar_matrix, dt);
-   double h0[3][3] = {
-      {lvec1.x, lvec1.y, lvec1.z}, {lvec2.x, lvec2.y, lvec2.z}, {lvec3.x, lvec3.y, lvec3.z}};
+   double h0[3][3] = {{lvec1.x, lvec1.y, lvec1.z}, {lvec2.x, lvec2.y, lvec2.z}, {lvec3.x, lvec3.y, lvec3.z}};
    matmul3(h0, scal);
    lvec1.x = h0[0][0], lvec1.y = h0[0][1], lvec1.z = h0[0][2];
    lvec2.x = h0[1][0], lvec2.y = h0[1][1], lvec2.z = h0[1][2];

src/md/thermostat.cpp

@@ -8,6 +8,16 @@
 
 #include <cmath>
 
+namespace tinker {
+namespace {
+class DummyThermostat : public BasicThermostat
+{
+public:
+   void control2(time_prec, bool) override {}
+};
+}
+}
+
 namespace tinker {
 void BasicThermostat::printBasic(FILE* o)
 {
@@ -23,9 +33,9 @@ void BasicThermostat::printDetail(FILE* o) {}
 
 void BasicThermostat::control1(time_prec) {}
 
-void BasicThermostat::control2(time_prec, bool save)
+void BasicThermostat::control2(time_prec, bool calcEkin)
 {
-   if (save) {
+   if (calcEkin) {
       T_prec temp;
       kinetic(temp);
    }
@@ -40,10 +50,14 @@ BasicThermostat* BasicThermostat::create(ThermostatEnum te)
       break;
    case ThermostatEnum::NHC:
       t = new NhcDevice(5, 5, static_cast<double>(mdstuf::nfree), //
-         NhcDevice::kineticAtomic,                                //
+         &eksum, NhcDevice::kineticAtomic,                        //
          NhcDevice::scaleVelocityAtomic,                          //
          std::string("NHC Atomic Temperature"));
       break;
+   case ThermostatEnum::m_LP2022:
+   case ThermostatEnum::m_NHC2006:
+      t = new DummyThermostat;
+      break;
    default:
       t = new BasicThermostat;
       break;
@@ -69,7 +83,7 @@ void BussiThermostat::control2(time_prec dt, bool)
 }
 
 namespace tinker {
-void NhcDevice::controlImpl(time_prec dt)
+void NhcDevice::controlImpl(time_prec dt, bool calcEkin)
 {
    const int nc = nhc_nc;
    constexpr int ns = nhc_nsy;
@@ -82,7 +96,11 @@ void NhcDevice::controlImpl(time_prec dt)
    double kbt = units::gasconst * bath::kelvin;
 
    double dtc = dt / nc;
-   double eksum0 = f_kin();
+   double eksum0;
+   if (calcEkin or m_kin_ptr == nullptr)
+      eksum0 = f_kin();
+   else
+      eksum0 = *m_kin_ptr;
    double velsc0 = 1.0;
    for (int k = 0; k < nc; ++k) {
       for (int j = 0; j < ns; ++j) {
@@ -129,12 +147,13 @@ void NhcDevice::controlImpl(time_prec dt)
    scale_vel(velsc0);
 }
 
-NhcDevice::NhcDevice(
-   int nhclen, int nc, double dfree, double (*kin)(), void (*scale)(double), std::string str)
+NhcDevice::NhcDevice(int nhclen, int nc, double dfree, double* kin_ptr, double (*kin)(), void (*scale)(double),
+   std::string str)
    : BasicThermostat()
    , nnose(nhclen)
    , nhc_nc(nc)
    , g0(dfree)
+   , m_kin_ptr(kin_ptr)
    , f_kin(kin)
    , scale_vel(scale)
    , name(str)
@@ -152,6 +171,8 @@ NhcDevice::NhcDevice(
       vnh[i] = 0;
    }
    qnh[0] *= g0;
+
+   f_kin();
 }
 
 void NhcDevice::printDetail(FILE* o)
@@ -172,12 +193,13 @@ void NhcDevice::printDetail(FILE* o)
 
 void NhcDevice::control1(time_prec dt)
 {
-   this->controlImpl(dt);
+   bool calcEkin = false;
+   this->controlImpl(dt, calcEkin);
 }
 
-void NhcDevice::control2(time_prec dt, bool)
+void NhcDevice::control2(time_prec dt, bool calcEkin)
 {
-   this->controlImpl(dt);
+   this->controlImpl(dt, calcEkin);
 }
 
 double NhcDevice::kineticAtomic()
@@ -209,16 +231,16 @@ Nhc06Thermostat::Nhc06Thermostat()
    // tpart
    if (atomic) {
       dofT = mdstuf::nfree;
-      m_tpart = new NhcDevice(nhclen, nc, dofT, NhcDevice::kineticAtomic,
-         NhcDevice::scaleVelocityAtomic, "NHC Atomic Temperature");
+      m_tpart = new NhcDevice(nhclen, nc, dofT, &eksum, NhcDevice::kineticAtomic, NhcDevice::scaleVelocityAtomic,
+         "NHC Atomic Temperature");
    } else {
       dofT = 3.0 * (rattle_dmol.nmol - 1);
-      m_tpart = new NhcDevice(nhclen, nc, dofT, Nhc06Thermostat::kineticRattleGroup,
+      m_tpart = new NhcDevice(nhclen, nc, dofT, &hc_eksum, Nhc06Thermostat::kineticRattleGroup,
          Nhc06Thermostat::scaleVelocityRattleGroup, "NHC Group Temperature");
    }
 
    // tbaro
-   m_tbaro = new NhcDevice(nhclen, nc, dofVbar(), Nhc06Thermostat::kineticVbar,
+   m_tbaro = new NhcDevice(nhclen, nc, dofVbar(), nullptr, Nhc06Thermostat::kineticVbar,
       Nhc06Thermostat::scaleVelocityVbar, "NHC Barostat Temperature");
 }
 
@@ -241,11 +263,17 @@ void Nhc06Thermostat::control1(time_prec dt)
    m_tpart->control1(dt);
 }
 
-void Nhc06Thermostat::control2(time_prec dt, bool save)
+void Nhc06Thermostat::control2(time_prec dt, bool)
 {
-   m_tpart->control2(dt, save);
-   if (applyBaro)
-      m_tbaro->control2(dt, false);
+   bool calcEkin;
+   if (applyBaro) {
+      calcEkin = false; // in nhc06 barostat, ek will be calculated in the barostat
+      m_tpart->control2(dt, calcEkin);
+      m_tbaro->control2(dt, calcEkin);
+   } else {
+      calcEkin = true;
+      m_tpart->control2(dt, calcEkin);
+   }
 }
 
 double Nhc06Thermostat::kineticRattleGroup()