Enable splitting of listed interaction calculation
authorBerk Hess <hess@kth.se>
Wed, 2 Sep 2020 15:27:23 +0000 (15:27 +0000)
committerPaul Bauer <paul.bauer.q@gmail.com>
Wed, 2 Sep 2020 15:27:23 +0000 (15:27 +0000)
This replaces the ListedForces unique_ptr in t_forcerec by a vector
of unique_ptr and adds a selection mechanism for bonded interactions.

Also moved the t_fcdata storage out of ListedForces into t_forcerec
because we need only one instance.

17 files changed:
src/gromacs/domdec/mdsetup.cpp
src/gromacs/listed_forces/listed_forces.cpp
src/gromacs/listed_forces/listed_forces.h
src/gromacs/mdlib/force.cpp
src/gromacs/mdlib/forcerec.cpp
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdrun/md.cpp
src/gromacs/mdrun/mimic.cpp
src/gromacs/mdrun/minimize.cpp
src/gromacs/mdrun/rerun.cpp
src/gromacs/mdrun/runner.cpp
src/gromacs/mdtypes/forcerec.h
src/gromacs/mdtypes/simulation_workload.h
src/gromacs/modularsimulator/modularsimulator.cpp
src/gromacs/modularsimulator/simulatoralgorithm.cpp
src/gromacs/topology/ifunc.cpp
src/gromacs/topology/ifunc.h

index 6a8559217f633a7300ebf1a9162cba2a4b3b87d6..90f28d18a77bc64e8ea5fcc0c6f6c39caf8e47a0 100644 (file)
@@ -127,7 +127,10 @@ void mdAlgorithmsSetupAtomData(const t_commrec*     cr,
         make_local_shells(cr, mdatoms, shellfc);
     }
 
-    fr->listedForces->setup(top->idef, fr->natoms_force, fr->gpuBonded != nullptr);
+    for (auto& listedForces : fr->listedForces)
+    {
+        listedForces.setup(top->idef, fr->natoms_force, fr->gpuBonded != nullptr);
+    }
 
     if (EEL_PME(fr->ic->eeltype) && (cr->duty & DUTY_PME))
     {
index c35a06a895fdf25fb27a2df3a8b90a25f60e6294..9548ba44338c6ab72491e053ec5c203f43f78c18 100644 (file)
 #include "manage_threading.h"
 #include "utilities.h"
 
-ListedForces::ListedForces(const int numEnergyGroups, const int numThreads, FILE* fplog) :
+ListedForces::ListedForces(const gmx_ffparams_t&      ffparams,
+                           const int                  numEnergyGroups,
+                           const int                  numThreads,
+                           const InteractionSelection interactionSelection,
+                           FILE*                      fplog) :
+    idefSelection_(ffparams),
     threading_(std::make_unique<bonded_threading_t>(numThreads, numEnergyGroups, fplog)),
-    fcdata_(std::make_unique<t_fcdata>())
+    interactionSelection_(interactionSelection)
 {
 }
 
+ListedForces::ListedForces(ListedForces&& o) noexcept = default;
+
 ListedForces::~ListedForces() = default;
 
-void ListedForces::setup(const InteractionDefinitions& idef, const int numAtomsForce, const bool useGpu)
+//! Copies the selection interactions from \p idefSrc to \p idef
+static void selectInteractions(InteractionDefinitions*                  idef,
+                               const InteractionDefinitions&            idefSrc,
+                               const ListedForces::InteractionSelection interactionSelection)
 {
-    idef_ = &idef;
+    const bool selectPairs =
+            interactionSelection.test(static_cast<int>(ListedForces::InteractionGroup::Pairs));
+    const bool selectDihedrals =
+            interactionSelection.test(static_cast<int>(ListedForces::InteractionGroup::Dihedrals));
+    const bool selectRest =
+            interactionSelection.test(static_cast<int>(ListedForces::InteractionGroup::Rest));
+
+    for (int ftype = 0; ftype < F_NRE; ftype++)
+    {
+        const t_interaction_function& ifunc = interaction_function[ftype];
+        if (ifunc.flags & IF_BOND)
+        {
+            bool assign = false;
+            if (ifunc.flags & IF_PAIR)
+            {
+                assign = selectPairs;
+            }
+            else if (ifunc.flags & IF_DIHEDRAL)
+            {
+                assign = selectDihedrals;
+            }
+            else
+            {
+                assign = selectRest;
+            }
+            if (assign)
+            {
+                idef->il[ftype] = idefSrc.il[ftype];
+            }
+            else
+            {
+                idef->il[ftype].clear();
+            }
+        }
+    }
+}
+
+void ListedForces::setup(const InteractionDefinitions& domainIdef, const int numAtomsForce, const bool useGpu)
+{
+    if (interactionSelection_.all())
+    {
+        // Avoid the overhead of copying all interaction lists by simply setting the reference to the domain idef
+        idef_ = &domainIdef;
+    }
+    else
+    {
+        idef_ = &idefSelection_;
+
+        selectInteractions(&idefSelection_, domainIdef, interactionSelection_);
+
+        idefSelection_.ilsort = domainIdef.ilsort;
+    }
 
     setup_bonded_threading(threading_.get(), numAtomsForce, useGpu, *idef_);
 
@@ -478,13 +539,12 @@ static void calcBondedForces(const InteractionDefinitions& idef,
     }
 }
 
-bool ListedForces::haveRestraints() const
+bool ListedForces::haveRestraints(const t_fcdata& fcdata) const
 {
-    GMX_ASSERT(fcdata_, "Need valid fcdata");
-    GMX_ASSERT(fcdata_->orires && fcdata_->disres, "NMR restraints objects should be set up");
+    GMX_ASSERT(fcdata.orires && fcdata.disres, "NMR restraints objects should be set up");
 
-    return (!idef_->il[F_POSRES].empty() || !idef_->il[F_FBPOSRES].empty()
-            || fcdata_->orires->nr > 0 || fcdata_->disres->nres > 0);
+    return (!idef_->il[F_POSRES].empty() || !idef_->il[F_FBPOSRES].empty() || fcdata.orires->nr > 0
+            || fcdata.disres->nres > 0);
 }
 
 bool ListedForces::haveCpuBondeds() const
@@ -492,9 +552,9 @@ bool ListedForces::haveCpuBondeds() const
     return threading_->haveBondeds;
 }
 
