SUBROUTINE LSTMAKE

      include 'global.inc'
	CUT = (ctof + 1.)*(ctof + 1.)

* *** *******************************************************
*     maintain the atom in the box

      do i = 1,npart
       if(xc(i).lt.0.0) xc(i) = xc(i) + xpore
       if(xc(i).gt.xpore) xc(i) = xc(i) - xpore     
       if(yc(i).lt.0.0) yc(i) = yc(i) + ypore
       if(yc(i).gt.ypore) yc(i) = yc(i) - ypore     
      enddo

* *** ******************************************************
*     calculate the distance between molecules L1 and L2
*     and compare with the coordination radius CTOF
*     for choosing neighbors of L1 molecule

      DO L1 = 1,npart  !L1
        IT1 = 0 
        IT2 = 0 
        DO L2 = 1,npart 	 !L2
           IF (L2.ne.L1) then

            ZCL = ZC(L2) - ZC(L1)
            YCL = YC(L2) - YC(L1)
            XCL = XC(L2) - XC(L1)

            ZP = ZCL
			 
              DO Ly = -1,1	  !Ly

              YP = YCL + float(Ly)*ypore   ! periodic along y
			 
                DO Lx = -1,1	     !Lx
 
                 XP = XCL + float(Lx)*xpore   ! periodic along x 
 
                 DS = XP*XP + YP*YP + ZP*ZP

                 IF(DS.le.CUT) then
                    IT1 = IT1 +1
                    LIST1(L1,IT1,4) = L2
                    LIST1(L1,IT1,3) = 0 
                    LIST1(L1,IT1,2) = Ly 
                    LIST1(L1,IT1,1) = Lx 


                 ENDIF
	          enddo	!Lx
              enddo !Ly

  	     endif
       enddo !L2

       NABL1(L1) = IT1

      enddo !L1

      END

* *** *********************************************************
* ***

      subroutine pair(DS,PSP,ia)

      include 'global.inc'

      real psp,ds,dx
      integer ia     ! index of current contact 
      integer num

      NUM = INT( DS * 10.)
      DX = (DS - R2(NUM,ia))
      PSP = CS(1,NUM,ia) + DX*(CS(2,NUM,ia) 
     &    + DX*(CS(3,NUM,ia) + DX*CS(4,NUM,ia)))

      end

! *** -----------------------------

      subroutine wall_spline(DS,PSP,ia)

      include 'global.inc'

      real psp,ds,dx
      integer ia     ! index of current contact 
      integer num

      NUM = INT( DS * 10.)
      DX = (DS - R2w(NUM,ia))
      if(NUM.gt.4000)then
        NUM=4000
        DX=0.0
      endif
      PSP = CSw(1,NUM,ia) + DX*(CSw(2,NUM,ia) 
     &    + DX*(CSw(3,NUM,ia) + DX*CSw(4,NUM,ia)))

      end

* *** ********************************************************
* *** 

      subroutine potgr(xin,yin,zin,EG)

* *** *************************************************
*
*     calculates the wall interaction of IL particle (at position x,y,z)
*
*     input NE1 must be 0 if average potential (no corrugation) !!!!!!!
*
* *** *************************************************

      include 'global.inc'

	real DS1, DS2, EG, EG1, EG2

	  call reduce (x,y,z,xin,yin,zin)	  ! transform into positions	
      									  ! equivalent of the grid
	DS1 = z*z
	call wall_spline(DS1, EG1,1)
	DS2 = (zpore-z)*(zpore-z)
	call wall_spline(DS2, EG2,1)

!	print *, EG1, EG2, z
	EG = EG1 + EG2

      end

!*************************************************************
!** 
      subroutine reduce (x,y,z,xin,yin,zin)

! ***   Here any reduction due to periodic substrate
!       must be applied

!          No reduction needed for simple slits!!!

	x = xin
	y = yin
	z = zin

	end

!*************************************************************

      subroutine ener1 (x,y,z,ene,IL)

! *** *************************************************

!     calculates the pair interaction of IL-th particle (at position x,y,z)
!
! *** *************************************************

      include 'global.inc'

