src/gromacs/gmxana/gmx_bar.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  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 "gromacs/utility/smalloc.h"
  48 #include "gromacs/fileio/futil.h"
  49 #include "gromacs/commandline/pargs.h"
  50 #include "macros.h"
  51 #include "gromacs/fileio/enxio.h"
  52 #include "physics.h"
  53 #include "gmx_fatal.h"
  54 #include "xvgr.h"
  55 #include "gmx_ana.h"
  56 #include "gromacs/math/utilities.h"
  57 #include "gromacs/utility/cstringutil.h"
  58 #include "names.h"
  59 #include "mdebin.h"
  60
  61
  62 /* Structure for the names of lambda vector components */
  63 typedef struct lambda_components_t
  64 {
  65     char **names;  /* Array of strings with names for the lambda vector
  66                       components */
  67     int    N;      /* The number of components */
  68     int    Nalloc; /* The number of allocated components */
  69 } lambda_components_t;
  70
  71 /* Structure for a lambda vector or a dhdl derivative direction */
  72 typedef struct lambda_vec_t
  73 {
  74     double                    *val;   /* The lambda vector component values. Only valid if
  75                                          dhdl == -1 */
  76     int                        dhdl;  /* The coordinate index for the derivative described by this
  77                                          structure, or -1 */
  78     const lambda_components_t *lc;    /* the associated lambda_components
  79                                          structure */
  80     int                        index; /* The state number (init-lambda-state) of this lambda
  81                                          vector, if known. If not, it is set to -1 */
  82 } lambda_vec_t;
  83
  84 /* the dhdl.xvg data from a simulation */
  85 typedef struct xvg_t
  86 {
  87     const char   *filename;
  88     int           ftp;           /* file type */
  89     int           nset;          /* number of lambdas, including dhdl */
  90     int          *np;            /* number of data points (du or hists) per lambda */
  91     int           np_alloc;      /* number of points (du or hists) allocated */
  92     double        temp;          /* temperature */
  93     lambda_vec_t *lambda;        /* the lambdas (of first index for y). */
  94     double       *t;             /* the times (of second index for y) */
  95     double      **y;             /* the dU values. y[0] holds the derivative, while
  96                                     further ones contain the energy differences between
  97                                     the native lambda and the 'foreign' lambdas. */
  98     lambda_vec_t  native_lambda; /* the native lambda */
  99
 100     struct xvg_t *next, *prev;   /*location in the global linked list of xvg_ts*/
 101 } xvg_t;
 102
 103
 104 typedef struct hist_t
 105 {
 106     unsigned int   *bin[2];                 /* the (forward + reverse) histogram values */
 107     double          dx[2];                  /* the histogram spacing. The reverse
 108                                                dx is the negative of the forward dx.*/
 109     gmx_int64_t     x0[2];                  /* the (forward + reverse) histogram start
 110                                                    point(s) as int */
 111
 112     int             nbin[2];                /* the (forward+reverse) number of bins */
 113     gmx_int64_t     sum;                    /* the total number of counts. Must be
 114                                                    the same for forward + reverse.  */
 115     int             nhist;                  /* number of hist datas (forward or reverse) */
 116
 117     double          start_time, delta_time; /* start time, end time of histogram */
 118 } hist_t;
 119
 120
 121 /* an aggregate of samples for partial free energy calculation */
 122 typedef struct samples_t
 123 {
 124     lambda_vec_t *native_lambda;  /* pointer to native lambda vector */
 125     lambda_vec_t *foreign_lambda; /* pointer to foreign lambda vector */
 126     double        temp;           /* the temperature */
 127     gmx_bool      derivative;     /* whether this sample is a derivative */
 128
 129     /* The samples come either as either delta U lists: */
 130     int     ndu;                    /* the number of delta U samples */
 131     double *du;                     /* the delta u's */
 132     double *t;                      /* the times associated with those samples, or: */
 133     double  start_time, delta_time; /*start time and delta time for linear time*/
 134
 135     /* or as histograms: */
 136     hist_t *hist; /* a histogram */
 137
 138     /* allocation data: (not NULL for data 'owned' by this struct) */
 139     double         *du_alloc, *t_alloc;  /* allocated delta u arrays  */
 140     size_t          ndu_alloc, nt_alloc; /* pre-allocated sizes */
 141     hist_t         *hist_alloc;          /* allocated hist */
 142
 143     gmx_int64_t     ntot;                /* total number of samples */
 144     const char     *filename;            /* the file name this sample comes from */
 145 } samples_t;
 146
 147 /* a sample range (start to end for du-style data, or boolean
 148     for both du-style data and histograms */
 149 typedef struct sample_range_t
 150 {
 151     int        start, end; /* start and end index for du style data */
 152     gmx_bool   use;        /* whether to use this sample */
 153
 154     samples_t *s;          /* the samples this range belongs to */
 155 } sample_range_t;
 156
 157
 158 /* a collection of samples for a partial free energy calculation
 159     (i.e. the collection of samples from one native lambda to one
 160     foreign lambda) */
 161 typedef struct sample_coll_t
 162 {
 163     lambda_vec_t   *native_lambda;     /* these should be the same for all samples
 164                                           in the histogram */
 165     lambda_vec_t   *foreign_lambda;    /* collection */
 166     double          temp;              /* the temperature */
 167
 168     int             nsamples;          /* the number of samples */
 169     samples_t     **s;                 /* the samples themselves */
 170     sample_range_t *r;                 /* the sample ranges */
 171     int             nsamples_alloc;    /* number of allocated samples */
 172
 173     gmx_int64_t     ntot;              /* total number of samples in the ranges of
 174                                               this collection */
 175
 176     struct sample_coll_t *next, *prev; /* next and previous in the list */
 177 } sample_coll_t;
 178
 179 /* all the samples associated with a lambda point */
 180 typedef struct lambda_data_t
 181 {
 182     lambda_vec_t         *lambda;      /* the native lambda (at start time if dynamic) */
 183     double                temp;        /* temperature */
 184
 185     sample_coll_t        *sc;          /* the samples */
 186
 187     sample_coll_t         sc_head;     /*the pre-allocated list head for the linked list.*/
 188
 189     struct lambda_data_t *next, *prev; /* the next and prev in the list */
 190 } lambda_data_t;
 191
 192 /* Top-level data structure of simulation data */
 193 typedef struct sim_data_t
 194 {
 195     lambda_data_t      *lb;      /* a lambda data linked list */
 196     lambda_data_t       lb_head; /* The head element of the linked list */
 197
 198     lambda_components_t lc;      /* the allowed components of the lambda
 199                                     vectors */
 200 } sim_data_t;
 201
 202 /* Top-level data structure with calculated values. */
 203 typedef struct {
 204     sample_coll_t *a, *b;            /* the simulation data */
 205
 206     double         dg;               /* the free energy difference */
 207     double         dg_err;           /* the free energy difference */
 208
 209     double         dg_disc_err;      /* discretization error */
 210     double         dg_histrange_err; /* histogram range error */
 211
 212     double         sa;               /* relative entropy of b in state a */
 213     double         sa_err;           /* error in sa */
 214     double         sb;               /* relative entropy of a in state b */
 215     double         sb_err;           /* error in sb */
 216
 217     double         dg_stddev;        /* expected dg stddev per sample */
 218     double         dg_stddev_err;    /* error in dg_stddev */
 219 } barres_t;
 220
 221
 222 /* Initialize a lambda_components structure */
 223 static void lambda_components_init(lambda_components_t *lc)
 224 {
 225     lc->N      = 0;
 226     lc->Nalloc = 2;
 227     snew(lc->names, lc->Nalloc);
 228 }
 229
 230 /* Add a component to a lambda_components structure */
 231 static void lambda_components_add(lambda_components_t *lc,
 232                                   const char *name, size_t name_length)
 233 {
 234     while (lc->N + 1 > lc->Nalloc)
 235     {
 236         lc->Nalloc = (lc->Nalloc == 0) ? 2 : 2*lc->Nalloc;
 237         srenew(lc->names, lc->Nalloc);
 238     }
 239     snew(lc->names[lc->N], name_length+1);
 240     strncpy(lc->names[lc->N], name, name_length);
 241     lc->N++;
 242 }
 243
 244 /* check whether a component with index 'index' matches the given name, or
 245    is also NULL. Returns TRUE if this is the case.
 246    the string name does not need to end */
 247 static gmx_bool lambda_components_check(const lambda_components_t *lc,
 248                                         int                        index,
 249                                         const char                *name,
 250                                         size_t                     name_length)
 251 {
 252     size_t len;
 253     if (index >= lc->N)
 254     {
 255         return FALSE;
 256     }
 257     if (name == NULL && lc->names[index] == NULL)
 258     {
 259         return TRUE;
 260     }
 261     if ((name == NULL) != (lc->names[index] == NULL))
 262     {
 263         return FALSE;
 264     }
 265     len = strlen(lc->names[index]);
 266     if (len != name_length)
 267     {
 268         return FALSE;
 269     }
 270     if (strncmp(lc->names[index], name, name_length) == 0)
 271     {
 272         return TRUE;
 273     }
 274     return FALSE;
 275 }
 276
 277 /* Find the index of a given lambda component name, or -1 if not found */
 278 static int lambda_components_find(const lambda_components_t *lc,
 279                                   const char                *name,
 280                                   size_t                     name_length)
 281 {
 282     int i;
 283
 284     for (i = 0; i < lc->N; i++)
 285     {
 286         if (strncmp(lc->names[i], name, name_length) == 0)
 287         {
 288             return i;
 289         }
 290     }
 291     return -1;
 292 }
 293
 294
 295
 296 /* initialize a lambda vector */
 297 static void lambda_vec_init(lambda_vec_t *lv, const lambda_components_t *lc)
 298 {
 299     snew(lv->val, lc->N);
 300     lv->index = -1;
 301     lv->dhdl  = -1;
 302     lv->lc    = lc;
 303 }
 304
 305 static void lambda_vec_destroy(lambda_vec_t *lv)
 306 {
 307     sfree(lv->val);
 308 }
 309
 310 static void lambda_vec_copy(lambda_vec_t *lv, const lambda_vec_t *orig)
 311 {
 312     int i;
 313
 314     lambda_vec_init(lv, orig->lc);
 315     lv->dhdl  = orig->dhdl;
 316     lv->index = orig->index;
 317     for (i = 0; i < lv->lc->N; i++)
 318     {
 319         lv->val[i] = orig->val[i];
 320     }
 321 }
 322
 323 /* write a lambda vec to a preallocated string */
 324 static void lambda_vec_print(const lambda_vec_t *lv, char *str, gmx_bool named)
 325 {
 326     int    i;
 327     size_t np;
 328
 329     str[0] = 0; /* reset the string */
 330     if (lv->dhdl < 0)
 331     {
 332         if (named)
 333         {
 334             str += sprintf(str, "delta H to ");
 335         }
 336         if (lv->lc->N > 1)
 337         {
 338             str += sprintf(str, "(");
 339         }
 340         for (i = 0; i < lv->lc->N; i++)
 341         {
 342             str += sprintf(str, "%g", lv->val[i]);
 343             if (i < lv->lc->N-1)
 344             {
 345                 str += sprintf(str, ", ");
 346             }
 347         }
 348         if (lv->lc->N > 1)
 349         {
 350             str += sprintf(str, ")");
 351         }
 352     }
 353     else
 354     {
 355         /* this lambda vector describes a derivative */
 356         str += sprintf(str, "dH/dl");
 357         if (strlen(lv->lc->names[lv->dhdl]) > 0)
 358         {
 359             str += sprintf(str, " (%s)", lv->lc->names[lv->dhdl]);
 360         }
 361     }
 362 }
 363
 364 /* write a shortened version of the lambda vec to a preallocated string */
 365 static void lambda_vec_print_short(const lambda_vec_t *lv, char *str)
 366 {
 367     int    i;
 368     size_t np;
 369
 370     if (lv->index >= 0)
 371     {
 372         sprintf(str, "%6d", lv->index);
 373     }
 374     else
 375     {
 376         if (lv->dhdl < 0)
 377         {
 378             sprintf(str, "%6.3f", lv->val[0]);
 379         }
 380         else
 381         {
 382             sprintf(str, "dH/dl[%d]", lv->dhdl);
 383         }
 384     }
 385 }
 386
 387 /* write an intermediate version of two lambda vecs to a preallocated string */
 388 static void lambda_vec_print_intermediate(const lambda_vec_t *a,
 389                                           const lambda_vec_t *b, char *str)
 390 {
 391     int    i;
 392     size_t np;
 393
 394     str[0] = 0;
 395     if ( (a->index >= 0) && (b->index >= 0) )
 396     {
 397         sprintf(str, "%6.3f", ((double)a->index+(double)b->index)/2.);
 398     }
 399     else
 400     {
 401         if ( (a->dhdl < 0) && (b->dhdl < 0) )
 402         {
 403             sprintf(str, "%6.3f", (a->val[0]+b->val[0])/2.);
 404         }
 405     }
 406 }
 407
 408
 409
 410 /* calculate the difference in lambda vectors: c = a-b.
 411    c must be initialized already, and a and b must describe non-derivative
 412    lambda points */
 413 static void lambda_vec_diff(const lambda_vec_t *a, const lambda_vec_t *b,
 414                             lambda_vec_t *c)
 415 {
 416     int i;
 417
 418     if ( (a->dhdl > 0)  || (b->dhdl > 0) )
 419     {
 420         gmx_fatal(FARGS,
 421                   "Trying to calculate the difference between derivatives instead of lambda points");
 422     }
 423     if ((a->lc != b->lc) || (a->lc != c->lc) )
 424     {
 425         gmx_fatal(FARGS,
 426                   "Trying to calculate the difference lambdas with differing basis set");
 427     }
 428     for (i = 0; i < a->lc->N; i++)
 429     {
 430         c->val[i] = a->val[i] - b->val[i];
 431     }
 432 }
 433
 434 /* calculate and return the absolute difference in lambda vectors: c = |a-b|.
 435    a and b must describe non-derivative lambda points */
 436 static double lambda_vec_abs_diff(const lambda_vec_t *a, const lambda_vec_t *b)
 437 {
 438     int    i;
 439     double ret = 0.;
 440
 441     if ( (a->dhdl > 0)  || (b->dhdl > 0) )
 442     {
 443         gmx_fatal(FARGS,
 444                   "Trying to calculate the difference between derivatives instead of lambda points");
 445     }
 446     if (a->lc != b->lc)
 447     {
 448         gmx_fatal(FARGS,
 449                   "Trying to calculate the difference lambdas with differing basis set");
 450     }
 451     for (i = 0; i < a->lc->N; i++)
 452     {
 453         double df = a->val[i] - b->val[i];
 454         ret += df*df;
 455     }
 456     return sqrt(ret);
 457 }
 458
 459
 460 /* check whether two lambda vectors are the same */
 461 static gmx_bool lambda_vec_same(const lambda_vec_t *a, const lambda_vec_t *b)
 462 {
 463     int i;
 464
 465     if (a->lc != b->lc)
 466     {
 467         return FALSE;
 468     }
 469     if (a->dhdl < 0)
 470     {
 471         for (i = 0; i < a->lc->N; i++)
 472         {
 473             if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
 474             {
 475                 return FALSE;
 476             }
 477         }
 478         return TRUE;
 479     }
 480     else
 481     {
 482         /* they're derivatives, so we check whether the indices match */
 483         return (a->dhdl == b->dhdl);
 484     }
 485 }
 486
 487 /* Compare the sort order of two foreign lambda vectors
 488
 489     returns 1 if a is 'bigger' than b,
 490     returns 0 if they're the same,
 491     returns -1 if a is 'smaller' than b.*/
 492 static gmx_bool lambda_vec_cmp_foreign(const lambda_vec_t *a,
 493                                        const lambda_vec_t *b)
 494 {
 495     int      i;
 496     double   norm_a    = 0, norm_b = 0;
 497     gmx_bool different = FALSE;
 498
 499     if (a->lc != b->lc)
 500     {
 501         gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
 502     }
 503     /* if either one has an index we sort based on that */
 504     if ((a->index >= 0) || (b->index >= 0))
 505     {
 506         if (a->index == b->index)
 507         {
 508             return 0;
 509         }
 510         return (a->index > b->index) ? 1 : -1;
 511     }
 512     if (a->dhdl >= 0 || b->dhdl >= 0)
 513     {
 514         /* lambda vectors that are derivatives always sort higher than those
 515            without derivatives */
 516         if ((a->dhdl >= 0)  != (b->dhdl >= 0) )
 517         {
 518             return (a->dhdl >= 0) ? 1 : -1;
 519         }
 520         return a->dhdl > b->dhdl;
 521     }
 522
 523     /* neither has an index, so we can only sort on the lambda components,
 524        which is only valid if there is one component */
 525     for (i = 0; i < a->lc->N; i++)
 526     {
 527         if (!gmx_within_tol(a->val[i], b->val[i], 10*GMX_REAL_EPS))
 528         {
 529             different = TRUE;
 530         }
 531         norm_a += a->val[i]*a->val[i];
 532         norm_b += b->val[i]*b->val[i];
 533     }
 534     if (!different)
 535     {
 536         return 0;
 537     }
 538     return norm_a > norm_b;
 539 }
 540
 541 /* Compare the sort order of two native lambda vectors
 542
 543     returns 1 if a is 'bigger' than b,
 544     returns 0 if they're the same,
 545     returns -1 if a is 'smaller' than b.*/
 546 static gmx_bool lambda_vec_cmp_native(const lambda_vec_t *a,
 547                                       const lambda_vec_t *b)
 548 {
 549     int i;
 550
 551     if (a->lc != b->lc)
 552     {
 553         gmx_fatal(FARGS, "Can't compare lambdas with differing basis sets");
 554     }
 555     /* if either one has an index we sort based on that */
 556     if ((a->index >= 0) || (b->index >= 0))
 557     {
 558         if (a->index == b->index)
 559         {
 560             return 0;
 561         }
 562         return (a->index > b->index) ? 1 : -1;
 563     }
 564     /* neither has an index, so we can only sort on the lambda components,
 565        which is only valid if there is one component */
 566     if (a->lc->N > 1)
 567     {
 568         gmx_fatal(FARGS,
 569                   "Can't compare lambdas with no index and > 1 component");
 570     }
 571     if (a->dhdl >= 0 || b->dhdl >= 0)
 572     {
 573         gmx_fatal(FARGS,
 574                   "Can't compare native lambdas that are derivatives");
 575     }
 576     if (gmx_within_tol(a->val[0], b->val[0], 10*GMX_REAL_EPS))
 577     {
 578         return 0;
 579     }
 580     return a->val[0] > b->val[0] ? 1 : -1;
 581 }
 582
 583
 584
 585
 586 static void hist_init(hist_t *h, int nhist, int *nbin)
 587 {
 588     int i;
 589     if (nhist > 2)
 590     {
 591         gmx_fatal(FARGS, "histogram with more than two sets of data!");
 592     }
 593     for (i = 0; i < nhist; i++)
 594     {
 595         snew(h->bin[i], nbin[i]);
 596         h->x0[i]      = 0;
 597         h->nbin[i]    = nbin[i];
 598         h->start_time = h->delta_time = 0;
 599         h->dx[i]      = 0;
 600     }
 601     h->sum   = 0;
 602     h->nhist = nhist;
 603 }
 604
 605 static void hist_destroy(hist_t *h)
 606 {
 607     sfree(h->bin);
 608 }
 609
 610
 611 static void xvg_init(xvg_t *ba)
 612 {
 613     ba->filename = NULL;
 614     ba->nset     = 0;
 615     ba->np_alloc = 0;
 616     ba->np       = NULL;
 617     ba->y        = NULL;
 618 }
 619
 620 static void samples_init(samples_t *s, lambda_vec_t *native_lambda,
 621                          lambda_vec_t *foreign_lambda, double temp,
 622                          gmx_bool derivative, const char *filename)
 623 {
 624     s->native_lambda  = native_lambda;
 625     s->foreign_lambda = foreign_lambda;
 626     s->temp           = temp;
 627     s->derivative     = derivative;
 628
 629     s->ndu        = 0;
 630     s->du         = NULL;
 631     s->t          = NULL;
 632     s->start_time = s->delta_time = 0;
 633     s->hist       = NULL;
 634     s->du_alloc   = NULL;
 635     s->t_alloc    = NULL;
 636     s->hist_alloc = NULL;
 637     s->ndu_alloc  = 0;
 638     s->nt_alloc   = 0;
 639
 640     s->ntot     = 0;
 641     s->filename = filename;
 642 }
 643
 644 static void sample_range_init(sample_range_t *r, samples_t *s)
 645 {
 646     r->start = 0;
 647     r->end   = s->ndu;
 648     r->use   = TRUE;
 649     r->s     = NULL;
 650 }
 651
 652 static void sample_coll_init(sample_coll_t *sc, lambda_vec_t *native_lambda,
 653                              lambda_vec_t *foreign_lambda, double temp)
 654 {
 655     sc->native_lambda  = native_lambda;
 656     sc->foreign_lambda = foreign_lambda;
 657     sc->temp           = temp;
 658
 659     sc->nsamples       = 0;
 660     sc->s              = NULL;
 661     sc->r              = NULL;
 662     sc->nsamples_alloc = 0;
 663
 664     sc->ntot = 0;
 665     sc->next = sc->prev = NULL;
 666 }
 667
 668 static void sample_coll_destroy(sample_coll_t *sc)
 669 {
 670     /* don't free the samples themselves */
 671     sfree(sc->r);
 672     sfree(sc->s);
 673 }
 674
 675
 676 static void lambda_data_init(lambda_data_t *l, lambda_vec_t *native_lambda,
 677                              double temp)
 678 {
 679     l->lambda = native_lambda;
 680     l->temp   = temp;
 681
 682     l->next = NULL;
 683     l->prev = NULL;
 684
 685     l->sc = &(l->sc_head);
 686
 687     sample_coll_init(l->sc, native_lambda, NULL, 0.);
 688     l->sc->next = l->sc;
 689     l->sc->prev = l->sc;
 690 }
 691
 692 static void barres_init(barres_t *br)
 693 {
 694     br->dg            = 0;
 695     br->dg_err        = 0;
 696     br->sa            = 0;
 697     br->sa_err        = 0;
 698     br->sb            = 0;
 699     br->sb_err        = 0;
 700     br->dg_stddev     = 0;
 701     br->dg_stddev_err = 0;
 702
 703     br->a = NULL;
 704     br->b = NULL;
 705 }
 706
 707
 708 /* calculate the total number of samples in a sample collection */
 709 static void sample_coll_calc_ntot(sample_coll_t *sc)
 710 {
 711     int i;
 712
 713     sc->ntot = 0;
 714     for (i = 0; i < sc->nsamples; i++)
 715     {
 716         if (sc->r[i].use)
 717         {
 718             if (sc->s[i]->hist)
 719             {
 720                 sc->ntot += sc->s[i]->ntot;
 721             }
 722             else
 723             {
 724                 sc->ntot += sc->r[i].end - sc->r[i].start;
 725             }
 726         }
 727     }
 728 }
 729
 730
 731 /* find the barsamples_t associated with a lambda that corresponds to
 732    a specific foreign lambda */
 733 static sample_coll_t *lambda_data_find_sample_coll(lambda_data_t *l,
 734                                                    lambda_vec_t  *foreign_lambda)
 735 {
 736     sample_coll_t *sc = l->sc->next;
 737
 738     while (sc != l->sc)
 739     {
 740         if (lambda_vec_same(sc->foreign_lambda, foreign_lambda))
 741         {
 742             return sc;
 743         }
 744         sc = sc->next;
 745     }
 746
 747     return NULL;
 748 }
 749
 750 /* insert li into an ordered list of lambda_colls */
 751 static void lambda_data_insert_sample_coll(lambda_data_t *l, sample_coll_t *sc)
 752 {
 753     sample_coll_t *scn = l->sc->next;
 754     while ( (scn != l->sc) )
 755     {
 756         if (lambda_vec_cmp_foreign(scn->foreign_lambda, sc->foreign_lambda) > 0)
 757         {
 758             break;
 759         }
 760         scn = scn->next;
 761     }
 762     /* now insert it before the found scn */
 763     sc->next        = scn;
 764     sc->prev        = scn->prev;
 765     scn->prev->next = sc;
 766     scn->prev       = sc;
 767 }
 768
 769 /* insert li into an ordered list of lambdas */
 770 static void lambda_data_insert_lambda(lambda_data_t *head, lambda_data_t *li)
 771 {
 772     lambda_data_t *lc = head->next;
 773     while (lc != head)
 774     {
 775         if (lambda_vec_cmp_native(lc->lambda, li->lambda) > 0)
 776         {
 777             break;
 778         }
 779         lc = lc->next;
 780     }
 781     /* now insert ourselves before the found lc */
 782     li->next       = lc;
 783     li->prev       = lc->prev;
 784     lc->prev->next = li;
 785     lc->prev       = li;
 786 }
 787
 788 /* insert a sample and a sample_range into a sample_coll. The
 789     samples are stored as a pointer, the range is copied. */
 790 static void sample_coll_insert_sample(sample_coll_t *sc, samples_t *s,
 791                                       sample_range_t *r)
 792 {
 793     /* first check if it belongs here */
 794     if (sc->temp != s->temp)
 795     {
 796         gmx_fatal(FARGS, "Temperatures in files %s and %s are not the same!",
 797                   s->filename, sc->next->s[0]->filename);
 798     }
 799     if (!lambda_vec_same(sc->native_lambda, s->native_lambda))
 800     {
 801         gmx_fatal(FARGS, "Native lambda in files %s and %s are not the same (and they should be)!",
 802                   s->filename, sc->next->s[0]->filename);
 803     }
 804     if (!lambda_vec_same(sc->foreign_lambda, s->foreign_lambda))
 805     {
 806         gmx_fatal(FARGS, "Foreign lambda in files %s and %s are not the same (and they should be)!",
 807                   s->filename, sc->next->s[0]->filename);
 808     }
 809
 810     /* check if there's room */
 811     if ( (sc->nsamples + 1) > sc->nsamples_alloc)
 812     {
 813         sc->nsamples_alloc = max(2*sc->nsamples_alloc, 2);
 814         srenew(sc->s, sc->nsamples_alloc);
 815         srenew(sc->r, sc->nsamples_alloc);
 816     }
 817     sc->s[sc->nsamples] = s;
 818     sc->r[sc->nsamples] = *r;
 819     sc->nsamples++;
 820
 821     sample_coll_calc_ntot(sc);
 822 }
 823
 824 /* insert a sample into a lambda_list, creating the right sample_coll if
 825    neccesary */
 826 static void lambda_data_list_insert_sample(lambda_data_t *head, samples_t *s)
 827 {
 828     gmx_bool       found = FALSE;
 829     sample_coll_t *sc;
 830     sample_range_t r;
 831
 832     lambda_data_t *l = head->next;
 833
 834     /* first search for the right lambda_data_t */
 835     while (l != head)
 836     {
 837         if (lambda_vec_same(l->lambda, s->native_lambda) )
 838         {
 839             found = TRUE;
 840             break;
 841         }
 842         l = l->next;
 843     }
 844
 845     if (!found)
 846     {
 847         snew(l, 1);                                     /* allocate a new one */
 848         lambda_data_init(l, s->native_lambda, s->temp); /* initialize it */
 849         lambda_data_insert_lambda(head, l);             /* add it to the list */
 850     }
 851
 852     /* now look for a sample collection */
 853     sc = lambda_data_find_sample_coll(l, s->foreign_lambda);
 854     if (!sc)
 855     {
 856         snew(sc, 1); /* allocate a new one */
 857         sample_coll_init(sc, s->native_lambda, s->foreign_lambda, s->temp);
 858         lambda_data_insert_sample_coll(l, sc);
 859     }
 860
 861     /* now insert the samples into the sample coll */
 862     sample_range_init(&r, s);
 863     sample_coll_insert_sample(sc, s, &r);
 864 }
 865
 866
 867 /* make a histogram out of a sample collection */
 868 static void sample_coll_make_hist(sample_coll_t *sc, int **bin,
 869                                   int *nbin_alloc, int *nbin,
 870                                   double *dx, double *xmin, int nbin_default)
 871 {
 872     int      i, j, k;
 873     gmx_bool dx_set   = FALSE;
 874     gmx_bool xmin_set = FALSE;
 875
 876     gmx_bool xmax_set      = FALSE;
 877     gmx_bool xmax_set_hard = FALSE; /* whether the xmax is bounded by the
 878                                        limits of a histogram */
 879     double   xmax          = -1;
 880
 881     /* first determine dx and xmin; try the histograms */
 882     for (i = 0; i < sc->nsamples; i++)
 883     {
 884         if (sc->s[i]->hist)
 885         {
 886             hist_t *hist = sc->s[i]->hist;
 887             for (k = 0; k < hist->nhist; k++)
 888             {
 889                 double hdx      = hist->dx[k];
 890                 double xmax_now = (hist->x0[k]+hist->nbin[k])*hdx;
 891
 892                 /* we use the biggest dx*/
 893                 if ( (!dx_set) || hist->dx[0] > *dx)
 894                 {
 895                     dx_set = TRUE;
 896                     *dx    = hist->dx[0];
 897                 }
 898                 if ( (!xmin_set) || (hist->x0[k]*hdx) < *xmin)
 899                 {
 900                     xmin_set = TRUE;
 901                     *xmin    = (hist->x0[k]*hdx);
 902                 }
 903
 904                 if ( (!xmax_set) || (xmax_now > xmax && !xmax_set_hard) )
 905                 {
 906                     xmax_set = TRUE;
 907                     xmax     = xmax_now;
 908                     if (hist->bin[k][hist->nbin[k]-1] != 0)
 909                     {
 910                         xmax_set_hard = TRUE;
 911                     }
 912                 }
 913                 if (hist->bin[k][hist->nbin[k]-1] != 0 && (xmax_now < xmax) )
 914                 {
 915                     xmax_set_hard = TRUE;
 916                     xmax          = xmax_now;
 917                 }
 918             }
 919         }
 920     }
 921     /* and the delta us */
 922     for (i = 0; i < sc->nsamples; i++)
 923     {
 924         if (sc->s[i]->ndu > 0)
 925         {
 926             /* determine min and max */
 927             int    starti  = sc->r[i].start;
 928             int    endi    = sc->r[i].end;
 929             double du_xmin = sc->s[i]->du[starti];
 930             double du_xmax = sc->s[i]->du[starti];
 931             for (j = starti+1; j < endi; j++)
 932             {
 933                 if (sc->s[i]->du[j] < du_xmin)
 934                 {
 935                     du_xmin = sc->s[i]->du[j];
 936                 }
 937                 if (sc->s[i]->du[j] > du_xmax)
 938                 {
 939                     du_xmax = sc->s[i]->du[j];
 940                 }
 941             }
 942
 943             /* and now change the limits */
 944             if ( (!xmin_set) || (du_xmin < *xmin) )
 945             {
 946                 xmin_set = TRUE;
 947                 *xmin    = du_xmin;
 948             }
 949             if ( (!xmax_set) || ((du_xmax > xmax) &&  !xmax_set_hard) )
 950             {
 951                 xmax_set = TRUE;
 952                 xmax     = du_xmax;
 953             }
 954         }
 955     }
 956
 957     if (!xmax_set || !xmin_set)
 958     {
 959         *nbin = 0;
 960         return;
 961     }
 962
 963
 964     if (!dx_set)
 965     {
 966         *nbin = nbin_default;
 967         *dx   = (xmax-(*xmin))/((*nbin)-2); /* -2 because we want the last bin to
 968                                                be 0, and we count from 0 */
 969     }
 970     else
 971     {
 972         *nbin = (xmax-(*xmin))/(*dx);
 973     }
 974
 975     if (*nbin > *nbin_alloc)
 976     {
 977         *nbin_alloc = *nbin;
 978         srenew(*bin, *nbin_alloc);
 979     }
 980
 981     /* reset the histogram */
 982     for (i = 0; i < (*nbin); i++)
 983     {
 984         (*bin)[i] = 0;
 985     }
 986
 987     /* now add the actual data */
 988     for (i = 0; i < sc->nsamples; i++)
 989     {
 990         if (sc->s[i]->hist)
 991         {
 992             hist_t *hist = sc->s[i]->hist;
 993             for (k = 0; k < hist->nhist; k++)
 994             {
 995                 double hdx       = hist->dx[k];
 996                 double xmin_hist = hist->x0[k]*hdx;
 997                 for (j = 0; j < hist->nbin[k]; j++)
 998                 {
 999                     /* calculate the bin corresponding to the middle of the
1000                        original bin */
1001                     double x     = hdx*(j+0.5) + xmin_hist;
1002                     int    binnr = (int)((x-(*xmin))/(*dx));
1003
1004                     if (binnr >= *nbin || binnr < 0)
1005                     {
1006                         binnr = (*nbin)-1;
1007                     }
1008
1009                     (*bin)[binnr] += hist->bin[k][j];
1010                 }
1011             }
1012         }
1013         else
1014         {
1015             int starti = sc->r[i].start;
1016             int endi   = sc->r[i].end;
1017             for (j = starti; j < endi; j++)
1018             {
1019                 int binnr = (int)((sc->s[i]->du[j]-(*xmin))/(*dx));
1020                 if (binnr >= *nbin || binnr < 0)
1021                 {
1022                     binnr = (*nbin)-1;
1023                 }
1024
1025                 (*bin)[binnr]++;
1026             }
1027         }
1028     }
1029 }
1030
1031 /* write a collection of histograms to a file */
1032 void sim_data_histogram(sim_data_t *sd, const char *filename,
1033                         int nbin_default, const output_env_t oenv)
1034 {
1035     char           label_x[STRLEN];
1036     const char    *dhdl    = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1037     const char    *title   = "N(\\DeltaH)";
1038     const char    *label_y = "Samples";
1039     FILE          *fp;
1040     lambda_data_t *bl;
1041     int            nsets     = 0;
1042     char         **setnames  = NULL;
1043     gmx_bool       first_set = FALSE;
1044     /* histogram data: */
1045     int           *hist       = NULL;
1046     int            nbin       = 0;
1047     int            nbin_alloc = 0;
1048     double         dx         = 0;
1049     double         min        = 0;
1050     int            i;
1051     lambda_data_t *bl_head = sd->lb;
1052
1053     printf("\nWriting histogram to %s\n", filename);
1054     sprintf(label_x, "\\DeltaH (%s)", unit_energy);
1055
1056     fp = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1057
1058     /* first get all the set names */
1059     bl = bl_head->next;
1060     /* iterate over all lambdas */
1061     while (bl != bl_head)
1062     {
1063         sample_coll_t *sc = bl->sc->next;
1064
1065         /* iterate over all samples */
1066         while (sc != bl->sc)
1067         {
1068             char buf[STRLEN], buf2[STRLEN];
1069
1070             nsets++;
1071             srenew(setnames, nsets);
1072             snew(setnames[nsets-1], STRLEN);
1073             if (sc->foreign_lambda->dhdl < 0)
1074             {
1075                 lambda_vec_print(sc->native_lambda, buf, FALSE);
1076                 lambda_vec_print(sc->foreign_lambda, buf2, FALSE);
1077                 sprintf(setnames[nsets-1], "N(%s(%s=%s) | %s=%s)",
1078                         deltag, lambda, buf2, lambda, buf);
1079             }
1080             else
1081             {
1082                 lambda_vec_print(sc->native_lambda, buf, FALSE);
1083                 sprintf(setnames[nsets-1], "N(%s | %s=%s)",
1084                         dhdl, lambda, buf);
1085             }
1086             sc = sc->next;
1087         }
1088
1089         bl = bl->next;
1090     }
1091     xvgr_legend(fp, nsets, (const char**)setnames, oenv);
1092
1093
1094     /* now make the histograms */
1095     bl = bl_head->next;
1096     /* iterate over all lambdas */
1097     while (bl != bl_head)
1098     {
1099         sample_coll_t *sc = bl->sc->next;
1100
1101         /* iterate over all samples */
1102         while (sc != bl->sc)
1103         {
1104             if (!first_set)
1105             {
1106                 xvgr_new_dataset(fp, 0, 0, NULL, oenv);
1107             }
1108
1109             sample_coll_make_hist(sc, &hist, &nbin_alloc, &nbin, &dx, &min,
1110                                   nbin_default);
1111
1112             for (i = 0; i < nbin; i++)
1113             {
1114                 double xmin = i*dx + min;
1115                 double xmax = (i+1)*dx + min;
1116
1117                 fprintf(fp, "%g %d\n%g %d\n", xmin, hist[i], xmax, hist[i]);
1118             }
1119
1120             first_set = FALSE;
1121             sc        = sc->next;
1122         }
1123
1124         bl = bl->next;
1125     }
1126
1127     if (hist)
1128     {
1129         sfree(hist);
1130     }
1131
1132     xvgrclose(fp);
1133 }
1134
1135 /* create a collection (array) of barres_t object given a ordered linked list
1136    of barlamda_t sample collections */
1137 static barres_t *barres_list_create(sim_data_t *sd, int *nres,
1138                                     gmx_bool use_dhdl)
1139 {
1140     lambda_data_t *bl;
1141     int            nlambda = 0;
1142     barres_t      *res;
1143     int            i;
1144     gmx_bool       dhdl    = FALSE;
1145     gmx_bool       first   = TRUE;
1146     lambda_data_t *bl_head = sd->lb;
1147
1148     /* first count the lambdas */
1149     bl = bl_head->next;
1150     while (bl != bl_head)
1151     {
1152         nlambda++;
1153         bl = bl->next;
1154     }
1155     snew(res, nlambda-1);
1156
1157     /* next put the right samples in the res */
1158     *nres = 0;
1159     bl    = bl_head->next->next; /* we start with the second one. */
1160     while (bl != bl_head)
1161     {
1162         sample_coll_t *sc, *scprev;
1163         barres_t      *br = &(res[*nres]);
1164         /* there is always a previous one. we search for that as a foreign
1165            lambda: */
1166         scprev = lambda_data_find_sample_coll(bl->prev, bl->lambda);
1167         sc     = lambda_data_find_sample_coll(bl, bl->prev->lambda);
1168
1169         barres_init(br);
1170
1171         if (use_dhdl)
1172         {
1173             /* we use dhdl */
1174
1175             scprev = lambda_data_find_sample_coll(bl->prev, bl->prev->lambda);
1176             sc     = lambda_data_find_sample_coll(bl, bl->lambda);
1177
1178             if (first)
1179             {
1180                 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");
1181                 dhdl = TRUE;
1182             }
1183             if (!dhdl)
1184             {
1185                 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");
1186             }
1187         }
1188         else if (!scprev && !sc)
1189         {
1190             gmx_fatal(FARGS, "There is no path from lambda=%f -> %f that is covered by foreign lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl by calculating the averages of dH/dl\nwith g_analyze and integrating them.\nAlternatively, use the -extp option if (and only if) the Hamiltonian\ndepends linearly on lambda, which is NOT normally the case.\n", bl->prev->lambda, bl->lambda);
1191         }
1192
1193         /* normal delta H */
1194         if (!scprev)
1195         {
1196             gmx_fatal(FARGS, "Could not find a set for foreign lambda = %f\nin the files for lambda = %f", bl->lambda, bl->prev->lambda);
1197         }
1198         if (!sc)
1199         {
1200             gmx_fatal(FARGS, "Could not find a set for foreign lambda = %f\nin the files for lambda = %f", bl->prev->lambda, bl->lambda);
1201         }
1202         br->a = scprev;
1203         br->b = sc;
1204
1205         first = FALSE;
1206         (*nres)++;
1207         bl = bl->next;
1208     }
1209     return res;
1210 }
1211
1212 /* estimate the maximum discretization error */
1213 static double barres_list_max_disc_err(barres_t *res, int nres)
1214 {
1215     int    i, j;
1216     double disc_err = 0.;
1217     double delta_lambda;
1218
1219     for (i = 0; i < nres; i++)
1220     {
1221         barres_t *br = &(res[i]);
1222
1223         delta_lambda = lambda_vec_abs_diff(br->b->native_lambda,
1224                                            br->a->native_lambda);
1225
1226         for (j = 0; j < br->a->nsamples; j++)
1227         {
1228             if (br->a->s[j]->hist)
1229             {
1230                 double Wfac = 1.;
1231                 if (br->a->s[j]->derivative)
1232                 {
1233                     Wfac =  delta_lambda;
1234                 }
1235
1236                 disc_err = max(disc_err, Wfac*br->a->s[j]->hist->dx[0]);
1237             }
1238         }
1239         for (j = 0; j < br->b->nsamples; j++)
1240         {
1241             if (br->b->s[j]->hist)
1242             {
1243                 double Wfac = 1.;
1244                 if (br->b->s[j]->derivative)
1245                 {
1246                     Wfac =  delta_lambda;
1247                 }
1248                 disc_err = max(disc_err, Wfac*br->b->s[j]->hist->dx[0]);
1249             }
1250         }
1251     }
1252     return disc_err;
1253 }
1254
1255
1256 /* impose start and end times on a sample collection, updating sample_ranges */
1257 static void sample_coll_impose_times(sample_coll_t *sc, double begin_t,
1258                                      double end_t)
1259 {
1260     int i;
1261     for (i = 0; i < sc->nsamples; i++)
1262     {
1263         samples_t      *s = sc->s[i];
1264         sample_range_t *r = &(sc->r[i]);
1265         if (s->hist)
1266         {
1267             double end_time = s->hist->delta_time*s->hist->sum +
1268                 s->hist->start_time;
1269             if (s->hist->start_time < begin_t || end_time > end_t)
1270             {
1271                 r->use = FALSE;
1272             }
1273         }
1274         else
1275         {
1276             if (!s->t)
1277             {
1278                 double end_time;
1279                 if (s->start_time < begin_t)
1280                 {
1281                     r->start = (int)((begin_t - s->start_time)/s->delta_time);
1282                 }
1283                 end_time = s->delta_time*s->ndu + s->start_time;
1284                 if (end_time > end_t)
1285                 {
1286                     r->end = (int)((end_t - s->start_time)/s->delta_time);
1287                 }
1288             }
1289             else
1290             {
1291                 int j;
1292                 for (j = 0; j < s->ndu; j++)
1293                 {
1294                     if (s->t[j] < begin_t)
1295                     {
1296                         r->start = j;
1297                     }
1298
1299                     if (s->t[j] >= end_t)
1300                     {
1301                         r->end = j;
1302                         break;
1303                     }
1304                 }
1305             }
1306             if (r->start > r->end)
1307             {
1308                 r->use = FALSE;
1309             }
1310         }
1311     }
1312     sample_coll_calc_ntot(sc);
1313 }
1314
1315 static void sim_data_impose_times(sim_data_t *sd, double begin, double end)
1316 {
1317     double         first_t, last_t;
1318     double         begin_t, end_t;
1319     lambda_data_t *lc;
1320     lambda_data_t *head = sd->lb;
1321     int            j;
1322
1323     if (begin <= 0 && end < 0)
1324     {
1325         return;
1326     }
1327
1328     /* first determine the global start and end times */
1329     first_t = -1;
1330     last_t  = -1;
1331     lc      = head->next;
1332     while (lc != head)
1333     {
1334         sample_coll_t *sc = lc->sc->next;
1335         while (sc != lc->sc)
1336         {
1337             for (j = 0; j < sc->nsamples; j++)
1338             {
1339                 double start_t, end_t;
1340
1341                 start_t = sc->s[j]->start_time;
1342                 end_t   =   sc->s[j]->start_time;
1343                 if (sc->s[j]->hist)
1344                 {
1345                     end_t += sc->s[j]->delta_time*sc->s[j]->hist->sum;
1346                 }
1347                 else
1348                 {
1349                     if (sc->s[j]->t)
1350                     {
1351                         end_t = sc->s[j]->t[sc->s[j]->ndu-1];
1352                     }
1353                     else
1354                     {
1355                         end_t += sc->s[j]->delta_time*sc->s[j]->ndu;
1356                     }
1357                 }
1358
1359                 if (start_t < first_t || first_t < 0)
1360                 {
1361                     first_t = start_t;
1362                 }
1363                 if (end_t > last_t)
1364                 {
1365                     last_t = end_t;
1366                 }
1367             }
1368             sc = sc->next;
1369         }
1370         lc = lc->next;
1371     }
1372
1373     /* calculate the actual times */
1374     if (begin > 0)
1375     {
1376         begin_t = begin;
1377     }
1378     else
1379     {
1380         begin_t = first_t;
1381     }
1382
1383     if (end > 0)
1384     {
1385         end_t = end;
1386     }
1387     else
1388     {
1389         end_t = last_t;
1390     }
1391     printf("\n   Samples in time interval: %.3f - %.3f\n", first_t, last_t);
1392
1393     if (begin_t > end_t)
1394     {
1395         return;
1396     }
1397     printf("Removing samples outside of: %.3f - %.3f\n", begin_t, end_t);
1398
1399     /* then impose them */
1400     lc = head->next;
1401     while (lc != head)
1402     {
1403         sample_coll_t *sc = lc->sc->next;
1404         while (sc != lc->sc)
1405         {
1406             sample_coll_impose_times(sc, begin_t, end_t);
1407             sc = sc->next;
1408         }
1409         lc = lc->next;
1410     }
1411 }
1412
1413
1414 /* create subsample i out of ni from an existing sample_coll */
1415 static gmx_bool sample_coll_create_subsample(sample_coll_t  *sc,
1416                                              sample_coll_t *sc_orig,
1417                                              int i, int ni)
1418 {
1419     int             j;
1420     int             hist_start, hist_end;
1421
1422     gmx_int64_t     ntot_start;
1423     gmx_int64_t     ntot_end;
1424     gmx_int64_t     ntot_so_far;
1425
1426     *sc = *sc_orig; /* just copy all fields */
1427
1428     /* allocate proprietary memory */
1429     snew(sc->s, sc_orig->nsamples);
1430     snew(sc->r, sc_orig->nsamples);
1431
1432     /* copy the samples */
1433     for (j = 0; j < sc_orig->nsamples; j++)
1434     {
1435         sc->s[j] = sc_orig->s[j];
1436         sc->r[j] = sc_orig->r[j]; /* copy the ranges too */
1437     }
1438
1439     /* now fix start and end fields */
1440     /* the casts avoid possible overflows */
1441     ntot_start  = (gmx_int64_t)(sc_orig->ntot*(double)i/(double)ni);
1442     ntot_end    = (gmx_int64_t)(sc_orig->ntot*(double)(i+1)/(double)ni);
1443     ntot_so_far = 0;
1444     for (j = 0; j < sc->nsamples; j++)
1445     {
1446         gmx_int64_t ntot_add;
1447         gmx_int64_t new_start, new_end;
1448
1449         if (sc->r[j].use)
1450         {
1451             if (sc->s[j]->hist)
1452             {
1453                 ntot_add = sc->s[j]->hist->sum;
1454             }
1455             else
1456             {
1457                 ntot_add = sc->r[j].end - sc->r[j].start;
1458             }
1459         }
1460         else
1461         {
1462             ntot_add = 0;
1463         }
1464
1465         if (!sc->s[j]->hist)
1466         {
1467             if (ntot_so_far < ntot_start)
1468             {
1469                 /* adjust starting point */
1470                 new_start = sc->r[j].start + (ntot_start - ntot_so_far);
1471             }
1472             else
1473             {
1474                 new_start = sc->r[j].start;
1475             }
1476             /* adjust end point */
1477             new_end = sc->r[j].start + (ntot_end - ntot_so_far);
1478             if (new_end > sc->r[j].end)
1479             {
1480                 new_end = sc->r[j].end;
1481             }
1482
1483             /* check if we're in range at all */
1484             if ( (new_end < new_start) || (new_start > sc->r[j].end) )
1485             {
1486                 new_start = 0;
1487                 new_end   = 0;
1488             }
1489             /* and write the new range */
1490             sc->r[j].start = (int)new_start;
1491             sc->r[j].end   = (int)new_end;
1492         }
1493         else
1494         {
1495             if (sc->r[j].use)
1496             {
1497                 double overlap;
1498                 double ntot_start_norm, ntot_end_norm;
1499                 /* calculate the amount of overlap of the
1500                    desired range (ntot_start -- ntot_end) onto
1501                    the histogram range (ntot_so_far -- ntot_so_far+ntot_add)*/
1502
1503                 /* first calculate normalized bounds
1504                    (where 0 is the start of the hist range, and 1 the end) */
1505                 ntot_start_norm = (ntot_start-ntot_so_far)/(double)ntot_add;
1506                 ntot_end_norm   = (ntot_end-ntot_so_far)/(double)ntot_add;
1507
1508                 /* now fix the boundaries */
1509                 ntot_start_norm = min(1, max(0., ntot_start_norm));
1510                 ntot_end_norm   = max(0, min(1., ntot_end_norm));
1511
1512                 /* and calculate the overlap */
1513                 overlap = ntot_end_norm - ntot_start_norm;
1514
1515                 if (overlap > 0.95) /* we allow for 5% slack */
1516                 {
1517                     sc->r[j].use = TRUE;
1518                 }
1519                 else if (overlap < 0.05)
1520                 {
1521                     sc->r[j].use = FALSE;
1522                 }
1523                 else
1524                 {
1525                     return FALSE;
1526                 }
1527             }
1528         }
1529         ntot_so_far += ntot_add;
1530     }
1531     sample_coll_calc_ntot(sc);
1532
1533     return TRUE;
1534 }
1535
1536 /* calculate minimum and maximum work values in sample collection */
1537 static void sample_coll_min_max(sample_coll_t *sc, double Wfac,
1538                                 double *Wmin, double *Wmax)
1539 {
1540     int i, j;
1541
1542     *Wmin = FLT_MAX;
1543     *Wmax = -FLT_MAX;
1544
1545     for (i = 0; i < sc->nsamples; i++)
1546     {
1547         samples_t      *s = sc->s[i];
1548         sample_range_t *r = &(sc->r[i]);
1549         if (r->use)
1550         {
1551             if (!s->hist)
1552             {
1553                 for (j = r->start; j < r->end; j++)
1554                 {
1555                     *Wmin = min(*Wmin, s->du[j]*Wfac);
1556                     *Wmax = max(*Wmax, s->du[j]*Wfac);
1557                 }
1558             }
1559             else
1560             {
1561                 int    hd = 0; /* determine the histogram direction: */
1562                 double dx;
1563                 if ( (s->hist->nhist > 1) && (Wfac < 0) )
1564                 {
1565                     hd = 1;
1566                 }
1567                 dx = s->hist->dx[hd];
1568
1569                 for (j = s->hist->nbin[hd]-1; j >= 0; j--)
1570                 {
1571                     *Wmin = min(*Wmin, Wfac*(s->hist->x0[hd])*dx);
1572                     *Wmax = max(*Wmax, Wfac*(s->hist->x0[hd])*dx);
1573                     /* look for the highest value bin with values */
1574                     if (s->hist->bin[hd][j] > 0)
1575                     {
1576                         *Wmin = min(*Wmin, Wfac*(j+s->hist->x0[hd]+1)*dx);
1577                         *Wmax = max(*Wmax, Wfac*(j+s->hist->x0[hd]+1)*dx);
1578                         break;
1579                     }
1580                 }
1581             }
1582         }
1583     }
1584 }
1585
1586 /* Initialize a sim_data structure */
1587 static void sim_data_init(sim_data_t *sd)
1588 {
1589     /* make linked list */
1590     sd->lb       = &(sd->lb_head);
1591     sd->lb->next = sd->lb;
1592     sd->lb->prev = sd->lb;
1593
1594     lambda_components_init(&(sd->lc));
1595 }
1596
1597
1598 static double calc_bar_sum(int n, const double *W, double Wfac, double sbMmDG)
1599 {
1600     int    i;
1601     double sum;
1602
1603     sum = 0;
1604
1605     for (i = 0; i < n; i++)
1606     {
1607         sum += 1./(1. + exp(Wfac*W[i] + sbMmDG));
1608     }
1609
1610     return sum;
1611 }
1612
1613 /* calculate the BAR average given a histogram
1614
1615     if type== 0, calculate the best estimate for the average,
1616     if type==-1, calculate the minimum possible value given the histogram
1617     if type== 1, calculate the maximum possible value given the histogram */
1618 static double calc_bar_sum_hist(const hist_t *hist, double Wfac, double sbMmDG,
1619                                 int type)
1620 {
1621     double sum = 0.;
1622     int    i;
1623     int    max;
1624     /* normalization factor multiplied with bin width and
1625        number of samples (we normalize through M): */
1626     double normdx = 1.;
1627     int    hd     = 0; /* determine the histogram direction: */
1628     double dx;
1629
1630     if ( (hist->nhist > 1) && (Wfac < 0) )
1631     {
1632         hd = 1;
1633     }
1634     dx  = hist->dx[hd];
1635     max = hist->nbin[hd]-1;
1636     if (type == 1)
1637     {
1638         max = hist->nbin[hd]; /* we also add whatever was out of range */
1639     }
1640
1641     for (i = 0; i < max; i++)
1642     {
1643         double x    = Wfac*((i+hist->x0[hd])+0.5)*dx; /* bin middle */
1644         double pxdx = hist->bin[0][i]*normdx;         /* p(x)dx */
1645
1646         sum += pxdx/(1. + exp(x + sbMmDG));
1647     }
1648
1649     return sum;
1650 }
1651
1652 static double calc_bar_lowlevel(sample_coll_t *ca, sample_coll_t *cb,
1653                                 double temp, double tol, int type)
1654 {
1655     double kT, beta, M;
1656     double DG;
1657     int    i, j;
1658     double Wfac1, Wfac2, Wmin, Wmax;
1659     double DG0, DG1, DG2, dDG1;
1660     double sum1, sum2;
1661     double n1, n2; /* numbers of samples as doubles */
1662
1663     kT   = BOLTZ*temp;
1664     beta = 1/kT;
1665
1666     /* count the numbers of samples */
1667     n1 = ca->ntot;
1668     n2 = cb->ntot;
1669
1670     M = log(n1/n2);
1671
1672     /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1673     if (ca->foreign_lambda->dhdl < 0)
1674     {
1675         /* this is the case when the delta U were calculated directly
1676            (i.e. we're not scaling dhdl) */
1677         Wfac1 = beta;
1678         Wfac2 = beta;
1679     }
1680     else
1681     {
1682         /* we're using dhdl, so delta_lambda needs to be a
1683            multiplication factor.  */
1684         /*double delta_lambda=cb->native_lambda-ca->native_lambda;*/
1685         double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1686                                                   ca->native_lambda);
1687         if (cb->native_lambda->lc->N > 1)
1688         {
1689             gmx_fatal(FARGS,
1690                       "Can't (yet) do multi-component dhdl interpolation");
1691         }
1692
1693         Wfac1 =  beta*delta_lambda;
1694         Wfac2 = -beta*delta_lambda;
1695     }
1696
1697     if (beta < 1)
1698     {
1699         /* We print the output both in kT and kJ/mol.
1700          * Here we determine DG in kT, so when beta < 1
1701          * the precision has to be increased.
1702          */
1703         tol *= beta;
1704     }
1705
1706     /* Calculate minimum and maximum work to give an initial estimate of
1707      * delta G  as their average.
1708      */
1709     {
1710         double Wmin1, Wmin2, Wmax1, Wmax2;
1711         sample_coll_min_max(ca, Wfac1, &Wmin1, &Wmax1);
1712         sample_coll_min_max(cb, Wfac2, &Wmin2, &Wmax2);
1713
1714         Wmin = min(Wmin1, Wmin2);
1715         Wmax = max(Wmax1, Wmax2);
1716     }
1717
1718     DG0 = Wmin;
1719     DG2 = Wmax;
1720
1721     if (debug)
1722     {
1723         fprintf(debug, "DG %9.5f %9.5f\n", DG0, DG2);
1724     }
1725     /* We approximate by bisection: given our initial estimates
1726        we keep checking whether the halfway point is greater or
1727        smaller than what we get out of the BAR averages.
1728
1729        For the comparison we can use twice the tolerance. */
1730     while (DG2 - DG0 > 2*tol)
1731     {
1732         DG1 = 0.5*(DG0 + DG2);
1733
1734         /* calculate the BAR averages */
1735         dDG1 = 0.;
1736
1737         for (i = 0; i < ca->nsamples; i++)
1738         {
1739             samples_t      *s = ca->s[i];
1740             sample_range_t *r = &(ca->r[i]);
1741             if (r->use)
1742             {
1743                 if (s->hist)
1744                 {
1745                     dDG1 += calc_bar_sum_hist(s->hist, Wfac1, (M-DG1), type);
1746                 }
1747                 else
1748                 {
1749                     dDG1 += calc_bar_sum(r->end - r->start, s->du + r->start,
1750                                          Wfac1, (M-DG1));
1751                 }
1752             }
1753         }
1754         for (i = 0; i < cb->nsamples; i++)
1755         {
1756             samples_t      *s = cb->s[i];
1757             sample_range_t *r = &(cb->r[i]);
1758             if (r->use)
1759             {
1760                 if (s->hist)
1761                 {
1762                     dDG1 -= calc_bar_sum_hist(s->hist, Wfac2, -(M-DG1), type);
1763                 }
1764                 else
1765                 {
1766                     dDG1 -= calc_bar_sum(r->end - r->start, s->du + r->start,
1767                                          Wfac2, -(M-DG1));
1768                 }
1769             }
1770         }
1771
1772         if (dDG1 < 0)
1773         {
1774             DG0 = DG1;
1775         }
1776         else
1777         {
1778             DG2 = DG1;
1779         }
1780         if (debug)
1781         {
1782             fprintf(debug, "DG %9.5f %9.5f\n", DG0, DG2);
1783         }
1784     }
1785
1786     return 0.5*(DG0 + DG2);
1787 }
1788
1789 static void calc_rel_entropy(sample_coll_t *ca, sample_coll_t *cb,
1790                              double temp, double dg, double *sa, double *sb)
1791 {
1792     int    i, j;
1793     double W_ab = 0.;
1794     double W_ba = 0.;
1795     double kT, beta;
1796     double Wfac1, Wfac2;
1797     double n1, n2;
1798
1799     kT   = BOLTZ*temp;
1800     beta = 1/kT;
1801
1802     /* count the numbers of samples */
1803     n1 = ca->ntot;
1804     n2 = cb->ntot;
1805
1806     /* to ensure the work values are the same as during the delta_G */
1807     /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1808     if (ca->foreign_lambda->dhdl < 0)
1809     {
1810         /* this is the case when the delta U were calculated directly
1811            (i.e. we're not scaling dhdl) */
1812         Wfac1 = beta;
1813         Wfac2 = beta;
1814     }
1815     else
1816     {
1817         /* we're using dhdl, so delta_lambda needs to be a
1818            multiplication factor.  */
1819         double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1820                                                   ca->native_lambda);
1821         Wfac1 =  beta*delta_lambda;
1822         Wfac2 = -beta*delta_lambda;
1823     }
1824
1825     /* first calculate the average work in both directions */
1826     for (i = 0; i < ca->nsamples; i++)
1827     {
1828         samples_t      *s = ca->s[i];
1829         sample_range_t *r = &(ca->r[i]);
1830         if (r->use)
1831         {
1832             if (!s->hist)
1833             {
1834                 for (j = r->start; j < r->end; j++)
1835                 {
1836                     W_ab += Wfac1*s->du[j];
1837                 }
1838             }
1839             else
1840             {
1841                 /* normalization factor multiplied with bin width and
1842                    number of samples (we normalize through M): */
1843                 double normdx = 1.;
1844                 double dx;
1845                 int    hd = 0; /* histogram direction */
1846                 if ( (s->hist->nhist > 1) && (Wfac1 < 0) )
1847                 {
1848                     hd = 1;
1849                 }
1850                 dx = s->hist->dx[hd];
1851
1852                 for (j = 0; j < s->hist->nbin[0]; j++)
1853                 {
1854                     double x    = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1855                     double pxdx = s->hist->bin[0][j]*normdx;         /* p(x)dx */
1856                     W_ab += pxdx*x;
1857                 }
1858             }
1859         }
1860     }
1861     W_ab /= n1;
1862
1863     for (i = 0; i < cb->nsamples; i++)
1864     {
1865         samples_t      *s = cb->s[i];
1866         sample_range_t *r = &(cb->r[i]);
1867         if (r->use)
1868         {
1869             if (!s->hist)
1870             {
1871                 for (j = r->start; j < r->end; j++)
1872                 {
1873                     W_ba += Wfac1*s->du[j];
1874                 }
1875             }
1876             else
1877             {
1878                 /* normalization factor multiplied with bin width and
1879                    number of samples (we normalize through M): */
1880                 double normdx = 1.;
1881                 double dx;
1882                 int    hd = 0; /* histogram direction */
1883                 if ( (s->hist->nhist > 1) && (Wfac2 < 0) )
1884                 {
1885                     hd = 1;
1886                 }
1887                 dx = s->hist->dx[hd];
1888
1889                 for (j = 0; j < s->hist->nbin[0]; j++)
1890                 {
1891                     double x    = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1892                     double pxdx = s->hist->bin[0][j]*normdx;         /* p(x)dx */
1893                     W_ba += pxdx*x;
1894                 }
1895             }
1896         }
1897     }
1898     W_ba /= n2;
1899
1900     /* then calculate the relative entropies */
1901     *sa = (W_ab - dg);
1902     *sb = (W_ba + dg);
1903 }
1904
1905 static void calc_dg_stddev(sample_coll_t *ca, sample_coll_t *cb,
1906                            double temp, double dg, double *stddev)
1907 {
1908     int    i, j;
1909     double M;
1910     double sigmafact = 0.;
1911     double kT, beta;
1912     double Wfac1, Wfac2;
1913     double n1, n2;
1914
1915     kT   = BOLTZ*temp;
1916     beta = 1/kT;
1917
1918     /* count the numbers of samples */
1919     n1 = ca->ntot;
1920     n2 = cb->ntot;
1921
1922     /* to ensure the work values are the same as during the delta_G */
1923     /*if (!lambda_vec_same(ca->native_lambda, ca->foreign_lambda))*/
1924     if (ca->foreign_lambda->dhdl < 0)
1925     {
1926         /* this is the case when the delta U were calculated directly
1927            (i.e. we're not scaling dhdl) */
1928         Wfac1 = beta;
1929         Wfac2 = beta;
1930     }
1931     else
1932     {
1933         /* we're using dhdl, so delta_lambda needs to be a
1934            multiplication factor.  */
1935         double delta_lambda = lambda_vec_abs_diff(cb->native_lambda,
1936                                                   ca->native_lambda);
1937         Wfac1 =  beta*delta_lambda;
1938         Wfac2 = -beta*delta_lambda;
1939     }
1940
1941     M = log(n1/n2);
1942
1943
1944     /* calculate average in both directions */
1945     for (i = 0; i < ca->nsamples; i++)
1946     {
1947         samples_t      *s = ca->s[i];
1948         sample_range_t *r = &(ca->r[i]);
1949         if (r->use)
1950         {
1951             if (!s->hist)
1952             {
1953                 for (j = r->start; j < r->end; j++)
1954                 {
1955                     sigmafact += 1./(2. + 2.*cosh((M + Wfac1*s->du[j] - dg)));
1956                 }
1957             }
1958             else
1959             {
1960                 /* normalization factor multiplied with bin width and
1961                    number of samples (we normalize through M): */
1962                 double normdx = 1.;
1963                 double dx;
1964                 int    hd = 0; /* histogram direction */
1965                 if ( (s->hist->nhist > 1) && (Wfac1 < 0) )
1966                 {
1967                     hd = 1;
1968                 }
1969                 dx = s->hist->dx[hd];
1970
1971                 for (j = 0; j < s->hist->nbin[0]; j++)
1972                 {
1973                     double x    = Wfac1*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
1974                     double pxdx = s->hist->bin[0][j]*normdx;         /* p(x)dx */
1975
1976                     sigmafact += pxdx/(2. + 2.*cosh((M + x - dg)));
1977                 }
1978             }
1979         }
1980     }
1981     for (i = 0; i < cb->nsamples; i++)
1982     {
1983         samples_t      *s = cb->s[i];
1984         sample_range_t *r = &(cb->r[i]);
1985         if (r->use)
1986         {
1987             if (!s->hist)
1988             {
1989                 for (j = r->start; j < r->end; j++)
1990                 {
1991                     sigmafact += 1./(2. + 2.*cosh((M - Wfac2*s->du[j] - dg)));
1992                 }
1993             }
1994             else
1995             {
1996                 /* normalization factor multiplied with bin width and
1997                    number of samples (we normalize through M): */
1998                 double normdx = 1.;
1999                 double dx;
2000                 int    hd = 0; /* histogram direction */
2001                 if ( (s->hist->nhist > 1) && (Wfac2 < 0) )
2002                 {
2003                     hd = 1;
2004                 }
2005                 dx = s->hist->dx[hd];
2006
2007                 for (j = 0; j < s->hist->nbin[0]; j++)
2008                 {
2009                     double x    = Wfac2*((j+s->hist->x0[0])+0.5)*dx; /*bin ctr*/
2010                     double pxdx = s->hist->bin[0][j]*normdx;         /* p(x)dx */
2011
2012                     sigmafact += pxdx/(2. + 2.*cosh((M - x - dg)));
2013                 }
2014             }
2015         }
2016     }
2017
2018     sigmafact /= (n1 + n2);
2019
2020
2021     /* Eq. 10 from
2022        Shirts, Bair, Hooker & Pande, Phys. Rev. Lett 91, 140601 (2003): */
2023     *stddev = sqrt(((1./sigmafact) - ( (n1+n2)/n1 + (n1+n2)/n2 )));
2024 }
2025
2026
2027
2028 static void calc_bar(barres_t *br, double tol,
2029                      int npee_min, int npee_max, gmx_bool *bEE,
2030                      double *partsum)
2031 {
2032     int      npee, p;
2033     double   dg_sig2, sa_sig2, sb_sig2, stddev_sig2; /* intermediate variance values
2034                                                         for calculated quantities */
2035     int      nsample1, nsample2;
2036     double   temp = br->a->temp;
2037     int      i, j;
2038     double   dg_min, dg_max;
2039     gmx_bool have_hist = FALSE;
2040
2041     br->dg = calc_bar_lowlevel(br->a, br->b, temp, tol, 0);
2042
2043     br->dg_disc_err      = 0.;
2044     br->dg_histrange_err = 0.;
2045
2046     /* check if there are histograms */
2047     for (i = 0; i < br->a->nsamples; i++)
2048     {
2049         if (br->a->r[i].use && br->a->s[i]->hist)
2050         {
2051             have_hist = TRUE;
2052             break;
2053         }
2054     }
2055     if (!have_hist)
2056     {
2057         for (i = 0; i < br->b->nsamples; i++)
2058         {
2059             if (br->b->r[i].use && br->b->s[i]->hist)
2060             {
2061                 have_hist = TRUE;
2062                 break;
2063             }
2064         }
2065     }
2066
2067     /* calculate histogram-specific errors */
2068     if (have_hist)
2069     {
2070         dg_min = calc_bar_lowlevel(br->a, br->b, temp, tol, -1);
2071         dg_max = calc_bar_lowlevel(br->a, br->b, temp, tol, 1);
2072
2073         if (fabs(dg_max - dg_min) > GMX_REAL_EPS*10)
2074         {
2075             /* the histogram range  error is the biggest of the differences
2076                between the best estimate and the extremes */
2077             br->dg_histrange_err = fabs(dg_max - dg_min);
2078         }
2079         br->dg_disc_err = 0.;
2080         for (i = 0; i < br->a->nsamples; i++)
2081         {
2082             if (br->a->s[i]->hist)
2083             {
2084                 br->dg_disc_err = max(br->dg_disc_err, br->a->s[i]->hist->dx[0]);
2085             }
2086         }
2087         for (i = 0; i < br->b->nsamples; i++)
2088         {
2089             if (br->b->s[i]->hist)
2090             {
2091                 br->dg_disc_err = max(br->dg_disc_err, br->b->s[i]->hist->dx[0]);
2092             }
2093         }
2094     }
2095     calc_rel_entropy(br->a, br->b, temp, br->dg, &(br->sa), &(br->sb));
2096
2097     calc_dg_stddev(br->a, br->b, temp, br->dg, &(br->dg_stddev) );
2098
2099     dg_sig2     = 0;
2100     sa_sig2     = 0;
2101     sb_sig2     = 0;
2102     stddev_sig2 = 0;
2103
2104     *bEE = TRUE;
2105     {
2106         sample_coll_t ca, cb;
2107
2108         /* initialize the samples */
2109         sample_coll_init(&ca, br->a->native_lambda, br->a->foreign_lambda,
2110                          br->a->temp);
2111         sample_coll_init(&cb, br->b->native_lambda, br->b->foreign_lambda,
2112                          br->b->temp);
2113
2114         for (npee = npee_min; npee <= npee_max; npee++)
2115         {
2116             double dgs      = 0;
2117             double dgs2     = 0;
2118             double dsa      = 0;
2119             double dsb      = 0;
2120             double dsa2     = 0;
2121             double dsb2     = 0;
2122             double dstddev  = 0;
2123             double dstddev2 = 0;
2124
2125
2126             for (p = 0; p < npee; p++)
2127             {
2128                 double   dgp;
2129                 double   stddevc;
2130                 double   sac, sbc;
2131                 gmx_bool cac, cbc;
2132
2133                 cac = sample_coll_create_subsample(&ca, br->a, p, npee);
2134                 cbc = sample_coll_create_subsample(&cb, br->b, p, npee);
2135
2136                 if (!cac || !cbc)
2137                 {
2138                     printf("WARNING: histogram number incompatible with block number for averaging: can't do error estimate\n");
2139                     *bEE = FALSE;
2140                     if (cac)
2141                     {
2142                         sample_coll_destroy(&ca);
2143                     }
2144                     if (cbc)
2145                     {
2146                         sample_coll_destroy(&cb);
2147                     }
2148                     return;
2149                 }
2150
2151                 dgp   = calc_bar_lowlevel(&ca, &cb, temp, tol, 0);
2152                 dgs  += dgp;
2153                 dgs2 += dgp*dgp;
2154
2155                 partsum[npee*(npee_max+1)+p] += dgp;
2156
2157                 calc_rel_entropy(&ca, &cb, temp, dgp, &sac, &sbc);
2158                 dsa  += sac;
2159                 dsa2 += sac*sac;
2160                 dsb  += sbc;
2161                 dsb2 += sbc*sbc;
2162                 calc_dg_stddev(&ca, &cb, temp, dgp, &stddevc );
2163
2164                 dstddev  += stddevc;
2165                 dstddev2 += stddevc*stddevc;
2166
2167                 sample_coll_destroy(&ca);
2168                 sample_coll_destroy(&cb);
2169             }
2170             dgs     /= npee;
2171             dgs2    /= npee;
2172             dg_sig2 += (dgs2-dgs*dgs)/(npee-1);
2173
2174             dsa     /= npee;
2175             dsa2    /= npee;
2176             dsb     /= npee;
2177             dsb2    /= npee;
2178             sa_sig2 += (dsa2-dsa*dsa)/(npee-1);
2179             sb_sig2 += (dsb2-dsb*dsb)/(npee-1);
2180
2181             dstddev     /= npee;
2182             dstddev2    /= npee;
2183             stddev_sig2 += (dstddev2-dstddev*dstddev)/(npee-1);
2184         }
2185         br->dg_err        = sqrt(dg_sig2/(npee_max - npee_min + 1));
2186         br->sa_err        = sqrt(sa_sig2/(npee_max - npee_min + 1));
2187         br->sb_err        = sqrt(sb_sig2/(npee_max - npee_min + 1));
2188         br->dg_stddev_err = sqrt(stddev_sig2/(npee_max - npee_min + 1));
2189     }
2190 }
2191
2192
2193 static double bar_err(int nbmin, int nbmax, const double *partsum)
2194 {
2195     int    nb, b;
2196     double svar, s, s2, dg;
2197
2198     svar = 0;
2199     for (nb = nbmin; nb <= nbmax; nb++)
2200     {
2201         s  = 0;
2202         s2 = 0;
2203         for (b = 0; b < nb; b++)
2204         {
2205             dg  = partsum[nb*(nbmax+1)+b];
2206             s  += dg;
2207             s2 += dg*dg;
2208         }
2209         s    /= nb;
2210         s2   /= nb;
2211         svar += (s2 - s*s)/(nb - 1);
2212     }
2213
2214     return sqrt(svar/(nbmax + 1 - nbmin));
2215 }
2216
2217
2218 /* Seek the end of an identifier (consecutive non-spaces), followed by
2219    an optional number of spaces or '='-signs. Returns a pointer to the
2220    first non-space value found after that. Returns NULL if the string
2221    ends before that.
2222  */
2223 static const char *find_value(const char *str)
2224 {
2225     gmx_bool name_end_found = FALSE;
2226
2227     /* if the string is a NULL pointer, return a NULL pointer. */
2228     if (str == NULL)
2229     {
2230         return NULL;
2231     }
2232     while (*str != '\0')
2233     {
2234         /* first find the end of the name */
2235         if (!name_end_found)
2236         {
2237             if (isspace(*str) || (*str == '=') )
2238             {
2239                 name_end_found = TRUE;
2240             }
2241         }
2242         else
2243         {
2244             if (!( isspace(*str) || (*str == '=') ))
2245             {
2246                 return str;
2247             }
2248         }
2249         str++;
2250     }
2251     return NULL;
2252 }
2253
2254
2255
2256 /* read a vector-notation description of a lambda vector */
2257 static gmx_bool read_lambda_compvec(const char                *str,
2258                                     lambda_vec_t              *lv,
2259                                     const lambda_components_t *lc_in,
2260                                     lambda_components_t       *lc_out,
2261                                     const char               **end,
2262                                     const char                *fn)
2263 {
2264     gmx_bool    initialize_lc = FALSE; /* whether to initialize the lambda
2265                                           components, or to check them */
2266     gmx_bool    start_reached = FALSE; /* whether the start of component names
2267                                           has been reached */
2268     gmx_bool    vector        = FALSE; /* whether there are multiple components */
2269     int         n             = 0;     /* current component number */
2270     const char *val_start     = NULL;  /* start of the component name, or NULL
2271                                           if not in a value */
2272     char       *strtod_end;
2273     gmx_bool    OK = TRUE;
2274
2275     if (end)
2276     {
2277         *end = str;
2278     }
2279
2280
2281     if (lc_out && lc_out->N == 0)
2282     {
2283         initialize_lc = TRUE;
2284     }
2285
2286     if (lc_in == NULL)
2287     {
2288         lc_in = lc_out;
2289     }
2290
2291     while (1)
2292     {
2293         if (!start_reached)
2294         {
2295             if (isalnum(*str))
2296             {
2297                 vector        = FALSE;
2298                 start_reached = TRUE;
2299                 val_start     = str;
2300             }
2301             else if (*str == '(')
2302             {
2303                 vector        = TRUE;
2304                 start_reached = TRUE;
2305             }
2306             else if (!isspace(*str))
2307             {
2308                 gmx_fatal(FARGS, "Error in lambda components in %s", fn);
2309             }
2310         }
2311         else
2312         {
2313             if (val_start)
2314             {
2315                 if (isspace(*str) || *str == ')' || *str == ',' || *str == '\0')
2316                 {
2317                     /* end of value */
2318                     if (lv == NULL)
2319                     {
2320                         if (initialize_lc)
2321                         {
2322                             lambda_components_add(lc_out, val_start,
2323                                                   (str-val_start));
2324                         }
2325                         else
2326                         {
2327                             if (!lambda_components_check(lc_out, n, val_start,
2328                                                          (str-val_start)))
2329                             {
2330                                 return FALSE;
2331                             }
2332                         }
2333                     }
2334                     else
2335                     {
2336                         /* add a vector component to lv */
2337                         lv->val[n] = strtod(val_start, &strtod_end);
2338                         if (val_start == strtod_end)
2339                         {
2340                             gmx_fatal(FARGS,
2341                                       "Error reading lambda vector in %s", fn);
2342                         }
2343                     }
2344                     /* reset for the next identifier */
2345                     val_start = NULL;
2346                     n++;
2347                     if (!vector)
2348                     {
2349                         return OK;
2350                     }
2351                 }
2352             }
2353             else if (isalnum(*str))
2354             {
2355                 val_start = str;
2356             }
2357             if (*str == ')')
2358             {
2359                 str++;
2360                 if (end)
2361                 {
2362                     *end = str;
2363                 }
2364                 if (!vector)
2365                 {
2366                     gmx_fatal(FARGS, "Error in lambda components in %s", fn);
2367                 }
2368                 else
2369                 {
2370                     if (n == lc_in->N)
2371                     {
2372                         return OK;
2373                     }
2374                     else if (lv == NULL)
2375                     {
2376                         return FALSE;
2377                     }
2378                     else
2379                     {
2380                         gmx_fatal(FARGS, "Incomplete lambda vector data in %s",
2381                                   fn);
2382                         return FALSE;
2383                     }
2384
2385                 }
2386             }
2387         }
2388         if (*str == '\0')
2389         {
2390             break;
2391         }
2392         str++;
2393         if (end)
2394         {
2395             *end = str;
2396         }
2397     }
2398     if (vector)
2399     {
2400         gmx_fatal(FARGS, "Incomplete lambda components data in %s", fn);
2401         return FALSE;
2402     }
2403     return OK;
2404 }
2405
2406 /* read and check the component names from a string */
2407 static gmx_bool read_lambda_components(const char          *str,
2408                                        lambda_components_t *lc,
2409                                        const char         **end,
2410                                        const char          *fn)
2411 {
2412     return read_lambda_compvec(str, NULL, NULL, lc, end, fn);
2413 }
2414
2415 /* read an initialized lambda vector from a string */
2416 static gmx_bool read_lambda_vector(const char   *str,
2417                                    lambda_vec_t *lv,
2418                                    const char  **end,
2419                                    const char   *fn)
2420 {
2421     return read_lambda_compvec(str, lv, lv->lc, NULL, end, fn);
2422 }
2423
2424
2425
2426 /* deduce lambda value from legend.
2427     fn = the file name
2428     legend = the legend string
2429     ba = the xvg data
2430     lam = the initialized lambda vector
2431    returns whether to use the data in this set.
2432  */
2433 static gmx_bool legend2lambda(const char   *fn,
2434                               const char   *legend,
2435                               lambda_vec_t *lam)
2436 {
2437     double        lambda = 0;
2438     const char   *ptr    = NULL, *ptr2 = NULL;
2439     gmx_bool      ok     = FALSE;
2440     gmx_bool      bdhdl  = FALSE;
2441     const char   *tostr  = " to ";
2442
2443     if (legend == NULL)
2444     {
2445         gmx_fatal(FARGS, "There is no legend in file '%s', can not deduce lambda", fn);
2446     }
2447
2448     /* look for the last 'to': */
2449     ptr2 = legend;
2450     do
2451     {
2452         ptr2 = strstr(ptr2, tostr);
2453         if (ptr2 != NULL)
2454         {
2455             ptr = ptr2;
2456             ptr2++;
2457         }
2458     }
2459     while (ptr2 != NULL && *ptr2 != '\0');
2460
2461     if (ptr)
2462     {
2463         ptr += strlen(tostr)-1; /* and advance past that 'to' */
2464     }
2465     else
2466     {
2467         /* look for the = sign */
2468         ptr = strrchr(legend, '=');
2469         if (!ptr)
2470         {
2471             /* otherwise look for the last space */
2472             ptr = strrchr(legend, ' ');
2473         }
2474     }
2475
2476     if (strstr(legend, "dH"))
2477     {
2478         ok    = TRUE;
2479         bdhdl = TRUE;
2480     }
2481     else if (strchr(legend, 'D') != NULL && strchr(legend, 'H') != NULL)
2482     {
2483         ok    = TRUE;
2484         bdhdl = FALSE;
2485     }
2486     else /*if (strstr(legend, "pV"))*/
2487     {
2488         return FALSE;
2489     }
2490     if (!ptr)
2491     {
2492         ok = FALSE;
2493     }
2494
2495     if (!ok)
2496     {
2497         gmx_fatal(FARGS, "There is no proper lambda legend in file '%s', can not deduce lambda", fn);
2498     }
2499     if (!bdhdl)
2500     {
2501         ptr = find_value(ptr);
2502         if (!ptr || !read_lambda_vector(ptr, lam, NULL, fn))
2503         {
2504             gmx_fatal(FARGS, "lambda vector '%s' %s faulty", legend, fn);
2505         }
2506     }
2507     else
2508     {
2509         int         dhdl_index;
2510         const char *end;
2511         char        buf[STRLEN];
2512
2513         ptr = strrchr(legend, '=');
2514         end = ptr;
2515         if (ptr)
2516         {
2517             /* there must be a component name */
2518             ptr--;
2519             if (ptr < legend)
2520             {
2521                 gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2522             }
2523             /* now backtrack to the start of the identifier */
2524             while (isspace(*ptr))
2525             {
2526                 end = ptr;
2527                 ptr--;
2528                 if (ptr < legend)
2529                 {
2530                     gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2531                 }
2532             }
2533             while (!isspace(*ptr))
2534             {
2535                 ptr--;
2536                 if (ptr < legend)
2537                 {
2538                     gmx_fatal(FARGS, "dhdl legend '%s' %s faulty", legend, fn);
2539                 }
2540             }
2541             ptr++;
2542             strncpy(buf, ptr, (end-ptr));
2543             buf[(end-ptr)] = '\0';
2544             dhdl_index     = lambda_components_find(lam->lc, ptr, (end-ptr));
2545             if (dhdl_index < 0)
2546             {
2547                 char buf[STRLEN];
2548                 strncpy(buf, ptr, (end-ptr));
2549                 buf[(end-ptr)] = '\0';
2550                 gmx_fatal(FARGS,
2551                           "Did not find lambda component for '%s' in %s",
2552                           buf, fn);
2553             }
2554         }
2555         else
2556         {
2557             if (lam->lc->N > 1)
2558             {
2559                 gmx_fatal(FARGS,
2560                           "dhdl without component name with >1 lambda component in %s",
2561                           fn);
2562             }
2563             dhdl_index = 0;
2564         }
2565         lam->dhdl = dhdl_index;
2566     }
2567     return TRUE;
2568 }
2569
2570 static gmx_bool subtitle2lambda(const char *subtitle, xvg_t *ba, const char *fn,
2571                                 lambda_components_t *lc)
2572 {
2573     gmx_bool    bFound;
2574     const char *ptr;
2575     char       *end;
2576     double      native_lambda;
2577
2578     bFound = FALSE;
2579
2580     /* first check for a state string */
2581     ptr = strstr(subtitle, "state");
2582     if (ptr)
2583     {
2584         int         index = -1;
2585         const char *val_end;
2586
2587         /* the new 4.6 style lambda vectors */
2588         ptr = find_value(ptr);
2589         if (ptr)
2590         {
2591             index = strtol(ptr, &end, 10);
2592             if (ptr == end)
2593             {
2594                 gmx_fatal(FARGS, "Incomplete state data in %s", fn);
2595                 return FALSE;
2596             }
2597             ptr = end;
2598         }
2599         else
2600         {
2601             gmx_fatal(FARGS, "Incomplete state data in %s", fn);
2602             return FALSE;
2603         }
2604         /* now find the lambda vector component names */
2605         while (*ptr != '(' && !isalnum(*ptr))
2606         {
2607             ptr++;
2608             if (*ptr == '\0')
2609             {
2610                 gmx_fatal(FARGS,
2611                           "Incomplete lambda vector component data in %s", fn);
2612                 return FALSE;
2613             }
2614         }
2615         val_end = ptr;
2616         if (!read_lambda_components(ptr, lc, &val_end, fn))
2617         {
2618             gmx_fatal(FARGS,
2619                       "lambda vector components in %s don't match those previously read",
2620                       fn);
2621         }
2622         ptr = find_value(val_end);
2623         if (!ptr)
2624         {
2625             gmx_fatal(FARGS, "Incomplete state data in %s", fn);
2626             return FALSE;
2627         }
2628         lambda_vec_init(&(ba->native_lambda), lc);
2629         if (!read_lambda_vector(ptr, &(ba->native_lambda), NULL, fn))
2630         {
2631             gmx_fatal(FARGS, "lambda vector in %s faulty", fn);
2632         }
2633         ba->native_lambda.index = index;
2634         bFound                  = TRUE;
2635     }
2636     else
2637     {
2638         /* compatibility mode: check for lambda in other ways. */
2639         /* plain text lambda string */
2640         ptr = strstr(subtitle, "lambda");
2641         if (ptr == NULL)
2642         {
2643             /* xmgrace formatted lambda string */
2644             ptr = strstr(subtitle, "\\xl\\f{}");
2645         }
2646         if (ptr == NULL)
2647         {
2648             /* xmgr formatted lambda string */
2649             ptr = strstr(subtitle, "\\8l\\4");
2650         }
2651         if (ptr != NULL)
2652         {
2653             ptr = strstr(ptr, "=");
2654         }
2655         if (ptr != NULL)
2656         {
2657             bFound = (sscanf(ptr+1, "%lf", &(native_lambda)) == 1);
2658             /* add the lambda component name as an empty string */
2659             if (lc->N > 0)
2660             {
2661                 if (!lambda_components_check(lc, 0, "", 0))
2662                 {
2663                     gmx_fatal(FARGS,
2664                               "lambda vector components in %s don't match those previously read",
2665                               fn);
2666                 }
2667             }
2668             else
2669             {
2670                 lambda_components_add(lc, "", 0);
2671             }
2672             lambda_vec_init(&(ba->native_lambda), lc);
2673             ba->native_lambda.val[0] = native_lambda;
2674         }
2675     }
2676
2677     return bFound;
2678 }
2679
2680 static void filename2lambda(const char *fn)
2681 {
2682     double        lambda;
2683     const char   *ptr, *digitptr;
2684     char         *endptr;
2685     int           dirsep;
2686     ptr = fn;
2687     /* go to the end of the path string and search backward to find the last
2688        directory in the path which has to contain the value of lambda
2689      */
2690     while (ptr[1] != '\0')
2691     {
2692         ptr++;
2693     }
2694     /* searching backward to find the second directory separator */
2695     dirsep   = 0;
2696     digitptr = NULL;
2697     while (ptr >= fn)
2698     {
2699         if (ptr[0] != DIR_SEPARATOR && ptr[1] == DIR_SEPARATOR)
2700         {
2701             if (dirsep == 1)
2702             {
2703                 break;
2704             }
2705             dirsep++;
2706         }
2707         /* save the last position of a digit between the last two
2708            separators = in the last dirname */
2709         if (dirsep > 0 && isdigit(*ptr))
2710         {
2711             digitptr = ptr;
2712         }
2713         ptr--;
2714     }
2715     if (!digitptr)
2716     {
2717         gmx_fatal(FARGS, "While trying to read the lambda value from the file path:"
2718                   " last directory in the path '%s' does not contain a number", fn);
2719     }
2720     if (digitptr[-1] == '-')
2721     {
2722         digitptr--;
2723     }
2724     lambda = strtod(digitptr, &endptr);
2725     if (endptr == digitptr)
2726     {
2727         gmx_fatal(FARGS, "Malformed number in file path '%s'", fn);
2728     }
2729 }
2730
2731 static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
2732                                   lambda_components_t *lc)
2733 {
2734     int          i;
2735     char        *subtitle, **legend, *ptr;
2736     int          np;
2737     gmx_bool     native_lambda_read = FALSE;
2738     char         buf[STRLEN];
2739     lambda_vec_t lv;
2740
2741     xvg_init(ba);
2742
2743     ba->filename = fn;
2744
2745     np = read_xvg_legend(fn, &ba->y, &ba->nset, &subtitle, &legend);
2746     if (!ba->y)
2747     {
2748         gmx_fatal(FARGS, "File %s contains no usable data.", fn);
2749     }
2750     /* Reorder the data */
2751     ba->t  = ba->y[0];
2752     for (i = 1; i < ba->nset; i++)
2753     {
2754         ba->y[i-1] = ba->y[i];
2755     }
2756     ba->nset--;
2757
2758     snew(ba->np, ba->nset);
2759     for (i = 0; i < ba->nset; i++)
2760     {
2761         ba->np[i] = np;
2762     }
2763
2764     ba->temp = -1;
2765     if (subtitle != NULL)
2766     {
2767         /* try to extract temperature */
2768         ptr = strstr(subtitle, "T =");
2769         if (ptr != NULL)
2770         {
2771             ptr += 3;
2772             if (sscanf(ptr, "%lf", &ba->temp) == 1)
2773             {
2774                 if (ba->temp <= 0)
2775                 {
2776                     gmx_fatal(FARGS, "Found temperature of %f in file '%s'",
2777                               ba->temp, fn);
2778                 }
2779             }
2780         }
2781     }
2782     if (ba->temp < 0)
2783     {
2784         if (*temp <= 0)
2785         {
2786             gmx_fatal(FARGS, "Did not find a temperature in the subtitle in file '%s', use the -temp option of [TT]gmx bar[tt]", fn);
2787         }
2788         ba->temp = *temp;
2789     }
2790
2791     /* Try to deduce lambda from the subtitle */
2792     if (subtitle)
2793     {
2794         if (subtitle2lambda(subtitle, ba, fn, lc))
2795         {
2796             native_lambda_read = TRUE;
2797         }
2798     }
2799     snew(ba->lambda, ba->nset);
2800     if (legend == NULL)
2801     {
2802         /* Check if we have a single set, no legend, nset=1 means t and dH/dl */
2803         if (ba->nset == 1)
2804         {
2805             if (!native_lambda_read)
2806             {
2807                 /* Deduce lambda from the file name */
2808                 filename2lambda(fn);
2809                 native_lambda_read = TRUE;
2810             }
2811             ba->lambda[0] = ba->native_lambda;
2812         }
2813         else
2814         {
2815             gmx_fatal(FARGS, "File %s contains multiple sets but no legends, can not determine the lambda values", fn);
2816         }
2817     }
2818     else
2819     {
2820         for (i = 0; i < ba->nset; )
2821         {
2822             gmx_bool use = FALSE;
2823             /* Read lambda from the legend */
2824             lambda_vec_init( &(ba->lambda[i]), lc );
2825             lambda_vec_copy( &(ba->lambda[i]), &(ba->native_lambda));
2826             use = legend2lambda(fn, legend[i], &(ba->lambda[i]));
2827             if (use)
2828             {
2829                 lambda_vec_print(&(ba->lambda[i]), buf, FALSE);
2830                 i++;
2831             }
2832             else
2833             {
2834                 int j;
2835                 printf("%s: Ignoring set '%s'.\n", fn, legend[i]);
2836                 for (j = i+1; j < ba->nset; j++)
2837                 {
2838                     ba->y[j-1]  = ba->y[j];
2839                     legend[j-1] = legend[j];
2840                 }
2841                 ba->nset--;
2842             }
2843         }
2844     }
2845
2846     if (!native_lambda_read)
2847     {
2848         gmx_fatal(FARGS, "File %s contains multiple sets but no indication of the native lambda", fn);
2849     }
2850
2851     if (legend != NULL)
2852     {
2853         for (i = 0; i < ba->nset-1; i++)
2854         {
2855             sfree(legend[i]);
2856         }
2857         sfree(legend);
2858     }
2859 }
2860
2861 static void read_bar_xvg(char *fn, real *temp, sim_data_t *sd)
2862 {
2863     xvg_t     *barsim;
2864     samples_t *s;
2865     int        i;
2866     double    *lambda;
2867
2868     snew(barsim, 1);
2869
2870     read_bar_xvg_lowlevel(fn, temp, barsim, &(sd->lc));
2871
2872     if (barsim->nset < 1)
2873     {
2874         gmx_fatal(FARGS, "File '%s' contains fewer than two columns", fn);
2875     }
2876
2877     if (!gmx_within_tol(*temp, barsim->temp, GMX_FLOAT_EPS) && (*temp > 0) )
2878     {
2879         gmx_fatal(FARGS, "Temperature in file %s different from earlier files or setting\n", fn);
2880     }
2881     *temp = barsim->temp;
2882
2883     /* now create a series of samples_t */
2884     snew(s, barsim->nset);
2885     for (i = 0; i < barsim->nset; i++)
2886     {
2887         samples_init(s+i, &(barsim->native_lambda), &(barsim->lambda[i]),
2888                      barsim->temp, lambda_vec_same(&(barsim->native_lambda),
2889                                                    &(barsim->lambda[i])),
2890                      fn);
2891         s[i].du  = barsim->y[i];
2892         s[i].ndu = barsim->np[i];
2893         s[i].t   = barsim->t;
2894
2895         lambda_data_list_insert_sample(sd->lb, s+i);
2896     }
2897     {
2898         char buf[STRLEN];
2899
2900         lambda_vec_print(s[0].native_lambda, buf, FALSE);
2901         printf("%s: %.1f - %.1f; lambda = %s\n    dH/dl & foreign lambdas:\n",
2902                fn, s[0].t[0], s[0].t[s[0].ndu-1], buf);
2903         for (i = 0; i < barsim->nset; i++)
2904         {
2905             lambda_vec_print(s[i].foreign_lambda, buf, TRUE);
2906             printf("        %s (%d pts)\n", buf, s[i].ndu);
2907         }
2908     }
2909     printf("\n\n");
2910 }
2911
2912 static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
2913                                  double start_time, double delta_time,
2914                                  lambda_vec_t *native_lambda, double temp,
2915                                  double *last_t, const char *filename)
2916 {
2917     int           i, j;
2918     gmx_bool      allocated;
2919     double        old_foreign_lambda;
2920     lambda_vec_t *foreign_lambda;
2921     int           type;
2922     samples_t    *s; /* convenience pointer */
2923     int           startj;
2924
2925     /* check the block types etc. */
2926     if ( (blk->nsub < 3) ||
2927          (blk->sub[0].type != xdr_datatype_int) ||
2928          (blk->sub[1].type != xdr_datatype_double) ||
2929          (
2930              (blk->sub[2].type != xdr_datatype_float) &&
2931              (blk->sub[2].type != xdr_datatype_double)
2932          ) ||
2933          (blk->sub[0].nr < 1) ||
2934          (blk->sub[1].nr < 1) )
2935     {
2936         gmx_fatal(FARGS,
2937                   "Unexpected/corrupted block data in file %s around time %f.",
2938                   filename, start_time);
2939     }
2940
2941     snew(foreign_lambda, 1);
2942     lambda_vec_init(foreign_lambda, native_lambda->lc);
2943     lambda_vec_copy(foreign_lambda, native_lambda);
2944     type = blk->sub[0].ival[0];
2945     if (type == dhbtDH)
2946     {
2947         for (i = 0; i < native_lambda->lc->N; i++)
2948         {
2949             foreign_lambda->val[i] = blk->sub[1].dval[i];
2950         }
2951     }
2952     else
2953     {
2954         if (blk->sub[0].nr > 1)
2955         {
2956             foreign_lambda->dhdl = blk->sub[0].ival[1];
2957         }
2958         else
2959         {
2960             foreign_lambda->dhdl = 0;
2961         }
2962     }
2963
2964     if (!*smp)
2965     {
2966         /* initialize the samples structure if it's empty. */
2967         snew(*smp, 1);
2968         samples_init(*smp, native_lambda, foreign_lambda, temp,
2969                      type == dhbtDHDL, filename);
2970         (*smp)->start_time = start_time;
2971         (*smp)->delta_time = delta_time;
2972     }
2973
2974     /* set convenience pointer */
2975     s = *smp;
2976
2977     /* now double check */
2978     if (!lambda_vec_same(s->foreign_lambda, foreign_lambda) )
2979     {
2980         char buf[STRLEN], buf2[STRLEN];
2981         lambda_vec_print(foreign_lambda, buf, FALSE);
2982         lambda_vec_print(s->foreign_lambda, buf2, FALSE);
2983         fprintf(stderr, "Got foreign lambda=%s, expected: %s\n", buf, buf2);
2984         gmx_fatal(FARGS, "Corrupted data in file %s around t=%f.",
2985                   filename, start_time);
2986     }
2987
2988     /* make room for the data */
2989     if (s->ndu_alloc < (size_t)(s->ndu + blk->sub[2].nr) )
2990     {
2991         s->ndu_alloc += (s->ndu_alloc < (size_t)blk->sub[2].nr) ?
2992             blk->sub[2].nr*2 : s->ndu_alloc;
2993         srenew(s->du_alloc, s->ndu_alloc);
2994         s->du = s->du_alloc;
2995     }
2996     startj   = s->ndu;
2997     s->ndu  += blk->sub[2].nr;
2998     s->ntot += blk->sub[2].nr;
2999     *ndu     = blk->sub[2].nr;
3000
3001     /* and copy the data*/
3002     for (j = 0; j < blk->sub[2].nr; j++)
3003     {
3004         if (blk->sub[2].type == xdr_datatype_float)
3005         {
3006             s->du[startj+j] = blk->sub[2].fval[j];
3007         }
3008         else
3009         {
3010             s->du[startj+j] = blk->sub[2].dval[j];
3011         }
3012     }
3013     if (start_time + blk->sub[2].nr*delta_time > *last_t)
3014     {
3015         *last_t = start_time + blk->sub[2].nr*delta_time;
3016     }
3017 }
3018
3019 static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
3020                                       double start_time, double delta_time,
3021                                       lambda_vec_t *native_lambda, double temp,
3022                                       double *last_t, const char *filename)
3023 {
3024     int           i, j;
3025     samples_t    *s;
3026     int           nhist;
3027     double        old_foreign_lambda;
3028     lambda_vec_t *foreign_lambda;
3029     int           type;
3030     int           nbins[2];
3031
3032     /* check the block types etc. */
3033     if ( (blk->nsub < 2) ||
3034          (blk->sub[0].type != xdr_datatype_double) ||
3035          (blk->sub[1].type != xdr_datatype_int64) ||
3036          (blk->sub[0].nr < 2)  ||
3037          (blk->sub[1].nr < 2) )
3038     {
3039         gmx_fatal(FARGS,
3040                   "Unexpected/corrupted block data in file %s around time %f",
3041                   filename, start_time);
3042     }
3043
3044     nhist = blk->nsub-2;
3045     if (nhist == 0)
3046     {
3047         return NULL;
3048     }
3049     if (nhist > 2)
3050     {
3051         gmx_fatal(FARGS,
3052                   "Unexpected/corrupted block data in file %s around time %f",
3053                   filename, start_time);
3054     }
3055
3056     snew(s, 1);
3057     *nsamples = 1;
3058
3059     snew(foreign_lambda, 1);
3060     lambda_vec_init(foreign_lambda, native_lambda->lc);
3061     lambda_vec_copy(foreign_lambda, native_lambda);
3062     type = (int)(blk->sub[1].lval[1]);
3063     if (type == dhbtDH)
3064     {
3065         double old_foreign_lambda;
3066
3067         old_foreign_lambda = blk->sub[0].dval[0];
3068         if (old_foreign_lambda >= 0)
3069         {
3070             foreign_lambda->val[0] = old_foreign_lambda;
3071             if (foreign_lambda->lc->N > 1)
3072             {
3073                 gmx_fatal(FARGS,
3074                           "Single-component lambda in multi-component file %s",
3075                           filename);
3076             }
3077         }
3078         else
3079         {
3080             for (i = 0; i < native_lambda->lc->N; i++)
3081             {
3082                 foreign_lambda->val[i] = blk->sub[0].dval[i+2];
3083             }
3084         }
3085     }
3086     else
3087     {
3088         if (foreign_lambda->lc->N > 1)
3089         {
3090             if (blk->sub[1].nr < 3 + nhist)
3091             {
3092                 gmx_fatal(FARGS,
3093                           "Missing derivative coord in multi-component file %s",
3094                           filename);
3095             }
3096             foreign_lambda->dhdl = blk->sub[1].lval[2 + nhist];
3097         }
3098         else
3099         {
3100             foreign_lambda->dhdl = 0;
3101         }
3102     }
3103
3104     samples_init(s, native_lambda, foreign_lambda, temp, type == dhbtDHDL,
3105                  filename);
3106     snew(s->hist, 1);
3107
3108     for (i = 0; i < nhist; i++)
3109     {
3110         nbins[i] = blk->sub[i+2].nr;
3111     }
3112
3113     hist_init(s->hist, nhist, nbins);
3114
3115     for (i = 0; i < nhist; i++)
3116     {
3117         s->hist->x0[i] = blk->sub[1].lval[2+i];
3118         s->hist->dx[i] = blk->sub[0].dval[1];
3119         if (i == 1)
3120         {
3121             s->hist->dx[i] = -s->hist->dx[i];
3122         }
3123     }
3124
3125     s->hist->start_time = start_time;
3126     s->hist->delta_time = delta_time;
3127     s->start_time       = start_time;
3128     s->delta_time       = delta_time;
3129
3130     for (i = 0; i < nhist; i++)
3131     {
3132         int             nbin;
3133         gmx_int64_t     sum = 0;
3134
3135         for (j = 0; j < s->hist->nbin[i]; j++)
3136         {
3137             int binv = (int)(blk->sub[i+2].ival[j]);
3138
3139             s->hist->bin[i][j] = binv;
3140             sum               += binv;
3141
3142         }
3143         if (i == 0)
3144         {
3145             s->ntot      = sum;
3146             s->hist->sum = sum;
3147         }
3148         else
3149         {
3150             if (s->ntot != sum)
3151             {
3152                 gmx_fatal(FARGS, "Histogram counts don't match in %s",
3153                           filename);
3154             }
3155         }
3156     }
3157
3158     if (start_time + s->hist->sum*delta_time > *last_t)
3159     {
3160         *last_t = start_time + s->hist->sum*delta_time;
3161     }
3162     return s;
3163 }
3164
3165
3166 static void read_barsim_edr(char *fn, real *temp, sim_data_t *sd)
3167 {
3168     int            i, j;
3169     ener_file_t    fp;
3170     t_enxframe    *fr;
3171     int            nre;
3172     gmx_enxnm_t   *enm           = NULL;
3173     double         first_t       = -1;
3174     double         last_t        = -1;
3175     samples_t    **samples_rawdh = NULL; /* contains samples for raw delta_h  */
3176     int           *nhists        = NULL; /* array to keep count & print at end */
3177     int           *npts          = NULL; /* array to keep count & print at end */
3178     lambda_vec_t **lambdas       = NULL; /* array to keep count & print at end */
3179     lambda_vec_t  *native_lambda;
3180     double         end_time;             /* the end time of the last batch of samples */
3181     int            nsamples = 0;
3182     lambda_vec_t   start_lambda;
3183
3184     fp = open_enx(fn, "r");
3185     do_enxnms(fp, &nre, &enm);
3186     snew(fr, 1);
3187
3188     snew(native_lambda, 1);
3189     start_lambda.lc = NULL;
3190
3191     while (do_enx(fp, fr))
3192     {
3193         /* count the data blocks */
3194         int    nblocks_raw  = 0;
3195         int    nblocks_hist = 0;
3196         int    nlam         = 0;
3197         int    k;
3198         /* DHCOLL block information: */
3199         double start_time = 0, delta_time = 0, old_start_lambda = 0, delta_lambda = 0;
3200         double rtemp      = 0;
3201
3202         /* count the blocks and handle collection information: */
3203         for (i = 0; i < fr->nblock; i++)
3204         {
3205             if (fr->block[i].id == enxDHHIST)
3206             {
3207                 nblocks_hist++;
3208             }
3209             if (fr->block[i].id == enxDH)
3210             {
3211                 nblocks_raw++;
3212             }
3213             if (fr->block[i].id == enxDHCOLL)
3214             {
3215                 nlam++;
3216                 if ( (fr->block[i].nsub < 1) ||
3217                      (fr->block[i].sub[0].type != xdr_datatype_double) ||
3218                      (fr->block[i].sub[0].nr < 5))
3219                 {
3220                     gmx_fatal(FARGS, "Unexpected block data in file %s", fn);
3221                 }
3222
3223                 /* read the data from the DHCOLL block */
3224                 rtemp            =        fr->block[i].sub[0].dval[0];
3225                 start_time       =   fr->block[i].sub[0].dval[1];
3226                 delta_time       =   fr->block[i].sub[0].dval[2];
3227                 old_start_lambda = fr->block[i].sub[0].dval[3];
3228                 delta_lambda     = fr->block[i].sub[0].dval[4];
3229
3230                 if (delta_lambda > 0)
3231                 {
3232                     gmx_fatal(FARGS, "Lambda values not constant in %s: can't apply BAR method", fn);
3233                 }
3234                 if ( ( *temp != rtemp) && (*temp > 0) )
3235                 {
3236                     gmx_fatal(FARGS, "Temperature in file %s different from earlier files or setting\n", fn);
3237                 }
3238                 *temp = rtemp;
3239
3240                 if (old_start_lambda >= 0)
3241                 {
3242                     if (sd->lc.N > 0)
3243                     {
3244                         if (!lambda_components_check(&(sd->lc), 0, "", 0))
3245                         {
3246                             gmx_fatal(FARGS,
3247                                       "lambda vector components in %s don't match those previously read",
3248                                       fn);
3249                         }
3250                     }
3251                     else
3252                     {
3253                         lambda_components_add(&(sd->lc), "", 0);
3254                     }
3255                     if (!start_lambda.lc)
3256                     {
3257                         lambda_vec_init(&start_lambda, &(sd->lc));
3258                     }
3259                     start_lambda.val[0] = old_start_lambda;
3260                 }
3261                 else
3262                 {
3263                     /* read lambda vector */
3264                     int      n_lambda_vec;
3265                     gmx_bool check = (sd->lc.N > 0);
3266                     if (fr->block[i].nsub < 2)
3267                     {
3268                         gmx_fatal(FARGS,
3269                                   "No lambda vector, but start_lambda=%f\n",
3270                                   old_start_lambda);
3271                     }
3272                     n_lambda_vec = fr->block[i].sub[1].ival[1];
3273                     for (j = 0; j < n_lambda_vec; j++)
3274                     {
3275                         const char *name =
3276                             efpt_singular_names[fr->block[i].sub[1].ival[1+j]];
3277                         if (check)
3278                         {
3279                             /* check the components */
3280                             lambda_components_check(&(sd->lc), j, name,
3281                                                     strlen(name));
3282                         }
3283                         else
3284                         {
3285                             lambda_components_add(&(sd->lc), name,
3286                                                   strlen(name));
3287                         }
3288                     }
3289                     lambda_vec_init(&start_lambda, &(sd->lc));
3290                     start_lambda.index = fr->block[i].sub[1].ival[0];
3291                     for (j = 0; j < n_lambda_vec; j++)
3292                     {
3293                         start_lambda.val[j] = fr->block[i].sub[0].dval[5+j];
3294                     }
3295                 }
3296                 if (first_t < 0)
3297                 {
3298                     first_t = start_time;
3299                 }
3300             }
3301         }
3302
3303         if (nlam != 1)
3304         {
3305             gmx_fatal(FARGS, "Did not find delta H information in file %s", fn);
3306         }
3307         if (nblocks_raw > 0 && nblocks_hist > 0)
3308         {
3309             gmx_fatal(FARGS, "Can't handle both raw delta U data and histograms in the same file %s", fn);
3310         }
3311
3312         if (nsamples > 0)
3313         {
3314             /* check the native lambda */
3315             if (!lambda_vec_same(&start_lambda, native_lambda) )
3316             {
3317                 gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %f, and becomes %f at time %f",
3318                           fn, native_lambda, start_lambda, start_time);
3319             }
3320             /* check the number of samples against the previous number */
3321             if ( ((nblocks_raw+nblocks_hist) != nsamples) || (nlam != 1 ) )
3322             {
3323                 gmx_fatal(FARGS, "Unexpected block count in %s: was %d, now %d\n",
3324                           fn, nsamples+1, nblocks_raw+nblocks_hist+nlam);
3325             }
3326             /* check whether last iterations's end time matches with
3327                the currrent start time */
3328             if ( (fabs(last_t - start_time) > 2*delta_time)  && last_t >= 0)
3329             {
3330                 /* it didn't. We need to store our samples and reallocate */
3331                 for (i = 0; i < nsamples; i++)
3332                 {
3333                     if (samples_rawdh[i])
3334                     {
3335                         /* insert it into the existing list */
3336                         lambda_data_list_insert_sample(sd->lb,
3337                                                        samples_rawdh[i]);
3338                         /* and make sure we'll allocate a new one this time
3339                            around */
3340                         samples_rawdh[i] = NULL;
3341                     }
3342                 }
3343             }
3344         }
3345         else
3346         {
3347             /* this is the first round; allocate the associated data
3348                structures */
3349             /*native_lambda=start_lambda;*/
3350             lambda_vec_init(native_lambda, &(sd->lc));
3351             lambda_vec_copy(native_lambda, &start_lambda);
3352             nsamples = nblocks_raw+nblocks_hist;
3353             snew(nhists, nsamples);
3354             snew(npts, nsamples);
3355             snew(lambdas, nsamples);
3356             snew(samples_rawdh, nsamples);
3357             for (i = 0; i < nsamples; i++)
3358             {
3359                 nhists[i]        = 0;
3360                 npts[i]          = 0;
3361                 lambdas[i]       = NULL;
3362                 samples_rawdh[i] = NULL; /* init to NULL so we know which
3363                                             ones contain values */
3364             }
3365         }
3366
3367         /* and read them */
3368         k = 0; /* counter for the lambdas, etc. arrays */
3369         for (i = 0; i < fr->nblock; i++)
3370         {
3371             if (fr->block[i].id == enxDH)
3372             {
3373                 int type = (fr->block[i].sub[0].ival[0]);
3374                 if (type == dhbtDH || type == dhbtDHDL)
3375                 {
3376                     int ndu;
3377                     read_edr_rawdh_block(&(samples_rawdh[k]),
3378                                          &ndu,
3379                                          &(fr->block[i]),
3380                                          start_time, delta_time,
3381                                          native_lambda, rtemp,
3382                                          &last_t, fn);
3383                     npts[k] += ndu;
3384                     if (samples_rawdh[k])
3385                     {
3386                         lambdas[k] = samples_rawdh[k]->foreign_lambda;
3387                     }
3388                     k++;
3389                 }
3390             }
3391             else if (fr->block[i].id == enxDHHIST)
3392             {
3393                 int type = (int)(fr->block[i].sub[1].lval[1]);
3394                 if (type == dhbtDH || type == dhbtDHDL)
3395                 {
3396                     int        j;
3397                     int        nb = 0;
3398                     samples_t *s; /* this is where the data will go */
3399                     s = read_edr_hist_block(&nb, &(fr->block[i]),
3400                                             start_time, delta_time,
3401                                             native_lambda, rtemp,
3402                                             &last_t, fn);
3403                     nhists[k] += nb;
3404                     if (nb > 0)
3405                     {
3406                         lambdas[k] = s->foreign_lambda;
3407                     }
3408                     k++;
3409                     /* and insert the new sample immediately */
3410                     for (j = 0; j < nb; j++)
3411                     {
3412                         lambda_data_list_insert_sample(sd->lb, s+j);
3413                     }
3414                 }
3415             }
3416         }
3417     }
3418     /* Now store all our extant sample collections */
3419     for (i = 0; i < nsamples; i++)
3420     {
3421         if (samples_rawdh[i])
3422         {
3423             /* insert it into the existing list */
3424             lambda_data_list_insert_sample(sd->lb, samples_rawdh[i]);
3425         }
3426     }
3427
3428
3429     {
3430         char buf[STRLEN];
3431         printf("\n");
3432         lambda_vec_print(native_lambda, buf, FALSE);
3433         printf("%s: %.1f - %.1f; lambda = %s\n    foreign lambdas:\n",
3434                fn, first_t, last_t, buf);
3435         for (i = 0; i < nsamples; i++)
3436         {
3437             if (lambdas[i])
3438             {
3439                 lambda_vec_print(lambdas[i], buf, TRUE);
3440                 if (nhists[i] > 0)
3441                 {
3442                     printf("        %s (%d hists)\n", buf, nhists[i]);
3443                 }
3444                 else
3445                 {
3446                     printf("        %s (%d pts)\n", buf, npts[i]);
3447                 }
3448             }
3449         }
3450     }
3451     printf("\n\n");
3452     sfree(npts);
3453     sfree(nhists);
3454     sfree(lambdas);
3455 }
3456
3457
3458 int gmx_bar(int argc, char *argv[])
3459 {
3460     static const char *desc[] = {
3461         "[THISMODULE] calculates free energy difference estimates through ",
3462         "Bennett's acceptance ratio method (BAR). It also automatically",
3463         "adds series of individual free energies obtained with BAR into",
3464         "a combined free energy estimate.[PAR]",
3465
3466         "Every individual BAR free energy difference relies on two ",
3467         "simulations at different states: say state A and state B, as",
3468         "controlled by a parameter, [GRK]lambda[grk] (see the [TT].mdp[tt] parameter",
3469         "[TT]init_lambda[tt]). The BAR method calculates a ratio of weighted",
3470         "average of the Hamiltonian difference of state B given state A and",
3471         "vice versa.",
3472         "The energy differences to the other state must be calculated",
3473         "explicitly during the simulation. This can be done with",
3474         "the [TT].mdp[tt] option [TT]foreign_lambda[tt].[PAR]",
3475
3476         "Input option [TT]-f[tt] expects multiple [TT]dhdl.xvg[tt] files. ",
3477         "Two types of input files are supported:[BR]",
3478         "[TT]*[tt]  Files with more than one [IT]y[it]-value. ",
3479         "The files should have columns ",
3480         "with dH/d[GRK]lambda[grk] and [GRK]Delta[grk][GRK]lambda[grk]. ",
3481         "The [GRK]lambda[grk] values are inferred ",
3482         "from the legends: [GRK]lambda[grk] of the simulation from the legend of ",
3483         "dH/d[GRK]lambda[grk] and the foreign [GRK]lambda[grk] values from the ",
3484         "legends of Delta H",
3485         "[BR]",
3486         "[TT]*[tt]  Files with only one [IT]y[it]-value. Using the",
3487         "[TT]-extp[tt] option for these files, it is assumed",
3488         "that the [IT]y[it]-value is dH/d[GRK]lambda[grk] and that the ",
3489         "Hamiltonian depends linearly on [GRK]lambda[grk]. ",
3490         "The [GRK]lambda[grk] value of the simulation is inferred from the ",
3491         "subtitle (if present), otherwise from a number in the subdirectory ",
3492         "in the file name.[PAR]",
3493
3494         "The [GRK]lambda[grk] of the simulation is parsed from ",
3495         "[TT]dhdl.xvg[tt] file's legend containing the string 'dH', the ",
3496         "foreign [GRK]lambda[grk] values from the legend containing the ",
3497         "capitalized letters 'D' and 'H'. The temperature is parsed from ",
3498         "the legend line containing 'T ='.[PAR]",
3499
3500         "The input option [TT]-g[tt] expects multiple [TT].edr[tt] files. ",
3501         "These can contain either lists of energy differences (see the ",
3502         "[TT].mdp[tt] option [TT]separate_dhdl_file[tt]), or a series of ",
3503         "histograms (see the [TT].mdp[tt] options [TT]dh_hist_size[tt] and ",
3504         "[TT]dh_hist_spacing[tt]).", "The temperature and [GRK]lambda[grk] ",
3505         "values are automatically deduced from the [TT]ener.edr[tt] file.[PAR]",
3506
3507         "In addition to the [TT].mdp[tt] option [TT]foreign_lambda[tt], ",
3508         "the energy difference can also be extrapolated from the ",
3509         "dH/d[GRK]lambda[grk] values. This is done with the[TT]-extp[tt]",
3510         "option, which assumes that the system's Hamiltonian depends linearly",
3511         "on [GRK]lambda[grk], which is not normally the case.[PAR]",
3512
3513         "The free energy estimates are determined using BAR with bisection, ",
3514         "with the precision of the output set with [TT]-prec[tt]. ",
3515         "An error estimate taking into account time correlations ",
3516         "is made by splitting the data into blocks and determining ",
3517         "the free energy differences over those blocks and assuming ",
3518         "the blocks are independent. ",
3519         "The final error estimate is determined from the average variance ",
3520         "over 5 blocks. A range of block numbers for error estimation can ",
3521         "be provided with the options [TT]-nbmin[tt] and [TT]-nbmax[tt].[PAR]",
3522
3523         "[THISMODULE] tries to aggregate samples with the same 'native' and ",
3524         "'foreign' [GRK]lambda[grk] values, but always assumes independent ",
3525         "samples. [BB]Note[bb] that when aggregating energy ",
3526         "differences/derivatives with different sampling intervals, this is ",
3527         "almost certainly not correct. Usually subsequent energies are ",
3528         "correlated and different time intervals mean different degrees ",
3529         "of correlation between samples.[PAR]",
3530
3531         "The results are split in two parts: the last part contains the final ",
3532         "results in kJ/mol, together with the error estimate for each part ",
3533         "and the total. The first part contains detailed free energy ",
3534         "difference estimates and phase space overlap measures in units of ",
3535         "kT (together with their computed error estimate). The printed ",
3536         "values are:[BR]",
3537         "[TT]*[tt]  lam_A: the [GRK]lambda[grk] values for point A.[BR]",
3538         "[TT]*[tt]  lam_B: the [GRK]lambda[grk] values for point B.[BR]",
3539         "[TT]*[tt]     DG: the free energy estimate.[BR]",
3540         "[TT]*[tt]    s_A: an estimate of the relative entropy of B in A.[BR]",
3541         "[TT]*[tt]    s_B: an estimate of the relative entropy of A in B.[BR]",
3542         "[TT]*[tt]  stdev: an estimate expected per-sample standard deviation.[PAR]",
3543
3544         "The relative entropy of both states in each other's ensemble can be ",
3545         "interpreted as a measure of phase space overlap: ",
3546         "the relative entropy s_A of the work samples of lambda_B in the ",
3547         "ensemble of lambda_A (and vice versa for s_B), is a ",
3548         "measure of the 'distance' between Boltzmann distributions of ",
3549         "the two states, that goes to zero for identical distributions. See ",
3550         "Wu & Kofke, J. Chem. Phys. 123 084109 (2005) for more information.",
3551         "[PAR]",
3552         "The estimate of the expected per-sample standard deviation, as given ",
3553         "in Bennett's original BAR paper: Bennett, J. Comp. Phys. 22, p 245 (1976).",
3554         "Eq. 10 therein gives an estimate of the quality of sampling (not directly",
3555         "of the actual statistical error, because it assumes independent samples).[PAR]",
3556
3557         "To get a visual estimate of the phase space overlap, use the ",
3558         "[TT]-oh[tt] option to write series of histograms, together with the ",
3559         "[TT]-nbin[tt] option.[PAR]"
3560     };
3561     static real        begin    = 0, end = -1, temp = -1;
3562     int                nd       = 2, nbmin = 5, nbmax = 5;
3563     int                nbin     = 100;
3564     gmx_bool           use_dhdl = FALSE;
3565     gmx_bool           calc_s, calc_v;
3566     t_pargs            pa[] = {
3567         { "-b",    FALSE, etREAL, {&begin},  "Begin time for BAR" },
3568         { "-e",    FALSE, etREAL, {&end},    "End time for BAR" },
3569         { "-temp", FALSE, etREAL, {&temp},   "Temperature (K)" },
3570         { "-prec", FALSE, etINT,  {&nd},     "The number of digits after the decimal point" },
3571         { "-nbmin",  FALSE, etINT,  {&nbmin}, "Minimum number of blocks for error estimation" },
3572         { "-nbmax",  FALSE, etINT,  {&nbmax}, "Maximum number of blocks for error estimation" },
3573         { "-nbin",  FALSE, etINT, {&nbin}, "Number of bins for histogram output"},
3574         { "-extp",  FALSE, etBOOL, {&use_dhdl}, "Whether to linearly extrapolate dH/dl values to use as energies"}
3575     };
3576
3577     t_filenm           fnm[] = {
3578         { efXVG, "-f",  "dhdl",   ffOPTRDMULT },
3579         { efEDR, "-g",  "ener",   ffOPTRDMULT },
3580         { efXVG, "-o",  "bar",    ffOPTWR },
3581         { efXVG, "-oi", "barint", ffOPTWR },
3582         { efXVG, "-oh", "histogram", ffOPTWR }
3583     };
3584 #define NFILE asize(fnm)
3585
3586     int          f, i, j;
3587     int          nf = 0;    /* file counter */
3588     int          nbs;
3589     int          nfile_tot; /* total number of input files */
3590     int          nxvgfile = 0;
3591     int          nedrfile = 0;
3592     char       **fxvgnms;
3593     char       **fedrnms;
3594     sim_data_t   sim_data; /* the simulation data */
3595     barres_t    *results;  /* the results */
3596     int          nresults; /* number of results in results array */
3597
3598     double      *partsum;
3599     double       prec, dg_tot, dg, sig, dg_tot_max, dg_tot_min;
3600     FILE        *fpb, *fpi;
3601     char         dgformat[20], xvg2format[STRLEN], xvg3format[STRLEN];
3602     char         buf[STRLEN], buf2[STRLEN];
3603     char         ktformat[STRLEN], sktformat[STRLEN];
3604     char         kteformat[STRLEN], skteformat[STRLEN];
3605     output_env_t oenv;
3606     double       kT, beta;
3607     gmx_bool     result_OK = TRUE, bEE = TRUE;
3608
3609     gmx_bool     disc_err          = FALSE;
3610     double       sum_disc_err      = 0.; /* discretization error */
3611     gmx_bool     histrange_err     = FALSE;
3612     double       sum_histrange_err = 0.; /* histogram range error */
3613     double       stat_err          = 0.; /* statistical error */
3614
3615     if (!parse_common_args(&argc, argv,
3616                            PCA_CAN_VIEW,
3617                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
3618     {
3619         return 0;
3620     }
3621
3622     if (opt2bSet("-f", NFILE, fnm))
3623     {
3624         nxvgfile = opt2fns(&fxvgnms, "-f", NFILE, fnm);
3625     }
3626     if (opt2bSet("-g", NFILE, fnm))
3627     {
3628         nedrfile = opt2fns(&fedrnms, "-g", NFILE, fnm);
3629     }
3630
3631     sim_data_init(&sim_data);
3632 #if 0
3633     /* make linked list */
3634     lb = &lambda_head;
3635     lambda_data_init(lb, 0, 0);
3636     lb->next = lb;
3637     lb->prev = lb;
3638 #endif
3639
3640
3641     nfile_tot = nxvgfile + nedrfile;
3642
3643     if (nfile_tot == 0)
3644     {
3645         gmx_fatal(FARGS, "No input files!");
3646     }
3647
3648     if (nd < 0)
3649     {
3650         gmx_fatal(FARGS, "Can not have negative number of digits");
3651     }
3652     prec = pow(10, -nd);
3653
3654     snew(partsum, (nbmax+1)*(nbmax+1));
3655     nf = 0;
3656
3657     /* read in all files. First xvg files */
3658     for (f = 0; f < nxvgfile; f++)
3659     {
3660         read_bar_xvg(fxvgnms[f], &temp, &sim_data);
3661         nf++;
3662     }
3663     /* then .edr files */
3664     for (f = 0; f < nedrfile; f++)
3665     {
3666         read_barsim_edr(fedrnms[f], &temp, &sim_data);;
3667         nf++;
3668     }
3669
3670     /* fix the times to allow for equilibration */
3671     sim_data_impose_times(&sim_data, begin, end);
3672
3673     if (opt2bSet("-oh", NFILE, fnm))
3674     {
3675         sim_data_histogram(&sim_data, opt2fn("-oh", NFILE, fnm), nbin, oenv);
3676     }
3677
3678     /* assemble the output structures from the lambdas */
3679     results = barres_list_create(&sim_data, &nresults, use_dhdl);
3680
3681     sum_disc_err = barres_list_max_disc_err(results, nresults);
3682
3683     if (nresults == 0)
3684     {
3685         printf("\nNo results to calculate.\n");
3686         return 0;
3687     }
3688
3689     if (sum_disc_err > prec)
3690     {
3691         prec = sum_disc_err;
3692         nd   = ceil(-log10(prec));
3693         printf("WARNING: setting the precision to %g because that is the minimum\n         reasonable number, given the expected discretization error.\n", prec);
3694     }
3695
3696
3697     /*sprintf(lamformat,"%%6.3f");*/
3698     sprintf( dgformat, "%%%d.%df", 3+nd, nd);
3699     /* the format strings of the results in kT */
3700     sprintf( ktformat, "%%%d.%df", 5+nd, nd);
3701     sprintf( sktformat, "%%%ds", 6+nd);
3702     /* the format strings of the errors in kT */
3703     sprintf( kteformat, "%%%d.%df", 3+nd, nd);
3704     sprintf( skteformat, "%%%ds", 4+nd);
3705     sprintf(xvg2format, "%s %s\n", "%s", dgformat);
3706     sprintf(xvg3format, "%s %s %s\n", "%s", dgformat, dgformat);
3707
3708
3709
3710     fpb = NULL;
3711     if (opt2bSet("-o", NFILE, fnm))
3712     {
3713         sprintf(buf, "%s (%s)", "\\DeltaG", "kT");
3714         fpb = xvgropen_type(opt2fn("-o", NFILE, fnm), "Free energy differences",
3715                             "\\lambda", buf, exvggtXYDY, oenv);
3716     }
3717
3718     fpi = NULL;
3719     if (opt2bSet("-oi", NFILE, fnm))
3720     {
3721         sprintf(buf, "%s (%s)", "\\DeltaG", "kT");
3722         fpi = xvgropen(opt2fn("-oi", NFILE, fnm), "Free energy integral",
3723                        "\\lambda", buf, oenv);
3724     }
3725
3726
3727
3728     if (nbmin > nbmax)
3729     {
3730         nbmin = nbmax;
3731     }
3732
3733     /* first calculate results */
3734     bEE      = TRUE;
3735     disc_err = FALSE;
3736     for (f = 0; f < nresults; f++)
3737     {
3738         /* Determine the free energy difference with a factor of 10
3739          * more accuracy than requested for printing.
3740          */
3741         calc_bar(&(results[f]), 0.1*prec, nbmin, nbmax,
3742                  &bEE, partsum);
3743
3744         if (results[f].dg_disc_err > prec/10.)
3745         {
3746             disc_err = TRUE;
3747         }
3748         if (results[f].dg_histrange_err > prec/10.)
3749         {
3750             histrange_err = TRUE;
3751         }
3752     }
3753
3754     /* print results in kT */
3755     kT   = BOLTZ*temp;
3756     beta = 1/kT;
3757
3758     printf("\nTemperature: %g K\n", temp);
3759
3760     printf("\nDetailed results in kT (see help for explanation):\n\n");
3761     printf("%6s ", " lam_A");
3762     printf("%6s ", " lam_B");
3763     printf(sktformat,  "DG ");
3764     if (bEE)
3765     {
3766         printf(skteformat, "+/- ");
3767     }
3768     if (disc_err)
3769     {
3770         printf(skteformat, "disc ");
3771     }
3772     if (histrange_err)
3773     {
3774         printf(skteformat, "range ");
3775     }
3776     printf(sktformat,  "s_A ");
3777     if (bEE)
3778     {
3779         printf(skteformat, "+/- " );
3780     }
3781     printf(sktformat,  "s_B ");
3782     if (bEE)
3783     {
3784         printf(skteformat, "+/- " );
3785     }
3786     printf(sktformat,  "stdev ");
3787     if (bEE)
3788     {
3789         printf(skteformat, "+/- ");
3790     }
3791     printf("\n");
3792     for (f = 0; f < nresults; f++)
3793     {
3794         lambda_vec_print_short(results[f].a->native_lambda, buf);
3795         printf("%s ", buf);
3796         lambda_vec_print_short(results[f].b->native_lambda, buf);
3797         printf("%s ", buf);
3798         printf(ktformat,  results[f].dg);
3799         printf(" ");
3800         if (bEE)
3801         {
3802             printf(kteformat, results[f].dg_err);
3803             printf(" ");
3804         }
3805         if (disc_err)
3806         {
3807             printf(kteformat, results[f].dg_disc_err);
3808             printf(" ");
3809         }
3810         if (histrange_err)
3811         {
3812             printf(kteformat, results[f].dg_histrange_err);
3813             printf(" ");
3814         }
3815         printf(ktformat,  results[f].sa);
3816         printf(" ");
3817         if (bEE)
3818         {
3819             printf(kteformat, results[f].sa_err);
3820             printf(" ");
3821         }
3822         printf(ktformat,  results[f].sb);
3823         printf(" ");
3824         if (bEE)
3825         {
3826             printf(kteformat, results[f].sb_err);
3827             printf(" ");
3828         }
3829         printf(ktformat,  results[f].dg_stddev);
3830         printf(" ");
3831         if (bEE)
3832         {
3833             printf(kteformat, results[f].dg_stddev_err);
3834         }
3835         printf("\n");
3836
3837         /* Check for negative relative entropy with a 95% certainty. */
3838         if (results[f].sa < -2*results[f].sa_err ||
3839             results[f].sb < -2*results[f].sb_err)
3840         {
3841             result_OK = FALSE;
3842         }
3843     }
3844
3845     if (!result_OK)
3846     {
3847         printf("\nWARNING: Some of these results violate the Second Law of "
3848                "Thermodynamics: \n"
3849                "         This is can be the result of severe undersampling, or "
3850                "(more likely)\n"
3851                "         there is something wrong with the simulations.\n");
3852     }
3853
3854
3855     /* final results in kJ/mol */
3856     printf("\n\nFinal results in kJ/mol:\n\n");
3857     dg_tot  = 0;
3858     for (f = 0; f < nresults; f++)
3859     {
3860
3861         if (fpi != NULL)
3862         {
3863             lambda_vec_print_short(results[f].a->native_lambda, buf);
3864             fprintf(fpi, xvg2format, buf, dg_tot);
3865         }
3866
3867
3868         if (fpb != NULL)
3869         {
3870             lambda_vec_print_intermediate(results[f].a->native_lambda,
3871                                           results[f].b->native_lambda,
3872                                           buf);
3873
3874             fprintf(fpb, xvg3format, buf, results[f].dg, results[f].dg_err);
3875         }
3876
3877         printf("point ");
3878         lambda_vec_print_short(results[f].a->native_lambda, buf);
3879         lambda_vec_print_short(results[f].b->native_lambda, buf2);
3880         printf("%s - %s", buf, buf2);
3881         printf(",   DG ");
3882
3883         printf(dgformat, results[f].dg*kT);
3884         if (bEE)
3885         {
3886             printf(" +/- ");
3887             printf(dgformat, results[f].dg_err*kT);
3888         }
3889         if (histrange_err)
3890         {
3891             printf(" (max. range err. = ");
3892             printf(dgformat, results[f].dg_histrange_err*kT);
3893             printf(")");
3894             sum_histrange_err += results[f].dg_histrange_err*kT;
3895         }
3896
3897         printf("\n");
3898         dg_tot += results[f].dg;
3899     }
3900     printf("\n");
3901     printf("total ");
3902     lambda_vec_print_short(results[0].a->native_lambda, buf);
3903     lambda_vec_print_short(results[nresults-1].b->native_lambda, buf2);
3904     printf("%s - %s", buf, buf2);
3905     printf(",   DG ");
3906
3907     printf(dgformat, dg_tot*kT);
3908     if (bEE)
3909     {
3910         stat_err = bar_err(nbmin, nbmax, partsum)*kT;
3911         printf(" +/- ");
3912         printf(dgformat, max(max(stat_err, sum_disc_err), sum_histrange_err));
3913     }
3914     printf("\n");
3915     if (disc_err)
3916     {
3917         printf("\nmaximum discretization error = ");
3918         printf(dgformat, sum_disc_err);
3919         if (bEE && stat_err < sum_disc_err)
3920         {
3921             printf("WARNING: discretization error (%g) is larger than statistical error.\n       Decrease histogram spacing for more accurate results\n", stat_err);
3922         }
3923     }
3924     if (histrange_err)
3925     {
3926         printf("\nmaximum histogram range error = ");
3927         printf(dgformat, sum_histrange_err);
3928         if (bEE && stat_err < sum_histrange_err)
3929         {
3930             printf("WARNING: histogram range error (%g) is larger than statistical error.\n       Increase histogram range for more accurate results\n", stat_err);
3931         }
3932
3933     }
3934     printf("\n");
3935
3936
3937     if (fpi != NULL)
3938     {
3939         lambda_vec_print_short(results[nresults-1].b->native_lambda, buf);
3940         fprintf(fpi, xvg2format, buf, dg_tot);
3941         gmx_ffclose(fpi);
3942     }
3943
3944     do_view(oenv, opt2fn_null("-o", NFILE, fnm), "-xydy");
3945     do_view(oenv, opt2fn_null("-oi", NFILE, fnm), "-xydy");
3946
3947     return 0;
3948 }