-bool ListedForces::haveCpuListedForces() const
+bool ListedForces::haveCpuListedForces(const t_fcdata& fcdata) const
 {
-    return haveCpuBondeds() || haveRestraints();
+    return haveCpuBondeds() || haveRestraints(fcdata);
 }
 
 namespace
@@ -625,6 +685,7 @@ void ListedForces::calculate(struct gmx_wallcycle*          wcycle,
                              const gmx_multisim_t*          ms,
                              const rvec                     x[],
                              gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
+                             t_fcdata*                      fcdata,
                              history_t*                     hist,
                              gmx::ForceOutputs*             forceOutputs,
                              const t_forcerec*              fr,
@@ -636,16 +697,15 @@ void ListedForces::calculate(struct gmx_wallcycle*          wcycle,
                              int*                           global_atom_index,
                              const gmx::StepWorkload&       stepWork)
 {
-    if (!stepWork.computeListedForces)
+    if (interactionSelection_.none() || !stepWork.computeListedForces)
     {
         return;
     }
 
-    const InteractionDefinitions& idef   = *idef_;
-    t_fcdata&                     fcdata = *fcdata_;
+    const InteractionDefinitions& idef = *idef_;
 
     t_pbc pbc_full; /* Full PBC is needed for position restraints */
-    if (haveRestraints())
+    if (haveRestraints(*fcdata))
     {
         if (!idef.il[F_POSRES].empty() || !idef.il[F_FBPOSRES].empty())
         {
@@ -671,24 +731,24 @@ void ListedForces::calculate(struct gmx_wallcycle*          wcycle,
         }
 
         /* Do pre force calculation stuff which might require communication */
-        if (fcdata.orires->nr > 0)
+        if (fcdata->orires->nr > 0)
         {
             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, xWholeMolecules, x, fr->bMolPBC ? pbc : nullptr, fcdata.orires, hist);
+                    md, xWholeMolecules, x, fr->bMolPBC ? pbc : nullptr, fcdata->orires, hist);
         }
-        if (fcdata.disres->nres > 0)
+        if (fcdata->disres->nres > 0)
         {
             calc_disres_R_6(cr, ms, idef.il[F_DISRES].size(), idef.il[F_DISRES].iatoms.data(), x,
-                            fr->bMolPBC ? pbc : nullptr, fcdata.disres, hist);
+                            fr->bMolPBC ? pbc : nullptr, fcdata->disres, hist);
         }
 
         wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
     }
 
     calc_listed(wcycle, idef, threading_.get(), x, forceOutputs, fr, pbc, enerd, nrnb, lambda, md,
-                &fcdata, global_atom_index, stepWork);
+                fcdata, global_atom_index, stepWork);
 
     /* Check if we have to determine energy differences
      * at foreign lambda's.
@@ -718,7 +778,7 @@ void ListedForces::calculate(struct gmx_wallcycle*          wcycle,
                 }
                 calc_listed_lambda(idef, threading_.get(), x, fr, pbc, forceBufferLambda_,
                                    shiftForceBufferLambda_, &(enerd->foreign_grpp), enerd->foreign_term,
-                                   dvdl, nrnb, lam_i, md, &fcdata, global_atom_index);
+                                   dvdl, nrnb, lam_i, md, fcdata, global_atom_index);
                 sum_epot(enerd->foreign_grpp, enerd->foreign_term);
                 const double dvdlSum = std::accumulate(std::begin(dvdl), std::end(dvdl), 0.);
                 std::fill(std::begin(dvdl), std::end(dvdl), 0.0);
index d821973fc3eb62950cebc04ecb194ca905a034c6..51df97c71a1fb20403ee2efde6a366c776e63202 100644 (file)
 
 #include <memory>
 
+#include <bitset>
+
 #include "gromacs/math/vectypes.h"
+#include "gromacs/topology/idef.h"
 #include "gromacs/topology/ifunc.h"
 #include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/basedefinitions.h"
@@ -84,7 +87,6 @@ struct gmx_grppairener_t;
 struct gmx_localtop_t;
 struct gmx_multisim_t;
 class history_t;
-class InteractionDefinitions;
 struct t_commrec;
 struct t_fcdata;
 struct t_forcerec;
@@ -120,28 +122,59 @@ BondedFunction bondedFunction(int ftype);
 
 /*! \libinternal
  * \brief Class for calculating listed interactions, uses OpenMP parallelization
+ *
+ * Listed interactions can be divided over multiple instances of ListedForces
+ * using the selection flags passed to the constructor.
  */
 class ListedForces
 {
 public:
+    //! Enum for selecting groups of listed interaction types
+    enum class InteractionGroup : int
+    {
+        Pairs,     //!< Pair interactions
+        Dihedrals, //!< Dihedrals, including cmap
+        Rest,      //!< All listed interactions that are not any of the above
+        Count      //!< The number of items above
+    };
+
+    //! Type for specifying selections of groups of interaction types
+    using InteractionSelection = std::bitset<static_cast<int>(InteractionGroup::Count)>;
+
+    //! Returns a selection with all listed interaction types selected
+    static InteractionSelection interactionSelectionAll()
+    {
+        InteractionSelection is;
+        return is.flip();
+    }
+
     /*! \brief Constructor
      *
+     * \param[in] ffparams         The force field parameters
      * \param[in] numEnergyGroups  The number of energy groups, used for storage of pair energies
      * \param[in] numThreads       The number of threads used for computed listed interactions
+     * \param[in] interactionSelection  Select of interaction groups through bits set
      * \param[in] fplog            Log file for printing env.var. override, can be nullptr
      */
-    ListedForces(int numEnergyGroups, int numThreads, FILE* fplog);
+    ListedForces(const gmx_ffparams_t& ffparams,
+                 int                   numEnergyGroups,
+                 int                   numThreads,
+                 InteractionSelection  interactionSelection,
+                 FILE*                 fplog);
+
+    //! Move constructor, default, but in the source file to hide implementation classes
+    ListedForces(ListedForces&& o) noexcept;
 
     //! Destructor which is actually default but in the source file to hide implementation classes
     ~ListedForces();
 
     /*! \brief Copy the listed interactions from \p idef and set up the thread parallelization
      *
-     * \param[in] idef           The idef with all listed interactions to be computed on this rank
+     * \param[in] domainIdef     Interaction definitions for all listed interactions to be computed on this domain/rank
      * \param[in] numAtomsForce  Force are, potentially, computed for atoms 0 to \p numAtomsForce
      * \param[in] useGpu         Whether a GPU is used to compute (part of) the listed interactions
      */
-    void setup(const InteractionDefinitions& idef, int numAtomsForce, bool useGpu);
+    void setup(const InteractionDefinitions& domainIdef, int numAtomsForce, bool useGpu);
 
     /*! \brief Do all aspects of energy and force calculations for mdrun
      * on the set of listed interactions
@@ -156,6 +189,7 @@ public:
                    const gmx_multisim_t*          ms,
                    const rvec                     x[],
                    gmx::ArrayRef<const gmx::RVec> xWholeMolecules,
+                   t_fcdata*                      fcdata,
                    history_t*                     hist,
                    gmx::ForceOutputs*             forceOutputs,
                    const t_forcerec*              fr,
@@ -175,30 +209,26 @@ public:
      * NOTE: the current implementation returns true if there are position restraints
      * or any bonded interactions computed on the CPU.
      */
-    bool haveCpuListedForces() const;
+    bool haveCpuListedForces(const t_fcdata& fcdata) const;
 
     //! Returns true if there are position, distance or orientation restraints
-    bool haveRestraints() const;
-
-    //! Returns a reference to the force calculation data
-    const t_fcdata& fcdata() const { return *fcdata_; }
-
-    //! Returns a reference to the force calculation data
-    t_fcdata& fcdata() { return *fcdata_; }
+    bool haveRestraints(const t_fcdata& fcdata) const;
 
 private:
-    //! The interaction definitions
+    //! Pointer to the interaction definitions
     InteractionDefinitions const* idef_ = nullptr;
+    //! Interaction defintions used for storing selections
+    InteractionDefinitions idefSelection_;
     //! Thread parallelization setup, unique_ptr to avoid declaring bonded_threading_t
     std::unique_ptr<bonded_threading_t> threading_;
-    //! Data for bonded tables and NMR restraining, needs to be refactored
-    std::unique_ptr<t_fcdata> fcdata_;
+    //! Tells which interactions to select for computation
+    const InteractionSelection interactionSelection_;
     //! Force buffer for free-energy forces
     std::vector<real> forceBufferLambda_;
     //! Shift force buffer for free-energy forces
     std::vector<gmx::RVec> shiftForceBufferLambda_;
 
-    GMX_DISALLOW_COPY_MOVE_AND_ASSIGN(ListedForces);
+    GMX_DISALLOW_COPY_AND_ASSIGN(ListedForces);
 };
 
 #endif
