A check for perturbed listed pairs beyond rlist

Excluded atoms pairs should be within the pairlist cut-off, as
otherwise long-range corrections are not applied. This is problematic
for free-energy calculations where intra-molecular interactions are
not coupled and replaced by pair interactions. Now a fatal error
is generated when perturbed pair interactions are beyond rlist.

Also added PME grid correction for LJ excluded pairs beyond rvdw,
which was missing from #3403.

Added missing release note for the fix of #3403.

Fixes #3403
Fixes #3809

diff --git a/api/nblib/gmxsetup.cpp b/api/nblib/gmxsetup.cpp
index 8501ff062a6ce9c22bf7cf02679a88fd60bc443c..f6b320950005b64388051bec54cf6a7a7475a90f 100644 (file)
--- a/api/nblib/gmxsetup.cpp
+++ b/api/nblib/gmxsetup.cpp
@@ -280,7 +280,7 @@ void NbvSetupUtil::setupInteractionConst(const NBKernelOptions& options)
         }
         gmxForceCalculator_->interactionConst_->coulombEwaldTables =
                 std::make_unique<EwaldCorrectionTables>();
-        init_interaction_const_tables(nullptr, gmxForceCalculator_->interactionConst_.get(), 0);
+        init_interaction_const_tables(nullptr, gmxForceCalculator_->interactionConst_.get(), 0, 0);
     }
 }
 

diff --git a/docs/release-notes/2021/major/bugs-fixed.rst b/docs/release-notes/2021/major/bugs-fixed.rst
index a38d2f26ba2dfbfa6e0785159e7ae0082f4b8144..5fec2401a9589292747a9c5a0df8cc60c85492c8 100644 (file)
--- a/docs/release-notes/2021/major/bugs-fixed.rst
+++ b/docs/release-notes/2021/major/bugs-fixed.rst
@@ -25,6 +25,18 @@ H atoms in particular.
 
 :issue:`3469`
 
+Correct excluded perturbed interactions beyond the non-bonded cut-off distance
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+With free-energy calculations without coupling of intermolecular interactions,
+non-bonded pair interactions at distance longer than the cut-off distance can
+be excluded. These interactions would still have PME long-range contributions.
+The contributions are now removed. In addition, mdrun will stop with a fatal
+error when interactions beyond the pair-list cut-off are present.
+
+:issue:`3403`
+:issue:`3808`
+
 Corrected AWH initial histogram size
 """"""""""""""""""""""""""""""""""""
 

diff --git a/src/gromacs/ewald/pme_load_balancing.cpp b/src/gromacs/ewald/pme_load_balancing.cpp
index 73b40cc43aa3471d02e851b59d1b947731378327..542c730156912408abe24c004ffb4ec1b242c431 100644 (file)
--- a/src/gromacs/ewald/pme_load_balancing.cpp
+++ b/src/gromacs/ewald/pme_load_balancing.cpp
@@ -833,7 +833,7 @@ static void pme_load_balance(pme_load_balancing_t*          pme_lb,
     }
 
     /* We always re-initialize the tables whether they are used or not */
-    init_interaction_const_tables(nullptr, ic, ir.tabext);
+    init_interaction_const_tables(nullptr, ic, set->rlistOuter, ir.tabext);
 
     Nbnxm::gpu_pme_loadbal_update_param(nbv, ic);
 

diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
index 53f3dbf9d7f90f59b7942dd73300829ab9dcf70a..6573583c7e50bb06fc1e0ee8594a4f66fc7bb7de 100644 (file)
--- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
+++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
@@ -375,6 +375,10 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
     real* gmx_restrict f      = &(forceWithShiftForces->force()[0][0]);
     real* gmx_restrict fshift = &(forceWithShiftForces->shiftForces()[0][0]);
 
+    const real rlistSquared = gmx::square(fr->rlist);
+
+    int numExcludedPairsBeyondRlist = 0;
+
     for (int n = 0; n < nri; n++)
     {
         int npair_within_cutoff = 0;
@@ -431,6 +435,11 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
             }
             npair_within_cutoff++;
 
+            if (rsq > rlistSquared)
+            {
+                numExcludedPairsBeyondRlist++;
+            }
+
             if (rsq > 0)
             {
                 /* Note that unlike in the nbnxn kernels, we do not need
@@ -722,7 +731,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
                 }
             }
 
-            if (vdwInteractionTypeIsEwald && r < rvdw)
+            if (vdwInteractionTypeIsEwald && (r < rvdw || !bPairIncluded))
             {
                 /* See comment in the preamble. When using LJ-Ewald interactions
                  * (unless we use a switch modifier) we subtract the reciprocal-space
@@ -839,6 +848,19 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
      */
 #pragma omp atomic
     inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri * 12 + nlist->jindex[nri] * 150);
