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