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