      SUBROUTINE  AKIMA (NP,ND,XD,YD,NI,XI, YI)
C
C PREVIOUSLY THIS ROUTINE WAS CALLED UVIP3P AND WAS IN SINGLE PRECISION
C IT WAS DOWNLOADED FROM GAMS AT NIST ON OCT 21 1998 
C IT PROVIDES SLIGHTLY BETTER FITTING THAN THE IMSL ROUTINE DCSAKM
C THIS IS A SUBJECTIVE COMMENT - IT DOES HOWEVER PREVENT
C OSCILLATIONS BETWEEN SPARCE DATA POINTS
C
C Univariate Interpolation (Improved Akima Method)
C
C Hiroshi Akima
C U.S. Department of Commerce, NTIA/ITS
C Version of 89/07/04
C
C This subroutine performs univariate interpolation.  It is based
C on the improved A method developed by Hiroshi Akima, 'A method
C of univariate interpolation that has the accuracy of a third-
C degree polynomial,' ACM TOMS, vol. xx, pp. xxx-xxx, 19xx.  (The
C equation numbers referred to in the comments below are those in
C the paper.)
C
C In this method, the interpolating function is a piecewise
C function composed of a set of polynomials applicable to
C successive intervals of the given data points.  This method
C uses third-degree polynomials as the default, but the user has
C an option to use higher-degree polynomial to reduce undulations
C in resulting curves.
C
C This method has the accuracy of a third-degree polynomial if
C the degree of the polynomials for the interpolating function is
C set to three.
C
C The input arguments are
C   NP = degree of the polynomials for the interpolating
C        function,
C   ND = number of input data points
C        (must be equal to 2 or greater),
C   XD = array of dimension ND, containing the abscissas of
C        the input data points
C        (must be in a monotonic increasing order),
C   YD = array of dimension ND, containing the ordinates of
C        the input data points,
C   NI = number of points for which interpolation is desired
C        (must be equal to 1 or greater),
C   XI = array of dimension NI, containing the abscissas of
C        the desired points.
C
C The output argument is
C   YI = array of dimension NI, where the ordinates of the
C        desired points are to be stored.
C
C If an integer value smaller than 3 is given to the NP argument,
C this subroutine assumes NP = 3.
C
C The XI array elements need not be monotonic, but this
C subroutine interpolates faster if the XI array elements are
C given in a monotonic order.
C
C If the XI array element is less than XD(1) or greater than
C XD(ND), this subroutine linearly interpolates the YI value.
C
C
C Specification statement
      implicit real*8 (a-h,o-z)
      DIMENSION   XD(*),YD(*),XI(*), YI(*)
C Error check
   10 IF (ND.LE.1)   GO TO 90
      IF (NI.LE.0)   GO TO 91
      DO 11  ID=2,ND
        IF (XD(ID).LE.XD(ID-1))     GO TO 92
   11 CONTINUE
C Branches off special cases.
      IF (ND.LE.4)   GO TO 50
C General case  --  Five data points of more
C Calculates some local variables.
   20 NP0=MAX(3,NP)
      NPM1=NP0-1
      RENPM1=NPM1
      RENNM2=NP0*(NP0-2)
C Main calculation for the general case
C First (outermost) DO-loop with respect to the desired points
   30 DO 39  II=1,NI
        IF (II.EQ.1)      IINTPV=-1
        XII=XI(II)
C Locates the interval that includes the desired point by binary
C search.
        IF (XII.LE.XD(1))  THEN
          IINT=0
        ELSE IF (XII.LT.XD(ND))  THEN
          IDMN=1
          IDMX=ND
          IDMD=(IDMN+IDMX)/2
   31     IF (XII.GE.XD(IDMD))  THEN
            IDMN=IDMD
          ELSE
            IDMX=IDMD
          END IF
          IDMD=(IDMN+IDMX)/2
          IF (IDMD.GT.IDMN)    GO TO 31
          IINT=IDMD
        ELSE
          IINT=ND
        END IF
C End of locating the interval of interest
C Interpolation or extrapolation in one of the three subcases
        IF (IINT.LE.0)  THEN
C Subcase 1  --  Linear extrapolation when the abscissa of the
C                desired point is equal to that of the first data
C                point or less.
C Estimates the first derivative when the interval is not the
C same as the one for the previous desired point.  --
C cf. Equation (8)
          IF (IINT.NE.IINTPV)  THEN
            IINTPV=IINT
            X0=XD(1)
            X1=XD(2)-X0
            X2=XD(3)-X0
            X3=XD(4)-X0
            Y0=YD(1)
            Y1=YD(2)-Y0
            Y2=YD(3)-Y0
            Y3=YD(4)-Y0
            DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1)
            A1=(((X2*X3)**2)*(X3-X2)*Y1
     1         +((X3*X1)**2)*(X1-X3)*Y2
     2         +((X1*X2)**2)*(X2-X1)*Y3)/DLT
          END IF
C Evaluates the YI value.
          YI(II)=Y0+A1*(XII-X0)
C End of Subcase 1
        ELSE IF (IINT.GE.ND)  THEN
C Subcase 2  --  Linear extrapolation when the abscissa of the
C                desired point is equal to that of the last data
C                point or greater.
C Estimates the first derivative when the interval is not the
C same as the one for the previous desired point.  --
C cf. Equation (8)
          IF (IINT.NE.IINTPV)  THEN
            IINTPV=IINT
            X0=XD(ND)
            X1=XD(ND-1)-X0
            X2=XD(ND-2)-X0
            X3=XD(ND-3)-X0
            Y0=YD(ND)
            Y1=YD(ND-1)-Y0
            Y2=YD(ND-2)-Y0
            Y3=YD(ND-3)-Y0
            DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1)
            A1=(((X2*X3)**2)*(X3-X2)*Y1
     1         +((X3*X1)**2)*(X1-X3)*Y2
     2         +((X1*X2)**2)*(X2-X1)*Y3)/DLT
          END IF
C Evaluates the YI value.
          YI(II)=Y0+A1*(XII-X0)
