src/gromacs/gmxana/gmx_wham.cpp

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