Use WholeMoleculeTransform for orires and epsilon-surface
authorBerk Hess <hess@kth.se>
Thu, 5 Mar 2020 11:47:59 +0000 (12:47 +0100)
committerBerk Hess <hess@kth.se>
Wed, 18 Mar 2020 22:04:44 +0000 (23:04 +0100)
Remove direct use of the graph from the orientation restraint
and Ewald epsilon-surface term code by passing in whole molecules
created by the WholeMolecules class.
Since this was the only remaining use of graph in do_force(),
it can now be removed.

Also enabled the epsilon-surface mdrun-test.

Fixes #3441

Change-Id: Idfec6508c6dfd9e1a656cf23613ede3793794901

13 files changed:
src/gromacs/listed_forces/listed_forces.cpp
src/gromacs/listed_forces/listed_forces.h
src/gromacs/listed_forces/orires.cpp
src/gromacs/listed_forces/orires.h
src/gromacs/mdlib/calcmu.cpp
src/gromacs/mdlib/calcmu.h
src/gromacs/mdlib/force.cpp
src/gromacs/mdlib/force.h
src/gromacs/mdlib/forcerec.cpp
src/gromacs/mdlib/sim_util.cpp
src/gromacs/modularsimulator/modularsimulator.cpp
src/programs/mdrun/tests/ewaldsurfaceterm.cpp
src/programs/mdrun/tests/refdata/EwaldSurfaceTerm_EwaldSurfaceTermTest_WithinTolerances_2.xml [new file with mode: 0644]

index fbfc56927a0c780c8086427a0c9d1571fd6cd61e..b04cfaf93b203dd3426ca447beb9b3e7dc3a3fb6 100644 (file)
@@ -474,11 +474,20 @@ bool haveCpuListedForces(const t_forcerec& fr, const InteractionDefinitions& ide
     return haveCpuBondeds(fr) || haveRestraints(idef, fcd);
 }
 
