Avoid numerical overflow with overlapping atoms

The verlet kernels did not allow overlapping atoms, even if they were
not interacting (in contrast to the group kernels). Fixed by clamping
the interaction distance so it can not become smaller than ~6e-4
in single and ~1e-18 in double, and when this number is later
multiplied by zero parameters it will not influence forces. The
clamping should never affect normal interactions; we would previously
crash for distances that were this small.
On Haswell, RF and PME kernels get 3% and 1% slower, respectively.
On CUDA, RF and PME kernels get 1% and 2% faster, respectively.

Fixes #1958.

Change-Id: I83b88f0e9ca34dc151a8b907f334a95a1a4301cc

diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
index f23eefbf378ce20f4640701ea1e8630b19bef7f2..97d1e7e0811f6b9fd95212f3434a59fb9ca6766e 100644 (file)
--- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
+++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.cpp
@@ -390,6 +390,12 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict    nlist,
 
             if (rsq > 0)
             {
+                /* Note that unlike in the nbnxn kernels, we do not need
+                 * to clamp the value of rsq before taking the invsqrt
+                 * to avoid NaN in the LJ calculation, since here we do
+                 * not calculate LJ interactions when C6 and C12 are zero.
+                 */
+
                 rinv         = gmx::invsqrt(rsq);
                 r            = rsq*rinv;
             }

