-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * 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.
*
- * If you want to redistribute modifications, 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 www.gromacs.org.
+ * 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.
*
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * 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.
*
- * For more info, check our website at http://www.gromacs.org
+ * 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.
*
- * And Hey:
- * GROwing Monsters And Cloning Shrimps
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "gromacs/legacyheaders/tables.h"
#include <math.h>
-#include "maths.h"
-#include "typedefs.h"
-#include "names.h"
-#include "smalloc.h"
-#include "gmx_fatal.h"
-#include "futil.h"
-#include "xvgr.h"
-#include "vec.h"
-#include "main.h"
-#include "network.h"
-#include "physics.h"
-#include "force.h"
-#include "gmxfio.h"
-#include "macros.h"
-#include "tables.h"
+
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
/* All the possible (implemented) table functions */
enum {
etabEwaldSwitch,
etabEwaldUser,
etabEwaldUserSwitch,
+ etabLJ6Ewald,
etabLJ6Switch,
etabLJ12Switch,
etabCOULSwitch,
{ "Ewald-Switch", TRUE },
{ "Ewald-User", TRUE },
{ "Ewald-User-Switch", TRUE },
+ { "LJ6Ewald", FALSE },
{ "LJ6Switch", FALSE },
{ "LJ12Switch", FALSE },
{ "COULSwitch", TRUE },
{ "LJ12-Encad shift", FALSE },
{ "COUL-Encad shift", TRUE },
{ "EXPMIN", FALSE },
- { "USER", FALSE }
+ { "USER", FALSE },
};
/* Index in the table that says which function to use */
#define pow4(x) ((x)*(x)*(x)*(x))
#define pow5(x) ((x)*(x)*(x)*(x)*(x))
-
-static double v_ewald_lr(double beta, double r)
+double v_q_ewald_lr(double beta, double r)
{
if (r == 0)
{
}
}
-void table_spline3_fill_ewald_lr(real *table_f,
- real *table_v,
- real *table_fdv0,
- int ntab,
- real dx,
- real beta)
+double v_lj_ewald_lr(double beta, double r)
+{
+ double br, br2, br4, r6, factor;
+ if (r == 0)
+ {
+ return pow(beta, 6)/6;
+ }
+ else
+ {
+ br = beta*r;
+ br2 = br*br;
+ br4 = br2*br2;
+ r6 = pow(r, 6.0);
+ factor = (1.0 - exp(-br2)*(1 + br2 + 0.5*br4))/r6;
+ return factor;
+ }
+}
+
+void table_spline3_fill_ewald_lr(real *table_f,
+ real *table_v,
+ real *table_fdv0,
+ int ntab,
+ double dx,
+ real beta,
+ real_space_grid_contribution_computer v_lr)
{
real tab_max;
int i, i_inrange;
double v_r0, v_r1, v_inrange, vi, a0, a1, a2dx;
double x_r0;
+ /* This function is called using either v_ewald_lr or v_lj_ewald_lr as a function argument
+ * depending on wether we should create electrostatic or Lennard-Jones Ewald tables.
+ */
+
if (ntab < 2)
{
gmx_fatal(FARGS, "Can not make a spline table with less than 2 points");
{
x_r0 = i*dx;
- v_r0 = v_ewald_lr(beta, x_r0);
+ v_r0 = (*v_lr)(beta, x_r0);
if (!bOutOfRange)
{
}
/* Get the potential at table point i-1 */
- v_r1 = v_ewald_lr(beta, (i-1)*dx);
+ v_r1 = (*v_lr)(beta, (i-1)*dx);
if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max)
{
/* Calculate the average second derivative times dx over interval i-1 to i.
* Using the function values at the end points and in the middle.
*/
- a2dx = (v_r0 + v_r1 - 2*v_ewald_lr(beta, x_r0-0.5*dx))/(0.25*dx);
+ a2dx = (v_r0 + v_r1 - 2*(*v_lr)(beta, x_r0-0.5*dx))/(0.25*dx);
/* Set the derivative of the spline to match the difference in potential
* over the interval plus the average effect of the quadratic term.
* This is the essential step for minimizing the error in the force.
}
}
+/* 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: 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);
+ }
- /* 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;
+ 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);
- return max(sc_f, sc_e);
+ 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 sc;
}
/* Calculate the potential and force for an r value
}
}
-static void init_table(FILE *fp, int n, int nx0,
+static void init_table(int n, int nx0,
double tabscale, t_tabledata *td, gmx_bool bAlloc)
{
int i;
for (k = 0; (k < ntab); k++)
{
- init_table(fp, nx, nx0, tabscale, &(td[k]), TRUE);
+ init_table(nx, nx0, tabscale, &(td[k]), TRUE);
for (i = 0; (i < nx); i++)
{
td[k].x[i] = yy[0][i];
sfree(td->f);
}
-static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr)
+static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
+ gmx_bool b14only)
{
/* Fill the table according to the formulas in the manual.
* In principle, we only need the potential and the second
int i;
double reppow, p;
double r1, rc, r12, r13;
- double r, r2, r6, rc6;
+ double r, r2, r6, rc2, rc6, rc12;
double expr, Vtab, Ftab;
/* Parameters for David's function */
double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
/* Parameters for the switching function */
double ksw, swi, swi1;
/* Temporary parameters */
- gmx_bool bSwitch, bShift;
- double ewc = fr->ewaldcoeff;
+ gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
+ double ewc = fr->ewaldcoeff_q;
+ double ewclj = fr->ewaldcoeff_lj;
+ double Vcut = 0;
- bSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
- (tp == etabCOULSwitch) ||
- (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch));
- bShift = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
- (tp == etabShift));
+ if (b14only)
+ {
+ bPotentialSwitch = FALSE;
+ bForceSwitch = FALSE;
+ bPotentialShift = FALSE;
+ }
+ else
+ {
+ bPotentialSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
+ (tp == etabCOULSwitch) ||
+ (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSWITCH)));
+ bForceSwitch = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
+ (tp == etabShift) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodFORCESWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodFORCESWITCH)));
+ bPotentialShift = ((tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSHIFT)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSHIFT)));
+ }
reppow = fr->reppow;
r1 = fr->rvdw_switch;
rc = fr->rvdw;
}
- if (bSwitch)
+ if (bPotentialSwitch)
{
ksw = 1.0/(pow5(rc-r1));
}
{
ksw = 0.0;
}
- if (bShift)
+ if (bForceSwitch)
{
if (tp == etabShift)
{
fp = xvgropen("switch.xvg", "switch", "r", "s");
#endif
+ if (bPotentialShift)
+ {
+ rc2 = rc*rc;
+ rc6 = 1.0/(rc2*rc2*rc2);
+ if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
+ {
+ rc12 = rc6*rc6;
+ }
+ else
+ {
+ rc12 = pow(rc, -reppow);
+ }
+
+ switch (tp)
+ {
+ case etabLJ6:
+ /* Dispersion */
+ Vcut = -rc6;
+ break;
+ case etabLJ6Ewald:
+ Vcut = -rc6*exp(-ewclj*ewclj*rc2)*(1 + ewclj*ewclj*rc2 + pow4(ewclj)*rc2*rc2/2);
+ break;
+ case etabLJ12:
+ /* Repulsion */
+ Vcut = rc12;
+ break;
+ case etabCOUL:
+ Vcut = 1.0/rc;
+ break;
+ case etabEwald:
+ case etabEwaldSwitch:
+ Vtab = gmx_erfc(ewc*rc)/rc;
+ break;
+ case etabEwaldUser:
+ /* Only calculate minus the reciprocal space contribution */
+ Vtab = -gmx_erf(ewc*rc)/rc;
+ break;
+ case etabRF:
+ case etabRF_ZERO:
+ /* No need for preventing the usage of modifiers with RF */
+ Vcut = 0.0;
+ break;
+ case etabEXPMIN:
+ Vcut = exp(-rc);
+ break;
+ default:
+ gmx_fatal(FARGS, "Cannot apply new potential-shift modifier to interaction type '%s' yet. (%s,%d)",
+ tprops[tp].name, __FILE__, __LINE__);
+ }
+ }
+
for (i = td->nx0; (i < td->nx); i++)
{
r = td->x[i];
}
Vtab = 0.0;
Ftab = 0.0;
- if (bSwitch)
+ if (bPotentialSwitch)
{
/* swi is function, swi1 1st derivative and swi2 2nd derivative */
/* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
break;
case etabEwaldUser:
case etabEwaldUserSwitch:
- /* Only calculate minus the reciprocal space contribution */
+ /* Only calculate the negative of the reciprocal space contribution */
Vtab = -gmx_erf(ewc*r)/r;
Ftab = -gmx_erf(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
break;
+ case etabLJ6Ewald:
+ Vtab = -r6*exp(-ewclj*ewclj*r2)*(1 + ewclj*ewclj*r2 + pow4(ewclj)*r2*r2/2);
+ Ftab = 6.0*Vtab/r - r6*exp(-ewclj*ewclj*r2)*pow5(ewclj)*ewclj*r2*r2*r;
+ break;
case etabRF:
case etabRF_ZERO:
Vtab = 1.0/r + fr->k_rf*r2 - fr->c_rf;
gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)",
tp, __FILE__, __LINE__);
}
- if (bShift)
+ if (bForceSwitch)
{
/* Normal coulomb with cut-off correction for potential */
if (r < rc)
Ftab += A*r12 + B*r13;
}
}
+ else
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ }
+ if (bPotentialShift)
+ {
+ if (r < rc)
+ {
+ Vtab -= Vcut;
+ }
+ else
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
}
if (ETAB_USER(tp))
Ftab += td->f[i];
}
- if ((r > r1) && bSwitch)
+ if (bPotentialSwitch)
{
- Ftab = Ftab*swi - Vtab*swi1;
- Vtab = Vtab*swi;
+ if (r >= rc)
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ else if (r > r1)
+ {
+ Ftab = Ftab*swi - Vtab*swi1;
+ Vtab = Vtab*swi;
+ }
}
-
/* Convert to single precision when we store to mem */
td->v[i] = Vtab;
td->f[i] = Ftab;
tabsel[etiLJ6] = etabLJ6Encad;
tabsel[etiLJ12] = etabLJ12Encad;
break;
+ case evdwPME:
+ tabsel[etiLJ6] = etabLJ6Ewald;
+ tabsel[etiLJ12] = etabLJ12;
+ break;
default:
gmx_fatal(FARGS, "Invalid vdwtype %d in %s line %d", vdwtype,
__FILE__, __LINE__);
}
+
+ if (!b14only && fr->vdw_modifier != eintmodNONE)
+ {
+ if (fr->vdw_modifier != eintmodPOTSHIFT &&
+ fr->vdwtype != evdwCUT)
+ {
+ gmx_incons("Potential modifiers other than potential-shift are only implemented for LJ cut-off");
+ }
+
+ /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
+ * to the original interaction forms when we fill the table, so we only check cutoffs here.
+ */
+ if (fr->vdwtype == evdwCUT)
+ {
+ switch (fr->vdw_modifier)
+ {
+ case eintmodNONE:
+ case eintmodPOTSHIFT:
+ case eintmodEXACTCUTOFF:
+ /* No modification */
+ break;
+ case eintmodPOTSWITCH:
+ tabsel[etiLJ6] = etabLJ6Switch;
+ tabsel[etiLJ12] = etabLJ12Switch;
+ break;
+ case eintmodFORCESWITCH:
+ tabsel[etiLJ6] = etabLJ6Shift;
+ tabsel[etiLJ12] = etabLJ12Shift;
+ break;
+ default:
+ gmx_incons("Unsupported vdw_modifier");
+ }
+ }
+ }
}
}
{
if (tabsel[k] != etabUSER)
{
- init_table(out, nx, nx0,
+ init_table(nx, nx0,
(tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale,
&(td[k]), !bReadTab);
- fill_table(&(td[k]), tabsel[k], fr);
+ fill_table(&(td[k]), tabsel[k], fr, b14only);
if (out)
{
fprintf(out, "%s table with %d data points for %s%s.\n"
return table;
}
-t_forcetable make_gb_table(FILE *out, const output_env_t oenv,
- const t_forcerec *fr,
- const char *fn,
- real rtab)
+t_forcetable make_gb_table(const output_env_t oenv,
+ const t_forcerec *fr)
{
const char *fns[3] = { "gbctab.xvg", "gbdtab.xvg", "gbrtab.xvg" };
const char *fns14[3] = { "gbctab14.xvg", "gbdtab14.xvg", "gbrtab14.xvg" };
snew_aligned(table.data, 4*nx, 32);
- init_table(out, nx, nx0, table.scale, &(td[0]), !bReadTab);
+ init_table(nx, nx0, table.scale, &(td[0]), !bReadTab);
/* Local implementation so we don't have to use the etabGB
* enum above, which will cause problems later when
for (i = nx0; i < nx; i++)
{
- Vtab = 0.0;
- Ftab = 0.0;
r = td->x[i];
r2 = r*r;
expterm = exp(-0.25*r2);
fprintf(fp, "%15.10e %15.10e %15.10e\n", x0, y0, yp);
}
- ffclose(fp);
+ gmx_ffclose(fp);
}
done_tabledata(&(td[0]));
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff