*******************************************************************************
C                             ENGTOT_LAYR1.F
C
C Calculate the kinetic and total energies.
C Uses ratint,locate.
*******************************************************************************
      subroutine energy(x,y,z,r1,rgrlat,rgroff,rgrsc,rgroffsc,
     x     du00,ugrlat,
     x     ugroff,dvgrsc,dvgroffsc,nd,ngrlat,ngroff,ngrsc,ngroffsc,
     x     eng_tot,eng_pot)

      include 'pimc_sqr3sc.dcl'
      dimension r1(nr1),du00(nr1)
      dimension x(2*m,nat),y(2*m,nat),z(2*m,nat)
      dimension ydum1(nait),ydum2(nait),ydum3(nait)
      dimension rgrlat(nr1),ugrlat(nr1)
      dimension rgroff(nr1),ugroff(nr1)
      dimension rgrsc(nr1),dvgrsc(nr1)
      dimension rgroffsc(nr1),dvgroffsc(nr1)

      pi=4.*atan(1.)
      eps=10.94                      ! Aziz potential parameters.
      rm=2.970
      a=1.9221529e5
      beta=-1.89296514
      alpha=10.73520708
      c6=1.34920045
      c8=0.41365922
      c10=0.17078164
      d=1.4135

      rcut_gr=rgrlat(ngrlat)

      delr1=r1(2)-r1(1)
      delrgr=rgrlat(2)-rgrlat(1)
      delrgr2=rgroff(2)-rgroff(1)
      delrgrsc=rgrsc(2)-rgrsc(1)
      delrgroffsc=rgroffsc(2)-rgroffsc(1)

*///////////////////////////////////////////////////////////////////////////
C Long-range correction to the he-he total and kinetic energies.  
C The active layer densities are about .085.  
*///////////////////////////////////////////////////////////////////////////

      den_comp1=.120             ! Density of each active layer

C Check to see if first active layer is complete.  If not, use
C appropriate density.
      
      den_act1=nat/xdim/ydim
      den_act2=0.
      if(den_act1.gt.den_comp1) then
         den_tot=den_act1
         den_act1=den_comp1
         den_act2=den_tot-den_act1             ! No atoms at this level
      end if

      rlow=(rcut-r1(1))/delr1
      ilow=int(rlow)+1
      indx=min(max(ilow-(nait-1)/2,1),nd+1-nait)
      call ratint(r1(indx),du00(indx),nait,rcut,ucut,err)
      scale=ucut*(rcut/rm)**6

      fac_act1=pi*den_act1*rm**6*c6*scale/4./rcut**4
      fac_act2=pi*den_act2*rm**6*c6*scale/4./rcut**4
      
      ucor=fac_act1+fac_act2

      vcor=-pi*(den_act1+den_act2)
     x     *eps*rm**6*c6/4./rcut**4

*///////////////////////////////////////////////////////////////////////////
C Sum particle-particle and particle substrate energies.
*///////////////////////////////////////////////////////////////////////////

      dudtsumf=0.0
      sumho=0.0                    ! 'Harmonic Osc.' contribution to KE

      do islice=1,m
         vsum_gr=0.0
         vsum_act=0.0
         dudtsum_gr=0.0
         dudtsum_act=0.0

         isp1=islice+1
         do i=1,nat
*-------------------------Spring contribution to k.e.-----------------------

            call dstpbc(x(islice,i),x(isp1,i),xdim,x1,sgn)
            call dstpbc(y(islice,i),y(isp1,i),ydim,y1,sgn)
            z1=z(islice,i)-z(isp1,i)
            sumho=sumho+x1*x1+y1*y1+z1*z1