C End of Subcase 2
        ELSE
C Subcase 3  --  Interpolation when the abscissa of the desired
C                point is  between those of the first and last
C                data points.
C Calculates the coefficients of the third-degree polynomial (for
C NP.LE.3) or the factors for the higher-degree polynomials (for
C NP.GT.3), when the interval is not the same as the one for the
C previous desired point.
          IF (IINT.NE.IINTPV)  THEN
            IINTPV=IINT
C The second DO-loop with respect to the two endpoints of the
C interval
            DO 37  IEPT=1,2
C Calculates the estimate of the first derivative at an endpoint.
C Initial setting for calculation
              ID0=IINT+IEPT-1
              X0=XD(ID0)
              Y0=YD(ID0)
              SMPEF=0.0d0
              SMWTF=0.0d0
              SMPEI=0.0d0
              SMWTI=0.0d0
C The third (innermost) DO-loop with respect to the four primary
C estimate of the first derivative
              DO 36  IPE=1,4
C Selects point numbers of four consecutive data points for
C calculating the primary estimate of the first derivative.
                IF (IPE.EQ.1)  THEN
                  ID1=ID0-3
                  ID2=ID0-2
                  ID3=ID0-1
                ELSE IF (IPE.EQ.2)  THEN
                  ID1=ID0+1
                ELSE IF (IPE.EQ.3)  THEN
                  ID2=ID0+2
                ELSE
                  ID3=ID0+3
                END IF
C Checks if any point number falls outside the legitimate range
C (between 1 and ND).  Skips calculation of the primary estimate
C if any does.
                IF (ID1.LT.1.OR.ID2.LT.1.OR.ID3.LT.1.OR.
     1              ID1.GT.ND.OR.ID2.GT.ND.OR.ID3.GT.ND)
     2               GO TO 36
C Calculates the primary estimate of the first derivative  --
C cf. Equation (8)
                X1=XD(ID1)-X0
                X2=XD(ID2)-X0
                X3=XD(ID3)-X0
                Y1=YD(ID1)-Y0
                Y2=YD(ID2)-Y0
                Y3=YD(ID3)-Y0
                DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1)
                PE=(((X2*X3)**2)*(X3-X2)*Y1
     1             +((X3*X1)**2)*(X1-X3)*Y2
     2             +((X1*X2)**2)*(X2-X1)*Y3)/DLT
C Calculates the volatility factor, VOL, and distance factor,
C SXX, for the primary estimate.  --  cf. Equations (9) and (11)
                SX=X1+X2+X3
                SY=Y1+Y2+Y3
                SXX=X1*X1+X2*X2+X3*X3
                SXY=X1*Y1+X2*Y2+X3*Y3
                DNM=4.0d0*SXX-SX*SX
                B0=(SXX*SY-SX*SXY)/DNM
                B1=(4.0d0*SXY-SX*SY)/DNM
                DY0=-B0
                DY1=Y1-(B0+B1*X1)
                DY2=Y2-(B0+B1*X2)
                DY3=Y3-(B0+B1*X3)
                VOL=DY0*DY0+DY1*DY1+DY2*DY2+DY3*DY3
C Calculates the EPSLN value, which is used to decide whether or
C not the volatility factor, VOL, is essentially zero.
                EPSLN=(YD(ID0)**2+YD(ID1)**2
     1                +YD(ID2)**2+YD(ID3)**2)*1.0D-12
C Accumulates the weighted primary estimates.  --
C cf. Equations (13) and (14)
                IF (VOL.GT.EPSLN)  THEN
C - For finite weight.
                  WT=1.0d0/(VOL*SXX)
                  SMPEF=SMPEF+PE*WT
                  SMWTF=SMWTF+WT
                ELSE
C - For infinite weight.
                  SMPEI=SMPEI+PE
                  SMWTI=SMWTI+1.0d0
                END IF
   36         CONTINUE
C End of the third DO-loop
C Calculates the final estimate of the first derivative.  --
C cf. Equation (14)
              IF (SMWTI.LT.0.5d0)  THEN
C - When no infinite weights exist.
                YP=SMPEF/SMWTF
              ELSE
C - When infinite weights exist.
                YP=SMPEI/SMWTI
              END IF
              IF (IEPT.EQ.1)  THEN
                YP0=YP
              ELSE
                YP1=YP
              END IF
C End of the calculation of the estimate of the first derivative
C at an endpoint
   37       CONTINUE
C End of the second DO-loop
            IF (NP0.LE.3)  THEN
C Calculates the coefficients of the third-degree polynomial
C (when NP.LE.3).  --  cf. Equation (4)
              DX=XD(IINT+1)-XD(IINT)
              DY=YD(IINT+1)-YD(IINT)
              A0=YD(IINT)
              A1=YP0
              YP1=YP1-YP0
              YP0=YP0-DY/DX
              A2=-(3.0d0*YP0+YP1)/DX
              A3= (2.0d0*YP0+YP1)/(DX*DX)
            ELSE
C Calculates the factors for the higher-degree polynomials
C (when NP.GT.3).  --  cf. Equation (20)
              DX=XD(IINT+1)-XD(IINT)
              DY=YD(IINT+1)-YD(IINT)
              T0=YP0*DX-DY
              T1=YP1*DX-DY
              AA0= (T0+RENPM1*T1)/RENNM2
              AA1=-(RENPM1*T0+T1)/RENNM2
            END IF
          END IF
C End of the calculation of the coefficients of the third-degree
C polynomial (when NP.LE.3) or the factors for the higher-degree
C polynomials (when NP.GT.3), when the interval is not the same
C as the one for the previous desired point.
C Evaluates the YI value.
          IF (NP0.LE.3)  THEN
C - With a third-degree polynomial.  --  cf. Equation (3)
            XX=XII-XD(IINT)
            YI(II)=A0+XX*(A1+XX*(A2+XX*A3))
          ELSE
C - With a higher-degree polynomial.  --  cf. Equation (19)
            U=(XII-XD(IINT))/DX
            UC=1.0d0-U
            V=AA0*((U**NP0)-U)+AA1*((UC**NP0)-UC)
            YI(II)=YD(IINT)+DY*U+V
          END IF
