Avoid energy drift with verlet-buffer-tolerance

Added a warning with weak coupling and SD when tau_t and
verlet-buffer-tolerance as so large that the temperature might
be off by more than 0.2%.
Turned off nstlist tuning in mdrun with nstlist=1.
Also fixed Andersen checks in mdrun that were useless, since they
were done before temperature coupling groups were set.

Change-Id: I41bd316449c9e03e21d33dd82ebeab034d019fc6

diff --git a/share/html/online/mdp_opt.html b/share/html/online/mdp_opt.html
index 64736bad2e53d41aaa278fc38b808188f9942c2a..7ab2d5df163d2070f0e8ea1216dfa353720df7a8 100644 (file)
--- a/share/html/online/mdp_opt.html
+++ b/share/html/online/mdp_opt.html
@@ -403,7 +403,7 @@ the neighbor list is made only once.
 With energy minimization the neighborlist will be updated for every
 energy evaluation when <b>nstlist</b><tt>&gt;0</tt>.
 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.
+<b>nstlist</b> is actually a minimum value and <tt>mdrun</tt> might increase it, unless it is set to 1.
 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

diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c
index 521a6e3268d0faf49b30d22a78d389bc6b0d8c51..2d206d01dfa03044cfb03d33e9084c932f186427 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.c
+++ b/src/gromacs/gmxpreprocess/readir.c
@@ -235,7 +235,10 @@ static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
 
 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
               warninp_t wi)
-/* Check internal consistency */
+/* Check internal consistency.
+ * NOTE: index groups are not set here yet, don't check things
+ * like temperature coupling group options here, but in triple_check
+ */
 {
     /* Strange macro: first one fills the err_buf, and then one can check
      * the condition, which will print the message and increase the error
@@ -955,14 +958,6 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
         CHECK(!(EI_VV(ir->eI)));
 
-        for (i = 0; i < ir->opts.ngtc; i++)
-        {
-            sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
-            CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
-            sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
-                    i, ir->opts.tau_t[i]);
-            CHECK(ir->opts.tau_t[i] < 0);
-        }
         if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
         {
             sprintf(warn_buf, "Center of mass removal not necessary for %s.  All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.", etcoupl_names[ir->etc]);
@@ -971,14 +966,8 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
 
         sprintf(err_buf, "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step", ir->nstcomm, etcoupl_names[ir->etc]);
         CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
-
-        for (i = 0; i < ir->opts.ngtc; i++)
-        {
-            int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
-            sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
-            CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
-        }
     }
+
     if (ir->etc == etcBERENDSEN)
     {
         sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
@@ -3950,7 +3939,7 @@ check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
                   warninp_t wi)
 {
-    char                      err_buf[256];
+    char                      err_buf[STRLEN];
     int                       i, m, c, nmol, npct;
     gmx_bool                  bCharge, bAcc;
     real                      gdt_max, *mgrp, mt;
@@ -3963,6 +3952,71 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
 
     set_warning_line(wi, mdparin, -1);
 
+    if (ir->cutoff_scheme == ecutsVERLET &&
+        ir->verletbuf_tol > 0 &&
+        ir->nstlist > 1 &&
+        ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
+         (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
+    {
+        /* Check if a too small Verlet buffer might potentially
+         * cause more drift than the thermostat can couple off.
+         */
+        /* Temperature error fraction for warning and suggestion */
+        const real T_error_warn    = 0.002;
+        const real T_error_suggest = 0.001;
+        /* For safety: 2 DOF per atom (typical with constraints) */
+        const real nrdf_at         = 2;
+        real       T, tau, max_T_error;
+        int        i;
+
+        T   = 0;
+        tau = 0;
+        for (i = 0; i < ir->opts.ngtc; i++)
+        {
+            T   = max(T, ir->opts.ref_t[i]);
+            tau = max(tau, ir->opts.tau_t[i]);
+        }
+        if (T > 0)
+        {
+            /* This is a worst case estimate of the temperature error,
+             * assuming perfect buffer estimation and no cancelation
+             * of errors. The factor 0.5 is because energy distributes
+             * equally over Ekin and Epot.
+             */
+            max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
+            if (max_T_error > T_error_warn)
+            {
+                sprintf(warn_buf, "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to %.0e or decrease tau_t.",
+                        ir->verletbuf_tol, T, tau,
+                        100*max_T_error,
+                        100*T_error_suggest,
+                        ir->verletbuf_tol*T_error_suggest/max_T_error);
+                warning(wi, warn_buf);
+            }
+        }
+    }
+
+    if (ETC_ANDERSEN(ir->etc))
+    {
+        int i;
+
+        for (i = 0; i < ir->opts.ngtc; i++)
+        {
+            sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
+            CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
+            sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
+                    i, ir->opts.tau_t[i]);
+            CHECK(ir->opts.tau_t[i] < 0);
+        }
+
+        for (i = 0; i < ir->opts.ngtc; i++)
+        {
+            int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
+            sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
+            CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
+        }
+    }
+
     if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
         ir->comm_mode == ecmNO &&
         !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&

diff --git a/src/programs/mdrun/runner.c b/src/programs/mdrun/runner.c
index a0d7b7f8a4eb3f811424b7b5416f01511a03fa71..bf6fd312e787196968de98aefe4dd1b863af4a1d 100644 (file)
--- a/src/programs/mdrun/runner.c
+++ b/src/programs/mdrun/runner.c
@@ -523,6 +523,14 @@ static void increase_nstlist(FILE *fp, t_commrec *cr,
 
     if (nstlist_cmdline <= 0)
     {
+        if (ir->nstlist == 1)
+        {
+            /* The user probably set nstlist=1 for a reason,
+             * don't mess with the settings.
+             */
+            return;
+        }
+
         if (fp != NULL && bGPU && ir->nstlist < nstlist_try[0])
         {
             fprintf(fp, nstl_gpu, ir->nstlist);