********************************************************************************
*                                                                              *
*  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 ### zz1 ###
c
c     Part of this routine has been taken from GRASP2 package (ver. 1.00,  
c     Copyright (C) 1992 by F.A. Parpia, I.P. Grant, and C.F. Fischer)
c
c     Evaluate the nuclear potential for Fermi models.  
c                                                                      
      function zz1(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/

      if (z1.eq.0.d0) then
         zz1=0.d0
         return
      endif

      fmtoau=1.0d-13/ainfcm

c     potential energy is V(i,j)=-zz1(i,j)/r1 -z2/r2
c     set atomic weight 

      atw=z1atmass
      at3=atw**(one/three)
      rrmsfm = 0.836d 00*at3+0.570d 00
      tfm = 2.30d 00

c     change units from fm into bohr
      rrms = rrmsfm*fmtoau
      t = tfm*fmtoau  
      a = t/(4.0d 00*log (3.0d 00))         
      facto1 = rrms**2-(7.0d 00/5.0d 00)*(pii**2)*(a**2)
      c = sqrt (5.0d 00/3.0d 00) * sqrt (facto1) 

      abc = a/c
      tabc = two*abc
      abc2 = abc*abc
      thabc2 = three*abc2
      abc3 = abc2*abc
      cba = c/a
      pi2 = pii*pii
      hpiac2 = half*pi2*abc2
      six = two*three
      h3 = half*three
      h3php = h3+hpiac2                                              
      call es (-cba,s2mcba,s3mcba)
      sabc3 = six*abc3
      dmsas = -sabc3*s3mcba
      en = one + abc2*pi2 + dmsas
      zbn = z1/en

      ri = r*(vxi(i)+veta(j))/2.d0
      rmc = ri-c
      rmcba = rmc/a
      rbc = ri/c
      if (rbc .le. one) then
         call es (rmcba,s2rcba,s3rcba)
         zz1 = zbn*( dmsas + sabc3*s3rcba
     &        +rbc*( h3php-thabc2*s2rcba
     &        -half*rbc*rbc) )
      else
         call es (-rmcba,s2rcba,s3rcba)
         zz1 = z1 * ( one
     &        +thabc2 * ( rbc *s2rcba
     &        +tabc*s3rcba ) / en )
      endif
      
      if (idbg(451).ne.0) then
         if (abs(zz1-z1).gt.1.d-8) then
            write(*,*) i,j,ri,zz1
         endif
      endif
      end