2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
56 #include "gmx_fatal.h"
63 /* Suppress Cygwin compiler warnings from using newlib version of
70 /* Structure for the names of lambda vector components */
71 typedef struct lambda_components_t
73 char **names; /* Array of strings with names for the lambda vector
75 int N; /* The number of components */
76 int Nalloc; /* The number of allocated components */
77 } lambda_components_t;
79 /* Structure for a lambda vector or a dhdl derivative direction */
80 typedef struct lambda_vec_t
82 double *val; /* The lambda vector component values. Only valid if
84 int dhdl; /* The coordinate index for the derivative described by this
86 const lambda_components_t *lc; /* the associated lambda_components
88 int index; /* The state number (init-lambda-state) of this lambda
89 vector, if known. If not, it is set to -1 */
92 /* the dhdl.xvg data from a simulation */
96 int ftp; /* file type */
97 int nset; /* number of lambdas, including dhdl */
98 int *np; /* number of data points (du or hists) per lambda */
99 int np_alloc; /* number of points (du or hists) allocated */
100 double temp; /* temperature */
101 lambda_vec_t *lambda; /* the lambdas (of first index for y). */
102 double *t; /* the times (of second index for y) */
103 double **y; /* the dU values. y[0] holds the derivative, while
104 further ones contain the energy differences between
105 the native lambda and the 'foreign' lambdas. */
106 lambda_vec_t native_lambda; /* the native lambda */
108 struct xvg_t *next, *prev; /*location in the global linked list of xvg_ts*/
112 typedef struct hist_t
114 unsigned int *bin[2]; /* the (forward + reverse) histogram values */
115 double dx[2]; /* the histogram spacing. The reverse
116 dx is the negative of the forward dx.*/
117 gmx_large_int_t x0[2]; /* the (forward + reverse) histogram start
120 int nbin[2]; /* the (forward+reverse) number of bins */
121 gmx_large_int_t sum; /* the total number of counts. Must be
122 the same for forward + reverse. */
123 int nhist; /* number of hist datas (forward or reverse) */
125 double start_time, delta_time; /* start time, end time of histogram */
129 /* an aggregate of samples for partial free energy calculation */
130 typedef struct samples_t
132 lambda_vec_t *native_lambda; /* pointer to native lambda vector */
133 lambda_vec_t *foreign_lambda; /* pointer to foreign lambda vector */
134 double temp; /* the temperature */
135 gmx_bool derivative; /* whether this sample is a derivative */
137 /* The samples come either as either delta U lists: */
138 int ndu; /* the number of delta U samples */
139 double *du; /* the delta u's */
140 double *t; /* the times associated with those samples, or: */
141 double start_time,delta_time;/*start time and delta time for linear time*/
143 /* or as histograms: */
144 hist_t *hist; /* a histogram */
146 /* allocation data: (not NULL for data 'owned' by this struct) */
147 double *du_alloc, *t_alloc; /* allocated delta u arrays */
148 size_t ndu_alloc, nt_alloc; /* pre-allocated sizes */
149 hist_t *hist_alloc; /* allocated hist */
151 gmx_large_int_t ntot; /* total number of samples */
152 const char *filename; /* the file name this sample comes from */
155 /* a sample range (start to end for du-style data, or boolean
156 for both du-style data and histograms */
157 typedef struct sample_range_t
159 int start, end; /* start and end index for du style data */
160 gmx_bool use; /* whether to use this sample */
162 samples_t *s; /* the samples this range belongs to */
166 /* a collection of samples for a partial free energy calculation
167 (i.e. the collection of samples from one native lambda to one
169 typedef struct sample_coll_t
171 lambda_vec_t *native_lambda; /* these should be the same for all samples
173 lambda_vec_t *foreign_lambda; /* collection */
174 double temp; /* the temperature */
176 int nsamples; /* the number of samples */
177 samples_t **s; /* the samples themselves */
178 sample_range_t *r; /* the sample ranges */
179 int nsamples_alloc; /* number of allocated samples */
181 gmx_large_int_t ntot; /* total number of samples in the ranges of
184 struct sample_coll_t *next, *prev; /* next and previous in the list */
187 /* all the samples associated with a lambda point */
188 typedef struct lambda_data_t
190 lambda_vec_t *lambda; /* the native lambda (at start time if dynamic) */
191 double temp; /* temperature */
193 sample_coll_t *sc; /* the samples */
195 sample_coll_t sc_head; /*the pre-allocated list head for the linked list.*/
197 struct lambda_data_t *next, *prev; /* the next and prev in the list */
200 /* Top-level data structure of simulation data */
201 typedef struct sim_data_t
203 lambda_data_t *lb; /* a lambda data linked list */
204 lambda_data_t lb_head; /* The head element of the linked list */
206 lambda_components_t lc; /* the allowed components of the lambda
210 /* Top-level data structure with calculated values. */
212 sample_coll_t *a, *b; /* the simulation data */
214 double dg; /* the free energy difference */
215 double dg_err; /* the free energy difference */
217 double dg_disc_err; /* discretization error */
218 double dg_histrange_err; /* histogram range error */
220 double sa; /* relative entropy of b in state a */
221 double sa_err; /* error in sa */
222 double sb; /* relative entropy of a in state b */
223 double sb_err; /* error in sb */
225 double dg_stddev; /* expected dg stddev per sample */
226 double dg_stddev_err; /* error in dg_stddev */
230 /* Initialize a lambda_components structure */
231 static void lambda_components_init(lambda_components_t *lc)
235 snew(lc->names, lc->Nalloc);
238 /* Add a component to a lambda_components structure */
239 static void lambda_components_add(lambda_components_t *lc,
240 const char *name, size_t name_length)
242 while (lc->N + 1 > lc->Nalloc)
244 lc->Nalloc = (lc->Nalloc == 0) ? 2 : 2*lc->Nalloc;
245 srealloc( lc->names, lc->Nalloc );
247 snew(lc->names[lc->N], name_length+1);
248 strncpy(lc->names[lc->N], name, name_length);
252 /* check whether a component with index 'index' matches the given name, or
253 is also NULL. Returns TRUE if this is the case.
254 the string name does not need to end */
255 static gmx_bool lambda_components_check(const lambda_components_t *lc,
263 if (name == NULL && lc->names[index] == NULL)
265 if ((name == NULL) != (lc->names[index] == NULL))
267 len = strlen(lc->names[index]);
268 if (len != name_length)
270 if (strncmp(lc->names[index], name, name_length)== 0)
275 /* Find the index of a given lambda component name, or -1 if not found */
276 static int lambda_components_find(const lambda_components_t *lc,
284 if (strncmp(lc->names[i], name, name_length) == 0)
292 /* initialize a lambda vector */
293 static void lambda_vec_init(lambda_vec_t *lv, const lambda_components_t *lc)
295 snew(lv->val, lc->N);
301 static void lambda_vec_destroy(lambda_vec_t *lv)
306 static void lambda_vec_copy(lambda_vec_t *lv, const lambda_vec_t *orig)
310 lambda_vec_init(lv, orig->lc);
312 lv->index=orig->index;
313 for(i=0;i<lv->lc->N;i++)
315 lv->val[i] = orig->val[i];
319 /* write a lambda vec to a preallocated string */
320 static void lambda_vec_print(const lambda_vec_t *lv, char *str, gmx_bool named)
325 str[0]=0; /* reset the string */
330 str += sprintf(str, "delta H to ");
334 str += sprintf(str, "(");
336 for(i=0;i<lv->lc->N;i++)
338 str += sprintf(str, "%g", lv->val[i]);
341 str += sprintf(str, ", ");
346 str += sprintf(str, ")");
351 /* this lambda vector describes a derivative */
352 str += sprintf(str, "dH/dl");
353 if (strlen(lv->lc->names[lv->dhdl]) > 0)
355 str += sprintf(str, " (%s)", lv->lc->names[lv->dhdl]);
360 /* write a shortened version of the lambda vec to a preallocated string */
361 static void lambda_vec_print_short(const lambda_vec_t *lv, char *str)
368 sprintf(str, "%6d", lv->index);
374 sprintf(str, "%6.3f", lv->val[0]);
378 sprintf(str, "dH/dl[%d]", lv->dhdl);
383 /* write an intermediate version of two lambda vecs to a preallocated string */
384 static void lambda_vec_print_intermediate(const lambda_vec_t *a,
385 const lambda_vec_t *b, char *str)
391 if ( (a->index >= 0) && (b->index>=0) )
393 sprintf(str, "%6.3f", ((double)a->index+(double)b->index)/2.);
397 if ( (a->dhdl < 0) && (b->dhdl < 0) )
399 sprintf(str, "%6.3f", (a->val[0]+b->val[0])/2.);
406 /* calculate the difference in lambda vectors: c = a-b.
407 c must be initialized already, and a and b must describe non-derivative
409 static void lambda_vec_diff(const lambda_vec_t *a, const lambda_vec_t *b,
414 if ( (a->dhdl > 0) || (b->dhdl > 0) )
417 "Trying to calculate the difference between derivatives instead of lambda points");
419 if ((a->lc != b->lc) || (a->lc != c->lc) )
422 "Trying to calculate the difference lambdas with differing basis set");
424 for(i=0;i<a->lc->N;i++)
426 c->val[i] = a->val[i] - b->val[i];
430 /* calculate and return the absolute difference in lambda vectors: c = |a-b|.
431 a and b must describe non-derivative lambda points */
432 static double lambda_vec_abs_diff(const lambda_vec_t *a, const lambda_vec_t *b)
437 if ( (a->dhdl > 0) || (b->dhdl > 0) )
440 "Trying to calculate the difference between derivatives instead of lambda points");
445 "Trying to calculate the difference lambdas with differing basis set");
447 for(i=0;i<a->lc->N;i++)
449 double df=a->val[i] - b->val[i];
456 /* check whether two lambda vectors are the same */
457 static gmx_bool lambda_vec_same(const lambda_vec_t *a, const lambda_vec_t *b)
465 for(i=0;i<a->lc->N;i++)
467 if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
476 /* they're derivatives, so we check whether the indices match */
477 return (a->dhdl == b->dhdl);
481 /* Compare the sort order of two foreign lambda vectors
483 returns 1 if a is 'bigger' than b,
484 returns 0 if they're the same,
485 returns -1 if a is 'smaller' than b.*/
486 static gmx_bool lambda_vec_cmp_foreign(const lambda_vec_t *a,
487 const lambda_vec_t *b)
490 double norm_a=0, norm_b=0;
491 gmx_bool different=FALSE;
495 gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
497 /* if either one has an index we sort based on that */
498 if ((a->index >= 0) || (b->index >= 0))
500 if (a->index == b->index)
502 return (a->index > b->index) ? 1 : -1;
504 if (a->dhdl >= 0 || b->dhdl >= 0)
506 /* lambda vectors that are derivatives always sort higher than those
507 without derivatives */
508 if ((a->dhdl >= 0) != (b->dhdl >= 0) )
510 return (a->dhdl >= 0) ? 1 : -1 ;
512 return a->dhdl > b->dhdl;
515 /* neither has an index, so we can only sort on the lambda components,
516 which is only valid if there is one component */
517 for(i=0;i<a->lc->N;i++)
519 if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
523 norm_a += a->val[i]*a->val[i];
524 norm_b += b->val[i]*b->val[i];
530 return norm_a > norm_b;
533 /* Compare the sort order of two native lambda vectors
535 returns 1 if a is 'bigger' than b,
536 returns 0 if they're the same,
537 returns -1 if a is 'smaller' than b.*/
538 static gmx_bool lambda_vec_cmp_native(const lambda_vec_t *a,
539 const lambda_vec_t *b)
545 gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
547 /* if either one has an index we sort based on that */
548 if ((a->index >= 0) || (b->index >= 0))
550 if (a->index == b->index)
552 return (a->index > b->index) ? 1 : -1;
554 /* neither has an index, so we can only sort on the lambda components,
555 which is only valid if there is one component */
559 "Can't compare lambdas with no index and > 1 component");
561 if (a->dhdl >= 0 || b->dhdl >= 0)
564 "Can't compare native lambdas that are derivatives");
566 if (gmx_within_tol(a->val[0], b->val[0], 10*GMX_REAL_EPS))
570 return a->val[0] > b->val[0] ? 1 : -1;
576 static void hist_init(hist_t *h, int nhist, int *nbin)
581 gmx_fatal(FARGS, "histogram with more than two sets of data!");
585 snew(h->bin[i], nbin[i]);
588 h->start_time=h->delta_time=0;
595 static void hist_destroy(hist_t *h)
601 static void xvg_init(xvg_t *ba)
610 static void samples_init(samples_t *s, lambda_vec_t *native_lambda,
611 lambda_vec_t *foreign_lambda, double temp,
612 gmx_bool derivative, const char *filename)
614 s->native_lambda=native_lambda;
615 s->foreign_lambda=foreign_lambda;
617 s->derivative=derivative;
622 s->start_time = s->delta_time = 0;
631 s->filename=filename;
634 static void sample_range_init(sample_range_t *r, samples_t *s)
642 static void sample_coll_init(sample_coll_t *sc, lambda_vec_t *native_lambda,
643 lambda_vec_t *foreign_lambda, double temp)
645 sc->native_lambda = native_lambda;
646 sc->foreign_lambda = foreign_lambda;
652 sc->nsamples_alloc=0;
655 sc->next=sc->prev=NULL;
658 static void sample_coll_destroy(sample_coll_t *sc)
660 /* don't free the samples themselves */
666 static void lambda_data_init(lambda_data_t *l, lambda_vec_t *native_lambda,
669 l->lambda=native_lambda;
677 sample_coll_init(l->sc, native_lambda, NULL, 0.);
682 static void barres_init(barres_t *br)
698 /* calculate the total number of samples in a sample collection */
699 static void sample_coll_calc_ntot(sample_coll_t *sc)
704 for(i=0;i<sc->nsamples;i++)
710 sc->ntot += sc->s[i]->ntot;
714 sc->ntot += sc->r[i].end - sc->r[i].start;
721 /* find the barsamples_t associated with a lambda that corresponds to
722 a specific foreign lambda */
723 static sample_coll_t *lambda_data_find_sample_coll(lambda_data_t *l,
724 lambda_vec_t *foreign_lambda)
726 sample_coll_t *sc=l->sc->next;
730 if (lambda_vec_same(sc->foreign_lambda,foreign_lambda))
740 /* insert li into an ordered list of lambda_colls */
741 static void lambda_data_insert_sample_coll(lambda_data_t *l, sample_coll_t *sc)
743 sample_coll_t *scn=l->sc->next;
744 while ( (scn!=l->sc) )
746 if (lambda_vec_cmp_foreign(scn->foreign_lambda, sc->foreign_lambda) > 0)
750 /* now insert it before the found scn */
757 /* insert li into an ordered list of lambdas */
758 static void lambda_data_insert_lambda(lambda_data_t *head, lambda_data_t *li)
760 lambda_data_t *lc=head->next;
763 if (lambda_vec_cmp_native(lc->lambda, li->lambda) > 0)
767 /* now insert ourselves before the found lc */
774 /* insert a sample and a sample_range into a sample_coll. The
775 samples are stored as a pointer, the range is copied. */
776 static void sample_coll_insert_sample(sample_coll_t *sc, samples_t *s,
779 /* first check if it belongs here */
780 if (sc->temp != s->temp)
782 gmx_fatal(FARGS, "Temperatures in files %s and %s are not the same!",
783 s->filename, sc->next->s[0]->filename);
785 if (!lambda_vec_same(sc->native_lambda, s->native_lambda))
787 gmx_fatal(FARGS, "Native lambda in files %s and %s are not the same (and they should be)!",
788 s->filename, sc->next->s[0]->filename);
790 if (!lambda_vec_same(sc->foreign_lambda,s->foreign_lambda))
792 gmx_fatal(FARGS, "Foreign lambda in files %s and %s are not the same (and they should be)!",
793 s->filename, sc->next->s[0]->filename);
796 /* check if there's room */
797 if ( (sc->nsamples + 1) > sc->nsamples_alloc)
799 sc->nsamples_alloc = max(2*sc->nsamples_alloc, 2);
800 srenew(sc->s, sc->nsamples_alloc);
801 srenew(sc->r, sc->nsamples_alloc);
803 sc->s[sc->nsamples]=s;
804 sc->r[sc->nsamples]=*r;
807 sample_coll_calc_ntot(sc);
810 /* insert a sample into a lambda_list, creating the right sample_coll if
812 static void lambda_data_list_insert_sample(lambda_data_t *head, samples_t *s)
814 gmx_bool found=FALSE;
818 lambda_data_t *l=head->next;
820 /* first search for the right lambda_data_t */
823 if (lambda_vec_same(l->lambda, s->native_lambda) )
833 snew(l, 1); /* allocate a new one */
834 lambda_data_init(l, s->native_lambda, s->temp); /* initialize it */
835 lambda_data_insert_lambda(head, l); /* add it to the list */
838 /* now look for a sample collection */
839 sc=lambda_data_find_sample_coll(l, s->foreign_lambda);
842 snew(sc, 1); /* allocate a new one */
843 sample_coll_init(sc, s->native_lambda, s->foreign_lambda, s->temp);
844 lambda_data_insert_sample_coll(l, sc);
847 /* now insert the samples into the sample coll */
848 sample_range_init(&r, s);
849 sample_coll_insert_sample(sc, s, &r);
853 /* make a histogram out of a sample collection */
854 static void sample_coll_make_hist(sample_coll_t *sc, int **bin,
855 int *nbin_alloc, int *nbin,
856 double *dx, double *xmin, int nbin_default)
859 gmx_bool dx_set=FALSE;
860 gmx_bool xmin_set=FALSE;
862 gmx_bool xmax_set=FALSE;
863 gmx_bool xmax_set_hard=FALSE; /* whether the xmax is bounded by the
864 limits of a histogram */
867 /* first determine dx and xmin; try the histograms */
868 for(i=0;i<sc->nsamples;i++)
872 hist_t *hist=sc->s[i]->hist;
873 for(k=0;k<hist->nhist;k++)
875 double hdx=hist->dx[k];
876 double xmax_now=(hist->x0[k]+hist->nbin[k])*hdx;
878 /* we use the biggest dx*/
879 if ( (!dx_set) || hist->dx[0] > *dx)
884 if ( (!xmin_set) || (hist->x0[k]*hdx) < *xmin)
887 *xmin = (hist->x0[k]*hdx);
890 if ( (!xmax_set) || (xmax_now>xmax && !xmax_set_hard) )
894 if (hist->bin[k][hist->nbin[k]-1] != 0)
897 if ( hist->bin[k][hist->nbin[k]-1]!=0 && (xmax_now < xmax) )
905 /* and the delta us */
906 for(i=0;i<sc->nsamples;i++)
908 if (sc->s[i]->ndu > 0)
910 /* determine min and max */
911 int starti=sc->r[i].start;
912 int endi=sc->r[i].end;
913 double du_xmin=sc->s[i]->du[starti];
914 double du_xmax=sc->s[i]->du[starti];
915 for(j=starti+1;j<endi;j++)
917 if (sc->s[i]->du[j] < du_xmin)
918 du_xmin = sc->s[i]->du[j];
919 if (sc->s[i]->du[j] > du_xmax)
920 du_xmax = sc->s[i]->du[j];
923 /* and now change the limits */
924 if ( (!xmin_set) || (du_xmin < *xmin) )
929 if ( (!xmax_set) || ((du_xmax > xmax) && !xmax_set_hard) )
937 if (!xmax_set || !xmin_set)
947 *dx=(xmax-(*xmin))/((*nbin)-2); /* -2 because we want the last bin to
948 be 0, and we count from 0 */
952 *nbin=(xmax-(*xmin))/(*dx);
955 if (*nbin > *nbin_alloc)
958 srenew(*bin, *nbin_alloc);
961 /* reset the histogram */
962 for(i=0;i<(*nbin);i++)
967 /* now add the actual data */
968 for(i=0;i<sc->nsamples;i++)
972 hist_t *hist=sc->s[i]->hist;
973 for(k=0;k<hist->nhist;k++)
975 double hdx = hist->dx[k];
976 double xmin_hist=hist->x0[k]*hdx;
977 for(j=0;j<hist->nbin[k];j++)
979 /* calculate the bin corresponding to the middle of the
981 double x=hdx*(j+0.5) + xmin_hist;
982 int binnr=(int)((x-(*xmin))/(*dx));
984 if (binnr >= *nbin || binnr<0)
987 (*bin)[binnr] += hist->bin[k][j];
993 int starti=sc->r[i].start;
994 int endi=sc->r[i].end;
995 for(j=starti;j<endi;j++)
997 int binnr=(int)((sc->s[i]->du[j]-(*xmin))/(*dx));
998 if (binnr >= *nbin || binnr<0)
1007 /* write a collection of histograms to a file */
1008 void sim_data_histogram(sim_data_t *sd, const char *filename,
1009 int nbin_default, const output_env_t oenv)
1011 char label_x[STRLEN];
1012 const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
1013 const char *title="N(\\DeltaH)";
1014 const char *label_y="Samples";
1018 char **setnames=NULL;
1019 gmx_bool first_set=FALSE;
1020 /* histogram data: */
1027 lambda_data_t *bl_head = sd->lb;
1029 printf("\nWriting histogram to %s\n", filename);
1030 sprintf(label_x, "\\DeltaH (%s)", unit_energy);
1032 fp=xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1034 /* first get all the set names */
1036 /* iterate over all lambdas */
1039 sample_coll_t *sc=bl->sc->next;
1041 /* iterate over all samples */
1044 char buf[STRLEN], buf2[STRLEN];
1047 srenew(setnames, nsets);
1048 snew(setnames[nsets-1], STRLEN);
1049 if (sc->foreign_lambda->dhdl < 0)
1051 lambda_vec_print(sc->native_lambda, buf, FALSE);
1052 lambda_vec_print(sc->foreign_lambda, buf2, FALSE);
1053 sprintf(setnames[nsets-1], "N(%s(%s=%s) | %s=%s)",
1054 deltag, lambda, buf2, lambda, buf);
1058 lambda_vec_print(sc->native_lambda, buf, FALSE);
1059 sprintf(setnames[nsets-1], "N(%s | %s=%s)",
1067 xvgr_legend(fp, nsets, (const char**)setnames, oenv);
1070 /* now make the histograms */
1072 /* iterate over all lambdas */
1075 sample_coll_t *sc=bl->sc->next;
1077 /* iterate over all samples */
1082 xvgr_new_dataset(fp, 0, 0, NULL, oenv);
1085 sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &min,
1090 double xmin=i*dx + min;
1091 double xmax=(i+1)*dx + min;
1093 fprintf(fp, "%g %d\n%g %d\n", xmin, hist[i], xmax, hist[i]);
1109 /* create a collection (array) of barres_t object given a ordered linked list
1110 of barlamda_t sample collections */
1111 static barres_t *barres_list_create(sim_data_t *sd, int *nres,
1118 gmx_bool dhdl=FALSE;
1119 gmx_bool first=TRUE;
1120 lambda_data_t *bl_head=sd->lb;
1122 /* first count the lambdas */
1129 snew(res, nlambda-1);
1131 /* next put the right samples in the res */
1133 bl=bl_head->next->next; /* we start with the second one. */
1136 sample_coll_t *sc, *scprev;
1137 barres_t *br=&(res[*nres]);
1138 /* there is always a previous one. we search for that as a foreign
1140 scprev=lambda_data_find_sample_coll(bl->prev, bl->lambda);
1141 sc=lambda_data_find_sample_coll(bl, bl->prev->lambda);
1149 scprev=lambda_data_find_sample_coll(bl->prev, bl->prev->lambda);
1150 sc=lambda_data_find_sample_coll(bl, bl->lambda);
1154 printf("\nWARNING: Using the derivative data (dH/dlambda) to extrapolate delta H values.\nThis will only work if the Hamiltonian is linear in lambda.\n");
1159 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");
1162 else if (!scprev && !sc)
1164 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);
1167 /* normal delta H */
1170 gmx_fatal(FARGS,"Could not find a set for foreign lambda = %g\nin the files for lambda = %g",bl->lambda,bl->prev->lambda);
1174 gmx_fatal(FARGS,"Could not find a set for foreign lambda = %g\nin the files for lambda = %g",bl->prev->lambda,bl->lambda);
1186 /* estimate the maximum discretization error */
1187 static double barres_list_max_disc_err(barres_t *res, int nres)
1191 double delta_lambda;
1195 barres_t *br=&(res[i]);
1197 delta_lambda=lambda_vec_abs_diff(br->b->native_lambda,
1198 br->a->native_lambda);
1200 for(j=0;j<br->a->nsamples;j++)
1202 if (br->a->s[j]->hist)
1205 if (br->a->s[j]->derivative)
1206 Wfac = delta_lambda;
1208 disc_err=max(disc_err, Wfac*br->a->s[j]->hist->dx[0]);
1211 for(j=0;j<br->b->nsamples;j++)
1213 if (br->b->s[j]->hist)
1216 if (br->b->s[j]->derivative)
1217 Wfac = delta_lambda;
1218 disc_err=max(disc_err, Wfac*br->b->s[j]->hist->dx[0]);
1226 /* impose start and end times on a sample collection, updating sample_ranges */
1227 static void sample_coll_impose_times(sample_coll_t *sc, double begin_t,
1231 for(i=0;i<sc->nsamples;i++)
1233 samples_t *s=sc->s[i];
1234 sample_range_t *r=&(sc->r[i]);
1237 double end_time=s->hist->delta_time*s->hist->sum +
1238 s->hist->start_time;
1239 if (s->hist->start_time < begin_t || end_time > end_t)
1249 if (s->start_time < begin_t)
1251 r->start=(int)((begin_t - s->start_time)/s->delta_time);
1253 end_time=s->delta_time*s->ndu + s->start_time;
1254 if (end_time > end_t)
1256 r->end=(int)((end_t - s->start_time)/s->delta_time);
1262 for(j=0;j<s->ndu;j++)
1264 if (s->t[j] <begin_t)
1269 if (s->t[j] >= end_t)
1276 if (r->start > r->end)
1282 sample_coll_calc_ntot(sc);
1285 static void sim_data_impose_times(sim_data_t *sd, double begin, double end)
1287 double first_t, last_t;
1288 double begin_t, end_t;
1290 lambda_data_t *head=sd->lb;
1293 if (begin<=0 && end<0)
1298 /* first determine the global start and end times */
1304 sample_coll_t *sc=lc->sc->next;
1307 for(j=0;j<sc->nsamples;j++)
1309 double start_t,end_t;
1311 start_t = sc->s[j]->start_time;
1312 end_t = sc->s[j]->start_time;
1315 end_t += sc->s[j]->delta_time*sc->s[j]->hist->sum;
1321 end_t = sc->s[j]->t[sc->s[j]->ndu-1];
1325 end_t += sc->s[j]->delta_time*sc->s[j]->ndu;
1329 if (start_t < first_t || first_t<0)
1343 /* calculate the actual times */
1361 printf("\n Samples in time interval: %.3f - %.3f\n", first_t, last_t);
1363 if (begin_t > end_t)
1367 printf("Removing samples outside of: %.3f - %.3f\n", begin_t, end_t);
1369 /* then impose them */
1373 sample_coll_t *sc=lc->sc->next;
1376 sample_coll_impose_times(sc, begin_t, end_t);
1384 /* create subsample i out of ni from an existing sample_coll */
1385 static gmx_bool sample_coll_create_subsample(sample_coll_t *sc,
1386 sample_coll_t *sc_orig,
1390 int hist_start, hist_end;
1392 gmx_large_int_t ntot_start;
1393 gmx_large_int_t ntot_end;
1394 gmx_large_int_t ntot_so_far;
1396 *sc = *sc_orig; /* just copy all fields */
1398 /* allocate proprietary memory */
1399 snew(sc->s, sc_orig->nsamples);
1400 snew(sc->r, sc_orig->nsamples);
1402 /* copy the samples */
1403 for(j=0;j<sc_orig->nsamples;j++)
1405 sc->s[j] = sc_orig->s[j];
1406 sc->r[j] = sc_orig->r[j]; /* copy the ranges too */
1409 /* now fix start and end fields */
1410 /* the casts avoid possible overflows */
1411 ntot_start=(gmx_large_int_t)(sc_orig->ntot*(double)i/(double)ni);
1412 ntot_end =(gmx_large_int_t)(sc_orig->ntot*(double)(i+1)/(double)ni);
1414 for(j=0;j<sc->nsamples;j++)
1416 gmx_large_int_t ntot_add;
1417 gmx_large_int_t new_start, new_end;
1423 ntot_add = sc->s[j]->hist->sum;
1427 ntot_add = sc->r[j].end - sc->r[j].start;
1435 if (!sc->s[j]->hist)
1437 if (ntot_so_far < ntot_start)
1439 /* adjust starting point */
1440 new_start = sc->r[j].start + (ntot_start - ntot_so_far);
1444 new_start = sc->r[j].start;
1446 /* adjust end point */
1447 new_end = sc->r[j].start + (ntot_end - ntot_so_far);
1448 if (new_end > sc->r[j].end)
1450 new_end=sc->r[j].end;
1453 /* check if we're in range at all */
1454 if ( (new_end < new_start) || (new_start > sc->r[j].end) )
1459 /* and write the new range */
1460 sc->r[j].start=(int)new_start;
1461 sc->r[j].end=(int)new_end;
1468 double ntot_start_norm, ntot_end_norm;
1469 /* calculate the amount of overlap of the
1470 desired range (ntot_start -- ntot_end) onto
1471 the histogram range (ntot_so_far -- ntot_so_far+ntot_add)*/
1473 /* first calculate normalized bounds
1474 (where 0 is the start of the hist range, and 1 the end) */
1475 ntot_start_norm = (ntot_start-ntot_so_far)/(double)ntot_add;
1476 ntot_end_norm = (ntot_end-ntot_so_far)/(double)ntot_add;
1478 /* now fix the boundaries */
1479 ntot_start_norm = min(1, max(0., ntot_start_norm));
1480 ntot_end_norm = max(0, min(1., ntot_end_norm));
1482 /* and calculate the overlap */
1483 overlap = ntot_end_norm - ntot_start_norm;
1485 if (overlap > 0.95) /* we allow for 5% slack */
1487 sc->r[j].use = TRUE;
1489 else if (overlap < 0.05)
1491 sc->r[j].use = FALSE;
1499 ntot_so_far += ntot_add;
1501 sample_coll_calc_ntot(sc);
1506 /* calculate minimum and maximum work values in sample collection */
1507 static void sample_coll_min_max(sample_coll_t *sc, double Wfac,
1508 double *Wmin, double *Wmax)
1515 for(i=0;i<sc->nsamples;i++)
1517 samples_t *s=sc->s[i];
1518 sample_range_t *r=&(sc->r[i]);
1523 for(j=r->start; j<r->end; j++)
1525 *Wmin = min(*Wmin,s->du[j]*Wfac);
1526 *Wmax = max(*Wmax,s->du[j]*Wfac);
1531 int hd=0; /* determine the histogram direction: */
1533 if ( (s->hist->nhist>1) && (Wfac<0) )
1539 for(j=s->hist->nbin[hd]-1;j>=0;j--)
1541 *Wmin=min(*Wmin,Wfac*(s->hist->x0[hd])*dx);
1542 *Wmax=max(*Wmax,Wfac*(s->hist->x0[hd])*dx);
1543 /* look for the highest value bin with values */
1544 if (s->hist->bin[hd][j]>0)
1546 *Wmin=min(*Wmin,Wfac*(j+s->hist->x0[hd]+1)*dx);
1547 *Wmax=max(*Wmax,Wfac*(j+s->hist->x0[hd]+1)*dx);
1556 /* Initialize a sim_data structure */
1557 static void sim_data_init(sim_data_t *sd)
1559 /* make linked list */
1560 sd->lb = &(sd->lb_head);
1561 sd->lb->next = sd->lb;
1562 sd->lb->prev = sd->lb;
1564 lambda_components_init(&(sd->lc));
1568 static double calc_bar_sum(int n,const double *W,double Wfac,double sbMmDG)
1577 sum += 1./(1. + exp(Wfac*W[i] + sbMmDG));
1583 /* calculate the BAR average given a histogram
1585 if type== 0, calculate the best estimate for the average,
1586 if type==-1, calculate the minimum possible value given the histogram
1587 if type== 1, calculate the maximum possible value given the histogram */
1588 static double calc_bar_sum_hist(const hist_t *hist, double Wfac, double sbMmDG,
1594 /* normalization factor multiplied with bin width and
1595 number of samples (we normalize through M): */
1597 int hd=0; /* determine the histogram direction: */
1600 if ( (hist->nhist>1) && (Wfac<0) )
1605 max=hist->nbin[hd]-1;
1608 max=hist->nbin[hd]; /* we also add whatever was out of range */
1613 double x=Wfac*((i+hist->x0[hd])+0.5)*dx; /* bin middle */
1614 double pxdx=hist->bin[0][i]*normdx; /* p(x)dx */
1616 sum += pxdx/(1. + exp(x + sbMmDG));
1622 static double calc_bar_lowlevel(sample_coll_t *ca, sample_coll_t *cb,
1623 double temp, double tol, int type)
1628 double Wfac1,Wfac2,Wmin,Wmax;
1629 double DG0,DG1,DG2,dDG1;
1631 double n1, n2; /* numbers of samples as doubles */
1636 /* count the numbers of samples */
1642 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1643 if (ca->foreign_lambda->dhdl<0)
1645 /* this is the case when the delta U were calculated directly
1646 (i.e. we're not scaling dhdl) */
1652 /* we're using dhdl, so delta_lambda needs to be a
1653 multiplication factor. */
1654 /*double delta_lambda=cb->native_lambda-ca->native_lambda;*/
1655 double delta_lambda=lambda_vec_abs_diff(cb->native_lambda,
1657 if (cb->native_lambda->lc->N > 1)
1660 "Can't (yet) do multi-component dhdl interpolation");
1663 Wfac1 = beta*delta_lambda;
1664 Wfac2 = -beta*delta_lambda;
1669 /* We print the output both in kT and kJ/mol.
1670 * Here we determine DG in kT, so when beta < 1
1671 * the precision has to be increased.
1676 /* Calculate minimum and maximum work to give an initial estimate of
1677 * delta G as their average.
1680 double Wmin1, Wmin2, Wmax1, Wmax2;
1681 sample_coll_min_max(ca, Wfac1, &Wmin1, &Wmax1);
1682 sample_coll_min_max(cb, Wfac2, &Wmin2, &Wmax2);
1684 Wmin=min(Wmin1, Wmin2);
1685 Wmax=max(Wmax1, Wmax2);
1693 fprintf(debug,"DG %9.5f %9.5f\n",DG0,DG2);
1695 /* We approximate by bisection: given our initial estimates
1696 we keep checking whether the halfway point is greater or
1697 smaller than what we get out of the BAR averages.
1699 For the comparison we can use twice the tolerance. */
1700 while (DG2 - DG0 > 2*tol)
1702 DG1 = 0.5*(DG0 + DG2);
1704 /*printf("Wfac1=%g, Wfac2=%g, beta=%g, DG1=%g\n",Wfac1,Wfac2,beta,
1707 /* calculate the BAR averages */
1710 for(i=0;i<ca->nsamples;i++)
1712 samples_t *s=ca->s[i];
1713 sample_range_t *r=&(ca->r[i]);
1718 dDG1 += calc_bar_sum_hist(s->hist, Wfac1, (M-DG1), type);
1722 dDG1 += calc_bar_sum(r->end - r->start, s->du + r->start,
1727 for(i=0;i<cb->nsamples;i++)
1729 samples_t *s=cb->s[i];
1730 sample_range_t *r=&(cb->r[i]);
1735 dDG1 -= calc_bar_sum_hist(s->hist, Wfac2, -(M-DG1), type);
1739 dDG1 -= calc_bar_sum(r->end - r->start, s->du + r->start,
1755 fprintf(debug,"DG %9.5f %9.5f\n",DG0,DG2);
1759 return 0.5*(DG0 + DG2);
1762 static void calc_rel_entropy(sample_coll_t *ca, sample_coll_t *cb,
1763 double temp, double dg, double *sa, double *sb)
1769 double Wfac1, Wfac2;
1775 /* count the numbers of samples */
1779 /* to ensure the work values are the same as during the delta_G */
1780 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1781 if (ca->foreign_lambda->dhdl<0)
1783 /* this is the case when the delta U were calculated directly
1784 (i.e. we're not scaling dhdl) */
1790 /* we're using dhdl, so delta_lambda needs to be a
1791 multiplication factor. */
1792 double delta_lambda=lambda_vec_abs_diff(cb->native_lambda,
1794 Wfac1 = beta*delta_lambda;
1795 Wfac2 = -beta*delta_lambda;
1798 /* first calculate the average work in both directions */
1799 for(i=0;i<ca->nsamples;i++)
1801 samples_t *s=ca->s[i];
1802 sample_range_t *r=&(ca->r[i]);
1807 for(j=r->start;j<r->end;j++)
1808 W_ab += Wfac1*s->du[j];
1812 /* normalization factor multiplied with bin width and
1813 number of samples (we normalize through M): */
1816 int hd=0; /* histogram direction */
1817 if ( (s->hist->nhist>1) && (Wfac1<0) )
1823 for(j=0;j<s->hist->nbin[0];j++)
1825 double x=Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1826 double pxdx=s->hist->bin[0][j]*normdx; /* p(x)dx */
1834 for(i=0;i<cb->nsamples;i++)
1836 samples_t *s=cb->s[i];
1837 sample_range_t *r=&(cb->r[i]);
1842 for(j=r->start;j<r->end;j++)
1843 W_ba += Wfac1*s->du[j];
1847 /* normalization factor multiplied with bin width and
1848 number of samples (we normalize through M): */
1851 int hd=0; /* histogram direction */
1852 if ( (s->hist->nhist>1) && (Wfac2<0) )
1858 for(j=0;j<s->hist->nbin[0];j++)
1860 double x=Wfac1*((j+s->hist->x0[0])+0.5)*dx;/*bin ctr*/
1861 double pxdx=s->hist->bin[0][j]*normdx; /* p(x)dx */
1869 /* then calculate the relative entropies */
1874 static void calc_dg_stddev(sample_coll_t *ca, sample_coll_t *cb,
1875 double temp, double dg, double *stddev)
1879 double sigmafact=0.;
1881 double Wfac1, Wfac2;
1887 /* count the numbers of samples */
1891 /* to ensure the work values are the same as during the delta_G */
1892 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1893 if (ca->foreign_lambda->dhdl<0)
1895 /* this is the case when the delta U were calculated directly
1896 (i.e. we're not scaling dhdl) */
1902 /* we're using dhdl, so delta_lambda needs to be a
1903 multiplication factor. */
1904 double delta_lambda=lambda_vec_abs_diff(cb->native_lambda,
1906 Wfac1 = beta*delta_lambda;
1907 Wfac2 = -beta*delta_lambda;
1913 /* calculate average in both directions */
1914 for(i=0;i<ca->nsamples;i++)
1916 samples_t *s=ca->s[i];
1917 sample_range_t *r=&(ca->r[i]);
1922 for(j=r->start;j<r->end;j++)
1924 sigmafact += 1./(2. + 2.*cosh((M + Wfac1*s->du[j] - dg)));
1929 /* normalization factor multiplied with bin width and
1930 number of samples (we normalize through M): */
1933 int hd=0; /* histogram direction */
1934 if ( (s->hist->nhist>1) && (Wfac1<0) )
1940 for(j=0;j<s->hist->nbin[0];j++)
1942 double x=Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1943 double pxdx=s->hist->bin[0][j]*normdx; /* p(x)dx */
1945 sigmafact += pxdx/(2. + 2.*cosh((M + x - dg)));
1950 for(i=0;i<cb->nsamples;i++)
1952 samples_t *s=cb->s[i];
1953 sample_range_t *r=&(cb->r[i]);
1958 for(j=r->start;j<r->end;j++)
1960 sigmafact += 1./(2. + 2.*cosh((M - Wfac2*s->du[j] - dg)));
1965 /* normalization factor multiplied with bin width and
1966 number of samples (we normalize through M): */
1969 int hd=0; /* histogram direction */
1970 if ( (s->hist->nhist>1) && (Wfac2<0) )
1976 for(j=0;j<s->hist->nbin[0];j++)
1978 double x=Wfac2*((j+s->hist->x0[0])+0.5)*dx;/*bin ctr*/
1979 double pxdx=s->hist->bin[0][j]*normdx; /* p(x)dx */
1981 sigmafact += pxdx/(2. + 2.*cosh((M - x - dg)));
1987 sigmafact /= (n1 + n2);
1991 Shirts, Bair, Hooker & Pande, Phys. Rev. Lett 91, 140601 (2003): */
1992 *stddev = sqrt(((1./sigmafact) - ( (n1+n2)/n1 + (n1+n2)/n2 )));
1997 static void calc_bar(barres_t *br, double tol,
1998 int npee_min, int npee_max, gmx_bool *bEE,
2002 double dg_sig2,sa_sig2,sb_sig2,stddev_sig2; /* intermediate variance values
2003 for calculated quantities */
2004 int nsample1, nsample2;
2005 double temp=br->a->temp;
2007 double dg_min, dg_max;
2008 gmx_bool have_hist=FALSE;
2010 br->dg=calc_bar_lowlevel(br->a, br->b, temp, tol, 0);
2012 br->dg_disc_err = 0.;
2013 br->dg_histrange_err = 0.;
2015 /* check if there are histograms */
2016 for(i=0;i<br->a->nsamples;i++)
2018 if (br->a->r[i].use && br->a->s[i]->hist)
2026 for(i=0;i<br->b->nsamples;i++)
2028 if (br->b->r[i].use && br->b->s[i]->hist)
2036 /* calculate histogram-specific errors */
2039 dg_min=calc_bar_lowlevel(br->a, br->b, temp, tol, -1);
2040 dg_max=calc_bar_lowlevel(br->a, br->b, temp, tol, 1);
2042 if (fabs(dg_max - dg_min) > GMX_REAL_EPS*10)
2044 /* the histogram range error is the biggest of the differences
2045 between the best estimate and the extremes */
2046 br->dg_histrange_err = fabs(dg_max - dg_min);
2048 br->dg_disc_err = 0.;
2049 for(i=0;i<br->a->nsamples;i++)
2051 if (br->a->s[i]->hist)
2052 br->dg_disc_err=max(br->dg_disc_err, br->a->s[i]->hist->dx[0]);
2054 for(i=0;i<br->b->nsamples;i++)
2056 if (br->b->s[i]->hist)
2057 br->dg_disc_err=max(br->dg_disc_err, br->b->s[i]->hist->dx[0]);
2060 calc_rel_entropy(br->a, br->b, temp, br->dg, &(br->sa), &(br->sb));
2062 calc_dg_stddev(br->a, br->b, temp, br->dg, &(br->dg_stddev) );
2071 sample_coll_t ca, cb;
2073 /* initialize the samples */
2074 sample_coll_init(&ca, br->a->native_lambda, br->a->foreign_lambda,
2076 sample_coll_init(&cb, br->b->native_lambda, br->b->foreign_lambda,
2079 for(npee=npee_min; npee<=npee_max; npee++)
2088 double dstddev2 = 0;
2091 for(p=0; p<npee; p++)
2098 cac=sample_coll_create_subsample(&ca, br->a, p, npee);
2099 cbc=sample_coll_create_subsample(&cb, br->b, p, npee);
2103 printf("WARNING: histogram number incompatible with block number for averaging: can't do error estimate\n");
2106 sample_coll_destroy(&ca);
2108 sample_coll_destroy(&cb);
2112 dgp=calc_bar_lowlevel(&ca, &cb, temp, tol, 0);
2116 partsum[npee*(npee_max+1)+p] += dgp;
2118 calc_rel_entropy(&ca, &cb, temp, dgp, &sac, &sbc);
2123 calc_dg_stddev(&ca, &cb, temp, dgp, &stddevc );
2126 dstddev2 += stddevc*stddevc;
2128 sample_coll_destroy(&ca);
2129 sample_coll_destroy(&cb);
2133 dg_sig2 += (dgs2-dgs*dgs)/(npee-1);
2139 sa_sig2 += (dsa2-dsa*dsa)/(npee-1);
2140 sb_sig2 += (dsb2-dsb*dsb)/(npee-1);
2144 stddev_sig2 += (dstddev2-dstddev*dstddev)/(npee-1);
2146 br->dg_err = sqrt(dg_sig2/(npee_max - npee_min + 1));
2147 br->sa_err = sqrt(sa_sig2/(npee_max - npee_min + 1));
2148 br->sb_err = sqrt(sb_sig2/(npee_max - npee_min + 1));
2149 br->dg_stddev_err = sqrt(stddev_sig2/(npee_max - npee_min + 1));
2154 static double bar_err(int nbmin, int nbmax, const double *partsum)
2157 double svar,s,s2,dg;
2160 for(nb=nbmin; nb<=nbmax; nb++)
2166 dg = partsum[nb*(nbmax+1)+b];
2172 svar += (s2 - s*s)/(nb - 1);
2175 return sqrt(svar/(nbmax + 1 - nbmin));
2179 /* Seek the end of an identifier (consecutive non-spaces), followed by
2180 an optional number of spaces or '='-signs. Returns a pointer to the
2181 first non-space value found after that. Returns NULL if the string
2184 static const char *find_value(const char *str)
2186 gmx_bool name_end_found=FALSE;
2188 /* if the string is a NULL pointer, return a NULL pointer. */
2195 /* first find the end of the name */
2196 if (! name_end_found)
2198 if ( isspace(*str) || (*str == '=') )
2200 name_end_found=TRUE;
2205 if (! ( isspace(*str) || (*str == '=') ))
2217 /* read a vector-notation description of a lambda vector */
2218 static gmx_bool read_lambda_compvec(const char *str,
2220 const lambda_components_t *lc_in,
2221 lambda_components_t *lc_out,
2225 gmx_bool initialize_lc = FALSE; /* whether to initialize the lambda
2226 components, or to check them */
2227 gmx_bool start_reached=FALSE; /* whether the start of component names
2229 gmx_bool vector=FALSE; /* whether there are multiple components */
2230 int n = 0; /* current component number */
2231 const char *val_start=NULL; /* start of the component name, or NULL
2232 if not in a value */
2240 if (lc_out && lc_out->N == 0)
2242 initialize_lc = TRUE;
2263 else if (! isspace(*str))
2265 gmx_fatal(FARGS, "Error in lambda components in %s", fn);
2272 if (isspace(*str) || *str==')' || *str==',' || *str=='\0')
2279 lambda_components_add(lc_out, val_start,
2284 if (!lambda_components_check(lc_out, n, val_start,
2293 /* add a vector component to lv */
2294 lv->val[n] = strtod(val_start, &strtod_end);
2295 if (val_start == strtod_end)
2298 "Error reading lambda vector in %s", fn);
2301 /* reset for the next identifier */
2310 else if (isalnum(*str))
2321 gmx_fatal(FARGS, "Error in lambda components in %s", fn);
2329 else if (lv == NULL)
2335 gmx_fatal(FARGS, "Incomplete lambda vector data in %s",
2355 gmx_fatal(FARGS, "Incomplete lambda components data in %s", fn);
2361 /* read and check the component names from a string */
2362 static gmx_bool read_lambda_components(const char *str,
2363 lambda_components_t *lc,
2367 return read_lambda_compvec(str, NULL, NULL, lc, end, fn);
2370 /* read an initialized lambda vector from a string */
2371 static gmx_bool read_lambda_vector(const char *str,
2376 return read_lambda_compvec(str, lv, lv->lc, NULL, end, fn);
2381 /* deduce lambda value from legend.
2383 legend = the legend string
2385 lam = the initialized lambda vector
2386 returns whether to use the data in this set.
2388 static gmx_bool legend2lambda(const char *fn,
2394 const char *ptr=NULL, *ptr2=NULL;
2396 gmx_bool bdhdl=FALSE;
2397 const char *tostr=" to ";
2401 gmx_fatal(FARGS,"There is no legend in file '%s', can not deduce lambda",fn);
2404 /* look for the last 'to': */
2408 ptr2 = strstr(ptr2, tostr);
2415 while(ptr2!=NULL && *ptr2!= '\0' );
2419 ptr += strlen(tostr)-1; /* and advance past that 'to' */
2423 /* look for the = sign */
2424 ptr = strrchr(legend,'=');
2427 /* otherwise look for the last space */
2428 ptr = strrchr(legend,' ');
2432 if (strstr(legend,"dH"))
2437 else if (strchr(legend,'D') != NULL && strchr(legend,'H') != NULL)
2442 else /*if (strstr(legend, "pV"))*/
2453 gmx_fatal(FARGS,"There is no proper lambda legend in file '%s', can not deduce lambda",fn);
2457 ptr=find_value(ptr);
2458 if (!ptr || !read_lambda_vector(ptr, lam, NULL, fn))
2460 gmx_fatal(FARGS, "lambda vector '%s' %s faulty", legend, fn);
2469 ptr = strrchr(legend, '=');
2473 /* there must be a component name */
2477 gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2479 /* now backtrack to the start of the identifier */
2480 while(isspace(*ptr))
2486 gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2489 while(!isspace(*ptr))
2494 gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2498 strncpy(buf, ptr, (end-ptr));
2499 buf[(end-ptr)]='\0';
2500 dhdl_index=lambda_components_find(lam->lc, ptr, (end-ptr));
2504 strncpy(buf, ptr, (end-ptr));
2505 buf[(end-ptr)]='\0';
2507 "Did not find lambda component for '%s' in %s",
2516 "dhdl without component name with >1 lambda component in %s",
2521 lam->dhdl = dhdl_index;
2526 static gmx_bool subtitle2lambda(const char *subtitle, xvg_t *ba, const char *fn,
2527 lambda_components_t *lc)
2532 double native_lambda;
2536 /* first check for a state string */
2537 ptr = strstr(subtitle, "state");
2541 const char *val_end;
2543 /* the new 4.6 style lambda vectors */
2544 ptr = find_value(ptr);
2547 index = strtol(ptr, &end, 10);
2550 gmx_fatal(FARGS,"Incomplete state data in %s", fn);
2557 gmx_fatal(FARGS,"Incomplete state data in %s", fn);
2560 /* now find the lambda vector component names */
2561 while(*ptr!='(' && ! isalnum(*ptr))
2567 "Incomplete lambda vector component data in %s", fn);
2572 if (!read_lambda_components(ptr, lc, &val_end, fn))
2575 "lambda vector components in %s don't match those previously read",
2578 ptr=find_value(val_end);
2581 gmx_fatal(FARGS,"Incomplete state data in %s", fn);
2584 lambda_vec_init(&(ba->native_lambda), lc);
2585 if (!read_lambda_vector(ptr, &(ba->native_lambda), NULL, fn))
2587 gmx_fatal(FARGS, "lambda vector in %s faulty", fn);
2589 ba->native_lambda.index=index;
2594 /* compatibility mode: check for lambda in other ways. */
2595 /* plain text lambda string */
2596 ptr = strstr(subtitle,"lambda");
2599 /* xmgrace formatted lambda string */
2600 ptr = strstr(subtitle,"\\xl\\f{}");
2604 /* xmgr formatted lambda string */
2605 ptr = strstr(subtitle,"\\8l\\4");
2609 ptr = strstr(ptr,"=");
2613 bFound = (sscanf(ptr+1,"%lf",&(native_lambda)) == 1);
2614 /* add the lambda component name as an empty string */
2617 if (!lambda_components_check(lc, 0, "", 0))
2620 "lambda vector components in %s don't match those previously read",
2626 lambda_components_add(lc, "", 0);
2628 lambda_vec_init(&(ba->native_lambda), lc);
2629 ba->native_lambda.val[0]=native_lambda;
2636 static void filename2lambda(const char *fn, xvg_t *ba)
2639 const char *ptr,*digitptr;
2643 /* go to the end of the path string and search backward to find the last
2644 directory in the path which has to contain the value of lambda
2646 while (ptr[1] != '\0')
2650 /* searching backward to find the second directory separator */
2655 if (ptr[0] != DIR_SEPARATOR && ptr[1] == DIR_SEPARATOR)
2657 if (dirsep == 1) break;
2660 /* save the last position of a digit between the last two
2661 separators = in the last dirname */
2662 if (dirsep > 0 && isdigit(*ptr))
2670 gmx_fatal(FARGS,"While trying to read the lambda value from the file path:"
2671 " last directory in the path '%s' does not contain a number",fn);
2673 if (digitptr[-1] == '-')
2677 lambda = strtod(digitptr,&endptr);
2678 if (endptr == digitptr)
2680 gmx_fatal(FARGS,"Malformed number in file path '%s'",fn);
2684 static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
2685 lambda_components_t *lc)
2688 char *subtitle,**legend,*ptr;
2690 gmx_bool native_lambda_read=FALSE;
2698 np = read_xvg_legend(fn,&ba->y,&ba->nset,&subtitle,&legend);
2701 gmx_fatal(FARGS,"File %s contains no usable data.",fn);
2703 /* Reorder the data */
2705 for(i=1; i<ba->nset; i++)
2707 ba->y[i-1] = ba->y[i];
2711 snew(ba->np,ba->nset);
2712 for(i=0;i<ba->nset;i++)
2716 if (subtitle != NULL)
2718 /* try to extract temperature */
2719 ptr = strstr(subtitle,"T =");
2723 if (sscanf(ptr,"%lf",&ba->temp) == 1)
2727 gmx_fatal(FARGS,"Found temperature of %g in file '%s'",
2737 gmx_fatal(FARGS,"Did not find a temperature in the subtitle in file '%s', use the -temp option of [TT]g_bar[tt]",fn);
2742 /* Try to deduce lambda from the subtitle */
2745 if (subtitle2lambda(subtitle,ba, fn, lc))
2747 native_lambda_read=TRUE;
2750 snew(ba->lambda,ba->nset);
2753 /* Check if we have a single set, no legend, nset=1 means t and dH/dl */
2756 if (!native_lambda_read)
2758 /* Deduce lambda from the file name */
2759 filename2lambda(fn, ba);
2760 native_lambda_read=TRUE;
2762 ba->lambda[0] = ba->native_lambda;
2766 gmx_fatal(FARGS,"File %s contains multiple sets but no legends, can not determine the lambda values",fn);
2771 for(i=0; i<ba->nset; )
2774 /* Read lambda from the legend */
2775 lambda_vec_init( &(ba->lambda[i]), lc );
2776 lambda_vec_copy( &(ba->lambda[i]), &(ba->native_lambda));
2777 use=legend2lambda(fn,legend[i], ba, &(ba->lambda[i]));
2780 lambda_vec_print(&(ba->lambda[i]), buf, FALSE);
2786 printf("%s: Ignoring set '%s'.\n", fn, legend[i]);
2787 for(j=i+1; j<ba->nset; j++)
2789 ba->y[j-1] = ba->y[j];
2790 legend[j-1] = legend[j];
2797 if (!native_lambda_read)
2799 gmx_fatal(FARGS,"File %s contains multiple sets but no indication of the native lambda",fn);
2804 for(i=0; i<ba->nset-1; i++)
2812 static void read_bar_xvg(char *fn, real *temp, sim_data_t *sd)
2821 read_bar_xvg_lowlevel(fn, temp, barsim, &(sd->lc));
2823 if (barsim->nset <1 )
2825 gmx_fatal(FARGS,"File '%s' contains fewer than two columns", fn);
2828 if ( !gmx_within_tol(*temp,barsim->temp,GMX_FLOAT_EPS) && (*temp > 0) )
2830 gmx_fatal(FARGS,"Temperature in file %s different from earlier files or setting\n", fn);
2834 /* now create a series of samples_t */
2835 snew(s, barsim->nset);
2836 for(i=0;i<barsim->nset;i++)
2838 samples_init(s+i, &(barsim->native_lambda), &(barsim->lambda[i]),
2839 barsim->temp, lambda_vec_same(&(barsim->native_lambda),
2840 &(barsim->lambda[i])),
2842 s[i].du=barsim->y[i];
2843 s[i].ndu=barsim->np[i];
2846 lambda_data_list_insert_sample(sd->lb, s+i);
2851 lambda_vec_print(s[0].native_lambda, buf, FALSE);
2852 printf("%s: %.1f - %.1f; lambda = %s\n dH/dl & foreign lambdas:\n",
2853 fn, s[0].t[0], s[0].t[s[0].ndu-1], buf);
2854 for(i=0;i<barsim->nset;i++)
2856 lambda_vec_print(s[i].foreign_lambda, buf, TRUE);
2857 printf(" %s (%d pts)\n", buf, s[i].ndu);
2863 static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
2864 double start_time, double delta_time,
2865 lambda_vec_t *native_lambda, double temp,
2866 double *last_t, const char *filename)
2870 double old_foreign_lambda;
2871 lambda_vec_t *foreign_lambda;
2873 samples_t *s; /* convenience pointer */
2876 /* check the block types etc. */
2877 if ( (blk->nsub < 3) ||
2878 (blk->sub[0].type != xdr_datatype_int) ||
2879 (blk->sub[1].type != xdr_datatype_double) ||
2881 (blk->sub[2].type != xdr_datatype_float) &&
2882 (blk->sub[2].type != xdr_datatype_double)
2884 (blk->sub[0].nr < 1) ||
2885 (blk->sub[1].nr < 1) )
2888 "Unexpected/corrupted block data in file %s around time %g.",
2889 filename, start_time);
2892 snew(foreign_lambda, 1);
2893 lambda_vec_init(foreign_lambda, native_lambda->lc);
2894 lambda_vec_copy(foreign_lambda, native_lambda);
2895 type = blk->sub[0].ival[0];
2898 for(i=0;i<native_lambda->lc->N;i++)
2900 foreign_lambda->val[i] = blk->sub[1].dval[i];
2905 if (blk->sub[0].nr > 1)
2907 foreign_lambda->dhdl = blk->sub[0].ival[1];
2910 foreign_lambda->dhdl = 0;
2915 /* initialize the samples structure if it's empty. */
2917 samples_init(*smp, native_lambda, foreign_lambda, temp,
2918 type==dhbtDHDL, filename);
2919 (*smp)->start_time=start_time;
2920 (*smp)->delta_time=delta_time;
2923 /* set convenience pointer */
2926 /* now double check */
2927 if ( ! lambda_vec_same(s->foreign_lambda, foreign_lambda) )
2929 char buf[STRLEN], buf2[STRLEN];
2930 lambda_vec_print(foreign_lambda, buf, FALSE);
2931 lambda_vec_print(s->foreign_lambda, buf2, FALSE);
2932 fprintf(stderr, "Got foreign lambda=%s, expected: %s\n", buf, buf2);
2933 gmx_fatal(FARGS, "Corrupted data in file %s around t=%g.",
2934 filename, start_time);
2937 /* make room for the data */
2938 if (s->ndu_alloc < (size_t)(s->ndu + blk->sub[2].nr) )
2940 s->ndu_alloc += (s->ndu_alloc < (size_t)blk->sub[2].nr) ?
2941 blk->sub[2].nr*2 : s->ndu_alloc;
2942 srenew(s->du_alloc, s->ndu_alloc);
2946 s->ndu += blk->sub[2].nr;
2947 s->ntot += blk->sub[2].nr;
2948 *ndu = blk->sub[2].nr;
2950 /* and copy the data*/
2951 for(j=0;j<blk->sub[2].nr;j++)
2953 if (blk->sub[2].type == xdr_datatype_float)
2955 s->du[startj+j] = blk->sub[2].fval[j];
2959 s->du[startj+j] = blk->sub[2].dval[j];
2962 if (start_time + blk->sub[2].nr*delta_time > *last_t)
2964 *last_t = start_time + blk->sub[2].nr*delta_time;
2968 static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
2969 double start_time, double delta_time,
2970 lambda_vec_t *native_lambda, double temp,
2971 double *last_t, const char *filename)
2976 double old_foreign_lambda;
2977 lambda_vec_t *foreign_lambda;
2981 /* check the block types etc. */
2982 if ( (blk->nsub < 2) ||
2983 (blk->sub[0].type != xdr_datatype_double) ||
2984 (blk->sub[1].type != xdr_datatype_large_int) ||
2985 (blk->sub[0].nr < 2) ||
2986 (blk->sub[1].nr < 2) )
2989 "Unexpected/corrupted block data in file %s around time %g",
2990 filename, start_time);
3001 "Unexpected/corrupted block data in file %s around time %g",
3002 filename, start_time);
3008 snew(foreign_lambda, 1);
3009 lambda_vec_init(foreign_lambda, native_lambda->lc);
3010 lambda_vec_copy(foreign_lambda, native_lambda);
3011 type = (int)(blk->sub[1].lval[1]);
3014 double old_foreign_lambda;
3016 old_foreign_lambda=blk->sub[0].dval[0];
3017 if (old_foreign_lambda >= 0)
3019 foreign_lambda->val[0]=old_foreign_lambda;
3020 if (foreign_lambda->lc->N > 1)
3023 "Single-component lambda in multi-component file %s",
3029 for(i=0;i<native_lambda->lc->N;i++)
3031 foreign_lambda->val[i] = blk->sub[0].dval[i+2];
3037 if (foreign_lambda->lc->N > 1)
3039 if (blk->sub[1].nr < 3 + nhist)
3042 "Missing derivative coord in multi-component file %s",
3045 foreign_lambda->dhdl = blk->sub[1].lval[2 + nhist];
3049 foreign_lambda->dhdl = 0;
3053 samples_init(s, native_lambda, foreign_lambda, temp, type==dhbtDHDL,
3057 for(i=0;i<nhist;i++)
3059 nbins[i] = blk->sub[i+2].nr;
3062 hist_init(s->hist, nhist, nbins);
3064 for(i=0;i<nhist;i++)
3066 s->hist->x0[i] = blk->sub[1].lval[2+i];
3067 s->hist->dx[i] = blk->sub[0].dval[1];
3070 s->hist->dx[i] = - s->hist->dx[i];
3074 s->hist->start_time = start_time;
3075 s->hist->delta_time = delta_time;
3076 s->start_time = start_time;
3077 s->delta_time = delta_time;
3079 for(i=0;i<nhist;i++)
3082 gmx_large_int_t sum=0;
3084 for(j=0;j<s->hist->nbin[i];j++)
3086 int binv=(int)(blk->sub[i+2].ival[j]);
3088 s->hist->bin[i][j] = binv;
3101 gmx_fatal(FARGS, "Histogram counts don't match in %s",
3107 if (start_time + s->hist->sum*delta_time > *last_t)
3109 *last_t = start_time + s->hist->sum*delta_time;
3115 static void read_barsim_edr(char *fn, real *temp, sim_data_t *sd)
3121 gmx_enxnm_t *enm=NULL;
3124 samples_t **samples_rawdh=NULL; /* contains samples for raw delta_h */
3125 int *nhists=NULL; /* array to keep count & print at end */
3126 int *npts=NULL; /* array to keep count & print at end */
3127 lambda_vec_t **lambdas=NULL; /* array to keep count & print at end */
3128 lambda_vec_t *native_lambda;
3129 double end_time; /* the end time of the last batch of samples */
3131 lambda_vec_t start_lambda;
3133 fp = open_enx(fn,"r");
3134 do_enxnms(fp,&nre,&enm);
3137 snew(native_lambda, 1);
3138 start_lambda.lc = NULL;
3140 while(do_enx(fp, fr))
3142 /* count the data blocks */
3147 /* DHCOLL block information: */
3148 double start_time=0, delta_time=0, old_start_lambda=0, delta_lambda=0;
3151 /* count the blocks and handle collection information: */
3152 for(i=0;i<fr->nblock;i++)
3154 if (fr->block[i].id == enxDHHIST )
3156 if ( fr->block[i].id == enxDH )
3158 if (fr->block[i].id == enxDHCOLL )
3161 if ( (fr->block[i].nsub < 1) ||
3162 (fr->block[i].sub[0].type != xdr_datatype_double) ||
3163 (fr->block[i].sub[0].nr < 5))
3165 gmx_fatal(FARGS, "Unexpected block data in file %s", fn);
3168 /* read the data from the DHCOLL block */
3169 rtemp = fr->block[i].sub[0].dval[0];
3170 start_time = fr->block[i].sub[0].dval[1];
3171 delta_time = fr->block[i].sub[0].dval[2];
3172 old_start_lambda = fr->block[i].sub[0].dval[3];
3173 delta_lambda = fr->block[i].sub[0].dval[4];
3177 gmx_fatal(FARGS, "Lambda values not constant in %s: can't apply BAR method", fn);
3179 if ( ( *temp != rtemp) && (*temp > 0) )
3181 gmx_fatal(FARGS,"Temperature in file %s different from earlier files or setting\n", fn);
3185 if (old_start_lambda >= 0)
3189 if (!lambda_components_check(&(sd->lc), 0, "", 0))
3192 "lambda vector components in %s don't match those previously read",
3198 lambda_components_add(&(sd->lc), "", 0);
3200 if (!start_lambda.lc)
3202 lambda_vec_init(&start_lambda, &(sd->lc));
3204 start_lambda.val[0]=old_start_lambda;
3208 /* read lambda vector */
3210 gmx_bool check=(sd->lc.N > 0);
3211 if (fr->block[i].nsub < 2)
3214 "No lambda vector, but start_lambda=%g\n",
3217 n_lambda_vec=fr->block[i].sub[1].ival[1];
3218 for(j=0;j<n_lambda_vec;j++)
3221 efpt_singular_names[fr->block[i].sub[1].ival[1+j]];
3224 /* check the components */
3225 lambda_components_check(&(sd->lc), j, name,
3230 lambda_components_add(&(sd->lc), name,
3234 lambda_vec_init(&start_lambda, &(sd->lc));
3235 start_lambda.index=fr->block[i].sub[1].ival[0];
3236 for(j=0;j<n_lambda_vec;j++)
3238 start_lambda.val[j]=fr->block[i].sub[0].dval[5+j];
3248 gmx_fatal(FARGS, "Did not find delta H information in file %s", fn);
3250 if (nblocks_raw > 0 && nblocks_hist > 0 )
3252 gmx_fatal(FARGS, "Can't handle both raw delta U data and histograms in the same file %s", fn);
3257 /* check the native lambda */
3258 if (!lambda_vec_same(&start_lambda, native_lambda) )
3260 gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %g, and becomes %g at time %g",
3261 fn, native_lambda, start_lambda, start_time);
3263 /* check the number of samples against the previous number */
3264 if ( ((nblocks_raw+nblocks_hist)!=nsamples) || (nlam!=1 ) )
3266 gmx_fatal(FARGS, "Unexpected block count in %s: was %d, now %d\n",
3267 fn, nsamples+1, nblocks_raw+nblocks_hist+nlam);
3269 /* check whether last iterations's end time matches with
3270 the currrent start time */
3271 if ( (fabs(last_t - start_time) > 2*delta_time) && last_t>=0)
3273 /* it didn't. We need to store our samples and reallocate */
3274 for(i=0;i<nsamples;i++)
3276 if (samples_rawdh[i])
3278 /* insert it into the existing list */
3279 lambda_data_list_insert_sample(sd->lb,
3281 /* and make sure we'll allocate a new one this time
3283 samples_rawdh[i]=NULL;
3290 /* this is the first round; allocate the associated data
3292 /*native_lambda=start_lambda;*/
3293 lambda_vec_init(native_lambda, &(sd->lc));
3294 lambda_vec_copy(native_lambda, &start_lambda);
3295 nsamples=nblocks_raw+nblocks_hist;
3296 snew(nhists, nsamples);
3297 snew(npts, nsamples);
3298 snew(lambdas, nsamples);
3299 snew(samples_rawdh, nsamples);
3300 for(i=0;i<nsamples;i++)
3305 samples_rawdh[i]=NULL; /* init to NULL so we know which
3306 ones contain values */
3311 k=0; /* counter for the lambdas, etc. arrays */
3312 for(i=0;i<fr->nblock;i++)
3314 if (fr->block[i].id == enxDH)
3316 int type=(fr->block[i].sub[0].ival[0]);
3317 if (type == dhbtDH || type == dhbtDHDL)
3320 read_edr_rawdh_block(&(samples_rawdh[k]),
3323 start_time, delta_time,
3324 native_lambda, rtemp,
3327 if (samples_rawdh[k])
3329 lambdas[k]=samples_rawdh[k]->foreign_lambda;
3334 else if (fr->block[i].id == enxDHHIST)
3336 int type=(int)(fr->block[i].sub[1].lval[1]);
3337 if (type == dhbtDH || type == dhbtDHDL)
3341 samples_t *s; /* this is where the data will go */
3342 s=read_edr_hist_block(&nb, &(fr->block[i]),
3343 start_time, delta_time,
3344 native_lambda, rtemp,
3349 lambdas[k]=s->foreign_lambda;
3352 /* and insert the new sample immediately */
3355 lambda_data_list_insert_sample(sd->lb, s+j);
3361 /* Now store all our extant sample collections */
3362 for(i=0;i<nsamples;i++)
3364 if (samples_rawdh[i])
3366 /* insert it into the existing list */
3367 lambda_data_list_insert_sample(sd->lb, samples_rawdh[i]);
3375 lambda_vec_print(native_lambda, buf, FALSE);
3376 printf("%s: %.1f - %.1f; lambda = %s\n foreign lambdas:\n",
3377 fn, first_t, last_t, buf);
3378 for(i=0;i<nsamples;i++)
3382 lambda_vec_print(lambdas[i], buf, TRUE);
3385 printf(" %s (%d hists)\n", buf, nhists[i]);
3389 printf(" %s (%d pts)\n", buf, npts[i]);
3401 int gmx_bar(int argc,char *argv[])
3403 static const char *desc[] = {
3404 "[TT]g_bar[tt] calculates free energy difference estimates through ",
3405 "Bennett's acceptance ratio method (BAR). It also automatically",
3406 "adds series of individual free energies obtained with BAR into",
3407 "a combined free energy estimate.[PAR]",
3409 "Every individual BAR free energy difference relies on two ",
3410 "simulations at different states: say state A and state B, as",
3411 "controlled by a parameter, [GRK]lambda[grk] (see the [TT].mdp[tt] parameter",
3412 "[TT]init_lambda[tt]). The BAR method calculates a ratio of weighted",
3413 "average of the Hamiltonian difference of state B given state A and",
3415 "The energy differences to the other state must be calculated",
3416 "explicitly during the simulation. This can be done with",
3417 "the [TT].mdp[tt] option [TT]foreign_lambda[tt].[PAR]",
3419 "Input option [TT]-f[tt] expects multiple [TT]dhdl.xvg[tt] files. ",
3420 "Two types of input files are supported:[BR]",
3421 "[TT]*[tt] Files with more than one [IT]y[it]-value. ",
3422 "The files should have columns ",
3423 "with dH/d[GRK]lambda[grk] and [GRK]Delta[grk][GRK]lambda[grk]. ",
3424 "The [GRK]lambda[grk] values are inferred ",
3425 "from the legends: [GRK]lambda[grk] of the simulation from the legend of ",
3426 "dH/d[GRK]lambda[grk] and the foreign [GRK]lambda[grk] values from the ",
3427 "legends of Delta H",
3429 "[TT]*[tt] Files with only one [IT]y[it]-value. Using the",
3430 "[TT]-extp[tt] option for these files, it is assumed",
3431 "that the [IT]y[it]-value is dH/d[GRK]lambda[grk] and that the ",
3432 "Hamiltonian depends linearly on [GRK]lambda[grk]. ",
3433 "The [GRK]lambda[grk] value of the simulation is inferred from the ",
3434 "subtitle (if present), otherwise from a number in the subdirectory ",
3435 "in the file name.[PAR]",
3437 "The [GRK]lambda[grk] of the simulation is parsed from ",
3438 "[TT]dhdl.xvg[tt] file's legend containing the string 'dH', the ",
3439 "foreign [GRK]lambda[grk] values from the legend containing the ",
3440 "capitalized letters 'D' and 'H'. The temperature is parsed from ",
3441 "the legend line containing 'T ='.[PAR]",
3443 "The input option [TT]-g[tt] expects multiple [TT].edr[tt] files. ",
3444 "These can contain either lists of energy differences (see the ",
3445 "[TT].mdp[tt] option [TT]separate_dhdl_file[tt]), or a series of ",
3446 "histograms (see the [TT].mdp[tt] options [TT]dh_hist_size[tt] and ",
3447 "[TT]dh_hist_spacing[tt]).", "The temperature and [GRK]lambda[grk] ",
3448 "values are automatically deduced from the [TT]ener.edr[tt] file.[PAR]",
3450 "In addition to the [TT].mdp[tt] option [TT]foreign_lambda[tt], ",
3451 "the energy difference can also be extrapolated from the ",
3452 "dH/d[GRK]lambda[grk] values. This is done with the[TT]-extp[tt]",
3453 "option, which assumes that the system's Hamiltonian depends linearly",
3454 "on [GRK]lambda[grk], which is not normally the case.[PAR]",
3456 "The free energy estimates are determined using BAR with bisection, ",
3457 "with the precision of the output set with [TT]-prec[tt]. ",
3458 "An error estimate taking into account time correlations ",
3459 "is made by splitting the data into blocks and determining ",
3460 "the free energy differences over those blocks and assuming ",
3461 "the blocks are independent. ",
3462 "The final error estimate is determined from the average variance ",
3463 "over 5 blocks. A range of block numbers for error estimation can ",
3464 "be provided with the options [TT]-nbmin[tt] and [TT]-nbmax[tt].[PAR]",
3466 "[TT]g_bar[tt] tries to aggregate samples with the same 'native' and ",
3467 "'foreign' [GRK]lambda[grk] values, but always assumes independent ",
3468 "samples. [BB]Note[bb] that when aggregating energy ",
3469 "differences/derivatives with different sampling intervals, this is ",
3470 "almost certainly not correct. Usually subsequent energies are ",
3471 "correlated and different time intervals mean different degrees ",
3472 "of correlation between samples.[PAR]",
3474 "The results are split in two parts: the last part contains the final ",
3475 "results in kJ/mol, together with the error estimate for each part ",
3476 "and the total. The first part contains detailed free energy ",
3477 "difference estimates and phase space overlap measures in units of ",
3478 "kT (together with their computed error estimate). The printed ",
3480 "[TT]*[tt] lam_A: the [GRK]lambda[grk] values for point A.[BR]",
3481 "[TT]*[tt] lam_B: the [GRK]lambda[grk] values for point B.[BR]",
3482 "[TT]*[tt] DG: the free energy estimate.[BR]",
3483 "[TT]*[tt] s_A: an estimate of the relative entropy of B in A.[BR]",
3484 "[TT]*[tt] s_A: an estimate of the relative entropy of A in B.[BR]",
3485 "[TT]*[tt] stdev: an estimate expected per-sample standard deviation.[PAR]",
3487 "The relative entropy of both states in each other's ensemble can be ",
3488 "interpreted as a measure of phase space overlap: ",
3489 "the relative entropy s_A of the work samples of lambda_B in the ",
3490 "ensemble of lambda_A (and vice versa for s_B), is a ",
3491 "measure of the 'distance' between Boltzmann distributions of ",
3492 "the two states, that goes to zero for identical distributions. See ",
3493 "Wu & Kofke, J. Chem. Phys. 123 084109 (2005) for more information.",
3495 "The estimate of the expected per-sample standard deviation, as given ",
3496 "in Bennett's original BAR paper: Bennett, J. Comp. Phys. 22, p 245 (1976).",
3497 "Eq. 10 therein gives an estimate of the quality of sampling (not directly",
3498 "of the actual statistical error, because it assumes independent samples).[PAR]",
3500 "To get a visual estimate of the phase space overlap, use the ",
3501 "[TT]-oh[tt] option to write series of histograms, together with the ",
3502 "[TT]-nbin[tt] option.[PAR]"
3504 static real begin=0,end=-1,temp=-1;
3505 int nd=2,nbmin=5,nbmax=5;
3507 gmx_bool use_dhdl=FALSE;
3508 gmx_bool calc_s,calc_v;
3510 { "-b", FALSE, etREAL, {&begin}, "Begin time for BAR" },
3511 { "-e", FALSE, etREAL, {&end}, "End time for BAR" },
3512 { "-temp", FALSE, etREAL, {&temp}, "Temperature (K)" },
3513 { "-prec", FALSE, etINT, {&nd}, "The number of digits after the decimal point" },
3514 { "-nbmin", FALSE, etINT, {&nbmin}, "Minimum number of blocks for error estimation" },
3515 { "-nbmax", FALSE, etINT, {&nbmax}, "Maximum number of blocks for error estimation" },
3516 { "-nbin", FALSE, etINT, {&nbin}, "Number of bins for histogram output"},
3517 { "-extp", FALSE, etBOOL, {&use_dhdl}, "Whether to linearly extrapolate dH/dl values to use as energies"}
3521 { efXVG, "-f", "dhdl", ffOPTRDMULT },
3522 { efEDR, "-g", "ener", ffOPTRDMULT },
3523 { efXVG, "-o", "bar", ffOPTWR },
3524 { efXVG, "-oi", "barint", ffOPTWR },
3525 { efXVG, "-oh", "histogram", ffOPTWR }
3527 #define NFILE asize(fnm)
3530 int nf=0; /* file counter */
3532 int nfile_tot; /* total number of input files */
3537 sim_data_t sim_data; /* the simulation data */
3538 barres_t *results; /* the results */
3539 int nresults; /* number of results in results array */
3542 double prec,dg_tot,dg,sig, dg_tot_max, dg_tot_min;
3544 char dgformat[20],xvg2format[STRLEN],xvg3format[STRLEN];
3545 char buf[STRLEN], buf2[STRLEN];
3546 char ktformat[STRLEN], sktformat[STRLEN];
3547 char kteformat[STRLEN], skteformat[STRLEN];
3550 gmx_bool result_OK=TRUE,bEE=TRUE;
3552 gmx_bool disc_err=FALSE;
3553 double sum_disc_err=0.; /* discretization error */
3554 gmx_bool histrange_err=FALSE;
3555 double sum_histrange_err=0.; /* histogram range error */
3556 double stat_err=0.; /* statistical error */
3558 CopyRight(stderr,argv[0]);
3559 parse_common_args(&argc,argv,
3561 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
3563 if (opt2bSet("-f", NFILE, fnm))
3565 nxvgfile = opt2fns(&fxvgnms,"-f",NFILE,fnm);
3567 if (opt2bSet("-g", NFILE, fnm))
3569 nedrfile = opt2fns(&fedrnms,"-g",NFILE,fnm);
3572 sim_data_init(&sim_data);
3574 /* make linked list */
3576 lambda_data_init(lb, 0, 0);
3582 nfile_tot = nxvgfile + nedrfile;
3586 gmx_fatal(FARGS,"No input files!");
3591 gmx_fatal(FARGS,"Can not have negative number of digits");
3595 snew(partsum,(nbmax+1)*(nbmax+1));
3598 /* read in all files. First xvg files */
3599 for(f=0; f<nxvgfile; f++)
3601 read_bar_xvg(fxvgnms[f],&temp,&sim_data);
3604 /* then .edr files */
3605 for(f=0; f<nedrfile; f++)
3607 read_barsim_edr(fedrnms[f],&temp,&sim_data);;
3611 /* fix the times to allow for equilibration */
3612 sim_data_impose_times(&sim_data, begin, end);
3614 if (opt2bSet("-oh",NFILE,fnm))
3616 sim_data_histogram(&sim_data, opt2fn("-oh",NFILE,fnm), nbin, oenv);
3619 /* assemble the output structures from the lambdas */
3620 results=barres_list_create(&sim_data, &nresults, use_dhdl);
3622 sum_disc_err=barres_list_max_disc_err(results, nresults);
3626 printf("\nNo results to calculate.\n");
3630 if (sum_disc_err > prec)
3633 nd = ceil(-log10(prec));
3634 printf("WARNING: setting the precision to %g because that is the minimum\n reasonable number, given the expected discretization error.\n", prec);
3638 /*sprintf(lamformat,"%%6.3f");*/
3639 sprintf( dgformat,"%%%d.%df",3+nd,nd);
3640 /* the format strings of the results in kT */
3641 sprintf( ktformat,"%%%d.%df",5+nd,nd);
3642 sprintf( sktformat,"%%%ds",6+nd);
3643 /* the format strings of the errors in kT */
3644 sprintf( kteformat,"%%%d.%df",3+nd,nd);
3645 sprintf( skteformat,"%%%ds",4+nd);
3646 sprintf(xvg2format,"%s %s\n","%s",dgformat);
3647 sprintf(xvg3format,"%s %s %s\n","%s",dgformat,dgformat);
3652 if (opt2bSet("-o",NFILE,fnm))
3654 sprintf(buf,"%s (%s)","\\DeltaG","kT");
3655 fpb = xvgropen_type(opt2fn("-o",NFILE,fnm),"Free energy differences",
3656 "\\lambda",buf,exvggtXYDY,oenv);
3660 if (opt2bSet("-oi",NFILE,fnm))
3662 sprintf(buf,"%s (%s)","\\DeltaG","kT");
3663 fpi = xvgropen(opt2fn("-oi",NFILE,fnm),"Free energy integral",
3664 "\\lambda",buf,oenv);
3672 /* first calculate results */
3675 for(f=0; f<nresults; f++)
3677 /* Determine the free energy difference with a factor of 10
3678 * more accuracy than requested for printing.
3680 calc_bar(&(results[f]), 0.1*prec, nbmin, nbmax,
3683 if (results[f].dg_disc_err > prec/10.)
3685 if (results[f].dg_histrange_err > prec/10.)
3689 /* print results in kT */
3693 printf("\nTemperature: %g K\n", temp);
3695 printf("\nDetailed results in kT (see help for explanation):\n\n");
3696 printf("%6s ", " lam_A");
3697 printf("%6s ", " lam_B");
3698 printf(sktformat, "DG ");
3700 printf(skteformat, "+/- ");
3702 printf(skteformat, "disc ");
3704 printf(skteformat, "range ");
3705 printf(sktformat, "s_A ");
3707 printf(skteformat, "+/- " );
3708 printf(sktformat, "s_B ");
3710 printf(skteformat, "+/- " );
3711 printf(sktformat, "stdev ");
3713 printf(skteformat, "+/- ");
3715 for(f=0; f<nresults; f++)
3717 lambda_vec_print_short(results[f].a->native_lambda, buf);
3719 lambda_vec_print_short(results[f].b->native_lambda, buf);
3721 printf(ktformat, results[f].dg);
3725 printf(kteformat, results[f].dg_err);
3730 printf(kteformat, results[f].dg_disc_err);
3735 printf(kteformat, results[f].dg_histrange_err);
3738 printf(ktformat, results[f].sa);
3742 printf(kteformat, results[f].sa_err);
3745 printf(ktformat, results[f].sb);
3749 printf(kteformat, results[f].sb_err);
3752 printf(ktformat, results[f].dg_stddev);
3756 printf(kteformat, results[f].dg_stddev_err);
3760 /* Check for negative relative entropy with a 95% certainty. */
3761 if (results[f].sa < -2*results[f].sa_err ||
3762 results[f].sb < -2*results[f].sb_err)
3770 printf("\nWARNING: Some of these results violate the Second Law of "
3771 "Thermodynamics: \n"
3772 " This is can be the result of severe undersampling, or "
3774 " there is something wrong with the simulations.\n");
3778 /* final results in kJ/mol */
3779 printf("\n\nFinal results in kJ/mol:\n\n");
3781 for(f=0; f<nresults; f++)
3786 lambda_vec_print_short(results[f].a->native_lambda, buf);
3787 fprintf(fpi, xvg2format, buf, dg_tot);
3793 lambda_vec_print_intermediate(results[f].a->native_lambda,
3794 results[f].b->native_lambda,
3797 fprintf(fpb, xvg3format, buf, results[f].dg,results[f].dg_err);
3801 lambda_vec_print_short(results[f].a->native_lambda, buf);
3802 lambda_vec_print_short(results[f].b->native_lambda, buf2);
3803 printf("%s - %s", buf, buf2);
3806 printf(dgformat,results[f].dg*kT);
3810 printf(dgformat,results[f].dg_err*kT);
3814 printf(" (max. range err. = ");
3815 printf(dgformat,results[f].dg_histrange_err*kT);
3817 sum_histrange_err += results[f].dg_histrange_err*kT;
3821 dg_tot += results[f].dg;
3825 lambda_vec_print_short(results[0].a->native_lambda, buf);
3826 lambda_vec_print_short(results[nresults-1].b->native_lambda, buf2);
3827 printf("%s - %s", buf, buf2);
3830 printf(dgformat,dg_tot*kT);
3833 stat_err=bar_err(nbmin,nbmax,partsum)*kT;
3835 printf(dgformat, max(max(stat_err, sum_disc_err), sum_histrange_err));
3840 printf("\nmaximum discretization error = ");
3841 printf(dgformat,sum_disc_err);
3842 if (bEE && stat_err < sum_disc_err)
3844 printf("WARNING: discretization error (%g) is larger than statistical error.\n Decrease histogram spacing for more accurate results\n", stat_err);
3849 printf("\nmaximum histogram range error = ");
3850 printf(dgformat,sum_histrange_err);
3851 if (bEE && stat_err < sum_histrange_err)
3853 printf("WARNING: histogram range error (%g) is larger than statistical error.\n Increase histogram range for more accurate results\n", stat_err);
3862 lambda_vec_print_short(results[nresults-1].b->native_lambda, buf);
3863 fprintf(fpi, xvg2format, buf, dg_tot);
3867 do_view(oenv,opt2fn_null("-o",NFILE,fnm),"-xydy");
3868 do_view(oenv,opt2fn_null("-oi",NFILE,fnm),"-xydy");