C End of Subcase 3
        END IF
   39 CONTINUE
C End of the first DO-loop
C End of general case
      RETURN
C Special cases  --  Four data points or less
C Preliminary processing for the special cases
   50 X0=XD(1)
      Y0=YD(1)
      X1=XD(2)-X0
      Y1=YD(2)-Y0
      IF (ND.EQ.2)   GO TO 60
      X2=XD(3)-X0
      Y2=YD(3)-Y0
      IF (ND.EQ.3)   GO TO 70
      X3=XD(4)-X0
      Y3=YD(4)-Y0
      GO TO 80
C Special Case 1  --  Two data points
C (Linear interpolation and extrapolation)
   60 A1=Y1/X1
      DO 61  II=1,NI
        YI(II)=Y0+A1*(XI(II)-X0)
   61 CONTINUE
C End of Special Case 1
      RETURN
C Special Case 2  --  Three data points
C (Quadratic interpolation and linear extrapolation)
   70 DLT=X1*X2*(X2-X1)
      A1=(X2*X2*Y1-X1*X1*Y2)/DLT
      A2=(X1*Y2-X2*Y1)/DLT
      A12=2.0d0*A2*X2+A1
      DO 71  II=1,NI
        XX=XI(II)-X0
        IF (XX.LE.0.0d0)  THEN
          YI(II)=Y0+A1*XX
        ELSE IF (XX.LT.X2) THEN
          YI(II)=Y0+XX*(A1+XX*A2)
        ELSE
          YI(II)=Y0+Y2+A12*(XX-X2)
        END IF
   71 CONTINUE
C End of Special Case 2
      RETURN
C Special Case 3  --  Four data points
C (Cubic interpolation and linear extrapolation)
   80 DLT=X1*X2*X3*(X2-X1)*(X3-X2)*(X3-X1)
      A1=(((X2*X3)**2)*(X3-X2)*Y1
     1   +((X3*X1)**2)*(X1-X3)*Y2
     2   +((X1*X2)**2)*(X2-X1)*Y3)/DLT
      A2=(X2*X3*(X2*X2-X3*X3)*Y1
     1   +X3*X1*(X3*X3-X1*X1)*Y2
     2   +X1*X2*(X1*X1-X2*X2)*Y3)/DLT
      A3=(X2*X3*(X3-X2)*Y1
     1   +X3*X1*(X1-X3)*Y2
     2   +X1*X2*(X2-X1)*Y3)/DLT
      A13=(3.0d0*A3*X3+2.0d0*A2)*X3+A1
      DO 81  II=1,NI
        XX=XI(II)-X0
        IF (XX.LE.0.0d0)  THEN
          YI(II)=Y0+A1*XX
        ELSE IF (XX.LT.X3) THEN
          YI(II)=Y0+XX*(A1+XX*(A2+XX*A3))
        ELSE
          YI(II)=Y0+Y3+A13*(XX-X3)
        END IF
   81 CONTINUE
C End of Special Case 3
      RETURN
C Error exit
   90 WRITE (*,99090) ND
      GO TO 99
   91 WRITE (*,99091) NI
      GO TO 99
   92 WRITE (*,99092) ID,XD(ID-1),XD(ID)
   99 WRITE (*,99099)
      RETURN
C Format statements for error messages
99090 FORMAT (1X/ ' ***   Insufficient data points.'
     1  7X,'ND =',I3)
99091 FORMAT (1X/ ' ***   No desired points.'
     1  7X,'NI =',I3)
99092 FORMAT (1X/ ' ***   Two data points identical or out of ',
     1  'sequence.'/
     2  7X,'ID, XD(ID-1), XD(ID) =',I5,2F10.3)
99099 FORMAT (' Error detected in the akima subroutine'/)
      END
      !-------------------------------------------------------------------
      ! subroutine: mpot
      !
      ! package   : li3 Potentials
      !
      ! Language  : Fortran 77
      !
      ! author    : V. Venturi (vanessa.venturi@nist.gov)
      !
      ! revision  : F. Colavecchia (flavioc@lanl.gov)
      !
      ! purpose   :  Get the potential of the singlet and triplet Li2,
      !             once the data and parameters are read at syspot3;
      !             according to eq. (8) of  JCP 118 5484 (2003).
      !
      ! input     :   np -> degree of polynomial used in the akima interpolation. 
      !            ndata -> number of ab initio+rkr+exchange data points
      !            rdata -> distance in a.u. for each data point
      !            vdata -> potentials at each data point.
      !            c6    -> data for Van der Waals form, eq. (3)
      !            c8    -> data for Van der Waals form, eq. (3)
      !            c10   -> data for Van der Waals form, eq. (3) 
      !            cA    -> data for exchange form, eq. (6)
      !            beta  -> data for exchange form, eq. (6)
      !            gamma -> data for exchange form, eq. (6)
      !            a     -> data for short range form, eq. (7)
      !            b     -> data for short range form, eq. (7)
      !            nv    -> singlet (nv=1) or triplet (nv=2) potential
      !            r     -> internuclear distance in a.u.                   
      !
      ! output    :    v -> li2 potential in a.u.
      !
      !
      ! notes     : The original code included RKR data from table 3
      !             (Linton et al J. Mol. Spect. 196, 20 (1999)).
      !             Now we make use of new RKR data provided by A. Ross
      !             (ross@in2p3.fr), read from rkr.dat. See also
      !             readme_li2.txt for the values of De, Do, a, rmin.
      !             Comments noted by CFC were included by F. Colavecchia
      !
      !-------------------------------------------------------------------
      subroutine mpot(np,ndata,rdata,vdata,c6,c8,c10,
     >                CA,beta,gamma,a,b,nv,r,v)
C
      double precision rdata,vdata,c6,c8,c10,CA,beta,gamma,a,b,r,v
      integer  np,ndata,nv
      dimension rdata(ndata),vdata(ndata)
C
C ... Local Variables
C  
      double precision vex,vdisp,rfinal

      rfinal = rdata(ndata)
