      FUNCTION POTFUN  ( R, C )
      IMPLICIT REAL*8  ( A-H, O-Z )
      REAL*8             A( 6 )
      PARAMETER        ( ME = 6, ML = 6, MT = 2 * ME + 4 * ML + 4 )
      REAL*8             X( MT )
*
*     Potential energy function for the Ne - CO (rigid) complex obtained
*     from the fit to the ab initio results of S.M. Cybulski
*     r(C-O) was held constant at 1.128 A.
*     This routine produces the S2 (default) or S1 (by changing
*     comments in the DATA blocks below) Ne-CO potential surfaces
*     from the paper of McBane
*     and Cybulski, JCP  110 (24): 11734-11741 (1999).
*
*     On input:
*     =========
*     R - intermolecular distance (in angstroms)
*     C - cos( theta ), where theta is the angle between the bond axis
*         of a diatomic and a vector from Ar to the center of mass
*         of a diatomic.  theta=0 corresponds to the Ne-O-C linear geometry.
*
*     On return:
*     ==========
*     POTFUN - interaction energy in mE_h (milihartrees)


      DATA               X  /
c     'Ne-CO, CCSD(T) surface 2, 124 points'
     1  0.12131043D+02,-0.12488996D+00,-0.78833852D-01,-0.77755921D-02,
     2  0.38796877D+00,-0.73207269D-01,-0.31482161D+01,-0.89744451D-01,
     3  0.20725815D+00, 0.86820100D-04,-0.12426411D+00, 0.23044662D-01,
     4  0.40824590D+00,-0.55709087D+00, 0.58290514D+00, 0.25816728D+00,
     5  0.23434471D+00, 0.25747327D+00,-0.56581497D-01, 0.74644290D+00,
     6 -0.17880629D+00,-0.67439160D-01,-0.15587397D+00,-0.20270330D+00,
     7 -0.12372032D-01,-0.21445825D+00, 0.29742644D-01, 0.10092801D-01,
     8  0.43035936D-01, 0.54128863D-01, 0.45848082D-03, 0.17017634D-01,
     9 -0.46443096D-02,-0.18225895D-02,-0.42605078D-02,-0.47387062D-02,
     A -0.59708526D+03,-0.11431630D+03,-0.12066476D+04,-0.25836209D+03/
c  'Ne-CO, CCSD(T) surface 1, 124 points'
c    1   0.12287790D+02  0.49775350D+00  0.30776674D+00 -0.18777375D-01
c    2   0.16814082D+00 -0.39170014D-02 -0.31289714D+01 -0.27673456D+00
c    3   0.54625293D-01  0.49449239D-02 -0.48797266D-01 -0.44483726D-04
c    4   0.32521871D+00 -0.85981745D+00  0.24178347D+00 -0.76120427D-01
c    5   0.45580148D-01  0.80746309D-01 -0.76652404D-01  0.82705136D+00
c    6  -0.64946618D-01  0.11538292D+00 -0.16199725D-01 -0.59986433D-01
c    7   0.75035729D-02 -0.20301692D+00  0.34967388D-01 -0.21004145D-01
c    8   0.59992043D-02  0.15809871D-01 -0.18646740D-02  0.13947276D-01
c    9  -0.75765071D-02  0.41662687D-04 -0.66898802D-03 -0.13462145D-02
c    A  -0.58323782D+03 -0.81398215D+02 -0.13254086D+04 -0.43232382D+03/

*
      PARAMETER        ( ONE = 1.0D0, THREE = 3.0D0, F15 = 15.0D0,
     *                   ZERO = 0.0D0, FIVE = 5.0D0, THIRTY = 30.0D0,
     *                   SEVNTY = 70.0D0, HALF = 0.5D0, T35 = 35.0D0,
     *                   S63 = 63.0D0, THRESH = 1.0D-8 )
*
*     Start by evaluating the required Legendre polynomials
*
      T = -C
      T2 = T * T
      T3 = T2 * T
      T4 = T3 * T
      T5 = T4 * T
      T6 = T5 * T
      A( 1 ) = ONE
      A( 2 ) = T
      A( 3 ) = HALF * ( THREE * T2 - ONE )
      A( 4 ) = HALF * ( FIVE * T3 - THREE * T )
      A( 5 ) = 0.125D0 * ( T35 * T4 - THIRTY * T2 + THREE )
      A( 6 ) = 0.125D0 * ( S63 * T5 - SEVNTY * T3 + F15 * T )
*
*     VSH( R, THETA ) part:
*
      D = ZERO
      B = ZERO
      DO 100 I = 1, ME
         D = D + X( I ) * A( I )
         B = B + X( I + ME ) * A( I )
 100  CONTINUE
      NE = 2 * ME
      R2 = R * R
      R3 = R2 * R
      G = ZERO
      DO 200 I = 1, ML
         IE = I + NE
         G = G + ( X( IE ) + X( IE + ML ) * R + X( IE + 2 * ML )
     *         * R2 + X( IE + 3 * ML ) * R3 ) * A( I )
 200  CONTINUE
      VSH = G * EXP( D + B * R )
*
*     VAS( R, THETA ) part:
*
      BR = B * R
      S = ONE
      T = ONE
      DBR = DABS( BR )
      DO 300 K = 1, 6
         T = T * DBR / DFLOAT( K )
         S = S + T
 300  CONTINUE
      F6 = ONE - DEXP( BR ) * S
      T = T * DBR / DFLOAT( 7 )
      S = S + T
      F7 = ONE - DEXP( BR ) * S
*
*     If F6 or f7 are smaller than 1.0D-8, they should be recalculated
*     and the segment below will do the trick.
*
      IF( DABS( F6 ) .LT. THRESH ) THEN
         F6 = ZERO
         DO 400 I = 7, 1000
            T = T * BR / DFLOAT( I )
            F6 = F6 + T
            IF( ( T / F6 ) .LT. THRESH ) GO TO 500
 400     CONTINUE
         WRITE( 6, * ) 'No convergence for F6.'
         STOP
 500     CONTINUE
         F6 = F6 * DEXP( BR )
      END IF
*
      IF( DABS( F7 ) .LT. THRESH ) THEN
         F7 = ZERO
         DO 600 I = 8, 1000
            T = T * BR / DFLOAT( I )
            F7 = F7 + T
            IF( ( T / F7 ) .LT. THRESH ) GO TO 700
 600     CONTINUE
         WRITE( 6, * ) 'NO CONVERGENCE IN D'
 700     CONTINUE
         F7 = F7 * DEXP( BR )
      END IF
*
      R6 = R3 * R3
      R7 = R6 * R
      NE = NE + 4 * ML
      VAS = F6 * ( X( 1 + NE ) + X( 2 + NE ) * A( 3 ) ) / R6
     *    + F7 * ( X( 3 + NE ) * A( 2 ) + X( 4 + NE ) * A( 4 ) ) / R7
      POTFUN = VSH + VAS
      RETURN
      END



