Merge branch 'release-4-6' into master

[alexxy/gromacs.git] / src / gromacs / mdlib / shakef.c

diff --git a/src/gromacs/mdlib/shakef.c b/src/gromacs/mdlib/shakef.c
index 57a6a8d38691bbaafa6257e1a79098f41c7a6b21..659fc32889790448ce61da6406afa4850402dfa0 100644 (file)
--- a/src/gromacs/mdlib/shakef.c
+++ b/src/gromacs/mdlib/shakef.c
@@ -436,13 +436,28 @@ gmx_bool bshakef(FILE *log, gmx_shakedata_t shaked,
     {
         if (ir->efep != efepNO)
         {
+            real bondA,bondB;
             dt_2 = 1/sqr(ir->delta_t);
             dvdl = 0;
             for (i = 0; i < ncons; i++)
             {
                 type  = idef->il[F_CONSTR].iatoms[3*i];
-                dvdl += lagr[i]*dt_2*
-                    (idef->iparams[type].constr.dB-idef->iparams[type].constr.dA);
+
+                /* dh/dl contribution from constraint force is  dh/dr (constraint force) dot dr/dl */
+                /* constraint force is -\sum_i lagr_i* d(constraint)/dr, with constrant = r^2-d^2  */
+                /* constraint force is -\sum_i lagr_i* 2 r  */
+                /* so dh/dl = -\sum_i lagr_i* 2 r * dr/dl */
+                /* However, by comparison with lincs and with
+                   comparison with a full thermodynamics cycle (see
+                   redmine issue #1255), this is off by a factor of
+                   two -- the 2r should apparently just be r.  Further
+                   investigation should be done at some point to
+                   understand why and see if there is something deeper
+                   we are missing */
+
+                bondA = idef->iparams[type].constr.dA;
+                bondB = idef->iparams[type].constr.dB;
+                dvdl += lagr[i] * dt_2 * ((1.0-lambda)*bondA + lambda*bondB) * (bondB-bondA);
             }
             *dvdlambda += dvdl;
         }