******************************************************************************
C                       PIMCDISTRIB_GR.F
C Calculate potential, kinetic, and total energies.
******************************************************************************
C Subroutines and functions:
      include 'metrop_sqr3sc.f'
      include 'getdata_sqr3sc.f'
      include 'initialize_sqr3sc.f'
      include 'permid_sqr3sc.f'         
c      include 'perm1b_sqr3sc.f'
      include 'mgauss_sqr3sc.f'
c      include 'samdist_sqr3sc.f'
      include 'samdist_sqr3sc75.f'
      include 'distrib_sqr3sc.f'
      include 'zdist_sqr3sc.f'
      include 'dstpbc_sqr3sc.f'
      include 'ratint.f'
      include 'gasdev.f'
      include 'ran1.f'
      include 'ran2.f'
      include 'pimc_sqr3sc.dcl'
******************************************************************************

      dimension x(2*m,nat),y(2*m,nat),z(2*m,nat)
      dimension xp(2*m,nat),yp(2*m,nat),zp(2*m,nat)
      dimension r1(nr1),u00(nr1)
      dimension dutfix(lev,nr1),ddutfix(lev,nr1)
      dimension rgrlat(nr1),ugrlat(nr1)
      dimension rgroff(nr1),ugroff(nr1)
      dimension dhistogram(nhistx,nhisty),zdist(nhistz)
      
      integer image(2*m,nat)

      call getdata(r1,rgrlat,rgroff,
     x     dutfix,ddutfix,u00,ugrlat,ugroff,nd,ngrlat,ngroff)

      open(99,file='data/pthermsqr3t200m150dn36.dat',
     x     status='old')
c      open(100,file='data/pthermsqr3t200m150dn36.dat',
c     x     status='unknown')
      open(100,file='junk.dat',status='unknown')

      open(101,file='data/distribt200m100d636n36c75.dat',
     x     status='unknown')
c      open(102,file='data/zdistt200m20d6x3n1.dat',
c     x     status='unknown')

*//////////////////////////////////////////////////////////////////////////
C Read positions.
*//////////////////////////////////////////////////////////////////////////

      do i=1,2*m
         do j=1,nat
            read(99,*) x(i,j),y(i,j),z(i,j),image(i,j)
            xp(i,j)=x(i,j)
            yp(i,j)=y(i,j)
            zp(i,j)=z(i,j)
         end do
      end do
      idum=-abs(int(x(1,1)*1.e5))

*//////////////////////////////////////////////////////////////////////////
C Calculate observables.
*//////////////////////////////////////////////////////////////////////////
      
      nacc=0
      sum_final=0.0
      sum2_final=0.0
      xb=float(nbin)
      xn=float(niter)
      fac=1./(4.*float(m)*tau*float(nat)*alambda)
      do iter=1,niter
         do ibin=1,nbin
            call metrop(x,y,z,xp,yp,zp,r1,rgrlat,rgroff,
     x           dutfix,ddutfix,u00,ugrlat,ugroff,
     x           image,idum,nd,ngrlat,ngroff,icount2,icount3,
     x           icount4,nacc)
         end do
         call distrib(x,y,dhistogram)
c         call zdistobsv(z,zdist)
      end do
      
      fac=1./float(m)/float(niter)
      delx=xdim/float(nhistx)
      dely=ydim/float(nhisty)

      do i=1,2*m
         do j=1,nat
            write(100,*) x(i,j),y(i,j),z(i,j),image(i,j)
         end do
      end do
      
      delx=xdim/float(nhistx)
      dely=ydim/float(nhisty)
       do i=1,nhistx
         xhist=float(i-1)*delx+delx/2.
         do j=1,nhisty
            yhist=float(j-1)*dely+dely/2.
            write(101,*) xhist,yhist,dhistogram(i,j)*fac
        end do
      end do

c      delz=zdim/float(nhistz)
c       do i=1,nhistz
c         zhist=float(i-1)*delz+delz/2.
c         write(102,*) zhist,zdist(i)*fac
c      end do
      end

