src/gromacs/gmxana/gmx_bar.cpp

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