REAL(Kind=WP_Kind) FUNCTION VNeH(r)
USE Numeric_Kinds_Module
!  combining rule potential for NeH.  See K. T. Tang and J. P. Toennies,
!  Chem. Phys. 156, 413 (1991)
!  atomic units for both vneh and r.
IMPLICIT NONE
CHARACTER(LEN=21), PARAMETER:: ProcName='vneh    '
INTEGER n, twon, icall
REAL(Kind=WP_Kind) r, a, b, c, f, rs, sum, br, ex, termo, term, one, a2, b2
DIMENSION c(18), f(18), rs(9), sum(0:18)
data icall /0/
SAVE
!  ------1--------2---------3---------4---------5---------6---------7--
!  determine parameters on the first CALL.  All are in au.     
IF(icall==0)THEN
   a=11.78d0
   b=1.873d0
   c(6) = 5.71d0
   c(8) = 92.59d0
   c(10) = 2007.0d0
   one = 1.d0
   a2 = 11.71258660d0
   b2 = 1.887834483d0
   DO 1 n=6,8
      twon = 2*n
      c(twon) = (c(twon-2)/c(twon-4))**3 * c(twon-6)
1        continue
   icall = 1
ENDIF
!
!  if r is too large just set v=0.
vneh = 0.d0
IF(r>1.d03) RETURN
!  if r is too small, tang-toennies form is unstable.  Use following.
IF(r==0.0d0)THEN
   vneh = a2
   RETURN
ENDIF
IF(r<0.5d0)THEN
   vneh = a2*exp(-b2*r)
   RETURN
ENDIF
!  begin calculation at the given r.
br = b*r
ex = exp(-br)         
rs(1) = r*r
DO 3 n=2,8
   rs(n) = rs(1)*rs(n-1)
3     continue
sum(0) = one
termo = one
DO 4 n=1,16
   term = br*termo/n
   sum(n) = sum(n-1) + term
   f(n) = one - ex*sum(n)
   termo = term
4     continue
!  construct attractive, damped vdw potential
vneh = 0.d0
DO 5 n=3,8
   twon = 2*n
   vneh = vneh -f(twon)*c(twon)/rs(n)
5     continue
!  add scf term to get final potential
vneh = vneh + a*ex
!
RETURN
END FUNCTION VNeH