********************************************************************************
*                                                                              * 
*  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 ### zz1g ###
c 
c     Evaluate the nuclear potential for Gaussian models.  
c                                                                
      function zz1g(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 (z1.eq.0.d0) then
         zz1g=0.d0
         eta1=1d35
         return
      endif
      
c     set atomic weight 
      
      atw=z1atmass
      at3=atw**(one/three)
      rrmsfm = 0.836d 00*at3+0.570d 00
      
c     change units from fm into bohr
      
      rrms = rrmsfm*fmtoau
c
c     exponent of the Gaussian distribution $\rho_0 exp(-\eta r^2)$
c     rrms=(3/(2\eta)^{1/2)
      
      eta1=3d0/(2d0*rrms*rrms)

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

c     set the normalization constant $Z=rho_0 (pi/eta)^{3/2}$

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

      return
      end