Add lambda state to dhdl.xvg with AWH fep
[alexxy/gromacs.git] / src / gromacs / mdlib / tests / energyoutput.cpp
index dae26029590b2056f3f0aa21573f8dce4752469b..51aed47c64def29dd2bdbc518983c654fada1b64 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,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.
@@ -68,7 +68,7 @@
 #include "gromacs/mdtypes/state.h"
 #include "gromacs/topology/topology.h"
 #include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/mdmodulenotification.h"
+#include "gromacs/utility/mdmodulesnotifiers.h"
 #include "gromacs/utility/stringutil.h"
 #include "gromacs/utility/textreader.h"
 #include "gromacs/utility/unique_cptr.h"
@@ -98,11 +98,11 @@ void fcloseWrapper(FILE* fp)
 struct EnergyOutputTestParameters
 {
     //! Thermostat (enum)
-    int temperatureCouplingScheme;
+    TemperatureCoupling temperatureCouplingScheme;
     //! Barostat (enum)
-    int pressureCouplingScheme;
+    PressureCoupling pressureCouplingScheme;
     //! Integrator
-    int integrator;
+    IntegrationAlgorithm integrator;
     //! Number of saved energy frames (to test averages output).
     int numFrames;
     //! If output should be initialized as a rerun.
@@ -116,17 +116,19 @@ struct EnergyOutputTestParameters
  * Only several combinations of the parameters are used. Using all possible combinations will
  * require ~10 MB of test data and ~2 sec to run the tests.
  */
-const EnergyOutputTestParameters parametersSets[] = { { etcNO, epcNO, eiMD, 1, false, false },
-                                                      { etcNO, epcNO, eiMD, 1, true, false },
-                                                      { etcNO, epcNO, eiMD, 1, false, true },
-                                                      { etcNO, epcNO, eiMD, 0, false, false },
-                                                      { etcNO, epcNO, eiMD, 10, false, false },
-                                                      { etcVRESCALE, epcNO, eiMD, 1, false, false },
-                                                      { etcNOSEHOOVER, epcNO, eiMD, 1, false, false },
-                                                      { etcNO, epcPARRINELLORAHMAN, eiMD, 1, false, false },
-                                                      { etcNO, epcMTTK, eiMD, 1, false, false },
-                                                      { etcNO, epcNO, eiVV, 1, false, false },
-                                                      { etcNO, epcMTTK, eiVV, 1, false, false } };
+const EnergyOutputTestParameters parametersSets[] = {
+    { TemperatureCoupling::No, PressureCoupling::No, IntegrationAlgorithm::MD, 1, false, false },
+    { TemperatureCoupling::No, PressureCoupling::No, IntegrationAlgorithm::MD, 1, true, false },
+    { TemperatureCoupling::No, PressureCoupling::No, IntegrationAlgorithm::MD, 1, false, true },
+    { TemperatureCoupling::No, PressureCoupling::No, IntegrationAlgorithm::MD, 0, false, false },
+    { TemperatureCoupling::No, PressureCoupling::No, IntegrationAlgorithm::MD, 10, false, false },
+    { TemperatureCoupling::VRescale, PressureCoupling::No, IntegrationAlgorithm::MD, 1, false, false },
+    { TemperatureCoupling::NoseHoover, PressureCoupling::No, IntegrationAlgorithm::MD, 1, false, false },
+    { TemperatureCoupling::No, PressureCoupling::ParrinelloRahman, IntegrationAlgorithm::MD, 1, false, false },
+    { TemperatureCoupling::No, PressureCoupling::Mttk, IntegrationAlgorithm::MD, 1, false, false },
+    { TemperatureCoupling::No, PressureCoupling::No, IntegrationAlgorithm::VV, 1, false, false },
+    { TemperatureCoupling::No, PressureCoupling::Mttk, IntegrationAlgorithm::VV, 1, false, false }
+};
 
 /*! \brief Test fixture to test energy output.
  *
@@ -136,6 +138,9 @@ const EnergyOutputTestParameters parametersSets[] = { { etcNO, epcNO, eiMD, 1, f
  */
 class EnergyOutputTest : public ::testing::TestWithParam<EnergyOutputTestParameters>
 {
+    int  numTempCouplingGroups_ = 3;
+    real cosAccel_              = 1.0;
+
 public:
     //! File manager
     TestFileManager fileManager_;
@@ -145,8 +150,6 @@ public:
     t_inputrec inputrec_;
     //! Topology
     gmx_mtop_t mtop_;
-    //! MD atoms
-    t_mdatoms mdatoms_;
     //! Simulation time
     double time_;
     //! Total mass
@@ -159,10 +162,6 @@ public:
     t_state state_;
     //! PBC box
     matrix box_;
-    //! Virial from constraints
-    tensor constraintsVirial_;
-    //! Virial from force computation
-    tensor forceVirial_;
     //! Total virial
     tensor totalVirial_;
     //! Pressure
@@ -175,8 +174,6 @@ public:
     std::vector<char*> groupNameHandles_;
     //! Total dipole momentum
     rvec muTotal_;
-    //! Distance and orientation restraints data
-    t_fcdata fcd_;
     //! Communication record
     t_commrec cr_;
     //! Constraints object (for constraints RMSD output in case of LINCS)
@@ -195,6 +192,7 @@ public:
     TestReferenceChecker checker_;
 
     EnergyOutputTest() :
+        ekindata_(numTempCouplingGroups_, cosAccel_, 1),
         logFilename_(fileManager_.getTemporaryFilePath(".log")),
         edrFilename_(fileManager_.getTemporaryFilePath(".edr")),
         log_(std::fopen(logFilename_.c_str(), "w")),
@@ -204,27 +202,25 @@ public:
         // Input record
         inputrec_.delta_t = 0.001;
 
-        // F_EQM
-        inputrec_.bQMMM = true;
         // F_RF_EXCL will not be tested - group scheme is not supported any more
-        inputrec_.cutoff_scheme = ecutsVERLET;
+        inputrec_.cutoff_scheme = CutoffScheme::Verlet;
         // F_COUL_RECIP
-        inputrec_.coulombtype = eelPME;
+        inputrec_.coulombtype = CoulombInteractionType::Pme;
         // F_LJ_RECIP
-        inputrec_.vdwtype = evdwPME;
+        inputrec_.vdwtype = VanDerWaalsType::Pme;
 
         // F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED, F_DVDL_RESTRAINT, F_DKDL and F_DVDL
-        inputrec_.efep                                  = efepYES;
-        inputrec_.fepvals->separate_dvdl[efptCOUL]      = true;
-        inputrec_.fepvals->separate_dvdl[efptVDW]       = true;
-        inputrec_.fepvals->separate_dvdl[efptBONDED]    = true;
-        inputrec_.fepvals->separate_dvdl[efptRESTRAINT] = true;
-        inputrec_.fepvals->separate_dvdl[efptMASS]      = true;
-        inputrec_.fepvals->separate_dvdl[efptCOUL]      = true;
-        inputrec_.fepvals->separate_dvdl[efptFEP]       = true;
+        inputrec_.efep = FreeEnergyPerturbationType::Yes;
+        inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Coul]      = true;
+        inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Vdw]       = true;
+        inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Bonded]    = true;
+        inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Restraint] = true;
+        inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Mass]      = true;
+        inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Coul]      = true;
+        inputrec_.fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep]       = true;
 
         // F_DISPCORR and F_PDISPCORR
