src/tools/gmx_bar.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  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.
  15
  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.
  20  *
  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.
  27  *
  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.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Green Red Orange Magenta Azure Cyan Skyblue
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <string.h>
  42 #include <ctype.h>
  43 #include <math.h>
  44 #include <float.h>
  45
  46 #include "sysstuff.h"
  47 #include "typedefs.h"
  48 #include "smalloc.h"
  49 #include "futil.h"
  50 #include "statutil.h"
  51 #include "copyrite.h"
  52 #include "macros.h"
  53 #include "enxio.h"
  54 #include "physics.h"
  55 #include "gmx_fatal.h"
  56 #include "xvgr.h"
  57 #include "gmx_ana.h"
  58 #include "maths.h"
  59 #include "string2.h"
  60
  61 /* the dhdl.xvg data from a simulation (actually obsolete, but still
  62     here for reading the dhdl.xvg file*/
  63 typedef struct xvg_t
  64 {
  65     char   *filename;
  66     int    ftp;     /* file type */
  67     int    nset;    /* number of lambdas, including dhdl */
  68     int *np;        /* number of data points (du or hists) per lambda */
  69     int  np_alloc;  /* number of points (du or hists) allocated */
  70     double temp;    /* temperature */
  71     double *lambda; /* the lambdas (of first index for y). */
  72     double *t;      /* the times (of second index for y) */
  73     double **y;     /* the dU values. y[0] holds the derivative, while
  74                        further ones contain the energy differences between
  75                        the native lambda and the 'foreign' lambdas. */
  76
  77     double native_lambda; /* the native lambda */
  78
  79     struct xvg_t *next, *prev; /*location in the global linked list of xvg_ts*/
  80 } xvg_t;
  81
  82
  83
  84 typedef struct hist_t
  85 {
  86     unsigned int *bin[2];       /* the (forward + reverse) histogram values */
  87     double dx[2];               /* the histogram spacing. The reverse
  88                                    dx is the negative of the forward dx.*/
  89     gmx_large_int_t x0[2];      /* the (forward + reverse) histogram start
  90                                    point(s) as int */
  91
  92     int nbin[2];                /* the (forward+reverse) number of bins */
  93     gmx_large_int_t sum;        /* the total number of counts. Must be
  94                                    the same for forward + reverse.  */
  95     int nhist;                  /* number of hist datas (forward or reverse) */
  96
  97     double start_time, delta_time;  /* start time, end time of histogram */
  98 } hist_t;
  99
 100
 101 /* an aggregate of samples for partial free energy calculation */
 102 typedef struct samples_t
 103 {
 104     double native_lambda;
 105     double foreign_lambda;
 106     double temp; /* the temperature */
 107     gmx_bool derivative; /* whether this sample is a derivative */
 108
 109     /* The samples come either as either delta U lists: */
 110     int ndu; /* the number of delta U samples */
 111     double *du; /* the delta u's */
 112     double *t; /* the times associated with those samples, or: */
 113     double start_time,delta_time;/*start time and delta time for linear time*/
 114
 115     /* or as histograms: */
 116     hist_t *hist; /* a histogram */
 117
 118     /* allocation data: (not NULL for data 'owned' by this struct) */
 119     double *du_alloc, *t_alloc; /* allocated delta u arrays  */
 120     size_t ndu_alloc, nt_alloc; /* pre-allocated sizes */
 121     hist_t *hist_alloc; /* allocated hist */
 122
 123     gmx_large_int_t ntot; /* total number of samples */
 124     const char *filename; /* the file name this sample comes from */
 125 } samples_t;
 126
 127 /* a sample range (start to end for du-style data, or boolean
 128     for both du-style data and histograms */
 129 typedef struct sample_range_t
 130 {
 131     int start, end; /* start and end index for du style data */
 132     gmx_bool use; /* whether to use this sample */
 133
 134     samples_t *s; /* the samples this range belongs to */
 135 } sample_range_t;
 136
 137
 138 /* a collection of samples for a partial free energy calculation
 139     (i.e. the collection of samples from one native lambda to one
 140     foreign lambda) */
 141 typedef struct sample_coll_t
 142 {
 143     double native_lambda;  /* these should be the same for all samples in the histogram?*/
 144     double foreign_lambda; /* collection */
 145     double temp; /* the temperature */
 146
 147     int nsamples; /* the number of samples */
 148     samples_t **s; /* the samples themselves */
 149     sample_range_t *r; /* the sample ranges */
 150     int nsamples_alloc; /* number of allocated samples */
 151
 152     gmx_large_int_t ntot; /* total number of samples in the ranges of
 153                              this collection */
 154
 155     struct sample_coll_t *next, *prev; /* next and previous in the list */
 156 } sample_coll_t;
 157
 158 /* all the samples associated with a lambda point */
 159 typedef struct lambda_t
 160 {
 161     double lambda; /* the native lambda (at start time if dynamic) */
 162     double temp; /* temperature */
 163
 164     sample_coll_t *sc; /* the samples */
 165
 166     sample_coll_t sc_head; /*the pre-allocated list head for the linked list.*/
 167
 168     struct lambda_t *next, *prev; /* the next and prev in the list */
 169 } lambda_t;
 170
 171
 172 /* calculated values. */
 173 typedef struct {
 174     sample_coll_t *a, *b; /* the simulation data */
 175
 176     double dg; /* the free energy difference */
 177     double dg_err; /* the free energy difference */
 178
 179     double dg_disc_err; /* discretization error */
 180     double dg_histrange_err; /* histogram range error */
 181
 182     double sa; /* relative entropy of b in state a */
 183     double sa_err; /* error in sa */
 184     double sb; /* relative entropy of a in state b */
 185     double sb_err; /* error in sb */
 186
 187     double dg_stddev; /* expected dg stddev per sample */
 188     double dg_stddev_err; /* error in dg_stddev */
 189 } barres_t;
 190
 191
 192
 193
 194 static void hist_init(hist_t *h, int nhist, int *nbin)
 195 {
 196     int i;
 197     if (nhist>2)
 198     {
 199         gmx_fatal(FARGS, "histogram with more than two sets of data!");
 200     }
 201     for(i=0;i<nhist;i++)
 202     {
 203         snew(h->bin[i], nbin[i]);
 204         h->x0[i]=0;
 205         h->nbin[i]=nbin[i];
 206         h->start_time=h->delta_time=0;
 207         h->dx[i]=0;
 208     }
 209     h->sum=0;
 210     h->nhist=nhist;
 211 }
 212
 213 static void hist_destroy(hist_t *h)
 214 {
 215     sfree(h->bin);
 216 }
 217
 218
 219 static void xvg_init(xvg_t *ba)
 220 {
 221     ba->filename=NULL;
 222     ba->nset=0;
 223     ba->np_alloc=0;
 224     ba->np=NULL;
 225     ba->y=NULL;
 226 }
 227
 228 static void samples_init(samples_t *s, double native_lambda,
 229                          double foreign_lambda, double temp,
 230                          gmx_bool derivative, const char *filename)
 231 {
 232     s->native_lambda=native_lambda;
 233     s->foreign_lambda=foreign_lambda;
 234     s->temp=temp;
 235     s->derivative=derivative;
 236
 237     s->ndu=0;
 238     s->du=NULL;
 239     s->t=NULL;
 240     s->start_time = s->delta_time = 0;
 241     s->hist=NULL;
 242     s->du_alloc=NULL;
 243     s->t_alloc=NULL;
 244     s->hist_alloc=NULL;
 245     s->ndu_alloc=0;
 246     s->nt_alloc=0;
 247
 248     s->ntot=0;
 249     s->filename=filename;
 250 }
 251
 252 static void sample_range_init(sample_range_t *r, samples_t *s)
 253 {
 254     r->start=0;
 255     r->end=s->ndu;
 256     r->use=TRUE;
 257     r->s=NULL;
 258 }
 259
 260 static void sample_coll_init(sample_coll_t *sc, double native_lambda,
 261                              double foreign_lambda, double temp)
 262 {
 263     sc->native_lambda = native_lambda;
 264     sc->foreign_lambda = foreign_lambda;
 265     sc->temp = temp;
 266
 267     sc->nsamples=0;
 268     sc->s=NULL;
 269     sc->r=NULL;
 270     sc->nsamples_alloc=0;
 271
 272     sc->ntot=0;
 273     sc->next=sc->prev=NULL;
 274 }
 275
 276 static void sample_coll_destroy(sample_coll_t *sc)
 277 {
 278     /* don't free the samples themselves */
 279     sfree(sc->r);
 280     sfree(sc->s);
 281 }
 282
 283
 284 static void lambda_init(lambda_t *l, double native_lambda, double temp)
 285 {
 286     l->lambda=native_lambda;
 287     l->temp=temp;
 288
 289     l->next=NULL;
 290     l->prev=NULL;
 291
 292     l->sc=&(l->sc_head);
 293
 294     sample_coll_init(l->sc, native_lambda, 0., 0.);
 295     l->sc->next=l->sc;
 296     l->sc->prev=l->sc;
 297 }
 298
 299 static void barres_init(barres_t *br)
 300 {
 301     br->dg=0;
 302     br->dg_err=0;
 303     br->sa=0;
 304     br->sa_err=0;
 305     br->sb=0;
 306     br->sb_err=0;
 307     br->dg_stddev=0;
 308     br->dg_stddev_err=0;
 309
 310     br->a=NULL;
 311     br->b=NULL;
 312 }
 313
 314
 315
 316
 317 static gmx_bool lambda_same(double lambda1, double lambda2)
 318 {
 319     return gmx_within_tol(lambda1, lambda2, 10*GMX_REAL_EPS);
 320 }
 321
 322 /* calculate the total number of samples in a sample collection */
 323 static void sample_coll_calc_ntot(sample_coll_t *sc)
 324 {
 325     int i;
 326
 327     sc->ntot=0;
 328     for(i=0;i<sc->nsamples;i++)
 329     {
 330         if (sc->r[i].use)
 331         {
 332             if (sc->s[i]->hist)
 333             {
 334                 sc->ntot += sc->s[i]->ntot;
 335             }
 336             else
 337             {
 338                 sc->ntot += sc->r[i].end - sc->r[i].start;
 339             }
 340         }
 341     }
 342 }
 343
 344
 345 /* find the barsamples_t associated with a lambda that corresponds to
 346    a specific foreign lambda */
 347 static sample_coll_t *lambda_find_sample_coll(lambda_t *l,
 348                                               double foreign_lambda)
 349 {
 350     sample_coll_t *sc=l->sc->next;
 351
 352     while(sc != l->sc)
 353     {
 354         if (lambda_same(sc->foreign_lambda,foreign_lambda))
 355         {
 356             return sc;
 357         }
 358         sc=sc->next;
 359     }
 360
 361     return NULL;
 362 }
 363
 364 /* insert li into an ordered list of lambda_colls */
 365 static void lambda_insert_sample_coll(lambda_t *l, sample_coll_t *sc)
 366 {
 367     sample_coll_t *scn=l->sc->next;
 368     while ( (scn!=l->sc) )
 369     {
 370         if (scn->foreign_lambda > sc->foreign_lambda)
 371             break;
 372         scn=scn->next;
 373     }
 374     /* now insert it before the found scn */
 375     sc->next=scn;
 376     sc->prev=scn->prev;
 377     scn->prev->next=sc;
 378     scn->prev=sc;
 379 }
 380
 381 /* insert li into an ordered list of lambdas */
 382 static void lambda_insert_lambda(lambda_t *head, lambda_t *li)
 383 {
 384     lambda_t *lc=head->next;
 385     while (lc!=head)
 386     {
 387         if (lc->lambda > li->lambda)
 388             break;
 389         lc=lc->next;
 390     }
 391     /* now insert ourselves before the found lc */
 392     li->next=lc;
 393     li->prev=lc->prev;
 394     lc->prev->next=li;
 395     lc->prev=li;
 396 }
 397
 398 /* insert a sample and a sample_range into a sample_coll. The
 399     samples are stored as a pointer, the range is copied. */
 400 static void sample_coll_insert_sample(sample_coll_t *sc, samples_t *s,
 401                                       sample_range_t *r)
 402 {
 403     /* first check if it belongs here */
 404     if (sc->temp != s->temp)
 405     {
 406         gmx_fatal(FARGS, "Temperatures in files %s and %s are not the same!",
 407                    s->filename, sc->next->s[0]->filename);
 408     }
 409     if (sc->native_lambda != s->native_lambda)
 410     {
 411         gmx_fatal(FARGS, "Native lambda in files %s and %s are not the same (and they should be)!",
 412                    s->filename, sc->next->s[0]->filename);
 413     }
 414     if (sc->foreign_lambda != s->foreign_lambda)
 415     {
 416         gmx_fatal(FARGS, "Foreign lambda in files %s and %s are not the same (and they should be)!",
 417                    s->filename, sc->next->s[0]->filename);
 418     }
 419
 420     /* check if there's room */
 421     if ( (sc->nsamples + 1) > sc->nsamples_alloc)
 422     {
 423         sc->nsamples_alloc = max(2*sc->nsamples_alloc, 2);
 424         srenew(sc->s, sc->nsamples_alloc);
 425         srenew(sc->r, sc->nsamples_alloc);
 426     }
 427     sc->s[sc->nsamples]=s;
 428     sc->r[sc->nsamples]=*r;
 429     sc->nsamples++;
 430
 431     sample_coll_calc_ntot(sc);
 432 }
 433
 434 /* insert a sample into a lambda_list, creating the right sample_coll if
 435    neccesary */
 436 static void lambda_list_insert_sample(lambda_t *head, samples_t *s)
 437 {
 438     gmx_bool found=FALSE;
 439     sample_coll_t *sc;
 440     sample_range_t r;
 441
 442     lambda_t *l=head->next;
 443
 444     /* first search for the right lambda_t */
 445     while(l != head)
 446     {
 447         if (lambda_same(l->lambda, s->native_lambda) )
 448         {
 449             found=TRUE;
 450             break;
 451         }
 452         l=l->next;
 453     }
 454
 455     if (!found)
 456     {
 457         snew(l, 1); /* allocate a new one */
 458         lambda_init(l, s->native_lambda, s->temp); /* initialize it */
 459         lambda_insert_lambda(head, l); /* add it to the list */
 460     }
 461
 462     /* now look for a sample collection */
 463     sc=lambda_find_sample_coll(l, s->foreign_lambda);
 464     if (!sc)
 465     {
 466         snew(sc, 1); /* allocate a new one */
 467         sample_coll_init(sc, s->native_lambda, s->foreign_lambda, s->temp);
 468         lambda_insert_sample_coll(l, sc);
 469     }
 470
 471     /* now insert the samples into the sample coll */
 472     sample_range_init(&r, s);
 473     sample_coll_insert_sample(sc, s, &r);
 474 }
 475
 476
 477 /* make a histogram out of a sample collection */
 478 static void sample_coll_make_hist(sample_coll_t *sc, int **bin,
 479                                   int *nbin_alloc, int *nbin,
 480                                   double *dx, double *xmin, int nbin_default)
 481 {
 482     int i,j,k;
 483     gmx_bool dx_set=FALSE;
 484     gmx_bool xmin_set=FALSE;
 485
 486     gmx_bool xmax_set=FALSE;
 487     gmx_bool xmax_set_hard=FALSE; /* whether the xmax is bounded by the
 488                                      limits of a histogram */
 489     double xmax=-1;
 490
 491     /* first determine dx and xmin; try the histograms */
 492     for(i=0;i<sc->nsamples;i++)
 493     {
 494         if (sc->s[i]->hist)
 495         {
 496             hist_t *hist=sc->s[i]->hist;
 497             for(k=0;k<hist->nhist;k++)
 498             {
 499                 double hdx=hist->dx[k];
 500                 double xmax_now=(hist->x0[k]+hist->nbin[k])*hdx;
 501
 502                 /* we use the biggest dx*/
 503                 if ( (!dx_set) || hist->dx[0] > *dx)
 504                 {
 505                     dx_set=TRUE;
 506                     *dx = hist->dx[0];
 507                 }
 508                 if ( (!xmin_set) || (hist->x0[k]*hdx) < *xmin)
 509                 {
 510                     xmin_set=TRUE;
 511                     *xmin = (hist->x0[k]*hdx);
 512                 }
 513
 514                 if ( (!xmax_set) || (xmax_now>xmax && !xmax_set_hard) )
 515                 {
 516                     xmax_set=TRUE;
 517                     xmax = xmax_now;
 518                     if (hist->bin[k][hist->nbin[k]-1] != 0)
 519                         xmax_set_hard=TRUE;
 520                 }
 521                 if ( hist->bin[k][hist->nbin[k]-1]!=0 && (xmax_now < xmax) )
 522                 {
 523                     xmax_set_hard=TRUE;
 524                     xmax = xmax_now;
 525                 }
 526             }
 527         }
 528     }
 529     /* and the delta us */
 530     for(i=0;i<sc->nsamples;i++)
 531     {
 532         if (sc->s[i]->ndu > 0)
 533         {
 534             /* determine min and max */
 535             int starti=sc->r[i].start;
 536             int endi=sc->r[i].end;
 537             double du_xmin=sc->s[i]->du[starti];
 538             double du_xmax=sc->s[i]->du[starti];
 539             for(j=starti+1;j<endi;j++)
 540             {
 541                 if (sc->s[i]->du[j] < du_xmin)
 542                     du_xmin = sc->s[i]->du[j];
 543                 if (sc->s[i]->du[j] > du_xmax)
 544                     du_xmax = sc->s[i]->du[j];
 545             }
 546
 547             /* and now change the limits */
 548             if ( (!xmin_set) || (du_xmin < *xmin) )
 549             {
 550                 xmin_set=TRUE;
 551                 *xmin=du_xmin;
 552             }
 553             if ( (!xmax_set) || ((du_xmax > xmax) &&  !xmax_set_hard) )
 554             {
 555                 xmax_set=TRUE;
 556                 xmax=du_xmax;
 557             }
 558         }
 559     }
 560
 561     if (!xmax_set || !xmin_set)
 562     {
 563         *nbin=0;
 564         return;
 565     }
 566
 567
 568     if (!dx_set)
 569     {
 570         *nbin=nbin_default;
 571         *dx=(xmax-(*xmin))/((*nbin)-2); /* -2 because we want the last bin to
 572                                            be 0, and we count from 0 */
 573     }
 574     else
 575     {
 576         *nbin=(xmax-(*xmin))/(*dx);
 577     }
 578
 579     if (*nbin > *nbin_alloc)
 580     {
 581         *nbin_alloc=*nbin;
 582         srenew(*bin, *nbin_alloc);
 583     }
 584
 585     /* reset the histogram */
 586     for(i=0;i<(*nbin);i++)
 587     {
 588         (*bin)[i] = 0;
 589     }
 590
 591     /* now add the actual data */
 592     for(i=0;i<sc->nsamples;i++)
 593     {
 594         if (sc->s[i]->hist)
 595         {
 596             hist_t *hist=sc->s[i]->hist;
 597             for(k=0;k<hist->nhist;k++)
 598             {
 599                 double hdx = hist->dx[k];
 600                 double xmin_hist=hist->x0[k]*hdx;
 601                 for(j=0;j<hist->nbin[k];j++)
 602                 {
 603                     /* calculate the bin corresponding to the middle of the
 604                        original bin */
 605                     double x=hdx*(j+0.5) + xmin_hist;
 606                     int binnr=(int)((x-(*xmin))/(*dx));
 607
 608                     if (binnr >= *nbin || binnr<0)
 609                         binnr = (*nbin)-1;
 610
 611                     (*bin)[binnr] += hist->bin[k][j];
 612                 }
 613             }
 614         }
 615         else
 616         {
 617             int starti=sc->r[i].start;
 618             int endi=sc->r[i].end;
 619             for(j=starti;j<endi;j++)
 620             {
 621                 int binnr=(int)((sc->s[i]->du[j]-(*xmin))/(*dx));
 622                 if (binnr >= *nbin || binnr<0)
 623                     binnr = (*nbin)-1;
 624
 625                 (*bin)[binnr] ++;
 626             }
 627         }
 628     }
 629 }
 630
 631 /* write a collection of histograms to a file */
 632 void lambdas_histogram(lambda_t *bl_head, const char *filename,
 633                        int nbin_default, const output_env_t oenv)
 634 {
 635     char label_x[STRLEN];
 636     const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
 637     const char *title="N(\\DeltaH)";
 638     const char *label_y="Samples";
 639     FILE *fp;
 640     lambda_t *bl;
 641     int nsets=0;
 642     char **setnames=NULL;
 643     gmx_bool first_set=FALSE;
 644     /* histogram data: */
 645     int *hist=NULL;
 646     int nbin=0;
 647     int nbin_alloc=0;
 648     double dx=0;
 649     double min=0;
 650     int i;
 651
 652     printf("\nWriting histogram to %s\n", filename);
 653     sprintf(label_x, "\\DeltaH (%s)", unit_energy);
 654
 655     fp=xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
 656
 657     /* first get all the set names */
 658     bl=bl_head->next;
 659     /* iterate over all lambdas */
 660     while(bl!=bl_head)
 661     {
 662         sample_coll_t *sc=bl->sc->next;
 663
 664         /* iterate over all samples */
 665         while(sc!=bl->sc)
 666         {
 667             nsets++;
 668             srenew(setnames, nsets);
 669             snew(setnames[nsets-1], STRLEN);
 670             if (!lambda_same(sc->foreign_lambda, sc->native_lambda))
 671             {
 672                 sprintf(setnames[nsets-1], "N(%s(%s=%g) | %s=%g)",
 673                         deltag, lambda, sc->foreign_lambda, lambda,
 674                         sc->native_lambda);
 675             }
 676             else
 677             {
 678                 sprintf(setnames[nsets-1], "N(%s | %s=%g)",
 679                         dhdl, lambda, sc->native_lambda);
 680             }
 681             sc=sc->next;
 682         }
 683
 684         bl=bl->next;
 685     }
 686     xvgr_legend(fp, nsets, (const char**)setnames, oenv);
 687
 688
 689     /* now make the histograms */
 690     bl=bl_head->next;
 691     /* iterate over all lambdas */
 692     while(bl!=bl_head)
 693     {
 694         sample_coll_t *sc=bl->sc->next;
 695
 696         /* iterate over all samples */
 697         while(sc!=bl->sc)
 698         {
 699             if (!first_set)
 700             {
 701                 xvgr_new_dataset(fp, 0, 0, NULL, oenv);
 702             }
 703
 704             sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &min,
 705                                   nbin_default);
 706
 707             for(i=0;i<nbin;i++)
 708             {
 709                 double xmin=i*dx + min;
 710                 double xmax=(i+1)*dx + min;
 711
 712                 fprintf(fp, "%g %d\n%g %d\n", xmin, hist[i], xmax, hist[i]);
 713             }
 714
 715             first_set=FALSE;
 716             sc=sc->next;
 717         }
 718
 719         bl=bl->next;
 720     }
 721
 722     if(hist)
 723         sfree(hist);
 724
 725     xvgrclose(fp);
 726 }
 727
 728 /* create a collection (array) of barres_t object given a ordered linked list
 729    of barlamda_t sample collections */
 730 static barres_t *barres_list_create(lambda_t *bl_head, int *nres,
 731                                     gmx_bool use_dhdl)
 732 {
 733     lambda_t *bl;
 734     int nlambda=0;
 735     barres_t *res;
 736     int i;
 737     gmx_bool dhdl=FALSE;
 738     gmx_bool first=TRUE;
 739
 740     /* first count the lambdas */
 741     bl=bl_head->next;
 742     while(bl!=bl_head)
 743     {
 744         nlambda++;
 745         bl=bl->next;
 746     }
 747     snew(res, nlambda-1);
 748
 749     /* next put the right samples in the res */
 750     *nres=0;
 751     bl=bl_head->next->next; /* we start with the second one. */
 752     while(bl!=bl_head)
 753     {
 754         sample_coll_t *sc, *scprev;
 755         barres_t *br=&(res[*nres]);
 756         /* there is always a previous one. we search for that as a foreign
 757            lambda: */
 758         scprev=lambda_find_sample_coll(bl->prev, bl->lambda);
 759         sc=lambda_find_sample_coll(bl, bl->prev->lambda);
 760
 761         barres_init(br);
 762
 763         if (use_dhdl)
 764         {
 765             /* we use dhdl */
 766
 767             scprev=lambda_find_sample_coll(bl->prev, bl->prev->lambda);
 768             sc=lambda_find_sample_coll(bl, bl->lambda);
 769
 770             if (first)
 771             {
 772                 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");
 773                 dhdl=TRUE;
 774             }
 775             if (!dhdl)
 776             {
 777                 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");
 778             }
 779         }
 780         else if (!scprev && !sc)
 781         {
 782             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);
 783         }
 784
 785         /* normal delta H */
 786         if (!scprev)
 787         {
 788             gmx_fatal(FARGS,"Could not find a set for foreign lambda = %g\nin the files for lambda = %g",bl->lambda,bl->prev->lambda);
 789         }
 790         if (!sc)
 791         {
 792             gmx_fatal(FARGS,"Could not find a set for foreign lambda = %g\nin the files for lambda = %g",bl->prev->lambda,bl->lambda);
 793         }
 794         br->a = scprev;
 795         br->b = sc;
 796
 797         first=FALSE;
 798         (*nres)++;
 799         bl=bl->next;
 800     }
 801     return res;
 802 }
 803
 804 /* estimate the maximum discretization error */
 805 static double barres_list_max_disc_err(barres_t *res, int nres)
 806 {
 807     int i,j;
 808     double disc_err=0.;
 809     double delta_lambda;
 810
 811     for(i=0;i<nres;i++)
 812     {
 813         barres_t *br=&(res[i]);
 814
 815         delta_lambda=fabs(br->b->native_lambda-br->a->native_lambda);
 816
 817         for(j=0;j<br->a->nsamples;j++)
 818         {
 819             if (br->a->s[j]->hist)
 820             {
 821                 double Wfac=1.;
 822                 if (br->a->s[j]->derivative)
 823                     Wfac =  delta_lambda;
 824
 825                 disc_err=max(disc_err, Wfac*br->a->s[j]->hist->dx[0]);
 826             }
 827         }
 828         for(j=0;j<br->b->nsamples;j++)
 829         {
 830             if (br->b->s[j]->hist)
 831             {
 832                 double Wfac=1.;
 833                 if (br->b->s[j]->derivative)
 834                     Wfac =  delta_lambda;
 835                 disc_err=max(disc_err, Wfac*br->b->s[j]->hist->dx[0]);
 836             }
 837         }
 838     }
 839     return disc_err;
 840 }
 841
 842
 843 /* impose start and end times on a sample collection, updating sample_ranges */
 844 static void sample_coll_impose_times(sample_coll_t *sc, double begin_t,
 845                                      double end_t)
 846 {
 847     int i;
 848     for(i=0;i<sc->nsamples;i++)
 849     {
 850         samples_t *s=sc->s[i];
 851         sample_range_t *r=&(sc->r[i]);
 852         if (s->hist)
 853         {
 854             double end_time=s->hist->delta_time*s->hist->sum +
 855                             s->hist->start_time;
 856             if (s->hist->start_time < begin_t || end_time > end_t)
 857             {
 858                 r->use=FALSE;
 859             }
 860         }
 861         else
 862         {
 863             if (!s->t)
 864             {
 865                 double end_time;
 866                 if (s->start_time < begin_t)
 867                 {
 868                     r->start=(int)((begin_t - s->start_time)/s->delta_time);
 869                 }
 870                 end_time=s->delta_time*s->ndu + s->start_time;
 871                 if (end_time > end_t)
 872                 {
 873                     r->end=(int)((end_t - s->start_time)/s->delta_time);
 874                 }
 875             }
 876             else
 877             {
 878                 int j;
 879                 for(j=0;j<s->ndu;j++)
 880                 {
 881                     if (s->t[j] <begin_t)
 882                     {
 883                         r->start = j;
 884                     }
 885
 886                     if (s->t[j] >= end_t)
 887                     {
 888                         r->end = j;
 889                         break;
 890                     }
 891                 }
 892             }
 893             if (r->start > r->end)
 894             {
 895                 r->use=FALSE;
 896             }
 897         }
 898     }
 899     sample_coll_calc_ntot(sc);
 900 }
 901
 902 static void lambdas_impose_times(lambda_t *head, double begin, double end)
 903 {
 904     double first_t, last_t;
 905     double begin_t, end_t;
 906     lambda_t *lc;
 907     int j;
 908
 909     if (begin<=0 && end<0)
 910     {
 911         return;
 912     }
 913
 914     /* first determine the global start and end times */
 915     first_t = -1;
 916     last_t = -1;
 917     lc=head->next;
 918     while(lc!=head)
 919     {
 920         sample_coll_t *sc=lc->sc->next;
 921         while(sc != lc->sc)
 922         {
 923             for(j=0;j<sc->nsamples;j++)
 924             {
 925                 double start_t,end_t;
 926
 927                 start_t = sc->s[j]->start_time;
 928                 end_t =   sc->s[j]->start_time;
 929                 if (sc->s[j]->hist)
 930                 {
 931                     end_t += sc->s[j]->delta_time*sc->s[j]->hist->sum;
 932                 }
 933                 else
 934                 {
 935                     if (sc->s[j]->t)
 936                     {
 937                         end_t = sc->s[j]->t[sc->s[j]->ndu-1];
 938                     }
 939                     else
 940                     {
 941                         end_t += sc->s[j]->delta_time*sc->s[j]->ndu;
 942                     }
 943                 }
 944
 945                 if (start_t < first_t || first_t<0)
 946                 {
 947                     first_t=start_t;
 948                 }
 949                 if (end_t > last_t)
 950                 {
 951                     last_t=end_t;
 952                 }
 953             }
 954             sc=sc->next;
 955         }
 956         lc=lc->next;
 957     }
 958
 959     /* calculate the actual times */
 960     if (begin > 0 )
 961     {
 962         begin_t = begin;
 963     }
 964     else
 965     {
 966         begin_t = first_t;
 967     }
 968
 969     if (end >0 )
 970     {
 971         end_t = end;
 972     }
 973     else
 974     {
 975         end_t = last_t;
 976     }
 977     printf("\n   Samples in time interval: %.3f - %.3f\n", first_t, last_t);
 978
 979     if (begin_t > end_t)
 980     {
 981         return;
 982     }
 983     printf("Removing samples outside of: %.3f - %.3f\n", begin_t, end_t);
 984
 985     /* then impose them */
 986     lc=head->next;
 987     while(lc!=head)
 988     {
 989         sample_coll_t *sc=lc->sc->next;
 990         while(sc != lc->sc)
 991         {
 992             sample_coll_impose_times(sc, begin_t, end_t);
 993             sc=sc->next;
 994         }
 995         lc=lc->next;
 996     }
 997 }
 998
 999
