*dk potfh2
      subroutine potfh2(v,rij)
c
c   $RCSfile: potfh2.f,v $ $Revision: 1.4 $
c   $Date: 89/08/04 14:22:14 $
c   $State: Stable $
c
c   fh2 new surface no. 5 or 5a
c  quantum version.  does not calc all derivatives.
c   this version makes ab and ac sato parameters a function of phi.
c   phi is angle of bc-to-a vector with bc axis.
c   bc is a homonuclear diatom.  the bc sato parameter is made a
c   function of the r(bc) distance and the bend angle phi
c    a three center energy term is also added to correct the hfh
c    barrier height
c  parameters are set on first call.
c        coordinates, potential energy, and derivatives are passed
c           through common /fh2pot/ rab,rbc,rac,e,dab,dbc,dac
c
      implicit real*8 (a-h,o-z)
      dimension rij(3)
      dimension d(3),re(3),beta(3),     zpo(3),op3z(3),zp3(3),tzp3(3),
     1 top3z(3),do4z(3),b(3),e3(3),de1(3),de3(3)
      dimension x(3),coul(3),exch(3)
      common /coscom/ cosphi,ratio,x2,y2,z2,rc3,x1,y,z1
      common /phicom/ deg1,z0,fder,deg
      common /fh2pot/r(3),e,dedr(3)
      common /z3com/ z(3),lchi
      common/derivs/ds1dr1,ds1dr2,ds1dr3
      data r2/1.41421356d0/
c
c   potential energy surface parameters
c        energies in kcal/mol, lengths in angstoms
      data (d(i),i=1,3)/141.196d0,109.449d0,141.196d0/
      data (re(i),i=1,3)/0.9170d0,0.7419d0,0.9170d0/
      data (beta(i),i=1,3)/2.2187d0,1.9420d0,2.2187d0/
c  the initial values of the hf Sato parameters are 0.175 for Surface 5,
c  and 0.170 for Surface 5a
      data (z(i),i=1,3)/0.170d0,0.104d0,0.170d0/
      data icall/0/
c  the internuclear distances are:
c   rij(1) = r(HH),  rij(2) = r(FH),  rij(3) = r(FH)
c  shift to truhlar's notation for r
      r(2)=rij(1)
      r(3)=rij(2)
      r(1)=rij(3)
c  clip, if necessary
      v=100.d0
      if(min(r(1),r(2),r(3)).lt.0.01d0) return
      if(icall.ne.0) go to 101
      nc = 0
c   input ab,ac sato parameters apply only to small phi.
c      write (6,602) d,re,beta,z
  602 format (/36h potential energy surface parameters//13h sato-polanyi
     1 //5h bond,20x,2hab,8x,2hbc,8x,2hac//15h diss. energies,5x,
     2 3f10.5//12h equilibrium,8x,3f10.5//11h morse beta,9x,3f10.5//
     3 16h sato parameters,4x,3f10.5//)
c   the next two lines set data:
      lchi=0
c     lchi should be set ne.0 for debugging.  it causes lots of print.
c   note:  lchi controls printing in pot and fphi; set 0 for
c      production.
      deg1 = 10.0d0
c      write (6,612) deg1
  612 format (1h0,47hthe above sato parameters apply only for phi.lt  ,
     1   f10.3,1x,7hdegrees)
      if (z(1).ne.z(3)) stop 88
      z0=z(1)
      do  10 i = 1,3
c   convert to atomic units
      d(i)=d(i)/627.5095d0
      re(i) = re(i)/0.52917706d0
      beta(i) = beta(i)*0.52917706d0
c   compute useful constants
      zpo(i) = 1.0d0 + z(i)
      op3z(i) = 1.0d0 + 3.0d0*z(i)
      top3z(i) = 2.0d0*op3z(i)
      zp3(i) = z(i) + 3.0d0
      tzp3(i) = 2.0d0*zp3(i)
      do4z(i) = d(i)/4.0d0/zpo(i)
   10 b(i) = beta(i)*do4z(i)*2.0d0
      icall=1
  101 continue
      call phical (phi)
      if (lchi.ne.0) write (6,700) r(1),r(2),r(3),phi
  700 format (1h0,3f8.3,49x,f10.4)
      z(1)=fphi(phi)
      z(3) = z(1)
c
c     surface #5 differs in this next section from the previous
c     4 surfaces.  the hh sato parameter is made a function of
c     the hh distance and the bend angle phi.
c
      a = 1.26d0
      b1 = 1.60d0
      ta = tanh(a*(r(2) - b1))
      z(2) = .0395d0 + .201d0*(.5d0*(ta+1.0d0))
c
c     compute the change in z resulting from the bend
      cphi = cosphi
      sphi = sqrt(1.0d0 - cphi**2)
      dist = r(1)**2 + r(2)**2 - r(3)**2
      ctang = dist/(2.0d0*r(1)*r(2))
      arg2 = (r(1)**2 + (r(2)/2.0d0)**2 - r(1)*r(2)*ctang)
      if(arg2 .lt. 0.0d0) arg2 = 1.0d-28
      rfh2 = sqrt(arg2)
      rfh22 = rfh2**2
      g = 0.2d0
      g2 = g**2
      rat = rfh22/(g2+rfh22)
      dem = g2 + rfh22
      drat = 2.0d0*rfh2*(g2/dem**2)
c
      gscale = (z(2) - .120d0)*9.0033d0
      angter = -.101168d0*(sphi**2) - .46183d0*(sphi**4)
      z(2) = z(2) + angter*gscale*rat
c
      do 11 i=1,3
       zpo(i) = 1.0d0 + z(i)
       op3z(i) = 1.0d0 + 3.0d0*z(i)
       top3z(i) = 2.0d0*op3z(i)
       zp3(i) = z(i) + 3.0d0
       tzp3(i) = 2.0d0*zp3(i)
   11  do4z(i) = d(i)/4.0d0/zpo(i)
c
      s = 0.0d0
      do 21 i = 1,3
       x(i) = exp(-beta(i)*(r(i)-re(i)))
       coul(i) = do4z(i)*(zp3(i)*x(i)-top3z(i))*x(i)
       exch(i) = do4z(i)*(op3z(i)*x(i)-tzp3(i))*x(i)
   21 s = s+coul(i)
      rad = sqrt((exch(1)-exch(2))**2+(exch(2)-exch(3))**2+
     1      (exch(3)-exch(1))**2)
      e = -rad/r2
      e = e+s
      e = e + d(1)
c
c   this next section is the 3-centered term:
c
      con = .270063d0
      alp = 0.18d0
      aprime = .078d0
      par = .95d0
      atid = 2.14d0
      ones = exp(-alp*(r(1)-r(3))**2)
      fth = exp(-atid*(r(1)-r(3))**4)
      e3c = con*((1.0d0-par)*ones + par*fth)
c   this is a scaling factor to give the correct asym. limits
      dist = r(1)**2 + r(3)**2 - r(2)**2
      ctheta = dist/(2.0d0*r(1)*r(3))
      escale = (.48d0 + .52d0*ctheta**2)*exp(-aprime*(r(1)+r(3))**2)
      e3c = escale * e3c
      e = e3c + e
      v=e
      return
      end