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