Clean up non-PME part of ewald module
authorMark Abraham <mark.j.abraham@gmail.com>
Sun, 12 Oct 2014 16:35:16 +0000 (18:35 +0200)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Fri, 5 Dec 2014 20:19:30 +0000 (21:19 +0100)
ewald/ewald and ewald/ewald-util code had several different kinds of
stuff in it. Separated code that does Ewald long-range and
Ewald-family charge correction from the code specific to the group
scheme.

Moved general-purpose routines that are used in a few other places to
calculate Ewald splitting parameters to the math module.

Minimized header dependencies.

Removed unused things: FLBS FLBSZ

Removed use of typedef for existing opaque struct ewald_tab, per
policy in #1490. Renamed to gmx_ewald_tab_t.

Change-Id: I1394bbd02aa92e6581d011e52c5bee12406a0144

13 files changed:
src/gromacs/ewald/ewald.c
src/gromacs/ewald/ewald.h [new file with mode: 0644]
src/gromacs/ewald/long-range-correction.c [moved from src/gromacs/ewald/ewald-util.c with 90% similarity]
src/gromacs/ewald/long-range-correction.h [moved from src/gromacs/ewald/ewald-util.h with 67% similarity]
src/gromacs/gmxana/gmx_pme_error.cpp
src/gromacs/gmxpreprocess/calc_verletbuf.c
src/gromacs/legacyheaders/types/forcerec.h
src/gromacs/math/CMakeLists.txt
src/gromacs/math/calculate-ewald-splitting-coefficient.cpp [new file with mode: 0644]
src/gromacs/math/calculate-ewald-splitting-coefficient.h [new file with mode: 0644]
src/gromacs/mdlib/force.cpp
src/gromacs/mdlib/forcerec.cpp
src/programs/mdrun/runner.cpp

index dee78cbf74ebc2278491034b0438ff6ba26dd94d..eda483e35d2ab0dac086a909c934d3e31905fa5c 100644 (file)
  */
 #include "gmxpre.h"
 
+#include "ewald.h"
+
 #include <math.h>
 #include <stdio.h>
-#include <stdlib.h>
 
-#include "gromacs/ewald/ewald-util.h"
 #include "gromacs/legacyheaders/macros.h"
-#include "gromacs/legacyheaders/typedefs.h"
 #include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/legacyheaders/types/inputrec.h"
 #include "gromacs/math/gmxcomplex.h"
-#include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/futil.h"
+#include "gromacs/math/vectypes.h"
 #include "gromacs/utility/smalloc.h"
 
-#define TOL 2e-5
-
-struct ewald_tab
+struct gmx_ewald_tab_t
 {
     int        nx, ny, nz, kmax;
     cvec     **eir;
     t_complex *tab_xy, *tab_qxyz;
 };
 
+void init_ewald_tab(struct gmx_ewald_tab_t **et, const t_inputrec *ir, FILE *fp)
+{
+    int n;
 
+    snew(*et, 1);
+    if (fp)
+    {
+        fprintf(fp, "Will do ordinary reciprocal space Ewald sum.\n");
+    }
 
-/* TODO: fix thread-safety */
+    (*et)->nx       = ir->nkx+1;
+    (*et)->ny       = ir->nky+1;
+    (*et)->nz       = ir->nkz+1;
+    (*et)->kmax     = max((*et)->nx, max((*et)->ny, (*et)->nz));
+    (*et)->eir      = NULL;
+    (*et)->tab_xy   = NULL;
+    (*et)->tab_qxyz = NULL;
+}
 
 /* the other routines are in complex.h */
 static t_complex conjmul(t_complex a, t_complex b)
@@ -75,9 +86,6 @@ static t_complex conjmul(t_complex a, t_complex b)
     return c;
 }
 
-
-
-
 static void tabulate_eir(int natom, rvec x[], int kmax, cvec **eir, rvec lll)
 {
     int  i, j, m;
@@ -111,27 +119,6 @@ static void tabulate_eir(int natom, rvec x[], int kmax, cvec **eir, rvec lll)
     }
 }
 