-        inputrec_.eDispCorr = edispcEner;
+        inputrec_.eDispCorr = DispersionCorrectionType::Ener;
         inputrec_.bRot      = true;
 
         // F_ECONSERVED
@@ -233,13 +229,10 @@ public:
         inputrec_.ref_p[ZZ][YY] = 0.0;
 
         // Dipole (mu)
-        inputrec_.ewald_geometry = eewg3DC;
+        inputrec_.ewald_geometry = EwaldGeometry::ThreeDC;
 
-        // GMX_CONSTRAINTVIR environment variable should also be
-        // set to print constraints and force virials separately.
-        gmxSetenv("GMX_CONSTRAINTVIR", "true", 1);
         // To print constrain RMSD, constraints algorithm should be set to LINCS.
-        inputrec_.eConstrAlg = econtLINCS;
+        inputrec_.eConstrAlg = ConstraintAlgorithm::Lincs;
 
         mtop_.bIntermolecularInteractions = false;
 
@@ -314,20 +307,16 @@ public:
             mtop_.groups.groupNames.emplace_back(&handle);
         }
 
-        mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].resize(3);
+        mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput].resize(numTempCouplingGroups_);
         mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput][0] = 0;
         mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput][1] = 1;
         mtop_.groups.groups[SimulationAtomGroupType::EnergyOutput][2] = 2;
 
-        mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].resize(3);
+        mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].resize(numTempCouplingGroups_);
         mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling][0] = 0;
         mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling][1] = 1;
         mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling][2] = 2;
 
-        mtop_.groups.groups[SimulationAtomGroupType::Acceleration].resize(2);
-        mtop_.groups.groups[SimulationAtomGroupType::Acceleration][0] = 0;
-        mtop_.groups.groups[SimulationAtomGroupType::Acceleration][1] = 2;
-
         // Nose-Hoover chains
         inputrec_.bPrintNHChains     = true;
         inputrec_.opts.nhchainlength = 2;
@@ -348,23 +337,17 @@ public:
 
         // Kinetic energy and related data
         ekindata_.tcstat.resize(mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].size());
-        ekindata_.grpstat.resize(mtop_.groups.groups[SimulationAtomGroupType::Acceleration].size());
 
         // This is needed so that the ebin space will be allocated
-        inputrec_.cos_accel = 1.0;
-        // This is to keep the destructor happy (otherwise sfree() segfaults)
-        ekindata_.nthreads = 0;
-        snew(ekindata_.ekin_work_alloc, 1);
-        snew(ekindata_.ekin_work, 1);
-        snew(ekindata_.dekindl_work, 1);
+        inputrec_.cos_accel = cosAccel_;
 
         // Group options for annealing output
-        inputrec_.opts.ngtc = 3;
+        inputrec_.opts.ngtc = numTempCouplingGroups_;
         snew(inputrec_.opts.ref_t, inputrec_.opts.ngtc);
         snew(inputrec_.opts.annealing, inputrec_.opts.ngtc);
-        inputrec_.opts.annealing[0] = eannNO;
-        inputrec_.opts.annealing[1] = eannSINGLE;
-        inputrec_.opts.annealing[2] = eannPERIODIC;
+        inputrec_.opts.annealing[0] = SimulatedAnnealing::No;
+        inputrec_.opts.annealing[1] = SimulatedAnnealing::Single;
+        inputrec_.opts.annealing[2] = SimulatedAnnealing::Periodic;
 
         // This is to keep done_inputrec happy (otherwise sfree() segfaults)
         snew(inputrec_.opts.anneal_time, inputrec_.opts.ngtc);
@@ -378,12 +361,8 @@ public:
         // TODO EnergyOutput should not take Constraints object
         // TODO This object will always return zero as RMSD value.
         //      It is more relevant to have non-zero value for testing.
-        constraints_ = makeConstraints(mtop_, inputrec_, nullptr, false, nullptr, mdatoms_, &cr_,
-                                       nullptr, nullptr, nullptr, false);
-
-        // No position/orientation restraints
-        fcd_.disres.npair = 0;
-        fcd_.orires.nr    = 0;
+        constraints_ = makeConstraints(
+                mtop_, inputrec_, nullptr, false, nullptr, &cr_, false, nullptr, nullptr, nullptr, false, nullptr);
     }
 
     /*! \brief Helper function to generate synthetic data to output
@@ -431,9 +410,9 @@ public:
         // Group pairs
         for (int i = 0; i < enerdata_->grpp.nener; i++)
         {
-            for (int k = 0; k < egNR; k++)
+            for (int k = 0; k < static_cast<int>(NonBondedEnergyTerms::Count); k++)
             {
-                enerdata_->grpp.ener[k][i] = (*testValue += 0.1);
+                enerdata_->grpp.energyGroupPairTerms[k][i] = (*testValue += 0.1);
             }
         }
 
@@ -443,12 +422,9 @@ public:
             tcstat.T      = (*testValue += 0.1);
             tcstat.lambda = (*testValue += 0.1);
         }
-        for (auto& grpstat : ekindata_.grpstat)
-        {
-            grpstat.u[XX] = (*testValue += 0.1);
-            grpstat.u[YY] = (*testValue += 0.1);
-            grpstat.u[ZZ] = (*testValue += 0.1);
-        }
+        // Removing constant acceleration removed a total increment of 0.6
+        // To avoid unnecessary changes in reference data, we keep the increment
+        (*testValue += 0.6);
 
         // This conditional is to check whether the ebin was allocated.
         // Otherwise it will print cosacc data into the first bin.
@@ -478,25 +454,9 @@ public:
         box_[ZZ][YY] = (*testValue += 0.1);
         box_[ZZ][ZZ] = (*testValue += 0.1);
 
-        constraintsVirial_[XX][XX] = (*testValue += 0.1);
-        constraintsVirial_[XX][YY] = (*testValue += 0.1);
-        constraintsVirial_[XX][ZZ] = (*testValue += 0.1);
-        constraintsVirial_[YY][XX] = (*testValue += 0.1);
-        constraintsVirial_[YY][YY] = (*testValue += 0.1);
-        constraintsVirial_[YY][ZZ] = (*testValue += 0.1);
-        constraintsVirial_[ZZ][XX] = (*testValue += 0.1);
-        constraintsVirial_[ZZ][YY] = (*testValue += 0.1);
-        constraintsVirial_[ZZ][ZZ] = (*testValue += 0.1);
-
-        forceVirial_[XX][XX] = (*testValue += 0.1);
-        forceVirial_[XX][YY] = (*testValue += 0.1);
-        forceVirial_[XX][ZZ] = (*testValue += 0.1);
-        forceVirial_[YY][XX] = (*testValue += 0.1);
-        forceVirial_[YY][YY] = (*testValue += 0.1);
-        forceVirial_[YY][ZZ] = (*testValue += 0.1);
-        forceVirial_[ZZ][XX] = (*testValue += 0.1);
-        forceVirial_[ZZ][YY] = (*testValue += 0.1);
-        forceVirial_[ZZ][ZZ] = (*testValue += 0.1);
+        // Removing GMX_CONSTRVIR removed a total increment of 1.8
+        // To avoid unnecessary changes in reference data, we keep the increment
+        (*testValue += 1.8);
 
         totalVirial_[XX][XX] = (*testValue += 0.1);
         totalVirial_[XX][YY] = (*testValue += 0.1);
@@ -632,23 +592,45 @@ TEST_P(EnergyOutputTest, CheckOutput)
         inputrec_.ref_p[YY][XX] = 1.0;
     }
 
-    MdModulesNotifier             mdModulesNotifier;
-    std::unique_ptr<EnergyOutput> energyOutput = std::make_unique<EnergyOutput>(
-            energyFile_, &mtop_, &inputrec_, nullptr, nullptr, parameters.isRerun,
-            StartingBehavior::NewSimulation, mdModulesNotifier);
+    MDModulesNotifiers            mdModulesNotifiers;
+    std::unique_ptr<EnergyOutput> energyOutput =
+            std::make_unique<EnergyOutput>(energyFile_,
+                                           mtop_,
+                                           inputrec_,
+                                           nullptr,
+                                           nullptr,
+                                           parameters.isRerun,
+                                           StartingBehavior::NewSimulation,
+                                           false,
+                                           mdModulesNotifiers);
 
     // Add synthetic data for a single step
     double testValue = 10.0;
     for (int frame = 0; frame < parameters.numFrames; frame++)
     {
         setStepData(&testValue);
-        energyOutput->addDataAtEnergyStep(false, true, time_, tmass_, enerdata_.get(), &state_, nullptr,
-                                          nullptr, box_, constraintsVirial_, forceVirial_, totalVirial_,
-                                          pressure_, &ekindata_, muTotal_, constraints_.get());
+        energyOutput->addDataAtEnergyStep(false,
+                                          true,
+                                          time_,
+                                          tmass_,
+                                          enerdata_.get(),
+                                          nullptr,
+                                          box_,
+                                          PTCouplingArrays({ state_.boxv,
+                                                             state_.nosehoover_xi,
+                                                             state_.nosehoover_vxi,
+                                                             state_.nhpres_xi,
+                                                             state_.nhpres_vxi }),
+                                          state_.fep_state,
+                                          totalVirial_,
+                                          pressure_,
+                                          &ekindata_,
+                                          muTotal_,
+                                          constraints_.get());
 
         energyOutput->printAnnealingTemperatures(log_, &mtop_.groups, &inputrec_.opts);
-        energyOutput->printStepToEnergyFile(energyFile_, true, false, false, log_, 100 * frame,
-                                            time_, nullptr, nullptr);
+        energyOutput->printStepToEnergyFile(
+                energyFile_, true, false, false, log_, 100 * frame, time_, nullptr, nullptr);
         time_ += 1.0;
     }
 
@@ -674,7 +656,7 @@ TEST_P(EnergyOutputTest, CheckOutput)
     checker_.checkString(TextReader::readFileToString(logFilename_), "log");
 }
 
-INSTANTIATE_TEST_CASE_P(WithParameters, EnergyOutputTest, ::testing::ValuesIn(parametersSets));
+INSTANTIATE_TEST_SUITE_P(WithParameters, EnergyOutputTest, ::testing::ValuesIn(parametersSets));
 
 } // namespace
 } // namespace test