+namespace
+{
+
+/*! \brief Calculates all listed force interactions.
+ *
+ * Note that pbc_full is used only for position restraints, and is
+ * not initialized if there are none.
+ */
 void calc_listed(const t_commrec*              cr,
                  const gmx_multisim_t*         ms,
                  struct gmx_wallcycle*         wcycle,
                  const InteractionDefinitions& idef,
                  const rvec                    x[],
+                 ArrayRef<const gmx::RVec>     xWholeMolecules,
                  history_t*                    hist,
                  gmx::ForceOutputs*            forceOutputs,
                  const t_forcerec*             fr,
@@ -527,17 +536,10 @@ void calc_listed(const t_commrec*              cr,
         /* Do pre force calculation stuff which might require communication */
         if (fcd->orires.nr > 0)
         {
-            /* This assertion is to ensure we have whole molecules.
-             * Unfortunately we do not have an mdrun state variable that tells
-             * us if molecules in x are not broken over PBC, so we have to make
-             * do with checking graph!=nullptr, which should tell us if we made
-             * molecules whole before calling the current function.
-             */
-            GMX_RELEASE_ASSERT(fr->pbcType == PbcType::No || g != nullptr,
-                               "With orientation restraints molecules should be whole");
+            GMX_ASSERT(!xWholeMolecules.empty(), "Need whole molecules for orienation restraints");
             enerd->term[F_ORIRESDEV] =
                     calc_orires_dev(ms, idef.il[F_ORIRES].size(), idef.il[F_ORIRES].iatoms.data(),
-                                    idef.iparams.data(), md, x, pbc_null, fcd, hist);
+                                    idef.iparams.data(), md, xWholeMolecules, x, pbc_null, fcd, hist);
         }
         if (fcd->disres.nres > 0)
         {
@@ -582,6 +584,11 @@ void calc_listed(const t_commrec*              cr,
     }
 }
 
+/*! \brief As calc_listed(), but only determines the potential energy
+ * for the perturbed interactions.
+ *
+ * The shift forces in fr are not affected.
+ */
 void calc_listed_lambda(const InteractionDefinitions& idef,
                         const rvec                    x[],
                         const t_forcerec*             fr,
@@ -646,25 +653,28 @@ void calc_listed_lambda(const InteractionDefinitions& idef,
     sfree(f);
 }
 
-void do_force_listed(struct gmx_wallcycle*         wcycle,
-                     const matrix                  box,
-                     const t_lambda*               fepvals,
-                     const t_commrec*              cr,
-                     const gmx_multisim_t*         ms,
-                     const InteractionDefinitions& idef,
-                     const rvec                    x[],
-                     history_t*                    hist,
-                     gmx::ForceOutputs*            forceOutputs,
-                     const t_forcerec*             fr,
-                     const struct t_pbc*           pbc,
-                     const struct t_graph*         graph,
-                     gmx_enerdata_t*               enerd,
-                     t_nrnb*                       nrnb,
-                     const real*                   lambda,
-                     const t_mdatoms*              md,
-                     t_fcdata*                     fcd,
-                     int*                          global_atom_index,
-                     const gmx::StepWorkload&      stepWork)
+} // namespace
+
+void do_force_listed(struct gmx_wallcycle*          wcycle,
+                     const matrix                   box,
+                     const t_lambda*                fepvals,
+                     const t_commrec*               cr,
+                     const gmx_multisim_t*          ms,
+                     const InteractionDefinitions&  idef,
+                     const rvec                     x[],
+                     gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
+                     history_t*                     hist,
+                     gmx::ForceOutputs*             forceOutputs,
+                     const t_forcerec*              fr,
+                     const struct t_pbc*            pbc,
+                     const struct t_graph*          graph,
+                     gmx_enerdata_t*                enerd,
+                     t_nrnb*                        nrnb,
+                     const real*                    lambda,
+                     const t_mdatoms*               md,
+                     t_fcdata*                      fcd,
+                     int*                           global_atom_index,
+                     const gmx::StepWorkload&       stepWork)
 {
     t_pbc pbc_full; /* Full PBC is needed for position restraints */
 
@@ -678,8 +688,8 @@ void do_force_listed(struct gmx_wallcycle*         wcycle,
         /* Not enough flops to bother counting */
         set_pbc(&pbc_full, fr->pbcType, box);
     }
-    calc_listed(cr, ms, wcycle, idef, x, hist, forceOutputs, fr, pbc, &pbc_full, graph, enerd, nrnb,
-                lambda, md, fcd, global_atom_index, stepWork);
+    calc_listed(cr, ms, wcycle, idef, x, xWholeMolecules, hist, forceOutputs, fr, pbc, &pbc_full,
+                graph, enerd, nrnb, lambda, md, fcd, global_atom_index, stepWork);
 
     /* Check if we have to determine energy differences
      * at foreign lambda's.
index 479a4d8fff90384c7e1d9c4a28d9212715c14347..c25e9973c141040837a9729fac7eeb493c2d6de7 100644 (file)
@@ -111,68 +111,32 @@ using BondedFunction = real (*)(int              nbonds,
 //! Getter for finding a callable CPU function to compute an \c ftype interaction.
 BondedFunction bondedFunction(int ftype);
 
-/*! \brief Calculates all listed force interactions.
- *
- * Note that pbc_full is used only for position restraints, and is
- * not initialized if there are none. */
-void calc_listed(const t_commrec*              cr,
-                 const gmx_multisim_t*         ms,
-                 struct gmx_wallcycle*         wcycle,
-                 const InteractionDefinitions& idef,
-                 const rvec                    x[],
-                 history_t*                    hist,
-                 gmx::ForceOutputs*            forceOutputs,
-                 const t_forcerec*             fr,
-                 const struct t_pbc*           pbc,
-                 const struct t_pbc*           pbc_full,
-                 const struct t_graph*         g,
-                 gmx_enerdata_t*               enerd,
-                 t_nrnb*                       nrnb,
-                 const real*                   lambda,
-                 const t_mdatoms*              md,
-                 struct t_fcdata*              fcd,
-                 int*                          ddgatindex,
-                 const gmx::StepWorkload&      stepWork);
-
-/*! \brief As calc_listed(), but only determines the potential energy
- * for the perturbed interactions.
- *
- * The shift forces in fr are not affected. */
-void calc_listed_lambda(const InteractionDefinitions& idef,
-                        const rvec                    x[],
-                        const t_forcerec*             fr,
-                        const struct t_pbc*           pbc,
-                        const struct t_graph*         g,
-                        gmx_grppairener_t*            grpp,
-                        real*                         epot,
-                        gmx::ArrayRef<real>           dvdl,
-                        t_nrnb*                       nrnb,
-                        const real*                   lambda,
-                        const t_mdatoms*              md,
-                        struct t_fcdata*              fcd,
-                        int*                          global_atom_index);
-
 /*! \brief Do all aspects of energy and force calculations for mdrun
- * on the set of listed interactions */
-void do_force_listed(struct gmx_wallcycle*         wcycle,
-                     const matrix                  box,
-                     const t_lambda*               fepvals,
-                     const t_commrec*              cr,
-                     const gmx_multisim_t*         ms,
-                     const InteractionDefinitions& idef,
-                     const rvec                    x[],
-                     history_t*                    hist,
-                     gmx::ForceOutputs*            forceOutputs,
-                     const t_forcerec*             fr,
-                     const struct t_pbc*           pbc,
-                     const struct t_graph*         graph,
-                     gmx_enerdata_t*               enerd,
-                     t_nrnb*                       nrnb,
-                     const real*                   lambda,
-                     const t_mdatoms*              md,
-                     struct t_fcdata*              fcd,
-                     int*                          global_atom_index,
-                     const gmx::StepWorkload&      stepWork);
+ * on the set of listed interactions
+ *
+ * xWholeMolecules only needs to contain whole molecules when orientation
+ * restraints need to be computed and can be empty otherwise.
+ */
+void do_force_listed(struct gmx_wallcycle*          wcycle,
+                     const matrix                   box,
+                     const t_lambda*                fepvals,
+                     const t_commrec*               cr,
+                     const gmx_multisim_t*          ms,
+                     const InteractionDefinitions&  idef,
+                     const rvec                     x[],
+                     gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
+                     history_t*                     hist,
+                     gmx::ForceOutputs*             forceOutputs,
+                     const t_forcerec*              fr,
+                     const struct t_pbc*            pbc,
+                     const struct t_graph*          graph,
+                     gmx_enerdata_t*                enerd,
+                     t_nrnb*                        nrnb,
+                     const real*                    lambda,
+                     const t_mdatoms*               md,
+                     struct t_fcdata*               fcd,
+                     int*                           global_atom_index,
+                     const gmx::StepWorkload&       stepWork);
 
 /*! \brief Returns true if there are position, distance or orientation restraints. */
 bool haveRestraints(const InteractionDefinitions& idef, const t_fcdata& fcd);
index 1f9665b62fd040c05c71881b221ad77fc13a0d79..fe375477ab7dd50b607cd6d592219e17473a3e60 100644 (file)
 #include "gromacs/topology/ifunc.h"
 #include "gromacs/topology/mtop_util.h"
 #include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/pleasecite.h"
 #include "gromacs/utility/smalloc.h"
 
+using gmx::ArrayRef;
+using gmx::RVec;
+
 // TODO This implementation of ensemble orientation restraints is nasty because
 // a user can't just do multi-sim with single-sim orientation restraints.
 
@@ -382,6 +386,7 @@ real calc_orires_dev(const gmx_multisim_t* ms,
                      const t_iatom         forceatoms[],
                      const t_iparams       ip[],
                      const t_mdatoms*      md,
+                     ArrayRef<const RVec>  xWholeMolecules,
                      const rvec            x[],
                      const t_pbc*          pbc,
                      t_fcdata*             fcd,
@@ -445,7 +450,7 @@ real calc_orires_dev(const gmx_multisim_t* ms,
     {
         if (md->cORF[i] == 0)
         {
-            copy_rvec(x[i], xtmp[j]);
+            copy_rvec(xWholeMolecules[i], xtmp[j]);
             mref[j] = md->massT[i];
             for (int d = 0; d < DIM; d++)
             {
index 44988ddba17c3b9af117dd8a57442757dd422ed8..f29557986b90842a86d531bc0b4ceadd6954db97 100644 (file)
@@ -59,6 +59,12 @@ struct t_fcdata;
 struct t_oriresdata;
 class t_state;
 
+namespace gmx
+{
+template<typename>
+class ArrayRef;
+} // namespace gmx
+
 /*! \brief
  * Decides whether orientation restraints can work, and initializes
  * all the orientation restraint stuff in *od (and assumes *od is
@@ -80,15 +86,16 @@ void init_orires(FILE*                 fplog,
  *
  * Returns the weighted RMS deviation of the orientation restraints.
  */
-real calc_orires_dev(const gmx_multisim_t* ms,
-                     int                   nfa,
-                     const t_iatom         fa[],
-                     const t_iparams       ip[],
-                     const t_mdatoms*      md,
-                     const rvec            x[],
-                     const t_pbc*          pbc,
-                     t_fcdata*             fcd,
-                     history_t*            hist);
+real calc_orires_dev(const gmx_multisim_t*          ms,
+                     int                            nfa,
+                     const t_iatom                  fa[],
+                     const t_iparams                ip[],
+                     const t_mdatoms*               md,
+                     gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
+                     const rvec                     x[],
+                     const t_pbc*                   pbc,
+                     t_fcdata*                      fcd,
+                     history_t*                     hist);
 
 /*! \brief
  * Diagonalizes the order tensor(s) of the orienation restraints.
index ef6b77786f0a77f38f28045b7f5d4762365d467a..006fe7fffa580f16c8e585707c9b7aa2dd6071de 100644 (file)
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/gmx_omp_nthreads.h"
 
-void calc_mu(int                      start,
-             int                      homenr,
-             gmx::ArrayRef<gmx::RVec> x,
-             const real               q[],
-             const real               qB[],
-             int                      nChargePerturbed,
-             dvec                     mu,
-             dvec                     mu_B)
+void calc_mu(int                            start,
+             int                            homenr,
+             gmx::ArrayRef<const gmx::RVec> x,
+             const real                     q[],
+             const real                     qB[],
+             int                            nChargePerturbed,
+             dvec                           mu,
+             dvec                           mu_B)
 {
     int    end, m;
     double mu_x, mu_y, mu_z;
index 822fe43e30efd9ab95f39c132f7d1841a3626f72..b96f12e86a6d5f995a951945777a837294808dd9 100644 (file)
 #include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/basedefinitions.h"
 
-void calc_mu(int                      start,
-             int                      homenr,
-             gmx::ArrayRef<gmx::RVec> x,
-             const real               q[],
-             const real               qB[],
-             int                      nChargePerturbed,
-             dvec                     mu,
-             dvec                     mu_B);
+void calc_mu(int                            start,
+             int                            homenr,
+             gmx::ArrayRef<const gmx::RVec> x,
+             const real                     q[],
+             const real                     qB[],
+             int                            nChargePerturbed,
+             dvec                           mu,
+             dvec                           mu_B);
 
 gmx_bool read_mu(FILE* fp, rvec mu, real* vol);
 /* Return true on succes */
index 9feb4003e66261bff9f76f0889a35fdce4082e54..fc407a0fe35635cbcfa62291c190696dc47ec3d8 100644 (file)
@@ -108,6 +108,7 @@ void do_force_lowlevel(t_forcerec*                         fr,
                        gmx_wallcycle_t                     wcycle,
                        const t_mdatoms*                    md,
                        gmx::ArrayRefWithPadding<gmx::RVec> coordinates,
+                       gmx::ArrayRef<const gmx::RVec>      xWholeMolecules,
                        history_t*                          hist,
                        gmx::ForceOutputs*                  forceOutputs,
                        gmx_enerdata_t*                     enerd,
@@ -180,8 +181,8 @@ void do_force_lowlevel(t_forcerec*                         fr,
             set_pbc_dd(&pbc, fr->pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : nullptr, TRUE, box);
         }
 
-        do_force_listed(wcycle, box, ir->fepvals, cr, ms, idef, x, hist, forceOutputs, fr, &pbc,
-                        graph, enerd, nrnb, lambda, md, fcd,
+        do_force_listed(wcycle, box, ir->fepvals, cr, ms, idef, x, xWholeMolecules, hist,
+                        forceOutputs, fr, &pbc, graph, enerd, nrnb, lambda, md, fcd,
                         DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr, stepWork);
     }
 
index 17a90bb3bfd020db33af31df2cae6a72f5a3c4de..1c36420304d39375f34898d2d947aad8d48489e8 100644 (file)
@@ -115,6 +115,11 @@ void do_force(FILE*                               log,
  */
 
 
+/* Compute listed forces, Ewald, PME corrections add when (when used).
+ *
+ * xWholeMolecules only needs to contain whole molecules when orientation
+ * restraints need to be computed and can be empty otherwise.
+ */
 void do_force_lowlevel(t_forcerec*                         fr,
                        const t_inputrec*                   ir,
                        const InteractionDefinitions&       idef,
@@ -124,6 +129,7 @@ void do_force_lowlevel(t_forcerec*                         fr,
                        gmx_wallcycle*                      wcycle,
                        const t_mdatoms*                    md,
                        gmx::ArrayRefWithPadding<gmx::RVec> coordinates,
+                       gmx::ArrayRef<const gmx::RVec>      xWholeMolecules,
                        history_t*                          hist,
                        gmx::ForceOutputs*                  forceOutputs,
                        gmx_enerdata_t*                     enerd,
index e7a04a2d28061b82c254446c8ac7551edffc4493..0d1828fe49ab8c51e612a2dc802e3d059ef846d6 100644 (file)
@@ -918,48 +918,6 @@ static void init_interaction_const(FILE*                 fp,
     *interaction_const = ic;
 }
 
-bool areMoleculesDistributedOverPbc(const t_inputrec& ir, const gmx_mtop_t& mtop, const gmx::MDLogger& mdlog)
-{
-    bool areMoleculesDistributedOverPbc = false;
-
-    if (ir.pbcType != PbcType::No)
-    {
-        /* Not making molecules whole is faster in most cases,
-         * but with orientation restraints or non-tinfoil boundary
-         * conditions we need whole molecules.
-         */
-        const bool useEwaldSurfaceCorrection = (EEL_PME_EWALD(ir.coulombtype) && ir.epsilon_surface != 0);
-        areMoleculesDistributedOverPbc =
-                (gmx_mtop_ftype_count(mtop, F_ORIRES) == 0 && !useEwaldSurfaceCorrection);
-
-        if (getenv("GMX_USE_GRAPH") != nullptr)
-        {
-            areMoleculesDistributedOverPbc = false;
-
-            GMX_LOG(mdlog.warning)
-                    .asParagraph()
-                    .appendText(
-                            "GMX_USE_GRAPH is set, using the graph for bonded "
-                            "interactions");
-
-            if (mtop.bIntermolecularInteractions)
-            {
-                GMX_LOG(mdlog.warning)
-                        .asParagraph()
-                        .appendText(
-                                "WARNING: Molecules linked by intermolecular interactions "
-                                "have to reside in the same periodic image, otherwise "
-                                "artifacts will occur!");
-            }
-        }
-
-        GMX_RELEASE_ASSERT(areMoleculesDistributedOverPbc || !mtop.bIntermolecularInteractions,
-                           "We need to use PBC within molecules with inter-molecular interactions");
-    }
-
-    return areMoleculesDistributedOverPbc;
-}
-
 void init_forcerec(FILE*                            fp,
                    const gmx::MDLogger&             mdlog,
                    t_forcerec*                      fr,
@@ -1087,25 +1045,30 @@ void init_forcerec(FILE*                            fp,
     {
         const bool useEwaldSurfaceCorrection =
                 (EEL_PME_EWALD(ir->coulombtype) && ir->epsilon_surface != 0);
+        const bool haveOrientationRestraints = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
         if (!DOMAINDECOMP(cr))
         {
-            fr->bMolPBC = areMoleculesDistributedOverPbc(*ir, *mtop, mdlog);
-            // The assert below is equivalent to fcd->orires.nr != gmx_mtop_ftype_count(mtop, F_ORIRES)
-            GMX_RELEASE_ASSERT(!fr->bMolPBC || fcd->orires.nr == 0,
-                               "Molecules broken over PBC exist in a simulation including "
-                               "orientation restraints. "
-                               "This likely means that the global topology and the force constant "
-                               "data have gotten out of sync.");
+            fr->bMolPBC = true;
+
+            if (useEwaldSurfaceCorrection || haveOrientationRestraints)
+            {
+                fr->wholeMoleculeTransform =
+                        std::make_unique<gmx::WholeMoleculeTransform>(*mtop, ir->pbcType);
+            }
         }
         else
         {
             fr->bMolPBC = dd_bonded_molpbc(cr->dd, fr->pbcType);
 
-            if (useEwaldSurfaceCorrection && !dd_moleculesAreAlwaysWhole(*cr->dd))
+            /* With Ewald surface correction it is useful to support e.g. running water
+             * in parallel with update groups.
+             * With orientation restraints there is no sensible use case supported with DD.
+             */
+            if ((useEwaldSurfaceCorrection && !dd_moleculesAreAlwaysWhole(*cr->dd)) || haveOrientationRestraints)
             {
                 gmx_fatal(FARGS,
-                          "You requested dipole correction (epsilon_surface > 0), but molecules "
-                          "are broken "
+                          "You requested Ewald surface correction or orientation restraints, "
+                          "but molecules are broken "
                           "over periodic boundary conditions by the domain decomposition. "
                           "Run without domain decomposition instead.");
             }
@@ -1113,8 +1076,7 @@ void init_forcerec(FILE*                            fp,
 
         if (useEwaldSurfaceCorrection)
         {
-            GMX_RELEASE_ASSERT((!DOMAINDECOMP(cr) && !fr->bMolPBC)
-                                       || (DOMAINDECOMP(cr) && dd_moleculesAreAlwaysWhole(*cr->dd)),
+            GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr) || dd_moleculesAreAlwaysWhole(*cr->dd),
                                "Molecules can not be broken by PBC with epsilon_surface > 0");
         }
     }
index d2332cb5b44d13040e41b4980ee151d42790b0ba..8cda9b4e9d00e755ac6fb5e105fbb4f635362079 100644 (file)
@@ -1405,6 +1405,12 @@ void do_force(FILE*                               fplog,
         wallcycle_stop(wcycle, ewcLAUNCH_GPU);
     }
 
+    gmx::ArrayRef<const gmx::RVec> xWholeMolecules;
+    if (fr->wholeMoleculeTransform)
+    {
+        xWholeMolecules = fr->wholeMoleculeTransform->wholeMoleculeCoordinates(x.unpaddedArrayRef(), box);
+    }
+
     DipoleData dipoleData;
 
     if (simulationWork.computeMuTot)
@@ -1414,7 +1420,9 @@ void do_force(FILE*                               fplog,
         /* Calculate total (local) dipole moment in a temporary common array.
          * This makes it possible to sum them over nodes faster.
          */
-        calc_mu(start, mdatoms->homenr, x.unpaddedArrayRef(), mdatoms->chargeA, mdatoms->chargeB,
+        gmx::ArrayRef<const gmx::RVec> xRef =
+                (xWholeMolecules.empty() ? x.unpaddedArrayRef() : xWholeMolecules);
+        calc_mu(start, mdatoms->homenr, xRef, mdatoms->chargeA, mdatoms->chargeB,
                 mdatoms->nChargePerturbed, dipoleData.muStaging[0], dipoleData.muStaging[1]);
 
         reduceAndUpdateMuTot(&dipoleData, cr, (fr->efep != efepNO), lambda, muTotal, ddBalanceRegionHandler);
@@ -1531,9 +1539,9 @@ void do_force(FILE*                               fplog,
         stateGpu->waitCoordinatesReadyOnHost(AtomLocality::NonLocal);
     }
     /* Compute the bonded and non-bonded energies and optionally forces */
-    do_force_lowlevel(fr, inputrec, top->idef, cr, ms, nrnb, wcycle, mdatoms, x, hist, &forceOut,
-                      enerd, fcd, box, lambda.data(), graph, as_rvec_array(dipoleData.muStateAB),
-                      stepWork, ddBalanceRegionHandler);
+    do_force_lowlevel(fr, inputrec, top->idef, cr, ms, nrnb, wcycle, mdatoms, x, xWholeMolecules,
+                      hist, &forceOut, enerd, fcd, box, lambda.data(), graph,
+                      as_rvec_array(dipoleData.muStateAB), stepWork, ddBalanceRegionHandler);
 
     wallcycle_stop(wcycle, ewcFORCE);
 
index dfedd40457b5635dda9c5837a46fd9cf3221d359..da948918e76d4b2e788b17f80353610d084fa225 100644 (file)
@@ -934,10 +934,6 @@ bool ModularSimulator::isInputCompatible(bool                             exitOn
             isInputCompatible
             && conditionalAssert(!doMembed,
                                  "Membrane embedding is not supported by the modular simulator.");
-    const bool useGraph = !areMoleculesDistributedOverPbc(*inputrec, globalTopology, MDLogger());
-    isInputCompatible =
-            isInputCompatible
-            && conditionalAssert(!useGraph, "Graph is not supported by the modular simulator.");
     // TODO: Change this to the boolean passed when we merge the user interface change for the GPU update.
     isInputCompatible =
             isInputCompatible
index 0f7735683a9f2a7c3e2ecc9ae3848a82816df619..367d0abb427789d551e6928b69bcf9b1fd0bb880 100644 (file)
@@ -151,8 +151,8 @@ TEST_P(EwaldSurfaceTermTest, WithinTolerances)
         CommandLine mdrunCaller;
         ASSERT_EQ(0, runner_.callMdrun(mdrunCaller));
         EnergyTermsToCompare energyTermsToCompare{
-            { { interaction_function[F_EPOT].longname, relativeToleranceAsFloatingPoint(1, 1e-4) },
-              { interaction_function[F_ETOT].longname, relativeToleranceAsFloatingPoint(1, 1e-4) } }
+            { { interaction_function[F_EPOT].longname, absoluteTolerance(1e-3) },
+              { interaction_function[F_ETOT].longname, absoluteTolerance(1e-3) } }
         };
         TestReferenceData refData;
         auto checker = refData.rootChecker().checkCompound("Simulation", simulationName);
@@ -178,9 +178,7 @@ TEST_P(EwaldSurfaceTermTest, WithinTolerances)
 
 //! Containers of systems to test.
 //! \{
-// Enable when the epsilon-surface bug is fixed
-//    std::vector<std::string> surfaceTerm = { "3DC", "epsilon-surface-constraint", "epsilon-surface" };
-std::vector<std::string> surfaceTerm = { "3DC", "epsilon-surface-constraint" };
+std::vector<std::string> surfaceTerm = { "3DC", "epsilon-surface-constraint", "epsilon-surface" };
 //! \}
 
 INSTANTIATE_TEST_CASE_P(EwaldSurfaceTerm, EwaldSurfaceTermTest, ::testing::ValuesIn(surfaceTerm));
diff --git a/src/programs/mdrun/tests/refdata/EwaldSurfaceTerm_EwaldSurfaceTermTest_WithinTolerances_2.xml b/src/programs/mdrun/tests/refdata/EwaldSurfaceTerm_EwaldSurfaceTermTest_WithinTolerances_2.xml
new file mode 100644 (file)
index 0000000..1ece247
--- /dev/null
@@ -0,0 +1,142 @@
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+  <Simulation Name="epsilon-surface">
+    <Energy Name="Total Energy">
+      <Real Name="Time 0.000000 Step 0 in frame 0">20.32147</Real>
+      <Real Name="Time 0.010000 Step 4 in frame 1">20.300804</Real>
+      <Real Name="Time 0.020000 Step 8 in frame 2">20.259245</Real>
+      <Real Name="Time 0.030000 Step 12 in frame 3">20.223288</Real>
+      <Real Name="Time 0.040000 Step 16 in frame 4">20.215561</Real>
+      <Real Name="Time 0.050000 Step 20 in frame 5">20.240704</Real>
+    </Energy>
+    <Energy Name="Potential">
+      <Real Name="Time 0.000000 Step 0 in frame 0">20.214478</Real>
+      <Real Name="Time 0.010000 Step 4 in frame 1">16.217239</Real>
+      <Real Name="Time 0.020000 Step 8 in frame 2">8.2757568</Real>
+      <Real Name="Time 0.030000 Step 12 in frame 3">1.3982239</Real>
+      <Real Name="Time 0.040000 Step 16 in frame 4">-0.07736969</Real>
+      <Real Name="Time 0.050000 Step 20 in frame 5">4.7798805</Real>
+    </Energy>
+    <Vector Name="Time 0.000000 Step 0 F[0]">
+      <Real Name="X">1.6004291</Real>
+      <Real Name="Y">522.04138</Real>
+      <Real Name="Z">370.92786</Real>
+    </Vector>
+    <Vector Name="Time 0.000000 Step 0 F[1]">
+      <Real Name="X">-1.900637</Real>
+      <Real Name="Y">-522.46527</Real>
+      <Real Name="Z">-371.22949</Real>
+    </Vector>
+    <Vector Name="Time 0.000000 Step 0 F[2]">
+      <Real Name="X">20.091843</Real>
+      <Real Name="Y">-1.0637188</Real>
+      <Real Name="Z">18.672668</Real>
+    </Vector>
+    <Vector Name="Time 0.000000 Step 0 F[3]">
+      <Real Name="X">-19.791595</Real>
+      <Real Name="Y">1.4878402</Real>
+      <Real Name="Z">-18.371048</Real>
+    </Vector>
+    <Vector Name="Time 0.010000 Step 4 F[0]">
+      <Real Name="X">1.5716085</Real>
+      <Real Name="Y">468.56262</Real>
+      <Real Name="Z">332.87283</Real>
+    </Vector>
+    <Vector Name="Time 0.010000 Step 4 F[1]">
+      <Real Name="X">-1.8655734</Real>
+      <Real Name="Y">-468.97778</Real>
+      <Real Name="Z">-333.16809</Real>
+    </Vector>
+    <Vector Name="Time 0.010000 Step 4 F[2]">
+      <Real Name="X">18.114136</Real>
+      <Real Name="Y">-1.0417742</Real>
+      <Real Name="Z">16.720505</Real>
+    </Vector>
+    <Vector Name="Time 0.010000 Step 4 F[3]">
+      <Real Name="X">-17.82019</Real>
+      <Real Name="Y">1.4571649</Real>
+      <Real Name="Z">-16.425293</Real>
+    </Vector>
+    <Vector Name="Time 0.020000 Step 8 F[0]">
+      <Real Name="X">1.5227978</Real>
+      <Real Name="Y">338.07782</Real>
+      <Real Name="Z">239.98503</Real>
+    </Vector>
+    <Vector Name="Time 0.020000 Step 8 F[1]">
+      <Real Name="X">-1.801703</Real>
+      <Real Name="Y">-338.47165</Real>
+      <Real Name="Z">-240.26477</Real>
+    </Vector>
+    <Vector Name="Time 0.020000 Step 8 F[2]">
+      <Real Name="X">13.286667</Real>
+      <Real Name="Y">-0.99115908</Real>
+      <Real Name="Z">11.956177</Real>
+    </Vector>
+    <Vector Name="Time 0.020000 Step 8 F[3]">
+      <Real Name="X">-13.007904</Real>
+      <Real Name="Y">1.3853663</Real>
+      <Real Name="Z">-11.676315</Real>
+    </Vector>
+    <Vector Name="Time 0.030000 Step 12 F[0]">
+      <Real Name="X">1.5251596</Real>
+      <Real Name="Y">152.05093</Real>
+      <Real Name="Z">107.45719</Real>
+    </Vector>
+    <Vector Name="Time 0.030000 Step 12 F[1]">
+      <Real Name="X">-1.7829351</Real>
+      <Real Name="Y">-152.41469</Real>
+      <Real Name="Z">-107.71528</Real>
+    </Vector>
+    <Vector Name="Time 0.030000 Step 12 F[2]">
+      <Real Name="X">6.4033661</Real>
+      <Real Name="Y">-0.92074591</Real>
+      <Real Name="Z">5.1584015</Real>
+    </Vector>
+    <Vector Name="Time 0.030000 Step 12 F[3]">
+      <Real Name="X">-6.1460266</Real>
+      <Real Name="Y">1.2851006</Real>
+      <Real Name="Z">-4.9001007</Real>
+    </Vector>
+    <Vector Name="Time 0.040000 Step 16 F[0]">
+      <Real Name="X">1.6834296</Real>
+      <Real Name="Y">-58.880402</Real>
+      <Real Name="Z">-43.009636</Real>
+    </Vector>
+    <Vector Name="Time 0.040000 Step 16 F[1]">
+      <Real Name="X">-1.9177186</Real>
+      <Real Name="Y">58.550491</Real>
+      <Real Name="Z">42.775757</Real>
+    </Vector>
+    <Vector Name="Time 0.040000 Step 16 F[2]">
+      <Real Name="X">-1.4052582</Real>
+      <Real Name="Y">-0.84537828</Real>
+      <Real Name="Z">-2.5614777</Real>
+    </Vector>
+    <Vector Name="Time 0.040000 Step 16 F[3]">
+      <Real Name="X">1.6388092</Real>
+      <Real Name="Y">1.1762809</Real>
+      <Real Name="Z">2.7957458</Real>
+    </Vector>
+    <Vector Name="Time 0.050000 Step 20 F[0]">
+      <Real Name="X">2.1072018</Real>
+      <Real Name="Y">-259.95499</Real>
+      <Real Name="Z">-186.73514</Real>
+    </Vector>
+    <Vector Name="Time 0.050000 Step 20 F[1]">
+      <Real Name="X">-2.3196166</Real>
+      <Real Name="Y">259.6571</Real>
+      <Real Name="Z">186.52396</Real>
+    </Vector>
+    <Vector Name="Time 0.050000 Step 20 F[2]">
+      <Real Name="X">-8.8654175</Real>
+      <Real Name="Y">-0.77974892</Real>
+      <Real Name="Z">-9.9438782</Real>
+    </Vector>
+    <Vector Name="Time 0.050000 Step 20 F[3]">
+      <Real Name="X">9.0768127</Real>
+      <Real Name="Y">1.0791473</Real>
+      <Real Name="Z">10.155777</Real>
+    </Vector>
+  </Simulation>
+</ReferenceData>