diff --git a/src/gromacs/mdlib/nbnxn_consts.h b/src/gromacs/mdlib/nbnxn_consts.h
index 20e0190b2456eb96d761033236c09eb4b142cb52..f7e00d7c62afa996cdb2769ca5db54e5ffb1ad8e 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_consts.h
+++ b/src/gromacs/mdlib/nbnxn_consts.h
@@ -48,14 +48,18 @@ extern "C" {
 
 #define NBNXN_CPU_CLUSTER_I_SIZE_2LOG  2
 
-/* To avoid NaN when excluded atoms are at zero distance, we add a small
- * number to r^2. NBNXN_AVOID_SING_R2_INC^-3 should fit in real.
- */
-#if !GMX_DOUBLE
-#define NBNXN_AVOID_SING_R2_INC  1.0e-12f
+// Lower limit for square interaction distances in nonbonded kernels.
+// For smaller values we will overflow when calculating r^-1 or r^-12, but
+// to keep it simple we always apply the limit from the tougher r^-12 condition.
+#if GMX_DOUBLE
+// Some double precision SIMD architectures use single precision in the first
+// step, so although the double precision criterion would allow smaller rsq,
+// we need to stay in single precision with some margin for the N-R iterations.
+#define NBNXN_MIN_RSQ         1.0e-36
 #else
-/* The double prec. x86 SIMD kernels use a single prec. invsqrt, so > 1e-38 */
-#define NBNXN_AVOID_SING_R2_INC  1.0e-36
+// The worst intermediate value we might evaluate is r^-12, which
+// means we should ensure r^2 stays above pow(GMX_FLOAT_MAX,-1.0/6.0)*1.01 (some margin)
+#define NBNXN_MIN_RSQ         3.82e-07f  // r > 6.2e-4
 #endif
 
 

diff --git a/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel.cuh b/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel.cuh
index 737c93f5ab9e56a1c91f37dd6ee87f42b21c787d..c63eda7f43bbb1f761d06018c4570a26482ea6a9 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel.cuh
+++ b/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel.cuh
@@ -441,8 +441,8 @@ __global__ void NB_KERNEL_FUNC_NAME(nbnxn_kernel, _F_cuda)
 #endif                          /* LJ_COMB_GEOM */
 #endif                          /* LJ_COMB */
 
-                                /* avoid NaN for excluded pairs at r=0 */
-                                r2      += (1.0f - int_bit) * NBNXN_AVOID_SING_R2_INC;
+                                // Ensure distance do not become so small that r^-12 overflows
+                                r2      = max(r2, NBNXN_MIN_RSQ);
 
                                 inv_r   = rsqrt(r2);
                                 inv_r2  = inv_r * inv_r;

diff --git a/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_fermi.cuh b/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_fermi.cuh
index 9321e16cda33a688c07abe8c5b8060379b8cd19d..0623603903a27cc36f9b88916859241d6c2937ce 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_fermi.cuh
+++ b/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_fermi.cuh
@@ -388,8 +388,8 @@ __global__ void NB_KERNEL_FUNC_NAME(nbnxn_kernel, _F_cuda)
 #endif                          /* LJ_COMB_GEOM */
 #endif                          /* LJ_COMB */
 
-                                /* avoid NaN for excluded pairs at r=0 */
-                                r2      += (1.0f - int_bit) * NBNXN_AVOID_SING_R2_INC;
+                                // Ensure distance do not become so small that r^-12 overflows
+                                r2      = max(r2, NBNXN_MIN_RSQ);
 
                                 inv_r   = rsqrt(r2);
                                 inv_r2  = inv_r * inv_r;

diff --git a/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_gpu_ref.cpp b/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_gpu_ref.cpp
index 04c28a52f73cbb2c8dc2629449ce35173fea57fc..92946aeea5a130fd80765f7af7cacc5632703ab3 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_gpu_ref.cpp
+++ b/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_gpu_ref.cpp
@@ -40,6 +40,8 @@
 
 #include <cmath>
 
+#include <algorithm>
+
 #include "gromacs/math/functions.h"
 #include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
@@ -256,8 +258,8 @@ nbnxn_kernel_gpu_ref(const nbnxn_pairlist_t     *nbl,
                                     npair++;
                                 }
 
-                                /* avoid NaN for excluded pairs at r=0 */
-                                rsq             += (1.0 - int_bit)*NBNXN_AVOID_SING_R2_INC;
+                                // Ensure distance do not become so small that r^-12 overflows
+                                rsq              = std::max(rsq, NBNXN_MIN_RSQ);
 
                                 rinv             = gmx::invsqrt(rsq);
                                 rinvsq           = rinv*rinv;

diff --git a/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref.cpp b/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref.cpp
index 42ba187fea92707abc1b435389d9104e5000fed4..fe8d744083284b7f28c432e5c4bc832fe2c97b60 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref.cpp
+++ b/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref.cpp
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -42,6 +42,8 @@
 
 #include <cmath>
 
+#include <algorithm>
+
 #include "gromacs/math/functions.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/force.h"

diff --git a/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h b/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h
index a5c9ca103fe6f8afd98c8a3886330a3c0919510f..846a281037a762a79ffd16910a90d619da019aba 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h
+++ b/src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref_inner.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -128,13 +128,9 @@
             skipmask = (rsq >= rcut2) ? 0 : skipmask;
             /* 9 flops for r^2 + cut-off check */
 
-#ifdef CHECK_EXCLS
-            /* Excluded atoms are allowed to be on top of each other.
-             * To avoid overflow of rinv, rinvsq and rinvsix
-             * we add a small number to rsq for excluded pairs only.
-             */
-            rsq += (1 - interact)*NBNXN_AVOID_SING_R2_INC;
-#endif
+            // Ensure the distances do not fall below the limit where r^-12 overflows.
+            // This should never happen for normal interactions.
+            rsq = std::max(rsq, NBNXN_MIN_RSQ);
 
 #ifdef COUNT_PAIRS
             npair++;

diff --git a/src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h b/src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h
index a7838d9fa7027bde5bcb83ebf106c982ec8f86e2..6d3c057d4088750375b883b6f479260cecf1c7d9 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h
+++ b/src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_inner.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -295,11 +295,10 @@
     }
 #endif
 
-#ifdef CHECK_EXCLS
-    /* For excluded pairs add a small number to avoid r^-6 = NaN */
-    rsq_S0      = rsq_S0 + selectByNotMask(avoid_sing_S, interact_S0);
-    rsq_S2      = rsq_S2 + selectByNotMask(avoid_sing_S, interact_S2);
-#endif
+    // Ensure the distances do not fall below the limit where r^-12 overflows.
+    // This should never happen for normal interactions.
+    rsq_S0      = max(rsq_S0, minRsq_S);
+    rsq_S2      = max(rsq_S2, minRsq_S);
 
     /* Calculate 1/r */
     rinv_S0     = invsqrt(rsq_S0);

