subroutine gmais(iai,iaj,rij,epif2,depif2)
      implicit none
!--------------------------------------------------------------------
!     Function: computes the potential interaction function for water
!               as described in Chem. Phys. Lett. 330 (2000) p118-124
!     Date: December 15, 2000
!     Author: Gerald Monard
!--------------------------------------------------------------------
#include "constants.h"
      integer iai                  ! input: the atomic number of i
      integer iaj                  ! input: the atomic number of j
      double precision rij         ! input: the distance between atoms i and j
      double precision epif2      ! output: the interaction energy
      double precision depif2

      ! local variable
      double precision a0conv
      double precision r
      double precision alpha(3), beta(3), gamma(3)
      integer i
      integer imin
      integer imax

      ! description of the function
      epif2 = 0.0d0
      depif2 = 0.0d0
      imin = min(iai,iaj)
      imax = max(iai,iaj)
      a0conv = 1.0d0/Bohr2Ang
      r = rij*a0conv

      if (imin.eq.1.and.imax.eq.1) then
         ! H-H interaction
         alpha(1) =  0.000362d0
         alpha(2) =  0.009138d0
         alpha(3) =  0.007175d0
         beta(1)  =  0.385490d0
         beta(2)  =  0.227530d0
         beta(3)  =  3.013020d0
         gamma(1) =  6.23090d0
         gamma(2) =  2.21168d0
         gamma(3) =  2.22572d0
      else if (imin.eq.1.and.imax.eq.8) then
         ! O-H interaction
         alpha(1) = -0.071576d0
         alpha(2) =  0.015255d0
         alpha(3) =  0.029575d0
         beta(1)  =  0.192055d0
         beta(2)  =  1.284460d0
         beta(3)  =  0.345024d0
         gamma(1) =  1.44174d0
         gamma(2) =  2.19381d0
         gamma(3) =  2.58602d0
      else if (imin.eq.8.and.imax.eq.8) then
         ! O-O interaction
         alpha(1) = -6.797745d0
         alpha(2) =  6.912401d0
         alpha(3) =  0.074600d0
         beta(1)  =  0.326159d0
         beta(2)  =  0.320157d0
         beta(3)  =  1.268150d0
         gamma(1) =  1.21137d0
         gamma(2) =  1.20945d0
         gamma(3) =  2.57953d0
      else 
         write(*,'("error: bad use of pif2 function")')
         write(*,'("error: pif2 call with atomic number",i5," and",i5)')
     .         imin, imax
         stop
      endif 

      do i = 1, 3
         epif2 = epif2 + alpha(i)*dexp(-beta(i)*(gamma(i)-r)**2)
         depif2 = depif2 + 2.0d0*alpha(i)*beta(i)*(gamma(i)-r)
     .                           *dexp(-beta(i)*(gamma(i)-r)**2)
      end do
      epif2 = epif2 * hartree2eV      ! converts hartree to eV
      depif2 = depif2 * hartree2eV     ! converts hartree to eV
      return
      end