C
      if(r.ge.0.0d0 .and. r.lt.rdata(1)) then
          !
          !     Extrapolate short range data using the form: aexp(-br)
          !	
          v=a*exp(-b*r)         
      elseif(r.ge.rdata(1) .and. r.lt.rfinal) then
          !
          !     Interpolate short range data 
          !	
          call akima(np,ndata,rdata,vdata,1,r,v)
          v=v/r**6 
      elseif(r.ge.rfinal) then
          !
          !     Calculate long range form
          !
          vex=CA*exp(-beta*r)*r**gamma 
          vdisp=-c6/r**6-c8/r**8-c10/r**10
          if(nv.eq.1) then
            v=vdisp-vex
          elseif(nv.eq.2) then
              v=vdisp+vex
          endif  
          
      endif 
C
      return
      end 
      !-------------------------------------------------------------------
      ! function  : correction_ss
      !
      ! package   : li3 Potentials
      !
      ! Language  : Fortran 77
      !
      ! author    : V. Venturi (vanessa.venturi@nist.gov)
      !
      ! revision  : F. Colavecchia (flavioc@lanl.gov)
      !
      ! purpose   : Harmonic correction for the potential.
      !
      !
      ! input     :    r    -> interatomic distance in bohr.
      !               shift -> harmonic constant of the correction
      !                re   -> distance at the minimum of the potential
      !
      ! output    :    v -> correction potential
      !
      !-------------------------------------------------------------------      
      double precision function correction_ss(r,shift,re)
      implicit none
      double precision r,shift,c6,re

      correction_ss=0.0d0
      if (r.lt.re) correction_ss=shift*(r-re)**2

      return
      end
      !-------------------------------------------------------------------
      ! function  : sort
      !
      ! package   : li3 Potentials
      !
      ! Language  : Fortran 77
      !
      ! author    : V. Venturi (vanessa.venturi@nist.gov)
      !
      ! purpose   : sorting vector a and b in ascending order
      !-------------------------------------------------------------------      
      subroutine sort(a,b,n)
      implicit none
      integer n,i,j
      real*8 a(n),atem,b(n),btem
      do 10 j=1,n
         do 20 i=1,n-j
            if (a(i).gt.a(i+1)) then
               atem=a(i+1)
               a(i+1)=a(i)
               a(i)=atem
               btem=b(i+1)
               b(i+1)=b(i)
               b(i)=btem
            end if
20       continue
10    continue
      return
      end
      !-------------------------------------------------------------------
      ! function  : sort
      !
      ! package   : li3 Potentials
      !
      ! Language  : Fortran 77
      !
      ! author    : V. Venturi (vanessa.venturi@nist.gov)
      !
      ! purpose   : sorting vector a,b and c in ascending order
      !-------------------------------------------------------------------      
      subroutine sort2(a,b,c,n)
      implicit none
      integer n,i,j
      real*8 a(n),atem,b(n),btem,c(n),ctem
      do 10 j=1,n
         do 20 i=1,n-j
            if (a(i).gt.a(i+1)) then
               atem=a(i+1)
               a(i+1)=a(i)
               a(i)=atem
               btem=b(i+1)
               b(i+1)=b(i)
               b(i)=btem
               ctem=c(i+1)
               c(i+1)=c(i)
               c(i)=ctem
            end if
20       continue
10    continue
      return
      end
    !-------------------------------------------------------------------
    ! function  : syspot2
    !
    ! package   : li3 Potentials
    !
    ! Language  : Fortran 77
    !
    ! author    : V. Venturi (vanessa.venturi@nist.gov)
    !
    ! revision  : F. Colavecchia (flavioc@lanl.gov) 
    !
    ! purpose   :  Get the potential of the singlet and triplet Li2.
    !
    !
    ! input     :    r    -> interatomic distance in bohr.
    !               nv    -> singlet (nv=1) or triplet (nv=2) potential
    !              lsetup -> sets up the reading of data or calculation
    !                        of potential
    !
    ! output    :    v -> li2 potential in a.u.  
    !             
    !
    ! modules   :
    !
    ! common    :
    !
    ! notes     : The original code included RKR data from table 3
    !             (Linton et al J. Mol. Spect. 196, 20 (1999)).
    !             Now we make use of new RKR data provided by A. Ross
    !             (ross@in2p3.fr), read from rkr.dat. See also
    !             readme_li2.txt for the values of De, Do, a, rmin.
    !             Comments noted by CFC were included by F. Colavecchia
    !
    !-------------------------------------------------------------------
      DOUBLE PRECISION FUNCTION SYS_SS_pot_3(R,V,NV,LSETUP)
      IMPLICIT NONE
      DOUBLE PRECISION R,V
      INTEGER NV
      LOGICAL LSETUP
C
C ... Local Variables
C
      double precision c6,c8,c10,CA,beta,gamma,as,bs,at,bt,
     >  re1,re3,shift1,shift3,rdatas,vdatas,rdatat,vdatat
      save c6,c8,c10,CA,beta,gamma,as,bs,at,bt,
     >  re1,re3,shift1,shift3,rdatas,vdatas,rdatat,vdatat
      double precision hartree,bohra
      save  hartree,bohra
      double precision correction_ss
      double precision rs1,vs1,rs2,vs2,rs3,vs3,rt0,vt0,rt1,vt1,rt2,vt2,
     >  rt3,vt3,res,des,ret,det,rdat,vdats,vdatt
      double precision  vs,vt,vav,vdiff,vex,vdisp,amp,vshift,
     >  rfit,rjoin,rfinal,r0,rl,ri,r1,r2,v1,v2,delta
      
      integer np,ns,nt,ndatas,ndatat
      save np,ns,nt,ndatas,ndatat
      integer  i,j,k,ni,nsl,ntl,ndat,nt2off
C
      integer  ns1,ns2,ns3,nt1,nt2,nt3,nt0
      integer nmax
      parameter (nmax=500)
      
      parameter (ns1=3,ns2=42,ns3=41,nt1=11,nt0=11,nt2=27,nt3=26)
      
C
      real*8 vdatatp(nmax),vdatasp(nmax)
      dimension rs1(ns1),vs1(ns1),rs2(ns2),vs2(ns2),rs3(ns3),vs3(ns3),
     >  rt1(nt1),vt1(nt1),rt2(nt2),vt2(nt2),rt3(nt3),vt3(nt3)
      dimension rdatas(nmax),vdatas(nmax),rdatat(nmax),vdatat(nmax),
     >  rdat(nmax),vdats(nmax),vdatt(nmax)
      dimension vt0(nt0),rt0(nt0)
      logical do_setup
      data do_setup /.true./
      save  do_setup
      logical order,norepeat
C
C     output flags for potential routine:
C     output_data  - output raw potential data
C     output_fitex - exchange fitting constant A as a function of R 
C     output_v     - output final potentials to be used in calculations
C     output_vmesh - output interpolated mesh of final potentials +
C                    correction included
C
      logical  output_data, output_fitex, output_v, output_vmesh
      data output_data /.false./
      data output_fitex /.false./
      data output_v /.false./
      data output_vmesh /.false./
C
C=====Data for ^7Li2====================================================
C
C     Physical constants used in potential routines
C
      DATA HARTREE /219474.630680D0/, BOHRA /0.5291772490D0/
C
C     np is the degree of polynomial used in the akima spline interpolation
C
      np=3
C
C     rfit= point at which the exchange form is fit to the data and at
C     which is calculated the shift in the data required to match the
C     short and long range parts 
C     rjoin= point at which the short range data is smoothly connected
C     to the long range form
C     rfinal= last point for short range data set used in the calculations
C
      DATA rfit /15.5D0/, rjoin /17.D0/,rfinal /20.0d0/
C
C --- Long range: 
C     Yan, Babb, Dalgarno and Drake, Phys. Rev. A 54, 2824 (1996)
C     (in atomic units)
C
      DATA C6 /1393.39D0/, C8 /83425.8D0/, C10 /7372100.0D0/
C
C --- Long range exchange: 
C     Cote, Dalgarno, and Jamieson, Phys. Rev. A 50, 399 (1994)
C     (in atomic units)
C
      DATA  beta /1.25902D0/, gamma /4.55988D0/

C
C --- Short-range singlet potential data
C     ----------------------------------
C     De determined by Zemke and Stwalley, J. Phys. Chem. 97, 2053 (1993)
C     (in units of cm-1)   
C
      DATA Des /8516.61D0/
C
C     Section with E>0 : small R
C     1st 2 pts from Konowalow and Olson J. Chem. Phys. 71, 450 (1979)
C     3rd pt from Schmidt-Mink et al Chem. Phys. 92, 263 (1985)
C     (in atomic units)
C
      DATA rs1 / 2.75, 3.00, 3.25 /
C
      DATA vs1 / 0.061866D0, 0.036411D0, 0.013546D0 /
C
C     Section with E<0 : RKR potential
C     from Barakat et al Chem. Phys. 102, 215 (1986)  
C     note: misprints in data corrected, see Cote et al 
C           Phys. Rev A 50 399 (1994)
C     (in units of angstroms and cm-1)
C
      DATA res / 2.673247D0 /
C
      DATA rs2 / 1.823201D0, 1.823462D0, 1.824057D0,
     &           1.825118D0, 1.826727D0, 1.828918D0,
     &           1.831708D0, 1.834867D0, 1.838851D0,
     &           1.843462D0, 1.848639D0, 1.854376D0,
     &           1.860716D0, 1.867431D0, 1.874904D0,
     &           1.882864D0, 1.891368D0, 1.900497D0,
     &           1.910173D0, 1.920282D0, 1.931173D0,
     &           1.942717D0, 1.954984D0, 1.967947D0,
     &           1.981710D0, 1.996354D0, 2.011906D0,
     &           2.028431D0, 2.046073D0, 2.064963D0,
     &           2.085263D0, 2.107177D0, 2.130967D0,
     &           2.156977D0, 2.185675D0, 2.217724D0,
     &           2.254114D0, 2.296439D0, 2.347580D0,
     &           2.413867D0, 2.517207D0, 2.673247D0 /
C
      DATA vs2 / 8514.7777D0, 8507.8421D0, 8492.0440D0,
     &           8463.9645D0, 8421.6125D0, 8364.3068D0,
     &           8292.0294D0, 8205.2323D0, 8104.4730D0,
     &           7990.4162D0, 7863.7083D0, 7724.9765D0,
     &           7574.8736D0, 7413.8431D0, 7242.5556D0,
     &           7061.4199D0, 6870.8931D0, 6671.3979D0,
     &           6463.3140D0, 6246.9482D0, 6022.6578D0,
     &           5790.7056D0, 5551.3992D0, 5304.9322D0,
     &           5051.5343D0, 4791.4274D0, 4524.7756D0,
     &           4251.7309D0, 3972.4624D0, 3687.1094D0,
     &           3395.7978D0, 3090.6412D0, 2795.7419D0,
     &           2487.1914D0, 2173.0721D0, 1853.4573D0,
     &           1528.4128D0, 1197.9974D0,  862.2642D0,
     &            521.2611D0,  175.0320D0,    0.0000D0  /
C
      DATA rs3 / 12.639382D0, 9.863785D0, 8.495553D0,
     &           7.657218D0, 7.090991D0, 6.675227D0,
     &           6.352341D0, 6.089904D0, 5.870023D0,
     &           5.680896D0, 5.514781D0, 5.366746D0,
CX   &           5.580896D0, 5.514781D0, 5.366746D0,
     &           5.232506D0, 5.109851D0, 4.996566D0,
     &           4.890860D0, 4.791595D0, 4.697788D0,
CX   &           4.490860D0, 4.791595D0, 4.697788D0,
     &           4.608426D0, 4.522846D0, 4.440679D0,
     &           4.361370D0, 4.284400D0, 4.209387D0,
     &           4.135091D0, 4.064160D0, 3.993263D0,
     &           3.923163D0, 3.853594D0, 3.784295D0,
     &           3.714965D0, 3.645290D0, 3.574908D0,
     &           3.503385D0, 3.430172D0, 3.354534D0,
     &           3.275416D0, 3.191166D0, 3.098852D0,
     &           2.992096D0, 2.848951D0 /
C
      DATA vs3 / 8514.7777D0, 8507.8421D0, 8492.0440D0,
     &           8463.9645D0, 8421.6125D0, 8364.3068D0,
     &           8292.0294D0, 8205.2323D0, 8104.4730D0,
     &           7990.4162D0, 7863.7083D0, 7724.9765D0,
     &           7574.8736D0, 7413.8431D0, 7242.5556D0,
     &           7061.4199D0, 6870.8931D0, 6671.3979D0,
     &           6463.3140D0, 6246.9482D0, 6022.6578D0,
     &           5790.7056D0, 5551.3992D0, 5304.9322D0,
     &           5051.5343D0, 4791.4274D0, 4524.7756D0,
     &           4251.7309D0, 3972.4624D0, 3687.1094D0,
     &           3395.7978D0, 3090.6412D0, 2795.7419D0,
     &           2487.1914D0, 2173.0721D0, 1853.4573D0,
     &           1528.4128D0, 1197.9974D0,  862.2642D0,
     &            521.2611D0,  175.0320D0 /
C
C --- Short-range triplet potential data
C     ----------------------------------
C     De determined by Linton et al, J. Mol. Spect. 196, 20 (1999)
C     (in units of cm-1)
C
      DATA Det /333.70D0/
C
C     Section with E>0 : small R
C     1st pt from Konowalow and Olson J. Chem. Phys. 81, 4534 (1984)
C     7 pts from Schmidt-Mink et al Chem. Phys. 92, 263 (1985)
C     (in atomic units)
C
c      DATA rt1 / 3.00D0 , 3.25D0,  3.50D0,
c     &           4.00D0 , 4.50D0,  5.00D0,
c     &           5.50D0 , 6.00D0/
C
c      DATA vt1 / 0.104471D0, 0.082074D0, 0.063725D0,
c     &           0.036778D0, 0.020006D0, 0.010053D0,
c     &           0.004370D0, 0.001258D0 /  
C     
C    11 pts from Halls et al, Chem. Phys. Lett.  339, 427 (2001)        
C          distances in angstroms
c
      DATA rt1 / 2.00d0, 2.125d0, 2.25D0,  2.375D0,
     &           2.50d0, 2.625d0, 2.75d0,
     &           2.875d0, 3.00d0,  3.125d0, 3.25d0/

     
      DATA vt1 / 0.0472751972d0,
     &       0.0361788060d0, 0.0273239836d0, 0.0203412490d0,
     &       0.0148869076d0, 0.0106605276d0, 0.0074111640d0,
     &       0.0049337380d0, 0.0030628654d0, 0.0016656470d0,
     &       0.0006356372d0 /
C
C     Section with E<0 : RKR potential
C     from Linton et al J. Mol. Spect. 196, 20 (1999)  
C     (in units of angstroms and cm-1)
C
      DATA ret / 4.173D0 /
C
CFC      DATA rt2 / 3.356D0, 3.358D0, 3.365D0,
CFC     &           3.377D0, 3.395D0, 3.422D0,
CFC     &           3.458D0, 3.505D0, 3.571D0,
CFC     &           3.668D0, 3.846D0, 4.159D0,4.173D0 /
C
CFC      DATA vt2 / 333.269D0, 330.170D0, 322.155D0,
CFC     &           308.098D0, 287.665D0, 260.837D0,
CFC     &           227.679D0, 188.240D0, 142.523D0,
CFC     &           90.453D0,  31.857D0, 0.024D0,  0.000D0 /
C
CFC      DATA rt3 / 16.052D0, 11.392D0, 9.441D0,
CFC     &            8.297D0,  7.501D0,  6.885D0,
CFC     &            6.373D0,  5.922D0,  5.503D0,
CFC     &            5.092D0,  4.630D0 , 4.185D0/
C
CFC      DATA vt3 / 333.269D0, 330.170D0, 322.155D0,
CFC     &           308.098D0, 287.665D0, 260.837D0,
CFC     &           227.679D0, 188.240D0, 142.523D0,
CFC     &            90.453D0,  31.857D0,   0.024D0 /
C
C=============Initial setup=============================================
C
      IF(LSETUP) THEN 
          if (do_setup) then
          shift1=0.0d0
          shift3=0.13d-5
          c6=1393.39d0
CFC
!
!   Read new RKR points
!   
          open(unit=1200,file='./Potentials/ALLNEW/rkr.dat',
     +status='old')
          do i=1,nt2
              read(1200,*) rt2(i),vt2(i)
              !write(*,*) i,rt2(i),vt2(i)
          end do
          do i=1,nt3
              read(1200,*) rt3(i),vt3(i)
              !write(*,*) i,rt3(i),vt3(i)              
          end do
          close(1200)
          nt2off = 4   ! skip the first nt2off values of rkr data
       
C
C --- Calculate total number of potential data points
C
          ns=ns1+ns2+ns3 
          nt=nt1+nt2+nt3-nt2off
          if (ns.gt.nmax.or.nt.gt.nmax) then
              write(*,1000) ns,nt,nmax
          endif
C
C --- Order each set of data, in ascending order, into an array    
C     and smooth data by multiplying by R^6 for interpolation.
C     Convert data to atomic units.
C
          re1=res/bohra 
          re3=ret/bohra 
C
          do i=1,ns1      
              rdatas(i)=rs1(i) 
              vdatas(i)=vs1(i)
              vdatas(i)=vdatas(i)*rdatas(i)**6
          enddo
          do i=1,ns2
              j=ns1+i
              rdatas(j)=rs2(i)/bohra 
              vdatas(j)=(vs2(i)-des)/hartree
              vdatas(j)=vdatas(j)*rdatas(j)**6
          enddo
          do i=1,ns3
              j=ns1+ns2+i
              rdatas(j)=rs3(i)/bohra 
              vdatas(j)=(vs3(i)-des)/hartree
              vdatas(j)=vdatas(j)*rdatas(j)**6
          enddo
C   
          order=.true.     
          do i=2,ns
              if (rdatas(i).eq.rdatas(i-1)) then
                  write(*,1010) i-1,i
                  stop
              elseif (rdatas(i).lt.rdatas(i-1)) then
                  order=.false.
              endif
          enddo
          if (.not. order) then
              call sort (rdatas,vdatas,ns)
          endif
C
c          do i=1,nt0
c              rdatat(i)=rt0(i)/bohra 
c              vdatat(i)=vt0(i)
c              vdatat(i)=vdatat(i)*rdatat(i)**6
c          enddo
       
          do i=1,nt1
              j=i
              rdatat(j)=rt1(i)/bohra 
              vdatat(j)=vt1(i)
              vdatat(j)=vdatat(j)*rdatat(j)**6
          enddo
C
C         Note: 1st and 2nd points of the RKR data had to be omitted  
C               to prevent oscillations in the spline fit between
C               the two sets of data 
CFC
CFC      Now we omit the first nt2off points of the left side
CFC      of RKR data.
CFC
          do i=nt2off+1,nt2
              j=nt1-nt2off+i
              rdatat(j)=rt2(i)/bohra 
              vdatat(j)=(vt2(i)-det)/hartree
              vdatat(j)=vdatat(j)*rdatat(j)**6
          enddo
          do i=1,nt3
              j=nt1+nt2-nt2off+i
              rdatat(j)=rt3(i)/bohra 
              vdatat(j)=(vt3(i)-det)/hartree
              vdatat(j)=vdatat(j)*rdatat(j)**6
          enddo
C
          order=.false.     
          do i=2,nt
              if (rdatat(i).eq.rdatat(i-1)) then
                  write(*,*) 'Error in RKR data'
                  write(*,1010) i-1,i
                  write(*,'(2f20.8)') rdatat(i-1),vdatat(i-1)
                  write(*,'(2f20.8)') rdatat(i),vdatat(i)                  
                  stop
              elseif (rdatat(i).lt.rdatat(i-1)) then
                  order=.false.
              endif
          enddo
          if (.not. order) then
              call sort (rdatat,vdatat,nt)
          endif
C
          if (output_data) then
              do i=1,ns 
                  write(50,1100) rdatas(i),vdatas(i)/rdatas(i)**6
              enddo  
              do i=1,nt 
                  write(51,1100) rdatat(i),vdatat(i)/rdatat(i)**6
              enddo  
C
              do i=1,ns
                  rdat(i)=rdatas(i)
                  vdats(i)=vdatas(i)
                  do j=1,nt
                      if (rdatas(i).eq.rdatat(j)) then
                          vdatt(i)=vdatat(j)
                      else
                          call akima(np,nt,rdatat,vdatat,1,rdat(i),vt)
                          vdatt(i)=vt
                      endif
                  enddo
              enddo
C
              k=ns
              do i=1,nt
                  norepeat=.true.
                  do j=1,ns                
                      if (rdatat(i).eq.rdatas(j)) then
                          norepeat=.false.
                      endif
                  enddo
                  if (norepeat) then  
                      k=k+1 
                      rdat(k)=rdatat(i)
                      vdatt(k)=vdatat(i)
                      call akima(np,ns,rdatas,vdatas,1,rdat(k),vs)
                      vdats(k)=vs
                  endif
              enddo
              ndat=k
              call sort2 (rdat,vdats,vdatt,ndat)
C
              do i=1,ndat
                  write(52,1110) rdat(i),vdats(i),vdatt(i)
              enddo
C
          endif 
C 
C --- Fit Smirnov and Chibisov exchange form: A exp(-beta R) R^gamma 
C     to the exchange energy determined from the potential data to
C     obtain the amplitude A
C     
C     Calculate A at the chosen fitting point 
C
C      Interpolate potential data using akima spline interpolation
C      
          call akima(np,ns,rdatas,vdatas,1,rfit,vs)
          call akima(np,nt,rdatat,vdatat,1,rfit,vt)
C
C      Determine fitting constant A
C
          vt=vt/rfit**6
          vs=vs/rfit**6 
          vex=log((vt-vs)/2.0d0)
          CA=vex+beta*rfit-gamma*log(rfit)
          CA=exp(CA)
          write(*,*) CA
C
C --- Calculate the average of vs and vt, which should correspond 
C     with vdisp 
C
          vav=(vt+vs)/2.0d0
          vdisp=-c6/rfit**6-c8/rfit**8-c10/rfit**10
          vshift=vdisp-vav
          write(*,*) vshift
          write(*,*) vdisp
          write(*,*) exp(vex)
C
C --- Shift short range data by calculated difference to match the 
C     long range form
C
         ! vshift = 0d0
          nsl=0
          do i=1,ns
          !    vdatas(i)=vdatas(i)+vshift*rdatas(i)**6
              if(nsl.eq.0.and.rdatas(i).ge.rfit) nsl=i
          enddo
          ntl=0
          do i=1,nt
          !    vdatat(i)=vdatat(i)+vshift*rdatat(i)**6
              if(ntl.eq.0.and.rdatat(i).ge.rfit) ntl=i
          enddo
          if(nsl.eq.0.or.ntl.eq.0) then
              write(*,1020) rfit
              stop
          endif
C
C --- Connect short range data onto long range form using akima spline
C     to smooth across the connection 
C
!          call akima(np,ns,rdatas,vdatas,1,rfit,vs)
!          call akima(np,nt,rdatat,vdatat,1,rfit,vt)
          rdatas(nsl)=rfit      
          rdatat(ntl)=rfit
          vdatas(nsl)=(-(vt-vs)/2.0d0+vdisp)*rfit**6
          vdatat(ntl)=( (vt-vs)/2.0d0+vdisp)*rfit**6
          write(*,*) 'ntl =',ntl
          write(*,*) 'nt  =',nt          
