Fix check-pointing for density guided simulations

Previously, checkpointing for density-guided simulations did not take
into account that the checkpointing occurs in the middle of the md-loop.

Now, the state is stored before the force-calculations, so that
simulations are restarted correctly.

refs #2282
fixes #3215

Change-Id: Ide4ebf94f2801acbe476a5cd348af66d0a6f28ac

diff --git a/src/gromacs/applied_forces/densityfitting.cpp b/src/gromacs/applied_forces/densityfitting.cpp
index 2f34d1281897680dfb270fc633a13b7e56655d43..28ead0fcfa1b5398876e0b3766a5f9f373683d18 100644 (file)
--- a/src/gromacs/applied_forces/densityfitting.cpp
+++ b/src/gromacs/applied_forces/densityfitting.cpp
@@ -379,12 +379,15 @@ public:
      * \param[in] checkpointWriting enables writing to the Key-Value-Tree
      *                              that is used for storing the checkpoint
      *                              information
+     *
+     * \note The provided state to checkpoint has to change if checkpointing
+     *       is moved before the force provider call in the MD-loop.
      */
     void writeCheckpointData(MdModulesWriteCheckpointData checkpointWriting)
     {
         if (densityFittingOptions_.active())
         {
-            const DensityFittingForceProviderState& state = forceProvider_->state();
+            const DensityFittingForceProviderState& state = forceProvider_->stateToCheckpoint();
             checkpointWriting.builder_.addValue<std::int64_t>(
                     DensityFittingModuleInfo::name_ + "-stepsSinceLastCalculation",
                     state.stepsSinceLastCalculation_);

diff --git a/src/gromacs/applied_forces/densityfittingforceprovider.cpp b/src/gromacs/applied_forces/densityfittingforceprovider.cpp
index 24e40526fb07dbb731f800622aa32150b7da694e..106831562de37c10b3d9e98faf16e22a6ffdee99 100644 (file)
--- a/src/gromacs/applied_forces/densityfittingforceprovider.cpp
+++ b/src/gromacs/applied_forces/densityfittingforceprovider.cpp
@@ -102,12 +102,13 @@ public:
     ~Impl();
     void calculateForces(const ForceProviderInput& forceProviderInput,
                          ForceProviderOutput*      forceProviderOutput);
-
-    DensityFittingForceProviderState state();
+    const DensityFittingForceProviderState& stateToCheckpoint();
 
 private:
+    DensityFittingForceProviderState state();
     const DensityFittingParameters&  parameters_;
     DensityFittingForceProviderState state_;
+    DensityFittingForceProviderState stateToCheckpoint_;
     LocalAtomSet                     localAtomSet_;
 
     GaussianSpreadKernelParameters::Shape spreadKernel_;
@@ -176,6 +177,8 @@ DensityFittingForceProvider::Impl::Impl(const DensityFittingParameters&
 void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput& forceProviderInput,
                                                         ForceProviderOutput* forceProviderOutput)
 {
+    // TODO change if checkpointing moves to the start of the md loop
+    stateToCheckpoint_ = state();
     // do nothing but count number of steps when not in density fitting step
     if (state_.stepsSinceLastCalculation_ % parameters_.calculationIntervalInSteps_ != 0)
     {
@@ -298,6 +301,10 @@ DensityFittingForceProviderState DensityFittingForceProvider::Impl::state()
     return state_;
 }
 
+const DensityFittingForceProviderState& DensityFittingForceProvider::Impl::stateToCheckpoint()
+{
+    return stateToCheckpoint_;
+}
 /********************************************************************
  * DensityFittingForceProvider
  */
@@ -321,9 +328,9 @@ void DensityFittingForceProvider::calculateForces(const ForceProviderInput& forc
     impl_->calculateForces(forceProviderInput, forceProviderOutput);
 }
 
-DensityFittingForceProviderState DensityFittingForceProvider::state()
+const DensityFittingForceProviderState& DensityFittingForceProvider::stateToCheckpoint()
 {
-    return impl_->state();
+    return impl_->stateToCheckpoint();
 }
 
 } // namespace gmx

diff --git a/src/gromacs/applied_forces/densityfittingforceprovider.h b/src/gromacs/applied_forces/densityfittingforceprovider.h
index 3c979865c5bea522a2f36f0cbce0ab6abaca00cf..16a121ba3432d21ee32dc30e9393bac563b80678 100644 (file)
--- a/src/gromacs/applied_forces/densityfittingforceprovider.h
+++ b/src/gromacs/applied_forces/densityfittingforceprovider.h
@@ -93,8 +93,11 @@ public:
     void calculateForces(const ForceProviderInput& forceProviderInput,
                          ForceProviderOutput*      forceProviderOutput) override;
 
-    //! Return the state of the forceprovider.
-    DensityFittingForceProviderState state();
+    /*! \brief Return the state of the forceprovider to be checkpointed
+     * TODO update this routine if checkpointing is moved to the beginning of
+     *      the md loop
+     */
+    const DensityFittingForceProviderState& stateToCheckpoint();
 
 private:
     class Impl;

diff --git a/src/programs/mdrun/tests/densityfittingmodule.cpp b/src/programs/mdrun/tests/densityfittingmodule.cpp
index fa098d82c90409ae19b2babc4bf68017d3a2598c..366a5342e9ecd375a3a3908e06aa4851a9dc6ae5 100644 (file)
--- a/src/programs/mdrun/tests/densityfittingmodule.cpp
+++ b/src/programs/mdrun/tests/densityfittingmodule.cpp
@@ -127,6 +127,10 @@ public:
     const std::string mdpEnergyAndDensityfittingIntervalMismatch_ = formatString(
             "nstcalcenergy = 7\n"
             "density-guided-simulation-nst = 3\n");
+    //! Set mdp values so that we skip every other step
+    const std::string mdpSkipDensityfittingEveryOtherStep_ = formatString(
+            "nstenergy = 2\n"
+            "density-guided-simulation-nst = 2\n");
     //! The command line to call mdrun
     CommandLine commandLineForMdrun_;
 };
@@ -171,5 +175,33 @@ TEST_F(DensityFittingTest, GromppErrorWhenEnergyEvaluationFrequencyMismatch)
                               ".*is not a multiple of density-guided-simulation-nst.*");
 }
 
+/* Fit a subset of three of twelve argon atoms into a reference density
+ * whose origin is offset from the simulation box origin. Stop the simulation,
+ * then restart.
+ *
+ * All density fitting mdp parameters are set to defaults
+ */
+TEST_F(DensityFittingTest, CheckpointWorks)
+{
+    runner_.useStringAsMdpFile(mdpMdDensfitYesUnsetValues + mdpSkipDensityfittingEveryOtherStep_);
+    runner_.cptFileName_ = fileManager_.getTemporaryFilePath(".cpt");
+    commandLineForMdrun_.addOption("-cpo", runner_.cptFileName_);
+
+    ASSERT_EQ(0, runner_.callGrompp());
+    ASSERT_EQ(0, runner_.callMdrun(commandLineForMdrun_));
+
+    // checkMdrun(expectedEnergyTermMagnitude);
+
+    CommandLine commandLineForRestart;
+    commandLineForRestart.addOption("-cpi", runner_.cptFileName_);
+    commandLineForRestart.addOption("-noappend");
+    runner_.nsteps_ = 4;
+    ASSERT_EQ(0, runner_.callMdrun(commandLineForRestart));
+
+    const real expectedEnergyTermMagnitude = -3378.825928;
+    checkMdrun(expectedEnergyTermMagnitude);
+}
+
+
 } // namespace test
 } // namespace gmx

diff --git a/src/programs/mdrun/tests/refdata/DensityFittingTest_CheckpointWorks.xml b/src/programs/mdrun/tests/refdata/DensityFittingTest_CheckpointWorks.xml
new file mode 100644 (file)
index 0000000..b1b3ab3
--- /dev/null
+++ b/src/programs/mdrun/tests/refdata/DensityFittingTest_CheckpointWorks.xml
@@ -0,0 +1,12 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <Energy Name="Potential">
+    <Real Name="Time 0.000000 Step 0 in frame 0">-3379.9202</Real>
+    <Real Name="Time 0.002000 Step 2 in frame 1">-3390.9541</Real>
+  </Energy>
+  <Energy Name="Density fitting">
+    <Real Name="Time 0.000000 Step 0 in frame 0">-3378.9026</Real>
+    <Real Name="Time 0.002000 Step 2 in frame 1">-3389.9324</Real>
+  </Energy>
+</ReferenceData>