Add lambda state to dhdl.xvg with AWH fep
[alexxy/gromacs.git] / src / gromacs / mdlib / tests / energyoutput.cpp
index fac276b9f65285b5f61a53fe232234ea81356dc4..51aed47c64def29dd2bdbc518983c654fada1b64 100644 (file)
@@ -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"
@@ -138,6 +138,9 @@ const EnergyOutputTestParameters parametersSets[] = {
  */
 class EnergyOutputTest : public ::testing::TestWithParam<EnergyOutputTestParameters>
 {
+    int  numTempCouplingGroups_ = 3;
+    real cosAccel_              = 1.0;
+
 public:
     //! File manager
     TestFileManager fileManager_;
@@ -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
@@ -193,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")),
@@ -202,8 +202,6 @@ 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 = CutoffScheme::Verlet;
         // F_COUL_RECIP
@@ -233,9 +231,6 @@ public:
         // Dipole (mu)
         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 = ConstraintAlgorithm::Lincs;
 
@@ -312,12 +307,12 @@ 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;
@@ -344,15 +339,10 @@ public:
         ekindata_.tcstat.resize(mtop_.groups.groups[SimulationAtomGroupType::TemperatureCoupling].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] = SimulatedAnnealing::No;
@@ -372,7 +362,7 @@ public:
         // 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, &cr_, nullptr, nullptr, nullptr, false);
+                mtop_, inputrec_, nullptr, false, nullptr, &cr_, false, nullptr, nullptr, nullptr, false, nullptr);
     }
 
     /*! \brief Helper function to generate synthetic data to output
@@ -464,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);
@@ -618,7 +592,7 @@ TEST_P(EnergyOutputTest, CheckOutput)
         inputrec_.ref_p[YY][XX] = 1.0;
     }
 
-    MdModulesNotifier             mdModulesNotifier;
+    MDModulesNotifiers            mdModulesNotifiers;
     std::unique_ptr<EnergyOutput> energyOutput =
             std::make_unique<EnergyOutput>(energyFile_,
                                            mtop_,
@@ -628,7 +602,7 @@ TEST_P(EnergyOutputTest, CheckOutput)
                                            parameters.isRerun,
                                            StartingBehavior::NewSimulation,
                                            false,
-                                           mdModulesNotifier);
+                                           mdModulesNotifiers);
 
     // Add synthetic data for a single step
     double testValue = 10.0;
@@ -641,7 +615,6 @@ TEST_P(EnergyOutputTest, CheckOutput)
                                           tmass_,
                                           enerdata_.get(),
                                           nullptr,
-                                          nullptr,
                                           box_,
                                           PTCouplingArrays({ state_.boxv,
                                                              state_.nosehoover_xi,
@@ -649,8 +622,6 @@ TEST_P(EnergyOutputTest, CheckOutput)
                                                              state_.nhpres_xi,
                                                              state_.nhpres_vxi }),
                                           state_.fep_state,
-                                          constraintsVirial_,
-                                          forceVirial_,
                                           totalVirial_,
                                           pressure_,
                                           &ekindata_,
@@ -685,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