Old mdp files are still processed correctly.
Also let mdrun tune nstlist with CPU-only runs. Added mdrun option
-nstlist to set nstlist manually.
Change-Id: Iac2c2b480c1f7f4233e57af6740424715c1bd1f1
The Verlet scheme specific settings in the {\tt mdp} file are:
\begin{verbatim}
-cutoff-scheme = Verlet
-verlet-buffer-drift = 0.005
+cutoff-scheme = Verlet
+verlet-buffer-tolerance = 0.005
\end{verbatim}
The Verlet buffer size is determined from the latter option, which is
-by default set to 0.005 kJ/mol/ps energy drift per atom. Note that the
-total energy drift in the system is affected by many factors and it is
-usually much smaller than this default setting for the estimate. For
-constant energy (NVE) simulations, this drift should be set to -1 and
-a buffer has to be set manually by specifying {\tt rlist} $>$ {\tt
+by default set to 0.005 kJ/mol/ps pair energy error per atom. Note that
+errors in pair energies cancel and the effect on the total energy drift
+is usually at least an order of magnitude smaller than the tolerance.
+Furthermore, the drift of the total energy is affected by many other
+factors, the constraint contribution is often the dominating one.
+For constant energy (NVE) simulations, this drift should be set to -1
+and a buffer has to be set manually by specifying {\tt rlist} $>$ {\tt
rcoulomb}. The simplest way to get a reasonable buffer size is to
use an NVT {\tt mdp} file with the target temperature set to what you
expect in your NVE simulation, and transfer the buffer size printed by
\ifthenelse{\equal{\gmxlite}{1}}{}{
\subsubsection{Energy drift and pair-list buffering}
-For a canonical ensemble, the average energy drift caused by the
+For a canonical ensemble, the average energy error caused by the
finite Verlet buffer size can be determined from the atomic
displacements and the shape of the potential at the cut-off.
%Since we are interested in the small drift regime, we will assume
interact with other particles over time $t$ and therefore the real
displacement distribution is much narrower. Given a non-bonded
interaction cut-off distance of $r_c$ and a pair-list cut-off
-$r_\ell=r_c+r_b$, we can then write the average energy drift after
+$r_\ell=r_c+r_b$, we can then write the average energy error after
time $t$ for pair interactions between one particle of type 1
surrounded by particles of type 2 with number density $\rho_2$, when
the inter particle distance changes from $r_0$ to $r_t$, as:
where $G$ is a Gaussian distribution with 0 mean and unit variance and
$E(x)=\frac{1}{2}\mathrm{erfc}(x/\sqrt{2})$. We always want to achieve
-small energy drift, so $\sigma$ will be small compared to both $r_c$
+small energy error, so $\sigma$ will be small compared to both $r_c$
and $r_\ell$, thus the approximations in the equations above are good,
-since the Gaussian distribution decays rapidly. The energy drift needs
+since the Gaussian distribution decays rapidly. The energy error needs
to be averaged over all particle pair types and weighted with the
-particle counts. In {\gromacs} we don't allow cancellation of drift
+particle counts. In {\gromacs} we don't allow cancellation of error
between pair types, so we average the absolute values. To obtain the
-average energy drift per unit time, it needs to be divided by the
+average energy error per unit time, it needs to be divided by the
neighbor-list life time $t = ({\tt nstlist} - 1)\times{\tt dt}$. This
function can not be inverted analytically, so we use bisection to
obtain the buffer size $r_b$ for a target drift. Again we note that
-in practice the drift we usually be much smaller than this estimate,
+in practice the error we usually be much smaller than this estimate,
as in the condensed phase particle displacements will be much smaller
than for freely moving particles, which is the assumption used here.
When (bond) constraints are present, some particles will have fewer
-degrees of freedom. This will reduce the energy drift. The
+degrees of freedom. This will reduce the energy errors. The
displacement in an arbitrary direction of a particle with 2 degrees of
freedom is not Gaussian, but rather follows the complementary error
function: \beq
\frac{\sqrt{\pi}}{2\sqrt{2}\sigma}\,\mathrm{erfc}\left(\frac{|r|}{\sqrt{2}\,\sigma}\right)
\eeq where $\sigma^2$ is again $k_B T/m$. This distribution can no
-longer be integrated analytically to obtain the energy drift. But we
+longer be integrated analytically to obtain the energy error. But we
can generate a tight upper bound using a scaled and shifted Gaussian
distribution (not shown). This Gaussian distribution can then be used
-to calculate the energy drift as described above. We consider
+to calculate the energy error as described above. We consider
particles constrained, i.e. having 2 degrees of freedom or fewer, when
they are connected by constraints to particles with a total mass of at
least 1.5 times the mass of the particles itself. For a particle with
freedom.
There is one important implementation detail that reduces the energy
-drift caused by the finite Verlet buffer list size. The derivation
+errors caused by the finite Verlet buffer list size. The derivation
above assumes a particle pair-list. However, the {\gromacs}
implementation uses a cluster pair-list for efficiency. The pair list
consists of pairs of clusters of 4 particles in most cases, also
reasonable assumptions. The fraction decreases with increasing
pair-list range, meaning that a smaller buffer can be used. For
typical all-atom simulations with a cut-off of 0.9 nm this fraction is
-around 0.9, which gives a reduction in the energy drift of a factor of
+around 0.9, which gives a reduction in the energy errors of a factor of
10. This reduction is taken into account during the automatic Verlet
buffer calculation and results in a smaller buffer size.
a time step of 2 fs and a pair-list update period of 10 steps
(pair-list life time: 18 fs). PME was used with {\tt ewald-rtol} set
to 10$^{-5}$; this parameter affects the shape of the potential at
- the cut-off. Drift estimates due to finite Verlet buffer size are
+ the cut-off. Error estimates due to finite Verlet buffer size are
shown for a $1 \times 1$ atom pair list and $4 \times 4$ atom pair
list without and with (dashed line) cancellation of positive and
- negative drift. Real energy drift is shown for double- and
+ negative errors. Real energy drift is shown for double- and
single-precision simulations. Single-precision rounding errors in
the SETTLE constraint algorithm cause the drift to become negative
at large buffer size. Note that at zero buffer size, the real drift
- is small because the positive (H-H) and negative (O-H) drift
- cancels.}
+ is small because positive (H-H) and negative (O-H) energy errors
+ cancel.}
\label{fig:verletdrift}
\end{figure}
-In \figref{verletdrift} one can see that for water with a pair-list
-life time of 18 fs, the drift estimate is a factor of 6 higher than
-the real drift, or alternatively the buffer estimate is 0.024 nm too
-large. This is because the protons don't move freely over 18 fs, but
-rather vibrate.
+In \figref{verletdrift} one can see that for small buffer sizes the drift
+of the total energy is much smaller than the pair energy error tolerance,
+due to cancellation of errors. For larger buffer size, the error estimate
+is a factor of 6 higher than drift of the total energy, or alternatively
+the buffer estimate is 0.024 nm too large. This is because the protons
+don't move freely over 18 fs, but rather vibrate.
%At a buffer size of zero there is cancellation of
%drift due to repulsive (H-H) and attractive (O-H) interactions.
<li><a HREF="#shellmd"><b>shell molecular dynamics</b></a>(emtol,niter,fcstep)
<li><a HREF="#tpi"><b>test particle insertion</b></a>(rtpi)
<li><A HREF="#out"><b>output control</b></A> (nstxout, nstvout, nstfout, nstlog, nstcalcenergy, nstenergy, nstxtcout, xtc-precision, xtc-grps, energygrps)
-<li><A HREF="#nl"><b>neighbor searching</b></A> (cutoff-scheme, nstlist, nstcalclr, ns-type, pbc, periodic-molecules, verlet-buffer-drift, rlist, rlistlong)
+<li><A HREF="#nl"><b>neighbor searching</b></A> (cutoff-scheme, nstlist, nstcalclr, ns-type, pbc, periodic-molecules, verlet-buffer-tolerance, rlist, rlistlong)
<li><A HREF="#el"><b>electrostatics</b></A> (coulombtype, coulomb-modifier, rcoulomb-switch, rcoulomb, epsilon-r, epsilon-rf)
<li><A HREF="#vdw"><b>VdW</b></A> (vdwtype, vdw-modifier, rvdw-switch, rvdw, DispCorr)
<li><A HREF="#table"><b>tables</b></A> (table-extension, energygrp-table)
<dt><b>Verlet</b></dt>
<dd>Generate a pair list with buffering. The buffer size is automatically set
-based on <b>verlet-buffer-drift</b>, unless this is set to -1, in which case
+based on <b>verlet-buffer-tolerance</b>, unless this is set to -1, in which case
<b>rlist</b> will be used. This option has an explicit, exact cut-off at
<b>rvdw</b>=<b>rcoulomb</b>. Currently only cut-off, reaction-field,
PME electrostatics and plain LJ are supported. Some <tt>mdrun</tt> functionality
the neighbor list is made only once.
With energy minimization the neighborlist will be updated for every
energy evaluation when <b>nstlist</b><tt>>0</tt>.
-With non-bonded force calculation on the GPU, a value of 20 or more gives
-the best performance.</dd>
+With <b>cutoff-scheme=Verlet</b> and <b>verlet-buffer-tolerance</b> set,
+<b>nstlist</b> is actually a minimum value and <tt>mdrun</tt> might increase it.
+With parallel simulations and/or non-bonded force calculation on the GPU,
+a value of 20 or 40 often gives the best performance.
+With <b>cutoff-scheme=Group</b> and non-exact cut-off's, <b>nstlist</b> will
+affect the accuracy of your simulation and it can not be chosen freely.
+</dd>
<dt><b>0</b></dt>
<dd>The neighbor list is only constructed once and never updated.
This is mainly useful for vacuum simulations in which all particles
and molecules are not made whole in the output</dd>
</dl></dd>
-<dt><b>verlet-buffer-drift: (0.005) [kJ/mol/ps]</b></dt>
-<dd>Useful only with <b>cutoff-scheme</b>=<b>Verlet</b>. This sets the target energy drift
-per particle caused by the Verlet buffer, which indirectly sets <b>rlist</b>.
+<dt><b>verlet-buffer-tolerance: (0.005) [kJ/mol/ps]</b></dt>
+<dd>Useful only with <b>cutoff-scheme</b>=<b>Verlet</b>. This sets the maximum
+allowed error for pair interactions per particle caused by the Verlet buffer,
+which indirectly sets <b>rlist</b>.
As both <b>nstlist</b> and the Verlet buffer size are fixed
(for performance reasons), particle pairs not in the pair list can occasionally
get within the cut-off distance during <b>nstlist</b>-1 nsteps. This
-generates energy drift. In a constant-temperature ensemble, the drift can be
+causes very small jumps in the energy. In a constant-temperature ensemble,
+these very small energy jumps can be
estimated for a given cut-off and <b>rlist</b>. The estimate assumes a
-homogeneous particle distribution, hence the drift might be slightly
+homogeneous particle distribution, hence the errors might be slightly
underestimated for multi-phase systems. For longer pair-list life-time
-(<b>nstlist</b>-1)*dt the drift is overestimated, because the interactions
+(<b>nstlist</b>-1)*dt the buffer is overestimated, because the interactions
between particles are ignored. Combined with cancellation of errors,
-the actual energy drift is usually one to two orders of magnitude smaller.
+the actual drift of the total energy is usually one to two orders of magnitude
+smaller.
Note that the generated buffer size takes into account that
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,
-use <b>verlet-buffer-drift</b>=-1 and set <b>rlist</b> manually.</dd>
+use <b>verlet-buffer-tolerance</b>=-1 and set <b>rlist</b> manually.</dd>
<dt><b>rlist: (1) [nm]</b></dt>
<dd>Cut-off distance for the short-range neighbor list.
With <b>cutoff-scheme</b>=<b>Verlet</b>, this is by default set by the
-<b>verlet-buffer-drift</b> option and the value of <b>rlist</b> is ignored.</dd>
+<b>verlet-buffer-tolerance</b> option and the value of <b>rlist</b> is ignored.</dd>
<dt><b>rlistlong: (-1) [nm]</b></dt>
<dd>Cut-off distance for the long-range neighbor list.
<A HREF="#user">userreal4</A><br>
<A HREF="#vdw">vdwtype</A><br>
<A HREF="#vdw">vdw-modifier</A><br>
-<A HREF="#nl">verlet-buffer-drift</A><br>
+<A HREF="#nl">verlet-buffer-tolerance</A><br>
<A HREF="#out">xtc-grps</A><br>
<A HREF="#out">xtc-precision</A><br>
<A HREF="#sa">zero-temp-time</A><br>
}
if (file_version >= 81)
{
- gmx_fio_do_real(fio, ir->verletbuf_drift);
+ gmx_fio_do_real(fio, ir->verletbuf_tol);
}
else
{
- ir->verletbuf_drift = 0;
+ ir->verletbuf_tol = 0;
}
gmx_fio_do_real(fio, ir->rlist);
if (file_version >= 67)
{
pr_rvec(fp, indent, "posres-comB", ir->posres_comB, DIM, TRUE);
}
- PR("verlet-buffer-drift", ir->verletbuf_drift);
+ PR("verlet-buffer-tolerance", ir->verletbuf_tol);
PR("rlist", ir->rlist);
PR("rlistlong", ir->rlistlong);
PR("nstcalclr", ir->nstcalclr);
}
void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
- const t_inputrec *ir, real drift_target,
+ const t_inputrec *ir,
const verletbuf_list_setup_t *list_setup,
int *n_nonlin_vsite,
real *rlist)
drift);
}
- if (fabs(drift) > drift_target)
+ if (fabs(drift) > ir->verletbuf_tol)
{
ib0 = ib;
}
/* 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.
- * The target is a maximum energy drift.
+ * 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.
* Returns the pair-list cut-off.
*/
void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol,
- const t_inputrec *ir, real drift_target,
+ const t_inputrec *ir,
const verletbuf_list_setup_t *list_setup,
int *n_nonlin_vsite,
real *rlist);
warning_error(wi, "rvdw should be >= 0");
}
if (ir->rlist < 0 &&
- !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_drift > 0))
+ !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
{
warning_error(wi, "rlist should be >= 0");
}
rc_max = max(ir->rvdw, ir->rcoulomb);
- if (ir->verletbuf_drift <= 0)
+ if (ir->verletbuf_tol <= 0)
{
- if (ir->verletbuf_drift == 0)
+ if (ir->verletbuf_tol == 0)
{
- warning_error(wi, "Can not have an energy drift of exactly 0");
+ warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
}
if (ir->rlist < rc_max)
{
if (ir->rlist > rc_max)
{
- warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-drift > 0. Will set rlist using verlet-buffer-drift.");
+ warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
}
if (ir->nstlist == 1)
{
if (EI_MD(ir->eI) && ir->etc == etcNO)
{
- warning_error(wi, "Temperature coupling is required for calculating rlist using the energy drift with verlet-buffer-drift > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-drift = -1.");
+ 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-drift > 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-drift = -1.");
+ 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.");
}
/* Set rlist temporarily so we can continue processing */
ir->rlist = rc_max;
/* replace the following commands with the clearer new versions*/
REPL_TYPE("unconstrained-start", "continuation");
REPL_TYPE("foreign-lambda", "fep-lambdas");
+ REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
CTYPE ("Preprocessor information: use cpp syntax.");
CTYPE ("Periodic boundary conditions: xyz, no, xy");
EETYPE("pbc", ir->ePBC, epbc_names);
EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
- CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
+ CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
CTYPE ("a value of -1 means: use rlist");
- RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
+ RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
CTYPE ("nblist cut-off");
RTYPE ("rlist", ir->rlist, 1.0);
CTYPE ("long-range cut-off for switched potentials");
gmx_bool bCompact, int nstglobalcomm, ivec ddxyz, int dd_node_order,
real rdd, real rconstr, const char *dddlb_opt, real dlb_scale,
const char *ddcsx, const char *ddcsy, const char *ddcsz,
- const char *nbpu_opt,
+ const char *nbpu_opt, int nstlist_cmdline,
gmx_large_int_t nsteps_cmdline, int nstepout, int resetstep,
int nmultisim, int repl_ex_nst, int repl_ex_nex,
int repl_ex_seed, real pforce, real cpt_period, real max_hours,
rvec posres_com; /* The COM of the posres atoms */
rvec posres_comB; /* The B-state COM of the posres atoms */
int andersen_seed; /* Random seed for Andersen thermostat (obsolete) */
- real verletbuf_drift; /* Max. drift (kJ/mol/ps/atom) for list buffer */
+ real verletbuf_tol; /* Per atom pair energy drift tolerance (kJ/mol/ps/atom) for list buffer */
real rlist; /* short range pairlist cut-off (nm) */
real rlistlong; /* long range pairlist cut-off (nm) */
int nstcalclr; /* Frequency of evaluating direct space long-range interactions */
}
#endif
-static real nbnxn_rlist_inc_nonloc_fac = 0.6;
+/* Clusters at the cut-off only increase rlist by 60% of their size */
+static real nbnxn_rlist_inc_outside_fac = 0.6;
/* Due to the cluster size the effective pair-list is longer than
* that of a simple atom pair-list. This function gives the extra distance.
*/
-real nbnxn_get_rlist_effective_inc(int cluster_size, real atom_density)
+real nbnxn_get_rlist_effective_inc(int cluster_size_j, real atom_density)
{
- return ((0.5 + nbnxn_rlist_inc_nonloc_fac)*sqr(((cluster_size) - 1.0)/(cluster_size))*pow((cluster_size)/(atom_density), 1.0/3.0));
+ int cluster_size_i;
+ real vol_inc_i, vol_inc_j;
+
+ /* We should get this from the setup, but currently it's the same for
+ * all setups, including GPUs.
+ */
+ cluster_size_i = NBNXN_CPU_CLUSTER_I_SIZE;
+
+ vol_inc_i = (cluster_size_i - 1)/atom_density;
+ vol_inc_j = (cluster_size_j - 1)/atom_density;
+
+ return nbnxn_rlist_inc_outside_fac*pow(vol_inc_i + vol_inc_j, 1.0/3.0);
}
/* Estimates the interaction volume^2 for non-local interactions */
xy_diag2 = ls[XX]*ls[XX] + ls[YY]*ls[YY] + ls[ZZ]*ls[ZZ];
/* The formulas below are a heuristic estimate of the average nsj per si*/
- r_eff_sup = rlist + nbnxn_rlist_inc_nonloc_fac*sqr((grid->na_c - 1.0)/grid->na_c)*sqrt(xy_diag2/3);
+ r_eff_sup = rlist + nbnxn_rlist_inc_outside_fac*sqr((grid->na_c - 1.0)/grid->na_c)*sqrt(xy_diag2/3);
if (!nbs->DomDec || nbs->zones->n == 1)
{
static void set_verlet_buffer(const gmx_mtop_t *mtop,
t_inputrec *ir,
matrix box,
- real verletbuf_drift,
warninp_t wi)
{
real ref_T;
{
if (ir->opts.ref_t[i] < 0)
{
- warning(wi, "Some atom groups do not use temperature coupling. This cannot be accounted for in the energy drift estimation for the Verlet buffer size. The energy drift and the Verlet buffer might be underestimated.");
+ 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.");
}
else
{
}
}
- printf("Determining Verlet buffer for an energy drift of %g kJ/mol/ps at %g K\n", verletbuf_drift, ref_T);
+ 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++)
{
/* 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, verletbuf_drift,
+ calc_verlet_buffer_size(mtop, det(box), ir,
&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, verletbuf_drift,
+ calc_verlet_buffer_size(mtop, det(box), ir,
&ls, &n_nonlin_vsite, &ir->rlist);
if (n_nonlin_vsite > 0)
{
- sprintf(warn_buf, "There are %d non-linear virtual site constructions. Their contribution to the energy drift is approximated. In most cases this does not affect the energy drift significantly.", n_nonlin_vsite);
+ sprintf(warn_buf, "There are %d non-linear virtual site constructions. Their contribution to the energy error is approximated. In most cases this does not affect the error significantly.", n_nonlin_vsite);
warning_note(wi, warn_buf);
}
if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
{
- gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-drift.", ir->rlistlong, sqrt(max_cutoff2(ir->ePBC, box)));
+ gmx_fatal(FARGS, "The pair-list cut-off (%g nm) is longer than half the shortest box vector or longer than the smallest box diagonal element (%g nm). Increase the box size or decrease nstlist or increase verlet-buffer-tolerance.", ir->rlistlong, sqrt(max_cutoff2(ir->ePBC, box)));
}
}
bGenVel ? state.v : NULL,
wi);
- if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_drift > 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)
{
- set_verlet_buffer(sys, ir, state.box, ir->verletbuf_drift, wi);
+ set_verlet_buffer(sys, ir, state.box, wi);
}
}
cmp_int(fp, "refcoord_scaling", -1, ir1->refcoord_scaling, ir2->refcoord_scaling);
cmp_rvec(fp, "inputrec->posres_com", -1, ir1->posres_com, ir2->posres_com, ftol, abstol);
cmp_rvec(fp, "inputrec->posres_comB", -1, ir1->posres_comB, ir2->posres_comB, ftol, abstol);
- cmp_real(fp, "inputrec->verletbuf_drift", -1, ir1->verletbuf_drift, ir2->verletbuf_drift, ftol, abstol);
+ cmp_real(fp, "inputrec->verletbuf_tol", -1, ir1->verletbuf_tol, ir2->verletbuf_tol, ftol, abstol);
cmp_real(fp, "inputrec->rlist", -1, ir1->rlist, ir2->rlist, ftol, abstol);
cmp_real(fp, "inputrec->rlistlong", -1, ir1->rlistlong, ir2->rlistlong, ftol, abstol);
cmp_int(fp, "inputrec->nstcalclr", -1, ir1->nstcalclr, ir2->nstcalclr);
"multiple times, e.g. \"[TT]0011[tt]\" for four ranks sharing two GPUs in this node.",
"This works within a single simulation, or a multi-simulation, with any form of MPI.",
"[PAR]",
+ "With the Verlet cut-off scheme and verlet-buffer-tolerance set,",
+ "the pair-list update interval nstlist can be chosen freely with",
+ "the option [TT]-nstlist[tt]. [TT]mdrun[tt] will then adjust",
+ "the pair-list cut-off to maintain accuracy.",
+ "By default [TT]mdrun[tt] will try to increase nstlist to improve",
+ "the performance. For CPU runs nstlist might increase to 20, for GPU",
+ "runs up till 40. But for medium to high parallelization or with",
+ "fast GPUs, a (user supplied) larger nstlist value can give much",
+ "better performance.",
+ "[PAR]",
"When using PME with separate PME nodes or with a GPU, the two major",
"compute tasks, the non-bonded force calculation and the PME calculation",
"run on different compute resources. If this load is not balanced,",
gmx_bool bReproducible = FALSE;
int npme = -1;
+ int nstlist = 0;
int nmultisim = 0;
int nstglobalcomm = -1;
int repl_ex_nst = 0;
"Global communication frequency" },
{ "-nb", FALSE, etENUM, {&nbpu_opt},
"Calculate non-bonded interactions on" },
+ { "-nstlist", FALSE, etINT, {&nstlist},
+ "Set nstlist when using a Verlet buffer tolerance (0 is guess)" },
{ "-tunepme", FALSE, etBOOL, {&bTunePME},
"Optimize PME load between PP/PME nodes or GPU/CPU" },
{ "-testverlet", FALSE, etBOOL, {&bTestVerlet},
rc = mdrunner(&hw_opt, fplog, cr, NFILE, fnm, oenv, bVerbose, bCompact,
nstglobalcomm, ddxyz, dd_node_order, rdd, rconstr,
dddlb_opt[0], dlb_scale, ddcsx, ddcsy, ddcsz,
- nbpu_opt[0],
+ nbpu_opt[0], nstlist,
nsteps, nstepout, resetstep,
nmultisim, repl_ex_nst, repl_ex_nex, repl_ex_seed,
pforce, cpt_period, max_hours, deviceOptions, Flags);
const char *ddcsy;
const char *ddcsz;
const char *nbpu_opt;
+ int nstlist_cmdline;
gmx_large_int_t nsteps_cmdline;
int nstepout;
int resetstep;
mc.ddxyz, mc.dd_node_order, mc.rdd,
mc.rconstr, mc.dddlb_opt, mc.dlb_scale,
mc.ddcsx, mc.ddcsy, mc.ddcsz,
- mc.nbpu_opt,
+ mc.nbpu_opt, mc.nstlist_cmdline,
mc.nsteps_cmdline, mc.nstepout, mc.resetstep,
mc.nmultisim, mc.repl_ex_nst, mc.repl_ex_nex, mc.repl_ex_seed, mc.pforce,
mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.Flags);
ivec ddxyz, int dd_node_order, real rdd, real rconstr,
const char *dddlb_opt, real dlb_scale,
const char *ddcsx, const char *ddcsy, const char *ddcsz,
- const char *nbpu_opt,
+ const char *nbpu_opt, int nstlist_cmdline,
gmx_large_int_t nsteps_cmdline,
int nstepout, int resetstep,
int nmultisim, int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
mda->ddcsy = ddcsy;
mda->ddcsz = ddcsz;
mda->nbpu_opt = nbpu_opt;
+ mda->nstlist_cmdline= nstlist_cmdline;
mda->nsteps_cmdline = nsteps_cmdline;
mda->nstepout = nstepout;
mda->resetstep = resetstep;
#endif /* GMX_THREAD_MPI */
-/* Environment variable for setting nstlist */
-static const char* NSTLIST_ENVVAR = "GMX_NSTLIST";
-/* Try to increase nstlist when using a GPU with nstlist less than this */
-static const int NSTLIST_GPU_ENOUGH = 20;
-/* Increase nstlist until the non-bonded cost increases more than this factor */
-static const float NBNXN_GPU_LIST_OK_FAC = 1.20;
-/* Don't increase nstlist beyond a non-bonded cost increases of this factor.
+/* We determine the extra cost of the non-bonded kernels compared to
+ * a reference nstlist value of 10 (which is the default in grompp).
+ */
+static const int nbnxn_reference_nstlist = 10;
+/* The values to try when switching */
+const int nstlist_try[] = { 20, 25, 40 };
+#define NNSTL sizeof(nstlist_try)/sizeof(nstlist_try[0])
+/* Increase nstlist until the non-bonded cost increases more than listfac_ok,
+ * but never more than listfac_max.
* A standard (protein+)water system at 300K with PME ewald_rtol=1e-5
* needs 1.28 at rcoulomb=0.9 and 1.24 at rcoulomb=1.0 to get to nstlist=40.
+ * Note that both CPU and GPU factors are conservative. Performance should
+ * not go down due to this tuning, except with a relatively slow GPU.
+ * On the other hand, at medium/high parallelization or with fast GPUs
+ * nstlist will not be increased enough to reach optimal performance.
*/
-static const float NBNXN_GPU_LIST_MAX_FAC = 1.30;
-
-/* Try to increase nstlist when running on a GPU */
+/* CPU: pair-search is about a factor 1.5 slower than the non-bonded kernel */
+static const float nbnxn_cpu_listfac_ok = 1.05;
+static const float nbnxn_cpu_listfac_max = 1.09;
+/* GPU: pair-search is a factor 1.5-3 slower than the non-bonded kernel */
+static const float nbnxn_gpu_listfac_ok = 1.20;
+static const float nbnxn_gpu_listfac_max = 1.30;
+
+/* Try to increase nstlist when using the Verlet cut-off scheme */
static void increase_nstlist(FILE *fp, t_commrec *cr,
- t_inputrec *ir, const gmx_mtop_t *mtop, matrix box)
+ t_inputrec *ir, int nstlist_cmdline,
+ const gmx_mtop_t *mtop, matrix box,
+ gmx_bool bGPU)
{
- char *env;
+ float listfac_ok, listfac_max;
int nstlist_orig, nstlist_prev;
verletbuf_list_setup_t ls;
real rlist_nstlist10, rlist_inc, rlist_ok, rlist_max;
real rlist_new, rlist_prev;
- int i;
+ int nstlist_ind = 0;
t_state state_tmp;
gmx_bool bBox, bDD, bCont;
- const char *nstl_fmt = "\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 *vbd_err = "Can not increase nstlist for GPU run because verlet-buffer-drift is not set or used";
- const char *box_err = "Can not increase nstlist for GPU run because the box is too small";
- const char *dd_err = "Can not increase nstlist for GPU run because of domain decomposition limitations";
+ 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 *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";
char buf[STRLEN];
- /* Number of + nstlist alternative values to try when switching */
- const int nstl[] = { 20, 25, 40 };
-#define NNSTL sizeof(nstl)/sizeof(nstl[0])
-
- env = getenv(NSTLIST_ENVVAR);
- if (env == NULL)
+ if (nstlist_cmdline <= 0)
{
- if (fp != NULL)
+ if (fp != NULL && bGPU && ir->nstlist < nstlist_try[0])
+ {
+ fprintf(fp, nstl_gpu, ir->nstlist);
+ }
+ nstlist_ind = 0;
+ while (nstlist_ind < NNSTL && ir->nstlist >= nstlist_try[nstlist_ind])
+ {
+ nstlist_ind++;
+ }
+ if (nstlist_ind == NNSTL)
{
- fprintf(fp, nstl_fmt, ir->nstlist);
+ /* There are no larger nstlist value to try */
+ return;
}
}
- if (ir->verletbuf_drift == 0)
+ 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");
}
- if (ir->verletbuf_drift < 0)
+ if (ir->verletbuf_tol < 0)
{
if (MASTER(cr))
{
return;
}
+ if (bGPU)
+ {
+ listfac_ok = nbnxn_gpu_listfac_ok;
+ listfac_max = nbnxn_gpu_listfac_max;
+ }
+ else
+ {
+ listfac_ok = nbnxn_cpu_listfac_ok;
+ listfac_max = nbnxn_cpu_listfac_max;
+ }
+
nstlist_orig = ir->nstlist;
- if (env != NULL)
+ if (nstlist_cmdline > 0)
{
- sprintf(buf, "Getting nstlist from environment variable GMX_NSTLIST=%s", env);
- if (MASTER(cr))
+ if (fp)
{
- fprintf(stderr, "%s\n", buf);
+ sprintf(buf, "Getting nstlist=%d from command line option",
+ nstlist_cmdline);
}
- if (fp != NULL)
- {
- fprintf(fp, "%s\n", buf);
- }
- sscanf(env, "%d", &ir->nstlist);
+ ir->nstlist = nstlist_cmdline;
}
- verletbuf_get_list_setup(TRUE, &ls);
+ verletbuf_get_list_setup(bGPU, &ls);
/* Allow rlist to make the list a given factor larger than the list
* would be with nstlist=10.
*/
nstlist_prev = ir->nstlist;
ir->nstlist = 10;
- calc_verlet_buffer_size(mtop, det(box), ir, ir->verletbuf_drift, &ls,
- NULL, &rlist_nstlist10);
+ calc_verlet_buffer_size(mtop, det(box), ir, &ls, NULL, &rlist_nstlist10);
ir->nstlist = nstlist_prev;
/* Determine the pair list size increase due to zero interactions */
- rlist_inc = nbnxn_get_rlist_effective_inc(NBNXN_GPU_CLUSTER_SIZE, mtop->natoms/det(box));
- rlist_ok = (rlist_nstlist10 + rlist_inc)*pow(NBNXN_GPU_LIST_OK_FAC, 1.0/3.0) - rlist_inc;
- rlist_max = (rlist_nstlist10 + rlist_inc)*pow(NBNXN_GPU_LIST_MAX_FAC, 1.0/3.0) - rlist_inc;
+ rlist_inc = nbnxn_get_rlist_effective_inc(ls.cluster_size_j,
+ mtop->natoms/det(box));
+ rlist_ok = (rlist_nstlist10 + rlist_inc)*pow(listfac_ok, 1.0/3.0) - rlist_inc;
+ rlist_max = (rlist_nstlist10 + rlist_inc)*pow(listfac_max, 1.0/3.0) - rlist_inc;
if (debug)
{
- fprintf(debug, "GPU nstlist tuning: rlist_inc %.3f rlist_max %.3f\n",
- rlist_inc, rlist_max);
+ fprintf(debug, "nstlist tuning: rlist_inc %.3f rlist_ok %.3f rlist_max %.3f\n",
+ rlist_inc, rlist_ok, rlist_max);
}
- i = 0;
nstlist_prev = nstlist_orig;
rlist_prev = ir->rlist;
do
{
- if (env == NULL)
+ if (nstlist_cmdline <= 0)
{
- ir->nstlist = nstl[i];
+ ir->nstlist = nstlist_try[nstlist_ind];
}
/* Set the pair-list buffer size in ir */
- calc_verlet_buffer_size(mtop, det(box), ir, ir->verletbuf_drift, &ls,
- NULL, &rlist_new);
+ calc_verlet_buffer_size(mtop, det(box), ir, &ls, NULL, &rlist_new);
/* Does rlist fit in the box? */
bBox = (sqr(rlist_new) < max_cutoff2(ir->ePBC, box));
bDD = change_dd_cutoff(cr, &state_tmp, ir, rlist_new);
}
+ if (debug)
+ {
+ fprintf(debug, "nstlist %d rlist %.3f bBox %d bDD %d\n",
+ ir->nstlist, rlist_new, bBox, bDD);
+ }
+
bCont = FALSE;
- if (env == NULL)
+ if (nstlist_cmdline <= 0)
{
if (bBox && bDD && rlist_new <= rlist_max)
{
/* Increase nstlist */
nstlist_prev = ir->nstlist;
rlist_prev = rlist_new;
- bCont = (i+1 < NNSTL && rlist_new < rlist_ok);
+ bCont = (nstlist_ind+1 < NNSTL && rlist_new < rlist_ok);
}
else
{
}
}
- i++;
+ nstlist_ind++;
}
while (bCont);
static void prepare_verlet_scheme(FILE *fplog,
t_commrec *cr,
t_inputrec *ir,
+ int nstlist_cmdline,
const gmx_mtop_t *mtop,
matrix box,
gmx_bool bUseGPU)
{
- if (ir->verletbuf_drift > 0)
+ if (ir->verletbuf_tol > 0)
{
/* 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,
- ir->verletbuf_drift, &ls,
- NULL, &rlist_new);
+ calc_verlet_buffer_size(mtop, det(box), ir, &ls, NULL, &rlist_new);
+
if (rlist_new != ir->rlist)
{
if (fplog != NULL)
}
}
- /* With GPU or emulation we should check nstlist for performance */
- if ((EI_DYNAMICS(ir->eI) &&
- bUseGPU &&
- ir->nstlist < NSTLIST_GPU_ENOUGH) ||
- getenv(NSTLIST_ENVVAR) != NULL)
+ if (nstlist_cmdline > 0 && (!EI_DYNAMICS(ir->eI) || ir->verletbuf_tol <= 0))
+ {
+ gmx_fatal(FARGS, "Can not set nstlist without %s",
+ !EI_DYNAMICS(ir->eI) ? "dynamics" : "verlet-buffer-tolerance");
+ }
+
+ if (EI_DYNAMICS(ir->eI))
{
- /* Choose a better nstlist */
- increase_nstlist(fplog, cr, ir, mtop, box);
+ /* Set or try nstlist values */
+ increase_nstlist(fplog, cr, ir, nstlist_cmdline, mtop, box, bUseGPU);
}
}
md_print_warn(NULL, fplog, "%s\n", conv_mesg);
- ir->cutoff_scheme = ecutsVERLET;
- ir->verletbuf_drift = 0.005;
+ ir->cutoff_scheme = ecutsVERLET;
+ ir->verletbuf_tol = 0.005;
if (ir->rcoulomb != ir->rvdw)
{
}
}
- /* We set the target energy drift to a small number.
+ /* We set the pair energy error tolerance to a small number.
* Note that this is only for testing. For production the user
* should think about this and set the mdp options.
*/
- ir->verletbuf_drift = 1e-4;
+ ir->verletbuf_tol = 1e-4;
}
if (inputrec2nboundeddim(ir) != 3)
verletbuf_list_setup_t ls;
verletbuf_get_list_setup(FALSE, &ls);
- calc_verlet_buffer_size(mtop, box_vol, ir, ir->verletbuf_drift, &ls,
- NULL, &ir->rlist);
+ calc_verlet_buffer_size(mtop, box_vol, ir, &ls, NULL, &ir->rlist);
}
else
{
- ir->verletbuf_drift = -1;
+ ir->verletbuf_tol = -1;
ir->rlist = 1.05*max(ir->rvdw, ir->rcoulomb);
}
ivec ddxyz, int dd_node_order, real rdd, real rconstr,
const char *dddlb_opt, real dlb_scale,
const char *ddcsx, const char *ddcsy, const char *ddcsz,
- const char *nbpu_opt,
+ const char *nbpu_opt, int nstlist_cmdline,
gmx_large_int_t nsteps_cmdline, int nstepout, int resetstep,
int gmx_unused nmultisim, int repl_ex_nst, int repl_ex_nex,
int repl_ex_seed, real pforce, real cpt_period, real max_hours,
getenv("GMX_EMULATE_GPU") != NULL);
prepare_verlet_scheme(fplog, cr,
- inputrec, mtop, state->box,
+ inputrec, nstlist_cmdline, mtop, state->box,
bUseGPU);
}
- else if (hwinfo->gpu_info.ncuda_dev_compatible > 0)
+ else
{
- md_print_warn(cr, fplog,
- "NOTE: GPU(s) found, but the current simulation can not use GPUs\n"
- " To use a GPU, set the mdp option: cutoff-scheme = Verlet\n"
- " (for quick performance testing you can use the -testverlet option)\n");
+ if (nstlist_cmdline > 0)
+ {
+ gmx_fatal(FARGS, "Can not set nstlist with the group cut-off scheme");
+ }
+
+ if (hwinfo->gpu_info.ncuda_dev_compatible > 0)
+ {
+ md_print_warn(cr, fplog,
+ "NOTE: GPU(s) found, but the current simulation can not use GPUs\n"
+ " To use a GPU, set the mdp option: cutoff-scheme = Verlet\n"
+ " (for quick performance testing you can use the -testverlet option)\n");
+ }
if (bForceUseGPU)
{
gmx_fatal(FARGS, "GPU requested, but can't be used without cutoff-scheme=Verlet");
}
- }
+
#ifdef GMX_TARGET_BGQ
- else
- {
md_print_warn(cr, fplog,
"NOTE: There is no SIMD implementation of the group scheme kernels on\n"
" BlueGene/Q. You will observe better performance from using the\n"
" Verlet cut-off scheme.\n");
- }
#endif
+ }
}
/* Check for externally set OpenMP affinity and turn off internal
oenv, bVerbose, bCompact, nstglobalcomm,
ddxyz, dd_node_order, rdd, rconstr,
dddlb_opt, dlb_scale, ddcsx, ddcsy, ddcsz,
- nbpu_opt,
+ nbpu_opt, nstlist_cmdline,
nsteps_cmdline, nstepout, resetstep, nmultisim,
repl_ex_nst, repl_ex_nex, repl_ex_seed, pforce,
cpt_period, max_hours, deviceOptions,
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff