Fix bug resetting mdatoms masses to lambda=0 state

When DD is used, the mdatoms masses were reset to their lambda=0 state
after every DD step. This is fixed by making the DomDecHelper aware of
the FEP element.

Note that the use of MDAtoms should likely be revisited to allow for a
more elegant solution, but this would likely require changes which are
less local than the proposed solution here. Further development is
therefore deferred to #3700.

diff --git a/docs/release-notes/2020/2020.4.rst b/docs/release-notes/2020/2020.4.rst
index 6479820ee8eb4bb9ca93b1f5cc6590ef2c4d92e3..d810a6740c72793ae481332ec7869060643cfe90 100644 (file)
--- a/docs/release-notes/2020/2020.4.rst
+++ b/docs/release-notes/2020/2020.4.rst
@@ -42,6 +42,12 @@ an assertion failure during PME tuning.
 
 :issue:`3677`
 
+Bug fix for FEP calculations with modular simulator and domain decomposition
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+When using the modular simulator, domain decomposition and free energy
+calculations with perturbed masses, the simulation would always be
+performed using the masses at lambda=0 instead of the actual lambda value.
+
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/modularsimulator/domdechelper.cpp b/src/gromacs/modularsimulator/domdechelper.cpp
index e731c73bd733dc1c7b761c7bb553f71137fbe344..d256c59967970f6f7539d0f363f0d8717f994eb5 100644 (file)
--- a/src/gromacs/modularsimulator/domdechelper.cpp
+++ b/src/gromacs/modularsimulator/domdechelper.cpp
@@ -50,6 +50,7 @@
 #include "gromacs/mdtypes/state.h"
 #include "gromacs/pbcutil/pbc.h"
 
+#include "freeenergyperturbationelement.h"
 #include "statepropagatordata.h"
 #include "topologyholder.h"
 
@@ -58,6 +59,7 @@ namespace gmx
 DomDecHelper::DomDecHelper(bool                               isVerbose,
                            int                                verbosePrintInterval,
                            StatePropagatorData*               statePropagatorData,
+                           FreeEnergyPerturbationElement*     freeEnergyPerturbationElement,
                            TopologyHolder*                    topologyHolder,
                            CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback,
                            int                                nstglobalcomm,
@@ -78,6 +80,7 @@ DomDecHelper::DomDecHelper(bool                               isVerbose,
     verbosePrintInterval_(verbosePrintInterval),
     nstglobalcomm_(nstglobalcomm),
     statePropagatorData_(statePropagatorData),
+    freeEnergyPerturbationElement_(freeEnergyPerturbationElement),
     topologyHolder_(topologyHolder),
     checkBondedInteractionsCallback_(std::move(checkBondedInteractionsCallback)),
     fplog_(fplog),
@@ -98,10 +101,6 @@ DomDecHelper::DomDecHelper(bool                               isVerbose,
 
 void DomDecHelper::setup()
 {
-    std::unique_ptr<t_state> localState   = statePropagatorData_->localState();
-    t_state*                 globalState  = statePropagatorData_->globalState();
-    PaddedHostVector<RVec>*  forcePointer = statePropagatorData_->forcePointer();
-
     // constant choices for this call to dd_partition_system
     const bool     verbose       = false;
     const bool     isMasterState = true;
@@ -109,14 +108,8 @@ void DomDecHelper::setup()
     gmx_wallcycle* wcycle        = nullptr;
 
     // Distribute the charge groups over the nodes from the master node
-    dd_partition_system(fplog_, mdlog_, inputrec_->init_step, cr_, isMasterState, nstglobalcomm,
-                        globalState, topologyHolder_->globalTopology(), inputrec_, imdSession_,
-                        pull_work_, localState.get(), forcePointer, mdAtoms_,
-                        topologyHolder_->localTopology_.get(), fr_, vsite_, constr_, nrnb_, wcycle,
-                        verbose);
-    topologyHolder_->updateLocalTopology();
-    (*checkBondedInteractionsCallback_)();
-    statePropagatorData_->setLocalState(std::move(localState));
+    partitionSystem(verbose, isMasterState, nstglobalcomm, wcycle,
+                    statePropagatorData_->localState(), statePropagatorData_->globalState());
 }
 
 void DomDecHelper::run(Step step, Time gmx_unused time)
@@ -125,9 +118,8 @@ void DomDecHelper::run(Step step, Time gmx_unused time)
     {
         return;
     }
-    std::unique_ptr<t_state> localState   = statePropagatorData_->localState();
-    t_state*                 globalState  = statePropagatorData_->globalState();
-    PaddedHostVector<RVec>*  forcePointer = statePropagatorData_->forcePointer();
+    std::unique_ptr<t_state> localState  = statePropagatorData_->localState();
+    t_state*                 globalState = statePropagatorData_->globalState();
 
     // constant choices for this call to dd_partition_system
     const bool verbose = isVerbose_ && (step % verbosePrintInterval_ == 0 || step == inputrec_->init_step);
@@ -150,13 +142,31 @@ void DomDecHelper::run(Step step, Time gmx_unused time)
     }
 
     // Distribute the charge groups over the nodes from the master node
-    dd_partition_system(fplog_, mdlog_, step, cr_, isMasterState, nstglobalcomm_, globalState,
-                        topologyHolder_->globalTopology(), inputrec_, imdSession_, pull_work_,
-                        localState.get(), forcePointer, mdAtoms_, topologyHolder_->localTopology_.get(),
-                        fr_, vsite_, constr_, nrnb_, wcycle_, verbose);
+    partitionSystem(verbose, isMasterState, nstglobalcomm_, wcycle_, std::move(localState), globalState);
+}
+
+void DomDecHelper::partitionSystem(bool                     verbose,
+                                   bool                     isMasterState,
+                                   int                      nstglobalcomm,
+                                   gmx_wallcycle*           wcycle,
+                                   std::unique_ptr<t_state> localState,
+                                   t_state*                 globalState)
+{
+    PaddedHostVector<RVec>* forcePointer = statePropagatorData_->forcePointer();
+
+    // Distribute the charge groups over the nodes from the master node
+    dd_partition_system(fplog_, mdlog_, inputrec_->init_step, cr_, isMasterState, nstglobalcomm,
+                        globalState, topologyHolder_->globalTopology(), inputrec_, imdSession_,
+                        pull_work_, localState.get(), forcePointer, mdAtoms_,
+                        topologyHolder_->localTopology_.get(), fr_, vsite_, constr_, nrnb_, wcycle,
+                        verbose);
     topologyHolder_->updateLocalTopology();
     (*checkBondedInteractionsCallback_)();
     statePropagatorData_->setLocalState(std::move(localState));
+    if (freeEnergyPerturbationElement_)
+    {
+        freeEnergyPerturbationElement_->updateMDAtoms();
+    }
 }
 
 SignallerCallbackPtr DomDecHelper::registerNSCallback()

diff --git a/src/gromacs/modularsimulator/domdechelper.h b/src/gromacs/modularsimulator/domdechelper.h
index 181698f3cc6de26a2f56499bb5bd3c2a96212717..89c3d5ccd0166d12b413b14be6b7e10258f15727 100644 (file)
--- a/src/gromacs/modularsimulator/domdechelper.h
+++ b/src/gromacs/modularsimulator/domdechelper.h
@@ -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, 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.
@@ -56,6 +56,7 @@ struct t_nrnb;
 namespace gmx
 {
 class Constraints;
+class FreeEnergyPerturbationElement;
 class ImdSession;
 class MDAtoms;
 class MDLogger;
@@ -90,6 +91,7 @@ public:
     DomDecHelper(bool                               isVerbose,
                  int                                verbosePrintInterval,
                  StatePropagatorData*               statePropagatorData,
+                 FreeEnergyPerturbationElement*     freeEnergyPerturbationElement,
                  TopologyHolder*                    topologyHolder,
                  CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback,
                  int                                nstglobalcomm,
@@ -135,11 +137,21 @@ private:
 
     //! Pointer to the micro state
     StatePropagatorData* statePropagatorData_;
+    //! Pointer to the free energy element
+    FreeEnergyPerturbationElement* freeEnergyPerturbationElement_;
     //! Pointer to the topology
     TopologyHolder* topologyHolder_;
     //! Pointer to the ComputeGlobalsHelper object - to ask for # of bonded interaction checking
     CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback_;
 
+    //! Helper function unifying the DD partitioning calls in setup() and run()
+    void partitionSystem(bool                     verbose,
+                         bool                     isMasterState,
+                         int                      nstglobalcomm,
+                         gmx_wallcycle*           wcycle,
+                         std::unique_ptr<t_state> localState,
+                         t_state*                 globalState);
+
     // Access to ISimulator data
     //! Handles logging.
     FILE* fplog_;

diff --git a/src/gromacs/modularsimulator/freeenergyperturbationelement.cpp b/src/gromacs/modularsimulator/freeenergyperturbationelement.cpp
index e8f1d56bba05e0301d50f24551ad1e59070afe4a..f44920fe12aa31e0799e72c7ebea3a3b54de636d 100644 (file)
--- a/src/gromacs/modularsimulator/freeenergyperturbationelement.cpp
+++ b/src/gromacs/modularsimulator/freeenergyperturbationelement.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, 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.
@@ -46,6 +46,7 @@
 #include "gromacs/mdlib/md_support.h"
 #include "gromacs/mdlib/mdatoms.h"
 #include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/mdtypes/state.h"
 
 namespace gmx
@@ -65,7 +66,7 @@ FreeEnergyPerturbationElement::FreeEnergyPerturbationElement(FILE*             f
     lambda_.fill(0);
     lambda0_.fill(0);
     initialize_lambdas(fplog_, *inputrec_, true, &currentFEPState_, lambda_, lambda0_.data());
-    update_mdatoms(mdAtoms_->mdatoms(), lambda_[efptMASS]);
+    updateMDAtoms();
 }
 
 void FreeEnergyPerturbationElement::scheduleTask(Step step,
@@ -83,7 +84,7 @@ void FreeEnergyPerturbationElement::updateLambdas(Step step)
 {
     // at beginning of step (if lambdas change...)
     setCurrentLambdasLocal(step, inputrec_->fepvals, lambda0_.data(), lambda_, currentFEPState_);
-    update_mdatoms(mdAtoms_->mdatoms(), lambda_[efptMASS]);
+    updateMDAtoms();
 }
 
 ArrayRef<real> FreeEnergyPerturbationElement::lambdaView()
@@ -101,6 +102,11 @@ int FreeEnergyPerturbationElement::currentFEPState()
     return currentFEPState_;
 }
 
+void FreeEnergyPerturbationElement::updateMDAtoms()
+{
+    update_mdatoms(mdAtoms_->mdatoms(), lambda_[efptMASS]);
+}
+
 void FreeEnergyPerturbationElement::writeCheckpoint(t_state* localState, t_state gmx_unused* globalState)
 {
     localState->fep_state = currentFEPState_;

diff --git a/src/gromacs/modularsimulator/freeenergyperturbationelement.h b/src/gromacs/modularsimulator/freeenergyperturbationelement.h
index f1566ffd764e30c1574f4dcd42a1932c10207768..d1d95a5e3a1adcc83293c77df94a87264d9077ca 100644 (file)
--- a/src/gromacs/modularsimulator/freeenergyperturbationelement.h
+++ b/src/gromacs/modularsimulator/freeenergyperturbationelement.h
@@ -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, 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.
@@ -77,6 +77,8 @@ public:
     ArrayRef<const real> constLambdaView();
     //! Get the current FEP state
     int currentFEPState();
+    //! Update MDAtoms (public because it's called by DomDec - see #3700)
+    void updateMDAtoms();
 
     //! Update lambda and mdatoms
     void scheduleTask(Step step, Time time, const RegisterRunFunctionPtr& registerRunFunction) override;

diff --git a/src/gromacs/modularsimulator/modularsimulator.cpp b/src/gromacs/modularsimulator/modularsimulator.cpp
index 8036261e0e1149c7c1c22dc1b2e551ae34220051..2d785ab50fa68c0645b8fc56fa0c76b47ae76ce5 100644 (file)
--- a/src/gromacs/modularsimulator/modularsimulator.cpp
+++ b/src/gromacs/modularsimulator/modularsimulator.cpp
@@ -450,8 +450,9 @@ void ModularSimulator::constructElementsAndSignallers()
                    "interactions.");
         domDecHelper_ = std::make_unique<DomDecHelper>(
                 mdrunOptions.verbose, mdrunOptions.verboseStepPrintInterval, statePropagatorDataPtr,
-                topologyHolder_.get(), std::move(checkBondedInteractionsCallback), nstglobalcomm_, fplog,
-                cr, mdlog, constr, inputrec, mdAtoms, nrnb, wcycle, fr, vsite, imdSession, pull_work);
+                freeEnergyPerturbationElementPtr, topologyHolder_.get(),
+                std::move(checkBondedInteractionsCallback), nstglobalcomm_, fplog, cr, mdlog,
+                constr, inputrec, mdAtoms, nrnb, wcycle, fr, vsite, imdSession, pull_work);
         neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(domDecHelper_.get()));
     }