1 ! Weighed Histogram Analysis Method (WHAM) code
2 ! Written by A. Liwo based on the work of Kumar et al.,
3 ! J.Comput.Chem., 13, 1011 (1992)
5 integer MaxN,MaxR,MaxIt
6 parameter (MaxN=1000,MaxR=10,MaxIt=1000)
7 double precision finorm_max,fimin
8 parameter (finorm_max=1.0d0,fimin=1.0d-6)
9 integer i,j,k,m,ind,iter,nR,t,tmin,tmax,n(MaxR),nint,nbin,nn
10 integer start,end,iharm,ts,te
11 double precision h(0:MaxN,MaxR),v(0:MaxN,MaxR),f(MaxR)
12 double precision r0(MaxR),Kh(MaxR),snk(MaxR),
13 & fi(MaxR),expf(MaxR),htot(0:MaxN),w(0:MaxN),
14 & delta,dd,dd1,dd2,hh,beta,dmin,denom,finorm,stot,avefi,pom
15 character*80 filein,fileout,statfile
16 double precision pi,vol,d0,gnmr1,q
17 character*8 tchar,qchar
18 double precision tt0,qq0,slope1,slope2
25 pi = 4.0d0*datan(1.0d0)
28 DO while (beta.ne.TT0)
35 ! Read in the histograms from subsequent simulations.
36 read (1,*) nR,delta,beta
37 write (2,*) 'nr=',nr,' delta=',delta,' beta=',beta
44 read (1,*) nbin,r0(1),Kh(1)
45 write (2,*) 'r0=',r0(1),' kh=',kh(1)
49 ind=int((q+1.0e-8)/delta)
50 h(ind,1)=h(ind,1)+nn*1.0d0
51 snk(1)=snk(1)+nn*1.0d0
52 if (ind.gt.tmax) tmax=ind
54 if (beta.eq.TT0 .and. r0(1).eq.qq0) exit
60 beta=1.0d0/(beta*1.987D-3)
61 write (2,*) 'beta=',beta,' nint=',nint
66 htot(t)=htot(t)+h(t,k)
75 write (2,*) 'tmax=',tmax
77 write (2,'(i5,20f10.0)') t,(h(t,m),m=1,nr),htot(t)
79 write (2,'(/a5,20f10.0)') 'Total',(snk(m),m=1,nr),stot
87 ! Simple iteration to calculate free energies corresponding to all simulation
95 ! Compute the Boltzmann factor corresponing to restrain potentials in different
104 & dexp(-beta*Kh(k)*(dd+pom*(i+0.5d0)-r0(k))**2)
110 ! Compute the weights of the bins using current free-energy values.
114 denom=denom+snk(m)*expf(m)*v(t,m)
118 write (2,'(i5,e20.4)') (t,w(t),t=0,tmax)
120 ! Compute new free-energy values corresponding to the righ-hand side of the
121 ! equation and their derivatives.
125 pom=htot(t)*w(t)*v(t,i)
135 write (2,'(8f10.5)') (fi(i),i=1,nR)
140 write (2,'(8f10.5)') (fi(i),i=1,nR)
141 write (2,'(8f10.5)') (f(i),i=1,nR)
143 ! Compute the norm of free-energy increments.
146 finorm=finorm+dabs(fi(i)-f(i))
150 write (2,*) 'Iteration',iter,' finorm',finorm
152 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
153 if (finorm.lt.fimin) then
154 write (2,*) 'Iteration converged'
161 ! Handle middle zeros in htot; remove if they are next to first (last) bin,
162 ! assign a small default otherwise
164 do while (htot(tmin).eq.0 .or. htot(tmin+1).eq.0)
167 do while (htot(tmax).eq.0 .or. htot(tmax-1).eq.0)
172 if (htot(t).eq.0) then
176 do while (htot(te+1).eq.0)
179 slope1=(dlog(htot(ts-1))-dlog(htot(ts-2)))
180 slope2=(dlog(htot(te+2))-dlog(htot(te+1)))
181 write (2,*) "ts",ts," te",te," slope1",slope1,"slope2",slope2
183 if (htot(te+2).eq.0 .or. t-ts.lt.te-t) then
184 htot(t)=dexp(dlog(htot(ts-1))+slope1*(t-ts+1))
185 write (2,*) "htot<",htot(t)
186 else if (t-ts.gt.te-t) then
187 htot(t)=dexp(dlog(htot(te+1))+slope2*(t-te-1))
188 write (2,*) "htot>",htot(t)
190 htot(t)=dexp((dlog(htot(ts-1))+slope1*(t-ts+1)+
191 & dlog(htot(te+1))+slope2*(t-te-1))/2)
192 write (2,*) "htot=",htot(t)
198 ! Now, put together the histograms from all simulations, in order to get the
199 ! unbiased total histogram.
203 c htot(t)=htot(t)+h(t,i)
205 write (2,*) "t, htot, w",t,htot(t),w(t)
209 c write (2,*) "tmin",tmin," tmax",tmax
210 c do t=tmin+1,tmin+(tmax-tmin)/2
211 c write (2,*) t,htot(t-1),htot(t),htot(t+1)
212 c if (htot(t-1).gt.0 .and. htot(t).eq.0 .and. htot(t+1).gt.0)
214 c write (2,*) "tmin",tmin
216 c write (2,*) "tmin",tmin
217 c do t=tmax-1,tmin+(tmax-tmin)/2,-1
218 c write (2,*) t,htot(t-1),htot(t),htot(t+1)
219 c if (htot(t-1).gt.0 .and. htot(t).eq.0 .and. htot(t+1).gt.0)
222 c write (2,*) "tmax",tmax
223 write (2,'(a)') 'Final histogram'
225 dd=dmin+(t+0.5d0)*delta
226 write (2,'(f6.4,2e20.10)') dmin+t*delta,htot(t)
232 c----------------------------------------------------------------------------
233 double precision function gnmr1(y,ymin,ymax)
235 double precision y,ymin,ymax
236 double precision wykl /4.0d0/
238 gnmr1=(ymin-y)**wykl/wykl
239 else if (y.gt.ymax) then
240 gnmr1=(y-ymax)**wykl/wykl