src/gromacs/gmxana/gmx_wham.cpp

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