1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
55 #include "gmx_fatal.h"
64 /* Structure for the names of lambda vector components */
65 typedef struct lambda_components_t
67 char **names; /* Array of strings with names for the lambda vector
69 int N; /* The number of components */
70 int Nalloc; /* The number of allocated components */
71 } lambda_components_t;
73 /* Structure for a lambda vector or a dhdl derivative direction */
74 typedef struct lambda_vec_t
76 double *val; /* The lambda vector component values. Only valid if
78 int dhdl; /* The coordinate index for the derivative described by this
80 const lambda_components_t *lc; /* the associated lambda_components
82 int index; /* The state number (init-lambda-state) of this lambda
83 vector, if known. If not, it is set to -1 */
86 /* the dhdl.xvg data from a simulation */
90 int ftp; /* file type */
91 int nset; /* number of lambdas, including dhdl */
92 int *np; /* number of data points (du or hists) per lambda */
93 int np_alloc; /* number of points (du or hists) allocated */
94 double temp; /* temperature */
95 lambda_vec_t *lambda; /* the lambdas (of first index for y). */
96 double *t; /* the times (of second index for y) */
97 double **y; /* the dU values. y[0] holds the derivative, while
98 further ones contain the energy differences between
99 the native lambda and the 'foreign' lambdas. */
100 lambda_vec_t native_lambda; /* the native lambda */
102 struct xvg_t *next, *prev; /*location in the global linked list of xvg_ts*/
106 typedef struct hist_t
108 unsigned int *bin[2]; /* the (forward + reverse) histogram values */
109 double dx[2]; /* the histogram spacing. The reverse
110 dx is the negative of the forward dx.*/
111 gmx_large_int_t x0[2]; /* the (forward + reverse) histogram start
114 int nbin[2]; /* the (forward+reverse) number of bins */
115 gmx_large_int_t sum; /* the total number of counts. Must be
116 the same for forward + reverse. */
117 int nhist; /* number of hist datas (forward or reverse) */
119 double start_time, delta_time; /* start time, end time of histogram */
123 /* an aggregate of samples for partial free energy calculation */
124 typedef struct samples_t
126 lambda_vec_t *native_lambda; /* pointer to native lambda vector */
127 lambda_vec_t *foreign_lambda; /* pointer to foreign lambda vector */
128 double temp; /* the temperature */
129 gmx_bool derivative; /* whether this sample is a derivative */
131 /* The samples come either as either delta U lists: */
132 int ndu; /* the number of delta U samples */
133 double *du; /* the delta u's */
134 double *t; /* the times associated with those samples, or: */
135 double start_time, delta_time; /*start time and delta time for linear time*/
137 /* or as histograms: */
138 hist_t *hist; /* a histogram */
140 /* allocation data: (not NULL for data 'owned' by this struct) */
141 double *du_alloc, *t_alloc; /* allocated delta u arrays */
142 size_t ndu_alloc, nt_alloc; /* pre-allocated sizes */
143 hist_t *hist_alloc; /* allocated hist */
145 gmx_large_int_t ntot; /* total number of samples */
146 const char *filename; /* the file name this sample comes from */
149 /* a sample range (start to end for du-style data, or boolean
150 for both du-style data and histograms */
151 typedef struct sample_range_t
153 int start, end; /* start and end index for du style data */
154 gmx_bool use; /* whether to use this sample */
156 samples_t *s; /* the samples this range belongs to */
160 /* a collection of samples for a partial free energy calculation
161 (i.e. the collection of samples from one native lambda to one
163 typedef struct sample_coll_t
165 lambda_vec_t *native_lambda; /* these should be the same for all samples
167 lambda_vec_t *foreign_lambda; /* collection */
168 double temp; /* the temperature */
170 int nsamples; /* the number of samples */
171 samples_t **s; /* the samples themselves */
172 sample_range_t *r; /* the sample ranges */
173 int nsamples_alloc; /* number of allocated samples */
175 gmx_large_int_t ntot; /* total number of samples in the ranges of
178 struct sample_coll_t *next, *prev; /* next and previous in the list */
181 /* all the samples associated with a lambda point */
182 typedef struct lambda_data_t
184 lambda_vec_t *lambda; /* the native lambda (at start time if dynamic) */
185 double temp; /* temperature */
187 sample_coll_t *sc; /* the samples */
189 sample_coll_t sc_head; /*the pre-allocated list head for the linked list.*/
191 struct lambda_data_t *next, *prev; /* the next and prev in the list */
194 /* Top-level data structure of simulation data */
195 typedef struct sim_data_t
197 lambda_data_t *lb; /* a lambda data linked list */
198 lambda_data_t lb_head; /* The head element of the linked list */
200 lambda_components_t lc; /* the allowed components of the lambda
204 /* Top-level data structure with calculated values. */
206 sample_coll_t *a, *b; /* the simulation data */
208 double dg; /* the free energy difference */
209 double dg_err; /* the free energy difference */
211 double dg_disc_err; /* discretization error */
212 double dg_histrange_err; /* histogram range error */
214 double sa; /* relative entropy of b in state a */
215 double sa_err; /* error in sa */
216 double sb; /* relative entropy of a in state b */
217 double sb_err; /* error in sb */
219 double dg_stddev; /* expected dg stddev per sample */
220 double dg_stddev_err; /* error in dg_stddev */
224 /* Initialize a lambda_components structure */
225 static void lambda_components_init(lambda_components_t *lc)
229 snew(lc->names, lc->Nalloc);
232 /* Add a component to a lambda_components structure */
233 static void lambda_components_add(lambda_components_t *lc,
234 const char *name, size_t name_length)
236 while (lc->N + 1 > lc->Nalloc)
238 lc->Nalloc = (lc->Nalloc == 0) ? 2 : 2*lc->Nalloc;
239 srealloc( lc->names, lc->Nalloc );
241 snew(lc->names[lc->N], name_length+1);
242 strncpy(lc->names[lc->N], name, name_length);
246 /* check whether a component with index 'index' matches the given name, or
247 is also NULL. Returns TRUE if this is the case.
248 the string name does not need to end */
249 static gmx_bool lambda_components_check(const lambda_components_t *lc,
259 if (name == NULL && lc->names[index] == NULL)
263 if ((name == NULL) != (lc->names[index] == NULL))
267 len = strlen(lc->names[index]);
268 if (len != name_length)
272 if (strncmp(lc->names[index], name, name_length) == 0)
279 /* Find the index of a given lambda component name, or -1 if not found */
280 static int lambda_components_find(const lambda_components_t *lc,
286 for (i = 0; i < lc->N; i++)
288 if (strncmp(lc->names[i], name, name_length) == 0)
298 /* initialize a lambda vector */
299 static void lambda_vec_init(lambda_vec_t *lv, const lambda_components_t *lc)
301 snew(lv->val, lc->N);
307 static void lambda_vec_destroy(lambda_vec_t *lv)
312 static void lambda_vec_copy(lambda_vec_t *lv, const lambda_vec_t *orig)
316 lambda_vec_init(lv, orig->lc);
317 lv->dhdl = orig->dhdl;
318 lv->index = orig->index;
319 for (i = 0; i < lv->lc->N; i++)
321 lv->val[i] = orig->val[i];
325 /* write a lambda vec to a preallocated string */
326 static void lambda_vec_print(const lambda_vec_t *lv, char *str, gmx_bool named)
331 str[0] = 0; /* reset the string */
336 str += sprintf(str, "delta H to ");
340 str += sprintf(str, "(");
342 for (i = 0; i < lv->lc->N; i++)
344 str += sprintf(str, "%g", lv->val[i]);
347 str += sprintf(str, ", ");
352 str += sprintf(str, ")");
357 /* this lambda vector describes a derivative */
358 str += sprintf(str, "dH/dl");
359 if (strlen(lv->lc->names[lv->dhdl]) > 0)
361 str += sprintf(str, " (%s)", lv->lc->names[lv->dhdl]);
366 /* write a shortened version of the lambda vec to a preallocated string */
367 static void lambda_vec_print_short(const lambda_vec_t *lv, char *str)
374 sprintf(str, "%6d", lv->index);
380 sprintf(str, "%6.3f", lv->val[0]);
384 sprintf(str, "dH/dl[%d]", lv->dhdl);
389 /* write an intermediate version of two lambda vecs to a preallocated string */
390 static void lambda_vec_print_intermediate(const lambda_vec_t *a,
391 const lambda_vec_t *b, char *str)
397 if ( (a->index >= 0) && (b->index >= 0) )
399 sprintf(str, "%6.3f", ((double)a->index+(double)b->index)/2.);
403 if ( (a->dhdl < 0) && (b->dhdl < 0) )
405 sprintf(str, "%6.3f", (a->val[0]+b->val[0])/2.);
412 /* calculate the difference in lambda vectors: c = a-b.
413 c must be initialized already, and a and b must describe non-derivative
415 static void lambda_vec_diff(const lambda_vec_t *a, const lambda_vec_t *b,
420 if ( (a->dhdl > 0) || (b->dhdl > 0) )
423 "Trying to calculate the difference between derivatives instead of lambda points");
425 if ((a->lc != b->lc) || (a->lc != c->lc) )
428 "Trying to calculate the difference lambdas with differing basis set");
430 for (i = 0; i < a->lc->N; i++)
432 c->val[i] = a->val[i] - b->val[i];
436 /* calculate and return the absolute difference in lambda vectors: c = |a-b|.
437 a and b must describe non-derivative lambda points */
438 static double lambda_vec_abs_diff(const lambda_vec_t *a, const lambda_vec_t *b)
443 if ( (a->dhdl > 0) || (b->dhdl > 0) )
446 "Trying to calculate the difference between derivatives instead of lambda points");
451 "Trying to calculate the difference lambdas with differing basis set");
453 for (i = 0; i < a->lc->N; i++)
455 double df = a->val[i] - b->val[i];
462 /* check whether two lambda vectors are the same */
463 static gmx_bool lambda_vec_same(const lambda_vec_t *a, const lambda_vec_t *b)
473 for (i = 0; i < a->lc->N; i++)
475 if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
484 /* they're derivatives, so we check whether the indices match */
485 return (a->dhdl == b->dhdl);
489 /* Compare the sort order of two foreign lambda vectors
491 returns 1 if a is 'bigger' than b,
492 returns 0 if they're the same,
493 returns -1 if a is 'smaller' than b.*/
494 static gmx_bool lambda_vec_cmp_foreign(const lambda_vec_t *a,
495 const lambda_vec_t *b)
498 double norm_a = 0, norm_b = 0;
499 gmx_bool different = FALSE;
503 gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
505 /* if either one has an index we sort based on that */
506 if ((a->index >= 0) || (b->index >= 0))
508 if (a->index == b->index)
512 return (a->index > b->index) ? 1 : -1;
514 if (a->dhdl >= 0 || b->dhdl >= 0)
516 /* lambda vectors that are derivatives always sort higher than those
517 without derivatives */
518 if ((a->dhdl >= 0) != (b->dhdl >= 0) )
520 return (a->dhdl >= 0) ? 1 : -1;
522 return a->dhdl > b->dhdl;
525 /* neither has an index, so we can only sort on the lambda components,
526 which is only valid if there is one component */
527 for (i = 0; i < a->lc->N; i++)
529 if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
533 norm_a += a->val[i]*a->val[i];
534 norm_b += b->val[i]*b->val[i];
540 return norm_a > norm_b;
543 /* Compare the sort order of two native lambda vectors
545 returns 1 if a is 'bigger' than b,
546 returns 0 if they're the same,
547 returns -1 if a is 'smaller' than b.*/
548 static gmx_bool lambda_vec_cmp_native(const lambda_vec_t *a,
549 const lambda_vec_t *b)
555 gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
557 /* if either one has an index we sort based on that */
558 if ((a->index >= 0) || (b->index >= 0))
560 if (a->index == b->index)
564 return (a->index > b->index) ? 1 : -1;
566 /* neither has an index, so we can only sort on the lambda components,
567 which is only valid if there is one component */
571 "Can't compare lambdas with no index and > 1 component");
573 if (a->dhdl >= 0 || b->dhdl >= 0)
576 "Can't compare native lambdas that are derivatives");
578 if (gmx_within_tol(a->val[0], b->val[0], 10*GMX_REAL_EPS))
582 return a->val[0] > b->val[0] ? 1 : -1;
588 static void hist_init(hist_t *h, int nhist, int *nbin)
593 gmx_fatal(FARGS, "histogram with more than two sets of data!");
595 for (i = 0; i < nhist; i++)
597 snew(h->bin[i], nbin[i]);
599 h->nbin[i] = nbin[i];
600 h->start_time = h->delta_time = 0;
607 static void hist_destroy(hist_t *h)
613 static void xvg_init(xvg_t *ba)
622 static void samples_init(samples_t *s, lambda_vec_t *native_lambda,
623 lambda_vec_t *foreign_lambda, double temp,
624 gmx_bool derivative, const char *filename)
626 s->native_lambda = native_lambda;
627 s->foreign_lambda = foreign_lambda;
629 s->derivative = derivative;
634 s->start_time = s->delta_time = 0;
638 s->hist_alloc = NULL;
643 s->filename = filename;
646 static void sample_range_init(sample_range_t *r, samples_t *s)
654 static void sample_coll_init(sample_coll_t *sc, lambda_vec_t *native_lambda,
655 lambda_vec_t *foreign_lambda, double temp)
657 sc->native_lambda = native_lambda;
658 sc->foreign_lambda = foreign_lambda;
664 sc->nsamples_alloc = 0;
667 sc->next = sc->prev = NULL;
670 static void sample_coll_destroy(sample_coll_t *sc)
672 /* don't free the samples themselves */
678 static void lambda_data_init(lambda_data_t *l, lambda_vec_t *native_lambda,
681 l->lambda = native_lambda;
687 l->sc = &(l->sc_head);
689 sample_coll_init(l->sc, native_lambda, NULL, 0.);
694 static void barres_init(barres_t *br)
703 br->dg_stddev_err = 0;
710 /* calculate the total number of samples in a sample collection */
711 static void sample_coll_calc_ntot(sample_coll_t *sc)
716 for (i = 0; i < sc->nsamples; i++)
722 sc->ntot += sc->s[i]->ntot;
726 sc->ntot += sc->r[i].end - sc->r[i].start;
733 /* find the barsamples_t associated with a lambda that corresponds to
734 a specific foreign lambda */
735 static sample_coll_t *lambda_data_find_sample_coll(lambda_data_t *l,
736 lambda_vec_t *foreign_lambda)
738 sample_coll_t *sc = l->sc->next;
742 if (lambda_vec_same(sc->foreign_lambda, foreign_lambda))
752 /* insert li into an ordered list of lambda_colls */
753 static void lambda_data_insert_sample_coll(lambda_data_t *l, sample_coll_t *sc)
755 sample_coll_t *scn = l->sc->next;
756 while ( (scn != l->sc) )
758 if (lambda_vec_cmp_foreign(scn->foreign_lambda, sc->foreign_lambda) > 0)
764 /* now insert it before the found scn */
766 sc->prev = scn->prev;
767 scn->prev->next = sc;
771 /* insert li into an ordered list of lambdas */
772 static void lambda_data_insert_lambda(lambda_data_t *head, lambda_data_t *li)
774 lambda_data_t *lc = head->next;
777 if (lambda_vec_cmp_native(lc->lambda, li->lambda) > 0)
783 /* now insert ourselves before the found lc */
790 /* insert a sample and a sample_range into a sample_coll. The
791 samples are stored as a pointer, the range is copied. */
792 static void sample_coll_insert_sample(sample_coll_t *sc, samples_t *s,
795 /* first check if it belongs here */
796 if (sc->temp != s->temp)
798 gmx_fatal(FARGS, "Temperatures in files %s and %s are not the same!",
799 s->filename, sc->next->s[0]->filename);
801 if (!lambda_vec_same(sc->native_lambda, s->native_lambda))
803 gmx_fatal(FARGS, "Native lambda in files %s and %s are not the same (and they should be)!",
804 s->filename, sc->next->s[0]->filename);
806 if (!lambda_vec_same(sc->foreign_lambda, s->foreign_lambda))
808 gmx_fatal(FARGS, "Foreign lambda in files %s and %s are not the same (and they should be)!",
809 s->filename, sc->next->s[0]->filename);
812 /* check if there's room */
813 if ( (sc->nsamples + 1) > sc->nsamples_alloc)
815 sc->nsamples_alloc = max(2*sc->nsamples_alloc, 2);
816 srenew(sc->s, sc->nsamples_alloc);
817 srenew(sc->r, sc->nsamples_alloc);
819 sc->s[sc->nsamples] = s;
820 sc->r[sc->nsamples] = *r;
823 sample_coll_calc_ntot(sc);
826 /* insert a sample into a lambda_list, creating the right sample_coll if
828 static void lambda_data_list_insert_sample(lambda_data_t *head, samples_t *s)
830 gmx_bool found = FALSE;
834 lambda_data_t *l = head->next;
836 /* first search for the right lambda_data_t */
839 if (lambda_vec_same(l->lambda, s->native_lambda) )
849 snew(l, 1); /* allocate a new one */
850 lambda_data_init(l, s->native_lambda, s->temp); /* initialize it */
851 lambda_data_insert_lambda(head, l); /* add it to the list */
854 /* now look for a sample collection */
855 sc = lambda_data_find_sample_coll(l, s->foreign_lambda);
858 snew(sc, 1); /* allocate a new one */
859 sample_coll_init(sc, s->native_lambda, s->foreign_lambda, s->temp);
860 lambda_data_insert_sample_coll(l, sc);
863 /* now insert the samples into the sample coll */
864 sample_range_init(&r, s);
865 sample_coll_insert_sample(sc, s, &r);
869 /* make a histogram out of a sample collection */
870 static void sample_coll_make_hist(sample_coll_t *sc, int **bin,
871 int *nbin_alloc, int *nbin,
872 double *dx, double *xmin, int nbin_default)
875 gmx_bool dx_set = FALSE;
876 gmx_bool xmin_set = FALSE;
878 gmx_bool xmax_set = FALSE;
879 gmx_bool xmax_set_hard = FALSE; /* whether the xmax is bounded by the
880 limits of a histogram */
883 /* first determine dx and xmin; try the histograms */
884 for (i = 0; i < sc->nsamples; i++)
888 hist_t *hist = sc->s[i]->hist;
889 for (k = 0; k < hist->nhist; k++)
891 double hdx = hist->dx[k];
892 double xmax_now = (hist->x0[k]+hist->nbin[k])*hdx;
894 /* we use the biggest dx*/
895 if ( (!dx_set) || hist->dx[0] > *dx)
900 if ( (!xmin_set) || (hist->x0[k]*hdx) < *xmin)
903 *xmin = (hist->x0[k]*hdx);
906 if ( (!xmax_set) || (xmax_now > xmax && !xmax_set_hard) )
910 if (hist->bin[k][hist->nbin[k]-1] != 0)
912 xmax_set_hard = TRUE;
915 if (hist->bin[k][hist->nbin[k]-1] != 0 && (xmax_now < xmax) )
917 xmax_set_hard = TRUE;
923 /* and the delta us */
924 for (i = 0; i < sc->nsamples; i++)
926 if (sc->s[i]->ndu > 0)
928 /* determine min and max */
929 int starti = sc->r[i].start;
930 int endi = sc->r[i].end;
931 double du_xmin = sc->s[i]->du[starti];
932 double du_xmax = sc->s[i]->du[starti];
933 for (j = starti+1; j < endi; j++)
935 if (sc->s[i]->du[j] < du_xmin)
937 du_xmin = sc->s[i]->du[j];
939 if (sc->s[i]->du[j] > du_xmax)
941 du_xmax = sc->s[i]->du[j];
945 /* and now change the limits */
946 if ( (!xmin_set) || (du_xmin < *xmin) )
951 if ( (!xmax_set) || ((du_xmax > xmax) && !xmax_set_hard) )
959 if (!xmax_set || !xmin_set)
968 *nbin = nbin_default;
969 *dx = (xmax-(*xmin))/((*nbin)-2); /* -2 because we want the last bin to
970 be 0, and we count from 0 */
974 *nbin = (xmax-(*xmin))/(*dx);
977 if (*nbin > *nbin_alloc)
980 srenew(*bin, *nbin_alloc);
983 /* reset the histogram */
984 for (i = 0; i < (*nbin); i++)
989 /* now add the actual data */
990 for (i = 0; i < sc->nsamples; i++)
994 hist_t *hist = sc->s[i]->hist;
995 for (k = 0; k < hist->nhist; k++)
997 double hdx = hist->dx[k];
998 double xmin_hist = hist->x0[k]*hdx;
999 for (j = 0; j < hist->nbin[k]; j++)
1001 /* calculate the bin corresponding to the middle of the
1003 double x = hdx*(j+0.5) + xmin_hist;
1004 int binnr = (int)((x-(*xmin))/(*dx));
1006 if (binnr >= *nbin || binnr < 0)
1011 (*bin)[binnr] += hist->bin[k][j];
1017 int starti = sc->r[i].start;
1018 int endi = sc->r[i].end;
1019 for (j = starti; j < endi; j++)
1021 int binnr = (int)((sc->s[i]->du[j]-(*xmin))/(*dx));
1022 if (binnr >= *nbin || binnr < 0)
1033 /* write a collection of histograms to a file */
1034 void sim_data_histogram(sim_data_t *sd, const char *filename,
1035 int nbin_default, const output_env_t oenv)
1037 char label_x[STRLEN];
1038 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1039 const char *title = "N(\\DeltaH)";
1040 const char *label_y = "Samples";
1044 char **setnames = NULL;
1045 gmx_bool first_set = FALSE;
1046 /* histogram data: */
1053 lambda_data_t *bl_head = sd->lb;
1055 printf("\nWriting histogram to %s\n", filename);
1056 sprintf(label_x, "\\DeltaH (%s)", unit_energy);
1058 fp = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1060 /* first get all the set names */
1062 /* iterate over all lambdas */
1063 while (bl != bl_head)
1065 sample_coll_t *sc = bl->sc->next;
1067 /* iterate over all samples */
1068 while (sc != bl->sc)
1070 char buf[STRLEN], buf2[STRLEN];
1073 srenew(setnames, nsets);
1074 snew(setnames[nsets-1], STRLEN);
1075 if (sc->foreign_lambda->dhdl < 0)
1077 lambda_vec_print(sc->native_lambda, buf, FALSE);
1078 lambda_vec_print(sc->foreign_lambda, buf2, FALSE);
1079 sprintf(setnames[nsets-1], "N(%s(%s=%s) | %s=%s)",
1080 deltag, lambda, buf2, lambda, buf);
1084 lambda_vec_print(sc->native_lambda, buf, FALSE);
1085 sprintf(setnames[nsets-1], "N(%s | %s=%s)",
1093 xvgr_legend(fp, nsets, (const char**)setnames, oenv);
1096 /* now make the histograms */
1098 /* iterate over all lambdas */
1099 while (bl != bl_head)
1101 sample_coll_t *sc = bl->sc->next;
1103 /* iterate over all samples */
1104 while (sc != bl->sc)
1108 xvgr_new_dataset(fp, 0, 0, NULL, oenv);
1111 sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &min,
1114 for (i = 0; i < nbin; i++)
1116 double xmin = i*dx + min;
1117 double xmax = (i+1)*dx + min;
1119 fprintf(fp, "%g %d\n%g %d\n", xmin, hist[i], xmax, hist[i]);
1137 /* create a collection (array) of barres_t object given a ordered linked list
1138 of barlamda_t sample collections */
1139 static barres_t *barres_list_create(sim_data_t *sd, int *nres,
1146 gmx_bool dhdl = FALSE;
1147 gmx_bool first = TRUE;
1148 lambda_data_t *bl_head = sd->lb;
1150 /* first count the lambdas */
1152 while (bl != bl_head)
1157 snew(res, nlambda-1);
1159 /* next put the right samples in the res */
1161 bl = bl_head->next->next; /* we start with the second one. */
1162 while (bl != bl_head)
1164 sample_coll_t *sc, *scprev;
1165 barres_t *br = &(res[*nres]);
1166 /* there is always a previous one. we search for that as a foreign
1168 scprev = lambda_data_find_sample_coll(bl->prev, bl->lambda);
1169 sc = lambda_data_find_sample_coll(bl, bl->prev->lambda);
1177 scprev = lambda_data_find_sample_coll(bl->prev, bl->prev->lambda);
1178 sc = lambda_data_find_sample_coll(bl, bl->lambda);
1182 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");
1187 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");
1190 else if (!scprev && !sc)
1192 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);
1195 /* normal delta H */
1198 gmx_fatal(FARGS, "Could not find a set for foreign lambda = %g\nin the files for lambda = %g", bl->lambda, bl->prev->lambda);
1202 gmx_fatal(FARGS, "Could not find a set for foreign lambda = %g\nin the files for lambda = %g", bl->prev->lambda, bl->lambda);
1214 /* estimate the maximum discretization error */
1215 static double barres_list_max_disc_err(barres_t *res, int nres)
1218 double disc_err = 0.;
1219 double delta_lambda;
1221 for (i = 0; i < nres; i++)
1223 barres_t *br = &(res[i]);
1225 delta_lambda = lambda_vec_abs_diff(br->b->native_lambda,
1226 br->a->native_lambda);
1228 for (j = 0; j < br->a->nsamples; j++)
1230 if (br->a->s[j]->hist)
1233 if (br->a->s[j]->derivative)
1235 Wfac = delta_lambda;
1238 disc_err = max(disc_err, Wfac*br->a->s[j]->hist->dx[0]);
1241 for (j = 0; j < br->b->nsamples; j++)
1243 if (br->b->s[j]->hist)
1246 if (br->b->s[j]->derivative)
1248 Wfac = delta_lambda;
1250 disc_err = max(disc_err, Wfac*br->b->s[j]->hist->dx[0]);
1258 /* impose start and end times on a sample collection, updating sample_ranges */
1259 static void sample_coll_impose_times(sample_coll_t *sc, double begin_t,
1263 for (i = 0; i < sc->nsamples; i++)
1265 samples_t *s = sc->s[i];
1266 sample_range_t *r = &(sc->r[i]);
1269 double end_time = s->hist->delta_time*s->hist->sum +
1270 s->hist->start_time;
1271 if (s->hist->start_time < begin_t || end_time > end_t)
1281 if (s->start_time < begin_t)
1283 r->start = (int)((begin_t - s->start_time)/s->delta_time);
1285 end_time = s->delta_time*s->ndu + s->start_time;
1286 if (end_time > end_t)
1288 r->end = (int)((end_t - s->start_time)/s->delta_time);
1294 for (j = 0; j < s->ndu; j++)
1296 if (s->t[j] < begin_t)
1301 if (s->t[j] >= end_t)
1308 if (r->start > r->end)
1314 sample_coll_calc_ntot(sc);
1317 static void sim_data_impose_times(sim_data_t *sd, double begin, double end)
1319 double first_t, last_t;
1320 double begin_t, end_t;
1322 lambda_data_t *head = sd->lb;
1325 if (begin <= 0 && end < 0)
1330 /* first determine the global start and end times */
1336 sample_coll_t *sc = lc->sc->next;
1337 while (sc != lc->sc)
1339 for (j = 0; j < sc->nsamples; j++)
1341 double start_t, end_t;
1343 start_t = sc->s[j]->start_time;
1344 end_t = sc->s[j]->start_time;
1347 end_t += sc->s[j]->delta_time*sc->s[j]->hist->sum;
1353 end_t = sc->s[j]->t[sc->s[j]->ndu-1];
1357 end_t += sc->s[j]->delta_time*sc->s[j]->ndu;
1361 if (start_t < first_t || first_t < 0)
1375 /* calculate the actual times */
1393 printf("\n Samples in time interval: %.3f - %.3f\n", first_t, last_t);
1395 if (begin_t > end_t)
1399 printf("Removing samples outside of: %.3f - %.3f\n", begin_t, end_t);
1401 /* then impose them */
1405 sample_coll_t *sc = lc->sc->next;
1406 while (sc != lc->sc)
1408 sample_coll_impose_times(sc, begin_t, end_t);
1416 /* create subsample i out of ni from an existing sample_coll */
1417 static gmx_bool sample_coll_create_subsample(sample_coll_t *sc,
1418 sample_coll_t *sc_orig,
1422 int hist_start, hist_end;
1424 gmx_large_int_t ntot_start;
1425 gmx_large_int_t ntot_end;
1426 gmx_large_int_t ntot_so_far;
1428 *sc = *sc_orig; /* just copy all fields */
1430 /* allocate proprietary memory */
1431 snew(sc->s, sc_orig->nsamples);
1432 snew(sc->r, sc_orig->nsamples);
1434 /* copy the samples */
1435 for (j = 0; j < sc_orig->nsamples; j++)
1437 sc->s[j] = sc_orig->s[j];
1438 sc->r[j] = sc_orig->r[j]; /* copy the ranges too */
1441 /* now fix start and end fields */
1442 /* the casts avoid possible overflows */
1443 ntot_start = (gmx_large_int_t)(sc_orig->ntot*(double)i/(double)ni);
1444 ntot_end = (gmx_large_int_t)(sc_orig->ntot*(double)(i+1)/(double)ni);
1446 for (j = 0; j < sc->nsamples; j++)
1448 gmx_large_int_t ntot_add;
1449 gmx_large_int_t new_start, new_end;
1455 ntot_add = sc->s[j]->hist->sum;
1459 ntot_add = sc->r[j].end - sc->r[j].start;
1467 if (!sc->s[j]->hist)
1469 if (ntot_so_far < ntot_start)
1471 /* adjust starting point */
1472 new_start = sc->r[j].start + (ntot_start - ntot_so_far);
1476 new_start = sc->r[j].start;
1478 /* adjust end point */
1479 new_end = sc->r[j].start + (ntot_end - ntot_so_far);
1480 if (new_end > sc->r[j].end)
1482 new_end = sc->r[j].end;
1485 /* check if we're in range at all */
1486 if ( (new_end < new_start) || (new_start > sc->r[j].end) )
1491 /* and write the new range */
1492 sc->r[j].start = (int)new_start;
1493 sc->r[j].end = (int)new_end;
1500 double ntot_start_norm, ntot_end_norm;
1501 /* calculate the amount of overlap of the
1502 desired range (ntot_start -- ntot_end) onto
1503 the histogram range (ntot_so_far -- ntot_so_far+ntot_add)*/
1505 /* first calculate normalized bounds
1506 (where 0 is the start of the hist range, and 1 the end) */
1507 ntot_start_norm = (ntot_start-ntot_so_far)/(double)ntot_add;
1508 ntot_end_norm = (ntot_end-ntot_so_far)/(double)ntot_add;
1510 /* now fix the boundaries */
1511 ntot_start_norm = min(1, max(0., ntot_start_norm));
1512 ntot_end_norm = max(0, min(1., ntot_end_norm));
1514 /* and calculate the overlap */
1515 overlap = ntot_end_norm - ntot_start_norm;
1517 if (overlap > 0.95) /* we allow for 5% slack */
1519 sc->r[j].use = TRUE;
1521 else if (overlap < 0.05)
1523 sc->r[j].use = FALSE;
1531 ntot_so_far += ntot_add;
1533 sample_coll_calc_ntot(sc);
1538 /* calculate minimum and maximum work values in sample collection */
1539 static void sample_coll_min_max(sample_coll_t *sc, double Wfac,
1540 double *Wmin, double *Wmax)
1547 for (i = 0; i < sc->nsamples; i++)
1549 samples_t *s = sc->s[i];
1550 sample_range_t *r = &(sc->r[i]);
1555 for (j = r->start; j < r->end; j++)
1557 *Wmin = min(*Wmin, s->du[j]*Wfac);
1558 *Wmax = max(*Wmax, s->du[j]*Wfac);
1563 int hd = 0; /* determine the histogram direction: */
1565 if ( (s->hist->nhist > 1) && (Wfac < 0) )
1569 dx = s->hist->dx[hd];
1571 for (j = s->hist->nbin[hd]-1; j >= 0; j--)
1573 *Wmin = min(*Wmin, Wfac*(s->hist->x0[hd])*dx);
1574 *Wmax = max(*Wmax, Wfac*(s->hist->x0[hd])*dx);
1575 /* look for the highest value bin with values */
1576 if (s->hist->bin[hd][j] > 0)
1578 *Wmin = min(*Wmin, Wfac*(j+s->hist->x0[hd]+1)*dx);
1579 *Wmax = max(*Wmax, Wfac*(j+s->hist->x0[hd]+1)*dx);
1588 /* Initialize a sim_data structure */
1589 static void sim_data_init(sim_data_t *sd)
1591 /* make linked list */
1592 sd->lb = &(sd->lb_head);
1593 sd->lb->next = sd->lb;
1594 sd->lb->prev = sd->lb;
1596 lambda_components_init(&(sd->lc));
1600 static double calc_bar_sum(int n, const double *W, double Wfac, double sbMmDG)
1607 for (i = 0; i < n; i++)
1609 sum += 1./(1. + exp(Wfac*W[i] + sbMmDG));
1615 /* calculate the BAR average given a histogram
1617 if type== 0, calculate the best estimate for the average,
1618 if type==-1, calculate the minimum possible value given the histogram
1619 if type== 1, calculate the maximum possible value given the histogram */
1620 static double calc_bar_sum_hist(const hist_t *hist, double Wfac, double sbMmDG,
1626 /* normalization factor multiplied with bin width and
1627 number of samples (we normalize through M): */
1629 int hd = 0; /* determine the histogram direction: */
1632 if ( (hist->nhist > 1) && (Wfac < 0) )
1637 max = hist->nbin[hd]-1;
1640 max = hist->nbin[hd]; /* we also add whatever was out of range */
1643 for (i = 0; i < max; i++)
1645 double x = Wfac*((i+hist->x0[hd])+0.5)*dx; /* bin middle */
1646 double pxdx = hist->bin[0][i]*normdx; /* p(x)dx */
1648 sum += pxdx/(1. + exp(x + sbMmDG));
1654 static double calc_bar_lowlevel(sample_coll_t *ca, sample_coll_t *cb,
1655 double temp, double tol, int type)
1660 double Wfac1, Wfac2, Wmin, Wmax;
1661 double DG0, DG1, DG2, dDG1;
1663 double n1, n2; /* numbers of samples as doubles */
1668 /* count the numbers of samples */
1674 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1675 if (ca->foreign_lambda->dhdl < 0)
1677 /* this is the case when the delta U were calculated directly
1678 (i.e. we're not scaling dhdl) */
1684 /* we're using dhdl, so delta_lambda needs to be a
1685 multiplication factor. */
1686 /*double delta_lambda=cb->native_lambda-ca->native_lambda;*/
1687 double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1689 if (cb->native_lambda->lc->N > 1)
1692 "Can't (yet) do multi-component dhdl interpolation");
1695 Wfac1 = beta*delta_lambda;
1696 Wfac2 = -beta*delta_lambda;
1701 /* We print the output both in kT and kJ/mol.
1702 * Here we determine DG in kT, so when beta < 1
1703 * the precision has to be increased.
1708 /* Calculate minimum and maximum work to give an initial estimate of
1709 * delta G as their average.
1712 double Wmin1, Wmin2, Wmax1, Wmax2;
1713 sample_coll_min_max(ca, Wfac1, &Wmin1, &Wmax1);
1714 sample_coll_min_max(cb, Wfac2, &Wmin2, &Wmax2);
1716 Wmin = min(Wmin1, Wmin2);
1717 Wmax = max(Wmax1, Wmax2);
1725 fprintf(debug, "DG %9.5f %9.5f\n", DG0, DG2);
1727 /* We approximate by bisection: given our initial estimates
1728 we keep checking whether the halfway point is greater or
1729 smaller than what we get out of the BAR averages.
1731 For the comparison we can use twice the tolerance. */
1732 while (DG2 - DG0 > 2*tol)
1734 DG1 = 0.5*(DG0 + DG2);
1736 /*printf("Wfac1=%g, Wfac2=%g, beta=%g, DG1=%g\n",Wfac1,Wfac2,beta,
1739 /* calculate the BAR averages */
1742 for (i = 0; i < ca->nsamples; i++)
1744 samples_t *s = ca->s[i];
1745 sample_range_t *r = &(ca->r[i]);
1750 dDG1 += calc_bar_sum_hist(s->hist, Wfac1, (M-DG1), type);
1754 dDG1 += calc_bar_sum(r->end - r->start, s->du + r->start,
1759 for (i = 0; i < cb->nsamples; i++)
1761 samples_t *s = cb->s[i];
1762 sample_range_t *r = &(cb->r[i]);
1767 dDG1 -= calc_bar_sum_hist(s->hist, Wfac2, -(M-DG1), type);
1771 dDG1 -= calc_bar_sum(r->end - r->start, s->du + r->start,
1787 fprintf(debug, "DG %9.5f %9.5f\n", DG0, DG2);
1791 return 0.5*(DG0 + DG2);
1794 static void calc_rel_entropy(sample_coll_t *ca, sample_coll_t *cb,
1795 double temp, double dg, double *sa, double *sb)
1801 double Wfac1, Wfac2;
1807 /* count the numbers of samples */
1811 /* to ensure the work values are the same as during the delta_G */
1812 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1813 if (ca->foreign_lambda->dhdl < 0)
1815 /* this is the case when the delta U were calculated directly
1816 (i.e. we're not scaling dhdl) */
1822 /* we're using dhdl, so delta_lambda needs to be a
1823 multiplication factor. */
1824 double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1826 Wfac1 = beta*delta_lambda;
1827 Wfac2 = -beta*delta_lambda;
1830 /* first calculate the average work in both directions */
1831 for (i = 0; i < ca->nsamples; i++)
1833 samples_t *s = ca->s[i];
1834 sample_range_t *r = &(ca->r[i]);
1839 for (j = r->start; j < r->end; j++)
1841 W_ab += Wfac1*s->du[j];
1846 /* normalization factor multiplied with bin width and
1847 number of samples (we normalize through M): */
1850 int hd = 0; /* histogram direction */
1851 if ( (s->hist->nhist > 1) && (Wfac1 < 0) )
1855 dx = s->hist->dx[hd];
1857 for (j = 0; j < s->hist->nbin[0]; j++)
1859 double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1860 double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
1868 for (i = 0; i < cb->nsamples; i++)
1870 samples_t *s = cb->s[i];
1871 sample_range_t *r = &(cb->r[i]);
1876 for (j = r->start; j < r->end; j++)
1878 W_ba += Wfac1*s->du[j];
1883 /* normalization factor multiplied with bin width and
1884 number of samples (we normalize through M): */
1887 int hd = 0; /* histogram direction */
1888 if ( (s->hist->nhist > 1) && (Wfac2 < 0) )
1892 dx = s->hist->dx[hd];
1894 for (j = 0; j < s->hist->nbin[0]; j++)
1896 double x = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1897 double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
1905 /* then calculate the relative entropies */
1910 static void calc_dg_stddev(sample_coll_t *ca, sample_coll_t *cb,
1911 double temp, double dg, double *stddev)
1915 double sigmafact = 0.;
1917 double Wfac1, Wfac2;
1923 /* count the numbers of samples */
1927 /* to ensure the work values are the same as during the delta_G */
1928 /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1929 if (ca->foreign_lambda->dhdl < 0)
1931 /* this is the case when the delta U were calculated directly
1932 (i.e. we're not scaling dhdl) */
1938 /* we're using dhdl, so delta_lambda needs to be a
1939 multiplication factor. */
1940 double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1942 Wfac1 = beta*delta_lambda;
1943 Wfac2 = -beta*delta_lambda;
1949 /* calculate average in both directions */
1950 for (i = 0; i < ca->nsamples; i++)
1952 samples_t *s = ca->s[i];
1953 sample_range_t *r = &(ca->r[i]);
1958 for (j = r->start; j < r->end; j++)
1960 sigmafact += 1./(2. + 2.*cosh((M + Wfac1*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) && (Wfac1 < 0) )
1974 dx = s->hist->dx[hd];
1976 for (j = 0; j < s->hist->nbin[0]; j++)
1978 double x = Wfac1*((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)));
1986 for (i = 0; i < cb->nsamples; i++)
1988 samples_t *s = cb->s[i];
1989 sample_range_t *r = &(cb->r[i]);
1994 for (j = r->start; j < r->end; j++)
1996 sigmafact += 1./(2. + 2.*cosh((M - Wfac2*s->du[j] - dg)));
2001 /* normalization factor multiplied with bin width and
2002 number of samples (we normalize through M): */
2005 int hd = 0; /* histogram direction */
2006 if ( (s->hist->nhist > 1) && (Wfac2 < 0) )
2010 dx = s->hist->dx[hd];
2012 for (j = 0; j < s->hist->nbin[0]; j++)
2014 double x = Wfac2*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
2015 double pxdx = s->hist->bin[0][j]*normdx; /* p(x)dx */
2017 sigmafact += pxdx/(2. + 2.*cosh((M - x - dg)));
2023 sigmafact /= (n1 + n2);
2027 Shirts, Bair, Hooker & Pande, Phys. Rev. Lett 91, 140601 (2003): */
2028 *stddev = sqrt(((1./sigmafact) - ( (n1+n2)/n1 + (n1+n2)/n2 )));
2033 static void calc_bar(barres_t *br, double tol,
2034 int npee_min, int npee_max, gmx_bool *bEE,
2038 double dg_sig2, sa_sig2, sb_sig2, stddev_sig2; /* intermediate variance values
2039 for calculated quantities */
2040 int nsample1, nsample2;
2041 double temp = br->a->temp;
2043 double dg_min, dg_max;
2044 gmx_bool have_hist = FALSE;
2046 br->dg = calc_bar_lowlevel(br->a, br->b, temp, tol, 0);
2048 br->dg_disc_err = 0.;
2049 br->dg_histrange_err = 0.;
2051 /* check if there are histograms */
2052 for (i = 0; i < br->a->nsamples; i++)
2054 if (br->a->r[i].use && br->a->s[i]->hist)
2062 for (i = 0; i < br->b->nsamples; i++)
2064 if (br->b->r[i].use && br->b->s[i]->hist)
2072 /* calculate histogram-specific errors */
2075 dg_min = calc_bar_lowlevel(br->a, br->b, temp, tol, -1);
2076 dg_max = calc_bar_lowlevel(br->a, br->b, temp, tol, 1);
2078 if (fabs(dg_max - dg_min) > GMX_REAL_EPS*10)
2080 /* the histogram range error is the biggest of the differences
2081 between the best estimate and the extremes */
2082 br->dg_histrange_err = fabs(dg_max - dg_min);
2084 br->dg_disc_err = 0.;
2085 for (i = 0; i < br->a->nsamples; i++)
2087 if (br->a->s[i]->hist)
2089 br->dg_disc_err = max(br->dg_disc_err, br->a->s[i]->hist->dx[0]);
2092 for (i = 0; i < br->b->nsamples; i++)
2094 if (br->b->s[i]->hist)
2096 br->dg_disc_err = max(br->dg_disc_err, br->b->s[i]->hist->dx[0]);
2100 calc_rel_entropy(br->a, br->b, temp, br->dg, &(br->sa), &(br->sb));
2102 calc_dg_stddev(br->a, br->b, temp, br->dg, &(br->dg_stddev) );
2111 sample_coll_t ca, cb;
2113 /* initialize the samples */
2114 sample_coll_init(&ca, br->a->native_lambda, br->a->foreign_lambda,
2116 sample_coll_init(&cb, br->b->native_lambda, br->b->foreign_lambda,
2119 for (npee = npee_min; npee <= npee_max; npee++)
2128 double dstddev2 = 0;
2131 for (p = 0; p < npee; p++)
2138 cac = sample_coll_create_subsample(&ca, br->a, p, npee);
2139 cbc = sample_coll_create_subsample(&cb, br->b, p, npee);
2143 printf("WARNING: histogram number incompatible with block number for averaging: can't do error estimate\n");
2147 sample_coll_destroy(&ca);
2151 sample_coll_destroy(&cb);
2156 dgp = calc_bar_lowlevel(&ca, &cb, temp, tol, 0);
2160 partsum[npee*(npee_max+1)+p] += dgp;
2162 calc_rel_entropy(&ca, &cb, temp, dgp, &sac, &sbc);
2167 calc_dg_stddev(&ca, &cb, temp, dgp, &stddevc );
2170 dstddev2 += stddevc*stddevc;
2172 sample_coll_destroy(&ca);
2173 sample_coll_destroy(&cb);
2177 dg_sig2 += (dgs2-dgs*dgs)/(npee-1);
2183 sa_sig2 += (dsa2-dsa*dsa)/(npee-1);
2184 sb_sig2 += (dsb2-dsb*dsb)/(npee-1);
2188 stddev_sig2 += (dstddev2-dstddev*dstddev)/(npee-1);
2190 br->dg_err = sqrt(dg_sig2/(npee_max - npee_min + 1));
2191 br->sa_err = sqrt(sa_sig2/(npee_max - npee_min + 1));
2192 br->sb_err = sqrt(sb_sig2/(npee_max - npee_min + 1));
2193 br->dg_stddev_err = sqrt(stddev_sig2/(npee_max - npee_min + 1));
2198 static double bar_err(int nbmin, int nbmax, const double *partsum)
2201 double svar, s, s2, dg;
2204 for (nb = nbmin; nb <= nbmax; nb++)
2208 for (b = 0; b < nb; b++)
2210 dg = partsum[nb*(nbmax+1)+b];
2216 svar += (s2 - s*s)/(nb - 1);
2219 return sqrt(svar/(nbmax + 1 - nbmin));
2223 /* Seek the end of an identifier (consecutive non-spaces), followed by
2224 an optional number of spaces or '='-signs. Returns a pointer to the
2225 first non-space value found after that. Returns NULL if the string
2228 static const char *find_value(const char *str)
2230 gmx_bool name_end_found = FALSE;
2232 /* if the string is a NULL pointer, return a NULL pointer. */
2237 while (*str != '\0')
2239 /* first find the end of the name */
2240 if (!name_end_found)
2242 if (isspace(*str) || (*str == '=') )
2244 name_end_found = TRUE;
2249 if (!( isspace(*str) || (*str == '=') ))
2261 /* read a vector-notation description of a lambda vector */
2262 static gmx_bool read_lambda_compvec(const char *str,
2264 const lambda_components_t *lc_in,
2265 lambda_components_t *lc_out,
2269 gmx_bool initialize_lc = FALSE; /* whether to initialize the lambda
2270 components, or to check them */
2271 gmx_bool start_reached = FALSE; /* whether the start of component names
2273 gmx_bool vector = FALSE; /* whether there are multiple components */
2274 int n = 0; /* current component number */
2275 const char *val_start = NULL; /* start of the component name, or NULL
2276 if not in a value */
2286 if (lc_out && lc_out->N == 0)
2288 initialize_lc = TRUE;
2303 start_reached = TRUE;
2306 else if (*str == '(')
2309 start_reached = TRUE;
2311 else if (!isspace(*str))
2313 gmx_fatal(FARGS, "Error in lambda components in %s", fn);
2320 if (isspace(*str) || *str == ')' || *str == ',' || *str == '\0')
2327 lambda_components_add(lc_out, val_start,
2332 if (!lambda_components_check(lc_out, n, val_start,
2341 /* add a vector component to lv */
2342 lv->val[n] = strtod(val_start, &strtod_end);
2343 if (val_start == strtod_end)
2346 "Error reading lambda vector in %s", fn);
2349 /* reset for the next identifier */
2358 else if (isalnum(*str))
2371 gmx_fatal(FARGS, "Error in lambda components in %s", fn);
2379 else if (lv == NULL)
2385 gmx_fatal(FARGS, "Incomplete lambda vector data in %s",
2405 gmx_fatal(FARGS, "Incomplete lambda components data in %s", fn);
2411 /* read and check the component names from a string */
2412 static gmx_bool read_lambda_components(const char *str,
2413 lambda_components_t *lc,
2417 return read_lambda_compvec(str, NULL, NULL, lc, end, fn);
2420 /* read an initialized lambda vector from a string */
2421 static gmx_bool read_lambda_vector(const char *str,
2426 return read_lambda_compvec(str, lv, lv->lc, NULL, end, fn);
2431 /* deduce lambda value from legend.
2433 legend = the legend string
2435 lam = the initialized lambda vector
2436 returns whether to use the data in this set.
2438 static gmx_bool legend2lambda(const char *fn,
2444 const char *ptr = NULL, *ptr2 = NULL;
2445 gmx_bool ok = FALSE;
2446 gmx_bool bdhdl = FALSE;
2447 const char *tostr = " to ";
2451 gmx_fatal(FARGS, "There is no legend in file '%s', can not deduce lambda", fn);
2454 /* look for the last 'to': */
2458 ptr2 = strstr(ptr2, tostr);
2465 while (ptr2 != NULL && *ptr2 != '\0');
2469 ptr += strlen(tostr)-1; /* and advance past that 'to' */
2473 /* look for the = sign */
2474 ptr = strrchr(legend, '=');
2477 /* otherwise look for the last space */
2478 ptr = strrchr(legend, ' ');
2482 if (strstr(legend, "dH"))
2487 else if (strchr(legend, 'D') != NULL && strchr(legend, 'H') != NULL)
2492 else /*if (strstr(legend, "pV"))*/
2503 gmx_fatal(FARGS, "There is no proper lambda legend in file '%s', can not deduce lambda", fn);
2507 ptr = find_value(ptr);
2508 if (!ptr || !read_lambda_vector(ptr, lam, NULL, fn))
2510 gmx_fatal(FARGS, "lambda vector '%s' %s faulty", legend, fn);
2519 ptr = strrchr(legend, '=');
2523 /* there must be a component name */
2527 gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2529 /* now backtrack to the start of the identifier */
2530 while (isspace(*ptr))
2536 gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2539 while (!isspace(*ptr))
2544 gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2548 strncpy(buf, ptr, (end-ptr));
2549 buf[(end-ptr)] = '\0';
2550 dhdl_index = lambda_components_find(lam->lc, ptr, (end-ptr));
2554 strncpy(buf, ptr, (end-ptr));
2555 buf[(end-ptr)] = '\0';
2557 "Did not find lambda component for '%s' in %s",
2566 "dhdl without component name with >1 lambda component in %s",
2571 lam->dhdl = dhdl_index;
2576 static gmx_bool subtitle2lambda(const char *subtitle, xvg_t *ba, const char *fn,
2577 lambda_components_t *lc)
2582 double native_lambda;
2586 /* first check for a state string */
2587 ptr = strstr(subtitle, "state");
2591 const char *val_end;
2593 /* the new 4.6 style lambda vectors */
2594 ptr = find_value(ptr);
2597 index = strtol(ptr, &end, 10);
2600 gmx_fatal(FARGS, "Incomplete state data in %s", fn);
2607 gmx_fatal(FARGS, "Incomplete state data in %s", fn);
2610 /* now find the lambda vector component names */
2611 while (*ptr != '(' && !isalnum(*ptr))
2617 "Incomplete lambda vector component data in %s", fn);
2622 if (!read_lambda_components(ptr, lc, &val_end, fn))
2625 "lambda vector components in %s don't match those previously read",
2628 ptr = find_value(val_end);
2631 gmx_fatal(FARGS, "Incomplete state data in %s", fn);
2634 lambda_vec_init(&(ba->native_lambda), lc);
2635 if (!read_lambda_vector(ptr, &(ba->native_lambda), NULL, fn))
2637 gmx_fatal(FARGS, "lambda vector in %s faulty", fn);
2639 ba->native_lambda.index = index;
2644 /* compatibility mode: check for lambda in other ways. */
2645 /* plain text lambda string */
2646 ptr = strstr(subtitle, "lambda");
2649 /* xmgrace formatted lambda string */
2650 ptr = strstr(subtitle, "\\xl\\f{}");
2654 /* xmgr formatted lambda string */
2655 ptr = strstr(subtitle, "\\8l\\4");
2659 ptr = strstr(ptr, "=");
2663 bFound = (sscanf(ptr+1, "%lf", &(native_lambda)) == 1);
2664 /* add the lambda component name as an empty string */
2667 if (!lambda_components_check(lc, 0, "", 0))
2670 "lambda vector components in %s don't match those previously read",
2676 lambda_components_add(lc, "", 0);
2678 lambda_vec_init(&(ba->native_lambda), lc);
2679 ba->native_lambda.val[0] = native_lambda;
2686 static void filename2lambda(const char *fn, xvg_t *ba)
2689 const char *ptr, *digitptr;
2693 /* go to the end of the path string and search backward to find the last
2694 directory in the path which has to contain the value of lambda
2696 while (ptr[1] != '\0')
2700 /* searching backward to find the second directory separator */
2705 if (ptr[0] != DIR_SEPARATOR && ptr[1] == DIR_SEPARATOR)
2713 /* save the last position of a digit between the last two
2714 separators = in the last dirname */
2715 if (dirsep > 0 && isdigit(*ptr))
2723 gmx_fatal(FARGS, "While trying to read the lambda value from the file path:"
2724 " last directory in the path '%s' does not contain a number", fn);
2726 if (digitptr[-1] == '-')
2730 lambda = strtod(digitptr, &endptr);
2731 if (endptr == digitptr)
2733 gmx_fatal(FARGS, "Malformed number in file path '%s'", fn);
2737 static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
2738 lambda_components_t *lc)
2741 char *subtitle, **legend, *ptr;
2743 gmx_bool native_lambda_read = FALSE;
2751 np = read_xvg_legend(fn, &ba->y, &ba->nset, &subtitle, &legend);
2754 gmx_fatal(FARGS, "File %s contains no usable data.", fn);
2756 /* Reorder the data */
2758 for (i = 1; i < ba->nset; i++)
2760 ba->y[i-1] = ba->y[i];
2764 snew(ba->np, ba->nset);
2765 for (i = 0; i < ba->nset; i++)
2771 if (subtitle != NULL)
2773 /* try to extract temperature */
2774 ptr = strstr(subtitle, "T =");
2778 if (sscanf(ptr, "%lf", &ba->temp) == 1)
2782 gmx_fatal(FARGS, "Found temperature of %g in file '%s'",
2792 gmx_fatal(FARGS, "Did not find a temperature in the subtitle in file '%s', use the -temp option of [TT]g_bar[tt]", fn);
2797 /* Try to deduce lambda from the subtitle */
2800 if (subtitle2lambda(subtitle, ba, fn, lc))
2802 native_lambda_read = TRUE;
2805 snew(ba->lambda, ba->nset);
2808 /* Check if we have a single set, no legend, nset=1 means t and dH/dl */
2811 if (!native_lambda_read)
2813 /* Deduce lambda from the file name */
2814 filename2lambda(fn, ba);
2815 native_lambda_read = TRUE;
2817 ba->lambda[0] = ba->native_lambda;
2821 gmx_fatal(FARGS, "File %s contains multiple sets but no legends, can not determine the lambda values", fn);
2826 for (i = 0; i < ba->nset; )
2828 gmx_bool use = FALSE;
2829 /* Read lambda from the legend */
2830 lambda_vec_init( &(ba->lambda[i]), lc );
2831 lambda_vec_copy( &(ba->lambda[i]), &(ba->native_lambda));
2832 use = legend2lambda(fn, legend[i], ba, &(ba->lambda[i]));
2835 lambda_vec_print(&(ba->lambda[i]), buf, FALSE);
2841 printf("%s: Ignoring set '%s'.\n", fn, legend[i]);
2842 for (j = i+1; j < ba->nset; j++)
2844 ba->y[j-1] = ba->y[j];
2845 legend[j-1] = legend[j];
2852 if (!native_lambda_read)
2854 gmx_fatal(FARGS, "File %s contains multiple sets but no indication of the native lambda", fn);
2859 for (i = 0; i < ba->nset-1; i++)
2867 static void read_bar_xvg(char *fn, real *temp, sim_data_t *sd)
2876 read_bar_xvg_lowlevel(fn, temp, barsim, &(sd->lc));
2878 if (barsim->nset < 1)
2880 gmx_fatal(FARGS, "File '%s' contains fewer than two columns", fn);
2883 if (!gmx_within_tol(*temp, barsim->temp, GMX_FLOAT_EPS) && (*temp > 0) )
2885 gmx_fatal(FARGS, "Temperature in file %s different from earlier files or setting\n", fn);
2887 *temp = barsim->temp;
2889 /* now create a series of samples_t */
2890 snew(s, barsim->nset);
2891 for (i = 0; i < barsim->nset; i++)
2893 samples_init(s+i, &(barsim->native_lambda), &(barsim->lambda[i]),
2894 barsim->temp, lambda_vec_same(&(barsim->native_lambda),
2895 &(barsim->lambda[i])),
2897 s[i].du = barsim->y[i];
2898 s[i].ndu = barsim->np[i];
2901 lambda_data_list_insert_sample(sd->lb, s+i);
2906 lambda_vec_print(s[0].native_lambda, buf, FALSE);
2907 printf("%s: %.1f - %.1f; lambda = %s\n dH/dl & foreign lambdas:\n",
2908 fn, s[0].t[0], s[0].t[s[0].ndu-1], buf);
2909 for (i = 0; i < barsim->nset; i++)
2911 lambda_vec_print(s[i].foreign_lambda, buf, TRUE);
2912 printf(" %s (%d pts)\n", buf, s[i].ndu);
2918 static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
2919 double start_time, double delta_time,
2920 lambda_vec_t *native_lambda, double temp,
2921 double *last_t, const char *filename)
2925 double old_foreign_lambda;
2926 lambda_vec_t *foreign_lambda;
2928 samples_t *s; /* convenience pointer */
2931 /* check the block types etc. */
2932 if ( (blk->nsub < 3) ||
2933 (blk->sub[0].type != xdr_datatype_int) ||
2934 (blk->sub[1].type != xdr_datatype_double) ||
2936 (blk->sub[2].type != xdr_datatype_float) &&
2937 (blk->sub[2].type != xdr_datatype_double)
2939 (blk->sub[0].nr < 1) ||
2940 (blk->sub[1].nr < 1) )
2943 "Unexpected/corrupted block data in file %s around time %g.",
2944 filename, start_time);
2947 snew(foreign_lambda, 1);
2948 lambda_vec_init(foreign_lambda, native_lambda->lc);
2949 lambda_vec_copy(foreign_lambda, native_lambda);
2950 type = blk->sub[0].ival[0];
2953 for (i = 0; i < native_lambda->lc->N; i++)
2955 foreign_lambda->val[i] = blk->sub[1].dval[i];
2960 if (blk->sub[0].nr > 1)
2962 foreign_lambda->dhdl = blk->sub[0].ival[1];
2966 foreign_lambda->dhdl = 0;
2972 /* initialize the samples structure if it's empty. */
2974 samples_init(*smp, native_lambda, foreign_lambda, temp,
2975 type == dhbtDHDL, filename);
2976 (*smp)->start_time = start_time;
2977 (*smp)->delta_time = delta_time;
2980 /* set convenience pointer */
2983 /* now double check */
2984 if (!lambda_vec_same(s->foreign_lambda, foreign_lambda) )
2986 char buf[STRLEN], buf2[STRLEN];
2987 lambda_vec_print(foreign_lambda, buf, FALSE);
2988 lambda_vec_print(s->foreign_lambda, buf2, FALSE);
2989 fprintf(stderr, "Got foreign lambda=%s, expected: %s\n", buf, buf2);
2990 gmx_fatal(FARGS, "Corrupted data in file %s around t=%g.",
2991 filename, start_time);
2994 /* make room for the data */
2995 if (s->ndu_alloc < (size_t)(s->ndu + blk->sub[2].nr) )
2997 s->ndu_alloc += (s->ndu_alloc < (size_t)blk->sub[2].nr) ?
2998 blk->sub[2].nr*2 : s->ndu_alloc;
2999 srenew(s->du_alloc, s->ndu_alloc);
3000 s->du = s->du_alloc;
3003 s->ndu += blk->sub[2].nr;
3004 s->ntot += blk->sub[2].nr;
3005 *ndu = blk->sub[2].nr;
3007 /* and copy the data*/
3008 for (j = 0; j < blk->sub[2].nr; j++)
3010 if (blk->sub[2].type == xdr_datatype_float)
3012 s->du[startj+j] = blk->sub[2].fval[j];
3016 s->du[startj+j] = blk->sub[2].dval[j];
3019 if (start_time + blk->sub[2].nr*delta_time > *last_t)
3021 *last_t = start_time + blk->sub[2].nr*delta_time;
3025 static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
3026 double start_time, double delta_time,
3027 lambda_vec_t *native_lambda, double temp,
3028 double *last_t, const char *filename)
3033 double old_foreign_lambda;
3034 lambda_vec_t *foreign_lambda;
3038 /* check the block types etc. */
3039 if ( (blk->nsub < 2) ||
3040 (blk->sub[0].type != xdr_datatype_double) ||
3041 (blk->sub[1].type != xdr_datatype_large_int) ||
3042 (blk->sub[0].nr < 2) ||
3043 (blk->sub[1].nr < 2) )
3046 "Unexpected/corrupted block data in file %s around time %g",
3047 filename, start_time);
3050 nhist = blk->nsub-2;
3058 "Unexpected/corrupted block data in file %s around time %g",
3059 filename, start_time);
3065 snew(foreign_lambda, 1);
3066 lambda_vec_init(foreign_lambda, native_lambda->lc);
3067 lambda_vec_copy(foreign_lambda, native_lambda);
3068 type = (int)(blk->sub[1].lval[1]);
3071 double old_foreign_lambda;
3073 old_foreign_lambda = blk->sub[0].dval[0];
3074 if (old_foreign_lambda >= 0)
3076 foreign_lambda->val[0] = old_foreign_lambda;
3077 if (foreign_lambda->lc->N > 1)
3080 "Single-component lambda in multi-component file %s",
3086 for (i = 0; i < native_lambda->lc->N; i++)
3088 foreign_lambda->val[i] = blk->sub[0].dval[i+2];
3094 if (foreign_lambda->lc->N > 1)
3096 if (blk->sub[1].nr < 3 + nhist)
3099 "Missing derivative coord in multi-component file %s",
3102 foreign_lambda->dhdl = blk->sub[1].lval[2 + nhist];
3106 foreign_lambda->dhdl = 0;
3110 samples_init(s, native_lambda, foreign_lambda, temp, type == dhbtDHDL,
3114 for (i = 0; i < nhist; i++)
3116 nbins[i] = blk->sub[i+2].nr;
3119 hist_init(s->hist, nhist, nbins);
3121 for (i = 0; i < nhist; i++)
3123 s->hist->x0[i] = blk->sub[1].lval[2+i];
3124 s->hist->dx[i] = blk->sub[0].dval[1];
3127 s->hist->dx[i] = -s->hist->dx[i];
3131 s->hist->start_time = start_time;
3132 s->hist->delta_time = delta_time;
3133 s->start_time = start_time;
3134 s->delta_time = delta_time;
3136 for (i = 0; i < nhist; i++)
3139 gmx_large_int_t sum = 0;
3141 for (j = 0; j < s->hist->nbin[i]; j++)
3143 int binv = (int)(blk->sub[i+2].ival[j]);
3145 s->hist->bin[i][j] = binv;
3158 gmx_fatal(FARGS, "Histogram counts don't match in %s",
3164 if (start_time + s->hist->sum*delta_time > *last_t)
3166 *last_t = start_time + s->hist->sum*delta_time;
3172 static void read_barsim_edr(char *fn, real *temp, sim_data_t *sd)
3178 gmx_enxnm_t *enm = NULL;
3179 double first_t = -1;
3181 samples_t **samples_rawdh = NULL; /* contains samples for raw delta_h */
3182 int *nhists = NULL; /* array to keep count & print at end */
3183 int *npts = NULL; /* array to keep count & print at end */
3184 lambda_vec_t **lambdas = NULL; /* array to keep count & print at end */
3185 lambda_vec_t *native_lambda;
3186 double end_time; /* the end time of the last batch of samples */
3188 lambda_vec_t start_lambda;
3190 fp = open_enx(fn, "r");
3191 do_enxnms(fp, &nre, &enm);
3194 snew(native_lambda, 1);
3195 start_lambda.lc = NULL;
3197 while (do_enx(fp, fr))
3199 /* count the data blocks */
3200 int nblocks_raw = 0;
3201 int nblocks_hist = 0;
3204 /* DHCOLL block information: */
3205 double start_time = 0, delta_time = 0, old_start_lambda = 0, delta_lambda = 0;
3208 /* count the blocks and handle collection information: */
3209 for (i = 0; i < fr->nblock; i++)
3211 if (fr->block[i].id == enxDHHIST)
3215 if (fr->block[i].id == enxDH)
3219 if (fr->block[i].id == enxDHCOLL)
3222 if ( (fr->block[i].nsub < 1) ||
3223 (fr->block[i].sub[0].type != xdr_datatype_double) ||
3224 (fr->block[i].sub[0].nr < 5))
3226 gmx_fatal(FARGS, "Unexpected block data in file %s", fn);
3229 /* read the data from the DHCOLL block */
3230 rtemp = fr->block[i].sub[0].dval[0];
3231 start_time = fr->block[i].sub[0].dval[1];
3232 delta_time = fr->block[i].sub[0].dval[2];
3233 old_start_lambda = fr->block[i].sub[0].dval[3];
3234 delta_lambda = fr->block[i].sub[0].dval[4];
3236 if (delta_lambda > 0)
3238 gmx_fatal(FARGS, "Lambda values not constant in %s: can't apply BAR method", fn);
3240 if ( ( *temp != rtemp) && (*temp > 0) )
3242 gmx_fatal(FARGS, "Temperature in file %s different from earlier files or setting\n", fn);
3246 if (old_start_lambda >= 0)
3250 if (!lambda_components_check(&(sd->lc), 0, "", 0))
3253 "lambda vector components in %s don't match those previously read",
3259 lambda_components_add(&(sd->lc), "", 0);
3261 if (!start_lambda.lc)
3263 lambda_vec_init(&start_lambda, &(sd->lc));
3265 start_lambda.val[0] = old_start_lambda;
3269 /* read lambda vector */
3271 gmx_bool check = (sd->lc.N > 0);
3272 if (fr->block[i].nsub < 2)
3275 "No lambda vector, but start_lambda=%g\n",
3278 n_lambda_vec = fr->block[i].sub[1].ival[1];
3279 for (j = 0; j < n_lambda_vec; j++)
3282 efpt_singular_names[fr->block[i].sub[1].ival[1+j]];
3285 /* check the components */
3286 lambda_components_check(&(sd->lc), j, name,
3291 lambda_components_add(&(sd->lc), name,
3295 lambda_vec_init(&start_lambda, &(sd->lc));
3296 start_lambda.index = fr->block[i].sub[1].ival[0];
3297 for (j = 0; j < n_lambda_vec; j++)
3299 start_lambda.val[j] = fr->block[i].sub[0].dval[5+j];
3304 first_t = start_time;
3311 gmx_fatal(FARGS, "Did not find delta H information in file %s", fn);
3313 if (nblocks_raw > 0 && nblocks_hist > 0)
3315 gmx_fatal(FARGS, "Can't handle both raw delta U data and histograms in the same file %s", fn);
3320 /* check the native lambda */
3321 if (!lambda_vec_same(&start_lambda, native_lambda) )
3323 gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %g, and becomes %g at time %g",
3324 fn, native_lambda, start_lambda, start_time);
3326 /* check the number of samples against the previous number */
3327 if ( ((nblocks_raw+nblocks_hist) != nsamples) || (nlam != 1 ) )
3329 gmx_fatal(FARGS, "Unexpected block count in %s: was %d, now %d\n",
3330 fn, nsamples+1, nblocks_raw+nblocks_hist+nlam);
3332 /* check whether last iterations's end time matches with
3333 the currrent start time */
3334 if ( (fabs(last_t - start_time) > 2*delta_time) && last_t >= 0)
3336 /* it didn't. We need to store our samples and reallocate */
3337 for (i = 0; i < nsamples; i++)
3339 if (samples_rawdh[i])
3341 /* insert it into the existing list */
3342 lambda_data_list_insert_sample(sd->lb,
3344 /* and make sure we'll allocate a new one this time
3346 samples_rawdh[i] = NULL;
3353 /* this is the first round; allocate the associated data
3355 /*native_lambda=start_lambda;*/
3356 lambda_vec_init(native_lambda, &(sd->lc));
3357 lambda_vec_copy(native_lambda, &start_lambda);
3358 nsamples = nblocks_raw+nblocks_hist;
3359 snew(nhists, nsamples);
3360 snew(npts, nsamples);
3361 snew(lambdas, nsamples);
3362 snew(samples_rawdh, nsamples);
3363 for (i = 0; i < nsamples; i++)
3368 samples_rawdh[i] = NULL; /* init to NULL so we know which
3369 ones contain values */
3374 k = 0; /* counter for the lambdas, etc. arrays */
3375 for (i = 0; i < fr->nblock; i++)
3377 if (fr->block[i].id == enxDH)
3379 int type = (fr->block[i].sub[0].ival[0]);
3380 if (type == dhbtDH || type == dhbtDHDL)
3383 read_edr_rawdh_block(&(samples_rawdh[k]),
3386 start_time, delta_time,
3387 native_lambda, rtemp,
3390 if (samples_rawdh[k])
3392 lambdas[k] = samples_rawdh[k]->foreign_lambda;
3397 else if (fr->block[i].id == enxDHHIST)
3399 int type = (int)(fr->block[i].sub[1].lval[1]);
3400 if (type == dhbtDH || type == dhbtDHDL)
3404 samples_t *s; /* this is where the data will go */
3405 s = read_edr_hist_block(&nb, &(fr->block[i]),
3406 start_time, delta_time,
3407 native_lambda, rtemp,
3412 lambdas[k] = s->foreign_lambda;
3415 /* and insert the new sample immediately */
3416 for (j = 0; j < nb; j++)
3418 lambda_data_list_insert_sample(sd->lb, s+j);
3424 /* Now store all our extant sample collections */
3425 for (i = 0; i < nsamples; i++)
3427 if (samples_rawdh[i])
3429 /* insert it into the existing list */
3430 lambda_data_list_insert_sample(sd->lb, samples_rawdh[i]);
3438 lambda_vec_print(native_lambda, buf, FALSE);
3439 printf("%s: %.1f - %.1f; lambda = %s\n foreign lambdas:\n",
3440 fn, first_t, last_t, buf);
3441 for (i = 0; i < nsamples; i++)
3445 lambda_vec_print(lambdas[i], buf, TRUE);
3448 printf(" %s (%d hists)\n", buf, nhists[i]);
3452 printf(" %s (%d pts)\n", buf, npts[i]);
3464 int gmx_bar(int argc, char *argv[])
3466 static const char *desc[] = {
3467 "[TT]g_bar[tt] calculates free energy difference estimates through ",
3468 "Bennett's acceptance ratio method (BAR). It also automatically",
3469 "adds series of individual free energies obtained with BAR into",
3470 "a combined free energy estimate.[PAR]",
3472 "Every individual BAR free energy difference relies on two ",
3473 "simulations at different states: say state A and state B, as",
3474 "controlled by a parameter, [GRK]lambda[grk] (see the [TT].mdp[tt] parameter",
3475 "[TT]init_lambda[tt]). The BAR method calculates a ratio of weighted",
3476 "average of the Hamiltonian difference of state B given state A and",
3478 "The energy differences to the other state must be calculated",
3479 "explicitly during the simulation. This can be done with",
3480 "the [TT].mdp[tt] option [TT]foreign_lambda[tt].[PAR]",
3482 "Input option [TT]-f[tt] expects multiple [TT]dhdl.xvg[tt] files. ",
3483 "Two types of input files are supported:[BR]",
3484 "[TT]*[tt] Files with more than one [IT]y[it]-value. ",
3485 "The files should have columns ",
3486 "with dH/d[GRK]lambda[grk] and [GRK]Delta[grk][GRK]lambda[grk]. ",
3487 "The [GRK]lambda[grk] values are inferred ",
3488 "from the legends: [GRK]lambda[grk] of the simulation from the legend of ",
3489 "dH/d[GRK]lambda[grk] and the foreign [GRK]lambda[grk] values from the ",
3490 "legends of Delta H",
3492 "[TT]*[tt] Files with only one [IT]y[it]-value. Using the",
3493 "[TT]-extp[tt] option for these files, it is assumed",
3494 "that the [IT]y[it]-value is dH/d[GRK]lambda[grk] and that the ",
3495 "Hamiltonian depends linearly on [GRK]lambda[grk]. ",
3496 "The [GRK]lambda[grk] value of the simulation is inferred from the ",
3497 "subtitle (if present), otherwise from a number in the subdirectory ",
3498 "in the file name.[PAR]",
3500 "The [GRK]lambda[grk] of the simulation is parsed from ",
3501 "[TT]dhdl.xvg[tt] file's legend containing the string 'dH', the ",
3502 "foreign [GRK]lambda[grk] values from the legend containing the ",
3503 "capitalized letters 'D' and 'H'. The temperature is parsed from ",
3504 "the legend line containing 'T ='.[PAR]",
3506 "The input option [TT]-g[tt] expects multiple [TT].edr[tt] files. ",
3507 "These can contain either lists of energy differences (see the ",
3508 "[TT].mdp[tt] option [TT]separate_dhdl_file[tt]), or a series of ",
3509 "histograms (see the [TT].mdp[tt] options [TT]dh_hist_size[tt] and ",
3510 "[TT]dh_hist_spacing[tt]).", "The temperature and [GRK]lambda[grk] ",
3511 "values are automatically deduced from the [TT]ener.edr[tt] file.[PAR]",
3513 "In addition to the [TT].mdp[tt] option [TT]foreign_lambda[tt], ",
3514 "the energy difference can also be extrapolated from the ",
3515 "dH/d[GRK]lambda[grk] values. This is done with the[TT]-extp[tt]",
3516 "option, which assumes that the system's Hamiltonian depends linearly",
3517 "on [GRK]lambda[grk], which is not normally the case.[PAR]",
3519 "The free energy estimates are determined using BAR with bisection, ",
3520 "with the precision of the output set with [TT]-prec[tt]. ",
3521 "An error estimate taking into account time correlations ",
3522 "is made by splitting the data into blocks and determining ",
3523 "the free energy differences over those blocks and assuming ",
3524 "the blocks are independent. ",
3525 "The final error estimate is determined from the average variance ",
3526 "over 5 blocks. A range of block numbers for error estimation can ",
3527 "be provided with the options [TT]-nbmin[tt] and [TT]-nbmax[tt].[PAR]",
3529 "[TT]g_bar[tt] tries to aggregate samples with the same 'native' and ",
3530 "'foreign' [GRK]lambda[grk] values, but always assumes independent ",
3531 "samples. [BB]Note[bb] that when aggregating energy ",
3532 "differences/derivatives with different sampling intervals, this is ",
3533 "almost certainly not correct. Usually subsequent energies are ",
3534 "correlated and different time intervals mean different degrees ",
3535 "of correlation between samples.[PAR]",
3537 "The results are split in two parts: the last part contains the final ",
3538 "results in kJ/mol, together with the error estimate for each part ",
3539 "and the total. The first part contains detailed free energy ",
3540 "difference estimates and phase space overlap measures in units of ",
3541 "kT (together with their computed error estimate). The printed ",
3543 "[TT]*[tt] lam_A: the [GRK]lambda[grk] values for point A.[BR]",
3544 "[TT]*[tt] lam_B: the [GRK]lambda[grk] values for point B.[BR]",
3545 "[TT]*[tt] DG: the free energy estimate.[BR]",
3546 "[TT]*[tt] s_A: an estimate of the relative entropy of B in A.[BR]",
3547 "[TT]*[tt] s_A: an estimate of the relative entropy of A in B.[BR]",
3548 "[TT]*[tt] stdev: an estimate expected per-sample standard deviation.[PAR]",
3550 "The relative entropy of both states in each other's ensemble can be ",
3551 "interpreted as a measure of phase space overlap: ",
3552 "the relative entropy s_A of the work samples of lambda_B in the ",
3553 "ensemble of lambda_A (and vice versa for s_B), is a ",
3554 "measure of the 'distance' between Boltzmann distributions of ",
3555 "the two states, that goes to zero for identical distributions. See ",
3556 "Wu & Kofke, J. Chem. Phys. 123 084109 (2005) for more information.",
3558 "The estimate of the expected per-sample standard deviation, as given ",
3559 "in Bennett's original BAR paper: Bennett, J. Comp. Phys. 22, p 245 (1976).",
3560 "Eq. 10 therein gives an estimate of the quality of sampling (not directly",
3561 "of the actual statistical error, because it assumes independent samples).[PAR]",
3563 "To get a visual estimate of the phase space overlap, use the ",
3564 "[TT]-oh[tt] option to write series of histograms, together with the ",
3565 "[TT]-nbin[tt] option.[PAR]"
3567 static real begin = 0, end = -1, temp = -1;
3568 int nd = 2, nbmin = 5, nbmax = 5;
3570 gmx_bool use_dhdl = FALSE;
3571 gmx_bool calc_s, calc_v;
3573 { "-b", FALSE, etREAL, {&begin}, "Begin time for BAR" },
3574 { "-e", FALSE, etREAL, {&end}, "End time for BAR" },
3575 { "-temp", FALSE, etREAL, {&temp}, "Temperature (K)" },
3576 { "-prec", FALSE, etINT, {&nd}, "The number of digits after the decimal point" },
3577 { "-nbmin", FALSE, etINT, {&nbmin}, "Minimum number of blocks for error estimation" },
3578 { "-nbmax", FALSE, etINT, {&nbmax}, "Maximum number of blocks for error estimation" },
3579 { "-nbin", FALSE, etINT, {&nbin}, "Number of bins for histogram output"},
3580 { "-extp", FALSE, etBOOL, {&use_dhdl}, "Whether to linearly extrapolate dH/dl values to use as energies"}
3584 { efXVG, "-f", "dhdl", ffOPTRDMULT },
3585 { efEDR, "-g", "ener", ffOPTRDMULT },
3586 { efXVG, "-o", "bar", ffOPTWR },
3587 { efXVG, "-oi", "barint", ffOPTWR },
3588 { efXVG, "-oh", "histogram", ffOPTWR }
3590 #define NFILE asize(fnm)
3593 int nf = 0; /* file counter */
3595 int nfile_tot; /* total number of input files */
3600 sim_data_t sim_data; /* the simulation data */
3601 barres_t *results; /* the results */
3602 int nresults; /* number of results in results array */
3605 double prec, dg_tot, dg, sig, dg_tot_max, dg_tot_min;
3607 char dgformat[20], xvg2format[STRLEN], xvg3format[STRLEN];
3608 char buf[STRLEN], buf2[STRLEN];
3609 char ktformat[STRLEN], sktformat[STRLEN];
3610 char kteformat[STRLEN], skteformat[STRLEN];
3613 gmx_bool result_OK = TRUE, bEE = TRUE;
3615 gmx_bool disc_err = FALSE;
3616 double sum_disc_err = 0.; /* discretization error */
3617 gmx_bool histrange_err = FALSE;
3618 double sum_histrange_err = 0.; /* histogram range error */
3619 double stat_err = 0.; /* statistical error */
3621 parse_common_args(&argc, argv,
3623 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
3625 if (opt2bSet("-f", NFILE, fnm))
3627 nxvgfile = opt2fns(&fxvgnms, "-f", NFILE, fnm);
3629 if (opt2bSet("-g", NFILE, fnm))
3631 nedrfile = opt2fns(&fedrnms, "-g", NFILE, fnm);
3634 sim_data_init(&sim_data);
3636 /* make linked list */
3638 lambda_data_init(lb, 0, 0);
3644 nfile_tot = nxvgfile + nedrfile;
3648 gmx_fatal(FARGS, "No input files!");
3653 gmx_fatal(FARGS, "Can not have negative number of digits");
3655 prec = pow(10, -nd);
3657 snew(partsum, (nbmax+1)*(nbmax+1));
3660 /* read in all files. First xvg files */
3661 for (f = 0; f < nxvgfile; f++)
3663 read_bar_xvg(fxvgnms[f], &temp, &sim_data);
3666 /* then .edr files */
3667 for (f = 0; f < nedrfile; f++)
3669 read_barsim_edr(fedrnms[f], &temp, &sim_data);;
3673 /* fix the times to allow for equilibration */
3674 sim_data_impose_times(&sim_data, begin, end);
3676 if (opt2bSet("-oh", NFILE, fnm))
3678 sim_data_histogram(&sim_data, opt2fn("-oh", NFILE, fnm), nbin, oenv);
3681 /* assemble the output structures from the lambdas */
3682 results = barres_list_create(&sim_data, &nresults, use_dhdl);
3684 sum_disc_err = barres_list_max_disc_err(results, nresults);
3688 printf("\nNo results to calculate.\n");
3692 if (sum_disc_err > prec)
3694 prec = sum_disc_err;
3695 nd = ceil(-log10(prec));
3696 printf("WARNING: setting the precision to %g because that is the minimum\n reasonable number, given the expected discretization error.\n", prec);
3700 /*sprintf(lamformat,"%%6.3f");*/
3701 sprintf( dgformat, "%%%d.%df", 3+nd, nd);
3702 /* the format strings of the results in kT */
3703 sprintf( ktformat, "%%%d.%df", 5+nd, nd);
3704 sprintf( sktformat, "%%%ds", 6+nd);
3705 /* the format strings of the errors in kT */
3706 sprintf( kteformat, "%%%d.%df", 3+nd, nd);
3707 sprintf( skteformat, "%%%ds", 4+nd);
3708 sprintf(xvg2format, "%s %s\n", "%s", dgformat);
3709 sprintf(xvg3format, "%s %s %s\n", "%s", dgformat, dgformat);
3714 if (opt2bSet("-o", NFILE, fnm))
3716 sprintf(buf, "%s (%s)", "\\DeltaG", "kT");
3717 fpb = xvgropen_type(opt2fn("-o", NFILE, fnm), "Free energy differences",
3718 "\\lambda", buf, exvggtXYDY, oenv);
3722 if (opt2bSet("-oi", NFILE, fnm))
3724 sprintf(buf, "%s (%s)", "\\DeltaG", "kT");
3725 fpi = xvgropen(opt2fn("-oi", NFILE, fnm), "Free energy integral",
3726 "\\lambda", buf, oenv);
3736 /* first calculate results */
3739 for (f = 0; f < nresults; f++)
3741 /* Determine the free energy difference with a factor of 10
3742 * more accuracy than requested for printing.
3744 calc_bar(&(results[f]), 0.1*prec, nbmin, nbmax,
3747 if (results[f].dg_disc_err > prec/10.)
3751 if (results[f].dg_histrange_err > prec/10.)
3753 histrange_err = TRUE;
3757 /* print results in kT */
3761 printf("\nTemperature: %g K\n", temp);
3763 printf("\nDetailed results in kT (see help for explanation):\n\n");
3764 printf("%6s ", " lam_A");
3765 printf("%6s ", " lam_B");
3766 printf(sktformat, "DG ");
3769 printf(skteformat, "+/- ");
3773 printf(skteformat, "disc ");
3777 printf(skteformat, "range ");
3779 printf(sktformat, "s_A ");
3782 printf(skteformat, "+/- " );
3784 printf(sktformat, "s_B ");
3787 printf(skteformat, "+/- " );
3789 printf(sktformat, "stdev ");
3792 printf(skteformat, "+/- ");
3795 for (f = 0; f < nresults; f++)
3797 lambda_vec_print_short(results[f].a->native_lambda, buf);
3799 lambda_vec_print_short(results[f].b->native_lambda, buf);
3801 printf(ktformat, results[f].dg);
3805 printf(kteformat, results[f].dg_err);
3810 printf(kteformat, results[f].dg_disc_err);
3815 printf(kteformat, results[f].dg_histrange_err);
3818 printf(ktformat, results[f].sa);
3822 printf(kteformat, results[f].sa_err);
3825 printf(ktformat, results[f].sb);
3829 printf(kteformat, results[f].sb_err);
3832 printf(ktformat, results[f].dg_stddev);
3836 printf(kteformat, results[f].dg_stddev_err);
3840 /* Check for negative relative entropy with a 95% certainty. */
3841 if (results[f].sa < -2*results[f].sa_err ||
3842 results[f].sb < -2*results[f].sb_err)
3850 printf("\nWARNING: Some of these results violate the Second Law of "
3851 "Thermodynamics: \n"
3852 " This is can be the result of severe undersampling, or "
3854 " there is something wrong with the simulations.\n");
3858 /* final results in kJ/mol */
3859 printf("\n\nFinal results in kJ/mol:\n\n");
3861 for (f = 0; f < nresults; f++)
3866 lambda_vec_print_short(results[f].a->native_lambda, buf);
3867 fprintf(fpi, xvg2format, buf, dg_tot);
3873 lambda_vec_print_intermediate(results[f].a->native_lambda,
3874 results[f].b->native_lambda,
3877 fprintf(fpb, xvg3format, buf, results[f].dg, results[f].dg_err);
3881 lambda_vec_print_short(results[f].a->native_lambda, buf);
3882 lambda_vec_print_short(results[f].b->native_lambda, buf2);
3883 printf("%s - %s", buf, buf2);
3886 printf(dgformat, results[f].dg*kT);
3890 printf(dgformat, results[f].dg_err*kT);
3894 printf(" (max. range err. = ");
3895 printf(dgformat, results[f].dg_histrange_err*kT);
3897 sum_histrange_err += results[f].dg_histrange_err*kT;
3901 dg_tot += results[f].dg;
3905 lambda_vec_print_short(results[0].a->native_lambda, buf);
3906 lambda_vec_print_short(results[nresults-1].b->native_lambda, buf2);
3907 printf("%s - %s", buf, buf2);
3910 printf(dgformat, dg_tot*kT);
3913 stat_err = bar_err(nbmin, nbmax, partsum)*kT;
3915 printf(dgformat, max(max(stat_err, sum_disc_err), sum_histrange_err));
3920 printf("\nmaximum discretization error = ");
3921 printf(dgformat, sum_disc_err);
3922 if (bEE && stat_err < sum_disc_err)
3924 printf("WARNING: discretization error (%g) is larger than statistical error.\n Decrease histogram spacing for more accurate results\n", stat_err);
3929 printf("\nmaximum histogram range error = ");
3930 printf(dgformat, sum_histrange_err);
3931 if (bEE && stat_err < sum_histrange_err)
3933 printf("WARNING: histogram range error (%g) is larger than statistical error.\n Increase histogram range for more accurate results\n", stat_err);
3942 lambda_vec_print_short(results[nresults-1].b->native_lambda, buf);
3943 fprintf(fpi, xvg2format, buf, dg_tot);
3947 do_view(oenv, opt2fn_null("-o", NFILE, fnm), "-xydy");
3948 do_view(oenv, opt2fn_null("-oi", NFILE, fnm), "-xydy");