Automated Verlet rlist for NVE simulations

[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / calc_verletbuf.c

diff --git a/src/gromacs/gmxpreprocess/calc_verletbuf.c b/src/gromacs/gmxpreprocess/calc_verletbuf.c
index f16efb0646124932ec15fc23c91348b7142e110d..7c28034583338e39f53e9626b15ce7f827b32920 100644 (file)
--- 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;