index 6805a794549d574783cf42fdfc86f71d70627980..1aa00226b63916e5e0daaf1d2ff26a1a28d1db5d 100644 (file)
@@ -138,9 +138,9 @@ void do_force_lowlevel(t_forcerec*                          fr,
         t_pbc pbc;
 
         /* Check whether we need to take into account PBC in listed interactions. */
-        ListedForces& listedForces = *fr->listedForces;
-        const auto    needPbcForListedForces =
-                fr->bMolPBC && stepWork.computeListedForces && listedForces.haveCpuListedForces();
+        ListedForces& listedForces           = fr->listedForces[0];
+        const auto    needPbcForListedForces = fr->bMolPBC && stepWork.computeListedForces
+                                            && listedForces.haveCpuListedForces(*fr->fcdata);
         if (needPbcForListedForces)
         {
             /* Since all atoms are in the rectangular or triclinic unit-cell,
@@ -149,8 +149,8 @@ void do_force_lowlevel(t_forcerec*                          fr,
             set_pbc_dd(&pbc, fr->pbcType, DOMAINDECOMP(cr) ? cr->dd->numCells : nullptr, TRUE, box);
         }
 
-        listedForces.calculate(wcycle, box, ir->fepvals, cr, ms, x, xWholeMolecules, hist,
-                               forceOutputs, fr, &pbc, enerd, nrnb, lambda, md,
+        listedForces.calculate(wcycle, box, ir->fepvals, cr, ms, x, xWholeMolecules, fr->fcdata.get(),
+                               hist, forceOutputs, fr, &pbc, enerd, nrnb, lambda, md,
                                DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr, stepWork);
     }
 
index 7e7ab9702d55f10dfefa515616443c6feffce856..54446643e3454171724a363db248fc381750dbd8 100644 (file)
@@ -1259,18 +1259,16 @@ void init_forcerec(FILE*                            fp,
         make_wall_tables(fp, ir, tabfn, &mtop->groups, fr);
     }
 
-    /* Initialize the thread working data for bonded interactions */
-    fr->listedForces = std::make_unique<ListedForces>(
-            mtop->groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
-            gmx_omp_nthreads_get(emntBonded), fp);
+    fr->fcdata = std::make_unique<t_fcdata>();
 
     if (!tabbfnm.empty())
     {
-        t_fcdata& fcdata = fr->listedForces->fcdata();
+        t_fcdata& fcdata = *fr->fcdata;
         // Need to catch std::bad_alloc
         // TODO Don't need to catch this here, when merging with master branch
         try
         {
+            //! \todo move these tables into a separate struct and store reference in ListedForces
             fcdata.bondtab  = make_bonded_tables(fp, F_TABBONDS, F_TABBONDSNC, mtop, tabbfnm, "b");
             fcdata.angletab = make_bonded_tables(fp, F_TABANGLES, -1, mtop, tabbfnm, "a");
             fcdata.dihtab   = make_bonded_tables(fp, F_TABDIHS, -1, mtop, tabbfnm, "d");
@@ -1287,6 +1285,11 @@ void init_forcerec(FILE*                            fp,
         }
     }
 
+    /* Initialize the thread working data for bonded interactions */
+    fr->listedForces.emplace_back(
+            mtop->ffparams, mtop->groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
+            gmx_omp_nthreads_get(emntBonded), ListedForces::interactionSelectionAll(), fp);
+
     // QM/MM initialization if requested
     if (ir->bQMMM)
     {
index dda01b711fd0cf8420644e157c3240f5fd4351bf..7b0c7c7e9ccbca53da78415ea1d30c03d4500c90 100644 (file)
@@ -851,10 +851,20 @@ static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec&
     // Note that haveSpecialForces is constant over the whole run
     domainWork.haveSpecialForces =
             haveSpecialForces(inputrec, *fr.forceProviders, pull_work, stepWork.computeForces, ed);
-    domainWork.haveCpuBondedWork  = fr.listedForces->haveCpuBondeds();
-    domainWork.haveGpuBondedWork  = ((fr.gpuBonded != nullptr) && fr.gpuBonded->haveInteractions());
-    domainWork.haveRestraintsWork = fr.listedForces->haveRestraints();
-    domainWork.haveCpuListedForceWork = fr.listedForces->haveCpuListedForces();
+    domainWork.haveCpuListedForceWork = false;
+    domainWork.haveCpuBondedWork      = false;
+    for (const auto& listedForces : fr.listedForces)
+    {
+        if (listedForces.haveCpuListedForces(*fr.fcdata))
+        {
+            domainWork.haveCpuListedForceWork = true;
+        }
+        if (listedForces.haveCpuBondeds())
+        {
+            domainWork.haveCpuBondedWork = true;
+        }
+    }
+    domainWork.haveGpuBondedWork = ((fr.gpuBonded != nullptr) && fr.gpuBonded->haveInteractions());
     // Note that haveFreeEnergyWork is constant over the whole run
     domainWork.haveFreeEnergyWork = (fr.efep != efepNO && mdatoms.nPerturbed != 0);
     // We assume we have local force work if there are CPU
index ac08b56afc0c24752e833e6f335b619a6110c8a0..e924f43e4deb4a4e22c399516eb5b78d66360846 100644 (file)
@@ -260,7 +260,7 @@ void gmx::LegacySimulator::do_md()
     const bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
     const bool useReplicaExchange   = (replExParams.exchangeInterval > 0);
 
-    const t_fcdata& fcdata = fr->listedForces->fcdata();
+    const t_fcdata& fcdata = *fr->fcdata;
 
     bool simulationsShareState = false;
     int  nstSignalComm         = nstglobalcomm;
@@ -1530,7 +1530,7 @@ void gmx::LegacySimulator::do_md()
             {
                 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
                                                    do_log ? fplog : nullptr, step, t,
-                                                   &fr->listedForces->fcdata(), awh.get());
+                                                   fr->fcdata.get(), awh.get());
             }
 
             if (ir->bPull)
index 29d846b8515ac559ec732ed34b35034e8a2adc79..188522b43f949459f64aaa2645d5daeee3288758 100644 (file)
@@ -527,8 +527,7 @@ void gmx::LegacySimulator::do_mimic()
 
             EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
             energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
-                                               do_log ? fplog : nullptr, step, t,
-                                               &fr->listedForces->fcdata(), awh);
+                                               do_log ? fplog : nullptr, step, t, fr->fcdata.get(), awh);
 
             if (do_per_step(step, ir->nstlog))
             {
index 7e6b030e3c9b177fa57966e13247a53da7d64017..e4f94ab0ea451877d2b5b858f4bae508ff0bc20f 100644 (file)
@@ -1128,7 +1128,7 @@ void LegacySimulator::do_cg()
 
         EnergyOutput::printHeader(fplog, step, step);
         energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step,
-                                           step, &fr->listedForces->fcdata(), nullptr);
+                                           step, fr->fcdata.get(), nullptr);
     }
 
     /* Estimate/guess the initial stepsize */
