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