source/unres/src-HCD-5D/WHAM_q-single.f

   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)
   4       implicit none
   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
  19
  20       call getarg(1,tchar)
  21       call getarg(2,qchar)
  22       read (tchar,*) TT0
  23       read (qchar,*) qq0
  24
  25       pi = 4.0d0*datan(1.0d0)
  26
  27
  28       DO while (beta.ne.TT0)
  29
  30       do i=1,MaxR
  31         do j=0,MaxN
  32           h(j,i)=0.0d0
  33         enddo
  34       enddo
  35 ! Read in the histograms from subsequent simulations.
  36       read (1,*) nR,delta,beta
  37       write (2,*) 'nr=',nr,' delta=',delta,' beta=',beta
  38       nint=1000*delta
  39       do i=1,nR
  40         do j=0,MaxN
  41           h(j,1)=0.0d0
  42         enddo
  43         tmax=0
  44         read (1,*) nbin,r0(1),Kh(1)
  45         write (2,*) 'r0=',r0(1),' kh=',kh(1)
  46         snk(1)=0
  47         do j=1,nbin
  48           read (1,*) q,nn
  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
  53         enddo ! j
  54         if (beta.eq.TT0 .and. r0(1).eq.qq0) exit
  55       enddo ! i
  56
  57       ENDDO
  58
  59       nR = 1
  60       beta=1.0d0/(beta*1.987D-3)
  61       write (2,*) 'beta=',beta,' nint=',nint
  62
  63       do t=0,tmax
  64         htot(t)=0.0d0
  65         do k=1,nR
  66           htot(t)=htot(t)+h(t,k)
  67         enddo
  68       enddo
  69
  70       stot=0.0d0
  71       do t=0,tmax
  72         stot=stot+htot(t)
  73       enddo
  74
  75       write (2,*) 'tmax=',tmax
  76       do t=0,tmax
  77         write (2,'(i5,20f10.0)') t,(h(t,m),m=1,nr),htot(t)
  78       enddo
  79       write (2,'(/a5,20f10.0)') 'Total',(snk(m),m=1,nr),stot
  80       write (2,'(a)')
  81
  82       do i=1,nR
  83 !        f(i)=dlog(i+0.0d0)
  84         f(i)=0.0d0
  85       enddo
  86
  87 ! Simple iteration to calculate free energies corresponding to all simulation
  88 ! runs.
  89       do iter=1,maxit
  90
  91         do k=1,nR
  92           expf(k)=dexp(f(k))
  93         enddo
  94
  95 ! Compute the Boltzmann factor corresponing to restrain potentials in different
  96 ! simulations.
  97         pom=delta/nint
  98         do k=1,nR
  99           do t=0,tmax
 100              dd = dmin+delta*t
 101              v(t,k)=0.0d0
 102              do i=0,nint-1
 103                v(t,k)=v(t,k)+
 104      &            dexp(-beta*Kh(k)*(dd+pom*(i+0.5d0)-r0(k))**2)
 105              enddo
 106              v(t,k)=v(t,k)/nint
 107           enddo
 108         enddo
 109
 110 ! Compute the weights of the bins using current free-energy values.
 111         do t=0,tmax
 112           denom=0.0d0
 113           do m=1,nR
 114             denom=denom+snk(m)*expf(m)*v(t,m)
 115           enddo
 116           w(t)=1.0d0/denom
 117         enddo
 118         write (2,'(i5,e20.4)') (t,w(t),t=0,tmax)
 119
 120 ! Compute new free-energy values corresponding to the righ-hand side of the
 121 ! equation and their derivatives.
 122         do i=1,nR
 123           fi(i)=0.0d0
 124           do t=0,tmax
 125             pom=htot(t)*w(t)*v(t,i)
 126             fi(i)=fi(i)+pom
 127           enddo ! t
 128         enddo ! i
 129
 130         avefi=0.0d0
 131         do i=1,nR
 132           fi(i)=-dlog(fi(i))
 133           avefi=avefi+fi(i)
 134         enddo
 135         write (2,'(8f10.5)') (fi(i),i=1,nR)
 136         avefi=avefi/nR
 137         do i=1,nR
 138           fi(i)=fi(i)-avefi
 139         enddo
 140         write (2,'(8f10.5)') (fi(i),i=1,nR)
 141         write (2,'(8f10.5)') (f(i),i=1,nR)
 142
 143 ! Compute the norm of free-energy increments.
 144         finorm=0.0d0
 145         do i=1,nR
 146           finorm=finorm+dabs(fi(i)-f(i))
 147           f(i)=fi(i)
 148         enddo
 149
 150         write (2,*) 'Iteration',iter,' finorm',finorm
 151
 152 ! Exit, if the increment norm is smaller than pre-assigned tolerance.
 153         if (finorm.lt.fimin) then
 154           write (2,*) 'Iteration converged'
 155           goto 20
 156         endif
 157
 158       enddo ! iter
 159
 160    20 continue
 161 ! Handle middle zeros in htot; remove if they are next to first (last) bin,
 162 ! assign a small default otherwise
 163       tmin=0
 164       do while (htot(tmin).eq.0 .or. htot(tmin+1).eq.0)
 165         tmin=tmin+1
 166       enddo
 167       do while (htot(tmax).eq.0 .or. htot(tmax-1).eq.0)
 168         tmax=tmax-1
 169       enddo
 170       t=tmin+2
 171       do while (t.lt.tmax)
 172         if (htot(t).eq.0) then
 173           write (2,*) "t=",t
 174           ts=t
 175           te=t
 176           do while (htot(te+1).eq.0)
 177             te=te+1
 178           enddo
 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
 182           do t=ts,te
 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)
 189             else
 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)
 193             endif
 194           enddo
 195         endif
 196         t=t+1
 197       enddo
 198 ! Now, put together the histograms from all simulations, in order to get the
 199 ! unbiased total histogram.
 200       do t=tmin,tmax
 201 c        htot(t)=0.0d0
 202 c        do i=1,nR
 203 c          htot(t)=htot(t)+h(t,i)
 204 c        enddo
 205         write (2,*) "t, htot, w",t,htot(t),w(t)
 206         htot(t)=htot(t)*w(t)
 207       enddo
 208
 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)
 213 c     &     tmin=t+1
 214 c        write (2,*) "tmin",tmin
 215 c      enddo
 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)
 220 c     &     tmax=t-1
 221 c      enddo
 222 c      write (2,*) "tmax",tmax
 223       write (2,'(a)') 'Final histogram'
 224       do t=tmin,tmax
 225         dd=dmin+(t+0.5d0)*delta
 226         write (2,'(f6.4,2e20.10)') dmin+t*delta,htot(t)
 227       enddo
 228
 229       stop 111
 230
 231       end
 232 c----------------------------------------------------------------------------
 233       double precision function gnmr1(y,ymin,ymax)
 234       implicit none
 235       double precision y,ymin,ymax
 236       double precision wykl /4.0d0/
 237       if (y.lt.ymin) then
 238         gnmr1=(ymin-y)**wykl/wykl
 239       else if (y.gt.ymax) then
 240         gnmr1=(y-ymax)**wykl/wykl
 241       else
 242         gnmr1=0.0d0
 243       endif
 244       return
 245       end