Merge release-4-6 into master
[alexxy/gromacs.git] / src / tools / gmx_bar.c
index a7eb79b98c4a3b355b1e839aada455b204af0853..bebdcb37f273d6cefc69acd3330551d7e9e2f4db 100644 (file)
@@ -727,7 +727,8 @@ void lambdas_histogram(lambda_t *bl_head, const char *filename,
 
 /* create a collection (array) of barres_t object given a ordered linked list 
    of barlamda_t sample collections */
-static barres_t *barres_list_create(lambda_t *bl_head, int *nres)
+static barres_t *barres_list_create(lambda_t *bl_head, int *nres,
+                                    gmx_bool use_dhdl)
 {
     lambda_t *bl;
     int nlambda=0;
@@ -759,7 +760,7 @@ static barres_t *barres_list_create(lambda_t *bl_head, int *nres)
 
         barres_init(br);
 
-        if (!scprev && !sc)
+        if (use_dhdl)
         {
             /* we use dhdl */
 
@@ -775,16 +776,20 @@ static barres_t *barres_list_create(lambda_t *bl_head, int *nres)
             {
                 gmx_fatal(FARGS,"Some dhdl files contain only one value (dH/dl), while others \ncontain multiple values (dH/dl and/or Delta H), will not proceed \nbecause of possible inconsistencies.\n");
             }
+        } 
+        else if (!scprev && !sc)
+        {
+            gmx_fatal(FARGS,"There is no path from lambda=%g -> %g that is covered by foreign lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl by calculating the averages of dH/dl\nwith g_analyze and integrating them.\nAlternatively, use the -extp option if (and only if) the Hamiltonian\ndepends linearly on lambda, which is NOT normally the case.\n", bl->prev->lambda, bl->lambda);
         }
         
         /* normal delta H */
         if (!scprev)
         {
-            gmx_fatal(FARGS,"Could not find a set for lambda = %g in the files for lambda = %g",bl->lambda,bl->prev->lambda);
+            gmx_fatal(FARGS,"Could not find a set for foreign lambda = %g\nin the files for lambda = %g",bl->lambda,bl->prev->lambda);
         }
         if (!sc)
         {
-            gmx_fatal(FARGS,"Could not find a set for lambda = %g in the files for lambda = %g",bl->prev->lambda,bl->lambda);
+            gmx_fatal(FARGS,"Could not find a set for foreign lambda = %g\nin the files for lambda = %g",bl->prev->lambda,bl->lambda);
         }
         br->a = scprev;
         br->b = sc;
@@ -2467,36 +2472,47 @@ int gmx_bar(int argc,char *argv[])
         "controlled by a parameter, [GRK]lambda[grk] (see the [TT].mdp[tt] parameter",
         "[TT]init_lambda[tt]). The BAR method calculates a ratio of weighted",
         "average of the Hamiltonian difference of state B given state A and",
-        "vice versa. If the Hamiltonian does not depend linearly on [GRK]lambda[grk]",
-        "(in which case we can extrapolate the derivative of the Hamiltonian",
-        "with respect to [GRK]lambda[grk], as is the default when [TT]free_energy[tt] is on),",
-        "the energy differences to the other state need to be calculated",
-        "explicitly during the simulation. This can be controlled with",
+        "vice versa.",
+        "The energy differences to the other state must be calculated",
+        "explicitly during the simulation. This can be done with",
         "the [TT].mdp[tt] option [TT]foreign_lambda[tt].[PAR]",
 
         "Input option [TT]-f[tt] expects multiple [TT]dhdl.xvg[tt] files. ",
         "Two types of input files are supported:[BR]",
-        "[TT]*[tt]  Files with only one [IT]y[it]-value, for such files it is assumed ",
-        "   that the [IT]y[it]-value is dH/d[GRK]lambda[grk] and that the Hamiltonian depends ",
-        "   linearly on [GRK]lambda[grk]. The [GRK]lambda[grk] value of the simulation is inferred ",
-        "   from the subtitle (if present), otherwise from a number in the",
-        "   subdirectory in the file name.",
+        "[TT]*[tt]  Files with more than one [IT]y[it]-value. ",
+        "The files should have columns ",
+        "with dH/d[GRK]lambda[grk] and [GRK]Delta[grk][GRK]lambda[grk]. ",
+        "The [GRK]lambda[grk] values are inferred ",
+        "from the legends: [GRK]lambda[grk] of the simulation from the legend of ",
+        "dH/d[GRK]lambda[grk] and the foreign [GRK]lambda[grk] values from the ",
+        "legends of Delta H",
         "[BR]",
-        "[TT]*[tt]  Files with more than one [IT]y[it]-value. The files should have columns ",
-        "   with dH/d[GRK]lambda[grk] and [GRK]Delta[grk][GRK]lambda[grk]. The [GRK]lambda[grk] values are inferred ",
-        "   from the legends: [GRK]lambda[grk] of the simulation from the legend of dH/d[GRK]lambda[grk] ",
-        "   and the foreign [GRK]lambda[grk] values from the legends of Delta H.[PAR]",
-        "The [GRK]lambda[grk] of the simulation is parsed from [TT]dhdl.xvg[tt] file's legend ",
-        "containing the string 'dH', the foreign [GRK]lambda[grk] values from the legend ",
-        "containing the capitalized letters 'D' and 'H'. The temperature ",
-        "is parsed from the legend line containing 'T ='.[PAR]",
+        "[TT]*[tt]  Files with only one [IT]y[it]-value. Using the",
+        "[TT]-extp[tt] option for these files, it is assumed",
+        "that the [IT]y[it]-value is dH/d[GRK]lambda[grk] and that the ",
+        "Hamiltonian depends linearly on [GRK]lambda[grk]. ",
+        "The [GRK]lambda[grk] value of the simulation is inferred from the ",
+        "subtitle (if present), otherwise from a number in the subdirectory ",
+        "in the file name.[PAR]",
+
+        "The [GRK]lambda[grk] of the simulation is parsed from ",
+        "[TT]dhdl.xvg[tt] file's legend containing the string 'dH', the ",
+        "foreign [GRK]lambda[grk] values from the legend containing the ",
+        "capitalized letters 'D' and 'H'. The temperature is parsed from ",
+        "the legend line containing 'T ='.[PAR]",
 
         "The input option [TT]-g[tt] expects multiple [TT].edr[tt] files. ",
-        "These can contain either lists of energy differences (see the",
-        "[TT].mdp[tt] option [TT]separate_dhdl_file[tt]), or a series of histograms",
-        "(see the [TT].mdp[tt] options [TT]dh_hist_size[tt] and [TT]dh_hist_spacing[tt]).",
-        "The temperature and [GRK]lambda[grk] values are automatically deduced from",
-        "the [TT]ener.edr[tt] file.[PAR]"
+        "These can contain either lists of energy differences (see the ",
+        "[TT].mdp[tt] option [TT]separate_dhdl_file[tt]), or a series of ",
+        "histograms (see the [TT].mdp[tt] options [TT]dh_hist_size[tt] and ",
+        "[TT]dh_hist_spacing[tt]).", "The temperature and [GRK]lambda[grk] ",
+        "values are automatically deduced from the [TT]ener.edr[tt] file.[PAR]",
+
+        "In addition to the [TT].mdp[tt] option [TT]foreign_lambda[tt], ",
+        "the energy difference can also be extrapolated from the ",
+        "dH/d[GRK]lambda[grk] values. This is done with the[TT]-extp[tt]",
+        "option, which assumes that the system's Hamiltonian depends linearly",
+        "on [GRK]lambda[grk], which is not normally the case.[PAR]",
 
         "The free energy estimates are determined using BAR with bisection, ",
         "with the precision of the output set with [TT]-prec[tt]. ",
@@ -2508,12 +2524,13 @@ int gmx_bar(int argc,char *argv[])
         "over 5 blocks. A range of block numbers for error estimation can ",
         "be provided with the options [TT]-nbmin[tt] and [TT]-nbmax[tt].[PAR]",
 
-        "[TT]g_bar[tt] tries to aggregate samples with the same 'native' and 'foreign'",
-        "[GRK]lambda[grk] values, but always assumes independent samples. [BB]Note[bb] that",
-        "when aggregating energy differences/derivatives with different",
-        "sampling intervals, this is almost certainly not correct. Usually",
-        "subsequent energies are correlated and different time intervals mean",
-        "different degrees of correlation between samples.[PAR]",
+        "[TT]g_bar[tt] tries to aggregate samples with the same 'native' and ",
+        "'foreign' [GRK]lambda[grk] values, but always assumes independent ",
+        "samples. [BB]Note[bb] that when aggregating energy ",
+        "differences/derivatives with different sampling intervals, this is ",
+        "almost certainly not correct. Usually subsequent energies are ",
+        "correlated and different time intervals mean different degrees ",
+        "of correlation between samples.[PAR]",
 
         "The results are split in two parts: the last part contains the final ",
         "results in kJ/mol, together with the error estimate for each part ",
@@ -2548,6 +2565,7 @@ int gmx_bar(int argc,char *argv[])
     static real begin=0,end=-1,temp=-1;
     int nd=2,nbmin=5,nbmax=5;
     int nbin=100;
+    gmx_bool use_dhdl=FALSE;
     gmx_bool calc_s,calc_v;
     t_pargs pa[] = {
         { "-b",    FALSE, etREAL, {&begin},  "Begin time for BAR" },
@@ -2556,7 +2574,8 @@ int gmx_bar(int argc,char *argv[])
         { "-prec", FALSE, etINT,  {&nd},     "The number of digits after the decimal point" },
         { "-nbmin",  FALSE, etINT,  {&nbmin}, "Minimum number of blocks for error estimation" },
         { "-nbmax",  FALSE, etINT,  {&nbmax}, "Maximum number of blocks for error estimation" },
-        { "-nbin",  FALSE, etINT, {&nbin}, "Number of bins for histogram output"}
+        { "-nbin",  FALSE, etINT, {&nbin}, "Number of bins for histogram output"},
+        { "-extp",  FALSE, etBOOL, {&use_dhdl}, "Whether to linearly extrapolate dH/dl values to use as energies"}
     };
     
     t_filenm   fnm[] = {
@@ -2657,7 +2676,7 @@ int gmx_bar(int argc,char *argv[])
     }
    
     /* assemble the output structures from the lambdas */
-    results=barres_list_create(lb, &nresults);
+    results=barres_list_create(lb, &nresults, use_dhdl);
 
     sum_disc_err=barres_list_max_disc_err(results, nresults);