*-----------------------------Graphite Contribution-------------------------
            
            dutmp1=0.0
            dutmp2=0.0
            duoff1=0.0
            duoff2=0.0
            dvtmp1=0.0
            dvtmp2=0.0
            distgr1=z(islice,i)
            distgr2=z(islice+1,i)

            if(distgr1.le.rcut_gr) then 
               ! Get value for laterally averaged potential V_0.
               rlow=(distgr1-rgrlat(1))/delrgr
               ilow=int(rlow)+1
               indx=min(max(ilow-(nait-1)/2,1),ngrlat+1-nait)
               call ratint(rgrlat(indx),ugrlat(indx),nait,distgr1,
     x              dutmp1,dy)
               dutmp1=dutmp1/tau

               ! Get value for z-derivative of V_0.
               rlsc=(distgr1-rgrsc(1))/delrgrsc
               ilsc=int(rlsc)+1
               xi=rlsc-ilsc+1.
               uk=dvgrsc(ilsc)
               uk1=dvgrsc(ilsc+1)
               uk2=dvgrsc(ilsc+2)
               t1=uk+(uk1-uk)*xi
               t2=uk1+(uk2-uk1)*(xi-1.)
               dvtmp1=t1+(t2-t1)*xi*.5
            end if


            x1=x(islice,i)
            y1=y(islice,i)
            ! Calculate the corrugation.
            corrug1=cos(2.*pi/3.+4.*pi*x1/sqrt(3.)/a_gr)
     x           + cos(4.*pi/3.+4.*pi*x1/sqrt(3.)/a_gr)
     x           + cos(2.*pi/3.-2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           + cos(4.*pi/3.-2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           + cos(4.*pi/3.+2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           + cos(8.*pi/3.+2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)

            ! Calculate the x-derivative of the corrugation.
            corrug1dx=-4.*pi/sqrt(3.)/a_gr
     x                *sin(2.*pi/3.+4.*pi*x1/sqrt(3.)/a_gr)
     x           -4.*pi/sqrt(3.)/a_gr
     x                *sin(4.*pi/3.+4.*pi*x1/sqrt(3.)/a_gr)
     x           +2.*pi/sqrt(3.)/a_gr
     x                *sin(2.*pi/3.-2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           +2.*pi/sqrt(3.)/a_gr
     x                *sin(4.*pi/3.-2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           -2.*pi/sqrt(3.)/a_gr
     x                *sin(4.*pi/3.+2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           -2.*pi/sqrt(3.)/a_gr
     x                *sin(8.*pi/3.+2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)

            ! Calculate the y-derivative of the corrugation.
            corrug1dy=-2.*pi/a_gr*sin(2.*pi/3.-2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           -2.*pi/a_gr*sin(4.*pi/3.-2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           -2.*pi/a_gr*sin(4.*pi/3.+2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)
     x           -2.*pi/a_gr*sin(8.*pi/3.+2.*pi*x1/sqrt(3.)/a_gr
     x                                  +2.*pi*y1/a_gr)

            ! Calculate the z-dependent expansion coefficient.
            utmp1=0.0
            rlow=(distgr1-rgroff(1))/delrgr2
            ilow=int(rlow)+1
            indx=min(max(ilow-(nait-1)/2,1),ngroff+1-nait)
            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.)
            utmp1=t1+(t2-t1)*xi*.5
            utmp1=utmp1/tau
            duoff1=2.*utmp1*corrug1   ! This is the contribution of the
                                     ! corrugation terms to the energy.

            ! Calculate the z-derivative of the expansion coefficient.
            dvofftmp1=0.0
            rlow=(distgr1-rgroffsc(1))/delrgroffsc
            ilow=int(rlow)+1
            xi=rlow-ilow+1.
            uk=dvgroffsc(ilow)
            uk1=dvgroffsc(ilow+1)
            uk2=dvgroffsc(ilow+2)
            t1=uk+(uk1-uk)*xi
            t2=uk1+(uk2-uk1)*(xi-1.)
            dvofftmp1=t1+(t2-t1)*xi*.5  ! This is the derivative.

            sc_cor11=(dvtmp1+2.*dvofftmp1*corrug1)
     x           *(dvtmp1+2.*dvofftmp1*corrug1)
     x           *alambda*tau*tau/4.
            sc_cor12=4.*utmp1*utmp1*corrug1dx*corrug1dx
     x           *alambda*tau*tau/4.
            sc_cor13=4.*utmp1*utmp1*corrug1dy*corrug1dy
     x           *alambda*tau*tau/4.

            if(distgr2.le.rcut_gr) then 
               rlow=(distgr2-rgrlat(1))/delrgr
               ilow=int(rlow)+1
               indx=min(max(ilow-(nait-1)/2,1),ngrlat+1-nait)
               call ratint(rgrlat(indx),ugrlat(indx),nait,distgr2,
     x              dutmp2,dy)
               dutmp2=dutmp2/tau

               rlsc=(distgr2-rgrsc(1))/delrgrsc
               ilsc=int(rlsc)+1
               xi=rlsc-ilsc+1.
               uk=dvgrsc(ilsc)
               uk1=dvgrsc(ilsc+1)
               uk2=dvgrsc(ilsc+2)
               t1=uk+(uk1-uk)*xi
               t2=uk1+(uk2-uk1)*(xi-1.)
               dvtmp2=t1+(t2-t1)*xi*.5
            end if
                   
            x2=x(islice+1,i)
            y2=y(islice+1,i)
            corrug2=cos(2.*pi/3.+4.*pi*x2/sqrt(3.)/a_gr)
     x           + cos(4.*pi/3.+4.*pi*x2/sqrt(3.)/a_gr)
     x           + cos(2.*pi/3.-2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           + cos(4.*pi/3.-2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           + cos(4.*pi/3.+2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           + cos(8.*pi/3.+2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
            
             corrug2dx=-4.*pi/sqrt(3.)/a_gr
     x                *sin(2.*pi/3.+4.*pi*x2/sqrt(3.)/a_gr)
     x           -4.*pi/sqrt(3.)/a_gr
     x                *sin(4.*pi/3.+4.*pi*x2/sqrt(3.)/a_gr)
     x           +2.*pi/sqrt(3.)/a_gr
     x                *sin(2.*pi/3.-2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           +2.*pi/sqrt(3.)/a_gr
     x                *sin(4.*pi/3.-2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           -2.*pi/sqrt(3.)/a_gr
     x                *sin(4.*pi/3.+2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           -2.*pi/sqrt(3.)/a_gr
     x                *sin(8.*pi/3.+2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)

             corrug2dy=-2.*pi/a_gr*sin(2.*pi/3.-2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           -2.*pi/a_gr*sin(4.*pi/3.-2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           -2.*pi/a_gr*sin(4.*pi/3.+2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
     x           -2.*pi/a_gr*sin(8.*pi/3.+2.*pi*x2/sqrt(3.)/a_gr
     x                                  +2.*pi*y2/a_gr)
            
            utmp2=0.0
            rlow=(distgr2-rgroff(1))/delrgr2
            ilow=int(rlow)+1
            indx=min(max(ilow-(nait-1)/2,1),ngroff+1-nait)
            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.)
            utmp2=t1+(t2-t1)*xi*.5
            utmp2=utmp2/tau
            duoff2=2.*utmp2*corrug2  

            dvofftmp2=0.0
            rlow=(distgr2-rgroffsc(1))/delrgroffsc
            ilow=int(rlow)+1
            xi=rlow-ilow+1.
            uk=dvgroffsc(ilow)
            uk1=dvgroffsc(ilow+1)
            uk2=dvgroffsc(ilow+2)
            t1=uk+(uk1-uk)*xi
            t2=uk1+(uk2-uk1)*(xi-1.)
            dvofftmp2=t1+(t2-t1)*xi*.5

            sc_cor21=(dvtmp2+2.*dvofftmp2*corrug2)
     x           *(dvtmp2+2.*dvofftmp2*corrug2)
     x           *alambda*tau*tau/4.
            sc_cor22=4.*utmp2*utmp2*corrug2dx*corrug2dx
     x           *alambda*tau*tau/4.
            sc_cor23=4.*utmp2*utmp2*corrug2dy*corrug2dy
     x           *alambda*tau*tau/4.

            sc_cor1=sc_cor11+sc_cor12+sc_cor13
            sc_cor2=sc_cor21+sc_cor22+sc_cor23

            dudtsum_gr=dudtsum_gr+(dutmp1+dutmp2+duoff1+duoff2)/2.
     x           +(sc_cor1+sc_cor2)/2.

*--------------------------Active Helium Contribution-----------------------
            do j=i+1,nat
               if(i.eq.j) goto 10
               dutmp1=0.0
               dutmp2=0.0
               
               call dstpbc(x(islice,i),x(islice,j),xdim,dx2,sgn)
               call dstpbc(y(islice,i),y(islice,j),ydim,dy2,sgn)
               dz2=z(islice,i)-z(islice,j)
               dist1=sqrt(dx2*dx2+dy2*dy2+dz2*dz2)

               call dstpbc(x(islice+1,i),x(islice+1,j),xdim,dx3,sgn)
               call dstpbc(y(islice+1,i),y(islice+1,j),ydim,dy3,sgn)
               dz3=z(islice+1,i)-z(islice+1,j)
               dist2=sqrt(dx3*dx3+dy3*dy3+dz3*dz3)

               if(dist1.le.rcut) then
                  if(dist1.le.2.1) dist1=2.1
                  rlow=(dist1-r1(1))/delr1
                  ilow=int(rlow)+1
                  indx=min(max(ilow-(nait-1)/2,1),nd+1-nait)
                  call ratint(r1(indx),du00(indx),nait,dist1,
     x                 dutmp1,dy)
                  dutmp1=dutmp1
               end if
               if(dist2.le.rcut) then
                  if(dist2.le.2.1) dist2=2.1
                  rlow=(dist2-r1(1))/delr1
                  ilow=int(rlow)+1
                  indx=min(max(ilow-(nait-1)/2,1),nd+1-nait)
                  call ratint(r1(indx),du00(indx),nait,dist2,
     x                 dutmp2,dy)
                  dutmp2=dutmp2
               end if

               dudtsum_act=dudtsum_act+(dutmp1+dutmp2)/2.
 10            continue
            end do
         end do
         dudtsumf=dudtsumf+dudtsum_act+dudtsum_gr
      end do
      fac=1./float(m)/float(nat)
      sumho=sumho/4./alambda/tau/tau*fac
      
c      print*,1.5/tau,-sumho,dudtsumf*fac,ucor
      eng_tot=1.5/tau-sumho+dudtsumf*fac+ucor
      eng_pot=dudtsumf*fac+ucor

      return
      end
