subroutine energyandforces(nat,alat, rxyz0, fxyz, etot, count)
!*****************************************************************************************
! This subroutine evaluates the Stillinger Weber Silicon potential with linear scaling
! COPYRIGHT
!    Copyright (C) 2009 AIST, UNIBAS
!    This file is distributed under the terms of the
!    GNU General Public License, see 
!    http://www.gnu.org/copyleft/gpl.txt .
!
! Implementation: Original version was written by Tetsuya Morishita, AIST Tsukuba (JP) 
!                 Improved by M. Amsler, S. Goedecker, Basel University (CH), 2009
!
! Note:
!
!     aa is the parameter A  given on page 5263 of PRB 31, 5262 (1985).
!     bb is the parameter B  given on page 5263 of PRB 31, 5262 (1985).
!     ra is the parameter a  given on page 5263 of PRB 31, 5262 (1985).
!     gam and ramda are the parameters gamma and lambda 
!     for the 3-body term, respectively (see Eq. (2.5) in the paper).
!    
! Input:
!     nat, integer: the number of atoms
!     alat,  real(8), dim(3)   : the three edges of the orthoromic simulation cell, periodic boundaries are applied
!                                and atoms outside the cell will be brought back into the box
!     rxyz, real(8), dim(3,nat) : the xyz cartesian components of the atomic positions in Angstroem    
!                                
! Output:
!     fxyz,  real(8), dim(3,nat): the xyz cartesian forces in on the corresponding atomic components n eV/A  
!     etot, real(8)            : total potential energy, 2-body and 3-body, in eV
!     count, real(8)           : increased by 1.d0 at each call of this subroutine,
!                                needs to be initialized to 0.d0 before calling this routine for the first time
!
! Other variables:
!     p:  the 2-body potential energy
!     p3: the 3-body potential energy 
!     fx(i) , fy(i) , fz(i)  are the 2-body forces on atom i.
!     fx3(i), fy3(i), fz3(i) are the 3-body forces on atom i.
!     fxyz(3,nat) contains both 2-body and 3-body forces  
!
!     All units follow the description in PRB 31, 5262 (1985).
!*****************************************************************************************

        implicit real(8) (a-h,o-z)
        real*8                :: rxyz0(3,nat)
	real*8                :: fxyz(3,nat),etot
	real*8                :: x(nat),y(nat),z(nat), fx(nat), fy(nat), fz(nat), fx3(nat), fy3(nat), fz3(nat)
	real*8                :: lx,ly,lz, gam, ramda, aa, ra, raa, bb, xij, yij, zij, vl,f,pe, p, p3
	real*8                :: xjj, yjj, zjj,rij,xik,yik,zik, xkk,ykk,zkk,force
	real*8                :: rik, xjk, yjk, zjk, rjk, cosjik, cosijk, cosikj, cosjik3, cosijk3, cosikj3
	real*8                :: rija, rika, rjka, refi, refj, refk
	real*8                :: hi, hj, hk, hixij0, hixij1, hixij, hixik0,  hixik, HKXIK1
	real*8                :: hjxij0, hjxij1, hjxij, hkxik0, hixik1, hkxik, hjxjk0, hjxjk1, hjxjk, hkxkj0, hkxkj1, hkxkj, hiyij0, hiyij1
	real*8                :: hiyij, hiyik0, hiyik1, hiyik, hjyij0, hjyij1, hjyij
	real*8                :: hkyik0, hkyik1, hkyik, hjyjk0, hjyjk1, hjyjk, hkykj0, hkykj1, hkykj
	real*8                :: HIZIJ0,HIZIJ1,HIZIJ,HIZIK0,HIZIK1,HIZIK,HJZIJ0,HJZIJ1,HJZIJ,HKZIK0,HKZIK1,HKZIK
	real*8                :: HJZJK0,HJZJK1,HJZJK,HKZKJ0,HKZKJ1,HKZKJ,PV3
	real*8                :: c, d, uc
	real*8                :: xbar, ybar, zbar, xmax, ymax, zmax, xmin, ymin, zmin,xx, yy,zz
        real*8                :: alat(3), count,sigma,eps,esigma,isigma,c4,crainv
        real*8                :: invrij, invrik, invrjk, invrija, invrika, invrjka
        real(8), ALLOCATABLE, DIMENSION(:,:)     :: rxyz
        real(8), ALLOCATABLE, DIMENSION(:,:)     :: rel
        real(8), ALLOCATABLE, DIMENSION(:,:)     :: txyz
        real(8), ALLOCATABLE, DIMENSION(:,:)     :: f2ij,f3ij,f3ik
        integer, ALLOCATABLE, DIMENSION(:,:)     :: lsta
        integer, ALLOCATABLE, DIMENSION(:)       :: lstb
        integer, ALLOCATABLE, DIMENSION(:)       :: lay
        integer, ALLOCATABLE, DIMENSION(:,:,:,:) :: icell
        integer                                  :: nat,nnbrx
        integer                                  :: i, ii, j, nij, k
        parameter(gam=1.2d0,ramda=21.0d0,aa=7.049556277d0)
        parameter(ra=1.8d0,bb=0.6022245584d0)
        parameter(uc=4.5d0)
        parameter(sigma=2.0951d0,eps=2.167239428587d0)

        count=count+1.d0
        cut=sigma*ra*2.d0
        isigma=1.d0/sigma
        esigma=eps*isigma

