Adding every N-steps option to densityfitting
authorChristian Blau <cblau@gwdg.de>
Wed, 4 Sep 2019 13:26:19 +0000 (15:26 +0200)
committerPaul Bauer <paul.bauer.q@gmail.com>
Mon, 9 Sep 2019 05:58:50 +0000 (07:58 +0200)
Adds an option to the density fitting module to apply forces only every
N steps.

refs #2282

Change-Id: I8d264734e1d9a82ff9717d270552e38a28b34729

src/gromacs/applied_forces/densityfitting.cpp
src/gromacs/applied_forces/densityfittingforceprovider.cpp
src/gromacs/applied_forces/densityfittingforceprovider.h
src/gromacs/applied_forces/densityfittingoptions.cpp
src/gromacs/applied_forces/densityfittingoptions.h
src/gromacs/applied_forces/densityfittingparameters.h
src/gromacs/applied_forces/tests/densityfitting.cpp
src/gromacs/applied_forces/tests/densityfittingoptions.cpp
src/gromacs/fileio/checkpoint.cpp
src/gromacs/fileio/checkpoint.h

index 0d122c2e1f2c3b17a263e946d1ad745819fd96f5..1867f339ee8d79eb35b25754e2619e96d1beb5aa 100644 (file)
 #include <memory>
 
 #include "gromacs/domdec/localatomsetmanager.h"
+#include "gromacs/fileio/checkpoint.h"
 #include "gromacs/fileio/mrcdensitymap.h"
 #include "gromacs/math/multidimarray.h"
+#include "gromacs/mdlib/broadcaststructs.h"
+#include "gromacs/mdtypes/commrec.h"
 #include "gromacs/mdtypes/imdmodule.h"
 #include "gromacs/utility/classhelpers.h"
 #include "gromacs/utility/exceptions.h"
@@ -212,7 +215,13 @@ class DensityFitting final : public IMDModule
          *   - constructing local atom sets in the simulation parameter setup
          *     by taking a LocalAtomSetManager * as parameter
          *   - the type of periodic boundary conditions that are used
-         *     by taking a PeriodicBoundaryConditionTypeEnum as parameter
+         *     by taking a PeriodicBoundaryConditionType as parameter
+         *   - the writing of checkpoint data
+         *     by taking a MdModulesWriteCheckpointData as parameter
+         *   - the reading of checkpoint data
+         *     by taking a MdModulesCheckpointReadingDataOnMaster as parameter
+         *   - the broadcasting of checkpoint data
+         *     by taking MdModulesCheckpointReadingBroadcast as parameter
          */
         explicit DensityFitting(MdModulesNotifier *notifier)
         {
@@ -256,6 +265,24 @@ class DensityFitting final : public IMDModule
                         this->setEnergyOutputRequest(energyOutputRequest);
                     };
             notifier->notifier_.subscribe(requestEnergyOutput);
+
+            // writing checkpoint data
+            const auto checkpointDataWriting = [this](MdModulesWriteCheckpointData checkpointData) {
+                    this->writeCheckpointData(checkpointData);
+                };
+            notifier->notifier_.subscribe(checkpointDataWriting);
+
+            // reading checkpoint data
+            const auto checkpointDataReading = [this](MdModulesCheckpointReadingDataOnMaster checkpointData) {
+                    this->readCheckpointDataOnMaster(checkpointData);
+                };
+            notifier->notifier_.subscribe(checkpointDataReading);
+
+            // broadcasting checkpoint data
+            const auto checkpointDataBroadcast = [this](MdModulesCheckpointReadingBroadcast checkpointData) {
+                    this->broadcastCheckpointData(checkpointData);
+                };
+            notifier->notifier_.subscribe(checkpointDataBroadcast);
         }
 
         //! From IMDModule
@@ -273,7 +300,8 @@ class DensityFitting final : public IMDModule
                             densityFittingSimulationParameters_.referenceDensity(),
                             densityFittingSimulationParameters_.transformationToDensityLattice(),
                             densityFittingSimulationParameters_.localAtomSet(),
-                            densityFittingSimulationParameters_.periodicBoundaryConditionType());
+                            densityFittingSimulationParameters_.periodicBoundaryConditionType(),
+                            densityFittingState_);
                 forceProviders->addForceProvider(forceProvider_.get());
             }
         }
@@ -305,6 +333,55 @@ class DensityFitting final : public IMDModule
             energyOutputRequest->energyOutputToDensityFitting_ = densityFittingOptions_.active();
         }
 
+        /*! \brief Write internal density fitting data to checkpoint file.
+         * \param[in] checkpointWriting enables writing to the Key-Value-Tree
+         *                              that is used for storing the checkpoint
+         *                              information
+         */
+        void writeCheckpointData(MdModulesWriteCheckpointData checkpointWriting)
+        {
+            if (densityFittingOptions_.active())
+            {
+                const DensityFittingForceProviderState &state =
+                    forceProvider_->state();
+                checkpointWriting.builder_.addValue<std::int64_t>(
+                        DensityFittingModuleInfo::name_ + "-stepsSinceLastCalculation",
+                        state.stepsSinceLastCalculation_);
+            }
+        }
+
+        /*! \brief Read the internal parameters from the checkpoint file on master
+         * \param[in] checkpointReading holding the checkpoint information
+         */
+        void readCheckpointDataOnMaster(MdModulesCheckpointReadingDataOnMaster checkpointReading)
+        {
+            if (densityFittingOptions_.active())
+            {
+                if (checkpointReading.checkpointedData_.keyExists(
+                            DensityFittingModuleInfo::name_ + "-stepsSinceLastCalculation"))
+                {
+                    densityFittingState_.stepsSinceLastCalculation_ =
+                        checkpointReading.checkpointedData_[
+                            DensityFittingModuleInfo::name_ + "-stepsSinceLastCalculation"].cast<std::int64_t>();
+                }
+            }
+        }
+
+        /*! \brief Broadcast the internal parameters.
+         * \param[in] checkpointBroadcast containing the communication record to
+         *                                broadcast the checkpoint information
+         */
+        void broadcastCheckpointData(MdModulesCheckpointReadingBroadcast checkpointBroadcast)
+        {
+            if (densityFittingOptions_.active())
+            {
+                if (PAR(&(checkpointBroadcast.cr_)))
+                {
+                    block_bc(&(checkpointBroadcast.cr_), densityFittingState_.stepsSinceLastCalculation_);
+                }
+            }
+        }
+
     private:
         //! The output provider
         DensityFittingOutputProvider                 densityFittingOutputProvider_;
@@ -315,7 +392,9 @@ class DensityFitting final : public IMDModule
         /*! \brief Parameters for density fitting that become available at
          * simulation setup time.
          */
-        DensityFittingSimulationParameterSetup       densityFittingSimulationParameters_;
+        DensityFittingSimulationParameterSetup        densityFittingSimulationParameters_;
+        //! The internal parameters of density fitting force provider
+        DensityFittingForceProviderState              densityFittingState_;
 
         GMX_DISALLOW_COPY_AND_ASSIGN(DensityFitting);
 };
index 57ccd175cc86dbb0523b07346f7be15c248f05df..60c6ce1001e51a535b731ce04ceeab023fa24638 100644 (file)
@@ -97,11 +97,16 @@ class DensityFittingForceProvider::Impl
              basic_mdspan<const float, dynamicExtents3D> referenceDensity,
              const TranslateAndScale &transformationToDensityLattice,
              const LocalAtomSet &localAtomSet,
-             int pbcType);
+             int pbcType,
+             const DensityFittingForceProviderState &state);
         ~Impl();
         void calculateForces(const ForceProviderInput &forceProviderInput, ForceProviderOutput *forceProviderOutput);
+
+        DensityFittingForceProviderState state();
+
     private:
         const DensityFittingParameters       &parameters_;
+        DensityFittingForceProviderState      state_;
         LocalAtomSet                          localAtomSet_;
 
         GaussianSpreadKernelParameters::Shape spreadKernel_;
@@ -115,6 +120,8 @@ class DensityFittingForceProvider::Impl
         TranslateAndScale                     transformationToDensityLattice_;
         RVec                                  referenceDensityCenter_;
         int                                   pbcType_;
+        real                                  forceConstantScale_;
+
 };
 
 DensityFittingForceProvider::Impl::~Impl() = default;
@@ -123,8 +130,10 @@ DensityFittingForceProvider::Impl::Impl(const DensityFittingParameters &paramete
                                         basic_mdspan<const float, dynamicExtents3D> referenceDensity,
                                         const TranslateAndScale &transformationToDensityLattice,
                                         const LocalAtomSet &localAtomSet,
-                                        int pbcType) :
+                                        int pbcType,
+                                        const DensityFittingForceProviderState &state) :
     parameters_(parameters),
