From: Berk Hess <hess@kth.se>
Date: Sun, 29 Jun 2014 09:31:14 +0000 (+0200)
Subject: Take LJ-PME into account for table spacing
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=ca3084144ed75f1ec67fe56312e56cbe45da6080;p=alexxy%2Fgromacs.git

Take LJ-PME into account for table spacing

The table spacing for the grid correction for Ewald interactions
only took Coulomb and not LJ into account. Note that with default
settings LJ tables can use a coarser spacing than Coulomb.

Change-Id: Ib0f831a2177bcf10ad9e03c6764cc155e89abdbb
---

diff --git a/src/gromacs/legacyheaders/tables.h b/src/gromacs/legacyheaders/tables.h
index b16abe9141..6f8695aea9 100644
--- a/src/gromacs/legacyheaders/tables.h
+++ b/src/gromacs/legacyheaders/tables.h
@@ -36,6 +36,7 @@
 #ifndef _tables_h
 #define _tables_h
 #include "types/simple.h"
+#include "types/interaction_const.h"
 
 #ifdef __cplusplus
 extern "C" {
@@ -64,7 +65,7 @@ void table_spline3_fill_ewald_lr(real                                 *table_F,
  * The force can then be interpolated linearly.
  */
 
-real ewald_spline3_table_scale(real ewaldcoeff, real rc);
+real ewald_spline3_table_scale(const interaction_const_t *ic);
 /* Return the scaling for the Ewald quadratic spline tables. */
 
 double v_q_ewald_lr(double beta, double r);
diff --git a/src/gromacs/mdlib/forcerec.c b/src/gromacs/mdlib/forcerec.c
index 9eea665b9a..fa66ca664b 100644
--- a/src/gromacs/mdlib/forcerec.c
+++ b/src/gromacs/mdlib/forcerec.c
@@ -1870,11 +1870,10 @@ static void init_ewald_f_table(interaction_const_t *ic,
 
     if (bUsesSimpleTables)
     {
-        /* With a spacing of 0.0005 we are at the force summation accuracy
-         * for the SSE kernels for "normal" atomistic simulations.
+        /* Get the Ewald table spacing based on Coulomb and/or LJ
+         * Ewald coefficients and rtol.
          */
-        ic->tabq_scale = ewald_spline3_table_scale(ic->ewaldcoeff_q,
-                                                   ic->rcoulomb);
+        ic->tabq_scale = ewald_spline3_table_scale(ic);
 
         maxr           = (rtab > ic->rcoulomb) ? rtab : ic->rcoulomb;
         ic->tabq_size  = (int)(maxr*ic->tabq_scale) + 2;
diff --git a/src/gromacs/mdlib/tables.c b/src/gromacs/mdlib/tables.c
index feb36f60f4..b5b628b3ea 100644
--- a/src/gromacs/mdlib/tables.c
+++ b/src/gromacs/mdlib/tables.c
@@ -310,26 +310,86 @@ void table_spline3_fill_ewald_lr(real                                 *table_f,
     }
 }
 
+/* Returns the spacing for a function using the maximum of
+ * the third derivative, x_scale (unit 1/length)
+ * and function tolerance.
+ */
+static double spline3_table_scale(double third_deriv_max,
+                                  double x_scale,
+                                  double func_tol)
+{
+    double deriv_tol;
+    double sc_deriv, sc_func;
+
+    /* Force tolerance: single precision accuracy */
+    deriv_tol = GMX_FLOAT_EPS;
+    sc_deriv  = sqrt(third_deriv_max/(6*4*deriv_tol*x_scale))*x_scale;
+
+    /* Don't try to be more accurate on energy than the precision */
+    func_tol  = max(func_tol, GMX_REAL_EPS);
+    sc_func   = pow(third_deriv_max/(6*12*sqrt(3)*func_tol), 1.0/3.0)*x_scale;
+
+    return max(sc_deriv, sc_func);
+}
+
 /* The scale (1/spacing) for third order spline interpolation
  * of the Ewald mesh contribution which needs to be subtracted
  * from the non-bonded interactions.
+ * Since there is currently only one spacing for Coulomb and LJ,
+ * the finest spacing is used if both Ewald types are present.
+ *
+ * Note that we could also implement completely separate tables
+ * for Coulomb and LJ Ewald, each with their own spacing.
+ * The current setup with the same spacing can provide slightly
+ * faster kernels with both Coulomb and LJ Ewald, especially
+ * when interleaving both tables (currently not implemented).
  */
-real ewald_spline3_table_scale(real ewaldcoeff, real rc)
+real ewald_spline3_table_scale(const interaction_const_t *ic)
 {
-    double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
-    double ftol, etol;
-    double sc_f, sc_e;
+    real sc;
 
-    /* Force tolerance: single precision accuracy */
-    ftol = GMX_FLOAT_EPS;
-    sc_f = sqrt(erf_x_d3/(6*4*ftol*ewaldcoeff))*ewaldcoeff;
+    sc = 0;
+
+    if (ic->eeltype == eelEWALD || EEL_PME(ic->eeltype))
+    {
+        double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
+        double etol;
+        real   sc_q;
 
-    /* Energy tolerance: 10x more accurate than the cut-off jump */
-    etol = 0.1*gmx_erfc(ewaldcoeff*rc);
-    etol = max(etol, GMX_REAL_EPS);
-    sc_e = pow(erf_x_d3/(6*12*sqrt(3)*etol), 1.0/3.0)*ewaldcoeff;
+        /* Energy tolerance: 0.1 times the cut-off jump */
+        etol  = 0.1*gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+
+        sc_q  = spline3_table_scale(erf_x_d3, ic->ewaldcoeff_q, etol);
+
+        if (debug)
+        {
+            fprintf(debug, "Ewald Coulomb quadratic spline table spacing: %f 1/nm\n", 1/sc_q);
+        }
+
+        sc    = max(sc, sc_q);
+    }
+
+    if (EVDW_PME(ic->vdwtype))
+    {
+        double func_d3 = 0.42888; /* max of (x^-6 (1 - exp(-x^2)(1+x^2+x^4/2)))''' */
+        double xrc2, etol;
+        real   sc_lj;
+
+        /* Energy tolerance: 0.1 times the cut-off jump */
+        xrc2  = sqr(ic->ewaldcoeff_lj*ic->rvdw);
+        etol  = 0.1*exp(-xrc2)*(1 + xrc2 + xrc2*xrc2/2.0);
+
+        sc_lj = spline3_table_scale(func_d3, ic->ewaldcoeff_lj, etol);
+
+        if (debug)
+        {
+            fprintf(debug, "Ewald LJ quadratic spline table spacing: %f 1/nm\n", 1/sc_lj);
+        }
+
+        sc = max(sc, sc_lj);
+    }
 
-    return max(sc_f, sc_e);
+    return sc;
 }
 
 /* Calculate the potential and force for an r value