      function potntl(ar)
c
      implicit real*8 ( a - h, o - z )
c
c      ref T. Joseph and N. Satyamurthy, jcp 86,704,(1987) (ci-sm)
c     r(1) = r(h-h) = rbc
c     r(2) = r(he-h) = en-to-end = rac
c     r(3) = r(he-h) = rab
c
      real*8 ar(3),rr(3),vv(3),a2(3),a3(3),r0(3),
     .  c(5,5,5), dd(3), a1(3), re(3)
      common /leps/de(3),beta(3),req(3),dummy(3),vbase
c
cc    external hehhdat
c
      data gm1,gm2 /1.3749d0,2.2225d0/
      data a2,a3,r0 /2.9402d0,1.5746d0,2.9402d0,
     .              2.6185d0,1.2451d0,2.6185d0,
     .              0.9525d0,0.9525d0,0.9525d0/
      data c       /6.3465d0, -3.6338d0, -2.9693d0, 4.7870d0, -0.4599d0,
     .              9.3453d0, -0.0874d0, -9.7453d0, 4.2826d0, 0.0d0,
     .              5.1718d0, 5.8826d0, -4.0011d0, 0.0d0, 0.0d0,
     .              0.6141d0, 6.4306d0, 0.0d0, 0.0d0, 0.0d0,
     .              -1.9547d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
c
     .              -3.6338d0, 6.9382d0, -3.5791d0, -3.8079d0, 0.0d0,
     .              -0.0874d0, 4.1563d0, 0.2004d0, 0.0d0, 0.0d0,
     .              5.8826d0, -0.1249d0, 0.0d0, 0.0d0, 0.0d0,
     .              6.4306d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
     .              0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
c
     .              -2.9693d0, -3.5791d0, 7.9695d0, 0.0d0, 0.0d0,
     .              -9.7453d0, 0.2004d0, 0.0d0, 0.0d0, 0.0d0,
     .              -4.0011d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
     .              10*0.0d0,
c
     .              4.7870d0, -3.8079d0, 0.0d0, 0.0d0, 0.0d0,
     .              4.2826d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
     .              15*0.0d0,
c
     .              -0.4599d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
     .              20*0.0d0/
c
      data ifirst / 2 /
      if (ifirst.ne.1) go to 5
      print 105,gm1,gm2
  105 format (' gm1 =' f10.5,5x,'gm2 =' f10.5)
      print 110,d
  110 format (' d array'/3f10.5)
      print 120,a1
  120 format (' a1 array'/3f10.5)
      print 130,re
  130 format (' re array'/3f10.5)
      print 140,r0
  140 format (' r0 array'/3f10.5)
      print 150,a2
  150 format (' a2 array'/3f10.5)
      print 160,a3
  160 format (' a3 array'/3f10.5)
      print 170,c
  170 format (' c array'/(5f10.5))
      ifirst = 2
    5 continue
      if (ifirst.ne.2) go to 8
c       energy in eV and distance in Angstrom
      dd(2) = 2.7865d0
      dd(1) = 2.0279d0
      dd(3) = 2.0279d0
      a1(2) = 2.6180d0
      a1(1) = 4.6093d0
      a1(3) = 4.6093d0
      re(2) = 1.0536d0
      re(1) = 0.7760d0
      re(3) = 0.7760d0
      ifirst = 3
    8 continue
c      switch to the papers arrangement of r
c       (and convert from bohr to Angstrom)
c
      rr(1) = ar(2)*0.529177249d0
      rr(2) = ar(1)*0.529177249d0
      rr(3) = ar(3)*0.529177249d0
c
      do 10 i=1,3
      rho=rr(i)-re(i)
      texp=exp(-a1(i)*rho)
      vv(i)=-dd(i)*(1.d0+rho*(a1(i)+rho*(a2(i)+rho*a3(i))))*texp
   10 continue
      sum=0.d0
      rho1=rr(1)-r0(1)
      rho2=rr(2)-r0(2)
      rho3=rr(3)-r0(3)
      do 40 i=1,5
      rho1i=rho1**(i-1)
      do 30 j=1,5
      rho2j=rho2**(j-1)
      do 20 k=1,5
      if (i+j+k.gt.7) go to 20
      rho3k=rho3**(k-1)
      sum=sum+c(i,j,k)*rho1i*rho2j*rho3k
   20 continue
   30 continue
   40 continue
      atanh=tanh(gm1*0.5d0*(rho1+rho3))
      btanh=tanh(gm2*rho2*0.5d0)
      ta=(1.d0-atanh)*(1.d0-btanh)
      potntl=vv(1)+vv(2)+vv(3)+ta*sum + 2.7865d0
cc    if (crlmp) potntl = potntl + brlm/(rs**2 + ds**2)
      potntl=potntl/27.2113961d0
      return
      end
