1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*- */
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
50 #include "gmx_random.h"
60 #define WHAM_MAXFILELEN 2048
62 /* enum for energy units */
63 enum { enSel, en_kJ, en_kCal, en_kT, enNr };
64 /* enum for type of input files (pdos, tpr, or pullf) */
65 enum { whamin_unknown, whamin_tpr, whamin_pullxf, whamin_pdo };
66 /* enum for bootstrapping method (
67 - bootstrap complete histograms with continuous weights (Bayesian bootstrap)
68 - bootstrap complete histograms
69 - bootstrap trajectories from given umbrella histograms
70 - bootstrap trajectories from Gaussian with mu/sigam computed from
71 the respective histogram
73 ********************************************************************
74 FOR MORE DETAILS ON THE BOOTSTRAP METHODS (INCLUDING EXAMPLES), SEE
75 JS Hub, BL de Groot, D van der Spoel
76 [TT]g_wham[tt] - A free weighted histogram analysis implementation including
77 robust error and autocorrelation estimates,
78 J Chem Theory Comput, accepted (2010)
79 ********************************************************************
81 enum { bsMethod_unknown, bsMethod_BayesianHist, bsMethod_hist,
82 bsMethod_traj, bsMethod_trajGauss };
87 /* umbrella with pull code of gromacs 4.x */
88 int npullgrps; /* nr of pull groups in tpr file */
89 int pull_geometry; /* such as distance, position */
90 ivec pull_dim; /* pull dimension with geometry distance */
91 int pull_ndim; /* nr of pull_dim != 0 */
92 real *k; /* force constants in tpr file */
93 rvec *init_dist; /* reference displacements */
94 real *umbInitDist; /* reference displacement in umbrella direction */
96 /* From here, old pdo stuff */
102 char PullName[4][256];
104 double UmbCons[4][3];
109 int nPull; /* nr of pull groups in this pdo or pullf/x file */
110 double **Histo,**cum; /* nPull histograms and nPull cumulative distr. funct */
111 int nBin; /* nr of bins. identical to opt->bins */
112 double *k; /* force constants for the nPull groups */
113 double *pos; /* umbrella positions for the nPull groups */
114 double *z; /* z=(-Fi/kT) for the nPull groups. These values are
115 iteratively computed during wham */
116 double *N, *Ntot; /* nr of data points in nPull histograms. N and Ntot
117 only differ if bHistEq==TRUE */
119 double *g,*tau,*tausmooth; /* g = 1 + 2*tau[int]/dt where tau is the integrated
120 autocorrelation time. Compare, e.g.
121 Ferrenberg/Swendsen, PRL 63:1195 (1989)
122 Kumar et al, J Comp Chem 13, 1011-1021 (1992), eq. 28 */
124 double dt; /* timestep in the input data. Can be adapted with
126 gmx_bool **bContrib; /* TRUE, if any data point of the histogram is within min
127 and max, otherwise FALSE. */
128 real **ztime; /* input data z(t) as a function of time. Required to
130 real *forceAv; /* average force estimated from average displacement, fAv=dzAv*k
131 Used for integration to guess the potential. */
132 real *aver,*sigma; /* average and stddev of histograms */
133 double *bsWeight; /* for bootstrapping complete histograms with continuous weights */
140 const char *fnTpr,*fnPullf;
141 const char *fnPdo,*fnPullx; /* file names of input */
142 gmx_bool bTpr,bPullf,bPdo,bPullx;/* input file types given? */
143 real tmin, tmax, dt; /* only read input within tmin and tmax with dt */
145 gmx_bool bInitPotByIntegration; /* before WHAM, guess potential by force integration. Yields
146 1.5 to 2 times faster convergence */
147 int stepUpdateContrib; /* update contribution table every ... iterations. Accelerates
150 /* BASIC WHAM OPTIONS */
151 int bins; /* nr of bins, min, max, and dz of profile */
153 real Temperature,Tolerance; /* temperature, converged when probability changes less
155 gmx_bool bCycl; /* generate cyclic (periodic) PMF */
158 gmx_bool bLog; /* energy output (instead of probability) for profile */
159 int unit; /* unit for PMF output kJ/mol or kT or kCal/mol */
160 gmx_bool bSym; /* symmetrize PMF around z=0 after WHAM, useful for
162 real zProf0; /* after wham, set prof to zero at this z-position
163 When bootstrapping, set zProf0 to a "stable" reference
165 gmx_bool bProf0Set; /* setting profile to 0 at zProf0? */
167 gmx_bool bBoundsOnly,bHistOnly; /* determine min and max, or write histograms and exit */
168 gmx_bool bAuto; /* determine min and max automatically but do not exit */
170 gmx_bool verbose; /* more noisy wham mode */
171 int stepchange; /* print maximum change in prof after how many interations */
172 output_env_t oenv; /* xvgr options */
174 /* AUTOCORRELATION STUFF */
175 gmx_bool bTauIntGiven,bCalcTauInt;/* IACT given or should be calculated? */
176 real sigSmoothIact; /* sigma of Gaussian to smooth ACTs */
177 gmx_bool bAllowReduceIact; /* Allow to reduce ACTs during smoothing. Otherwise
178 ACT are only increased during smoothing */
179 real acTrestart; /* when computing ACT, time between restarting points */
180 gmx_bool bHistEq; /* Enforce the same weight for each umbella window, that is
181 calculate with the same number of data points for
182 each window. That can be reasonable, if the histograms
183 have different length, but due to autocorrelation,
184 a longer simulation should not have larger weightin wham. */
186 /* BOOTSTRAPPING STUFF */
187 int nBootStrap; /* nr of bootstraps (50 is usually enough) */
188 int bsMethod; /* if == bsMethod_hist, consider complete histograms as independent
189 data points and, hence, only mix complete histograms.
190 if == bsMethod_BayesianHist, consider complete histograms
191 as independent data points, but assign random weights
192 to the histograms during the bootstrapping ("Bayesian bootstrap")
193 In case of long correlations (e.g., inside a channel), these
194 will yield a more realistic error.
195 if == bsMethod_traj(Gauss), generate synthetic histograms
197 histogram by generating an autocorrelated random sequence
198 that is distributed according to the respective given
199 histogram. With bsMethod_trajGauss, bootstrap from a Gaussian
200 (instead of from the umbrella histogram) to generate a new
203 real tauBootStrap; /* autocorrelation time (ACT) used to generate synthetic
204 histograms. If ==0, use calculated ACF */
205 int histBootStrapBlockLength; /* when mixing histograms, mix only histograms withing blocks
206 long the reaction coordinate xi. Avoids gaps along xi. */
207 int bsSeed; /* random seed for bootstrapping */
208 gmx_bool bs_verbose; /* Write cumulative distribution functions (CDFs) of histograms
209 and write the generated histograms for each bootstrap */
211 /* tabulated umbrella potential stuff */
213 double *tabX,*tabY,tabMin,tabMax,tabDz;
216 gmx_rng_t rng; /* gromacs random number generator */
220 t_UmbrellaWindow * initUmbrellaWindows(int nwin)
222 t_UmbrellaWindow *win;
225 for (i=0; i<nwin; i++)
227 win[i].Histo = win[i].cum = 0;
228 win[i].k = win[i].pos = win[i].z =0;
229 win[i].N = win[i].Ntot = 0;
230 win[i].g = win[i].tau = win[i].tausmooth = 0;
234 win[i].aver = win[i].sigma = 0;
240 void freeUmbrellaWindows(t_UmbrellaWindow *win, int nwin)
243 for (i=0; i<nwin; i++)
246 for (j=0;j<win[i].nPull;j++)
247 sfree(win[i].Histo[j]);
249 for (j=0;j<win[i].nPull;j++)
250 sfree(win[i].cum[j]);
252 for (j=0;j<win[i].nPull;j++)
253 sfree(win[i].bContrib[j]);
263 sfree(win[i].tausmooth);
264 sfree(win[i].bContrib);
266 sfree(win[i].forceAv);
269 sfree(win[i].bsWeight);
274 /* Read and setup tabulated umbrella potential */
275 void setup_tab(const char *fn,t_UmbrellaOptions *opt)
280 printf("Setting up tabulated potential from file %s\n",fn);
281 nl=read_xvg(fn,&y,&ny);
284 gmx_fatal(FARGS,"Found %d columns in %s. Expected 2.\n",ny,fn);
286 opt->tabMax=y[0][nl-1];
287 opt->tabDz=(opt->tabMax-opt->tabMin)/(nl-1);
289 gmx_fatal(FARGS,"The tabulated potential in %s must be provided in \n"
290 "ascending z-direction",fn);
292 if (fabs(y[0][i+1]-y[0][i]-opt->tabDz) > opt->tabDz/1e6)
293 gmx_fatal(FARGS,"z-values in %s are not equally spaced.\n",ny,fn);
298 opt->tabX[i]=y[0][i];
299 opt->tabY[i]=y[1][i];
301 printf("Found equally spaced tabulated potential from %g to %g, spacing %g\n",
302 opt->tabMin,opt->tabMax,opt->tabDz);
305 void read_pdo_header(FILE * file,t_UmbrellaHeader * header, t_UmbrellaOptions *opt)
308 char Buffer0[256], Buffer1[256], Buffer2[256], Buffer3[256], Buffer4[256];
312 if (fgets(line,2048,file) == NULL)
313 gmx_fatal(FARGS,"Error reading header from pdo file\n");
314 sscanf(line,"%s%s%s",Buffer0,Buffer1,Buffer2);
315 if(strcmp(Buffer1,"UMBRELLA"))
316 gmx_fatal(FARGS,"This does not appear to be a valid pdo file. Found %s, expected %s\n"
317 "(Found in first line: `%s')\n",
318 Buffer1, "UMBRELLA",line);
319 if(strcmp(Buffer2,"3.0"))
320 gmx_fatal(FARGS,"This does not appear to be a version 3.0 pdo file");
323 if (fgets(line,2048,file) == NULL)
324 gmx_fatal(FARGS,"Error reading header from pdo file\n");
325 sscanf(line,"%s%s%s%d%d%d",Buffer0,Buffer1,Buffer2,
326 &(header->Dims[0]),&(header->Dims[1]),&(header->Dims[2]));
328 /* printf("%d %d %d\n", header->Dims[0],header->Dims[1],header->Dims[2]); */
330 header->nDim = header->Dims[0] + header->Dims[1] + header->Dims[2];
332 gmx_fatal(FARGS,"Currently only supports one dimension");
335 if (fgets(line,2048,file) == NULL)
336 gmx_fatal(FARGS,"Error reading header from pdo file\n");
337 sscanf(line,"%s%s%d",Buffer0,Buffer1,&(header->nSkip));
340 if (fgets(line,2048,file) == NULL)
341 gmx_fatal(FARGS,"Error reading header from pdo file\n");
342 sscanf(line,"%s%s%s%s",Buffer0,Buffer1,Buffer2,header->Reference);
345 if (fgets(line,2048,file) == NULL)
346 gmx_fatal(FARGS,"Error reading header from pdo file\n");
347 sscanf(line,"%s%s%s%s%s%d",Buffer0,Buffer1,Buffer2,Buffer3,Buffer4,&(header->nPull));
350 printf("\tFound nPull=%d , nSkip=%d, ref=%s\n",header->nPull,header->nSkip,
353 for(i=0;i<header->nPull;++i)
355 if (fgets(line,2048,file) == NULL)
356 gmx_fatal(FARGS,"Error reading header from pdo file\n");
357 sscanf(line,"%*s%*s%*s%s%*s%*s%lf%*s%*s%lf",header->PullName[i]
358 ,&(header->UmbPos[i][0]),&(header->UmbCons[i][0]));
360 printf("\tpullgroup %d, pullname = %s, UmbPos = %g, UmbConst = %g\n",
361 i,header->PullName[i],header->UmbPos[i][0],header->UmbCons[i][0]);
364 if (fgets(line,2048,file) == NULL)
365 gmx_fatal(FARGS,"Cannot read from file\n");
366 sscanf(line,"%s",Buffer3);
367 if (strcmp(Buffer3,"#####") != 0)
368 gmx_fatal(FARGS,"Expected '#####', found %s. Hick.\n",Buffer3);
372 static char *fgets3(FILE *fp,char ptr[],int *len)
377 if (fgets(ptr,*len-1,fp) == NULL)
380 while ((strchr(ptr,'\n') == NULL) && (!feof(fp)))
382 /* This line is longer than len characters, let's increase len! */
386 if (fgets(p-1,STRLEN,fp) == NULL)
390 if (ptr[slen-1] == '\n')
396 void read_pdo_data(FILE * file, t_UmbrellaHeader * header,
397 int fileno, t_UmbrellaWindow * win,
398 t_UmbrellaOptions *opt,
399 gmx_bool bGetMinMax,real *mintmp,real *maxtmp)
401 int i,inttemp,bins,count,ntot;
402 real min,max,minfound=1e20,maxfound=-1e20;
403 double temp,time,time0=0,dt;
405 t_UmbrellaWindow * window=0;
406 gmx_bool timeok,dt_ok=1;
407 char *tmpbuf=0,fmt[256],fmtign[256];
408 int len=STRLEN,dstep=1;
409 const int blocklen=4096;
419 /* Need to alocate memory and set up structure */
420 window->nPull=header->nPull;
423 snew(window->Histo,window->nPull);
424 snew(window->z,window->nPull);
425 snew(window->k,window->nPull);
426 snew(window->pos,window->nPull);
427 snew(window->N, window->nPull);
428 snew(window->Ntot, window->nPull);
429 snew(window->g, window->nPull);
430 snew(window->bsWeight, window->nPull);
434 if (opt->bCalcTauInt)
435 snew(window->ztime, window->nPull);
438 snew(lennow,window->nPull);
440 for(i=0;i<window->nPull;++i)
443 window->bsWeight[i]=1.;
444 snew(window->Histo[i],bins);
445 window->k[i]=header->UmbCons[i][0];
446 window->pos[i]=header->UmbPos[i][0];
450 if (opt->bCalcTauInt)
454 /* Done with setup */
460 min=max=bins=0; /* Get rid of warnings */
465 while ( (ptr=fgets3(file,tmpbuf,&len)) != NULL)
469 if (ptr[0] == '#' || strlen(ptr)<2)
472 /* Initiate format string */
474 strcat(fmtign,"%*s");
476 sscanf(ptr,"%lf",&time); /* printf("Time %f\n",time); */
477 /* Round time to fs */
478 time=1.0/1000*( (int) (time*1000+0.5) );
480 /* get time step of pdo file */
488 dstep=(int)(opt->dt/dt+0.5);
497 dt_ok=((count-1)%dstep == 0);
498 timeok=(dt_ok && time >= opt->tmin && time <= opt->tmax);
500 printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
501 time,opt->tmin, opt->tmax, dt_ok,timeok); */
505 for(i=0;i<header->nPull;++i)
508 strcat(fmt,"%lf"); /* Creating a format stings such as "%*s...%*s%lf" */
509 strcat(fmtign,"%*s"); /* ignoring one more entry in the next loop */
510 if(sscanf(ptr,fmt,&temp))
512 temp+=header->UmbPos[i][0];
521 if (opt->bCalcTauInt)
523 /* save time series for autocorrelation analysis */
524 ntot=window->Ntot[i];
528 srenew(window->ztime[i],lennow[i]);
530 window->ztime[i][ntot]=temp;
543 else if (inttemp >= bins)
547 if(inttemp >= 0 && inttemp < bins)
549 window->Histo[i][inttemp]+=1.;
560 printf("time %f larger than tmax %f, stop reading pdo file\n",time,opt->tmax);
575 void enforceEqualWeights(t_UmbrellaWindow * window,int nWindows)
580 NEnforced=window[0].Ntot[0];
581 printf("\nFound -hist-eq. Enforcing equal weights for all histograms, \ni.e. doing a "
582 "non-weighted histogram analysis method. Ndata = %d\n",NEnforced);
583 /* enforce all histograms to have the same weight as the very first histogram */
585 for(j=0;j<nWindows;++j)
587 for(k=0;k<window[j].nPull;++k)
589 ratio=1.0*NEnforced/window[j].Ntot[k];
590 for(i=0;i<window[0].nBin;++i)
592 window[j].Histo[k][i]*=ratio;
594 window[j].N[k]=(int)(ratio*window[j].N[k] + 0.5);
599 /* Simple linear interpolation between two given tabulated points */
600 double tabulated_pot(double dist, t_UmbrellaOptions *opt)
605 jl=floor((dist-opt->tabMin)/opt->tabDz);
607 if (jl<0 || ju>=opt->tabNbins)
609 gmx_fatal(FARGS,"Distance %f out of bounds of tabulated potential (jl=%d, ju=%d).\n"
610 "Provide an extended table.",dist,jl,ju);
614 dz=dist-opt->tabX[jl];
615 dp=(pu-pl)*dz/opt->tabDz;
620 /* Don't worry, that routine does not mean we compute the PMF in limited precision.
621 After rapid convergence (using only substiantal contributions), we always switch to
623 void setup_acc_wham(double *profile,t_UmbrellaWindow * window,int nWindows,
624 t_UmbrellaOptions *opt)
626 int i,j,k,nGrptot=0,nContrib=0,nTot=0;
627 double U,min=opt->min,dz=opt->dz,temp,ztot_half,distance,ztot,contrib1,contrib2;
628 gmx_bool bAnyContrib;
630 static double wham_contrib_lim;
634 for(i=0;i<nWindows;++i)
636 nGrptot+=window[i].nPull;
638 wham_contrib_lim=opt->Tolerance/nGrptot;
641 ztot=opt->max-opt->min;
644 for(i=0;i<nWindows;++i)
646 if ( ! window[i].bContrib)
648 snew(window[i].bContrib,window[i].nPull);
650 for(j=0;j<window[i].nPull;++j)
652 if ( ! window[i].bContrib[j])
653 snew(window[i].bContrib[j],opt->bins);
655 for(k=0;k<opt->bins;++k)
657 temp=(1.0*k+0.5)*dz+min;
658 distance = temp - window[i].pos[j]; /* distance to umbrella center */
660 { /* in cyclic wham: */
661 if (distance > ztot_half) /* |distance| < ztot_half */
663 else if (distance < -ztot_half)
666 /* Note: there are two contributions to bin k in the wham equations:
667 i) N[j]*exp(- U/(8.314e-3*opt->Temperature) + window[i].z[j])
668 ii) exp(- U/(8.314e-3*opt->Temperature))
669 where U is the umbrella potential
670 If any of these number is larger wham_contrib_lim, I set contrib=TRUE
674 U=0.5*window[i].k[j]*sqr(distance); /* harmonic potential assumed. */
676 U=tabulated_pot(distance,opt); /* Use tabulated potential */
678 contrib1=profile[k]*exp(- U/(8.314e-3*opt->Temperature));
679 contrib2=window[i].N[j]*exp(- U/(8.314e-3*opt->Temperature) + window[i].z[j]);
680 window[i].bContrib[j][k] = (contrib1 > wham_contrib_lim || contrib2 > wham_contrib_lim);
681 bAnyContrib = (bAnyContrib | window[i].bContrib[j][k]);
682 if (window[i].bContrib[j][k])
686 /* If this histo is far outside min and max all bContrib may be FALSE,
687 causing a floating point exception later on. To avoid that, switch
690 for(k=0;k<opt->bins;++k)
691 window[i].bContrib[j][k]=TRUE;
695 printf("Initialized rapid wham stuff (contrib tolerance %g)\n"
696 "Evaluating only %d of %d expressions.\n\n",wham_contrib_lim,nContrib, nTot);
699 printf("Updated rapid wham stuff. (evaluating only %d of %d contributions)\n",
705 void calc_profile(double *profile,t_UmbrellaWindow * window, int nWindows,
706 t_UmbrellaOptions *opt, gmx_bool bExact)
709 double num,ztot_half,ztot,distance,min=opt->min,dz=opt->dz;
710 double denom,U=0,temp=0,invg;
712 ztot=opt->max-opt->min;
715 for(i=0;i<opt->bins;++i)
718 for(j=0;j<nWindows;++j)
720 for(k=0;k<window[j].nPull;++k)
722 invg = 1.0/window[j].g[k] * window[j].bsWeight[k];
723 temp = (1.0*i+0.5)*dz+min;
724 num += invg*window[j].Histo[k][i];
726 if (! (bExact || window[j].bContrib[k][i]))
728 distance = temp - window[j].pos[k]; /* distance to umbrella center */
730 { /* in cyclic wham: */
731 if (distance > ztot_half) /* |distance| < ztot_half */
733 else if (distance < -ztot_half)
738 U=0.5*window[j].k[k]*sqr(distance); /* harmonic potential assumed. */
740 U=tabulated_pot(distance,opt); /* Use tabulated potential */
741 denom+=invg*window[j].N[k]*exp(- U/(8.314e-3*opt->Temperature) + window[j].z[k]);
744 profile[i]=num/denom;
749 double calc_z(double * profile,t_UmbrellaWindow * window, int nWindows,
750 t_UmbrellaOptions *opt, gmx_bool bExact)
753 double U=0,min=opt->min,dz=opt->dz,temp,ztot_half,distance,ztot,totalMax;
754 double MAX=-1e20, total=0;
756 ztot=opt->max-opt->min;
759 for(i=0;i<nWindows;++i)
761 for(j=0;j<window[i].nPull;++j)
764 for(k=0;k<window[i].nBin;++k)
766 if (! (bExact || window[i].bContrib[j][k]))
768 temp=(1.0*k+0.5)*dz+min;
769 distance = temp - window[i].pos[j]; /* distance to umbrella center */
771 { /* in cyclic wham: */
772 if (distance > ztot_half) /* |distance| < ztot_half */
774 else if (distance < -ztot_half)
779 U=0.5*window[i].k[j]*sqr(distance); /* harmonic potential assumed. */
781 U=tabulated_pot(distance,opt); /* Use tabulated potential */
783 total+=profile[k]*exp(-U/(8.314e-3*opt->Temperature));
785 /* Avoid floating point exception if window is far outside min and max */
790 temp = fabs(total - window[i].z[j]);
796 window[i].z[j] = total;
802 void symmetrizeProfile(double* profile,t_UmbrellaOptions *opt)
805 double *prof2,bins=opt->bins,min=opt->min,max=opt->max,dz=opt->dz,zsym,deltaz,profsym;
808 if (min>0. || max<0.)
809 gmx_fatal(FARGS,"Cannot symmetrize profile around z=0 with min=%f and max=%f\n",
818 /* bin left of zsym */
819 j=floor((zsym-min)/dz-0.5);
820 if (j>=0 && (j+1)<bins)
822 /* interpolate profile linearly between bins j and j+1 */
825 profsym=profile[j] + (profile[j+1]-profile[j])/dz*deltaz;
826 /* average between left and right */
827 prof2[i]=0.5*(profsym+profile[i]);
835 memcpy(profile,prof2,bins*sizeof(double));
839 void prof_normalization_and_unit(double * profile, t_UmbrellaOptions *opt)
842 double unit_factor=1., R_MolarGasConst, diff;
844 R_MolarGasConst=8.314472e-3; /* in kJ/(mol*K) */
847 /* Not log? Nothing to do! */
851 /* Get profile in units of kT, kJ/mol, or kCal/mol */
852 if (opt->unit == en_kT)
854 else if (opt->unit == en_kJ)
855 unit_factor=R_MolarGasConst*opt->Temperature;
856 else if (opt->unit == en_kCal)
857 unit_factor=R_MolarGasConst*opt->Temperature/4.1868;
859 gmx_fatal(FARGS,"Sorry, I don't know this energy unit.");
863 profile[i]=-log(profile[i])*unit_factor;
865 /* shift to zero at z=opt->zProf0 */
872 /* Get bin with shortest distance to opt->zProf0
873 (-0.5 from bin position and +0.5 from rounding cancel) */
874 imin=(int)((opt->zProf0-opt->min)/opt->dz);
887 void getRandomIntArray(int nPull,int blockLength,int* randomArray,gmx_rng_t rng)
889 int ipull,blockBase,nr,ipullRandom;
894 for (ipull=0; ipull<nPull; ipull++)
896 blockBase=(ipull/blockLength)*blockLength;
898 { /* make sure nothing bad happens in the last block */
899 nr=(int)(gmx_rng_uniform_real(rng)*blockLength);
900 ipullRandom = blockBase + nr;
901 } while (ipullRandom >= nPull);
902 if (ipullRandom<0 || ipullRandom>=nPull)
903 gmx_fatal(FARGS,"Ups, random iWin = %d, nPull = %d, nr = %d, "
904 "blockLength = %d, blockBase = %d\n",
905 ipullRandom,nPull,nr,blockLength,blockBase);
906 randomArray[ipull]=ipullRandom;
908 /*for (ipull=0; ipull<nPull; ipull++)
909 printf("%d ",randomArray[ipull]); printf("\n"); */
912 void copy_pullgrp_to_synthwindow(t_UmbrellaWindow *synthWindow,
913 t_UmbrellaWindow *thisWindow,int pullid)
915 synthWindow->N [0]=thisWindow->N [pullid];
916 synthWindow->Histo [0]=thisWindow->Histo [pullid];
917 synthWindow->pos [0]=thisWindow->pos [pullid];
918 synthWindow->z [0]=thisWindow->z [pullid];
919 synthWindow->k [0]=thisWindow->k [pullid];
920 synthWindow->bContrib[0]=thisWindow->bContrib [pullid];
921 synthWindow->g [0]=thisWindow->g [pullid];
922 synthWindow->bsWeight[0]=thisWindow->bsWeight [pullid];
925 /* Calculate cumulative distribution function of of all histograms. They
926 allow to create random number sequences
927 which are distributed according to the histograms. Required to generate
928 the "synthetic" histograms for the Bootstrap method */
929 void calc_cumulatives(t_UmbrellaWindow *window,int nWindows,
930 t_UmbrellaOptions *opt,const char *fnhist)
939 snew(fn,strlen(fnhist)+10);
940 snew(buf,strlen(fnhist)+10);
941 sprintf(fn,"%s_cumul.xvg",strncpy(buf,fnhist,strlen(fnhist)-4));
942 fp=xvgropen(fn,"CDFs of umbrella windows","z","CDF",opt->oenv);
946 for (i=0; i<nWindows; i++)
948 snew(window[i].cum,window[i].nPull);
949 for (j=0; j<window[i].nPull; j++)
951 snew(window[i].cum[j],nbin+1);
952 window[i].cum[j][0]=0.;
953 for (k=1; k<=nbin; k++)
954 window[i].cum[j][k] = window[i].cum[j][k-1]+window[i].Histo[j][k-1];
956 /* normalize CDFs. Ensure cum[nbin]==1 */
957 last = window[i].cum[j][nbin];
958 for (k=0; k<=nbin; k++)
959 window[i].cum[j][k] /= last;
963 printf("Cumulative distriubtion functions of all histograms created.\n");
966 for (k=0; k<=nbin; k++)
968 fprintf(fp,"%g\t",opt->min+k*opt->dz);
969 for (i=0; i<nWindows; i++)
970 for (j=0; j<window[i].nPull; j++)
971 fprintf(fp,"%g\t",window[i].cum[j][k]);
974 printf("Wrote cumulative distribution functions to %s\n",fn);
982 /* Return j such that xx[j] <= x < xx[j+1] */
983 void searchCumulative(double xx[], int n, double x, int *j)
1005 void create_synthetic_histo(t_UmbrellaWindow *synthWindow, t_UmbrellaWindow *thisWindow,
1006 int pullid,t_UmbrellaOptions *opt)
1008 int N,i,nbins,r_index,ibin;
1009 double r,tausteps=0.0,a,ap,dt,x,invsqrt2,g,y,sig=0.,z,mu=0.;
1012 N=thisWindow->N[pullid];
1014 nbins=thisWindow->nBin;
1016 /* tau = autocorrelation time */
1017 if (opt->tauBootStrap>0.0)
1018 tausteps=opt->tauBootStrap/dt;
1019 else if (opt->bTauIntGiven || opt->bCalcTauInt)
1021 /* calc tausteps from g=1+2tausteps */
1022 g=thisWindow->g[pullid];
1028 "When generating hypothetical trajctories from given umbrella histograms,\n"
1029 "autocorrelation times (ACTs) are required. Otherwise the statistical error\n"
1030 "cannot be predicted. You have 3 options:\n"
1031 "1) Make g_wham estimate the ACTs (options -ac and -acsig).\n"
1032 "2) Calculate the ACTs by yourself (e.g. with g_analyze) and provide them\n");
1034 " with option -iiact for all umbrella windows.\n"
1035 "3) If all ACTs are identical and know, you can define them with -bs-tau.\n"
1036 " Use option (3) only if you are sure what you're doing, you may severely\n"
1037 " underestimate the error if a too small ACT is given.\n");
1038 gmx_fatal(FARGS,errstr);
1041 synthWindow->N [0]=N;
1042 synthWindow->pos [0]=thisWindow->pos[pullid];
1043 synthWindow->z [0]=thisWindow->z[pullid];
1044 synthWindow->k [0]=thisWindow->k[pullid];
1045 synthWindow->bContrib[0]=thisWindow->bContrib[pullid];
1046 synthWindow->g [0]=thisWindow->g [pullid];
1047 synthWindow->bsWeight[0]=thisWindow->bsWeight[pullid];
1049 for (i=0;i<nbins;i++)
1050 synthWindow->Histo[0][i]=0.;
1052 if (opt->bsMethod==bsMethod_trajGauss)
1054 sig = thisWindow->sigma [pullid];
1055 mu = thisWindow->aver [pullid];
1058 /* Genrate autocorrelated Gaussian random variable with autocorrelation time tau
1060 If x and y are random numbers from N(0,1) (Gaussian with average 0 and sigma=1),
1062 z = a*x + sqrt(1-a^2)*y
1063 is also from N(0,1), and cov(z,x) = a. Thus, by gerenating a sequence
1064 x' = a*x + sqrt(1-a^2)*y, the sequnce x(t) is from N(0,1) and has an autocorrelation
1066 C(t) = exp(-t/tau) with tau=-1/ln(a)
1068 Then, use error function to turn the Gaussian random variable into a uniformly
1069 distributed one in [0,1]. Eventually, use cumulative distribution function of
1070 histogram to get random variables distributed according to histogram.
1071 Note: The ACT of the flat distribution and of the generated histogram is not
1072 100% exactly tau, but near tau (my test was 3.8 instead of 4).
1074 a=exp(-1.0/tausteps);
1076 invsqrt2=1./sqrt(2.0);
1078 /* init random sequence */
1079 x=gmx_rng_gaussian_table(opt->rng);
1081 if (opt->bsMethod==bsMethod_traj)
1083 /* bootstrap points from the umbrella histograms */
1086 y=gmx_rng_gaussian_table(opt->rng);
1088 /* get flat distribution in [0,1] using cumulative distribution function of Gauusian
1089 Note: CDF(Gaussian) = 0.5*{1+erf[x/sqrt(2)]}
1091 r=0.5*(1+gmx_erf(x*invsqrt2));
1092 searchCumulative(thisWindow->cum[pullid],nbins+1 ,r,&r_index);
1093 synthWindow->Histo[0][r_index]+=1.;
1096 else if (opt->bsMethod==bsMethod_trajGauss)
1098 /* bootstrap points from a Gaussian with the same average and sigma
1099 as the respective umbrella histogram. The idea was, that -given
1100 limited sampling- the bootstrapped histograms are otherwise biased
1101 from the limited sampling of the US histos. However, bootstrapping from
1102 the Gaussian seems to yield a similar estimate. */
1106 y=gmx_rng_gaussian_table(opt->rng);
1109 ibin=floor((z-opt->min)/opt->dz);
1113 while ( (ibin+=nbins) < 0);
1114 else if (ibin>=nbins)
1115 while ( (ibin-=nbins) >= nbins);
1118 if (ibin>=0 && ibin<nbins)
1120 synthWindow->Histo[0][ibin]+=1.;
1127 gmx_fatal(FARGS,"Unknown bsMethod (id %d). That should not happen.\n",opt->bsMethod);
1132 void print_histograms(const char *fnhist, t_UmbrellaWindow * window, int nWindows,
1133 int bs_index,t_UmbrellaOptions *opt)
1135 char *fn=0,*buf=0,title[256];
1142 strcpy(title,"Umbrella histograms");
1146 snew(fn,strlen(fnhist)+10);
1147 snew(buf,strlen(fnhist)+1);
1148 sprintf(fn,"%s_bs%d.xvg",strncpy(buf,fnhist,strlen(fnhist)-4),bs_index);
1149 sprintf(title,"Umbrella histograms. Bootstrap #%d",bs_index);
1152 fp=xvgropen(fn,title,"z","count",opt->oenv);
1155 /* Write histograms */
1158 fprintf(fp,"%e\t",(double)(l+0.5)*opt->dz+opt->min);
1159 for(i=0;i<nWindows;++i)
1161 for(j=0;j<window[i].nPull;++j)
1163 fprintf(fp,"%e\t",window[i].Histo[j][l]);
1170 printf("Wrote %s\n",fn);
1178 int func_wham_is_larger(const void *a, const void *b)
1192 void setRandomBsWeights(t_UmbrellaWindow *synthwin,int nAllPull, t_UmbrellaOptions *opt)
1199 /* generate ordered random numbers between 0 and nAllPull */
1200 for (i=0; i<nAllPull-1; i++)
1202 r[i] = gmx_rng_uniform_real(opt->rng) * nAllPull;
1204 qsort((void *)r,nAllPull-1, sizeof(double), &func_wham_is_larger);
1205 r[nAllPull-1]=1.0*nAllPull;
1207 synthwin[0].bsWeight[0]=r[0];
1208 for (i=1; i<nAllPull; i++)
1210 synthwin[i].bsWeight[0]=r[i]-r[i-1];
1213 /* avoid to have zero weight by adding a tiny value */
1214 for (i=0; i<nAllPull; i++)
1215 if (synthwin[i].bsWeight[0] < 1e-5)
1216 synthwin[i].bsWeight[0] = 1e-5;
1221 void do_bootstrapping(const char *fnres, const char* fnprof, const char *fnhist,
1222 char* ylabel, double *profile,
1223 t_UmbrellaWindow * window, int nWindows, t_UmbrellaOptions *opt)
1225 t_UmbrellaWindow * synthWindow;
1226 double *bsProfile,*bsProfiles_av, *bsProfiles_av2,maxchange=1e20,tmp,stddev;
1227 int i,j,*randomArray=0,winid,pullid,ib;
1228 int iAllPull,nAllPull,*allPull_winId,*allPull_pullId;
1230 gmx_bool bExact=FALSE;
1232 /* init random generator */
1233 if (opt->bsSeed==-1)
1234 opt->rng=gmx_rng_init(gmx_rng_make_seed());
1236 opt->rng=gmx_rng_init(opt->bsSeed);
1238 snew(bsProfile, opt->bins);
1239 snew(bsProfiles_av, opt->bins);
1240 snew(bsProfiles_av2,opt->bins);
1242 /* Create array of all pull groups. Note that different windows
1243 may have different nr of pull groups
1244 First: Get total nr of pull groups */
1246 for (i=0;i<nWindows;i++)
1247 nAllPull+=window[i].nPull;
1248 snew(allPull_winId,nAllPull);
1249 snew(allPull_pullId,nAllPull);
1251 /* Setup one array of all pull groups */
1252 for (i=0;i<nWindows;i++)
1254 for (j=0;j<window[i].nPull;j++)
1256 allPull_winId[iAllPull]=i;
1257 allPull_pullId[iAllPull]=j;
1262 /* setup stuff for synthetic windows */
1263 snew(synthWindow,nAllPull);
1264 for (i=0;i<nAllPull;i++)
1266 synthWindow[i].nPull=1;
1267 synthWindow[i].nBin=opt->bins;
1268 snew(synthWindow[i].Histo,1);
1269 if (opt->bsMethod == bsMethod_traj || opt->bsMethod == bsMethod_trajGauss)
1270 snew(synthWindow[i].Histo[0],opt->bins);
1271 snew(synthWindow[i].N,1);
1272 snew(synthWindow[i].pos,1);
1273 snew(synthWindow[i].z,1);
1274 snew(synthWindow[i].k,1);
1275 snew(synthWindow[i].bContrib,1);
1276 snew(synthWindow[i].g,1);
1277 snew(synthWindow[i].bsWeight,1);
1280 switch(opt->bsMethod)
1283 snew(randomArray,nAllPull);
1284 printf("\n\nWhen computing statistical errors by bootstrapping entire histograms:\n");
1285 please_cite(stdout,"Hub2006");
1287 case bsMethod_BayesianHist :
1288 /* just copy all histogams into synthWindow array */
1289 for (i=0;i<nAllPull;i++)
1291 winid =allPull_winId [i];
1292 pullid=allPull_pullId[i];
1293 copy_pullgrp_to_synthwindow(synthWindow+i,window+winid,pullid);
1297 case bsMethod_trajGauss:
1298 calc_cumulatives(window,nWindows,opt,fnhist);
1301 gmx_fatal(FARGS,"Unknown bootstrap method. That should not have happened.\n");
1304 /* do bootstrapping */
1305 fp=xvgropen(fnprof,"Boot strap profiles","z",ylabel,opt->oenv);
1306 for (ib=0;ib<opt->nBootStrap;ib++)
1308 printf(" *******************************************\n"
1309 " ******** Start bootstrap nr %d ************\n"
1310 " *******************************************\n",ib+1);
1312 switch(opt->bsMethod)
1315 /* bootstrap complete histograms from given histograms */
1316 getRandomIntArray(nAllPull,opt->histBootStrapBlockLength,randomArray,opt->rng);
1317 for (i=0;i<nAllPull;i++){
1318 winid =allPull_winId [randomArray[i]];
1319 pullid=allPull_pullId[randomArray[i]];
1320 copy_pullgrp_to_synthwindow(synthWindow+i,window+winid,pullid);
1323 case bsMethod_BayesianHist:
1324 /* keep histos, but assign random weights ("Bayesian bootstrap") */
1325 setRandomBsWeights(synthWindow,nAllPull,opt);
1328 case bsMethod_trajGauss:
1329 /* create new histos from given histos, that is generate new hypothetical
1331 for (i=0;i<nAllPull;i++)
1333 winid=allPull_winId[i];
1334 pullid=allPull_pullId[i];
1335 create_synthetic_histo(synthWindow+i,window+winid,pullid,opt);
1340 /* write histos in case of verbose output */
1341 if (opt->bs_verbose)
1342 print_histograms(fnhist,synthWindow,nAllPull,ib,opt);
1348 memcpy(bsProfile,profile,opt->bins*sizeof(double)); /* use profile as guess */
1351 if ( (i%opt->stepUpdateContrib) == 0)
1352 setup_acc_wham(bsProfile,synthWindow,nAllPull,opt);
1353 if (maxchange<opt->Tolerance)
1355 if (((i%opt->stepchange) == 0 || i==1) && !i==0)
1356 printf("\t%4d) Maximum change %e\n",i,maxchange);
1357 calc_profile(bsProfile,synthWindow,nAllPull,opt,bExact);
1359 } while( (maxchange=calc_z(bsProfile, synthWindow, nAllPull, opt,bExact)) > opt->Tolerance || !bExact);
1360 printf("\tConverged in %d iterations. Final maximum change %g\n",i,maxchange);
1363 prof_normalization_and_unit(bsProfile,opt);
1365 /* symmetrize profile around z=0 */
1367 symmetrizeProfile(bsProfile,opt);
1369 /* save stuff to get average and stddev */
1370 for (i=0;i<opt->bins;i++)
1373 bsProfiles_av[i]+=tmp;
1374 bsProfiles_av2[i]+=tmp*tmp;
1375 fprintf(fp,"%e\t%e\n",(i+0.5)*opt->dz+opt->min,tmp);
1381 /* write average and stddev */
1382 fp=xvgropen(fnres,"Average and stddev from bootstrapping","z",ylabel,opt->oenv);
1383 fprintf(fp,"@TYPE xydy\n");
1384 for (i=0;i<opt->bins;i++)
1386 bsProfiles_av [i]/=opt->nBootStrap;
1387 bsProfiles_av2[i]/=opt->nBootStrap;
1388 tmp=bsProfiles_av2[i]-sqr(bsProfiles_av[i]);
1389 stddev=(tmp>=0.) ? sqrt(tmp) : 0.; /* Catch rouding errors */
1390 fprintf(fp,"%e\t%e\t%e\n",(i+0.5)*opt->dz+opt->min,bsProfiles_av [i],stddev);
1393 printf("Wrote boot strap result to %s\n",fnres);
1396 int whaminFileType(char *fn)
1400 if (strcmp(fn+len-3,"tpr")==0)
1402 else if (strcmp(fn+len-3,"xvg")==0 || strcmp(fn+len-6,"xvg.gz")==0)
1403 return whamin_pullxf;
1404 else if (strcmp(fn+len-3,"pdo")==0 || strcmp(fn+len-6,"pdo.gz")==0)
1407 gmx_fatal(FARGS,"Unknown file type of %s. Should be tpr, xvg, or pdo.\n",fn);
1408 return whamin_unknown;
1411 void read_wham_in(const char *fn,char ***filenamesRet, int *nfilesRet,
1412 t_UmbrellaOptions *opt)
1414 char **filename=0,tmp[STRLEN];
1415 int nread,sizenow,i,block=1;
1421 while (fscanf(fp,"%s",tmp) != EOF)
1423 if (strlen(tmp)>=WHAM_MAXFILELEN)
1424 gmx_fatal(FARGS,"Filename too long. Only %d characters allowed\n",WHAM_MAXFILELEN);
1428 srenew(filename,sizenow);
1429 for (i=sizenow-block;i<sizenow;i++)
1430 snew(filename[i],WHAM_MAXFILELEN);
1432 strcpy(filename[nread],tmp);
1434 printf("Found file %s in %s\n",filename[nread],fn);
1437 *filenamesRet=filename;
1442 FILE *open_pdo_pipe(const char *fn, t_UmbrellaOptions *opt,gmx_bool *bPipeOpen)
1444 char Buffer[1024],gunzip[1024],*Path=0;
1446 static gmx_bool bFirst=1;
1448 /* gzipped pdo file? */
1449 if ((strcmp(fn+strlen(fn)-3,".gz")==0))
1451 /* search gunzip executable */
1452 if(!(Path=getenv("GMX_PATH_GZIP")))
1454 if (gmx_fexist("/bin/gunzip"))
1455 sprintf(gunzip,"%s","/bin/gunzip");
1456 else if (gmx_fexist("/usr/bin/gunzip"))
1457 sprintf(gunzip,"%s","/usr/bin/gunzip");
1459 gmx_fatal(FARGS,"Cannot find executable gunzip in /bin or /usr/bin.\n"
1460 "You may want to define the path to gunzip "
1461 "with the environment variable GMX_PATH_GZIP.",gunzip);
1465 sprintf(gunzip,"%s/gunzip",Path);
1466 if (!gmx_fexist(gunzip))
1467 gmx_fatal(FARGS,"Cannot find executable %s. Please define the path to gunzip"
1468 " in the environmental varialbe GMX_PATH_GZIP.",gunzip);
1472 printf("Using gunzig executable %s\n",gunzip);
1475 if (!gmx_fexist(fn))
1477 gmx_fatal(FARGS,"File %s does not exist.\n",fn);
1479 sprintf(Buffer,"%s -c < %s",gunzip,fn);
1481 printf("Executing command '%s'\n",Buffer);
1483 if((pipe=popen(Buffer,"r"))==NULL)
1485 gmx_fatal(FARGS,"Unable to open pipe to `%s'\n",Buffer);
1488 gmx_fatal(FARGS,"Cannot open a compressed file on platform without pipe support");
1493 pipe=ffopen(fn,"r");
1500 void pdo_close_file(FILE *fp)
1509 /* Reading pdo files */
1510 void read_pdo_files(char **fn, int nfiles, t_UmbrellaHeader* header,
1511 t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
1514 real mintmp,maxtmp,done=0.;
1517 /* char Buffer0[1000]; */
1520 gmx_fatal(FARGS,"No files found. Hick.");
1522 /* if min and max are not given, get min and max from the input files */
1525 printf("Automatic determination of boundaries from %d pdo files...\n",nfiles);
1528 for(i=0;i<nfiles;++i)
1530 file=open_pdo_pipe(fn[i],opt,&bPipeOpen);
1531 /*fgets(Buffer0,999,file);
1532 fprintf(stderr,"First line '%s'\n",Buffer0); */
1533 done=100.0*(i+1)/nfiles;
1534 printf("\rOpening %s ... [%2.0f%%]",fn[i],done); fflush(stdout);
1537 read_pdo_header(file,header,opt);
1538 /* here only determine min and max of this window */
1539 read_pdo_data(file,header,i,NULL,opt,TRUE,&mintmp,&maxtmp);
1540 if (maxtmp>opt->max)
1542 if (mintmp<opt->min)
1545 pdo_close_file(file);
1550 printf("\nDetermined boundaries to %f and %f\n\n",opt->min,opt->max);
1551 if (opt->bBoundsOnly)
1553 printf("Found option -boundsonly, now exiting.\n");
1557 /* store stepsize in profile */
1558 opt->dz=(opt->max-opt->min)/opt->bins;
1560 /* Having min and max, we read in all files */
1561 /* Loop over all files */
1562 for(i=0;i<nfiles;++i)
1564 done=100.0*(i+1)/nfiles;
1565 printf("\rOpening %s ... [%2.0f%%]",fn[i],done); fflush(stdout);
1568 file=open_pdo_pipe(fn[i],opt,&bPipeOpen);
1569 read_pdo_header(file,header,opt);
1570 /* load data into window */
1571 read_pdo_data(file,header,i,window,opt,FALSE,NULL,NULL);
1572 if ((window+i)->Ntot[0] == 0.0)
1573 fprintf(stderr,"\nWARNING, no data points read from file %s (check -b option)\n", fn[i]);
1575 pdo_close_file(file);
1580 for(i=0;i<nfiles;++i)
1585 #define int2YN(a) (((a)==0)?("N"):("Y"))
1587 void read_tpr_header(const char *fn,t_UmbrellaHeader* header,t_UmbrellaOptions *opt)
1594 /* printf("Reading %s \n",fn); */
1595 read_tpx_state(fn,&ir,&state,NULL,NULL);
1597 if (ir.ePull != epullUMBRELLA)
1598 gmx_fatal(FARGS,"This is not a tpr of an umbrella simulation. Found pull type \"%s\" "
1599 " (ir.ePull = %d)\n", epull_names[ir.ePull],ir.ePull);
1601 /* nr of pull groups */
1604 gmx_fatal(FARGS,"This is not a tpr of umbrella simulation. Found only %d pull groups\n",ngrp);
1606 header->npullgrps=ir.pull->ngrp;
1607 header->pull_geometry=ir.pull->eGeom;
1608 copy_ivec(ir.pull->dim,header->pull_dim);
1609 header->pull_ndim=header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2];
1610 if (header->pull_geometry==epullgPOS && header->pull_ndim>1)
1612 gmx_fatal(FARGS,"Found pull geometry 'position' and more than 1 pull dimension (%d).\n"
1613 "Hence, the pull potential does not correspond to a one-dimensional umbrella potential.\n"
1614 "If you have some special umbrella setup you may want to write your own pdo files\n"
1615 "and feed them into g_wham. Check g_wham -h !\n",header->pull_ndim);
1617 snew(header->k,ngrp);
1618 snew(header->init_dist,ngrp);
1619 snew(header->umbInitDist,ngrp);
1621 /* only z-direction with epullgCYL? */
1622 if (header->pull_geometry == epullgCYL)
1624 if (header->pull_dim[XX] || header->pull_dim[YY] || (!header->pull_dim[ZZ]))
1625 gmx_fatal(FARGS,"With pull geometry 'cylinder', expected pulling in Z direction only.\n"
1626 "However, found dimensions [%s %s %s]\n",
1627 int2YN(header->pull_dim[XX]),int2YN(header->pull_dim[YY]),
1628 int2YN(header->pull_dim[ZZ]));
1631 for (i=0;i<ngrp;i++)
1633 header->k[i]=ir.pull->grp[i+1].k;
1634 if (header->k[i]==0.0)
1635 gmx_fatal(FARGS,"Pull group %d has force constant of of 0.0 in %s.\n"
1636 "That doesn't seem to be an Umbrella tpr.\n",
1638 copy_rvec(ir.pull->grp[i+1].init,header->init_dist[i]);
1640 /* initial distance to reference */
1641 switch(header->pull_geometry)
1645 if (header->pull_dim[d])
1646 header->umbInitDist[i]=header->init_dist[i][d];
1649 /* umbrella distance stored in init_dist[i][0] for geometry cylinder (not in ...[i][ZZ]) */
1653 header->umbInitDist[i]=header->init_dist[i][0];
1656 gmx_fatal(FARGS,"Pull geometry %s not supported\n",epullg_names[header->pull_geometry]);
1660 if (opt->verbose || first)
1662 printf("File %s, %d groups, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n",
1663 fn,header->npullgrps,epullg_names[header->pull_geometry],
1664 int2YN(header->pull_dim[0]),int2YN(header->pull_dim[1]),int2YN(header->pull_dim[2]),
1666 for (i=0;i<ngrp;i++)
1667 printf("\tgrp %d) k = %-5g position = %g\n",i,header->k[i],header->umbInitDist[i]);
1669 if (!opt->verbose && first)
1670 printf("\tUse option -v to see this output for all input tpr files\n");
1676 double dist_ndim(double **dx,int ndim,int line)
1680 for (i=0;i<ndim;i++)
1681 r2+=sqr(dx[i][line]);
1685 void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header,
1686 t_UmbrellaWindow * window,
1687 t_UmbrellaOptions *opt,
1688 gmx_bool bGetMinMax,real *mintmp,real *maxtmp)
1690 double **y=0,pos=0.,t,force,time0=0.,dt;
1691 int ny,nt,bins,ibin,i,g,dstep=1,nColPerGrp,nColRefOnce,nColRefEachGrp,nColExpect,ntot;
1692 real min,max,minfound=1e20,maxfound=-1e20;
1693 gmx_bool dt_ok,timeok,bHaveForce;
1694 const char *quantity;
1695 const int blocklen=4096;
1699 in force output pullf.xvg:
1700 No reference, one column per pull group
1701 in position output pullx.xvg (not cylinder)
1702 ndim reference, ndim columns per pull group
1703 in position output pullx.xvg (in geometry cylinder):
1704 ndim*2 columns per pull group (ndim for ref, ndim for group)
1707 nColPerGrp = opt->bPullx ? header->pull_ndim : 1;
1708 quantity = opt->bPullx ? "position" : "force";
1712 if (header->pull_geometry == epullgCYL)
1714 /* Geometry cylinder -> reference group before each pull group */
1715 nColRefEachGrp=header->pull_ndim;
1720 /* Geometry NOT cylinder -> reference group only once after time column */
1722 nColRefOnce=header->pull_ndim;
1725 else /* read forces, no reference groups */
1731 nColExpect = 1 + nColRefOnce + header->npullgrps*(nColRefEachGrp+nColPerGrp);
1732 bHaveForce = opt->bPullf;
1734 /* With geometry "distance" or "distance_periodic", only force reading is supported so far.
1735 That avoids the somewhat tedious extraction of the right columns from the pullx files
1736 to compute the distances projection on the vector. Sorry for the laziness. */
1737 if ( (header->pull_geometry==epullgDIR || header->pull_geometry==epullgDIRPBC)
1740 gmx_fatal(FARGS,"With pull geometries \"direction\" and \"direction_periodic\", only pull force "
1741 "reading \n(option -if) is supported at present, "
1742 "not pull position reading (options -ix).\nMake sure mdrun writes the pull "
1743 "forces (pullf.xvg files)\nand provide them to g_wham with option -if.",
1744 epullg_names[header->pull_geometry]);
1747 nt=read_xvg(fn,&y,&ny);
1749 /* Check consistency */
1752 gmx_fatal(FARGS,"Empty pull %s file %s\n",quantity,fn);
1754 if (ny != nColExpect)
1756 gmx_fatal(FARGS,"Found %d pull groups in %s,\n but %d data columns in %s (expected %d)\n"
1757 "\nMaybe you confused options -ix and -if ?\n",
1758 header->npullgrps,fntpr,ny-1,fn,nColExpect-1);
1762 printf("Found %d times and %d %s sets %s\n",nt,(ny-1)/nColPerGrp,quantity,fn);
1771 window->dt=y[0][1]-y[0][0];
1773 else if (opt->nBootStrap && opt->tauBootStrap!=0.0)
1775 fprintf(stderr,"\n *** WARNING, Could not determine time step in %s\n",fn);
1778 /* Need to alocate memory and set up structure */
1779 window->nPull=header->npullgrps;
1782 snew(window->Histo,window->nPull);
1783 snew(window->z,window->nPull);
1784 snew(window->k,window->nPull);
1785 snew(window->pos,window->nPull);
1786 snew(window->N, window->nPull);
1787 snew(window->Ntot, window->nPull);
1788 snew(window->g, window->nPull);
1789 snew(window->bsWeight, window->nPull);
1792 if (opt->bCalcTauInt)
1793 snew(window->ztime,window->nPull);
1796 snew(lennow,window->nPull);
1798 for(g=0;g<window->nPull;++g)
1801 window->bsWeight[g]=1.;
1802 snew(window->Histo[g],bins);
1803 window->k[g]=header->k[g];
1807 window->pos[g]=header->umbInitDist[g];
1808 if (opt->bCalcTauInt)
1809 window->ztime[g]=NULL;
1814 { /* only determine min and max */
1817 min=max=bins=0; /* Get rid of warnings */
1822 /* Do you want that time frame? */
1823 t=1.0/1000*( (int) ((y[0][i]*1000) + 0.5)); /* round time to fs */
1825 /* get time step of pdo file and get dstep from opt->dt */
1835 dstep=(int)(opt->dt/dt+0.5);
1840 window->dt=dt*dstep;
1843 dt_ok=(i%dstep == 0);
1844 timeok=(dt_ok && t >= opt->tmin && t <= opt->tmax);
1846 printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
1847 t,opt->tmin, opt->tmax, dt_ok,timeok); */
1851 for(g=0;g<header->npullgrps;++g)
1855 /* y has 1 time column y[0] and one column per force y[1],...,y[nGrps] */
1857 pos= -force/header->k[g] + header->umbInitDist[g];
1861 switch (header->pull_geometry)
1864 /* y has 1 time column y[0] and nColPerGrps columns per pull group;
1865 Distance to reference: */
1866 /* pos=dist_ndim(y+1+nColRef+g*nColPerGrp,header->pull_ndim,i); gmx 4.0 */
1867 pos=dist_ndim(y + 1 + nColRefOnce + g*nColPerGrp,header->pull_ndim,i);
1871 Time ref[ndim] group1[ndim] group2[ndim] ... */
1874 Time ref1[ndim] group1[ndim] ref2[ndim] group2[ndim] ... */
1876 /* * with geometry==position, we have the reference once (nColRefOnce==ndim), but
1877 no extra reference group columns before each group (nColRefEachGrp==0)
1879 * with geometry==cylinder, we have no initial ref group column (nColRefOnce==0),
1880 but ndim ref group colums before every group (nColRefEachGrp==ndim)
1881 Distance to reference: */
1882 pos=y[1 + nColRefOnce + nColRefEachGrp + g*(nColRefEachGrp+nColPerGrp)][i];
1885 gmx_fatal(FARGS,"Bad error, this error should have been catched before. Ups.\n");
1889 /* printf("grp %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
1899 if (opt->bCalcTauInt && !bGetMinMax)
1901 /* save time series for autocorrelation analysis */
1902 ntot=window->Ntot[g];
1903 /* printf("i %d, ntot %d, lennow[g] = %d\n",i,ntot,lennow[g]); */
1904 if (ntot>=lennow[g])
1906 lennow[g]+=blocklen;
1907 srenew(window->ztime[g],lennow[g]);
1909 window->ztime[g][ntot]=pos;
1912 ibin=(int) floor((pos-min)/(max-min)*bins);
1916 while ( (ibin+=bins) < 0);
1917 else if (ibin>=bins)
1918 while ( (ibin-=bins) >= bins);
1920 if(ibin >= 0 && ibin < bins)
1922 window->Histo[g][ibin]+=1.;
1929 else if (t>opt->tmax)
1932 printf("time %f larger than tmax %f, stop reading this pullx/pullf file\n",t,opt->tmax);
1947 void read_tpr_pullxf_files(char **fnTprs,char **fnPull,int nfiles,
1948 t_UmbrellaHeader* header,
1949 t_UmbrellaWindow *window, t_UmbrellaOptions *opt)
1954 printf("Reading %d tpr and pullf files\n",nfiles/2);
1956 /* min and max not given? */
1959 printf("Automatic determination of boundaries...\n");
1962 for (i=0;i<nfiles; i++)
1964 if (whaminFileType(fnTprs[i]) != whamin_tpr)
1965 gmx_fatal(FARGS,"Expected the %d'th file in input file to be a tpr file\n",i);
1966 read_tpr_header(fnTprs[i],header,opt);
1967 if (whaminFileType(fnPull[i]) != whamin_pullxf)
1968 gmx_fatal(FARGS,"Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n",i);
1969 read_pull_xf(fnPull[i],fnTprs[i],header,NULL,opt,TRUE,&mintmp,&maxtmp);
1970 if (maxtmp>opt->max)
1972 if (mintmp<opt->min)
1975 printf("\nDetermined boundaries to %f and %f\n\n",opt->min,opt->max);
1976 if (opt->bBoundsOnly)
1978 printf("Found option -boundsonly, now exiting.\n");
1982 /* store stepsize in profile */
1983 opt->dz=(opt->max-opt->min)/opt->bins;
1985 for (i=0;i<nfiles; i++)
1987 if (whaminFileType(fnTprs[i]) != whamin_tpr)
1988 gmx_fatal(FARGS,"Expected the %d'th file in input file to be a tpr file\n",i);
1989 read_tpr_header(fnTprs[i],header,opt);
1990 if (whaminFileType(fnPull[i]) != whamin_pullxf)
1991 gmx_fatal(FARGS,"Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n",i);
1992 read_pull_xf(fnPull[i],fnTprs[i],header,window+i,opt,FALSE,NULL,NULL);
1993 if (window[i].Ntot[0] == 0.0)
1994 fprintf(stderr,"\nWARNING, no data points read from file %s (check -b option)\n", fnPull[i]);
1997 for (i=0;i<nfiles; i++)
2006 /* Note: Here we consider tau[int] := int_0^inf ACF(t) as the integrated autocorrelation time.
2007 The factor `g := 1 + 2*tau[int]` subsequently enters the uncertainty.
2009 void readIntegratedAutocorrelationTimes(t_UmbrellaWindow *window,int nwins,t_UmbrellaOptions *opt,
2015 printf("Readging Integrated autocorrelation times from %s ...\n",fn);
2016 nlines=read_xvg(fn,&iact,&ny);
2018 gmx_fatal(FARGS,"Found %d lines with integrated autocorrelation times in %s.\nExpected %d",
2020 for (i=0;i<nlines;i++){
2021 if (window[i].nPull != ny)
2022 gmx_fatal(FARGS,"You are providing autocorrelation times with option -iiact and the\n"
2023 "number of pull groups is different in different simulations. That is not\n"
2024 "supported yet. Sorry.\n");
2025 for (ig=0;ig<window[i].nPull;ig++){
2026 /* compare Kumar et al, J Comp Chem 13, 1011-1021 (1992) */
2027 window[i].g[ig]=1+2*iact[ig][i]/window[i].dt;
2029 if (iact[ig][i] <= 0.0)
2030 fprintf(stderr,"\nWARNING, IACT = %f (window %d, group %d)\n", iact[ig][i],i,ig);
2036 /* Smooth autocorreltion times along the reaction coordinate. This is useful
2037 if the ACT is subject to high uncertainty in case if limited sampling. Note
2038 that -in case of limited sampling- the ACT may be severely underestimated.
2039 Note: the g=1+2tau are overwritten.
2040 if opt->bAllowReduceIact==FALSE, the ACTs are never reduced, only increased
2043 void smoothIact(t_UmbrellaWindow *window,int nwins,t_UmbrellaOptions *opt)
2046 double pos,dpos2,siglim,siglim2,gaufact,invtwosig2,w,weight,tausm;
2048 /* only evaluate within +- 3sigma of the Gausian */
2049 siglim=3.0*opt->sigSmoothIact;
2050 siglim2=dsqr(siglim);
2051 /* pre-factor of Gaussian */
2052 gaufact=1.0/(sqrt(2*M_PI)*opt->sigSmoothIact);
2053 invtwosig2=0.5/dsqr(opt->sigSmoothIact);
2055 for (i=0;i<nwins;i++)
2057 snew(window[i].tausmooth,window[i].nPull);
2058 for (ig=0;ig<window[i].nPull;ig++)
2062 pos=window[i].pos[ig];
2063 for (j=0;j<nwins;j++)
2065 for (jg=0;jg<window[j].nPull;jg++)
2067 dpos2=dsqr(window[j].pos[jg]-pos);
2069 w=gaufact*exp(-dpos2*invtwosig2);
2071 tausm+=w*window[j].tau[jg];
2072 /*printf("Weight %g dpos2=%g pos=%g gaufact=%g invtwosig2=%g\n",
2073 w,dpos2,pos,gaufact,invtwosig2); */
2078 if (opt->bAllowReduceIact || tausm>window[i].tau[ig])
2079 window[i].tausmooth[ig]=tausm;
2081 window[i].tausmooth[ig]=window[i].tau[ig];
2082 window[i].g[ig] = 1+2*tausm/window[i].dt;
2087 /* try to compute the autocorrelation time for each umbrealla window */
2088 #define WHAM_AC_ZERO_LIMIT 0.05
2089 void calcIntegratedAutocorrelationTimes(t_UmbrellaWindow *window,int nwins,
2090 t_UmbrellaOptions *opt, const char *fn)
2092 int i,ig,ncorr,ntot,j,k,*count,restart;
2093 real *corr,c0,dt,timemax,tmp;
2094 real *ztime,av,tausteps;
2098 fpcorr=xvgropen("hist_autocorr.xvg","Autocorrelation functions of umbrella windows",
2099 "time [ps]","autocorrelation function",opt->oenv);
2102 for (i=0;i<nwins;i++)
2104 printf("\rEstimating integrated autocorreltion times ... [%2.0f%%] ...",100.*(i+1)/nwins);
2106 ntot=window[i].Ntot[0];
2108 /* using half the maximum time as length of autocorrelation function */
2111 gmx_fatal(FARGS,"Tryig to estimtate autocorrelation time from only %d"
2112 " points. Provide more pull data!",ntot);
2114 /* snew(corrSq,ncorr); */
2118 snew(window[i].tau,window[i].nPull);
2119 restart=(int)(opt->acTrestart/dt+0.5);
2123 for (ig=0;ig<window[i].nPull;ig++)
2125 if (ntot != window[i].Ntot[ig])
2126 gmx_fatal(FARGS,"Encountered different nr of frames in different pull groups.\n"
2127 "That should not happen. (%d and %d)\n", ntot,window[i].Ntot[ig]);
2128 ztime=window[i].ztime[ig];
2130 /* calc autocorrelation function C(t) = < [z(tau)-<z>]*[z(tau+t)-<z>]> */
2131 for(j=0, av=0; (j<ntot); j++)
2134 for(k=0; (k<ncorr); k++)
2139 for(j=0; (j<ntot); j+=restart)
2140 for(k=0; (k<ncorr) && (j+k < ntot); k++)
2142 tmp=(ztime[j]-av)*(ztime[j+k]-av);
2144 /* corrSq[k] += tmp*tmp; */
2147 /* divide by nr of frames for each time displacement */
2148 for(k=0; (k<ncorr); k++)
2150 /* count probably = (ncorr-k+(restart-1))/restart; */
2151 corr[k] = corr[k]/count[k];
2152 /* variance of autocorrelation function */
2153 /* corrSq[k]=corrSq[k]/count[k]; */
2155 /* normalize such that corr[0] == 0 */
2157 for(k=0; (k<ncorr); k++)
2160 /* corrSq[k]*=c0*c0; */
2163 /* write ACFs in verbose mode */
2166 for(k=0; (k<ncorr); k++)
2167 fprintf(fpcorr,"%g %g\n",k*dt,corr[k]);
2168 fprintf(fpcorr,"&\n");
2171 /* esimate integrated correlation time, fitting is too unstable */
2172 tausteps = 0.5*corr[0];
2173 /* consider corr below WHAM_AC_ZERO_LIMIT as noise */
2174 for(j=1; (j<ncorr) && (corr[j]>WHAM_AC_ZERO_LIMIT); j++)
2175 tausteps += corr[j];
2177 /* g = 1+2*tau, see. Ferrenberg/Swendsen, PRL 63:1195 (1989) or
2178 Kumar et al, eq. 28 ff. */
2179 window[i].tau[ig] = tausteps*dt;
2180 window[i].g[ig] = 1+2*tausteps;
2181 /* printf("win %d, group %d, estimated correlation time = %g ps\n",i,ig,window[i].tau[ig]); */
2190 /* plot IACT along reaction coordinate */
2191 fp=xvgropen(fn,"Integrated autocorrelation times","z","IACT [ps]",opt->oenv);
2192 fprintf(fp,"@ s0 symbol 1\n@ s0 symbol size 0.5\n@ s0 line linestyle 0\n");
2193 fprintf(fp,"# WIN tau(gr1) tau(gr2) ...\n");
2194 for (i=0;i<nwins;i++)
2196 fprintf(fp,"# %3d ",i);
2197 for (ig=0;ig<window[i].nPull;ig++)
2198 fprintf(fp," %11g",window[i].tau[ig]);
2201 for (i=0;i<nwins;i++)
2202 for (ig=0;ig<window[i].nPull;ig++)
2203 fprintf(fp,"%8g %8g\n",window[i].pos[ig],window[i].tau[ig]);
2204 if (opt->sigSmoothIact > 0.0)
2206 printf("Smoothing autocorrelation times along reaction coordinate with Gaussian of sig = %g\n",
2207 opt->sigSmoothIact);
2208 /* smooth IACT along reaction coordinate and overwrite g=1+2tau */
2209 smoothIact(window,nwins,opt);
2210 fprintf(fp,"&\n@ s1 symbol 1\n@ s1 symbol size 0.5\n@ s1 line linestyle 0\n");
2211 fprintf(fp,"@ s1 symbol color 2\n");
2212 for (i=0;i<nwins;i++)
2213 for (ig=0;ig<window[i].nPull;ig++)
2214 fprintf(fp,"%8g %8g\n",window[i].pos[ig],window[i].tausmooth[ig]);
2217 printf("Wrote %s\n",fn);
2220 /* compute average and sigma of each umbrella window */
2221 void averageSigma(t_UmbrellaWindow *window,int nwins,t_UmbrellaOptions *opt)
2224 real av,sum2,sig,diff,*ztime,nSamplesIndep;
2226 for (i=0;i<nwins;i++)
2228 snew(window[i].aver, window[i].nPull);
2229 snew(window[i].sigma,window[i].nPull);
2231 ntot=window[i].Ntot[0];
2232 for (ig=0;ig<window[i].nPull;ig++)
2234 ztime=window[i].ztime[ig];
2235 for (k=0, av=0.; k<ntot; k++)
2238 for (k=0, sum2=0.; k<ntot; k++)
2243 sig=sqrt(sum2/ntot);
2244 window[i].aver[ig]=av;
2246 /* Note: This estimate for sigma is biased from the limited sampling.
2247 Correct sigma by n/(n-1) where n = number of independent
2248 samples. Only possible if IACT is known.
2252 nSamplesIndep=window[i].N[ig]/(window[i].tau[ig]/window[i].dt);
2253 window[i].sigma[ig]=sig * nSamplesIndep/(nSamplesIndep-1);
2256 window[i].sigma[ig]=sig;
2257 printf("win %d, aver = %f sig = %f\n",i,av,window[i].sigma[ig]);
2263 /* Use histograms to compute average force on pull group.
2264 In addition, compute the sigma of the histogram.
2266 void computeAverageForce(t_UmbrellaWindow *window,int nWindows,t_UmbrellaOptions *opt)
2268 int i,j,bins=opt->bins,k;
2269 double dz,min=opt->min,max=opt->max,displAv,displAv2,temp,distance,ztot,ztot_half,w,weight;
2276 /* Compute average displacement from histograms */
2277 for(j=0;j<nWindows;++j)
2279 snew(window[j].forceAv,window[j].nPull);
2280 for(k=0;k<window[j].nPull;++k)
2285 for(i=0;i<opt->bins;++i)
2287 temp=(1.0*i+0.5)*dz+min;
2288 distance = temp - window[j].pos[k];
2290 { /* in cyclic wham: */
2291 if (distance > ztot_half) /* |distance| < ztot_half */
2293 else if (distance < -ztot_half)
2296 w=window[j].Histo[k][i]/window[j].g[k];
2297 displAv += w*distance;
2298 displAv2 += w*sqr(distance);
2300 /* Are we near min or max? We are getting wron forces from the histgrams since
2301 the histigrams are zero outside [min,max). Therefore, assume that the position
2302 on the other side of the histomgram center is equally likely. */
2305 posmirrored=window[j].pos[k]-distance;
2306 if (posmirrored>=max || posmirrored<min)
2308 displAv += -w*distance;
2309 displAv2 += w*sqr(-distance);
2317 /* average force from average displacement */
2318 window[j].forceAv[k] = displAv*window[j].k[k];
2319 /* sigma from average square displacement */
2320 /* window[j].sigma [k] = sqrt(displAv2); */
2321 /* printf("Win %d, sigma = %f\n",j,sqrt(displAv2)); */
2326 /* Check if the complete reaction coordinate is covered by the histograms */
2327 void checkReactionCoordinateCovered(t_UmbrellaWindow *window,int nwins,
2328 t_UmbrellaOptions *opt)
2330 int i,ig,j,bins=opt->bins,bBoundary;
2331 real avcount=0,z,relcount,*count;
2332 snew(count,opt->bins);
2334 for(j=0;j<opt->bins;++j)
2336 for (i=0;i<nwins;i++){
2337 for (ig=0;ig<window[i].nPull;ig++)
2338 count[j]+=window[i].Histo[ig][j];
2340 avcount+=1.0*count[j];
2345 relcount=count[j]/avcount;
2346 z=(j+0.5)*opt->dz+opt->min;
2347 bBoundary=( j<bins/20 || (bins-j)>bins/20 );
2348 /* check for bins with no data */
2350 fprintf(stderr, "\nWARNING, no data point in bin %d (z=%g) !\n"
2351 "You may not get a reasonable profile. Check your histograms!\n",j,z);
2352 /* and check for poor sampling */
2353 else if (relcount<0.005 && !bBoundary)
2354 fprintf(stderr, "Warning, poor sampling bin %d (z=%g). Check your histograms!\n",j,z);
2360 void guessPotByIntegration(t_UmbrellaWindow *window,int nWindows,t_UmbrellaOptions *opt,
2363 int i,j,ig,bins=opt->bins,nHist,winmin,groupmin;
2364 double dz,min=opt->min,*pot,pos,hispos,dist,diff,fAv,distmin,*f;
2367 dz=(opt->max-min)/bins;
2369 printf("Getting initial potential by integration.\n");
2371 /* Compute average displacement from histograms */
2372 computeAverageForce(window,nWindows,opt);
2374 /* Get force for each bin from all histograms in this bin, or, alternatively,
2375 if no histograms are inside this bin, from the closest histogram */
2378 for(j=0;j<opt->bins;++j)
2380 pos=(1.0*j+0.5)*dz+min;
2385 for (i=0;i<nWindows;i++)
2387 for (ig=0;ig<window[i].nPull;ig++)
2389 hispos=window[i].pos[ig];
2390 dist=fabs(hispos-pos);
2391 /* average force within bin */
2395 fAv+=window[i].forceAv[ig];
2397 /* at the same time, rememer closest histogram */
2405 /* if no histogram found in this bin, use closest histogram */
2409 fAv=window[winmin].forceAv[groupmin];
2413 for(j=1;j<opt->bins;++j)
2414 pot[j] = pot[j-1] - 0.5*dz*(f[j-1]+f[j]);
2416 /* cyclic wham: linearly correct possible offset */
2419 diff=(pot[bins-1]-pot[0])/(bins-1);
2420 for(j=1;j<opt->bins;++j)
2425 fp=xvgropen("pmfintegrated.xvg","PMF from force integration","z","PMF [kJ/mol]",opt->oenv);
2426 for(j=0;j<opt->bins;++j)
2427 fprintf(fp,"%g %g\n",(j+0.5)*dz+opt->min,pot[j]);
2429 printf("verbose mode: wrote %s with PMF from interated forces\n","pmfintegrated.xvg");
2432 /* get initial z=exp(-F[i]/kT) from integrated potential, where F[i] denote the free
2433 energy offsets which are usually determined by wham
2434 First: turn pot into probabilities:
2436 for(j=0;j<opt->bins;++j)
2437 pot[j]=exp(-pot[j]/(8.314e-3*opt->Temperature));
2438 calc_z(pot,window,nWindows,opt,TRUE);
2445 int gmx_wham(int argc,char *argv[])
2447 const char *desc[] = {
2448 "This is an analysis program that implements the Weighted",
2449 "Histogram Analysis Method (WHAM). It is intended to analyze",
2450 "output files generated by umbrella sampling simulations to ",
2451 "compute a potential of mean force (PMF). [PAR] ",
2452 "At present, three input modes are supported.[BR]",
2453 "[TT]*[tt] With option [TT]-it[tt], the user provides a file which contains the",
2454 " file names of the umbrella simulation run-input files ([TT].tpr[tt] files),",
2455 " AND, with option [TT]-ix[tt], a file which contains file names of",
2456 " the pullx [TT]mdrun[tt] output files. The [TT].tpr[tt] and pullx files must",
2457 " be in corresponding order, i.e. the first [TT].tpr[tt] created the",
2458 " first pullx, etc.[BR]",
2459 "[TT]*[tt] Same as the previous input mode, except that the the user",
2460 " provides the pull force output file names ([TT]pullf.xvg[tt]) with option [TT]-if[tt].",
2461 " From the pull force the position in the umbrella potential is",
2462 " computed. This does not work with tabulated umbrella potentials.[BR]"
2463 "[TT]*[tt] With option [TT]-ip[tt], the user provides file names of (gzipped) [TT].pdo[tt] files, i.e.",
2464 " the GROMACS 3.3 umbrella output files. If you have some unusual"
2465 " reaction coordinate you may also generate your own [TT].pdo[tt] files and",
2466 " feed them with the [TT]-ip[tt] option into to [TT]g_wham[tt]. The [TT].pdo[tt] file header",
2467 " must be similar to the following:[PAR]",
2468 "[TT]# UMBRELLA 3.0[BR]",
2469 "# Component selection: 0 0 1[BR]",
2471 "# Ref. Group 'TestAtom'[BR]",
2472 "# Nr. of pull groups 2[BR]",
2473 "# Group 1 'GR1' Umb. Pos. 5.0 Umb. Cons. 1000.0[BR]",
2474 "# Group 2 'GR2' Umb. Pos. 2.0 Umb. Cons. 500.0[BR]",
2476 "The number of pull groups, umbrella positions, force constants, and names ",
2477 "may (of course) differ. Following the header, a time column and ",
2478 "a data column for each pull group follows (i.e. the displacement",
2479 "with respect to the umbrella center). Up to four pull groups are possible ",
2480 "per [TT].pdo[tt] file at present.[PAR]",
2481 "By default, the output files are[BR]",
2482 " [TT]-o[tt] PMF output file[BR]",
2483 " [TT]-hist[tt] Histograms output file[BR]",
2484 "Always check whether the histograms sufficiently overlap.[PAR]",
2485 "The umbrella potential is assumed to be harmonic and the force constants are ",
2486 "read from the [TT].tpr[tt] or [TT].pdo[tt] files. If a non-harmonic umbrella force was applied ",
2487 "a tabulated potential can be provided with [TT]-tab[tt].[PAR]",
2488 "WHAM OPTIONS[BR]------------[BR]",
2489 " [TT]-bins[tt] Number of bins used in analysis[BR]",
2490 " [TT]-temp[tt] Temperature in the simulations[BR]",
2491 " [TT]-tol[tt] Stop iteration if profile (probability) changed less than tolerance[BR]",
2492 " [TT]-auto[tt] Automatic determination of boundaries[BR]",
2493 " [TT]-min,-max[tt] Boundaries of the profile [BR]",
2494 "The data points that are used to compute the profile",
2495 "can be restricted with options [TT]-b[tt], [TT]-e[tt], and [TT]-dt[tt]. ",
2496 "Adjust [TT]-b[tt] to ensure sufficient equilibration in each ",
2497 "umbrella window.[PAR]",
2498 "With [TT]-log[tt] (default) the profile is written in energy units, otherwise ",
2499 "(with [TT]-nolog[tt]) as probability. The unit can be specified with [TT]-unit[tt]. ",
2500 "With energy output, the energy in the first bin is defined to be zero. ",
2501 "If you want the free energy at a different ",
2502 "position to be zero, set [TT]-zprof0[tt] (useful with bootstrapping, see below).[PAR]",
2503 "For cyclic or periodic reaction coordinates (dihedral angle, channel PMF",
2504 "without osmotic gradient), the option [TT]-cycl[tt] is useful. [TT]g_wham[tt] will make use of the ",
2505 "periodicity of the system and generate a periodic PMF. The first and the last bin of the",
2506 "reaction coordinate will assumed be be neighbors.[PAR]",
2507 "Option [TT]-sym[tt] symmetrizes the profile around z=0 before output, ",
2508 "which may be useful for, e.g. membranes.[PAR]",
2509 "AUTOCORRELATIONS[BR]----------------[BR]",
2510 "With [TT]-ac[tt], [TT]g_wham[tt] estimates the integrated autocorrelation ",
2511 "time (IACT) [GRK]tau[grk] for each umbrella window and weights the respective ",
2512 "window with 1/[1+2*[GRK]tau[grk]/dt]. The IACTs are written ",
2513 "to the file defined with [TT]-oiact[tt]. In verbose mode, all ",
2514 "autocorrelation functions (ACFs) are written to [TT]hist_autocorr.xvg[tt]. ",
2515 "Because the IACTs can be severely underestimated in case of limited ",
2516 "sampling, option [TT]-acsig[tt] allows one to smooth the IACTs along the ",
2517 "reaction coordinate with a Gaussian ([GRK]sigma[grk] provided with [TT]-acsig[tt], ",
2518 "see output in [TT]iact.xvg[tt]). Note that the IACTs are estimated by simple ",
2519 "integration of the ACFs while the ACFs are larger 0.05.",
2520 "If you prefer to compute the IACTs by a more sophisticated (but possibly ",
2521 "less robust) method such as fitting to a double exponential, you can ",
2522 "compute the IACTs with [TT]g_analyze[tt] and provide them to [TT]g_wham[tt] with the file ",
2523 "[TT]iact-in.dat[tt] (option [TT]-iiact[tt]), which should contain one line per ",
2524 "input file ([TT].pdo[tt] or pullx/f file) and one column per pull group in the respective file.[PAR]",
2525 "ERROR ANALYSIS[BR]--------------[BR]",
2526 "Statistical errors may be estimated with bootstrap analysis. Use it with care, ",
2527 "otherwise the statistical error may be substantially underestimated. ",
2528 "More background and examples for the bootstrap technique can be found in ",
2529 "Hub, de Groot and Van der Spoel, JCTC (2010) 6: 3713-3720.[BR]",
2530 "[TT]-nBootstrap[tt] defines the number of bootstraps (use, e.g., 100). ",
2531 "Four bootstrapping methods are supported and ",
2532 "selected with [TT]-bs-method[tt].[BR]",
2533 " (1) [TT]b-hist[tt] Default: complete histograms are considered as independent ",
2534 "data points, and the bootstrap is carried out by assigning random weights to the ",
2535 "histograms (\"Bayesian bootstrap\"). Note that each point along the reaction coordinate",
2536 "must be covered by multiple independent histograms (e.g. 10 histograms), otherwise the ",
2537 "statistical error is underestimated.[BR]",
2538 " (2) [TT]hist[tt] Complete histograms are considered as independent data points. ",
2539 "For each bootstrap, N histograms are randomly chosen from the N given histograms ",
2540 "(allowing duplication, i.e. sampling with replacement).",
2541 "To avoid gaps without data along the reaction coordinate blocks of histograms ",
2542 "([TT]-histbs-block[tt]) may be defined. In that case, the given histograms are ",
2543 "divided into blocks and only histograms within each block are mixed. Note that ",
2544 "the histograms within each block must be representative for all possible histograms, ",
2545 "otherwise the statistical error is underestimated.[BR]",
2546 " (3) [TT]traj[tt] The given histograms are used to generate new random trajectories,",
2547 "such that the generated data points are distributed according the given histograms ",
2548 "and properly autocorrelated. The autocorrelation time (ACT) for each window must be ",
2549 "known, so use [TT]-ac[tt] or provide the ACT with [TT]-iiact[tt]. If the ACT of all ",
2550 "windows are identical (and known), you can also provide them with [TT]-bs-tau[tt]. ",
2551 "Note that this method may severely underestimate the error in case of limited sampling, ",
2552 "that is if individual histograms do not represent the complete phase space at ",
2553 "the respective positions.[BR]",
2554 " (4) [TT]traj-gauss[tt] The same as method [TT]traj[tt], but the trajectories are ",
2555 "not bootstrapped from the umbrella histograms but from Gaussians with the average ",
2556 "and width of the umbrella histograms. That method yields similar error estimates ",
2557 "like method [TT]traj[tt].[PAR]"
2558 "Bootstrapping output:[BR]",
2559 " [TT]-bsres[tt] Average profile and standard deviations[BR]",
2560 " [TT]-bsprof[tt] All bootstrapping profiles[BR]",
2561 "With [TT]-vbs[tt] (verbose bootstrapping), the histograms of each bootstrap are written, ",
2562 "and, with bootstrap method [TT]traj[tt], the cumulative distribution functions of ",
2566 const char *en_unit[]={NULL,"kJ","kCal","kT",NULL};
2567 const char *en_unit_label[]={"","E (kJ mol\\S-1\\N)","E (kcal mol\\S-1\\N)","E (kT)",NULL};
2568 const char *en_bsMethod[]={ NULL,"b-hist", "hist", "traj", "traj-gauss", NULL };
2570 static t_UmbrellaOptions opt;
2573 { "-min", FALSE, etREAL, {&opt.min},
2574 "Minimum coordinate in profile"},
2575 { "-max", FALSE, etREAL, {&opt.max},
2576 "Maximum coordinate in profile"},
2577 { "-auto", FALSE, etBOOL, {&opt.bAuto},
2578 "Determine min and max automatically"},
2579 { "-bins",FALSE, etINT, {&opt.bins},
2580 "Number of bins in profile"},
2581 { "-temp", FALSE, etREAL, {&opt.Temperature},
2583 { "-tol", FALSE, etREAL, {&opt.Tolerance},
2585 { "-v", FALSE, etBOOL, {&opt.verbose},
2587 { "-b", FALSE, etREAL, {&opt.tmin},
2588 "First time to analyse (ps)"},
2589 { "-e", FALSE, etREAL, {&opt.tmax},
2590 "Last time to analyse (ps)"},
2591 { "-dt", FALSE, etREAL, {&opt.dt},
2592 "Analyse only every dt ps"},
2593 { "-histonly", FALSE, etBOOL, {&opt.bHistOnly},
2594 "Write histograms and exit"},
2595 { "-boundsonly", FALSE, etBOOL, {&opt.bBoundsOnly},
2596 "Determine min and max and exit (with [TT]-auto[tt])"},
2597 { "-log", FALSE, etBOOL, {&opt.bLog},
2598 "Calculate the log of the profile before printing"},
2599 { "-unit", FALSE, etENUM, {en_unit},
2600 "Energy unit in case of log output" },
2601 { "-zprof0", FALSE, etREAL, {&opt.zProf0},
2602 "Define profile to 0.0 at this position (with [TT]-log[tt])"},
2603 { "-cycl", FALSE, etBOOL, {&opt.bCycl},
2604 "Create cyclic/periodic profile. Assumes min and max are the same point."},
2605 { "-sym", FALSE, etBOOL, {&opt.bSym},
2606 "Symmetrize profile around z=0"},
2607 { "-hist-eq", FALSE, etBOOL, {&opt.bHistEq},
2608 "HIDDENEnforce equal weight for all histograms. (Non-Weighed-HAM)"},
2609 { "-ac", FALSE, etBOOL, {&opt.bCalcTauInt},
2610 "Calculate integrated autocorrelation times and use in wham"},
2611 { "-acsig", FALSE, etREAL, {&opt.sigSmoothIact},
2612 "Smooth autocorrelation times along reaction coordinate with Gaussian of this [GRK]sigma[grk]"},
2613 { "-ac-trestart", FALSE, etREAL, {&opt.acTrestart},
2614 "When computing autocorrelation functions, restart computing every .. (ps)"},
2615 { "-acred", FALSE, etBOOL, {&opt.bAllowReduceIact},
2616 "HIDDENWhen smoothing the ACTs, allow to reduce ACTs. Otherwise, only increase ACTs "
2617 "during smoothing"},
2618 { "-nBootstrap", FALSE, etINT, {&opt.nBootStrap},
2619 "nr of bootstraps to estimate statistical uncertainty (e.g., 200)" },
2620 { "-bs-method", FALSE, etENUM, {en_bsMethod},
2621 "Bootstrap method" },
2622 { "-bs-tau", FALSE, etREAL, {&opt.tauBootStrap},
2623 "Autocorrelation time (ACT) assumed for all histograms. Use option [TT]-ac[tt] if ACT is unknown."},
2624 { "-bs-seed", FALSE, etINT, {&opt.bsSeed},
2625 "Seed for bootstrapping. (-1 = use time)"},
2626 { "-histbs-block", FALSE, etINT, {&opt.histBootStrapBlockLength},
2627 "When mixing histograms only mix within blocks of [TT]-histbs-block[tt]."},
2628 { "-vbs", FALSE, etBOOL, {&opt.bs_verbose},
2629 "Verbose bootstrapping. Print the CDFs and a histogram file for each bootstrap."},
2630 { "-stepout", FALSE, etINT, {&opt.stepchange},
2631 "HIDDENWrite maximum change every ... (set to 1 with [TT]-v[tt])"},
2632 { "-updateContr", FALSE, etINT, {&opt.stepUpdateContrib},
2633 "HIDDENUpdate table with significan contributions to WHAM every ... iterations"},
2637 { efDAT, "-ix","pullx-files",ffOPTRD}, /* wham input: pullf.xvg's and tprs */
2638 { efDAT, "-if","pullf-files",ffOPTRD}, /* wham input: pullf.xvg's and tprs */
2639 { efDAT, "-it","tpr-files",ffOPTRD}, /* wham input: tprs */
2640 { efDAT, "-ip","pdo-files",ffOPTRD}, /* wham input: pdo files (gmx3 style) */
2641 { efXVG, "-o", "profile", ffWRITE }, /* output file for profile */
2642 { efXVG, "-hist","histo", ffWRITE}, /* output file for histograms */
2643 { efXVG, "-oiact","iact",ffOPTWR}, /* writing integrated autocorrelation times */
2644 { efDAT, "-iiact","iact-in",ffOPTRD}, /* reading integrated autocorrelation times */
2645 { efXVG, "-bsres","bsResult", ffOPTWR}, /* average and errors of bootstrap analysis */
2646 { efXVG, "-bsprof","bsProfs", ffOPTWR}, /* output file for bootstrap profiles */
2647 { efDAT, "-tab","umb-pot",ffOPTRD}, /* Tabulated umbrella potential (if not harmonic) */
2650 int i,j,l,nfiles,nwins,nfiles2;
2651 t_UmbrellaHeader header;
2652 t_UmbrellaWindow * window=NULL;
2653 double *profile,maxchange=1e20;
2654 gmx_bool bMinSet,bMaxSet,bAutoSet,bExact=FALSE;
2655 char **fninTpr,**fninPull,**fninPdo;
2657 FILE *histout,*profout;
2658 char ylabel[256],title[256];
2660 #define NFILE asize(fnm)
2662 CopyRight(stderr,argv[0]);
2666 opt.bHistOnly=FALSE;
2675 /* bootstrapping stuff */
2677 opt.bsMethod=bsMethod_hist;
2678 opt.tauBootStrap=0.0;
2680 opt.histBootStrapBlockLength=8;
2681 opt.bs_verbose=FALSE;
2686 opt.Temperature=298;
2688 opt.bBoundsOnly=FALSE;
2690 opt.bCalcTauInt=FALSE;
2691 opt.sigSmoothIact=0.0;
2692 opt.bAllowReduceIact=TRUE;
2693 opt.bInitPotByIntegration=TRUE;
2696 opt.stepUpdateContrib=100;
2698 parse_common_args(&argc,argv,PCA_BE_NICE,
2699 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&opt.oenv);
2701 opt.unit=nenum(en_unit);
2702 opt.bsMethod=nenum(en_bsMethod);
2704 opt.bProf0Set=opt2parg_bSet("-zprof0", asize(pa), pa);
2706 opt.bTab=opt2bSet("-tab",NFILE,fnm);
2707 opt.bPdo=opt2bSet("-ip",NFILE,fnm);
2708 opt.bTpr=opt2bSet("-it",NFILE,fnm);
2709 opt.bPullx=opt2bSet("-ix",NFILE,fnm);
2710 opt.bPullf=opt2bSet("-if",NFILE,fnm);
2711 opt.bTauIntGiven=opt2bSet("-iiact",NFILE,fnm);
2712 if (opt.bTab && opt.bPullf)
2713 gmx_fatal(FARGS,"Force input does not work with tabulated potentials. "
2714 "Provide pullx.xvg or pdo files!");
2716 #define WHAMBOOLXOR(a,b) ( ((!(a))&&(b)) || ((a)&&(!(b))))
2717 if (!opt.bPdo && !WHAMBOOLXOR(opt.bPullx,opt.bPullf))
2718 gmx_fatal(FARGS,"Give either pullx (-ix) OR pullf (-if) data. Not both.");
2719 if ( !opt.bPdo && !(opt.bTpr || opt.bPullf || opt.bPullx))
2720 gmx_fatal(FARGS,"g_wham supports three input modes, pullx, pullf, or pdo file input."
2721 "\n\n Check g_wham -h !");
2723 opt.fnPdo=opt2fn("-ip",NFILE,fnm);
2724 opt.fnTpr=opt2fn("-it",NFILE,fnm);
2725 opt.fnPullf=opt2fn("-if",NFILE,fnm);
2726 opt.fnPullx=opt2fn("-ix",NFILE,fnm);
2728 bMinSet = opt2parg_bSet("-min", asize(pa), pa);
2729 bMaxSet = opt2parg_bSet("-max", asize(pa), pa);
2730 bAutoSet = opt2parg_bSet("-auto", asize(pa), pa);
2731 if ( (bMinSet || bMaxSet) && bAutoSet)
2732 gmx_fatal(FARGS,"With -auto, do not give -min or -max\n");
2734 if ( (bMinSet && !bMaxSet) || (!bMinSet && bMaxSet))
2735 gmx_fatal(FARGS,"When giving -min, you must give -max (and vice versa), too\n");
2737 if (bMinSet && opt.bAuto)
2739 printf("Note: min and max given, switching off -auto.\n");
2743 if (opt.bTauIntGiven && opt.bCalcTauInt)
2744 gmx_fatal(FARGS,"Either read (option -iiact) or calculate (option -ac) the\n"
2745 "the autocorrelation times. Not both.");
2747 if (opt.tauBootStrap>0.0 && opt2parg_bSet("-ac",asize(pa), pa))
2748 gmx_fatal(FARGS,"Either compute autocorrelation times (ACTs) (option -ac) or "
2749 "provide it with -bs-tau for bootstrapping. Not Both.\n");
2750 if (opt.tauBootStrap>0.0 && opt2bSet("-iiact",NFILE,fnm))
2751 gmx_fatal(FARGS,"Either provide autocorrelation times (ACTs) with file iact-in.dat "
2752 "(option -iiact) or define all ACTs with -bs-tau for bootstrapping\n. Not Both.");
2755 /* Reading gmx4 pull output and tpr files */
2756 if (opt.bTpr || opt.bPullf || opt.bPullx)
2758 read_wham_in(opt.fnTpr,&fninTpr,&nfiles,&opt);
2760 fnPull=opt.bPullf ? opt.fnPullf : opt.fnPullx;
2761 read_wham_in(fnPull,&fninPull,&nfiles2,&opt);
2762 printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n",
2763 nfiles,nfiles2,opt.bPullf ? "force" : "position",opt.fnTpr,fnPull);
2764 if (nfiles!=nfiles2)
2765 gmx_fatal(FARGS,"Found %d file names in %s, but %d in %s\n",nfiles,
2766 opt.fnTpr,nfiles2,fnPull);
2767 window=initUmbrellaWindows(nfiles);
2768 read_tpr_pullxf_files(fninTpr,fninPull,nfiles, &header, window, &opt);
2771 { /* reading pdo files */
2772 read_wham_in(opt.fnPdo,&fninPdo,&nfiles,&opt);
2773 printf("Found %d pdo files in %s\n",nfiles,opt.fnPdo);
2774 window=initUmbrellaWindows(nfiles);
2775 read_pdo_files(fninPdo,nfiles, &header, window, &opt);
2779 /* enforce equal weight for all histograms? */
2781 enforceEqualWeights(window,nwins);
2783 /* write histograms */
2784 histout=xvgropen(opt2fn("-hist",NFILE,fnm),"Umbrella histograms",
2785 "z","count",opt.oenv);
2786 for(l=0;l<opt.bins;++l)
2788 fprintf(histout,"%e\t",(double)(l+0.5)/opt.bins*(opt.max-opt.min)+opt.min);
2789 for(i=0;i<nwins;++i)
2791 for(j=0;j<window[i].nPull;++j)
2793 fprintf(histout,"%e\t",window[i].Histo[j][l]);
2796 fprintf(histout,"\n");
2799 printf("Wrote %s\n",opt2fn("-hist",NFILE,fnm));
2802 printf("Wrote histograms to %s, now exiting.\n",opt2fn("-hist",NFILE,fnm));
2806 /* Using tabulated umbrella potential */
2808 setup_tab(opt2fn("-tab",NFILE,fnm),&opt);
2810 /* Integrated autocorrelation times provided ? */
2811 if (opt.bTauIntGiven)
2812 readIntegratedAutocorrelationTimes(window,nwins,&opt,opt2fn("-iiact",NFILE,fnm));
2814 /* Compute integrated autocorrelation times */
2815 if (opt.bCalcTauInt)
2816 calcIntegratedAutocorrelationTimes(window,nwins,&opt,opt2fn("-oiact",NFILE,fnm));
2818 /* calc average and sigma for each histogram
2819 (maybe required for bootstrapping. If not, this is fast anyhow) */
2820 if (opt.nBootStrap && opt.bsMethod==bsMethod_trajGauss)
2821 averageSigma(window,nwins,&opt);
2823 /* Get initial potential by simple integration */
2824 if (opt.bInitPotByIntegration)
2825 guessPotByIntegration(window,nwins,&opt,0);
2827 /* Check if complete reaction coordinate is covered */
2828 checkReactionCoordinateCovered(window,nwins,&opt);
2830 /* Calculate profile */
2831 snew(profile,opt.bins);
2837 if ( (i%opt.stepUpdateContrib) == 0)
2838 setup_acc_wham(profile,window,nwins,&opt);
2839 if (maxchange<opt.Tolerance)
2842 /* if (opt.verbose) */
2843 printf("Switched to exact iteration in iteration %d\n",i);
2845 calc_profile(profile,window,nwins,&opt,bExact);
2846 if (((i%opt.stepchange) == 0 || i==1) && !i==0)
2847 printf("\t%4d) Maximum change %e\n",i,maxchange);
2849 } while ( (maxchange=calc_z(profile, window, nwins, &opt,bExact)) > opt.Tolerance || !bExact);
2850 printf("Converged in %d iterations. Final maximum change %g\n",i,maxchange);
2852 /* calc error from Kumar's formula */
2853 /* Unclear how the error propagates along reaction coordinate, therefore
2855 /* calc_error_kumar(profile,window, nwins,&opt); */
2857 /* Write profile in energy units? */
2860 prof_normalization_and_unit(profile,&opt);
2861 strcpy(ylabel,en_unit_label[opt.unit]);
2862 strcpy(title,"Umbrella potential");
2866 strcpy(ylabel,"Density of states");
2867 strcpy(title,"Density of states");
2870 /* symmetrize profile around z=0? */
2872 symmetrizeProfile(profile,&opt);
2874 /* write profile or density of states */
2875 profout=xvgropen(opt2fn("-o",NFILE,fnm),title,"z",ylabel,opt.oenv);
2876 for(i=0;i<opt.bins;++i)
2877 fprintf(profout,"%e\t%e\n",(double)(i+0.5)/opt.bins*(opt.max-opt.min)+opt.min,profile[i]);
2879 printf("Wrote %s\n",opt2fn("-o",NFILE,fnm));
2881 /* Bootstrap Method */
2883 do_bootstrapping(opt2fn("-bsres",NFILE,fnm),opt2fn("-bsprof",NFILE,fnm),
2884 opt2fn("-hist",NFILE,fnm),
2885 ylabel, profile, window, nwins, &opt);
2888 freeUmbrellaWindows(window,nfiles);
2890 printf("\nIn case you use results from g_wham for a publication, please cite:\n");
2891 please_cite(stdout,"Hub2010");