Allow printing potential energies to dhdl.xvg file

Previously, only total energies were printed out to the dhdl.xvg
file.  In many cases, it is more useful to print out just
the potential energies.  For example, when doing mdrun -rerun,
and comparing to output run without rerun, the velocities may not
be identical, which makes it difficult to compare the total
energies between runs when computing free energies and other
observables. dhdl-print-energy can now be 'no', 'total',
or 'potential'. The alternative 'kinetic' should generally not
be needed, since that can be done analytically in all cases.

Backwards compatible in the .tpr since gmx_bool is an int. The
old false value will be interpreted as 'no', and the 'yes' value
as 'total' energy, which will reproduce the old behavior.

Fixes #1329

Change-Id: I16828ef07c46bcfc61fe03744cbf251c00160636

diff --git a/share/html/online/mdp_opt.html b/share/html/online/mdp_opt.html
index 246f49cd0b5f14ca856ba2834ae8b32e7360ed15..310a6c9a3dad436b56456b9a91664f8dc58c1f2a 100644 (file)
--- a/share/html/online/mdp_opt.html
+++ b/share/html/online/mdp_opt.html
@@ -1630,7 +1630,8 @@ the molecule definition in the topology.</dd>
 <dt><b>dhdl-derivatives: (yes)</b></dt>
 <dd>If yes (the default), the derivatives of the Hamiltonian with respect to lambda at each <b>nstdhdl</b> step are written out. These values are needed for interpolation of linear energy differences with <tt>g_bar</tt> (although the same can also be achieved with the right <b>foreign lambda</b> setting, that may not be as flexible), or with thermodynamic integration</dd>
 <dt><b>dhdl-print-energy: (no)</b></dt>
-<dd> Include the total energy in the dhdl file.  This information is needed for later analysis if the states of interest in the free e energy calculation are at different temperatures.  If all are at the same temperature, this information is not needed.</dd>
+<dd> Include either the total or the potential energy in the dhdl file.  Options are 'no', 'potential', or 'total'. This information is needed for later free energy analysis if the states of interest are at different temperatures.  If all states are at the same temperature, this information is not needed. 'potential' is useful in case one is using <tt>mdrun -rerun</tt> to generate the <tt>dhdl.xvg</tt> file.  When rerunning from an existing trajectory, the kinetic energy will often not be correct, and thus one must compute the residual free energy from the potential alone, with the kinetic energy component computed analytically.
+</dd>
 <dt><b>separate-dhdl-file: (yes)</b></dt>
 <dd><dl compact>
 <dt><b>yes</b></dt>

diff --git a/src/gromacs/fileio/tpxio.c b/src/gromacs/fileio/tpxio.c
index 62e777ff55f11c5a897c3159ea2d301b4a8005a5..de2e2640c9b87b7515926aca732141101b0d75de 100644 (file)
--- a/src/gromacs/fileio/tpxio.c
+++ b/src/gromacs/fileio/tpxio.c
@@ -606,11 +606,11 @@ static void do_fepvals(t_fileio *fio, t_lambda *fepvals, gmx_bool bRead, int fil
     }
     if (file_version >= 79)
     {
-        gmx_fio_do_int(fio, fepvals->bPrintEnergy);
+        gmx_fio_do_int(fio, fepvals->edHdLPrintEnergy);
     }
     else
     {
-        fepvals->bPrintEnergy = FALSE;
+        fepvals->edHdLPrintEnergy = edHdLPrintEnergyNO;
     }
 
     /* handle lambda_neighbors */

diff --git a/src/gromacs/gmxlib/mvdata.c b/src/gromacs/gmxlib/mvdata.c
index 62b88991556b53c4ed3f7f836941cec56736f78b..b9d3a3a57b0960ddfea63e97334c7e2b3d3ee46d 100644 (file)
--- a/src/gromacs/gmxlib/mvdata.c
+++ b/src/gromacs/gmxlib/mvdata.c
@@ -574,7 +574,7 @@ static void bc_fepvals(const t_commrec *cr, t_lambda *fep)
     block_bc(cr, fep->init_lambda);
     block_bc(cr, fep->init_fep_state);
     block_bc(cr, fep->delta_lambda);
-    block_bc(cr, fep->bPrintEnergy);
+    block_bc(cr, fep->edHdLPrintEnergy);
     block_bc(cr, fep->n_lambda);
     if (fep->n_lambda > 0)
     {

diff --git a/src/gromacs/gmxlib/names.c b/src/gromacs/gmxlib/names.c
index 35ecbab8c0ea798d64db3c5fd12359f649404b7b..649db23f1bb32f078d22e2a37cfc012f1cc1864c 100644 (file)
--- a/src/gromacs/gmxlib/names.c
+++ b/src/gromacs/gmxlib/names.c
@@ -163,6 +163,10 @@ const char *efpt_singular_names[efptNR+1] = {
     "fep-lambda", "mass-lambda", "coul-lambda", "vdw-lambda", "bonded-lambda", "restraint-lambda", "temperature-lambda", NULL
 };
 
+const char *edHdLPrintEnergy_names[edHdLPrintEnergyNR+1] = {
+    "no", "total", "potential", "yes", NULL
+};
+
 const char *elamstats_names[elamstatsNR+1] = {
     "no", "metropolis-transition", "barker-transition", "minvar", "wang-landau", "weighted-wang-landau", NULL
 };

diff --git a/src/gromacs/gmxlib/txtdump.c b/src/gromacs/gmxlib/txtdump.c
index 182a41949e00b572251f8c6f55647f1aeaa894f9..7d94359da578ef303bf076eb1ce1bdfd2f212c88 100644 (file)
--- a/src/gromacs/gmxlib/txtdump.c
+++ b/src/gromacs/gmxlib/txtdump.c
@@ -738,7 +738,7 @@ static void pr_fepvals(FILE *fp, int indent, t_lambda *fep, gmx_bool bMDPformat)
         }
     }
     PI("calc-lambda-neighbors", fep->lambda_neighbors);
-    PS("dhdl-print-energy", EBOOL(fep->bPrintEnergy));
+    PS("dhdl-print-energy", edHdLPrintEnergy_names[fep->edHdLPrintEnergy]);
     PR("sc-alpha", fep->sc_alpha);
     PI("sc-power", fep->sc_power);
     PR("sc-r-power", fep->sc_r_power);

diff --git a/src/gromacs/gmxpreprocess/compute_io.c b/src/gromacs/gmxpreprocess/compute_io.c
index 67ef12f210e639a35d5d8d60df09ea90bf43ce02..52b62a10477d4c6da657c49c8848b26e74792dcd 100644 (file)
--- a/src/gromacs/gmxpreprocess/compute_io.c
+++ b/src/gromacs/gmxpreprocess/compute_io.c
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014, 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.
@@ -110,7 +110,7 @@ double compute_io(t_inputrec *ir, int natoms, gmx_groups_t *groups,
                 nchars += 5;   /* alchemical state */
             }
 
-            if (ir->fepvals->bPrintEnergy)
+            if (ir->fepvals->edHdLPrintEnergy != edHdLPrintEnergyNO)
             {
                 nchars += 12; /* energy for dhdl */
             }