@@ -1546,7 +1546,7 @@ void LegacySimulator::do_cg()
             }
             energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
                                                do_log ? fplog : nullptr, step, step,
-                                               &fr->listedForces->fcdata(), nullptr);
+                                               fr->fcdata.get(), nullptr);
         }
 
         /* Send energies and positions to the IMD client if bIMD is TRUE. */
@@ -1590,7 +1590,7 @@ void LegacySimulator::do_cg()
             /* Write final energy file entries */
             energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
                                                !do_log ? fplog : nullptr, step, step,
-                                               &fr->listedForces->fcdata(), nullptr);
+                                               fr->fcdata.get(), nullptr);
         }
     }
 
@@ -1778,7 +1778,7 @@ void LegacySimulator::do_lbfgs()
 
         EnergyOutput::printHeader(fplog, step, step);
         energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step,
-                                           step, &fr->listedForces->fcdata(), nullptr);
+                                           step, fr->fcdata.get(), nullptr);
     }
 
     /* Set the initial step.
@@ -2271,7 +2271,7 @@ void LegacySimulator::do_lbfgs()
             }
             energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
                                                do_log ? fplog : nullptr, step, step,
-                                               &fr->listedForces->fcdata(), nullptr);
+                                               fr->fcdata.get(), nullptr);
         }
 
         /* Send x and E to IMD client, if bIMD is TRUE. */
