fixed the reported temperature with Brownian dynamics

Reduced the (meaningless) masses with BD by a factor 2.
Since v with BD is defined as dx/dt, the reported Ekin and T are
now correct (and if there are too high, dt is too large).
Fixes #797

Change-Id: Iee8613ff9637deb58600ad03289f7bcce7f32df2

diff --git a/src/kernel/md.c b/src/kernel/md.c
index 1f6004f1db1396f4e26dad8999356d50d79efbb9..1022231573f91902ff86939130ecb3d827479fdb 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -162,7 +162,7 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
     gmx_bool        bAppend;
     gmx_bool        bResetCountersHalfMaxH=FALSE;
     gmx_bool        bVV,bIterations,bFirstIterate,bTemp,bPres,bTrotter;
-    real        temp0,mu_aver=0,dvdl;
+    real        mu_aver=0,dvdl;
     int         a0,a1,gnx=0,ii;
     atom_id     *grpindex=NULL;
     char        *grpname;
@@ -500,7 +500,6 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
         enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step,
                                      and there is no previous step */
     }
-    temp0 = enerd->term[F_TEMP];
     
     /* if using an iterative algorithm, we need to create a working directory for the state. */
     if (bIterations) 
@@ -528,7 +527,10 @@ double do_md(FILE *fplog,t_commrec *cr,int nfile,const t_filenm fnm[],
                     "RMS relative constraint deviation after constraining: %.2e\n",
                     constr_rmsd(constr,FALSE));
         }
-        fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
+        if (EI_STATE_VELOCITY(ir->eI))
+        {
+            fprintf(fplog,"Initial temperature: %g K\n",enerd->term[F_TEMP]);
+        }
         if (bRerunMD)
         {
             fprintf(stderr,"starting md rerun '%s', reading coordinates from"

diff --git a/src/mdlib/mdatom.c b/src/mdlib/mdatom.c
index 0db216b5c099234f241024af29f21037192a2be6..a3ad311e223b88749980ed129d3d22ccaff22768 100644 (file)
--- a/src/mdlib/mdatom.c
+++ b/src/mdlib/mdatom.c
@@ -1,4 +1,5 @@
-/*
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ *
  * 
  *                This source code is part of
  * 
@@ -196,25 +197,42 @@ void atoms2md(gmx_mtop_t *mtop,t_inputrec *ir,
     if (md->cFREEZE) {
       md->cFREEZE[i] = ggrpnr(groups,egcFREEZE,ag);
     }
-    if (EI_ENERGY_MINIMIZATION(ir->eI)) {
-      mA = 1.0;
-      mB = 1.0;
-    } else if (ir->eI == eiBD) {
-      /* Make the mass proportional to the friction coefficient for BD.
-       * This is necessary for the constraint algorithms.
-       */
-      if (ir->bd_fric) {
-       mA = ir->bd_fric*ir->delta_t;
-       mB = ir->bd_fric*ir->delta_t;
-      } else {
-       fac = ir->delta_t/opts->tau_t[md->cTC ? groups->grpnr[egcTC][ag] : 0];
-       mA = atom->m*fac;
-       mB = atom->mB*fac;
-      }
-    } else {
-      mA = atom->m;
-      mB = atom->mB;
-    }
+        if (EI_ENERGY_MINIMIZATION(ir->eI))
+        {
+            /* Displacement is proportional to F, masses used for constraints */
+            mA = 1.0;
+            mB = 1.0;
+        }
+        else if (ir->eI == eiBD)
+        {
+            /* With BD the physical masses are irrelevant.
+             * To keep the code simple we use most of the normal MD code path
+             * for BD. Thus for constraining the masses should be proportional
+             * to the friction coefficient. We set the absolute value such that
+             * m/2<(dx/dt)^2> = m/2*2kT/fric*dt = kT/2 => m=fric*dt/2
+             * Then if we set the (meaningless) velocity to v=dx/dt, we get the
+             * correct kinetic energy and temperature using the usual code path.
+             * Thus with BD v*dt will give the displacement and the reported
+             * temperature can signal bad integration (too large time step).
+             */
+            if (ir->bd_fric > 0)
+            {
+                mA = 0.5*ir->bd_fric*ir->delta_t;
+                mB = 0.5*ir->bd_fric*ir->delta_t;
+            }
+            else
+            {
+                /* The friction coefficient is mass/tau_t */
+                fac = ir->delta_t/opts->tau_t[md->cTC ? groups->grpnr[egcTC][ag] : 0];
+                mA = 0.5*atom->m*fac;
+                mB = 0.5*atom->mB*fac;
+            }
+        }
+        else
+        {
+            mA = atom->m;
+            mB = atom->mB;
+        }
     if (md->nMassPerturbed) {
       md->massA[i]     = mA;
       md->massB[i]     = mB;

diff --git a/src/mdlib/update.c b/src/mdlib/update.c
index b4228ea5552a0f58d4b28886023dcca9fa54696f..96ed94565a556604823129d7a1447d42a645ceac 100644 (file)
--- a/src/mdlib/update.c
+++ b/src/mdlib/update.c
@@ -749,9 +749,9 @@ static void do_update_bd(int start,int nrend,double dt,
                 } 
                 else 
                 {
-                    /* NOTE: invmass = 1/(mass*friction_constant*dt) */
-                    vn = invmass[n]*f[n][d]*dt 
-                        + sqrt(invmass[n])*rf[gt]*gmx_rng_gaussian_table(gaussrand);
+                    /* NOTE: invmass = 2/(mass*friction_constant*dt) */
+                    vn = 0.5*invmass[n]*f[n][d]*dt 
+                        + sqrt(0.5*invmass[n])*rf[gt]*gmx_rng_gaussian_table(gaussrand);
                 }
 
                 v[n][d]      = vn;