1000 /* create subsample i out of ni from an existing sample_coll */
1001 static gmx_bool sample_coll_create_subsample(sample_coll_t  *sc,
1002                                              sample_coll_t *sc_orig,
1003                                              int i, int ni)
1004 {
1005     int j;
1006     int hist_start, hist_end;
1007
1008     gmx_large_int_t ntot_start;
1009     gmx_large_int_t ntot_end;
1010     gmx_large_int_t ntot_so_far;
1011
1012     *sc = *sc_orig; /* just copy all fields */
1013
1014     /* allocate proprietary memory */
1015     snew(sc->s, sc_orig->nsamples);
1016     snew(sc->r, sc_orig->nsamples);
1017
1018     /* copy the samples */
1019     for(j=0;j<sc_orig->nsamples;j++)
1020     {
1021         sc->s[j] = sc_orig->s[j];
1022         sc->r[j] = sc_orig->r[j]; /* copy the ranges too */
1023     }
1024
1025     /* now fix start and end fields */
1026     /* the casts avoid possible overflows */
1027     ntot_start=(gmx_large_int_t)(sc_orig->ntot*(double)i/(double)ni);
1028     ntot_end  =(gmx_large_int_t)(sc_orig->ntot*(double)(i+1)/(double)ni);
1029     ntot_so_far = 0;
1030     for(j=0;j<sc->nsamples;j++)
1031     {
1032         gmx_large_int_t ntot_add;
1033         gmx_large_int_t new_start, new_end;
1034
1035         if (sc->r[j].use)
1036         {
1037             if (sc->s[j]->hist)
1038             {
1039                 ntot_add = sc->s[j]->hist->sum;
1040             }
1041             else
1042             {
1043                 ntot_add = sc->r[j].end - sc->r[j].start;
1044             }
1045         }
1046         else
1047         {
1048             ntot_add = 0;
1049         }
1050
1051         if (!sc->s[j]->hist)
1052         {
1053             if (ntot_so_far < ntot_start)
1054             {
1055                 /* adjust starting point */
1056                 new_start = sc->r[j].start + (ntot_start - ntot_so_far);
1057             }
1058             else
1059             {
1060                 new_start = sc->r[j].start;
1061             }
1062             /* adjust end point */
1063             new_end = sc->r[j].start + (ntot_end - ntot_so_far);
1064             if (new_end > sc->r[j].end)
1065             {
1066                 new_end=sc->r[j].end;
1067             }
1068
1069             /* check if we're in range at all */
1070             if ( (new_end < new_start) || (new_start > sc->r[j].end) )
1071             {
1072                 new_start=0;
1073                 new_end=0;
1074             }
1075             /* and write the new range */
1076             sc->r[j].start=(int)new_start;
1077             sc->r[j].end=(int)new_end;
1078         }
1079         else
1080         {
1081             if (sc->r[j].use)
1082             {
1083                 double overlap;
1084                 double ntot_start_norm, ntot_end_norm;
1085                 /* calculate the amount of overlap of the
1086                    desired range (ntot_start -- ntot_end) onto
1087                    the histogram range (ntot_so_far -- ntot_so_far+ntot_add)*/
1088
1089                 /* first calculate normalized bounds
1090                    (where 0 is the start of the hist range, and 1 the end) */
1091                 ntot_start_norm = (ntot_start-ntot_so_far)/(double)ntot_add;
1092                 ntot_end_norm = (ntot_end-ntot_so_far)/(double)ntot_add;
1093
1094                 /* now fix the boundaries */
1095                 ntot_start_norm = min(1, max(0., ntot_start_norm));
1096                 ntot_end_norm = max(0, min(1., ntot_end_norm));
1097
1098                 /* and calculate the overlap */
1099                 overlap = ntot_end_norm - ntot_start_norm;
1100
1101                 if (overlap > 0.95) /* we allow for 5% slack */
1102                 {
1103                     sc->r[j].use = TRUE;
1104                 }
1105                 else if (overlap < 0.05)
1106                 {
1107                     sc->r[j].use = FALSE;
1108                 }
1109                 else
1110                 {
1111                     return FALSE;
1112                 }
1113             }
1114         }
1115         ntot_so_far += ntot_add;
1116     }
1117     sample_coll_calc_ntot(sc);
1118
1119     return TRUE;
1120 }
1121
1122 /* calculate minimum and maximum work values in sample collection */
1123 static void sample_coll_min_max(sample_coll_t *sc, double Wfac,
1124                                 double *Wmin, double *Wmax)
1125 {
1126     int i,j;
1127
1128     *Wmin=FLT_MAX;
1129     *Wmax=-FLT_MAX;
1130
1131     for(i=0;i<sc->nsamples;i++)
1132     {
1133         samples_t *s=sc->s[i];
1134         sample_range_t *r=&(sc->r[i]);
1135         if (r->use)
1136         {
1137             if (!s->hist)
1138             {
1139                 for(j=r->start; j<r->end; j++)
1140                 {
1141                     *Wmin = min(*Wmin,s->du[j]*Wfac);
1142                     *Wmax = max(*Wmax,s->du[j]*Wfac);
1143                 }
1144             }
1145             else
1146             {
1147                 int hd=0; /* determine the histogram direction: */
1148                 double dx;
1149                 if ( (s->hist->nhist>1) && (Wfac<0) )
1150                 {
1151                     hd=1;
1152                 }
1153                 dx=s->hist->dx[hd];
1154
1155                 for(j=s->hist->nbin[hd]-1;j>=0;j--)
1156                 {
1157                     *Wmin=min(*Wmin,Wfac*(s->hist->x0[hd])*dx);
1158                     *Wmax=max(*Wmax,Wfac*(s->hist->x0[hd])*dx);
1159                     /* look for the highest value bin with values */
1160                     if (s->hist->bin[hd][j]>0)
1161                     {
1162                         *Wmin=min(*Wmin,Wfac*(j+s->hist->x0[hd]+1)*dx);
1163                         *Wmax=max(*Wmax,Wfac*(j+s->hist->x0[hd]+1)*dx);
1164                         break;
1165                     }
1166                 }
1167             }
1168         }
1169     }
1170 }
1171
1172
1173 static double calc_bar_sum(int n,const double *W,double Wfac,double sbMmDG)
1174 {
1175     int    i;
1176     double sum;
1177
1178     sum = 0;
1179
1180     for(i=0; i<n; i++)
1181     {
1182         sum += 1./(1. + exp(Wfac*W[i] + sbMmDG));
1183     }
1184
1185     return sum;
1186 }
1187
1188 /* calculate the BAR average given a histogram
1189
1190     if type== 0, calculate the best estimate for the average,
1191     if type==-1, calculate the minimum possible value given the histogram
1192     if type== 1, calculate the maximum possible value given the histogram */
1193 static double calc_bar_sum_hist(const hist_t *hist, double Wfac, double sbMmDG,
1194                                 int type)
1195 {
1196     double sum=0.;
1197     int i;
1198     int max;
1199     /* normalization factor multiplied with bin width and
1200        number of samples (we normalize through M): */
1201     double normdx = 1.;
1202     int hd=0; /* determine the histogram direction: */
1203     double dx;
1204
1205     if ( (hist->nhist>1) && (Wfac<0) )
1206     {
1207         hd=1;
1208     }
1209     dx=hist->dx[hd];
1210     max=hist->nbin[hd]-1;
1211     if (type==1)
1212     {
1213         max=hist->nbin[hd]; /* we also add whatever was out of range */
1214     }
1215
1216     for(i=0;i<max;i++)
1217     {
1218         double x=Wfac*((i+hist->x0[hd])+0.5)*dx; /* bin middle */
1219         double pxdx=hist->bin[0][i]*normdx; /* p(x)dx */
1220
1221         sum += pxdx/(1. + exp(x + sbMmDG));
1222     }
1223
1224     return sum;
1225 }
1226
1227 static double calc_bar_lowlevel(sample_coll_t *ca, sample_coll_t *cb,
1228                                 double temp, double tol, int type)
1229 {
1230     double kT,beta,M;
1231     double DG;
1232     int i,j;
1233     double Wfac1,Wfac2,Wmin,Wmax;
1234     double DG0,DG1,DG2,dDG1;
1235     double sum1,sum2;
1236     double n1, n2; /* numbers of samples as doubles */
1237
1238     kT   = BOLTZ*temp;
1239     beta = 1/kT;
1240
1241     /* count the numbers of samples */
1242     n1 = ca->ntot;
1243     n2 = cb->ntot;
1244
1245     M = log(n1/n2);
1246
1247     if (!lambda_same(ca->native_lambda, ca->foreign_lambda))
1248     {
1249         /* this is the case when the delta U were calculated directly
1250            (i.e. we're not scaling dhdl) */
1251         Wfac1 = beta;
1252         Wfac2 = beta;
1253     }
1254     else
1255     {
1256         /* we're using dhdl, so delta_lambda needs to be a
1257            multiplication factor.  */
1258         double delta_lambda=cb->native_lambda-ca->native_lambda;
1259         Wfac1 =  beta*delta_lambda;
1260         Wfac2 = -beta*delta_lambda;
1261     }
1262
1263     if (beta < 1)
1264     {
1265         /* We print the output both in kT and kJ/mol.
1266          * Here we determine DG in kT, so when beta < 1
1267          * the precision has to be increased.
1268          */
1269         tol *= beta;
1270     }
1271
1272     /* Calculate minimum and maximum work to give an initial estimate of
1273      * delta G  as their average.
1274      */
1275     {
1276         double Wmin1, Wmin2, Wmax1, Wmax2;
1277         sample_coll_min_max(ca, Wfac1, &Wmin1, &Wmax1);
1278         sample_coll_min_max(cb, Wfac2, &Wmin2, &Wmax2);
1279
1280         Wmin=min(Wmin1, Wmin2);
1281         Wmax=max(Wmax1, Wmax2);
1282     }
1283
1284     DG0 = Wmin;
1285     DG2 = Wmax;
1286
1287     if (debug)
1288     {
1289         fprintf(debug,"DG %9.5f %9.5f\n",DG0,DG2);
1290     }
1291     /* We approximate by bisection: given our initial estimates
1292        we keep checking whether the halfway point is greater or
1293        smaller than what we get out of the BAR averages.
1294
1295        For the comparison we can use twice the tolerance. */
1296     while (DG2 - DG0 > 2*tol)
1297     {
1298         DG1 = 0.5*(DG0 + DG2);
1299
1300         /*printf("Wfac1=%g, Wfac2=%g, beta=%g, DG1=%g\n",Wfac1,Wfac2,beta,
1301           DG1);*/
1302
1303         /* calculate the BAR averages */
1304         dDG1=0.;
1305
1306         for(i=0;i<ca->nsamples;i++)
1307         {
1308             samples_t *s=ca->s[i];
1309             sample_range_t *r=&(ca->r[i]);
1310             if (r->use)
1311             {
1312                 if (s->hist)
1313                 {
1314                     dDG1 += calc_bar_sum_hist(s->hist, Wfac1, (M-DG1), type);
1315                 }
1316                 else
1317                 {
1318                     dDG1 += calc_bar_sum(r->end - r->start, s->du + r->start,
1319                                          Wfac1, (M-DG1));
1320                 }
1321             }
1322         }
1323         for(i=0;i<cb->nsamples;i++)
1324         {
1325             samples_t *s=cb->s[i];
1326             sample_range_t *r=&(cb->r[i]);
1327             if (r->use)
1328             {
1329                 if (s->hist)
1330                 {
1331                     dDG1 -= calc_bar_sum_hist(s->hist, Wfac2, -(M-DG1), type);
1332                 }
1333                 else
1334                 {
1335                     dDG1 -= calc_bar_sum(r->end - r->start, s->du + r->start,
1336                                          Wfac2, -(M-DG1));
1337                 }
1338             }
1339         }
1340
1341         if (dDG1 < 0)
1342         {
1343             DG0 = DG1;
1344         }
1345         else
1346         {
1347             DG2 = DG1;
1348         }
1349         if (debug)
1350         {
1351             fprintf(debug,"DG %9.5f %9.5f\n",DG0,DG2);
1352         }
1353     }
1354
1355     return 0.5*(DG0 + DG2);
1356 }
1357
1358 static void calc_rel_entropy(sample_coll_t *ca, sample_coll_t *cb,
1359                              double temp, double dg, double *sa, double *sb)
1360 {
1361     int i,j;
1362     double W_ab=0.;
1363     double W_ba=0.;
1364     double kT, beta;
1365     double Wfac1, Wfac2;
1366     double n1,n2;
1367
1368     kT   = BOLTZ*temp;
1369     beta = 1/kT;
1370
1371     /* count the numbers of samples */
1372     n1 = ca->ntot;
1373     n2 = cb->ntot;
1374
1375     /* to ensure the work values are the same as during the delta_G */
1376     if (!lambda_same(ca->native_lambda, ca->foreign_lambda))
1377     {
1378         /* this is the case when the delta U were calculated directly
1379            (i.e. we're not scaling dhdl) */
1380         Wfac1 = beta;
1381         Wfac2 = beta;
1382     }
1383     else
1384     {
1385         /* we're using dhdl, so delta_lambda needs to be a
1386            multiplication factor.  */
1387         double delta_lambda=cb->native_lambda-ca->native_lambda;
1388         Wfac1 =  beta*delta_lambda;
1389         Wfac2 = -beta*delta_lambda;
1390     }
1391
1392     /* first calculate the average work in both directions */
1393     for(i=0;i<ca->nsamples;i++)
1394     {
1395         samples_t *s=ca->s[i];
1396         sample_range_t *r=&(ca->r[i]);
1397         if (r->use)
1398         {
1399             if (!s->hist)
1400             {
1401                 for(j=r->start;j<r->end;j++)
1402                     W_ab += Wfac1*s->du[j];
1403             }
1404             else
1405             {
1406                 /* normalization factor multiplied with bin width and
1407                    number of samples (we normalize through M): */
1408                 double normdx = 1.;
1409                 double dx;
1410                 int hd=0; /* histogram direction */
1411                 if ( (s->hist->nhist>1) && (Wfac1<0) )
1412                 {
1413                     hd=1;
1414                 }
1415                 dx=s->hist->dx[hd];
1416
1417                 for(j=0;j<s->hist->nbin[0];j++)
1418                 {
1419                     double x=Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1420                     double pxdx=s->hist->bin[0][j]*normdx; /* p(x)dx */
1421                     W_ab += pxdx*x;
1422                 }
1423             }
1424         }
1425     }
1426     W_ab/=n1;
1427
1428     for(i=0;i<cb->nsamples;i++)
1429     {
1430         samples_t *s=cb->s[i];
1431         sample_range_t *r=&(cb->r[i]);
1432         if (r->use)
1433         {
1434             if (!s->hist)
1435             {
1436                 for(j=r->start;j<r->end;j++)
1437                     W_ba += Wfac1*s->du[j];
1438             }
1439             else
1440             {
1441                 /* normalization factor multiplied with bin width and
1442                    number of samples (we normalize through M): */
1443                 double normdx = 1.;
1444                 double dx;
1445                 int hd=0; /* histogram direction */
1446                 if ( (s->hist->nhist>1) && (Wfac2<0) )
1447                 {
1448                     hd=1;
1449                 }
1450                 dx=s->hist->dx[hd];
1451
1452                 for(j=0;j<s->hist->nbin[0];j++)
1453                 {
1454                     double x=Wfac1*((j+s->hist->x0[0])+0.5)*dx;/*bin ctr*/
1455                     double pxdx=s->hist->bin[0][j]*normdx; /* p(x)dx */
1456                     W_ba += pxdx*x;
1457                 }
1458             }
1459         }
1460     }
1461     W_ba/=n2;
1462
1463     /* then calculate the relative entropies */
1464     *sa = (W_ab - dg);
1465     *sb = (W_ba + dg);
1466 }
1467
1468 static void calc_dg_stddev(sample_coll_t *ca, sample_coll_t *cb,
1469                            double temp, double dg, double *stddev)
1470 {
1471     int i,j;
1472     double M;
1473     double sigmafact=0.;
1474     double kT, beta;
1475     double Wfac1, Wfac2;
1476     double n1, n2;
1477
1478     kT   = BOLTZ*temp;
1479     beta = 1/kT;
1480
1481     /* count the numbers of samples */
1482     n1 = ca->ntot;
1483     n2 = cb->ntot;
1484
1485     /* to ensure the work values are the same as during the delta_G */
1486     if (!lambda_same(ca->native_lambda, ca->foreign_lambda))
1487     {
1488         /* this is the case when the delta U were calculated directly
1489            (i.e. we're not scaling dhdl) */
1490         Wfac1 = beta;
1491         Wfac2 = beta;
1492     }
1493     else
1494     {
1495         /* we're using dhdl, so delta_lambda needs to be a
1496            multiplication factor.  */
1497         double delta_lambda=cb->native_lambda-ca->native_lambda;
1498         Wfac1 =  beta*delta_lambda;
1499         Wfac2 = -beta*delta_lambda;
1500     }
1501
1502     M = log(n1/n2);
1503
1504
1505     /* calculate average in both directions */
1506     for(i=0;i<ca->nsamples;i++)
1507     {
1508         samples_t *s=ca->s[i];
1509         sample_range_t *r=&(ca->r[i]);
1510         if (r->use)
1511         {
1512             if (!s->hist)
1513             {
1514                 for(j=r->start;j<r->end;j++)
1515                 {
1516                     sigmafact += 1./(2. + 2.*cosh((M + Wfac1*s->du[j] - dg)));
1517                 }
1518             }
1519             else
1520             {
1521                 /* normalization factor multiplied with bin width and
1522                    number of samples (we normalize through M): */
1523                 double normdx = 1.;
1524                 double dx;
1525                 int hd=0; /* histogram direction */
1526                 if ( (s->hist->nhist>1) && (Wfac1<0) )
1527                 {
1528                     hd=1;
1529                 }
1530                 dx=s->hist->dx[hd];
1531
1532                 for(j=0;j<s->hist->nbin[0];j++)
1533                 {
1534                     double x=Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1535                     double pxdx=s->hist->bin[0][j]*normdx; /* p(x)dx */
1536
1537                     sigmafact += pxdx/(2. + 2.*cosh((M + x - dg)));
1538                 }
1539             }
1540         }
1541     }
1542     for(i=0;i<cb->nsamples;i++)
1543     {
1544         samples_t *s=cb->s[i];
1545         sample_range_t *r=&(cb->r[i]);
1546         if (r->use)
1547         {
1548             if (!s->hist)
1549             {
1550                 for(j=r->start;j<r->end;j++)
1551                 {
1552                     sigmafact += 1./(2. + 2.*cosh((M - Wfac2*s->du[j] - dg)));
1553                 }
1554             }
1555             else
1556             {
1557                 /* normalization factor multiplied with bin width and
1558                    number of samples (we normalize through M): */
1559                 double normdx = 1.;
1560                 double dx;
1561                 int hd=0; /* histogram direction */
1562                 if ( (s->hist->nhist>1) && (Wfac2<0) )
1563                 {
1564                     hd=1;
1565                 }
1566                 dx=s->hist->dx[hd];
1567
1568                 for(j=0;j<s->hist->nbin[0];j++)
1569                 {
1570                     double x=Wfac2*((j+s->hist->x0[0])+0.5)*dx;/*bin ctr*/
1571                     double pxdx=s->hist->bin[0][j]*normdx; /* p(x)dx */
1572
1573                     sigmafact += pxdx/(2. + 2.*cosh((M - x - dg)));
1574                 }
1575             }
1576         }
1577     }
1578
1579     sigmafact /= (n1 + n2);
1580
1581
1582     /* Eq. 10 from
1583        Shirts, Bair, Hooker & Pande, Phys. Rev. Lett 91, 140601 (2003): */
1584     *stddev = sqrt(((1./sigmafact) - ( (n1+n2)/n1 + (n1+n2)/n2 )));
1585 }
1586
1587
1588
1589 static void calc_bar(barres_t *br, double tol,
1590                      int npee_min, int npee_max, gmx_bool *bEE,
1591                      double *partsum)
1592 {
1593     int npee,p;
1594     double dg_sig2,sa_sig2,sb_sig2,stddev_sig2; /* intermediate variance values
1595                                                    for calculated quantities */
1596     int nsample1, nsample2;
1597     double temp=br->a->temp;
1598     int i,j;
1599     double dg_min, dg_max;
1600     gmx_bool have_hist=FALSE;
1601
1602     br->dg=calc_bar_lowlevel(br->a, br->b, temp, tol, 0);
1603
1604     br->dg_disc_err = 0.;
1605     br->dg_histrange_err = 0.;
1606
1607     /* check if there are histograms */
1608     for(i=0;i<br->a->nsamples;i++)
1609     {
1610         if (br->a->r[i].use && br->a->s[i]->hist)
1611         {
1612             have_hist=TRUE;
1613             break;
1614         }
1615     }
1616     if (!have_hist)
1617     {
1618         for(i=0;i<br->b->nsamples;i++)
1619         {
1620             if (br->b->r[i].use && br->b->s[i]->hist)
1621             {
1622                 have_hist=TRUE;
1623                 break;
1624             }
1625         }
1626     }
1627
1628     /* calculate histogram-specific errors */
1629     if (have_hist)
1630     {
1631         dg_min=calc_bar_lowlevel(br->a, br->b, temp, tol, -1);
1632         dg_max=calc_bar_lowlevel(br->a, br->b, temp, tol, 1);
1633
1634         if (fabs(dg_max - dg_min) > GMX_REAL_EPS*10)
1635         {
1636             /* the histogram range  error is the biggest of the differences
1637                between the best estimate and the extremes */
1638             br->dg_histrange_err = fabs(dg_max - dg_min);
1639         }
1640         br->dg_disc_err = 0.;
1641         for(i=0;i<br->a->nsamples;i++)
1642         {
1643             if (br->a->s[i]->hist)
1644                 br->dg_disc_err=max(br->dg_disc_err, br->a->s[i]->hist->dx[0]);
1645         }
1646         for(i=0;i<br->b->nsamples;i++)
1647         {
1648             if (br->b->s[i]->hist)
1649                 br->dg_disc_err=max(br->dg_disc_err, br->b->s[i]->hist->dx[0]);
1650         }
1651     }
1652     calc_rel_entropy(br->a, br->b, temp, br->dg, &(br->sa), &(br->sb));
1653
1654     calc_dg_stddev(br->a, br->b, temp, br->dg, &(br->dg_stddev) );
1655
1656     dg_sig2 = 0;
1657     sa_sig2 = 0;
1658     sb_sig2 = 0;
1659     stddev_sig2 = 0;
1660
1661     *bEE=TRUE;
1662     {
1663         sample_coll_t ca, cb;
1664
1665         /* initialize the samples */
1666         sample_coll_init(&ca, br->a->native_lambda, br->a->foreign_lambda,
1667                         br->a->temp);
1668         sample_coll_init(&cb, br->b->native_lambda, br->b->foreign_lambda,
1669                         br->b->temp);
1670
1671         for(npee=npee_min; npee<=npee_max; npee++)
1672         {
1673             double dgs      = 0;
1674             double dgs2     = 0;
1675             double dsa      = 0;
1676             double dsb      = 0;
1677             double dsa2     = 0;
1678             double dsb2     = 0;
1679             double dstddev  = 0;
1680             double dstddev2 = 0;
1681
1682
1683             for(p=0; p<npee; p++)
1684             {
1685                 double dgp;
1686                 double stddevc;
1687                 double sac, sbc;
1688                 gmx_bool cac, cbc;
1689
1690                 cac=sample_coll_create_subsample(&ca, br->a, p, npee);
1691                 cbc=sample_coll_create_subsample(&cb, br->b, p, npee);
1692
1693                 if (!cac || !cbc)
1694                 {
1695                     printf("WARNING: histogram number incompatible with block number for averaging: can't do error estimate\n");
1696                     *bEE=FALSE;
1697                     if (cac)
1698                         sample_coll_destroy(&ca);
1699                     if (cbc)
1700                         sample_coll_destroy(&cb);
1701                     return;
1702                 }
1703
1704                 dgp=calc_bar_lowlevel(&ca, &cb, temp, tol, 0);
1705                 dgs  += dgp;
1706                 dgs2 += dgp*dgp;
1707
1708                 partsum[npee*(npee_max+1)+p] += dgp;
1709
1710                 calc_rel_entropy(&ca, &cb, temp, dgp, &sac, &sbc);
1711                 dsa  += sac;
1712                 dsa2 += sac*sac;
1713                 dsb  += sbc;
1714                 dsb2 += sbc*sbc;
1715                 calc_dg_stddev(&ca, &cb, temp, dgp, &stddevc );
1716
1717                 dstddev  += stddevc;
1718                 dstddev2 += stddevc*stddevc;
1719
1720                 sample_coll_destroy(&ca);
1721                 sample_coll_destroy(&cb);
1722             }
1723             dgs  /= npee;
1724             dgs2 /= npee;
1725             dg_sig2 += (dgs2-dgs*dgs)/(npee-1);
1726
1727             dsa  /= npee;
1728             dsa2 /= npee;
1729             dsb  /= npee;
1730             dsb2 /= npee;
1731             sa_sig2 += (dsa2-dsa*dsa)/(npee-1);
1732             sb_sig2 += (dsb2-dsb*dsb)/(npee-1);
1733
1734             dstddev  /= npee;
1735             dstddev2 /= npee;
1736             stddev_sig2 += (dstddev2-dstddev*dstddev)/(npee-1);
1737         }
1738         br->dg_err = sqrt(dg_sig2/(npee_max - npee_min + 1));
1739         br->sa_err = sqrt(sa_sig2/(npee_max - npee_min + 1));
1740         br->sb_err = sqrt(sb_sig2/(npee_max - npee_min + 1));
1741         br->dg_stddev_err = sqrt(stddev_sig2/(npee_max - npee_min + 1));
1742     }
1743 }
1744
1745
1746 static double bar_err(int nbmin, int nbmax, const double *partsum)
1747 {
1748     int nb,b;
1749     double svar,s,s2,dg;
1750
1751     svar = 0;
1752     for(nb=nbmin; nb<=nbmax; nb++)
1753     {
1754         s  = 0;
1755         s2 = 0;
1756         for(b=0; b<nb; b++)
1757         {
1758             dg  = partsum[nb*(nbmax+1)+b];
1759             s  += dg;
1760             s2 += dg*dg;
1761         }
1762         s  /= nb;
1763         s2 /= nb;
1764         svar += (s2 - s*s)/(nb - 1);
1765     }
1766
1767     return sqrt(svar/(nbmax + 1 - nbmin));
1768 }
1769
1770 /* deduce lambda value from legend.
1771 input:
1772     bdhdl = if true, value may be a derivative.
1773 output:
1774     bdhdl = whether the legend was for a derivative.
1775     */
1776 static double legend2lambda(char *fn,const char *legend,gmx_bool *bdhdl)
1777 {
1778     double lambda=0;
1779     const char   *ptr;
1780     gmx_bool ok=FALSE;
1781
1782     if (legend == NULL)
1783     {
1784         gmx_fatal(FARGS,"There is no legend in file '%s', can not deduce lambda",fn);
1785     }
1786     ptr = strrchr(legend,' ');
1787
1788     if (strstr(legend,"dH"))
1789     {
1790         if (! (*bdhdl))
1791         {
1792             ok=FALSE;
1793         }
1794         else
1795         {
1796             ok=TRUE;
1797         }
1798     }
1799     else
1800     {
1801         if (strchr(legend,'D') != NULL && strchr(legend,'H') != NULL)
1802         {
1803             ok=TRUE;
1804             *bdhdl=FALSE;
1805         }
1806     }
1807     if (!ptr)
1808     {
1809         ok=FALSE;
1810     }
1811
1812     if (!ok)
1813     {
1814         printf("%s\n", legend);
1815         gmx_fatal(FARGS,"There is no proper lambda legend in file '%s', can not deduce lambda",fn);
1816     }
1817     if (sscanf(ptr,"%lf",&lambda) != 1)
1818     {
1819         gmx_fatal(FARGS,"There is no proper lambda legend in file '%s', can not deduce lambda",fn);
1820     }
1821
1822     return lambda;
1823 }
1824
1825 static gmx_bool subtitle2lambda(const char *subtitle,double *lambda)
1826 {
1827     gmx_bool bFound;
1828     char *ptr;
1829
1830     bFound = FALSE;
1831
1832     /* plain text lambda string */
1833     ptr = strstr(subtitle,"lambda");
1834     if (ptr == NULL)
1835     {
1836         /* xmgrace formatted lambda string */
1837         ptr = strstr(subtitle,"\\xl\\f{}");
1838     }
1839     if (ptr == NULL)
1840     {
1841         /* xmgr formatted lambda string */
1842         ptr = strstr(subtitle,"\\8l\\4");
1843     }
1844     if (ptr != NULL)
1845     {
1846         ptr = strstr(ptr,"=");
1847     }
1848     if (ptr != NULL)
1849     {
1850         bFound = (sscanf(ptr+1,"%lf",lambda) == 1);
1851     }
1852
1853     return bFound;
1854 }
1855
1856 static double filename2lambda(char *fn)
1857 {
1858     double lambda;
1859     char   *ptr,*endptr,*digitptr;
1860     int     dirsep;
1861     ptr = fn;
1862     /* go to the end of the path string and search backward to find the last
1863        directory in the path which has to contain the value of lambda
1864      */
1865     while (ptr[1] != '\0')
1866     {
1867         ptr++;
1868     }
1869     /* searching backward to find the second directory separator */
1870     dirsep = 0;
1871     digitptr = NULL;
1872     while (ptr >= fn)
1873     {
1874         if (ptr[0] != DIR_SEPARATOR && ptr[1] == DIR_SEPARATOR)
1875         {
1876             if (dirsep == 1) break;
1877             dirsep++;
1878         }
1879         /* save the last position of a digit between the last two
1880            separators = in the last dirname */
1881         if (dirsep > 0 && isdigit(*ptr))
1882         {
1883             digitptr = ptr;
1884         }
1885         ptr--;
1886     }
1887     if (!digitptr)
1888     {
1889         gmx_fatal(FARGS,"While trying to read the lambda value from the file path:"
1890                     " last directory in the path '%s' does not contain a number",fn);
1891     }
1892     if (digitptr[-1] == '-')
1893     {
1894         digitptr--;
1895     }
1896     lambda = strtod(digitptr,&endptr);
1897     if (endptr == digitptr)
1898     {
1899         gmx_fatal(FARGS,"Malformed number in file path '%s'",fn);
1900     }
1901
1902     return lambda;
1903 }
1904
1905 static void read_bar_xvg_lowlevel(char *fn, real *temp, xvg_t *ba)
1906 {
1907     int  i;
1908     char *subtitle,**legend,*ptr;
1909     int np;
1910     gmx_bool native_lambda_read=FALSE;
1911
1912     xvg_init(ba);
1913
1914     ba->filename = fn;
1915
1916     np = read_xvg_legend(fn,&ba->y,&ba->nset,&subtitle,&legend);
1917     if (!ba->y)
1918     {
1919         gmx_fatal(FARGS,"File %s contains no usable data.",fn);
1920     }
1921     ba->t  = ba->y[0];
1922
1923     snew(ba->np,ba->nset);
1924     for(i=0;i<ba->nset;i++)
1925         ba->np[i]=np;
1926
1927
1928     ba->temp = -1;
1929     if (subtitle != NULL)
1930     {
1931         ptr = strstr(subtitle,"T =");
1932         if (ptr != NULL)
1933         {
1934             ptr += 3;
1935             if (sscanf(ptr,"%lf",&ba->temp) == 1)
1936             {
1937                 if (ba->temp <= 0)
1938                 {
1939                     gmx_fatal(FARGS,"Found temperature of %g in file '%s'",
1940                               ba->temp,fn);
1941                 }
1942             }
1943         }
1944     }
1945     if (ba->temp < 0)
1946     {
1947         if (*temp <= 0)
1948         {
1949             gmx_fatal(FARGS,"Did not find a temperature in the subtitle in file '%s', use the -temp option of [TT]g_bar[tt]",fn);
1950         }
1951         ba->temp = *temp;
1952     }
1953
1954     /* Try to deduce lambda from the subtitle */
1955     if (subtitle)
1956     {
1957         if (subtitle2lambda(subtitle,&(ba->native_lambda)))
1958         {
1959             native_lambda_read=TRUE;
1960         }
1961     }
1962     snew(ba->lambda,ba->nset-1);
1963     if (legend == NULL)
1964     {
1965         /* Check if we have a single set, no legend, nset=2 means t and dH/dl */
1966         if (ba->nset == 2)
1967         {
1968             if (!native_lambda_read)
1969             {
1970                 /* Deduce lambda from the file name */
1971                 ba->native_lambda = filename2lambda(fn);
1972                 native_lambda_read=TRUE;
1973             }
1974             ba->lambda[0] = ba->native_lambda;
1975         }
1976         else
1977         {
1978             gmx_fatal(FARGS,"File %s contains multiple sets but no legends, can not determine the lambda values",fn);
1979         }
1980     }
1981     else
1982     {
1983         for(i=0; i<ba->nset-1; i++)
1984         {
1985             gmx_bool is_dhdl=(i==0);
1986             /* Read lambda from the legend */
1987             ba->lambda[i] = legend2lambda(fn,legend[i], &is_dhdl);
1988
1989             if (is_dhdl && !native_lambda_read)
1990             {
1991                 ba->native_lambda = ba->lambda[i];
1992                 native_lambda_read=TRUE;
1993             }
1994         }
1995     }
1996
1997     if (!native_lambda_read)
1998     {
1999         gmx_fatal(FARGS,"File %s contains multiple sets but no indication of the native lambda",fn);
2000     }
2001
2002     /* Reorder the data */
2003     for(i=1; i<ba->nset; i++)
2004     {
2005         ba->y[i-1] = ba->y[i];
2006     }
2007     if (legend != NULL)
2008     {
2009         for(i=0; i<ba->nset-1; i++)
2010         {
2011             sfree(legend[i]);
2012         }
2013         sfree(legend);
2014     }
2015     ba->nset--;
2016 }
2017
2018 static void read_bar_xvg(char *fn, real *temp, lambda_t *lambda_head)
2019 {
2020     xvg_t *barsim;
2021     samples_t *s;
2022     int i;
2023     double *lambda;
2024
2025     snew(barsim,1);
2026
2027     read_bar_xvg_lowlevel(fn, temp, barsim);
2028
2029     if (barsim->nset <1 )
2030     {
2031         gmx_fatal(FARGS,"File '%s' contains fewer than two columns", fn);
2032     }
2033
2034     if ( !gmx_within_tol(*temp,barsim->temp,GMX_FLOAT_EPS) && (*temp > 0) )
2035     {
2036         gmx_fatal(FARGS,"Temperature in file %s different from earlier files or setting\n", fn);
2037     }
2038     *temp=barsim->temp;
2039
2040     /* now create a series of samples_t */
2041     snew(s, barsim->nset);
2042     for(i=0;i<barsim->nset;i++)
2043     {
2044         samples_init(s+i, barsim->native_lambda, barsim->lambda[i],
2045                      barsim->temp, lambda_same(barsim->native_lambda,
2046                                                barsim->lambda[i]),
2047                      fn);
2048         s[i].du=barsim->y[i];
2049         s[i].ndu=barsim->np[i];
2050         s[i].t=barsim->t;
2051
2052         lambda_list_insert_sample(lambda_head, s+i);
2053     }
2054     printf("%s: %.1f - %.1f; lambda = %.3f\n    foreign lambdas:",
2055            fn, s[0].t[0], s[0].t[s[0].ndu-1], s[0].native_lambda);
2056     for(i=0;i<barsim->nset;i++)
2057     {
2058         printf(" %.3f (%d pts)", s[i].foreign_lambda, s[i].ndu);
2059     }
2060     printf("\n\n");
2061 }
2062
2063 static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
2064                                  double start_time, double delta_time,
2065                                  double native_lambda, double temp,
2066                                  double *last_t, const char *filename)
2067 {
2068     int j;
2069     gmx_bool allocated;
2070     double foreign_lambda;
2071     int derivative;
2072     samples_t *s; /* convenience pointer */
2073     int startj;
2074
2075     /* check the block types etc. */
2076     if ( (blk->nsub < 3) ||
2077          (blk->sub[0].type != xdr_datatype_int) ||
2078          (blk->sub[1].type != xdr_datatype_double) ||
2079          (
2080           (blk->sub[2].type != xdr_datatype_float) &&
2081           (blk->sub[2].type != xdr_datatype_double)
2082          ) ||
2083          (blk->sub[0].nr < 1) ||
2084          (blk->sub[1].nr < 1) )
2085     {
2086         gmx_fatal(FARGS,
2087                   "Unexpected/corrupted block data in file %s around time %g.",
2088                   filename, start_time);
2089     }
2090
2091     derivative = blk->sub[0].ival[0];
2092     foreign_lambda = blk->sub[1].dval[0];
2093
2094     if (! *smp)
2095     {
2096         /* initialize the samples structure if it's empty. */
2097         snew(*smp, 1);
2098         samples_init(*smp, native_lambda, foreign_lambda, temp,
2099                      derivative!=0, filename);
2100         (*smp)->start_time=start_time;
2101         (*smp)->delta_time=delta_time;
2102     }
2103
2104     /* set convenience pointer */
2105     s=*smp;
2106
2107     /* now double check */
2108     if ( ! lambda_same(s->foreign_lambda, foreign_lambda) ||
2109          (  (derivative!=0) != (s->derivative!=0) ) )
2110     {
2111         fprintf(stderr, "Got foreign lambda=%g, expected: %g\n",
2112                 foreign_lambda, s->foreign_lambda);
2113         fprintf(stderr, "Got derivative=%d, expected: %d\n",
2114                 derivative, s->derivative);
2115         gmx_fatal(FARGS, "Corrupted data in file %s around t=%g.",
2116                   filename, start_time);
2117     }
2118
2119     /* make room for the data */
2120     if (s->ndu_alloc < (size_t)(s->ndu + blk->sub[2].nr) )
2121     {
2122         s->ndu_alloc += (s->ndu_alloc < (size_t)blk->sub[2].nr) ?
2123                             blk->sub[2].nr*2 : s->ndu_alloc;
2124         srenew(s->du_alloc, s->ndu_alloc);
2125         s->du=s->du_alloc;
2126     }
2127     startj = s->ndu;
2128     s->ndu += blk->sub[2].nr;
2129     s->ntot += blk->sub[2].nr;
2130     *ndu = blk->sub[2].nr;
2131
2132     /* and copy the data*/
2133     for(j=0;j<blk->sub[2].nr;j++)
2134     {
2135         if (blk->sub[2].type == xdr_datatype_float)
2136         {
2137             s->du[startj+j] = blk->sub[2].fval[j];
2138         }
2139         else
2140         {
2141             s->du[startj+j] = blk->sub[2].dval[j];
2142         }
2143     }
2144     if (start_time + blk->sub[2].nr*delta_time > *last_t)
2145     {
2146         *last_t = start_time + blk->sub[2].nr*delta_time;
2147     }
2148 }
2149
2150 static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
2151                                       double start_time, double delta_time,
2152                                       double native_lambda, double temp,
2153                                       double *last_t, const char *filename)
2154 {
2155     int i,j;
2156     samples_t *s;
2157     int nhist;
2158     double foreign_lambda;
2159     int derivative;
2160     int nbins[2];
2161
2162     /* check the block types etc. */
2163     if ( (blk->nsub < 2) ||
2164          (blk->sub[0].type != xdr_datatype_double) ||
2165          (blk->sub[1].type != xdr_datatype_large_int) ||
2166          (blk->sub[0].nr < 2)  ||
2167          (blk->sub[1].nr < 2) )
2168     {
2169         gmx_fatal(FARGS,
2170                   "Unexpected/corrupted block data in file %s around time %g",
2171                   filename, start_time);
2172     }
2173
2174     nhist=blk->nsub-2;
2175     if (nhist == 0)
2176     {
2177         return NULL;
2178     }
2179     if (nhist > 2)
2180     {
2181         gmx_fatal(FARGS,
2182                   "Unexpected/corrupted block data in file %s around time %g",
2183                   filename, start_time);
2184     }
2185
2186     snew(s, 1);
2187     *nsamples=1;
2188
2189     foreign_lambda=blk->sub[0].dval[0];
2190     derivative=(int)(blk->sub[1].lval[1]);
2191     if (derivative)
2192         foreign_lambda=native_lambda;
2193
2194     samples_init(s, native_lambda, foreign_lambda, temp,
2195                  derivative!=0, filename);
2196     snew(s->hist, 1);
2197
2198     for(i=0;i<nhist;i++)
2199     {
2200         nbins[i] = blk->sub[i+2].nr;
2201     }
2202
2203     hist_init(s->hist, nhist, nbins);
2204
2205     for(i=0;i<nhist;i++)
2206     {
2207         s->hist->x0[i]=blk->sub[1].lval[2+i];
2208         s->hist->dx[i] = blk->sub[0].dval[1];
2209         if (i==1)
2210             s->hist->dx[i] = - s->hist->dx[i];
2211     }
2212
2213     s->hist->start_time = start_time;
2214     s->hist->delta_time = delta_time;
2215     s->start_time = start_time;
2216     s->delta_time = delta_time;
2217
2218     for(i=0;i<nhist;i++)
2219     {
2220         int nbin;
2221         gmx_large_int_t sum=0;
2222
2223         for(j=0;j<s->hist->nbin[i];j++)
2224         {
2225             int binv=(int)(blk->sub[i+2].ival[j]);
2226
2227             s->hist->bin[i][j] = binv;
2228             sum += binv;
2229
2230         }
2231         if (i==0)
2232         {
2233             s->ntot = sum;
2234             s->hist->sum = sum;
2235         }
2236         else
2237         {
2238             if (s->ntot != sum)
2239             {
2240                 gmx_fatal(FARGS, "Histogram counts don't match in %s",
2241                           filename);
2242             }
2243         }
2244     }
2245
2246     if (start_time + s->hist->sum*delta_time > *last_t)
2247     {
2248         *last_t = start_time + s->hist->sum*delta_time;
2249     }
2250     return s;
2251 }
2252
2253
2254 static void read_barsim_edr(char *fn, real *temp, lambda_t *lambda_head)
2255 {
2256     int i;
2257     ener_file_t fp;
2258     t_enxframe *fr;
2259     int nre;
2260     gmx_enxnm_t *enm=NULL;
2261     double first_t=-1;
2262     double last_t=-1;
2263     samples_t **samples_rawdh=NULL; /* contains samples for raw delta_h  */
2264     int *nhists=NULL;       /* array to keep count & print at end */
2265     int *npts=NULL;         /* array to keep count & print at end */
2266     double *lambdas=NULL;   /* array to keep count & print at end */
2267     double native_lambda=-1;
2268     double end_time;        /* the end time of the last batch of samples */
2269     int nsamples=0;
2270
2271     fp = open_enx(fn,"r");
2272     do_enxnms(fp,&nre,&enm);
2273     snew(fr, 1);
2274
2275     while(do_enx(fp, fr))
2276     {
2277         /* count the data blocks */
2278         int nblocks_raw=0;
2279         int nblocks_hist=0;
2280         int nlam=0;
2281         int k;
2282         /* DHCOLL block information: */
2283         double start_time=0, delta_time=0, start_lambda=0, delta_lambda=0;
2284         double rtemp=0;
2285
2286         /* count the blocks and handle collection information: */
2287         for(i=0;i<fr->nblock;i++)
2288         {
2289             if (fr->block[i].id == enxDHHIST )
2290                 nblocks_hist++;
2291             if ( fr->block[i].id == enxDH )
2292                 nblocks_raw++;
2293             if (fr->block[i].id == enxDHCOLL )
2294             {
2295                 nlam++;
2296                 if ( (fr->block[i].nsub < 1) ||
2297                      (fr->block[i].sub[0].type != xdr_datatype_double) ||
2298                      (fr->block[i].sub[0].nr < 5))
2299                 {
2300                     gmx_fatal(FARGS, "Unexpected block data in file %s", fn);
2301                 }
2302
2303                 /* read the data from the DHCOLL block */
2304                 rtemp =        fr->block[i].sub[0].dval[0];
2305                 start_time =   fr->block[i].sub[0].dval[1];
2306                 delta_time =   fr->block[i].sub[0].dval[2];
2307                 start_lambda = fr->block[i].sub[0].dval[3];
2308                 delta_lambda = fr->block[i].sub[0].dval[4];
2309
2310                 if (delta_lambda>0)
2311                 {
2312                     gmx_fatal(FARGS, "Lambda values not constant in %s: can't apply BAR method", fn);
2313                 }
2314                 if ( ( *temp != rtemp) && (*temp > 0) )
2315                 {
2316                     gmx_fatal(FARGS,"Temperature in file %s different from earlier files or setting\n", fn);
2317                 }
2318                 *temp=rtemp;
2319
2320                 if (first_t < 0)
2321                     first_t=start_time;
2322             }
2323         }
2324
2325         if (nlam != 1)
2326         {
2327             gmx_fatal(FARGS, "Did not find a delta h information in file %s" , fn);
2328         }
2329         if (nblocks_raw > 0 && nblocks_hist > 0 )
2330         {
2331             gmx_fatal(FARGS, "Can't handle both raw delta U data and histograms in the same file %s", fn);
2332         }
2333
2334         if (nsamples > 0)
2335         {
2336             /* check the native lambda */
2337             if (!lambda_same(start_lambda, native_lambda) )
2338             {
2339                 gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %g, and becomes %g at time %g",
2340                           fn, native_lambda, start_lambda, start_time);
2341             }
2342             /* check the number of samples against the previous number */
2343             if ( ((nblocks_raw+nblocks_hist)!=nsamples) || (nlam!=1 ) )
2344             {
2345                 gmx_fatal(FARGS, "Unexpected block count in %s: was %d, now %d\n",
2346                           fn, nsamples+1, nblocks_raw+nblocks_hist+nlam);
2347             }
2348             /* check whether last iterations's end time matches with
2349                the currrent start time */
2350             if ( (fabs(last_t - start_time) > 2*delta_time)  && last_t>=0)
2351             {
2352                 /* it didn't. We need to store our samples and reallocate */
2353                 for(i=0;i<nsamples;i++)
2354                 {
2355                     if (samples_rawdh[i])
2356                     {
2357                         /* insert it into the existing list */
2358                         lambda_list_insert_sample(lambda_head,
2359                                                   samples_rawdh[i]);
2360                         /* and make sure we'll allocate a new one this time
2361                            around */
2362                         samples_rawdh[i]=NULL;
2363                     }
2364                 }
2365             }
2366         }
2367         else
2368         {
2369             /* this is the first round; allocate the associated data
2370                structures */
2371             native_lambda=start_lambda;
2372             nsamples=nblocks_raw+nblocks_hist;
2373             snew(nhists, nsamples);
2374             snew(npts, nsamples);
2375             snew(lambdas, nsamples);
2376             snew(samples_rawdh, nsamples);
2377             for(i=0;i<nsamples;i++)
2378             {
2379                 nhists[i]=0;
2380                 npts[i]=0;
2381                 lambdas[i]=-1;
2382                 samples_rawdh[i]=NULL; /* init to NULL so we know which
2383                                           ones contain values */
2384             }
2385         }
2386
2387         /* and read them */
2388         k=0; /* counter for the lambdas, etc. arrays */
2389         for(i=0;i<fr->nblock;i++)
2390         {
2391             if (fr->block[i].id == enxDH)
2392             {
2393                 int ndu;
2394                 read_edr_rawdh_block(&(samples_rawdh[k]),
2395                                      &ndu,
2396                                      &(fr->block[i]),
2397                                      start_time, delta_time,
2398                                      start_lambda, rtemp,
2399                                      &last_t, fn);
2400                 npts[k] += ndu;
2401                 if (samples_rawdh[k])
2402                 {
2403                     lambdas[k]=samples_rawdh[k]->foreign_lambda;
2404                 }
2405                 k++;
2406             }
2407             else if (fr->block[i].id == enxDHHIST)
2408             {
2409                 int j;
2410                 int nb=0;
2411                 samples_t *s; /* this is where the data will go */
2412                 s=read_edr_hist_block(&nb, &(fr->block[i]),
2413                                       start_time, delta_time,
2414                                       start_lambda, rtemp,
2415                                       &last_t, fn);
2416                 nhists[k] += nb;
2417                 if (nb>0)
2418                 {
2419                     lambdas[k]= s->foreign_lambda;
2420                 }
2421                 k++;
2422                 /* and insert the new sample immediately */
2423                 for(j=0;j<nb;j++)
2424                 {
2425                     lambda_list_insert_sample(lambda_head, s+j);
2426                 }
2427             }
2428         }
2429     }
2430     /* Now store all our extant sample collections */
2431     for(i=0;i<nsamples;i++)
2432     {
2433         if (samples_rawdh[i])
2434         {
2435             /* insert it into the existing list */
2436             lambda_list_insert_sample(lambda_head, samples_rawdh[i]);
2437         }
2438     }
2439
2440
2441     fprintf(stderr, "\n");
2442     printf("%s: %.1f - %.1f; lambda = %.3f\n    foreign lambdas:",
2443            fn, first_t, last_t, native_lambda);
2444     for(i=0;i<nsamples;i++)
2445     {
2446         if (nhists[i] > 0)
2447         {
2448             printf(" %.3f (%d hists)", lambdas[i], nhists[i]);
2449         }
2450         else
2451         {
2452             printf(" %.3f (%d pts)", lambdas[i], npts[i]);
2453         }
2454     }
2455     printf("\n\n");
2456     sfree(npts);
2457     sfree(nhists);
2458     sfree(lambdas);
2459 }
2460
2461
2462 int gmx_bar(int argc,char *argv[])
2463 {
2464     static const char *desc[] = {
2465         "[TT]g_bar[tt] calculates free energy difference estimates through ",
2466         "Bennett's acceptance ratio method (BAR). It also automatically",
2467         "adds series of individual free energies obtained with BAR into",
2468         "a combined free energy estimate.[PAR]",
2469
2470         "Every individual BAR free energy difference relies on two ",
2471         "simulations at different states: say state A and state B, as",
2472         "controlled by a parameter, [GRK]lambda[grk] (see the [TT].mdp[tt] parameter",
2473         "[TT]init_lambda[tt]). The BAR method calculates a ratio of weighted",
2474         "average of the Hamiltonian difference of state B given state A and",
2475         "vice versa.",
2476         "The energy differences to the other state must be calculated",
2477         "explicitly during the simulation. This can be done with",
2478         "the [TT].mdp[tt] option [TT]foreign_lambda[tt].[PAR]",
2479
2480         "Input option [TT]-f[tt] expects multiple [TT]dhdl.xvg[tt] files. ",
2481         "Two types of input files are supported:[BR]",
2482         "[TT]*[tt]  Files with more than one [IT]y[it]-value. ",
2483         "The files should have columns ",
2484         "with dH/d[GRK]lambda[grk] and [GRK]Delta[grk][GRK]lambda[grk]. ",
2485         "The [GRK]lambda[grk] values are inferred ",
2486         "from the legends: [GRK]lambda[grk] of the simulation from the legend of ",
2487         "dH/d[GRK]lambda[grk] and the foreign [GRK]lambda[grk] values from the ",
2488         "legends of Delta H",
2489         "[BR]",
2490         "[TT]*[tt]  Files with only one [IT]y[it]-value. Using the",
2491         "[TT]-extp[tt] option for these files, it is assumed",
2492         "that the [IT]y[it]-value is dH/d[GRK]lambda[grk] and that the ",
2493         "Hamiltonian depends linearly on [GRK]lambda[grk]. ",
2494         "The [GRK]lambda[grk] value of the simulation is inferred from the ",
2495         "subtitle (if present), otherwise from a number in the subdirectory ",
2496         "in the file name.[PAR]",
2497
2498         "The [GRK]lambda[grk] of the simulation is parsed from ",
2499         "[TT]dhdl.xvg[tt] file's legend containing the string 'dH', the ",
2500         "foreign [GRK]lambda[grk] values from the legend containing the ",
2501         "capitalized letters 'D' and 'H'. The temperature is parsed from ",
2502         "the legend line containing 'T ='.[PAR]",
2503
2504         "The input option [TT]-g[tt] expects multiple [TT].edr[tt] files. ",
2505         "These can contain either lists of energy differences (see the ",
2506         "[TT].mdp[tt] option [TT]separate_dhdl_file[tt]), or a series of ",
2507         "histograms (see the [TT].mdp[tt] options [TT]dh_hist_size[tt] and ",
2508         "[TT]dh_hist_spacing[tt]).", "The temperature and [GRK]lambda[grk] ",
2509         "values are automatically deduced from the [TT]ener.edr[tt] file.[PAR]",
2510
2511         "In addition to the [TT].mdp[tt] option [TT]foreign_lambda[tt], ",
2512         "the energy difference can also be extrapolated from the ",
2513         "dH/d[GRK]lambda[grk] values. This is done with the[TT]-extp[tt]",
2514         "option, which assumes that the system's Hamiltonian depends linearly",
2515         "on [GRK]lambda[grk], which is not normally the case.[PAR]",
2516
2517         "The free energy estimates are determined using BAR with bisection, ",
2518         "with the precision of the output set with [TT]-prec[tt]. ",
2519         "An error estimate taking into account time correlations ",
2520         "is made by splitting the data into blocks and determining ",
2521         "the free energy differences over those blocks and assuming ",
2522         "the blocks are independent. ",
2523         "The final error estimate is determined from the average variance ",
2524         "over 5 blocks. A range of block numbers for error estimation can ",
2525         "be provided with the options [TT]-nbmin[tt] and [TT]-nbmax[tt].[PAR]",
2526
2527         "[TT]g_bar[tt] tries to aggregate samples with the same 'native' and ",
2528         "'foreign' [GRK]lambda[grk] values, but always assumes independent ",
2529         "samples. [BB]Note[bb] that when aggregating energy ",
2530         "differences/derivatives with different sampling intervals, this is ",
2531         "almost certainly not correct. Usually subsequent energies are ",
2532         "correlated and different time intervals mean different degrees ",
2533         "of correlation between samples.[PAR]",
2534
2535         "The results are split in two parts: the last part contains the final ",
2536         "results in kJ/mol, together with the error estimate for each part ",
2537         "and the total. The first part contains detailed free energy ",
2538         "difference estimates and phase space overlap measures in units of ",
2539         "kT (together with their computed error estimate). The printed ",
2540         "values are:[BR]",
2541         "[TT]*[tt]  lam_A: the [GRK]lambda[grk] values for point A.[BR]",
2542         "[TT]*[tt]  lam_B: the [GRK]lambda[grk] values for point B.[BR]",
2543         "[TT]*[tt]     DG: the free energy estimate.[BR]",
2544         "[TT]*[tt]    s_A: an estimate of the relative entropy of B in A.[BR]",
2545         "[TT]*[tt]    s_A: an estimate of the relative entropy of A in B.[BR]",
2546         "[TT]*[tt]  stdev: an estimate expected per-sample standard deviation.[PAR]",
2547
2548         "The relative entropy of both states in each other's ensemble can be ",
2549         "interpreted as a measure of phase space overlap: ",
2550         "the relative entropy s_A of the work samples of lambda_B in the ",
2551         "ensemble of lambda_A (and vice versa for s_B), is a ",
2552         "measure of the 'distance' between Boltzmann distributions of ",
2553         "the two states, that goes to zero for identical distributions. See ",
2554         "Wu & Kofke, J. Chem. Phys. 123 084109 (2005) for more information.",
2555         "[PAR]",
2556         "The estimate of the expected per-sample standard deviation, as given ",
2557         "in Bennett's original BAR paper: Bennett, J. Comp. Phys. 22, p 245 (1976).",
2558         "Eq. 10 therein gives an estimate of the quality of sampling (not directly",
2559         "of the actual statistical error, because it assumes independent samples).[PAR]",
2560
2561         "To get a visual estimate of the phase space overlap, use the ",
2562         "[TT]-oh[tt] option to write series of histograms, together with the ",
2563         "[TT]-nbin[tt] option.[PAR]"
2564     };
2565     static real begin=0,end=-1,temp=-1;
2566     int nd=2,nbmin=5,nbmax=5;
2567     int nbin=100;
2568     gmx_bool use_dhdl=FALSE;
2569     gmx_bool calc_s,calc_v;
2570     t_pargs pa[] = {
2571         { "-b",    FALSE, etREAL, {&begin},  "Begin time for BAR" },
2572         { "-e",    FALSE, etREAL, {&end},    "End time for BAR" },
2573         { "-temp", FALSE, etREAL, {&temp},   "Temperature (K)" },
2574         { "-prec", FALSE, etINT,  {&nd},     "The number of digits after the decimal point" },
2575         { "-nbmin",  FALSE, etINT,  {&nbmin}, "Minimum number of blocks for error estimation" },
2576         { "-nbmax",  FALSE, etINT,  {&nbmax}, "Maximum number of blocks for error estimation" },
2577         { "-nbin",  FALSE, etINT, {&nbin}, "Number of bins for histogram output"},
2578         { "-extp",  FALSE, etBOOL, {&use_dhdl}, "Whether to linearly extrapolate dH/dl values to use as energies"}
2579     };
2580
2581     t_filenm   fnm[] = {
2582         { efXVG, "-f",  "dhdl",   ffOPTRDMULT },
2583         { efEDR, "-g",  "ener",   ffOPTRDMULT },
2584         { efXVG, "-o",  "bar",    ffOPTWR },
2585         { efXVG, "-oi", "barint", ffOPTWR },
2586         { efXVG, "-oh", "histogram", ffOPTWR }
2587     };
2588 #define NFILE asize(fnm)
2589
2590     int      f,i,j;
2591     int      nf=0; /* file counter */
2592     int      nbs;
2593     int      nfile_tot; /* total number of input files */
2594     int      nxvgfile=0;
2595     int      nedrfile=0;
2596     char     **fxvgnms;
2597     char     **fedrnms;
2598     lambda_t *lb;    /* the pre-processed lambda data (linked list head) */
2599     lambda_t lambda_head; /* the head element */
2600     barres_t *results;  /* the results */
2601     int    nresults;  /* number of results in results array */
2602
2603     double   *partsum;
2604     double   prec,dg_tot,dg,sig, dg_tot_max, dg_tot_min;
2605     FILE     *fpb,*fpi;
2606     char     lamformat[20];
2607     char     dgformat[20],xvg2format[STRLEN],xvg3format[STRLEN],buf[STRLEN];
2608     char     ktformat[STRLEN], sktformat[STRLEN];
2609     char     kteformat[STRLEN], skteformat[STRLEN];
2610     output_env_t oenv;
2611     double   kT, beta;
2612     gmx_bool     result_OK=TRUE,bEE=TRUE;
2613
2614     gmx_bool     disc_err=FALSE;
2615     double   sum_disc_err=0.; /* discretization error */
2616     gmx_bool     histrange_err=FALSE;
2617     double   sum_histrange_err=0.; /* histogram range error */
2618     double   stat_err=0.; /* statistical error */
2619
2620     parse_common_args(&argc,argv,
2621                       PCA_CAN_VIEW,
2622                       NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
2623
2624     if (opt2bSet("-f", NFILE, fnm))
2625     {
2626         nxvgfile = opt2fns(&fxvgnms,"-f",NFILE,fnm);
2627     }
2628     if (opt2bSet("-g", NFILE, fnm))
2629     {
2630         nedrfile = opt2fns(&fedrnms,"-g",NFILE,fnm);
2631     }
2632
2633     /* make linked list */
2634     lb=&lambda_head;
2635     lambda_init(lb, 0, 0);
2636     lb->next=lb;
2637     lb->prev=lb;
2638
2639
2640     nfile_tot = nxvgfile + nedrfile;
2641
2642     if (nfile_tot == 0)
2643     {
2644         gmx_fatal(FARGS,"No input files!");
2645     }
2646
2647     if (nd < 0)
2648     {
2649         gmx_fatal(FARGS,"Can not have negative number of digits");
2650     }
2651     prec = pow(10,-nd);
2652
2653     snew(partsum,(nbmax+1)*(nbmax+1));
2654     nf = 0;
2655
2656     /* read in all files. First xvg files */
2657     for(f=0; f<nxvgfile; f++)
2658     {
2659         read_bar_xvg(fxvgnms[f],&temp,lb);
2660         nf++;
2661     }
2662     /* then .edr files */
2663     for(f=0; f<nedrfile; f++)
2664     {
2665         read_barsim_edr(fedrnms[f],&temp,lb);;
2666         nf++;
2667     }
2668
2669     /* fix the times to allow for equilibration */
2670     lambdas_impose_times(lb, begin, end);
2671
2672     if (opt2bSet("-oh",NFILE,fnm))
2673     {
2674         lambdas_histogram(lb, opt2fn("-oh",NFILE,fnm), nbin, oenv);
2675     }
2676
2677     /* assemble the output structures from the lambdas */
2678     results=barres_list_create(lb, &nresults, use_dhdl);
2679
2680     sum_disc_err=barres_list_max_disc_err(results, nresults);
2681
2682     if (nresults == 0)
2683     {
2684         printf("\nNo results to calculate.\n");
2685         return 0;
2686     }
2687
2688     if (sum_disc_err > prec)
2689     {
2690         prec=sum_disc_err;
2691         nd = ceil(-log10(prec));
2692         printf("WARNING: setting the precision to %g because that is the minimum\n         reasonable number, given the expected discretization error.\n", prec);
2693     }
2694
2695
2696     sprintf(lamformat,"%%6.3f");
2697     sprintf( dgformat,"%%%d.%df",3+nd,nd);
2698     /* the format strings of the results in kT */
2699     sprintf( ktformat,"%%%d.%df",5+nd,nd);
2700     sprintf( sktformat,"%%%ds",6+nd);
2701     /* the format strings of the errors in kT */
2702     sprintf( kteformat,"%%%d.%df",3+nd,nd);
2703     sprintf( skteformat,"%%%ds",4+nd);
2704     sprintf(xvg2format,"%s %s\n","%g",dgformat);
2705     sprintf(xvg3format,"%s %s %s\n","%g",dgformat,dgformat);
2706
2707
2708
2709     fpb = NULL;
2710     if (opt2bSet("-o",NFILE,fnm))
2711     {
2712         sprintf(buf,"%s (%s)","\\DeltaG","kT");
2713         fpb = xvgropen_type(opt2fn("-o",NFILE,fnm),"Free energy differences",
2714                             "\\lambda",buf,exvggtXYDY,oenv);
2715     }
2716
2717     fpi = NULL;
2718     if (opt2bSet("-oi",NFILE,fnm))
2719     {
2720         sprintf(buf,"%s (%s)","\\DeltaG","kT");
2721         fpi = xvgropen(opt2fn("-oi",NFILE,fnm),"Free energy integral",
2722                       "\\lambda",buf,oenv);
2723     }
2724
2725
2726
2727     if (nbmin > nbmax)
2728         nbmin=nbmax;
2729
2730     /* first calculate results */
2731     bEE = TRUE;
2732     disc_err = FALSE;
2733     for(f=0; f<nresults; f++)
2734     {
2735         /* Determine the free energy difference with a factor of 10
2736          * more accuracy than requested for printing.
2737          */
2738         calc_bar(&(results[f]), 0.1*prec, nbmin, nbmax,
2739                  &bEE, partsum);
2740
2741         if (results[f].dg_disc_err > prec/10.)
2742             disc_err=TRUE;
2743         if (results[f].dg_histrange_err > prec/10.)
2744             histrange_err=TRUE;
2745     }
2746
2747     /* print results in kT */
2748     kT   = BOLTZ*temp;
2749     beta = 1/kT;
2750
2751     printf("\nTemperature: %g K\n", temp);
2752
2753     printf("\nDetailed results in kT (see help for explanation):\n\n");
2754     printf("%6s ", " lam_A");
2755     printf("%6s ", " lam_B");
2756     printf(sktformat,  "DG ");
2757     if (bEE)
2758         printf(skteformat, "+/- ");
2759     if (disc_err)
2760         printf(skteformat, "disc ");
2761     if (histrange_err)
2762         printf(skteformat, "range ");
2763     printf(sktformat,  "s_A ");
2764     if (bEE)
2765         printf(skteformat, "+/- " );
2766     printf(sktformat,  "s_B ");
2767     if (bEE)
2768         printf(skteformat, "+/- " );
2769     printf(sktformat,  "stdev ");
2770     if (bEE)
2771         printf(skteformat, "+/- ");
2772     printf("\n");
2773     for(f=0; f<nresults; f++)
2774     {
2775         printf(lamformat, results[f].a->native_lambda);
2776         printf(" ");
2777         printf(lamformat, results[f].b->native_lambda);
2778         printf(" ");
2779         printf(ktformat,  results[f].dg);
2780         printf(" ");
2781         if (bEE)
2782         {
2783             printf(kteformat, results[f].dg_err);
2784             printf(" ");
2785         }
2786         if (disc_err)
2787         {
2788             printf(kteformat, results[f].dg_disc_err);
2789             printf(" ");
2790         }
2791         if (histrange_err)
2792         {
2793             printf(kteformat, results[f].dg_histrange_err);
2794             printf(" ");
2795         }
2796         printf(ktformat,  results[f].sa);
2797         printf(" ");
2798         if (bEE)
2799         {
2800             printf(kteformat, results[f].sa_err);
2801             printf(" ");
2802         }
2803         printf(ktformat,  results[f].sb);
2804         printf(" ");
2805         if (bEE)
2806         {
2807             printf(kteformat, results[f].sb_err);
2808             printf(" ");
2809         }
2810         printf(ktformat,  results[f].dg_stddev);
2811         printf(" ");
2812         if (bEE)
2813         {
2814             printf(kteformat, results[f].dg_stddev_err);
2815         }
2816         printf("\n");
2817
2818         /* Check for negative relative entropy with a 95% certainty. */
2819         if (results[f].sa < -2*results[f].sa_err ||
2820             results[f].sb < -2*results[f].sb_err)
2821         {
2822             result_OK=FALSE;
2823         }
2824     }
2825
2826     if (!result_OK)
2827     {
2828         printf("\nWARNING: Some of these results violate the Second Law of "
2829                "Thermodynamics: \n"
2830                "         This is can be the result of severe undersampling, or "
2831                "(more likely)\n"
2832                "         there is something wrong with the simulations.\n");
2833     }
2834
2835
2836     /* final results in kJ/mol */
2837     printf("\n\nFinal results in kJ/mol:\n\n");
2838     dg_tot  = 0;
2839     for(f=0; f<nresults; f++)
2840     {
2841
2842         if (fpi != NULL)
2843         {
2844             fprintf(fpi, xvg2format, results[f].a->native_lambda, dg_tot);
2845         }
2846
2847
2848         if (fpb != NULL)
2849         {
2850             fprintf(fpb, xvg3format,
2851                     0.5*(results[f].a->native_lambda +
2852                          results[f].b->native_lambda),
2853                     results[f].dg,results[f].dg_err);
2854         }
2855
2856         /*printf("lambda %4.2f - %4.2f, DG ", results[f].lambda_a,
2857                                               results[f].lambda_b);*/
2858         printf("lambda ");
2859         printf(lamformat, results[f].a->native_lambda);
2860         printf(" - ");
2861         printf(lamformat, results[f].b->native_lambda);
2862         printf(",   DG ");
2863
2864         printf(dgformat,results[f].dg*kT);
2865         if (bEE)
2866         {
2867             printf(" +/- ");
2868             printf(dgformat,results[f].dg_err*kT);
2869         }
2870         if (histrange_err)
2871         {
2872             printf(" (max. range err. = ");
2873             printf(dgformat,results[f].dg_histrange_err*kT);
2874             printf(")");
2875             sum_histrange_err += results[f].dg_histrange_err*kT;
2876         }
2877
2878         printf("\n");
2879         dg_tot += results[f].dg;
2880     }
2881     printf("\n");
2882     printf("total  ");
2883     printf(lamformat, results[0].a->native_lambda);
2884     printf(" - ");
2885     printf(lamformat, results[nresults-1].b->native_lambda);
2886     printf(",   DG ");
2887
2888     printf(dgformat,dg_tot*kT);
2889     if (bEE)
2890     {
2891         stat_err=bar_err(nbmin,nbmax,partsum)*kT;
2892         printf(" +/- ");
2893         printf(dgformat, max(max(stat_err, sum_disc_err), sum_histrange_err));
2894     }
2895     printf("\n");
2896     if (disc_err)
2897     {
2898         printf("\nmaximum discretization error = ");
2899         printf(dgformat,sum_disc_err);
2900         if (bEE && stat_err < sum_disc_err)
2901         {
2902             printf("WARNING: discretization error (%g) is larger than statistical error.\n       Decrease histogram spacing for more accurate results\n", stat_err);
2903         }
2904     }
2905     if (histrange_err)
2906     {
2907         printf("\nmaximum histogram range error = ");
2908         printf(dgformat,sum_histrange_err);
2909         if (bEE && stat_err < sum_histrange_err)
2910         {
2911             printf("WARNING: histogram range error (%g) is larger than statistical error.\n       Increase histogram range for more accurate results\n", stat_err);
2912         }
2913
2914     }
2915     printf("\n");
2916
2917
2918     if (fpi != NULL)
2919     {
2920         fprintf(fpi, xvg2format,
2921                 results[nresults-1].b->native_lambda, dg_tot);
2922         ffclose(fpi);
2923     }
2924
2925     do_view(oenv,opt2fn_null("-o",NFILE,fnm),"-xydy");
2926     do_view(oenv,opt2fn_null("-oi",NFILE,fnm),"-xydy");
2927
2928     thanx(stderr);
2929
2930     return 0;
2931 }