Fixed minor issue with NPT conserved energy
authorMark Abraham <mark.j.abraham@gmail.com>
Sat, 20 Sep 2014 10:48:41 +0000 (12:48 +0200)
committerMark Abraham <mark.j.abraham@gmail.com>
Thu, 9 Oct 2014 14:28:31 +0000 (16:28 +0200)
The value of the number of degrees of freedom in a group is a real,
but it was being truncated to an integer before promoting it back to a
real. This is wrong, whether or not the number can reasonably be
non-integral. This only affected the conserved energy quantity,
not the kinetic or total energy used for coupling and output.

Change-Id: I5fcb6428dc0e333cabd79262cc2a9cffa5fec03a

src/gromacs/mdlib/coupling.c

index 41806dc44b4338b324cf84d733cc81e2527d4505..b156c61fad22009667945749571563f0af1b46bd 100644 (file)
@@ -1236,7 +1236,8 @@ int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool b
 
 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
 {
-    int     i, j, nd, ndj, bmass, qmass, ngtcall;
+    int     i, j, bmass, qmass, ngtcall;
+    real    nd, ndj;
     real    ener_npt, reft, eta, kT, tau;
     double *ivxi, *ixi;
     double *iQinv;
@@ -1319,7 +1320,7 @@ real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
             reft = max(ir->opts.ref_t[i], 0);
             kT   = BOLTZ * reft;
 
-            if (nd > 0)
+            if (nd > 0.0)
             {
                 if (IR_NVT_TROTTER(ir))
                 {
@@ -1336,7 +1337,7 @@ real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
                             }
                             else
                             {
-                                ndj = 1;
+                                ndj = 1.0;
                             }
                             ener_npt += ndj*ixi[j]*kT;
                         }