src/gromacs/gmxana/gmx_wham.cpp

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