!      do i = 1,Npart
!         UTT(i,IL) = 0.0   ! important for those that are NOT neighbors
!	enddo

      MAX = NABL1(IL)

	ene = 0.0
      if(max.gt.0) then

        DO LL= 1,MAX
	     L = FLOAT(LIST1(IL,LL,4))

           ZR = ZC(L) - Z 
           YR = YC(L) - Y 
           XR = XC(L) - X  

           ZP = ZR  
           YP = YR + FLOAT(LIST1(IL,LL,2))*ypore
           XP = XR + FLOAT(LIST1(IL,LL,1))*xpore

           DS = XP * XP + YP * YP + ZP * ZP

           IF (DS.le.CTOFS) then
	 
              call pair(DS,psp,1)
	        if (psp.gt.1.e6) psp = 1.e6

!              UTT(L,IL) = UTT(L,IL) + PSP
              ene = ene + psp

	   endif

        enddo

      endif

      end

**************************************************************

      subroutine ener0 (x,y,z,ene)

* *** *************************************************
*
*     calculates the pair interaction of new particle 
*     (at position x,y,z) - NO LISTMAKE used !!!
*
* *** *************************************************

      include 'global.inc'

	ene = 0.0

        DO L= 1,npart
!	  UTT(L,Npart+1) = 0.0

           ZR = ZC(L) - Z 
           YR = YC(L) - Y 
           XR = XC(L) - X  

           ZP = ZR 

           do LL = -1,1

             YP = YR + FLOAT(LL)*ypore

             do LK = -1,1

               XP = XR + FLOAT(LK)*xpore

               DS = XP * XP + YP * YP + ZP * ZP

                IF (DS.le.CTOFS) then
	            psp = 1.e6
	            if (DS.gt.2.25) then
	 
                    call pair(DS,psp,1)	   								 
 
	            endif
 	 
!                UTT(L,Npart+1) = UTT(L,Npart+1) + PSP       
                ene = ene + psp

	          endif
	      enddo
         enddo
       enddo


      end

!**************************************************************

      SUBROUTINE ENERGY_tot  

*        compute the entire enthalpy
*             H = U + EQ + PV
* 
*         HTOT - total energy of the layer
*
*        UGTOT - total energy of wall interaction
*
*        UTT(L,l1) - temporary pair interaction
*
*        EGT(L1) - temporary wall interaction of L1 particle
*
*        ENMOLT(L1) - temporary total interaction of L1 particle
*
* *** ********************************************************
	
	include 'global.inc'

      call lstmake

      UGTOT = 0.0 
      UTOT = 0.0

      do L1 = 1,npart

         call period_check(xc(l1),yc(l1),xcb,ycb,iper)
         CALL POTgr(xcb,ycb,zc(l1),EG)
!	              print *, 'In ENERGY_tot : int_gr finished: ', EG

         call ener1(xcb,ycb,zc(l1),ene,l1)
!	              print *, 'In ENERGY_tot : inter particle: ',ene

	   xc(l1) = xcb
	   yc(l1) = ycb
	  
         EGT(L1) = EG 
         ENMOLT(L1) = ene + EG 

         UGTOT = UGTOT + EG
         UTOT = UTOT + EG + ene/2.0 

      enddo

      VTOT = zport*yport*xport
      HTOT = UTOT

!      HTOT = UTOT + PRESS * VTOT   !!!!!!!!!!!!

      END

!**************************************************************

      subroutine update

      include 'global.inc'

	call lstmake

	call ENERGY_tot

       ENTH = HTOT
       VOL = VTOT
	 UGT = UGTOT

       DO J = 1,N 

         EGG(J)=EGT(J)
	   ENMOL(J) = ENMOLT(J)

!	   UU(J,J) = 0.0

!         DO I = J+1,N
!            UU(I,J) = UTT(I,J)
!	      UU(j,i) = UU(i,j)
!	   enddo

!        if(ib.ge.21) write(67,'(2f12.5)')UGT/npart, ENTH/npart
       enddo

	end