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