program GaussMC

     implicit none
     
     integer iunit

     integer imove,iblock
     integer ntry,nacc
     integer, parameter:: nbins=100,nblock=10,nmove=1e5

     real*8 mu,sigma2,xnew,xold
     real*8 beta,potnew,potold,deltapot
     real*8 xranf,urnd,delta,acceptance
     real*8 mean,mean2,variance
 
     real*8 blmean,blmean2,blvariance


     integer bin
     real*8 minbin,maxbin,binwidth
     real*8 ghist(nbins),nghist
     real*8 gaussian(nbins),ngaussian

     logical ife

     open(17,file='acceptance.dat',status='unknown')
     open(18,file='position.xyz',status='unknown')
     open(19,file='avgs.dat',status='unknown')
     open(20,file='histogram.dat',status='unknown')
     open(21,file='blavgs.dat',status='unknown')


     iunit=0
     inquire(file='rmaseed,dat',exist=ife)
     if(ife) iunit=1
     call rmaset(iunit)


! Set Gaussians parameters, "inverse temperature", and MC displacement
     mu=1.
     sigma2=5.
     beta=1./sigma2
     delta=3.*sigma2

     minbin=mu-4.*sigma2
     maxbin=mu+4.*sigma2
     binwidth=(maxbin-minbin)/dble(nbins)
     
     
     
! Initialize position and monitoring parameters
     xold=mu
     xnew=xold

     ntry=0
     nacc=0
     acceptance=0.
     mean=0.
     mean2=0.
     variance=0.



     do iblock=1,nblock

        blmean=0.
        blmean2=0.
        blvariance=0.


! Compute "potential"

        potold=0.5*(xold-mu)**2
        potnew=potold

        ! Begin MC loop
        
        do imove=1,nmove
           
           ntry=ntry+1
           
! Generate trial position
           call rmar(urnd)
           xnew=xold+(urnd-0.5)*delta

! Compute new "potential"
           potnew = 0.5*(xnew-mu)**2

! Compute variation in the "potential"
           deltapot=potnew-potold

! Select random number for Metropolis test
           call rmar(urnd)

! Perform Metropolis test and update variables accordingly
           if(urnd.lt.exp(-beta*deltapot))then
              xold=xnew
              potold=potnew
              nacc=nacc+1
           end if

! Compute acceptance and running averages
           acceptance=dble(nacc)/dble(ntry)
           
           mean=(dble(ntry-1)*mean+xold)/dble(ntry)
           mean2=(dble(ntry-1)*mean2+xold**2)/dble(ntry)
           variance=dble(ntry)*(mean2-mean**2)/dble(max(1,ntry-1))
           

           bin=nint(1+(xold-minbin)/binwidth)
           if(bin.lt.1)bin=1
           if(bin.gt.nbins)bin=nbins
           ghist(bin)=ghist(bin)+1
           
           if(imove.eq.nmove)then
              write(17,222) ntry, acceptance
              write(18,555) xnew, 0., 0.
              write(19,222) imove, mean, mean2, variance
           end if
        
           
        end do

! Accumulate block averages
           blmean=(dble(nblock-1)*blmean+mean)/dble(nblock)
           blmean2=(dble(nblock-1)*blmean2+mean2)/dble(nblock)
           blvariance=dble(nblock)*(blmean2-blmean**2)/dble(max(1,nblock-1))

           write(21,222) iblock,blmean,blmean2,blvariance


     end do

! Create histogram of the reconstructed probability density
        do bin=1,nbins
           gaussian(bin)=exp(-0.5*dble((minbin+(bin-1)*binwidth)-mu)**2/sigma2)
        end do

        nghist=sum(ghist)*dble(binwidth)
        ngaussian=sum(gaussian)*dble(binwidth)

        do bin=1,nbins
           write(20,333) minbin+(bin-1)*binwidth, ghist(bin)/nghist,gaussian(bin)/ngaussian
        end do



222 format(I20,4E20.10)     
333 format(3E20.10)    
555 format('Ar',3(1x, e15.6)) 

     

      close(17)
      close(18)
      close(19)
      close(20)
      close(21)

   end program GaussMC



!------------------------------------------------------------------------!   
! Better random number generator                                         !
!------------------------------------------------------------------------!
      subroutine rmaset(iunit)
!     initializing routine for rmar, must be called before 
!     beginning to generate pseudorandom numbers by means of the
!     subroutine rmar
!     ranges : 0<=ij<=31328 and 0<=kl<=30081.
!     implicit real*8 (a-h,o-z)
      implicit none
      real*8 u,c,cd,cm,s,t
      integer i,j,k,m,iunit,n,kl,ij,ii,jj

      common/raset1/u(97),c,cd,cm,i,j
      data ij,kl/1802,9373/
!
      if(iunit.eq.0)                            then
       i=mod(ij/177, 177)+2
       j=mod(ij,     177)+2
       k=mod(kl/169, 178)+1
       m=mod(kl,     169)
       print '(a,2i7,4i4)',' marsaglia initialized (seed=0): ',ij,kl,i,j,k,m
!
       do ii=1,97
        s=0.0d00
        t=0.5d00
        do jj=1,24
         n=mod(mod(i*j,179)*k, 179)
         i=j
         j=k
         k=n
         m=mod(53*m+1, 169)
         if(mod(m*n,64).ge.32) s=s+t
         t=0.5d00*t
        enddo
        u(ii)=s
       enddo
!
       c =  362436.d00/16777216.d00
       cd= 7654321.d00/16777216.d00
       cm=16777213.d00/16777216.d00
!
       i=97
       j=33
!
      else
      
       print*,'marsaglia initialized, seed=',iunit
       open(41,file='rmaseed.dat',status='old',form='unformatted')
       rewind 41
       read(41) u,c,cd,cm,i,j
       close(41)
      end if
!
      return
      end

!--------------------------------------------------------------------------------------------!
      subroutine rmar(xr)          
!     pseudo random number generator 
!     proposed by marsaglia, zaman and tsang
!     implicit real*8 (a-h,o-z)
      implicit none
      real*8 xr,u,c,cd,cm,uni
      integer i,j
      common/raset1/u(97),c,cd,cm,i,j
!
      uni=u(i)-u(j)
      if(uni.lt. 0.0e00) uni=uni+1.0e00
      u(i)=uni           
      i=i-1              
      if(i.eq.0) i=97    
      j=j-1              
      if(j.eq.0) j=97   
      c=c-cd            
      if(c  .lt. 0.0e00)   c=c+cm
      uni=uni-c
      if(uni.lt. 0) uni=uni+1
      xr=uni
!
      return
      end