!-------------------------------------------------------------------
    !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_ab_initio.dat
    !                     li3_coef.dat
    !
    !-------------------------------------------------------------------
      subroutine interpolate(s1,s2,s3,pot)
      implicit none
      real*8 s1,s2,s3,pot

      integer j,Num,co,i
      logical firstcall
      real*8 x(1200),y(1200),z(1200),g(1200)
      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 
      r = 1.08d0
      Num=1120
      if(firstcall) then
          firstcall=.false.
          !
          !     Read the data of interpolation files
          !
          open(unit=1200,file='li3_ab_initio.dat',status='old')
          open(unit=1300,file='li3_coef.dat',status='old')
          do i=1,Num
              read(1200,*) rho,x(i),y(i),z(i),e2
              read(1300,*) j,g(j),dummy
          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

c      write(*,200) ds1,ds2,ds3,d2s1s1,d2s2s2,d2s3s3,d2s1s2,
c     +     d2s1s3,d2s2s3
 200  format(9d9.3)
      end