Add missing mass term to deltaH

The mass contribution was missing from foreign Hamiltonian difference
output from mdrun used for BAR.

Refs #2703

Change-Id: I415c835e9e1a8473243041be004013af7141a1cc

diff --git a/docs/release-notes/2018/2018.4.rst b/docs/release-notes/2018/2018.4.rst
index 189d6608395237bc7296753fb2f311ce386ca415..e754d0b84a385719f31f6a6dda3a3aedc976b9f2 100644 (file)
--- a/docs/release-notes/2018/2018.4.rst
+++ b/docs/release-notes/2018/2018.4.rst
@@ -29,6 +29,14 @@ such as BAR.
 
 :issue: `2703`
 
+Add mass contribution to foreign Hamiltonian differences
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+For free energy calculations with perturbed masses, the kinetic energy
+contribution was missing from the foreign Hamiltonian values.
+
+:issue: `2703`
+
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp
index 9d93859af17c853439e26e146324afb761f1a228..c939c4e8a7f13553ec52fde621a8812dc35727b7 100644 (file)
--- a/src/gromacs/mdlib/force.cpp
+++ b/src/gromacs/mdlib/force.cpp
@@ -724,7 +724,6 @@ void sum_epot(gmx_grppairener_t *grpp, real *epot)
 void sum_dhdl(gmx_enerdata_t *enerd, gmx::ArrayRef<const real> lambda, t_lambda *fepvals)
 {
     int    index;
-    double dlam;
 
     enerd->dvdl_lin[efptVDW] += enerd->term[F_DVDL_VDW];  /* include dispersion correction */
     enerd->term[F_DVDL]       = 0.0;
@@ -792,11 +791,14 @@ void sum_dhdl(gmx_enerdata_t *enerd, gmx::ArrayRef<const real> lambda, t_lambda
            current lambda, because the contributions to the current
            lambda are automatically zeroed */
 
+        double &enerpart_lambda = enerd->enerpart_lambda[i + 1];
+
         for (size_t j = 0; j < lambda.size(); j++)
         {
             /* Note that this loop is over all dhdl components, not just the separated ones */
-            dlam = (fepvals->all_lambda[j][i] - lambda[j]);
-            enerd->enerpart_lambda[i+1] += dlam*enerd->dvdl_lin[j];
+            const double dlam  = fepvals->all_lambda[j][i] - lambda[j];
+
+            enerpart_lambda   += dlam*enerd->dvdl_lin[j];
 
             /* Constraints can not be evaluated at foreign lambdas, so we add
              * a linear extrapolation. This is an approximation, but usually
@@ -805,14 +807,19 @@ void sum_dhdl(gmx_enerdata_t *enerd, gmx::ArrayRef<const real> lambda, t_lambda
             if ((j == efptBONDED && fepvals->separate_dvdl[efptBONDED]) ||
                 (j == efptFEP && !fepvals->separate_dvdl[efptBONDED]))
             {
-                enerd->enerpart_lambda[i + 1] += dlam*enerd->term[F_DVDL_CONSTR];
+                enerpart_lambda += dlam*enerd->term[F_DVDL_CONSTR];
+            }
+
+            if (j == efptMASS)
+            {
+                enerpart_lambda += dlam*enerd->term[F_DKDL];
             }
 
             if (debug)
             {
                 fprintf(debug, "enerdiff lam %g: (%15s), non-linear %f linear %f*%f\n",
                         fepvals->all_lambda[j][i], efpt_names[j],
-                        (enerd->enerpart_lambda[i+1] - enerd->enerpart_lambda[0]),
+                        enerpart_lambda - enerd->enerpart_lambda[0],
                         dlam, enerd->dvdl_lin[j]);
             }
         }