********************************************************************************
*                                                                              * 
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### zz2g ###
c
c     Evaluate the nuclear potential for Gaussian models.  
                                                                
      function zz2g(i,j)
      implicit integer*4 (i-n)
      implicit real*8 (A-H, O-Z)
      include 'commons8.inc'
c     ainfcm: bohr radius in cm
      data ainfcm/0.529 177 249  d-08/

      fmtoau=1.0d-13/ainfcm

      if (z2.eq.0.d0) then
         zz2g=0.d0
         eta2=1d35
         return
      endif

c     set atomic weight 

      atw=z2atmass
      at3=atw**(one/three)
      rrmsfm = 0.836d 00*at3+0.570d 00

c     change units from fm into bohr

      rrms = rrmsfm*fmtoau

c     exponent of the Gaussian distribution $\rho_0 exp(-\eta r^2)$
c     rrms=(3/(2\eta)^{1/2)
      
      eta2=3d0/(2d0*rrms*rrms)
      
c     set the normalization constant $Z=rho_0 (pi/eta)^{3/2}$
c      rho0=z2/(pii/eta2)**1.5d0

      ri = r*(vxi(i)-veta(j))/2.d0

c      x=eta2*ri*ri
c      dg15= dgamit(1.5d0,x)*x**1.5d0*dgamma(1.5d0)  
c      zz2g = 2.0d0*pii*rho0*(dg15/eta2**1.5d0+ri*exp(-x)/eta2)
c
c     making use of F77 intrinsic function 
c     erf(x)=2/(Pi)^(1/2)*int_0^x exp(-t^2)dt
c     if available
c
      zz2g=z2*erf(sqrt(eta2)*ri)

      return
      end