Move non-bonded initialization out of the forcerec
authorArtem Zhmurov <zhmurov@gmail.com>
Tue, 4 Feb 2020 10:51:45 +0000 (11:51 +0100)
committerArtem Zhmurov <zhmurov@gmail.com>
Wed, 5 Feb 2020 20:37:09 +0000 (21:37 +0100)
Forcerec should not be responsible for initializing modules.

Change-Id: I666075fe03441c815b191f1eb3809e294c62d541

src/gromacs/mdlib/forcerec.cpp
src/gromacs/mdlib/forcerec.h
src/gromacs/mdrun/runner.cpp
src/gromacs/nbnxm/nbnxm.h
src/gromacs/nbnxm/nbnxm_setup.cpp
src/gromacs/topology/mtop_util.cpp
src/gromacs/topology/mtop_util.h

index 1aa3402e78804208734ebb2168d01985f8d24304..cadf5e6bc1de809f7216018e9eb8afb795e9e85e 100644 (file)
@@ -188,7 +188,7 @@ enum
     acSETTLE
 };
 
-static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t* mtop, const t_forcerec* fr, gmx_bool* bFEP_NonBonded)
+static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t* mtop, const t_forcerec* fr)
 {
     gmx_bool* type_VDW;
     int*      a_con;
@@ -204,8 +204,6 @@ static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t* mtop, const t_f
         }
     }
 
-    *bFEP_NonBonded = FALSE;
-
     std::vector<cginfo_mb_t> cginfoPerMolblock;
     int                      a_offset = 0;
     for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
@@ -315,7 +313,6 @@ static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t* mtop, const t_f
                 if (fr->efep != efepNO && PERTURBED(atom))
                 {
                     SET_CGINFO_FEP(atomInfo);
-                    *bFEP_NonBonded = TRUE;
                 }
             }
         }
