program wham
      implicit none
      integer*4 NH0,NINT0
      parameter (NH0=100,NINT0=1000)
      integer*4 NW,i,j,istart,iend,NNT,ii
      real*8 TEMP,dis1,dis2,dis3,pmin,vi,
     1roi(NH0,NINT0),wi(NH0,NINT0),fi(NH0),xi(NINT0),ro(NINT0),
     1an,beta,sum,rob,gg,fiold
      character*10 so
c
      TEMP=298.0d0
      NW=50
      NNT=50
    
      write(6,600)
600   format('     WHAM analysis for Tinker',/,
     1       '       Input: HAM.OPT')
      open(9,file='HAM.OPT')
1     read(9,900,end=99,err=99)so
900   format(a10)
      write(6,900)so
      if(so(1:4).eq.'NINT')read(9,*)NNT
      if(so(1:4).eq.'TEMP')read(9,*)TEMP
      if(so(1:2).eq.'NW'  )read(9,*)NW
      if(so(1:4).eq.'END')goto 99
      goto 1
99    close(9)
      write(6,*)'Options read'

c
c     List of options
c     ---------------
c     dynamic the dynamic command, with system path (e.g. /bin/dynamic)
c     minimize the minimize command, with system path
c     xmin lower coordinate limit
c     xmax upper coordinate limit
c     NEQUI number of equilibrium steps
c     NPROD number of production steps
c     NENS  ensemble type (2-NVT, 4-NPT)
c     NFS   MD integration step in fs
c     NPS   MD time between dumps in ps --- only one dump per simulation
c     NINT histogram integration step -- not used here
c     NW   number of weighted histograms
c     MTOL minimization gradient limit
c     PRES pressure (for NPT)
c     PPER periodic period
c     LPER use periodic coordinate
c     TEMP temperature
c
      if(NH0.lt.NW)call report('too many histograms')
      if(NINT0.lt.NNT)call report('too many points')
c
c     beta=1/(kT)
      beta=1.0d0/(TEMP*8.314d0/4.184d3)

c     read in all histograms
      do 2 i=1,NW
      write(so,'(i10)')i
      do 21 istart=1,10
21    if(so(istart:istart).ne.' ')goto 22
22    do 23 iend=10,1,-1
23    if(so(iend:iend).ne.' ')goto 24
24    open(9,file=so(istart:iend)//'.prn')
      do 2 j=1,NNT
      read(9,*)xi(j),ii,vi
      wi(i,j)=exp(-beta*vi)
2     roi(i,j)=dble(ii)
      close(9)
      write(6,*)NW,' histograms read'

c     initialize free energies   fi = exp(-Fi/kT)
      do 4 i=1,NW
4     fi(i)=1.0d0
c
c     this should probably make biased distribution functions:
      do 5 i=1,NW
      an=0.0d0
      do 6 j=1,NNT
6     an=an+roi(i,j)
      do 5 j=1,NNT
5     roi(i,j)=roi(i,j)/an
      write(6,*)'Histograms have been normalized'
      
c     density calculation
90000 do 10 j=1,NNT 
      ro(j)=0.0d0
      sum=0.0d0
      do 11 i=1,NW
11    sum=sum+wi(i,j)/fi(i)
      do 10 i=1,NW
c
c     biased density
      rob=roi(i,j)
10    if(sum.ne.0.0d0)ro(j)=ro(j)+rob/sum
c
c     ensure positive density
      do 14 j=1,NNT
14    ro(j)=max(1.0d-100,ro(j))

      gg=0.0d0
      do 8 i=1,NW
      fiold=fi(i)
      fi(i)=0.0d0
      do 9 j=1,NNT
9     fi(i)=fi(i)+wi(i,j)*ro(j)
8     gg=gg+(fiold-fi(i))**2/fiold**2

      write(6,*)gg
      if(gg.gt.1.0d-6)goto 90000

      write(6,*)'converged'

c     write down the PMF:
      pmin=-(1.0d0/beta)*log(ro(1))
      do 15 j=1,NNT
15    if(-(1.0d0/beta)*log(ro(j)).lt.pmin)pmin=-(1.0d0/beta)*log(ro(j))
      open(9,file='pmf.prn')
      do 12 j=1,NNT
12    write(9,1212)xi(j),-(1.0d0/beta)*log(ro(j))-pmin
1212  format(2g12.4)
      close(9)
      write(6,*)'PMF written into pmf.prn'

c     write down the density:
      open(9,file='ro.prn')
      do 16 j=1,NNT
16    write(9,1212)xi(j),ro(j)
      close(9)
      write(6,*)'Density written into ro.prn'

      end

      subroutine report(s)
      character*(*) s
      write(6,60)
      write(6,'(a)')s
      write(6,60)
60    format(70(1h*))
      stop
      end