+
+    if (numExcludedPairsBeyondRlist > 0)
+    {
+        gmx_fatal(FARGS,
+                  "There are %d perturbed non-bonded pair interactions beyond the pair-list cutoff "
+                  "of %g nm, which is not supported. This can happen because the system is "
+                  "unstable or because intra-molecular interactions at long distances are "
+                  "excluded. If the "
+                  "latter is the case, you can try to increase nstlist or rlist to avoid this."
+                  "The error is likely triggered by the use of couple-intramol=no "
+                  "and the maximal distance in the decoupled molecule exceeding rlist.",
+                  numExcludedPairsBeyondRlist, fr->rlist);
+    }
 }
 
 typedef void (*KernelFunction)(const t_nblist* gmx_restrict nlist,

diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp
index 0a9f0dc89f6d760a31fcef32f4cdacd9a7982daf..8d7eb271086c8cf129c422602cb7d6263aecdb47 100644 (file)
--- a/src/gromacs/mdlib/forcerec.cpp
+++ b/src/gromacs/mdlib/forcerec.cpp
@@ -724,7 +724,8 @@ static void initVdwEwaldParameters(FILE* fp, const t_inputrec* ir, interaction_c
  * both accuracy requirements, when relevant.
  */
 static void init_ewald_f_table(const interaction_const_t& ic,
-                               const real                 tableExtensionLength,
+                               const real                 rlist,
+                               const real                 tabext,
                                EwaldCorrectionTables*     coulombTables,
                                EwaldCorrectionTables*     vdwTables)
 {
@@ -739,7 +740,7 @@ static void init_ewald_f_table(const interaction_const_t& ic,
     const bool havePerturbedNonbondeds = (ic.softCoreParameters != nullptr);
 
     real tableLen = ic.rcoulomb;
-    if (useCoulombTable && havePerturbedNonbondeds && tableExtensionLength > 0.0)
+    if ((useCoulombTable || useVdwTable) && havePerturbedNonbondeds && rlist + tabext > 0.0)
     {
         /* TODO: Ideally this should also check if couple-intramol == no, but that isn't
          * stored in ir. Grompp puts that info into an opts structure that doesn't make it into the tpr.
@@ -747,7 +748,7 @@ static void init_ewald_f_table(const interaction_const_t& ic,
          * couple-intramol == no. Meanwhile, always having larger tables should only affect
          * memory consumption, not speed (barring cache issues).
          */
-        tableLen = ic.rcoulomb + tableExtensionLength;
+        tableLen = rlist + tabext;
     }
     const int tableSize = static_cast<int>(tableLen * tableScale) + 2;
 
@@ -763,11 +764,11 @@ static void init_ewald_f_table(const interaction_const_t& ic,
     }
 }
 
-void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real tableExtensionLength)
+void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real rlist, const real tableExtensionLength)
 {
     if (EEL_PME_EWALD(ic->eeltype) || EVDW_PME(ic->vdwtype))
     {
-        init_ewald_f_table(*ic, tableExtensionLength, ic->coulombEwaldTables.get(),
+        init_ewald_f_table(*ic, rlist, tableExtensionLength, ic->coulombEwaldTables.get(),
                            ic->vdwEwaldTables.get());
         if (fp != nullptr)
         {
@@ -1081,7 +1082,7 @@ void init_forcerec(FILE*                            fp,
 
     /* fr->ic is used both by verlet and group kernels (to some extent) now */
     init_interaction_const(fp, &fr->ic, ir, mtop, systemHasNetCharge);
-    init_interaction_const_tables(fp, fr->ic, ir->tabext);
+    init_interaction_const_tables(fp, fr->ic, fr->rlist, ir->tabext);
 
     const interaction_const_t* ic = fr->ic;
 

diff --git a/src/gromacs/mdlib/forcerec.h b/src/gromacs/mdlib/forcerec.h
index aed645f1867def785f18a22fd570c5953069ca4b..1a01ff141330ba927ce925138fa8a53c39e7d0fa 100644 (file)
--- a/src/gromacs/mdlib/forcerec.h
+++ b/src/gromacs/mdlib/forcerec.h
@@ -82,9 +82,10 @@ void forcerec_set_ranges(t_forcerec* fr, int natoms_force, int natoms_force_cons
  * Initializes the tables in the interaction constant data structure.
  * \param[in] fp                     File for debugging output
  * \param[in] ic                     Structure holding the table constant
+ * \param[in] rlist                  Length of the neighbour list
  * \param[in] tableExtensionLength   Length by which to extend the tables. Taken from the input record.
  */
-void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, real tableExtensionLength);
+void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, real rlist, real tableExtensionLength);
 
 /*! \brief Initialize forcerec structure.
  *

diff --git a/src/gromacs/nbnxm/benchmark/bench_setup.cpp b/src/gromacs/nbnxm/benchmark/bench_setup.cpp
index 9fffb257ea21f4f4db102c79b1b4d7e1d225fba8..af0b2926995fea593ab5f4921f8a60867986f600 100644 (file)
--- a/src/gromacs/nbnxm/benchmark/bench_setup.cpp
+++ b/src/gromacs/nbnxm/benchmark/bench_setup.cpp
@@ -165,7 +165,7 @@ static interaction_const_t setupInteractionConst(const KernelBenchOptions& optio
         GMX_RELEASE_ASSERT(options.ewaldcoeff_q > 0, "Ewald coefficient should be > 0");
         ic.ewaldcoeff_q       = options.ewaldcoeff_q;
         ic.coulombEwaldTables = std::make_unique<EwaldCorrectionTables>();
-        init_interaction_const_tables(nullptr, &ic, 0);
+        init_interaction_const_tables(nullptr, &ic, 0, 0);
     }
 
     return ic;