subroutine gwalid(iai,iaj,rij,epif3,depif3)
      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
!     Modified:
!       * March 10, 2010 (GM): add supplementary parameters
!                              from Harb's thesis
!       * July 15, 2013 (GM): correction from Harb's thesis: 2nd set
!                             of parameters
!--------------------------------------------------------------------
#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 epif3      ! output: the interaction energy
      double precision depif3

      ! local variable
      double precision a0conv
      double precision r
      double precision alpha, beta, chi, delta, epsilon
      integer i
      integer imin
      integer imax

      ! description of the function
      epif3 = 0.0d0
      depif3 = 0.0d0
      imin = min(iai,iaj)
      imax = max(iai,iaj)
      a0conv = 1.0d0/Bohr2Ang ! converts angstroms to atomic units
      r = rij*a0conv

      if (imin.eq.1.and.imax.eq.1) then
         ! H-H interaction
         alpha   =    2.47949d0
         beta    =    2.896120d0
         chi     =   47.262d0
         delta   = -505.73d0
         epsilon = 1549.81d0
       else if (imin.eq.1.and.imax.eq.6) then
         ! C-H interaction
! PIF version: 1st (no C-N and no N-N interaction, part 1 of Walid Harb's thesis)
!        alpha   =    0.42691d0
!        beta    =    0.765030d0
!        chi     = -377.328d0
!        delta   = 4668.31d0
!        epsilon =-7005.27 d0
! PIF version: 2nd (C-N and N-N interaction + modification of H-H interaction: part 2 of Walid Harb's thesis)
         alpha   =    0.53126d0
         beta    =    0.830561d0
         chi     = -353.470d0
         delta   = 4643.40d0
         epsilon =-7083.63d0
      else if (imin.eq.1.and.imax.eq.7) then
         ! N-H interaction  
         alpha   =    33.11381d0
         beta    =     1.782993d0
         chi     =  -172.487d0
         delta   =   158.43d0
         epsilon =  2423.97d0    
      else if (imin.eq.1.and.imax.eq.8) then
         ! O-H interaction
         alpha   =   29.32517d0
         beta    =    2.092648d0
         chi     =  -55.779d0
         delta   =   44.53d0
         epsilon =  313.44d0
!      else if (imin.eq.6.and.imax.eq.6) then
!        ! C-C interaction
!        alpha   = 
!        beta    = 
!        chi     = 
!        delta   = 
!        epsilon = 
       else if (imin.eq.6.and.imax.eq.7) then
         ! C-N interaction
         alpha   =     0.94697
         beta    =     0.981346
         chi     =     2.085
         delta   =     6.12        
         epsilon =   102.30
       else if (imin.eq.6.and.imax.eq.8) then
         ! O-C interaction
         alpha   =     0.87884d0
         beta    =     0.886032d0
         chi     =     1.879d0
         delta   =    33.95d0
         epsilon =   843.67d0
       else if (imin.eq.7.and.imax.eq.7) then
         ! N-N interaction
         alpha   =     1.43533
         beta    =     0.841558
         chi     =     1.941
         delta   =    35.45 
         epsilon =   887.74
       else if (imin.eq.7.and.imax.eq.8) then
         ! O-N interaction
         alpha   =     1.94702d0
         beta    =     1.016194d0
         chi     =     1.941d0
         delta   =    35.45d0
         epsilon =   887.74d0
      else if (imin.eq.8.and.imax.eq.8) then
         ! O-O interaction
         alpha   =    75.15466d0
         beta    =     1.512063d0
         chi     =   338.656d0
         delta   = -32185.68d0
         epsilon =274634.10d0
      else 
         write(*,'("error: bad use of pif2 function")')
         write(*,'("error: pif2 call with atomic number",i5," and",i5)')
     .         imin, imax
         stop
      endif 

  
        epif3  =   alpha*dexp(-beta*r) + chi/r**6 + delta/r**8
     &           + epsilon/r**10
        depif3 = - alpha*beta*dexp(-beta*r) - 6.0d0*chi/r**7
     &           - 8.0d0*delta/r**9 - 10.0d0*epsilon/r**11
     .                      
   
      epif3 = epif3 * hartree2eV      ! converts hartree to eV
      depif3 = depif3 * hartree2eV     ! converts hartree to eV
      return
      end