******************************************************************************
C                            PIMCOBSV_GRSOL.F
C Calculate the energy and superfluid density (using winding numbers).
******************************************************************************
C Subroutines and functions:
      include 'metrop_sqr3sc.f'
      include 'getdatadt_sqr3sc.f'       !!!!!!Check!!!!!!!!!!!
      include 'commen_sqr3sc.f'
      include 'permid_sqr3sc.f'
      include 'mgauss_sqr3sc.f'       
      include 'samdist_sqr3sc.f'       !!!!!!!!!Check!!!!!!!!!!!
      include 'ssfpeak_sqr3sc.f'
      include 'engtot_sqr3sc.f'       !!!!!!!!!Check!!!!!!!!!!!
      include 'distrib_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),du00(nr1)
      dimension dutfix(lev,nr1),ddutfix(lev,nr1)
      dimension rgrlat(nr1),rgroff(nr1),rgrsc(nr1),rgroffsc(nr1)
      dimension ugrlat(nr1),dvgrsc(nr1),dvgroffsc(nr1)
      dimension ugroff(nr1)
      dimension ssf(0:nssfx,0:nssfy),ssf_sum(0:nssfx,0:nssfy)
      dimension dhistogram(nhistx,nhisty),zdist(nhistz)
      integer image(2*m,nat)

      call getdata(r1,rgrlat,rgroff,rgrsc,rgroffsc,
     x     dutfix,ddutfix,u00,du00,
     x     ugrlat,ugroff,dvgrsc,dvgroffsc,nd,ngrlat,ngroff,ngrsc,
     x     ngroffsc)

      open(89,file='data/pthermsqr3t200m140d358n36.dat',
     x     status='old')
      open(90,file='data/pthermsqr3t200m140d358n36.dat',
     x     status='unknown')
c      open(100,file='junk636m100.dat',status='unknown')

      itherm=1            ! 0 to start anew, 1 to start from old data

*//////////////////////////////////////////////////////////////////////////////
C Initialize particle positions.  For ITHERM=0, use subroutine SOLID.F.
C For ITHERM=1, read positions from old data data file.
*//////////////////////////////////////////////////////////////////////////////
      
      if(itherm.eq.0) then
         call initialize(x,y,z,xp,yp,zp)
         do i=1,2*m
            do j=1,nat
               image(i,j)=j
            end do
         end do
         idum=-abs(int(x(1,1)*1.e5))
      else
         do i=1,2*m
            do j=1,nat
               read(89,*) 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))
      end if
      
*//////////////////////////////////////////////////////////////////////////////
C Calculate observables.  Note the superfluid density must be found with
C respect to all the helium atoms, both frozen and dynamic.
*//////////////////////////////////////////////////////////////////////////////
      
      icount2=0
      icount3=0
      icount4=0
      nacc=0
      sum=0.
      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,
     x           rgroff,dutfix,ddutfix,u00,ugrlat,
     x           ugroff,image,idum,nd,ngrlat,ngroff,icount2,
     x           icount3,icount4,nacc)
         end do
         call energy(x,y,z,r1,rgrlat,rgroff,rgrsc,rgroffsc,
     x        du00,ugrlat,
     x        ugroff,dvgrsc,dvgroffsc,nd,ngrlat,ngroff,ngrsc,
     x        ngroffsc,eng_tot,eng_pot)
c         call ssfobsv(x,y,z,ssfpeak1,ssfpeak2)
         print*,eng_tot,eng_pot,ssfpeak1,ssfpeak2
      end do

      do i=1,2*m
         do j=1,nat
            write(90,*) x(i,j),y(i,j),z(i,j),image(i,j)
         end do
      end do
      end


