************************************************************************
C                         SAMDIST_LAYR1.F
C The interactions are found by using a 3 point interpolation formula
C on the tabulated data.  See Allen and Tildesley, page 144.
************************************************************************
      subroutine samdist(x,y,z,r1,rgrlat,rgroff,u00,ugrlat,
     x     ugroff,iatomg,nsl1,nsl2,ioffd,islice,nd,ngrlat,
     x     ngroff,ipc,natp,move,result)

      include 'pimc_sqr3sc.dcl'
      dimension r1(nr1),r2(nr1)
      dimension u00(nr1),dist(3),uint(3)
      dimension xx(3),yy(3),uoff(3)
      dimension x(2*m,nat),y(2*m,nat),z(2*m,nat)
      dimension rgrlat(nr1),rgroff(nr1)
      dimension ugrlat(nr1),ugroff(nr1)
      dimension ydum1(nait),ydum2(nait)

      integer im(npol),in(npol),iatomg(nperm,4)

      rm=2.9673                       ! Aziz potential parameters 
      c6=1.3732412                         
      c8=0.4253785
      c10=0.178100
      c1=c6*rm**6
      c2=c8*rm**8
      c3=c10*rm**10

      pi=4.*atan(1.)                ! Other useful stuff.  Note we must
      nexp=lev-ioffd                 ! calculate values on the two end
      nslm=(nsl1+nsl2)/2              ! slices and their bisecting slice.
      nsld=(nsl1-nsl2)/2
      dtau=float(nsld)*tau

      vec1=0.0                       ! Zero these accumulators
      vec2=0.0
      usum1=0.0
      usum2=0.0
      result=0.0
      uchk1=0.0
      uchk2=0.0

      rcut_gr=rgrlat(ngrlat-5)   ! Cutoff value for He-graphite interaction
      ucut_gr=ugrlat(ngrlat-5)
      delr1=r1(2)-r1(1)
      delrgr=rgrlat(2)-rgrlat(1)
      delrgr2=rgroff(2)-rgroff(1)

      fac1=4.*alambda*dtau
      iarg1=islice+nsl1              ! Upper slice
      iarg2=islice+nsl2              ! Lower slice
      iargm=islice+nslm              ! Middle (bisecting) slice

*---------------------------------------------------------------------------
C Calculate sums for spring and interaction terms.
*---------------------------------------------------------------------------

      do i=1,natp
                                ! Interslice displacements
c         print*,ipc,iatomg(ipc,i)
         iag=iatomg(ipc,i) 
         tx1=abs(x(iargm,iag)-x(iarg2,iag))  
         if(tx1.gt.xdim/2.) tx1=xdim-tx1
         ty1=abs(y(iargm,iag)-y(iarg2,iag))
         if(ty1.gt.ydim/2.) ty1=ydim-ty1
         tz1=abs(z(iargm,iag)-z(iarg2,iag))
         
	 tx2=abs(x(iarg1,iag)-x(iargm,iag))
         if(tx2.gt.xdim/2.) tx2=xdim-tx2
 	 ty2=abs(y(iarg1,iag)-y(iargm,iag))
         if(ty2.gt.ydim/2.) ty2=ydim-ty2
	 tz2=abs(z(iarg1,iag)-z(iargm,iag))

	 vec1=vec1+tx1*tx1+ty1*ty1+tz1*tz1
	 vec2=vec2+tx2*tx2+ty2*ty2+tz2*tz2

*/////////////////////////////////////////////////////////////////////////
C Graphite Contribution.                                                 /  
*-------------------------------------------------------------------------
C Because of the distance between the substrate and the active layers,   /
C There's no long range correction for the graphite interaction.         /
C Just calculate the interaction out to whatever distance is needed      /
C beforehand--at least 10 Angstroms.
*/////////////////////////////////////////////////////////////////////////

         xx(1)=x(iarg2,iag)
         xx(2)=x(iargm,iag)
         xx(3)=x(iarg1,iag)
         yy(1)=y(iarg2,iag)
         yy(2)=y(iargm,iag)
         yy(3)=y(iarg1,iag)
         dist(1)=z(iarg2,iag)
         dist(2)=z(iargm,iag)
         dist(3)=z(iarg1,iag)

         do ix=1,3
            uint(ix)=0.0               ! First do the lateral part.
            utmp=0.0
