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);
1599 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
1603 /* write average and stddev */
1604 fp = xvgropen(fnres, "Average and stddev from bootstrapping", "z", ylabel, opt->oenv);
1605 if(output_env_get_print_xvgr_codes(opt->oenv))
1607 fprintf(fp, "@TYPE xydy\n");
1609 for (i = 0; i < opt->bins; i++)
1611 bsProfiles_av [i] /= opt->nBootStrap;
1612 bsProfiles_av2[i] /= opt->nBootStrap;
1613 tmp = bsProfiles_av2[i]-sqr(bsProfiles_av[i]);
1614 stddev = (tmp >= 0.) ? sqrt(tmp) : 0.; /* Catch rouding errors */
1615 fprintf(fp, "%e\t%e\t%e\n", (i+0.5)*opt->dz+opt->min, bsProfiles_av [i], stddev);
1618 printf("Wrote boot strap result to %s\n", fnres);
1621 int whaminFileType(char *fn)
1625 if (strcmp(fn+len-3, "tpr") == 0)
1629 else if (strcmp(fn+len-3, "xvg") == 0 || strcmp(fn+len-6, "xvg.gz") == 0)
1631 return whamin_pullxf;
1633 else if (strcmp(fn+len-3, "pdo") == 0 || strcmp(fn+len-6, "pdo.gz") == 0)
1639 gmx_fatal(FARGS, "Unknown file type of %s. Should be tpr, xvg, or pdo.\n", fn);
1641 return whamin_unknown;
1644 void read_wham_in(const char *fn, char ***filenamesRet, int *nfilesRet,
1645 t_UmbrellaOptions *opt)
1647 char **filename = 0, tmp[STRLEN];
1648 int nread, sizenow, i, block = 1;
1651 fp = ffopen(fn, "r");
1654 while (fscanf(fp, "%s", tmp) != EOF)
1656 if (strlen(tmp) >= WHAM_MAXFILELEN)
1658 gmx_fatal(FARGS, "Filename too long. Only %d characters allowed\n", WHAM_MAXFILELEN);
1660 if (nread >= sizenow)
1663 srenew(filename, sizenow);
1664 for (i = sizenow-block; i < sizenow; i++)
1666 snew(filename[i], WHAM_MAXFILELEN);
1669 strcpy(filename[nread], tmp);
1672 printf("Found file %s in %s\n", filename[nread], fn);
1676 *filenamesRet = filename;
1681 FILE *open_pdo_pipe(const char *fn, t_UmbrellaOptions *opt, gmx_bool *bPipeOpen)
1683 char Buffer[1024], gunzip[1024], *Path = 0;
1685 static gmx_bool bFirst = 1;
1687 /* gzipped pdo file? */
1688 if ((strcmp(fn+strlen(fn)-3, ".gz") == 0))
1690 /* search gunzip executable */
1691 if (!(Path = getenv("GMX_PATH_GZIP")))
1693 if (gmx_fexist("/bin/gunzip"))
1695 sprintf(gunzip, "%s", "/bin/gunzip");
1697 else if (gmx_fexist("/usr/bin/gunzip"))
1699 sprintf(gunzip, "%s", "/usr/bin/gunzip");
1703 gmx_fatal(FARGS, "Cannot find executable gunzip in /bin or /usr/bin.\n"
1704 "You may want to define the path to gunzip "
1705 "with the environment variable GMX_PATH_GZIP.", gunzip);
1710 sprintf(gunzip, "%s/gunzip", Path);
1711 if (!gmx_fexist(gunzip))
1713 gmx_fatal(FARGS, "Cannot find executable %s. Please define the path to gunzip"
1714 " in the environmental varialbe GMX_PATH_GZIP.", gunzip);
1719 printf("Using gunzig executable %s\n", gunzip);
1722 if (!gmx_fexist(fn))
1724 gmx_fatal(FARGS, "File %s does not exist.\n", fn);
1726 sprintf(Buffer, "%s -c < %s", gunzip, fn);
1729 printf("Executing command '%s'\n", Buffer);
1732 if ((pipe = popen(Buffer, "r")) == NULL)
1734 gmx_fatal(FARGS, "Unable to open pipe to `%s'\n", Buffer);
1737 gmx_fatal(FARGS, "Cannot open a compressed file on platform without pipe support");
1743 pipe = ffopen(fn, "r");
1750 void pdo_close_file(FILE *fp)
1759 /* Reading pdo files */
1760 void read_pdo_files(char **fn, int nfiles, t_UmbrellaHeader* header,
1761 t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
1764 real mintmp, maxtmp, done = 0.;
1767 /* char Buffer0[1000]; */
1771 gmx_fatal(FARGS, "No files found. Hick.");
1774 /* if min and max are not given, get min and max from the input files */
1777 printf("Automatic determination of boundaries from %d pdo files...\n", nfiles);
1780 for (i = 0; i < nfiles; ++i)
1782 file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
1783 /*fgets(Buffer0,999,file);
1784 fprintf(stderr,"First line '%s'\n",Buffer0); */
1785 done = 100.0*(i+1)/nfiles;
1786 printf("\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
1791 read_pdo_header(file, header, opt);
1792 /* here only determine min and max of this window */
1793 read_pdo_data(file, header, i, NULL, opt, TRUE, &mintmp, &maxtmp);
1794 if (maxtmp > opt->max)
1798 if (mintmp < opt->min)
1804 pdo_close_file(file);
1812 printf("\nDetermined boundaries to %f and %f\n\n", opt->min, opt->max);
1813 if (opt->bBoundsOnly)
1815 printf("Found option -boundsonly, now exiting.\n");
1819 /* store stepsize in profile */
1820 opt->dz = (opt->max-opt->min)/opt->bins;
1822 /* Having min and max, we read in all files */
1823 /* Loop over all files */
1824 for (i = 0; i < nfiles; ++i)
1826 done = 100.0*(i+1)/nfiles;
1827 printf("\rOpening %s ... [%2.0f%%]", fn[i], done); fflush(stdout);
1832 file = open_pdo_pipe(fn[i], opt, &bPipeOpen);
1833 read_pdo_header(file, header, opt);
1834 /* load data into window */
1835 read_pdo_data(file, header, i, window, opt, FALSE, NULL, NULL);
1836 if ((window+i)->Ntot[0] == 0.0)
1838 fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fn[i]);
1842 pdo_close_file(file);
1850 for (i = 0; i < nfiles; ++i)
1857 #define int2YN(a) (((a) == 0) ? ("N") : ("Y"))
1859 void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions *opt)
1864 static int first = 1;
1866 /* printf("Reading %s \n",fn); */
1867 read_tpx_state(fn, &ir, &state, NULL, NULL);
1869 if (ir.ePull != epullUMBRELLA)
1871 gmx_fatal(FARGS, "This is not a tpr of an umbrella simulation. Found pull type \"%s\" "
1872 " (ir.ePull = %d)\n", epull_names[ir.ePull], ir.ePull);
1875 /* nr of pull groups */
1876 ngrp = ir.pull->ngrp;
1879 gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull groups\n", ngrp);
1882 header->npullgrps = ir.pull->ngrp;
1883 header->pull_geometry = ir.pull->eGeom;
1884 copy_ivec(ir.pull->dim, header->pull_dim);
1885 header->pull_ndim = header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2];
1886 if (header->pull_geometry == epullgPOS && header->pull_ndim > 1)
1888 gmx_fatal(FARGS, "Found pull geometry 'position' and more than 1 pull dimension (%d).\n"
1889 "Hence, the pull potential does not correspond to a one-dimensional umbrella potential.\n"
1890 "If you have some special umbrella setup you may want to write your own pdo files\n"
1891 "and feed them into g_wham. Check g_wham -h !\n", header->pull_ndim);
1893 snew(header->k, ngrp);
1894 snew(header->init_dist, ngrp);
1895 snew(header->umbInitDist, ngrp);
1897 /* only z-direction with epullgCYL? */
1898 if (header->pull_geometry == epullgCYL)
1900 if (header->pull_dim[XX] || header->pull_dim[YY] || (!header->pull_dim[ZZ]))
1902 gmx_fatal(FARGS, "With pull geometry 'cylinder', expected pulling in Z direction only.\n"
1903 "However, found dimensions [%s %s %s]\n",
1904 int2YN(header->pull_dim[XX]), int2YN(header->pull_dim[YY]),
1905 int2YN(header->pull_dim[ZZ]));
1909 for (i = 0; i < ngrp; i++)
1911 header->k[i] = ir.pull->grp[i+1].k;
1912 if (header->k[i] == 0.0)
1914 gmx_fatal(FARGS, "Pull group %d has force constant of of 0.0 in %s.\n"
1915 "That doesn't seem to be an Umbrella tpr.\n",
1918 copy_rvec(ir.pull->grp[i+1].init, header->init_dist[i]);
1920 /* initial distance to reference */
1921 switch (header->pull_geometry)
1924 for (d = 0; d < DIM; d++)
1926 if (header->pull_dim[d])
1928 header->umbInitDist[i] = header->init_dist[i][d];
1933 /* umbrella distance stored in init_dist[i][0] for geometry cylinder (not in ...[i][ZZ]) */
1937 header->umbInitDist[i] = header->init_dist[i][0];
1940 gmx_fatal(FARGS, "Pull geometry %s not supported\n", epullg_names[header->pull_geometry]);
1944 if (opt->verbose || first)
1946 printf("File %s, %d groups, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n",
1947 fn, header->npullgrps, epullg_names[header->pull_geometry],
1948 int2YN(header->pull_dim[0]), int2YN(header->pull_dim[1]), int2YN(header->pull_dim[2]),
1950 for (i = 0; i < ngrp; i++)
1952 printf("\tgrp %d) k = %-5g position = %g\n", i, header->k[i], header->umbInitDist[i]);
1955 if (!opt->verbose && first)
1957 printf("\tUse option -v to see this output for all input tpr files\n");
1964 double dist_ndim(double **dx, int ndim, int line)
1968 for (i = 0; i < ndim; i++)
1970 r2 += sqr(dx[i][line]);
1975 void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header,
1976 t_UmbrellaWindow * window,
1977 t_UmbrellaOptions *opt,
1978 gmx_bool bGetMinMax, real *mintmp, real *maxtmp)
1980 double **y = 0, pos = 0., t, force, time0 = 0., dt;
1981 int ny, nt, bins, ibin, i, g, dstep = 1, nColPerGrp, nColRefOnce, nColRefEachGrp, nColExpect, ntot;
1982 real min, max, minfound = 1e20, maxfound = -1e20;
1983 gmx_bool dt_ok, timeok, bHaveForce;
1984 const char *quantity;
1985 const int blocklen = 4096;
1989 in force output pullf.xvg:
1990 No reference, one column per pull group
1991 in position output pullx.xvg (not cylinder)
1992 ndim reference, ndim columns per pull group
1993 in position output pullx.xvg (in geometry cylinder):
1994 ndim*2 columns per pull group (ndim for ref, ndim for group)
1997 nColPerGrp = opt->bPullx ? header->pull_ndim : 1;
1998 quantity = opt->bPullx ? "position" : "force";
2002 if (header->pull_geometry == epullgCYL)
2004 /* Geometry cylinder -> reference group before each pull group */
2005 nColRefEachGrp = header->pull_ndim;
2010 /* Geometry NOT cylinder -> reference group only once after time column */
2012 nColRefOnce = header->pull_ndim;
2015 else /* read forces, no reference groups */
2021 nColExpect = 1 + nColRefOnce + header->npullgrps*(nColRefEachGrp+nColPerGrp);
2022 bHaveForce = opt->bPullf;
2024 /* With geometry "distance" or "distance_periodic", only force reading is supported so far.
2025 That avoids the somewhat tedious extraction of the right columns from the pullx files
2026 to compute the distances projection on the vector. Sorry for the laziness. */
2027 if ( (header->pull_geometry == epullgDIR || header->pull_geometry == epullgDIRPBC)
2030 gmx_fatal(FARGS, "With pull geometries \"direction\" and \"direction_periodic\", only pull force "
2031 "reading \n(option -if) is supported at present, "
2032 "not pull position reading (options -ix).\nMake sure mdrun writes the pull "
2033 "forces (pullf.xvg files)\nand provide them to g_wham with option -if.",
2034 epullg_names[header->pull_geometry]);
2037 nt = read_xvg(fn, &y, &ny);
2039 /* Check consistency */
2042 gmx_fatal(FARGS, "Empty pull %s file %s\n", quantity, fn);
2044 if (ny != nColExpect)
2046 gmx_fatal(FARGS, "Found %d pull groups in %s,\n but %d data columns in %s (expected %d)\n"
2047 "\nMaybe you confused options -ix and -if ?\n",
2048 header->npullgrps, fntpr, ny-1, fn, nColExpect-1);
2053 printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerGrp, quantity, fn);
2063 window->dt = y[0][1]-y[0][0];
2065 else if (opt->nBootStrap && opt->tauBootStrap != 0.0)
2067 fprintf(stderr, "\n *** WARNING, Could not determine time step in %s\n", fn);
2070 /* Need to alocate memory and set up structure */
2071 window->nPull = header->npullgrps;
2072 window->nBin = bins;
2074 snew(window->Histo, window->nPull);
2075 snew(window->z, window->nPull);
2076 snew(window->k, window->nPull);
2077 snew(window->pos, window->nPull);
2078 snew(window->N, window->nPull);
2079 snew(window->Ntot, window->nPull);
2080 snew(window->g, window->nPull);
2081 snew(window->bsWeight, window->nPull);
2082 window->bContrib = 0;
2084 if (opt->bCalcTauInt)
2086 snew(window->ztime, window->nPull);
2090 window->ztime = NULL;
2092 snew(lennow, window->nPull);
2094 for (g = 0; g < window->nPull; ++g)
2097 window->bsWeight[g] = 1.;
2098 snew(window->Histo[g], bins);
2099 window->k[g] = header->k[g];
2101 window->Ntot[g] = 0;
2103 window->pos[g] = header->umbInitDist[g];
2104 if (opt->bCalcTauInt)
2106 window->ztime[g] = NULL;
2112 { /* only determine min and max */
2115 min = max = bins = 0; /* Get rid of warnings */
2118 for (i = 0; i < nt; i++)
2120 /* Do you want that time frame? */
2121 t = 1.0/1000*( (int) ((y[0][i]*1000) + 0.5)); /* round time to fs */
2123 /* get time step of pdo file and get dstep from opt->dt */
2133 dstep = (int)(opt->dt/dt+0.5);
2141 window->dt = dt*dstep;
2145 dt_ok = (i%dstep == 0);
2146 timeok = (dt_ok && t >= opt->tmin && t <= opt->tmax);
2148 printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
2149 t,opt->tmin, opt->tmax, dt_ok,timeok); */
2153 for (g = 0; g < header->npullgrps; ++g)
2157 /* y has 1 time column y[0] and one column per force y[1],...,y[nGrps] */
2159 pos = -force/header->k[g] + header->umbInitDist[g];
2163 switch (header->pull_geometry)
2166 /* y has 1 time column y[0] and nColPerGrps columns per pull group;
2167 Distance to reference: */
2168 /* pos=dist_ndim(y+1+nColRef+g*nColPerGrp,header->pull_ndim,i); gmx 4.0 */
2169 pos = dist_ndim(y + 1 + nColRefOnce + g*nColPerGrp, header->pull_ndim, i);
2173 Time ref[ndim] group1[ndim] group2[ndim] ... */
2176 Time ref1[ndim] group1[ndim] ref2[ndim] group2[ndim] ... */
2178 /* * with geometry==position, we have the reference once (nColRefOnce==ndim), but
2179 no extra reference group columns before each group (nColRefEachGrp==0)
2181 * with geometry==cylinder, we have no initial ref group column (nColRefOnce==0),
2182 but ndim ref group colums before every group (nColRefEachGrp==ndim)
2183 Distance to reference: */
2184 pos = y[1 + nColRefOnce + nColRefEachGrp + g*(nColRefEachGrp+nColPerGrp)][i];
2187 gmx_fatal(FARGS, "Bad error, this error should have been catched before. Ups.\n");
2191 /* printf("grp %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
2205 if (opt->bCalcTauInt && !bGetMinMax)
2207 /* save time series for autocorrelation analysis */
2208 ntot = window->Ntot[g];
2209 /* printf("i %d, ntot %d, lennow[g] = %d\n",i,ntot,lennow[g]); */
2210 if (ntot >= lennow[g])
2212 lennow[g] += blocklen;
2213 srenew(window->ztime[g], lennow[g]);
2215 window->ztime[g][ntot] = pos;
2218 ibin = (int) floor((pos-min)/(max-min)*bins);
2223 while ( (ibin += bins) < 0)
2228 else if (ibin >= bins)
2230 while ( (ibin -= bins) >= bins)
2236 if (ibin >= 0 && ibin < bins)
2238 window->Histo[g][ibin] += 1.;
2245 else if (t > opt->tmax)
2249 printf("time %f larger than tmax %f, stop reading this pullx/pullf file\n", t, opt->tmax);
2261 for (i = 0; i < ny; i++)
2267 void read_tpr_pullxf_files(char **fnTprs, char **fnPull, int nfiles,
2268 t_UmbrellaHeader* header,
2269 t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
2272 real mintmp, maxtmp;
2274 printf("Reading %d tpr and pullf files\n", nfiles/2);
2276 /* min and max not given? */
2279 printf("Automatic determination of boundaries...\n");
2282 for (i = 0; i < nfiles; i++)
2284 if (whaminFileType(fnTprs[i]) != whamin_tpr)
2286 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a tpr file\n", i);
2288 read_tpr_header(fnTprs[i], header, opt);
2289 if (whaminFileType(fnPull[i]) != whamin_pullxf)
2291 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
2293 read_pull_xf(fnPull[i], fnTprs[i], header, NULL, opt, TRUE, &mintmp, &maxtmp);
2294 if (maxtmp > opt->max)
2298 if (mintmp < opt->min)
2303 printf("\nDetermined boundaries to %f and %f\n\n", opt->min, opt->max);
2304 if (opt->bBoundsOnly)
2306 printf("Found option -boundsonly, now exiting.\n");
2310 /* store stepsize in profile */
2311 opt->dz = (opt->max-opt->min)/opt->bins;
2313 for (i = 0; i < nfiles; i++)
2315 if (whaminFileType(fnTprs[i]) != whamin_tpr)
2317 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a tpr file\n", i);
2319 read_tpr_header(fnTprs[i], header, opt);
2320 if (whaminFileType(fnPull[i]) != whamin_pullxf)
2322 gmx_fatal(FARGS, "Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n", i);
2324 read_pull_xf(fnPull[i], fnTprs[i], header, window+i, opt, FALSE, NULL, NULL);
2325 if (window[i].Ntot[0] == 0.0)
2327 fprintf(stderr, "\nWARNING, no data points read from file %s (check -b option)\n", fnPull[i]);
2331 for (i = 0; i < nfiles; i++)
2340 /* Note: Here we consider tau[int] := int_0^inf ACF(t) as the integrated autocorrelation time.
2341 The factor `g := 1 + 2*tau[int]` subsequently enters the uncertainty.
2343 void readIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt,
2346 int nlines, ny, i, ig;
2349 printf("Readging Integrated autocorrelation times from %s ...\n", fn);
2350 nlines = read_xvg(fn, &iact, &ny);
2351 if (nlines != nwins)
2353 gmx_fatal(FARGS, "Found %d lines with integrated autocorrelation times in %s.\nExpected %d",
2356 for (i = 0; i < nlines; i++)
2358 if (window[i].nPull != ny)
2360 gmx_fatal(FARGS, "You are providing autocorrelation times with option -iiact and the\n"
2361 "number of pull groups is different in different simulations. That is not\n"
2362 "supported yet. Sorry.\n");
2364 for (ig = 0; ig < window[i].nPull; ig++)
2366 /* compare Kumar et al, J Comp Chem 13, 1011-1021 (1992) */
2367 window[i].g[ig] = 1+2*iact[ig][i]/window[i].dt;
2369 if (iact[ig][i] <= 0.0)
2371 fprintf(stderr, "\nWARNING, IACT = %f (window %d, group %d)\n", iact[ig][i], i, ig);
2378 /* Smooth autocorreltion times along the reaction coordinate. This is useful
2379 if the ACT is subject to high uncertainty in case if limited sampling. Note
2380 that -in case of limited sampling- the ACT may be severely underestimated.
2381 Note: the g=1+2tau are overwritten.
2382 if opt->bAllowReduceIact==FALSE, the ACTs are never reduced, only increased
2385 void smoothIact(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt)
2388 double pos, dpos2, siglim, siglim2, gaufact, invtwosig2, w, weight, tausm;
2390 /* only evaluate within +- 3sigma of the Gausian */
2391 siglim = 3.0*opt->sigSmoothIact;
2392 siglim2 = dsqr(siglim);
2393 /* pre-factor of Gaussian */
2394 gaufact = 1.0/(sqrt(2*M_PI)*opt->sigSmoothIact);
2395 invtwosig2 = 0.5/dsqr(opt->sigSmoothIact);
2397 for (i = 0; i < nwins; i++)
2399 snew(window[i].tausmooth, window[i].nPull);
2400 for (ig = 0; ig < window[i].nPull; ig++)
2404 pos = window[i].pos[ig];
2405 for (j = 0; j < nwins; j++)
2407 for (jg = 0; jg < window[j].nPull; jg++)
2409 dpos2 = dsqr(window[j].pos[jg]-pos);
2410 if (dpos2 < siglim2)
2412 w = gaufact*exp(-dpos2*invtwosig2);
2414 tausm += w*window[j].tau[jg];
2415 /*printf("Weight %g dpos2=%g pos=%g gaufact=%g invtwosig2=%g\n",
2416 w,dpos2,pos,gaufact,invtwosig2); */
2421 if (opt->bAllowReduceIact || tausm > window[i].tau[ig])
2423 window[i].tausmooth[ig] = tausm;
2427 window[i].tausmooth[ig] = window[i].tau[ig];
2429 window[i].g[ig] = 1+2*tausm/window[i].dt;
2434 /* try to compute the autocorrelation time for each umbrealla window */
2435 #define WHAM_AC_ZERO_LIMIT 0.05
2436 void calcIntegratedAutocorrelationTimes(t_UmbrellaWindow *window, int nwins,
2437 t_UmbrellaOptions *opt, const char *fn)
2439 int i, ig, ncorr, ntot, j, k, *count, restart;
2440 real *corr, c0, dt, timemax, tmp;
2441 real *ztime, av, tausteps;
2442 FILE *fp, *fpcorr = 0;
2446 fpcorr = xvgropen("hist_autocorr.xvg", "Autocorrelation functions of umbrella windows",
2447 "time [ps]", "autocorrelation function", opt->oenv);
2451 for (i = 0; i < nwins; i++)
2453 printf("\rEstimating integrated autocorreltion times ... [%2.0f%%] ...", 100.*(i+1)/nwins);
2455 ntot = window[i].Ntot[0];
2457 /* using half the maximum time as length of autocorrelation function */
2461 gmx_fatal(FARGS, "Tryig to estimtate autocorrelation time from only %d"
2462 " points. Provide more pull data!", ntot);
2465 /* snew(corrSq,ncorr); */
2469 snew(window[i].tau, window[i].nPull);
2470 restart = (int)(opt->acTrestart/dt+0.5);
2476 for (ig = 0; ig < window[i].nPull; ig++)
2478 if (ntot != window[i].Ntot[ig])
2480 gmx_fatal(FARGS, "Encountered different nr of frames in different pull groups.\n"
2481 "That should not happen. (%d and %d)\n", ntot, window[i].Ntot[ig]);
2483 ztime = window[i].ztime[ig];
2485 /* calc autocorrelation function C(t) = < [z(tau)-<z>]*[z(tau+t)-<z>]> */
2486 for (j = 0, av = 0; (j < ntot); j++)
2491 for (k = 0; (k < ncorr); k++)
2496 for (j = 0; (j < ntot); j += restart)
2498 for (k = 0; (k < ncorr) && (j+k < ntot); k++)
2500 tmp = (ztime[j]-av)*(ztime[j+k]-av);
2502 /* corrSq[k] += tmp*tmp; */
2506 /* divide by nr of frames for each time displacement */
2507 for (k = 0; (k < ncorr); k++)
2509 /* count probably = (ncorr-k+(restart-1))/restart; */
2510 corr[k] = corr[k]/count[k];
2511 /* variance of autocorrelation function */
2512 /* corrSq[k]=corrSq[k]/count[k]; */
2514 /* normalize such that corr[0] == 0 */
2516 for (k = 0; (k < ncorr); k++)
2519 /* corrSq[k]*=c0*c0; */
2522 /* write ACFs in verbose mode */
2525 for (k = 0; (k < ncorr); k++)
2527 fprintf(fpcorr, "%g %g\n", k*dt, corr[k]);
2529 fprintf(fpcorr, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
2532 /* esimate integrated correlation time, fitting is too unstable */
2533 tausteps = 0.5*corr[0];
2534 /* consider corr below WHAM_AC_ZERO_LIMIT as noise */
2535 for (j = 1; (j < ncorr) && (corr[j] > WHAM_AC_ZERO_LIMIT); j++)
2537 tausteps += corr[j];
2540 /* g = 1+2*tau, see. Ferrenberg/Swendsen, PRL 63:1195 (1989) or
2541 Kumar et al, eq. 28 ff. */
2542 window[i].tau[ig] = tausteps*dt;
2543 window[i].g[ig] = 1+2*tausteps;
2544 /* printf("win %d, group %d, estimated correlation time = %g ps\n",i,ig,window[i].tau[ig]); */
2555 /* plot IACT along reaction coordinate */
2556 fp = xvgropen(fn, "Integrated autocorrelation times", "z", "IACT [ps]", opt->oenv);
2557 if(output_env_get_print_xvgr_codes(opt->oenv))
2559 fprintf(fp, "@ s0 symbol 1\n@ s0 symbol size 0.5\n@ s0 line linestyle 0\n");
2560 fprintf(fp, "# WIN tau(gr1) tau(gr2) ...\n");
2561 for (i = 0; i < nwins; i++)
2563 fprintf(fp, "# %3d ", i);
2564 for (ig = 0; ig < window[i].nPull; ig++)
2566 fprintf(fp, " %11g", window[i].tau[ig]);
2571 for (i = 0; i < nwins; i++)
2573 for (ig = 0; ig < window[i].nPull; ig++)
2575 fprintf(fp, "%8g %8g\n", window[i].pos[ig], window[i].tau[ig]);
2578 if (opt->sigSmoothIact > 0.0)
2580 printf("Smoothing autocorrelation times along reaction coordinate with Gaussian of sig = %g\n",
2581 opt->sigSmoothIact);
2582 /* smooth IACT along reaction coordinate and overwrite g=1+2tau */
2583 smoothIact(window, nwins, opt);
2584 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(opt->oenv) ? "&" : "");
2585 if(output_env_get_print_xvgr_codes(opt->oenv))
2587 fprintf(fp, "@ s1 symbol 1\n@ s1 symbol size 0.5\n@ s1 line linestyle 0\n");
2588 fprintf(fp, "@ s1 symbol color 2\n");
2590 for (i = 0; i < nwins; i++)
2592 for (ig = 0; ig < window[i].nPull; ig++)
2594 fprintf(fp, "%8g %8g\n", window[i].pos[ig], window[i].tausmooth[ig]);
2599 printf("Wrote %s\n", fn);
2602 /* compute average and sigma of each umbrella window */
2603 void averageSigma(t_UmbrellaWindow *window, int nwins, t_UmbrellaOptions *opt)
2606 real av, sum2, sig, diff, *ztime, nSamplesIndep;
2608 for (i = 0; i < nwins; i++)
2610 snew(window[i].aver, window[i].nPull);
2611 snew(window[i].sigma, window[i].nPull);
2613 ntot = window[i].Ntot[0];
2614 for (ig = 0; ig < window[i].nPull; ig++)
2616 ztime = window[i].ztime[ig];
2617 for (k = 0, av = 0.; k < ntot; k++)
2622 for (k = 0, sum2 = 0.; k < ntot; k++)
2627 sig = sqrt(sum2/ntot);
2628 window[i].aver[ig] = av;
2630 /* Note: This estimate for sigma is biased from the limited sampling.
2631 Correct sigma by n/(n-1) where n = number of independent
2632 samples. Only possible if IACT is known.
2636 nSamplesIndep = window[i].N[ig]/(window[i].tau[ig]/window[i].dt);
2637 window[i].sigma[ig] = sig * nSamplesIndep/(nSamplesIndep-1);
2641 window[i].sigma[ig] = sig;
2643 printf("win %d, aver = %f sig = %f\n", i, av, window[i].sigma[ig]);
2649 /* Use histograms to compute average force on pull group.
2650 In addition, compute the sigma of the histogram.
2652 void computeAverageForce(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt)
2654 int i, j, bins = opt->bins, k;
2655 double dz, min = opt->min, max = opt->max, displAv, displAv2, temp, distance, ztot, ztot_half, w, weight;
2658 dz = (max-min)/bins;
2659 ztot = opt->max-min;
2662 /* Compute average displacement from histograms */
2663 for (j = 0; j < nWindows; ++j)
2665 snew(window[j].forceAv, window[j].nPull);
2666 for (k = 0; k < window[j].nPull; ++k)
2671 for (i = 0; i < opt->bins; ++i)
2673 temp = (1.0*i+0.5)*dz+min;
2674 distance = temp - window[j].pos[k];
2676 { /* in cyclic wham: */
2677 if (distance > ztot_half) /* |distance| < ztot_half */
2681 else if (distance < -ztot_half)
2686 w = window[j].Histo[k][i]/window[j].g[k];
2687 displAv += w*distance;
2688 displAv2 += w*sqr(distance);
2690 /* Are we near min or max? We are getting wron forces from the histgrams since
2691 the histigrams are zero outside [min,max). Therefore, assume that the position
2692 on the other side of the histomgram center is equally likely. */
2695 posmirrored = window[j].pos[k]-distance;
2696 if (posmirrored >= max || posmirrored < min)
2698 displAv += -w*distance;
2699 displAv2 += w*sqr(-distance);
2707 /* average force from average displacement */
2708 window[j].forceAv[k] = displAv*window[j].k[k];
2709 /* sigma from average square displacement */
2710 /* window[j].sigma [k] = sqrt(displAv2); */
2711 /* printf("Win %d, sigma = %f\n",j,sqrt(displAv2)); */
2716 /* Check if the complete reaction coordinate is covered by the histograms */
2717 void checkReactionCoordinateCovered(t_UmbrellaWindow *window, int nwins,
2718 t_UmbrellaOptions *opt)
2720 int i, ig, j, bins = opt->bins, bBoundary;
2721 real avcount = 0, z, relcount, *count;
2722 snew(count, opt->bins);
2724 for (j = 0; j < opt->bins; ++j)
2726 for (i = 0; i < nwins; i++)
2728 for (ig = 0; ig < window[i].nPull; ig++)
2730 count[j] += window[i].Histo[ig][j];
2733 avcount += 1.0*count[j];
2736 for (j = 0; j < bins; ++j)
2738 relcount = count[j]/avcount;
2739 z = (j+0.5)*opt->dz+opt->min;
2740 bBoundary = ( j<bins/20 || (bins-j)>bins/20 );
2741 /* check for bins with no data */
2744 fprintf(stderr, "\nWARNING, no data point in bin %d (z=%g) !\n"
2745 "You may not get a reasonable profile. Check your histograms!\n", j, z);
2747 /* and check for poor sampling */
2748 else if (relcount < 0.005 && !bBoundary)
2750 fprintf(stderr, "Warning, poor sampling bin %d (z=%g). Check your histograms!\n", j, z);
2757 void guessPotByIntegration(t_UmbrellaWindow *window, int nWindows, t_UmbrellaOptions *opt,
2760 int i, j, ig, bins = opt->bins, nHist, winmin, groupmin;
2761 double dz, min = opt->min, *pot, pos, hispos, dist, diff, fAv, distmin, *f;
2764 dz = (opt->max-min)/bins;
2766 printf("Getting initial potential by integration.\n");
2768 /* Compute average displacement from histograms */
2769 computeAverageForce(window, nWindows, opt);
2771 /* Get force for each bin from all histograms in this bin, or, alternatively,
2772 if no histograms are inside this bin, from the closest histogram */
2775 for (j = 0; j < opt->bins; ++j)
2777 pos = (1.0*j+0.5)*dz+min;
2781 groupmin = winmin = 0;
2782 for (i = 0; i < nWindows; i++)
2784 for (ig = 0; ig < window[i].nPull; ig++)
2786 hispos = window[i].pos[ig];
2787 dist = fabs(hispos-pos);
2788 /* average force within bin */
2792 fAv += window[i].forceAv[ig];
2794 /* at the same time, rememer closest histogram */
2803 /* if no histogram found in this bin, use closest histogram */
2810 fAv = window[winmin].forceAv[groupmin];
2814 for (j = 1; j < opt->bins; ++j)
2816 pot[j] = pot[j-1] - 0.5*dz*(f[j-1]+f[j]);
2819 /* cyclic wham: linearly correct possible offset */
2822 diff = (pot[bins-1]-pot[0])/(bins-1);
2823 for (j = 1; j < opt->bins; ++j)
2830 fp = xvgropen("pmfintegrated.xvg", "PMF from force integration", "z", "PMF [kJ/mol]", opt->oenv);
2831 for (j = 0; j < opt->bins; ++j)
2833 fprintf(fp, "%g %g\n", (j+0.5)*dz+opt->min, pot[j]);
2836 printf("verbose mode: wrote %s with PMF from interated forces\n", "pmfintegrated.xvg");
2839 /* get initial z=exp(-F[i]/kT) from integrated potential, where F[i] denote the free
2840 energy offsets which are usually determined by wham
2841 First: turn pot into probabilities:
2843 for (j = 0; j < opt->bins; ++j)
2845 pot[j] = exp(-pot[j]/(8.314e-3*opt->Temperature));
2847 calc_z(pot, window, nWindows, opt, TRUE);
2854 int gmx_wham(int argc, char *argv[])
2856 const char *desc[] = {
2857 "This is an analysis program that implements the Weighted",
2858 "Histogram Analysis Method (WHAM). It is intended to analyze",
2859 "output files generated by umbrella sampling simulations to ",
2860 "compute a potential of mean force (PMF). [PAR] ",
2861 "At present, three input modes are supported.[BR]",
2862 "[TT]*[tt] With option [TT]-it[tt], the user provides a file which contains the",
2863 " file names of the umbrella simulation run-input files ([TT].tpr[tt] files),",
2864 " AND, with option [TT]-ix[tt], a file which contains file names of",
2865 " the pullx [TT]mdrun[tt] output files. The [TT].tpr[tt] and pullx files must",
2866 " be in corresponding order, i.e. the first [TT].tpr[tt] created the",
2867 " first pullx, etc.[BR]",
2868 "[TT]*[tt] Same as the previous input mode, except that the the user",
2869 " provides the pull force output file names ([TT]pullf.xvg[tt]) with option [TT]-if[tt].",
2870 " From the pull force the position in the umbrella potential is",
2871 " computed. This does not work with tabulated umbrella potentials.[BR]"
2872 "[TT]*[tt] With option [TT]-ip[tt], the user provides file names of (gzipped) [TT].pdo[tt] files, i.e.",
2873 " the GROMACS 3.3 umbrella output files. If you have some unusual"
2874 " reaction coordinate you may also generate your own [TT].pdo[tt] files and",
2875 " feed them with the [TT]-ip[tt] option into to [TT]g_wham[tt]. The [TT].pdo[tt] file header",
2876 " must be similar to the following:[PAR]",
2877 "[TT]# UMBRELLA 3.0[BR]",
2878 "# Component selection: 0 0 1[BR]",
2880 "# Ref. Group 'TestAtom'[BR]",
2881 "# Nr. of pull groups 2[BR]",
2882 "# Group 1 'GR1' Umb. Pos. 5.0 Umb. Cons. 1000.0[BR]",
2883 "# Group 2 'GR2' Umb. Pos. 2.0 Umb. Cons. 500.0[BR]",
2885 "The number of pull groups, umbrella positions, force constants, and names ",
2886 "may (of course) differ. Following the header, a time column and ",
2887 "a data column for each pull group follows (i.e. the displacement",
2888 "with respect to the umbrella center). Up to four pull groups are possible ",
2889 "per [TT].pdo[tt] file at present.[PAR]",
2890 "By default, the output files are[BR]",
2891 " [TT]-o[tt] PMF output file[BR]",
2892 " [TT]-hist[tt] Histograms output file[BR]",
2893 "Always check whether the histograms sufficiently overlap.[PAR]",
2894 "The umbrella potential is assumed to be harmonic and the force constants are ",
2895 "read from the [TT].tpr[tt] or [TT].pdo[tt] files. If a non-harmonic umbrella force was applied ",
2896 "a tabulated potential can be provided with [TT]-tab[tt].[PAR]",
2897 "WHAM OPTIONS[BR]------------[BR]",
2898 " [TT]-bins[tt] Number of bins used in analysis[BR]",
2899 " [TT]-temp[tt] Temperature in the simulations[BR]",
2900 " [TT]-tol[tt] Stop iteration if profile (probability) changed less than tolerance[BR]",
2901 " [TT]-auto[tt] Automatic determination of boundaries[BR]",
2902 " [TT]-min,-max[tt] Boundaries of the profile [BR]",
2903 "The data points that are used to compute the profile",
2904 "can be restricted with options [TT]-b[tt], [TT]-e[tt], and [TT]-dt[tt]. ",
2905 "Adjust [TT]-b[tt] to ensure sufficient equilibration in each ",
2906 "umbrella window.[PAR]",
2907 "With [TT]-log[tt] (default) the profile is written in energy units, otherwise ",
2908 "(with [TT]-nolog[tt]) as probability. The unit can be specified with [TT]-unit[tt]. ",
2909 "With energy output, the energy in the first bin is defined to be zero. ",
2910 "If you want the free energy at a different ",
2911 "position to be zero, set [TT]-zprof0[tt] (useful with bootstrapping, see below).[PAR]",
2912 "For cyclic or periodic reaction coordinates (dihedral angle, channel PMF",
2913 "without osmotic gradient), the option [TT]-cycl[tt] is useful. [TT]g_wham[tt] will make use of the ",
2914 "periodicity of the system and generate a periodic PMF. The first and the last bin of the",
2915 "reaction coordinate will assumed be be neighbors.[PAR]",
2916 "Option [TT]-sym[tt] symmetrizes the profile around z=0 before output, ",
2917 "which may be useful for, e.g. membranes.[PAR]",
2918 "AUTOCORRELATIONS[BR]----------------[BR]",
2919 "With [TT]-ac[tt], [TT]g_wham[tt] estimates the integrated autocorrelation ",
2920 "time (IACT) [GRK]tau[grk] for each umbrella window and weights the respective ",
2921 "window with 1/[1+2*[GRK]tau[grk]/dt]. The IACTs are written ",
2922 "to the file defined with [TT]-oiact[tt]. In verbose mode, all ",
2923 "autocorrelation functions (ACFs) are written to [TT]hist_autocorr.xvg[tt]. ",
2924 "Because the IACTs can be severely underestimated in case of limited ",
2925 "sampling, option [TT]-acsig[tt] allows one to smooth the IACTs along the ",
2926 "reaction coordinate with a Gaussian ([GRK]sigma[grk] provided with [TT]-acsig[tt], ",
2927 "see output in [TT]iact.xvg[tt]). Note that the IACTs are estimated by simple ",
2928 "integration of the ACFs while the ACFs are larger 0.05.",
2929 "If you prefer to compute the IACTs by a more sophisticated (but possibly ",
2930 "less robust) method such as fitting to a double exponential, you can ",
2931 "compute the IACTs with [TT]g_analyze[tt] and provide them to [TT]g_wham[tt] with the file ",
2932 "[TT]iact-in.dat[tt] (option [TT]-iiact[tt]), which should contain one line per ",
2933 "input file ([TT].pdo[tt] or pullx/f file) and one column per pull group in the respective file.[PAR]",
2934 "ERROR ANALYSIS[BR]--------------[BR]",
2935 "Statistical errors may be estimated with bootstrap analysis. Use it with care, ",
2936 "otherwise the statistical error may be substantially underestimated. ",
2937 "More background and examples for the bootstrap technique can be found in ",
2938 "Hub, de Groot and Van der Spoel, JCTC (2010) 6: 3713-3720.[BR]",
2939 "[TT]-nBootstrap[tt] defines the number of bootstraps (use, e.g., 100). ",
2940 "Four bootstrapping methods are supported and ",
2941 "selected with [TT]-bs-method[tt].[BR]",
2942 " (1) [TT]b-hist[tt] Default: complete histograms are considered as independent ",
2943 "data points, and the bootstrap is carried out by assigning random weights to the ",
2944 "histograms (\"Bayesian bootstrap\"). Note that each point along the reaction coordinate",
2945 "must be covered by multiple independent histograms (e.g. 10 histograms), otherwise the ",
2946 "statistical error is underestimated.[BR]",
2947 " (2) [TT]hist[tt] Complete histograms are considered as independent data points. ",
2948 "For each bootstrap, N histograms are randomly chosen from the N given histograms ",
2949 "(allowing duplication, i.e. sampling with replacement).",
2950 "To avoid gaps without data along the reaction coordinate blocks of histograms ",
2951 "([TT]-histbs-block[tt]) may be defined. In that case, the given histograms are ",
2952 "divided into blocks and only histograms within each block are mixed. Note that ",
2953 "the histograms within each block must be representative for all possible histograms, ",
2954 "otherwise the statistical error is underestimated.[BR]",
2955 " (3) [TT]traj[tt] The given histograms are used to generate new random trajectories,",
2956 "such that the generated data points are distributed according the given histograms ",
2957 "and properly autocorrelated. The autocorrelation time (ACT) for each window must be ",
2958 "known, so use [TT]-ac[tt] or provide the ACT with [TT]-iiact[tt]. If the ACT of all ",
2959 "windows are identical (and known), you can also provide them with [TT]-bs-tau[tt]. ",
2960 "Note that this method may severely underestimate the error in case of limited sampling, ",
2961 "that is if individual histograms do not represent the complete phase space at ",
2962 "the respective positions.[BR]",
2963 " (4) [TT]traj-gauss[tt] The same as method [TT]traj[tt], but the trajectories are ",
2964 "not bootstrapped from the umbrella histograms but from Gaussians with the average ",
2965 "and width of the umbrella histograms. That method yields similar error estimates ",
2966 "like method [TT]traj[tt].[PAR]"
2967 "Bootstrapping output:[BR]",
2968 " [TT]-bsres[tt] Average profile and standard deviations[BR]",
2969 " [TT]-bsprof[tt] All bootstrapping profiles[BR]",
2970 "With [TT]-vbs[tt] (verbose bootstrapping), the histograms of each bootstrap are written, ",
2971 "and, with bootstrap method [TT]traj[tt], the cumulative distribution functions of ",
2975 const char *en_unit[] = {NULL, "kJ", "kCal", "kT", NULL};
2976 const char *en_unit_label[] = {"", "E (kJ mol\\S-1\\N)", "E (kcal mol\\S-1\\N)", "E (kT)", NULL};
2977 const char *en_bsMethod[] = { NULL, "b-hist", "hist", "traj", "traj-gauss", NULL };
2979 static t_UmbrellaOptions opt;
2982 { "-min", FALSE, etREAL, {&opt.min},
2983 "Minimum coordinate in profile"},
2984 { "-max", FALSE, etREAL, {&opt.max},
2985 "Maximum coordinate in profile"},
2986 { "-auto", FALSE, etBOOL, {&opt.bAuto},
2987 "Determine min and max automatically"},
2988 { "-bins", FALSE, etINT, {&opt.bins},
2989 "Number of bins in profile"},
2990 { "-temp", FALSE, etREAL, {&opt.Temperature},
2992 { "-tol", FALSE, etREAL, {&opt.Tolerance},
2994 { "-v", FALSE, etBOOL, {&opt.verbose},
2996 { "-b", FALSE, etREAL, {&opt.tmin},
2997 "First time to analyse (ps)"},
2998 { "-e", FALSE, etREAL, {&opt.tmax},
2999 "Last time to analyse (ps)"},
3000 { "-dt", FALSE, etREAL, {&opt.dt},
3001 "Analyse only every dt ps"},
3002 { "-histonly", FALSE, etBOOL, {&opt.bHistOnly},
3003 "Write histograms and exit"},
3004 { "-boundsonly", FALSE, etBOOL, {&opt.bBoundsOnly},
3005 "Determine min and max and exit (with [TT]-auto[tt])"},
3006 { "-log", FALSE, etBOOL, {&opt.bLog},
3007 "Calculate the log of the profile before printing"},
3008 { "-unit", FALSE, etENUM, {en_unit},
3009 "Energy unit in case of log output" },
3010 { "-zprof0", FALSE, etREAL, {&opt.zProf0},
3011 "Define profile to 0.0 at this position (with [TT]-log[tt])"},
3012 { "-cycl", FALSE, etBOOL, {&opt.bCycl},
3013 "Create cyclic/periodic profile. Assumes min and max are the same point."},
3014 { "-sym", FALSE, etBOOL, {&opt.bSym},
3015 "Symmetrize profile around z=0"},
3016 { "-hist-eq", FALSE, etBOOL, {&opt.bHistEq},
3017 "HIDDENEnforce equal weight for all histograms. (Non-Weighed-HAM)"},
3018 { "-ac", FALSE, etBOOL, {&opt.bCalcTauInt},
3019 "Calculate integrated autocorrelation times and use in wham"},
3020 { "-acsig", FALSE, etREAL, {&opt.sigSmoothIact},
3021 "Smooth autocorrelation times along reaction coordinate with Gaussian of this [GRK]sigma[grk]"},
3022 { "-ac-trestart", FALSE, etREAL, {&opt.acTrestart},
3023 "When computing autocorrelation functions, restart computing every .. (ps)"},
3024 { "-acred", FALSE, etBOOL, {&opt.bAllowReduceIact},
3025 "HIDDENWhen smoothing the ACTs, allow to reduce ACTs. Otherwise, only increase ACTs "
3026 "during smoothing"},
3027 { "-nBootstrap", FALSE, etINT, {&opt.nBootStrap},
3028 "nr of bootstraps to estimate statistical uncertainty (e.g., 200)" },
3029 { "-bs-method", FALSE, etENUM, {en_bsMethod},
3030 "Bootstrap method" },
3031 { "-bs-tau", FALSE, etREAL, {&opt.tauBootStrap},
3032 "Autocorrelation time (ACT) assumed for all histograms. Use option [TT]-ac[tt] if ACT is unknown."},
3033 { "-bs-seed", FALSE, etINT, {&opt.bsSeed},
3034 "Seed for bootstrapping. (-1 = use time)"},
3035 { "-histbs-block", FALSE, etINT, {&opt.histBootStrapBlockLength},
3036 "When mixing histograms only mix within blocks of [TT]-histbs-block[tt]."},
3037 { "-vbs", FALSE, etBOOL, {&opt.bs_verbose},
3038 "Verbose bootstrapping. Print the CDFs and a histogram file for each bootstrap."},
3039 { "-stepout", FALSE, etINT, {&opt.stepchange},
3040 "HIDDENWrite maximum change every ... (set to 1 with [TT]-v[tt])"},
3041 { "-updateContr", FALSE, etINT, {&opt.stepUpdateContrib},
3042 "HIDDENUpdate table with significan contributions to WHAM every ... iterations"},
3046 { efDAT, "-ix", "pullx-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs */
3047 { efDAT, "-if", "pullf-files", ffOPTRD}, /* wham input: pullf.xvg's and tprs */
3048 { efDAT, "-it", "tpr-files", ffOPTRD}, /* wham input: tprs */
3049 { efDAT, "-ip", "pdo-files", ffOPTRD}, /* wham input: pdo files (gmx3 style) */
3050 { efXVG, "-o", "profile", ffWRITE }, /* output file for profile */
3051 { efXVG, "-hist", "histo", ffWRITE}, /* output file for histograms */
3052 { efXVG, "-oiact", "iact", ffOPTWR}, /* writing integrated autocorrelation times */
3053 { efDAT, "-iiact", "iact-in", ffOPTRD}, /* reading integrated autocorrelation times */
3054 { efXVG, "-bsres", "bsResult", ffOPTWR}, /* average and errors of bootstrap analysis */
3055 { efXVG, "-bsprof", "bsProfs", ffOPTWR}, /* output file for bootstrap profiles */
3056 { efDAT, "-tab", "umb-pot", ffOPTRD}, /* Tabulated umbrella potential (if not harmonic) */
3059 int i, j, l, nfiles, nwins, nfiles2;
3060 t_UmbrellaHeader header;
3061 t_UmbrellaWindow * window = NULL;
3062 double *profile, maxchange = 1e20;
3063 gmx_bool bMinSet, bMaxSet, bAutoSet, bExact = FALSE;
3064 char **fninTpr, **fninPull, **fninPdo;
3066 FILE *histout, *profout;
3067 char ylabel[256], title[256];
3069 #define NFILE asize(fnm)
3071 CopyRight(stderr, argv[0]);
3074 opt.verbose = FALSE;
3075 opt.bHistOnly = FALSE;
3084 /* bootstrapping stuff */
3086 opt.bsMethod = bsMethod_hist;
3087 opt.tauBootStrap = 0.0;
3089 opt.histBootStrapBlockLength = 8;
3090 opt.bs_verbose = FALSE;
3095 opt.Temperature = 298;
3096 opt.Tolerance = 1e-6;
3097 opt.bBoundsOnly = FALSE;
3099 opt.bCalcTauInt = FALSE;
3100 opt.sigSmoothIact = 0.0;
3101 opt.bAllowReduceIact = TRUE;
3102 opt.bInitPotByIntegration = TRUE;
3103 opt.acTrestart = 1.0;
3104 opt.stepchange = 100;
3105 opt.stepUpdateContrib = 100;
3107 parse_common_args(&argc, argv, PCA_BE_NICE,
3108 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &opt.oenv);
3110 opt.unit = nenum(en_unit);
3111 opt.bsMethod = nenum(en_bsMethod);
3113 opt.bProf0Set = opt2parg_bSet("-zprof0", asize(pa), pa);
3115 opt.bTab = opt2bSet("-tab", NFILE, fnm);
3116 opt.bPdo = opt2bSet("-ip", NFILE, fnm);
3117 opt.bTpr = opt2bSet("-it", NFILE, fnm);
3118 opt.bPullx = opt2bSet("-ix", NFILE, fnm);
3119 opt.bPullf = opt2bSet("-if", NFILE, fnm);
3120 opt.bTauIntGiven = opt2bSet("-iiact", NFILE, fnm);
3121 if (opt.bTab && opt.bPullf)
3123 gmx_fatal(FARGS, "Force input does not work with tabulated potentials. "
3124 "Provide pullx.xvg or pdo files!");
3127 #define WHAMBOOLXOR(a, b) ( ((!(a)) && (b)) || ((a) && (!(b))))
3128 if (!opt.bPdo && !WHAMBOOLXOR(opt.bPullx, opt.bPullf))
3130 gmx_fatal(FARGS, "Give either pullx (-ix) OR pullf (-if) data. Not both.");
3132 if (!opt.bPdo && !(opt.bTpr || opt.bPullf || opt.bPullx))
3134 gmx_fatal(FARGS, "g_wham supports three input modes, pullx, pullf, or pdo file input."
3135 "\n\n Check g_wham -h !");
3138 opt.fnPdo = opt2fn("-ip", NFILE, fnm);
3139 opt.fnTpr = opt2fn("-it", NFILE, fnm);
3140 opt.fnPullf = opt2fn("-if", NFILE, fnm);
3141 opt.fnPullx = opt2fn("-ix", NFILE, fnm);
3143 bMinSet = opt2parg_bSet("-min", asize(pa), pa);
3144 bMaxSet = opt2parg_bSet("-max", asize(pa), pa);
3145 bAutoSet = opt2parg_bSet("-auto", asize(pa), pa);
3146 if ( (bMinSet || bMaxSet) && bAutoSet)
3148 gmx_fatal(FARGS, "With -auto, do not give -min or -max\n");
3151 if ( (bMinSet && !bMaxSet) || (!bMinSet && bMaxSet))
3153 gmx_fatal(FARGS, "When giving -min, you must give -max (and vice versa), too\n");
3156 if (bMinSet && opt.bAuto)
3158 printf("Note: min and max given, switching off -auto.\n");
3162 if (opt.bTauIntGiven && opt.bCalcTauInt)
3164 gmx_fatal(FARGS, "Either read (option -iiact) or calculate (option -ac) the\n"
3165 "the autocorrelation times. Not both.");
3168 if (opt.tauBootStrap > 0.0 && opt2parg_bSet("-ac", asize(pa), pa))
3170 gmx_fatal(FARGS, "Either compute autocorrelation times (ACTs) (option -ac) or "
3171 "provide it with -bs-tau for bootstrapping. Not Both.\n");
3173 if (opt.tauBootStrap > 0.0 && opt2bSet("-iiact", NFILE, fnm))
3175 gmx_fatal(FARGS, "Either provide autocorrelation times (ACTs) with file iact-in.dat "
3176 "(option -iiact) or define all ACTs with -bs-tau for bootstrapping\n. Not Both.");
3180 /* Reading gmx4 pull output and tpr files */
3181 if (opt.bTpr || opt.bPullf || opt.bPullx)
3183 read_wham_in(opt.fnTpr, &fninTpr, &nfiles, &opt);
3185 fnPull = opt.bPullf ? opt.fnPullf : opt.fnPullx;
3186 read_wham_in(fnPull, &fninPull, &nfiles2, &opt);
3187 printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n",
3188 nfiles, nfiles2, opt.bPullf ? "force" : "position", opt.fnTpr, fnPull);
3189 if (nfiles != nfiles2)
3191 gmx_fatal(FARGS, "Found %d file names in %s, but %d in %s\n", nfiles,
3192 opt.fnTpr, nfiles2, fnPull);
3194 window = initUmbrellaWindows(nfiles);
3195 read_tpr_pullxf_files(fninTpr, fninPull, nfiles, &header, window, &opt);
3198 { /* reading pdo files */
3199 read_wham_in(opt.fnPdo, &fninPdo, &nfiles, &opt);
3200 printf("Found %d pdo files in %s\n", nfiles, opt.fnPdo);
3201 window = initUmbrellaWindows(nfiles);
3202 read_pdo_files(fninPdo, nfiles, &header, window, &opt);
3206 /* enforce equal weight for all histograms? */
3209 enforceEqualWeights(window, nwins);
3212 /* write histograms */
3213 histout = xvgropen(opt2fn("-hist", NFILE, fnm), "Umbrella histograms",
3214 "z", "count", opt.oenv);
3215 for (l = 0; l < opt.bins; ++l)
3217 fprintf(histout, "%e\t", (double)(l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
3218 for (i = 0; i < nwins; ++i)
3220 for (j = 0; j < window[i].nPull; ++j)
3222 fprintf(histout, "%e\t", window[i].Histo[j][l]);
3225 fprintf(histout, "\n");
3228 printf("Wrote %s\n", opt2fn("-hist", NFILE, fnm));
3231 printf("Wrote histograms to %s, now exiting.\n", opt2fn("-hist", NFILE, fnm));
3235 /* Using tabulated umbrella potential */
3238 setup_tab(opt2fn("-tab", NFILE, fnm), &opt);
3241 /* Integrated autocorrelation times provided ? */
3242 if (opt.bTauIntGiven)
3244 readIntegratedAutocorrelationTimes(window, nwins, &opt, opt2fn("-iiact", NFILE, fnm));
3247 /* Compute integrated autocorrelation times */
3248 if (opt.bCalcTauInt)
3250 calcIntegratedAutocorrelationTimes(window, nwins, &opt, opt2fn("-oiact", NFILE, fnm));
3253 /* calc average and sigma for each histogram
3254 (maybe required for bootstrapping. If not, this is fast anyhow) */
3255 if (opt.nBootStrap && opt.bsMethod == bsMethod_trajGauss)
3257 averageSigma(window, nwins, &opt);
3260 /* Get initial potential by simple integration */
3261 if (opt.bInitPotByIntegration)
3263 guessPotByIntegration(window, nwins, &opt, 0);
3266 /* Check if complete reaction coordinate is covered */
3267 checkReactionCoordinateCovered(window, nwins, &opt);
3269 /* Calculate profile */
3270 snew(profile, opt.bins);
3278 if ( (i%opt.stepUpdateContrib) == 0)
3280 setup_acc_wham(profile, window, nwins, &opt);
3282 if (maxchange < opt.Tolerance)
3285 /* if (opt.verbose) */
3286 printf("Switched to exact iteration in iteration %d\n", i);
3288 calc_profile(profile, window, nwins, &opt, bExact);
3289 if (((i%opt.stepchange) == 0 || i == 1) && !i == 0)
3291 printf("\t%4d) Maximum change %e\n", i, maxchange);
3295 while ( (maxchange = calc_z(profile, window, nwins, &opt, bExact)) > opt.Tolerance || !bExact);
3296 printf("Converged in %d iterations. Final maximum change %g\n", i, maxchange);
3298 /* calc error from Kumar's formula */
3299 /* Unclear how the error propagates along reaction coordinate, therefore
3301 /* calc_error_kumar(profile,window, nwins,&opt); */
3303 /* Write profile in energy units? */
3306 prof_normalization_and_unit(profile, &opt);
3307 strcpy(ylabel, en_unit_label[opt.unit]);
3308 strcpy(title, "Umbrella potential");
3312 strcpy(ylabel, "Density of states");
3313 strcpy(title, "Density of states");
3316 /* symmetrize profile around z=0? */
3319 symmetrizeProfile(profile, &opt);
3322 /* write profile or density of states */
3323 profout = xvgropen(opt2fn("-o", NFILE, fnm), title, "z", ylabel, opt.oenv);
3324 for (i = 0; i < opt.bins; ++i)
3326 fprintf(profout, "%e\t%e\n", (double)(i+0.5)/opt.bins*(opt.max-opt.min)+opt.min, profile[i]);
3329 printf("Wrote %s\n", opt2fn("-o", NFILE, fnm));
3331 /* Bootstrap Method */
3334 do_bootstrapping(opt2fn("-bsres", NFILE, fnm), opt2fn("-bsprof", NFILE, fnm),
3335 opt2fn("-hist", NFILE, fnm),
3336 ylabel, profile, window, nwins, &opt);
3340 freeUmbrellaWindows(window, nfiles);
3342 printf("\nIn case you use results from g_wham for a publication, please cite:\n");
3343 please_cite(stdout, "Hub2010");