@@ -2313,8 +2313,8 @@ void LegacySimulator::do_lbfgs()
     if (!do_ene || !do_log) /* Write final energy file entries */
     {
         energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
-                                           !do_log ? fplog : nullptr, step, step,
-                                           &fr->listedForces->fcdata(), nullptr);
+                                           !do_log ? fplog : nullptr, step, step, fr->fcdata.get(),
+                                           nullptr);
     }
 
     /* Print some stuff... */
@@ -2477,8 +2477,8 @@ void LegacySimulator::do_steep()
 
                 const bool do_dr = do_per_step(steps_accepted, inputrec->nstdisreout);
                 const bool do_or = do_per_step(steps_accepted, inputrec->nstorireout);
-                energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, do_dr, do_or, fplog,
-                                                   count, count, &fr->listedForces->fcdata(), nullptr);
+                energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, do_dr, do_or,
+                                                   fplog, count, count, fr->fcdata.get(), nullptr);
                 fflush(fplog);
             }
         }
index 8f85eb82278b10afc3f9900f4d97e1786cb9c038..12dcb98b7b49845d389be1122d11e95f6a7df52f 100644 (file)
@@ -607,8 +607,7 @@ void gmx::LegacySimulator::do_rerun()
 
             EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
             energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
-                                               do_log ? fplog : nullptr, step, t,
-                                               &fr->listedForces->fcdata(), awh);
+                                               do_log ? fplog : nullptr, step, t, fr->fcdata.get(), awh);
 
             if (do_per_step(step, ir->nstlog))
             {
index 1d1676206b678fbd144b9781167ad2131200905b..8bc7dc71f6f8bd9dbeacb8ed10669b040e5664fc 100644 (file)
@@ -1388,8 +1388,8 @@ int Mdrunner::mdrunner()
                       opt2fn("-tablep", filenames.size(), filenames.data()),
                       opt2fns("-tableb", filenames.size(), filenames.data()), pforce);
         // Dirty hack, for fixing disres and orires should be made mdmodules
