*deck vh2(kw)
      function vh2 ( r )
c
c     call with r in a.u.(bohr)
c     returns h2 potential energy in au if units=1.
c
      implicit real*8 (a-h,o-z)
      common / potpar / vbase, emass, units
      common / poth2 / c6,c8, rkw(87), ekw(87), wkw(87)
      dimension tab(3)
      data ncall / 0 /
      data conv1, conv2, conv3 / .52917706d-8, 27.2106d0,
     >     1.6021892d-12 /
      if ( ncall .gt. 0 ) go to 1
      call h2data
      do 100 i = 1, 87
      rkw(i) = rkw(i)/conv1
      ekw(i) = ekw(i)/(conv2*conv3)
      wkw(i) = wkw(i)/(conv2*conv3)*conv1*conv1
  100 continue
      c6 = c6/(conv2*conv3*conv1**6)
      c8 = c8/(conv2*conv3*conv1**8)
    1 continue
      iflag = 1
   10 continue
      ncall = ncall + 1
      if(r.gt.10.d0)go to 200
      if(r.lt.0.24262d0) go to 300
      call spl1d2(87,rkw,ekw,wkw,1,r,tab)
      vh2=tab(iflag)*units
c      if ( iflag .eq. 1 ) vh2 = vh2 - vbase
      return
  200 continue
      c6r6 = c6/r**6
      c8r8 = c8/r**8
      tab(1) = -units*(c6r6 + c8r8)
      tab(2) = units*(6.d0*c6r6 + 8.d0*c8r8)/r
      tab(3) = -units*(42.d0*c6r6 + 72.d0*c8r8)/r**2
      vh2 = tab(iflag)
      return
c
  300 continue
      if(r.lt.1.d-40)then
      vh2=1.d+40
      return
      endif
      vh2=units*(-1.90372d0 + 1.d0/r)
      return
c
      entry dvh2
      iflag = 2
      go to 10
c
      entry ddvh2
      iflag = 3
      go to 10
c
      end
