Fixing the constraint dhdl term in shake

Addresses and corrects new issues with #1255,
with additional comments explaining what else
may need to be done in the future to make sure
that the theory of the changes are completely
correctly founded.  See the redmine entry #1255
to see the numerical validation of the current choices.
Also adds possible missing constraint dHdl contribution.
Fixes #1268

Change-Id: Ibe3980c2e265cced49b03b966e02143c5cf9f07c

diff --git a/src/kernel/md.c b/src/kernel/md.c
index 24eb4d6cf69f512d1f5f98a753032a18d68e8f95..0ab2b9313fe62312153b4cdbc654c5592a1e16c2 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -1802,7 +1802,28 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 {
                     fprintf(fplog, sepdvdlformat, "Constraint dV/dl", 0.0, dvdl_constr);
                 }
-                enerd->term[F_DVDL_CONSTR] += dvdl_constr;
+                if (bVV)
+                {
+                    /* this factor or 2 correction is necessary
+                       because half of the constraint force is removed
+                       in the vv step, so we have to double it.  See
+                       the Redmine issue #1255.  It is not yet clear
+                       if the factor of 2 is exact, or just a very
+                       good approximation, and this will be
+                       investigated.  The next step is to see if this
+                       can be done adding a dhdl contribution from the
+                       rattle step, but this is somewhat more
+                       complicated with the current code. Will be
+                       investigated, hopefully for 4.6.3. However,
+                       this current solution is much better than
+                       having it completely wrong.
+                    */
+                    enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
+                }
+                else
+                {
+                    enerd->term[F_DVDL_CONSTR] += dvdl_constr;
+                }
             }
             else if (graph)
             {
@@ -1882,8 +1903,6 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             bFirstIterate = FALSE;
         }
 
-        /* only add constraint dvdl after constraints */
-        enerd->term[F_DVDL_CONSTR] += dvdl_constr;
         if (!bVV || bRerunMD)
         {
             /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */

diff --git a/src/mdlib/force.c b/src/mdlib/force.c
index 4666743c00830bc312ad98d9b289853bbad8e132..75da6bd137d3f9847956ba7763602eb5b0395975 100644 (file)
--- a/src/mdlib/force.c
+++ b/src/mdlib/force.c
@@ -861,7 +861,14 @@ void sum_dhdl(gmx_enerdata_t *enerd, real *lambda, t_lambda *fepvals)
      * which is a very good approximation (except for exotic settings).
      * (investigate how to overcome this post 4.6 - MRS)
      */
-    enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
+    if (fepvals->separate_dvdl[efptBONDED])
+    {
+        enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
+    }
+    else
+    {
+        enerd->term[F_DVDL] += enerd->term[F_DVDL_CONSTR];
+    }
     enerd->term[F_DVDL_CONSTR] = 0;
 
     for (i = 0; i < fepvals->n_lambda; i++)

diff --git a/src/mdlib/shakef.c b/src/mdlib/shakef.c
index fd71f315482b60207850ce774d4f7b7016874519..a73680878088b6e06806c126d982b04d2f5f070b 100644 (file)
--- a/src/mdlib/shakef.c
+++ b/src/mdlib/shakef.c
@@ -438,13 +438,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;
         }