Fixed problem with bond constraint dhdl being left out

At some point, the dhdl from the holonomic bond constraints got left out,
because the term got zeroed before being added to the bonded energy.

Also, I made some of the variable names involved more descriptive to
clarify the logic.

Removed velocity constraint calls, as don't contribute to dhdl.

Fixes redmine #1255

Change-Id: I2eebf357d1a922b7636dc015b09f286827283dd0

diff --git a/include/types/idef.h b/include/types/idef.h
index 91e4ce633511b6c503dd9d440e7451310a54fef8..5f4ce1e0b7223f75da4fef76f17eebec8690c8e7 100644 (file)
--- a/include/types/idef.h
+++ b/include/types/idef.h
@@ -136,7 +136,7 @@ enum {
     F_VTEMP_NOLONGERUSED,
     F_PDISPCORR,
     F_PRES,
-    F_DHDL_CON,
+    F_DVDL_CONSTR,
     F_DVDL,
     F_DKDL,
     F_DVDL_COUL,

diff --git a/src/gmxlib/tpxio.c b/src/gmxlib/tpxio.c
index d916097ac775ba2fdab355990f8127d446f48677..f79118c66e87bd97bcf0a533c210607b2b08686e 100644 (file)
--- a/src/gmxlib/tpxio.c
+++ b/src/gmxlib/tpxio.c
@@ -174,14 +174,13 @@ static const t_ftupd ftupd[] = {
     { 46, F_ECONSERVED        },
     { 69, F_VTEMP_NOLONGERUSED},
     { 66, F_PDISPCORR         },
-    { 54, F_DHDL_CON          },
+    { 54, F_DVDL_CONSTR       },
     { 76, F_ANHARM_POL        },
     { 79, F_DVDL_COUL         },
     { 79, F_DVDL_VDW,         },
     { 79, F_DVDL_BONDED,      },
     { 79, F_DVDL_RESTRAINT    },
     { 79, F_DVDL_TEMPERATURE  },
-    { 54, F_DHDL_CON          }
 };
 #define NFTUPD asize(ftupd)
 

diff --git a/src/kernel/md.c b/src/kernel/md.c
index 886f9cbcc6c06db2450e6133e9ef0740d2f4810d..24eb4d6cf69f512d1f5f98a753032a18d68e8f95 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -206,7 +206,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     gmx_bool          bResetCountersHalfMaxH = FALSE;
     gmx_bool          bVV, bIterativeCase, bFirstIterate, bTemp, bPres, bTrotter;
     gmx_bool          bUpdateDoLR;
-    real              mu_aver = 0, dvdl;
+    real              mu_aver = 0, dvdl_constr;
     int               a0, a1, gnx = 0, ii;
     atom_id          *grpindex = NULL;
     char             *grpname;
@@ -1273,9 +1273,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 bOK = TRUE;
                 if (!bRerunMD || rerun_fr.bV || bForceUpdate)     /* Why is rerun_fr.bV here?  Unclear. */
                 {
-                    dvdl = 0;
-
-                    update_constraints(fplog, step, &dvdl, ir, ekind, mdatoms,
+                    update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
                                        state, fr->bMolPBC, graph, f,
                                        &top->idef, shake_vir, NULL,
                                        cr, nrnb, wcycle, upd, constr,
@@ -1370,7 +1368,6 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
             if (bTrotter && !bInitStep)
             {
-                enerd->term[F_DVDL_BONDED] += dvdl;        /* only add after iterations */
                 copy_mat(shake_vir, state->svir_prev);
                 copy_mat(force_vir, state->fvir_prev);
                 if (IR_NVT_TROTTER(ir) && ir->eI == eiVV)
@@ -1386,12 +1383,6 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 copy_rvecn(cbuf, state->v, 0, state->natoms);
             }
 
-            if (fr->bSepDVDL && fplog && do_log)
-            {
-                fprintf(fplog, sepdvdlformat, "Constraint", 0.0, dvdl);
-            }
-            enerd->term[F_DVDL_BONDED] += dvdl;
-
             GMX_MPE_LOG(ev_timestep1);
         }
 
@@ -1660,7 +1651,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
                 if (constr && bIfRandomize)
                 {
-                    update_constraints(fplog, step, &dvdl, ir, ekind, mdatoms,
+                    update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
                                        state, fr->bMolPBC, graph, f,
                                        &top->idef, tmp_vir, NULL,
                                        cr, nrnb, wcycle, upd, constr,
@@ -1698,10 +1689,11 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             copy_mat(state->box, lastbox);
 
             bOK = TRUE;
+            dvdl_constr = 0;
+
             if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate))
             {
                 wallcycle_start(wcycle, ewcUPDATE);
-                dvdl = 0;
                 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
                 if (bTrotter)
                 {
@@ -1761,7 +1753,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                               ekind, M, wcycle, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef);
                 wallcycle_stop(wcycle, ewcUPDATE);
 
-                update_constraints(fplog, step, &dvdl, ir, ekind, mdatoms, state,
+                update_constraints(fplog, step, &dvdl_constr, ir, ekind, mdatoms, state,
                                    fr->bMolPBC, graph, f,
                                    &top->idef, shake_vir, force_vir,
                                    cr, nrnb, wcycle, upd, constr,
@@ -1794,7 +1786,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                      * to numerical errors, or are they important
                      * physically? I'm thinking they are just errors, but not completely sure.
                      * For now, will call without actually constraining, constr=NULL*/
-                    update_constraints(fplog, step, &dvdl, ir, ekind, mdatoms,
+                    update_constraints(fplog, step, NULL, ir, ekind, mdatoms,
                                        state, fr->bMolPBC, graph, f,
                                        &top->idef, tmp_vir, force_vir,
                                        cr, nrnb, wcycle, upd, NULL,
@@ -1808,9 +1800,9 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
                 if (fr->bSepDVDL && fplog && do_log)
                 {
-                    fprintf(fplog, sepdvdlformat, "Constraint dV/dl", 0.0, dvdl);
+                    fprintf(fplog, sepdvdlformat, "Constraint dV/dl", 0.0, dvdl_constr);
                 }
-                enerd->term[F_DVDL_BONDED] += dvdl;
+                enerd->term[F_DVDL_CONSTR] += dvdl_constr;
             }
             else if (graph)
             {
@@ -1891,7 +1883,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
 
         /* only add constraint dvdl after constraints */
-        enerd->term[F_DVDL_BONDED] += dvdl;
+        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 c331ecb552bcaaef896a00046c73434c416b5bc7..4666743c00830bc312ad98d9b289853bbad8e132 100644 (file)
--- a/src/mdlib/force.c
+++ b/src/mdlib/force.c
@@ -861,6 +861,8 @@ 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];
+    enerd->term[F_DVDL_CONSTR] = 0;
 
     for (i = 0; i < fepvals->n_lambda; i++)
     {                                         /* note we are iterating over fepvals here!

diff --git a/src/mdlib/minimize.c b/src/mdlib/minimize.c
index 3b76d6c68e72902cd0d324af1cc84fb29725dc39..250e64ae6bca2f205383f3871787430732589155 100644 (file)
--- a/src/mdlib/minimize.c
+++ b/src/mdlib/minimize.c
@@ -317,7 +317,7 @@ void init_em(FILE *fplog, const char *title,
              gmx_mdoutf_t **outf, t_mdebin **mdebin)
 {
     int  start, homenr, i;
-    real dvdlambda;
+    real dvdl_constr;
 
     if (fplog)
     {
@@ -433,11 +433,11 @@ void init_em(FILE *fplog, const char *title,
         if (!ir->bContinuation)
         {
             /* Constrain the starting coordinates */
-            dvdlambda = 0;
+            dvdl_constr = 0;
             constrain(PAR(cr) ? NULL : fplog, TRUE, TRUE, constr, &(*top)->idef,
                       ir, NULL, cr, -1, 0, mdatoms,
                       ems->s.x, ems->s.x, NULL, fr->bMolPBC, ems->s.box,
-                      ems->s.lambda[efptFEP], &dvdlambda,
+                      ems->s.lambda[efptFEP], &dvdl_constr,
                       NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
         }
     }
@@ -552,7 +552,7 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
     int      i;
     int      start, end;
     rvec    *x1, *x2;
-    real     dvdlambda;
+    real     dvdl_constr;
 
     s1 = &ems1->s;
     s2 = &ems2->s;
@@ -650,11 +650,11 @@ static void do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
     if (constr)
     {
         wallcycle_start(wcycle, ewcCONSTR);
-        dvdlambda = 0;
+        dvdl_constr = 0;
         constrain(NULL, TRUE, TRUE, constr, &top->idef,
                   ir, NULL, cr, count, 0, md,
                   s1->x, s2->x, NULL, bMolPBC, s2->box,
-                  s2->lambda[efptBONDED], &dvdlambda,
+                  s2->lambda[efptBONDED], &dvdl_constr,
                   NULL, NULL, nrnb, econqCoord, FALSE, 0, 0);
         wallcycle_stop(wcycle, ewcCONSTR);
     }
@@ -695,7 +695,7 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
     gmx_bool bNS;
     int      nabnsb;
     tensor   force_vir, shake_vir, ekin;
-    real     dvdlambda, prescorr, enercorr, dvdlcorr;
+    real     dvdl_constr, prescorr, enercorr, dvdlcorr;
     real     terminate = 0;
 
     /* Set the time to the initial time, the time does not change during EM */
@@ -790,17 +790,17 @@ static void evaluate_energy(FILE *fplog, gmx_bool bVerbose, t_commrec *cr,
     {
         /* Project out the constraint components of the force */
         wallcycle_start(wcycle, ewcCONSTR);
-        dvdlambda = 0;
+        dvdl_constr = 0;
         constrain(NULL, FALSE, FALSE, constr, &top->idef,
                   inputrec, NULL, cr, count, 0, mdatoms,
                   ems->s.x, ems->f, ems->f, fr->bMolPBC, ems->s.box,
-                  ems->s.lambda[efptBONDED], &dvdlambda,
+                  ems->s.lambda[efptBONDED], &dvdl_constr,
                   NULL, &shake_vir, nrnb, econqForceDispl, FALSE, 0, 0);
         if (fr->bSepDVDL && fplog)
         {
-            fprintf(fplog, sepdvdlformat, "Constraints", t, dvdlambda);
+            fprintf(fplog, sepdvdlformat, "Constraints", t, dvdl_constr);
         }
-        enerd->term[F_DVDL_BONDED] += dvdlambda;
+        enerd->term[F_DVDL_CONSTR] += dvdl_constr;
         m_add(force_vir, shake_vir, vir);
         wallcycle_stop(wcycle, ewcCONSTR);
     }
@@ -2375,7 +2375,7 @@ double do_steep(FILE *fplog, t_commrec *cr,
     gmx_global_stat_t gstat;
     t_graph          *graph;
     real              stepsize, constepsize;
-    real              ustep, dvdlambda, fnormn;
+    real              ustep, fnormn;
     gmx_mdoutf_t     *outf;
     t_mdebin         *mdebin;
     gmx_bool          bDone, bAbort, do_x, do_f;