*****************************************************************************
C                            PTHERM_SQR3.F
C Thermalize positions for the first layer atoms using semiclassical
C interactions between the helium atoms and the helium atoms and the 
C substrate.  Substrate corrugations are included.
C Records thermalized postitions in external file.  
******************************************************************************
C Subroutines and functions:
      include 'metrop_sqr3sc.f'
      include 'getdata_sqr3sc.f'    
      include 'commen_sqr3sc.f'
      include 'permid_sqr3sc.f'     
      include 'mgauss_sqr3sc.f' 
      include 'samdist_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),rgroff(nr1)
      dimension ugrlat(nr1),ugroff(nr1)
      integer image(2*m,nat)

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

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

      itherm=0            ! 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(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))
      end if
      
*//////////////////////////////////////////////////////////////////////////////
C Thermalize the configuration. 
*//////////////////////////////////////////////////////////////////////////////
      
      icount2=0
      icount3=0
      icount4=0
      nacc=0
      sum=0.
      do itherm=1,ntherm
         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
      end do
      fac=float(ntherm)*float(m)*float(nmov)
      print*, float(nacc)/fac
      
      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
      end
