From b9fc004924a97f3b1b762515edbfcb7cdaecd243 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Thu, 9 Jan 2014 13:58:12 +0100 Subject: [PATCH] Automated Verlet rlist for NVE simulations NVE simulations with the Verlet scheme will now automatically have the list buffer set using verlet-buffer-tolerance and the initial temperature, unless this is zero in which case a 10% buffer is used. A warning is added for potentially large drift over the simulation. This change makes empty mdp files valid again. Change-Id: Iae2d3d8ba212dd9e57783d9fb3d84f7cfc1efdd0 --- share/html/online/mdp_opt.html | 5 +- src/gromacs/gmxpreprocess/calc_verletbuf.c | 30 +++- src/gromacs/gmxpreprocess/calc_verletbuf.h | 13 ++ src/gromacs/gmxpreprocess/grompp.c | 155 +++++++++++++++++---- src/gromacs/gmxpreprocess/readir.c | 8 +- src/programs/mdrun/runner.c | 39 +++++- 6 files changed, 205 insertions(+), 45 deletions(-) diff --git a/share/html/online/mdp_opt.html b/share/html/online/mdp_opt.html index a759d4df72..f8b4319974 100644 --- a/share/html/online/mdp_opt.html +++ b/share/html/online/mdp_opt.html @@ -531,7 +531,10 @@ 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, +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 verlet-buffer-tolerance=-1 and set rlist manually.
rlist: (1) [nm]
diff --git a/src/gromacs/gmxpreprocess/calc_verletbuf.c b/src/gromacs/gmxpreprocess/calc_verletbuf.c index f16efb0646..7c28034583 100644 --- a/src/gromacs/gmxpreprocess/calc_verletbuf.c +++ b/src/gromacs/gmxpreprocess/calc_verletbuf.c @@ -1,7 +1,7 @@ /* * 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. @@ -712,6 +712,7 @@ static real surface_frac(int cluster_size, real particle_distance, real rlist) 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) @@ -732,6 +733,29 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, 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; @@ -868,7 +892,7 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, * 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 */ @@ -899,7 +923,7 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, } 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; diff --git a/src/gromacs/gmxpreprocess/calc_verletbuf.h b/src/gromacs/gmxpreprocess/calc_verletbuf.h index e2bbb48ebe..65ef492350 100644 --- a/src/gromacs/gmxpreprocess/calc_verletbuf.h +++ b/src/gromacs/gmxpreprocess/calc_verletbuf.h @@ -49,6 +49,17 @@ typedef struct } 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. @@ -60,6 +71,7 @@ void verletbuf_get_list_setup(gmx_bool bGPU, /* 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. @@ -67,6 +79,7 @@ void verletbuf_get_list_setup(gmx_bool bGPU, */ 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); diff --git a/src/gromacs/gmxpreprocess/grompp.c b/src/gromacs/gmxpreprocess/grompp.c index 7d867b4b5f..d604b1fb29 100644 --- a/src/gromacs/gmxpreprocess/grompp.c +++ b/src/gromacs/gmxpreprocess/grompp.c @@ -1254,52 +1254,89 @@ static void check_gbsa_params(gpp_atomtype_t atype) } -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) @@ -1773,11 +1810,71 @@ int gmx_grompp(int argc, char *argv[]) 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); + } } } diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c index a5e2c307fc..a0c9f806ec 100644 --- a/src/gromacs/gmxpreprocess/readir.c +++ b/src/gromacs/gmxpreprocess/readir.c @@ -62,6 +62,7 @@ #include "mtop_util.h" #include "chargegroup.h" #include "inputrec.h" +#include "calc_verletbuf.h" #define MAXPTR 254 #define NOGID 255 @@ -380,11 +381,6 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts, { 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."); @@ -395,7 +391,7 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts, 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; } } } diff --git a/src/programs/mdrun/runner.c b/src/programs/mdrun/runner.c index fec0b0c096..2a29f9cda3 100644 --- a/src/programs/mdrun/runner.c +++ b/src/programs/mdrun/runner.c @@ -508,6 +508,7 @@ static void increase_nstlist(FILE *fp, t_commrec *cr, 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"; @@ -531,6 +532,20 @@ static void increase_nstlist(FILE *fp, t_commrec *cr, } } + 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"); @@ -579,7 +594,8 @@ static void increase_nstlist(FILE *fp, t_commrec *cr, */ 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 */ @@ -603,7 +619,7 @@ static void increase_nstlist(FILE *fp, t_commrec *cr, } /* 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)); @@ -685,7 +701,8 @@ static void prepare_verlet_scheme(FILE *fplog, 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; @@ -697,7 +714,7 @@ static void prepare_verlet_scheme(FILE *fplog, */ 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) { @@ -790,12 +807,22 @@ static void convert_to_verlet_scheme(FILE *fplog, 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); -- 2.22.0