+    state_(state),
     localAtomSet_(localAtomSet),
     spreadKernel_(makeSpreadKernel(parameters_.gaussianTransformSpreadingWidth_,
                                    parameters_.gaussianTransformSpreadingRangeInMultiplesOfWidth_,
@@ -135,7 +144,8 @@ DensityFittingForceProvider::Impl::Impl(const DensityFittingParameters &paramete
     transformedCoordinates_(localAtomSet_.numAtomsLocal()),
     amplitudeLookup_(parameters_.amplitudeLookupMethod_),
     transformationToDensityLattice_(transformationToDensityLattice),
-    pbcType_(pbcType)
+    pbcType_(pbcType),
+    forceConstantScale_(parameters_.calculationIntervalInSteps_)
 {
     referenceDensityCenter_  = {
         real(referenceDensity.extent(XX))/2,
@@ -149,6 +159,15 @@ DensityFittingForceProvider::Impl::Impl(const DensityFittingParameters &paramete
 void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput &forceProviderInput,
                                                         ForceProviderOutput      *forceProviderOutput)
 {
+    // do nothing but count number of steps when not in density fitting step
+    if (state_.stepsSinceLastCalculation_ % parameters_.calculationIntervalInSteps_ != 0)
+    {
+        ++(state_.stepsSinceLastCalculation_);
+        return;
+    }
+
+    state_.stepsSinceLastCalculation_ = 0;
+
     // do nothing if there are no density fitting atoms on this node
     if (localAtomSet_.numAtomsLocal() == 0)
     {
@@ -217,7 +236,7 @@ void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput
     for (const auto localAtomIndex : localAtomSet_.localIndex())
     {
         forceProviderOutput->forceWithVirial_.force_[localAtomIndex] +=
-            parameters_.forceConstant_ * *densityForceIterator;
+            forceConstantScale_ * parameters_.forceConstant_ * *densityForceIterator;
         ++densityForceIterator;
     }
 
@@ -227,6 +246,12 @@ void DensityFittingForceProvider::Impl::calculateForces(const ForceProviderInput
     forceProviderOutput->enerd_.term[F_DENSITYFITTING] += energy;
 }
 
+DensityFittingForceProviderState
+DensityFittingForceProvider::Impl::state()
+{
+    return state_;
+}
+
 /********************************************************************
  * DensityFittingForceProvider
  */
@@ -237,8 +262,9 @@ DensityFittingForceProvider::DensityFittingForceProvider(const DensityFittingPar
                                                          basic_mdspan<const float, dynamicExtents3D> referenceDensity,
                                                          const TranslateAndScale &transformationToDensityLattice,
                                                          const LocalAtomSet &localAtomSet,
-                                                         int pbcType)
-    : impl_(new Impl(parameters, referenceDensity, transformationToDensityLattice, localAtomSet, pbcType))
+                                                         int pbcType,
+                                                         const DensityFittingForceProviderState &state)
+    : impl_(new Impl(parameters, referenceDensity, transformationToDensityLattice, localAtomSet, pbcType, state))
 {}
 
 void DensityFittingForceProvider::calculateForces(const ForceProviderInput  &forceProviderInput,
@@ -247,4 +273,9 @@ void DensityFittingForceProvider::calculateForces(const ForceProviderInput  &for
     impl_->calculateForces(forceProviderInput, forceProviderOutput);
 }
 
+DensityFittingForceProviderState DensityFittingForceProvider::state()
+{
+    return impl_->state();
+}
+
 } // namespace gmx
index a5dd1b285e92373f6923433e9a615710f06b91ea..70b8dbf71e53bde677bd2f2adeeef7b60ba417a0 100644 (file)
@@ -55,6 +55,17 @@ namespace gmx
 
 struct DensityFittingParameters;
 
+/*! \internal
+ * \brief Parameters defining the internal density fitting force provider state.
+ */
+struct DensityFittingForceProviderState
+{
+    /*! \brief The steps since the last force calculation.
+     *  Used if density fitting is to be calculated every N steps.
+     */
+    std::int64_t stepsSinceLastCalculation_ = 0;
+};
+
 /*! \internal \brief
  * Implements IForceProvider for density-fitting forces.
  */
@@ -66,14 +77,19 @@ class DensityFittingForceProvider final : public IForceProvider
                                     basic_mdspan<const float, dynamicExtents3D> referenceDensity,
                                     const TranslateAndScale &transformationToDensityLattice,
                                     const LocalAtomSet &localAtomSet,
-                                    int pbcType);
+                                    int pbcType,
+                                    const DensityFittingForceProviderState &state);
         ~DensityFittingForceProvider();
-        /* \brief Calculate forces that maximise goodness-of-fit with a reference density map
+        /*!\brief Calculate forces that maximise goodness-of-fit with a reference density map.
          * \param[in] forceProviderInput input for force provider
          * \param[out] forceProviderOutput output for force provider
          */
         void calculateForces(const ForceProviderInput &forceProviderInput,
                              ForceProviderOutput      *forceProviderOutput) override;
+
+        //! Return the state of the forceprovider.
+        DensityFittingForceProviderState state();
+
     private:
         class Impl;
         PrivateImplPointer<Impl> impl_;
index ccdb0043de503bc6a878d130f0fe527efc34efa8..7317e4ca4fc0f116d7330eb8aa1d1d4d588434f4 100644 (file)
@@ -120,6 +120,7 @@ void DensityFittingOptions::initMdpTransform(IKeyValueTreeTransformRules * rules
     densityfittingMdpTransformFromString<real>(rules, &fromStdString<real>, c_gaussianTransformSpreadingWidthTag_);
     densityfittingMdpTransformFromString<real>(rules, &fromStdString<real>, c_gaussianTransformSpreadingRangeInMultiplesOfWidthTag_);
     densityfittingMdpTransformFromString<std::string>(rules, stringIdentityTransform, c_referenceDensityFileNameTag_);
+    densityfittingMdpTransformFromString<std::int64_t>(rules, &fromStdString<std::int64_t>, c_everyNStepsTag_);
 }
 
 void DensityFittingOptions::buildMdpOutput(KeyValueTreeObjectBuilder *builder) const
@@ -140,6 +141,7 @@ void DensityFittingOptions::buildMdpOutput(KeyValueTreeObjectBuilder *builder) c
         addDensityFittingMdpOutputValue(builder, parameters_.gaussianTransformSpreadingWidth_, c_gaussianTransformSpreadingWidthTag_);
         addDensityFittingMdpOutputValue(builder, parameters_.gaussianTransformSpreadingRangeInMultiplesOfWidth_, c_gaussianTransformSpreadingRangeInMultiplesOfWidthTag_);
         addDensityFittingMdpOutputValue(builder, referenceDensityFileName_, c_referenceDensityFileNameTag_);
+        addDensityFittingMdpOutputValue(builder, parameters_.calculationIntervalInSteps_, c_everyNStepsTag_);
     }
 }
 
@@ -162,6 +164,7 @@ void DensityFittingOptions::initMdpOptions(IOptionsContainerWithSections *option
     section.addOption(RealOption(c_gaussianTransformSpreadingWidthTag_.c_str()).store(&parameters_.gaussianTransformSpreadingWidth_));
     section.addOption(RealOption(c_gaussianTransformSpreadingRangeInMultiplesOfWidthTag_.c_str()).store(&parameters_.gaussianTransformSpreadingRangeInMultiplesOfWidth_));
     section.addOption(StringOption(c_referenceDensityFileNameTag_.c_str()).store(&referenceDensityFileName_));
+    section.addOption(Int64Option(c_everyNStepsTag_.c_str()).store(&parameters_.calculationIntervalInSteps_));
 }
 
 bool DensityFittingOptions::active() const
@@ -171,6 +174,12 @@ bool DensityFittingOptions::active() const
 
 const DensityFittingParameters &DensityFittingOptions::buildParameters()
 {
+    // the options modules does not know unsigned integers so any input of this
+    // kind is rectified here
+    if (parameters_.calculationIntervalInSteps_ < 1)
+    {
+        parameters_.calculationIntervalInSteps_ = 1;
+    }
     return parameters_;
 }
 
index fa5ed45e7c89bea4c037bacb05c98212b18cb782..047d42223523848c5fa9cfae0e0dbbf0895abcd8 100644 (file)
@@ -123,6 +123,8 @@ class DensityFittingOptions final : public IMdpOptionProvider
         const std::string c_referenceDensityFileNameTag_ = "reference-density-filename";
         std::string       referenceDensityFileName_      = "reference.mrc";
 
+        const std::string c_everyNStepsTag_ = "nst";
+
 
         DensityFittingParameters parameters_;
 };
index f585e26b465c2161b16965986a9a0adccc89331d..d1030d841f631cf8014e4acba3cec210996b3963 100644 (file)
@@ -73,6 +73,8 @@ struct DensityFittingParameters
     real                           gaussianTransformSpreadingWidth_ = 0.2;
     //! The spreading range for spreading atoms onto the grid in multiples of the spreading width
     real                           gaussianTransformSpreadingRangeInMultiplesOfWidth_ = 4.0;
+    //! Apply density fitting forces only every n-steps
+    std::int64_t                   calculationIntervalInSteps_ = 1;
 };
 
 /*!\brief Check if two structs holding density fitting parameters are equal.
index c0a188fd0bb49ba2743edc4f1c61e95c289984f7..dbbcc1505aa350f6df41bb085ca0f9fc25bc1d99 100644 (file)
@@ -115,7 +115,7 @@ class DensityFittingTest : public ::testing::Test
             densityFittingModule_->initForceProviders(&densityFittingForces_);
         }
 
-    private:
+    protected:
 
         KeyValueTreeBuilder        mdpValueBuilder_;
         MdModulesNotifier          notifier_;
@@ -145,7 +145,6 @@ TEST_F(DensityFittingTest, SingleAtom)
     EXPECT_ANY_THROW(intializeForceProviders());
 }
 
-
 } // namespace
 
 } // namespace gmx
index bd5d6d722d962c6880e78fd05a331d22c2c03f37..05e226e42bc3f3820b47182d52cb4554c606bcad 100644 (file)
@@ -197,6 +197,7 @@ TEST_F(DensityFittingOptionsTest, OutputDefaultValuesWhenActive)
         "density-guided-simulation-gaussian-transform-spreading-width = 0.2\n"
         "density-guided-simulation-gaussian-transform-spreading-range-in-multiples-of-width = 4\n"
         "density-guided-simulation-reference-density-filename = reference.mrc\n"
+        "density-guided-simulation-nst = 1\n"
         };
 
     EXPECT_EQ(expected, stream.toString());
index 6250e76ba0f4fb518e21f2f1b5e7a1afc920de5a..7cd4037ec22b3da0089e9540b4424eda18a006ab 100644 (file)
@@ -267,39 +267,6 @@ enum class StatePart
     pullHistory         //!< Pull history statistics (sums since last written output)
 };
 
-namespace gmx
-{
-
-struct MdModulesCheckpointReadingDataOnMaster
-{
-    //! The data of the MdModules that is stored in the checkpoint file
-    const KeyValueTreeObject &checkpointedData;
-    //! The version of the read ceckpoint file
-    int                       checkpointFileVersion_;
-};
-
-/*! \libinternal
- * \brief Provides the MdModules with the communication record to broadcast.
- */
-struct MdModulesCheckpointReadingBroadcast
-{
-    //! The communication record
-    const t_commrec &cr;
-    //! The version of the read file version
-    int              checkpointFileVersion_;
-};
-
-/*! \libinternal \brief Writing the MdModules data to a checkpoint file.
- */
-struct MdModulesWriteCheckpointData
-{
-    //! Builder for the Key-Value-Tree to store the MdModule checkpoint data
-    KeyValueTreeObjectBuilder builder;
-    //! The version of the read file version
-    int                       checkpointFileVersion_;
-};
-} // namespace gmx
-
 //! \brief Return the name of a checkpoint entry based on part and part entry
 static const char *entryName(StatePart part, int ecpt)
 {
index a99aff3e33e304ad84c4123ed559cd2403069f15..1941a6213f359cf49d526b55a8b7aa39267a805f 100644 (file)
@@ -44,6 +44,7 @@
 
 #include "gromacs/math/vectypes.h"
 #include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/keyvaluetreebuilder.h"
 
 class energyhistory_t;
 struct gmx_file_position_t;
@@ -53,10 +54,43 @@ struct t_fileio;
 struct t_inputrec;
 class t_state;
 struct t_trxframe;
+
 namespace gmx
 {
+
 struct MdModulesNotifier;
-}
+class KeyValueTreeObject;
+
+struct MdModulesCheckpointReadingDataOnMaster
+{
+    //! The data of the MdModules that is stored in the checkpoint file
+    const KeyValueTreeObject &checkpointedData_;
+    //! The version of the read ceckpoint file
+    int                       checkpointFileVersion_;
+};
+
+/*! \libinternal
+ * \brief Provides the MdModules with the communication record to broadcast.
+ */
+struct MdModulesCheckpointReadingBroadcast
+{
+    //! The communication record
+    const t_commrec &cr_;
+    //! The version of the read file version
+    int              checkpointFileVersion_;
+};
+
+/*! \libinternal \brief Writing the MdModules data to a checkpoint file.
+ */
+struct MdModulesWriteCheckpointData
+{
+    //! Builder for the Key-Value-Tree to store the MdModule checkpoint data
+    KeyValueTreeObjectBuilder builder_;
+    //! The version of the read file version
+    int                       checkpointFileVersion_;
+};
+
+} // namespace gmx
 
 /* the name of the environment variable to disable fsync failure checks with */
 #define GMX_IGNORE_FSYNC_FAILURE_ENV "GMX_IGNORE_FSYNC_FAILURE"