3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
61 /* enum for energy units */
62 enum { enSel
, en_kJ
, en_kCal
, en_kT
, enNr
};
63 /* enum for type of input files (pdos, tpr, or pullf) */
64 enum { whamin_unknown
, whamin_tpr
, whamin_pullxf
, whamin_pdo
};
65 /* enum for methods to make profile cyclic/periodic */
66 enum { enCycl
, enCycl_no
, enCycl_yes
, enCycl_weighted
, enCycl_nr
};
70 /* umbrella with pull code of gromacs 4 */
71 int npullgrps
; /* nr of pull groups in tpr file */
72 int pull_geometry
; /* such as distance, position */
73 ivec pull_dim
; /* pull dimension with geometry distance */
74 int pull_ndim
; /* nr of pull_dim != 0 */
75 real
*k
; /* force constants in tpr file */
76 rvec
*init_dist
; /* reference displacements */
77 real
*umbInitDist
; /* referebce displacement in umbrella direction */
79 /* From here, old pdo stuff */
85 char PullName
[4][256];
107 const char *fnTpr
,*fnPullf
,*fnPdo
,*fnPullx
;
108 bool bTpr
,bPullf
,bPdo
,bPullx
;
110 bool verbose
,bShift
,bAuto
,bBoundsOnly
;
113 real Temperature
,Tolerance
;
114 int nBootStrap
,histBootStrapBlockLength
;
115 real dtBootStrap
,zProfZero
,alpha
;
116 int bsSeed
,stepchange
;
117 bool bHistBootStrap
,bWeightedCycl
,bHistOutOnly
;
118 bool bAutobounds
,bNoprof
;
123 bool bProf0Set
,bs_verbose
;
125 double *tabX
,*tabY
,tabMin
,tabMax
,tabDz
;
130 /* Return j such that xx[j] <= x < xx[j+1] */
131 void searchOrderedTable(double xx
[], int n
, double x
, int *j
)
139 ascending
=(xx
[n
-1] > xx
[0]);
143 if ((x
>= xx
[jm
]) == ascending
)
149 else if (x
==xx
[n
-1]) *j
=n
-2;
154 /* Read and setup tabulated umbrella potential */
155 void setup_tab(const char *fn
,t_UmbrellaOptions
*opt
)
161 printf("Setting up tabulated potential from file %s\n",fn
);
162 nl
=read_xvg(fn
,&y
,&ny
);
165 gmx_fatal(FARGS
,"Found %d columns in %s. Expected 2.\n",ny
,fn
);
167 opt
->tabMax
=y
[0][nl
-1];
168 opt
->tabDz
=(opt
->tabMax
-opt
->tabMin
)/(nl
-1);
170 gmx_fatal(FARGS
,"The tabulated potential in %s must be provided in \n"
171 "ascending z-direction",fn
);
173 if (fabs(y
[0][i
+1]-y
[0][i
]-opt
->tabDz
) > opt
->tabDz
/1e6
)
174 gmx_fatal(FARGS
,"z-values in %s are not equally spaced.\n",ny
,fn
);
178 opt
->tabX
[i
]=y
[0][i
];
179 opt
->tabY
[i
]=y
[1][i
];
181 printf("Found equally spaced tabulated potential from %g to %g, spacing %g\n",
182 opt
->tabMin
,opt
->tabMax
,opt
->tabDz
);
186 void read_pdo_header(FILE * file
,t_UmbrellaHeader
* header
,
187 t_UmbrellaOptions
*opt
)
198 if(3 != fscanf(file
,"%s%s%s",Buffer0
,Buffer1
,Buffer2
))
200 gmx_fatal(FARGS
,"Error reading header from pdo file");
202 if(strcmp(Buffer1
,"UMBRELLA"))
203 gmx_fatal(FARGS
,"This does not appear to be a valid pdo file. Found %s, expected %s",
204 Buffer1
, "UMBRELLA");
205 if(strcmp(Buffer2
,"3.0"))
206 gmx_fatal(FARGS
,"This does not appear to be a version 3.0 pdo file");
209 if(6 != fscanf(file
,"%s%s%s%d%d%d",Buffer0
,Buffer1
,Buffer2
,
210 &(header
->Dims
[0]),&(header
->Dims
[1]),&(header
->Dims
[2])))
212 gmx_fatal(FARGS
,"Error reading dimensions in header from pdo file");
215 /* printf("%d %d %d\n", header->Dims[0],header->Dims[1],header->Dims[2]); */
217 header
->nDim
= header
->Dims
[0] + header
->Dims
[1] + header
->Dims
[2];
219 gmx_fatal(FARGS
,"Currently only supports one dimension");
222 if(3 != fscanf(file
,"%s%s%d",Buffer0
,Buffer1
,&(header
->nSkip
)))
224 gmx_fatal(FARGS
,"Error reading header from pdo file");
228 if(4 != fscanf(file
,"%s%s%s%s",Buffer0
,Buffer1
,Buffer2
,header
->Reference
))
230 gmx_fatal(FARGS
,"Error reading header from pdo file");
234 if(6 != fscanf(file
,"%s%s%s%s%s%d",Buffer0
,Buffer1
,Buffer2
,Buffer3
,Buffer4
,&(header
->nPull
)))
236 gmx_fatal(FARGS
,"Error reading header from pdo file");
240 printf("Found nPull=%d , nSkip=%d, ref=%s\n",header
->nPull
,header
->nSkip
,
243 for(i
=0;i
<header
->nPull
;++i
)
245 if(4 != fscanf(file
,"%s%s%s%s",Buffer0
,Buffer1
,Buffer2
,header
->PullName
[i
]))
247 gmx_fatal(FARGS
,"Error reading header from pdo file");
250 printf("pullgroup %d, pullname = %s\n",i
,header
->PullName
[i
]);
251 for(j
=0;j
<header
->nDim
;++j
)
253 if(6 != fscanf(file
,"%s%s%lf%s%s%lf",Buffer0
,Buffer1
,&(header
->UmbPos
[i
][j
]),
254 Buffer2
,Buffer3
,&(header
->UmbCons
[i
][j
])))
256 gmx_fatal(FARGS
,"Error reading header from pdo file");
261 /* We want to combine both halves of a profile into one */
262 if(header
->UmbPos
[i
][j
]<0)
264 header
->UmbPos
[i
][j
]= -header
->UmbPos
[i
][j
];
265 header
->Flipped
[i
]=TRUE
;
268 else header
->Flipped
[i
]=FALSE
;
269 /*printf("%f\t%f\n",header->UmbPos[i][j],header->UmbCons[i][j]);*/
273 if(1 != fscanf(file
,"%s",Buffer3
))
275 gmx_fatal(FARGS
,"Error reading header from pdo file");
278 if (strcmp(Buffer3
,"#####") != 0)
279 gmx_fatal(FARGS
,"Expected '#####', found %s. Hick.\n",Buffer3
);
283 static char *fgets3(FILE *fp
,char ptr
[],int *len
)
289 if (fgets(ptr
,*len
-1,fp
) == NULL
)
292 while ((strchr(ptr
,'\n') == NULL
) && (!feof(fp
)))
294 /* This line is longer than len characters, let's increase len! */
298 if (fgets(p
-1,STRLEN
,fp
) == NULL
)
302 if (ptr
[slen
-1] == '\n')
309 void read_pdo_data(FILE * file
, t_UmbrellaHeader
* header
,
310 int fileno
, t_UmbrellaWindow
* win
,
311 t_UmbrellaOptions
*opt
,
312 bool bGetMinMax
,real
*mintmp
,real
*maxtmp
)
314 int i
,inttemp
,bins
,count
;
315 real min
,max
,minfound
,maxfound
;
316 double temp
,time
,time0
=0,dt
;
318 t_UmbrellaWindow
* window
=0;
320 char *tmpbuf
,fmt
[256],fmtign
[256];
321 int len
=STRLEN
,dstep
=1;
333 /* Need to alocate memory and set up structure */
334 window
->nPull
=header
->nPull
;
337 snew(window
->Histo
,window
->nPull
);
338 snew(window
->z
,window
->nPull
);
339 snew(window
->k
,window
->nPull
);
340 snew(window
->pos
,window
->nPull
);
341 snew(window
->Flipped
,window
->nPull
);
342 snew(window
->N
, window
->nPull
);
343 snew(window
->Ntot
, window
->nPull
);
345 for(i
=0;i
<window
->nPull
;++i
)
348 snew(window
->Histo
[i
],bins
);
349 window
->k
[i
]=header
->UmbCons
[i
][0];
350 window
->pos
[i
]=header
->UmbPos
[i
][0];
351 window
->Flipped
[i
]=header
->Flipped
[i
];
355 /* Done with setup */
361 min
=max
=bins
=0; /* Get rid of warnings */
366 while ( (ptr
=fgets3(file
,tmpbuf
,&len
)) != NULL
)
370 if (ptr
[0] == '#' || strlen(ptr
)<2)
373 /* Initiate format string */
375 strcat(fmtign
,"%*s");
377 sscanf(ptr
,"%lf",&time
); /* printf("Time %f\n",time); */
378 /* Round time to fs */
379 time
=1.0/1000*( (int) (time
*1000+0.5) );
381 /* get time step of pdo file */
389 dstep
=(int)(opt
->dt
/dt
+0.5);
398 dt_ok
=((count
-1)%dstep
== 0);
399 timeok
=(dt_ok
&& time
>= opt
->tmin
&& time
<= opt
->tmax
);
401 printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
402 time,opt->tmin, opt->tmax, dt_ok,timeok); */
406 for(i
=0;i
<header
->nPull
;++i
)
409 strcat(fmt
,"%lf"); /* Creating a format stings such as "%*s...%*s%lf" */
410 strcat(fmtign
,"%*s"); /* ignoring one more entry in the next loop */
411 if(sscanf(ptr
,fmt
,&temp
))
415 if(header
->Flipped
[i
]) temp
=-temp
;
418 temp
+=header
->UmbPos
[i
][0];
433 if (opt
->cycl
==enCycl_yes
)
437 else if (inttemp
>= bins
)
441 if(inttemp
>= 0 && inttemp
< bins
)
443 window
->Histo
[i
][inttemp
]+=1;
454 printf("time %f larger than tmax %f, stop reading pdo file\n",time
,opt
->tmax
);
467 void enforceEqualWeights(t_UmbrellaWindow
* window
,int nWindows
)
473 NEnforced
=window
[0].Ntot
[0];
474 printf("\nFound -hist-eq. Enforcing equal weights for all histograms, \ni.e. doing a "
475 "non-weighted histogram analysis method. Ndata = %d\n",NEnforced
);
476 /* enforce all histograms to have the same weight as the very first histogram */
478 for(j
=0;j
<nWindows
;++j
)
479 for(k
=0;k
<window
[j
].nPull
;++k
)
481 ratio
=1.0*NEnforced
/window
[j
].Ntot
[k
];
482 for(i
=0;i
<window
[0].nBin
;++i
)
483 window
[j
].Histo
[k
][i
]*=ratio
;
484 window
[j
].N
[k
]=(int)(ratio
*window
[j
].N
[k
]+0.5);
489 /* Simple linear interpolation between two given tabulated points */
490 double tabulated_pot(double dist
, t_UmbrellaOptions
*opt
)
496 jl
=floor((dist
-opt
->tabMin
)/opt
->tabDz
);
498 if (jl
<0 || ju
>=opt
->tabNbins
)
499 gmx_fatal(FARGS
,"Distance %f out of bounds of tabulated potential (jl=%d, ju=%d).\n"
500 "Provide an extended table.",dist
,jl
,ju
);
503 dz
=dist
-opt
->tabX
[jl
];
504 dp
=(pu
-pl
)*dz
/opt
->tabDz
;
509 /* Don't worry, that routine does not mean we compute the PMF in limited precision.
510 After rapid convergence (using only substiantal contributions), we always switch to
512 #define WHAM_CONTRIB_LIM 1e-10
513 void setup_acc_wham(t_UmbrellaWindow
* window
,int nWindows
, t_UmbrellaOptions
*opt
)
516 double U
,min
=opt
->min
,dz
=opt
->dz
,temp
,ztot_half
,distance
,ztot
,contrib
;
520 ztot
=opt
->max
-opt
->min
;
523 for(i
=0;i
<nWindows
;++i
)
525 snew(window
[i
].bContrib
,window
[i
].nPull
);
526 for(j
=0;j
<window
[i
].nPull
;++j
)
528 snew(window
[i
].bContrib
[j
],opt
->bins
);
530 for(k
=0;k
<opt
->bins
;++k
)
532 temp
=(1.0*k
+0.5)*dz
+min
;
533 distance
= temp
- window
[i
].pos
[j
]; /* distance to umbrella center */
534 if (opt
->cycl
==enCycl_yes
)
535 { /* in cyclic wham: */
536 if (distance
> ztot_half
) /* |distance| < ztot_half */
538 else if (distance
< -ztot_half
)
542 U
=0.5*window
[i
].k
[j
]*sqr(distance
); /* harmonic potential assumed. */
544 U
=tabulated_pot(distance
,opt
); /* Use tabulated potential */
546 contrib
=exp(- U
/(8.314e-3*opt
->Temperature
));
547 window
[i
].bContrib
[j
][k
] = (contrib
> WHAM_CONTRIB_LIM
);
548 bAnyContrib
= (bAnyContrib
| window
[i
].bContrib
[j
][k
]);
550 /* If this histo is far outside min and max all bContrib may be FALSE,
551 * causing a floating point exception later on. To avoid that, switch
552 * them all to true.*/
554 for(k
=0;k
<opt
->bins
;++k
)
555 window
[i
].bContrib
[j
][k
]=TRUE
;
558 printf("Initialized rapid wham stuff.\n");
562 void calc_profile(double *profile
,t_UmbrellaWindow
* window
, int nWindows
, t_UmbrellaOptions
*opt
,
566 double num
,ztot_half
,ztot
,distance
,min
=opt
->min
,dz
=opt
->dz
;
567 double denom
,U
=0,temp
=0;
570 ztot
=opt
->max
-opt
->min
;
574 for(i
=0;i
<opt
->bins
;++i
)
577 for(j
=0;j
<nWindows
;++j
)
579 for(k
=0;k
<window
[j
].nPull
;++k
)
581 temp
=(1.0*i
+0.5)*dz
+min
;
582 num
+=window
[j
].Histo
[k
][i
];
584 if (! (bExact
|| window
[j
].bContrib
[k
][i
]))
586 distance
= temp
- window
[j
].pos
[k
]; /* distance to umbrella center */
587 if (opt
->cycl
==enCycl_yes
){ /* in cyclic wham: */
588 if (distance
> ztot_half
) /* |distance| < ztot_half */
590 else if (distance
< -ztot_half
)
595 U
=0.5*window
[j
].k
[k
]*sqr(distance
); /* harmonic potential assumed. */
597 U
=tabulated_pot(distance
,opt
); /* Use tabulated potential */
598 denom
+=window
[j
].N
[k
]*exp(- U
/(8.314e-3*opt
->Temperature
) + window
[j
].z
[k
]);
601 profile
[i
]=num
/denom
;
606 double calc_z(double * profile
,t_UmbrellaWindow
* window
, int nWindows
, t_UmbrellaOptions
*opt
,
610 double U
=0,min
=opt
->min
,dz
=opt
->dz
,temp
,ztot_half
,distance
,ztot
;
611 double MAX
=-1e20
, total
=0;
614 ztot
=opt
->max
-opt
->min
;
617 for(i
=0;i
<nWindows
;++i
)
619 for(j
=0;j
<window
[i
].nPull
;++j
)
622 for(k
=0;k
<window
[i
].nBin
;++k
)
624 if (! (bExact
|| window
[i
].bContrib
[j
][k
]))
626 temp
=(1.0*k
+0.5)*dz
+min
;
627 distance
= temp
- window
[i
].pos
[j
]; /* distance to umbrella center */
628 if (opt
->cycl
==enCycl_yes
)
629 { /* in cyclic wham: */
630 if (distance
> ztot_half
) /* |distance| < ztot_half */
632 else if (distance
< -ztot_half
)
637 U
=0.5*window
[i
].k
[j
]*sqr(distance
); /* harmonic potential assumed. */
639 U
=tabulated_pot(distance
,opt
); /* Use tabulated potential */
641 total
+=profile
[k
]*exp(-U
/(8.314e-3*opt
->Temperature
));
647 temp
= fabs(total
- window
[i
].z
[j
]);
648 if(temp
> MAX
) MAX
=temp
;
649 window
[i
].z
[j
] = total
;
656 void cyclicProfByWeightedCorr(double *profile
,t_UmbrellaWindow
*window
,
657 int nWindows
, t_UmbrellaOptions
* opt
,
658 bool bAppendCorr2File
, const char *fn
,
659 const output_env_t oenv
)
661 int i
,j
,k
,bins
=opt
->bins
;
663 double *weight
,sum
=0.,diff
,*histsum
,*corr
,sumCorr
=0.,dCorr
;
669 printf("\nEnforcing a periodic profile by a sampling wheighted correction.");
670 please_cite(stdout
,"Hub2008");
677 /* generate weights proportional to 1/(n(i)*n(i+1))^alpha
678 where n(i) is the total nur of data points in bin i from all histograms */
679 for(i
=0;i
<nWindows
;++i
)
680 for(j
=0;j
<window
[i
].nPull
;++j
)
682 histsum
[k
]+=window
[i
].Histo
[j
][k
];
684 for(k
=0,sum
=0.;k
<bins
-1;++k
)
686 weight
[k
]=1./pow(histsum
[k
]*histsum
[k
+1],opt
->alpha
);
689 for(k
=0;k
<bins
-1;++k
)
692 /* difference between last and first bin */
693 diff
=profile
[bins
-1]-profile
[0];
694 printf("Distributing %f between adjacent bins to enforce a cyclic profile\n",diff
);
696 for (i
=0;i
<bins
-1;i
++)
698 dCorr
=weight
[i
]*diff
;
703 for (i
=0;i
<bins
-1;i
++)
704 profile
[i
+1]-=corr
[i
];
706 if (bAppendCorr2File
)
708 fp
=xvgropen(fn
,"Corrections to enforce periodicity","z",
709 "\\f{12}D\\f{}G(z)",oenv
);
710 sprintf(buf
,"corrections propotional to 1/(n\\si\\Nn\\si+1\\N)\\S%.2f",
712 xvgr_subtitle(fp
,buf
,oenv
);
713 for (i
=0;i
<bins
-1;i
++)
714 fprintf(fp
,"%g %g\n",opt
->min
+opt
->dz
*(i
+1),-corr
[i
]);
725 void prof_normalization_and_unit(double * profile
, t_UmbrellaOptions
*opt
)
728 double unit_factor
=1., R_MolarGasConst
, diff
;
731 R_MolarGasConst
=8.314472e-3; /* in kJ/(mol*K) */
734 /* No log? Nothing to do! */
738 /* Get profile in units of kT, kJ/mol, or kCal/mol */
739 if (opt
->unit
== en_kT
)
741 else if (opt
->unit
== en_kJ
)
742 unit_factor
=R_MolarGasConst
*opt
->Temperature
;
743 else if (opt
->unit
== en_kCal
)
744 unit_factor
=R_MolarGasConst
*opt
->Temperature
/4.1868;
746 gmx_fatal(FARGS
,"Sorry I don't know this energy unit.");
750 profile
[i
]=-log(profile
[i
])*unit_factor
;
752 /* shift to zero at z=opt->zProf0 */
756 /* Get bin with shortest distance to opt->zProf0 */
757 imin
=(int)((opt
->zProf0
-opt
->min
)/opt
->dz
);
771 void getRandomIntArray(int nPull
,int blockLength
,int* randomArray
)
773 int ipull
,blockBase
,nr
,ipullRandom
;
779 for (ipull
=0; ipull
<nPull
; ipull
++)
781 blockBase
=(ipull
/blockLength
)*blockLength
;
782 do{ /* make sure nothing bad happens in the last block */
783 nr
=(int)((1.0*rand()/RAND_MAX
)*blockLength
);
784 ipullRandom
= blockBase
+ nr
;
785 } while (ipullRandom
>= nPull
);
786 if (ipullRandom
<0 || ipullRandom
>=nPull
)
787 gmx_fatal(FARGS
,"Ups, random iWin = %d, nPull = %d, nr = %d, "
788 "blockLength = %d, blockBase = %d\n",
789 ipullRandom
,nPull
,nr
,blockLength
,blockBase
);
790 randomArray
[ipull
]=ipullRandom
;
792 /*for (ipull=0; ipull<nPull; ipull++)
793 printf("%d ",randomArray[ipull]); printf("\n"); */
797 void copy_pullgrp_to_synthwindow(t_UmbrellaWindow
*synthWindow
,
798 t_UmbrellaWindow
*thisWindow
,int pullid
)
800 synthWindow
->N
[0]=thisWindow
->N
[pullid
];
801 synthWindow
->Histo
[0]=thisWindow
->Histo
[pullid
];
802 synthWindow
->pos
[0]=thisWindow
->pos
[pullid
];
803 synthWindow
->z
[0]=thisWindow
->z
[pullid
];
804 synthWindow
->k
[0]=thisWindow
->k
[pullid
];
805 synthWindow
->bContrib
[0]=thisWindow
->bContrib
[pullid
];
809 /* Calculate cummulative of all histograms. They allow to create random numbers
810 which are distributed according to the histograms. Required to generate
811 the "synthetic" histograms for the Bootstrap method */
812 void calc_cummulants(t_UmbrellaWindow
*window
,int nWindows
,
813 t_UmbrellaOptions
*opt
,const char *fnhist
,
814 const output_env_t oenv
)
824 snew(fn
,strlen(fnhist
)+10);
825 snew(buf
,strlen(fnhist
)+10);
826 sprintf(fn
,"%s_cumul.xvg",strncpy(buf
,fnhist
,strlen(fnhist
)-4));
827 fp
=xvgropen(fn
,"Cummulants of umbrella histograms","z","cummulant",
832 for (i
=0; i
<nWindows
; i
++)
834 snew(window
[i
].cum
,window
[i
].nPull
);
835 for (j
=0; j
<window
[i
].nPull
; j
++)
837 snew(window
[i
].cum
[j
],nbin
+1);
838 window
[i
].cum
[j
][0]=0.;
839 for (k
=1; k
<=nbin
; k
++)
840 window
[i
].cum
[j
][k
] = window
[i
].cum
[j
][k
-1]+window
[i
].Histo
[j
][k
-1];
842 /* normalize cummulant. Ensure cum[nbin]==1 */
843 last
= window
[i
].cum
[j
][nbin
];
844 for (k
=0; k
<=nbin
; k
++)
845 window
[i
].cum
[j
][k
] /= last
;
849 printf("Cumulants of all histograms created.\n");
852 for (k
=0; k
<=nbin
; k
++)
854 fprintf(fp
,"%g\t",opt
->min
+k
*opt
->dz
);
855 for (i
=0; i
<nWindows
; i
++)
856 for (j
=0; j
<window
[i
].nPull
; j
++)
857 fprintf(fp
,"%g\t",window
[i
].cum
[j
][k
]);
860 printf("Wrote cumulants to %s\n",fn
);
868 /* Return j such that xx[j] <= x < xx[j+1] */
869 void searchCummulant(double xx
[], int n
, double x
, int *j
)
893 void create_synthetic_histo(t_UmbrellaWindow
*synthWindow
,
894 t_UmbrellaWindow
*thisWindow
,
895 int pullid
,t_UmbrellaOptions
*opt
)
897 int nsynth
,N
,i
,nbins
,r_index
;
899 static bool bWarnout
=0;
902 N
=thisWindow
->N
[pullid
];
903 nbins
=thisWindow
->nBin
;
905 /* nsynth = nr of data points in synthetic histo */
906 if (opt
->dtBootStrap
==0.0)
910 nsynth
=(int)(thisWindow
->N
[pullid
]*thisWindow
->dt
/opt
->dtBootStrap
+0.5);
915 if (!bWarnout
&& nsynth
<10)
917 printf("\n++++ WARNING ++++\n\tOnly %d data points per synthetic histogram!\n"
918 "\tYou may want to consider option -bs-dt.\n\n",nsynth
);
922 synthWindow
->N
[0]=nsynth
;
923 synthWindow
->pos
[0]=thisWindow
->pos
[pullid
];
924 synthWindow
->z
[0]=thisWindow
->z
[pullid
];
925 synthWindow
->k
[0]=thisWindow
->k
[pullid
];
926 synthWindow
->bContrib
[0]=thisWindow
->bContrib
[pullid
];
928 for (i
=0;i
<nbins
;i
++)
929 synthWindow
->Histo
[0][i
]=0.;
931 for (i
=0;i
<nsynth
;i
++)
933 r
=1.0*rand()/RAND_MAX
;
934 searchCummulant(thisWindow
->cum
[pullid
],nbins
+1 ,r
,&r_index
);
935 synthWindow
->Histo
[0][r_index
]+=1.;
940 void print_histograms(const char *fnhist
, t_UmbrellaWindow
* window
,
941 int nWindows
, int bs_index
,t_UmbrellaOptions
*opt
,
942 const output_env_t oenv
)
945 char *buf
=0,title
[256];
953 strcpy(title
,"Umbrella histograms");
957 snew(fn
,strlen(fnhist
)+6);
958 snew(buf
,strlen(fnhist
)+1);
959 sprintf(fn
,"%s_bs%d.xvg",strncpy(buf
,fnhist
,strlen(fnhist
)-4),bs_index
);
960 sprintf(title
,"Umbrella histograms. Bootstrap #%d",bs_index
);
963 fp
=xvgropen(fn
,title
,"z","count",oenv
);
966 /* Write histograms */
969 fprintf(fp
,"%e\t",(double)(l
+0.5)*opt
->dz
+opt
->min
);
970 for(i
=0;i
<nWindows
;++i
)
972 for(j
=0;j
<window
[i
].nPull
;++j
)
974 fprintf(fp
,"%e\t",window
[i
].Histo
[j
][l
]);
981 printf("Wrote %s\n",fn
);
990 void do_bootstrapping(const char *fnres
, const char* fnprof
,
991 const char *fnhist
, char* ylabel
, double *profile
,
992 t_UmbrellaWindow
* window
, int nWindows
,
993 t_UmbrellaOptions
*opt
, const output_env_t oenv
)
995 t_UmbrellaWindow
* synthWindow
;
996 double *bsProfile
,*bsProfiles_av
, *bsProfiles_av2
,maxchange
=1e20
,tmp
,stddev
;
997 int i
,j
,*randomArray
=0,winid
,pullid
,ib
;
998 int iAllPull
,nAllPull
,*allPull_winId
,*allPull_pullId
;
1004 if (opt
->bsSeed
==-1)
1009 snew(bsProfile
, opt
->bins
);
1010 snew(bsProfiles_av
, opt
->bins
);
1011 snew(bsProfiles_av2
,opt
->bins
);
1013 /* Create array of all pull groups. Note that different windows
1014 may have different nr of pull groups
1015 First: Get total nr of pull groups */
1017 for (i
=0;i
<nWindows
;i
++)
1018 nAllPull
+=window
[i
].nPull
;
1019 snew(allPull_winId
,nAllPull
);
1020 snew(allPull_pullId
,nAllPull
);
1022 /* Setup one array of all pull groups */
1023 for (i
=0;i
<nWindows
;i
++)
1024 for (j
=0;j
<window
[i
].nPull
;j
++)
1026 allPull_winId
[iAllPull
]=i
;
1027 allPull_pullId
[iAllPull
]=j
;
1031 /* setup stuff for synthetic windows */
1032 snew(synthWindow
,nAllPull
);
1033 for (i
=0;i
<nAllPull
;i
++)
1035 synthWindow
[i
].nPull
=1;
1036 synthWindow
[i
].nBin
=opt
->bins
;
1037 snew(synthWindow
[i
].Histo
,1);
1038 if (!opt
->bHistBootStrap
)
1039 snew(synthWindow
[i
].Histo
[0],opt
->bins
);
1040 snew(synthWindow
[i
].N
,1);
1041 snew(synthWindow
[i
].pos
,1);
1042 snew(synthWindow
[i
].z
,1);
1043 snew(synthWindow
[i
].k
,1);
1044 snew(synthWindow
[i
].bContrib
,1);
1047 if (opt
->bHistBootStrap
)
1049 snew(randomArray
,nAllPull
);
1050 printf("\n\nWhen computing statistical errors by bootstrapping entire histograms:\n");
1051 please_cite(stdout
,"Hub2006");
1055 calc_cummulants(window
,nWindows
,opt
,fnhist
,oenv
);
1058 /* do bootstrapping */
1059 fp
=xvgropen(fnprof
,"Boot strap profiles","z",ylabel
,oenv
);
1060 for (ib
=0;ib
<opt
->nBootStrap
;ib
++)
1062 printf(" *******************************************\n"
1063 " ******** Start bootstrap nr %d ************\n"
1064 " *******************************************\n",ib
+1);
1066 if (opt
->bHistBootStrap
)
1068 /* only mix given histos */
1069 getRandomIntArray(nAllPull
,opt
->histBootStrapBlockLength
,randomArray
);
1070 for (i
=0;i
<nAllPull
;i
++)
1072 winid
=allPull_winId
[randomArray
[i
]];
1073 pullid
=allPull_pullId
[randomArray
[i
]];
1074 copy_pullgrp_to_synthwindow(synthWindow
+i
,window
+winid
,pullid
);
1079 /* create new histos from given histos */
1080 for (i
=0;i
<nAllPull
;i
++)
1082 winid
=allPull_winId
[i
];
1083 pullid
=allPull_pullId
[i
];
1084 create_synthetic_histo(synthWindow
+i
,window
+winid
,pullid
,opt
);
1088 /* print histos in case of verbose output */
1089 if (opt
->bs_verbose
)
1090 print_histograms(fnhist
,synthWindow
,nAllPull
,ib
,opt
,oenv
);
1096 memcpy(bsProfile
,profile
,opt
->bins
*sizeof(double)); /* use profile as guess */
1098 if (maxchange
<opt
->Tolerance
)
1100 if (((i
%opt
->stepchange
) == 0 || i
==1) && !i
==0)
1101 printf("\t%4d) Maximum change %e\n",i
,maxchange
);
1102 calc_profile(bsProfile
,synthWindow
,nAllPull
,opt
,bExact
);
1104 } while( (maxchange
=calc_z(bsProfile
, synthWindow
, nAllPull
, opt
,bExact
)) > opt
->Tolerance
|| !bExact
);
1105 printf("\tConverged in %d iterations. Final maximum change %g\n",i
,maxchange
);
1108 prof_normalization_and_unit(bsProfile
,opt
);
1109 /* Force cyclic profile by wheighted correction */
1110 if (opt
->cycl
==enCycl_weighted
)
1111 cyclicProfByWeightedCorr(bsProfile
,synthWindow
,nAllPull
,opt
,
1114 /* save stuff to get average and stddev */
1115 for (i
=0;i
<opt
->bins
;i
++)
1118 bsProfiles_av
[i
]+=tmp
;
1119 bsProfiles_av2
[i
]+=tmp
*tmp
;
1120 fprintf(fp
,"%e\t%e\n",(i
+0.5)*opt
->dz
+opt
->min
,tmp
);
1126 /* write average and stddev */
1127 fp
=ffopen(fnres
,"w");
1128 fprintf(fp
,"@ title \"%s\"\n","Average and stddev from bootstrapping");
1129 fprintf(fp
,"@ xaxis label \"%s\"\n","z");
1130 fprintf(fp
,"@ yaxis label \"%s\"\n",ylabel
);
1131 fprintf(fp
,"@TYPE xydy\n");
1132 for (i
=0;i
<opt
->bins
;i
++)
1134 bsProfiles_av
[i
]/=opt
->nBootStrap
;
1135 bsProfiles_av2
[i
]/=opt
->nBootStrap
;
1136 tmp
=bsProfiles_av2
[i
]-sqr(bsProfiles_av
[i
]);
1137 stddev
=(tmp
>=0.) ? sqrt(tmp
) : 0.; /* Catch rouding errors */
1138 fprintf(fp
,"%e\t%e\t%e\n",(i
+0.5)*opt
->dz
+opt
->min
,bsProfiles_av
[i
],stddev
);
1141 printf("Wrote boot strap result to %s\n",fnres
);
1145 int whaminFileType(const char *fn
)
1149 if (strcmp(fn
+len
-3,"tpr")==0)
1151 else if (strcmp(fn
+len
-3,"xvg")==0 || strcmp(fn
+len
-6,"xvg.gz")==0)
1152 return whamin_pullxf
;
1153 else if (strcmp(fn
+len
-3,"pdo")==0 || strcmp(fn
+len
-6,"pdo.gz")==0)
1156 gmx_fatal(FARGS
,"Unknown file type of %s. Should be tpr, xvg, or pdo.\n",fn
);
1157 return whamin_unknown
;
1161 void read_wham_in(const char *fn
,char ***filenamesRet
, int *nfilesRet
,
1162 t_UmbrellaOptions
*opt
)
1164 char **filename
,tmp
[STRLEN
];
1165 int nread
,sizenow
,i
,block
=10;
1167 #define MAXFILELEN 512
1172 snew(filename
,sizenow
);
1173 for (i
=0;i
<sizenow
;i
++)
1174 snew(filename
[i
],MAXFILELEN
);
1176 while (fscanf(fp
,"%s",tmp
) != EOF
)
1178 if (strlen(tmp
)>=MAXFILELEN
)
1179 gmx_fatal(FARGS
,"Filename too long. Only %d characters allowed\n",MAXFILELEN
);
1180 strcpy(filename
[nread
],tmp
);
1182 printf("Found file %s in %s\n",filename
[nread
],fn
);
1187 srenew(filename
,sizenow
);
1188 for (i
=sizenow
-block
;i
<sizenow
;i
++)
1189 snew(filename
[i
],MAXFILELEN
);
1192 *filenamesRet
=filename
;
1197 FILE *pdo_open_file(const char *fn
)
1199 char Buffer
[1024],gunzip
[1024],*Path
=0;
1202 if (!gmx_fexist(fn
))
1204 gmx_fatal(FARGS
,"File %s does not exist.\n",fn
);
1207 /* gzipped pdo file? */
1208 if (strcmp(fn
+strlen(fn
)-3,".gz")==0)
1211 if(!(Path
=getenv("GMX_PATH_GZIP")))
1212 sprintf(gunzip
,"%s","/bin/gunzip");
1214 sprintf(gunzip
,"%s/gunzip",Path
);
1215 if (!gmx_fexist(gunzip
))
1216 gmx_fatal(FARGS
,"Cannot find executable %s. You may want to define the path to gunzip "
1217 "with the environment variable GMX_PATH_GZIP.",gunzip
);
1218 sprintf(Buffer
,"%s -c < %s",gunzip
,fn
);
1219 if((fp
=popen(Buffer
,"r"))==NULL
)
1221 gmx_fatal(FARGS
,"Unable to open pipe to `%s'\n",Buffer
);
1224 gmx_fatal(FARGS
,"Cannot open a compressed file on platform without pipe support");
1229 if((fp
=ffopen(fn
,"r"))==NULL
)
1231 gmx_fatal(FARGS
,"Unable to open file %s\n",fn
);
1238 pdo_close_file(FILE *fp
)
1247 /* Reading pdo files */
1248 void read_pdo_files(char **fn
, int nfiles
, t_UmbrellaHeader
* header
,
1249 t_UmbrellaWindow
**window
, t_UmbrellaOptions
*opt
)
1257 gmx_fatal(FARGS
,"No files found. Hick.");
1259 /* if min and max are not given, get min and max from the input files */
1262 printf("Automatic determination of boundaries from %d pdo files...\n",nfiles
);
1265 for(i
=0;i
<nfiles
;++i
)
1267 file
=pdo_open_file(fn
[i
]);
1268 printf("\rOpening %s ...",fn
[i
]); fflush(stdout
);
1271 read_pdo_header(file
,header
,opt
);
1272 /* here only determine min and max of this window */
1273 read_pdo_data(file
,header
,i
,NULL
,opt
,TRUE
,&mintmp
,&maxtmp
);
1274 if (maxtmp
>opt
->max
)
1276 if (mintmp
<opt
->min
)
1278 pdo_close_file(file
);
1281 printf("\nDetermined boundaries to %f and %f\n\n",opt
->min
,opt
->max
);
1282 if (opt
->bBoundsOnly
)
1284 printf("Found option -boundsonly, now exiting.\n");
1288 /* store stepsize in profile */
1289 opt
->dz
=(opt
->max
-opt
->min
)/opt
->bins
;
1291 snew(*window
,nfiles
);
1293 /* Having min and max, we read in all files */
1294 /* Loop over all files */
1295 for(i
=0;i
<nfiles
;++i
)
1297 printf("\rOpening %s ...",fn
[i
]); fflush(stdout
);
1300 file
=pdo_open_file(fn
[i
]);
1301 /* read in the headers */
1302 read_pdo_header(file
,header
,opt
);
1303 /* load data into window */
1304 read_pdo_data(file
,header
,i
,*window
,opt
,FALSE
,NULL
,NULL
);
1305 pdo_close_file(file
);
1311 #define int2YN(a) (((a)==0)?("N"):("Y"))
1313 void read_tpr_header(const char *fn
,t_UmbrellaHeader
* header
,
1314 t_UmbrellaOptions
*opt
)
1322 /* printf("Reading %s \n",fn); */
1323 read_tpx_state(fn
,&ir
,&state
,NULL
,NULL
);
1325 if (ir
.ePull
!= epullUMBRELLA
)
1326 gmx_fatal(FARGS
,"This is not a tpr of an umbrella simulation. Found ir.ePull = %s\n",
1327 epull_names
[ir
.ePull
]);
1329 /* nr of pull groups */
1332 gmx_fatal(FARGS
,"This is not a tpr of umbrella simulation. Found only %d pull groups\n",ngrp
);
1334 header
->npullgrps
=ir
.pull
->ngrp
;
1335 header
->pull_geometry
=ir
.pull
->eGeom
;
1336 copy_ivec(ir
.pull
->dim
,header
->pull_dim
);
1337 header
->pull_ndim
=header
->pull_dim
[0]+header
->pull_dim
[1]+header
->pull_dim
[2];
1338 if (header
->pull_geometry
==epullgPOS
&& header
->pull_ndim
>1)
1339 gmx_fatal(FARGS
,"Found pull geometry 'position' and more than 1 pull dimension (%d).\n"
1340 "Hence, the pull potential does not correspond to a one-dimensional umbrella potential.\n"
1341 "If you have some special umbrella setup you may want to write your own pdo files\n"
1342 "and feed them into g_wham. Check g_wham -h !\n",header
->pull_ndim
);
1343 snew(header
->k
,ngrp
);
1344 snew(header
->init_dist
,ngrp
);
1345 snew(header
->umbInitDist
,ngrp
);
1347 for (i
=0;i
<ngrp
;i
++)
1349 header
->k
[i
]=ir
.pull
->grp
[i
+1].k
;
1350 if (header
->k
[i
]==0.0)
1351 gmx_fatal(FARGS
,"Pull group %d has force constant of of 0.0 in %s.\n"
1352 "That doesn't seem to be an Umbrella tpr.\n",
1354 copy_rvec(ir
.pull
->grp
[i
+1].init
,header
->init_dist
[i
]);
1355 header
->Flipped
[i
]=opt
->bFlipProf
;
1357 /* initial distance to reference */
1358 switch(header
->pull_geometry
)
1362 if (header
->pull_dim
[d
])
1363 header
->umbInitDist
[i
]=header
->init_dist
[i
][d
];
1366 header
->umbInitDist
[i
]=header
->init_dist
[i
][0];
1369 gmx_fatal(FARGS
,"Pull geometry %s not supported\n",epullg_names
[header
->pull_geometry
]);
1373 if (opt
->verbose
|| first
)
1375 printf("File %s, %d groups, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n",
1376 fn
,header
->npullgrps
,epullg_names
[header
->pull_geometry
],
1377 int2YN(header
->pull_dim
[0]),int2YN(header
->pull_dim
[1]),int2YN(header
->pull_dim
[2]),
1379 for (i
=0;i
<ngrp
;i
++)
1380 printf("\tgrp %d) k = %.3f inittial distance = %g\n",i
,header
->k
[i
],header
->umbInitDist
[i
]);
1386 double dist_ndim(double **dx
,int ndim
,int line
)
1392 for (i
=0;i
<ndim
;i
++)
1393 r2
+=sqr(dx
[i
][line
]);
1399 void read_pull_xf(const char *fn
, const char *fntpr
,
1400 t_UmbrellaHeader
* header
, t_UmbrellaWindow
* window
,
1401 t_UmbrellaOptions
*opt
, bool bGetMinMax
,real
*mintmp
,
1404 double **y
,pos
=0.,t
,force
,time0
=0.,dt
;
1405 int ny
,nt
,bins
,ibin
,i
,g
,dstep
=1,nColPerGrp
,nColRef
,nColExpect
;
1406 real min
,max
,minfound
,maxfound
;
1407 bool dt_ok
,timeok
,bHaveForce
;
1408 const char *quantity
;
1414 in force output pullf.xvg: No reference, one column per pull group
1415 in position output pullx.xvg: ndim reference, ndim columns per pull group
1417 nColPerGrp
= opt
->bPullx
? header
->pull_ndim
: 1;
1418 nColRef
= opt
->bPullx
? header
->pull_ndim
: 0;
1419 quantity
= opt
->bPullx
? "position" : "force";
1420 nColExpect
= 1 + nColRef
+ header
->npullgrps
*nColPerGrp
;
1421 bHaveForce
= opt
->bPullf
;
1423 nt
=read_xvg(fn
,&y
,&ny
);
1425 /* Check consistency */
1427 gmx_fatal(FARGS
,"Empty pull %s file %s\n",quantity
,fn
);
1428 if (ny
!= nColExpect
)
1429 gmx_fatal(FARGS
,"Found %d pull groups in %s,\n but %d data columns in %s (expected %d)\n",
1430 header
->npullgrps
,fntpr
,ny
-1,fn
,nColExpect
-1);
1433 printf("Found %d times and %d %s sets %s\n",nt
,(ny
-1)/nColPerGrp
,quantity
,fn
);
1441 window
->dt
=y
[0][1]-y
[0][0];
1442 else if (opt
->nBootStrap
&& opt
->dtBootStrap
!=0.0)
1443 fprintf(stderr
,"\n *** WARNING, Could not determine time step in %s\n",fn
);
1445 /* Need to alocate memory and set up structure */
1446 window
->nPull
=header
->npullgrps
;
1449 snew(window
->Histo
,window
->nPull
);
1450 snew(window
->z
,window
->nPull
);
1451 snew(window
->k
,window
->nPull
);
1452 snew(window
->pos
,window
->nPull
);
1453 snew(window
->Flipped
,window
->nPull
);
1454 snew(window
->N
, window
->nPull
);
1455 snew(window
->Ntot
, window
->nPull
);
1457 for(g
=0;g
<window
->nPull
;++g
)
1460 snew(window
->Histo
[g
],bins
);
1461 window
->k
[g
]=header
->k
[g
];
1462 window
->Flipped
[g
]=header
->Flipped
[g
];
1465 window
->pos
[g
]=header
->umbInitDist
[g
];
1470 /* only determine min and max */
1473 min
=max
=bins
=0; /* Get rid of warnings */
1476 if(header
->Flipped
[0])
1477 gmx_fatal(FARGS
,"Sorry, flipping not supported for gmx4 output\n");
1481 /* Do you want that time frame? */
1482 t
=1.0/1000*((int)(0.5+y
[0][i
]*1000)); /* round time to fs */
1484 /* get time step of pdo file and get dstep from opt->dt */
1492 dstep
=(int)(opt
->dt
/dt
+0.5);
1497 window
->dt
=dt
*dstep
;
1500 dt_ok
=(i
%dstep
== 0);
1501 timeok
=(dt_ok
&& t
>= opt
->tmin
&& t
<= opt
->tmax
);
1503 printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
1504 t,opt->tmin, opt->tmax, dt_ok,timeok); */
1508 for(g
=0;g
<header
->npullgrps
;++g
)
1512 /* y has 1 time column y[0] and one column per force y[1],...,y[nGrps] */
1514 pos
= -force
/header
->k
[g
] + header
->umbInitDist
[g
];
1518 switch (header
->pull_geometry
)
1521 /* y has 1 time column y[0] and nColPerGrps columns per pull group;
1522 Distance to reference: */
1523 pos
=dist_ndim(y
+1+nColRef
+g
*nColPerGrp
,header
->pull_ndim
,i
);
1526 /* with geometry==position, we have always one column per group;
1527 Distance to reference: */
1528 pos
=y
[1+nColRef
+g
][i
];
1531 gmx_fatal(FARGS
,"Bad error, this error should have been catched before. Ups.\n");
1535 /* printf("grp %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
1545 ibin
=(int) floor((pos
-min
)/(max
-min
)*bins
);
1546 if (opt
->cycl
==enCycl_yes
)
1549 while ( (ibin
+=bins
) < 0);
1550 else if (ibin
>=bins
)
1551 while ( (ibin
-=bins
) >= bins
);
1553 if(ibin
>= 0 && ibin
< bins
)
1555 window
->Histo
[g
][ibin
]+=1.;
1562 else if (t
>opt
->tmax
)
1565 printf("time %f larger than tmax %f, stop reading this pullx/pullf file\n",t
,opt
->tmax
);
1578 void read_tpr_pullxf_files(char **fnTprs
,char **fnPull
,
1579 int nfiles
, t_UmbrellaHeader
* header
,
1580 t_UmbrellaWindow
**window
, t_UmbrellaOptions
*opt
)
1586 printf("Reading %d tpr and pullf files\n",nfiles
);
1588 /* min and max not given? */
1591 printf("Automatic determination of boundaries...\n");
1594 for (i
=0;i
<nfiles
; i
++)
1596 if (whaminFileType(fnTprs
[i
]) != whamin_tpr
)
1597 gmx_fatal(FARGS
,"Expected the %d'th file in input file to be a tpr file\n",i
);
1598 read_tpr_header(fnTprs
[i
],header
,opt
);
1599 if (whaminFileType(fnPull
[i
]) != whamin_pullxf
)
1600 gmx_fatal(FARGS
,"Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n",i
);
1601 read_pull_xf(fnPull
[i
],fnTprs
[i
],header
,NULL
,opt
,TRUE
,&mintmp
,&maxtmp
);
1602 if (maxtmp
>opt
->max
)
1604 if (mintmp
<opt
->min
)
1607 printf("\nDetermined boundaries to %f and %f\n\n",opt
->min
,opt
->max
);
1608 if (opt
->bBoundsOnly
){
1609 printf("Found option -boundsonly, now exiting.\n");
1613 /* store stepsize in profile */
1614 opt
->dz
=(opt
->max
-opt
->min
)/opt
->bins
;
1616 snew(*window
,nfiles
);
1617 for (i
=0;i
<nfiles
; i
++)
1619 if (whaminFileType(fnTprs
[i
]) != whamin_tpr
)
1620 gmx_fatal(FARGS
,"Expected the %d'th file in input file to be a tpr file\n",i
);
1621 read_tpr_header(fnTprs
[i
],header
,opt
);
1622 if (whaminFileType(fnPull
[i
]) != whamin_pullxf
)
1623 gmx_fatal(FARGS
,"Expected the %d'th file in input file to be a xvg (pullx/pullf) file\n",i
);
1624 read_pull_xf(fnPull
[i
],fnTprs
[i
],header
,*window
+i
,opt
,FALSE
,NULL
,NULL
);
1629 int gmx_wham(int argc
,char *argv
[])
1631 const char *desc
[] = {
1632 "This is an analysis program that implements the Weighted",
1633 "Histogram Analysis Method (WHAM). It is intended to analyze",
1634 "output files generated by umbrella sampling simulations to ",
1635 "compute a potential of mean force (PMF). [PAR]",
1636 "At present, three input modes are supported:[BR]",
1637 "[TT]*[tt] With option [TT]-it[tt], the user provides a file which contains the",
1638 " filenames of the umbrella simulation run-input files (tpr files),",
1639 " AND, with option -ix, a file which contains filenames of",
1640 " the pullx mdrun output files. The tpr and pullx files must",
1641 " be in corresponding order, i.e. the first tpr created the",
1642 " first pullx, etc.[BR]",
1643 "[TT]*[tt] Same as the previous input mode, except that the the user",
1644 " provides the pull force ouput file names (pullf.xvg) with option -if.",
1645 " From the pull force the position in the ubrella potential is",
1646 " computed. This does not work with tabulated umbrella potentials.",
1647 "[TT]*[tt] With option [TT]-ip[tt], the user provides filenames of (gzipped) pdo files, i.e.",
1648 " the gromacs 3.3 umbrella output files. If you have some unusual",
1649 " reaction coordinate you may also generate your own pdo files and",
1650 " feed them with the -ip option into to g_wham. The pdo file header",
1651 " must be similar to the folowing:[BR]",
1652 "[TT]# UMBRELLA 3.0[BR]",
1653 "# Component selection: 0 0 1[BR]",
1655 "# Ref. Group 'TestAtom'[BR]",
1656 "# Nr. of pull groups 2[BR]",
1657 "# Group 1 'GR1' Umb. Pos. 5.0 Umb. Cons. 1000.0[BR]",
1658 "# Group 2 'GR2' Umb. Pos. 2.0 Umb. Cons. 500.0[BR]",
1660 " Nr of pull groups, umbrella positions, force constants, and names",
1661 " may (of course) differ. Following the header, a time column and",
1662 " a data columns for each pull group follow (i.e. the displacement",
1663 " with respect to the umbrella center). Up to four pull groups are possible",
1664 " at present.[PAR]",
1665 "By default, the output files are[BR]",
1666 " [TT]-o[tt] PMF output file[BR]",
1667 " [TT]-hist[tt] histograms output file[PAR]",
1668 "The umbrella potential is assumed to be harmonic and the force constants are ",
1669 "read from the tpr or pdo files. If a non-harmonic umbrella force was applied ",
1670 "a tabulated potential can be provied with -tab.[PAR]",
1671 "WHAM OPTIONS[PAR]",
1672 " [TT]-bins[tt] Nr of bins used in analysis[BR]",
1673 " [TT]-temp[tt] Temperature in the simulations[BR]",
1674 " [TT]-tol[tt] Stop iteration if profile (probability) changed less than tolerance[BR]",
1675 " [TT]-auto[tt] Automatic determination of boudndaries[BR]",
1676 " [TT]-min,-max[tt] Boundaries of the profile [BR]",
1677 "The data points which are used ",
1678 "to compute the profile can be restricted with options -b, -e, and -dt. ",
1679 "Play particularly with -b to ensure sufficient equilibration in each ",
1680 "umbrella window![PAR]",
1681 "With -log (default) the profile is written in energy units, otherwise (-nolog) as ",
1682 "probability. The unit can be specified with -unit. With energy output, ",
1683 "the energy in the first bin is defined to be zero. If you want the free energy at a different ",
1684 "position to be zero, choose with -zprof0 (useful with bootstrapping, see below).[PAR]",
1685 "For cyclic (or periodic) reaction coordinates (dihedral angle, channel PMF",
1686 "without osmotic gradient), -cycl is useful.[BR]",
1687 "[TT]-cycl yes[tt] min and max are assumed to",
1688 "be neighboring points and histogram points outside min and max are mapped into ",
1689 "the interval [min,max] (compare histogram output). [BR]",
1690 "[TT]-cycl weighted[tt] First, a non-cyclic profile is computed. Subsequently, ",
1691 "periodicity is enforced by adding corrections dG(i) between neighboring bins",
1692 "i and i+1. The correction is chosen proportional to 1/[n(i)*n(i+1)]^alpha, where",
1693 "n(i) denotes the total nr of data points in bin i as collected from all histograms.",
1694 "alpha is defined with -alpha. The corrections are written to the file defined by -wcorr.",
1695 " (Compare Hub and de Groot, PNAS 105:1198 (2008))[PAR]",
1696 "ERROR ANALYSIS[BR]",
1697 "Statistical errors may be estimated with bootstrap analysis. Use it with care, ",
1698 "otherwise the statistical error may be substantially undererstimated !![BR]",
1699 "[TT]-nBootstrap[tt] defines the nr of bootstraps. Two bootstrapping modes are supported.[BR]",
1700 "[TT]-histbs[tt] Complete histograms are considered as independent data points (default). For each",
1701 "bootstrap, N histograms are randomly chosen from the N given histograms (allowing duplication).",
1702 "To avoid gaps without data along the reaction coordinate blocks of histograms (-histbs-block)",
1703 "may be defined. In that case, the given histograms are divided into blocks and ",
1704 "only histograms within each block are mixed. Note that the histograms",
1705 "within each block must be representative for all possible histograms, otherwise the",
1706 "statistical error is undererstimated![BR]",
1707 "[TT]-nohistbs[tt] The given histograms are used to generate new random histograms,",
1708 "such that the generated data points are distributed according the given histograms. The number",
1709 "of points generated for each bootstrap histogram can be controlled with -bs-dt.",
1710 "Note that one data point should be generated for each *independent* point in the given",
1711 "histograms. With the long autocorrelations in MD simulations, this procedure may ",
1712 "easily understimate the error![BR]",
1713 "Bootstrapping output:[BR]",
1714 "[TT]-bsres[tt] Average profile and standard deviations[BR]",
1715 "[TT]-bsprof[tt] All bootstrapping profiles[BR]",
1716 "With [TT]-vbs[tt] (verbose bootstrapping), the histograms of each bootstrap are written, and, ",
1717 "with [TT]-nohistBS[tt], the cummulants of the histogram.",
1720 static t_UmbrellaOptions opt
;
1721 static bool bHistOnly
=FALSE
;
1723 const char *en_unit
[]={NULL
,"kJ","kCal","kT",NULL
};
1724 const char *en_unit_label
[]={"","E (kJ mol\\S-1\\N)","E (kcal mol\\S-1\\N)","E (kT)",};
1725 const char *en_cycl
[]={NULL
,"no","yes","weighted",NULL
};
1728 { "-min", FALSE
, etREAL
, {&opt
.min
},
1729 "Minimum coordinate in profile"},
1730 { "-max", FALSE
, etREAL
, {&opt
.max
},
1731 "Maximum coordinate in profile"},
1732 { "-auto", FALSE
, etBOOL
, {&opt
.bAuto
},
1733 "determine min and max automatically"},
1734 { "-bins",FALSE
, etINT
, {&opt
.bins
},
1735 "Number of bins in profile"},
1736 { "-temp", FALSE
, etREAL
, {&opt
.Temperature
},
1738 { "-tol", FALSE
, etREAL
, {&opt
.Tolerance
},
1740 { "-v", FALSE
, etBOOL
, {&opt
.verbose
},
1742 { "-b", FALSE
, etREAL
, {&opt
.tmin
},
1743 "first time to analyse (ps)"},
1744 { "-e", FALSE
, etREAL
, {&opt
.tmax
},
1745 "last time to analyse (ps)"},
1746 { "-dt", FALSE
, etREAL
, {&opt
.dt
},
1747 "Analyse only every dt ps"},
1748 { "-histonly", FALSE
, etBOOL
, {&bHistOnly
},
1749 "Write histograms and exit"},
1750 { "-boundsonly", FALSE
, etBOOL
, {&opt
.bBoundsOnly
},
1751 "Determine min and max and exit (with -auto)"},
1752 { "-log", FALSE
, etBOOL
, {&opt
.bLog
},
1753 "Calculate the log of the profile before printing"},
1754 { "-unit", FALSE
, etENUM
, {en_unit
},
1755 "energy unit in case of log output" },
1756 { "-zprof0", FALSE
, etREAL
, {&opt
.zProf0
},
1757 "Define profile to 0.0 at this position (with -log)"},
1758 { "-cycl", FALSE
, etENUM
, {en_cycl
},
1759 "Create cyclic/periodic profile. Assumes min and max are the same point."},
1760 { "-alpha", FALSE
, etREAL
, {&opt
.alpha
},
1761 "for '-cycl weighted', set parameter alpha"},
1762 { "-flip", FALSE
, etBOOL
, {&opt
.bFlipProf
},
1763 "Combine halves of profile (not supported)"},
1764 { "-hist-eq", FALSE
, etBOOL
, {&opt
.bHistEq
},
1765 "Enforce equal weight for all histograms. (Non-Weighed-HAM)"},
1766 { "-nBootstrap", FALSE
, etINT
, {&opt
.nBootStrap
},
1767 "nr of bootstraps to estimate statistical uncertainty" },
1768 { "-bs-dt", FALSE
, etREAL
, {&opt
.dtBootStrap
},
1769 "timestep for synthetic bootstrap histograms (ps). Ensure independent data points!"},
1770 { "-bs-seed", FALSE
, etINT
, {&opt
.bsSeed
},
1771 "seed for bootstrapping. (-1 = use time)"},
1772 { "-histbs", FALSE
, etBOOL
, {&opt
.bHistBootStrap
},
1773 "In bootstrapping, consider complete histograms as one data point. "
1774 "Accounts better for long autocorrelations."},
1775 { "-histbs-block", FALSE
, etINT
, {&opt
.histBootStrapBlockLength
},
1776 "when mixin histograms only mix within blocks of -histBS_block."},
1777 { "-vbs", FALSE
, etBOOL
, {&opt
.bs_verbose
},
1778 "verbose bootstrapping. Print the cummulants and a histogram file for each bootstrap."},
1782 { efDAT
, "-ix","pullx-files",ffOPTRD
}, /* wham input: pullf.xvg's and tprs */
1783 { efDAT
, "-if","pullf-files",ffOPTRD
}, /* wham input: pullf.xvg's and tprs */
1784 { efDAT
, "-it","tpr-files",ffOPTRD
}, /* wham input: tprs */
1785 { efDAT
, "-ip","pdo-files",ffOPTRD
}, /* wham input: pdo files (gmx3 style) */
1786 { efXVG
, "-o", "profile", ffWRITE
}, /* output file for profile */
1787 { efXVG
, "-hist","histo", ffWRITE
}, /* output file for histograms */
1788 { efXVG
, "-bsres","bsResult", ffOPTWR
}, /* average and errors of bootstrap analysis */
1789 { efXVG
, "-bsprof","bsProfs", ffOPTWR
}, /* output file for bootstrap profiles */
1790 { efDAT
, "-tab","umb-pot",ffOPTRD
}, /* Tabulated umbrella potential (if not harmonic) */
1791 { efXVG
, "-wcorr","cycl-corr",ffOPTRD
}, /* Corrections to profile in case -cycl weighted */
1794 int i
,j
,l
,nfiles
,nwins
,nfiles2
;
1795 t_UmbrellaHeader header
;
1796 t_UmbrellaWindow
* window
=NULL
;
1797 double *profile
,maxchange
=1e20
;
1798 bool bMinSet
,bMaxSet
,bAutoSet
,bExact
=FALSE
;
1799 char **fninTpr
,**fninPull
,**fninPdo
;
1801 FILE *histout
,*profout
;
1802 char ylabel
[256],title
[256];
1813 opt
.dtBootStrap
=0.0;
1815 opt
.bHistBootStrap
=TRUE
;
1816 opt
.histBootStrapBlockLength
=12;
1818 opt
.bWeightedCycl
=FALSE
;
1820 opt
.bHistOutOnly
=FALSE
;
1828 opt
.bHistBootStrap
=TRUE
;
1829 opt
.histBootStrapBlockLength
=8;
1830 opt
.bs_verbose
=FALSE
;
1831 opt
.Temperature
=298;
1832 opt
.bFlipProf
=FALSE
;
1835 opt
.bBoundsOnly
=FALSE
;
1838 #define NFILE asize(fnm)
1840 CopyRight(stderr
,argv
[0]);
1841 parse_common_args(&argc
,argv
,PCA_BE_NICE
,
1842 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
,&oenv
);
1844 opt
.unit
=nenum(en_unit
);
1845 opt
.cycl
=nenum(en_cycl
);
1847 opt
.bProf0Set
=opt2parg_bSet("-zprof0", asize(pa
), pa
);
1849 opt
.bTab
=opt2bSet("-tab",NFILE
,fnm
);
1850 opt
.bPdo
=opt2bSet("-ip",NFILE
,fnm
);
1851 opt
.bTpr
=opt2bSet("-it",NFILE
,fnm
);
1852 opt
.bPullx
=opt2bSet("-ix",NFILE
,fnm
);
1853 opt
.bPullf
=opt2bSet("-if",NFILE
,fnm
);
1854 if (opt
.bTab
&& opt
.bPullf
)
1855 gmx_fatal(FARGS
,"Force input does not work with tabulated potentials. "
1856 "Provide pullx.xvg or pdo files!");
1858 #define BOOLXOR(a,b) ( ((!(a))&&(b)) || ((a)&&(!(b))))
1859 if (!opt
.bPdo
&& !BOOLXOR(opt
.bPullx
,opt
.bPullf
))
1860 gmx_fatal(FARGS
,"Give either pullx (-ix) OR pullf (-if) data. Not both.");
1861 if ( !opt
.bPdo
&& !(opt
.bTpr
|| opt
.bPullf
|| opt
.bPullx
))
1862 gmx_fatal(FARGS
,"g_wham supports three input modes, pullx, pullf, or pdo file input."
1863 "\n\n Check g_wham -h !");
1865 opt
.fnPdo
=opt2fn("-ip",NFILE
,fnm
);
1866 opt
.fnTpr
=opt2fn("-it",NFILE
,fnm
);
1867 opt
.fnPullf
=opt2fn("-if",NFILE
,fnm
);
1868 opt
.fnPullx
=opt2fn("-ix",NFILE
,fnm
);
1870 bMinSet
= opt2parg_bSet("-min", asize(pa
), pa
);
1871 bMaxSet
= opt2parg_bSet("-max", asize(pa
), pa
);
1872 bAutoSet
= opt2parg_bSet("-auto", asize(pa
), pa
);
1873 if ( (bMinSet
|| bMaxSet
) && bAutoSet
)
1874 gmx_fatal(FARGS
,"With -auto, do not give -min or -max\n");
1876 if ( (bMinSet
&& !bMaxSet
) || (!bMinSet
&& bMaxSet
))
1877 gmx_fatal(FARGS
,"When giving -min, you must give -max (and vice versa), too\n");
1879 if (bMinSet
&& opt
.bAuto
)
1881 printf("Note: min and max given, switching off -auto.\n");
1885 /* Reading gmx4 pull output and tpr files */
1886 if (opt
.bTpr
|| opt
.bPullf
|| opt
.bPullx
)
1888 read_wham_in(opt
.fnTpr
,&fninTpr
,&nfiles
,&opt
);
1890 fnPull
=opt
.bPullf
? opt
.fnPullf
: opt
.fnPullx
;
1891 read_wham_in(fnPull
,&fninPull
,&nfiles2
,&opt
);
1892 printf("Found %d tpr and %d pull %s files in %s and %s, respectively\n",
1893 nfiles
,nfiles2
,opt
.bPullf
? "force" : "position",opt
.fnTpr
,fnPull
);
1894 if (nfiles
!=nfiles2
)
1895 gmx_fatal(FARGS
,"Found %d file names in %s, but %d in %s\n",nfiles
,
1896 opt
.fnTpr
,nfiles2
,fnPull
);
1897 read_tpr_pullxf_files(fninTpr
,fninPull
,nfiles
, &header
, &window
, &opt
);
1901 /* reading pdo files */
1902 read_wham_in(opt
.fnPdo
,&fninPdo
,&nfiles
,&opt
);
1903 printf("Found %d pdo files in %s\n",nfiles
,opt
.fnPdo
);
1904 read_pdo_files(fninPdo
,nfiles
, &header
, &window
, &opt
);
1908 /* enforce equal weight for all histograms? */
1910 enforceEqualWeights(window
,nwins
);
1912 /* write histograms */
1913 histout
=xvgropen(opt2fn("-hist",NFILE
,fnm
),"Umbrella histograms",
1915 for(l
=0;l
<opt
.bins
;++l
)
1917 fprintf(histout
,"%e\t",(double)(l
+0.5)/opt
.bins
*(opt
.max
-opt
.min
)+opt
.min
);
1918 for(i
=0;i
<nwins
;++i
)
1920 for(j
=0;j
<window
[i
].nPull
;++j
)
1922 fprintf(histout
,"%e\t",window
[i
].Histo
[j
][l
]);
1925 fprintf(histout
,"\n");
1930 printf("Wrote histograms to %s, now exiting.\n",opt2fn("-hist",NFILE
,fnm
));
1934 /* Using tabulated umbrella potential */
1936 setup_tab(opt2fn("-tab",NFILE
,fnm
),&opt
);
1938 setup_acc_wham(window
,nwins
,&opt
);
1940 /* Calculate profile */
1941 snew(profile
,opt
.bins
);
1942 opt
.stepchange
=(opt
.verbose
)? 1 : 100;
1945 if (maxchange
<opt
.Tolerance
)
1948 /* if (opt.verbose) */
1949 printf("Switched to exact iteration in iteration %d\n",i
);
1951 calc_profile(profile
,window
,nwins
,&opt
,bExact
);
1952 if (((i
%opt
.stepchange
) == 0 || i
==1) && !i
==0)
1953 printf("\t%4d) Maximum change %e\n",i
,maxchange
);
1955 } while( (maxchange
=calc_z(profile
, window
, nwins
, &opt
,bExact
)) > opt
.Tolerance
|| !bExact
);
1956 printf("Converged in %d iterations. Final maximum change %g\n",i
,maxchange
);
1958 /* Write profile in energy units? */
1961 prof_normalization_and_unit(profile
,&opt
);
1962 strcpy(ylabel
,en_unit_label
[opt
.unit
]);
1963 strcpy(title
,"Umbrella potential");
1967 strcpy(ylabel
,"Density of states");
1968 strcpy(title
,"Density of states");
1970 /* Force cyclic profile by wheighted correction? */
1971 if (opt
.cycl
==enCycl_weighted
)
1972 cyclicProfByWeightedCorr(profile
,window
,nwins
,&opt
,TRUE
,
1973 opt2fn("-wcorr",NFILE
,fnm
),oenv
);
1975 profout
=xvgropen(opt2fn("-o",NFILE
,fnm
),title
,"z",ylabel
,oenv
);
1976 for(i
=0;i
<opt
.bins
;++i
)
1977 fprintf(profout
,"%e\t%e\n",
1978 (double)(i
+0.5)/opt
.bins
*(opt
.max
-opt
.min
)+opt
.min
,profile
[i
]);
1980 printf("Wrote %s\n",opt2fn("-o",NFILE
,fnm
));
1982 /* Bootstrap Method */
1984 do_bootstrapping(opt2fn("-bsres",NFILE
,fnm
),opt2fn("-bsprof",NFILE
,fnm
),
1985 opt2fn("-hist",NFILE
,fnm
),
1986 ylabel
, profile
, window
, nwins
, &opt
,oenv
);