-        fr->listedForces->fcdata().disres = disresdata;
-        fr->listedForces->fcdata().orires = oriresdata;
+        fr->fcdata->disres = disresdata;
+        fr->fcdata->orires = oriresdata;
 
         // Save a handle to device stream manager to use elsewhere in the code
         // TODO: Forcerec is not a correct place to store it.
index ea81170069f9f5e5e35fb41e9e2fb51406c86d37..b48bd0910d71383784743b37edf5ce12517bf100 100644 (file)
@@ -56,6 +56,7 @@ struct bonded_threading_t;
 class DeviceContext;
 class DispersionCorrection;
 class ListedForces;
+struct t_fcdata;
 struct t_forcetable;
 struct t_QMMMrec;
 
@@ -298,8 +299,13 @@ struct t_forcerec
     real userreal3 = 0;
     real userreal4 = 0;
 
-    /* The listed forces calculation data */
-    std::unique_ptr<ListedForces> listedForces;
+    /* Data for special listed force calculations */
+    std::unique_ptr<t_fcdata> fcdata;
+
+    // The listed forces calculation data, 1 entry or multiple entries with multiple time stepping
+    std::vector<ListedForces> listedForces;
+
+    static constexpr size_t c_listedForcesFast = 0;
 
     /* TODO: Replace the pointer by an object once we got rid of C */
     gmx::GpuBonded* gpuBonded = nullptr;