-void init_ewald_tab(ewald_tab_t *et, const t_inputrec *ir, FILE *fp)
-{
-    int n;
-
-    snew(*et, 1);
-    if (fp)
-    {
-        fprintf(fp, "Will do ordinary reciprocal space Ewald sum.\n");
-    }
-
-    (*et)->nx       = ir->nkx+1;
-    (*et)->ny       = ir->nky+1;
-    (*et)->nz       = ir->nkz+1;
-    (*et)->kmax     = max((*et)->nx, max((*et)->ny, (*et)->nz));
-    (*et)->eir      = NULL;
-    (*et)->tab_xy   = NULL;
-    (*et)->tab_qxyz = NULL;
-}
-
-
-
 real do_ewald(t_inputrec *ir,
               rvec x[],        rvec f[],
               real chargeA[],  real chargeB[],
@@ -139,7 +126,7 @@ real do_ewald(t_inputrec *ir,
               t_commrec *cr,   int natoms,
               matrix lrvir,    real ewaldcoeff,
               real lambda,     real *dvdlambda,
-              ewald_tab_t et)
+              struct gmx_ewald_tab_t *et)
 {
     real     factor     = -1.0/(4*ewaldcoeff*ewaldcoeff);
     real     scaleRecip = 4.0*M_PI/(box[XX]*box[YY]*box[ZZ])*ONE_4PI_EPS0/ir->epsilon_r; /* 1/(Vol*e0) */
@@ -299,3 +286,45 @@ real do_ewald(t_inputrec *ir,
 
     return energy;
 }
+
+real ewald_charge_correction(t_commrec *cr, t_forcerec *fr, real lambda,
+                             matrix box,
+                             real *dvdlambda, tensor vir)
+
+{
+    real vol, fac, qs2A, qs2B, vc, enercorr;
+    int  d;
+
+    if (MASTER(cr))
+    {
+        /* Apply charge correction */
+        vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
+
+        fac = M_PI*ONE_4PI_EPS0/(fr->epsilon_r*2.0*vol*vol*sqr(fr->ewaldcoeff_q));
+
+        qs2A = fr->qsum[0]*fr->qsum[0];
+        qs2B = fr->qsum[1]*fr->qsum[1];
+
+        vc = (qs2A*(1 - lambda) + qs2B*lambda)*fac;
+
+        enercorr = -vol*vc;
+
+        *dvdlambda += -vol*(qs2B - qs2A)*fac;
+
+        for (d = 0; d < DIM; d++)
+        {
+            vir[d][d] += vc;
+        }
+
+        if (debug)
+        {
+            fprintf(debug, "Total charge correction: Vcharge=%g\n", enercorr);
+        }
+    }
+    else
+    {
+        enercorr = 0;
+    }
+
+    return enercorr;
+}
diff --git a/src/gromacs/ewald/ewald.h b/src/gromacs/ewald/ewald.h
new file mode 100644 (file)
index 0000000..227f786
--- /dev/null
@@ -0,0 +1,84 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+#ifndef GMX_EWALD_EWALD_H
+#define GMX_EWALD_EWALD_H
+
+#include <stdio.h>
+
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/legacyheaders/types/forcerec.h"
+#include "gromacs/legacyheaders/types/inputrec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/utility/real.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Forward declaration of type for managing Ewald tables */
+struct gmx_ewald_tab_t;
+
+/*! \brief Initialize the tables used in the Ewald long-ranged part */
+void
+init_ewald_tab(struct gmx_ewald_tab_t **et, const t_inputrec *ir,
+               FILE *fp);
+
+/*! \brief Do the long-ranged part of an Ewald calculation */
+real
+do_ewald(t_inputrec *ir,
+         rvec x[],        rvec f[],
+         real chargeA[],  real chargeB[],
+         rvec box,
+         t_commrec *cr,  int natoms,
+         matrix lrvir,   real ewaldcoeff,
+         real lambda,    real *dvdlambda,
+         struct gmx_ewald_tab_t *et);
+
+/*! \brief Calculate the correction to the Ewald sum, due to a net system
+ * charge.
+ *
+ * Should only be called on one thread. */
+real
+ewald_charge_correction(t_commrec *cr, t_forcerec *fr, real lambda, matrix box,
+                        real *dvdlambda, tensor vir);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
similarity index 90%
rename from src/gromacs/ewald/ewald-util.c
rename to src/gromacs/ewald/long-range-correction.c
index 83445f4642278826ffef19c7acb93c8c6aa427df..a71b031dc5970599868571a681c47c3cb267a1e2 100644 (file)
  */
 #include "gmxpre.h"
 
-#include "ewald-util.h"
+#include "long-range-correction.h"
 
 #include <math.h>
-#include <stdio.h>
 
 #include "gromacs/legacyheaders/macros.h"
 #include "gromacs/legacyheaders/names.h"
-#include "gromacs/legacyheaders/txtdump.h"
-#include "gromacs/legacyheaders/typedefs.h"
 #include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/legacyheaders/types/forcerec.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
-#include "gromacs/utility/futil.h"
-#include "gromacs/utility/smalloc.h"
-
-real calc_ewaldcoeff_q(real rc, real dtol)
-{
-    real x = 5, low, high;
-    int  n, i = 0;
-
-
-    do
-    {
-        i++;
-        x *= 2;
-    }
-    while (gmx_erfc(x*rc) > dtol);
-
-    n    = i+60; /* search tolerance is 2^-60 */
-    low  = 0;
-    high = x;
-    for (i = 0; i < n; i++)
-    {
-        x = (low+high)/2;
-        if (gmx_erfc(x*rc) > dtol)
-        {
-            low = x;
-        }
-        else
-        {
-            high = x;
-        }
-    }
-    return x;
-}
-
-static real ewald_function_lj(real x, real rc)
-{
-    real xrc, xrc2, xrc4, factor;
-    xrc  = x*rc;
-    xrc2 = xrc*xrc;
-    xrc4 = xrc2*xrc2;
-#ifdef GMX_DOUBLE
-    factor = exp(-xrc2)*(1 + xrc2 + xrc4/2.0);
-#else
-    factor = expf(-xrc2)*(1 + xrc2 + xrc4/2.0);
-#endif
-
-    return factor;
-}
-
-real calc_ewaldcoeff_lj(real rc, real dtol)
-{
-    real x = 5, low, high;
-    int  n, i = 0;
-
-    do
-    {
-        i++;
-        x *= 2.0;
-    }
-    while (ewald_function_lj(x, rc) > dtol);
-
-    n    = i + 60; /* search tolerance is 2^-60 */
-    low  = 0;
-    high = x;
-    for (i = 0; i < n; ++i)
-    {
-        x = (low + high) / 2.0;
-        if (ewald_function_lj(x, rc) > dtol)
-        {
-            low = x;
-        }
-        else
-        {
-            high = x;
-        }
-    }
-    return x;
-}
 
 void ewald_LRcorrection(int start, int end,
                         t_commrec *cr, int thread, t_forcerec *fr,
@@ -621,45 +541,3 @@ void ewald_LRcorrection(int start, int end,
         }
     }
 }