@@ -126,7 +126,7 @@ double compute_io(t_inputrec *ir, int natoms, gmx_groups_t *groups,
                 {
                     ndh_tot += 1;
                 }
-                if (ir->fepvals->bPrintEnergy)
+                if (ir->fepvals->edHdLPrintEnergy != edHdLPrintEnergyNO)
                 {
                     ndh_tot += 1;
                 }

diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c
index e26f584d90badf5a0f80d61384d4ed13d613d282..9e2dbf0311905cf816cebef8f7ff9cf722e966a0 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.c
+++ b/src/gromacs/gmxpreprocess/readir.c
@@ -2139,7 +2139,7 @@ void get_ir(const char *mdparin, const char *mdparout,
     STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
     ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
     STYPE ("init-lambda-weights", is->lambda_weights, NULL);
-    EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
+    EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
     RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
     ITYPE ("sc-power", fep->sc_power, 1);
     RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
@@ -2364,11 +2364,23 @@ void get_ir(const char *mdparin, const char *mdparout,
         }
     }
 
-    if (ir->bSimTemp)
+    if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
     {
-        fep->bPrintEnergy = TRUE;
-        /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
-           if the temperature is changing. */
+        fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
+        warning_note(wi, "Old option for dhdl-print-energy given: "
+                     "changing \"yes\" to \"total\"\n");
+    }
+
+    if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
+    {
+        /* always print out the energy to dhdl if we are doing
+           expanded ensemble, since we need the total energy for
+           analysis if the temperature is changing. In some
+           conditions one may only want the potential energy, so
+           we will allow that if the appropriate mdp setting has
+           been enabled. Otherwise, total it is:
+         */
+        fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
     }
 
     if ((ir->efep != efepNO) || ir->bSimTemp)

diff --git a/src/gromacs/legacyheaders/names.h b/src/gromacs/legacyheaders/names.h
index f79fa9eea43ab7bca33093002b695fc1615d7fba..424d69114dc15432c6bfaec66bbe5a1a7ee364e6 100644 (file)
--- a/src/gromacs/legacyheaders/names.h
+++ b/src/gromacs/legacyheaders/names.h
@@ -76,6 +76,7 @@ extern const char *esimtemp_names[esimtempNR+1];
 extern const char *efep_names[efepNR+1];
 extern const char *efpt_names[efptNR+1];
 extern const char *efpt_singular_names[efptNR+1];
+extern const char *edHdLPrintEnergy_names[edHdLPrintEnergyNR+1];
 extern const char *elamstats_names[elamstatsNR+1];
 extern const char *elmcmove_names[elmcmoveNR+1];
 extern const char *elmceq_names[elmceqNR+1];

diff --git a/src/gromacs/legacyheaders/types/enums.h b/src/gromacs/legacyheaders/types/enums.h
index d8eaab27f694f310ae2b4f479af9e08ead8a91c2..4b1161a4fa7d10e1b68a355778a64e9fd8c2ec65 100644 (file)
--- a/src/gromacs/legacyheaders/types/enums.h
+++ b/src/gromacs/legacyheaders/types/enums.h
@@ -191,6 +191,13 @@ enum {
     efptFEP, efptMASS, efptCOUL, efptVDW, efptBONDED, efptRESTRAINT, efptTEMPERATURE, efptNR
 };
 
+/* Printing the energy to the free energy dhdl file. YES is an alias to TOTAL, and
+ * will be converted in readir, so we never have to account for it in code.
+ */
+enum {
+    edHdLPrintEnergyNO, edHdLPrintEnergyTOTAL, edHdLPrintEnergyPOTENTIAL, edHdLPrintEnergyYES, edHdLPrintEnergyNR
+};
+
 /* How the lambda weights are calculated:
    elamstatsMETROPOLIS = using the metropolis criteria
    elamstatsBARKER = using the Barker critera for transition weights - also called unoptimized Bennett

diff --git a/src/gromacs/legacyheaders/types/inputrec.h b/src/gromacs/legacyheaders/types/inputrec.h
index 6c790977393aea73765b0fd3b87ae786044fd373..bd9ea9818d083479641acf0017e176b1bbd51e90 100644 (file)
--- a/src/gromacs/legacyheaders/types/inputrec.h
+++ b/src/gromacs/legacyheaders/types/inputrec.h
@@ -151,7 +151,7 @@ typedef struct {
                                        valid value if positive)   */
     int      init_fep_state;        /* the initial number of the state                   */
     double   delta_lambda;          /* change of lambda per time step (fraction of (0.1) */
-    gmx_bool bPrintEnergy;          /* Whether to print the energy in the dhdl           */
+    int      edHdLPrintEnergy;      /* print no, total or potential energies in dhdl    */
     int      n_lambda;              /* The number of foreign lambda points               */
     double **all_lambda;            /* The array of all lambda values                    */
     int      lambda_neighbors;      /* The number of neighboring lambda states to

diff --git a/src/gromacs/mdlib/mdebin.c b/src/gromacs/mdlib/mdebin.c
index 5ef8b7d79a98ca52491ca4e5782f34fc8ac3ea24..bb4809174c0fde7a3e5b9128a2b3467d865b02c6 100644 (file)
--- a/src/gromacs/mdlib/mdebin.c
+++ b/src/gromacs/mdlib/mdebin.c
@@ -825,7 +825,7 @@ extern FILE *open_dhdl(const char *filename, const t_inputrec *ir,
         nsets += 1;   /*add fep state for expanded ensemble */
     }
 
-    if (fep->bPrintEnergy)
+    if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)
     {
         nsets += 1;  /* add energy to the dhdl as well */
     }
@@ -850,9 +850,18 @@ extern FILE *open_dhdl(const char *filename, const t_inputrec *ir,
         s         += 1;
     }
 
