!-------------------------------------------------------------------
    !subroutine : interpolate
    !
    ! package   : Li3 Potential
    !
    ! Language  : Fortran 77
    !
    ! author    : M. Salazar (msalazar@uu.edu) 
    !             F. Colavecchia (flavioc@lanl.gov)
    !
    ! date      : 05/10/02       version:
    ! revision  :                version:
    !
    ! purpose   :  Compute the three body term of the spin aligned
    !              Li3 potential by  interpolating ab-initio data
    !
    ! input     :  s1  -> s1 symmetric triangular coordinate in atomic units  
    !              s2  -> s2 symmetric triangular coordinate in atomic units  
    !              s3  -> s3 symmetric triangular coordinate in atomic units  
    !
    ! output    :  pot -> three body term  
    !
    ! modules   :
    !
    ! common    :
    !
    ! notes     :  Needs: 'li3_sc_cyclic.dat'
    !                     'coef_li3_sc_cyclic.dat'
    !
    !-------------------------------------------------------------------
      subroutine interpolate(s1i,s2i,s3i,pot)
      implicit none
      real*8 s1i,s2i,s3i,pot
       
      real*8 s1,s2,s3

      integer j,Num,co,i
      logical firstcall
      real*8 x(6200),y(6200),z(6200),g(6200)
      real*8 r,rt,rho,e2,ds1,ds2,ds3
      real*8 d2s1s1,d2s2s2,d2s3s3
      real*8 d2s1s2,d2s1s3,d2s2s3
      real*8 dummy
      data firstcall /.true./
      save 
!
!     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
!
      s1 = s1i*0.529177d0
      s2 = s2i*0.529177d0
      s3 = s3i*0.529177d0
      r = 0.02d0
      Num=6124
      if(firstcall) then
          firstcall=.false.
          !
          !     Read the data of interpolation files
          !
          open(unit=1200,file='./Potentials/li3web/li3_sc_cyclic.dat'
     + ,status='old')
          open(unit=1300,file='./Potentials/li3web/coef_li3_sc_cyclic
     +.dat',status='old')
          do i=1,Num
              read(1200,*) x(i),y(i),z(i),dummy,e2
              read(1300,*) j,g(j)
          end do
          close(unit=1200)
          close(unit=1300)
      end if



      pot=0.0d0
      ds1=0.0d0
      ds2=0.0d0
      ds3=0.0d0
      d2s1s1=0.0d0
      d2s2s2=0.0d0
      d2s3s3=0.0d0
      d2s1s2=0.0d0
      d2s1s3=0.0d0
      d2s2s3=0.0d0

      do 100 j=1,Num
         rt=dsqrt((s1-x(j))**2 + (s2-y(j))**2 +
     +        (s3-z(j))**2 + r)

         pot=pot+g(j)*rt
         if(r.eq.0.0d0) then
            co=0
            if (dabs(s1-x(j)).lt.0.1d-4) then
               co=co+1
            endif
            if (dabs(s2-y(j)).lt.0.1d-4) then
               co=co+1
            endif
            if (dabs(s3-z(j)).lt.0.1d-4) then
               co=co+1
            endif
            if (co.eq.3) then
               goto 100
            endif
         endif
         ds1=ds1 + g(j)*(s1-x(j))/rt
         ds2=ds2 + g(j)*(s2-y(j))/rt
         ds3=ds3 + g(j)*(s3-z(j))/rt

         d2s1s1=d2s1s1 + g(j)*(1.0d0/rt - (s1-x(j))**2/rt**3)
         d2s2s2=d2s2s2 + g(j)*(1.0d0/rt - (s2-y(j))**2/rt**3)
         d2s3s3=d2s3s3 + g(j)*(1.0d0/rt - (s3-z(j))**2/rt**3)
         d2s1s2=d2s1s2 - g(j)*((s1-x(j))*(s2-y(j))/rt**3)
         d2s1s3=d2s1s3 - g(j)*((s1-x(j))*(s3-z(j))/rt**3)
         d2s2s3=d2s2s3 - g(j)*((s2-y(j))*(s3-z(j))/rt**3)
 100  continue

 200  format(9d9.3)
      end