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,2018, 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 // NOLINT(clang-analyzer-optin.performance.Padding)
 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 %d 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") != 0)
 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") != 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 = true;
 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<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(const 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(const 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(const 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, "%s", 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 //! Make random weights for histograms for the Bayesian bootstrap of complete histograms)
1552 static void setRandomBsWeights(t_UmbrellaWindow *synthwin, int nAllPull, t_UmbrellaOptions *opt)
1553 {
1554     int     i;
1555     double *r;
1556     gmx::UniformRealDistribution<real> dist(0, nAllPull);
1557
1558     snew(r, nAllPull);
1559
1560     /* generate ordered random numbers between 0 and nAllPull  */
1561     for (i = 0; i < nAllPull-1; i++)
1562     {
1563         r[i] = dist(opt->rng);
1564     }
1565     std::sort(r, r+nAllPull-1);
1566     r[nAllPull-1] = 1.0*nAllPull;
1567
1568     synthwin[0].bsWeight[0] = r[0];
1569     for (i = 1; i < nAllPull; i++)
1570     {
1571         synthwin[i].bsWeight[0] = r[i]-r[i-1];
1572     }
1573
1574     /* avoid to have zero weight by adding a tiny value */
1575     for (i = 0; i < nAllPull; i++)
1576     {
1577         if (synthwin[i].bsWeight[0] < 1e-5)
1578         {
1579             synthwin[i].bsWeight[0] = 1e-5;
1580         }
1581     }
1582
1583     sfree(r);
1584 }
1585
1586 //! The main bootstrapping routine
1587 static void do_bootstrapping(const char *fnres, const char* fnprof, const char *fnhist,
1588                              const char *xlabel, char* ylabel, double *profile,
1589                              t_UmbrellaWindow * window, int nWindows, t_UmbrellaOptions *opt)
1590 {
1591     t_UmbrellaWindow * synthWindow;
1592     double            *bsProfile, *bsProfiles_av, *bsProfiles_av2, maxchange = 1e20, tmp, stddev;
1593     int                i, j, *randomArray = nullptr, winid, pullid, ib;
1594     int                iAllPull, nAllPull, *allPull_winId, *allPull_pullId;
1595     FILE              *fp;
1596     gmx_bool           bExact = FALSE;
1597
1598     /* init random generator */
1599     if (opt->bsSeed == 0)
1600     {
1601         opt->bsSeed = static_cast<int>(gmx::makeRandomSeed());
1602     }
1603     opt->rng.seed(opt->bsSeed);
1604
1605     snew(bsProfile,     opt->bins);
1606     snew(bsProfiles_av, opt->bins);
1607     snew(bsProfiles_av2, opt->bins);
1608
1609     /* Create array of all pull groups. Note that different windows
1610        may have different nr of pull groups
1611        First: Get total nr of pull groups */
1612     nAllPull = 0;
1613     for (i = 0; i < nWindows; i++)
1614     {
1615         nAllPull += window[i].nPull;
1616     }
1617     snew(allPull_winId, nAllPull);
1618     snew(allPull_pullId, nAllPull);
1619     iAllPull = 0;
1620     /* Setup one array of all pull groups */
1621     for (i = 0; i < nWindows; i++)
1622     {
1623         for (j = 0; j < window[i].nPull; j++)
1624         {
1625             allPull_winId[iAllPull]  = i;
1626             allPull_pullId[iAllPull] = j;
1627             iAllPull++;
1628         }
1629     }
1630
1631     /* setup stuff for synthetic windows */
1632     snew(synthWindow, nAllPull);
1633     for (i = 0; i < nAllPull; i++)
1634     {
1635         synthWindow[i].nPull = 1;
1636         synthWindow[i].nBin  = opt->bins;
1637         snew(synthWindow[i].Histo, 1);
1638         if (opt->bsMethod == bsMethod_traj || opt->bsMethod == bsMethod_trajGauss)
1639         {
1640             snew(synthWindow[i].Histo[0], opt->bins);
1641         }
1642         snew(synthWindow[i].N, 1);
1643         snew(synthWindow[i].pos, 1);
1644         snew(synthWindow[i].z, 1);
1645         snew(synthWindow[i].k, 1);
1646         snew(synthWindow[i].bContrib, 1);
1647         snew(synthWindow[i].g, 1);
1648         snew(synthWindow[i].bsWeight, 1);
1649     }
1650
1651     switch (opt->bsMethod)
1652     {
1653         case bsMethod_hist:
1654             snew(randomArray, nAllPull);
1655             printf("\n\nWhen computing statistical errors by bootstrapping entire histograms:\n");
1656             please_cite(stdout, "Hub2006");
1657             break;
1658         case bsMethod_BayesianHist:
1659             /* just copy all histogams into synthWindow array */
1660             for (i = 0; i < nAllPull; i++)
1661             {
1662                 winid  = allPull_winId [i];
1663                 pullid = allPull_pullId[i];
1664                 copy_pullgrp_to_synthwindow(synthWindow+i, window+winid, pullid);
1665             }
1666             break;
1667         case bsMethod_traj:
1668         case bsMethod_trajGauss:
1669             calc_cumulatives(window, nWindows, opt, fnhist, xlabel);
1670             break;
1671         default:
1672             gmx_fatal(FARGS, "Unknown bootstrap method. That should not have happened.\n");
1673     }
1674
1675     /* do bootstrapping */
1676     fp = xvgropen(fnprof, "Bootstrap profiles", xlabel, ylabel, opt->oenv);
1677     for (ib = 0; ib < opt->nBootStrap; ib++)
1678     {
1679         printf("  *******************************************\n"
1680                "  ******** Start bootstrap nr %d ************\n"
1681                "  *******************************************\n", ib+1);
1682
1683         switch (opt->bsMethod)
1684         {
1685             case bsMethod_hist:
1686                 /* bootstrap complete histograms from given histograms */
1687                 getRandomIntArray(nAllPull, opt->histBootStrapBlockLength, randomArray, &opt->rng);
1688                 for (i = 0; i < nAllPull; i++)
1689                 {
1690                     winid  = allPull_winId [randomArray[i]];
1691                     pullid = allPull_pullId[randomArray[i]];
1692                     copy_pullgrp_to_synthwindow(synthWindow+i, window+winid, pullid);
1693                 }
1694                 break;
1695             case bsMethod_BayesianHist:
1696                 /* keep histos, but assign random weights ("Bayesian bootstrap") */
1697                 setRandomBsWeights(synthWindow, nAllPull, opt);
1698                 break;
1699             case bsMethod_traj:
1700             case bsMethod_trajGauss:
1701                 /* create new histos from given histos, that is generate new hypothetical
1702                    trajectories */
1703                 for (i = 0; i < nAllPull; i++)
1704                 {
1705                     winid  = allPull_winId[i];
1706                     pullid = allPull_pullId[i];
1707                     create_synthetic_histo(synthWindow+i, window+winid, pullid, opt);
1708                 }
1709                 break;
1710         }
1711
1712         /* write histos in case of verbose output */
1713         if (opt->bs_verbose)
1714         {
1715             print_histograms(fnhist, synthWindow, nAllPull, ib, opt, xlabel);
1716         }
1717
1718         /* do wham */
1719         i         = 0;
1720         bExact    = FALSE;
1721         maxchange = 1e20;
1722         std::memcpy(bsProfile, profile, opt->bins*sizeof(double)); /* use profile as guess */
1723         do
1724         {
1725             if ( (i%opt->stepUpdateContrib) == 0)
1726             {
1727                 setup_acc_wham(bsProfile, synthWindow, nAllPull, opt);
1728             }
1729             if (maxchange < opt->Tolerance)
1730             {
1731                 bExact = TRUE;
1732             }
1733             if (((i%opt->stepchange) == 0 || i == 1) && i != 0)
1734             {
1735                 printf("\t%4d) Maximum change %e\n", i, maxchange);
1736             }
1737             calc_profile(bsProfile, synthWindow, nAllPull, opt, bExact);
1738             i++;
1739         }
1740         while ( (maxchange = calc_z(bsProfile, synthWindow, nAllPull, opt, bExact)) > opt->Tolerance || !bExact);
1741         printf("\tConverged in %d iterations. Final maximum change %g\n", i, maxchange);
1742
1743         if (opt->bLog)
1744         {
1745             prof_normalization_and_unit(bsProfile, opt);
1746         }
1747
1748         /* symmetrize profile around z=0 */
1749         if (opt->bSym)
1750         {
1751             symmetrizeProfile(bsProfile, opt);
1752         }
1753
1754         /* save stuff to get average and stddev */
1755         for (i = 0; i < opt->bins; i++)
1756         {
1757             tmp                = bsProfile[i];
1758             bsProfiles_av[i]  += tmp;
1759             bsProfiles_av2[i] += tmp*tmp;
1760             fprintf(fp, "%e\t%e\n", (i+0.5)*opt->dz+opt->min, tmp);
1761         }
1762         fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
1763     }
1764     xvgrclose(fp);
1765
1766     /* write average and stddev */
1767     fp = xvgropen(fnres, "Average and stddev from bootstrapping", xlabel, ylabel, opt->oenv);
1768     if (output_env_get_print_xvgr_codes(opt->oenv))
1769     {
1770         fprintf(fp, "@TYPE xydy\n");
1771     }
1772     for (i = 0; i < opt->bins; i++)
1773     {
1774         bsProfiles_av [i] /= opt->nBootStrap;
1775         bsProfiles_av2[i] /= opt->nBootStrap;
1776         tmp                = bsProfiles_av2[i]-gmx::square(bsProfiles_av[i]);
1777         stddev             = (tmp >= 0.) ? std::sqrt(tmp) : 0.; /* Catch rouding errors */
1778         fprintf(fp, "%e\t%e\t%e\n", (i+0.5)*opt->dz+opt->min, bsProfiles_av [i], stddev);
1779     }
1780     xvgrclose(fp);
1781     printf("Wrote boot strap result to %s\n", fnres);
1782 }
1783
1784 //! Return type of input file based on file extension (xvg, pdo, or tpr)
1785 static int whaminFileType(char *fn)
1786 {
1787     int len;
1788     len = std::strlen(fn);
1789     if (std::strcmp(fn+len-3, "tpr") == 0)
1790     {
1791         return whamin_tpr;
1792     }
1793     else if (std::strcmp(fn+len-3, "xvg") == 0 || std::strcmp(fn+len-6, "xvg.gz") == 0)
1794     {
1795         return whamin_pullxf;
1796     }
1797     else if (std::strcmp(fn+len-3, "pdo") == 0 || std::strcmp(fn+len-6, "pdo.gz") == 0)
1798     {
1799         return whamin_pdo;
1800     }
1801     else
1802     {
1803         gmx_fatal(FARGS, "Unknown file type of %s. Should be tpr, xvg, or pdo.\n", fn);
1804     }
1805 }
1806
1807 //! Read the files names in pdo-files.dat, pullf/x-files.dat, tpr-files.dat
1808 static void read_wham_in(const char *fn, char ***filenamesRet, int *nfilesRet,
1809                          t_UmbrellaOptions *opt)
1810 {
1811     char **filename = nullptr, tmp[WHAM_MAXFILELEN+2];
1812     int    nread, sizenow, i, block = 1;
1813     FILE  *fp;
1814
1815     fp      = gmx_ffopen(fn, "r");
1816     nread   = 0;
1817     sizenow = 0;
1818     while (fgets(tmp, sizeof(tmp), fp) != nullptr)
1819     {
1820         if (std::strlen(tmp) >= WHAM_MAXFILELEN)
1821         {
1822             gmx_fatal(FARGS, "Filename too long in %s. Only %d characters allowed.\n", fn, WHAM_MAXFILELEN);
1823         }
1824         if (nread >= sizenow)
1825         {
1826             sizenow += block;
1827             srenew(filename, sizenow);
1828             for (i = sizenow-block; i < sizenow; i++)
1829             {
1830                 snew(filename[i], WHAM_MAXFILELEN);
1831             }
1832         }
1833         /* remove newline if there is one */
1834         if (tmp[std::strlen(tmp)-1] == '\n')
1835         {
1836             tmp[std::strlen(tmp)-1] = '\0';
1837         }
1838         std::strcpy(filename[nread], tmp);
1839         if (opt->verbose)
1840         {
1841             printf("Found file %s in %s\n", filename[nread], fn);
1842         }
1843         nread++;
1844     }
1845     *filenamesRet = filename;
1846     *nfilesRet    = nread;
1847 }
1848
1849 //! Open a file or a pipe to a gzipped file
1850 static FILE *open_pdo_pipe(const char *fn, t_UmbrellaOptions *opt, gmx_bool *bPipeOpen)
1851 {
1852     char            Buffer[1024], gunzip[1024], *Path = nullptr;
1853     FILE           *pipe   = nullptr;
1854     static gmx_bool bFirst = true;
1855
1856     /* gzipped pdo file? */
1857     if ((std::strcmp(fn+std::strlen(fn)-3, ".gz") == 0))
1858     {
1859         /* search gunzip executable */
1860         if (!(Path = getenv("GMX_PATH_GZIP")))
1861         {
1862             if (gmx_fexist("/bin/gunzip"))
1863             {
1864                 sprintf(gunzip, "%s", "/bin/gunzip");
1865             }
1866             else if (gmx_fexist("/usr/bin/gunzip"))
1867             {
1868                 sprintf(gunzip, "%s", "/usr/bin/gunzip");
1869             }
1870             else
1871             {
1872                 gmx_fatal(FARGS, "Cannot find executable gunzip in /bin or /usr/bin.\n"
1873                           "You may want to define the path to gunzip "
1874                           "with the environment variable GMX_PATH_GZIP.", gunzip);
1875             }
1876         }
1877         else
1878         {
1879             sprintf(gunzip, "%s/gunzip", Path);
1880             if (!gmx_fexist(gunzip))
1881             {
1882                 gmx_fatal(FARGS, "Cannot find executable %s. Please define the path to gunzip"
1883                           " in the environmental varialbe GMX_PATH_GZIP.", gunzip);
1884             }
1885         }
1886         if (bFirst)
1887         {
1888             printf("Using gunzip executable %s\n", gunzip);
1889             bFirst = false;
1890         }
1891         if (!gmx_fexist(fn))
1892         {
1893             gmx_fatal(FARGS, "File %s does not exist.\n", fn);
1894         }
1895         sprintf(Buffer, "%s -c < %s", gunzip, fn);
1896         if (opt->verbose)
1897         {
1898             printf("Executing command '%s'\n", Buffer);
1899         }
1900 #if HAVE_PIPES
1901         if ((pipe = popen(Buffer, "r")) == nullptr)
1902         {
1903             gmx_fatal(FARGS, "Unable to open pipe to `%s'\n", Buffer);
1904         }
1905 #else
1906         gmx_fatal(FARGS, "Cannot open a compressed file on platform without pipe support");
1907 #endif
1908         *bPipeOpen = TRUE;
1909     }
1910     else
1911     {
1912         pipe       = gmx_ffopen(fn, "r");
1913         *bPipeOpen = FALSE;
1914     }
1915
1916     return pipe;
1917 }
1918
1919 //! Close file or pipe
1920 static void pdo_close_file(FILE *fp)
1921 {
1922 #if HAVE_PIPES
1923     pclose(fp);
1924 #else
1925     gmx_ffclose(fp);
1926 #endif
1927 }
1928
1929 //! Reading all pdo files
1930 static void read_pdo_files(char **fn, int nfiles, t_UmbrellaHeader* header,
1931                            t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
1932 {
1933     FILE    *file;
1934     real     mintmp, maxtmp, done = 0.;
1935     int      i;
1936     gmx_bool bPipeOpen;
1937     /* char Buffer0[1000]; */
1938
1939     if (nfiles < 1)
1940     {
1941         gmx_fatal(FARGS, "No files found. Hick.");
1942     }
1943
1944     /* if min and max are not given, get min and max from the input files */
1945     if (opt->bAuto)
1946     {
1947         printf("Automatic determination of boundaries from %d pdo files...\n", nfiles);
1948         opt->min = 1e20;
1949         opt->max = -1e20;
1950         for (i = 0; i < nfiles; ++i)
1951         {
1952             file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
1953             /*fgets(Buffer0,999,file);
1954                fprintf(stderr,"First line '%s'\n",Buffer0); */
1955             done = 100.0*(i+1)/nfiles;
1956             fprintf(stdout, "\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
1957             if (opt->verbose)
1958             {
1959                 printf("\n");
1960             }
1961             read_pdo_header(file, header, opt);
1962             /* here only determine min and max of this window */
1963             read_pdo_data(file, header, i, nullptr, opt, TRUE, &mintmp, &maxtmp);
1964             if (maxtmp > opt->max)
1965             {
1966                 opt->max = maxtmp;
1967             }
1968             if (mintmp < opt->min)
1969             {
1970                 opt->min = mintmp;
1971             }
1972             if (bPipeOpen)
1973             {
1974                 pdo_close_file(file);
1975             }
1976             else
1977             {
1978                 gmx_ffclose(file);
1979             }
1980         }
1981         printf("\n");
1982         printf("\nDetermined boundaries to %f and %f\n\n", opt->min, opt->max);
1983         if (opt->bBoundsOnly)
1984         {
1985             printf("Found option -boundsonly, now exiting.\n");
1986             exit (0);
1987         }
1988     }
1989     /* store stepsize in profile */
1990     opt->dz = (opt->max-opt->min)/opt->bins;
1991
1992     /* Having min and max, we read in all files */
1993     /* Loop over all files */
1994     for (i = 0; i < nfiles; ++i)
1995     {
1996         done = 100.0*(i+1)/nfiles;
1997         fprintf(stdout, "\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
1998         if (opt->verbose)
1999         {
2000             printf("\n");
2001         }
2002         file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
2003         read_pdo_header(file, header, opt);
2004         /* load data into window */
2005         read_pdo_data(file, header, i, window, opt, FALSE, nullptr, nullptr);
2006         if ((window+i)->Ntot[0] == 0)
2007         {
2008             fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fn[i]);
2009         }
2010         if (bPipeOpen)
2011         {
2012             pdo_close_file(file);
2013         }
2014         else
2015         {
2016             gmx_ffclose(file);
2017         }
2018     }
2019     printf("\n");
2020     for (i = 0; i < nfiles; ++i)
2021     {
2022         sfree(fn[i]);
2023     }
2024     sfree(fn);
2025 }
2026
2027 //! translate 0/1 to N/Y to write pull dimensions
2028 #define int2YN(a) (((a) == 0) ? ("N") : ("Y"))
2029
2030 //! Read pull groups from a tpr file (including position, force const, geometry, number of groups)
2031 static void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions *opt, t_coordselection *coordsel)
2032 {
2033     t_inputrec      irInstance;
2034     t_inputrec     *ir = &irInstance;
2035     t_state         state;
2036     static int      first = 1;
2037
2038     /* printf("Reading %s \n",fn); */
2039     read_tpx_state(fn, ir, &state, nullptr);
2040
2041     if (!ir->bPull)
2042     {
2043         gmx_fatal(FARGS, "This is not a tpr with COM pulling");
2044     }
2045     if (ir->pull->ncoord == 0)
2046     {
2047         gmx_fatal(FARGS, "No pull coordinates found in %s", fn);
2048     }
2049
2050     /* Read overall pull info */
2051     header->npullcrds      = ir->pull->ncoord;
2052     header->bPrintCOM      = ir->pull->bPrintCOM;
2053     header->bPrintRefValue = ir->pull->bPrintRefValue;
2054     header->bPrintComp     = ir->pull->bPrintComp;
2055
2056     /* Read pull coordinates */
2057     snew(header->pcrd, header->npullcrds);
2058     for (int i = 0; i < ir->pull->ncoord; i++)
2059     {
2060         header->pcrd[i].pull_type     = ir->pull->coord[i].eType;
2061         header->pcrd[i].geometry      = ir->pull->coord[i].eGeom;
2062         header->pcrd[i].ngroup        = ir->pull->coord[i].ngroup;
2063         header->pcrd[i].k             = ir->pull->coord[i].k;
2064         header->pcrd[i].init_dist     = ir->pull->coord[i].init;
2065
2066         copy_ivec(ir->pull->coord[i].dim, header->pcrd[i].dim);
2067         header->pcrd[i].ndim         = header->pcrd[i].dim[XX] + header->pcrd[i].dim[YY] + header->pcrd[i].dim[ZZ];
2068
2069         std::strcpy(header->pcrd[i].coord_unit,
2070                     pull_coordinate_units(&ir->pull->coord[i]));
2071
2072         if (ir->efep != efepNO && ir->pull->coord[i].k != ir->pull->coord[i].kB)
2073         {
2074             gmx_fatal(FARGS, "Seems like you did free-energy perturbation, and you perturbed the force constant."
2075                       " This is not supported.\n");
2076         }
2077         if (coordsel && (coordsel->n != ir->pull->ncoord))
2078         {
2079             gmx_fatal(FARGS, "Found %d pull coordinates in %s, but %d columns in the respective line\n"
2080                       "coordinate selection file (option -is)\n", ir->pull->ncoord, fn, coordsel->n);
2081         }
2082     }
2083
2084     /* Check pull coords for consistency */
2085     int  geom          = -1;
2086     ivec thedim        = { 0, 0, 0 };
2087     bool geometryIsSet = false;
2088     for (int i = 0; i < ir->pull->ncoord; i++)
2089     {
2090         if (coordsel == nullptr || coordsel->bUse[i])
2091         {
2092             if (header->pcrd[i].pull_type != epullUMBRELLA)
2093             {
2094                 gmx_fatal(FARGS, "%s: Pull coordinate %d is of type \"%s\", expected \"umbrella\". Only umbrella coodinates can enter WHAM.\n"
2095                           "If you have umrella and non-umbrella coordinates, you can select the umbrella coordinates with gmx wham -is\n",
2096                           fn, i+1, epull_names[header->pcrd[i].pull_type]);
2097             }
2098             if (!geometryIsSet)
2099             {
2100                 geom = header->pcrd[i].geometry;
2101                 copy_ivec(header->pcrd[i].dim, thedim);
2102                 geometryIsSet = true;
2103             }
2104             if (geom != header->pcrd[i].geometry)
2105             {
2106                 gmx_fatal(FARGS, "%s: Your pull coordinates have different pull geometry (coordinate 1: %s, coordinate %d: %s)\n"
2107                           "If you want to use only some pull coordinates in WHAM, please select them with option gmx wham -is\n",
2108                           fn, epullg_names[geom], i+1, epullg_names[header->pcrd[i].geometry]);
2109             }
2110             if (thedim[XX] != header->pcrd[i].dim[XX] || thedim[YY] != header->pcrd[i].dim[YY] || thedim[ZZ] != header->pcrd[i].dim[ZZ])
2111             {
2112                 gmx_fatal(FARGS, "%s: Your pull coordinates have different pull dimensions (coordinate 1: %s %s %s, coordinate %d: %s %s %s)\n"
2113                           "If you want to use only some pull coordinates in WHAM, please select them with option gmx wham -is\n",
2114                           fn, int2YN(thedim[XX]), int2YN(thedim[YY]), int2YN(thedim[ZZ]), i+1,
2115                           int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]), int2YN(header->pcrd[i].dim[ZZ]));
2116             }
2117             if (header->pcrd[i].geometry == epullgCYL)
2118             {
2119                 if (header->pcrd[i].dim[XX] || header->pcrd[i].dim[YY] || (!header->pcrd[i].dim[ZZ]))
2120                 {
2121                     gmx_fatal(FARGS, "With pull geometry 'cylinder', expected pulling in Z direction only.\n"
2122                               "However, found dimensions [%s %s %s]\n",
2123                               int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]),
2124                               int2YN(header->pcrd[i].dim[ZZ]));
2125                 }
2126             }
2127             if (header->pcrd[i].k <= 0.0)
2128             {
2129                 gmx_fatal(FARGS, "%s: Pull coordinate %d has force constant of of %g in %s.\n"
2130                           "That doesn't seem to be an Umbrella tpr.\n",
2131                           fn, i+1, header->pcrd[i].k);
2132             }
2133         }
2134     }
2135
2136     if (opt->verbose || first)
2137     {
2138         printf("\nFile %s, %d coordinates, with these options:\n", fn, header->npullcrds);
2139         int maxlen = 0;
2140         for (int i = 0; i < ir->pull->ncoord; i++)
2141         {
2142             int lentmp = strlen(epullg_names[header->pcrd[i].geometry]);
2143             maxlen     = (lentmp > maxlen) ? lentmp : maxlen;
2144         }
2145         char fmt[STRLEN];
2146         sprintf(fmt, "\tGeometry %%-%ds  k = %%-8g  position = %%-8g  dimensions [%%s %%s %%s] (%%d dimensions). Used: %%s\n",
2147                 maxlen+1);
2148         for (int i = 0; i < ir->pull->ncoord; i++)
2149         {
2150             bool use = (coordsel == nullptr || coordsel->bUse[i]);
2151             printf(fmt,
2152                    epullg_names[header->pcrd[i].geometry], header->pcrd[i].k, header->pcrd[i].init_dist,
2153                    int2YN(header->pcrd[i].dim[XX]), int2YN(header->pcrd[i].dim[YY]), int2YN(header->pcrd[i].dim[ZZ]),
2154                    header->pcrd[i].ndim, use ? "Yes" : "No");
2155             printf("\tPull group coordinates of %d groups expected in pullx files.\n", ir->pull->bPrintCOM ? header->pcrd[i].ngroup : 0);
2156         }
2157         printf("\tReference value of the coordinate%s expected in pullx files.\n",
2158                header->bPrintRefValue ? "" : " not");
2159     }
2160     if (!opt->verbose && first)
2161     {
2162         printf("\tUse option -v to see this output for all input tpr files\n\n");
2163     }
2164
2165     first = 0;
2166 }
2167
2168 //! Read pullx.xvg or pullf.xvg
2169 static void read_pull_xf(const char *fn, t_UmbrellaHeader * header,
2170                          t_UmbrellaWindow * window,
2171                          t_UmbrellaOptions *opt,
2172                          gmx_bool bGetMinMax, real *mintmp, real *maxtmp,
2173                          t_coordselection *coordsel)
2174 {
2175     double        **y = nullptr, pos = 0., t, force, time0 = 0., dt;
2176     int             ny, nt, bins, ibin, i, g, gUsed, dstep = 1;
2177     int             nColExpect, ntot, column;
2178     real            min, max, minfound = 1e20, maxfound = -1e20;
2179     gmx_bool        dt_ok, timeok;
2180     const char     *quantity;
2181     const int       blocklen = 4096;
2182     int            *lennow   = nullptr;
2183     static gmx_bool bFirst   = TRUE;
2184
2185     /*
2186      * Data columns in pull output:
2187      *  - in force output pullf.xvg:
2188      *    No reference columns, one column per pull coordinate
2189      *
2190      *  - in position output pullx.xvg:
2191      *     * optionally, ndim columns for COMs of all groups (depending on on mdp options pull-print-com);
2192      *     * The displacement, always one column. Note: with pull-print-components = yes, the dx/dy/dz would
2193      *       be written separately into pullx file, but this is not supported and throws an error below;
2194      *     * optionally, the position of the reference coordinate (depending on pull-print-ref-value)
2195      */
2196
2197     if (header->bPrintComp && opt->bPullx)
2198     {
2199         gmx_fatal(FARGS, "gmx wham cannot read pullx files if the components of the coordinate was written\n"
2200                   "(mdp option pull-print-components). Provide the pull force files instead (with option -if).\n");
2201     }
2202
2203     int *nColThisCrd, *nColCOMCrd, *nColRefCrd;
2204     snew(nColThisCrd, header->npullcrds);
2205     snew(nColCOMCrd,  header->npullcrds);
2206     snew(nColRefCrd,  header->npullcrds);
2207
2208     if (!opt->bPullx)
2209     {
2210         /* pullf reading: simply one column per coordinate */
2211         for (g = 0; g < header->npullcrds; g++)
2212         {
2213             nColThisCrd[g] = 1;
2214             nColCOMCrd[g]  = 0;
2215             nColRefCrd[g]  = 0;
2216         }
2217     }
2218     else
2219     {
2220         /* pullx reading. Note explanation above. */
2221         for (g = 0; g < header->npullcrds; g++)
2222         {
2223             nColRefCrd[g]   = (header->bPrintRefValue ? 1 : 0);
2224             nColCOMCrd[g]   = (header->bPrintCOM ? header->pcrd[g].ndim*header->pcrd[g].ngroup : 0);
2225             nColThisCrd[g]  = 1 + nColCOMCrd[g] + nColRefCrd[g];
2226         }
2227     }
2228
2229     nColExpect = 1; /* time column */
2230     for (g = 0; g < header->npullcrds; g++)
2231     {
2232         nColExpect += nColThisCrd[g];
2233     }
2234
2235     /* read pullf or pullx. Could be optimized if min and max are given. */
2236     nt = read_xvg(fn, &y, &ny);
2237
2238     /* Check consistency */
2239     quantity  = opt->bPullx ? "position" : "force";
2240     if (nt < 1)
2241     {
2242         gmx_fatal(FARGS, "Empty pull %s file %s\n", quantity, fn);
2243     }
2244     if (bFirst || opt->verbose)
2245     {
2246         printf("\nReading pull %s file %s, expecting %d columns:\n", quantity, fn, nColExpect);
2247         for (i = 0; i < header->npullcrds; i++)
2248         {
2249             printf("\tColumns for pull coordinate %d\n", i+1);
2250             printf("\t\treaction coordinate:             %d\n"
2251                    "\t\tcenter-of-mass of groups:        %d\n"
2252                    "\t\treference position column:       %s\n",
2253                    1, nColCOMCrd[i], (header->bPrintRefValue ? "Yes" : "No"));
2254         }
2255         printf("\tFound %d times in %s\n", nt, fn);
2256         bFirst = FALSE;
2257     }
2258     if (nColExpect != ny)
2259     {
2260         gmx_fatal(FARGS, "Expected %d columns (including time column) in %s, but found %d."
2261                   " Maybe you confused options -if and -ix ?", nColExpect, fn, ny);
2262     }
2263
2264     if (!bGetMinMax)
2265     {
2266         assert(window);
2267         bins = opt->bins;
2268         min  = opt->min;
2269         max  = opt->max;
2270         if (nt > 1)
2271         {
2272             window->dt = y[0][1]-y[0][0];
2273         }
2274         else if (opt->nBootStrap && opt->tauBootStrap != 0.0)
2275         {
2276             fprintf(stderr, "\n *** WARNING, Could not determine time step in %s\n", fn);
2277         }
2278
2279         /* Need to alocate memory and set up structure for windows */
2280         if (coordsel)
2281         {
2282             /* Use only groups selected with option -is file */
2283             if (header->npullcrds != coordsel->n)
2284             {
2285                 gmx_fatal(FARGS, "tpr file contains %d pull groups, but expected %d from group selection file\n",
2286                           header->npullcrds, coordsel->n);
2287             }
2288             window->nPull = coordsel->nUse;
2289         }
2290         else
2291         {
2292             window->nPull = header->npullcrds;
2293         }
2294
2295         window->nBin = bins;
2296         snew(window->Histo,    window->nPull);
2297         snew(window->z,        window->nPull);
2298         snew(window->k,        window->nPull);
2299         snew(window->pos,      window->nPull);
2300         snew(window->N,        window->nPull);
2301         snew(window->Ntot,     window->nPull);
2302         snew(window->g,        window->nPull);
2303         snew(window->bsWeight, window->nPull);
2304         window->bContrib = nullptr;
2305
2306         if (opt->bCalcTauInt)
2307         {
2308             snew(window->ztime, window->nPull);
2309         }
2310         else
2311         {
2312             window->ztime = nullptr;
2313         }
2314         snew(lennow, window->nPull);
2315
2316         for (g = 0; g < window->nPull; ++g)
2317         {
2318             window->z        [g] = 1;
2319             window->bsWeight [g] = 1.;
2320             window->N        [g] = 0;
2321             window->Ntot     [g] = 0;
2322             window->g        [g] = 1.;
2323             snew(window->Histo[g], bins);
2324
2325             if (opt->bCalcTauInt)
2326             {
2327                 window->ztime[g] = nullptr;
2328             }
2329         }
2330
2331         /* Copying umbrella center and force const is more involved since not
2332            all pull groups from header (tpr file) may be used in window variable */
2333         for (g = 0, gUsed = 0; g < header->npullcrds; ++g)
2334         {
2335             if (coordsel && !coordsel->bUse[g])
2336             {
2337                 continue;
2338             }
2339             window->k  [gUsed] = header->pcrd[g].k;
2340             window->pos[gUsed] = header->pcrd[g].init_dist;
2341             gUsed++;
2342         }
2343     }
2344     else
2345     {   /* only determine min and max */
2346         minfound = 1e20;
2347         maxfound = -1e20;
2348         min      = max = bins = 0; /* Get rid of warnings */
2349     }
2350
2351
2352     for (i = 0; i < nt; i++)
2353     {
2354         /* Do you want that time frame? */
2355         t = 1.0/1000*( static_cast<int64_t> ((y[0][i]*1000) + 0.5)); /* round time to fs */
2356
2357         /* get time step of pdo file and get dstep from opt->dt */
2358         if (i == 0)
2359         {
2360             time0 = t;
2361         }
2362         else if (i == 1)
2363         {
2364             dt = t-time0;
2365             if (opt->dt > 0.0)
2366             {
2367                 dstep = static_cast<int>(opt->dt/dt+0.5);
2368                 if (dstep == 0)
2369                 {
2370                     dstep = 1;
2371                 }
2372             }
2373             if (!bGetMinMax)
2374             {
2375                 window->dt = dt*dstep;
2376             }
2377         }
2378
2379         dt_ok  = (i%dstep == 0);
2380         timeok = (dt_ok && t >= opt->tmin && t <= opt->tmax);
2381         /*if (opt->verbose)
2382            printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
2383            t,opt->tmin, opt->tmax, dt_ok,timeok); */
2384         if (timeok)
2385         {
2386             /* Note: if coordsel == NULL:
2387              *          all groups in pullf/x file are stored in this window, and gUsed == g
2388              *       if coordsel != NULL:
2389              *          only groups with coordsel.bUse[g]==TRUE are stored. gUsed is not always equal g
2390              */
2391             gUsed  = -1;
2392             for (g = 0; g < header->npullcrds; ++g)
2393             {
2394                 /* was this group selected for application in WHAM? */
2395                 if (coordsel && !coordsel->bUse[g])
2396                 {
2397                     continue;
2398                 }
2399                 gUsed++;
2400
2401                 if (opt->bPullf)
2402                 {
2403                     /* y has 1 time column y[0] and one column per force y[1],...,y[nCrds] */
2404                     force = y[g+1][i];
2405                     pos   = -force/header->pcrd[g].k + header->pcrd[g].init_dist;
2406                 }
2407                 else
2408                 {
2409                     /* Pick the correct column index.
2410                        Note that there is always exactly one displacement column.
2411                      */
2412                     column = 1;
2413                     for (int j = 0; j < g; j++)
2414                     {
2415                         column += nColThisCrd[j];
2416                     }
2417                     pos     = y[column][i];
2418                 }
2419
2420                 /* printf("crd %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
2421                 if (bGetMinMax)
2422                 {
2423                     if (pos < minfound)
2424                     {
2425                         minfound = pos;
2426                     }
2427                     if (pos > maxfound)
2428                     {
2429                         maxfound = pos;
2430                     }
2431                 }
2432                 else
2433                 {
2434                     if (gUsed >= window->nPull)
2435                     {
2436                         gmx_fatal(FARGS, "gUsed too large (%d, nPull=%d). This error should have been caught before.\n",
2437                                   gUsed, window->nPull);
2438                     }
2439
2440                     if (opt->bCalcTauInt && !bGetMinMax)
2441                     {
2442                         /* save time series for autocorrelation analysis */
2443                         ntot = window->Ntot[gUsed];
2444                         /* printf("i %d, ntot %d, lennow[g] = %d\n",i,ntot,lennow[g]); */
2445                         if (ntot >= lennow[gUsed])
2446                         {
2447                             lennow[gUsed] += blocklen;
2448                             srenew(window->ztime[gUsed], lennow[gUsed]);
2449                         }
2450                         window->ztime[gUsed][ntot] = pos;
2451                     }
2452
2453                     ibin = static_cast<int> (std::floor((pos-min)/(max-min)*bins));
2454                     if (opt->bCycl)
2455                     {
2456                         if (ibin < 0)
2457                         {
2458                             while ( (ibin += bins) < 0)
2459                             {
2460                                 ;
2461                             }
2462                         }
2463                         else if (ibin >= bins)
2464                         {
2465                             while ( (ibin -= bins) >= bins)
2466                             {
2467                                 ;
2468                             }
2469                         }
2470                     }
2471                     if (ibin >= 0 && ibin < bins)
2472                     {
2473                         window->Histo[gUsed][ibin] += 1.;
2474                         window->N[gUsed]++;
2475                     }
2476                     window->Ntot[gUsed]++;
2477                 }
2478             }
2479         }
2480         else if (t > opt->tmax)
2481         {
2482             if (opt->verbose)
2483             {
2484                 printf("time %f larger than tmax %f, stop reading this pullx/pullf file\n", t, opt->tmax);
2485             }
2486             break;
2487         }
2488     }
2489
2490     if (bGetMinMax)
2491     {
2492         *mintmp = minfound;
2493         *maxtmp = maxfound;
2494     }
2495     sfree(lennow);
2496     for (i = 0; i < ny; i++)
2497     {
2498         sfree(y[i]);
2499     }
2500 }
2501
2502 //! read pullf-files.dat or pullx-files.dat and tpr-files.dat
2503 static void read_tpr_pullxf_files(char **fnTprs, char **fnPull, int nfiles,
2504                                   t_UmbrellaHeader* header,
2505                                   t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
2506 {
2507     int  i;
2508     real mintmp, maxtmp;
2509
2510     printf("Reading %d tpr and pullf files\n", nfiles);
2511
2512     /* min and max not given? */
2513     if (opt->bAuto)
2514     {
2515         printf("Automatic determination of boundaries...\n");
2516         opt->min = 1e20;
2517         opt->max = -1e20;
2518         for (i = 0; i < nfiles; i++)
2519         {
2520             if (whaminFileType(fnTprs[i]) != whamin_tpr)
2521             {
2522                 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a tpr file\n", i);
2523             }
2524             read_tpr_header(fnTprs[i], header, opt, (opt->nCoordsel > 0) ? &opt->coordsel[i] : nullptr);
2525             if (whaminFileType(fnPull[i]) != whamin_pullxf)
2526             {
2527                 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
2528             }
2529             read_pull_xf(fnPull[i], header, nullptr, opt, TRUE, &mintmp, &maxtmp,
2530                          (opt->nCoordsel > 0) ? &opt->coordsel[i] : nullptr);
2531             if (maxtmp > opt->max)
2532             {
2533                 opt->max = maxtmp;
2534             }
2535             if (mintmp < opt->min)
2536             {
2537                 opt->min = mintmp;
2538             }
2539         }
2540         printf("\nDetermined boundaries to %f and %f\n\n", opt->min, opt->max);
2541         if (opt->bBoundsOnly)
2542         {
2543             printf("Found option -boundsonly, now exiting.\n");
2544             exit (0);
2545         }
2546     }
2547     /* store stepsize in profile */
2548     opt->dz = (opt->max-opt->min)/opt->bins;
2549
2550     bool foundData = false;
2551     for (i = 0; i < nfiles; i++)
2552     {
2553         if (whaminFileType(fnTprs[i]) != whamin_tpr)
2554         {
2555             gmx_fatal(FARGS, "Expected the %d'th file in input file to be a tpr file\n", i);
2556         }
2557         read_tpr_header(fnTprs[i], header, opt, (opt->nCoordsel > 0) ? &opt->coordsel[i] : nullptr);
2558         if (whaminFileType(fnPull[i]) != whamin_pullxf)
2559         {
2560             gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
2561         }
2562         read_pull_xf(fnPull[i], header, window+i, opt, FALSE, nullptr, nullptr,
2563                      (opt->nCoordsel > 0) ? &opt->coordsel[i] : nullptr);
2564         if (window[i].Ntot[0] == 0)
2565         {
2566             fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fnPull[i]);
2567         }
2568         else
2569         {
2570             foundData = true;
2571         }
2572     }
2573     if (!foundData)
2574     {
2575         gmx_fatal(FARGS, "No data points were found in pullf/pullx files. Maybe you need to specify the -b option?\n");
2576     }
2577
2578     for (i = 0; i < nfiles; i++)
2579     {
2580         sfree(fnTprs[i]);
2581         sfree(fnPull[i]);
2582     }
2583     sfree(fnTprs);
2584     sfree(fnPull);
2585 }
2586
2587 /*! \brief Read integrated autocorrelation times from input file (option -iiact)
2588  *
2589  * Note: Here we consider tau[int] := int_0^inf ACF(t) as the integrated autocorrelation times.
2590  * The factor `g := 1 + 2*tau[int]` subsequently enters the uncertainty.
2591  */
2592 static void readIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins, const char* fn)
2593 {
2594     int      nlines, ny, i, ig;
2595     double **iact;
2596
2597     printf("Readging Integrated autocorrelation times from %s ...\n", fn);
2598     nlines = read_xvg(fn, &iact, &ny);
2599     if (nlines != nwins)
2600     {
2601         gmx_fatal(FARGS, "Found %d lines with integrated autocorrelation times in %s.\nExpected %d",
2602                   nlines, fn, nwins);
2603     }
2604     for (i = 0; i < nlines; i++)
2605     {
2606         if (window[i].nPull != ny)
2607         {
2608             gmx_fatal(FARGS, "You are providing autocorrelation times with option -iiact and the\n"
2609                       "number of pull groups is different in different simulations. That is not\n"
2610                       "supported yet. Sorry.\n");
2611         }
2612         for (ig = 0; ig < window[i].nPull; ig++)
2613         {
2614             /* compare Kumar et al, J Comp Chem 13, 1011-1021 (1992) */
2615             window[i].g[ig] = 1+2*iact[ig][i]/window[i].dt;
2616
2617             if (iact[ig][i] <= 0.0)
2618             {
2619                 fprintf(stderr, "\nWARNING, IACT = %f (window %d, group %d)\n", iact[ig][i], i, ig);
2620             }
2621         }
2622     }
2623 }
2624
2625
2626 /*! \brief Smooth autocorreltion times along the reaction coordinate.
2627  *
2628  * This is useful
2629  * if the ACT is subject to high uncertainty in case if limited sampling. Note
2630  * that -in case of limited sampling- the ACT may be severely underestimated.
2631  * Note: the g=1+2tau are overwritten.
2632  * If opt->bAllowReduceIact==FALSE, the ACTs are never reduced, only increased
2633  * by the smoothing
2634  */
2635 static void smoothIact(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt)
2636 {
2637     int    i, ig, j, jg;
2638     double pos, dpos2, siglim, siglim2, gaufact, invtwosig2, w, weight, tausm;
2639
2640     /* only evaluate within +- 3sigma of the Gausian */
2641     siglim  = 3.0*opt->sigSmoothIact;
2642     siglim2 = gmx::square(siglim);
2643     /* pre-factor of Gaussian */
2644     gaufact    = 1.0/(std::sqrt(2*M_PI)*opt->sigSmoothIact);
2645     invtwosig2 = 0.5/gmx::square(opt->sigSmoothIact);
2646
2647     for (i = 0; i < nwins; i++)
2648     {
2649         snew(window[i].tausmooth, window[i].nPull);
2650         for (ig = 0; ig < window[i].nPull; ig++)
2651         {
2652             tausm  = 0.;
2653             weight = 0;
2654             pos    = window[i].pos[ig];
2655             for (j = 0; j < nwins; j++)
2656             {
2657                 for (jg = 0; jg < window[j].nPull; jg++)
2658                 {
2659                     dpos2 = gmx::square(window[j].pos[jg]-pos);
2660                     if (dpos2 < siglim2)
2661                     {
2662                         w       = gaufact*std::exp(-dpos2*invtwosig2);
2663                         weight += w;
2664                         tausm  += w*window[j].tau[jg];
2665                         /*printf("Weight %g dpos2=%g pos=%g gaufact=%g invtwosig2=%g\n",
2666                            w,dpos2,pos,gaufact,invtwosig2); */
2667                     }
2668                 }
2669             }
2670             tausm /= weight;
2671             if (opt->bAllowReduceIact || tausm > window[i].tau[ig])
2672             {
2673                 window[i].tausmooth[ig] = tausm;
2674             }
2675             else
2676             {
2677                 window[i].tausmooth[ig] = window[i].tau[ig];
2678             }
2679             window[i].g[ig] = 1+2*tausm/window[i].dt;
2680         }
2681     }
2682 }
2683
2684 //! Stop integrating autoccorelation function when ACF drops under this value
2685 #define WHAM_AC_ZERO_LIMIT 0.05
2686
2687 /*! \brief Try to compute the autocorrelation time for each umbrealla window
2688  */
2689 static void calcIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins,
2690                                                t_UmbrellaOptions *opt, const char *fn, const char *xlabel)
2691 {
2692     int   i, ig, ncorr, ntot, j, k, *count, restart;
2693     real *corr, c0, dt, tmp;
2694     real *ztime, av, tausteps;
2695     FILE *fp, *fpcorr = nullptr;
2696
2697     if (opt->verbose)
2698     {
2699         fpcorr = xvgropen("hist_autocorr.xvg", "Autocorrelation functions of umbrella windows",
2700                           "time [ps]", "autocorrelation function", opt->oenv);
2701     }
2702
2703     printf("\n");
2704     for (i = 0; i < nwins; i++)
2705     {
2706         fprintf(stdout, "\rEstimating integrated autocorrelation times ... [%2.0f%%] ...", 100.*(i+1)/nwins);
2707         fflush(stdout);
2708         ntot = window[i].Ntot[0];
2709
2710         /* using half the maximum time as length of autocorrelation function */
2711         ncorr = ntot/2;
2712         if (ntot < 10)
2713         {
2714             gmx_fatal(FARGS, "Tryig to estimtate autocorrelation time from only %d"
2715                       " points. Provide more pull data!", ntot);
2716         }
2717         snew(corr, ncorr);
2718         /* snew(corrSq,ncorr); */
2719         snew(count, ncorr);
2720         dt = window[i].dt;
2721         snew(window[i].tau, window[i].nPull);
2722         restart = static_cast<int>(opt->acTrestart/dt+0.5);
2723         if (restart == 0)
2724         {
2725             restart = 1;
2726         }
2727
2728         for (ig = 0; ig < window[i].nPull; ig++)
2729         {
2730             if (ntot != window[i].Ntot[ig])
2731             {
2732                 gmx_fatal(FARGS, "Encountered different nr of frames in different pull groups.\n"
2733                           "That should not happen. (%d and %d)\n", ntot, window[i].Ntot[ig]);
2734             }
2735             ztime = window[i].ztime[ig];
2736
2737             /* calc autocorrelation function C(t) = < [z(tau)-<z>]*[z(tau+t)-<z>]> */
2738             for (j = 0, av = 0; (j < ntot); j++)
2739             {
2740                 av += ztime[j];
2741             }
2742             av /= ntot;
2743             for (k = 0; (k < ncorr); k++)
2744             {
2745                 corr[k]  = 0.;
2746                 count[k] = 0;
2747             }
2748             for (j = 0; (j < ntot); j += restart)
2749             {
2750                 for (k = 0; (k < ncorr) && (j+k < ntot); k++)
2751                 {
2752                     tmp        = (ztime[j]-av)*(ztime[j+k]-av);
2753                     corr  [k] += tmp;
2754                     /* corrSq[k] += tmp*tmp; */
2755                     count[k]++;
2756                 }
2757             }
2758             /* divide by nr of frames for each time displacement */
2759             for (k = 0; (k < ncorr); k++)
2760             {
2761                 /* count probably = (ncorr-k+(restart-1))/restart; */
2762                 corr[k] = corr[k]/count[k];
2763                 /* variance of autocorrelation function */
2764                 /* corrSq[k]=corrSq[k]/count[k]; */
2765             }
2766             /* normalize such that corr[0] == 0 */
2767             c0 = 1./corr[0];
2768             for (k = 0; (k < ncorr); k++)
2769             {
2770                 corr[k] *= c0;
2771                 /* corrSq[k]*=c0*c0; */
2772             }
2773
2774             /* write ACFs in verbose mode */
2775             if (fpcorr)
2776             {
2777                 for (k = 0; (k < ncorr); k++)
2778                 {
2779                     fprintf(fpcorr, "%g  %g\n", k*dt, corr[k]);
2780                 }
2781                 fprintf(fpcorr, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
2782             }
2783
2784             /* esimate integrated correlation time, fitting is too unstable */
2785             tausteps = 0.5*corr[0];
2786             /* consider corr below WHAM_AC_ZERO_LIMIT as noise */
2787             for (j = 1; (j < ncorr) && (corr[j] > WHAM_AC_ZERO_LIMIT); j++)
2788             {
2789                 tausteps += corr[j];
2790             }
2791
2792             /* g = 1+2*tau, see. Ferrenberg/Swendsen, PRL 63:1195 (1989) or
2793                Kumar et al, eq. 28 ff. */
2794             window[i].tau[ig] = tausteps*dt;
2795             window[i].g[ig]   = 1+2*tausteps;
2796             /* printf("win %d, group %d, estimated correlation time = %g ps\n",i,ig,window[i].tau[ig]); */
2797         } /* ig loop */
2798         sfree(corr);
2799         sfree(count);
2800     }
2801     printf(" done\n");
2802     if (fpcorr)
2803     {
2804         xvgrclose(fpcorr);
2805     }
2806
2807     /* plot IACT along reaction coordinate */
2808     fp = xvgropen(fn, "Integrated autocorrelation times", xlabel, "IACT [ps]", opt->oenv);
2809     if (output_env_get_print_xvgr_codes(opt->oenv))
2810     {
2811         fprintf(fp, "@    s0 symbol 1\n@    s0 symbol size 0.5\n@    s0 line linestyle 0\n");
2812         fprintf(fp, "#  WIN   tau(gr1)  tau(gr2) ...\n");
2813         for (i = 0; i < nwins; i++)
2814         {
2815             fprintf(fp, "# %3d   ", i);
2816             for (ig = 0; ig < window[i].nPull; ig++)
2817             {
2818                 fprintf(fp, " %11g", window[i].tau[ig]);
2819             }
2820             fprintf(fp, "\n");
2821         }
2822     }
2823     for (i = 0; i < nwins; i++)
2824     {
2825         for (ig = 0; ig < window[i].nPull; ig++)
2826         {
2827             fprintf(fp, "%8g %8g\n", window[i].pos[ig], window[i].tau[ig]);
2828         }
2829     }
2830     if (opt->sigSmoothIact > 0.0)
2831     {
2832         printf("Smoothing autocorrelation times along reaction coordinate with Gaussian of sig = %g\n",
2833                opt->sigSmoothIact);
2834         /* smooth IACT along reaction coordinate and overwrite g=1+2tau */
2835         smoothIact(window, nwins, opt);
2836         fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
2837         if (output_env_get_print_xvgr_codes(opt->oenv))
2838         {
2839             fprintf(fp, "@    s1 symbol 1\n@    s1 symbol size 0.5\n@    s1 line linestyle 0\n");
2840             fprintf(fp, "@    s1 symbol color 2\n");
2841         }
2842         for (i = 0; i < nwins; i++)
2843         {
2844             for (ig = 0; ig < window[i].nPull; ig++)
2845             {
2846                 fprintf(fp, "%8g %8g\n", window[i].pos[ig], window[i].tausmooth[ig]);
2847             }
2848         }
2849     }
2850     xvgrclose(fp);
2851     printf("Wrote %s\n", fn);
2852 }
2853
2854 /*! \brief
2855  * compute average and sigma of each umbrella histogram
2856  */
2857 static void averageSigma(t_UmbrellaWindow *window, int nwins)
2858 {
2859     int  i, ig, ntot, k;
2860     real av, sum2, sig, diff, *ztime, nSamplesIndep;
2861
2862     for (i = 0; i < nwins; i++)
2863     {
2864         snew(window[i].aver, window[i].nPull);
2865         snew(window[i].sigma, window[i].nPull);
2866
2867         ntot = window[i].Ntot[0];
2868         for (ig = 0; ig < window[i].nPull; ig++)
2869         {
2870             ztime = window[i].ztime[ig];
2871             for (k = 0, av = 0.; k < ntot; k++)
2872             {
2873                 av += ztime[k];
2874             }
2875             av /= ntot;
2876             for (k = 0, sum2 = 0.; k < ntot; k++)
2877             {
2878                 diff  = ztime[k]-av;
2879                 sum2 += diff*diff;
2880             }
2881             sig                = std::sqrt(sum2/ntot);
2882             window[i].aver[ig] = av;
2883
2884             /* Note: This estimate for sigma is biased from the limited sampling.
2885                Correct sigma by n/(n-1) where n = number of independent
2886                samples. Only possible if IACT is known.
2887              */
2888             if (window[i].tau)
2889             {
2890                 nSamplesIndep       = window[i].N[ig]/(window[i].tau[ig]/window[i].dt);
2891                 window[i].sigma[ig] = sig * nSamplesIndep/(nSamplesIndep-1);
2892             }
2893             else
2894             {
2895                 window[i].sigma[ig] = sig;
2896             }
2897             printf("win %d, aver = %f  sig = %f\n", i, av, window[i].sigma[ig]);
2898         }
2899     }
2900 }
2901
2902
2903 /*! \brief
2904  * Use histograms to compute average force on pull group.
2905  */
2906 static void computeAverageForce(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt)
2907 {
2908     int    i, j, bins = opt->bins, k;
2909     double dz, min = opt->min, max = opt->max, displAv, temp, distance, ztot, ztot_half, w, weight;
2910     double posmirrored;
2911
2912     dz        = (max-min)/bins;
2913     ztot      = opt->max-min;
2914     ztot_half = ztot/2;
2915
2916     /* Compute average displacement from histograms */
2917     for (j = 0; j < nWindows; ++j)
2918     {
2919         snew(window[j].forceAv, window[j].nPull);
2920         for (k = 0; k < window[j].nPull; ++k)
2921         {
2922             displAv = 0.0;
2923             weight  = 0.0;
2924             for (i = 0; i < opt->bins; ++i)
2925             {
2926                 temp     = (1.0*i+0.5)*dz+min;
2927                 distance = temp - window[j].pos[k];
2928                 if (opt->bCycl)
2929                 {                                       /* in cyclic wham:             */
2930                     if (distance > ztot_half)           /*    |distance| < ztot_half   */
2931                     {
2932                         distance -= ztot;
2933                     }
2934                     else if (distance < -ztot_half)
2935                     {
2936                         distance += ztot;
2937                     }
2938                 }
2939                 w         = window[j].Histo[k][i]/window[j].g[k];
2940                 displAv  += w*distance;
2941                 weight   += w;
2942                 /* Are we near min or max? We are getting wrong forces from the histgrams since
2943                    the histograms are zero outside [min,max). Therefore, assume that the position
2944                    on the other side of the histomgram center is equally likely. */
2945                 if (!opt->bCycl)
2946                 {
2947                     posmirrored = window[j].pos[k]-distance;
2948                     if (posmirrored >= max || posmirrored < min)
2949                     {
2950                         displAv  += -w*distance;
2951                         weight   += w;
2952                     }
2953                 }
2954             }
2955             displAv  /= weight;
2956
2957             /* average force from average displacement */
2958             window[j].forceAv[k] = displAv*window[j].k[k];
2959             /* sigma from average square displacement */
2960             /* window[j].sigma  [k] = sqrt(displAv2); */
2961             /* printf("Win %d, sigma = %f\n",j,sqrt(displAv2)); */
2962         }
2963     }
2964 }
2965
2966 /*! \brief
2967  * Check if the complete reaction coordinate is covered by the histograms
2968  */
2969 static void  checkReactionCoordinateCovered(t_UmbrellaWindow *window, int nwins,
2970                                             t_UmbrellaOptions *opt)
2971 {
2972     int  i, ig, j, bins = opt->bins;
2973     bool bBoundary;
2974     real avcount = 0, z, relcount, *count;
2975     snew(count, opt->bins);
2976
2977     for (j = 0; j < opt->bins; ++j)
2978     {
2979         for (i = 0; i < nwins; i++)
2980         {
2981             for (ig = 0; ig < window[i].nPull; ig++)
2982             {
2983                 count[j] += window[i].Histo[ig][j];
2984             }
2985         }
2986         avcount += 1.0*count[j];
2987     }
2988     avcount /= bins;
2989     for (j = 0; j < bins; ++j)
2990     {
2991         relcount  = count[j]/avcount;
2992         z         = (j+0.5)*opt->dz+opt->min;
2993         bBoundary = j<bins/20 || (bins-j)>bins/20;
2994         /* check for bins with no data */
2995         if (count[j] == 0)
2996         {
2997             fprintf(stderr, "\nWARNING, no data point in bin %d (z=%g) !\n"
2998                     "You may not get a reasonable profile. Check your histograms!\n", j, z);
2999         }
3000         /* and check for poor sampling */
3001         else if (relcount < 0.005 && !bBoundary)
3002         {
3003             fprintf(stderr, "Warning, poor sampling bin %d (z=%g). Check your histograms!\n", j, z);
3004         }
3005     }
3006     sfree(count);
3007 }
3008
3009 /*! \brief Compute initial potential by integrating the average force
3010  *
3011  * This speeds up the convergence by roughly a factor of 2
3012  */
3013 static void guessPotByIntegration(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt, const char *xlabel)
3014 {
3015     int    i, j, ig, bins = opt->bins, nHist, winmin, groupmin;
3016     double dz, min = opt->min, *pot, pos, hispos, dist, diff, fAv, distmin, *f;
3017     FILE  *fp;
3018
3019     dz = (opt->max-min)/bins;
3020
3021     printf("Getting initial potential by integration.\n");
3022
3023     /* Compute average displacement from histograms */
3024     computeAverageForce(window, nWindows, opt);
3025
3026     /* Get force for each bin from all histograms in this bin, or, alternatively,
3027        if no histograms are inside this bin, from the closest histogram */
3028     snew(pot, bins);
3029     snew(f, bins);
3030     for (j = 0; j < opt->bins; ++j)
3031     {
3032         pos      = (1.0*j+0.5)*dz+min;
3033         nHist    = 0;
3034         fAv      = 0.;
3035         distmin  = 1e20;
3036         groupmin = winmin = 0;
3037         for (i = 0; i < nWindows; i++)
3038         {
3039             for (ig = 0; ig < window[i].nPull; ig++)
3040             {
3041                 hispos = window[i].pos[ig];
3042                 dist   = std::abs(hispos-pos);
3043                 /* average force within bin */
3044                 if (dist < dz/2)
3045                 {
3046                     nHist++;
3047                     fAv += window[i].forceAv[ig];
3048                 }
3049                 /* at the same time, remember closest histogram */
3050                 if (dist < distmin)
3051                 {
3052                     winmin   = i;
3053                     groupmin = ig;
3054                     distmin  = dist;
3055                 }
3056             }
3057         }
3058         /* if no histogram found in this bin, use closest histogram */
3059         if (nHist > 0)
3060         {
3061             fAv = fAv/nHist;
3062         }
3063         else
3064         {
3065             fAv = window[winmin].forceAv[groupmin];
3066         }
3067         f[j] = fAv;
3068     }
3069     for (j = 1; j < opt->bins; ++j)
3070     {
3071         pot[j] = pot[j-1] - 0.5*dz*(f[j-1]+f[j]);
3072     }
3073
3074     /* cyclic wham: linearly correct possible offset */
3075     if (opt->bCycl)
3076     {
3077         diff = (pot[bins-1]-pot[0])/(bins-1);
3078         for (j = 1; j < opt->bins; ++j)
3079         {
3080             pot[j] -= j*diff;
3081         }
3082     }
3083     if (opt->verbose)
3084     {
3085         fp = xvgropen("pmfintegrated.xvg", "PMF from force integration", xlabel, "PMF (kJ/mol)", opt->oenv);
3086         for (j = 0; j < opt->bins; ++j)
3087         {
3088             fprintf(fp, "%g  %g\n", (j+0.5)*dz+opt->min, pot[j]);
3089         }
3090         xvgrclose(fp);
3091         printf("verbose mode: wrote %s with PMF from interated forces\n", "pmfintegrated.xvg");
3092     }
3093
3094     /* get initial z=exp(-F[i]/kT) from integrated potential, where F[i] denote the free
3095        energy offsets which are usually determined by wham
3096        First: turn pot into probabilities:
3097      */
3098     for (j = 0; j < opt->bins; ++j)
3099     {
3100         pot[j] = std::exp(-pot[j]/(BOLTZ*opt->Temperature));
3101     }
3102     calc_z(pot, window, nWindows, opt, TRUE);
3103
3104     sfree(pot);
3105     sfree(f);
3106 }
3107
3108 //! Count number of words in an line
3109 static int wordcount(char *ptr)
3110 {
3111     int i, n, is[2];
3112     int cur = 0;
3113
3114     if (std::strlen(ptr) == 0)
3115     {
3116         return 0;
3117     }
3118     /* fprintf(stderr,"ptr='%s'\n",ptr); */
3119     n = 1;
3120     for (i = 0; (ptr[i] != '\0'); i++)
3121     {
3122         is[cur] = isspace(ptr[i]);
3123         if ((i > 0)  && (is[cur] && !is[1-cur]))
3124         {
3125             n++;
3126         }
3127         cur = 1-cur;
3128     }
3129     return n;
3130 }
3131
3132 /*! \brief Read input file for pull group selection (option -is)
3133  *
3134  * TO DO: ptr=fgets(...) is never freed (small memory leak)
3135  */
3136 static void readPullCoordSelection(t_UmbrellaOptions *opt, char **fnTpr, int nTpr)
3137 {
3138     FILE *fp;
3139     int   i, iline, n, len = STRLEN, temp;
3140     char *ptr = nullptr, *tmpbuf = nullptr;
3141     char  fmt[1024], fmtign[1024];
3142     int   block = 1, sizenow;
3143
3144     fp            = gmx_ffopen(opt->fnCoordSel, "r");
3145     opt->coordsel = nullptr;
3146
3147     snew(tmpbuf, len);
3148     sizenow = 0;
3149     iline   = 0;
3150     while ( (ptr = fgets3(fp, tmpbuf, &len)) != nullptr)
3151     {
3152         trim(ptr);
3153         n = wordcount(ptr);
3154
3155         if (iline >= sizenow)
3156         {
3157             sizenow += block;
3158             srenew(opt->coordsel, sizenow);
3159         }
3160         opt->coordsel[iline].n    = n;
3161         opt->coordsel[iline].nUse = 0;
3162         snew(opt->coordsel[iline].bUse, n);
3163
3164         fmtign[0] = '\0';
3165         for (i = 0; i < n; i++)
3166         {
3167             std::strcpy(fmt, fmtign);
3168             std::strcat(fmt, "%d");
3169             if (sscanf(ptr, fmt, &temp))
3170             {
3171                 opt->coordsel[iline].bUse[i] = (temp > 0);
3172                 if (opt->coordsel[iline].bUse[i])
3173                 {
3174                     opt->coordsel[iline].nUse++;
3175                 }
3176             }
3177             std::strcat(fmtign, "%*s");
3178         }
3179         iline++;
3180     }
3181     opt->nCoordsel = iline;
3182     if (nTpr != opt->nCoordsel)
3183     {
3184         gmx_fatal(FARGS, "Found %d tpr files but %d lines in %s\n", nTpr, opt->nCoordsel,
3185                   opt->fnCoordSel);
3186     }
3187
3188     printf("\nUse only these pull coordinates:\n");
3189     for (iline = 0; iline < nTpr; iline++)
3190     {
3191         printf("%s (%d of %d coordinates):", fnTpr[iline], opt->coordsel[iline].nUse, opt->coordsel[iline].n);
3192         for (i = 0; i < opt->coordsel[iline].n; i++)
3193         {
3194             if (opt->coordsel[iline].bUse[i])
3195             {
3196                 printf(" %d", i+1);
3197             }
3198         }
3199         printf("\n");
3200     }
3201     printf("\n");
3202
3203     sfree(tmpbuf);
3204 }
3205
3206 //! Boolean XOR
3207 #define WHAMBOOLXOR(a, b) ( ((!(a)) && (b)) || ((a) && (!(b))))
3208
3209 //! Number of elements in fnm (used for command line parsing)
3210 #define NFILE asize(fnm)
3211
3212 //! The main gmx wham routine
3213 int gmx_wham(int argc, char *argv[])
3214 {
3215     const char              *desc[] = {
3216         "[THISMODULE] is an analysis program that implements the Weighted",
3217         "Histogram Analysis Method (WHAM). It is intended to analyze",
3218         "output files generated by umbrella sampling simulations to ",
3219         "compute a potential of mean force (PMF).[PAR]",
3220         "",
3221         "[THISMODULE] is currently not fully up to date. It only supports pull setups",
3222         "where the first pull coordinate(s) is/are umbrella pull coordinates",
3223         "and, if multiple coordinates need to be analyzed, all used the same",
3224         "geometry and dimensions. In most cases this is not an issue.[PAR]",
3225         "At present, three input modes are supported.",
3226         "",
3227         "* With option [TT]-it[tt], the user provides a file which contains the",
3228         "  file names of the umbrella simulation run-input files ([REF].tpr[ref] files),",
3229         "  AND, with option [TT]-ix[tt], a file which contains file names of",
3230         "  the pullx [TT]mdrun[tt] output files. The [REF].tpr[ref] and pullx files must",
3231         "  be in corresponding order, i.e. the first [REF].tpr[ref] created the",
3232         "  first pullx, etc.",
3233         "* Same as the previous input mode, except that the the user",
3234         "  provides the pull force output file names ([TT]pullf.xvg[tt]) with option [TT]-if[tt].",
3235         "  From the pull force the position in the umbrella potential is",
3236         "  computed. This does not work with tabulated umbrella potentials.",
3237         "* With option [TT]-ip[tt], the user provides file names of (gzipped) [REF].pdo[ref] files, i.e.",
3238         "  the GROMACS 3.3 umbrella output files. If you have some unusual",
3239         "  reaction coordinate you may also generate your own [REF].pdo[ref] files and",
3240         "  feed them with the [TT]-ip[tt] option into to [THISMODULE]. The [REF].pdo[ref] file header",
3241         "  must be similar to the following::",
3242         "",
3243         "  # UMBRELLA      3.0",
3244         "  # Component selection: 0 0 1",
3245         "  # nSkip 1",
3246         "  # Ref. Group 'TestAtom'",
3247         "  # Nr. of pull groups 2",
3248         "  # Group 1 'GR1'  Umb. Pos. 5.0 Umb. Cons. 1000.0",
3249         "  # Group 2 'GR2'  Umb. Pos. 2.0 Umb. Cons. 500.0",
3250         "  #####",
3251         "",
3252         "  The number of pull groups, umbrella positions, force constants, and names ",
3253         "  may (of course) differ. Following the header, a time column and ",
3254         "  a data column for each pull group follows (i.e. the displacement",
3255         "  with respect to the umbrella center). Up to four pull groups are possible ",
3256         "  per [REF].pdo[ref] file at present.[PAR]",
3257         "By default, all pull coordinates found in all pullx/pullf files are used in WHAM. If only ",
3258         "some of the pull coordinates should be used, a pull coordinate selection file (option [TT]-is[tt]) can ",
3259         "be provided. The selection file must contain one line for each tpr file in tpr-files.dat.",
3260         "Each of these lines must contain one digit (0 or 1) for each pull coordinate in the tpr file. ",
3261         "Here, 1 indicates that the pull coordinate is used in WHAM, and 0 means it is omitted. Example:",
3262         "If you have three tpr files, each containing 4 pull coordinates, but only pull coordinates 1 and 2 should be ",
3263         "used, coordsel.dat looks like this::",
3264         "",
3265         "  1 1 0 0",
3266         "  1 1 0 0",
3267         "  1 1 0 0",
3268         "",
3269         "By default, the output files are::",
3270         "",
3271         "  [TT]-o[tt]      PMF output file",
3272         "  [TT]-hist[tt]   Histograms output file",
3273         "",
3274         "Always check whether the histograms sufficiently overlap.[PAR]",
3275         "The umbrella potential is assumed to be harmonic and the force constants are ",
3276         "read from the [REF].tpr[ref] or [REF].pdo[ref] files. If a non-harmonic umbrella force was applied ",
3277         "a tabulated potential can be provided with [TT]-tab[tt].",
3278         "",
3279         "WHAM options",
3280         "^^^^^^^^^^^^",
3281         "",
3282         "* [TT]-bins[tt]   Number of bins used in analysis",
3283         "* [TT]-temp[tt]   Temperature in the simulations",
3284         "* [TT]-tol[tt]    Stop iteration if profile (probability) changed less than tolerance",
3285         "* [TT]-auto[tt]   Automatic determination of boundaries",
3286         "* [TT]-min,-max[tt]   Boundaries of the profile",
3287         "",
3288         "The data points that are used to compute the profile",
3289         "can be restricted with options [TT]-b[tt], [TT]-e[tt], and [TT]-dt[tt]. ",
3290         "Adjust [TT]-b[tt] to ensure sufficient equilibration in each ",
3291         "umbrella window.[PAR]",
3292         "With [TT]-log[tt] (default) the profile is written in energy units, otherwise ",
3293         "(with [TT]-nolog[tt]) as probability. The unit can be specified with [TT]-unit[tt]. ",
3294         "With energy output, the energy in the first bin is defined to be zero. ",
3295         "If you want the free energy at a different ",
3296         "position to be zero, set [TT]-zprof0[tt] (useful with bootstrapping, see below).[PAR]",
3297         "For cyclic or periodic reaction coordinates (dihedral angle, channel PMF",
3298         "without osmotic gradient), the option [TT]-cycl[tt] is useful.",
3299         "[THISMODULE] will make use of the",
3300         "periodicity of the system and generate a periodic PMF. The first and the last bin of the",
3301         "reaction coordinate will assumed be be neighbors.[PAR]",
3302         "Option [TT]-sym[tt] symmetrizes the profile around z=0 before output, ",
3303         "which may be useful for, e.g. membranes.",
3304         "",
3305         "Parallelization",
3306         "^^^^^^^^^^^^^^^",
3307         "",
3308         "If available, the number of OpenMP threads used by gmx wham can be controlled by setting",
3309         "the [TT]OMP_NUM_THREADS[tt] environment variable.",
3310         "",
3311         "Autocorrelations",
3312         "^^^^^^^^^^^^^^^^",
3313         "",
3314         "With [TT]-ac[tt], [THISMODULE] estimates the integrated autocorrelation ",
3315         "time (IACT) [GRK]tau[grk] for each umbrella window and weights the respective ",
3316         "window with 1/[1+2*[GRK]tau[grk]/dt]. The IACTs are written ",
3317         "to the file defined with [TT]-oiact[tt]. In verbose mode, all ",
3318         "autocorrelation functions (ACFs) are written to [TT]hist_autocorr.xvg[tt]. ",
3319         "Because the IACTs can be severely underestimated in case of limited ",
3320         "sampling, option [TT]-acsig[tt] allows one to smooth the IACTs along the ",
3321         "reaction coordinate with a Gaussian ([GRK]sigma[grk] provided with [TT]-acsig[tt], ",
3322         "see output in [TT]iact.xvg[tt]). Note that the IACTs are estimated by simple ",
3323         "integration of the ACFs while the ACFs are larger 0.05.",
3324         "If you prefer to compute the IACTs by a more sophisticated (but possibly ",
3325         "less robust) method such as fitting to a double exponential, you can ",
3326         "compute the IACTs with [gmx-analyze] and provide them to [THISMODULE] with the file ",
3327         "[TT]iact-in.dat[tt] (option [TT]-iiact[tt]), which should contain one line per ",
3328         "input file ([REF].pdo[ref] or pullx/f file) and one column per pull coordinate in the respective file.",
3329         "",
3330         "Error analysis",
3331         "^^^^^^^^^^^^^^",
3332         "",
3333         "Statistical errors may be estimated with bootstrap analysis. Use it with care, ",
3334         "otherwise the statistical error may be substantially underestimated. ",
3335         "More background and examples for the bootstrap technique can be found in ",
3336         "Hub, de Groot and Van der Spoel, JCTC (2010) 6: 3713-3720.",
3337         "[TT]-nBootstrap[tt] defines the number of bootstraps (use, e.g., 100). ",
3338         "Four bootstrapping methods are supported and ",
3339         "selected with [TT]-bs-method[tt].",
3340         "",
3341         "* [TT]b-hist[tt]   Default: complete histograms are considered as independent ",
3342         "  data points, and the bootstrap is carried out by assigning random weights to the ",
3343         "  histograms (\"Bayesian bootstrap\"). Note that each point along the reaction coordinate",
3344         "  must be covered by multiple independent histograms (e.g. 10 histograms), otherwise the ",
3345         "  statistical error is underestimated.",
3346         "* [TT]hist[tt]    Complete histograms are considered as independent data points. ",
3347         "  For each bootstrap, N histograms are randomly chosen from the N given histograms ",
3348         "  (allowing duplication, i.e. sampling with replacement).",
3349         "  To avoid gaps without data along the reaction coordinate blocks of histograms ",
3350         "  ([TT]-histbs-block[tt]) may be defined. In that case, the given histograms are ",
3351         "  divided into blocks and only histograms within each block are mixed. Note that ",
3352         "  the histograms within each block must be representative for all possible histograms, ",
3353         "  otherwise the statistical error is underestimated.",
3354         "* [TT]traj[tt]  The given histograms are used to generate new random trajectories,",
3355         "  such that the generated data points are distributed according the given histograms ",
3356         "  and properly autocorrelated. The autocorrelation time (ACT) for each window must be ",
3357         "  known, so use [TT]-ac[tt] or provide the ACT with [TT]-iiact[tt]. If the ACT of all ",
3358         "  windows are identical (and known), you can also provide them with [TT]-bs-tau[tt]. ",
3359         "  Note that this method may severely underestimate the error in case of limited sampling, ",
3360         "  that is if individual histograms do not represent the complete phase space at ",
3361         "  the respective positions.",
3362         "* [TT]traj-gauss[tt]  The same as method [TT]traj[tt], but the trajectories are ",
3363         "  not bootstrapped from the umbrella histograms but from Gaussians with the average ",
3364         "  and width of the umbrella histograms. That method yields similar error estimates ",
3365         "  like method [TT]traj[tt].",
3366         "",
3367         "Bootstrapping output:",
3368         "",
3369         "* [TT]-bsres[tt]   Average profile and standard deviations",
3370         "* [TT]-bsprof[tt]  All bootstrapping profiles",
3371         "",
3372         "With [TT]-vbs[tt] (verbose bootstrapping), the histograms of each bootstrap are written, ",
3373         "and, with bootstrap method [TT]traj[tt], the cumulative distribution functions of ",
3374         "the histograms."
3375     };
3376
3377     const char              *en_unit[]       = {nullptr, "kJ", "kCal", "kT", nullptr};
3378     const char              *en_unit_label[] = {"", "E (kJ mol\\S-1\\N)", "E (kcal mol\\S-1\\N)", "E (kT)", nullptr};
3379     const char              *en_bsMethod[]   = { nullptr, "b-hist", "hist", "traj", "traj-gauss", nullptr };
3380     static t_UmbrellaOptions opt;
3381
3382     t_pargs                  pa[] = {
3383         { "-min", FALSE, etREAL, {&opt.min},
3384           "Minimum coordinate in profile"},
3385         { "-max", FALSE, etREAL, {&opt.max},
3386           "Maximum coordinate in profile"},
3387         { "-auto", FALSE, etBOOL, {&opt.bAuto},
3388           "Determine min and max automatically"},
3389         { "-bins", FALSE, etINT, {&opt.bins},
3390           "Number of bins in profile"},
3391         { "-temp", FALSE, etREAL, {&opt.Temperature},
3392           "Temperature"},
3393         { "-tol", FALSE, etREAL, {&opt.Tolerance},
3394           "Tolerance"},
3395         { "-v", FALSE, etBOOL, {&opt.verbose},
3396           "Verbose mode"},
3397         { "-b", FALSE, etREAL, {&opt.tmin},
3398           "First time to analyse (ps)"},
3399         { "-e", FALSE, etREAL, {&opt.tmax},
3400           "Last time to analyse (ps)"},
3401         { "-dt", FALSE, etREAL, {&opt.dt},
3402           "Analyse only every dt ps"},
3403         { "-histonly", FALSE, etBOOL, {&opt.bHistOnly},
3404           "Write histograms and exit"},
3405         { "-boundsonly", FALSE, etBOOL, {&opt.bBoundsOnly},
3406           "Determine min and max and exit (with [TT]-auto[tt])"},
3407         { "-log", FALSE, etBOOL, {&opt.bLog},
3408           "Calculate the log of the profile before printing"},
3409         { "-unit", FALSE,  etENUM, {en_unit},
3410           "Energy unit in case of log output" },
3411         { "-zprof0", FALSE, etREAL, {&opt.zProf0},
3412           "Define profile to 0.0 at this position (with [TT]-log[tt])"},
3413         { "-cycl", FALSE, etBOOL, {&opt.bCycl},
3414           "Create cyclic/periodic profile. Assumes min and max are the same point."},
3415         { "-sym", FALSE, etBOOL, {&opt.bSym},
3416           "Symmetrize profile around z=0"},
3417         { "-hist-eq", FALSE, etBOOL, {&opt.bHistEq},
3418           "HIDDENEnforce equal weight for all histograms. (Non-Weighed-HAM)"},
3419         { "-ac", FALSE, etBOOL, {&opt.bCalcTauInt},
3420           "Calculate integrated autocorrelation times and use in wham"},
3421         { "-acsig", FALSE, etREAL, {&opt.sigSmoothIact},
3422           "Smooth autocorrelation times along reaction coordinate with Gaussian of this [GRK]sigma[grk]"},
3423         { "-ac-trestart", FALSE, etREAL, {&opt.acTrestart},
3424           "When computing autocorrelation functions, restart computing every .. (ps)"},
3425         { "-acred", FALSE, etBOOL, {&opt.bAllowReduceIact},
3426           "HIDDENWhen smoothing the ACTs, allows one to reduce ACTs. Otherwise, only increase ACTs "
3427           "during smoothing"},
3428         { "-nBootstrap", FALSE,  etINT, {&opt.nBootStrap},
3429           "nr of bootstraps to estimate statistical uncertainty (e.g., 200)" },
3430         { "-bs-method", FALSE,  etENUM, {en_bsMethod},
3431           "Bootstrap method" },
3432         { "-bs-tau", FALSE, etREAL, {&opt.tauBootStrap},
3433           "Autocorrelation time (ACT) assumed for all histograms. Use option [TT]-ac[tt] if ACT is unknown."},
3434         { "-bs-seed", FALSE, etINT, {&opt.bsSeed},
3435           "Seed for bootstrapping. (-1 = use time)"},
3436         { "-histbs-block", FALSE, etINT, {&opt.histBootStrapBlockLength},
3437           "When mixing histograms only mix within blocks of [TT]-histbs-block[tt]."},
3438         { "-vbs", FALSE, etBOOL, {&opt.bs_verbose},
3439           "Verbose bootstrapping. Print the CDFs and a histogram file for each bootstrap."},
3440         { "-stepout", FALSE, etINT, {&opt.stepchange},
3441           "HIDDENWrite maximum change every ... (set to 1 with [TT]-v[tt])"},
3442         { "-updateContr", FALSE, etINT, {&opt.stepUpdateContrib},
3443           "HIDDENUpdate table with significan contributions to WHAM every ... iterations"},
3444     };
3445
3446     t_filenm                 fnm[] = {
3447         { efDAT, "-ix", "pullx-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs           */
3448         { efDAT, "-if", "pullf-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs           */
3449         { efDAT, "-it", "tpr-files", ffOPTRD},   /* wham input: tprs                           */
3450         { efDAT, "-ip", "pdo-files", ffOPTRD},   /* wham input: pdo files (gmx3 style)         */
3451         { efDAT, "-is", "coordsel", ffOPTRD},    /* input: select pull coords to use           */
3452         { efXVG, "-o", "profile", ffWRITE },     /* output file for profile                     */
3453         { efXVG, "-hist", "histo", ffWRITE},     /* output file for histograms                  */
3454         { efXVG, "-oiact", "iact", ffOPTWR},     /* writing integrated autocorrelation times    */
3455         { efDAT, "-iiact", "iact-in", ffOPTRD},  /* reading integrated autocorrelation times   */
3456         { efXVG, "-bsres", "bsResult", ffOPTWR}, /* average and errors of bootstrap analysis    */
3457         { efXVG, "-bsprof", "bsProfs", ffOPTWR}, /* output file for bootstrap profiles          */
3458         { efDAT, "-tab", "umb-pot", ffOPTRD},    /* Tabulated umbrella potential (if not harmonic) */
3459     };
3460
3461     int                      i, j, l, nfiles, nwins, nfiles2;
3462     t_UmbrellaHeader         header;
3463     t_UmbrellaWindow       * window = nullptr;
3464     double                  *profile, maxchange = 1e20;
3465     gmx_bool                 bMinSet, bMaxSet, bAutoSet, bExact = FALSE;
3466     char                   **fninTpr, **fninPull, **fninPdo;
3467     const char              *fnPull;
3468     FILE                    *histout, *profout;
3469     char                     xlabel[STRLEN], ylabel[256], title[256];
3470
3471     opt.bins      = 200;
3472     opt.verbose   = FALSE;
3473     opt.bHistOnly = FALSE;
3474     opt.bCycl     = FALSE;
3475     opt.tmin      = 50;
3476     opt.tmax      = 1e20;
3477     opt.dt        = 0.0;
3478     opt.min       = 0;
3479     opt.max       = 0;
3480     opt.bAuto     = TRUE;
3481     opt.nCoordsel = 0;
3482     opt.coordsel  = nullptr;
3483
3484     /* bootstrapping stuff */
3485     opt.nBootStrap               = 0;
3486     opt.bsMethod                 = bsMethod_hist;
3487     opt.tauBootStrap             = 0.0;
3488     opt.bsSeed                   = -1;
3489     opt.histBootStrapBlockLength = 8;
3490     opt.bs_verbose               = FALSE;
3491
3492     opt.bLog                  = TRUE;
3493     opt.unit                  = en_kJ;
3494     opt.zProf0                = 0.;
3495     opt.Temperature           = 298;
3496     opt.Tolerance             = 1e-6;
3497     opt.bBoundsOnly           = FALSE;
3498     opt.bSym                  = FALSE;
3499     opt.bCalcTauInt           = FALSE;
3500     opt.sigSmoothIact         = 0.0;
3501     opt.bAllowReduceIact      = TRUE;
3502     opt.bInitPotByIntegration = TRUE;
3503     opt.acTrestart            = 1.0;
3504     opt.stepchange            = 100;
3505     opt.stepUpdateContrib     = 100;
3506
3507     if (!parse_common_args(&argc, argv, 0,
3508                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &opt.oenv))
3509     {
3510         return 0;
3511     }
3512
3513     opt.unit     = nenum(en_unit);
3514     opt.bsMethod = nenum(en_bsMethod);
3515
3516     opt.bProf0Set = opt2parg_bSet("-zprof0",  asize(pa), pa);
3517
3518     opt.bTab         = opt2bSet("-tab", NFILE, fnm);
3519     opt.bPdo         = opt2bSet("-ip", NFILE, fnm);
3520     opt.bTpr         = opt2bSet("-it", NFILE, fnm);
3521     opt.bPullx       = opt2bSet("-ix", NFILE, fnm);
3522     opt.bPullf       = opt2bSet("-if", NFILE, fnm);
3523     opt.bTauIntGiven = opt2bSet("-iiact", NFILE, fnm);
3524     if  (opt.bTab && opt.bPullf)
3525     {
3526         gmx_fatal(FARGS, "Force input does not work with tabulated potentials. "
3527                   "Provide pullx.xvg or pdo files!");
3528     }
3529
3530     if (!opt.bPdo && !WHAMBOOLXOR(opt.bPullx, opt.bPullf))
3531     {
3532         gmx_fatal(FARGS, "Give either pullx (-ix) OR pullf (-if) data. Not both.");
3533     }
3534     if (!opt.bPdo && !(opt.bTpr || opt.bPullf || opt.bPullx))
3535     {
3536         gmx_fatal(FARGS, "gmx wham supports three input modes, pullx, pullf, or pdo file input."
3537                   "\n\n Check gmx wham -h !");
3538     }
3539
3540     opt.fnPdo      = opt2fn("-ip", NFILE, fnm);
3541     opt.fnTpr      = opt2fn("-it", NFILE, fnm);
3542     opt.fnPullf    = opt2fn("-if", NFILE, fnm);
3543     opt.fnPullx    = opt2fn("-ix", NFILE, fnm);
3544     opt.fnCoordSel = opt2fn_null("-is", NFILE, fnm);
3545
3546     bMinSet  = opt2parg_bSet("-min",  asize(pa), pa);
3547     bMaxSet  = opt2parg_bSet("-max",  asize(pa), pa);
3548     bAutoSet = opt2parg_bSet("-auto",  asize(pa), pa);
3549     if ( (bMinSet || bMaxSet) && bAutoSet)
3550     {
3551         gmx_fatal(FARGS, "With -auto, do not give -min or -max\n");
3552     }
3553
3554     if ( (bMinSet && !bMaxSet) || (!bMinSet && bMaxSet))
3555     {
3556         gmx_fatal(FARGS, "When giving -min, you must give -max (and vice versa), too\n");
3557     }
3558
3559     if (bMinSet && opt.bAuto)
3560     {
3561         printf("Note: min and max given, switching off -auto.\n");
3562         opt.bAuto = FALSE;
3563     }
3564
3565     if (opt.bTauIntGiven && opt.bCalcTauInt)
3566     {
3567         gmx_fatal(FARGS, "Either read (option -iiact) or calculate (option -ac) the\n"
3568                   "the autocorrelation times. Not both.");
3569     }
3570
3571     if (opt.tauBootStrap > 0.0 && opt2parg_bSet("-ac", asize(pa), pa))
3572     {
3573         gmx_fatal(FARGS, "Either compute autocorrelation times (ACTs) (option -ac) or "
3574                   "provide it with -bs-tau for bootstrapping. Not Both.\n");
3575     }
3576     if (opt.tauBootStrap > 0.0 && opt2bSet("-iiact", NFILE, fnm))
3577     {
3578         gmx_fatal(FARGS, "Either provide autocorrelation times (ACTs) with file iact-in.dat "
3579                   "(option -iiact) or define all ACTs with -bs-tau for bootstrapping\n. Not Both.");
3580     }
3581
3582     /* Reading gmx4/gmx5 pull output and tpr files */
3583     if (opt.bTpr || opt.bPullf || opt.bPullx)
3584     {
3585         read_wham_in(opt.fnTpr, &fninTpr, &nfiles, &opt);
3586
3587         fnPull = opt.bPullf ? opt.fnPullf : opt.fnPullx;
3588         read_wham_in(fnPull, &fninPull, &nfiles2, &opt);
3589         printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n",
3590                nfiles, nfiles2, opt.bPullf ? "force" : "position", opt.fnTpr, fnPull);
3591         if (nfiles != nfiles2)
3592         {
3593             gmx_fatal(FARGS, "Found %d file names in %s, but %d in %s\n", nfiles,
3594                       opt.fnTpr, nfiles2, fnPull);
3595         }
3596
3597         /* Read file that selects the pull group to be used */
3598         if (opt.fnCoordSel != nullptr)
3599         {
3600             readPullCoordSelection(&opt, fninTpr, nfiles);
3601         }
3602
3603         window = initUmbrellaWindows(nfiles);
3604         read_tpr_pullxf_files(fninTpr, fninPull, nfiles, &header, window, &opt);
3605     }
3606     else
3607     {   /* reading pdo files */
3608         if  (opt.fnCoordSel != nullptr)
3609         {
3610             gmx_fatal(FARGS, "Reading a -is file is not supported with PDO input files.\n"
3611                       "Use awk or a similar tool to pick the required pull groups from your PDO files\n");
3612         }
3613         read_wham_in(opt.fnPdo, &fninPdo, &nfiles, &opt);
3614         printf("Found %d pdo files in %s\n", nfiles, opt.fnPdo);
3615         window = initUmbrellaWindows(nfiles);
3616         read_pdo_files(fninPdo, nfiles, &header, window, &opt);
3617     }
3618
3619     /* It is currently assumed that all pull coordinates have the same geometry, so they also have the same coordinate units.
3620        We can therefore get the units for the xlabel from the first coordinate. */
3621     sprintf(xlabel, "\\xx\\f{} (%s)", header.pcrd[0].coord_unit);
3622
3623     nwins = nfiles;
3624
3625     /* enforce equal weight for all histograms? */
3626     if (opt.bHistEq)
3627     {
3628         enforceEqualWeights(window, nwins);
3629     }
3630
3631     /* write histograms */
3632     histout = xvgropen(opt2fn("-hist", NFILE, fnm), "Umbrella histograms",
3633                        xlabel, "count", opt.oenv);
3634     for (l = 0; l < opt.bins; ++l)
3635     {
3636         fprintf(histout, "%e\t", (l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
3637         for (i = 0; i < nwins; ++i)
3638         {
3639             for (j = 0; j < window[i].nPull; ++j)
3640             {
3641                 fprintf(histout, "%e\t", window[i].Histo[j][l]);
3642             }
3643         }
3644         fprintf(histout, "\n");
3645     }
3646     xvgrclose(histout);
3647     printf("Wrote %s\n", opt2fn("-hist", NFILE, fnm));
3648     if (opt.bHistOnly)
3649     {
3650         printf("Wrote histograms to %s, now exiting.\n", opt2fn("-hist", NFILE, fnm));
3651         return 0;
3652     }
3653
3654     /* Using tabulated umbrella potential */
3655     if (opt.bTab)
3656     {
3657         setup_tab(opt2fn("-tab", NFILE, fnm), &opt);
3658     }
3659
3660     /* Integrated autocorrelation times provided ? */
3661     if (opt.bTauIntGiven)
3662     {
3663         readIntegratedAutocorrelationTimes(window, nwins, opt2fn("-iiact", NFILE, fnm));
3664     }
3665
3666     /* Compute integrated autocorrelation times */
3667     if (opt.bCalcTauInt)
3668     {
3669         calcIntegratedAutocorrelationTimes(window, nwins, &opt, opt2fn("-oiact", NFILE, fnm), xlabel);
3670     }
3671
3672     /* calc average and sigma for each histogram
3673        (maybe required for bootstrapping. If not, this is fast anyhow) */
3674     if (opt.nBootStrap && opt.bsMethod == bsMethod_trajGauss)
3675     {
3676         averageSigma(window, nwins);
3677     }
3678
3679     /* Get initial potential by simple integration */
3680     if (opt.bInitPotByIntegration)
3681     {
3682         guessPotByIntegration(window, nwins, &opt, xlabel);
3683     }
3684
3685     /* Check if complete reaction coordinate is covered */
3686     checkReactionCoordinateCovered(window, nwins, &opt);
3687
3688     /* Calculate profile */
3689     snew(profile, opt.bins);
3690     if (opt.verbose)
3691     {
3692         opt.stepchange = 1;
3693     }
3694     i = 0;
3695     do
3696     {
3697         if ( (i%opt.stepUpdateContrib) == 0)
3698         {
3699             setup_acc_wham(profile, window, nwins, &opt);
3700         }
3701         if (maxchange < opt.Tolerance)
3702         {
3703             bExact = TRUE;
3704             /* if (opt.verbose) */
3705             printf("Switched to exact iteration in iteration %d\n", i);
3706         }
3707         calc_profile(profile, window, nwins, &opt, bExact);
3708         if (((i%opt.stepchange) == 0 || i == 1) && i != 0)
3709         {
3710             printf("\t%4d) Maximum change %e\n", i, maxchange);
3711         }
3712         i++;
3713     }
3714     while ( (maxchange = calc_z(profile, window, nwins, &opt, bExact)) > opt.Tolerance || !bExact);
3715     printf("Converged in %d iterations. Final maximum change %g\n", i, maxchange);
3716
3717     /* calc error from Kumar's formula */
3718     /* Unclear how the error propagates along reaction coordinate, therefore
3719        commented out  */
3720     /* calc_error_kumar(profile,window, nwins,&opt); */
3721
3722     /* Write profile in energy units? */
3723     if (opt.bLog)
3724     {
3725         prof_normalization_and_unit(profile, &opt);
3726         std::strcpy(ylabel, en_unit_label[opt.unit]);
3727         std::strcpy(title, "Umbrella potential");
3728     }
3729     else
3730     {
3731         std::strcpy(ylabel, "Density of states");
3732         std::strcpy(title, "Density of states");
3733     }
3734
3735     /* symmetrize profile around z=0? */
3736     if (opt.bSym)
3737     {
3738         symmetrizeProfile(profile, &opt);
3739     }
3740
3741     /* write profile or density of states */
3742     profout = xvgropen(opt2fn("-o", NFILE, fnm), title, xlabel, ylabel, opt.oenv);
3743     for (i = 0; i < opt.bins; ++i)
3744     {
3745         fprintf(profout, "%e\t%e\n", (i+0.5)/opt.bins*(opt.max-opt.min)+opt.min, profile[i]);
3746     }
3747     xvgrclose(profout);
3748     printf("Wrote %s\n", opt2fn("-o", NFILE, fnm));
3749
3750     /* Bootstrap Method */
3751     if (opt.nBootStrap)
3752     {
3753         do_bootstrapping(opt2fn("-bsres", NFILE, fnm), opt2fn("-bsprof", NFILE, fnm),
3754                          opt2fn("-hist", NFILE, fnm),
3755                          xlabel, ylabel, profile, window, nwins, &opt);
3756     }
3757
3758     sfree(profile);
3759     freeUmbrellaWindows(window, nfiles);
3760
3761     printf("\nIn case you use results from gmx wham for a publication, please cite:\n");
3762     please_cite(stdout, "Hub2010");
3763
3764     return 0;
3765 }