!-------------------------------------------------------------------
    ! 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:
    ! 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 perimeter less
    !              than pswitch, and a extrapolation for greater values
    !              of the perimeter. Asymptotic form corresponds
    !              to that of the axilrod-teller potential for Li3
    !              Needs: parameters.dat
    !                     li3_ab_initio.dat
    !                     li3_coef.dat
    !
    !
    !-------------------------------------------------------------------
      real*8 function vli3_3b(ri)
      implicit none
      logical firstcall,search,check,is_in,inter,debug
      integer i,num,j,ne
      integer r2s
      real*8 r(3),ri(3)
      real*8 x(6200),y(6200),z(6200),g(6200)
      real*8 xe(17)
      real*8 v,s1,s2,s3,ee,e2,rr,dummy
      real*8 pswitch,rho,vli3_at
      real*8 delper, s1switch, vinterp, vextrap, perimeter
      data firstcall /.true./
      data check /.false./      
      data debug /.false./
!      common /interpol/ inter       ! for test only.
      save x,y,z,g,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,
!     but the Li3 data we have are in Angtroms, so we
!     need to convert them
!
      vli3_3b = 0d0 !default value if the interpolation is an extrapolation
      r(1) = ri(1) 
      r(2) = ri(2) 
      r(3) = ri(3) 

      if(firstcall) then 
          firstcall=.false.
          !
          ! Read the pswitch value
          !
          !          open(unit=1200,file='rholim.dat',status='old')
          !          read(1200,*) pswitch
          !          close(1200)
      end if
!
!      print *,'olvier he0'
!     Override inter value.
!     This rho is valid only for X3 symmetries
!
      rho = dsqrt(dsqrt(3d0)/3d0*(r(1)**2+r(2)**2+r(3)**2))
      inter = .false.
      if (rho.le.pswitch) inter = .true.
      !
      !   get the s1,s2,s3
      !
      i = r2s(r(1),r(2),r(3),s1,s2,s3,.true.)
      !
      !   Optimal r value
      !
      rr=1.08d0
!      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.
!
      delper=0.075d0
      pswitch=23.4d0
      perimeter=r(1)+r(2)+r(3)
!     
      s1switch=0.5d0*(1.d0+tanh((perimeter-pswitch)/(pswitch*delper)))
      !write(*,*) s1switch
      call interpolate_old(s1,s2,s3,v)
      vinterp=v
      call extrapolate(s1,s2,s3,v)
      vextrap=v
      v=vinterp*(1.d0-s1switch)+(vextrap)*s1switch
      v=vinterp*(1.d0-s1switch)+(vextrap+vli3_at(r))*s1switch
!
!     Another switch to ATM
!
!      delper=0.1d0
!      pswitch=80d0
!      s1switch=0.5d0*(1.d0+tanh((perimeter-pswitch)/(pswitch*delper)))      
!      v=v*(1.d0-s1switch)+(v+vli3_at(r))*s1switch
!      write(*,100) r(1),r(2),r(3),s1,s2,s3,v
!      write(*,100) rho,perimeter,vinterp,vextrap,s1switch
! 100  format(7f11.6)
      vli3_3b = v
      return
      end