-
-real ewald_charge_correction(t_commrec *cr, t_forcerec *fr, real lambda,
-                             matrix box,
-                             real *dvdlambda, tensor vir)
-
-{
-    real vol, fac, qs2A, qs2B, vc, enercorr;
-    int  d;
-
-    if (MASTER(cr))
-    {
-        /* Apply charge correction */
-        vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
-
-        fac = M_PI*ONE_4PI_EPS0/(fr->epsilon_r*2.0*vol*vol*sqr(fr->ewaldcoeff_q));
-
-        qs2A = fr->qsum[0]*fr->qsum[0];
-        qs2B = fr->qsum[1]*fr->qsum[1];
-
-        vc = (qs2A*(1 - lambda) + qs2B*lambda)*fac;
-
-        enercorr = -vol*vc;
-
-        *dvdlambda += -vol*(qs2B - qs2A)*fac;
-
-        for (d = 0; d < DIM; d++)
-        {
-            vir[d][d] += vc;
-        }
-
-        if (debug)
-        {
-            fprintf(debug, "Total charge correction: Vcharge=%g\n", enercorr);
-        }
-    }
-    else
-    {
-        enercorr = 0;
-    }
-
-    return enercorr;
-}
similarity index 67%
rename from src/gromacs/ewald/ewald-util.h
rename to src/gromacs/ewald/long-range-correction.h
index e37709136368dfa9e31a24e4db7165e0379f4f3a..3404d75945cb134d6f10e7785fdad36e00499cb1 100644 (file)
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
-#ifndef GMX_EWALD_EWALD_UTIL_H
-#define GMX_EWALD_EWALD_UTIL_H
+#ifndef GMX_EWALD_LONG_RANGE_CORRECTION_H
+#define GMX_EWALD_LONG_RANGE_CORRECTION_H
 
