src/gromacs/gmxana/gmx_wham.cpp

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