function vneh(r)
c  combining rule potential for NeH.  See K. T. Tang and J. P. Toennies,
c  Chem. Phys. 156, 413 (1991)
c  atomic units for both vneh and r.
      implicit none
      character*8 procname
      parameter (procname='vneh    ')
      integer n, twon, icall
      real*8 vneh, 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
c  ------1--------2---------3---------4---------5---------6---------7--
c  determine parameters on the first call.  All are in au.     
      if (icall.eq.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
      end if
c
c  if r is too large just set v=0.
      vneh = 0.d0
      if(r.gt.1.d03) return
c  if r is too small, tang-toennies form is unstable.  Use following.
      if(r.eq.0.0d0)then
         vneh = a2
         return
      end if
      if(r.lt.0.5d0) then
         vneh = a2*exp(-b2*r)
         return
      end if
c  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
c  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
c  add scf term to get final potential
      vneh = vneh + a*ex
c
      return
      end