index 5dfe43c333c55c3519a5469f808af0cc2c38e88f..8f75ba7436b91c08bdfaf78f02e1f3ed2f17f3ca 100644 (file)
@@ -113,8 +113,6 @@ public:
     bool haveGpuBondedWork = false;
     //! Whether the current nstlist step-range has bonded work to run on he CPU.
     bool haveCpuBondedWork = false;
-    //! Whether the current nstlist step-range has restraints work to run on he CPU.
-    bool haveRestraintsWork = false;
     //! Whether the current nstlist step-range has listed forces work to run on he CPU.
     //  Note: currently this is haveCpuBondedWork | haveRestraintsWork
     bool haveCpuListedForceWork = false;
index dbb04d3fb5746b839d715b0e5a786d155a1d84b6..5604a8dec2b6ff0ec8278d752ae6381771e067b4 100644 (file)
@@ -358,10 +358,9 @@ ModularSimulator::ModularSimulator(std::unique_ptr<LegacySimulatorData> legacySi
 
 void ModularSimulator::checkInputForDisabledFunctionality()
 {
-    isInputCompatible(true, legacySimulatorData_->inputrec,
-                      legacySimulatorData_->mdrunOptions.rerun, *legacySimulatorData_->top_global,
-                      legacySimulatorData_->ms, legacySimulatorData_->replExParams,
-                      &legacySimulatorData_->fr->listedForces->fcdata(),
+    isInputCompatible(true, legacySimulatorData_->inputrec, legacySimulatorData_->mdrunOptions.rerun,
+                      *legacySimulatorData_->top_global, legacySimulatorData_->ms,
+                      legacySimulatorData_->replExParams, legacySimulatorData_->fr->fcdata.get(),
                       opt2bSet("-ei", legacySimulatorData_->nfile, legacySimulatorData_->fnm),
                       legacySimulatorData_->membed != nullptr);
     if (legacySimulatorData_->observablesHistory->edsamHistory)
index 25d5e018c8184e1f8fa05961784f941c5af5f594..44886b83ca1c2280eaa92ee2df97b08ca00774c2 100644 (file)
@@ -408,7 +408,7 @@ ModularSimulatorAlgorithmBuilder::ModularSimulatorAlgorithmBuilder(
             statePropagatorData_.get(), freeEnergyPerturbationData_.get(),
             legacySimulatorData->top_global, legacySimulatorData->inputrec, legacySimulatorData->mdAtoms,
             legacySimulatorData->enerd, legacySimulatorData->ekind, legacySimulatorData->constr,
-            legacySimulatorData->fplog, &legacySimulatorData->fr->listedForces->fcdata(),
+            legacySimulatorData->fplog, legacySimulatorData->fr->fcdata.get(),
             legacySimulatorData->mdModulesNotifier, MASTER(legacySimulatorData->cr),
             legacySimulatorData->observablesHistory, legacySimulatorData->startingBehavior);
 }
index e6b4028814a2a0069bb0ce823c60b66c0bd3ad55..0be9b8bb1ec0b5b0028c83ade77edef8209508f4 100644 (file)
@@ -45,9 +45,9 @@
         str, lstr, (nra), (nrpa), (nrpb), IF_BOND \
     }
 
-#define def_bondedz(str, lstr, nra, nrpa, nrpb)                \
-    {                                                          \
-        str, lstr, (nra), (nrpa), (nrpb), IF_BOND | IF_LIMZERO \
+#define def_pair(str, lstr, nra, nrpa, nrpb)                             \
+    {                                                                    \
+        str, lstr, (nra), (nrpa), (nrpb), IF_BOND | IF_PAIR | IF_LIMZERO \
     }
 
 #define def_bondedt(str, lstr, nra, nrpa, nrpb)                  \
         str, lstr, (nra), (nrpa), (nrpb), IF_BOND | IF_ATYPE \
     }
 
+#define def_dihedral(str, lstr, nra, nrpa, nrpb)                \
+    {                                                           \
+        str, lstr, (nra), (nrpa), (nrpb), IF_BOND | IF_DIHEDRAL \
+    }
+
+#define def_dihedral_tabulated(str, lstr, nra, nrpa, nrpb)                     \
+    {                                                                          \
+        str, lstr, (nra), (nrpa), (nrpb), IF_BOND | IF_DIHEDRAL | IF_TABULATED \
+    }
+
 #define def_bond(str, lstr, nra, nrpa, nrpb)                               \
     {                                                                      \
         str, lstr, (nra), (nrpa), (nrpb), IF_BOND | IF_CHEMBOND | IF_BTYPE \
@@ -118,15 +128,16 @@ const t_interaction_function interaction_function[F_NRE] = {
     def_bonded("CROSS_BOND_BOND", "Bond-Cross", 3, 3, 0),
     def_bonded("CROSS_BOND_ANGLE", "BA-Cross", 3, 4, 0), def_angle("UREY_BRADLEY", "U-B", 3, 4, 4),
     def_angle("QANGLES", "Quartic Angles", 3, 6, 0), def_bondedt("TABANGLES", "Tab. Angles", 3, 2, 2),
-    def_bonded("PDIHS", "Proper Dih.", 4, 3, 3), def_bonded("RBDIHS", "Ryckaert-Bell.", 4, 6, 6),
-    def_bonded("RESTRDIHS", "Restricted Dih.", 4, 2, 2), def_bonded("CBTDIHS", "CBT Dih.", 4, 6, 6),
-    def_bonded("FOURDIHS", "Fourier Dih.", 4, 4, 4), def_bonded("IDIHS", "Improper Dih.", 4, 2, 2),
-    def_bonded("PIDIHS", "Improper Dih.", 4, 3, 3), def_bondedt("TABDIHS", "Tab. Dih.", 4, 2, 2),
-    def_bonded("CMAP", "CMAP Dih.", 5, -1, -1), def_nofc("GB12", "GB 1-2 Pol. (unused)"),
+    def_dihedral("PDIHS", "Proper Dih.", 4, 3, 3), def_dihedral("RBDIHS", "Ryckaert-Bell.", 4, 6, 6),
+    def_dihedral("RESTRDIHS", "Restricted Dih.", 4, 2, 2),
+    def_dihedral("CBTDIHS", "CBT Dih.", 4, 6, 6), def_dihedral("FOURDIHS", "Fourier Dih.", 4, 4, 4),
+    def_dihedral("IDIHS", "Improper Dih.", 4, 2, 2), def_dihedral("PIDIHS", "Improper Dih.", 4, 3, 3),
+    def_dihedral_tabulated("TABDIHS", "Tab. Dih.", 4, 2, 2),
+    def_dihedral("CMAP", "CMAP Dih.", 5, -1, -1), def_nofc("GB12", "GB 1-2 Pol. (unused)"),
     def_nofc("GB13", "GB 1-3 Pol. (unused)"), def_nofc("GB14", "GB 1-4 Pol. (unused)"),
     def_nofc("GBPOL", "GB Polarization (unused)"), def_nofc("NPSOLVATION", "Nonpolar Sol. (unused)"),
-    def_bondedz("LJ14", "LJ-14", 2, 2, 2), def_nofc("COUL14", "Coulomb-14"),
-    def_bondedz("LJC14_Q", "LJC-14 q", 2, 5, 0), def_bondedz("LJC_NB", "LJC Pairs NB", 2, 4, 0),
+    def_pair("LJ14", "LJ-14", 2, 2, 2), def_nofc("COUL14", "Coulomb-14"),
+    def_pair("LJC14_Q", "LJC-14 q", 2, 5, 0), def_pair("LJC_NB", "LJC Pairs NB", 2, 4, 0),
     def_nb("LJ_SR", "LJ (SR)", 2, 2), def_nb("BHAM", "Buck.ham (SR)", 2, 3),
     def_nofc("LJ_LR", "LJ (unused)"), def_nofc("BHAM_LR", "B.ham (unused)"),
     def_nofc("DISPCORR", "Disper. corr."), def_nofc("COUL_SR", "Coulomb (SR)"),
index 584c57f0b76fb6c34276478de5798b2b9363a74e..0ac9c0c3275133e4c6d496701db8a085b20497ff 100644 (file)
@@ -75,8 +75,10 @@ constexpr unsigned int IF_CONSTRAINT = 1 << 2;
 constexpr unsigned int IF_CHEMBOND   = 1 << 3;
 constexpr unsigned int IF_BTYPE      = 1 << 4;
 constexpr unsigned int IF_ATYPE      = 1 << 5;
-constexpr unsigned int IF_TABULATED  = 1 << 6;
-constexpr unsigned int IF_LIMZERO    = 1 << 7;
+constexpr unsigned int IF_DIHEDRAL   = 1 << 6;
+constexpr unsigned int IF_PAIR       = 1 << 7;
+constexpr unsigned int IF_TABULATED  = 1 << 8;
+constexpr unsigned int IF_LIMZERO    = 1 << 9;
 /* These flags tell to some of the routines what can be done with this
  * item in the list.
  * With IF_BOND a bonded interaction will be calculated.