proposed fix for redmine bug #3403

Fixes #3403 by preventing the free energy kernel from skipping pairs
that are beyond the cutoff but also in the exclusion list generated by
couple-intramol=no. Forces calculation of the reciprocal-space Ewald
component for those pairs in the free energy kernel. Also extends the
size of the relevant tables to rcoulomb+tabext to avoid getting nonsense
energies/forces beyond the cutoff. Molecules being decoupled in this way
should be smaller than rcoulomb+tabext.

diff --git a/src/gromacs/ewald/pme_load_balancing.cpp b/src/gromacs/ewald/pme_load_balancing.cpp
index 4bbb44b91d616d6a73b70e798527909eabab949e..73b40cc43aa3471d02e851b59d1b947731378327 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);
+    init_interaction_const_tables(nullptr, ic, 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 8155e5beaed6b4253c7644c1fd56a7afbf4aa90d..a7ea1581761583bb85ab8a6956cc2ad3ae97e13d 100644 (file)
--- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
+++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
@@ -399,14 +399,20 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
             const RealType dz  = iz - x[j3 + 2];
             const RealType rsq = dx * dx + dy * dy + dz * dz;
             RealType       FscalC[NSTATES], FscalV[NSTATES];
+            /* Check if this pair on the exlusions list.*/
+            const bool bPairIncluded = nlist->excl_fep == nullptr || nlist->excl_fep[k];
 
-            if (rsq >= rcutoff_max2)
+            if (rsq >= rcutoff_max2 && bPairIncluded)
             {
                 /* We save significant time by skipping all code below.
                  * Note that with soft-core interactions, the actual cut-off
                  * check might be different. But since the soft-core distance
                  * is always larger than r, checking on r here is safe.
+                 * Exclusions outside the cutoff can not be skipped as
+                 * when using Ewald: the reciprocal-space
+                 * Ewald component still needs to be subtracted.
                  */
+
                 continue;
             }
             npair_within_cutoff++;
@@ -455,7 +461,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
             tj[STATE_A] = ntiA + 2 * typeA[jnr];
             tj[STATE_B] = ntiB + 2 * typeB[jnr];
 
-            if (nlist->excl_fep == nullptr || nlist->excl_fep[k])
+            if (bPairIncluded)
             {
                 c6[STATE_A] = nbfp[tj[STATE_A]];
                 c6[STATE_B] = nbfp[tj[STATE_B]];
@@ -613,7 +619,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
                         FscalC[i] *= rpinvC;
                         FscalV[i] *= rpinvV;
                     }
-                }
+                } // end for (int i = 0; i < NSTATES; i++)
 
                 /* Assemble A and B states */
                 for (int i = 0; i < NSTATES; i++)
@@ -637,7 +643,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
                         dvdl_vdw += Vvdw[i] * DLF[i];
                     }
                 }
-            }
+            } // end if (bPairIncluded)
             else if (icoul == GMX_NBKERNEL_ELEC_REACTIONFIELD)
             {
                 /* For excluded pairs, which are only in this pair list when
@@ -660,7 +666,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
                 }
             }
 
-            if (elecInteractionTypeIsEwald && r < rcoulomb)
+            if (elecInteractionTypeIsEwald && (r < rcoulomb || !bPairIncluded))
             {
                 /* See comment in the preamble. When using Ewald interactions
                  * (unless we use a switch modifier) we subtract the reciprocal-space
@@ -769,7 +775,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
 #pragma omp atomic
                 f[j3 + 2] -= tz;
             }
-        }
+        } // end for (int k = nj0; k < nj1; k++)
 
         /* The atomics below are expensive with many OpenMP threads.
          * Here unperturbed i-particles will usually only have a few
@@ -805,7 +811,7 @@ static void nb_free_energy_kernel(const t_nblist* gmx_restrict nlist,
                 Vv[ggid] += vvtot;
             }
         }
-    }
+    } // end for (int n = 0; n < nri; n++)
 
 #pragma omp atomic
     dvdl[efptCOUL] += dvdl_coul;

diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp
index 2f771943f6b04566c5469b4b764633cfd5080855..c8725816b3814318f76b905cd49b5e480971d2c4 100644 (file)
--- a/src/gromacs/mdlib/forcerec.cpp
+++ b/src/gromacs/mdlib/forcerec.cpp
@@ -712,6 +712,7 @@ 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,
                                EwaldCorrectionTables*     coulombTables,
                                EwaldCorrectionTables*     vdwTables)
 {
@@ -723,7 +724,18 @@ static void init_ewald_f_table(const interaction_const_t& ic,
      */
     const real tableScale = ewald_spline3_table_scale(ic, useCoulombTable, useVdwTable);
 
-    const int tableSize = static_cast<int>(ic.rcoulomb * tableScale) + 2;
+    real tableLen = ic.rcoulomb;
+    if (useCoulombTable && tableExtensionLength > 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.
+         * The alternative is to look through all the exclusions and check if they come from
+         * couple-intramol == no. Meanwhile, always having larger tables should only affect
+         * memory consumption, not speed (barring cache issues).
+         */
+        tableLen = ic.rcoulomb + tableExtensionLength;
+    }
+    const int tableSize = static_cast<int>(tableLen * tableScale) + 2;
 
     if (useCoulombTable)
     {
@@ -737,11 +749,11 @@ static void init_ewald_f_table(const interaction_const_t& ic,
     }
 }
 
-void init_interaction_const_tables(FILE* fp, interaction_const_t* ic)
+void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const real tableExtensionLength)
 {
     if (EEL_PME_EWALD(ic->eeltype))
     {
-        init_ewald_f_table(*ic, ic->coulombEwaldTables.get(), nullptr);
+        init_ewald_f_table(*ic, tableExtensionLength, ic->coulombEwaldTables.get(), nullptr);
         if (fp != nullptr)
         {
             fprintf(fp, "Initialized non-bonded Coulomb Ewald tables, spacing: %.2e size: %zu\n\n",
@@ -1089,7 +1101,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);
+    init_interaction_const_tables(fp, fr->ic, ir->tabext);
 
     const interaction_const_t* ic = fr->ic;
 

diff --git a/src/gromacs/mdlib/forcerec.h b/src/gromacs/mdlib/forcerec.h
index fea69c594f03732d140012eff569eb3393ffb5b2..a8343d62e67e0a3b29f9b37fda3b2b7e9c3950d7 100644 (file)
--- a/src/gromacs/mdlib/forcerec.h
+++ b/src/gromacs/mdlib/forcerec.h
@@ -82,10 +82,11 @@ void forcerec_set_ranges(t_forcerec* fr, int natoms_force, int natoms_force_cons
 /*! \brief Initiate table constants
  *
  * 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] fp                     File for debugging output
+ * \param[in] ic                     Structure holding the table constant
+ * \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);
+void init_interaction_const_tables(FILE* fp, interaction_const_t* ic, const 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 eb29b2092fd8a56409f4c0aab1e8ec273f2aa3b0..df4b03dfa9fa579c6ceb361b0df4acbb00276939 100644 (file)
--- a/src/gromacs/nbnxm/benchmark/bench_setup.cpp
+++ b/src/gromacs/nbnxm/benchmark/bench_setup.cpp
@@ -159,7 +159,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);
+        init_interaction_const_tables(nullptr, &ic, 0);
     }
 
     return ic;