SUBROUTINE h3pk (v, r1)
USE Numeric_Kind_Module
!
!   $RCSfile: h3pk.f,v $ $Revision: 1.3 $
!   $Date: 89/08/07 09:13:26 $
!   $State: Stable $
!
!-----------------------------------------------------------------------
! This is the Porter-Karplus HHH potential energy surface.
!  r1 is in bohr. v is in Hartree atomic units.
!-----------------------------------------------------------------------
INTEGER i, k
REAL(KIND=DP_Kind) a, alf, bet, c1, c2, c3, clip, d1, d3, delta, e1, e3, epsi, fj123
REAL(KIND=DP_Kind) fkap, flam, h, q, q1, r, rh, rz, s, u, v, z
REAL(KIND=DP_Kind) r1 (3), f1 (3), r2 (3), s2 (3), fj (3)
DATA u / 1.0d0 /, h / 0.5d0 /
DATA d1 / 4.746600d0 /, d3 / 1.966800d0 /, alf / 1.0443500d0 /, &
 bet / 1.00012200d0 /, delta / 28.200d0 /, epsi / - 17.500d0 /, &
 fkap / 0.6000d0 /, flam / 0.6500d0 /, rz / 1.4008300d0 /, rh / &
 1.d0 /
DATA clip / 8.7571556D-01 /
IF(dmin1 (r1 (1), r1 (2), r1 (3) ) .gt.0.25d0)THEN
   fj123 = epsi
   q = 0.0d0
   DO i = 1, 3
      IF(r1 (i) .eq.0.d0) RETURN
      r2 (i) = r1 (i) / rh
      f1 (i) = delta * (u + u / r2 (i) ) * exp ( - 2.0d0 * r2 (i) &
       )
   ENDDO
   s = u
   DO k = 1, 3
      r = r2 (k)
      z = (u + fkap * exp ( - flam * r) ) * r
      z = (u + z + z * z / 3d0) * exp ( - z)
      s2 (k) = z * z
      fj123 = fj123 * z
      s = s * z
      a = exp ( - bet * (r - rz) )
      e3 = d3 * (a * a + 2d0 * a)
      a = exp ( - alf * (r - rz) )
      e1 = d1 * (a * a - 2d0 * a)
      IF(K.eq.1)THEN
         a = f1 (2) + f1 (3)
      ELSEIF(K.eq.2)THEN
         a = f1 (1) + f1 (3)
      ELSEIF(K.eq.3)THEN
         a = f1 (1) + f1 (2)
      ELSE
         STOP "h3pk"
      ENDIF
      fj (k) = h * (e1 - e3) + s2 (k) * (h * (e1 + e3) + a)
      q1 = h * (e1 + e3 + s2 (k) * (e1 - e3) )
      q = q + q1
   ENDDO
   c3 = (q - fj123) **2 - h * ( (fj (1) - fj (2) ) **2 + (fj (2) &
    - fj (3) ) **2 + (fj (3) - fj (1) ) **2)
   c2 = - (q - fj123) * (u - s) + h * ( (fj (1) - fj (2) ) &
    * (s2 (1) - s2 (2) ) + (fj (2) - fj (3) ) * (s2 (2) - s2 (3) ) &
    + (fj (1) - fj (3) ) * (s2 (1) - s2 (3) ) )
   c1 = (u - s) **2 - h * ( (s2 (1) - s2 (2) ) **2 + (s2 (1) &
    - s2 (3) ) **2 + (s2 (2) - s2 (3) ) **2)
   v = - c2 / c1
   IF(c2**2.gt. (c1 * c3) ) v = v - sqrt (c2**2 - c1 * c3) &
    / c1
   v = v + 4.7466d0
   v = v / 27.2113961d0
!
ELSE
   v = clip*exp(5.0*(0.25d0-dmin1 (r1 (1), r1 (2), r1 (3) )))
ENDIF
RETURN
END SUBROUTINE h3pk