Reimplement constant acceleration groups

[alexxy/gromacs.git] / src / gromacs / mdlib / tests / leapfrogtestdata.cpp

diff --git a/src/gromacs/mdlib/tests/leapfrogtestdata.cpp b/src/gromacs/mdlib/tests/leapfrogtestdata.cpp
index 38ae157c642c15c0b593118fe8c8566e9be34d73..cf80cda412530db1c9f349f97979d534dc18b612 100644 (file)
--- a/src/gromacs/mdlib/tests/leapfrogtestdata.cpp
+++ b/src/gromacs/mdlib/tests/leapfrogtestdata.cpp
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2019, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -67,26 +67,33 @@ namespace gmx
 namespace test
 {
 
-LeapFrogTestData::LeapFrogTestData(int numAtoms, real timestep, const rvec v0, const rvec f0) :
+LeapFrogTestData::LeapFrogTestData(int        numAtoms,
+                                   real       timestep,
+                                   const rvec v0,
+                                   const rvec f0,
+                                   int        numTCoupleGroups,
+                                   int        nstpcouple) :
     numAtoms_(numAtoms),
     timestep_(timestep),
-    x0_(numAtoms_),
-    x_(numAtoms_),
-    xPrime_(numAtoms_),
-    v0_(numAtoms_),
-    v_(numAtoms_),
-    f_(numAtoms_),
-    inverseMasses_(numAtoms_),
-    inverseMassesPerDim_(numAtoms_)
+    x0_(numAtoms),
+    x_(numAtoms),
+    xPrime_(numAtoms),
+    v0_(numAtoms),
+    v_(numAtoms),
+    f_(numAtoms),
+    inverseMasses_(numAtoms),
+    inverseMassesPerDim_(numAtoms),
+    kineticEnergyData_(numTCoupleGroups == 0 ? 1 : numTCoupleGroups, 0.0, 1),
+    numTCoupleGroups_(numTCoupleGroups)
 {
     mdAtoms_.nr = numAtoms_;
 
     for (int i = 0; i < numAtoms_; i++)
     {
         // Typical PBC box size is tens of nanometers
-        x_[i][XX] = (i%21)*1.0;
-        x_[i][YY] = 6.5 + (i%13)*(-1.0);
-        x_[i][ZZ] = (i%32)*(0.0);
+        x_[i][XX] = (i % 21) * 1.0;
+        x_[i][YY] = 6.5 + (i % 13) * (-1.0);
+        x_[i][ZZ] = (i % 32) * (0.0);
 
         for (int d = 0; d < DIM; d++)
         {
@@ -100,7 +107,7 @@ LeapFrogTestData::LeapFrogTestData(int numAtoms, real timestep, const rvec v0, c
             v0_[i][d] = v_[i][d];
         }
         // Atom masses are ~1-100 g/mol
-        inverseMasses_[i] = 1.0/(1.0 + i%100);
+        inverseMasses_[i] = 1.0 / (1.0 + i % 100);
         for (int d = 0; d < DIM; d++)
         {
             inverseMassesPerDim_[i][d] = inverseMasses_[i];
@@ -109,11 +116,41 @@ LeapFrogTestData::LeapFrogTestData(int numAtoms, real timestep, const rvec v0, c
     mdAtoms_.invmass       = inverseMasses_.data();
     mdAtoms_.invMassPerDim = as_rvec_array(inverseMassesPerDim_.data());
 
-    // Data needed for current CPU-based implementation
-    inputRecord_.eI      = eiMD;
+    // Temperature coupling
+    snew(mdAtoms_.cTC, numAtoms_);
+
+    // To do temperature coupling at each step
+    inputRecord_.nsttcouple = 1;
+
+    if (numTCoupleGroups_ == 0)
+    {
+        inputRecord_.etc = TemperatureCoupling::No;
+        for (int i = 0; i < numAtoms_; i++)
+        {
+            mdAtoms_.cTC[i] = 0;
+        }
+        t_grp_tcstat temperatureCouplingGroupData;
+        temperatureCouplingGroupData.lambda = 1.0;
+        kineticEnergyData_.tcstat[0]        = temperatureCouplingGroupData;
+    }
+    else
+    {
+        inputRecord_.etc = TemperatureCoupling::Yes;
+        for (int i = 0; i < numAtoms_; i++)
+        {
+            mdAtoms_.cTC[i] = i % numTCoupleGroups_;
+        }
+        for (int i = 0; i < numTCoupleGroups_; i++)
+        {
+            real         tCoupleLambda = 1.0 - (i + 1.0) / 10.0;
+            t_grp_tcstat temperatureCouplingGroupData;
+            temperatureCouplingGroupData.lambda = tCoupleLambda;
+            kineticEnergyData_.tcstat[i]        = temperatureCouplingGroupData;
+        }
+    }
+
+    inputRecord_.eI      = IntegrationAlgorithm::MD;
     inputRecord_.delta_t = timestep_;
-    inputRecord_.etc     = etcNO;
-    inputRecord_.epc     = epcNO;
 
     state_.flags = 0;
 
@@ -129,40 +166,53 @@ LeapFrogTestData::LeapFrogTestData(int numAtoms, real timestep, const rvec v0, c
     state_.box[ZZ][YY] = 0.0;
     state_.box[ZZ][ZZ] = 10.0;
 
-    kineticEnergyData_.bNEMD            = false;
-    kineticEnergyData_.cosacc.cos_accel = 0.0;
-    t_grp_tcstat temperatureCouplingGroupData;
-    temperatureCouplingGroupData.lambda = 1;
-    kineticEnergyData_.tcstat.emplace_back(temperatureCouplingGroupData);
-
-    kineticEnergyData_.nthreads = 1;
-    snew(kineticEnergyData_.ekin_work_alloc, kineticEnergyData_.nthreads);
-    snew(kineticEnergyData_.ekin_work, kineticEnergyData_.nthreads);
-    snew(kineticEnergyData_.dekindl_work, kineticEnergyData_.nthreads);
-
     mdAtoms_.homenr                   = numAtoms_;
     mdAtoms_.haveVsites               = false;
     mdAtoms_.havePartiallyFrozenAtoms = false;
-    snew(mdAtoms_.cTC, numAtoms_);
-    for (int i = 0; i < numAtoms_; i++)
-    {
-        mdAtoms_.cTC[i] = 0;
-    }
+    mdAtoms_.cFREEZE                  = nullptr;
+    mdAtoms_.ptype                    = nullptr;
 
-    prVScalingMatrix_[XX][XX] = 1.0;
-    prVScalingMatrix_[XX][YY] = 0.0;
-    prVScalingMatrix_[XX][ZZ] = 0.0;
+    update_ = std::make_unique<Update>(inputRecord_, nullptr);
+    update_->updateAfterPartition(numAtoms,
+                                  gmx::ArrayRef<const unsigned short>(),
+                                  gmx::arrayRefFromArray(mdAtoms_.cTC, mdAtoms_.nr),
+                                  gmx::ArrayRef<const unsigned short>());
 
-    prVScalingMatrix_[YY][XX] = 0.0;
-    prVScalingMatrix_[YY][YY] = 1.0;
-    prVScalingMatrix_[YY][ZZ] = 0.0;
+    doPressureCouple_ = (nstpcouple != 0);
 
-    prVScalingMatrix_[ZZ][XX] = 0.0;
-    prVScalingMatrix_[ZZ][YY] = 0.0;
-    prVScalingMatrix_[ZZ][ZZ] = 1.0;
+    if (doPressureCouple_)
+    {
+        inputRecord_.epc        = PressureCoupling::ParrinelloRahman;
+        inputRecord_.nstpcouple = nstpcouple;
+        dtPressureCouple_       = inputRecord_.nstpcouple * inputRecord_.delta_t;
+
+        velocityScalingMatrix_[XX][XX] = 1.2;
+        velocityScalingMatrix_[XX][YY] = 0.0;
+        velocityScalingMatrix_[XX][ZZ] = 0.0;
+
+        velocityScalingMatrix_[YY][XX] = 0.0;
+        velocityScalingMatrix_[YY][YY] = 0.8;
+        velocityScalingMatrix_[YY][ZZ] = 0.0;
 
-    update_ = std::make_unique<Update>(&inputRecord_, nullptr);
-    update_->setNumAtoms(numAtoms);
+        velocityScalingMatrix_[ZZ][XX] = 0.0;
+        velocityScalingMatrix_[ZZ][YY] = 0.0;
+        velocityScalingMatrix_[ZZ][ZZ] = 0.9;
+    }
+    else
+    {
+        inputRecord_.epc               = PressureCoupling::No;
+        velocityScalingMatrix_[XX][XX] = 1.0;
+        velocityScalingMatrix_[XX][YY] = 0.0;
+        velocityScalingMatrix_[XX][ZZ] = 0.0;
+
+        velocityScalingMatrix_[YY][XX] = 0.0;
+        velocityScalingMatrix_[YY][YY] = 1.0;
+        velocityScalingMatrix_[YY][ZZ] = 0.0;
+
+        velocityScalingMatrix_[ZZ][XX] = 0.0;
+        velocityScalingMatrix_[ZZ][YY] = 0.0;
+        velocityScalingMatrix_[ZZ][ZZ] = 1.0;
+    }
 }
 
 LeapFrogTestData::~LeapFrogTestData()
@@ -170,5 +220,5 @@ LeapFrogTestData::~LeapFrogTestData()
     sfree(mdAtoms_.cTC);
 }
 
-}  // namespace test
-}  // namespace gmx
+} // namespace test
+} // namespace gmx