From 3ec1eb14d43512d8379f396bb57213c1c93ad33b Mon Sep 17 00:00:00 2001 From: Sander Pronk Date: Thu, 22 Nov 2012 18:28:02 +0100 Subject: [PATCH] Made the linear extrapolation in g_bar non-default behavior g_bar currently resorts to linear extrapolation of dH/dl if it doesn't find covering foreign lambda deltaH values. This is dangerous because the most likely case of this being used is when the user forgets to set foreign lambdas. This commit adds an option -extp to g_bar that the user must explictily set to enable the linear extrapolation. If not, and there is no foreign_lambda path, g_bar prints an error message with details of what to do. Change-Id: Ib39f964b4f43f79379f00a3331ea4dc88183ed63 --- src/tools/gmx_bar.c | 54 +++++++++++++++++++++++++++++++++-------------------- 1 file changed, 34 insertions(+), 20 deletions(-) diff --git a/src/tools/gmx_bar.c b/src/tools/gmx_bar.c index 6a49d389d7..25fd504f0d 100644 --- a/src/tools/gmx_bar.c +++ b/src/tools/gmx_bar.c @@ -732,7 +732,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; @@ -764,7 +765,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 */ @@ -780,16 +781,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; @@ -2472,25 +2477,26 @@ 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.", - "[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]", + "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 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 ", @@ -2503,6 +2509,12 @@ int gmx_bar(int argc,char *argv[]) "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 values. This is done with the[TT]-extp[tt] option," + "which assumes that the system's Hamiltonian depends linearly" + "on [GRK]lambda, 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]. ", "An error estimate taking into account time correlations ", @@ -2553,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" }, @@ -2561,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[] = { @@ -2662,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); -- 2.11.4.GIT