! Weighed Histogram Analysis Method (WHAM) code
! Written by A. Liwo based on the work of Kumar et al., 
! J.Comput.Chem., 13, 1011 (1992)
      implicit none
      integer MaxN,MaxR,MaxIt
      parameter (MaxN=1000,MaxR=10,MaxIt=1000)
      double precision finorm_max,fimin
      parameter (finorm_max=1.0d0,fimin=1.0d-6)
      integer i,j,k,m,ind,iter,nR,t,tmin,tmax,n(MaxR),nint,nbin,nn
      integer start,end,iharm,ts,te
      double precision h(0:MaxN,MaxR),v(0:MaxN,MaxR),f(MaxR)
      double precision r0(MaxR),Kh(MaxR),snk(MaxR),
     & fi(MaxR),expf(MaxR),htot(0:MaxN),w(0:MaxN),
     & delta,dd,dd1,dd2,hh,beta,dmin,denom,finorm,stot,avefi,pom
      character*80 filein,fileout,statfile
      double precision pi,vol,d0,gnmr1,q
      character*8 tchar,qchar
      double precision tt0,qq0,slope1,slope2

      call getarg(1,tchar)
      call getarg(2,qchar)
      read (tchar,*) TT0
      read (qchar,*) qq0
 
      pi = 4.0d0*datan(1.0d0)
      

      DO while (beta.ne.TT0) 

      do i=1,MaxR
        do j=0,MaxN
          h(j,i)=0.0d0
        enddo
      enddo
! Read in the histograms from subsequent simulations.
      read (1,*) nR,delta,beta
      write (2,*) 'nr=',nr,' delta=',delta,' beta=',beta
      nint=1000*delta
      do i=1,nR
        do j=0,MaxN
          h(j,1)=0.0d0
        enddo
        tmax=0
        read (1,*) nbin,r0(1),Kh(1)
        write (2,*) 'r0=',r0(1),' kh=',kh(1)
        snk(1)=0
        do j=1,nbin
          read (1,*) q,nn 
          ind=int((q+1.0e-8)/delta)
          h(ind,1)=h(ind,1)+nn*1.0d0
          snk(1)=snk(1)+nn*1.0d0
          if (ind.gt.tmax) tmax=ind
        enddo ! j
        if (beta.eq.TT0 .and. r0(1).eq.qq0) exit
      enddo ! i

      ENDDO
 
      nR = 1
      beta=1.0d0/(beta*1.987D-3)
      write (2,*) 'beta=',beta,' nint=',nint

      do t=0,tmax
        htot(t)=0.0d0
        do k=1,nR
          htot(t)=htot(t)+h(t,k)
        enddo
      enddo

      stot=0.0d0
      do t=0,tmax
        stot=stot+htot(t)
      enddo

      write (2,*) 'tmax=',tmax
      do t=0,tmax
        write (2,'(i5,20f10.0)') t,(h(t,m),m=1,nr),htot(t)
      enddo
      write (2,'(/a5,20f10.0)') 'Total',(snk(m),m=1,nr),stot
      write (2,'(a)')

      do i=1,nR
!        f(i)=dlog(i+0.0d0)
        f(i)=0.0d0
      enddo

! Simple iteration to calculate free energies corresponding to all simulation
! runs.
      do iter=1,maxit 
        
        do k=1,nR
          expf(k)=dexp(f(k))
        enddo

! Compute the Boltzmann factor corresponing to restrain potentials in different
! simulations.
        pom=delta/nint
        do k=1,nR
          do t=0,tmax
             dd = dmin+delta*t
             v(t,k)=0.0d0
             do i=0,nint-1
               v(t,k)=v(t,k)+
     &            dexp(-beta*Kh(k)*(dd+pom*(i+0.5d0)-r0(k))**2)
             enddo
             v(t,k)=v(t,k)/nint
          enddo
        enddo
        
! Compute the weights of the bins using current free-energy values.
        do t=0,tmax
          denom=0.0d0
          do m=1,nR
            denom=denom+snk(m)*expf(m)*v(t,m)
          enddo
          w(t)=1.0d0/denom
        enddo
        write (2,'(i5,e20.4)') (t,w(t),t=0,tmax)

! Compute new free-energy values corresponding to the righ-hand side of the 
! equation and their derivatives.
        do i=1,nR
          fi(i)=0.0d0
          do t=0,tmax
            pom=htot(t)*w(t)*v(t,i)
            fi(i)=fi(i)+pom
          enddo ! t
        enddo ! i

        avefi=0.0d0
        do i=1,nR
          fi(i)=-dlog(fi(i))
          avefi=avefi+fi(i)
        enddo
        write (2,'(8f10.5)') (fi(i),i=1,nR)
        avefi=avefi/nR
        do i=1,nR
          fi(i)=fi(i)-avefi
        enddo
        write (2,'(8f10.5)') (fi(i),i=1,nR)
        write (2,'(8f10.5)') (f(i),i=1,nR)

! Compute the norm of free-energy increments.
        finorm=0.0d0
        do i=1,nR
          finorm=finorm+dabs(fi(i)-f(i))
          f(i)=fi(i)
        enddo  

        write (2,*) 'Iteration',iter,' finorm',finorm

! Exit, if the increment norm is smaller than pre-assigned tolerance.
        if (finorm.lt.fimin) then
          write (2,*) 'Iteration converged'
          goto 20
        endif

      enddo ! iter

   20 continue
! Handle middle zeros in htot; remove if they are next to first (last) bin,
! assign a small default otherwise
      tmin=0
      do while (htot(tmin).eq.0 .or. htot(tmin+1).eq.0) 
        tmin=tmin+1
      enddo
      do while (htot(tmax).eq.0 .or. htot(tmax-1).eq.0)
        tmax=tmax-1
      enddo
      t=tmin+2
      do while (t.lt.tmax)
        if (htot(t).eq.0) then
          write (2,*) "t=",t
          ts=t
          te=t
          do while (htot(te+1).eq.0) 
            te=te+1
          enddo
          slope1=(dlog(htot(ts-1))-dlog(htot(ts-2)))
          slope2=(dlog(htot(te+2))-dlog(htot(te+1)))
          write (2,*) "ts",ts," te",te," slope1",slope1,"slope2",slope2
          do t=ts,te
            if (htot(te+2).eq.0 .or. t-ts.lt.te-t) then
              htot(t)=dexp(dlog(htot(ts-1))+slope1*(t-ts+1))
              write (2,*) "htot<",htot(t)
            else if (t-ts.gt.te-t) then
              htot(t)=dexp(dlog(htot(te+1))+slope2*(t-te-1))
              write (2,*) "htot>",htot(t)
            else
              htot(t)=dexp((dlog(htot(ts-1))+slope1*(t-ts+1)+
     &         dlog(htot(te+1))+slope2*(t-te-1))/2)
              write (2,*) "htot=",htot(t)
            endif
          enddo
        endif
        t=t+1
      enddo
! Now, put together the histograms from all simulations, in order to get the
! unbiased total histogram.
      do t=tmin,tmax
c        htot(t)=0.0d0
c        do i=1,nR
c          htot(t)=htot(t)+h(t,i)
c        enddo
        write (2,*) "t, htot, w",t,htot(t),w(t)
        htot(t)=htot(t)*w(t)
      enddo
      
c      write (2,*) "tmin",tmin," tmax",tmax
c      do t=tmin+1,tmin+(tmax-tmin)/2
c        write (2,*) t,htot(t-1),htot(t),htot(t+1)
c        if (htot(t-1).gt.0 .and. htot(t).eq.0 .and. htot(t+1).gt.0) 
c     &     tmin=t+1
c        write (2,*) "tmin",tmin
c      enddo 
c      write (2,*) "tmin",tmin
c      do t=tmax-1,tmin+(tmax-tmin)/2,-1
c        write (2,*) t,htot(t-1),htot(t),htot(t+1)
c        if (htot(t-1).gt.0 .and. htot(t).eq.0 .and. htot(t+1).gt.0) 
c     &     tmax=t-1
c      enddo
c      write (2,*) "tmax",tmax
      write (2,'(a)') 'Final histogram'
      do t=tmin,tmax
        dd=dmin+(t+0.5d0)*delta
        write (2,'(f6.4,2e20.10)') dmin+t*delta,htot(t)
      enddo

      stop 111

      end
c----------------------------------------------------------------------------
      double precision function gnmr1(y,ymin,ymax)
      implicit none
      double precision y,ymin,ymax
      double precision wykl /4.0d0/
      if (y.lt.ymin) then
        gnmr1=(ymin-y)**wykl/wykl
      else if (y.gt.ymax) then
        gnmr1=(y-ymax)**wykl/wykl
      else
        gnmr1=0.0d0
      endif
      return
      end