Generalize IForceProvider
authorTeemu Murtola <teemu.murtola@gmail.com>
Sun, 25 Jun 2017 14:24:01 +0000 (17:24 +0300)
committerErik Lindahl <erik.lindahl@gmail.com>
Thu, 29 Jun 2017 17:46:49 +0000 (19:46 +0200)
- Remove knowledge of individual modules from t_forcerec, and hide the
  number of modules, and whether they contribute to a virial or not,
  behind a generic interface.
- Improve hard-to-understand initialization of separate f_novirsum,
  resolving some TODOs.
- Add some parameters that probably make sense for other modules beyond
  the electric field one.

This also makes the requirement for all modules to have the same
parameters for the calculateForces() call explicit.

Change-Id: I4952515c4b707ba458fd267565fd000532ec281e

13 files changed:
docs/doxygen/lib/mdmodules.md
src/gromacs/applied-forces/electricfield.cpp
src/gromacs/applied-forces/tests/electricfield.cpp
src/gromacs/mdlib/forcerec.cpp
src/gromacs/mdlib/forcerec.h
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdrunutility/mdmodules.cpp
src/gromacs/mdrunutility/mdmodules.h
src/gromacs/mdtypes/forcerec.h
src/gromacs/mdtypes/iforceprovider.cpp [new file with mode: 0644]
src/gromacs/mdtypes/iforceprovider.h
src/gromacs/mdtypes/imdmodule.h
src/programs/mdrun/runner.cpp

index 7cbbd73ad78f592900466ec0737ef3d5b24965d4..761329d48fbdd3d35562ee15d368f507fe57b36c 100644 (file)
@@ -22,8 +22,8 @@ Structure of a module
 ---------------------
 
 Each module implements a factory that returns an instance of gmx::IMDModule.
-This interface has methods that in turn return instances of other interfaces:
-gmx::IMdpOptionProvider, gmx::IMDOutputProvider, and IForceProvider.
+This interface has methods that in turn refer to other interfaces:
+gmx::IMdpOptionProvider, gmx::IMDOutputProvider, and gmx::IForceProvider.
 The module also implements these interfaces (or a subset of them), and code
 outside the module only calls methods in these interfaces.
 
index 5a7e45587e063e9d5adfe53e07144581c9af765a..40da0f1a31e7302cdcf9efe9c2027bf8fa7947a6 100644 (file)
@@ -165,7 +165,13 @@ class ElectricField final : public IMDModule,
         // From IMDModule
         IMdpOptionProvider *mdpOptionProvider() override { return this; }
         IMDOutputProvider *outputProvider() override { return this; }
-        IForceProvider *forceProvider() override { return this; }
+        void initForceProviders(ForceProviders *forceProviders) override
+        {
+            if (isActive())
+            {
+                forceProviders->addForceProviderWithoutVirialContribution(this);
+            }
+        }
 
         // From IMdpOptionProvider
         void initMdpTransform(IKeyValueTreeTransformRules *transform) override;
@@ -177,12 +183,13 @@ class ElectricField final : public IMDModule,
         void finishOutput() override;
 
         // From IForceProvider
-        void initForcerec(t_forcerec *fr) override;
         //! \copydoc IForceProvider::calculateForces()
         void calculateForces(const t_commrec  *cr,
                              const t_mdatoms  *mdatoms,
-                             PaddedRVecVector *force,
-                             double            t) override;
+                             const matrix      box,
+                             double            t,
+                             const rvec       *x,
+                             ArrayRef<RVec>    force) override;
 
     private:
         //! Return whether or not to apply a field
@@ -357,15 +364,6 @@ void ElectricField::finishOutput()
     }
 }
 
-void ElectricField::initForcerec(t_forcerec *fr)
-{
-    if (isActive())
-    {
-        fr->bF_NoVirSum = TRUE;
-        fr->efield      = this;
-    }
-}
-
 real ElectricField::field(int dim, real t) const
 {
     return efield_[dim].evaluate(t);
@@ -386,12 +384,14 @@ void ElectricField::printComponents(double t) const
 
 void ElectricField::calculateForces(const t_commrec  *cr,
                                     const t_mdatoms  *mdatoms,
-                                    PaddedRVecVector *force,
-                                    double            t)
+                                    const matrix      /* box */,
+                                    double            t,
+                                    const rvec        * /* x */,
+                                    ArrayRef<RVec>    force)
 {
     if (isActive())
     {
-        rvec *f = as_rvec_array(force->data());
+        rvec *f = as_rvec_array(force.data());
 
         for (int m = 0; (m < DIM); m++)
         {
index b9c4734945f734dd394cef4d6b3af9d773c42e1f..d621f9012721b4e857a19cdab6961760ea8dc1c4 100644 (file)
@@ -53,6 +53,7 @@
 #include "gromacs/mdtypes/iforceprovider.h"
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/keyvaluetreebuilder.h"
 #include "gromacs/utility/keyvaluetreetransform.h"
 #include "gromacs/utility/real.h"
@@ -72,8 +73,6 @@ namespace
 class ElectricFieldTest : public ::testing::Test
 {
     public:
-        ElectricFieldTest() {}
-
         void test(int  dim,
                   real E0,
                   real omega,
@@ -108,13 +107,11 @@ class ElectricFieldTest : public ::testing::Test
             snew(md.chargeA, md.homenr);
             md.chargeA[0] = 1;
 
-            t_commrec  *cr       = init_commrec();
-            t_forcerec *forcerec = mk_forcerec();
-            module.forceProvider()->initForcerec(forcerec);
-            forcerec->efield->calculateForces(cr, &md, &f, 0);
+            t_commrec  *cr = init_commrec();
+            module.initForceProviders()->calculateForces(cr, &md, nullptr, 0, nullptr,
+                                                         gmx::EmptyArrayRef(), f);
             done_commrec(cr);
             EXPECT_REAL_EQ_TOL(f[0][dim], expectedValue, tolerance);
-            sfree(forcerec);
             sfree(md.chargeA);
         }
 };
index 8982c60a91d72e5ccc7dc1e91fe396a2e5e9b205..f2c3861e5778f6803ee7d34cdbb5f2b56a2be948 100644 (file)
@@ -2284,7 +2284,6 @@ void init_forcerec(FILE                *fp,
                    const gmx::MDLogger &mdlog,
                    t_forcerec          *fr,
                    t_fcdata            *fcd,
-                   IForceProvider      *forceProviders,
                    const t_inputrec    *ir,
                    const gmx_mtop_t    *mtop,
                    const t_commrec     *cr,
@@ -2799,14 +2798,10 @@ void init_forcerec(FILE                *fp,
     }
 
     fr->bF_NoVirSum = (EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype) ||
+                       fr->forceProviders->hasForcesWithoutVirialContribution() ||
                        gmx_mtop_ftype_count(mtop, F_POSRES) > 0 ||
                        gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0);
 
-    /* Initialization call after setting bF_NoVirSum,
-     * since it efield->initForcerec also sets this to true.
-     */
-    forceProviders->initForcerec(fr);
-
     if (fr->bF_NoVirSum)
     {
         fr->forceBufferNoVirialSummation = new PaddedRVecVector;
index 677a98abdadf7a4e151fbd96403634f6503e14d4..5eaafc4e49f25f82d7fd7bcf856309aa160b61ef 100644 (file)
@@ -44,7 +44,6 @@
 #include "gromacs/mdtypes/forcerec.h"
 #include "gromacs/timing/wallcycle.h"
 
-struct IForceProvider;
 struct t_commrec;
 struct t_fcdata;
 struct t_filenm;
@@ -100,7 +99,6 @@ void init_interaction_const_tables(FILE                   *fp,
  * \param[in]  mdlog       File for printing
  * \param[out] fr          The forcerec
  * \param[in]  fcd         Force constant data
- * \param[in]  forceProviders  Handle to modules providing forces
  * \param[in]  ir          Inputrec structure
  * \param[in]  mtop        Molecular topology
  * \param[in]  cr          Communication structures
@@ -116,7 +114,6 @@ void init_forcerec(FILE                   *fplog,
                    const gmx::MDLogger    &mdlog,
                    t_forcerec             *fr,
                    t_fcdata               *fcd,
-                   IForceProvider         *forceProviders,
                    const t_inputrec       *ir,
                    const gmx_mtop_t       *mtop,
                    const t_commrec        *cr,
index cdfc9bc719e13cb9943935b5e181d73be7dc44e0..44ff5d475a22fa6c70ab11f07635a00d20d00f71 100644 (file)
@@ -1387,11 +1387,13 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr,
 
     if (bDoForces)
     {
-        /* Compute forces due to electric field */
-        if (fr->efield != nullptr)
+        /* Collect forces from modules */
+        gmx::ArrayRef<gmx::RVec> fNoVirSum;
+        if (fr->bF_NoVirSum)
         {
-            fr->efield->calculateForces(cr, mdatoms, fr->f_novirsum, t);
+            fNoVirSum = *fr->f_novirsum;
         }
+        fr->forceProviders->calculateForces(cr, mdatoms, box, t, x, *force, fNoVirSum);
 
         /* If we have NoVirSum forces, but we do not calculate the virial,
          * we sum fr->f_novirsum=f later.
@@ -1745,11 +1747,13 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr,
 
     if (bDoForces)
     {
-        /* Compute forces due to electric field */
-        if (fr->efield != nullptr)
+        /* Collect forces from modules */
+        gmx::ArrayRef<gmx::RVec> fNoVirSum;
+        if (fr->bF_NoVirSum)
         {
-            fr->efield->calculateForces(cr, mdatoms, fr->f_novirsum, t);
+            fNoVirSum = *fr->f_novirsum;
         }
+        fr->forceProviders->calculateForces(cr, mdatoms, box, t, x, *force, fNoVirSum);
 
         /* Communicate the forces */
         if (DOMAINDECOMP(cr))
index 24f8450c5c78582cf9146ad60793b825571141c9..cc75cd8190abec10318ceea0ff58931d91221643 100644 (file)
@@ -53,7 +53,7 @@
 namespace gmx
 {
 
-class MDModules::Impl : public IMDOutputProvider, public IForceProvider
+class MDModules::Impl : public IMDOutputProvider
 {
     public:
 
@@ -81,20 +81,8 @@ class MDModules::Impl : public IMDOutputProvider, public IForceProvider
             field_->outputProvider()->finishOutput();
         }
 
-        // From IForceProvider
-        virtual void initForcerec(t_forcerec *fr)
-        {
-            field_->forceProvider()->initForcerec(fr);
-        }
-        virtual void calculateForces(const t_commrec  * /*cr*/,
-                                     const t_mdatoms  * /*mdatoms*/,
-                                     PaddedRVecVector * /*force*/,
-                                     double             /*t*/)
-        {
-            // not called currently
-        }
-
-        std::unique_ptr<IMDModule> field_;
+        std::unique_ptr<IMDModule>      field_;
+        std::unique_ptr<ForceProviders> forceProviders_;
 };
 
 MDModules::MDModules() : impl_(new Impl)
@@ -140,9 +128,13 @@ IMDOutputProvider *MDModules::outputProvider()
     return impl_.get();
 }
 
-IForceProvider *MDModules::forceProvider()
+ForceProviders *MDModules::initForceProviders()
 {
-    return impl_.get();
+    GMX_RELEASE_ASSERT(impl_->forceProviders_ == nullptr,
+                       "Force providers initialized multiple times");
+    impl_->forceProviders_.reset(new ForceProviders);
+    impl_->field_->initForceProviders(impl_->forceProviders_.get());
+    return impl_->forceProviders_.get();
 }
 
 } // namespace gmx
index b625dd26118e2a8fd9c4141084826eee300527e6..663d76ffe77a43a4055c8f10af37cde09aa2e281 100644 (file)
@@ -45,7 +45,8 @@
 
 #include "gromacs/utility/classhelpers.h"
 
-struct IForceProvider;
+struct ForceProviders;
+
 struct t_inputrec;
 
 namespace gmx
@@ -72,15 +73,11 @@ class KeyValueTreeObject;
  *
  * Currently, where the set of modules needs to be accessed, either a pointer
  * to MDModules is passed around, or an instance of IMDOutputProvider or
- * IForceProvider returned from MDModules.  The implementation of these
- * interfaces in MDModules calls the corresponding methods in the relevant
- * modules.  In the future, some additional logic may need to be introduced at
- * the call sites that can also influence the signature of the methods.  In
- * this case, a separate object may need to be introduced (e.g.,
- * ForceProvidersManager or similar) that can be passed around without
- * knowledge of the full MDModules.  t_forcerec also currently directly calls
- * individual modules through pointers to their interfaces, which should be
- * generalized in the future.
+ * ForceProviders returned from MDModules.  These objects returned from
+ * MDModules call the corresponding methods in the relevant modules.
+ * In the future, some additional logic may need to be introduced at
+ * the call sites that can also influence the signature of the methods,
+ * similar to what ForceProviders already does for force computation.
  *
  * The assignOptionsToModules() and adjustInputrecBasedOnModules() methods of
  * this class also take responsibility for wiring up the options (and their
@@ -128,9 +125,9 @@ class MDModules
          */
         IMDOutputProvider *outputProvider();
         /*! \brief
-         * Returns an interface for initializing modules providing forces.
+         * Returns an object for computing forces from the modules.
          */
-        IForceProvider *forceProvider();
+        ForceProviders *initForceProviders();
 
     private:
         class Impl;
index 4673b4c0ebd421a159571a6de711838221a32274..7cd11fad44f1b5c2da7f08bc706a1cf9b0d5227d 100644 (file)
@@ -47,7 +47,7 @@
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/real.h"
 
-struct IForceProvider;
+struct ForceProviders;
 
 /* Abstract type for PME that is defined only in the routine that use them. */
 struct gmx_genborn_t;
@@ -438,7 +438,7 @@ struct t_forcerec {
     int                         nthread_ewc;
     struct ewald_corr_thread_t *ewc_t;
 
-    struct IForceProvider      *efield;
+    struct ForceProviders      *forceProviders;
 };
 
 /* Important: Starting with Gromacs-4.6, the values of c6 and c12 in the nbfp array have
diff --git a/src/gromacs/mdtypes/iforceprovider.cpp b/src/gromacs/mdtypes/iforceprovider.cpp
new file mode 100644 (file)
index 0000000..c6f9c77
--- /dev/null
@@ -0,0 +1,99 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2017, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements classes from iforceprovider.h.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \ingroup module_mdtypes
+ */
+#include "gmxpre.h"
+
+#include "iforceprovider.h"
+
+#include <vector>
+
+#include "gromacs/utility/arrayref.h"
+
+using namespace gmx;
+
+class ForceProviders::Impl
+{
+    public:
+        std::vector<IForceProvider *> withVirialContribution_;
+        std::vector<IForceProvider *> withoutVirialContribution_;
+};
+
+ForceProviders::ForceProviders()
+    : impl_(new Impl)
+{
+}
+
+ForceProviders::~ForceProviders()
+{
+}
+
+void ForceProviders::addForceProvider(gmx::IForceProvider *provider)
+{
+    impl_->withVirialContribution_.push_back(provider);
+}
+
+void ForceProviders::addForceProviderWithoutVirialContribution(gmx::IForceProvider *provider)
+{
+    impl_->withoutVirialContribution_.push_back(provider);
+}
+
+bool ForceProviders::hasForcesWithoutVirialContribution() const
+{
+    return !impl_->withoutVirialContribution_.empty();
+}
+
+void ForceProviders::calculateForces(const t_commrec          *cr,
+                                     const t_mdatoms          *mdatoms,
+                                     const matrix              box,
+                                     double                    t,
+                                     const rvec               *x,
+                                     gmx::ArrayRef<gmx::RVec>  force,
+                                     gmx::ArrayRef<gmx::RVec>  f_novirsum) const
+{
+    for (auto provider : impl_->withVirialContribution_)
+    {
+        provider->calculateForces(cr, mdatoms, box, t, x, force);
+    }
+    for (auto provider : impl_->withoutVirialContribution_)
+    {
+        provider->calculateForces(cr, mdatoms, box, t, x, f_novirsum);
+    }
+}
index 3dd88fb98a01de68d4b1cc7370a07adf5a0e792c..cccc27650ff992e17b2cbbdf0e2d2ca9ff473dbf 100644 (file)
  */
 /*! \libinternal \file
  * \brief
- * Declares gmx::IForceProvider.
+ * Declares gmx::IForceProvider and ForceProviders.
  *
  * See \ref page_mdmodules for an overview of this and associated interfaces.
  *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
  * \inlibraryapi
  * \ingroup module_mdtypes
  */
 #ifndef GMX_MDTYPES_IFORCEPROVIDER_H
 #define GMX_MDTYPES_IFORCEPROVIDER_H
 
-#include "gromacs/math/paddedvector.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/utility/classhelpers.h"
 
 struct t_commrec;
 struct t_forcerec;
 struct t_mdatoms;
 
+namespace gmx
+{
+
+template <typename T>
+class ArrayRef;
+
 /*! \libinternal \brief
  * Interface for a component that provides forces during MD.
  *
- * This is typically part of a larger structure/class managing its own
- * data, such that it has the information on what to do stored locally.
+ * Modules implementing IMDModule generally implement this internally, and use
+ * IMDModule::initForceProviders() to register their implementation in
+ * ForceProviders.
  *
- * The interface is not very generic, as it has been written purely based on
- * extraction of existing functions related to electric field handling.
- * This needs to be generalized when more modules are moved to use the
- * interface.
+ * The interface most likely requires additional generalization for use in
+ * other modules than the current electric field implementation.
  *
  * \inlibraryapi
  * \ingroup module_mdtypes
  */
-struct IForceProvider
+class IForceProvider
 {
     public:
-        /*! \brief
-         * Sets relevant options in the forcerec structure.
-         *
-         * \param[inout] fr The forcerec structure
-         *
-         * \todo
-         * This should be replaced by a method that returns a set of
-         * flags/other options (either here, or where the IForceProvider
-         * instance is returned), and forcerec should be initialized based on
-         * that.
-         */
-        virtual void initForcerec(t_forcerec *fr) = 0;
-
         /*! \brief
          * Computes forces.
          *
          * \param[in]    cr      Communication record for parallel operations
          * \param[in]    mdatoms Atom information
-         * \param[inout] force   The forces
+         * \param[in]    box     The box
          * \param[in]    t       The actual time in the simulation (ps)
+         * \param[in]    x       The coordinates
+         * \param[inout] force   The forces
          */
-        virtual void calculateForces(const t_commrec  *cr,
-                                     const t_mdatoms  *mdatoms,
-                                     PaddedRVecVector *force,
-                                     double            t) = 0;
+        virtual void calculateForces(const t_commrec          *cr,
+                                     const t_mdatoms          *mdatoms,
+                                     const matrix              box,
+                                     double                    t,
+                                     const rvec               *x,
+                                     gmx::ArrayRef<gmx::RVec>  force) = 0;
 
     protected:
         ~IForceProvider() {}
 };
 
+} // namespace gmx
+
+/*! \libinternal \brief
+ * Evaluates forces from a collection of gmx::IForceProvider.
+ *
+ * This class is a `struct` outside the `gmx` namespace to make it possible to
+ * forward-declare it in forcerec.h, which still needs to compile when included
+ * from the C group kernels.
+ *
+ * \inlibraryapi
+ * \ingroup module_mdtypes
+ */
+struct ForceProviders
+{
+    public:
+        ForceProviders();
+        ~ForceProviders();
+
+        /*! \brief
+         * Adds a provider.
+         */
+        void addForceProvider(gmx::IForceProvider *provider);
+        /*! \brief
+         * Adds a provider whose forces should not contribute to the virial.
+         */
+        void addForceProviderWithoutVirialContribution(gmx::IForceProvider *provider);
+
+        //! Whether there are modules that do not contribute to the virial.
+        bool hasForcesWithoutVirialContribution() const;
+
+        //! Computes forces.
+        void calculateForces(const t_commrec          *cr,
+                             const t_mdatoms          *mdatoms,
+                             const matrix              box,
+                             double                    t,
+                             const rvec               *x,
+                             gmx::ArrayRef<gmx::RVec>  force,
+                             gmx::ArrayRef<gmx::RVec>  f_novirsum) const;
+
+    private:
+        class Impl;
+
+        gmx::PrivateImplPointer<Impl> impl_;
+};
+
 #endif
index 5c2fed976e0b4b7fd35abb1bd35d3432c1f9f5c8..94cba248261c4b622f7b3bebf392c7b83e29ac3e 100644 (file)
@@ -44,7 +44,7 @@
 #ifndef GMX_MDTYPES_IMDMODULE_H
 #define GMX_MDTYPES_IMDMODULE_H
 
-struct IForceProvider;
+struct ForceProviders;
 
 namespace gmx
 {
@@ -71,8 +71,8 @@ class IMDModule
         virtual IMdpOptionProvider *mdpOptionProvider() = 0;
         //! Returns an interface for handling output files during simulation.
         virtual IMDOutputProvider *outputProvider()     = 0;
-        //! Returns an interface for computing forces during simulation.
-        virtual IForceProvider *forceProvider()         = 0;
+        //! Initializes force providers from this module.
+        virtual void initForceProviders(ForceProviders *forceProviders) = 0;
 };
 
 } // namespace gmx
index 4aa7d11b81a00ad4a661eac2d54d6171592033dd..a9e19b56e774631012d57d44fb27f803e20d37c0 100644 (file)
@@ -1245,10 +1245,11 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         bcast_state(cr, state);
 
         /* Initiate forcerecord */
-        fr          = mk_forcerec();
-        fr->hwinfo  = hwinfo;
-        fr->gpu_opt = &hw_opt->gpu_opt;
-        init_forcerec(fplog, mdlog, fr, fcd, mdModules.forceProvider(),
+        fr                 = mk_forcerec();
+        fr->hwinfo         = hwinfo;
+        fr->gpu_opt        = &hw_opt->gpu_opt;
+        fr->forceProviders = mdModules.initForceProviders();
+        init_forcerec(fplog, mdlog, fr, fcd,
                       inputrec, mtop, cr, box,
                       opt2fn("-table", nfile, fnm),
                       opt2fn("-tablep", nfile, fnm),