      DOUBLE PRECISION FUNCTION FLI2T(R)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

C*********************************************************************C
C                                                                     C
C                     LI2  TRIPLET SIGMA U +                          C
C                                                                     C
C  TWICE THE REDUCED MASS (TWOMU)                                     C
C  FUNCTION IS CALLED ONCE WITH R=0.D0 TO GET SPLINE COEFFS BEFORE    C
C   MAIN RUN.  ONE R PER CALL.                                        C
C                                                                     C
C  POTENTIAL IS FROM KONAWALOW ET AL JCP 81 P.4534 1984               C
C  (see ZEMKE2.F for details)                                         C
C*********************************************************************C

      parameter (nmax=1000)
      PARAMETER (TWOMU=1.D0)
c     PARAMETER (TWOMU=6.015121D0*1822.88853D0)
c     PARAMETER (TWOMU=7.016003D0*1822.88853D0)
c     PARAMETER (TWOMU=2.d0*7.016003D0*6.015121D0*1822.88853D0/
c    &           (7.016003D0+6.015121D0))
      PARAMETER (HARTREE=219474.64D0, BOHR=0.529177249D0)
c Li6
c      PARAMETER (DET=333.74585D0,NCN=166)
c     parameter (det=333.74d0,ncn=166) ! approximate
c Li7 
      PARAMETER (DET=333.78857D0,NCN=166)
c     parameter (det=333.78d0,ncn=166) ! approximate
c Li6-Li7
c     parameter (det=333.76d0,ncn=166)

      parameter (csix=1393.39d0,ceight=83425.8d0,cten=7372100.d0)  
      PARAMETER (AMP=0.01288D0,ALPHA=4.558D0,BETA=1.259D0)
      DIMENSION RAT(8),PRT(8), rzs(157), ezs(157)
      dimension rk(7), pk1(7), pk2(7), pk3(7), pk4(7)
      DIMENSION RRT(nmax),PT(nmax),Y2(nmax),UU(nmax)

      OPEN(UNIT=33,FILE="zemke.pot",STATUS="UNKNOWN")

      IF (R.EQ.0.D0) THEN
 
C
C SECTION WITH E>0 : SMALL R
C FROM KONOWALOW (1ST PT) AND SCHMIDT-MINK (7 OTHER PTS)
C THESE ARE IN A.U.
C

      ddet = det 

C
      DATA PRT  / 0.104471D0,   0.082074D0,   0.063725D0,
     $           0.036778D0,   0.020006D0,   0.010053D0,
     $           0.004370D0,   0.001258D0 /

      DATA RAT /  3.00D0, 3.25D0, 3.50D0,
     $           4.00D0, 4.50D0, 5.00D0,
     $           5.50D0, 6.00D0 /


      data rk /  3.00D0, 3.50D0,
     $           4.00D0, 4.50D0, 5.00D0,
     $           5.50D0, 6.00D0 /

c data OK-OVC (Konowalow, Regan and Rosenkrantz)
 
      data pk1 / 0.105985d0,  0.065544d0,
     $           0.038288d0,  0.021133d0, 0.010850d0,
     $           0.004920d0,  0.001637d0 /

c data Corrected OVC (Konowalow, Regan and Rosenkrantz)

      data pk2 / 0.105607d0,  0.065326d0,
     $           0.038023d0,  0.020839d0, 0.010547d0,
     $           0.004625d0,  0.001363d0 /

c data KF (Konowalow, Regan and Rosenkrantz)

      data pk3 / 0.103335d0,  0.062060d0, 
     $           0.035123d0,  0.018623d0, 0.008962d0,
     $           0.003539d0,  0.000648d0 /

c data Composite (Konowalow, Regan and Rosenkrantz)

      data pk4 / 0.104471d0,  0.063693d0,
     $           0.036573d0,  0.019731d0, 0.009755d0,
     $           0.004082d0,  0.001005d0 /

C
C SECTION WITH E<0 : DRKR POTENTIAL
C FROM Zemke & Stwalley 
C THESE ARE IN SPECTSTROSCOPIC UNITS : CM-1 AND ANGSTROMS
C

      do 33 i = 1,157 
         read(33,*) rzs(i), ezs(i)
   33 continue

C
C CONSTRUCTION OF THE POTENTIAL
C

      do 202 i=1,7
         rrt(i) = rk(i)
         pt(i) = pk1(i)
  202 continue
      DO 210 I = 1, 157 
         RRT(I+7) = rzs(I)/BOHR
         PT(I+7) = (ezs(I)-DDET)/HARTREE
  210 CONTINUE
       
      rrt(165)=rrt(164)+(rrt(164)-rrt(163))
      rrt(166)=rrt(164)+2*(rrt(164)-rrt(163))
      z = rrt(165)
      C = CSIX/Z**6 + CEIGHT/Z**8 + CTEN/Z**10
      TAIL = AMP*(Z**ALPHA)*DEXP(-BETA*Z)
      pt(165) = -( C - TAIL )
      z = rrt(166)
      C = CSIX/Z**6 + CEIGHT/Z**8 + CTEN/Z**10
      TAIL = AMP*(Z**ALPHA)*DEXP(-BETA*Z)
      pt(166) = -( C - TAIL )

      U = RRT(NCN)
      SLOPE = (6*CSIX/(U**7)) + (8*CEIGHT/(U**9)) + (10*CTEN/(U**11))
     &        + (ALPHA/U - BETA)*AMP*(U**ALPHA)*DEXP(-BETA*U)
      YP1=1.D40
      YPN=SLOPE
      CALL SPLINE(RRT,PT,NCN,YP1,YPN,Y2,UU)
      u = RRT(1)
      CALL SPLINT(RRT,PT,Y2,NCN,u,S,1)
      BS = -S/PT(1)
      AS = PT(1)*DEXP(BS*RRT(1))

C
C THE NEXT THREE BRANCHES GIVE THE POTENTIAL DEPENDING OF R
C
      ELSE IF ((R.GT.0.D0).AND.(R.LT.RRT(1))) THEN
         FLI2T = TWOMU*AS*DEXP(-BS*R)
      ELSE IF ((R.GE.RRT(1)).AND.(R.LE.RRT(NCN))) THEN
         CALL SPLINT(RRT,PT,Y2,NCN,R,POTS,0)
         FLI2T = TWOMU*POTS
      ELSE
         C = CSIX/R**6 + CEIGHT/R**8 + CTEN/R**10
         TAIL = AMP*(R**ALPHA)*DEXP(-BETA*R)
         FLI2T = -TWOMU*( C - TAIL )
      END IF

      RETURN
      END