diff --git a/src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_outer.h b/src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_outer.h
index 2f4c91c871e1043922e338cdd0258f7466def678..c346e7a9de3561285352b732abb36eba67a0834b 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_outer.h
+++ b/src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_outer.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -147,7 +147,7 @@
 #endif
 #endif /* LJ_COMB_LB */
 
-    SimdReal  avoid_sing_S;
+    SimdReal  minRsq_S;
     SimdReal  rc2_S;
 #ifdef VDW_CUTOFF_CHECK
     SimdReal  rcvdw2_S;
@@ -305,7 +305,7 @@
     rcvdw2_S = SimdReal(ic->rvdw*ic->rvdw);
 #endif
 
-    avoid_sing_S = SimdReal(NBNXN_AVOID_SING_R2_INC);
+    minRsq_S            = SimdReal(NBNXN_MIN_RSQ);
 
     q                   = nbat->q;
     facel               = ic->epsfac;

diff --git a/src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_inner.h b/src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_inner.h
index 231419f2e0ab4d7072b9cdf3a36653f7f1038923..7d954abb841f4b671f8448977ec0841da979fbf9 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_inner.h
+++ b/src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_inner.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -379,13 +379,12 @@
     }
 #endif
 
-#ifdef CHECK_EXCLS
-    /* For excluded pairs add a small number to avoid r^-6 = NaN */
-    rsq_S0      = rsq_S0 + selectByNotMask(avoid_sing_S, interact_S0);
-    rsq_S1      = rsq_S1 + selectByNotMask(avoid_sing_S, interact_S1);
-    rsq_S2      = rsq_S2 + selectByNotMask(avoid_sing_S, interact_S2);
-    rsq_S3      = rsq_S3 + selectByNotMask(avoid_sing_S, interact_S3);
-#endif
+    // Ensure the distances do not fall below the limit where r^-12 overflows.
+    // This should never happen for normal interactions.
+    rsq_S0      = max(rsq_S0, minRsq_S);
+    rsq_S1      = max(rsq_S1, minRsq_S);
+    rsq_S2      = max(rsq_S2, minRsq_S);
+    rsq_S3      = max(rsq_S3, minRsq_S);
 
     /* Calculate 1/r */
 #if !GMX_DOUBLE

diff --git a/src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_outer.h b/src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_outer.h
index 3458e92366abc01f80a6ac7c76a4b58a8e68c081..86bf7e8897f1f3f7b65e8908cf9776baacf74ecf 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_outer.h
+++ b/src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_outer.h
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -154,7 +154,7 @@
 #endif
 #endif /* LJ_COMB_LB */
 
-    SimdReal  avoid_sing_S;
+    SimdReal  minRsq_S;
     SimdReal  rc2_S;
 #ifdef VDW_CUTOFF_CHECK
     SimdReal  rcvdw2_S;
@@ -326,7 +326,7 @@
     rcvdw2_S =  SimdReal(ic->rvdw*ic->rvdw);
 #endif
 
-    avoid_sing_S = SimdReal(NBNXN_AVOID_SING_R2_INC);
+    minRsq_S            = SimdReal(NBNXN_MIN_RSQ);
 
     q                   = nbat->q;
     facel               = ic->epsfac;