c            if(dist(ix).lt.rgrlat(1)) then
            if(dist(ix).lt.1.8) then
               move=2
               goto 998
            end if
            if(dist(ix).lt.rcut_gr) then
               rlow=(dist(ix)-rgrlat(1))/delrgr
               ilow=int(rlow)+1
               xi=rlow-ilow+1.

               uk=ugrlat(ilow)
               uk1=ugrlat(ilow+1)
               uk2=ugrlat(ilow+2)
               t1=uk+(uk1-uk)*xi
               t2=uk1+(uk2-uk1)*(xi-1.)
               utmp=t1+(t2-t1)*xi*.5

               uint(ix)=utmp*float(nsld)  
            end if
                                           ! Now add the corrugations
            corrug=cos(2.*pi/3.+4.*pi*xx(ix)/sqr3/a_gr)
     x           + cos(4.*pi/3.+4.*pi*xx(ix)/sqr3/a_gr)
     x           + cos(2.*pi/3.-2.*pi*xx(ix)/sqr3/a_gr
     x                                  +2.*pi*yy(ix)/a_gr)
     x           + cos(4.*pi/3.-2.*pi*xx(ix)/sqr3/a_gr
     x                                  +2.*pi*yy(ix)/a_gr)
     x           + cos(4.*pi/3.+2.*pi*xx(ix)/sqr3/a_gr
     x                                  +2.*pi*yy(ix)/a_gr)
     x           + cos(8.*pi/3.+2.*pi*xx(ix)/sqr3/a_gr
     x                                  +2.*pi*yy(ix)/a_gr)
            uoff(ix)=0.0
            utmp=0.0
            if(dist(ix).lt.1.8) then
               move=2
               goto 998
            end if

            if(dist(ix).lt.rgroff(ngroff-3)) then
               rlow=(dist(ix)-rgroff(1))/delrgr2
               ilow=int(rlow)+1
               xi=rlow-ilow+1.
               
               uk=ugroff(ilow)
               uk1=ugroff(ilow+1)
               uk2=ugroff(ilow+2)
               t1=uk+(uk1-uk)*xi
               t2=uk1+(uk2-uk1)*(xi-1.)
               utmp=t1+(t2-t1)*xi*.5
               
               uoff(ix)=2.*utmp*corrug*float(nsld)  
            end if
         end do

         usum1=usum1+(uint(1)+uint(2))/2.+(uoff(1)+uoff(2))/2.
         usum2=usum2+(uint(2)+uint(3))/2.+(uoff(2)+uoff(3))/2.

*/////////////////////////////////////////////////////////////////////////
C Helium film contributions
*/////////////////////////////////////////////////////////////////////////

         do j=1,nat
                                    ! Particle distances on indiviual slices
            if(j.eq.iag) goto 10
            xl=abs(x(iarg2,iag)-x(iarg2,j))
            if(xl.gt.xdim/2.) xl=xdim-xl
            xm=abs(x(iargm,iag)-x(iargm,j))
            if(xm.gt.xdim/2.) xm=xdim-xm
            xu=abs(x(iarg1,iag)-x(iarg1,j))
            if(xu.gt.xdim/2.) xu=xdim-xu

            yl=abs(y(iarg2,iag)-y(iarg2,j))
            if(yl.gt.ydim/2.) yl=ydim-yl
            ym=abs(y(iargm,iag)-y(iargm,j))
            if(ym.gt.ydim/2.) ym=ydim-ym
            yu=abs(y(iarg1,iag)-y(iarg1,j))
            if(yu.gt.ydim/2.) yu=ydim-yu

            zl=abs(z(iarg2,iag)-z(iarg2,j))
            zm=abs(z(iargm,iag)-z(iargm,j))
            zu=abs(z(iarg1,iag)-z(iarg1,j))

            dist(1)=sqrt(xl*xl+yl*yl+zl*zl)
            dist(2)=sqrt(xm*xm+ym*ym+zm*zm)
            dist(3)=sqrt(xu*xu+yu*yu+zu*zu)

            do id=1,3
               if(dist(id).lt.rmin) then
                  move=2
                  goto 998
               end if
            end do

*---------------------------Diagonal Elements-------------------------------
            do ix=1,3
               uint(ix)=0.0
               if(dist(ix).lt.rcut) then
                  rlow=(dist(ix)-r1(1))/delr1
                  ilow=int(rlow)+1
                  xi=rlow-ilow+1
                  uk=u00(ilow)
                  uk1=u00(ilow+1)
                  uk2=u00(ilow+2)
                  t1=uk+(uk1-uk)*xi
                  t2=uk1+(uk2-uk1)*(xi-1.)
                  utmp=t1+(t2-t1)*xi*.5
                  uint(ix)=utmp*float(nsld)  
               end if
            end do
*----------------------------------------------------------------------------
            usum1=usum1+(uint(1)+uint(2))/2.
            usum2=usum2+(uint(2)+uint(3))/2.
 10         continue
         end do
      end do

      explow=vec1/fac1+usum1   !+ucor_low
      expup=vec2/fac1+usum2    !+ucor_up
      result=explow+expup
c      print*,vec1/fac1,usum1
c      print*,explow,expup,result
 998  continue
      return
      end
