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