!-------------------------------------------------------------------
    ! function  : vli3_3b
    !
    ! package   : Li3 Potential
    !
    ! Language  : Fortran 77
    !
    ! author    : F. Colavecchia (flavioc@lanl.gov)
    !             G. A. Parker   (parker@nhn.ou.edu)
    !
    ! date      : 05/10/02       version: 0.1
    ! revision  :                version:
    !
    ! purpose   :  Compute the three body term of the spin aligned
    !              Li3 potential
    !
    ! input     :  r  -> interatomic distances in atomic units  
    !
    ! output    :  vli3_3b -> three body term  
    !
    ! modules   :
    !
    !
    ! common    :
    !
    !
    ! notes     :  Potential is build upon an interpolation of
    !              ab initio data for values of the scaled perimeter (s1) less
    !              than pswitch, and a extrapolation for greater values
    !              of the scaled perimeter. Asymptotic form corresponds
    !              to that of the Axilrod-Teller potential for Li3
    !              Needs: parameters.dat
    !                     li3_sc_cyclic.dat
    !                     coef_li3_sc_cyclic.dat
    !              in the inter and extrapolation routines.
    !
    !
    !
    !
    !-------------------------------------------------------------------
      real*8 function vli3_3b(ri)
      implicit none
      logical search,check,is_in
      integer i,num,j,ne
      integer r2s
      real*8 r(3),ri(3)      
      real*8 v,s1,s2,s3,ee,e2,rr,dummy
      real*8 pswitch,rho,vli3_at
      real*8 delper, s1switch, vinterp, vextrap
      data check /.false./
      save ne,pswitch,rr     
!      save firstcall,xe,ne,pswitch,rr
!
!     Compute the three body part of the Li3 potential
!
!     the requirement is that r(3) are in bohr atomic units
!     However, see interpolate.f
!
      vli3_3b = 0d0 
      r(1) = ri(1) 
      r(2) = ri(2) 
      r(3) = ri(3) 
      rho = dsqrt(dsqrt(3d0)/3d0*(r(1)**2+r(2)**2+r(3)**2))
      !
      !   get the s1,s2,s3
      !
      i = r2s(r(1),r(2),r(3),s1,s2,s3,.true.)
      !
      !   Optimal r value
      !
      rr=1.08d0
      !
      !     The triangular coordinates are defined as:
      !           
      !      s1=(r(1)+r(2)+r(3))/dsqrt(3.0d0)
      !      s2=(r(2)-r(3))/dsqrt(2.0d0)
      !      s3=(2.0d0*r(1)-r(2)-r(3))/dsqrt(6.0d0)
      !
      !     Now we see if the point (s1,s2,s3) lies in the interpolation
      !     region.
      !
      if(check) then
          is_in = search(s1,s2,s3)
          if(.not.is_in) return
      end if
      !
      !     Here we decide between inter or extrapolation. 
      !     See paper for details
      !
      delper=0.075d0
      pswitch=15.5d0

      s1switch=0.5d0*(1.d0+tanh((s1-pswitch)/(pswitch*delper)))
     
      call interpolate(s1,s2,s3,vinterp)
      call extrapolate(s1,s2,s3,vextrap)
      v=vinterp*(1.d0-s1switch)+(vextrap+vli3_at(r))*s1switch

      vli3_3b = v
      return
      end