!          vdatas(nsl)=vs
!          vdatat(ntl)=vt
C
          delta = 0.5d0
          ni=int((rfinal-rjoin)/delta)+1
!          ni = 2
          ndatas=nsl+ni
          ndatat=ntl+ni  
          ri=rjoin
          do i=1,ni             
              vex=CA*exp(-beta*ri)*ri**gamma
              vdisp=-c6/ri**6-c8/ri**8-c10/ri**10
C
              rdatas(nsl+i)=ri       
              rdatat(ntl+i)=ri
              vdatas(nsl+i)=(vdisp-vex)*ri**6       
              vdatat(ntl+i)=(vdisp+vex)*ri**6
C
              ri=ri+delta
          enddo
          rfinal=rdatas(ndatas)
C
C     Extrapolate short range data using the form: aexp(-br)
C  
          r1=rdatas(1)
          r2=rdatas(1)+1d-6
          call akima(np,ndatas,rdatas,vdatas,1,r1,v1)
          call akima(np,ndatas,rdatas,vdatas,1,r2,v2)
          if (v1*v2.lt.0.0d0) then
              write(*,1030) 
              stop
          endif
          v1=v1/r1**6
          v2=v2/r2**6 
          bs=log(v1/v2)/(r2-r1)
          as=v2*exp(bs*r2) 
C
          r1=rdatat(1)
          r2=rdatat(1)+1d-6
          !call akima(np,ndatat,rdatat,vdatat,1,r1,v1)
          call akima(np,ndatat,rdatat,vdatat,1,r2,v2)
          if (v1*v2.lt.0.0d0) then
              write(*,1030) 
              stop
          endif
          v1=vdatat(1)/r1**6
          v2=v2/r2**6
          bt=log(v1/v2)/(r2-r1)
          at=v2*exp(bt*r2) 
C
C     Output potential details/parameters
C
          write(*,1040) c6,c8,c10
          write(*,1050) rfit,CA,beta,gamma
          write(*,1060) vshift,vshift*hartree
          write(*,1070) rjoin,rfinal
          write(*,1080) as,bs,at,bt

      endif
C
      do_setup=.not.do_setup
      SYS_SS_pot_3=0.0d0
C
      RETURN 
      ENDIF    
C
C=======================================================================
C
      IF((NV.LE.0).OR.(NV.GT.2)) THEN
          WRITE(*,900)
          STOP
      END IF
C
      GO TO (10,20),NV
C
C=======================================================================
C
C --- SECTION FOR NV=1 -- Calculate singlet potential
C
   10 call mpot(np,ndatas,rdatas,vdatas,c6,c8,c10,
     >          CA,beta,gamma,as,bs,nv,r,v)
C
      v=v+correction_ss(r,shift1,re1)
      SYS_SS_pot_3=V
C
      RETURN
C
C=======================================================================
C
C --- SECTION FOR NV=2 -- Calculate triplet potential
C
   20 call mpot(np,ndatat,rdatat,vdatat,c6,c8,c10,
     >          CA,beta,gamma,at,bt,nv,r,v)
C
      v=v+correction_ss(r,shift3,re3)
      SYS_SS_pot_3=V
C
      RETURN
C
C=======================================================================
C
C --- FORMAT STATEMENTS ------------------------------------------------
C
  900 FORMAT(//,' Program HALTED in Potential Subroutine -- NV OUT ',
     X  'OF RANGE')
C
  990 FORMAT(//,'MOLECULAR POTENTIAL DATA')
 1000 format(//,'ERROR in Potential Subroutine -- no. of data points',
     X  'exceeds array size:',3i6)
 1010 format(//,'ERROR in Potential Subroutine -- identical data for',
     X  'pair:',2i6)
 1020 format(//,'ERROR in Potential Subroutine -- rfit (au)=',d14.6,
     X  'exceeds short range data')
 1030 format(//,'ERROR in Potential Subroutine -- data points used in',
     X  'short range extrapolation of data cross a zero')
 1040 format(//,'Long range dispersion coefficients (au):',/,
     X  ' C6=',d14.6,' C8=',d14.6,' C10=',d14.6)
 1050 format(//,'Long range exchange, Aexp(-beta R)R^gamma (au):',/,
     X  ' fit to short range data at R=',d14.6,/,
     X  ' A=',d14.6,' beta=',d14.6,' gamma=',d14.6)
 1060 format(/,'Shift applied to short range data to match to',
     X  ' long range form:',/,
     X  d14.6,' au',' or',d14.6,' cm-1')
 1070 format(//,'Long range form assumed at R (au)=',d14.6,/,
     X  'Short range data extended to R (au)=',d14.6)
 1080 format(//,'Short range extrapolation, aexp(-bR) (au)',/,
     X  'Singlet: a=',d14.6,' b=',d14.6,/,
     X  'Triplet: a=',d14.6,' b=',d14.6)
 1100 format(2d16.8)
 1110 format(3d16.8)
 1120 format(4d16.8)
C
      END
      !-------------------------------------------------------------------
      ! function  : vli2
      !
      ! package   : Li3 Potential
      !
      ! Language  : Fortran 77
      !
      ! author    : F. Colavecchia (flavioc@lanl.gov)
      !
      ! date      : 05/12/02
      ! revision  :
      !
      ! purpose   : wrapper for the sys_ss_pot_3 routine
      !
      ! input     :    r    -> interatomic distance in bohr.
      !
      ! output    :    vli2 -> li2 potential in a.u.  
      !
      !-------------------------------------------------------------------
      real*8 function vli2(r)
      implicit none
      real*8 r,v,sys_ss_pot_3
      integer nv
      logical lsetup
      data lsetup /.true./
      save lsetup           
      
      nv = 2   ! Triplet potential
      !nv = 1   ! Singlet potential
      v = 0d0 
      if(lsetup) then 
           vli2 = sys_ss_pot_3(r,v,nv,lsetup)
           lsetup = .false.
      end if
      vli2 = sys_ss_pot_3(r,v,nv,lsetup)

      if(r.eq.0.d0)vli2= 3.062396111034812d0

      return
      end  