! linear scaling calculation of verlet list, only serial 
        ll1=int(alat(1)/cut)
        if (ll1.lt.1) stop 'alat(1) too small'
        ll2=int(alat(2)/cut)
        if (ll2.lt.1) stop 'alat(2) too small'
        ll3=int(alat(3)/cut)
        if (ll3.lt.1) stop 'alat(3) too small'

        npr=1
        ncx=29
1234    ncx=ncx*2
        allocate(icell(0:ncx,-1:ll1,-1:ll2,-1:ll3))
        do 984,l3=-1,ll3
        do 984,l2=-1,ll2
        do 984,l1=-1,ll1
984     icell(0,l1,l2,l3)=0
        rlc1i=ll1/alat(1)
        rlc2i=ll2/alat(2)
        rlc3i=ll3/alat(3)

        do 983,iat=1,nat
        rxyz0(1,iat)=modulo(modulo(rxyz0(1,iat),alat(1)),alat(1))
        rxyz0(2,iat)=modulo(modulo(rxyz0(2,iat),alat(2)),alat(2))
        rxyz0(3,iat)=modulo(modulo(rxyz0(3,iat),alat(3)),alat(3))
        l1=int(rxyz0(1,iat)*rlc1i)
        l2=int(rxyz0(2,iat)*rlc2i)
        l3=int(rxyz0(3,iat)*rlc3i)

        ii=icell(0,l1,l2,l3)
        ii=ii+1
        icell(0,l1,l2,l3)=ii
        if (ii.gt.ncx) then
        write(10,*) count,'NCX too small',ncx
        deallocate(icell)
        goto 1234
        endif
        icell(ii,l1,l2,l3)=iat
983     continue


! duplicate all atoms within boundary layer
        laymx=ncx*(2*ll1*ll2+2*ll1*ll3+2*ll2*ll3+4*ll1+4*ll2+4*ll3+8)
        nn=nat+laymx
        allocate(rxyz(3,nn),lay(nn))
        do  iat=1,nat
        lay(iat)=iat
        rxyz(1,iat)=rxyz0(1,iat)
        rxyz(2,iat)=rxyz0(2,iat)
        rxyz(3,iat)=rxyz0(3,iat)
        enddo
        il=nat
! xy plane
        do l2=0,ll2-1
        do l1=0,ll1-1

        in=icell(0,l1,l2,0)
        icell(0,l1,l2,ll3)=in
        do ii=1,in
        i=icell(ii,l1,l2,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,l1,l2,ll3)=il
        rxyz(1,il)=rxyz(1,i)
        rxyz(2,il)=rxyz(2,i)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,l1,l2,ll3-1)
        icell(0,l1,l2,-1)=in
        do ii=1,in
        i=icell(ii,l1,l2,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,l1,l2,-1)=il
        rxyz(1,il)=rxyz(1,i)
        rxyz(2,il)=rxyz(2,i)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo

        enddo
        enddo


! yz plane
        do l3=0,ll3-1
        do l2=0,ll2-1

        in=icell(0,0,l2,l3)
        icell(0,ll1,l2,l3)=in
        do ii=1,in
        i=icell(ii,0,l2,l3)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,l2,l3)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)
        rxyz(3,il)=rxyz(3,i)
        enddo

        in=icell(0,ll1-1,l2,l3)
        icell(0,-1,l2,l3)=in
        do ii=1,in
        i=icell(ii,ll1-1,l2,l3)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,l2,l3)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)
        rxyz(3,il)=rxyz(3,i)
        enddo

        enddo
        enddo


! xz plane
        do l3=0,ll3-1
        do l1=0,ll1-1

        in=icell(0,l1,0,l3)
        icell(0,l1,ll2,l3)=in
        do ii=1,in
        i=icell(ii,l1,0,l3)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,l1,ll2,l3)=il
        rxyz(1,il)=rxyz(1,i)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)
        enddo

        in=icell(0,l1,ll2-1,l3)
        icell(0,l1,-1,l3)=in
        do ii=1,in
        i=icell(ii,l1,ll2-1,l3)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,l1,-1,l3)=il
        rxyz(1,il)=rxyz(1,i)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)
        enddo

        enddo
        enddo


! x axis
        do l1=0,ll1-1

        in=icell(0,l1,0,0)
        icell(0,l1,ll2,ll3)=in
        do ii=1,in
        i=icell(ii,l1,0,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,l1,ll2,ll3)=il
        rxyz(1,il)=rxyz(1,i)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,l1,0,ll3-1)
        icell(0,l1,ll2,-1)=in
        do ii=1,in
        i=icell(ii,l1,0,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,l1,ll2,-1)=il
        rxyz(1,il)=rxyz(1,i)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo

        in=icell(0,l1,ll2-1,0)
        icell(0,l1,-1,ll3)=in
        do ii=1,in
        i=icell(ii,l1,ll2-1,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,l1,-1,ll3)=il
        rxyz(1,il)=rxyz(1,i)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,l1,ll2-1,ll3-1)
        icell(0,l1,-1,-1)=in
        do ii=1,in
        i=icell(ii,l1,ll2-1,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,l1,-1,-1)=il
        rxyz(1,il)=rxyz(1,i)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo

        enddo


! y axis
        do l2=0,ll2-1

        in=icell(0,0,l2,0)
        icell(0,ll1,l2,ll3)=in
        do ii=1,in
        i=icell(ii,0,l2,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,l2,ll3)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,0,l2,ll3-1)
        icell(0,ll1,l2,-1)=in
        do ii=1,in
        i=icell(ii,0,l2,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,l2,-1)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo

        in=icell(0,ll1-1,l2,0)
        icell(0,-1,l2,ll3)=in
        do ii=1,in
        i=icell(ii,ll1-1,l2,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,l2,ll3)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,ll1-1,l2,ll3-1)
        icell(0,-1,l2,-1)=in
        do ii=1,in
        i=icell(ii,ll1-1,l2,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,l2,-1)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo

        enddo


! z axis
        do l3=0,ll3-1

        in=icell(0,0,0,l3)
        icell(0,ll1,ll2,l3)=in
        do ii=1,in
        i=icell(ii,0,0,l3)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,ll2,l3)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)
        enddo

        in=icell(0,ll1-1,0,l3)
        icell(0,-1,ll2,l3)=in
        do ii=1,in
        i=icell(ii,ll1-1,0,l3)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,ll2,l3)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)
        enddo

        in=icell(0,0,ll2-1,l3)
        icell(0,ll1,-1,l3)=in
        do ii=1,in
        i=icell(ii,0,ll2-1,l3)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,-1,l3)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)
        enddo

        in=icell(0,ll1-1,ll2-1,l3)
        icell(0,-1,-1,l3)=in
        do ii=1,in
        i=icell(ii,ll1-1,ll2-1,l3)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,-1,l3)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)
        enddo

        enddo


! corners
        in=icell(0,0,0,0)
        icell(0,ll1,ll2,ll3)=in
        do ii=1,in
        i=icell(ii,0,0,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,ll2,ll3)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,ll1-1,0,0)
        icell(0,-1,ll2,ll3)=in
        do ii=1,in
        i=icell(ii,ll1-1,0,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,ll2,ll3)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,0,ll2-1,0)
        icell(0,ll1,-1,ll3)=in
        do ii=1,in
        i=icell(ii,0,ll2-1,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,-1,ll3)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,ll1-1,ll2-1,0)
        icell(0,-1,-1,ll3)=in
        do ii=1,in
        i=icell(ii,ll1-1,ll2-1,0)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,-1,ll3)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)+alat(3)
        enddo

        in=icell(0,0,0,ll3-1)
        icell(0,ll1,ll2,-1)=in
        do ii=1,in
        i=icell(ii,0,0,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,ll2,-1)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo
        
        in=icell(0,ll1-1,0,ll3-1)
        icell(0,-1,ll2,-1)=in
        do ii=1,in
        i=icell(ii,ll1-1,0,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,ll2,-1)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)+alat(2)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo

        in=icell(0,0,ll2-1,ll3-1)
        icell(0,ll1,-1,-1)=in
        do ii=1,in
        i=icell(ii,0,ll2-1,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,ll1,-1,-1)=il
        rxyz(1,il)=rxyz(1,i)+alat(1)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo

        in=icell(0,ll1-1,ll2-1,ll3-1)
        icell(0,-1,-1,-1)=in
        do ii=1,in
        i=icell(ii,ll1-1,ll2-1,ll3-1)
        il=il+1
        if (il.gt.nn) stop 'enlarge laymx'
        lay(il)=i
        icell(ii,-1,-1,-1)=il
        rxyz(1,il)=rxyz(1,i)-alat(1)
        rxyz(2,il)=rxyz(2,i)-alat(2)
        rxyz(3,il)=rxyz(3,i)-alat(3)
        enddo



        allocate(lsta(2,nat))
        nnbrx=300
2345    nnbrx=3*nnbrx/2
        allocate(lstb(nnbrx*nat),rel(5,nnbrx*nat))

        indlstx=0


        npr=1
        iam=0

        cut2=cut**2
! assign contiguous portions of the arrays lstb and rel to the threads (this version only contains one thread)
        myspace=(nat*nnbrx)/npr  
        if (iam.eq.0) myspaceout=myspace
! Verlet list, relative positions
        indlst=0
      do 6000,l3=0,ll3-1
      do 6000,l2=0,ll2-1
      do 6000,l1=0,ll1-1
      do 6600,ii=1,icell(0,l1,l2,l3)
        iat=icell(ii,l1,l2,l3)
        if ( ((iat-1)*npr)/nat .eq. iam) then
        lsta(1,iat)=iam*myspace+indlst+1
        call sublstiat_l(iat,nn,ncx,ll1,ll2,ll3,l1,l2,l3,myspace, & 
             rxyz,icell,lstb(iam*myspace+1),lay,rel(1,iam*myspace+1),cut2,indlst)
        lsta(2,iat)=iam*myspace+indlst
        ipb=lsta(1,iat)
        ndat=lsta(2,iat)-lsta(1,iat)+1
        endif

6600    continue
6000    continue
        indlstx=max(indlstx,indlst)

           if (indlstx.ge.myspaceout) then
               write(10,*) count,'NNBRX too  small', nnbrx
               deallocate(lstb,rel)
               goto 2345
           endif

        npr=1
        iam=0
        npjx=300 ; npjkx=6000
!end of creating pairlist part------------------------------------------------------------

!start energy and force calculation-------------------------------------------------------
!set all variables to zero
        p=0.0d0
        p3=0.0d0
	pv3=0.0d0
        fx(:)=0.0d0
        fy(:)=0.0d0
        fz(:)=0.0d0
        fx3(:)=0.0d0
        fy3(:)=0.0d0
        fz3(:)=0.0d0
!-----------------------------------------------------------------------------------------
!     triple loop for the 2 and 3-body forces
!     do 20 i
!     do 30 j
!     do 40 k
!     the pairlists lsta and lstb are used for the perodic boundary conditions
!-----------------------------------------------------------------------------------------

      do 20 i=1,nat       !DO 20
      call subfeniat_l(i,nat,nnbrx,rel,cut2,p,p3,fx,fy,fz,fx3,fy3,fz3,lstb,lsta,isigma,sigma)
   20 continue
!-----------------------------------------------------------------------------------------
!*****if necessary,********
      do i=1,nat
        fx(i)=fx(i)+fx3(i)
        fy(i)=fy(i)+fy3(i)
        fz(i)=fz(i)+fz3(i)
      end do
!-----------------------------------------------------------------------------------------

        do i=1,nat
           fxyz(1,i)=fx(i)*esigma
           fxyz(2,i)=fy(i)*esigma
           fxyz(3,i)=fz(i)*esigma
         enddo
         etot=(p+p3)*eps
         deallocate(rxyz,icell,lay,lsta,lstb,rel)
end
!End of the force calculation-------------------------------------------------------------

    REAL*8  FUNCTION f(c)
        real*8 :: c,fm,ra, bb, aa,c4,crainv
        parameter(aa=7.049556277d0)
        parameter(ra=1.8d0,bb=0.6022245584d0)
        parameter(sigma=2.0951d0,eps=2.167239428587d0)

        if((c-ra).lt.0.d0) then
        crainv=1.0d0/(c-ra)
        c4=c*c*c*c
        f=aa*bb*4.0d0/(c4*c)*dexp(crainv)+aa*(bb/(c4)-1.0d0)*dexp(crainv)*crainv*crainv
        else
        f=0.d0
        endif

        return
        end function f

     REAL*8  function pe(d)
	real*8, parameter :: aa=7.049556277d0
	real*8, parameter :: bb=0.6022245584d0
        parameter(ra=1.8d0)
        parameter(sigma=2.0951d0,eps=2.167239428587d0)
        real*8::d

        if((d-ra).lt.0.d0) then  
          pe=aa*(bb/(d*d*d*d)-1.0d0)*dexp(1.0d0/(d-ra))
        else
          pe=0.d0
        endif
        return
      end function pe


!------------------------------------------------------------------------------------------
        subroutine sublstiat_l(iat,nn,ncx,ll1,ll2,ll3,l1,l2,l3,myspace, &
                   rxyz,icell,lstb,lay,rel,cut2,indlst)
! finds the neighbours of atom iat (specified by lsta and lstb) and and
! the relative position rel of iat with respect to these neighbours
        implicit real(8) (a-h,o-z)
        dimension rxyz(3,nn),lay(nn),icell(0:ncx,-1:ll1,-1:ll2,-1:ll3), &
                  lstb(0:myspace-1),rel(5,0:myspace-1)

        do 6363,k3=l3-1,l3+1
        do 6363,k2=l2-1,l2+1
        do 6363,k1=l1-1,l1+1
        do 6363,jj=1,icell(0,k1,k2,k3)
          jat=icell(jj,k1,k2,k3)
          if (jat.eq.iat) goto 6363
          xrel= rxyz(1,iat)-rxyz(1,jat)
          yrel= rxyz(2,iat)-rxyz(2,jat)
          zrel= rxyz(3,iat)-rxyz(3,jat)
          rr2=xrel**2 + yrel**2 + zrel**2
          if ( rr2 .le. cut2 ) then
           indlst=min(indlst,myspace-1)
           lstb(indlst)=lay(jat)
!        write(6,*) 'iat,indlst,lay(jat)',iat,indlst,lay(jat)
           tt=sqrt(rr2)
           tti=1.d0/tt
           rel(1,indlst)=xrel*tti
           rel(2,indlst)=yrel*tti
           rel(3,indlst)=zrel*tti
           rel(4,indlst)=tt
           rel(5,indlst)=tti
           indlst= indlst+1
          endif
6363    continue

        return
        end


subroutine subfeniat_l(i,nat,nnbrx,rel,cut2,p,p3,fx,fy,fz,fx3,fy3,fz3,lstb,lsta,isigma,sigma)
implicit none
        integer, intent(IN)   :: i,nat,nnbrx,lsta(2,nat),lstb(nnbrx*nat)
        integer               :: Ipb,Ipe,k,m,l,j,nij
	real*8, intent(IN)    :: rel(5,nnbrx*nat),cut2,isigma,sigma
        real*8, intent(INOUT) :: fx(nat),fy(nat),fz(nat),fx3(nat),fy3(nat),fz3(nat),p3,p
        real*8                :: hi, hj, hk, hixij0, hixij1, hixij, hixik0,  hixik, HKXIK1
	real*8                :: hjxij0, hjxij1, hjxij, hkxik0, hixik1, hkxik, hjxjk0, hjxjk1, hjxjk, hkxkj0, hkxkj1, hkxkj, hiyij0, hiyij1
	real*8                :: hiyij, hiyik0, hiyik1, hiyik, hjyij0, hjyij1, hjyij
	real*8                :: hkyik0, hkyik1, hkyik, hjyjk0, hjyjk1, hjyjk, hkykj0, hkykj1, hkykj
	real*8                :: HIZIJ0,HIZIJ1,HIZIJ,HIZIK0,HIZIK1,HIZIK,HJZIJ0,HJZIJ1,HJZIJ,HKZIK0,HKZIK1,HKZIK
	real*8                :: HJZJK0,HJZJK1,HJZJK,HKZKJ0,HKZKJ1,HKZKJ,PV3
	real*8                :: rij,xik,yik,zik,xij, yij, zij,f,pe
	real*8                :: rik, xjk, yjk, zjk, rjk, cosjik, cosijk, cosikj, cosjik3, cosijk3, cosikj3
	real*8                :: rija, rika, rjka, refi, refj, refk
        real*8                :: invrik,invrjk,invrij,gam,ramda,aa,ra,bb,uc,invrika,invrjka,invrija,force,c4,crainv

        parameter(gam=1.2d0,ramda=21.0d0,aa=7.049556277d0)
        parameter(ra=1.8d0,bb=0.6022245584d0)
        parameter(uc=4.5d0)
        
         Ipb = lsta(1,i)
         Ipe = lsta(2,i)

      do 30 l=Ipb,Ipe,+1 !DO 30
            j=lstb(l)
         if (j.le.i) go to 30
         nij=0
         rij=rel(4,l)*isigma
         invrij=rel(5,l)*sigma
         xij=rel(1,l)*rij
         yij=rel(2,l)*rij
         zij=rel(3,l)*rij


        if(rij.ge.2.d0*ra) go to 30
        if(rij.ge.ra) go to 31
        
        crainv=1.0d0/(rij-ra)
        c4=rij*rij*rij*rij
        force=aa*bb*4.0d0/(c4*rij)*dexp(crainv)+aa*(bb/(c4)-1.0d0)*dexp(crainv)*crainv*crainv

        fx(i)=force*xij*invrij+fx(i)
        fy(i)=force*yij*invrij+fy(i)
        fz(i)=force*zij*invrij+fz(i)

        fx(j)=-force*xij*invrij+fx(j)
        fy(j)=-force*yij*invrij+fy(j)
        fz(j)=-force*zij*invrij+fz(j)

        p=p+aa*(bb/(rij*rij*rij*rij)-1.0d0)*dexp(1.0d0/(rij-ra))

        nij=1

   31   continue

         do 40 m=Ipb,Ipe,+1 !DO 40
         k=lstb(m)
         if (k.le.j) goto 40
         invrik=rel(5,m)*sigma
         rik=rel(4,m)*isigma
         xik=rel(1,m)*rik
         yik=rel(2,m)*rik
         zik=rel(3,m)*rik

        if((rik.ge.ra).and.(nij.eq.0)) go to 40

        xjk=xik-xij
        yjk=yik-yij
        zjk=zik-zij

        rjk=sqrt(xjk*xjk+yjk*yjk+zjk*zjk)
        invrjk=1.d0/rjk

        if((rjk.ge.ra).and.(nij.eq.0)) go to 40
        cosjik=(xij*xik+yij*yik+zij*zik)*(invrij*invrik)
        cosijk=(-xij*xjk-yij*yjk-zij*zjk)*(invrij*invrjk)
        cosikj=(xik*xjk+yik*yjk+zik*zjk)*(invrik*invrjk)
        cosjik3=cosjik+1.0d0/3.0d0
        cosijk3=cosijk+1.0d0/3.0d0
        cosikj3=cosikj+1.0d0/3.0d0

        rija=rij-ra
        rika=rik-ra
        rjka=rjk-ra

        invrija=1.d0/rija
        invrika=1.d0/rika
        invrjka=1.d0/rjka

      if(rija.ge.0.0d0) then
      refi=0.0d0
      refj=0.0d0
      refk=ramda*exp(gam*invrika+gam*invrjka)
      else if((rija.lt.0.0d0).and.(rika.lt.0.0d0)) then
       if(rjka.lt.0.0d0) then
       refi=ramda*exp(gam*invrija+gam*invrika)
       refj=ramda*exp(gam*invrija+gam*invrjka)
       refk=ramda*exp(gam*invrika+gam*invrjka)
       else
       refi=ramda*exp(gam*invrija+gam*invrika)
       refj=0.0d0
       refk=0.0d0
       end if
      else if((rija.lt.0.0d0).and.(rjka.lt.0.0d0)) then
      refi=0.0d0
      refj=ramda*exp(gam*invrija+gam*invrjka)
      refk=0.0d0
else
goto 40
      refi=0.0d0
      refj=0.0d0
      refk=0.0d0
      end if

      hi=refi*cosjik3*cosjik3
      hj=refj*cosijk3*cosijk3
      hk=refk*cosikj3*cosikj3
      p3=p3+hi+hj+hk

      hixij0=2.0d0*(xik*invrik-xij*cosjik*invrij)
      hixij1=gam*xij*cosjik3*(invrija*invrija)
      hixij=refi*cosjik3*(hixij0-hixij1)*invrij
      hixik0=2.0d0*(xij*invrij-xik*cosjik*invrik)
      hixik1=gam*xik*cosjik3*(invrika*invrika)
      hixik=refi*cosjik3*(hixik0-hixik1)*invrik
      hjxij0=2.0d0*(-xjk*invrjk-xij*cosijk*invrij)
      hjxij1=gam*xij*cosijk3*(invrija*invrija)
      hjxij=refj*cosijk3*(hjxij0-hjxij1)*invrij
      hkxik0=2.0d0*(xjk*invrjk-xik*cosikj*invrik)
      hkxik1=gam*xik*cosikj3*(invrika*invrika)
      hkxik=refk*cosikj3*(hkxik0-hkxik1)*invrik
      hjxjk0=2.0d0*(-xij*invrij-xjk*cosijk*invrjk)
      hjxjk1=gam*xjk*cosijk3*(invrjka*invrjka)
      hjxjk=refj*cosijk3*(hjxjk0-hjxjk1)*invrjk
      hkxkj0=2.0d0*(-xik*invrik+xjk*cosikj*invrjk)
      hkxkj1=gam*xjk*cosikj3*(invrjka*invrjka)
      hkxkj=refk*cosikj3*(hkxkj0+hkxkj1)*invrjk


      hiyij0=2.0d0*(yik*invrik-yij*cosjik*invrij)
      hiyij1=gam*yij*cosjik3*(invrija*invrija)
      hiyij=refi*cosjik3*(hiyij0-hiyij1)*invrij
      hiyik0=2.0d0*(yij*invrij-yik*cosjik*invrik)
      hiyik1=gam*yik*cosjik3*(invrika*invrika)
      hiyik=refi*cosjik3*(hiyik0-hiyik1)*invrik
      hjyij0=2.0d0*(-yjk*invrjk-yij*cosijk*invrij)
      hjyij1=gam*yij*cosijk3*(invrija*invrija)
      hjyij=refj*cosijk3*(hjyij0-hjyij1)*invrij
      hkyik0=2.0d0*(yjk*invrjk-yik*cosikj*invrik)
      hkyik1=gam*yik*cosikj3*(invrika*invrika)
      hkyik=refk*cosikj3*(hkyik0-hkyik1)*invrik
      hjyjk0=2.0d0*(-yij*invrij-yjk*cosijk*invrjk)
      hjyjk1=gam*yjk*cosijk3*(invrjka*invrjka)
      hjyjk=refj*cosijk3*(hjyjk0-hjyjk1)*invrjk
      hkykj0=2.0d0*(-yik*invrik+yjk*cosikj*invrjk)
      hkykj1=gam*yjk*cosikj3*(invrjka*invrjka)
      hkykj=refk*cosikj3*(hkykj0+hkykj1)*invrjk

      hizij0=2.0d0*(zik*invrik-zij*cosjik*invrij)
      hizij1=gam*zij*cosjik3*(invrija*invrija)
      hizij=refi*cosjik3*(hizij0-hizij1)*invrij
      hizik0=2.0d0*(zij*invrij-zik*cosjik*invrik)
      hizik1=gam*zik*cosjik3*(invrika*invrika)
      hizik=refi*cosjik3*(hizik0-hizik1)*invrik
      hjzij0=2.0d0*(-zjk*invrjk-zij*cosijk*invrij)
      hjzij1=gam*zij*cosijk3*(invrija*invrija)
      hjzij=refj*cosijk3*(hjzij0-hjzij1)*invrij
      hkzik0=2.0d0*(zjk*invrjk-zik*cosikj*invrik)
      hkzik1=gam*zik*cosikj3*(invrika*invrika)
      hkzik=refk*cosikj3*(hkzik0-hkzik1)*invrik
      hjzjk0=2.0d0*(-zij*invrij-zjk*cosijk*invrjk)
      hjzjk1=gam*zjk*cosijk3*(invrjka*invrjka)
      hjzjk=refj*cosijk3*(hjzjk0-hjzjk1)*invrjk
      hkzkj0=2.0d0*(-zik*invrik+zjk*cosikj*invrjk)
      hkzkj1=gam*zjk*cosikj3*(invrjka*invrjka)
      hkzkj=refk*cosikj3*(hkzkj0+hkzkj1)*invrjk

        fx3(i)=fx3(i)-hixij-hixik-hjxij-hkxik
        fy3(i)=fy3(i)-hiyij-hiyik-hjyij-hkyik
        fz3(i)=fz3(i)-hizij-hizik-hjzij-hkzik

        fx3(j)=fx3(j)+hixij+hjxij-hjxjk+hkxkj
        fy3(j)=fy3(j)+hiyij+hjyij-hjyjk+hkykj
        fz3(j)=fz3(j)+hizij+hjzij-hjzjk+hkzkj

        fx3(k)=fx3(k)+hixik+hkxik-hkxkj+hjxjk
        fy3(k)=fy3(k)+hiyik+hkyik-hkykj+hjyjk
        fz3(k)=fz3(k)+hizik+hkzik-hkzkj+hjzjk
   40   continue
   30   continue
end subroutine

