/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2012,2013, by the GROMACS development team, led by
+ * Copyright (c) 2012,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.
* 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 <assert.h>
+#include <math.h>
+#include <stdlib.h>
#include <sys/types.h>
-#include <math.h>
-#include "typedefs.h"
-#include "physics.h"
-#include "smalloc.h"
-#include "gmx_fatal.h"
-#include "macros.h"
-#include "vec.h"
-#include "coulomb.h"
+
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/math/units.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/legacyheaders/coulomb.h"
#include "calc_verletbuf.h"
#include "../mdlib/nbnxn_consts.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
+
#ifdef GMX_NBNXN_SIMD
/* The include below sets the SIMD instruction type (precision+width)
* for all nbnxn SIMD search and non-bonded kernel code.
#ifdef GMX_NBNXN_HALF_WIDTH_SIMD
#define GMX_USE_HALF_WIDTH_SIMD_HERE
#endif
-#include "gromacs/simd/macros.h"
+#include "gromacs/simd/simd.h"
#endif
#ifndef GMX_NBNXN_SIMD
list_setup->cluster_size_j = NBNXN_CPU_CLUSTER_I_SIZE;
#else
- list_setup->cluster_size_j = GMX_SIMD_WIDTH_HERE;
+ list_setup->cluster_size_j = GMX_SIMD_REAL_WIDTH;
#ifdef GMX_NBNXN_SIMD_2XNN
/* We assume the smallest cluster size to be on the safe side */
list_setup->cluster_size_j /= 2;
for (i = 0; i < il->nr; i += 1+NRAL(ft))
{
const t_iparams *ip;
- real cam[5], inv_mass, m_aj;
+ real cam[5] = {0}, inv_mass, m_aj;
int a1, j, aj, coeff;
ip = &ffparams->iparams[il->iatoms[i]];
inv_mass += coeff*coeff/m_aj;
}
vsite_m[a1] = 1/inv_mass;
+ /* Correct for loop increment of i */
+ i += j - 1 - NRAL(ft);
break;
default:
/* Use the mass of the lightest constructing atom.
add_at(&att, &natt, &prop[a], nmol);
}
+ /* cppcheck-suppress uninitvar Fixed in cppcheck 1.65 */
sfree(vsite_m);
sfree(prop);
}
static real ener_drift(const verletbuf_atomtype_t *att, int natt,
const gmx_ffparams_t *ffp,
real kT_fac,
- real md_ljd, real md_ljr, real md_el, real dd_el,
+ real md1_ljd, real d2_ljd, real md3_ljd,
+ real md1_ljr, real d2_ljr, real md3_ljr,
+ real md1_el, real d2_el,
real r_buffer,
real rlist, real boxvol)
{
- double drift_tot, pot1, pot2, pot;
+ double drift_tot, pot1, pot2, pot3, pot;
int i, j;
real s2i_2d, s2i_3d, s2j_2d, s2j_3d, s2, s;
int ti, tj;
- real md, dd;
- real sc_fac, rsh;
+ real md1, d2, md3;
+ real sc_fac, rsh, rsh2;
double c_exp, c_erfc;
drift_tot = 0;
* pairs will partially cancel.
*/
/* -dV/dr at the cut-off for LJ + Coulomb */
- md =
- md_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
- md_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
- md_el*att[i].prop.q*att[j].prop.q;
-
- /* d2V/dr2 at the cut-off for Coulomb, we neglect LJ */
- dd = dd_el*att[i].prop.q*att[j].prop.q;
+ md1 =
+ md1_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
+ md1_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
+ md1_el*att[i].prop.q*att[j].prop.q;
+
+ /* d2V/dr2 at the cut-off for LJ + Coulomb */
+ d2 =
+ d2_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
+ d2_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12 +
+ d2_el*att[i].prop.q*att[j].prop.q;
+
+ /* -d3V/dr3 at the cut-off for LJ, we neglect Coulomb */
+ md3 =
+ md3_ljd*ffp->iparams[ti*ffp->atnr+tj].lj.c6 +
+ md3_ljr*ffp->iparams[ti*ffp->atnr+tj].lj.c12;
rsh = r_buffer;
sc_fac = 1.0;
c_exp = exp(-rsh*rsh/(2*s2))/sqrt(2*M_PI);
c_erfc = 0.5*gmx_erfc(rsh/(sqrt(2*s2)));
s = sqrt(s2);
+ rsh2 = rsh*rsh;
pot1 = sc_fac*
- md/2*((rsh*rsh + s2)*c_erfc - rsh*s*c_exp);
+ md1/2*((rsh2 + s2)*c_erfc - rsh*s*c_exp);
pot2 = sc_fac*
- dd/6*(s*(rsh*rsh + 2*s2)*c_exp - rsh*(rsh*rsh + 3*s2)*c_erfc);
- pot = pot1 + pot2;
+ d2/6*(s*(rsh2 + 2*s2)*c_exp - rsh*(rsh2 + 3*s2)*c_erfc);
+ pot3 =
+ md3/24*((rsh2*rsh2 + 6*rsh2*s2 + 3*s2*s2)*c_erfc - rsh*s*(rsh2 + 5*s2)*c_exp);
+ pot = pot1 + pot2 + pot3;
if (gmx_debug_at)
{
- fprintf(debug, "n %d %d d s %.3f %.3f %.3f %.3f con %d md %8.1e dd %8.1e pot1 %8.1e pot2 %8.1e pot %8.1e\n",
+ fprintf(debug, "n %d %d d s %.3f %.3f %.3f %.3f con %d -d1 %8.1e d2 %8.1e -d3 %8.1e pot1 %8.1e pot2 %8.1e pot3 %8.1e pot %8.1e\n",
att[i].n, att[j].n,
sqrt(s2i_2d), sqrt(s2i_3d),
sqrt(s2j_2d), sqrt(s2j_3d),
att[i].prop.bConstr+att[j].prop.bConstr,
- md, dd, pot1, pot2, pot);
+ md1, d2, md3,
+ pot1, pot2, pot3, pot);
}
/* Multiply by the number of atom pairs */
return area_rel/cluster_size;
}
+/* Returns the negative of the third derivative of a potential r^-p
+ * with a force-switch function, evaluated at the cut-off rc.
+ */
+static real md3_force_switch(real p, real rswitch, real rc)
+{
+ /* The switched force function is:
+ * p*r^-(p+1) + a*(r - rswitch)^2 + b*(r - rswitch)^3
+ */
+ real a, b;
+ real md3_pot, md3_sw;
+
+ a = -((p + 4)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*pow(rc-rswitch, 2));
+ b = ((p + 3)*rc - (p + 1)*rswitch)/(pow(rc, p+2)*pow(rc-rswitch, 3));
+
+ md3_pot = (p + 2)*(p + 1)*p*pow(rc, p+3);
+ md3_sw = 2*a + 6*b*(rc - rswitch);
+
+ return md3_pot + md3_sw;
+}
+
void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
const t_inputrec *ir,
+ real reference_temperature,
const verletbuf_list_setup_t *list_setup,
int *n_nonlin_vsite,
real *rlist)
verletbuf_atomtype_t *att = NULL;
int natt = -1, i;
double reppow;
- real md_ljd, md_ljr, md_el, dd_el;
+ real md1_ljd, d2_ljd, md3_ljd;
+ real md1_ljr, d2_ljr, md3_ljr;
+ real md1_el, d2_el;
real elfac;
real kT_fac, mass_min;
int ib0, ib1, ib;
real rb, rl;
real drift;
+ if (reference_temperature < 0)
+ {
+ if (EI_MD(ir->eI) && ir->etc == etcNO)
+ {
+ /* This case should be handled outside calc_verlet_buffer_size */
+ gmx_incons("calc_verlet_buffer_size called with an NVE ensemble and reference_temperature < 0");
+ }
+
+ /* We use the maximum temperature with multiple T-coupl groups.
+ * We could use a per particle temperature, but since particles
+ * interact, this might underestimate the buffer size.
+ */
+ reference_temperature = 0;
+ for (i = 0; i < ir->opts.ngtc; i++)
+ {
+ if (ir->opts.tau_t[i] >= 0)
+ {
+ reference_temperature = max(reference_temperature,
+ ir->opts.ref_t[i]);
+ }
+ }
+ }
+
/* Resolution of the buffer size */
resolution = 0.001;
fprintf(debug, "energy drift atom types: %d\n", natt);
}
- reppow = mtop->ffparams.reppow;
- md_ljd = 0;
- md_ljr = 0;
+ reppow = mtop->ffparams.reppow;
+ md1_ljd = 0;
+ d2_ljd = 0;
+ md3_ljd = 0;
+ md1_ljr = 0;
+ d2_ljr = 0;
+ md3_ljr = 0;
if (ir->vdwtype == evdwCUT)
{
- /* -dV/dr of -r^-6 and r^-repporw */
- md_ljd = -6*pow(ir->rvdw, -7.0);
- md_ljr = reppow*pow(ir->rvdw, -(reppow+1));
- /* The contribution of the second derivative is negligible */
+ real sw_range, md3_pswf;
+
+ switch (ir->vdw_modifier)
+ {
+ case eintmodNONE:
+ case eintmodPOTSHIFT:
+ /* -dV/dr of -r^-6 and r^-reppow */
+ md1_ljd = -6*pow(ir->rvdw, -7.0);
+ md1_ljr = reppow*pow(ir->rvdw, -(reppow+1));
+ /* The contribution of the higher derivatives is negligible */
+ break;
+ case eintmodFORCESWITCH:
+ /* At the cut-off: V=V'=V''=0, so we use only V''' */
+ md3_ljd = -md3_force_switch(6.0, ir->rvdw_switch, ir->rvdw);
+ md3_ljr = md3_force_switch(reppow, ir->rvdw_switch, ir->rvdw);
+ break;
+ case eintmodPOTSWITCH:
+ /* At the cut-off: V=V'=V''=0.
+ * V''' is given by the original potential times
+ * the third derivative of the switch function.
+ */
+ sw_range = ir->rvdw - ir->rvdw_switch;
+ md3_pswf = 60.0*pow(sw_range, -3.0);
+
+ md3_ljd = -pow(ir->rvdw, -6.0 )*md3_pswf;
+ md3_ljr = pow(ir->rvdw, -reppow)*md3_pswf;
+ break;
+ default:
+ gmx_incons("Unimplemented VdW modifier");
+ }
}
else if (EVDW_PME(ir->vdwtype))
{
br4 = br2*br2;
br6 = br4*br2;
/* -dV/dr of g(br)*r^-6 [where g(x) = exp(-x^2)(1+x^2+x^4/2), see LJ-PME equations in manual] and r^-reppow */
- md_ljd = -exp(-br2)*(br6 + 3.0*br4 + 6.0*br2 + 6.0)*pow(r, -7.0);
- md_ljr = reppow*pow(r, -(reppow+1));
- /* The contribution of the second derivative is negligible */
+ md1_ljd = -exp(-br2)*(br6 + 3.0*br4 + 6.0*br2 + 6.0)*pow(r, -7.0);
+ md1_ljr = reppow*pow(r, -(reppow+1));
+ /* The contribution of the higher derivatives is negligible */
}
else
{
elfac = ONE_4PI_EPS0/ir->epsilon_r;
/* Determine md=-dV/dr and dd=d^2V/dr^2 */
- md_el = 0;
- dd_el = 0;
+ md1_el = 0;
+ d2_el = 0;
if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
{
real eps_rf, k_rf;
if (eps_rf > 0)
{
- md_el = elfac*(pow(ir->rcoulomb, -2.0) - 2*k_rf*ir->rcoulomb);
+ md1_el = elfac*(pow(ir->rcoulomb, -2.0) - 2*k_rf*ir->rcoulomb);
}
- dd_el = elfac*(2*pow(ir->rcoulomb, -3.0) + 2*k_rf);
+ d2_el = elfac*(2*pow(ir->rcoulomb, -3.0) + 2*k_rf);
}
else if (EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD)
{
real b, rc, br;
- b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
- rc = ir->rcoulomb;
- br = b*rc;
- md_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + gmx_erfc(br)/(rc*rc));
- dd_el = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*gmx_erfc(br)/rc);
+ b = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
+ rc = ir->rcoulomb;
+ br = b*rc;
+ md1_el = elfac*(b*exp(-br*br)*M_2_SQRTPI/rc + gmx_erfc(br)/(rc*rc));
+ d2_el = elfac/(rc*rc)*(2*b*(1 + br*br)*exp(-br*br)*M_2_SQRTPI + 2*gmx_erfc(br)/rc);
}
else
{
* should be negligible (unless nstlist is extremely large, which
* you wouldn't do anyhow).
*/
- kT_fac = 2*BOLTZ*ir->opts.ref_t[0]*(ir->nstlist-1)*ir->delta_t;
+ kT_fac = 2*BOLTZ*reference_temperature*(ir->nstlist-1)*ir->delta_t;
if (ir->bd_fric > 0)
{
/* This is directly sigma^2 of the displacement */
}
else
{
- kT_fac = BOLTZ*ir->opts.ref_t[0]*sqr((ir->nstlist-1)*ir->delta_t);
+ kT_fac = BOLTZ*reference_temperature*sqr((ir->nstlist-1)*ir->delta_t);
}
mass_min = att[0].prop.mass;
if (debug)
{
- fprintf(debug, "md_ljd %e md_ljr %e\n", md_ljd, md_ljr);
- fprintf(debug, "md_el %e dd_el %e\n", md_el, dd_el);
+ fprintf(debug, "md1_ljd %9.2e d2_ljd %9.2e md3_ljd %9.2e\n", md1_ljd, d2_ljd, md3_ljd);
+ fprintf(debug, "md1_ljr %9.2e d2_ljr %9.2e md3_ljr %9.2e\n", md1_ljr, d2_ljr, md3_ljr);
+ fprintf(debug, "md1_el %9.2e d2_el %9.2e\n", md1_el, d2_el);
fprintf(debug, "sqrt(kT_fac) %f\n", sqrt(kT_fac));
fprintf(debug, "mass_min %f\n", mass_min);
}
*/
drift = ener_drift(att, natt, &mtop->ffparams,
kT_fac,
- md_ljd, md_ljr, md_el, dd_el, rb,
+ md1_ljd, d2_ljd, md3_ljd,
+ md1_ljr, d2_ljr, md3_ljr,
+ md1_el, d2_el,
+ rb,
rl, boxvol);
/* Correct for the fact that we are using a Ni x Nj particle pair list
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff