src/gromacs/gmxana/gmx_wham.cpp

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*- */
   2 /*
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Green Red Orange Magenta Azure Cyan Skyblue
  35  */
  36
  37 /*! \internal \file
  38  *  \brief Implementation of the Weighted Histogram Analysis Method (WHAM)
  39  *
  40  *  \author Jochen Hub <jhub@gwdg.de>
  41  */
  42
  43 #ifdef HAVE_CONFIG_H
  44 #include <config.h>
  45 #endif
  46 #include <stdio.h>
  47 #include <stdlib.h>
  48 #include <sstream>
  49
  50 #include "gromacs/commandline/pargs.h"
  51 #include "typedefs.h"
  52 #include "smalloc.h"
  53 #include "vec.h"
  54 #include "copyrite.h"
  55 #include "gromacs/fileio/tpxio.h"
  56 #include "names.h"
  57 #include "gmx_random.h"
  58 #include "gmx_ana.h"
  59 #include "macros.h"
  60
  61 #include "string2.h"
  62 #include "xvgr.h"
  63
  64 //! longest file names allowed in input files
  65 #define WHAM_MAXFILELEN 2048
  66
  67 /*! \brief
  68  * enum for energy units
  69  */
  70 enum {
  71     enSel, en_kJ, en_kCal, en_kT, enNr
  72 };
  73 /*! \brief
  74  * enum for type of input files (pdos, tpr, or pullf)
  75  */
  76 enum {
  77     whamin_unknown, whamin_tpr, whamin_pullxf, whamin_pdo
  78 };
  79
  80 /*! \brief enum for bootstrapping method
  81  *
  82  * These bootstrap methods are supported:
  83  *  - bootstrap complete histograms with continuous weights (Bayesian bootstrap)
  84  *    (bsMethod_BayesianHist)
  85  *  - bootstrap complete histograms (bsMethod_hist)
  86  *  - bootstrap trajectories from given umbrella histograms. This generates new
  87  *    "synthetic" histograms (bsMethod_traj)
  88  *  - bootstrap trajectories from Gaussian with mu/sigma computed from
  89  *    the respective histogram (bsMethod_trajGauss). This gives very similar
  90  *    results compared to bsMethod_traj.
  91  *
  92  *  ********************************************************************
  93  *  FOR MORE DETAILS ON THE BOOTSTRAP METHODS (INCLUDING EXAMPLES), SEE
  94  *  JS Hub, BL de Groot, D van der Spoel
  95  *  g_wham - A free weighted histogram analysis implementation including
  96  *  robust error and autocorrelation estimates,
  97  *  J Chem Theory Comput, 6(12), 3713-3720 (2010)
  98  *  ********************************************************************
  99  */
 100 enum {
 101     bsMethod_unknown, bsMethod_BayesianHist, bsMethod_hist,
 102     bsMethod_traj, bsMethod_trajGauss
 103 };
 104
 105
 106 //! Parameters of the umbrella potentials
 107 typedef struct
 108 {
 109     /*!
 110      * \name Using umbrella pull code since gromacs 4.x
 111      */
 112     /*!\{*/
 113     int   npullcrds;     //!< nr of pull coordinates in tpr file
 114     int   pull_geometry; //!< such as distance, direction
 115     ivec  pull_dim;      //!< pull dimension with geometry distance
 116     int   pull_ndim;     //!< nr of pull_dim != 0
 117     gmx_bool bPrintRef;  //!< Coordinates of reference groups written to pullx.xvg ?
 118     real *k;             //!< force constants in tpr file
 119     real *init_dist;     //!< reference displacements
 120     real *umbInitDist;   //!< reference displacement in umbrella direction
 121     /*!\}*/
 122     /*!
 123      * \name Using PDO files common until gromacs 3.x
 124      */
 125     /*!\{*/
 126     int    nSkip;
 127     char   Reference[256];
 128     int    nPull;
 129     int    nDim;
 130     ivec   Dims;
 131     char   PullName[4][256];
 132     double UmbPos[4][3];
 133     double UmbCons[4][3];
 134     /*!\}*/
 135 } t_UmbrellaHeader;
 136
 137 //! Data in the umbrella histograms
 138 typedef struct
 139 {
 140     int      nPull;       //!< nr of pull groups in this pdo or pullf/x file
 141     double **Histo;       //!< nPull histograms
 142     double **cum;         //!< nPull cumulative distribution functions
 143     int      nBin;        //!< nr of bins. identical to opt->bins
 144     double  *k;           //!< force constants for the nPull groups
 145     double  *pos;         //!< umbrella positions for the nPull groups
 146     double  *z;           //!< z=(-Fi/kT) for the nPull groups. These values are iteratively computed during wham
 147     int     *N;           //!< nr of data points in nPull histograms.
 148     int     *Ntot;        //!< also nr of data points. N and Ntot only differ if bHistEq==TRUE
 149
 150     /*! \brief  g = 1 + 2*tau[int]/dt where tau is the integrated autocorrelation time.
 151      *
 152      * Compare, e.g. Ferrenberg/Swendsen, PRL 63:1195 (1989),
 153      * Kumar et al, J Comp Chem 13, 1011-1021 (1992), eq. 28
 154      */
 155     double *g;
 156     double *tau;         //!< intetrated autocorrelation time (IACT)
 157     double *tausmooth;   //!< smoothed IACT
 158
 159     double  dt;          //!< timestep in the input data. Can be adapted with g_wham option -dt
 160
 161     /*! \brief TRUE, if any data point of the histogram is within min and max, otherwise FALSE */
 162     gmx_bool **bContrib;
 163     real     **ztime;     //!< input data z(t) as a function of time. Required to compute ACTs
 164
 165     /*! \brief average force estimated from average displacement, fAv=dzAv*k
 166      *
 167      *  Used for integration to guess the potential.
 168      */
 169     real   *forceAv;
 170     real   *aver;         //!< average of histograms
 171     real   *sigma;        //!< stddev of histograms
 172     double *bsWeight;     //!< for bootstrapping complete histograms with continuous weights
 173 } t_UmbrellaWindow;
 174
 175 //! Selection of pull groups to be used in WHAM (one structure for each tpr file)
 176 typedef struct
 177 {
 178     int       n;         //!< total nr of pull groups in this tpr file
 179     int       nUse;      //!< nr of pull groups used
 180     gmx_bool *bUse;      //!< boolean array of size n. =1 if used, =0 if not
 181 } t_groupselection;
 182
 183 //! Parameters of WHAM
 184 typedef struct
 185 {
 186     /*!
 187      * \name Input stuff
 188      */
 189     /*!\{*/
 190     const char       *fnTpr, *fnPullf, *fnGroupsel;
 191     const char       *fnPdo, *fnPullx;            //!< file names of input
 192     gmx_bool          bTpr, bPullf, bPdo, bPullx; //!< input file types given?
 193     real              tmin, tmax, dt;             //!< only read input within tmin and tmax with dt
 194
 195     gmx_bool          bInitPotByIntegration;      //!< before WHAM, guess potential by force integration. Yields 1.5 to 2 times faster convergence
 196     int               stepUpdateContrib;          //!< update contribution table every ... iterations. Accelerates WHAM.
 197     int               nGroupsel;                  //!< if >0: use only certain group in WHAM, if ==0: use all groups
 198     t_groupselection *groupsel;                   //!< for each tpr file: which pull groups to use in WHAM?
 199     /*!\}*/
 200     /*!
 201      * \name Basic WHAM options
 202      */
 203     /*!\{*/
 204     int      bins;                   //!< nr of bins, min, max, and dz of profile
 205     real     min, max, dz;
 206     real     Temperature, Tolerance; //!< temperature, converged when probability changes less than Tolerance
 207     gmx_bool bCycl;                  //!< generate cyclic (periodic) PMF
 208     /*!\}*/
 209     /*!
 210      * \name Output control
 211      */
 212     /*!\{*/
 213     gmx_bool bLog;                   //!< energy output (instead of probability) for profile
 214     int      unit;                   //!< unit for PMF output kJ/mol or kT or kCal/mol
 215     gmx_bool bSym;                   //!< symmetrize PMF around z=0 after WHAM, useful for membranes etc.
 216     /*! \brief after wham, set prof to zero at this z-position.
 217      * When bootstrapping, set zProf0 to a "stable" reference position.
 218      */
 219     real         zProf0;
 220     gmx_bool     bProf0Set;              //!< setting profile to 0 at zProf0?
 221
 222     gmx_bool     bBoundsOnly, bHistOnly; //!< determine min and max, or write histograms and exit
 223     gmx_bool     bAuto;                  //!< determine min and max automatically but do not exit
 224
 225     gmx_bool     verbose;                //!< more noisy wham mode
 226     int          stepchange;             //!< print maximum change in prof after how many interations
 227     output_env_t oenv;                   //!< xvgr options
 228     /*!\}*/
 229     /*!
 230      * \name Autocorrelation stuff
 231      */
 232     /*!\{*/
 233     gmx_bool bTauIntGiven, bCalcTauInt; //!< IACT given or should be calculated?
 234     real     sigSmoothIact;             //!< sigma of Gaussian to smooth ACTs
 235     gmx_bool bAllowReduceIact;          //!< Allow to reduce ACTs during smoothing. Otherwise ACT are only increased during smoothing
 236     real     acTrestart;                //!< when computing ACT, time between restarting points
 237
 238     /* \brief Enforce the same weight for each umbella window, that is
 239      *  calculate with the same number of data points for
 240      *  each window. That can be reasonable, if the histograms
 241      *  have different length, but due to autocorrelation,
 242      *  a longer simulation should not have larger weightin wham.
 243      */
 244     gmx_bool bHistEq;
 245     /*!\}*/
 246
 247     /*!
 248      * \name Bootstrapping stuff
 249      */
 250     /*!\{*/
 251     int nBootStrap;              //!< nr of bootstraps (50 is usually enough)
 252
 253     /* \brief bootstrap method
 254      *
 255      * if == bsMethod_hist, consider complete histograms as independent
 256      * data points and, hence, only mix complete histograms.
 257      * if == bsMethod_BayesianHist, consider complete histograms
 258      * as independent data points, but assign random weights
 259      * to the histograms during the bootstrapping ("Bayesian bootstrap")
 260      * In case of long correlations (e.g., inside a channel), these
 261      * will yield a more realistic error.
 262      * if == bsMethod_traj(Gauss), generate synthetic histograms
 263      * for each given
 264      * histogram by generating an autocorrelated random sequence
 265      * that is distributed according to the respective given
 266      * histogram. With bsMethod_trajGauss, bootstrap from a Gaussian
 267      * (instead of from the umbrella histogram) to generate a new
 268      * histogram.
 269      */
 270     int bsMethod;
 271
 272     /* \brief  autocorrelation time (ACT) used to generate synthetic histograms. If ==0, use calculated ACF */
 273     real tauBootStrap;
 274
 275     /* \brief when mixing histograms, mix only histograms withing blocks
 276               long the reaction coordinate xi. Avoids gaps along xi. */
 277     int histBootStrapBlockLength;
 278
 279     int bsSeed;                    //!< random seed for bootstrapping
 280
 281     /* \brief Write cumulative distribution functions (CDFs) of histograms
 282               and write the generated histograms for each bootstrap */
 283     gmx_bool bs_verbose;
 284     /*!\}*/
 285     /*!
 286      * \name tabulated umbrella potential stuff
 287      */
 288     /*!\{*/
 289     gmx_bool  bTab;
 290     double   *tabX, *tabY, tabMin, tabMax, tabDz;
 291     int       tabNbins;
 292     /*!\}*/
 293     gmx_rng_t rng;                  //!< gromacs random number generator
 294 } t_UmbrellaOptions;
 295
 296 //! Make an umbrella window (may contain several histograms)
 297 t_UmbrellaWindow * initUmbrellaWindows(int nwin)
 298 {
 299     t_UmbrellaWindow *win;
 300     int               i;
 301     snew(win, nwin);
 302     for (i = 0; i < nwin; i++)
 303     {
 304         win[i].Histo    = win[i].cum  = 0;
 305         win[i].k        = win[i].pos  = win[i].z = 0;
 306         win[i].N        = win[i].Ntot = 0;
 307         win[i].g        = win[i].tau  = win[i].tausmooth = 0;
 308         win[i].bContrib = 0;
 309         win[i].ztime    = 0;
 310         win[i].forceAv  = 0;
 311         win[i].aver     = win[i].sigma = 0;
 312         win[i].bsWeight = 0;
 313     }
 314     return win;
 315 }
 316
 317 //! Delete an umbrella window (may contain several histograms)
 318 void freeUmbrellaWindows(t_UmbrellaWindow *win, int nwin)
 319 {
 320     int i, j;
 321     for (i = 0; i < nwin; i++)
 322     {
 323         if (win[i].Histo)
 324         {
 325             for (j = 0; j < win[i].nPull; j++)
 326             {
 327                 sfree(win[i].Histo[j]);
 328             }
 329         }
 330         if (win[i].cum)
 331         {
 332             for (j = 0; j < win[i].nPull; j++)
 333             {
 334                 sfree(win[i].cum[j]);
 335             }
 336         }
 337         if (win[i].bContrib)
 338         {
 339             for (j = 0; j < win[i].nPull; j++)
 340             {
 341                 sfree(win[i].bContrib[j]);
 342             }
 343         }
 344         sfree(win[i].Histo);
 345         sfree(win[i].cum);
 346         sfree(win[i].k);
 347         sfree(win[i].pos);
 348         sfree(win[i].z);
 349         sfree(win[i].N);
 350         sfree(win[i].Ntot);
 351         sfree(win[i].g);
 352         sfree(win[i].tau);
 353         sfree(win[i].tausmooth);
 354         sfree(win[i].bContrib);
 355         sfree(win[i].ztime);
 356         sfree(win[i].forceAv);
 357         sfree(win[i].aver);
 358         sfree(win[i].sigma);
 359         sfree(win[i].bsWeight);
 360     }
 361     sfree(win);
 362 }
 363
 364 /*! \brief
 365  * Read and setup tabulated umbrella potential
 366  */
 367 void setup_tab(const char *fn, t_UmbrellaOptions *opt)
 368 {
 369     int      i, ny, nl;
 370     double **y;
 371
 372     printf("Setting up tabulated potential from file %s\n", fn);
 373     nl            = read_xvg(fn, &y, &ny);
 374     opt->tabNbins = nl;
 375     if (ny != 2)
 376     {
 377         gmx_fatal(FARGS, "Found %d columns in %s. Expected 2.\n", ny, fn);
 378     }
 379     opt->tabMin = y[0][0];
 380     opt->tabMax = y[0][nl-1];
 381     opt->tabDz  = (opt->tabMax-opt->tabMin)/(nl-1);
 382     if (opt->tabDz <= 0)
 383     {
 384         gmx_fatal(FARGS, "The tabulated potential in %s must be provided in \n"
 385                   "ascending z-direction", fn);
 386     }
 387     for (i = 0; i < nl-1; i++)
 388     {
 389         if  (fabs(y[0][i+1]-y[0][i]-opt->tabDz) > opt->tabDz/1e6)
 390         {
 391             gmx_fatal(FARGS, "z-values in %s are not equally spaced.\n", ny, fn);
 392         }
 393     }
 394     snew(opt->tabY, nl);
 395     snew(opt->tabX, nl);
 396     for (i = 0; i < nl; i++)
 397     {
 398         opt->tabX[i] = y[0][i];
 399         opt->tabY[i] = y[1][i];
 400     }
 401     printf("Found equally spaced tabulated potential from %g to %g, spacing %g\n",
 402            opt->tabMin, opt->tabMax, opt->tabDz);
 403 }
 404
 405 //! Read the header of an PDO file (position, force const, nr of groups)
 406 void read_pdo_header(FILE * file, t_UmbrellaHeader * header, t_UmbrellaOptions *opt)
 407 {
 408     char               line[2048];
 409     char               Buffer0[256], Buffer1[256], Buffer2[256], Buffer3[256], Buffer4[256];
 410     int                i;
 411     std::istringstream ist;
 412
 413     /*  line 1 */
 414     if (fgets(line, 2048, file) == NULL)
 415     {
 416         gmx_fatal(FARGS, "Error reading header from pdo file\n");
 417     }
 418     ist.str(line);
 419     ist >> Buffer0 >> Buffer1 >> Buffer2;
 420     if (strcmp(Buffer1, "UMBRELLA"))
 421     {
 422         gmx_fatal(FARGS, "This does not appear to be a valid pdo file. Found %s, expected %s\n"
 423                   "(Found in first line: `%s')\n",
 424                   Buffer1, "UMBRELLA", line);
 425     }
 426     if (strcmp(Buffer2, "3.0"))
 427     {
 428         gmx_fatal(FARGS, "This does not appear to be a version 3.0 pdo file");
 429     }
 430
 431     /*  line 2 */
 432     if (fgets(line, 2048, file) == NULL)
 433     {
 434         gmx_fatal(FARGS, "Error reading header from pdo file\n");
 435     }
 436     ist.str(line);
 437     ist >> Buffer0 >> Buffer1 >> Buffer2 >> header->Dims[0] >> header->Dims[1] >> header->Dims[2];
 438     /* printf("%d %d %d\n", header->Dims[0],header->Dims[1],header->Dims[2]); */
 439
 440     header->nDim = header->Dims[0] + header->Dims[1] + header->Dims[2];
 441     if (header->nDim != 1)
 442     {
 443         gmx_fatal(FARGS, "Currently only supports one dimension");
 444     }
 445
 446     /* line3 */
 447     if (fgets(line, 2048, file) == NULL)
 448     {
 449         gmx_fatal(FARGS, "Error reading header from pdo file\n");
 450     }
 451     ist.str(line);
 452     ist >> Buffer0 >> Buffer1 >> header->nSkip;
 453
 454     /* line 4 */
 455     if (fgets(line, 2048, file) == NULL)
 456     {
 457         gmx_fatal(FARGS, "Error reading header from pdo file\n");
 458     }
 459     ist.str(line);
 460     ist >> Buffer0 >> Buffer1 >> Buffer2 >> header->Reference;
 461
 462     /* line 5 */
 463     if (fgets(line, 2048, file) == NULL)
 464     {
 465         gmx_fatal(FARGS, "Error reading header from pdo file\n");
 466     }
 467     ist.str(line);
 468     ist >> Buffer0 >> Buffer1 >> Buffer2 >> Buffer3 >> Buffer4 >> header->nPull;
 469
 470     if (opt->verbose)
 471     {
 472         printf("\tFound nPull=%d , nSkip=%d, ref=%s\n", header->nPull, header->nSkip,
 473                header->Reference);
 474     }
 475
 476     for (i = 0; i < header->nPull; ++i)
 477     {
 478         if (fgets(line, 2048, file) == NULL)
 479         {
 480             gmx_fatal(FARGS, "Error reading header from pdo file\n");
 481         }
 482         ist.str(line);
 483         ist >> Buffer0 >> Buffer1 >> Buffer2 >> header->PullName[i];
 484         ist >> Buffer0 >> Buffer1 >> header->UmbPos[i][0];
 485         ist >> Buffer0 >> Buffer1 >> header->UmbCons[i][0];
 486
 487         if (opt->verbose)
 488         {
 489             printf("\tpullgroup %d, pullname = %s, UmbPos = %g, UmbConst = %g\n",
 490                    i, header->PullName[i], header->UmbPos[i][0], header->UmbCons[i][0]);
 491         }
 492     }
 493
 494     if (fgets(line, 2048, file) == NULL)
 495     {
 496         gmx_fatal(FARGS, "Cannot read from file\n");
 497     }
 498     ist.str(line);
 499     ist >> Buffer3;
 500     if (strcmp(Buffer3, "#####") != 0)
 501     {
 502         gmx_fatal(FARGS, "Expected '#####', found %s. Hick.\n", Buffer3);
 503     }
 504 }
 505
 506 //! smarter fgets
 507 static char *fgets3(FILE *fp, char ptr[], int *len)
 508 {
 509     char *p;
 510     int   slen;
 511
 512     if (fgets(ptr, *len-1, fp) == NULL)
 513     {
 514         return NULL;
 515     }
 516     p = ptr;
 517     while ((strchr(ptr, '\n') == NULL) && (!feof(fp)))
 518     {
 519         /* This line is longer than len characters, let's increase len! */
 520         *len += STRLEN;
 521         p    += STRLEN;
 522         srenew(ptr, *len);
 523         if (fgets(p-1, STRLEN, fp) == NULL)
 524         {
 525             break;
 526         }
 527     }
 528     slen = strlen(ptr);
 529     if (ptr[slen-1] == '\n')
 530     {
 531         ptr[slen-1] = '\0';
 532     }
 533
 534     return ptr;
 535 }
 536
 537 /*! \brief Read the data columns of and PDO file.
 538  *
 539  *  TO DO: Get rid of the scanf function to avoid the clang warning.
 540  *         At the moment, this warning is avoided by hiding the format string
 541  *         the variable fmtlf.
 542  */
 543 void read_pdo_data(FILE * file, t_UmbrellaHeader * header,
 544                    int fileno, t_UmbrellaWindow * win,
 545                    t_UmbrellaOptions *opt,
 546                    gmx_bool bGetMinMax, real *mintmp, real *maxtmp)
 547 {
 548     int                i, inttemp, bins, count, ntot;
 549     real               min, max, minfound = 1e20, maxfound = -1e20;
 550     double             temp, time, time0 = 0, dt;
 551     char              *ptr    = 0;
 552     t_UmbrellaWindow * window = 0;
 553     gmx_bool           timeok, dt_ok = 1;
 554     char              *tmpbuf   = 0, fmt[256], fmtign[256], fmtlf[5] = "%lf";
 555     int                len      = STRLEN, dstep = 1;
 556     const int          blocklen = 4096;
 557     int               *lennow   = 0;
 558
 559     if (!bGetMinMax)
 560     {
 561         bins = opt->bins;
 562         min  = opt->min;
 563         max  = opt->max;
 564
 565         window = win+fileno;
 566         /* Need to alocate memory and set up structure */
 567         window->nPull = header->nPull;
 568         window->nBin  = bins;
 569
 570         snew(window->Histo, window->nPull);
 571         snew(window->z, window->nPull);
 572         snew(window->k, window->nPull);
 573         snew(window->pos, window->nPull);
 574         snew(window->N, window->nPull);
 575         snew(window->Ntot, window->nPull);
 576         snew(window->g, window->nPull);
 577         snew(window->bsWeight, window->nPull);
 578
 579         window->bContrib = 0;
 580
 581         if (opt->bCalcTauInt)
 582         {
 583             snew(window->ztime, window->nPull);
 584         }
 585         else
 586         {
 587             window->ztime = 0;
 588         }
 589         snew(lennow, window->nPull);
 590
 591         for (i = 0; i < window->nPull; ++i)
 592         {
 593             window->z[i]        = 1;
 594             window->bsWeight[i] = 1.;
 595             snew(window->Histo[i], bins);
 596             window->k[i]    = header->UmbCons[i][0];
 597             window->pos[i]  = header->UmbPos[i][0];
 598             window->N[i]    = 0;
 599             window->Ntot[i] = 0;
 600             window->g[i]    = 1.;
 601             if (opt->bCalcTauInt)
 602             {
 603                 window->ztime[i] = 0;
 604             }
 605         }
 606
 607         /* Done with setup */
 608     }
 609     else
 610     {
 611         minfound = 1e20;
 612         maxfound = -1e20;
 613         min      = max = bins = 0; /* Get rid of warnings */
 614     }
 615
 616     count = 0;
 617     snew(tmpbuf, len);
 618     while ( (ptr = fgets3(file, tmpbuf, &len)) != NULL)
 619     {
 620         trim(ptr);
 621
 622         if (ptr[0] == '#' || strlen(ptr) < 2)
 623         {
 624             continue;
 625         }
 626
 627         /* Initiate format string */
 628         fmtign[0] = '\0';
 629         strcat(fmtign, "%*s");
 630
 631         sscanf(ptr, fmtlf, &time); /* printf("Time %f\n",time); */
 632         /* Round time to fs */
 633         time = 1.0/1000*( static_cast<int> (time*1000+0.5) );
 634
 635         /* get time step of pdo file */
 636         if (count == 0)
 637         {
 638             time0 = time;
 639         }
 640         else if (count == 1)
 641         {
 642             dt = time-time0;
 643             if (opt->dt > 0.0)
 644             {
 645                 dstep = static_cast<int>(opt->dt/dt+0.5);
 646                 if (dstep == 0)
 647                 {
 648                     dstep = 1;
 649                 }
 650             }
 651             if (!bGetMinMax)
 652             {
 653                 window->dt = dt*dstep;
 654             }
 655         }
 656         count++;
 657
 658         dt_ok  = ((count-1)%dstep == 0);
 659         timeok = (dt_ok && time >= opt->tmin && time <= opt->tmax);
 660         /* if (opt->verbose)
 661            printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
 662            time,opt->tmin, opt->tmax, dt_ok,timeok); */
 663
 664         if (timeok)
 665         {
 666             for (i = 0; i < header->nPull; ++i)
 667             {
 668                 strcpy(fmt, fmtign);
 669                 strcat(fmt, "%lf");      /* Creating a format stings such as "%*s...%*s%lf" */
 670                 strcat(fmtign, "%*s");   /* ignoring one more entry in the next loop */
 671                 if (sscanf(ptr, fmt, &temp))
 672                 {
 673                     temp += header->UmbPos[i][0];
 674                     if (bGetMinMax)
 675                     {
 676                         if (temp < minfound)
 677                         {
 678                             minfound = temp;
 679                         }
 680                         if (temp > maxfound)
 681                         {
 682                             maxfound = temp;
 683                         }
 684                     }
 685                     else
 686                     {
 687                         if (opt->bCalcTauInt)
 688                         {
 689                             /* save time series for autocorrelation analysis */
 690                             ntot = window->Ntot[i];
 691                             if (ntot >= lennow[i])
 692                             {
 693                                 lennow[i] += blocklen;
 694                                 srenew(window->ztime[i], lennow[i]);
 695                             }
 696                             window->ztime[i][ntot] = temp;
 697                         }
 698
 699                         temp -= min;
 700                         temp /= (max-min);
 701                         temp *= bins;
 702                         temp  = floor(temp);
 703
 704                         inttemp = static_cast<int> (temp);
 705                         if (opt->bCycl)
 706                         {
 707                             if (inttemp < 0)
 708                             {
 709                                 inttemp += bins;
 710                             }
 711                             else if (inttemp >= bins)
 712                             {
 713                                 inttemp -= bins;
 714                             }
 715                         }
 716
 717                         if (inttemp >= 0 && inttemp < bins)
 718                         {
 719                             window->Histo[i][inttemp] += 1.;
 720                             window->N[i]++;
 721                         }
 722                         window->Ntot[i]++;
 723                     }
 724                 }
 725             }
 726         }
 727         if (time > opt->tmax)
 728         {
 729             if (opt->verbose)
 730             {
 731                 printf("time %f larger than tmax %f, stop reading pdo file\n", time, opt->tmax);
 732             }
 733             break;
 734         }
 735     }
 736
 737     if (bGetMinMax)
 738     {
 739         *mintmp = minfound;
 740         *maxtmp = maxfound;
 741     }
 742
 743     sfree(lennow);
 744     sfree(tmpbuf);
 745 }
 746
 747 /*! \brief Set identical weights for all histograms
 748  *
 749  * Normally, the weight is given by the number data points in each
 750  * histogram, together with the autocorrelation time. This can be overwritten
 751  * by this routine (not recommended). Since we now support autocorrelations, it is better to set
 752  * an appropriate autocorrelation times instead of using this function.
 753  */
 754 void enforceEqualWeights(t_UmbrellaWindow * window, int nWindows)
 755 {
 756     int    i, k, j, NEnforced;
 757     double ratio;
 758
 759     NEnforced = window[0].Ntot[0];
 760     printf("\nFound -hist-eq. Enforcing equal weights for all histograms, \ni.e. doing a "
 761            "non-weighted histogram analysis method. Ndata = %d\n", NEnforced);
 762     /* enforce all histograms to have the same weight as the very first histogram */
 763
 764     for (j = 0; j < nWindows; ++j)
 765     {
 766         for (k = 0; k < window[j].nPull; ++k)
 767         {
 768             ratio = 1.0*NEnforced/window[j].Ntot[k];
 769             for (i = 0; i < window[0].nBin; ++i)
 770             {
 771                 window[j].Histo[k][i] *= ratio;
 772             }
 773             window[j].N[k] = static_cast<int>(ratio*window[j].N[k] + 0.5);
 774         }
 775     }
 776 }
 777
 778 /*! \brief Simple linear interpolation between two given tabulated points
 779  */
 780 double tabulated_pot(double dist, t_UmbrellaOptions *opt)
 781 {
 782     int    jl, ju;
 783     double pl, pu, dz, dp;
 784
 785     jl = static_cast<int> (floor((dist-opt->tabMin)/opt->tabDz));
 786     ju = jl+1;
 787     if (jl < 0 || ju >= opt->tabNbins)
 788     {
 789         gmx_fatal(FARGS, "Distance %f out of bounds of tabulated potential (jl=%d, ju=%d).\n"
 790                   "Provide an extended table.", dist, jl, ju);
 791     }
 792     pl = opt->tabY[jl];
 793     pu = opt->tabY[ju];
 794     dz = dist-opt->tabX[jl];
 795     dp = (pu-pl)*dz/opt->tabDz;
 796     return pl+dp;
 797 }
 798
 799
 800 /*! \brief
 801  * Check which bins substiantially contribute (accelerates WHAM)
 802  *
 803  * Don't worry, that routine does not mean we compute the PMF in limited precision.
 804  * After rapid convergence (using only substiantal contributions), we always switch to
 805  * full precision.
 806  */
 807 void setup_acc_wham(double *profile, t_UmbrellaWindow * window, int nWindows,
 808                     t_UmbrellaOptions *opt)
 809 {
 810     int           i, j, k, nGrptot = 0, nContrib = 0, nTot = 0;
 811     double        U, min = opt->min, dz = opt->dz, temp, ztot_half, distance, ztot, contrib1, contrib2;
 812     gmx_bool      bAnyContrib;
 813     static int    bFirst = 1;
 814     static double wham_contrib_lim;
 815
 816     if (bFirst)
 817     {
 818         for (i = 0; i < nWindows; ++i)
 819         {
 820             nGrptot += window[i].nPull;
 821         }
 822         wham_contrib_lim = opt->Tolerance/nGrptot;
 823     }
 824
 825     ztot      = opt->max-opt->min;
 826     ztot_half = ztot/2;
 827
 828     for (i = 0; i < nWindows; ++i)
 829     {
 830         if (!window[i].bContrib)
 831         {
 832             snew(window[i].bContrib, window[i].nPull);
 833         }
 834         for (j = 0; j < window[i].nPull; ++j)
 835         {
 836             if (!window[i].bContrib[j])
 837             {
 838                 snew(window[i].bContrib[j], opt->bins);
 839             }
 840             bAnyContrib = FALSE;
 841             for (k = 0; k < opt->bins; ++k)
 842             {
 843                 temp     = (1.0*k+0.5)*dz+min;
 844                 distance = temp - window[i].pos[j];   /* distance to umbrella center */
 845                 if (opt->bCycl)
 846                 {                                     /* in cyclic wham:             */
 847                     if (distance > ztot_half)         /*    |distance| < ztot_half   */
 848                     {
 849                         distance -= ztot;
 850                     }
 851                     else if (distance < -ztot_half)
 852                     {
 853                         distance += ztot;
 854                     }
 855                 }
 856                 /* Note: there are two contributions to bin k in the wham equations:
 857                    i)  N[j]*exp(- U/(8.314e-3*opt->Temperature) + window[i].z[j])
 858                    ii) exp(- U/(8.314e-3*opt->Temperature))
 859                    where U is the umbrella potential
 860                    If any of these number is larger wham_contrib_lim, I set contrib=TRUE
 861                  */
 862
 863                 if (!opt->bTab)
 864                 {
 865                     U = 0.5*window[i].k[j]*sqr(distance);       /* harmonic potential assumed. */
 866                 }
 867                 else
 868                 {
 869                     U = tabulated_pot(distance, opt);            /* Use tabulated potential     */
 870
 871                 }
 872                 contrib1                 = profile[k]*exp(-U/(8.314e-3*opt->Temperature));
 873                 contrib2                 = window[i].N[j]*exp(-U/(8.314e-3*opt->Temperature) + window[i].z[j]);
 874                 window[i].bContrib[j][k] = (contrib1 > wham_contrib_lim || contrib2 > wham_contrib_lim);
 875                 bAnyContrib              = (bAnyContrib | window[i].bContrib[j][k]);
 876                 if (window[i].bContrib[j][k])
 877                 {
 878                     nContrib++;
 879                 }
 880                 nTot++;
 881             }
 882             /* If this histo is far outside min and max all bContrib may be FALSE,
 883                causing a floating point exception later on. To avoid that, switch
 884                them all to true.*/
 885             if (!bAnyContrib)
 886             {
 887                 for (k = 0; k < opt->bins; ++k)
 888                 {
 889                     window[i].bContrib[j][k] = TRUE;
 890                 }
 891             }
 892         }
 893     }
 894     if (bFirst)
 895     {
 896         printf("Initialized rapid wham stuff (contrib tolerance %g)\n"
 897                "Evaluating only %d of %d expressions.\n\n", wham_contrib_lim, nContrib, nTot);
 898     }
 899
 900     if (opt->verbose)
 901     {
 902         printf("Updated rapid wham stuff. (evaluating only %d of %d contributions)\n",
 903                nContrib, nTot);
 904     }
 905     bFirst = 0;
 906 }
 907
 908 //! Compute the PMF (one of the two main WHAM routines)
 909 void calc_profile(double *profile, t_UmbrellaWindow * window, int nWindows,
 910                   t_UmbrellaOptions *opt, gmx_bool bExact)
 911 {
 912     int    i, k, j;
 913     double num, ztot_half, ztot, distance, min = opt->min, dz = opt->dz;
 914     double denom, U = 0, temp = 0, invg;
 915
 916     ztot      = opt->max-opt->min;
 917     ztot_half = ztot/2;
 918
 919     for (i = 0; i < opt->bins; ++i)
 920     {
 921         num = denom = 0.;
 922         for (j = 0; j < nWindows; ++j)
 923         {
 924             for (k = 0; k < window[j].nPull; ++k)
 925             {
 926                 invg = 1.0/window[j].g[k] * window[j].bsWeight[k];
 927                 temp = (1.0*i+0.5)*dz+min;
 928                 num += invg*window[j].Histo[k][i];
 929
 930                 if (!(bExact || window[j].bContrib[k][i]))
 931                 {
 932                     continue;
 933                 }
 934                 distance = temp - window[j].pos[k];   /* distance to umbrella center */
 935                 if (opt->bCycl)
 936                 {                                     /* in cyclic wham:             */
 937                     if (distance > ztot_half)         /*    |distance| < ztot_half   */
 938                     {
 939                         distance -= ztot;
 940                     }
 941                     else if (distance < -ztot_half)
 942                     {
 943                         distance += ztot;
 944                     }
 945                 }
 946
 947                 if (!opt->bTab)
 948                 {
 949                     U = 0.5*window[j].k[k]*sqr(distance);       /* harmonic potential assumed. */
 950                 }
 951                 else
 952                 {
 953                     U = tabulated_pot(distance, opt);            /* Use tabulated potential     */
 954                 }
 955                 denom += invg*window[j].N[k]*exp(-U/(8.314e-3*opt->Temperature) + window[j].z[k]);
 956             }
 957         }
 958         profile[i] = num/denom;
 959     }
 960 }
 961
 962 //! Compute the free energy offsets z (one of the two main WHAM routines)
 963 double calc_z(double * profile, t_UmbrellaWindow * window, int nWindows,
 964               t_UmbrellaOptions *opt, gmx_bool bExact)
 965 {
 966     int    i, j, k;
 967     double U   = 0, min = opt->min, dz = opt->dz, temp, ztot_half, distance, ztot;
 968     double MAX = -1e20, total = 0;
 969
 970     ztot      = opt->max-opt->min;
 971     ztot_half = ztot/2;
 972
 973     for (i = 0; i < nWindows; ++i)
 974     {
 975         for (j = 0; j < window[i].nPull; ++j)
 976         {
 977             total = 0;
 978             for (k = 0; k < window[i].nBin; ++k)
 979             {
 980                 if (!(bExact || window[i].bContrib[j][k]))
 981                 {
 982                     continue;
 983                 }
 984                 temp     = (1.0*k+0.5)*dz+min;
 985                 distance = temp - window[i].pos[j];   /* distance to umbrella center */
 986                 if (opt->bCycl)
 987                 {                                     /* in cyclic wham:             */
 988                     if (distance > ztot_half)         /*    |distance| < ztot_half   */
 989                     {
 990                         distance -= ztot;
 991                     }
 992                     else if (distance < -ztot_half)
 993                     {
 994                         distance += ztot;
 995                     }
 996                 }
 997
 998                 if (!opt->bTab)
 999                 {
1000                     U = 0.5*window[i].k[j]*sqr(distance);       /* harmonic potential assumed. */
1001                 }
1002                 else
1003                 {
1004                     U = tabulated_pot(distance, opt);            /* Use tabulated potential     */
1005
1006                 }
1007                 total += profile[k]*exp(-U/(8.314e-3*opt->Temperature));
1008             }
1009             /* Avoid floating point exception if window is far outside min and max */
1010             if (total != 0.0)
1011             {
1012                 total = -log(total);
1013             }
1014             else
1015             {
1016                 total = 1000.0;
1017             }
1018             temp = fabs(total - window[i].z[j]);
1019             if (temp > MAX)
1020             {
1021                 MAX = temp;
1022             }
1023             window[i].z[j] = total;
1024         }
1025     }
1026     return MAX;
1027 }
1028
1029 //! Make PMF symmetric around 0 (useful e.g. for membranes)
1030 void symmetrizeProfile(double* profile, t_UmbrellaOptions *opt)
1031 {
1032     int     i, j, bins = opt->bins;
1033     double *prof2, min = opt->min, max = opt->max, dz = opt->dz, zsym, deltaz, profsym;
1034     double  z, z1;
1035
1036     if (min > 0. || max < 0.)
1037     {
1038         gmx_fatal(FARGS, "Cannot symmetrize profile around z=0 with min=%f and max=%f\n",
1039                   opt->min, opt->max);
1040     }
1041
1042     snew(prof2, bins);
1043
1044     for (i = 0; i < bins; i++)
1045     {
1046         z    = min+(i+0.5)*dz;
1047         zsym = -z;
1048         /* bin left of zsym */
1049         j = static_cast<int> (floor((zsym-min)/dz-0.5));
1050         if (j >= 0 && (j+1) < bins)
1051         {
1052             /* interpolate profile linearly between bins j and j+1 */
1053             z1      = min+(j+0.5)*dz;
1054             deltaz  = zsym-z1;
1055             profsym = profile[j] + (profile[j+1]-profile[j])/dz*deltaz;
1056             /* average between left and right */
1057             prof2[i] = 0.5*(profsym+profile[i]);
1058         }
1059         else
1060         {
1061             prof2[i] = profile[i];
1062         }
1063     }
1064
1065     memcpy(profile, prof2, bins*sizeof(double));
1066     sfree(prof2);
1067 }
1068
1069 //! Set energy unit (kJ/mol,kT,kCal/mol) and set it to zero at opt->zProf0
1070 void prof_normalization_and_unit(double * profile, t_UmbrellaOptions *opt)
1071 {
1072     int    i, bins, imin;
1073     double unit_factor = 1., R_MolarGasConst, diff;
1074
1075     R_MolarGasConst = 8.314472e-3; /* in kJ/(mol*K) */
1076     bins            = opt->bins;
1077
1078     /* Not log? Nothing to do! */
1079     if (!opt->bLog)
1080     {
1081         return;
1082     }
1083
1084     /* Get profile in units of kT, kJ/mol, or kCal/mol */
1085     if (opt->unit == en_kT)
1086     {
1087         unit_factor = 1.0;
1088     }
1089     else if (opt->unit == en_kJ)
1090     {
1091         unit_factor = R_MolarGasConst*opt->Temperature;
1092     }
1093     else if (opt->unit == en_kCal)
1094     {
1095         unit_factor = R_MolarGasConst*opt->Temperature/4.1868;
1096     }
1097     else
1098     {
1099         gmx_fatal(FARGS, "Sorry, I don't know this energy unit.");
1100     }
1101
1102     for (i = 0; i < bins; i++)
1103     {
1104         if (profile[i] > 0.0)
1105         {
1106             profile[i] = -log(profile[i])*unit_factor;
1107         }
1108     }
1109
1110     /* shift to zero at z=opt->zProf0 */
1111     if (!opt->bProf0Set)
1112     {
1113         diff = profile[0];
1114     }
1115     else
1116     {
1117         /* Get bin with shortest distance to opt->zProf0
1118            (-0.5 from bin position and +0.5 from rounding cancel) */
1119         imin = static_cast<int>((opt->zProf0-opt->min)/opt->dz);
1120         if (imin < 0)
1121         {
1122             imin = 0;
1123         }
1124         else if (imin >= bins)
1125         {
1126             imin = bins-1;
1127         }
1128         diff = profile[imin];
1129     }
1130
1131     /* Shift to zero */
1132     for (i = 0; i < bins; i++)
1133     {
1134         profile[i] -= diff;
1135     }
1136 }
1137
1138 //! Make an array of random integers (used for bootstrapping)
1139 void getRandomIntArray(int nPull, int blockLength, int* randomArray, gmx_rng_t rng)
1140 {
1141     int ipull, blockBase, nr, ipullRandom;
1142
1143     if (blockLength == 0)
1144     {
1145         blockLength = nPull;
1146     }
1147
1148     for (ipull = 0; ipull < nPull; ipull++)
1149     {
1150         blockBase = (ipull/blockLength)*blockLength;
1151         do
1152         {      /* make sure nothing bad happens in the last block */
1153             nr          = static_cast<int>(gmx_rng_uniform_real(rng)*blockLength);
1154             ipullRandom = blockBase + nr;
1155         }
1156         while (ipullRandom >= nPull);
1157         if (ipullRandom < 0 || ipullRandom >= nPull)
1158         {
1159             gmx_fatal(FARGS, "Ups, random iWin = %d, nPull = %d, nr = %d, "
1160                       "blockLength = %d, blockBase = %d\n",
1161                       ipullRandom, nPull, nr, blockLength, blockBase);
1162         }
1163         randomArray[ipull] = ipullRandom;
1164     }
1165     /*for (ipull=0; ipull<nPull; ipull++)
1166        printf("%d ",randomArray[ipull]); printf("\n"); */
1167 }
1168
1169 /*! \brief Set pull group information of a synthetic histogram
1170  *
1171  * This is used when bootstapping new trajectories and thereby create new histogtrams,
1172  * but it is not required if we bootstrap complete histograms.
1173  */
1174 void copy_pullgrp_to_synthwindow(t_UmbrellaWindow *synthWindow,
1175                                  t_UmbrellaWindow *thisWindow, int pullid)
1176 {
1177     synthWindow->N       [0] = thisWindow->N        [pullid];
1178     synthWindow->Histo   [0] = thisWindow->Histo    [pullid];
1179     synthWindow->pos     [0] = thisWindow->pos      [pullid];
1180     synthWindow->z       [0] = thisWindow->z        [pullid];
1181     synthWindow->k       [0] = thisWindow->k        [pullid];
1182     synthWindow->bContrib[0] = thisWindow->bContrib [pullid];
1183     synthWindow->g       [0] = thisWindow->g        [pullid];
1184     synthWindow->bsWeight[0] = thisWindow->bsWeight [pullid];
1185 }
1186
1187 /*! \brief Calculate cumulative distribution function of of all histograms.
1188  *
1189  * This allow to create random number sequences
1190  * which are distributed according to the histograms. Required to generate
1191  * the "synthetic" histograms for the Bootstrap method
1192  */
1193 void calc_cumulatives(t_UmbrellaWindow *window, int nWindows,
1194                       t_UmbrellaOptions *opt, const char *fnhist)
1195 {
1196     int    i, j, k, nbin;
1197     double last;
1198     char  *fn = 0, *buf = 0;
1199     FILE  *fp = 0;
1200
1201     if (opt->bs_verbose)
1202     {
1203         snew(fn, strlen(fnhist)+10);
1204         snew(buf, strlen(fnhist)+10);
1205         sprintf(fn, "%s_cumul.xvg", strncpy(buf, fnhist, strlen(fnhist)-4));
1206         fp = xvgropen(fn, "CDFs of umbrella windows", "z", "CDF", opt->oenv);
1207     }
1208
1209     nbin = opt->bins;
1210     for (i = 0; i < nWindows; i++)
1211     {
1212         snew(window[i].cum, window[i].nPull);
1213         for (j = 0; j < window[i].nPull; j++)
1214         {
1215             snew(window[i].cum[j], nbin+1);
1216             window[i].cum[j][0] = 0.;
1217             for (k = 1; k <= nbin; k++)
1218             {
1219                 window[i].cum[j][k] = window[i].cum[j][k-1]+window[i].Histo[j][k-1];
1220             }
1221
1222             /* normalize CDFs. Ensure cum[nbin]==1 */
1223             last = window[i].cum[j][nbin];
1224             for (k = 0; k <= nbin; k++)
1225             {
1226                 window[i].cum[j][k] /= last;
1227             }
1228         }
1229     }
1230
1231     printf("Cumulative distriubtion functions of all histograms created.\n");
1232     if (opt->bs_verbose)
1233     {
1234         for (k = 0; k <= nbin; k++)
1235         {
1236             fprintf(fp, "%g\t", opt->min+k*opt->dz);
1237             for (i = 0; i < nWindows; i++)
1238             {
1239                 for (j = 0; j < window[i].nPull; j++)
1240                 {
1241                     fprintf(fp, "%g\t", window[i].cum[j][k]);
1242                 }
1243             }
1244             fprintf(fp, "\n");
1245         }
1246         printf("Wrote cumulative distribution functions to %s\n", fn);
1247         ffclose(fp);
1248         sfree(fn);
1249         sfree(buf);
1250     }
1251 }
1252
1253
1254 /*! \brief Return j such that xx[j] <= x < xx[j+1]
1255  *
1256  *  This is used to generate a random sequence distributed according to a histogram
1257  */
1258 void searchCumulative(double xx[], int n, double x, int *j)
1259 {
1260     int ju, jm, jl;
1261
1262     jl = -1;
1263     ju = n;
1264     while (ju-jl > 1)
1265     {
1266         jm = (ju+jl) >> 1;
1267         if (x >= xx[jm])
1268         {
1269             jl = jm;
1270         }
1271         else
1272         {
1273             ju = jm;
1274         }
1275     }
1276     if (x == xx[0])
1277     {
1278         *j = 0;
1279     }
1280     else if (x == xx[n-1])
1281     {
1282         *j = n-2;
1283     }
1284     else
1285     {
1286         *j = jl;
1287     }
1288 }
1289
1290 //! Bootstrap new trajectories and thereby generate new (bootstrapped) histograms
1291 void create_synthetic_histo(t_UmbrellaWindow *synthWindow, t_UmbrellaWindow *thisWindow,
1292                             int pullid, t_UmbrellaOptions *opt)
1293 {
1294     int    N, i, nbins, r_index, ibin;
1295     double r, tausteps = 0.0, a, ap, dt, x, invsqrt2, g, y, sig = 0., z, mu = 0.;
1296     char   errstr[1024];
1297
1298     N     = thisWindow->N[pullid];
1299     dt    = thisWindow->dt;
1300     nbins = thisWindow->nBin;
1301
1302     /* tau = autocorrelation time */
1303     if (opt->tauBootStrap > 0.0)
1304     {
1305         tausteps = opt->tauBootStrap/dt;
1306     }
1307     else if (opt->bTauIntGiven || opt->bCalcTauInt)
1308     {
1309         /* calc tausteps from g=1+2tausteps */
1310         g        = thisWindow->g[pullid];
1311         tausteps = (g-1)/2;
1312     }
1313     else
1314     {
1315         sprintf(errstr,
1316                 "When generating hypothetical trajctories from given umbrella histograms,\n"
1317                 "autocorrelation times (ACTs) are required. Otherwise the statistical error\n"
1318                 "cannot be predicted. You have 3 options:\n"
1319                 "1) Make g_wham estimate the ACTs (options -ac and -acsig).\n"
1320                 "2) Calculate the ACTs by yourself (e.g. with g_analyze) and provide them\n");
1321         strcat(errstr,
1322                "   with option -iiact for all umbrella windows.\n"
1323                "3) If all ACTs are identical and know, you can define them with -bs-tau.\n"
1324                "   Use option (3) only if you are sure what you're doing, you may severely\n"
1325                "   underestimate the error if a too small ACT is given.\n");
1326         gmx_fatal(FARGS, errstr);
1327     }
1328
1329     synthWindow->N       [0] = N;
1330     synthWindow->pos     [0] = thisWindow->pos[pullid];
1331     synthWindow->z       [0] = thisWindow->z[pullid];
1332     synthWindow->k       [0] = thisWindow->k[pullid];
1333     synthWindow->bContrib[0] = thisWindow->bContrib[pullid];
1334     synthWindow->g       [0] = thisWindow->g       [pullid];
1335     synthWindow->bsWeight[0] = thisWindow->bsWeight[pullid];
1336
1337     for (i = 0; i < nbins; i++)
1338     {
1339         synthWindow->Histo[0][i] = 0.;
1340     }
1341
1342     if (opt->bsMethod == bsMethod_trajGauss)
1343     {
1344         sig = thisWindow->sigma [pullid];
1345         mu  = thisWindow->aver  [pullid];
1346     }
1347
1348     /* Genrate autocorrelated Gaussian random variable with autocorrelation time tau
1349        Use the following:
1350        If x and y are random numbers from N(0,1) (Gaussian with average 0 and sigma=1),
1351        then
1352        z = a*x + sqrt(1-a^2)*y
1353        is also from N(0,1), and cov(z,x) = a. Thus, by gerenating a sequence
1354        x' = a*x + sqrt(1-a^2)*y, the sequnce x(t) is from N(0,1) and has an autocorrelation
1355        function
1356        C(t) = exp(-t/tau) with tau=-1/ln(a)
1357
1358        Then, use error function to turn the Gaussian random variable into a uniformly
1359        distributed one in [0,1]. Eventually, use cumulative distribution function of
1360        histogram to get random variables distributed according to histogram.
1361        Note: The ACT of the flat distribution and of the generated histogram is not
1362        100% exactly tau, but near tau (my test was 3.8 instead of 4).
1363      */
1364     a        = exp(-1.0/tausteps);
1365     ap       = sqrt(1-a*a);
1366     invsqrt2 = 1./sqrt(2.0);
1367
1368     /* init random sequence */
1369     x = gmx_rng_gaussian_table(opt->rng);
1370
1371     if (opt->bsMethod == bsMethod_traj)
1372     {
1373         /* bootstrap points from the umbrella histograms */
1374         for (i = 0; i < N; i++)
1375         {
1376             y = gmx_rng_gaussian_table(opt->rng);
1377             x = a*x+ap*y;
1378             /* get flat distribution in [0,1] using cumulative distribution function of Gauusian
1379                Note: CDF(Gaussian) = 0.5*{1+erf[x/sqrt(2)]}
1380              */
1381             r = 0.5*(1+gmx_erf(x*invsqrt2));
1382             searchCumulative(thisWindow->cum[pullid], nbins+1, r, &r_index);
1383             synthWindow->Histo[0][r_index] += 1.;
1384         }
1385     }
1386     else if (opt->bsMethod == bsMethod_trajGauss)
1387     {
1388         /* bootstrap points from a Gaussian with the same average and sigma
1389            as the respective umbrella histogram. The idea was, that -given
1390            limited sampling- the bootstrapped histograms are otherwise biased
1391            from the limited sampling of the US histos. However, bootstrapping from
1392            the Gaussian seems to yield a similar estimate. */
1393         i = 0;
1394         while (i < N)
1395         {
1396             y    = gmx_rng_gaussian_table(opt->rng);
1397             x    = a*x+ap*y;
1398             z    = x*sig+mu;
1399             ibin = static_cast<int> (floor((z-opt->min)/opt->dz));
1400             if (opt->bCycl)
1401             {
1402                 if (ibin < 0)
1403                 {
1404                     while ( (ibin += nbins) < 0)
1405                     {
1406                         ;
1407                     }
1408                 }
1409                 else if (ibin >= nbins)
1410                 {
1411                     while ( (ibin -= nbins) >= nbins)
1412                     {
1413                         ;
1414                     }
1415                 }
1416             }
1417
1418             if (ibin >= 0 && ibin < nbins)
1419             {
1420                 synthWindow->Histo[0][ibin] += 1.;
1421                 i++;
1422             }
1423         }
1424     }
1425     else
1426     {
1427         gmx_fatal(FARGS, "Unknown bsMethod (id %d). That should not happen.\n", opt->bsMethod);
1428     }
1429 }
1430
1431 /*! \brief Write all histograms to a file
1432  *
1433  * If bs_index>=0, a number is added to the output file name to allow the ouput of all
1434  * sets of bootstrapped histograms.
1435  */
1436 void print_histograms(const char *fnhist, t_UmbrellaWindow * window, int nWindows,
1437                       int bs_index, t_UmbrellaOptions *opt)
1438 {
1439     char *fn = 0, *buf = 0, title[256];
1440     FILE *fp;
1441     int   bins, l, i, j;
1442
1443     if (bs_index >= 0)
1444     {
1445         snew(fn, strlen(fnhist)+10);
1446         snew(buf, strlen(fnhist)+1);
1447         sprintf(fn, "%s_bs%d.xvg", strncpy(buf, fnhist, strlen(fnhist)-4), bs_index);
1448         sprintf(title, "Umbrella histograms. Bootstrap #%d", bs_index);
1449     }
1450     else
1451     {
1452         fn = strdup(fnhist);
1453         strcpy(title, "Umbrella histograms");
1454     }
1455
1456     fp   = xvgropen(fn, title, "z", "count", opt->oenv);
1457     bins = opt->bins;
1458
1459     /* Write histograms */
1460     for (l = 0; l < bins; ++l)
1461     {
1462         fprintf(fp, "%e\t", (double)(l+0.5)*opt->dz+opt->min);
1463         for (i = 0; i < nWindows; ++i)
1464         {
1465             for (j = 0; j < window[i].nPull; ++j)
1466             {
1467                 fprintf(fp, "%e\t", window[i].Histo[j][l]);
1468             }
1469         }
1470         fprintf(fp, "\n");
1471     }
1472
1473     ffclose(fp);
1474     printf("Wrote %s\n", fn);
1475     if (bs_index >= 0)
1476     {
1477         sfree(buf);
1478     }
1479     sfree(fn);
1480 }
1481
1482 //! Used for qsort to sort random numbers
1483 int func_wham_is_larger(const void *a, const void *b)
1484 {
1485     double *aa, *bb;
1486     aa = (double*)a;
1487     bb = (double*)b;
1488     if (*aa < *bb)
1489     {
1490         return -1;
1491     }
1492     else if (*aa > *bb)
1493     {
1494         return 1;
1495     }
1496     else
1497     {
1498         return 0;
1499     }
1500 }
1501
1502 //! Make random weights for histograms for the Bayesian bootstrap of complete histograms)
1503 void setRandomBsWeights(t_UmbrellaWindow *synthwin, int nAllPull, t_UmbrellaOptions *opt)
1504 {
1505     int     i;
1506     double *r;
1507
1508     snew(r, nAllPull);
1509
1510     /* generate ordered random numbers between 0 and nAllPull  */
1511     for (i = 0; i < nAllPull-1; i++)
1512     {
1513         r[i] = gmx_rng_uniform_real(opt->rng) * nAllPull;
1514     }
1515     qsort((void *)r, nAllPull-1, sizeof(double), &func_wham_is_larger);
1516     r[nAllPull-1] = 1.0*nAllPull;
1517
1518     synthwin[0].bsWeight[0] = r[0];
1519     for (i = 1; i < nAllPull; i++)
1520     {
1521         synthwin[i].bsWeight[0] = r[i]-r[i-1];
1522     }
1523
1524     /* avoid to have zero weight by adding a tiny value */
1525     for (i = 0; i < nAllPull; i++)
1526     {
1527         if (synthwin[i].bsWeight[0] < 1e-5)
1528         {
1529             synthwin[i].bsWeight[0] = 1e-5;
1530         }
1531     }
1532
1533     sfree(r);
1534 }
1535
1536 //! The main bootstrapping routine
1537 void do_bootstrapping(const char *fnres, const char* fnprof, const char *fnhist,
1538                       char* ylabel, double *profile,
1539                       t_UmbrellaWindow * window, int nWindows, t_UmbrellaOptions *opt)
1540 {
1541     t_UmbrellaWindow * synthWindow;
1542     double            *bsProfile, *bsProfiles_av, *bsProfiles_av2, maxchange = 1e20, tmp, stddev;
1543     int                i, j, *randomArray = 0, winid, pullid, ib;
1544     int                iAllPull, nAllPull, *allPull_winId, *allPull_pullId;
1545     FILE              *fp;
1546     gmx_bool           bExact = FALSE;
1547
1548     /* init random generator */
1549     if (opt->bsSeed == -1)
1550     {
1551         opt->rng = gmx_rng_init(gmx_rng_make_seed());
1552     }
1553     else
1554     {
1555         opt->rng = gmx_rng_init(opt->bsSeed);
1556     }
1557
1558     snew(bsProfile,     opt->bins);
1559     snew(bsProfiles_av, opt->bins);
1560     snew(bsProfiles_av2, opt->bins);
1561
1562     /* Create array of all pull groups. Note that different windows
1563        may have different nr of pull groups
1564        First: Get total nr of pull groups */
1565     nAllPull = 0;
1566     for (i = 0; i < nWindows; i++)
1567     {
1568         nAllPull += window[i].nPull;
1569     }
1570     snew(allPull_winId, nAllPull);
1571     snew(allPull_pullId, nAllPull);
1572     iAllPull = 0;
1573     /* Setup one array of all pull groups */
1574     for (i = 0; i < nWindows; i++)
1575     {
1576         for (j = 0; j < window[i].nPull; j++)
1577         {
1578             allPull_winId[iAllPull]  = i;
1579             allPull_pullId[iAllPull] = j;
1580             iAllPull++;
1581         }
1582     }
1583
1584     /* setup stuff for synthetic windows */
1585     snew(synthWindow, nAllPull);
1586     for (i = 0; i < nAllPull; i++)
1587     {
1588         synthWindow[i].nPull = 1;
1589         synthWindow[i].nBin  = opt->bins;
1590         snew(synthWindow[i].Histo, 1);
1591         if (opt->bsMethod == bsMethod_traj || opt->bsMethod == bsMethod_trajGauss)
1592         {
1593             snew(synthWindow[i].Histo[0], opt->bins);
1594         }
1595         snew(synthWindow[i].N, 1);
1596         snew(synthWindow[i].pos, 1);
1597         snew(synthWindow[i].z, 1);
1598         snew(synthWindow[i].k, 1);
1599         snew(synthWindow[i].bContrib, 1);
1600         snew(synthWindow[i].g, 1);
1601         snew(synthWindow[i].bsWeight, 1);
1602     }
1603
1604     switch (opt->bsMethod)
1605     {
1606         case bsMethod_hist:
1607             snew(randomArray, nAllPull);
1608             printf("\n\nWhen computing statistical errors by bootstrapping entire histograms:\n");
1609             please_cite(stdout, "Hub2006");
1610             break;
1611         case bsMethod_BayesianHist:
1612             /* just copy all histogams into synthWindow array */
1613             for (i = 0; i < nAllPull; i++)
1614             {
1615                 winid  = allPull_winId [i];
1616                 pullid = allPull_pullId[i];
1617                 copy_pullgrp_to_synthwindow(synthWindow+i, window+winid, pullid);
1618             }
1619             break;
1620         case bsMethod_traj:
1621         case bsMethod_trajGauss:
1622             calc_cumulatives(window, nWindows, opt, fnhist);
1623             break;
1624         default:
1625             gmx_fatal(FARGS, "Unknown bootstrap method. That should not have happened.\n");
1626     }
1627
1628     /* do bootstrapping */
1629     fp = xvgropen(fnprof, "Boot strap profiles", "z", ylabel, opt->oenv);
1630     for (ib = 0; ib < opt->nBootStrap; ib++)
1631     {
1632         printf("  *******************************************\n"
1633                "  ******** Start bootstrap nr %d ************\n"
1634                "  *******************************************\n", ib+1);
1635
1636         switch (opt->bsMethod)
1637         {
1638             case bsMethod_hist:
1639                 /* bootstrap complete histograms from given histograms */
1640                 getRandomIntArray(nAllPull, opt->histBootStrapBlockLength, randomArray, opt->rng);
1641                 for (i = 0; i < nAllPull; i++)
1642                 {
1643                     winid  = allPull_winId [randomArray[i]];
1644                     pullid = allPull_pullId[randomArray[i]];
1645                     copy_pullgrp_to_synthwindow(synthWindow+i, window+winid, pullid);
1646                 }
1647                 break;
1648             case bsMethod_BayesianHist:
1649                 /* keep histos, but assign random weights ("Bayesian bootstrap") */
1650                 setRandomBsWeights(synthWindow, nAllPull, opt);
1651                 break;
1652             case bsMethod_traj:
1653             case bsMethod_trajGauss:
1654                 /* create new histos from given histos, that is generate new hypothetical
1655                    trajectories */
1656                 for (i = 0; i < nAllPull; i++)
1657                 {
1658                     winid  = allPull_winId[i];
1659                     pullid = allPull_pullId[i];
1660                     create_synthetic_histo(synthWindow+i, window+winid, pullid, opt);
1661                 }
1662                 break;
1663         }
1664
1665         /* write histos in case of verbose output */
1666         if (opt->bs_verbose)
1667         {
1668             print_histograms(fnhist, synthWindow, nAllPull, ib, opt);
1669         }
1670
1671         /* do wham */
1672         i         = 0;
1673         bExact    = FALSE;
1674         maxchange = 1e20;
1675         memcpy(bsProfile, profile, opt->bins*sizeof(double)); /* use profile as guess */
1676         do
1677         {
1678             if ( (i%opt->stepUpdateContrib) == 0)
1679             {
1680                 setup_acc_wham(bsProfile, synthWindow, nAllPull, opt);
1681             }
1682             if (maxchange < opt->Tolerance)
1683             {
1684                 bExact = TRUE;
1685             }
1686             if (((i%opt->stepchange) == 0 || i == 1) && i != 0)
1687             {
1688                 printf("\t%4d) Maximum change %e\n", i, maxchange);
1689             }
1690             calc_profile(bsProfile, synthWindow, nAllPull, opt, bExact);
1691             i++;
1692         }
1693         while ( (maxchange = calc_z(bsProfile, synthWindow, nAllPull, opt, bExact)) > opt->Tolerance || !bExact);
1694         printf("\tConverged in %d iterations. Final maximum change %g\n", i, maxchange);
1695
1696         if (opt->bLog)
1697         {
1698             prof_normalization_and_unit(bsProfile, opt);
1699         }
1700
1701         /* symmetrize profile around z=0 */
1702         if (opt->bSym)
1703         {
1704             symmetrizeProfile(bsProfile, opt);
1705         }
1706
1707         /* save stuff to get average and stddev */
1708         for (i = 0; i < opt->bins; i++)
1709         {
1710             tmp                = bsProfile[i];
1711             bsProfiles_av[i]  += tmp;
1712             bsProfiles_av2[i] += tmp*tmp;
1713             fprintf(fp, "%e\t%e\n", (i+0.5)*opt->dz+opt->min, tmp);
1714         }
1715         fprintf(fp, "&\n");
1716     }
1717     ffclose(fp);
1718
1719     /* write average and stddev */
1720     fp = xvgropen(fnres, "Average and stddev from bootstrapping", "z", ylabel, opt->oenv);
1721     fprintf(fp, "@TYPE xydy\n");
1722     for (i = 0; i < opt->bins; i++)
1723     {
1724         bsProfiles_av [i] /= opt->nBootStrap;
1725         bsProfiles_av2[i] /= opt->nBootStrap;
1726         tmp                = bsProfiles_av2[i]-sqr(bsProfiles_av[i]);
1727         stddev             = (tmp >= 0.) ? sqrt(tmp) : 0.; /* Catch rouding errors */
1728         fprintf(fp, "%e\t%e\t%e\n", (i+0.5)*opt->dz+opt->min, bsProfiles_av [i], stddev);
1729     }
1730     ffclose(fp);
1731     printf("Wrote boot strap result to %s\n", fnres);
1732 }
1733
1734 //! Return type of input file based on file extension (xvg, pdo, or tpr)
1735 int whaminFileType(char *fn)
1736 {
1737     int len;
1738     len = strlen(fn);
1739     if (strcmp(fn+len-3, "tpr") == 0)
1740     {
1741         return whamin_tpr;
1742     }
1743     else if (strcmp(fn+len-3, "xvg") == 0 || strcmp(fn+len-6, "xvg.gz") == 0)
1744     {
1745         return whamin_pullxf;
1746     }
1747     else if (strcmp(fn+len-3, "pdo") == 0 || strcmp(fn+len-6, "pdo.gz") == 0)
1748     {
1749         return whamin_pdo;
1750     }
1751     else
1752     {
1753         gmx_fatal(FARGS, "Unknown file type of %s. Should be tpr, xvg, or pdo.\n", fn);
1754     }
1755     return whamin_unknown;
1756 }
1757
1758 //! Read the files names in pdo-files.dat, pullf/x-files.dat, tpr-files.dat
1759 void read_wham_in(const char *fn, char ***filenamesRet, int *nfilesRet,
1760                   t_UmbrellaOptions *opt)
1761 {
1762     char **filename = 0, tmp[WHAM_MAXFILELEN+2];
1763     int    nread, sizenow, i, block = 1;
1764     FILE  *fp;
1765
1766     fp      = ffopen(fn, "r");
1767     nread   = 0;
1768     sizenow = 0;
1769     while (fgets(tmp, sizeof(tmp), fp) != NULL)
1770     {
1771         if (strlen(tmp) >= WHAM_MAXFILELEN)
1772         {
1773             gmx_fatal(FARGS, "Filename too long in %s. Only %d characters allowed.\n", fn, WHAM_MAXFILELEN);
1774         }
1775         if (nread >= sizenow)
1776         {
1777             sizenow += block;
1778             srenew(filename, sizenow);
1779             for (i = sizenow-block; i < sizenow; i++)
1780             {
1781                 snew(filename[i], WHAM_MAXFILELEN);
1782             }
1783         }
1784         /* remove newline if there is one */
1785         if (tmp[strlen(tmp)-1] == '\n')
1786         {
1787             tmp[strlen(tmp)-1] = '\0';
1788         }
1789         strcpy(filename[nread], tmp);
1790         if (opt->verbose)
1791         {
1792             printf("Found file %s in %s\n", filename[nread], fn);
1793         }
1794         nread++;
1795     }
1796     *filenamesRet = filename;
1797     *nfilesRet    = nread;
1798 }
1799
1800 //! Open a file or a pipe to a gzipped file
1801 FILE *open_pdo_pipe(const char *fn, t_UmbrellaOptions *opt, gmx_bool *bPipeOpen)
1802 {
1803     char            Buffer[1024], gunzip[1024], *Path = 0;
1804     FILE           *pipe   = 0;
1805     static gmx_bool bFirst = 1;
1806
1807     /* gzipped pdo file? */
1808     if ((strcmp(fn+strlen(fn)-3, ".gz") == 0))
1809     {
1810         /* search gunzip executable */
1811         if (!(Path = getenv("GMX_PATH_GZIP")))
1812         {
1813             if (gmx_fexist("/bin/gunzip"))
1814             {
1815                 sprintf(gunzip, "%s", "/bin/gunzip");
1816             }
1817             else if (gmx_fexist("/usr/bin/gunzip"))
1818             {
1819                 sprintf(gunzip, "%s", "/usr/bin/gunzip");
1820             }
1821             else
1822             {
1823                 gmx_fatal(FARGS, "Cannot find executable gunzip in /bin or /usr/bin.\n"
1824                           "You may want to define the path to gunzip "
1825                           "with the environment variable GMX_PATH_GZIP.", gunzip);
1826             }
1827         }
1828         else
1829         {
1830             sprintf(gunzip, "%s/gunzip", Path);
1831             if (!gmx_fexist(gunzip))
1832             {
1833                 gmx_fatal(FARGS, "Cannot find executable %s. Please define the path to gunzip"
1834                           " in the environmental varialbe GMX_PATH_GZIP.", gunzip);
1835             }
1836         }
1837         if (bFirst)
1838         {
1839             printf("Using gunzig executable %s\n", gunzip);
1840             bFirst = 0;
1841         }
1842         if (!gmx_fexist(fn))
1843         {
1844             gmx_fatal(FARGS, "File %s does not exist.\n", fn);
1845         }
1846         sprintf(Buffer, "%s -c < %s", gunzip, fn);
1847         if (opt->verbose)
1848         {
1849             printf("Executing command '%s'\n", Buffer);
1850         }
1851 #ifdef HAVE_PIPES
1852         if ((pipe = popen(Buffer, "r")) == NULL)
1853         {
1854             gmx_fatal(FARGS, "Unable to open pipe to `%s'\n", Buffer);
1855         }
1856 #else
1857         gmx_fatal(FARGS, "Cannot open a compressed file on platform without pipe support");
1858 #endif
1859         *bPipeOpen = TRUE;
1860     }
1861     else
1862     {
1863         pipe       = ffopen(fn, "r");
1864         *bPipeOpen = FALSE;
1865     }
1866
1867     return pipe;
1868 }
1869
1870 //! Close file or pipe
1871 void pdo_close_file(FILE *fp)
1872 {
1873 #ifdef HAVE_PIPES
1874     pclose(fp);
1875 #else
1876     ffclose(fp);
1877 #endif
1878 }
1879
1880 //! Reading all pdo files
1881 void read_pdo_files(char **fn, int nfiles, t_UmbrellaHeader* header,
1882                     t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
1883 {
1884     FILE    *file;
1885     real     mintmp, maxtmp, done = 0.;
1886     int      i;
1887     gmx_bool bPipeOpen;
1888     /* char Buffer0[1000]; */
1889
1890     if (nfiles < 1)
1891     {
1892         gmx_fatal(FARGS, "No files found. Hick.");
1893     }
1894
1895     /* if min and max are not given, get min and max from the input files */
1896     if (opt->bAuto)
1897     {
1898         printf("Automatic determination of boundaries from %d pdo files...\n", nfiles);
1899         opt->min = 1e20;
1900         opt->max = -1e20;
1901         for (i = 0; i < nfiles; ++i)
1902         {
1903             file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
1904             /*fgets(Buffer0,999,file);
1905                fprintf(stderr,"First line '%s'\n",Buffer0); */
1906             done = 100.0*(i+1)/nfiles;
1907             printf("\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
1908             if (opt->verbose)
1909             {
1910                 printf("\n");
1911             }
1912             read_pdo_header(file, header, opt);
1913             /* here only determine min and max of this window */
1914             read_pdo_data(file, header, i, NULL, opt, TRUE, &mintmp, &maxtmp);
1915             if (maxtmp > opt->max)
1916             {
1917                 opt->max = maxtmp;
1918             }
1919             if (mintmp < opt->min)
1920             {
1921                 opt->min = mintmp;
1922             }
1923             if (bPipeOpen)
1924             {
1925                 pdo_close_file(file);
1926             }
1927             else
1928             {
1929                 ffclose(file);
1930             }
1931         }
1932         printf("\n");
1933         printf("\nDetermined boundaries to %f and %f\n\n", opt->min, opt->max);
1934         if (opt->bBoundsOnly)
1935         {
1936             printf("Found option -boundsonly, now exiting.\n");
1937             exit (0);
1938         }
1939     }
1940     /* store stepsize in profile */
1941     opt->dz = (opt->max-opt->min)/opt->bins;
1942
1943     /* Having min and max, we read in all files */
1944     /* Loop over all files */
1945     for (i = 0; i < nfiles; ++i)
1946     {
1947         done = 100.0*(i+1)/nfiles;
1948         printf("\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
1949         if (opt->verbose)
1950         {
1951             printf("\n");
1952         }
1953         file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
1954         read_pdo_header(file, header, opt);
1955         /* load data into window */
1956         read_pdo_data(file, header, i, window, opt, FALSE, NULL, NULL);
1957         if ((window+i)->Ntot[0] == 0)
1958         {
1959             fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fn[i]);
1960         }
1961         if (bPipeOpen)
1962         {
1963             pdo_close_file(file);
1964         }
1965         else
1966         {
1967             ffclose(file);
1968         }
1969     }
1970     printf("\n");
1971     for (i = 0; i < nfiles; ++i)
1972     {
1973         sfree(fn[i]);
1974     }
1975     sfree(fn);
1976 }
1977
1978 //! translate 0/1 to N/Y to write pull dimensions
1979 #define int2YN(a) (((a) == 0) ? ("N") : ("Y"))
1980
1981 //! Read pull groups from a tpr file (including position, force const, geometry, number of groups)
1982 void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions *opt)
1983 {
1984     t_inputrec  ir;
1985     int         i, ncrd;
1986     t_state     state;
1987     static int  first = 1;
1988
1989     /* printf("Reading %s \n",fn); */
1990     read_tpx_state(fn, &ir, &state, NULL, NULL);
1991
1992     if (ir.ePull != epullUMBRELLA)
1993     {
1994         gmx_fatal(FARGS, "This is not a tpr of an umbrella simulation. Found pull type \"%s\" "
1995                   " (ir.ePull = %d)\n", epull_names[ir.ePull], ir.ePull);
1996     }
1997
1998     /* nr of pull groups */
1999     ncrd = ir.pull->ncoord;
2000     if (ncrd < 1)
2001     {
2002         gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull coordinates\n", ncrd);
2003     }
2004
2005     header->npullcrds     = ir.pull->ncoord;
2006     header->pull_geometry = ir.pull->eGeom;
2007     header->bPrintRef     = ir.pull->bPrintRef;
2008     copy_ivec(ir.pull->dim, header->pull_dim);
2009     header->pull_ndim = header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2];
2010     snew(header->k, ncrd);
2011     snew(header->init_dist, ncrd);
2012     snew(header->umbInitDist, ncrd);
2013
2014     /* only z-direction with epullgCYL? */
2015     if (header->pull_geometry == epullgCYL)
2016     {
2017         if (header->pull_dim[XX] || header->pull_dim[YY] || (!header->pull_dim[ZZ]))
2018         {
2019             gmx_fatal(FARGS, "With pull geometry 'cylinder', expected pulling in Z direction only.\n"
2020                       "However, found dimensions [%s %s %s]\n",
2021                       int2YN(header->pull_dim[XX]), int2YN(header->pull_dim[YY]),
2022                       int2YN(header->pull_dim[ZZ]));
2023         }
2024     }
2025
2026     for (i = 0; i < ncrd; i++)
2027     {
2028         header->k[i] = ir.pull->coord[i].k;
2029         if (header->k[i] == 0.0)
2030         {
2031             gmx_fatal(FARGS, "Pull coordinate %d has force constant of of 0.0 in %s.\n"
2032                       "That doesn't seem to be an Umbrella tpr.\n",
2033                       i, fn);
2034         }
2035         header->init_dist[i] =  ir.pull->coord[i].init;
2036
2037         /* initial distance to reference */
2038         switch (header->pull_geometry)
2039         {
2040             case epullgCYL:
2041             /* umbrella distance stored in init_dist[i] for geometry cylinder (not in ...[i][ZZ]) */
2042             case epullgDIST:
2043             case epullgDIR:
2044             case epullgDIRPBC:
2045                 header->umbInitDist[i] = header->init_dist[i];
2046                 break;
2047             default:
2048                 gmx_fatal(FARGS, "Pull geometry %s not supported\n", epullg_names[header->pull_geometry]);
2049         }
2050     }
2051
2052     if (opt->verbose || first)
2053     {
2054         printf("File %s, %d coordinates, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n"
2055                "\tPull group coordinates%s expected in pullx files.\n",
2056                fn, header->npullcrds, epullg_names[header->pull_geometry],
2057                int2YN(header->pull_dim[0]), int2YN(header->pull_dim[1]), int2YN(header->pull_dim[2]),
2058                header->pull_ndim, (header->bPrintRef ? "" : " not"));
2059         for (i = 0; i < ncrd; i++)
2060         {
2061             printf("\tcrd %d) k = %-5g  position = %g\n", i, header->k[i], header->umbInitDist[i]);
2062         }
2063     }
2064     if (!opt->verbose && first)
2065     {
2066         printf("\tUse option -v to see this output for all input tpr files\n\n");
2067     }
2068
2069     first = 0;
2070 }
2071
2072 //! 2-norm in a ndim-dimensional space
2073 double dist_ndim(double **dx, int ndim, int line)
2074 {
2075     int    i;
2076     double r2 = 0.;
2077     for (i = 0; i < ndim; i++)
2078     {
2079         r2 += sqr(dx[i][line]);
2080     }
2081     return sqrt(r2);
2082 }
2083
2084 //! Read pullx.xvg or pullf.xvg
2085 void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header,
2086                   t_UmbrellaWindow * window,
2087                   t_UmbrellaOptions *opt,
2088                   gmx_bool bGetMinMax, real *mintmp, real *maxtmp,
2089                   t_groupselection *groupsel)
2090 {
2091     double        **y = 0, pos = 0., t, force, time0 = 0., dt;
2092     int             ny, nt, bins, ibin, i, g, gUsed, dstep = 1, nColPerCrd, nColRefEachCrd, nColExpect, ntot, column;
2093     real            min, max, minfound = 1e20, maxfound = -1e20;
2094     gmx_bool        dt_ok, timeok, bHaveForce;
2095     const char     *quantity;
2096     const int       blocklen = 4096;
2097     int            *lennow   = 0;
2098     static gmx_bool bFirst   = TRUE;
2099
2100     /*
2101      * Data columns in pull output:
2102      *  - in force output pullf.xvg:
2103      *    No reference columns, one column per pull coordinate
2104      *
2105      *  - in position output pullx.xvg
2106      *    bPrintRef == TRUE:  for each pull coordinate: ndim reference columns, and ndim dx columns
2107      *    bPrintRef == FALSE: for each pull coordinate: no   reference columns, but ndim dx columns
2108      */
2109
2110     nColPerCrd = opt->bPullx ? header->pull_ndim : 1;
2111     quantity   = opt->bPullx ? "position" : "force";
2112
2113     if (opt->bPullx && header->bPrintRef)
2114     {
2115         nColRefEachCrd = header->pull_ndim;
2116     }
2117     else
2118     {
2119         nColRefEachCrd = 0;
2120     }
2121
2122     nColExpect = 1 + header->npullcrds*(nColRefEachCrd+nColPerCrd);
2123     bHaveForce = opt->bPullf;
2124
2125     /* With geometry "distance" or "distance_periodic", only force reading is supported so far.
2126        That avoids the somewhat tedious extraction of the reaction coordinate from the pullx files.
2127        Sorry for the laziness, this is a To-do. */
2128     if  ( (header->pull_geometry == epullgDIR || header->pull_geometry == epullgDIRPBC)
2129           && opt->bPullx)
2130     {
2131         gmx_fatal(FARGS, "With pull geometries \"direction\" and \"direction_periodic\", only pull force "
2132                   "reading \n(option -if) is supported at present, "
2133                   "not pull position reading (options -ix).\nMake sure mdrun writes the pull "
2134                   "forces (pullf.xvg files)\nand provide them to g_wham with option -if.",
2135                   epullg_names[header->pull_geometry]);
2136     }
2137
2138     nt = read_xvg(fn, &y, &ny);
2139
2140     /* Check consistency */
2141     if (nt < 1)
2142     {
2143         gmx_fatal(FARGS, "Empty pull %s file %s\n", quantity, fn);
2144     }
2145     if (bFirst)
2146     {
2147         printf("Reading pull %s file with pull geometry %s and %d pull dimensions\n",
2148                bHaveForce ? "force" : "position", epullg_names[header->pull_geometry],
2149                header->pull_ndim);
2150         printf("Expecting these columns in pull file:\n"
2151                "\t%d reference columns for each individual pull coordinate\n"
2152                "\t%d data columns for each pull coordinate\n", nColRefEachCrd, nColPerCrd);
2153         printf("With %d pull groups, expect %d columns (including the time column)\n", header->npullcrds, nColExpect);
2154         bFirst = FALSE;
2155     }
2156     if (ny != nColExpect)
2157     {
2158         gmx_fatal(FARGS, "Found %d pull coodinates in %s,\n but %d data columns in %s (expected %d)\n"
2159                   "\nMaybe you confused options -ix and -if ?\n",
2160                   header->npullcrds, fntpr, ny-1, fn, nColExpect-1);
2161     }
2162
2163     if (opt->verbose)
2164     {
2165         printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerCrd, quantity, fn);
2166     }
2167
2168     if (!bGetMinMax)
2169     {
2170         bins = opt->bins;
2171         min  = opt->min;
2172         max  = opt->max;
2173         if (nt > 1)
2174         {
2175             window->dt = y[0][1]-y[0][0];
2176         }
2177         else if (opt->nBootStrap && opt->tauBootStrap != 0.0)
2178         {
2179             fprintf(stderr, "\n *** WARNING, Could not determine time step in %s\n", fn);
2180         }
2181
2182         /* Need to alocate memory and set up structure */
2183
2184         if (groupsel)
2185         {
2186             /* Use only groups selected with option -is file */
2187             if (header->npullcrds != groupsel->n)
2188             {
2189                 gmx_fatal(FARGS, "tpr file contains %d pull groups, but expected %d from group selection file\n",
2190                           header->npullcrds, groupsel->n);
2191             }
2192             window->nPull = groupsel->nUse;
2193         }
2194         else
2195         {
2196             window->nPull = header->npullcrds;
2197         }
2198
2199         window->nBin = bins;
2200         snew(window->Histo, window->nPull);
2201         snew(window->z, window->nPull);
2202         snew(window->k, window->nPull);
2203         snew(window->pos, window->nPull);
2204         snew(window->N, window->nPull);
2205         snew(window->Ntot, window->nPull);
2206         snew(window->g, window->nPull);
2207         snew(window->bsWeight, window->nPull);
2208         window->bContrib = 0;
2209
2210         if (opt->bCalcTauInt)
2211         {
2212             snew(window->ztime, window->nPull);
2213         }
2214         else
2215         {
2216             window->ztime = NULL;
2217         }
2218         snew(lennow, window->nPull);
2219
2220         for (g = 0; g < window->nPull; ++g)
2221         {
2222             window->z[g]        = 1;
2223             window->bsWeight[g] = 1.;
2224             snew(window->Histo[g], bins);
2225             window->N[g]    = 0;
2226             window->Ntot[g] = 0;
2227             window->g[g]    = 1.;
2228             if (opt->bCalcTauInt)
2229             {
2230                 window->ztime[g] = NULL;
2231             }
2232         }
2233
2234         /* Copying umbrella center and force const is more involved since not
2235            all pull groups from header (tpr file) may be used in window variable */
2236         for (g = 0, gUsed = 0; g < header->npullcrds; ++g)
2237         {
2238             if (groupsel && (groupsel->bUse[g] == FALSE))
2239             {
2240                 continue;
2241             }
2242             window->k[gUsed]   = header->k[g];
2243             window->pos[gUsed] = header->umbInitDist[g];
2244             gUsed++;
2245         }
2246     }
2247     else
2248     {   /* only determine min and max */
2249         minfound = 1e20;
2250         maxfound = -1e20;
2251         min      = max = bins = 0; /* Get rid of warnings */
2252     }
2253
2254
2255     for (i = 0; i < nt; i++)
2256     {
2257         /* Do you want that time frame? */
2258         t = 1.0/1000*( static_cast<int> ((y[0][i]*1000) + 0.5)); /* round time to fs */
2259
2260         /* get time step of pdo file and get dstep from opt->dt */
2261         if (i == 0)
2262         {
2263             time0 = t;
2264         }
2265         else if (i == 1)
2266         {
2267             dt = t-time0;
2268             if (opt->dt > 0.0)
2269             {
2270                 dstep = static_cast<int>(opt->dt/dt+0.5);
2271                 if (dstep == 0)
2272                 {
2273                     dstep = 1;
2274                 }
2275             }
2276             if (!bGetMinMax)
2277             {
2278                 window->dt = dt*dstep;
2279             }
2280         }
2281
2282         dt_ok  = (i%dstep == 0);
2283         timeok = (dt_ok && t >= opt->tmin && t <= opt->tmax);
2284         /*if (opt->verbose)
2285            printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
2286            t,opt->tmin, opt->tmax, dt_ok,timeok); */
2287         if (timeok)
2288         {
2289             /* Note: if groupsel == NULL:
2290              *          all groups in pullf/x file are stored in this window, and gUsed == g
2291              *       if groupsel != NULL:
2292              *          only groups with groupsel.bUse[g]==TRUE are stored. gUsed is not always equal g
2293              */
2294             gUsed = -1;
2295             for (g = 0; g < header->npullcrds; ++g)
2296             {
2297                 /* was this group selected for application in WHAM? */
2298                 if (groupsel && (groupsel->bUse[g] == FALSE))
2299                 {
2300                     continue;
2301                 }
2302
2303                 gUsed++;
2304
2305                 if (bHaveForce)
2306                 {
2307                     /* y has 1 time column y[0] and one column per force y[1],...,y[nCrds] */
2308                     force = y[g+1][i];
2309                     pos   = -force/header->k[g] + header->umbInitDist[g];
2310                 }
2311                 else
2312                 {
2313                     /* pick the right column from:
2314                      * time ref1[ndim] coord1[ndim] ref2[ndim] coord2[ndim] ...
2315                      */
2316                     column = 1 +  nColRefEachCrd + g*(nColRefEachCrd+nColPerCrd);
2317                     switch (header->pull_geometry)
2318                     {
2319                         case epullgDIST:
2320                             pos = dist_ndim(y + column, header->pull_ndim, i);
2321                             break;
2322                         case epullgCYL:
2323                             pos = y[column][i];
2324                             break;
2325                         default:
2326                             gmx_fatal(FARGS, "Bad error, this error should have been catched before. Ups.\n");
2327                     }
2328                 }
2329
2330                 /* printf("crd %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
2331                 if (bGetMinMax)
2332                 {
2333                     if (pos < minfound)
2334                     {
2335                         minfound = pos;
2336                     }
2337                     if (pos > maxfound)
2338                     {
2339                         maxfound = pos;
2340                     }
2341                 }
2342                 else
2343                 {
2344                     if (gUsed >= window->nPull)
2345                     {
2346                         gmx_fatal(FARGS, "gUsed too large (%d, nPull=%d). This error should have been catched before.\n",
2347                                   gUsed, window->nPull);
2348                     }
2349
2350                     if (opt->bCalcTauInt && !bGetMinMax)
2351                     {
2352                         /* save time series for autocorrelation analysis */
2353                         ntot = window->Ntot[gUsed];
2354                         /* printf("i %d, ntot %d, lennow[g] = %d\n",i,ntot,lennow[g]); */
2355                         if (ntot >= lennow[gUsed])
2356                         {
2357                             lennow[gUsed] += blocklen;
2358                             srenew(window->ztime[gUsed], lennow[gUsed]);
2359                         }
2360                         window->ztime[gUsed][ntot] = pos;
2361                     }
2362
2363                     ibin = static_cast<int> (floor((pos-min)/(max-min)*bins));
2364                     if (opt->bCycl)
2365                     {
2366                         if (ibin < 0)
2367                         {
2368                             while ( (ibin += bins) < 0)
2369                             {
2370                                 ;
2371                             }
2372                         }
2373                         else if (ibin >= bins)
2374                         {
2375                             while ( (ibin -= bins) >= bins)
2376                             {
2377                                 ;
2378                             }
2379                         }
2380                     }
2381                     if (ibin >= 0 && ibin < bins)
2382                     {
2383                         window->Histo[gUsed][ibin] += 1.;
2384                         window->N[gUsed]++;
2385                     }
2386                     window->Ntot[gUsed]++;
2387                 }
2388             }
2389         }
2390         else if (t > opt->tmax)
2391         {
2392             if (opt->verbose)
2393             {
2394                 printf("time %f larger than tmax %f, stop reading this pullx/pullf file\n", t, opt->tmax);
2395             }
2396             break;
2397         }
2398     }
2399
2400     if (bGetMinMax)
2401     {
2402         *mintmp = minfound;
2403         *maxtmp = maxfound;
2404     }
2405     sfree(lennow);
2406     for (i = 0; i < ny; i++)
2407     {
2408         sfree(y[i]);
2409     }
2410 }
2411
2412 //! read pullf-files.dat or pullx-files.dat and tpr-files.dat
2413 void read_tpr_pullxf_files(char **fnTprs, char **fnPull, int nfiles,
2414                            t_UmbrellaHeader* header,
2415                            t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
2416 {
2417     int  i;
2418     real mintmp, maxtmp;
2419
2420     printf("Reading %d tpr and pullf files\n", nfiles/2);
2421
2422     /* min and max not given? */
2423     if (opt->bAuto)
2424     {
2425         printf("Automatic determination of boundaries...\n");
2426         opt->min = 1e20;
2427         opt->max = -1e20;
2428         for (i = 0; i < nfiles; i++)
2429         {
2430             if (whaminFileType(fnTprs[i]) != whamin_tpr)
2431             {
2432                 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a tpr file\n", i);
2433             }
2434             read_tpr_header(fnTprs[i], header, opt);
2435             if (whaminFileType(fnPull[i]) != whamin_pullxf)
2436             {
2437                 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
2438             }
2439             read_pull_xf(fnPull[i], fnTprs[i], header, NULL, opt, TRUE, &mintmp, &maxtmp,
2440                          (opt->nGroupsel > 0) ? &opt->groupsel[i] : NULL);
2441             if (maxtmp > opt->max)
2442             {
2443                 opt->max = maxtmp;
2444             }
2445             if (mintmp < opt->min)
2446             {
2447                 opt->min = mintmp;
2448             }
2449         }
2450         printf("\nDetermined boundaries to %f and %f\n\n", opt->min, opt->max);
2451         if (opt->bBoundsOnly)
2452         {
2453             printf("Found option -boundsonly, now exiting.\n");
2454             exit (0);
2455         }
2456     }
2457     /* store stepsize in profile */
2458     opt->dz = (opt->max-opt->min)/opt->bins;
2459
2460     for (i = 0; i < nfiles; i++)
2461     {
2462         if (whaminFileType(fnTprs[i]) != whamin_tpr)
2463         {
2464             gmx_fatal(FARGS, "Expected the %d'th file in input file to be a tpr file\n", i);
2465         }
2466         read_tpr_header(fnTprs[i], header, opt);
2467         if (whaminFileType(fnPull[i]) != whamin_pullxf)
2468         {
2469             gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
2470         }
2471         read_pull_xf(fnPull[i], fnTprs[i], header, window+i, opt, FALSE, NULL, NULL,
2472                      (opt->nGroupsel > 0) ? &opt->groupsel[i] : NULL);
2473         if (window[i].Ntot[0] == 0)
2474         {
2475             fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fnPull[i]);
2476         }
2477     }
2478
2479     for (i = 0; i < nfiles; i++)
2480     {
2481         sfree(fnTprs[i]);
2482         sfree(fnPull[i]);
2483     }
2484     sfree(fnTprs);
2485     sfree(fnPull);
2486 }
2487
2488 /*! \brief Read integrated autocorrelation times from input file (option -iiact)
2489  *
2490  * Note: Here we consider tau[int] := int_0^inf ACF(t) as the integrated autocorrelation times.
2491  * The factor `g := 1 + 2*tau[int]` subsequently enters the uncertainty.
2492  */
2493 void readIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins, const char* fn)
2494 {
2495     int      nlines, ny, i, ig;
2496     double **iact;
2497
2498     printf("Readging Integrated autocorrelation times from %s ...\n", fn);
2499     nlines = read_xvg(fn, &iact, &ny);
2500     if (nlines != nwins)
2501     {
2502         gmx_fatal(FARGS, "Found %d lines with integrated autocorrelation times in %s.\nExpected %d",
2503                   nlines, fn, nwins);
2504     }
2505     for (i = 0; i < nlines; i++)
2506     {
2507         if (window[i].nPull != ny)
2508         {
2509             gmx_fatal(FARGS, "You are providing autocorrelation times with option -iiact and the\n"
2510                       "number of pull groups is different in different simulations. That is not\n"
2511                       "supported yet. Sorry.\n");
2512         }
2513         for (ig = 0; ig < window[i].nPull; ig++)
2514         {
2515             /* compare Kumar et al, J Comp Chem 13, 1011-1021 (1992) */
2516             window[i].g[ig] = 1+2*iact[ig][i]/window[i].dt;
2517
2518             if (iact[ig][i] <= 0.0)
2519             {
2520                 fprintf(stderr, "\nWARNING, IACT = %f (window %d, group %d)\n", iact[ig][i], i, ig);
2521             }
2522         }
2523     }
2524 }
2525
2526
2527 /*! \brief Smooth autocorreltion times along the reaction coordinate.
2528  *
2529  * This is useful
2530  * if the ACT is subject to high uncertainty in case if limited sampling. Note
2531  * that -in case of limited sampling- the ACT may be severely underestimated.
2532  * Note: the g=1+2tau are overwritten.
2533  * If opt->bAllowReduceIact==FALSE, the ACTs are never reduced, only increased
2534  * by the smoothing
2535  */
2536 void smoothIact(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt)
2537 {
2538     int    i, ig, j, jg;
2539     double pos, dpos2, siglim, siglim2, gaufact, invtwosig2, w, weight, tausm;
2540
2541     /* only evaluate within +- 3sigma of the Gausian */
2542     siglim  = 3.0*opt->sigSmoothIact;
2543     siglim2 = dsqr(siglim);
2544     /* pre-factor of Gaussian */
2545     gaufact    = 1.0/(sqrt(2*M_PI)*opt->sigSmoothIact);
2546     invtwosig2 = 0.5/dsqr(opt->sigSmoothIact);
2547
2548     for (i = 0; i < nwins; i++)
2549     {
2550         snew(window[i].tausmooth, window[i].nPull);
2551         for (ig = 0; ig < window[i].nPull; ig++)
2552         {
2553             tausm  = 0.;
2554             weight = 0;
2555             pos    = window[i].pos[ig];
2556             for (j = 0; j < nwins; j++)
2557             {
2558                 for (jg = 0; jg < window[j].nPull; jg++)
2559                 {
2560                     dpos2 = dsqr(window[j].pos[jg]-pos);
2561                     if (dpos2 < siglim2)
2562                     {
2563                         w       = gaufact*exp(-dpos2*invtwosig2);
2564                         weight += w;
2565                         tausm  += w*window[j].tau[jg];
2566                         /*printf("Weight %g dpos2=%g pos=%g gaufact=%g invtwosig2=%g\n",
2567                            w,dpos2,pos,gaufact,invtwosig2); */
2568                     }
2569                 }
2570             }
2571             tausm /= weight;
2572             if (opt->bAllowReduceIact || tausm > window[i].tau[ig])
2573             {
2574                 window[i].tausmooth[ig] = tausm;
2575             }
2576             else
2577             {
2578                 window[i].tausmooth[ig] = window[i].tau[ig];
2579             }
2580             window[i].g[ig] = 1+2*tausm/window[i].dt;
2581         }
2582     }
2583 }
2584
2585 //! Stop integrating autoccorelation function when ACF drops under this value
2586 #define WHAM_AC_ZERO_LIMIT 0.05
2587
2588 /*! \brief Try to compute the autocorrelation time for each umbrealla window
2589  */
2590 void calcIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins,
2591                                         t_UmbrellaOptions *opt, const char *fn)
2592 {
2593     int   i, ig, ncorr, ntot, j, k, *count, restart;
2594     real *corr, c0, dt, tmp;
2595     real *ztime, av, tausteps;
2596     FILE *fp, *fpcorr = 0;
2597
2598     if (opt->verbose)
2599     {
2600         fpcorr = xvgropen("hist_autocorr.xvg", "Autocorrelation functions of umbrella windows",
2601                           "time [ps]", "autocorrelation function", opt->oenv);
2602     }
2603
2604     printf("\n");
2605     for (i = 0; i < nwins; i++)
2606     {
2607         printf("\rEstimating integrated autocorreltion times ... [%2.0f%%] ...", 100.*(i+1)/nwins);
2608         fflush(stdout);
2609         ntot = window[i].Ntot[0];
2610
2611         /* using half the maximum time as length of autocorrelation function */
2612         ncorr = ntot/2;
2613         if (ntot < 10)
2614         {
2615             gmx_fatal(FARGS, "Tryig to estimtate autocorrelation time from only %d"
2616                       " points. Provide more pull data!", ntot);
2617         }
2618         snew(corr, ncorr);
2619         /* snew(corrSq,ncorr); */
2620         snew(count, ncorr);
2621         dt = window[i].dt;
2622         snew(window[i].tau, window[i].nPull);
2623         restart = static_cast<int>(opt->acTrestart/dt+0.5);
2624         if (restart == 0)
2625         {
2626             restart = 1;
2627         }
2628
2629         for (ig = 0; ig < window[i].nPull; ig++)
2630         {
2631             if (ntot != window[i].Ntot[ig])
2632             {
2633                 gmx_fatal(FARGS, "Encountered different nr of frames in different pull groups.\n"
2634                           "That should not happen. (%d and %d)\n", ntot, window[i].Ntot[ig]);
2635             }
2636             ztime = window[i].ztime[ig];
2637
2638             /* calc autocorrelation function C(t) = < [z(tau)-<z>]*[z(tau+t)-<z>]> */
2639             for (j = 0, av = 0; (j < ntot); j++)
2640             {
2641                 av += ztime[j];
2642             }
2643             av /= ntot;
2644             for (k = 0; (k < ncorr); k++)
2645             {
2646                 corr[k]  = 0.;
2647                 count[k] = 0;
2648             }
2649             for (j = 0; (j < ntot); j += restart)
2650             {
2651                 for (k = 0; (k < ncorr) && (j+k < ntot); k++)
2652                 {
2653                     tmp        = (ztime[j]-av)*(ztime[j+k]-av);
2654                     corr  [k] += tmp;
2655                     /* corrSq[k] += tmp*tmp; */
2656                     count[k]++;
2657                 }
2658             }
2659             /* divide by nr of frames for each time displacement */
2660             for (k = 0; (k < ncorr); k++)
2661             {
2662                 /* count probably = (ncorr-k+(restart-1))/restart; */
2663                 corr[k] = corr[k]/count[k];
2664                 /* variance of autocorrelation function */
2665                 /* corrSq[k]=corrSq[k]/count[k]; */
2666             }
2667             /* normalize such that corr[0] == 0 */
2668             c0 = 1./corr[0];
2669             for (k = 0; (k < ncorr); k++)
2670             {
2671                 corr[k] *= c0;
2672                 /* corrSq[k]*=c0*c0; */
2673             }
2674
2675             /* write ACFs in verbose mode */
2676             if (fpcorr)
2677             {
2678                 for (k = 0; (k < ncorr); k++)
2679                 {
2680                     fprintf(fpcorr, "%g  %g\n", k*dt, corr[k]);
2681                 }
2682                 fprintf(fpcorr, "&\n");
2683             }
2684
2685             /* esimate integrated correlation time, fitting is too unstable */
2686             tausteps = 0.5*corr[0];
2687             /* consider corr below WHAM_AC_ZERO_LIMIT as noise */
2688             for (j = 1; (j < ncorr) && (corr[j] > WHAM_AC_ZERO_LIMIT); j++)
2689             {
2690                 tausteps += corr[j];
2691             }
2692
2693             /* g = 1+2*tau, see. Ferrenberg/Swendsen, PRL 63:1195 (1989) or
2694                Kumar et al, eq. 28 ff. */
2695             window[i].tau[ig] = tausteps*dt;
2696             window[i].g[ig]   = 1+2*tausteps;
2697             /* printf("win %d, group %d, estimated correlation time = %g ps\n",i,ig,window[i].tau[ig]); */
2698         } /* ig loop */
2699         sfree(corr);
2700         sfree(count);
2701     }
2702     printf(" done\n");
2703     if (fpcorr)
2704     {
2705         ffclose(fpcorr);
2706     }
2707
2708     /* plot IACT along reaction coordinate */
2709     fp = xvgropen(fn, "Integrated autocorrelation times", "z", "IACT [ps]", opt->oenv);
2710     fprintf(fp, "@    s0 symbol 1\n@    s0 symbol size 0.5\n@    s0 line linestyle 0\n");
2711     fprintf(fp, "#  WIN   tau(gr1)  tau(gr2) ...\n");
2712     for (i = 0; i < nwins; i++)
2713     {
2714         fprintf(fp, "# %3d   ", i);
2715         for (ig = 0; ig < window[i].nPull; ig++)
2716         {
2717             fprintf(fp, " %11g", window[i].tau[ig]);
2718         }
2719         fprintf(fp, "\n");
2720     }
2721     for (i = 0; i < nwins; i++)
2722     {
2723         for (ig = 0; ig < window[i].nPull; ig++)
2724         {
2725             fprintf(fp, "%8g %8g\n", window[i].pos[ig], window[i].tau[ig]);
2726         }
2727     }
2728     if (opt->sigSmoothIact > 0.0)
2729     {
2730         printf("Smoothing autocorrelation times along reaction coordinate with Gaussian of sig = %g\n",
2731                opt->sigSmoothIact);
2732         /* smooth IACT along reaction coordinate and overwrite g=1+2tau */
2733         smoothIact(window, nwins, opt);
2734         fprintf(fp, "&\n@    s1 symbol 1\n@    s1 symbol size 0.5\n@    s1 line linestyle 0\n");
2735         fprintf(fp, "@    s1 symbol color 2\n");
2736         for (i = 0; i < nwins; i++)
2737         {
2738             for (ig = 0; ig < window[i].nPull; ig++)
2739             {
2740                 fprintf(fp, "%8g %8g\n", window[i].pos[ig], window[i].tausmooth[ig]);
2741             }
2742         }
2743     }
2744     ffclose(fp);
2745     printf("Wrote %s\n", fn);
2746 }
2747
2748 /*! \brief
2749  * compute average and sigma of each umbrella histogram
2750  */
2751 void averageSigma(t_UmbrellaWindow *window, int nwins)
2752 {
2753     int  i, ig, ntot, k;
2754     real av, sum2, sig, diff, *ztime, nSamplesIndep;
2755
2756     for (i = 0; i < nwins; i++)
2757     {
2758         snew(window[i].aver, window[i].nPull);
2759         snew(window[i].sigma, window[i].nPull);
2760
2761         ntot = window[i].Ntot[0];
2762         for (ig = 0; ig < window[i].nPull; ig++)
2763         {
2764             ztime = window[i].ztime[ig];
2765             for (k = 0, av = 0.; k < ntot; k++)
2766             {
2767                 av += ztime[k];
2768             }
2769             av /= ntot;
2770             for (k = 0, sum2 = 0.; k < ntot; k++)
2771             {
2772                 diff  = ztime[k]-av;
2773                 sum2 += diff*diff;
2774             }
2775             sig                = sqrt(sum2/ntot);
2776             window[i].aver[ig] = av;
2777
2778             /* Note: This estimate for sigma is biased from the limited sampling.
2779                Correct sigma by n/(n-1) where n = number of independent
2780                samples. Only possible if IACT is known.
2781              */
2782             if (window[i].tau)
2783             {
2784                 nSamplesIndep       = window[i].N[ig]/(window[i].tau[ig]/window[i].dt);
2785                 window[i].sigma[ig] = sig * nSamplesIndep/(nSamplesIndep-1);
2786             }
2787             else
2788             {
2789                 window[i].sigma[ig] = sig;
2790             }
2791             printf("win %d, aver = %f  sig = %f\n", i, av, window[i].sigma[ig]);
2792         }
2793     }
2794 }
2795
2796
2797 /*! \brief
2798  * Use histograms to compute average force on pull group.
2799  */
2800 void computeAverageForce(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt)
2801 {
2802     int    i, j, bins = opt->bins, k;
2803     double dz, min = opt->min, max = opt->max, displAv, temp, distance, ztot, ztot_half, w, weight;
2804     double posmirrored;
2805
2806     dz        = (max-min)/bins;
2807     ztot      = opt->max-min;
2808     ztot_half = ztot/2;
2809
2810     /* Compute average displacement from histograms */
2811     for (j = 0; j < nWindows; ++j)
2812     {
2813         snew(window[j].forceAv, window[j].nPull);
2814         for (k = 0; k < window[j].nPull; ++k)
2815         {
2816             displAv = 0.0;
2817             weight  = 0.0;
2818             for (i = 0; i < opt->bins; ++i)
2819             {
2820                 temp     = (1.0*i+0.5)*dz+min;
2821                 distance = temp - window[j].pos[k];
2822                 if (opt->bCycl)
2823                 {                                       /* in cyclic wham:             */
2824                     if (distance > ztot_half)           /*    |distance| < ztot_half   */
2825                     {
2826                         distance -= ztot;
2827                     }
2828                     else if (distance < -ztot_half)
2829                     {
2830                         distance += ztot;
2831                     }
2832                 }
2833                 w         = window[j].Histo[k][i]/window[j].g[k];
2834                 displAv  += w*distance;
2835                 weight   += w;
2836                 /* Are we near min or max? We are getting wrong forces from the histgrams since
2837                    the histograms are zero outside [min,max). Therefore, assume that the position
2838                    on the other side of the histomgram center is equally likely. */
2839                 if (!opt->bCycl)
2840                 {
2841                     posmirrored = window[j].pos[k]-distance;
2842                     if (posmirrored >= max || posmirrored < min)
2843                     {
2844                         displAv  += -w*distance;
2845                         weight   += w;
2846                     }
2847                 }
2848             }
2849             displAv  /= weight;
2850
2851             /* average force from average displacement */
2852             window[j].forceAv[k] = displAv*window[j].k[k];
2853             /* sigma from average square displacement */
2854             /* window[j].sigma  [k] = sqrt(displAv2); */
2855             /* printf("Win %d, sigma = %f\n",j,sqrt(displAv2)); */
2856         }
2857     }
2858 }
2859
2860 /*! \brief
2861  * Check if the complete reaction coordinate is covered by the histograms
2862  */
2863 void  checkReactionCoordinateCovered(t_UmbrellaWindow *window, int nwins,
2864                                      t_UmbrellaOptions *opt)
2865 {
2866     int  i, ig, j, bins = opt->bins, bBoundary;
2867     real avcount = 0, z, relcount, *count;
2868     snew(count, opt->bins);
2869
2870     for (j = 0; j < opt->bins; ++j)
2871     {
2872         for (i = 0; i < nwins; i++)
2873         {
2874             for (ig = 0; ig < window[i].nPull; ig++)
2875             {
2876                 count[j] += window[i].Histo[ig][j];
2877             }
2878         }
2879         avcount += 1.0*count[j];
2880     }
2881     avcount /= bins;
2882     for (j = 0; j < bins; ++j)
2883     {
2884         relcount  = count[j]/avcount;
2885         z         = (j+0.5)*opt->dz+opt->min;
2886         bBoundary = ( j<bins/20 || (bins-j)>bins/20 );
2887         /* check for bins with no data */
2888         if (count[j] == 0)
2889         {
2890             fprintf(stderr, "\nWARNING, no data point in bin %d (z=%g) !\n"
2891                     "You may not get a reasonable profile. Check your histograms!\n", j, z);
2892         }
2893         /* and check for poor sampling */
2894         else if (relcount < 0.005 && !bBoundary)
2895         {
2896             fprintf(stderr, "Warning, poor sampling bin %d (z=%g). Check your histograms!\n", j, z);
2897         }
2898     }
2899     sfree(count);
2900 }
2901
2902 /*! \brief Compute initial potential by integrating the average force
2903  *
2904  * This speeds up the convergence by roughly a factor of 2
2905  */
2906 void guessPotByIntegration(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt)
2907 {
2908     int    i, j, ig, bins = opt->bins, nHist, winmin, groupmin;
2909     double dz, min = opt->min, *pot, pos, hispos, dist, diff, fAv, distmin, *f;
2910     FILE  *fp;
2911
2912     dz = (opt->max-min)/bins;
2913
2914     printf("Getting initial potential by integration.\n");
2915
2916     /* Compute average displacement from histograms */
2917     computeAverageForce(window, nWindows, opt);
2918
2919     /* Get force for each bin from all histograms in this bin, or, alternatively,
2920        if no histograms are inside this bin, from the closest histogram */
2921     snew(pot, bins);
2922     snew(f, bins);
2923     for (j = 0; j < opt->bins; ++j)
2924     {
2925         pos      = (1.0*j+0.5)*dz+min;
2926         nHist    = 0;
2927         fAv      = 0.;
2928         distmin  = 1e20;
2929         groupmin = winmin = 0;
2930         for (i = 0; i < nWindows; i++)
2931         {
2932             for (ig = 0; ig < window[i].nPull; ig++)
2933             {
2934                 hispos = window[i].pos[ig];
2935                 dist   = fabs(hispos-pos);
2936                 /* average force within bin */
2937                 if (dist < dz/2)
2938                 {
2939                     nHist++;
2940                     fAv += window[i].forceAv[ig];
2941                 }
2942                 /* at the same time, rememer closest histogram */
2943                 if (dist < distmin)
2944                 {
2945                     winmin   = i;
2946                     groupmin = ig;
2947                     distmin  = dist;
2948                 }
2949             }
2950         }
2951         /* if no histogram found in this bin, use closest histogram */
2952         if (nHist > 0)
2953         {
2954             fAv = fAv/nHist;
2955         }
2956         else
2957         {
2958             fAv = window[winmin].forceAv[groupmin];
2959         }
2960         f[j] = fAv;
2961     }
2962     for (j = 1; j < opt->bins; ++j)
2963     {
2964         pot[j] = pot[j-1] - 0.5*dz*(f[j-1]+f[j]);
2965     }
2966
2967     /* cyclic wham: linearly correct possible offset */
2968     if (opt->bCycl)
2969     {
2970         diff = (pot[bins-1]-pot[0])/(bins-1);
2971         for (j = 1; j < opt->bins; ++j)
2972         {
2973             pot[j] -= j*diff;
2974         }
2975     }
2976     if (opt->verbose)
2977     {
2978         fp = xvgropen("pmfintegrated.xvg", "PMF from force integration", "z", "PMF [kJ/mol]", opt->oenv);
2979         for (j = 0; j < opt->bins; ++j)
2980         {
2981             fprintf(fp, "%g  %g\n", (j+0.5)*dz+opt->min, pot[j]);
2982         }
2983         ffclose(fp);
2984         printf("verbose mode: wrote %s with PMF from interated forces\n", "pmfintegrated.xvg");
2985     }
2986
2987     /* get initial z=exp(-F[i]/kT) from integrated potential, where F[i] denote the free
2988        energy offsets which are usually determined by wham
2989        First: turn pot into probabilities:
2990      */
2991     for (j = 0; j < opt->bins; ++j)
2992     {
2993         pot[j] = exp(-pot[j]/(8.314e-3*opt->Temperature));
2994     }
2995     calc_z(pot, window, nWindows, opt, TRUE);
2996
2997     sfree(pot);
2998     sfree(f);
2999 }
3000
3001 //! Count number of words in an line
3002 static int wordcount(char *ptr)
3003 {
3004     int i, n, is[2];
3005     int cur = 0;
3006
3007     if (strlen(ptr) == 0)
3008     {
3009         return 0;
3010     }
3011     /* fprintf(stderr,"ptr='%s'\n",ptr); */
3012     n = 1;
3013     for (i = 0; (ptr[i] != '\0'); i++)
3014     {
3015         is[cur] = isspace(ptr[i]);
3016         if ((i > 0)  && (is[cur] && !is[1-cur]))
3017         {
3018             n++;
3019         }
3020         cur = 1-cur;
3021     }
3022     return n;
3023 }
3024
3025 /*! \brief Read input file for pull group selection (option -is)
3026  *
3027  * TO DO: ptr=fgets(...) is never freed (small memory leak)
3028  */
3029 void readPullGroupSelection(t_UmbrellaOptions *opt, char **fnTpr, int nTpr)
3030 {
3031     FILE *fp;
3032     int   i, iline, n, len = STRLEN, temp;
3033     char *ptr = 0, *tmpbuf = 0;
3034     char  fmt[1024], fmtign[1024];
3035     int   block = 1, sizenow;
3036
3037     fp            = ffopen(opt->fnGroupsel, "r");
3038     opt->groupsel = NULL;
3039
3040     snew(tmpbuf, len);
3041     sizenow = 0;
3042     iline   = 0;
3043     while ( (ptr = fgets3(fp, tmpbuf, &len)) != NULL)
3044     {
3045         trim(ptr);
3046         n = wordcount(ptr);
3047
3048         if (iline >= sizenow)
3049         {
3050             sizenow += block;
3051             srenew(opt->groupsel, sizenow);
3052         }
3053         opt->groupsel[iline].n    = n;
3054         opt->groupsel[iline].nUse = 0;
3055         snew(opt->groupsel[iline].bUse, n);
3056
3057         fmtign[0] = '\0';
3058         for (i = 0; i < n; i++)
3059         {
3060             strcpy(fmt, fmtign);
3061             strcat(fmt, "%d");
3062             if (sscanf(ptr, fmt, &temp))
3063             {
3064                 opt->groupsel[iline].bUse[i] = (temp > 0);
3065                 if (opt->groupsel[iline].bUse[i])
3066                 {
3067                     opt->groupsel[iline].nUse++;
3068                 }
3069             }
3070             strcat(fmtign, "%*s");
3071         }
3072         iline++;
3073     }
3074     opt->nGroupsel = iline;
3075     if (nTpr != opt->nGroupsel)
3076     {
3077         gmx_fatal(FARGS, "Found %d tpr files but %d lines in %s\n", nTpr, opt->nGroupsel,
3078                   opt->fnGroupsel);
3079     }
3080
3081     printf("\nUse only these pull groups:\n");
3082     for (iline = 0; iline < nTpr; iline++)
3083     {
3084         printf("%s (%d of %d groups):", fnTpr[iline], opt->groupsel[iline].nUse, opt->groupsel[iline].n);
3085         for (i = 0; i < opt->groupsel[iline].n; i++)
3086         {
3087             if (opt->groupsel[iline].bUse[i])
3088             {
3089                 printf(" %d", i+1);
3090             }
3091         }
3092         printf("\n");
3093     }
3094     printf("\n");
3095
3096     sfree(tmpbuf);
3097 }
3098
3099 //! Boolean XOR
3100 #define WHAMBOOLXOR(a, b) ( ((!(a)) && (b)) || ((a) && (!(b))))
3101
3102 /*! Number of elements in fnm (used for command line parsing) */
3103 #define NFILE asize(fnm)
3104
3105 //! The main g_wham routine
3106 int gmx_wham(int argc, char *argv[])
3107 {
3108     const char              *desc[] = {
3109         "[THISMODULE] is an analysis program that implements the Weighted",
3110         "Histogram Analysis Method (WHAM). It is intended to analyze",
3111         "output files generated by umbrella sampling simulations to ",
3112         "compute a potential of mean force (PMF). [PAR] ",
3113         "At present, three input modes are supported.[BR]",
3114         "[TT]*[tt] With option [TT]-it[tt], the user provides a file which contains the",
3115         " file names of the umbrella simulation run-input files ([TT].tpr[tt] files),",
3116         " AND, with option [TT]-ix[tt], a file which contains file names of",
3117         " the pullx [TT]mdrun[tt] output files. The [TT].tpr[tt] and pullx files must",
3118         " be in corresponding order, i.e. the first [TT].tpr[tt] created the",
3119         " first pullx, etc.[BR]",
3120         "[TT]*[tt] Same as the previous input mode, except that the the user",
3121         " provides the pull force output file names ([TT]pullf.xvg[tt]) with option [TT]-if[tt].",
3122         " From the pull force the position in the umbrella potential is",
3123         " computed. This does not work with tabulated umbrella potentials.[BR]"
3124         "[TT]*[tt] With option [TT]-ip[tt], the user provides file names of (gzipped) [TT].pdo[tt] files, i.e.",
3125         " the GROMACS 3.3 umbrella output files. If you have some unusual"
3126         " reaction coordinate you may also generate your own [TT].pdo[tt] files and",
3127         " feed them with the [TT]-ip[tt] option into to [THISMODULE]. The [TT].pdo[tt] file header",
3128         " must be similar to the following:[PAR]",
3129         "[TT]# UMBRELLA      3.0[BR]",
3130         "# Component selection: 0 0 1[BR]",
3131         "# nSkip 1[BR]",
3132         "# Ref. Group 'TestAtom'[BR]",
3133         "# Nr. of pull groups 2[BR]",
3134         "# Group 1 'GR1'  Umb. Pos. 5.0 Umb. Cons. 1000.0[BR]",
3135         "# Group 2 'GR2'  Umb. Pos. 2.0 Umb. Cons. 500.0[BR]",
3136         "#####[tt][PAR]",
3137         "The number of pull groups, umbrella positions, force constants, and names ",
3138         "may (of course) differ. Following the header, a time column and ",
3139         "a data column for each pull group follows (i.e. the displacement",
3140         "with respect to the umbrella center). Up to four pull groups are possible ",
3141         "per [TT].pdo[tt] file at present.[PAR]",
3142         "By default, all pull groups found in all pullx/pullf files are used in WHAM. If only ",
3143         "some of the pull groups should be used, a pull group selection file (option [TT]-is[tt]) can ",
3144         "be provided. The selection file must contain one line for each tpr file in tpr-files.dat.",
3145         "Each of these lines must contain one digit (0 or 1) for each pull group in the tpr file. ",
3146         "Here, 1 indicates that the pull group is used in WHAM, and 0 means it is omitted. Example:",
3147         "If you have three tpr files, each containing 4 pull groups, but only pull group 1 and 2 should be ",
3148         "used, groupsel.dat looks like this:[BR]",
3149         "1 1 0 0[BR]",
3150         "1 1 0 0[BR]",
3151         "1 1 0 0[PAR]",
3152         "By default, the output files are[BR]",
3153         "  [TT]-o[tt]      PMF output file[BR]",
3154         "  [TT]-hist[tt]   Histograms output file[BR]",
3155         "Always check whether the histograms sufficiently overlap.[PAR]",
3156         "The umbrella potential is assumed to be harmonic and the force constants are ",
3157         "read from the [TT].tpr[tt] or [TT].pdo[tt] files. If a non-harmonic umbrella force was applied ",
3158         "a tabulated potential can be provided with [TT]-tab[tt].[PAR]",
3159         "WHAM OPTIONS[BR]------------[BR]",
3160         "  [TT]-bins[tt]   Number of bins used in analysis[BR]",
3161         "  [TT]-temp[tt]   Temperature in the simulations[BR]",
3162         "  [TT]-tol[tt]    Stop iteration if profile (probability) changed less than tolerance[BR]",
3163         "  [TT]-auto[tt]   Automatic determination of boundaries[BR]",
3164         "  [TT]-min,-max[tt]   Boundaries of the profile [BR]",
3165         "The data points that are used to compute the profile",
3166         "can be restricted with options [TT]-b[tt], [TT]-e[tt], and [TT]-dt[tt]. ",
3167         "Adjust [TT]-b[tt] to ensure sufficient equilibration in each ",
3168         "umbrella window.[PAR]",
3169         "With [TT]-log[tt] (default) the profile is written in energy units, otherwise ",
3170         "(with [TT]-nolog[tt]) as probability. The unit can be specified with [TT]-unit[tt]. ",
3171         "With energy output, the energy in the first bin is defined to be zero. ",
3172         "If you want the free energy at a different ",
3173         "position to be zero, set [TT]-zprof0[tt] (useful with bootstrapping, see below).[PAR]",
3174         "For cyclic or periodic reaction coordinates (dihedral angle, channel PMF",
3175         "without osmotic gradient), the option [TT]-cycl[tt] is useful.",
3176         "[THISMODULE] will make use of the",
3177         "periodicity of the system and generate a periodic PMF. The first and the last bin of the",
3178         "reaction coordinate will assumed be be neighbors.[PAR]",
3179         "Option [TT]-sym[tt] symmetrizes the profile around z=0 before output, ",
3180         "which may be useful for, e.g. membranes.[PAR]",
3181         "AUTOCORRELATIONS[BR]----------------[BR]",
3182         "With [TT]-ac[tt], [THISMODULE] estimates the integrated autocorrelation ",
3183         "time (IACT) [GRK]tau[grk] for each umbrella window and weights the respective ",
3184         "window with 1/[1+2*[GRK]tau[grk]/dt]. The IACTs are written ",
3185         "to the file defined with [TT]-oiact[tt]. In verbose mode, all ",
3186         "autocorrelation functions (ACFs) are written to [TT]hist_autocorr.xvg[tt]. ",
3187         "Because the IACTs can be severely underestimated in case of limited ",
3188         "sampling, option [TT]-acsig[tt] allows one to smooth the IACTs along the ",
3189         "reaction coordinate with a Gaussian ([GRK]sigma[grk] provided with [TT]-acsig[tt], ",
3190         "see output in [TT]iact.xvg[tt]). Note that the IACTs are estimated by simple ",
3191         "integration of the ACFs while the ACFs are larger 0.05.",
3192         "If you prefer to compute the IACTs by a more sophisticated (but possibly ",
3193         "less robust) method such as fitting to a double exponential, you can ",
3194         "compute the IACTs with [gmx-analyze] and provide them to [THISMODULE] with the file ",
3195         "[TT]iact-in.dat[tt] (option [TT]-iiact[tt]), which should contain one line per ",
3196         "input file ([TT].pdo[tt] or pullx/f file) and one column per pull group in the respective file.[PAR]",
3197         "ERROR ANALYSIS[BR]--------------[BR]",
3198         "Statistical errors may be estimated with bootstrap analysis. Use it with care, ",
3199         "otherwise the statistical error may be substantially underestimated. ",
3200         "More background and examples for the bootstrap technique can be found in ",
3201         "Hub, de Groot and Van der Spoel, JCTC (2010) 6: 3713-3720.[BR]",
3202         "[TT]-nBootstrap[tt] defines the number of bootstraps (use, e.g., 100). ",
3203         "Four bootstrapping methods are supported and ",
3204         "selected with [TT]-bs-method[tt].[BR]",
3205         "  (1) [TT]b-hist[tt]   Default: complete histograms are considered as independent ",
3206         "data points, and the bootstrap is carried out by assigning random weights to the ",
3207         "histograms (\"Bayesian bootstrap\"). Note that each point along the reaction coordinate",
3208         "must be covered by multiple independent histograms (e.g. 10 histograms), otherwise the ",
3209         "statistical error is underestimated.[BR]",
3210         "  (2) [TT]hist[tt]    Complete histograms are considered as independent data points. ",
3211         "For each bootstrap, N histograms are randomly chosen from the N given histograms ",
3212         "(allowing duplication, i.e. sampling with replacement).",
3213         "To avoid gaps without data along the reaction coordinate blocks of histograms ",
3214         "([TT]-histbs-block[tt]) may be defined. In that case, the given histograms are ",
3215         "divided into blocks and only histograms within each block are mixed. Note that ",
3216         "the histograms within each block must be representative for all possible histograms, ",
3217         "otherwise the statistical error is underestimated.[BR]",
3218         "  (3) [TT]traj[tt]  The given histograms are used to generate new random trajectories,",
3219         "such that the generated data points are distributed according the given histograms ",
3220         "and properly autocorrelated. The autocorrelation time (ACT) for each window must be ",
3221         "known, so use [TT]-ac[tt] or provide the ACT with [TT]-iiact[tt]. If the ACT of all ",
3222         "windows are identical (and known), you can also provide them with [TT]-bs-tau[tt]. ",
3223         "Note that this method may severely underestimate the error in case of limited sampling, ",
3224         "that is if individual histograms do not represent the complete phase space at ",
3225         "the respective positions.[BR]",
3226         "  (4) [TT]traj-gauss[tt]  The same as method [TT]traj[tt], but the trajectories are ",
3227         "not bootstrapped from the umbrella histograms but from Gaussians with the average ",
3228         "and width of the umbrella histograms. That method yields similar error estimates ",
3229         "like method [TT]traj[tt].[PAR]"
3230         "Bootstrapping output:[BR]",
3231         "  [TT]-bsres[tt]   Average profile and standard deviations[BR]",
3232         "  [TT]-bsprof[tt]  All bootstrapping profiles[BR]",
3233         "With [TT]-vbs[tt] (verbose bootstrapping), the histograms of each bootstrap are written, ",
3234         "and, with bootstrap method [TT]traj[tt], the cumulative distribution functions of ",
3235         "the histograms."
3236     };
3237
3238     const char              *en_unit[]       = {NULL, "kJ", "kCal", "kT", NULL};
3239     const char              *en_unit_label[] = {"", "E (kJ mol\\S-1\\N)", "E (kcal mol\\S-1\\N)", "E (kT)", NULL};
3240     const char              *en_bsMethod[]   = { NULL, "b-hist", "hist", "traj", "traj-gauss", NULL };
3241
3242     static t_UmbrellaOptions opt;
3243
3244     t_pargs                  pa[] = {
3245         { "-min", FALSE, etREAL, {&opt.min},
3246           "Minimum coordinate in profile"},
3247         { "-max", FALSE, etREAL, {&opt.max},
3248           "Maximum coordinate in profile"},
3249         { "-auto", FALSE, etBOOL, {&opt.bAuto},
3250           "Determine min and max automatically"},
3251         { "-bins", FALSE, etINT, {&opt.bins},
3252           "Number of bins in profile"},
3253         { "-temp", FALSE, etREAL, {&opt.Temperature},
3254           "Temperature"},
3255         { "-tol", FALSE, etREAL, {&opt.Tolerance},
3256           "Tolerance"},
3257         { "-v", FALSE, etBOOL, {&opt.verbose},
3258           "Verbose mode"},
3259         { "-b", FALSE, etREAL, {&opt.tmin},
3260           "First time to analyse (ps)"},
3261         { "-e", FALSE, etREAL, {&opt.tmax},
3262           "Last time to analyse (ps)"},
3263         { "-dt", FALSE, etREAL, {&opt.dt},
3264           "Analyse only every dt ps"},
3265         { "-histonly", FALSE, etBOOL, {&opt.bHistOnly},
3266           "Write histograms and exit"},
3267         { "-boundsonly", FALSE, etBOOL, {&opt.bBoundsOnly},
3268           "Determine min and max and exit (with [TT]-auto[tt])"},
3269         { "-log", FALSE, etBOOL, {&opt.bLog},
3270           "Calculate the log of the profile before printing"},
3271         { "-unit", FALSE,  etENUM, {en_unit},
3272           "Energy unit in case of log output" },
3273         { "-zprof0", FALSE, etREAL, {&opt.zProf0},
3274           "Define profile to 0.0 at this position (with [TT]-log[tt])"},
3275         { "-cycl", FALSE, etBOOL, {&opt.bCycl},
3276           "Create cyclic/periodic profile. Assumes min and max are the same point."},
3277         { "-sym", FALSE, etBOOL, {&opt.bSym},
3278           "Symmetrize profile around z=0"},
3279         { "-hist-eq", FALSE, etBOOL, {&opt.bHistEq},
3280           "HIDDENEnforce equal weight for all histograms. (Non-Weighed-HAM)"},
3281         { "-ac", FALSE, etBOOL, {&opt.bCalcTauInt},
3282           "Calculate integrated autocorrelation times and use in wham"},
3283         { "-acsig", FALSE, etREAL, {&opt.sigSmoothIact},
3284           "Smooth autocorrelation times along reaction coordinate with Gaussian of this [GRK]sigma[grk]"},
3285         { "-ac-trestart", FALSE, etREAL, {&opt.acTrestart},
3286           "When computing autocorrelation functions, restart computing every .. (ps)"},
3287         { "-acred", FALSE, etBOOL, {&opt.bAllowReduceIact},
3288           "HIDDENWhen smoothing the ACTs, allow to reduce ACTs. Otherwise, only increase ACTs "
3289           "during smoothing"},
3290         { "-nBootstrap", FALSE,  etINT, {&opt.nBootStrap},
3291           "nr of bootstraps to estimate statistical uncertainty (e.g., 200)" },
3292         { "-bs-method", FALSE,  etENUM, {en_bsMethod},
3293           "Bootstrap method" },
3294         { "-bs-tau", FALSE, etREAL, {&opt.tauBootStrap},
3295           "Autocorrelation time (ACT) assumed for all histograms. Use option [TT]-ac[tt] if ACT is unknown."},
3296         { "-bs-seed", FALSE, etINT, {&opt.bsSeed},
3297           "Seed for bootstrapping. (-1 = use time)"},
3298         { "-histbs-block", FALSE, etINT, {&opt.histBootStrapBlockLength},
3299           "When mixing histograms only mix within blocks of [TT]-histbs-block[tt]."},
3300         { "-vbs", FALSE, etBOOL, {&opt.bs_verbose},
3301           "Verbose bootstrapping. Print the CDFs and a histogram file for each bootstrap."},
3302         { "-stepout", FALSE, etINT, {&opt.stepchange},
3303           "HIDDENWrite maximum change every ... (set to 1 with [TT]-v[tt])"},
3304         { "-updateContr", FALSE, etINT, {&opt.stepUpdateContrib},
3305           "HIDDENUpdate table with significan contributions to WHAM every ... iterations"},
3306     };
3307
3308     t_filenm                 fnm[] = {
3309         { efDAT, "-ix", "pullx-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs           */
3310         { efDAT, "-if", "pullf-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs           */
3311         { efDAT, "-it", "tpr-files", ffOPTRD},   /* wham input: tprs                           */
3312         { efDAT, "-ip", "pdo-files", ffOPTRD},   /* wham input: pdo files (gmx3 style)         */
3313         { efDAT, "-is", "groupsel", ffOPTRD},    /* input: select pull groups to use           */
3314         { efXVG, "-o", "profile", ffWRITE },     /* output file for profile                     */
3315         { efXVG, "-hist", "histo", ffWRITE},     /* output file for histograms                  */
3316         { efXVG, "-oiact", "iact", ffOPTWR},     /* writing integrated autocorrelation times    */
3317         { efDAT, "-iiact", "iact-in", ffOPTRD},  /* reading integrated autocorrelation times   */
3318         { efXVG, "-bsres", "bsResult", ffOPTWR}, /* average and errors of bootstrap analysis    */
3319         { efXVG, "-bsprof", "bsProfs", ffOPTWR}, /* output file for bootstrap profiles          */
3320         { efDAT, "-tab", "umb-pot", ffOPTRD},    /* Tabulated umbrella potential (if not harmonic) */
3321     };
3322
3323     int                      i, j, l, nfiles, nwins, nfiles2;
3324     t_UmbrellaHeader         header;
3325     t_UmbrellaWindow       * window = NULL;
3326     double                  *profile, maxchange = 1e20;
3327     gmx_bool                 bMinSet, bMaxSet, bAutoSet, bExact = FALSE;
3328     char                   **fninTpr, **fninPull, **fninPdo;
3329     const char              *fnPull;
3330     FILE                    *histout, *profout;
3331     char                     ylabel[256], title[256];
3332
3333     opt.bins      = 200;
3334     opt.verbose   = FALSE;
3335     opt.bHistOnly = FALSE;
3336     opt.bCycl     = FALSE;
3337     opt.tmin      = 50;
3338     opt.tmax      = 1e20;
3339     opt.dt        = 0.0;
3340     opt.min       = 0;
3341     opt.max       = 0;
3342     opt.bAuto     = TRUE;
3343     opt.nGroupsel = 0;
3344
3345     /* bootstrapping stuff */
3346     opt.nBootStrap               = 0;
3347     opt.bsMethod                 = bsMethod_hist;
3348     opt.tauBootStrap             = 0.0;
3349     opt.bsSeed                   = -1;
3350     opt.histBootStrapBlockLength = 8;
3351     opt.bs_verbose               = FALSE;
3352
3353     opt.bLog                  = TRUE;
3354     opt.unit                  = en_kJ;
3355     opt.zProf0                = 0.;
3356     opt.Temperature           = 298;
3357     opt.Tolerance             = 1e-6;
3358     opt.bBoundsOnly           = FALSE;
3359     opt.bSym                  = FALSE;
3360     opt.bCalcTauInt           = FALSE;
3361     opt.sigSmoothIact         = 0.0;
3362     opt.bAllowReduceIact      = TRUE;
3363     opt.bInitPotByIntegration = TRUE;
3364     opt.acTrestart            = 1.0;
3365     opt.stepchange            = 100;
3366     opt.stepUpdateContrib     = 100;
3367
3368     if (!parse_common_args(&argc, argv, PCA_BE_NICE,
3369                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &opt.oenv))
3370     {
3371         return 0;
3372     }
3373
3374     opt.unit     = nenum(en_unit);
3375     opt.bsMethod = nenum(en_bsMethod);
3376
3377     opt.bProf0Set = opt2parg_bSet("-zprof0",  asize(pa), pa);
3378
3379     opt.bTab         = opt2bSet("-tab", NFILE, fnm);
3380     opt.bPdo         = opt2bSet("-ip", NFILE, fnm);
3381     opt.bTpr         = opt2bSet("-it", NFILE, fnm);
3382     opt.bPullx       = opt2bSet("-ix", NFILE, fnm);
3383     opt.bPullf       = opt2bSet("-if", NFILE, fnm);
3384     opt.bTauIntGiven = opt2bSet("-iiact", NFILE, fnm);
3385     if  (opt.bTab && opt.bPullf)
3386     {
3387         gmx_fatal(FARGS, "Force input does not work with tabulated potentials. "
3388                   "Provide pullx.xvg or pdo files!");
3389     }
3390
3391     if (!opt.bPdo && !WHAMBOOLXOR(opt.bPullx, opt.bPullf))
3392     {
3393         gmx_fatal(FARGS, "Give either pullx (-ix) OR pullf (-if) data. Not both.");
3394     }
3395     if (!opt.bPdo && !(opt.bTpr || opt.bPullf || opt.bPullx))
3396     {
3397         gmx_fatal(FARGS, "g_wham supports three input modes, pullx, pullf, or pdo file input."
3398                   "\n\n Check g_wham -h !");
3399     }
3400
3401     opt.fnPdo      = opt2fn("-ip", NFILE, fnm);
3402     opt.fnTpr      = opt2fn("-it", NFILE, fnm);
3403     opt.fnPullf    = opt2fn("-if", NFILE, fnm);
3404     opt.fnPullx    = opt2fn("-ix", NFILE, fnm);
3405     opt.fnGroupsel = opt2fn_null("-is", NFILE, fnm);
3406
3407     bMinSet  = opt2parg_bSet("-min",  asize(pa), pa);
3408     bMaxSet  = opt2parg_bSet("-max",  asize(pa), pa);
3409     bAutoSet = opt2parg_bSet("-auto",  asize(pa), pa);
3410     if ( (bMinSet || bMaxSet) && bAutoSet)
3411     {
3412         gmx_fatal(FARGS, "With -auto, do not give -min or -max\n");
3413     }
3414
3415     if ( (bMinSet && !bMaxSet) || (!bMinSet && bMaxSet))
3416     {
3417         gmx_fatal(FARGS, "When giving -min, you must give -max (and vice versa), too\n");
3418     }
3419
3420     if (bMinSet && opt.bAuto)
3421     {
3422         printf("Note: min and max given, switching off -auto.\n");
3423         opt.bAuto = FALSE;
3424     }
3425
3426     if (opt.bTauIntGiven && opt.bCalcTauInt)
3427     {
3428         gmx_fatal(FARGS, "Either read (option -iiact) or calculate (option -ac) the\n"
3429                   "the autocorrelation times. Not both.");
3430     }
3431
3432     if (opt.tauBootStrap > 0.0 && opt2parg_bSet("-ac", asize(pa), pa))
3433     {
3434         gmx_fatal(FARGS, "Either compute autocorrelation times (ACTs) (option -ac) or "
3435                   "provide it with -bs-tau for bootstrapping. Not Both.\n");
3436     }
3437     if (opt.tauBootStrap > 0.0 && opt2bSet("-iiact", NFILE, fnm))
3438     {
3439         gmx_fatal(FARGS, "Either provide autocorrelation times (ACTs) with file iact-in.dat "
3440                   "(option -iiact) or define all ACTs with -bs-tau for bootstrapping\n. Not Both.");
3441     }
3442
3443     /* Reading gmx4 pull output and tpr files */
3444     if (opt.bTpr || opt.bPullf || opt.bPullx)
3445     {
3446         read_wham_in(opt.fnTpr, &fninTpr, &nfiles, &opt);
3447
3448         fnPull = opt.bPullf ? opt.fnPullf : opt.fnPullx;
3449         read_wham_in(fnPull, &fninPull, &nfiles2, &opt);
3450         printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n",
3451                nfiles, nfiles2, opt.bPullf ? "force" : "position", opt.fnTpr, fnPull);
3452         if (nfiles != nfiles2)
3453         {
3454             gmx_fatal(FARGS, "Found %d file names in %s, but %d in %s\n", nfiles,
3455                       opt.fnTpr, nfiles2, fnPull);
3456         }
3457
3458         /* Read file that selects the pull group to be used */
3459         if (opt.fnGroupsel != NULL)
3460         {
3461             readPullGroupSelection(&opt, fninTpr, nfiles);
3462         }
3463
3464         window = initUmbrellaWindows(nfiles);
3465         read_tpr_pullxf_files(fninTpr, fninPull, nfiles, &header, window, &opt);
3466     }
3467     else
3468     {   /* reading pdo files */
3469         if  (opt.fnGroupsel != NULL)
3470         {
3471             gmx_fatal(FARGS, "Reading a -is file is not supported with PDO input files.\n"
3472                       "Use awk or a similar tool to pick the required pull groups from your PDO files\n");
3473         }
3474         read_wham_in(opt.fnPdo, &fninPdo, &nfiles, &opt);
3475         printf("Found %d pdo files in %s\n", nfiles, opt.fnPdo);
3476         window = initUmbrellaWindows(nfiles);
3477         read_pdo_files(fninPdo, nfiles, &header, window, &opt);
3478     }
3479     nwins = nfiles;
3480
3481     /* enforce equal weight for all histograms? */
3482     if (opt.bHistEq)
3483     {
3484         enforceEqualWeights(window, nwins);
3485     }
3486
3487     /* write histograms */
3488     histout = xvgropen(opt2fn("-hist", NFILE, fnm), "Umbrella histograms",
3489                        "z", "count", opt.oenv);
3490     for (l = 0; l < opt.bins; ++l)
3491     {
3492         fprintf(histout, "%e\t", (double)(l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
3493         for (i = 0; i < nwins; ++i)
3494         {
3495             for (j = 0; j < window[i].nPull; ++j)
3496             {
3497                 fprintf(histout, "%e\t", window[i].Histo[j][l]);
3498             }
3499         }
3500         fprintf(histout, "\n");
3501     }
3502     ffclose(histout);
3503     printf("Wrote %s\n", opt2fn("-hist", NFILE, fnm));
3504     if (opt.bHistOnly)
3505     {
3506         printf("Wrote histograms to %s, now exiting.\n", opt2fn("-hist", NFILE, fnm));
3507         return 0;
3508     }
3509
3510     /* Using tabulated umbrella potential */
3511     if (opt.bTab)
3512     {
3513         setup_tab(opt2fn("-tab", NFILE, fnm), &opt);
3514     }
3515
3516     /* Integrated autocorrelation times provided ? */
3517     if (opt.bTauIntGiven)
3518     {
3519         readIntegratedAutocorrelationTimes(window, nwins, opt2fn("-iiact", NFILE, fnm));
3520     }
3521
3522     /* Compute integrated autocorrelation times */
3523     if (opt.bCalcTauInt)
3524     {
3525         calcIntegratedAutocorrelationTimes(window, nwins, &opt, opt2fn("-oiact", NFILE, fnm));
3526     }
3527
3528     /* calc average and sigma for each histogram
3529        (maybe required for bootstrapping. If not, this is fast anyhow) */
3530     if (opt.nBootStrap && opt.bsMethod == bsMethod_trajGauss)
3531     {
3532         averageSigma(window, nwins);
3533     }
3534
3535     /* Get initial potential by simple integration */
3536     if (opt.bInitPotByIntegration)
3537     {
3538         guessPotByIntegration(window, nwins, &opt);
3539     }
3540
3541     /* Check if complete reaction coordinate is covered */
3542     checkReactionCoordinateCovered(window, nwins, &opt);
3543
3544     /* Calculate profile */
3545     snew(profile, opt.bins);
3546     if (opt.verbose)
3547     {
3548         opt.stepchange = 1;
3549     }
3550     i = 0;
3551     do
3552     {
3553         if ( (i%opt.stepUpdateContrib) == 0)
3554         {
3555             setup_acc_wham(profile, window, nwins, &opt);
3556         }
3557         if (maxchange < opt.Tolerance)
3558         {
3559             bExact = TRUE;
3560             /* if (opt.verbose) */
3561             printf("Switched to exact iteration in iteration %d\n", i);
3562         }
3563         calc_profile(profile, window, nwins, &opt, bExact);
3564         if (((i%opt.stepchange) == 0 || i == 1) && i != 0)
3565         {
3566             printf("\t%4d) Maximum change %e\n", i, maxchange);
3567         }
3568         i++;
3569     }
3570     while ( (maxchange = calc_z(profile, window, nwins, &opt, bExact)) > opt.Tolerance || !bExact);
3571     printf("Converged in %d iterations. Final maximum change %g\n", i, maxchange);
3572
3573     /* calc error from Kumar's formula */
3574     /* Unclear how the error propagates along reaction coordinate, therefore
3575        commented out  */
3576     /* calc_error_kumar(profile,window, nwins,&opt); */
3577
3578     /* Write profile in energy units? */
3579     if (opt.bLog)
3580     {
3581         prof_normalization_and_unit(profile, &opt);
3582         strcpy(ylabel, en_unit_label[opt.unit]);
3583         strcpy(title, "Umbrella potential");
3584     }
3585     else
3586     {
3587         strcpy(ylabel, "Density of states");
3588         strcpy(title, "Density of states");
3589     }
3590
3591     /* symmetrize profile around z=0? */
3592     if (opt.bSym)
3593     {
3594         symmetrizeProfile(profile, &opt);
3595     }
3596
3597     /* write profile or density of states */
3598     profout = xvgropen(opt2fn("-o", NFILE, fnm), title, "z", ylabel, opt.oenv);
3599     for (i = 0; i < opt.bins; ++i)
3600     {
3601         fprintf(profout, "%e\t%e\n", (double)(i+0.5)/opt.bins*(opt.max-opt.min)+opt.min, profile[i]);
3602     }
3603     ffclose(profout);
3604     printf("Wrote %s\n", opt2fn("-o", NFILE, fnm));
3605
3606     /* Bootstrap Method */
3607     if (opt.nBootStrap)
3608     {
3609         do_bootstrapping(opt2fn("-bsres", NFILE, fnm), opt2fn("-bsprof", NFILE, fnm),
3610                          opt2fn("-hist", NFILE, fnm),
3611                          ylabel, profile, window, nwins, &opt);
3612     }
3613
3614     sfree(profile);
3615     freeUmbrellaWindows(window, nfiles);
3616
3617     printf("\nIn case you use results from g_wham for a publication, please cite:\n");
3618     please_cite(stdout, "Hub2010");
3619
3620     return 0;
3621 }