-    if (fep->bPrintEnergy)
+    if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)
     {
-        sprintf(buf, "%s (%s)", "Energy", unit_energy);
+        switch (fep->edHdLPrintEnergy)
+        {
+            case edHdLPrintEnergyPOTENTIAL:
+                sprintf(buf, "%s (%s)", "Potential Energy", unit_energy);
+                break;
+            case edHdLPrintEnergyTOTAL:
+            case edHdLPrintEnergyYES:
+            default:
+                sprintf(buf, "%s (%s)", "Total Energy", unit_energy);
+        }
         setname[s] = strdup(buf);
         s         += 1;
     }
@@ -900,7 +909,7 @@ extern FILE *open_dhdl(const char *filename, const t_inputrec *ir,
             nsetsbegin = 0;
         }
 
-        if (fep->bPrintEnergy)
+        if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)
         {
             nsetsbegin += 1;
         }
@@ -1208,9 +1217,19 @@ void upd_mdebin(t_mdebin       *md,
                 fprintf(md->fp_dhdl, " %4d", state->fep_state);
             }
             /* total energy (for if the temperature changes */
-            if (fep->bPrintEnergy)
+
+            if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)
             {
-                store_energy = enerd->term[F_ETOT];
+                switch (fep->edHdLPrintEnergy)
+                {
+                    case edHdLPrintEnergyPOTENTIAL:
+                        store_energy = enerd->term[F_EPOT];
+                        break;
+                    case edHdLPrintEnergyTOTAL:
+                    case edHdLPrintEnergyYES:
+                    default:
+                        store_energy = enerd->term[F_ETOT];
+                }
                 fprintf(md->fp_dhdl, " %#.8g", store_energy);
             }
 

diff --git a/src/gromacs/mdlib/mdebin_bar.c b/src/gromacs/mdlib/mdebin_bar.c
index 9f057bd6fd4126dda7848c64f8149e4a3c01f354..b855e33a9d6b9f791c9c7a8a69de841911b5adfc 100644 (file)
--- a/src/gromacs/mdlib/mdebin_bar.c
+++ b/src/gromacs/mdlib/mdebin_bar.c
@@ -469,7 +469,7 @@ void mde_delta_h_coll_init(t_mde_delta_h_coll *dhc, const t_inputrec *ir)
                 bExpanded = TRUE;
             }
             /* whether to print energies */
-            if (ir->fepvals->bPrintEnergy)
+            if (ir->fepvals->edHdLPrintEnergy != edHdLPrintEnergyNO)
             {
                 dhc->ndh += 1;
                 bEnergy   = TRUE;

diff --git a/src/gromacs/tools/compare.c b/src/gromacs/tools/compare.c
index cbbf9f8ee473253c8c4c12adaea955b9dcfeb0f9..34775b0ef43e938c6d82013c2af8659db614a843 100644 (file)
--- a/src/gromacs/tools/compare.c
+++ b/src/gromacs/tools/compare.c
@@ -725,7 +725,7 @@ static void cmp_fepvals(FILE *fp, t_lambda *fep1, t_lambda *fep2, real ftol, rea
     cmp_int(fp, "inputrec->fepvals->sc_power", -1, fep1->sc_power, fep2->sc_power);
     cmp_real(fp, "inputrec->fepvals->sc_r_power", -1, fep1->sc_r_power, fep2->sc_r_power, ftol, abstol);
     cmp_real(fp, "inputrec->fepvals->sc_sigma", -1, fep1->sc_sigma, fep2->sc_sigma, ftol, abstol);
-    cmp_bool(fp, "inputrec->fepvals->bPrintEnergy", -1, fep1->bPrintEnergy, fep1->bPrintEnergy);
+    cmp_int(fp, "inputrec->fepvals->edHdLPrintEnergy", -1, fep1->edHdLPrintEnergy, fep1->edHdLPrintEnergy);
     cmp_bool(fp, "inputrec->fepvals->bScCoul", -1, fep1->bScCoul, fep1->bScCoul);
     cmp_int(fp, "inputrec->separate_dhdl_file", -1, fep1->separate_dhdl_file, fep2->separate_dhdl_file);
     cmp_int(fp, "inputrec->dhdl_derivatives", -1, fep1->dhdl_derivatives, fep2->dhdl_derivatives);