C***********************************************************************
C                                 LICNSRF
C                           
C Computes an analytic fit to the two-dimensional SCF potential energy
C surface of LiCN (the CN bond lenght is frozen at 2.186 a.u.).
C Formulas are taken from:
C     R. Essers, J. Tennyson and P.E.S. Wormer: Chem. Phys. Lett. 89,
C     223 (1982).
C Note that the parameters in that paper are not reported with the
C precision required to exactly reconstruct the surface. The parameters
C in the following routine, however, are the exact ones, taken from:
C     J. Tennyson: private communications.
C
C***********************************************************************


C***********************************************************************
C                                 LICNSRF 
C                           
C Computes an analytic fit to the two-dimensional SCF potential energy
C surface of LiCN. The potential fit is a Legendre-polynomial expansion:
C     V(r,theta) = Sum[V_k(r) P_k(cos(theta))], k=0,...,9.
C This routine computes V_k(r) only.
C
C Input parameters:
C    k: expansion term as defined above
C    r: distance (in a.u.) of the Li nucleus from the center of nuclear
C       mass of the CN diatom
C Output parameters:
C    error: error code having the following meaning:
C           0: everything was o.k.
C           1: illegal value for "r"
C Function value:
C    the value of V_k(r) (in a.u.) as defined above
C***********************************************************************

      real*8 function licnsrf (k,r,error)

      implicit none

      real*8  r,rmax,a0,r0,q(0:6),a(0:9),b(0:9),c(0:9),d(4:8,0:6)
      integer k,error,j

      data a0,r0,rmax / 1.515625d0, 1.900781d0, 9.0d0 /
      data q
     +   / -1.000000d0, -0.215135d0, -3.414573d0, -3.818815d0,
     +    -15.841520d0,-14.293740d0,-43.817190d0 /
      data a
     +   / 1.3832116000d0,  2.9579132000d0,  4.7420297000d0,
     +     1.8885299000d0,  4.4143354000d0,  4.0256496000d0,
     +     5.8425899000d0,  2.6168114000d0,  6.3446579000d0,
     +   -15.2022800000d0 /
      data b
     +   /-0.1400070600d0, -1.4797716000d0, -1.8119862000d0,
     +    -1.2875030000d0, -2.3229714000d0, -2.7753832000d0,
     +    -3.4808529000d0, -2.6555952000d0, -4.3449801000d0,
     +     6.5492537000d0 /
      data c
     +   /-0.2078921600d0,  0.0116130820d0,  0.0171806800d0,
     +    -0.0277284910d0,  0.0706927420d0,  0.1377197800d0,
     +     0.1863111400d0,  0.0058815504d0,  0.1529136800d0,
     +    -1.3025678000d0 /
      data (d(4,j),j=0,6)
     +   / -10.52710d0,   0.00000d0,  -3.17000d0,   0.00000d0,
     +       0.00000d0,   0.00000d0,   0.00000d0 /
      data (d(5,j),j=0,6)
     +   /   0.00000d0, -20.62328d0,   0.00000d0,   3.73208d0,
     +       0.00000d0,   0.00000d0,   0.00000d0 /
      data (d(6,j),j=0,6)
     +   / -57.49396d0,   0.00000d0,-106.81920d0,   0.00000d0,
     +      17.14139d0,   0.00000d0,   0.00000d0 /
      data (d(7,j),j=0,6)
     +   /   0.00000d0,-202.89720d0,   0.00000d0, -75.23207d0,
     +       0.00000d0, -28.45514d0,   0.00000d0 /
      data (d(8,j),j=0,6)
     +   /-458.20150d0,   0.00000d0,-353.73470d0,   0.00000d0,
     +    -112.64270d0,   0.00000d0,-108.27860d0 /

C --- INITIALISE VARIABLES ---

      licnsrf = 0.0d0

C --- CHECK WHETHER PARAMETERS LIE IN ALLOWED RANGE ---

      error = 0
      if (r .lt. r0) then
         error = 1
         return
      endif

C --- COMPUTE ELECTROSTATIC ENERGY ---

      if (k .le. 6) then
         licnsrf = licnsrf+q(k)/r**(k+1)
      endif

C --- ADD INDUCTION ENERGY ---

      if (k .le. 6) then
         do j = 4,8
            licnsrf = licnsrf+d(j,k)/r**j
         enddo
      endif

C --- MULTIPLY ENERGY WITH GAUSSIAN DAMPING FUNCTION ---

      licnsrf = licnsrf*(1.0d0-exp(-a0*(r-r0)**2))

C --- ADD SHORT-RANGE ENERGY ---

      if ((r .le. rmax) .and. (k .le. 9)) then
         licnsrf = licnsrf+exp(a(k)+(b(k)+c(k)*r)*r)
      endif
      
C --- RETURN TO CALLING ROUTINE ---

      return
      end