-#include <stdio.h>
-
-#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/legacyheaders/types/forcerec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/real.h"
 
 #ifdef __cplusplus
 extern "C" {
 #endif
 
-/* Ewald related stuff */
-
-void
-init_ewald_tab(ewald_tab_t *et, const t_inputrec *ir,
-               FILE *fp);
-/* initialize the ewald table (as found in the t_forcerec) */
-
-real
-calc_ewaldcoeff_q(real rc, real dtol);
-/* Determines the Ewald parameter, both for Ewald and PME */
-
-real calc_ewaldcoeff_lj(real rc, real dtol);
-/* Determines the Ewald parameters for LJ-PME */
-
-real
-do_ewald(t_inputrec *ir,
-         rvec x[],        rvec f[],
-         real chargeA[],  real chargeB[],
-         rvec box,
-         t_commrec *cr,  int natoms,
-         matrix lrvir,   real ewaldcoeff,
-         real lambda,    real *dvdlambda,
-         ewald_tab_t et);
-/* Do an Ewald calculation for the long range electrostatics. */
-
+/*! \brief Calculate long-range Ewald correction terms.
+ *
+ * For the group cutoff scheme (only), calculates the correction to
+ * the Ewald sums (electrostatic and/or LJ) due to pairs excluded from
+ * the long-ranged part.
+ *
+ * For both cutoff schemes, but only for Coulomb interactions,
+ * calculates correction for surface dipole terms. */
 void
 ewald_LRcorrection(int start, int end,
                    t_commrec *cr, int thread, t_forcerec *fr,
@@ -86,18 +72,6 @@ ewald_LRcorrection(int start, int end,
                    real *Vcorr_q, real *Vcorr_lj,
                    real lambda_q, real lambda_lj,
                    real *dvdlambda_q, real *dvdlambda_lj);
-/* Calculate the Long range correction to the Ewald sums,
- * electrostatic and/or LJ, due to excluded pairs and/or
- * surface dipole terms.
- */
-
-real
-ewald_charge_correction(t_commrec *cr, t_forcerec *fr, real lambda, matrix box,
-                        real *dvdlambda, tensor vir);
-/* Calculate the Long range correction to the Ewald sum,
- * due to a net system charge.
- * Should only be called on one thread.
- */
 
 #ifdef __cplusplus
 }
index 70ae7bf1d4513708d27b8e090a8ee044fee2d018..50e0c55f129e69a4bf6ee1df85be777562d5afe2 100644 (file)
@@ -39,7 +39,6 @@
 #include <algorithm>
 
 #include "gromacs/commandline/pargs.h"
-#include "gromacs/ewald/ewald-util.h"
 #include "gromacs/fileio/tpxio.h"
 #include "gromacs/gmxana/gmx_ana.h"
 #include "gromacs/legacyheaders/calcgrid.h"
@@ -52,6 +51,7 @@
 #include "gromacs/legacyheaders/readinp.h"
 #include "gromacs/legacyheaders/typedefs.h"
 #include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/calculate-ewald-splitting-coefficient.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/random/random.h"
index a2763fe1faae93d9bbb12fe98693d47c9bb26bf8..3ace4d48ab811bd135e36f95fc3ee50c2adb57c6 100644 (file)
@@ -42,9 +42,9 @@
 
 #include <sys/types.h>
 
-#include "gromacs/ewald/ewald-util.h"
 #include "gromacs/legacyheaders/macros.h"
 #include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/math/calculate-ewald-splitting-coefficient.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/nbnxn_consts.h"
index 151c68af6d7cb83fa36fd4f2b48ea68440217829..59cd9f4279d61196597b5d97cccad70202711039 100644 (file)
@@ -182,8 +182,8 @@ typedef struct {
 } cginfo_mb_t;
 
 
-/* ewald table type */
-typedef struct ewald_tab *ewald_tab_t;
+/* Forward declaration of type for managing Ewald tables */
+struct gmx_ewald_tab_t;
 
 typedef struct {
     rvec             *f;
@@ -367,10 +367,10 @@ typedef struct {
     tensor    vir_lj_recip;
 
     /* PME/Ewald stuff */
-    gmx_bool    bEwald;
-    real        ewaldcoeff_q;
-    real        ewaldcoeff_lj;
-    ewald_tab_t ewald_table;
+    gmx_bool                bEwald;
+    real                    ewaldcoeff_q;
+    real                    ewaldcoeff_lj;
+    struct gmx_ewald_tab_t *ewald_table;
 
     /* Virial Stuff */
     rvec *fshift;
index 6b3da2e9a655d76472553a35aac4864af89eb131..94c71d45a3d78d1cfbcd8f1e9b0bb5cdcda1c3bf 100644 (file)
@@ -37,8 +37,9 @@ set(LIBGROMACS_SOURCES ${LIBGROMACS_SOURCES} ${MATH_SOURCES} PARENT_SCOPE)
 
 gmx_install_headers(
     3dtransforms.h
-    gmxcomplex.h
+    calculate-ewald-splitting-coefficient.h
     do_fit.h
+    gmxcomplex.h
     units.h
     utilities.h
     vec.h
diff --git a/src/gromacs/math/calculate-ewald-splitting-coefficient.cpp b/src/gromacs/math/calculate-ewald-splitting-coefficient.cpp
new file mode 100644 (file)
index 0000000..f5b4b76
--- /dev/null
@@ -0,0 +1,117 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "calculate-ewald-splitting-coefficient.h"
+
+#include <cmath>
+
+#include "gromacs/math/utilities.h"
+#include "gromacs/utility/real.h"
+
+real calc_ewaldcoeff_q(real rc, real rtol)
+{
+    real beta = 5, low, high;
+    int  n, i = 0;
+
+    do
+    {
+        i++;
+        beta *= 2;
+    }
+    while (gmx_erfc(beta*rc) > rtol);
+
+    /* Do a binary search with tolerance 2^-60 */
+    n    = i+60;
+    low  = 0;
+    high = beta;
+    for (i = 0; i < n; i++)
+    {
+        beta = (low+high)/2;
+        if (gmx_erfc(beta*rc) > rtol)
+        {
+            low = beta;
+        }
+        else
+        {
+            high = beta;
+        }
+    }
+    return beta;
+}
+
+static real compute_lj_function(real beta, real rc)
+{
+    real xrc, xrc2, xrc4, result;
+    xrc    = beta*rc;
+    xrc2   = xrc*xrc;
+    xrc4   = xrc2*xrc2;
+    result = std::exp(-xrc2)*(1 + xrc2 + xrc4/2.0);
+
+    return result;
+}
+
+real calc_ewaldcoeff_lj(real rc, real rtol)
+{
+    real beta = 5, low, high;
+    int  n, i = 0;
+
+    do
+    {
+        i++;
+        beta *= 2.0;
+    }
+    while (compute_lj_function(beta, rc) > rtol);
+
+    /* Do a binary search with tolerance 2^-60 */
+    n    = i + 60;
+    low  = 0;
+    high = beta;
+    for (i = 0; i < n; ++i)
+    {
+        beta = (low + high) / 2.0;
+        if (compute_lj_function(beta, rc) > rtol)
+        {
+            low = beta;
+        }
+        else
+        {
+            high = beta;
+        }
+    }
+    return beta;
+}
diff --git a/src/gromacs/math/calculate-ewald-splitting-coefficient.h b/src/gromacs/math/calculate-ewald-splitting-coefficient.h
new file mode 100644 (file)
index 0000000..47ebb89
--- /dev/null
@@ -0,0 +1,90 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \file
+ *
+ * \brief Declares functions for computing Ewald splitting coefficients
+ *
+ * These belong in the maths module because they do simple maths and
+ * are used many parts of Gromacs.
+ *
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \inpublicapi
+ */
+#ifndef GMX_MATH_CALCULATE_EWALD_SPLITTING_COEFFICIENT_H
+#define GMX_MATH_CALCULATE_EWALD_SPLITTING_COEFFICIENT_H
+
+#include "gromacs/utility/real.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*! \brief Computes the Ewald splitting coefficient for Coulomb
+ *
+ * Returns a value of beta that satisfies rtol > erfc(beta * rc)
+ * (and is very close to equality). That value is used the same way in
+ * all Coulomb-based Ewald methods.
+ *
+ * \param[in] rc    Cutoff radius
+ * \param[in] rtol  Required maximum value of the short-ranged
+ *                  potential at the cutoff (ie. ewald-rtol)
+ * \return          The value of the splitting coefficient that
+ *                  produces the required dtol at rc.
+ */
+real
+calc_ewaldcoeff_q(real rc, real rtol);
+
+/*! \brief Computes the Ewald splitting coefficient for LJ
+ *
+ * Returns a value of beta that satisfies dtol > erfc(beta * rc) * (1
+ * + beta^2 * rc^2 + 0.5 * beta^4 * rc^4) (and is very close to
+ * equality), which is used in LJ-PME.
+ *
+ * \param[in] rc    Cutoff radius
+ * \param[in] rtol  Required maximum value of the short-ranged
+ *                  potential at the cutoff (ie. ewald-rtol-lj)
+ * \return          The value of the splitting coefficient that
+ *                  produces the required dtol at rc.
+ */
+real
+calc_ewaldcoeff_lj(real rc, real rtol);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
index 80a81d255cf2763f3eacd7044bfa7cb3821fc6ef..16393693189644140e2fab7f46aa16feb3e6efef 100644 (file)
@@ -45,7 +45,8 @@
 #include <string.h>
 
 #include "gromacs/domdec/domdec.h"
-#include "gromacs/ewald/ewald-util.h"
+#include "gromacs/ewald/ewald.h"
+#include "gromacs/ewald/long-range-correction.h"
 #include "gromacs/ewald/pme.h"
 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
 #include "gromacs/legacyheaders/macros.h"
index c25e2ba0b1fa4f269619b1767d9dc6f07647242a..30ee41c8b50a0e9f7364d7681c5901a6052e7a31 100644 (file)
@@ -43,7 +43,7 @@
 #include <string.h>
 
 #include "gromacs/domdec/domdec.h"
-#include "gromacs/ewald/ewald-util.h"
+#include "gromacs/ewald/ewald.h"
 #include "gromacs/legacyheaders/copyrite.h"
 #include "gromacs/legacyheaders/force.h"
 #include "gromacs/legacyheaders/gmx_detect_hardware.h"
@@ -64,6 +64,7 @@
 #include "gromacs/legacyheaders/typedefs.h"
 #include "gromacs/legacyheaders/types/commrec.h"
 #include "gromacs/legacyheaders/types/nbnxn_cuda_types_ext.h"
+#include "gromacs/math/calculate-ewald-splitting-coefficient.h"
 #include "gromacs/math/units.h"
 #include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
index c63b38474e8580db0a21d87da2dd288f3d305e63..a2e01a3e48e8f854cc0d846785edff8c27f1961f 100644 (file)
@@ -48,7 +48,6 @@
 
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/essentialdynamics/edsam.h"
-#include "gromacs/ewald/ewald-util.h"
 #include "gromacs/ewald/pme.h"
 #include "gromacs/fileio/tpxio.h"
 #include "gromacs/gmxpreprocess/calc_verletbuf.h"
@@ -73,6 +72,7 @@
 #include "gromacs/legacyheaders/sighandler.h"
 #include "gromacs/legacyheaders/txtdump.h"
 #include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/math/calculate-ewald-splitting-coefficient.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/nbnxn_consts.h"
 #include "gromacs/mdlib/nbnxn_search.h"