Improved the accuracy of the sd1 integrator

There was a small error in one of the coefficients of the sd1 integrator.
Correcting this makes sd1 as accurate as sd, unless constraints are used.

Change-Id: I2f71a80cc0487d92fd151c6ea4f243feb228626c

diff --git a/share/html/online/mdp_opt.html b/share/html/online/mdp_opt.html
index 803dff950ae09195b3c5e2ec66168325c6982541..8bb822db7d636f7046a94dd2efd6fe0269af9107 100644 (file)
--- a/share/html/online/mdp_opt.html
+++ b/share/html/online/mdp_opt.html
@@ -152,11 +152,10 @@ a Berendsen thermostat with the same <b>tau-t</b>.</dd>
 <dt><b>sd1</b></dt>
 <dd> An efficient leap-frog stochastic dynamics integrator.
 This integrator is equivalent to <b>sd</b>, except that it requires
-only one Gaussian random number and one constraint step.
-This integrator is less accurate.
-For water the relative error in the temperature with this integrator
-is 0.5 <b>delta-t</b>/<b>tau-t</b>, but for other systems it can be higher.
-Use with care and check the temperature.</dd>
+only one Gaussian random number and one constraint step and is therefore
+significantly faster. Without constraints the accuracy is the same as <b>sd</b>.
+With constraints the accuracy is significantly reduced, so then <b>sd</b>
+will often be preferred.</dd>
 <dt><b>bd</b></dt>
 <dd>An Euler integrator for Brownian or position Langevin dynamics, the
 velocity is the force divided by a friction coefficient 

diff --git a/src/kernel/readir.c b/src/kernel/readir.c
index 1adbdc6c626de5500111587d9246fcc2aecc99ab..d5ea1206f331c307555883ddd8958eceab45cea7 100644 (file)
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@ -2398,18 +2398,15 @@ void triple_check(const char *mdparin,t_inputrec *ir,gmx_mtop_t *sys,
            eel_names[eelGRF]);
     CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
   }
-    
-  if (ir->eI == eiSD1) {
-    gdt_max = 0;
-    for(i=0; (i<ir->opts.ngtc); i++)
-      gdt_max = max(gdt_max,ir->delta_t/ir->opts.tau_t[i]);
-    if (0.5*gdt_max > 0.0015) {
-      sprintf(warn_buf,"The relative error with integrator %s is 0.5*delta-t/tau-t = %g, you might want to switch to integrator %s\n",
-             ei_names[ir->eI],0.5*gdt_max,ei_names[eiSD2]);
-      warning_note(wi,warn_buf);
-    }
-  }
 
+    if (ir->eI == eiSD1 &&
+        (gmx_mtop_ftype_count(sys,F_CONSTR) > 0 ||
+         gmx_mtop_ftype_count(sys,F_SETTLE) > 0))
+    {
+        sprintf(warn_buf,"With constraints integrator %s is less accurate, consider using %s instead",ei_names[ir->eI],ei_names[eiSD2]);
+        warning_note(wi,warn_buf);
+    }
+    
   bAcc = FALSE;
   for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
     for(m=0; (m<DIM); m++) {

diff --git a/src/mdlib/update.c b/src/mdlib/update.c
index 861b48766d3be16fbd7cf5b5b28f95ab21a62d7c..2e86ae7f09b5a2fa93078e3be7385a8d2b31618f 100644 (file)
--- a/src/mdlib/update.c
+++ b/src/mdlib/update.c
@@ -550,7 +550,7 @@ static void do_update_sd1(gmx_stochd_t *sd,
   {
       kT = BOLTZ*ref_t[n];
       /* The mass is encounted for later, since this differs per atom */
-      sig[n].V  = sqrt(2*kT*(1 - sdc[n].em));
+      sig[n].V  = sqrt(kT*(1 - sdc[n].em*sdc[n].em));
   }
   
   for(n=start; n<start+homenr; n++)