@@ -936,18 +933,9 @@ void init_forcerec(FILE*                            fp,
                    const char*                      tabfn,
                    const char*                      tabpfn,
                    gmx::ArrayRef<const std::string> tabbfnm,
-                   const gmx_hw_info_t&             hardwareInfo,
-                   const gmx_device_info_t*         deviceInfo,
-                   const bool                       useGpuForBonded,
                    const bool                       pmeOnlyRankUsesGpu,
-                   real                             print_force,
-                   gmx_wallcycle*                   wcycle)
+                   real                             print_force)
 {
-    real     rtab;
-    char*    env;
-    double   dbl;
-    gmx_bool bFEP_NonBonded;
-
     /* By default we turn SIMD kernels on, but it might be turned off further down... */
     fr->use_simd_kernels = TRUE;
 
@@ -1011,10 +999,10 @@ void init_forcerec(FILE*                            fp,
     fr->sc_r_power    = ir->fepvals->sc_r_power;
     fr->sc_sigma6_def = gmx::power6(ir->fepvals->sc_sigma);
 
-    env = getenv("GMX_SCSIGMA_MIN");
+    char* env = getenv("GMX_SCSIGMA_MIN");
     if (env != nullptr)
     {
-        dbl = 0;
+        double dbl = 0;
         sscanf(env, "%20lf", &dbl);
         fr->sc_sigma6_min = gmx::power6(dbl);
         if (fp)
@@ -1318,7 +1306,7 @@ void init_forcerec(FILE*                            fp,
      * in that case grompp should already have checked that we do not need
      * normal tables and we only generate tables for 1-4 interactions.
      */
-    rtab = ir->rlist + ir->tabext;
+    real rtab = ir->rlist + ir->tabext;
 
     /* We want to use unmodified tables for 1-4 coulombic
      * interactions, so we must in general have an extra set of
@@ -1386,7 +1374,7 @@ void init_forcerec(FILE*                            fp,
     }
 
     /* Set all the static charge group info */
-    fr->cginfo_mb = init_cginfo_mb(mtop, fr, &bFEP_NonBonded);
+    fr->cginfo_mb = init_cginfo_mb(mtop, fr);
     if (!DOMAINDECOMP(cr))
     {
         fr->cginfo = cginfo_expand(mtop->molblock.size(), fr->cginfo_mb);
@@ -1406,38 +1394,6 @@ void init_forcerec(FILE*                            fp,
     fr->nthread_ewc = gmx_omp_nthreads_get(emntBonded);
     snew(fr->ewc_t, fr->nthread_ewc);
 
-    if (fr->cutoff_scheme == ecutsVERLET)
-    {
-        // We checked the cut-offs in grompp, but double-check here.
-        // We have PME+LJcutoff kernels for rcoulomb>rvdw.
-        if (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == eelCUT)
-        {
-            GMX_RELEASE_ASSERT(ir->rcoulomb >= ir->rvdw,
-                               "With Verlet lists and PME we should have rcoulomb>=rvdw");
-        }
-        else
-        {
-            GMX_RELEASE_ASSERT(
-                    ir->rcoulomb == ir->rvdw,
-                    "With Verlet lists and no PME rcoulomb and rvdw should be identical");
-        }
-
-        fr->nbv = Nbnxm::init_nb_verlet(mdlog, bFEP_NonBonded, ir, fr, cr, hardwareInfo, deviceInfo,
-                                        mtop, box, wcycle);
-
-        if (useGpuForBonded)
-        {
-            auto stream = havePPDomainDecomposition(cr)
-                                  ? Nbnxm::gpu_get_command_stream(
-                                            fr->nbv->gpu_nbv, gmx::InteractionLocality::NonLocal)
-                                  : Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv,
-                                                                  gmx::InteractionLocality::Local);
-            // TODO the heap allocation is only needed while
-            // t_forcerec lacks a constructor.
-            fr->gpuBonded = new gmx::GpuBonded(mtop->ffparams, stream, wcycle);
-        }
-    }
-
     if (ir->eDispCorr != edispcNO)
     {
         fr->dispersionCorrection = std::make_unique<DispersionCorrection>(
index e5504f48383ea626780de38d891512832f14dce7..53fc326e30ccfa737319aa8c3e3ca14d42655752 100644 (file)
@@ -100,12 +100,8 @@ void init_interaction_const_tables(FILE* fp, interaction_const_t* ic);
  * \param[in]  tabfn              Table potential file for non-bonded interactions
  * \param[in]  tabpfn             Table potential file for pair interactions
  * \param[in]  tabbfnm            Table potential files for bonded interactions
- * \param[in]  hardwareInfo       Information about hardware
- * \param[in]  deviceInfo         Info about GPU device to use for short-ranged work
- * \param[in]  useGpuForBonded    Whether bonded interactions will run on a GPU
  * \param[in]  pmeOnlyRankUsesGpu Whether there is a PME task on a GPU on a PME-only rank
  * \param[in]  print_force        Print forces for atoms with force >= print_force
- * \param[out] wcycle             Pointer to cycle counter object
  */
 void init_forcerec(FILE*                            fplog,
                    const gmx::MDLogger&             mdlog,
@@ -118,12 +114,8 @@ void init_forcerec(FILE*                            fplog,
                    const char*                      tabfn,
                    const char*                      tabpfn,
                    gmx::ArrayRef<const std::string> tabbfnm,
-                   const gmx_hw_info_t&             hardwareInfo,
-                   const gmx_device_info_t*         deviceInfo,
-                   bool                             useGpuForBonded,
                    bool                             pmeOnlyRankUsesGpu,
-                   real                             print_force,
-                   gmx_wallcycle*                   wcycle);
+                   real                             print_force);
 
 /*! \brief Divide exclusions over threads
  *
index 1c9b18e96c5412caee81d5f4f97d91a98906d69f..364d0bd9c29bc1f086b60d9307b274f6677b270b 100644 (file)
@@ -411,6 +411,18 @@ static void prepare_verlet_scheme(FILE*               fplog,
                                   bool                makeGpuPairList,
                                   const gmx::CpuInfo& cpuinfo)
 {
+    // We checked the cut-offs in grompp, but double-check here.
+    // We have PME+LJcutoff kernels for rcoulomb>rvdw.
+    if (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == eelCUT)
+    {
+        GMX_RELEASE_ASSERT(ir->rcoulomb >= ir->rvdw,
+                           "With Verlet lists and PME we should have rcoulomb>=rvdw");
+    }
+    else
+    {
+        GMX_RELEASE_ASSERT(ir->rcoulomb == ir->rvdw,
+                           "With Verlet lists and no PME rcoulomb and rvdw should be identical");
+    }
     /* For NVE simulations, we will retain the initial list buffer */
     if (EI_DYNAMICS(ir->eI) && ir->verletbuf_tol > 0 && !(EI_MD(ir->eI) && ir->etc == etcNO))
     {
@@ -1315,6 +1327,7 @@ int Mdrunner::mdrunner()
     const bool                   thisRankHasPmeGpuTask = gpuTaskAssignments.thisRankHasPmeGpuTask();
     std::unique_ptr<MDAtoms>     mdAtoms;
     std::unique_ptr<gmx_vsite_t> vsite;
+    std::unique_ptr<GpuBonded>   gpuBonded;
 
     t_nrnb nrnb;
     if (thisRankHasDuty(cr, DUTY_PP))
@@ -1329,9 +1342,21 @@ int Mdrunner::mdrunner()
         init_forcerec(fplog, mdlog, fr, fcd, inputrec, &mtop, cr, box,
                       opt2fn("-table", filenames.size(), filenames.data()),
                       opt2fn("-tablep", filenames.size(), filenames.data()),
-                      opt2fns("-tableb", filenames.size(), filenames.data()), *hwinfo,
-                      nonbondedDeviceInfo, useGpuForBonded,
-                      pmeRunMode == PmeRunMode::GPU && !thisRankHasDuty(cr, DUTY_PME), pforce, wcycle);
+                      opt2fns("-tableb", filenames.size(), filenames.data()),
+                      pmeRunMode == PmeRunMode::GPU && !thisRankHasDuty(cr, DUTY_PME), pforce);
+
+        fr->nbv = Nbnxm::init_nb_verlet(mdlog, inputrec, fr, cr, *hwinfo, nonbondedDeviceInfo,
+                                        &mtop, box, wcycle);
+        if (useGpuForBonded)
+        {
+            auto stream = havePPDomainDecomposition(cr)
+                                  ? Nbnxm::gpu_get_command_stream(
+                                            fr->nbv->gpu_nbv, gmx::InteractionLocality::NonLocal)
+                                  : Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv,
+                                                                  gmx::InteractionLocality::Local);
+            gpuBonded     = std::make_unique<GpuBonded>(mtop.ffparams, stream, wcycle);
+            fr->gpuBonded = gpuBonded.get();
+        }
 
         // TODO Move this to happen during domain decomposition setup,
         // once stream and event handling works well with that.
@@ -1641,6 +1666,7 @@ int Mdrunner::mdrunner()
     mdAtoms.reset(nullptr);
     globalState.reset(nullptr);
     mdModules_.reset(nullptr); // destruct force providers here as they might also use the GPU
+    gpuBonded.reset(nullptr);
     /* Free pinned buffers in *fr */
     delete fr;
     fr = nullptr;
index 3cd2a38ee4fdd01f96822218697403c08037072f..65b3c1a1ce0c312e410dc340e3add4664e41ea11 100644 (file)
@@ -404,7 +404,6 @@ namespace Nbnxm
 
 /*! \brief Creates an Nbnxm object */
 std::unique_ptr<nonbonded_verlet_t> init_nb_verlet(const gmx::MDLogger&     mdlog,
-                                                   gmx_bool                 bFEP_NonBonded,
                                                    const t_inputrec*        ir,
                                                    const t_forcerec*        fr,
                                                    const t_commrec*         cr,
index 35aea4ae3dcdfbd6ea5349e10ee58b0b34b6455e..c45d54decf96dde6ba8640999fa7d41b4e8f627c 100644 (file)
@@ -54,6 +54,7 @@
 #include "gromacs/nbnxm/nbnxm.h"
 #include "gromacs/nbnxm/pairlist_tuning.h"
 #include "gromacs/simd/simd.h"
+#include "gromacs/topology/mtop_util.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/logger.h"
 
@@ -357,7 +358,6 @@ static int getMinimumIlistCountForGpuBalancing(NbnxmGpu* nbnxmGpu)
 }
 
 std::unique_ptr<nonbonded_verlet_t> init_nb_verlet(const gmx::MDLogger&     mdlog,
-                                                   gmx_bool                 bFEP_NonBonded,
                                                    const t_inputrec*        ir,
                                                    const t_forcerec*        fr,
                                                    const t_commrec*         cr,
@@ -392,6 +392,7 @@ std::unique_ptr<nonbonded_verlet_t> init_nb_verlet(const gmx::MDLogger&     mdlo
 
     const bool haveMultipleDomains = (DOMAINDECOMP(cr) && cr->dd->nnodes > 1);
 
+    bool           bFEP_NonBonded = (fr->efep != efepNO) && haveFepPerturbedNBInteractions(mtop);
     PairlistParams pairlistParams(kernelSetup.kernelType, bFEP_NonBonded, ir->rlist,
                                   havePPDomainDecomposition(cr));
 
index 6f6a9d64b03c7a882eb919b9d97e80b30b65116b..e1f7ce822d809b58374ffff5e706e31291d13acd 100644 (file)
@@ -1139,3 +1139,24 @@ void convertAtomsToMtop(t_symtab* symtab, char** name, t_atoms* atoms, gmx_mtop_
 
     gmx_mtop_finalize(mtop);
 }
+
+bool haveFepPerturbedNBInteractions(const gmx_mtop_t* mtop)
+{
+    for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
+    {
+        const gmx_molblock_t& molb = mtop->molblock[mb];
+        const gmx_moltype_t&  molt = mtop->moltype[molb.type];
+        for (int m = 0; m < molb.nmol; m++)
+        {
+            for (int a = 0; a < molt.atoms.nr; a++)
+            {
+                const t_atom& atom = molt.atoms.atom[a];
+                if (PERTURBED(atom))
+                {
+                    return true;
+                }
+            }
+        }
+    }
+    return false;
+}
index 779d77f1aade98271a0574d07bb3efea06160483..2843f28c14bbd2c7cf5a53aa35ae2b51d0d497c3 100644 (file)
@@ -307,4 +307,11 @@ std::vector<int> get_atom_index(const gmx_mtop_t* mtop);
  */
 void convertAtomsToMtop(t_symtab* symtab, char** name, t_atoms* atoms, gmx_mtop_t* mtop);
 
+/*! \brief Checks if the non-bonded FEP should be performed in this run.
+ *
+ * \param[in]  mtop  Molecular topology.
+ * \returns Whether FEP non-bonded is requested.
+ */
+bool haveFepPerturbedNBInteractions(const gmx_mtop_t* mtop);
+
 #endif