subroutine gpif(iai,iaj,rij,epif,depif,lok)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     File     : pif.F
c     Author   : Gerald Monard
c     Date     : 01 fev 2000
c     Change   : divcon.1.0.C011
c     Function : compute Parametrizable Interaction Function for PM3
c                as described by Bernal-Uruchurtu et al.
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
#include "constants.h"
      real*8 epif,depif         ! result : pif function
      integer iai, iaj          ! atomic number of atom i and j
      real*8 rij,r              ! distance between i and j
      real*8 eorig              ! original interaction energy

      real*8 alpha, beta, gamma, delta, epsilon ! PIF parameters
      integer imin, imax        ! minimum and maximum atomic number
                                ! (they are sorted)
      real*8 r2, r6, r8, r10    ! local variables
      real*8 dr2, dr6, dr8, dr10! local variables
      real*8 a0conv
      logical lok               ! logical return: true if PIF computation was
                                ! succesful, false unless

      ! three types of interaction
      ! H-H, O-H and O-O

      epif = 0.0d0
      depif = 0.0d0
      imin = min(iai,iaj)
      imax = max(iai,iaj)
      a0conv = 1.0d0/Bohr2Ang
      r = rij*a0conv
      if (imin.eq.1) then
         if (imax.eq.1) then
            ! H-H interaction
            alpha   = .5197876E+00 
            beta    =  .2489980E+01
            gamma   =  .5010107E+02 
            delta   =  -.6763223E+03 
            epsilon = .2711768E+04 
         else if (imax.eq.8) then
            ! H-O interaction
            alpha   = .4936554E+02 
            beta    = .2276467E+01
            gamma   = -.3828750E+02
            delta   = -.1450515E+03
            epsilon =  .6518225E+03 
         else if (imax.eq.17) then
            ! H-Cl interaction
            alpha   =  27.42530
            beta    =   2.521105
            gamma   = -19.93522
            delta   =  29.06138
            epsilon = 280.3394
         else ! error
            goto 100
         endif
      else if (imin.eq.8) then
         if (imax.eq.8) then
            ! O-O interaction
            alpha   = .1245958E+02 
            beta    = .2502734E+01 
            gamma   =-.2662244E+03
            delta   = .1258753E+05
            epsilon =-.8948691E+05
         else if (imax.eq.17) then
            ! O-Cl interaction
            alpha   =    35.27732
            beta    =     1.292906
            gamma   =  -680.9553
            delta   =    19.13535
            epsilon = 57189.97
         else ! error
            goto 100
         endif
      else ! error
         goto 100
      endif
      r2 = 1.0d0/(r*r)
      dr2 = -2.0d0/(r**3)
      r6 = r2*r2*r2
      dr6 = -6.0d0/(r**7)
      r8 = r6*r2
      dr8 = -8.0d0/(r**9)
      r10 = r8*r2
      dr10 = -10.0d0/(r**11)
      epif = alpha*dexp(-beta*r) + gamma*r6 + delta*r8 + epsilon*r10
      depif = -beta*alpha*dexp(-beta*r) + gamma*dr6 
     %         + delta*dr8 + epsilon*dr10
      epif = epif * hartree2eV
      depif = depif * hartree2eV
c     write(*,*) 'pif',iai,iaj,rij,pif
      lok = .true.
      return

100   continue
c     write(*,'("WARNING WARNING !!!")')
c     write(*,'(" In pif.F : atomic number not recognized")')
c     write(*,'(" imin = ",i3)') imin
c     write(*,'(" imax = ",i3)') imax
      lok = .false.
      return
      end