Fixed dvdlambda for SHAKE + FE

Garbage values have been accumulated to dvdlambda whenever SHAKE was
used.

Refactoring in 513813b514 removed the lam variable. Previously, it was
used as a temporary iteration variable over SHAKE blocks, but that
refactoring to use the main variable changed the logic of the
following loop.

Fixes #2434

Change-Id: Ia9642f50358ab31ec98f8317f7a1dcda8622a9f1

diff --git a/docs/release-notes/2018/2018.1.rst b/docs/release-notes/2018/2018.1.rst
index 094a8cd17e82cd2a83412e7431ba8622a50fdbaf..13d0c7ed973b1b8c2161c860f2fb26266a68a97b 100644 (file)
--- a/docs/release-notes/2018/2018.1.rst
+++ b/docs/release-notes/2018/2018.1.rst
@@ -30,6 +30,13 @@ system.
 
 :issue:`2381`
 
+Fixed FEP calculations with SHAKE
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+All SHAKE + FEP calculations accumulated wrong values to dH/dl output,
+but in some cases the result will look the same.
+
+:issue:`2434`
+
 Fixed handling of mdp ``define`` statement assigning preprocessor values
 """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 Now .mdp files can configure the topology with values, as originally

diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h
index fb3c8e8e4f68e465850dd7124a32ae5098982e02..fe6cc79d965fb91a0489459e40341a4a382e7991 100644 (file)
--- a/src/gromacs/mdlib/constr.h
+++ b/src/gromacs/mdlib/constr.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018, 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.
@@ -96,7 +96,7 @@ gmx_bool bshakef(FILE           *log,          /* Log file                    */
                  rvec            x_s[],        /* Coords before update         */
                  rvec            prime[],      /* Output coords                */
                  t_nrnb         *nrnb,         /* Performance measure          */
-                 real           *lagr,         /* The Lagrange multipliers     */
+                 real * const    lagr,         /* The Lagrange multipliers     */
                  real            lambda,       /* FEP lambda                   */
                  real           *dvdlambda,    /* FEP force                    */
                  real            invdt,        /* 1/delta_t                    */

diff --git a/src/gromacs/mdlib/shakef.cpp b/src/gromacs/mdlib/shakef.cpp
index 19e3befc7b44fbe3494e9429b5b5411e5d321639..61f2dad9fed939225d6676918472dfc11284fa1e 100644 (file)
--- a/src/gromacs/mdlib/shakef.cpp
+++ b/src/gromacs/mdlib/shakef.cpp
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017,2018, 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.
@@ -393,12 +393,12 @@ static void check_cons(FILE *log, int nc, rvec x[], rvec prime[], rvec v[],
 gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked,
                  real invmass[], int nblocks, int sblock[],
                  t_idef *idef, t_inputrec *ir, rvec x_s[], rvec prime[],
-                 t_nrnb *nrnb, real *scaled_lagrange_multiplier, real lambda, real *dvdlambda,
+                 t_nrnb *nrnb, real * const scaled_lagrange_multiplier, real lambda, real *dvdlambda,
                  real invdt, rvec *v, gmx_bool bCalcVir, tensor vir_r_m_dr,
                  gmx_bool bDumpOnError, int econq)
 {
     t_iatom *iatoms;
-    real     dt_2, dvdl;
+    real    *lam, dt_2, dvdl;
     int      i, n0, ncon, blen, type, ll;
     int      tnit = 0, trij = 0;
 
@@ -413,7 +413,11 @@ gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked,
         scaled_lagrange_multiplier[ll] = 0;
     }
 
+    // TODO Rewrite this block so that it is obvious that i, iatoms
+    // and lam are all iteration variables. Is this easier if the
+    // sblock data structure is organized differently?
     iatoms = &(idef->il[F_CONSTR].iatoms[sblock[0]]);
+    lam    = scaled_lagrange_multiplier;
     for (i = 0; (i < nblocks); )
     {
         blen  = (sblock[i+1]-sblock[i]);
@@ -440,7 +444,7 @@ gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked,
         tnit                       += n0*blen;
         trij                       += blen;
         iatoms                     += 3*blen; /* Increment pointer! */
-        scaled_lagrange_multiplier += blen;
+        lam                        += blen;
         i++;
     }
     /* only for position part? */