diff --git a/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_jit_support.cpp b/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_jit_support.cpp
index f1a3b79b0e7d8d16039baca021b7e610db61de47..c6a63e79c15621cff1feb708cd895fee87c1e8f6 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_jit_support.cpp
+++ b/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_jit_support.cpp
@@ -197,13 +197,13 @@ nbnxn_gpu_compile_kernels(gmx_nbnxn_ocl_t *nb)
          * in the JIT compilation that happens at runtime.
          */
         extraDefines += gmx::formatString(
-                    " -DCENTRAL=%d -DNBNXN_GPU_NCLUSTER_PER_SUPERCLUSTER=%d -DNBNXN_GPU_CLUSTER_SIZE=%d -DNBNXN_GPU_JGROUP_SIZE=%d -DNBNXN_AVOID_SING_R2_INC=%s %s",
+                    " -DCENTRAL=%d -DNBNXN_GPU_NCLUSTER_PER_SUPERCLUSTER=%d -DNBNXN_GPU_CLUSTER_SIZE=%d -DNBNXN_GPU_JGROUP_SIZE=%d -DNBNXN_MIN_RSQ=%s %s",
                     CENTRAL,                                 /* Defined in ishift.h */
                     c_nbnxnGpuNumClusterPerSupercluster,     /* Defined in nbnxn_pairlist.h */
                     c_nbnxnGpuClusterSize,                   /* Defined in nbnxn_pairlist.h */
                     c_nbnxnGpuJgroupSize,                    /* Defined in nbnxn_pairlist.h */
-                    STRINGIFY_MACRO(NBNXN_AVOID_SING_R2_INC) /* Defined in nbnxn_consts.h */
-                                                             /* NBNXN_AVOID_SING_R2_INC passed as string to avoid
+                    STRINGIFY_MACRO(NBNXN_MIN_RSQ)           /* Defined in nbnxn_consts.h */
+                                                             /* NBNXN_MIN_RSQ passed as string to avoid
                                                                 floating point representation problems with sprintf */
                     , (nb->bPrefetchLjParam) ? "-DIATYPE_SHMEM" : ""
                     );

diff --git a/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_amd.clh b/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_amd.clh
index 056b7944afdeb7e5118298fe7b394f019a5462f6..d5a7fdf6fa5d0c5e1678d8a7b2c113e9780b8848 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_amd.clh
+++ b/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_amd.clh
@@ -440,8 +440,8 @@ __global float *restrict fshift,            /* stores float3 values */     /* OU
 #endif                          /* LJ_COMB_GEOM */
 #endif                          /* LJ_COMB */
 
-                                /* avoid NaN for excluded pairs at r=0 */
-                                r2      += (1.0f - int_bit) * NBNXN_AVOID_SING_R2_INC;
+                                // Ensure distance do not become so small that r^-12 overflows
+                                r2      = max(r2,NBNXN_MIN_RSQ);
 
                                 inv_r   = rsqrt(r2);
                                 inv_r2  = inv_r * inv_r;

diff --git a/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_nowarp.clh b/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_nowarp.clh
index 06be7b15f3ea606b44f9714115ac1e4b7933c2dd..d544a69f61eac6da17fb8ed00eba553b405c67d7 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_nowarp.clh
+++ b/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_nowarp.clh
@@ -443,8 +443,8 @@ __global float *restrict fshift,            /* stores float3 values */     /* OU
 #endif                          /* LJ_COMB_GEOM */
 #endif                          /* LJ_COMB */
 
-                                /* avoid NaN for excluded pairs at r=0 */
-                                r2      += (1.0f - int_bit) * NBNXN_AVOID_SING_R2_INC;
+                                // Ensure distance do not become so small that r^-12 overflows
+                                r2      = max(r2,NBNXN_MIN_RSQ);
 
                                 inv_r   = rsqrt(r2);
                                 inv_r2  = inv_r * inv_r;

diff --git a/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_nvidia.clh b/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_nvidia.clh
index 6e93c21e190eb9a4b9c96908beb10896908fd5f4..2fbcaee413a6ba0c180dc2bdac6594e3fc23dda3 100644 (file)
--- a/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_nvidia.clh
+++ b/src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_kernel_nvidia.clh
@@ -433,8 +433,8 @@ __global float *restrict fshift,            /* stores float3 values */     /* OU
 #endif                          /* LJ_COMB_GEOM */
 #endif                          /* LJ_COMB */
 
-                                /* avoid NaN for excluded pairs at r=0 */
-                                r2      += (1.0f - int_bit) * NBNXN_AVOID_SING_R2_INC;
+                                // Ensure distance do not become so small that r^-12 overflows
+                                r2      = max(r2,NBNXN_MIN_RSQ);
 
                                 inv_r   = rsqrt(r2);
                                 inv_r2  = inv_r * inv_r;