Grompp error for mismatching nst for energy calulation and densityfitting
authorChristian Blau <cblau@gwdg.de>
Fri, 20 Sep 2019 09:50:33 +0000 (11:50 +0200)
committerPaul Bauer <paul.bauer.q@gmail.com>
Mon, 23 Sep 2019 10:01:51 +0000 (12:01 +0200)
When using multiple-time-stepping with density-guided simulations, the
energy evaluation frequency of the density-guided simulations must
match the "nstcalcenergy" option.

Change-Id: Ifcde8cada8a91540d69c9f91079bfc91c4415ba5

docs/reference-manual/special/density-guided-simulation.rst
src/gromacs/applied_forces/densityfitting.cpp
src/gromacs/applied_forces/densityfittingoptions.cpp
src/gromacs/applied_forces/densityfittingoptions.h
src/gromacs/gmxpreprocess/grompp.cpp
src/gromacs/gmxpreprocess/readir.cpp
src/gromacs/gmxpreprocess/readir.h
src/gromacs/gmxpreprocess/tests/readir.cpp
src/gromacs/utility/mdmodulenotification.h
src/programs/mdrun/tests/densityfittingmodule.cpp

index 7362c70597e0fc5681567c122dcc6c5533a415e8..301da01260e2019e8047907ec270425814356cd9 100644 (file)
@@ -137,8 +137,8 @@ The cost of applying forces every integration step is reduced when applying the
 density-guided simulation forces only every :math:`N` steps. The applied force
 is scaled by :math:`N` to approximate the same effective Hamiltonian as when
 applying the forces every step, while maintaining time-reversibility and energy
-conservation. Note that for this setting, the energy output frequency should
-be a multiple of :math:`N`.
+conservation. Note that for this setting, the energy output frequency must be a
+multiple of :math:`N`.
 
 The maximal time-step should not exceed the fastest oscillation period of any
 atom within the map potential divided by :math:`\pi`. This oscillation period
index fb5477fa8282f03afc6ad330fd40acb5a26924c1..d067bcb727caeffaa53bedeea5a5ca615a3e8638 100644 (file)
@@ -263,6 +263,12 @@ class DensityFitting final : public IMDModule
                 };
             notifier->notifier_.subscribe(readInternalParametersFunction);
 
+            // Checking for consistency with all .mdp options
+            const auto checkEnergyCaluclationFrequencyFunction = [this](EnergyCalculationFrequencyErrors * energyCalculationFrequencyErrors) {
+                    densityFittingOptions_.checkEnergyCaluclationFrequency(energyCalculationFrequencyErrors);
+                };
+            notifier->notifier_.subscribe(checkEnergyCaluclationFrequencyFunction);
+
             // constructing local atom sets during simulation setup
             const auto setLocalAtomSetFunction = [this](LocalAtomSetManager *localAtomSetManager) {
                     this->constructLocalAtomSet(localAtomSetManager);
index 2543cd543c1692848da2a9fc7b982b6820e6620c..bf9f89aa127a98f38e5840478ca0a330f115e55f 100644 (file)
@@ -51,6 +51,7 @@
 #include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/keyvaluetreebuilder.h"
 #include "gromacs/utility/keyvaluetreetransform.h"
+#include "gromacs/utility/mdmodulenotification.h"
 #include "gromacs/utility/strconvert.h"
 
 #include "densityfittingamplitudelookup.h"
@@ -248,6 +249,17 @@ void DensityFittingOptions::readInternalParametersFromKvt(const KeyValueTreeObje
                    [](const KeyValueTreeValue &val) { return val.cast<std::int64_t>(); });
 }
 
+void DensityFittingOptions::checkEnergyCaluclationFrequency(EnergyCalculationFrequencyErrors * energyCalculationFrequencyErrors) const
+{
+    if (energyCalculationFrequencyErrors->energyCalculationIntervalInSteps() % parameters_.calculationIntervalInSteps_ != 0)
+    {
+        energyCalculationFrequencyErrors->addError("nstcalcenergy (" +
+                                                   toString(energyCalculationFrequencyErrors->energyCalculationIntervalInSteps())
+                                                   + ") is not a multiple of " + DensityFittingModuleInfo::name_ + "-" + c_everyNStepsTag_
+                                                   + " (" + toString(parameters_.calculationIntervalInSteps_) + ") .");
+    }
+}
+
 const std::string &DensityFittingOptions::referenceDensityFileName() const
 {
     return referenceDensityFileName_;
index d629a7905fac28c463d5e21e84398e04885ee1da..034f3d3a3354bfbe4840d7058a07bceae19c674b 100644 (file)
@@ -51,6 +51,7 @@
 namespace gmx
 {
 
+class EnergyCalculationFrequencyErrors;
 class IndexGroupsAndNames;
 class KeyValueTreeObject;
 class KeyValueTreeBuilder;
@@ -101,6 +102,9 @@ class DensityFittingOptions final : public IMdpOptionProvider
         //! Return the file name of the reference density
         const std::string &referenceDensityFileName() const;
 
+        //! Check if input parameters are consistent with other simulation parameters
+        void checkEnergyCaluclationFrequency(EnergyCalculationFrequencyErrors * energyCalculationFrequencyErrors) const;
+
     private:
         const std::string c_activeTag_ = "active";
 
index a3a2c3f05fb7d78c3c75cfb3239d8dfe694b2a45..ba57fa1f0d798492ee58b589176e59025d20ea10 100644 (file)
@@ -1882,7 +1882,7 @@ int gmx_grompp(int argc, char *argv[])
     {
         fprintf(stderr, "checking input for internal consistency...\n");
     }
-    check_ir(mdparin, ir, opts, wi);
+    check_ir(mdparin, mdModules.notifier(), ir, opts, wi);
 
     if (ir->ld_seed == -1)
     {
index 8ac8242b70edae5ef2768e3ee5ba1c63a92ba06c..339c7185806d32f4c3a15c83affd5cb2d495df68 100644 (file)
@@ -238,8 +238,8 @@ static void process_interaction_modifier(int *eintmod)
     }
 }
 
-void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
-              warninp_t wi)
+void check_ir(const char *mdparin, const gmx::MdModulesNotifier &mdModulesNotifier,
+              t_inputrec *ir, t_gromppopts *opts, warninp_t wi)
 /* Check internal consistency.
  * NOTE: index groups are not set here yet, don't check things
  * like temperature coupling group options here, but in triple_check
@@ -527,6 +527,17 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
             check_nst("nstcalcenergy", ir->nstcalcenergy,
                       "nstenergy", &ir->nstenergy, wi);
         }
+
+        // Inquire all MdModules, if their parameters match with the energy
+        // calculation frequency
+        gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
+        mdModulesNotifier.notifier_.notify(&energyCalculationFrequencyErrors);
+
+        // Emit all errors from the energy calculation frequency checks
+        for (const std::string &energyFrequencyErrorMessage : energyCalculationFrequencyErrors.errorMessages())
+        {
+            warning_error(wi, energyFrequencyErrorMessage);
+        }
     }
 
     if (ir->nsteps == 0 && !ir->bContinuation)
index 3be50585b575ec52eba40d31dfaa9816c2a45572..f5375791fbcef9c1e22ab20353d00e02cb346bb3 100644 (file)
@@ -93,8 +93,8 @@ void init_inputrec_strings();
 /*! \brief Clean up object that holds strings parsed from an .mdp file */
 void done_inputrec_strings();
 
-void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
-              warninp_t wi);
+void check_ir(const char *mdparin, const gmx::MdModulesNotifier &mdModulesNotifier,
+              t_inputrec *ir, t_gromppopts *opts, warninp_t wi);
 /* Validate inputrec data.
  * Fatal errors will be added to nerror.
  */
index a1c83c0522401de9c17fa37ae5f499c98454de68..a8a391ba09ce04c1567a44ee9017aea8a19863ed 100644 (file)
@@ -99,7 +99,7 @@ class GetIrTest : public ::testing::Test
             get_ir(inputMdpFilename.c_str(), outputMdpFilename.c_str(),
                    &mdModules_, &ir_, &opts_, WriteMdpHeader::no, wi_);
 
-            check_ir(inputMdpFilename.c_str(), &ir_, &opts_, wi_);
+            check_ir(inputMdpFilename.c_str(), mdModules_.notifier(), &ir_, &opts_, wi_);
             // Now check
             bool                 failure = warning_errors_exist(wi_);
             TestReferenceData    data;
index 7f669aa32c7642ce0edb8262fefc48fa2c5f536c..654007ecae860fb1d05b5741242c6a0625dbd1b6 100644 (file)
@@ -45,6 +45,7 @@
 #define GMX_MDRUNUTILITY_MDMODULENOTIFICATION_H
 
 #include <functional>
+#include <string>
 #include <vector>
 
 struct t_commrec;
@@ -209,11 +210,50 @@ struct MdModulesEnergyOutputToDensityFittingRequestChecker
     bool energyOutputToDensityFitting_ = false;
 };
 
+/*! \libinternal
+ * \brief Collect errors for the energy calculation frequency.
+ *
+ * Collect errors regarding energy calculation frequencies as strings that then
+ * may be used to issue errors.
+ *
+ * \note The mdp option "nstcalcenergy" is altered after reading the .mdp input
+ *       and only used in certain integrators, thus this class is to be used
+ *       only after all these operations are done.
+ */
+class EnergyCalculationFrequencyErrors
+{
+    public:
+        //! Construct by setting the energy calculation frequency
+        EnergyCalculationFrequencyErrors(int64_t energyCalculationIntervalInSteps) :
+            energyCalculationIntervalInSteps_(energyCalculationIntervalInSteps){}
+        //! Return the number of steps of an energy calculation interval
+        std::int64_t energyCalculationIntervalInSteps() const
+        {
+            return energyCalculationIntervalInSteps_;
+        }
+        //! Collect error messages
+        void addError(const std::string &errorMessage)
+        {
+            errorMessages_.push_back(errorMessage);
+        }
+        //! Return error messages
+        const std::vector<std::string> &errorMessages() const
+        {
+            return errorMessages_;
+        }
+    private:
+        //! The frequency of energy calculations
+        const std::int64_t       energyCalculationIntervalInSteps_;
+        //! The error messages
+        std::vector<std::string> errorMessages_;
+};
+
 struct MdModulesNotifier
 {
 //! Register callback function types for MdModule
     registerMdModuleNotification<
         const t_commrec &,
+        EnergyCalculationFrequencyErrors *,
         IndexGroupsAndNames,
         KeyValueTreeObjectBuilder,
         const KeyValueTreeObject &,
index 39797876eaba44f1fb44c5576c472316c538dc92..adf4e4d2418696d148e1e97479151ef56332b641 100644 (file)
@@ -44,6 +44,9 @@
 
 #include <string>
 
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
 #include "gromacs/topology/ifunc.h"
 #include "gromacs/utility/stringutil.h"
 
@@ -101,6 +104,16 @@ class DensityFittingTest : public MdrunTestFixture
                     "density-guided-simulation-reference-density-filename = %s\n",
                     TestFileManager::getInputFilePath("ellipsoid-density.mrc").c_str() );
 
+        //! Mdp values for md integrator with default density fitting parameters.
+        const std::string mdpMdDensfitYesUnsetValues = formatString(
+                    "integrator                       = md\n"
+                    "nsteps                           = 2\n"
+                    "cutoff-scheme                    = verlet\n"
+                    "density-guided-simulation-active = yes\n"
+                    "density-guided-simulation-group  = FirstThreeOfTwelve\n"
+                    "density-guided-simulation-reference-density-filename = %s\n",
+                    TestFileManager::getInputFilePath("ellipsoid-density.mrc").c_str() );
+
         //! Mdp values for steepest-decent energy minimization with density fitting values set to non-defaults.
         const std::string mdpDensiftAllDefaultsChanged_ = formatString(
                     "density-guided-simulation-similarity-measure = relative-entropy\n"
@@ -110,7 +123,11 @@ class DensityFittingTest : public MdrunTestFixture
                     "density-guided-simulation-gaussian-transform-spreading-range-in-multiples-of-width = 6\n"
                     "density-guided-simulation-normalize-densities = false\n"
                     );
-
+        //! Set mdp values so that energy calculation interval and density guided simulation interval mismatch.
+        const std::string mdpEnergyAndDensityfittingIntervalMismatch_ = formatString(
+                    "nstcalcenergy = 7\n"
+                    "density-guided-simulation-nst = 3\n"
+                    );
         //! The command line to call mdrun
         CommandLine       commandLineForMdrun_;
 };
@@ -144,6 +161,15 @@ TEST_F(DensityFittingTest, EnergyMinimizationEnergyCorrectForRelativeEntropy)
     checkMdrun(expectedEnergyTermMagnitude);
 }
 
+/* Test that grompp exits with error message if energy evaluation frequencies
+ * do not match.
+ */
+TEST_F(DensityFittingTest, GromppErrorWhenEnergyEvaluationFrequencyMismatch)
+{
+    runner_.useStringAsMdpFile(mdpMdDensfitYesUnsetValues + mdpEnergyAndDensityfittingIntervalMismatch_);
+
+    EXPECT_DEATH_IF_SUPPORTED(runner_.callGrompp(), ".*is not a multiple of density-guided-simulation-nst.*");
+}
 
 } // namespace test
 } // namespace gmx