1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*- */
3 * This file is part of the GROMACS molecular simulation package.
5 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
6 * Copyright (c) 2001-2004, The GROMACS development team,
7 * check out http://www.gromacs.org for more information.
8 * Copyright (c) 2012,2013, by the GROMACS development team, led by
9 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
10 * others, as listed in the AUTHORS file in the top-level source
11 * directory and at http://www.gromacs.org.
13 * GROMACS is free software; you can redistribute it and/or
14 * modify it under the terms of the GNU Lesser General Public License
15 * as published by the Free Software Foundation; either version 2.1
16 * of the License, or (at your option) any later version.
18 * GROMACS is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21 * Lesser General Public License for more details.
23 * You should have received a copy of the GNU Lesser General Public
24 * License along with GROMACS; if not, see
25 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
26 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
28 * If you want to redistribute modifications to GROMACS, please
29 * consider that scientific software is very special. Version
30 * control is crucial - bugs must be traceable. We will be happy to
31 * consider code for inclusion in the official distribution, but
32 * derived work must not be called official GROMACS. Details are found
33 * in the README & COPYING files - if they are missing, get the
34 * official version at http://www.gromacs.org.
36 * To help us fund GROMACS development, we humbly ask that you cite
37 * the research papers on the package. Check out http://www.gromacs.org.
53 #include "gmx_random.h"
63 #define WHAM_MAXFILELEN 2048
65 /* enum for energy units */
67 enSel, en_kJ, en_kCal, en_kT, enNr
69 /* enum for type of input files (pdos, tpr, or pullf) */
71 whamin_unknown, whamin_tpr, whamin_pullxf, whamin_pdo
73 /* enum for bootstrapping method (
74 - bootstrap complete histograms with continuous weights (Bayesian bootstrap)
75 - bootstrap complete histograms
76 - bootstrap trajectories from given umbrella histograms
77 - bootstrap trajectories from Gaussian with mu/sigam computed from
78 the respective histogram
80 ********************************************************************
81 FOR MORE DETAILS ON THE BOOTSTRAP METHODS (INCLUDING EXAMPLES), SEE
82 JS Hub, BL de Groot, D van der Spoel
83 [TT]g_wham[tt] - A free weighted histogram analysis implementation including
84 robust error and autocorrelation estimates,
85 J Chem Theory Comput, accepted (2010)
86 ********************************************************************
89 bsMethod_unknown, bsMethod_BayesianHist, bsMethod_hist,
90 bsMethod_traj, bsMethod_trajGauss
96 /* umbrella with pull code of gromacs 4.x */
97 int npullgrps; /* nr of pull groups in tpr file */
98 int pull_geometry; /* such as distance, position */
99 ivec pull_dim; /* pull dimension with geometry distance */
100 int pull_ndim; /* nr of pull_dim != 0 */
101 real *k; /* force constants in tpr file */
102 rvec *init_dist; /* reference displacements */
103 real *umbInitDist; /* reference displacement in umbrella direction */
105 /* From here, old pdo stuff */
111 char PullName[4][256];
113 double UmbCons[4][3];
118 int nPull; /* nr of pull groups in this pdo or pullf/x file */
119 double **Histo, **cum; /* nPull histograms and nPull cumulative distr. funct */
120 int nBin; /* nr of bins. identical to opt->bins */
121 double *k; /* force constants for the nPull groups */
122 double *pos; /* umbrella positions for the nPull groups */
123 double *z; /* z=(-Fi/kT) for the nPull groups. These values are
124 iteratively computed during wham */
125 double *N, *Ntot; /* nr of data points in nPull histograms. N and Ntot
126 only differ if bHistEq==TRUE */
128 double *g, *tau, *tausmooth; /* g = 1 + 2*tau[int]/dt where tau is the integrated
129 autocorrelation time. Compare, e.g.
130 Ferrenberg/Swendsen, PRL 63:1195 (1989)
131 Kumar et al, J Comp Chem 13, 1011-1021 (1992), eq. 28 */
133 double dt; /* timestep in the input data. Can be adapted with
135 gmx_bool **bContrib; /* TRUE, if any data point of the histogram is within min
136 and max, otherwise FALSE. */
137 real **ztime; /* input data z(t) as a function of time. Required to
139 real *forceAv; /* average force estimated from average displacement, fAv=dzAv*k
140 Used for integration to guess the potential. */
141 real *aver, *sigma; /* average and stddev of histograms */
142 double *bsWeight; /* for bootstrapping complete histograms with continuous weights */
149 const char *fnTpr, *fnPullf;
150 const char *fnPdo, *fnPullx; /* file names of input */
151 gmx_bool bTpr, bPullf, bPdo, bPullx; /* input file types given? */
152 real tmin, tmax, dt; /* only read input within tmin and tmax with dt */
154 gmx_bool bInitPotByIntegration; /* before WHAM, guess potential by force integration. Yields
155 1.5 to 2 times faster convergence */
156 int stepUpdateContrib; /* update contribution table every ... iterations. Accelerates
159 /* BASIC WHAM OPTIONS */
160 int bins; /* nr of bins, min, max, and dz of profile */
162 real Temperature, Tolerance; /* temperature, converged when probability changes less
164 gmx_bool bCycl; /* generate cyclic (periodic) PMF */
167 gmx_bool bLog; /* energy output (instead of probability) for profile */
168 int unit; /* unit for PMF output kJ/mol or kT or kCal/mol */
169 gmx_bool bSym; /* symmetrize PMF around z=0 after WHAM, useful for
171 real zProf0; /* after wham, set prof to zero at this z-position
172 When bootstrapping, set zProf0 to a "stable" reference
174 gmx_bool bProf0Set; /* setting profile to 0 at zProf0? */
176 gmx_bool bBoundsOnly, bHistOnly; /* determine min and max, or write histograms and exit */
177 gmx_bool bAuto; /* determine min and max automatically but do not exit */
179 gmx_bool verbose; /* more noisy wham mode */
180 int stepchange; /* print maximum change in prof after how many interations */
181 output_env_t oenv; /* xvgr options */
183 /* AUTOCORRELATION STUFF */
184 gmx_bool bTauIntGiven, bCalcTauInt; /* IACT given or should be calculated? */
185 real sigSmoothIact; /* sigma of Gaussian to smooth ACTs */
186 gmx_bool bAllowReduceIact; /* Allow to reduce ACTs during smoothing. Otherwise
187 ACT are only increased during smoothing */
188 real acTrestart; /* when computing ACT, time between restarting points */
189 gmx_bool bHistEq; /* Enforce the same weight for each umbella window, that is
190 calculate with the same number of data points for
191 each window. That can be reasonable, if the histograms
192 have different length, but due to autocorrelation,
193 a longer simulation should not have larger weightin wham. */
195 /* BOOTSTRAPPING STUFF */
196 int nBootStrap; /* nr of bootstraps (50 is usually enough) */
197 int bsMethod; /* if == bsMethod_hist, consider complete histograms as independent
198 data points and, hence, only mix complete histograms.
199 if == bsMethod_BayesianHist, consider complete histograms
200 as independent data points, but assign random weights
201 to the histograms during the bootstrapping ("Bayesian bootstrap")
202 In case of long correlations (e.g., inside a channel), these
203 will yield a more realistic error.
204 if == bsMethod_traj(Gauss), generate synthetic histograms
206 histogram by generating an autocorrelated random sequence
207 that is distributed according to the respective given
208 histogram. With bsMethod_trajGauss, bootstrap from a Gaussian
209 (instead of from the umbrella histogram) to generate a new
212 real tauBootStrap; /* autocorrelation time (ACT) used to generate synthetic
213 histograms. If ==0, use calculated ACF */
214 int histBootStrapBlockLength; /* when mixing histograms, mix only histograms withing blocks
215 long the reaction coordinate xi. Avoids gaps along xi. */
216 int bsSeed; /* random seed for bootstrapping */
217 gmx_bool bs_verbose; /* Write cumulative distribution functions (CDFs) of histograms
218 and write the generated histograms for each bootstrap */
220 /* tabulated umbrella potential stuff */
222 double *tabX, *tabY, tabMin, tabMax, tabDz;
225 gmx_rng_t rng; /* gromacs random number generator */
229 t_UmbrellaWindow * initUmbrellaWindows(int nwin)
231 t_UmbrellaWindow *win;
234 for (i = 0; i < nwin; i++)
236 win[i].Histo = win[i].cum = 0;
237 win[i].k = win[i].pos = win[i].z = 0;
238 win[i].N = win[i].Ntot = 0;
239 win[i].g = win[i].tau = win[i].tausmooth = 0;
243 win[i].aver = win[i].sigma = 0;
249 void freeUmbrellaWindows(t_UmbrellaWindow *win, int nwin)
252 for (i = 0; i < nwin; i++)
256 for (j = 0; j < win[i].nPull; j++)
258 sfree(win[i].Histo[j]);
263 for (j = 0; j < win[i].nPull; j++)
265 sfree(win[i].cum[j]);
270 for (j = 0; j < win[i].nPull; j++)
272 sfree(win[i].bContrib[j]);
284 sfree(win[i].tausmooth);
285 sfree(win[i].bContrib);
287 sfree(win[i].forceAv);
290 sfree(win[i].bsWeight);
295 /* Read and setup tabulated umbrella potential */
296 void setup_tab(const char *fn, t_UmbrellaOptions *opt)
301 printf("Setting up tabulated potential from file %s\n", fn);
302 nl = read_xvg(fn, &y, &ny);
306 gmx_fatal(FARGS, "Found %d columns in %s. Expected 2.\n", ny, fn);
308 opt->tabMin = y[0][0];
309 opt->tabMax = y[0][nl-1];
310 opt->tabDz = (opt->tabMax-opt->tabMin)/(nl-1);
313 gmx_fatal(FARGS, "The tabulated potential in %s must be provided in \n"
314 "ascending z-direction", fn);
316 for (i = 0; i < nl-1; i++)
318 if (fabs(y[0][i+1]-y[0][i]-opt->tabDz) > opt->tabDz/1e6)
320 gmx_fatal(FARGS, "z-values in %s are not equally spaced.\n", ny, fn);
325 for (i = 0; i < nl; i++)
327 opt->tabX[i] = y[0][i];
328 opt->tabY[i] = y[1][i];
330 printf("Found equally spaced tabulated potential from %g to %g, spacing %g\n",
331 opt->tabMin, opt->tabMax, opt->tabDz);
334 void read_pdo_header(FILE * file, t_UmbrellaHeader * header, t_UmbrellaOptions *opt)
337 char Buffer0[256], Buffer1[256], Buffer2[256], Buffer3[256], Buffer4[256];
341 if (fgets(line, 2048, file) == NULL)
343 gmx_fatal(FARGS, "Error reading header from pdo file\n");
345 sscanf(line, "%s%s%s", Buffer0, Buffer1, Buffer2);
346 if (strcmp(Buffer1, "UMBRELLA"))
348 gmx_fatal(FARGS, "This does not appear to be a valid pdo file. Found %s, expected %s\n"
349 "(Found in first line: `%s')\n",
350 Buffer1, "UMBRELLA", line);
352 if (strcmp(Buffer2, "3.0"))
354 gmx_fatal(FARGS, "This does not appear to be a version 3.0 pdo file");
358 if (fgets(line, 2048, file) == NULL)
360 gmx_fatal(FARGS, "Error reading header from pdo file\n");
362 sscanf(line, "%s%s%s%d%d%d", Buffer0, Buffer1, Buffer2,
363 &(header->Dims[0]), &(header->Dims[1]), &(header->Dims[2]));
365 /* printf("%d %d %d\n", header->Dims[0],header->Dims[1],header->Dims[2]); */
367 header->nDim = header->Dims[0] + header->Dims[1] + header->Dims[2];
368 if (header->nDim != 1)
370 gmx_fatal(FARGS, "Currently only supports one dimension");
374 if (fgets(line, 2048, file) == NULL)
376 gmx_fatal(FARGS, "Error reading header from pdo file\n");
378 sscanf(line, "%s%s%d", Buffer0, Buffer1, &(header->nSkip));
381 if (fgets(line, 2048, file) == NULL)
383 gmx_fatal(FARGS, "Error reading header from pdo file\n");
385 sscanf(line, "%s%s%s%s", Buffer0, Buffer1, Buffer2, header->Reference);
388 if (fgets(line, 2048, file) == NULL)
390 gmx_fatal(FARGS, "Error reading header from pdo file\n");
392 sscanf(line, "%s%s%s%s%s%d", Buffer0, Buffer1, Buffer2, Buffer3, Buffer4, &(header->nPull));
396 printf("\tFound nPull=%d , nSkip=%d, ref=%s\n", header->nPull, header->nSkip,
400 for (i = 0; i < header->nPull; ++i)
402 if (fgets(line, 2048, file) == NULL)
404 gmx_fatal(FARGS, "Error reading header from pdo file\n");
406 sscanf(line, "%*s%*s%*s%s%*s%*s%lf%*s%*s%lf", header->PullName[i],
407 &(header->UmbPos[i][0]), &(header->UmbCons[i][0]));
410 printf("\tpullgroup %d, pullname = %s, UmbPos = %g, UmbConst = %g\n",
411 i, header->PullName[i], header->UmbPos[i][0], header->UmbCons[i][0]);
415 if (fgets(line, 2048, file) == NULL)
417 gmx_fatal(FARGS, "Cannot read from file\n");
419 sscanf(line, "%s", Buffer3);
420 if (strcmp(Buffer3, "#####") != 0)
422 gmx_fatal(FARGS, "Expected '#####', found %s. Hick.\n", Buffer3);
427 static char *fgets3(FILE *fp, char ptr[], int *len)
432 if (fgets(ptr, *len-1, fp) == NULL)
437 while ((strchr(ptr, '\n') == NULL) && (!feof(fp)))
439 /* This line is longer than len characters, let's increase len! */
443 if (fgets(p-1, STRLEN, fp) == NULL)
449 if (ptr[slen-1] == '\n')
457 void read_pdo_data(FILE * file, t_UmbrellaHeader * header,
458 int fileno, t_UmbrellaWindow * win,
459 t_UmbrellaOptions *opt,
460 gmx_bool bGetMinMax, real *mintmp, real *maxtmp)
462 int i, inttemp, bins, count, ntot;
463 real min, max, minfound = 1e20, maxfound = -1e20;
464 double temp, time, time0 = 0, dt;
466 t_UmbrellaWindow * window = 0;
467 gmx_bool timeok, dt_ok = 1;
468 char *tmpbuf = 0, fmt[256], fmtign[256];
469 int len = STRLEN, dstep = 1;
470 const int blocklen = 4096;
480 /* Need to alocate memory and set up structure */
481 window->nPull = header->nPull;
484 snew(window->Histo, window->nPull);
485 snew(window->z, window->nPull);
486 snew(window->k, window->nPull);
487 snew(window->pos, window->nPull);
488 snew(window->N, window->nPull);
489 snew(window->Ntot, window->nPull);
490 snew(window->g, window->nPull);
491 snew(window->bsWeight, window->nPull);
493 window->bContrib = 0;
495 if (opt->bCalcTauInt)
497 snew(window->ztime, window->nPull);
503 snew(lennow, window->nPull);
505 for (i = 0; i < window->nPull; ++i)
508 window->bsWeight[i] = 1.;
509 snew(window->Histo[i], bins);
510 window->k[i] = header->UmbCons[i][0];
511 window->pos[i] = header->UmbPos[i][0];
515 if (opt->bCalcTauInt)
517 window->ztime[i] = 0;
521 /* Done with setup */
527 min = max = bins = 0; /* Get rid of warnings */
532 while ( (ptr = fgets3(file, tmpbuf, &len)) != NULL)
536 if (ptr[0] == '#' || strlen(ptr) < 2)
541 /* Initiate format string */
543 strcat(fmtign, "%*s");
545 sscanf(ptr, "%lf", &time); /* printf("Time %f\n",time); */
546 /* Round time to fs */
547 time = 1.0/1000*( (int) (time*1000+0.5) );
549 /* get time step of pdo file */
559 dstep = (int)(opt->dt/dt+0.5);
567 window->dt = dt*dstep;
572 dt_ok = ((count-1)%dstep == 0);
573 timeok = (dt_ok && time >= opt->tmin && time <= opt->tmax);
575 printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
576 time,opt->tmin, opt->tmax, dt_ok,timeok); */
580 for (i = 0; i < header->nPull; ++i)
583 strcat(fmt, "%lf"); /* Creating a format stings such as "%*s...%*s%lf" */
584 strcat(fmtign, "%*s"); /* ignoring one more entry in the next loop */
585 if (sscanf(ptr, fmt, &temp))
587 temp += header->UmbPos[i][0];
601 if (opt->bCalcTauInt)
603 /* save time series for autocorrelation analysis */
604 ntot = window->Ntot[i];
605 if (ntot >= lennow[i])
607 lennow[i] += blocklen;
608 srenew(window->ztime[i], lennow[i]);
610 window->ztime[i][ntot] = temp;
625 else if (inttemp >= bins)
631 if (inttemp >= 0 && inttemp < bins)
633 window->Histo[i][inttemp] += 1.;
641 if (time > opt->tmax)
645 printf("time %f larger than tmax %f, stop reading pdo file\n", time, opt->tmax);
661 void enforceEqualWeights(t_UmbrellaWindow * window, int nWindows)
663 int i, k, j, NEnforced;
666 NEnforced = window[0].Ntot[0];
667 printf("\nFound -hist-eq. Enforcing equal weights for all histograms, \ni.e. doing a "
668 "non-weighted histogram analysis method. Ndata = %d\n", NEnforced);
669 /* enforce all histograms to have the same weight as the very first histogram */
671 for (j = 0; j < nWindows; ++j)
673 for (k = 0; k < window[j].nPull; ++k)
675 ratio = 1.0*NEnforced/window[j].Ntot[k];
676 for (i = 0; i < window[0].nBin; ++i)
678 window[j].Histo[k][i] *= ratio;
680 window[j].N[k] = (int)(ratio*window[j].N[k] + 0.5);
685 /* Simple linear interpolation between two given tabulated points */
686 double tabulated_pot(double dist, t_UmbrellaOptions *opt)
689 double pl, pu, dz, dp;
691 jl = floor((dist-opt->tabMin)/opt->tabDz);
693 if (jl < 0 || ju >= opt->tabNbins)
695 gmx_fatal(FARGS, "Distance %f out of bounds of tabulated potential (jl=%d, ju=%d).\n"
696 "Provide an extended table.", dist, jl, ju);
700 dz = dist-opt->tabX[jl];
701 dp = (pu-pl)*dz/opt->tabDz;
706 /* Don't worry, that routine does not mean we compute the PMF in limited precision.
707 After rapid convergence (using only substiantal contributions), we always switch to
709 void setup_acc_wham(double *profile, t_UmbrellaWindow * window, int nWindows,
710 t_UmbrellaOptions *opt)
712 int i, j, k, nGrptot = 0, nContrib = 0, nTot = 0;
713 double U, min = opt->min, dz = opt->dz, temp, ztot_half, distance, ztot, contrib1, contrib2;
714 gmx_bool bAnyContrib;
715 static int bFirst = 1;
716 static double wham_contrib_lim;
720 for (i = 0; i < nWindows; ++i)
722 nGrptot += window[i].nPull;
724 wham_contrib_lim = opt->Tolerance/nGrptot;
727 ztot = opt->max-opt->min;
730 for (i = 0; i < nWindows; ++i)
732 if (!window[i].bContrib)
734 snew(window[i].bContrib, window[i].nPull);
736 for (j = 0; j < window[i].nPull; ++j)
738 if (!window[i].bContrib[j])
740 snew(window[i].bContrib[j], opt->bins);
743 for (k = 0; k < opt->bins; ++k)
745 temp = (1.0*k+0.5)*dz+min;
746 distance = temp - window[i].pos[j]; /* distance to umbrella center */
748 { /* in cyclic wham: */
749 if (distance > ztot_half) /* |distance| < ztot_half */
753 else if (distance < -ztot_half)
758 /* Note: there are two contributions to bin k in the wham equations:
759 i) N[j]*exp(- U/(8.314e-3*opt->Temperature) + window[i].z[j])
760 ii) exp(- U/(8.314e-3*opt->Temperature))
761 where U is the umbrella potential
762 If any of these number is larger wham_contrib_lim, I set contrib=TRUE
767 U = 0.5*window[i].k[j]*sqr(distance); /* harmonic potential assumed. */
771 U = tabulated_pot(distance, opt); /* Use tabulated potential */
774 contrib1 = profile[k]*exp(-U/(8.314e-3*opt->Temperature));
775 contrib2 = window[i].N[j]*exp(-U/(8.314e-3*opt->Temperature) + window[i].z[j]);
776 window[i].bContrib[j][k] = (contrib1 > wham_contrib_lim || contrib2 > wham_contrib_lim);
777 bAnyContrib = (bAnyContrib | window[i].bContrib[j][k]);
778 if (window[i].bContrib[j][k])
784 /* If this histo is far outside min and max all bContrib may be FALSE,
785 causing a floating point exception later on. To avoid that, switch
789 for (k = 0; k < opt->bins; ++k)
791 window[i].bContrib[j][k] = TRUE;
798 printf("Initialized rapid wham stuff (contrib tolerance %g)\n"
799 "Evaluating only %d of %d expressions.\n\n", wham_contrib_lim, nContrib, nTot);
804 printf("Updated rapid wham stuff. (evaluating only %d of %d contributions)\n",
811 void calc_profile(double *profile, t_UmbrellaWindow * window, int nWindows,
812 t_UmbrellaOptions *opt, gmx_bool bExact)
815 double num, ztot_half, ztot, distance, min = opt->min, dz = opt->dz;
816 double denom, U = 0, temp = 0, invg;
818 ztot = opt->max-opt->min;
821 for (i = 0; i < opt->bins; ++i)
824 for (j = 0; j < nWindows; ++j)
826 for (k = 0; k < window[j].nPull; ++k)
828 invg = 1.0/window[j].g[k] * window[j].bsWeight[k];
829 temp = (1.0*i+0.5)*dz+min;
830 num += invg*window[j].Histo[k][i];
832 if (!(bExact || window[j].bContrib[k][i]))
836 distance = temp - window[j].pos[k]; /* distance to umbrella center */
838 { /* in cyclic wham: */
839 if (distance > ztot_half) /* |distance| < ztot_half */
843 else if (distance < -ztot_half)
851 U = 0.5*window[j].k[k]*sqr(distance); /* harmonic potential assumed. */
855 U = tabulated_pot(distance, opt); /* Use tabulated potential */
857 denom += invg*window[j].N[k]*exp(-U/(8.314e-3*opt->Temperature) + window[j].z[k]);
860 profile[i] = num/denom;
865 double calc_z(double * profile, t_UmbrellaWindow * window, int nWindows,
866 t_UmbrellaOptions *opt, gmx_bool bExact)
868 int i, j, k, binMax = -1;
869 double U = 0, min = opt->min, dz = opt->dz, temp, ztot_half, distance, ztot, totalMax;
870 double MAX = -1e20, total = 0;
872 ztot = opt->max-opt->min;
875 for (i = 0; i < nWindows; ++i)
877 for (j = 0; j < window[i].nPull; ++j)
880 for (k = 0; k < window[i].nBin; ++k)
882 if (!(bExact || window[i].bContrib[j][k]))
886 temp = (1.0*k+0.5)*dz+min;
887 distance = temp - window[i].pos[j]; /* distance to umbrella center */
889 { /* in cyclic wham: */
890 if (distance > ztot_half) /* |distance| < ztot_half */
894 else if (distance < -ztot_half)
902 U = 0.5*window[i].k[j]*sqr(distance); /* harmonic potential assumed. */
906 U = tabulated_pot(distance, opt); /* Use tabulated potential */
909 total += profile[k]*exp(-U/(8.314e-3*opt->Temperature));
911 /* Avoid floating point exception if window is far outside min and max */
920 temp = fabs(total - window[i].z[j]);
927 window[i].z[j] = total;
933 void symmetrizeProfile(double* profile, t_UmbrellaOptions *opt)
936 double *prof2, bins = opt->bins, min = opt->min, max = opt->max, dz = opt->dz, zsym, deltaz, profsym;
939 if (min > 0. || max < 0.)
941 gmx_fatal(FARGS, "Cannot symmetrize profile around z=0 with min=%f and max=%f\n",
947 for (i = 0; i < bins; i++)
951 /* bin left of zsym */
952 j = floor((zsym-min)/dz-0.5);
953 if (j >= 0 && (j+1) < bins)
955 /* interpolate profile linearly between bins j and j+1 */
958 profsym = profile[j] + (profile[j+1]-profile[j])/dz*deltaz;
959 /* average between left and right */
960 prof2[i] = 0.5*(profsym+profile[i]);
964 prof2[i] = profile[i];
968 memcpy(profile, prof2, bins*sizeof(double));
972 void prof_normalization_and_unit(double * profile, t_UmbrellaOptions *opt)
975 double unit_factor = 1., R_MolarGasConst, diff;
977 R_MolarGasConst = 8.314472e-3; /* in kJ/(mol*K) */
980 /* Not log? Nothing to do! */
986 /* Get profile in units of kT, kJ/mol, or kCal/mol */
987 if (opt->unit == en_kT)
991 else if (opt->unit == en_kJ)
993 unit_factor = R_MolarGasConst*opt->Temperature;
995 else if (opt->unit == en_kCal)
997 unit_factor = R_MolarGasConst*opt->Temperature/4.1868;
1001 gmx_fatal(FARGS, "Sorry, I don't know this energy unit.");
1004 for (i = 0; i < bins; i++)
1006 if (profile[i] > 0.0)
1008 profile[i] = -log(profile[i])*unit_factor;
1012 /* shift to zero at z=opt->zProf0 */
1013 if (!opt->bProf0Set)
1019 /* Get bin with shortest distance to opt->zProf0
1020 (-0.5 from bin position and +0.5 from rounding cancel) */
1021 imin = (int)((opt->zProf0-opt->min)/opt->dz);
1026 else if (imin >= bins)
1030 diff = profile[imin];
1034 for (i = 0; i < bins; i++)
1040 void getRandomIntArray(int nPull, int blockLength, int* randomArray, gmx_rng_t rng)
1042 int ipull, blockBase, nr, ipullRandom;
1044 if (blockLength == 0)
1046 blockLength = nPull;
1049 for (ipull = 0; ipull < nPull; ipull++)
1051 blockBase = (ipull/blockLength)*blockLength;
1053 { /* make sure nothing bad happens in the last block */
1054 nr = (int)(gmx_rng_uniform_real(rng)*blockLength);
1055 ipullRandom = blockBase + nr;
1057 while (ipullRandom >= nPull);
1058 if (ipullRandom < 0 || ipullRandom >= nPull)
1060 gmx_fatal(FARGS, "Ups, random iWin = %d, nPull = %d, nr = %d, "
1061 "blockLength = %d, blockBase = %d\n",
1062 ipullRandom, nPull, nr, blockLength, blockBase);
1064 randomArray[ipull] = ipullRandom;
1066 /*for (ipull=0; ipull<nPull; ipull++)
1067 printf("%d ",randomArray[ipull]); printf("\n"); */
1070 void copy_pullgrp_to_synthwindow(t_UmbrellaWindow *synthWindow,
1071 t_UmbrellaWindow *thisWindow, int pullid)
1073 synthWindow->N [0] = thisWindow->N [pullid];
1074 synthWindow->Histo [0] = thisWindow->Histo [pullid];
1075 synthWindow->pos [0] = thisWindow->pos [pullid];
1076 synthWindow->z [0] = thisWindow->z [pullid];
1077 synthWindow->k [0] = thisWindow->k [pullid];
1078 synthWindow->bContrib[0] = thisWindow->bContrib [pullid];
1079 synthWindow->g [0] = thisWindow->g [pullid];
1080 synthWindow->bsWeight[0] = thisWindow->bsWeight [pullid];
1083 /* Calculate cumulative distribution function of of all histograms. They
1084 allow to create random number sequences
1085 which are distributed according to the histograms. Required to generate
1086 the "synthetic" histograms for the Bootstrap method */
1087 void calc_cumulatives(t_UmbrellaWindow *window, int nWindows,
1088 t_UmbrellaOptions *opt, const char *fnhist)
1092 char *fn = 0, *buf = 0;
1095 if (opt->bs_verbose)
1097 snew(fn, strlen(fnhist)+10);
1098 snew(buf, strlen(fnhist)+10);
1099 sprintf(fn, "%s_cumul.xvg", strncpy(buf, fnhist, strlen(fnhist)-4));
1100 fp = xvgropen(fn, "CDFs of umbrella windows", "z", "CDF", opt->oenv);
1104 for (i = 0; i < nWindows; i++)
1106 snew(window[i].cum, window[i].nPull);
1107 for (j = 0; j < window[i].nPull; j++)
1109 snew(window[i].cum[j], nbin+1);
1110 window[i].cum[j][0] = 0.;
1111 for (k = 1; k <= nbin; k++)
1113 window[i].cum[j][k] = window[i].cum[j][k-1]+window[i].Histo[j][k-1];
1116 /* normalize CDFs. Ensure cum[nbin]==1 */
1117 last = window[i].cum[j][nbin];
1118 for (k = 0; k <= nbin; k++)
1120 window[i].cum[j][k] /= last;
1125 printf("Cumulative distriubtion functions of all histograms created.\n");
1126 if (opt->bs_verbose)
1128 for (k = 0; k <= nbin; k++)
1130 fprintf(fp, "%g\t", opt->min+k*opt->dz);
1131 for (i = 0; i < nWindows; i++)
1133 for (j = 0; j < window[i].nPull; j++)
1135 fprintf(fp, "%g\t", window[i].cum[j][k]);
1140 printf("Wrote cumulative distribution functions to %s\n", fn);
1148 /* Return j such that xx[j] <= x < xx[j+1] */
1149 void searchCumulative(double xx[], int n, double x, int *j)
1171 else if (x == xx[n-1])
1181 void create_synthetic_histo(t_UmbrellaWindow *synthWindow, t_UmbrellaWindow *thisWindow,
1182 int pullid, t_UmbrellaOptions *opt)
1184 int N, i, nbins, r_index, ibin;
1185 double r, tausteps = 0.0, a, ap, dt, x, invsqrt2, g, y, sig = 0., z, mu = 0.;
1188 N = thisWindow->N[pullid];
1189 dt = thisWindow->dt;
1190 nbins = thisWindow->nBin;
1192 /* tau = autocorrelation time */
1193 if (opt->tauBootStrap > 0.0)
1195 tausteps = opt->tauBootStrap/dt;
1197 else if (opt->bTauIntGiven || opt->bCalcTauInt)
1199 /* calc tausteps from g=1+2tausteps */
1200 g = thisWindow->g[pullid];
1206 "When generating hypothetical trajctories from given umbrella histograms,\n"
1207 "autocorrelation times (ACTs) are required. Otherwise the statistical error\n"
1208 "cannot be predicted. You have 3 options:\n"
1209 "1) Make g_wham estimate the ACTs (options -ac and -acsig).\n"
1210 "2) Calculate the ACTs by yourself (e.g. with g_analyze) and provide them\n");
1212 " with option -iiact for all umbrella windows.\n"
1213 "3) If all ACTs are identical and know, you can define them with -bs-tau.\n"
1214 " Use option (3) only if you are sure what you're doing, you may severely\n"
1215 " underestimate the error if a too small ACT is given.\n");
1216 gmx_fatal(FARGS, errstr);
1219 synthWindow->N [0] = N;
1220 synthWindow->pos [0] = thisWindow->pos[pullid];
1221 synthWindow->z [0] = thisWindow->z[pullid];
1222 synthWindow->k [0] = thisWindow->k[pullid];
1223 synthWindow->bContrib[0] = thisWindow->bContrib[pullid];
1224 synthWindow->g [0] = thisWindow->g [pullid];
1225 synthWindow->bsWeight[0] = thisWindow->bsWeight[pullid];
1227 for (i = 0; i < nbins; i++)
1229 synthWindow->Histo[0][i] = 0.;
1232 if (opt->bsMethod == bsMethod_trajGauss)
1234 sig = thisWindow->sigma [pullid];
1235 mu = thisWindow->aver [pullid];
1238 /* Genrate autocorrelated Gaussian random variable with autocorrelation time tau
1240 If x and y are random numbers from N(0,1) (Gaussian with average 0 and sigma=1),
1242 z = a*x + sqrt(1-a^2)*y
1243 is also from N(0,1), and cov(z,x) = a. Thus, by gerenating a sequence
1244 x' = a*x + sqrt(1-a^2)*y, the sequnce x(t) is from N(0,1) and has an autocorrelation
1246 C(t) = exp(-t/tau) with tau=-1/ln(a)
1248 Then, use error function to turn the Gaussian random variable into a uniformly
1249 distributed one in [0,1]. Eventually, use cumulative distribution function of
1250 histogram to get random variables distributed according to histogram.
1251 Note: The ACT of the flat distribution and of the generated histogram is not
1252 100% exactly tau, but near tau (my test was 3.8 instead of 4).
1254 a = exp(-1.0/tausteps);
1256 invsqrt2 = 1./sqrt(2.0);
1258 /* init random sequence */
1259 x = gmx_rng_gaussian_table(opt->rng);
1261 if (opt->bsMethod == bsMethod_traj)
1263 /* bootstrap points from the umbrella histograms */
1264 for (i = 0; i < N; i++)
1266 y = gmx_rng_gaussian_table(opt->rng);
1268 /* get flat distribution in [0,1] using cumulative distribution function of Gauusian
1269 Note: CDF(Gaussian) = 0.5*{1+erf[x/sqrt(2)]}
1271 r = 0.5*(1+gmx_erf(x*invsqrt2));
1272 searchCumulative(thisWindow->cum[pullid], nbins+1, r, &r_index);
1273 synthWindow->Histo[0][r_index] += 1.;
1276 else if (opt->bsMethod == bsMethod_trajGauss)
1278 /* bootstrap points from a Gaussian with the same average and sigma
1279 as the respective umbrella histogram. The idea was, that -given
1280 limited sampling- the bootstrapped histograms are otherwise biased
1281 from the limited sampling of the US histos. However, bootstrapping from
1282 the Gaussian seems to yield a similar estimate. */
1286 y = gmx_rng_gaussian_table(opt->rng);
1289 ibin = floor((z-opt->min)/opt->dz);
1294 while ( (ibin += nbins) < 0)
1299 else if (ibin >= nbins)
1301 while ( (ibin -= nbins) >= nbins)
1308 if (ibin >= 0 && ibin < nbins)
1310 synthWindow->Histo[0][ibin] += 1.;
1317 gmx_fatal(FARGS, "Unknown bsMethod (id %d). That should not happen.\n", opt->bsMethod);
1322 void print_histograms(const char *fnhist, t_UmbrellaWindow * window, int nWindows,
1323 int bs_index, t_UmbrellaOptions *opt)
1325 char *fn = 0, *buf = 0, title[256];
1331 fn = strdup(fnhist);
1332 strcpy(title, "Umbrella histograms");
1336 snew(fn, strlen(fnhist)+10);
1337 snew(buf, strlen(fnhist)+1);
1338 sprintf(fn, "%s_bs%d.xvg", strncpy(buf, fnhist, strlen(fnhist)-4), bs_index);
1339 sprintf(title, "Umbrella histograms. Bootstrap #%d", bs_index);
1342 fp = xvgropen(fn, title, "z", "count", opt->oenv);
1345 /* Write histograms */
1346 for (l = 0; l < bins; ++l)
1348 fprintf(fp, "%e\t", (double)(l+0.5)*opt->dz+opt->min);
1349 for (i = 0; i < nWindows; ++i)
1351 for (j = 0; j < window[i].nPull; ++j)
1353 fprintf(fp, "%e\t", window[i].Histo[j][l]);
1360 printf("Wrote %s\n", fn);
1368 int func_wham_is_larger(const void *a, const void *b)
1388 void setRandomBsWeights(t_UmbrellaWindow *synthwin, int nAllPull, t_UmbrellaOptions *opt)
1395 /* generate ordered random numbers between 0 and nAllPull */
1396 for (i = 0; i < nAllPull-1; i++)
1398 r[i] = gmx_rng_uniform_real(opt->rng) * nAllPull;
1400 qsort((void *)r, nAllPull-1, sizeof(double), &func_wham_is_larger);
1401 r[nAllPull-1] = 1.0*nAllPull;
1403 synthwin[0].bsWeight[0] = r[0];
1404 for (i = 1; i < nAllPull; i++)
1406 synthwin[i].bsWeight[0] = r[i]-r[i-1];
1409 /* avoid to have zero weight by adding a tiny value */
1410 for (i = 0; i < nAllPull; i++)
1412 if (synthwin[i].bsWeight[0] < 1e-5)
1414 synthwin[i].bsWeight[0] = 1e-5;
1421 void do_bootstrapping(const char *fnres, const char* fnprof, const char *fnhist,
1422 char* ylabel, double *profile,
1423 t_UmbrellaWindow * window, int nWindows, t_UmbrellaOptions *opt)
1425 t_UmbrellaWindow * synthWindow;
1426 double *bsProfile, *bsProfiles_av, *bsProfiles_av2, maxchange = 1e20, tmp, stddev;
1427 int i, j, *randomArray = 0, winid, pullid, ib;
1428 int iAllPull, nAllPull, *allPull_winId, *allPull_pullId;
1430 gmx_bool bExact = FALSE;
1432 /* init random generator */
1433 if (opt->bsSeed == -1)
1435 opt->rng = gmx_rng_init(gmx_rng_make_seed());
1439 opt->rng = gmx_rng_init(opt->bsSeed);
1442 snew(bsProfile, opt->bins);
1443 snew(bsProfiles_av, opt->bins);
1444 snew(bsProfiles_av2, opt->bins);
1446 /* Create array of all pull groups. Note that different windows
1447 may have different nr of pull groups
1448 First: Get total nr of pull groups */
1450 for (i = 0; i < nWindows; i++)
1452 nAllPull += window[i].nPull;
1454 snew(allPull_winId, nAllPull);
1455 snew(allPull_pullId, nAllPull);
1457 /* Setup one array of all pull groups */
1458 for (i = 0; i < nWindows; i++)
1460 for (j = 0; j < window[i].nPull; j++)
1462 allPull_winId[iAllPull] = i;
1463 allPull_pullId[iAllPull] = j;
1468 /* setup stuff for synthetic windows */
1469 snew(synthWindow, nAllPull);
1470 for (i = 0; i < nAllPull; i++)
1472 synthWindow[i].nPull = 1;
1473 synthWindow[i].nBin = opt->bins;
1474 snew(synthWindow[i].Histo, 1);
1475 if (opt->bsMethod == bsMethod_traj || opt->bsMethod == bsMethod_trajGauss)
1477 snew(synthWindow[i].Histo[0], opt->bins);
1479 snew(synthWindow[i].N, 1);
1480 snew(synthWindow[i].pos, 1);
1481 snew(synthWindow[i].z, 1);
1482 snew(synthWindow[i].k, 1);
1483 snew(synthWindow[i].bContrib, 1);
1484 snew(synthWindow[i].g, 1);
1485 snew(synthWindow[i].bsWeight, 1);
1488 switch (opt->bsMethod)
1491 snew(randomArray, nAllPull);
1492 printf("\n\nWhen computing statistical errors by bootstrapping entire histograms:\n");
1493 please_cite(stdout, "Hub2006");
1495 case bsMethod_BayesianHist:
1496 /* just copy all histogams into synthWindow array */
1497 for (i = 0; i < nAllPull; i++)
1499 winid = allPull_winId [i];
1500 pullid = allPull_pullId[i];
1501 copy_pullgrp_to_synthwindow(synthWindow+i, window+winid, pullid);
1505 case bsMethod_trajGauss:
1506 calc_cumulatives(window, nWindows, opt, fnhist);
1509 gmx_fatal(FARGS, "Unknown bootstrap method. That should not have happened.\n");
1512 /* do bootstrapping */
1513 fp = xvgropen(fnprof, "Boot strap profiles", "z", ylabel, opt->oenv);
1514 for (ib = 0; ib < opt->nBootStrap; ib++)
1516 printf(" *******************************************\n"
1517 " ******** Start bootstrap nr %d ************\n"
1518 " *******************************************\n", ib+1);
1520 switch (opt->bsMethod)
1523 /* bootstrap complete histograms from given histograms */
1524 getRandomIntArray(nAllPull, opt->histBootStrapBlockLength, randomArray, opt->rng);
1525 for (i = 0; i < nAllPull; i++)
1527 winid = allPull_winId [randomArray[i]];
1528 pullid = allPull_pullId[randomArray[i]];
1529 copy_pullgrp_to_synthwindow(synthWindow+i, window+winid, pullid);
1532 case bsMethod_BayesianHist:
1533 /* keep histos, but assign random weights ("Bayesian bootstrap") */
1534 setRandomBsWeights(synthWindow, nAllPull, opt);
1537 case bsMethod_trajGauss:
1538 /* create new histos from given histos, that is generate new hypothetical
1540 for (i = 0; i < nAllPull; i++)
1542 winid = allPull_winId[i];
1543 pullid = allPull_pullId[i];
1544 create_synthetic_histo(synthWindow+i, window+winid, pullid, opt);
1549 /* write histos in case of verbose output */
1550 if (opt->bs_verbose)
1552 print_histograms(fnhist, synthWindow, nAllPull, ib, opt);
1559 memcpy(bsProfile, profile, opt->bins*sizeof(double)); /* use profile as guess */
1562 if ( (i%opt->stepUpdateContrib) == 0)
1564 setup_acc_wham(bsProfile, synthWindow, nAllPull, opt);
1566 if (maxchange < opt->Tolerance)
1570 if (((i%opt->stepchange) == 0 || i == 1) && !i == 0)
1572 printf("\t%4d) Maximum change %e\n", i, maxchange);
1574 calc_profile(bsProfile, synthWindow, nAllPull, opt, bExact);
1577 while ( (maxchange = calc_z(bsProfile, synthWindow, nAllPull, opt, bExact)) > opt->Tolerance || !bExact);
1578 printf("\tConverged in %d iterations. Final maximum change %g\n", i, maxchange);
1582 prof_normalization_and_unit(bsProfile, opt);
1585 /* symmetrize profile around z=0 */
1588 symmetrizeProfile(bsProfile, opt);
1591 /* save stuff to get average and stddev */
1592 for (i = 0; i < opt->bins; i++)
1595 bsProfiles_av[i] += tmp;
1596 bsProfiles_av2[i] += tmp*tmp;
1597 fprintf(fp, "%e\t%e\n", (i+0.5)*opt->dz+opt->min, tmp);
1603 /* write average and stddev */
1604 fp = xvgropen(fnres, "Average and stddev from bootstrapping", "z", ylabel, opt->oenv);
1605 fprintf(fp, "@TYPE xydy\n");
1606 for (i = 0; i < opt->bins; i++)
1608 bsProfiles_av [i] /= opt->nBootStrap;
1609 bsProfiles_av2[i] /= opt->nBootStrap;
1610 tmp = bsProfiles_av2[i]-sqr(bsProfiles_av[i]);
1611 stddev = (tmp >= 0.) ? sqrt(tmp) : 0.; /* Catch rouding errors */
1612 fprintf(fp, "%e\t%e\t%e\n", (i+0.5)*opt->dz+opt->min, bsProfiles_av [i], stddev);
1615 printf("Wrote boot strap result to %s\n", fnres);
1618 int whaminFileType(char *fn)
1622 if (strcmp(fn+len-3, "tpr") == 0)
1626 else if (strcmp(fn+len-3, "xvg") == 0 || strcmp(fn+len-6, "xvg.gz") == 0)
1628 return whamin_pullxf;
1630 else if (strcmp(fn+len-3, "pdo") == 0 || strcmp(fn+len-6, "pdo.gz") == 0)
1636 gmx_fatal(FARGS, "Unknown file type of %s. Should be tpr, xvg, or pdo.\n", fn);
1638 return whamin_unknown;
1641 void read_wham_in(const char *fn, char ***filenamesRet, int *nfilesRet,
1642 t_UmbrellaOptions *opt)
1644 char **filename = 0, tmp[STRLEN];
1645 int nread, sizenow, i, block = 1;
1648 fp = ffopen(fn, "r");
1651 while (fscanf(fp, "%s", tmp) != EOF)
1653 if (strlen(tmp) >= WHAM_MAXFILELEN)
1655 gmx_fatal(FARGS, "Filename too long. Only %d characters allowed\n", WHAM_MAXFILELEN);
1657 if (nread >= sizenow)
1660 srenew(filename, sizenow);
1661 for (i = sizenow-block; i < sizenow; i++)
1663 snew(filename[i], WHAM_MAXFILELEN);
1666 strcpy(filename[nread], tmp);
1669 printf("Found file %s in %s\n", filename[nread], fn);
1673 *filenamesRet = filename;
1678 FILE *open_pdo_pipe(const char *fn, t_UmbrellaOptions *opt, gmx_bool *bPipeOpen)
1680 char Buffer[1024], gunzip[1024], *Path = 0;
1682 static gmx_bool bFirst = 1;
1684 /* gzipped pdo file? */
1685 if ((strcmp(fn+strlen(fn)-3, ".gz") == 0))
1687 /* search gunzip executable */
1688 if (!(Path = getenv("GMX_PATH_GZIP")))
1690 if (gmx_fexist("/bin/gunzip"))
1692 sprintf(gunzip, "%s", "/bin/gunzip");
1694 else if (gmx_fexist("/usr/bin/gunzip"))
1696 sprintf(gunzip, "%s", "/usr/bin/gunzip");
1700 gmx_fatal(FARGS, "Cannot find executable gunzip in /bin or /usr/bin.\n"
1701 "You may want to define the path to gunzip "
1702 "with the environment variable GMX_PATH_GZIP.", gunzip);
1707 sprintf(gunzip, "%s/gunzip", Path);
1708 if (!gmx_fexist(gunzip))
1710 gmx_fatal(FARGS, "Cannot find executable %s. Please define the path to gunzip"
1711 " in the environmental varialbe GMX_PATH_GZIP.", gunzip);
1716 printf("Using gunzig executable %s\n", gunzip);
1719 if (!gmx_fexist(fn))
1721 gmx_fatal(FARGS, "File %s does not exist.\n", fn);
1723 sprintf(Buffer, "%s -c < %s", gunzip, fn);
1726 printf("Executing command '%s'\n", Buffer);
1729 if ((pipe = popen(Buffer, "r")) == NULL)
1731 gmx_fatal(FARGS, "Unable to open pipe to `%s'\n", Buffer);
1734 gmx_fatal(FARGS, "Cannot open a compressed file on platform without pipe support");
1740 pipe = ffopen(fn, "r");
1747 void pdo_close_file(FILE *fp)
1756 /* Reading pdo files */
1757 void read_pdo_files(char **fn, int nfiles, t_UmbrellaHeader* header,
1758 t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
1761 real mintmp, maxtmp, done = 0.;
1764 /* char Buffer0[1000]; */
1768 gmx_fatal(FARGS, "No files found. Hick.");
1771 /* if min and max are not given, get min and max from the input files */
1774 printf("Automatic determination of boundaries from %d pdo files...\n", nfiles);
1777 for (i = 0; i < nfiles; ++i)
1779 file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
1780 /*fgets(Buffer0,999,file);
1781 fprintf(stderr,"First line '%s'\n",Buffer0); */
1782 done = 100.0*(i+1)/nfiles;
1783 printf("\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
1788 read_pdo_header(file, header, opt);
1789 /* here only determine min and max of this window */
1790 read_pdo_data(file, header, i, NULL, opt, TRUE, &mintmp, &maxtmp);
1791 if (maxtmp > opt->max)
1795 if (mintmp < opt->min)
1801 pdo_close_file(file);
1809 printf("\nDetermined boundaries to %f and %f\n\n", opt->min, opt->max);
1810 if (opt->bBoundsOnly)
1812 printf("Found option -boundsonly, now exiting.\n");
1816 /* store stepsize in profile */
1817 opt->dz = (opt->max-opt->min)/opt->bins;
1819 /* Having min and max, we read in all files */
1820 /* Loop over all files */
1821 for (i = 0; i < nfiles; ++i)
1823 done = 100.0*(i+1)/nfiles;
1824 printf("\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
1829 file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
1830 read_pdo_header(file, header, opt);
1831 /* load data into window */
1832 read_pdo_data(file, header, i, window, opt, FALSE, NULL, NULL);
1833 if ((window+i)->Ntot[0] == 0.0)
1835 fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fn[i]);
1839 pdo_close_file(file);
1847 for (i = 0; i < nfiles; ++i)
1854 #define int2YN(a) (((a) == 0) ? ("N") : ("Y"))
1856 void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions *opt)
1861 static int first = 1;
1863 /* printf("Reading %s \n",fn); */
1864 read_tpx_state(fn, &ir, &state, NULL, NULL);
1866 if (ir.ePull != epullUMBRELLA)
1868 gmx_fatal(FARGS, "This is not a tpr of an umbrella simulation. Found pull type \"%s\" "
1869 " (ir.ePull = %d)\n", epull_names[ir.ePull], ir.ePull);
1872 /* nr of pull groups */
1873 ngrp = ir.pull->ngrp;
1876 gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull groups\n", ngrp);
1879 header->npullgrps = ir.pull->ngrp;
1880 header->pull_geometry = ir.pull->eGeom;
1881 copy_ivec(ir.pull->dim, header->pull_dim);
1882 header->pull_ndim = header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2];
1883 if (header->pull_geometry == epullgPOS && header->pull_ndim > 1)
1885 gmx_fatal(FARGS, "Found pull geometry 'position' and more than 1 pull dimension (%d).\n"
1886 "Hence, the pull potential does not correspond to a one-dimensional umbrella potential.\n"
1887 "If you have some special umbrella setup you may want to write your own pdo files\n"
1888 "and feed them into g_wham. Check g_wham -h !\n", header->pull_ndim);
1890 snew(header->k, ngrp);
1891 snew(header->init_dist, ngrp);
1892 snew(header->umbInitDist, ngrp);
1894 /* only z-direction with epullgCYL? */
1895 if (header->pull_geometry == epullgCYL)
1897 if (header->pull_dim[XX] || header->pull_dim[YY] || (!header->pull_dim[ZZ]))
1899 gmx_fatal(FARGS, "With pull geometry 'cylinder', expected pulling in Z direction only.\n"
1900 "However, found dimensions [%s %s %s]\n",
1901 int2YN(header->pull_dim[XX]), int2YN(header->pull_dim[YY]),
1902 int2YN(header->pull_dim[ZZ]));
1906 for (i = 0; i < ngrp; i++)
1908 header->k[i] = ir.pull->grp[i+1].k;
1909 if (header->k[i] == 0.0)
1911 gmx_fatal(FARGS, "Pull group %d has force constant of of 0.0 in %s.\n"
1912 "That doesn't seem to be an Umbrella tpr.\n",
1915 copy_rvec(ir.pull->grp[i+1].init, header->init_dist[i]);
1917 /* initial distance to reference */
1918 switch (header->pull_geometry)
1921 for (d = 0; d < DIM; d++)
1923 if (header->pull_dim[d])
1925 header->umbInitDist[i] = header->init_dist[i][d];
1930 /* umbrella distance stored in init_dist[i][0] for geometry cylinder (not in ...[i][ZZ]) */
1934 header->umbInitDist[i] = header->init_dist[i][0];
1937 gmx_fatal(FARGS, "Pull geometry %s not supported\n", epullg_names[header->pull_geometry]);
1941 if (opt->verbose || first)
1943 printf("File %s, %d groups, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n",
1944 fn, header->npullgrps, epullg_names[header->pull_geometry],
1945 int2YN(header->pull_dim[0]), int2YN(header->pull_dim[1]), int2YN(header->pull_dim[2]),
1947 for (i = 0; i < ngrp; i++)
1949 printf("\tgrp %d) k = %-5g position = %g\n", i, header->k[i], header->umbInitDist[i]);
1952 if (!opt->verbose && first)
1954 printf("\tUse option -v to see this output for all input tpr files\n");
1961 double dist_ndim(double **dx, int ndim, int line)
1965 for (i = 0; i < ndim; i++)
1967 r2 += sqr(dx[i][line]);
1972 void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header,
1973 t_UmbrellaWindow * window,
1974 t_UmbrellaOptions *opt,
1975 gmx_bool bGetMinMax, real *mintmp, real *maxtmp)
1977 double **y = 0, pos = 0., t, force, time0 = 0., dt;
1978 int ny, nt, bins, ibin, i, g, dstep = 1, nColPerGrp, nColRefOnce, nColRefEachGrp, nColExpect, ntot;
1979 real min, max, minfound = 1e20, maxfound = -1e20;
1980 gmx_bool dt_ok, timeok, bHaveForce;
1981 const char *quantity;
1982 const int blocklen = 4096;
1986 in force output pullf.xvg:
1987 No reference, one column per pull group
1988 in position output pullx.xvg (not cylinder)
1989 ndim reference, ndim columns per pull group
1990 in position output pullx.xvg (in geometry cylinder):
1991 ndim*2 columns per pull group (ndim for ref, ndim for group)
1994 nColPerGrp = opt->bPullx ? header->pull_ndim : 1;
1995 quantity = opt->bPullx ? "position" : "force";
1999 if (header->pull_geometry == epullgCYL)
2001 /* Geometry cylinder -> reference group before each pull group */
2002 nColRefEachGrp = header->pull_ndim;
2007 /* Geometry NOT cylinder -> reference group only once after time column */
2009 nColRefOnce = header->pull_ndim;
2012 else /* read forces, no reference groups */
2018 nColExpect = 1 + nColRefOnce + header->npullgrps*(nColRefEachGrp+nColPerGrp);
2019 bHaveForce = opt->bPullf;
2021 /* With geometry "distance" or "distance_periodic", only force reading is supported so far.
2022 That avoids the somewhat tedious extraction of the right columns from the pullx files
2023 to compute the distances projection on the vector. Sorry for the laziness. */
2024 if ( (header->pull_geometry == epullgDIR || header->pull_geometry == epullgDIRPBC)
2027 gmx_fatal(FARGS, "With pull geometries \"direction\" and \"direction_periodic\", only pull force "
2028 "reading \n(option -if) is supported at present, "
2029 "not pull position reading (options -ix).\nMake sure mdrun writes the pull "
2030 "forces (pullf.xvg files)\nand provide them to g_wham with option -if.",
2031 epullg_names[header->pull_geometry]);
2034 nt = read_xvg(fn, &y, &ny);
2036 /* Check consistency */
2039 gmx_fatal(FARGS, "Empty pull %s file %s\n", quantity, fn);
2041 if (ny != nColExpect)
2043 gmx_fatal(FARGS, "Found %d pull groups in %s,\n but %d data columns in %s (expected %d)\n"
2044 "\nMaybe you confused options -ix and -if ?\n",
2045 header->npullgrps, fntpr, ny-1, fn, nColExpect-1);
2050 printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerGrp, quantity, fn);
2060 window->dt = y[0][1]-y[0][0];
2062 else if (opt->nBootStrap && opt->tauBootStrap != 0.0)
2064 fprintf(stderr, "\n *** WARNING, Could not determine time step in %s\n", fn);
2067 /* Need to alocate memory and set up structure */
2068 window->nPull = header->npullgrps;
2069 window->nBin = bins;
2071 snew(window->Histo, window->nPull);
2072 snew(window->z, window->nPull);
2073 snew(window->k, window->nPull);
2074 snew(window->pos, window->nPull);
2075 snew(window->N, window->nPull);
2076 snew(window->Ntot, window->nPull);
2077 snew(window->g, window->nPull);
2078 snew(window->bsWeight, window->nPull);
2079 window->bContrib = 0;
2081 if (opt->bCalcTauInt)
2083 snew(window->ztime, window->nPull);
2087 window->ztime = NULL;
2089 snew(lennow, window->nPull);
2091 for (g = 0; g < window->nPull; ++g)
2094 window->bsWeight[g] = 1.;
2095 snew(window->Histo[g], bins);
2096 window->k[g] = header->k[g];
2098 window->Ntot[g] = 0;
2100 window->pos[g] = header->umbInitDist[g];
2101 if (opt->bCalcTauInt)
2103 window->ztime[g] = NULL;
2109 { /* only determine min and max */
2112 min = max = bins = 0; /* Get rid of warnings */
2115 for (i = 0; i < nt; i++)
2117 /* Do you want that time frame? */
2118 t = 1.0/1000*( (int) ((y[0][i]*1000) + 0.5)); /* round time to fs */
2120 /* get time step of pdo file and get dstep from opt->dt */
2130 dstep = (int)(opt->dt/dt+0.5);
2138 window->dt = dt*dstep;
2142 dt_ok = (i%dstep == 0);
2143 timeok = (dt_ok && t >= opt->tmin && t <= opt->tmax);
2145 printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
2146 t,opt->tmin, opt->tmax, dt_ok,timeok); */
2150 for (g = 0; g < header->npullgrps; ++g)
2154 /* y has 1 time column y[0] and one column per force y[1],...,y[nGrps] */
2156 pos = -force/header->k[g] + header->umbInitDist[g];
2160 switch (header->pull_geometry)
2163 /* y has 1 time column y[0] and nColPerGrps columns per pull group;
2164 Distance to reference: */
2165 /* pos=dist_ndim(y+1+nColRef+g*nColPerGrp,header->pull_ndim,i); gmx 4.0 */
2166 pos = dist_ndim(y + 1 + nColRefOnce + g*nColPerGrp, header->pull_ndim, i);
2170 Time ref[ndim] group1[ndim] group2[ndim] ... */
2173 Time ref1[ndim] group1[ndim] ref2[ndim] group2[ndim] ... */
2175 /* * with geometry==position, we have the reference once (nColRefOnce==ndim), but
2176 no extra reference group columns before each group (nColRefEachGrp==0)
2178 * with geometry==cylinder, we have no initial ref group column (nColRefOnce==0),
2179 but ndim ref group colums before every group (nColRefEachGrp==ndim)
2180 Distance to reference: */
2181 pos = y[1 + nColRefOnce + nColRefEachGrp + g*(nColRefEachGrp+nColPerGrp)][i];
2184 gmx_fatal(FARGS, "Bad error, this error should have been catched before. Ups.\n");
2188 /* printf("grp %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
2202 if (opt->bCalcTauInt && !bGetMinMax)
2204 /* save time series for autocorrelation analysis */
2205 ntot = window->Ntot[g];
2206 /* printf("i %d, ntot %d, lennow[g] = %d\n",i,ntot,lennow[g]); */
2207 if (ntot >= lennow[g])
2209 lennow[g] += blocklen;
2210 srenew(window->ztime[g], lennow[g]);
2212 window->ztime[g][ntot] = pos;
2215 ibin = (int) floor((pos-min)/(max-min)*bins);
2220 while ( (ibin += bins) < 0)
2225 else if (ibin >= bins)
2227 while ( (ibin -= bins) >= bins)
2233 if (ibin >= 0 && ibin < bins)
2235 window->Histo[g][ibin] += 1.;
2242 else if (t > opt->tmax)
2246 printf("time %f larger than tmax %f, stop reading this pullx/pullf file\n", t, opt->tmax);
2258 for (i = 0; i < ny; i++)
2264 void read_tpr_pullxf_files(char **fnTprs, char **fnPull, int nfiles,
2265 t_UmbrellaHeader* header,
2266 t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
2269 real mintmp, maxtmp;
2271 printf("Reading %d tpr and pullf files\n", nfiles/2);
2273 /* min and max not given? */
2276 printf("Automatic determination of boundaries...\n");
2279 for (i = 0; i < nfiles; i++)
2281 if (whaminFileType(fnTprs[i]) != whamin_tpr)
2283 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a tpr file\n", i);
2285 read_tpr_header(fnTprs[i], header, opt);
2286 if (whaminFileType(fnPull[i]) != whamin_pullxf)
2288 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
2290 read_pull_xf(fnPull[i], fnTprs[i], header, NULL, opt, TRUE, &mintmp, &maxtmp);
2291 if (maxtmp > opt->max)
2295 if (mintmp < opt->min)
2300 printf("\nDetermined boundaries to %f and %f\n\n", opt->min, opt->max);
2301 if (opt->bBoundsOnly)
2303 printf("Found option -boundsonly, now exiting.\n");
2307 /* store stepsize in profile */
2308 opt->dz = (opt->max-opt->min)/opt->bins;
2310 for (i = 0; i < nfiles; i++)
2312 if (whaminFileType(fnTprs[i]) != whamin_tpr)
2314 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a tpr file\n", i);
2316 read_tpr_header(fnTprs[i], header, opt);
2317 if (whaminFileType(fnPull[i]) != whamin_pullxf)
2319 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
2321 read_pull_xf(fnPull[i], fnTprs[i], header, window+i, opt, FALSE, NULL, NULL);
2322 if (window[i].Ntot[0] == 0.0)
2324 fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fnPull[i]);
2328 for (i = 0; i < nfiles; i++)
2337 /* Note: Here we consider tau[int] := int_0^inf ACF(t) as the integrated autocorrelation time.
2338 The factor `g := 1 + 2*tau[int]` subsequently enters the uncertainty.
2340 void readIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt,
2343 int nlines, ny, i, ig;
2346 printf("Readging Integrated autocorrelation times from %s ...\n", fn);
2347 nlines = read_xvg(fn, &iact, &ny);
2348 if (nlines != nwins)
2350 gmx_fatal(FARGS, "Found %d lines with integrated autocorrelation times in %s.\nExpected %d",
2353 for (i = 0; i < nlines; i++)
2355 if (window[i].nPull != ny)
2357 gmx_fatal(FARGS, "You are providing autocorrelation times with option -iiact and the\n"
2358 "number of pull groups is different in different simulations. That is not\n"
2359 "supported yet. Sorry.\n");
2361 for (ig = 0; ig < window[i].nPull; ig++)
2363 /* compare Kumar et al, J Comp Chem 13, 1011-1021 (1992) */
2364 window[i].g[ig] = 1+2*iact[ig][i]/window[i].dt;
2366 if (iact[ig][i] <= 0.0)
2368 fprintf(stderr, "\nWARNING, IACT = %f (window %d, group %d)\n", iact[ig][i], i, ig);
2375 /* Smooth autocorreltion times along the reaction coordinate. This is useful
2376 if the ACT is subject to high uncertainty in case if limited sampling. Note
2377 that -in case of limited sampling- the ACT may be severely underestimated.
2378 Note: the g=1+2tau are overwritten.
2379 if opt->bAllowReduceIact==FALSE, the ACTs are never reduced, only increased
2382 void smoothIact(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt)
2385 double pos, dpos2, siglim, siglim2, gaufact, invtwosig2, w, weight, tausm;
2387 /* only evaluate within +- 3sigma of the Gausian */
2388 siglim = 3.0*opt->sigSmoothIact;
2389 siglim2 = dsqr(siglim);
2390 /* pre-factor of Gaussian */
2391 gaufact = 1.0/(sqrt(2*M_PI)*opt->sigSmoothIact);
2392 invtwosig2 = 0.5/dsqr(opt->sigSmoothIact);
2394 for (i = 0; i < nwins; i++)
2396 snew(window[i].tausmooth, window[i].nPull);
2397 for (ig = 0; ig < window[i].nPull; ig++)
2401 pos = window[i].pos[ig];
2402 for (j = 0; j < nwins; j++)
2404 for (jg = 0; jg < window[j].nPull; jg++)
2406 dpos2 = dsqr(window[j].pos[jg]-pos);
2407 if (dpos2 < siglim2)
2409 w = gaufact*exp(-dpos2*invtwosig2);
2411 tausm += w*window[j].tau[jg];
2412 /*printf("Weight %g dpos2=%g pos=%g gaufact=%g invtwosig2=%g\n",
2413 w,dpos2,pos,gaufact,invtwosig2); */
2418 if (opt->bAllowReduceIact || tausm > window[i].tau[ig])
2420 window[i].tausmooth[ig] = tausm;
2424 window[i].tausmooth[ig] = window[i].tau[ig];
2426 window[i].g[ig] = 1+2*tausm/window[i].dt;
2431 /* try to compute the autocorrelation time for each umbrealla window */
2432 #define WHAM_AC_ZERO_LIMIT 0.05
2433 void calcIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins,
2434 t_UmbrellaOptions *opt, const char *fn)
2436 int i, ig, ncorr, ntot, j, k, *count, restart;
2437 real *corr, c0, dt, timemax, tmp;
2438 real *ztime, av, tausteps;
2439 FILE *fp, *fpcorr = 0;
2443 fpcorr = xvgropen("hist_autocorr.xvg", "Autocorrelation functions of umbrella windows",
2444 "time [ps]", "autocorrelation function", opt->oenv);
2448 for (i = 0; i < nwins; i++)
2450 printf("\rEstimating integrated autocorreltion times ... [%2.0f%%] ...", 100.*(i+1)/nwins);
2452 ntot = window[i].Ntot[0];
2454 /* using half the maximum time as length of autocorrelation function */
2458 gmx_fatal(FARGS, "Tryig to estimtate autocorrelation time from only %d"
2459 " points. Provide more pull data!", ntot);
2462 /* snew(corrSq,ncorr); */
2466 snew(window[i].tau, window[i].nPull);
2467 restart = (int)(opt->acTrestart/dt+0.5);
2473 for (ig = 0; ig < window[i].nPull; ig++)
2475 if (ntot != window[i].Ntot[ig])
2477 gmx_fatal(FARGS, "Encountered different nr of frames in different pull groups.\n"
2478 "That should not happen. (%d and %d)\n", ntot, window[i].Ntot[ig]);
2480 ztime = window[i].ztime[ig];
2482 /* calc autocorrelation function C(t) = < [z(tau)-<z>]*[z(tau+t)-<z>]> */
2483 for (j = 0, av = 0; (j < ntot); j++)
2488 for (k = 0; (k < ncorr); k++)
2493 for (j = 0; (j < ntot); j += restart)
2495 for (k = 0; (k < ncorr) && (j+k < ntot); k++)
2497 tmp = (ztime[j]-av)*(ztime[j+k]-av);
2499 /* corrSq[k] += tmp*tmp; */
2503 /* divide by nr of frames for each time displacement */
2504 for (k = 0; (k < ncorr); k++)
2506 /* count probably = (ncorr-k+(restart-1))/restart; */
2507 corr[k] = corr[k]/count[k];
2508 /* variance of autocorrelation function */
2509 /* corrSq[k]=corrSq[k]/count[k]; */
2511 /* normalize such that corr[0] == 0 */
2513 for (k = 0; (k < ncorr); k++)
2516 /* corrSq[k]*=c0*c0; */
2519 /* write ACFs in verbose mode */
2522 for (k = 0; (k < ncorr); k++)
2524 fprintf(fpcorr, "%g %g\n", k*dt, corr[k]);
2526 fprintf(fpcorr, "&\n");
2529 /* esimate integrated correlation time, fitting is too unstable */
2530 tausteps = 0.5*corr[0];
2531 /* consider corr below WHAM_AC_ZERO_LIMIT as noise */
2532 for (j = 1; (j < ncorr) && (corr[j] > WHAM_AC_ZERO_LIMIT); j++)
2534 tausteps += corr[j];
2537 /* g = 1+2*tau, see. Ferrenberg/Swendsen, PRL 63:1195 (1989) or
2538 Kumar et al, eq. 28 ff. */
2539 window[i].tau[ig] = tausteps*dt;
2540 window[i].g[ig] = 1+2*tausteps;
2541 /* printf("win %d, group %d, estimated correlation time = %g ps\n",i,ig,window[i].tau[ig]); */
2552 /* plot IACT along reaction coordinate */
2553 fp = xvgropen(fn, "Integrated autocorrelation times", "z", "IACT [ps]", opt->oenv);
2554 fprintf(fp, "@ s0 symbol 1\n@ s0 symbol size 0.5\n@ s0 line linestyle 0\n");
2555 fprintf(fp, "# WIN tau(gr1) tau(gr2) ...\n");
2556 for (i = 0; i < nwins; i++)
2558 fprintf(fp, "# %3d ", i);
2559 for (ig = 0; ig < window[i].nPull; ig++)
2561 fprintf(fp, " %11g", window[i].tau[ig]);
2565 for (i = 0; i < nwins; i++)
2567 for (ig = 0; ig < window[i].nPull; ig++)
2569 fprintf(fp, "%8g %8g\n", window[i].pos[ig], window[i].tau[ig]);
2572 if (opt->sigSmoothIact > 0.0)
2574 printf("Smoothing autocorrelation times along reaction coordinate with Gaussian of sig = %g\n",
2575 opt->sigSmoothIact);
2576 /* smooth IACT along reaction coordinate and overwrite g=1+2tau */
2577 smoothIact(window, nwins, opt);
2578 fprintf(fp, "&\n@ s1 symbol 1\n@ s1 symbol size 0.5\n@ s1 line linestyle 0\n");
2579 fprintf(fp, "@ s1 symbol color 2\n");
2580 for (i = 0; i < nwins; i++)
2582 for (ig = 0; ig < window[i].nPull; ig++)
2584 fprintf(fp, "%8g %8g\n", window[i].pos[ig], window[i].tausmooth[ig]);
2589 printf("Wrote %s\n", fn);
2592 /* compute average and sigma of each umbrella window */
2593 void averageSigma(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt)
2596 real av, sum2, sig, diff, *ztime, nSamplesIndep;
2598 for (i = 0; i < nwins; i++)
2600 snew(window[i].aver, window[i].nPull);
2601 snew(window[i].sigma, window[i].nPull);
2603 ntot = window[i].Ntot[0];
2604 for (ig = 0; ig < window[i].nPull; ig++)
2606 ztime = window[i].ztime[ig];
2607 for (k = 0, av = 0.; k < ntot; k++)
2612 for (k = 0, sum2 = 0.; k < ntot; k++)
2617 sig = sqrt(sum2/ntot);
2618 window[i].aver[ig] = av;
2620 /* Note: This estimate for sigma is biased from the limited sampling.
2621 Correct sigma by n/(n-1) where n = number of independent
2622 samples. Only possible if IACT is known.
2626 nSamplesIndep = window[i].N[ig]/(window[i].tau[ig]/window[i].dt);
2627 window[i].sigma[ig] = sig * nSamplesIndep/(nSamplesIndep-1);
2631 window[i].sigma[ig] = sig;
2633 printf("win %d, aver = %f sig = %f\n", i, av, window[i].sigma[ig]);
2639 /* Use histograms to compute average force on pull group.
2640 In addition, compute the sigma of the histogram.
2642 void computeAverageForce(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt)
2644 int i, j, bins = opt->bins, k;
2645 double dz, min = opt->min, max = opt->max, displAv, displAv2, temp, distance, ztot, ztot_half, w, weight;
2648 dz = (max-min)/bins;
2649 ztot = opt->max-min;
2652 /* Compute average displacement from histograms */
2653 for (j = 0; j < nWindows; ++j)
2655 snew(window[j].forceAv, window[j].nPull);
2656 for (k = 0; k < window[j].nPull; ++k)
2661 for (i = 0; i < opt->bins; ++i)
2663 temp = (1.0*i+0.5)*dz+min;
2664 distance = temp - window[j].pos[k];
2666 { /* in cyclic wham: */
2667 if (distance > ztot_half) /* |distance| < ztot_half */
2671 else if (distance < -ztot_half)
2676 w = window[j].Histo[k][i]/window[j].g[k];
2677 displAv += w*distance;
2678 displAv2 += w*sqr(distance);
2680 /* Are we near min or max? We are getting wron forces from the histgrams since
2681 the histigrams are zero outside [min,max). Therefore, assume that the position
2682 on the other side of the histomgram center is equally likely. */
2685 posmirrored = window[j].pos[k]-distance;
2686 if (posmirrored >= max || posmirrored < min)
2688 displAv += -w*distance;
2689 displAv2 += w*sqr(-distance);
2697 /* average force from average displacement */
2698 window[j].forceAv[k] = displAv*window[j].k[k];
2699 /* sigma from average square displacement */
2700 /* window[j].sigma [k] = sqrt(displAv2); */
2701 /* printf("Win %d, sigma = %f\n",j,sqrt(displAv2)); */
2706 /* Check if the complete reaction coordinate is covered by the histograms */
2707 void checkReactionCoordinateCovered(t_UmbrellaWindow *window, int nwins,
2708 t_UmbrellaOptions *opt)
2710 int i, ig, j, bins = opt->bins, bBoundary;
2711 real avcount = 0, z, relcount, *count;
2712 snew(count, opt->bins);
2714 for (j = 0; j < opt->bins; ++j)
2716 for (i = 0; i < nwins; i++)
2718 for (ig = 0; ig < window[i].nPull; ig++)
2720 count[j] += window[i].Histo[ig][j];
2723 avcount += 1.0*count[j];
2726 for (j = 0; j < bins; ++j)
2728 relcount = count[j]/avcount;
2729 z = (j+0.5)*opt->dz+opt->min;
2730 bBoundary = ( j<bins/20 || (bins-j)>bins/20 );
2731 /* check for bins with no data */
2734 fprintf(stderr, "\nWARNING, no data point in bin %d (z=%g) !\n"
2735 "You may not get a reasonable profile. Check your histograms!\n", j, z);
2737 /* and check for poor sampling */
2738 else if (relcount < 0.005 && !bBoundary)
2740 fprintf(stderr, "Warning, poor sampling bin %d (z=%g). Check your histograms!\n", j, z);
2747 void guessPotByIntegration(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt,
2750 int i, j, ig, bins = opt->bins, nHist, winmin, groupmin;
2751 double dz, min = opt->min, *pot, pos, hispos, dist, diff, fAv, distmin, *f;
2754 dz = (opt->max-min)/bins;
2756 printf("Getting initial potential by integration.\n");
2758 /* Compute average displacement from histograms */
2759 computeAverageForce(window, nWindows, opt);
2761 /* Get force for each bin from all histograms in this bin, or, alternatively,
2762 if no histograms are inside this bin, from the closest histogram */
2765 for (j = 0; j < opt->bins; ++j)
2767 pos = (1.0*j+0.5)*dz+min;
2771 groupmin = winmin = 0;
2772 for (i = 0; i < nWindows; i++)
2774 for (ig = 0; ig < window[i].nPull; ig++)
2776 hispos = window[i].pos[ig];
2777 dist = fabs(hispos-pos);
2778 /* average force within bin */
2782 fAv += window[i].forceAv[ig];
2784 /* at the same time, rememer closest histogram */
2793 /* if no histogram found in this bin, use closest histogram */
2800 fAv = window[winmin].forceAv[groupmin];
2804 for (j = 1; j < opt->bins; ++j)
2806 pot[j] = pot[j-1] - 0.5*dz*(f[j-1]+f[j]);
2809 /* cyclic wham: linearly correct possible offset */
2812 diff = (pot[bins-1]-pot[0])/(bins-1);
2813 for (j = 1; j < opt->bins; ++j)
2820 fp = xvgropen("pmfintegrated.xvg", "PMF from force integration", "z", "PMF [kJ/mol]", opt->oenv);
2821 for (j = 0; j < opt->bins; ++j)
2823 fprintf(fp, "%g %g\n", (j+0.5)*dz+opt->min, pot[j]);
2826 printf("verbose mode: wrote %s with PMF from interated forces\n", "pmfintegrated.xvg");
2829 /* get initial z=exp(-F[i]/kT) from integrated potential, where F[i] denote the free
2830 energy offsets which are usually determined by wham
2831 First: turn pot into probabilities:
2833 for (j = 0; j < opt->bins; ++j)
2835 pot[j] = exp(-pot[j]/(8.314e-3*opt->Temperature));
2837 calc_z(pot, window, nWindows, opt, TRUE);
2844 int gmx_wham(int argc, char *argv[])
2846 const char *desc[] = {
2847 "This is an analysis program that implements the Weighted",
2848 "Histogram Analysis Method (WHAM). It is intended to analyze",
2849 "output files generated by umbrella sampling simulations to ",
2850 "compute a potential of mean force (PMF). [PAR] ",
2851 "At present, three input modes are supported.[BR]",
2852 "[TT]*[tt] With option [TT]-it[tt], the user provides a file which contains the",
2853 " file names of the umbrella simulation run-input files ([TT].tpr[tt] files),",
2854 " AND, with option [TT]-ix[tt], a file which contains file names of",
2855 " the pullx [TT]mdrun[tt] output files. The [TT].tpr[tt] and pullx files must",
2856 " be in corresponding order, i.e. the first [TT].tpr[tt] created the",
2857 " first pullx, etc.[BR]",
2858 "[TT]*[tt] Same as the previous input mode, except that the the user",
2859 " provides the pull force output file names ([TT]pullf.xvg[tt]) with option [TT]-if[tt].",
2860 " From the pull force the position in the umbrella potential is",
2861 " computed. This does not work with tabulated umbrella potentials.[BR]"
2862 "[TT]*[tt] With option [TT]-ip[tt], the user provides file names of (gzipped) [TT].pdo[tt] files, i.e.",
2863 " the GROMACS 3.3 umbrella output files. If you have some unusual"
2864 " reaction coordinate you may also generate your own [TT].pdo[tt] files and",
2865 " feed them with the [TT]-ip[tt] option into to [TT]g_wham[tt]. The [TT].pdo[tt] file header",
2866 " must be similar to the following:[PAR]",
2867 "[TT]# UMBRELLA 3.0[BR]",
2868 "# Component selection: 0 0 1[BR]",
2870 "# Ref. Group 'TestAtom'[BR]",
2871 "# Nr. of pull groups 2[BR]",
2872 "# Group 1 'GR1' Umb. Pos. 5.0 Umb. Cons. 1000.0[BR]",
2873 "# Group 2 'GR2' Umb. Pos. 2.0 Umb. Cons. 500.0[BR]",
2875 "The number of pull groups, umbrella positions, force constants, and names ",
2876 "may (of course) differ. Following the header, a time column and ",
2877 "a data column for each pull group follows (i.e. the displacement",
2878 "with respect to the umbrella center). Up to four pull groups are possible ",
2879 "per [TT].pdo[tt] file at present.[PAR]",
2880 "By default, the output files are[BR]",
2881 " [TT]-o[tt] PMF output file[BR]",
2882 " [TT]-hist[tt] Histograms output file[BR]",
2883 "Always check whether the histograms sufficiently overlap.[PAR]",
2884 "The umbrella potential is assumed to be harmonic and the force constants are ",
2885 "read from the [TT].tpr[tt] or [TT].pdo[tt] files. If a non-harmonic umbrella force was applied ",
2886 "a tabulated potential can be provided with [TT]-tab[tt].[PAR]",
2887 "WHAM OPTIONS[BR]------------[BR]",
2888 " [TT]-bins[tt] Number of bins used in analysis[BR]",
2889 " [TT]-temp[tt] Temperature in the simulations[BR]",
2890 " [TT]-tol[tt] Stop iteration if profile (probability) changed less than tolerance[BR]",
2891 " [TT]-auto[tt] Automatic determination of boundaries[BR]",
2892 " [TT]-min,-max[tt] Boundaries of the profile [BR]",
2893 "The data points that are used to compute the profile",
2894 "can be restricted with options [TT]-b[tt], [TT]-e[tt], and [TT]-dt[tt]. ",
2895 "Adjust [TT]-b[tt] to ensure sufficient equilibration in each ",
2896 "umbrella window.[PAR]",
2897 "With [TT]-log[tt] (default) the profile is written in energy units, otherwise ",
2898 "(with [TT]-nolog[tt]) as probability. The unit can be specified with [TT]-unit[tt]. ",
2899 "With energy output, the energy in the first bin is defined to be zero. ",
2900 "If you want the free energy at a different ",
2901 "position to be zero, set [TT]-zprof0[tt] (useful with bootstrapping, see below).[PAR]",
2902 "For cyclic or periodic reaction coordinates (dihedral angle, channel PMF",
2903 "without osmotic gradient), the option [TT]-cycl[tt] is useful. [TT]g_wham[tt] will make use of the ",
2904 "periodicity of the system and generate a periodic PMF. The first and the last bin of the",
2905 "reaction coordinate will assumed be be neighbors.[PAR]",
2906 "Option [TT]-sym[tt] symmetrizes the profile around z=0 before output, ",
2907 "which may be useful for, e.g. membranes.[PAR]",
2908 "AUTOCORRELATIONS[BR]----------------[BR]",
2909 "With [TT]-ac[tt], [TT]g_wham[tt] estimates the integrated autocorrelation ",
2910 "time (IACT) [GRK]tau[grk] for each umbrella window and weights the respective ",
2911 "window with 1/[1+2*[GRK]tau[grk]/dt]. The IACTs are written ",
2912 "to the file defined with [TT]-oiact[tt]. In verbose mode, all ",
2913 "autocorrelation functions (ACFs) are written to [TT]hist_autocorr.xvg[tt]. ",
2914 "Because the IACTs can be severely underestimated in case of limited ",
2915 "sampling, option [TT]-acsig[tt] allows one to smooth the IACTs along the ",
2916 "reaction coordinate with a Gaussian ([GRK]sigma[grk] provided with [TT]-acsig[tt], ",
2917 "see output in [TT]iact.xvg[tt]). Note that the IACTs are estimated by simple ",
2918 "integration of the ACFs while the ACFs are larger 0.05.",
2919 "If you prefer to compute the IACTs by a more sophisticated (but possibly ",
2920 "less robust) method such as fitting to a double exponential, you can ",
2921 "compute the IACTs with [TT]g_analyze[tt] and provide them to [TT]g_wham[tt] with the file ",
2922 "[TT]iact-in.dat[tt] (option [TT]-iiact[tt]), which should contain one line per ",
2923 "input file ([TT].pdo[tt] or pullx/f file) and one column per pull group in the respective file.[PAR]",
2924 "ERROR ANALYSIS[BR]--------------[BR]",
2925 "Statistical errors may be estimated with bootstrap analysis. Use it with care, ",
2926 "otherwise the statistical error may be substantially underestimated. ",
2927 "More background and examples for the bootstrap technique can be found in ",
2928 "Hub, de Groot and Van der Spoel, JCTC (2010) 6: 3713-3720.[BR]",
2929 "[TT]-nBootstrap[tt] defines the number of bootstraps (use, e.g., 100). ",
2930 "Four bootstrapping methods are supported and ",
2931 "selected with [TT]-bs-method[tt].[BR]",
2932 " (1) [TT]b-hist[tt] Default: complete histograms are considered as independent ",
2933 "data points, and the bootstrap is carried out by assigning random weights to the ",
2934 "histograms (\"Bayesian bootstrap\"). Note that each point along the reaction coordinate",
2935 "must be covered by multiple independent histograms (e.g. 10 histograms), otherwise the ",
2936 "statistical error is underestimated.[BR]",
2937 " (2) [TT]hist[tt] Complete histograms are considered as independent data points. ",
2938 "For each bootstrap, N histograms are randomly chosen from the N given histograms ",
2939 "(allowing duplication, i.e. sampling with replacement).",
2940 "To avoid gaps without data along the reaction coordinate blocks of histograms ",
2941 "([TT]-histbs-block[tt]) may be defined. In that case, the given histograms are ",
2942 "divided into blocks and only histograms within each block are mixed. Note that ",
2943 "the histograms within each block must be representative for all possible histograms, ",
2944 "otherwise the statistical error is underestimated.[BR]",
2945 " (3) [TT]traj[tt] The given histograms are used to generate new random trajectories,",
2946 "such that the generated data points are distributed according the given histograms ",
2947 "and properly autocorrelated. The autocorrelation time (ACT) for each window must be ",
2948 "known, so use [TT]-ac[tt] or provide the ACT with [TT]-iiact[tt]. If the ACT of all ",
2949 "windows are identical (and known), you can also provide them with [TT]-bs-tau[tt]. ",
2950 "Note that this method may severely underestimate the error in case of limited sampling, ",
2951 "that is if individual histograms do not represent the complete phase space at ",
2952 "the respective positions.[BR]",
2953 " (4) [TT]traj-gauss[tt] The same as method [TT]traj[tt], but the trajectories are ",
2954 "not bootstrapped from the umbrella histograms but from Gaussians with the average ",
2955 "and width of the umbrella histograms. That method yields similar error estimates ",
2956 "like method [TT]traj[tt].[PAR]"
2957 "Bootstrapping output:[BR]",
2958 " [TT]-bsres[tt] Average profile and standard deviations[BR]",
2959 " [TT]-bsprof[tt] All bootstrapping profiles[BR]",
2960 "With [TT]-vbs[tt] (verbose bootstrapping), the histograms of each bootstrap are written, ",
2961 "and, with bootstrap method [TT]traj[tt], the cumulative distribution functions of ",
2965 const char *en_unit[] = {NULL, "kJ", "kCal", "kT", NULL};
2966 const char *en_unit_label[] = {"", "E (kJ mol\\S-1\\N)", "E (kcal mol\\S-1\\N)", "E (kT)", NULL};
2967 const char *en_bsMethod[] = { NULL, "b-hist", "hist", "traj", "traj-gauss", NULL };
2969 static t_UmbrellaOptions opt;
2972 { "-min", FALSE, etREAL, {&opt.min},
2973 "Minimum coordinate in profile"},
2974 { "-max", FALSE, etREAL, {&opt.max},
2975 "Maximum coordinate in profile"},
2976 { "-auto", FALSE, etBOOL, {&opt.bAuto},
2977 "Determine min and max automatically"},
2978 { "-bins", FALSE, etINT, {&opt.bins},
2979 "Number of bins in profile"},
2980 { "-temp", FALSE, etREAL, {&opt.Temperature},
2982 { "-tol", FALSE, etREAL, {&opt.Tolerance},
2984 { "-v", FALSE, etBOOL, {&opt.verbose},
2986 { "-b", FALSE, etREAL, {&opt.tmin},
2987 "First time to analyse (ps)"},
2988 { "-e", FALSE, etREAL, {&opt.tmax},
2989 "Last time to analyse (ps)"},
2990 { "-dt", FALSE, etREAL, {&opt.dt},
2991 "Analyse only every dt ps"},
2992 { "-histonly", FALSE, etBOOL, {&opt.bHistOnly},
2993 "Write histograms and exit"},
2994 { "-boundsonly", FALSE, etBOOL, {&opt.bBoundsOnly},
2995 "Determine min and max and exit (with [TT]-auto[tt])"},
2996 { "-log", FALSE, etBOOL, {&opt.bLog},
2997 "Calculate the log of the profile before printing"},
2998 { "-unit", FALSE, etENUM, {en_unit},
2999 "Energy unit in case of log output" },
3000 { "-zprof0", FALSE, etREAL, {&opt.zProf0},
3001 "Define profile to 0.0 at this position (with [TT]-log[tt])"},
3002 { "-cycl", FALSE, etBOOL, {&opt.bCycl},
3003 "Create cyclic/periodic profile. Assumes min and max are the same point."},
3004 { "-sym", FALSE, etBOOL, {&opt.bSym},
3005 "Symmetrize profile around z=0"},
3006 { "-hist-eq", FALSE, etBOOL, {&opt.bHistEq},
3007 "HIDDENEnforce equal weight for all histograms. (Non-Weighed-HAM)"},
3008 { "-ac", FALSE, etBOOL, {&opt.bCalcTauInt},
3009 "Calculate integrated autocorrelation times and use in wham"},
3010 { "-acsig", FALSE, etREAL, {&opt.sigSmoothIact},
3011 "Smooth autocorrelation times along reaction coordinate with Gaussian of this [GRK]sigma[grk]"},
3012 { "-ac-trestart", FALSE, etREAL, {&opt.acTrestart},
3013 "When computing autocorrelation functions, restart computing every .. (ps)"},
3014 { "-acred", FALSE, etBOOL, {&opt.bAllowReduceIact},
3015 "HIDDENWhen smoothing the ACTs, allow to reduce ACTs. Otherwise, only increase ACTs "
3016 "during smoothing"},
3017 { "-nBootstrap", FALSE, etINT, {&opt.nBootStrap},
3018 "nr of bootstraps to estimate statistical uncertainty (e.g., 200)" },
3019 { "-bs-method", FALSE, etENUM, {en_bsMethod},
3020 "Bootstrap method" },
3021 { "-bs-tau", FALSE, etREAL, {&opt.tauBootStrap},
3022 "Autocorrelation time (ACT) assumed for all histograms. Use option [TT]-ac[tt] if ACT is unknown."},
3023 { "-bs-seed", FALSE, etINT, {&opt.bsSeed},
3024 "Seed for bootstrapping. (-1 = use time)"},
3025 { "-histbs-block", FALSE, etINT, {&opt.histBootStrapBlockLength},
3026 "When mixing histograms only mix within blocks of [TT]-histbs-block[tt]."},
3027 { "-vbs", FALSE, etBOOL, {&opt.bs_verbose},
3028 "Verbose bootstrapping. Print the CDFs and a histogram file for each bootstrap."},
3029 { "-stepout", FALSE, etINT, {&opt.stepchange},
3030 "HIDDENWrite maximum change every ... (set to 1 with [TT]-v[tt])"},
3031 { "-updateContr", FALSE, etINT, {&opt.stepUpdateContrib},
3032 "HIDDENUpdate table with significan contributions to WHAM every ... iterations"},
3036 { efDAT, "-ix", "pullx-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs */
3037 { efDAT, "-if", "pullf-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs */
3038 { efDAT, "-it", "tpr-files", ffOPTRD}, /* wham input: tprs */
3039 { efDAT, "-ip", "pdo-files", ffOPTRD}, /* wham input: pdo files (gmx3 style) */
3040 { efXVG, "-o", "profile", ffWRITE }, /* output file for profile */
3041 { efXVG, "-hist", "histo", ffWRITE}, /* output file for histograms */
3042 { efXVG, "-oiact", "iact", ffOPTWR}, /* writing integrated autocorrelation times */
3043 { efDAT, "-iiact", "iact-in", ffOPTRD}, /* reading integrated autocorrelation times */
3044 { efXVG, "-bsres", "bsResult", ffOPTWR}, /* average and errors of bootstrap analysis */
3045 { efXVG, "-bsprof", "bsProfs", ffOPTWR}, /* output file for bootstrap profiles */
3046 { efDAT, "-tab", "umb-pot", ffOPTRD}, /* Tabulated umbrella potential (if not harmonic) */
3049 int i, j, l, nfiles, nwins, nfiles2;
3050 t_UmbrellaHeader header;
3051 t_UmbrellaWindow * window = NULL;
3052 double *profile, maxchange = 1e20;
3053 gmx_bool bMinSet, bMaxSet, bAutoSet, bExact = FALSE;
3054 char **fninTpr, **fninPull, **fninPdo;
3056 FILE *histout, *profout;
3057 char ylabel[256], title[256];
3059 #define NFILE asize(fnm)
3061 CopyRight(stderr, argv[0]);
3064 opt.verbose = FALSE;
3065 opt.bHistOnly = FALSE;
3074 /* bootstrapping stuff */
3076 opt.bsMethod = bsMethod_hist;
3077 opt.tauBootStrap = 0.0;
3079 opt.histBootStrapBlockLength = 8;
3080 opt.bs_verbose = FALSE;
3085 opt.Temperature = 298;
3086 opt.Tolerance = 1e-6;
3087 opt.bBoundsOnly = FALSE;
3089 opt.bCalcTauInt = FALSE;
3090 opt.sigSmoothIact = 0.0;
3091 opt.bAllowReduceIact = TRUE;
3092 opt.bInitPotByIntegration = TRUE;
3093 opt.acTrestart = 1.0;
3094 opt.stepchange = 100;
3095 opt.stepUpdateContrib = 100;
3097 parse_common_args(&argc, argv, PCA_BE_NICE,
3098 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &opt.oenv);
3100 opt.unit = nenum(en_unit);
3101 opt.bsMethod = nenum(en_bsMethod);
3103 opt.bProf0Set = opt2parg_bSet("-zprof0", asize(pa), pa);
3105 opt.bTab = opt2bSet("-tab", NFILE, fnm);
3106 opt.bPdo = opt2bSet("-ip", NFILE, fnm);
3107 opt.bTpr = opt2bSet("-it", NFILE, fnm);
3108 opt.bPullx = opt2bSet("-ix", NFILE, fnm);
3109 opt.bPullf = opt2bSet("-if", NFILE, fnm);
3110 opt.bTauIntGiven = opt2bSet("-iiact", NFILE, fnm);
3111 if (opt.bTab && opt.bPullf)
3113 gmx_fatal(FARGS, "Force input does not work with tabulated potentials. "
3114 "Provide pullx.xvg or pdo files!");
3117 #define WHAMBOOLXOR(a, b) ( ((!(a)) && (b)) || ((a) && (!(b))))
3118 if (!opt.bPdo && !WHAMBOOLXOR(opt.bPullx, opt.bPullf))
3120 gmx_fatal(FARGS, "Give either pullx (-ix) OR pullf (-if) data. Not both.");
3122 if (!opt.bPdo && !(opt.bTpr || opt.bPullf || opt.bPullx))
3124 gmx_fatal(FARGS, "g_wham supports three input modes, pullx, pullf, or pdo file input."
3125 "\n\n Check g_wham -h !");
3128 opt.fnPdo = opt2fn("-ip", NFILE, fnm);
3129 opt.fnTpr = opt2fn("-it", NFILE, fnm);
3130 opt.fnPullf = opt2fn("-if", NFILE, fnm);
3131 opt.fnPullx = opt2fn("-ix", NFILE, fnm);
3133 bMinSet = opt2parg_bSet("-min", asize(pa), pa);
3134 bMaxSet = opt2parg_bSet("-max", asize(pa), pa);
3135 bAutoSet = opt2parg_bSet("-auto", asize(pa), pa);
3136 if ( (bMinSet || bMaxSet) && bAutoSet)
3138 gmx_fatal(FARGS, "With -auto, do not give -min or -max\n");
3141 if ( (bMinSet && !bMaxSet) || (!bMinSet && bMaxSet))
3143 gmx_fatal(FARGS, "When giving -min, you must give -max (and vice versa), too\n");
3146 if (bMinSet && opt.bAuto)
3148 printf("Note: min and max given, switching off -auto.\n");
3152 if (opt.bTauIntGiven && opt.bCalcTauInt)
3154 gmx_fatal(FARGS, "Either read (option -iiact) or calculate (option -ac) the\n"
3155 "the autocorrelation times. Not both.");
3158 if (opt.tauBootStrap > 0.0 && opt2parg_bSet("-ac", asize(pa), pa))
3160 gmx_fatal(FARGS, "Either compute autocorrelation times (ACTs) (option -ac) or "
3161 "provide it with -bs-tau for bootstrapping. Not Both.\n");
3163 if (opt.tauBootStrap > 0.0 && opt2bSet("-iiact", NFILE, fnm))
3165 gmx_fatal(FARGS, "Either provide autocorrelation times (ACTs) with file iact-in.dat "
3166 "(option -iiact) or define all ACTs with -bs-tau for bootstrapping\n. Not Both.");
3170 /* Reading gmx4 pull output and tpr files */
3171 if (opt.bTpr || opt.bPullf || opt.bPullx)
3173 read_wham_in(opt.fnTpr, &fninTpr, &nfiles, &opt);
3175 fnPull = opt.bPullf ? opt.fnPullf : opt.fnPullx;
3176 read_wham_in(fnPull, &fninPull, &nfiles2, &opt);
3177 printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n",
3178 nfiles, nfiles2, opt.bPullf ? "force" : "position", opt.fnTpr, fnPull);
3179 if (nfiles != nfiles2)
3181 gmx_fatal(FARGS, "Found %d file names in %s, but %d in %s\n", nfiles,
3182 opt.fnTpr, nfiles2, fnPull);
3184 window = initUmbrellaWindows(nfiles);
3185 read_tpr_pullxf_files(fninTpr, fninPull, nfiles, &header, window, &opt);
3188 { /* reading pdo files */
3189 read_wham_in(opt.fnPdo, &fninPdo, &nfiles, &opt);
3190 printf("Found %d pdo files in %s\n", nfiles, opt.fnPdo);
3191 window = initUmbrellaWindows(nfiles);
3192 read_pdo_files(fninPdo, nfiles, &header, window, &opt);
3196 /* enforce equal weight for all histograms? */
3199 enforceEqualWeights(window, nwins);
3202 /* write histograms */
3203 histout = xvgropen(opt2fn("-hist", NFILE, fnm), "Umbrella histograms",
3204 "z", "count", opt.oenv);
3205 for (l = 0; l < opt.bins; ++l)
3207 fprintf(histout, "%e\t", (double)(l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
3208 for (i = 0; i < nwins; ++i)
3210 for (j = 0; j < window[i].nPull; ++j)
3212 fprintf(histout, "%e\t", window[i].Histo[j][l]);
3215 fprintf(histout, "\n");
3218 printf("Wrote %s\n", opt2fn("-hist", NFILE, fnm));
3221 printf("Wrote histograms to %s, now exiting.\n", opt2fn("-hist", NFILE, fnm));
3225 /* Using tabulated umbrella potential */
3228 setup_tab(opt2fn("-tab", NFILE, fnm), &opt);
3231 /* Integrated autocorrelation times provided ? */
3232 if (opt.bTauIntGiven)
3234 readIntegratedAutocorrelationTimes(window, nwins, &opt, opt2fn("-iiact", NFILE, fnm));
3237 /* Compute integrated autocorrelation times */
3238 if (opt.bCalcTauInt)
3240 calcIntegratedAutocorrelationTimes(window, nwins, &opt, opt2fn("-oiact", NFILE, fnm));
3243 /* calc average and sigma for each histogram
3244 (maybe required for bootstrapping. If not, this is fast anyhow) */
3245 if (opt.nBootStrap && opt.bsMethod == bsMethod_trajGauss)
3247 averageSigma(window, nwins, &opt);
3250 /* Get initial potential by simple integration */
3251 if (opt.bInitPotByIntegration)
3253 guessPotByIntegration(window, nwins, &opt, 0);
3256 /* Check if complete reaction coordinate is covered */
3257 checkReactionCoordinateCovered(window, nwins, &opt);
3259 /* Calculate profile */
3260 snew(profile, opt.bins);
3268 if ( (i%opt.stepUpdateContrib) == 0)
3270 setup_acc_wham(profile, window, nwins, &opt);
3272 if (maxchange < opt.Tolerance)
3275 /* if (opt.verbose) */
3276 printf("Switched to exact iteration in iteration %d\n", i);
3278 calc_profile(profile, window, nwins, &opt, bExact);
3279 if (((i%opt.stepchange) == 0 || i == 1) && !i == 0)
3281 printf("\t%4d) Maximum change %e\n", i, maxchange);
3285 while ( (maxchange = calc_z(profile, window, nwins, &opt, bExact)) > opt.Tolerance || !bExact);
3286 printf("Converged in %d iterations. Final maximum change %g\n", i, maxchange);
3288 /* calc error from Kumar's formula */
3289 /* Unclear how the error propagates along reaction coordinate, therefore
3291 /* calc_error_kumar(profile,window, nwins,&opt); */
3293 /* Write profile in energy units? */
3296 prof_normalization_and_unit(profile, &opt);
3297 strcpy(ylabel, en_unit_label[opt.unit]);
3298 strcpy(title, "Umbrella potential");
3302 strcpy(ylabel, "Density of states");
3303 strcpy(title, "Density of states");
3306 /* symmetrize profile around z=0? */
3309 symmetrizeProfile(profile, &opt);
3312 /* write profile or density of states */
3313 profout = xvgropen(opt2fn("-o", NFILE, fnm), title, "z", ylabel, opt.oenv);
3314 for (i = 0; i < opt.bins; ++i)
3316 fprintf(profout, "%e\t%e\n", (double)(i+0.5)/opt.bins*(opt.max-opt.min)+opt.min, profile[i]);
3319 printf("Wrote %s\n", opt2fn("-o", NFILE, fnm));
3321 /* Bootstrap Method */
3324 do_bootstrapping(opt2fn("-bsres", NFILE, fnm), opt2fn("-bsprof", NFILE, fnm),
3325 opt2fn("-hist", NFILE, fnm),
3326 ylabel, profile, window, nwins, &opt);
3330 freeUmbrellaWindows(window, nfiles);
3332 printf("\nIn case you use results from g_wham for a publication, please cite:\n");
3333 please_cite(stdout, "Hub2010");