the GROMACS pair-list setup leads to a reduction in the drift by
a factor 10, compared to a simple particle-pair based list.
Without dynamics (energy minimization etc.), the buffer is 5% of the cut-off.
-For dynamics without temperature coupling or to override the buffer size,
+For NVE simulations the initial temperature is used, unless this is zero,
+in which case a buffer of 10% is used. For NVE simulations the tolerance
+usually needs to be lowered to achieve proper energy conservation on
+the nanosecond time scale. To override the automated buffer setting,
use <b>verlet-buffer-tolerance</b>=-1 and set <b>rlist</b> manually.</dd>
<dt><b>rlist: (1) [nm]</b></dt>
/*
* 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.
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)
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;
* 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;
} verletbuf_list_setup_t;
+/* Add a 5% and 10% rlist buffer for simulations without dynamics (EM, NM, ...)
+ * and NVE simulations with zero initial temperature, respectively.
+ * 10% should be enough for any NVE simulation with PME and nstlist=10,
+ * for other settings it might not be enough, but then it's difficult
+ * to come up with any reasonable (not crazily expensive) value
+ * and grompp will notify the user when using the 10% buffer.
+ */
+static const real verlet_buffer_ratio_nodynamics = 0.05;
+static const real verlet_buffer_ratio_NVE_T0 = 0.10;
+
+
/* Sets the pair-list setup assumed for the current Gromacs configuration.
* The setup with smallest cluster sizes is return, such that the Verlet
* buffer size estimated with this setup will be conservative.
/* Calculate the non-bonded pair-list buffer size for the Verlet list
* based on the particle masses, temperature, LJ types, charges
* and constraints as well as the non-bonded force behavior at the cut-off.
+ * If reference_temperature < 0, the maximum coupling temperature will be used.
* The target is a maximum energy drift of ir->verletbuf_tol.
* Returns the number of non-linear virtual sites. For these it's difficult
* to determine their contribution to the drift exaclty, so we approximate.
*/
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);
}
-static void set_verlet_buffer(const gmx_mtop_t *mtop,
- t_inputrec *ir,
- matrix box,
- warninp_t wi)
+static real calc_temp(const gmx_mtop_t *mtop,
+ const t_inputrec *ir,
+ rvec *v)
{
- real ref_T;
- int i;
- verletbuf_list_setup_t ls;
- real rlist_1x1;
- int n_nonlin_vsite;
- char warn_buf[STRLEN];
+ double sum_mv2;
+ gmx_mtop_atomloop_all_t aloop;
+ t_atom *atom;
+ int a;
+ int nrdf, g;
+
+ sum_mv2 = 0;
+
+ aloop = gmx_mtop_atomloop_all_init(mtop);
+ while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
+ {
+ sum_mv2 += atom->m*norm2(v[a]);
+ }
- ref_T = 0;
+ nrdf = 0;
+ for (g = 0; g < ir->opts.ngtc; g++)
+ {
+ nrdf += ir->opts.nrdf[g];
+ }
+
+ return sum_mv2/(nrdf*BOLTZ);
+}
+
+static real get_max_reference_temp(const t_inputrec *ir,
+ warninp_t wi)
+{
+ real ref_t;
+ int i;
+ gmx_bool bNoCoupl;
+
+ ref_t = 0;
+ bNoCoupl = FALSE;
for (i = 0; i < ir->opts.ngtc; i++)
{
- if (ir->opts.ref_t[i] < 0)
+ if (ir->opts.tau_t[i] < 0)
{
- warning(wi, "Some atom groups do not use temperature coupling. This cannot be accounted for in the energy error estimation for the Verlet buffer size. The energy error and the Verlet buffer might be underestimated.");
+ bNoCoupl = TRUE;
}
else
{
- ref_T = max(ref_T, ir->opts.ref_t[i]);
+ ref_t = max(ref_t, ir->opts.ref_t[i]);
}
}
- printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, ref_T);
-
- for (i = 0; i < ir->opts.ngtc; i++)
+ if (bNoCoupl)
{
- if (ir->opts.ref_t[i] >= 0 && ir->opts.ref_t[i] != ref_T)
- {
- sprintf(warn_buf, "ref_T for group of %.1f DOFs is %g K, which is smaller than the maximum of %g K used for the buffer size calculation. The buffer size might be on the conservative (large) side.",
- ir->opts.nrdf[i], ir->opts.ref_t[i], ref_T);
- warning_note(wi, warn_buf);
- }
+ char buf[STRLEN];
+
+ sprintf(buf, "Some temperature coupling groups do not use temperature coupling. We will assume their temperature is not more than %.3f K. If their temperature is higher, the energy error and the Verlet buffer might be underestimated.",
+ ref_t);
+ warning(wi, buf);
}
+ return ref_t;
+}
+
+static void set_verlet_buffer(const gmx_mtop_t *mtop,
+ t_inputrec *ir,
+ real buffer_temp,
+ matrix box,
+ warninp_t wi)
+{
+ int i;
+ verletbuf_list_setup_t ls;
+ real rlist_1x1;
+ int n_nonlin_vsite;
+ char warn_buf[STRLEN];
+
+ printf("Determining Verlet buffer for a tolerance of %g kJ/mol/ps at %g K\n", ir->verletbuf_tol, buffer_temp);
+
/* Calculate the buffer size for simple atom vs atoms list */
ls.cluster_size_i = 1;
ls.cluster_size_j = 1;
- calc_verlet_buffer_size(mtop, det(box), ir,
+ calc_verlet_buffer_size(mtop, det(box), ir, buffer_temp,
&ls, &n_nonlin_vsite, &rlist_1x1);
/* Set the pair-list buffer size in ir */
verletbuf_get_list_setup(FALSE, &ls);
- calc_verlet_buffer_size(mtop, det(box), ir,
+ calc_verlet_buffer_size(mtop, det(box), ir, buffer_temp,
&ls, &n_nonlin_vsite, &ir->rlist);
if (n_nonlin_vsite > 0)
if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 &&
ir->nstlist > 1)
{
- if (EI_DYNAMICS(ir->eI) &&
- !(EI_MD(ir->eI) && ir->etc == etcNO) &&
- inputrec2nboundeddim(ir) == 3)
+ if (EI_DYNAMICS(ir->eI) && inputrec2nboundeddim(ir) == 3)
{
- set_verlet_buffer(sys, ir, state.box, wi);
+ real buffer_temp;
+
+ if (EI_MD(ir->eI) && ir->etc == etcNO)
+ {
+ if (bGenVel)
+ {
+ buffer_temp = opts->tempi;
+ }
+ else
+ {
+ buffer_temp = calc_temp(sys, ir, state.v);
+ }
+ if (buffer_temp > 0)
+ {
+ sprintf(warn_buf, "NVE simulation: will use the initial temperature of %.3f K for determining the Verlet buffer size", buffer_temp);
+ warning_note(wi, warn_buf);
+ }
+ else
+ {
+ sprintf(warn_buf, "NVE simulation with an initial temperature of zero: will use a Verlet buffer of %d%%. Check your energy drift!",
+ (int)(verlet_buffer_ratio_NVE_T0*100 + 0.5));
+ warning_note(wi, warn_buf);
+ }
+ }
+ else
+ {
+ buffer_temp = get_max_reference_temp(ir, wi);
+ }
+
+ if (EI_MD(ir->eI) && ir->etc == etcNO && buffer_temp == 0)
+ {
+ /* NVE with initial T=0: we add a fixed ratio to rlist.
+ * Since we don't actually use verletbuf_tol, we set it to -1
+ * so it can't be misused later.
+ */
+ ir->rlist *= 1.0 + verlet_buffer_ratio_NVE_T0;
+ ir->verletbuf_tol = -1;
+ }
+ else
+ {
+ /* We warn for NVE simulations with >1(.1)% drift tolerance */
+ const real drift_tol = 0.01;
+ real ener_runtime;
+
+ /* We use 2 DOF per atom = 2kT pot+kin energy, to be on
+ * the safe side with constraints (without constraints: 3 DOF).
+ */
+ ener_runtime = 2*BOLTZ*buffer_temp/(ir->nsteps*ir->delta_t);
+
+ if (EI_MD(ir->eI) && ir->etc == etcNO && ir->nstlist > 1 &&
+ ir->nsteps > 0 &&
+ ir->verletbuf_tol > 1.1*drift_tol*ener_runtime)
+ {
+ sprintf(warn_buf, "You are using a Verlet buffer tolerance of %g kJ/mol/ps for an NVE simulation of length %g ps, which can give a final drift of %d%%. For conserving energy to %d%%, you might need to set verlet-buffer-tolerance to %.1e.",
+ ir->verletbuf_tol, ir->nsteps*ir->delta_t,
+ (int)(ir->verletbuf_tol/ener_runtime*100 + 0.5),
+ (int)(100*drift_tol + 0.5),
+ drift_tol*ener_runtime);
+ warning_note(wi, warn_buf);
+ }
+
+ set_verlet_buffer(sys, ir, buffer_temp, state.box, wi);
+ }
}
}
#include "mtop_util.h"
#include "chargegroup.h"
#include "inputrec.h"
+#include "calc_verletbuf.h"
#define MAXPTR 254
#define NOGID 255
{
if (EI_DYNAMICS(ir->eI))
{
- if (EI_MD(ir->eI) && ir->etc == etcNO)
- {
- warning_error(wi, "Temperature coupling is required for calculating rlist using the energy tolerance with verlet-buffer-tolerance > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-tolerance = -1.");
- }
-
if (inputrec2nboundeddim(ir) < 3)
{
warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
else
{
/* Set the buffer to 5% of the cut-off */
- ir->rlist = 1.05*rc_max;
+ ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
}
}
}
t_state state_tmp;
gmx_bool bBox, bDD, bCont;
const char *nstl_gpu = "\nFor optimal performance with a GPU nstlist (now %d) should be larger.\nThe optimum depends on your CPU and GPU resources.\nYou might want to try several nstlist values.\n";
+ const char *nve_err = "Can not increase nstlist because an NVE ensemble is used";
const char *vbd_err = "Can not increase nstlist because verlet-buffer-tolerance is not set or used";
const char *box_err = "Can not increase nstlist because the box is too small";
const char *dd_err = "Can not increase nstlist because of domain decomposition limitations";
}
}
+ if (EI_MD(ir->eI) && ir->etc == etcNO)
+ {
+ if (MASTER(cr))
+ {
+ fprintf(stderr, "%s\n", nve_err);
+ }
+ if (fp != NULL)
+ {
+ fprintf(fp, "%s\n", nve_err);
+ }
+
+ return;
+ }
+
if (ir->verletbuf_tol == 0 && bGPU)
{
gmx_fatal(FARGS, "You are using an old tpr file with a GPU, please generate a new tpr file with an up to date version of grompp");
*/
nstlist_prev = ir->nstlist;
ir->nstlist = 10;
- calc_verlet_buffer_size(mtop, det(box), ir, &ls, NULL, &rlist_nstlist10);
+ calc_verlet_buffer_size(mtop, det(box), ir, -1, &ls, NULL,
+ &rlist_nstlist10);
ir->nstlist = nstlist_prev;
/* Determine the pair list size increase due to zero interactions */
}
/* Set the pair-list buffer size in ir */
- calc_verlet_buffer_size(mtop, det(box), ir, &ls, NULL, &rlist_new);
+ calc_verlet_buffer_size(mtop, det(box), ir, -1, &ls, NULL, &rlist_new);
/* Does rlist fit in the box? */
bBox = (sqr(rlist_new) < max_cutoff2(ir->ePBC, box));
matrix box,
gmx_bool bUseGPU)
{
- if (ir->verletbuf_tol > 0)
+ /* For NVE simulations, we will retain the initial list buffer */
+ if (ir->verletbuf_tol > 0 && !(EI_MD(ir->eI) && ir->etc == etcNO))
{
/* Update the Verlet buffer size for the current run setup */
verletbuf_list_setup_t ls;
*/
verletbuf_get_list_setup(bUseGPU, &ls);
- calc_verlet_buffer_size(mtop, det(box), ir, &ls, NULL, &rlist_new);
+ calc_verlet_buffer_size(mtop, det(box), ir, -1, &ls, NULL, &rlist_new);
if (rlist_new != ir->rlist)
{
verletbuf_list_setup_t ls;
verletbuf_get_list_setup(FALSE, &ls);
- calc_verlet_buffer_size(mtop, box_vol, ir, &ls, NULL, &ir->rlist);
+ calc_verlet_buffer_size(mtop, box_vol, ir, -1, &ls, NULL, &ir->rlist);
}
else
{
+ real rlist_fac;
+
+ if (EI_MD(ir->eI))
+ {
+ rlist_fac = 1 + verlet_buffer_ratio_NVE_T0;
+ }
+ else
+ {
+ rlist_fac = 1 + verlet_buffer_ratio_nodynamics;
+ }
ir->verletbuf_tol = -1;
- ir->rlist = 1.05*max(ir->rvdw, ir->rcoulomb);
+ ir->rlist = rlist_fac*max(ir->rvdw, ir->rcoulomb);
}
gmx_mtop_remove_chargegroups(mtop);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff