      subroutine gridxc(iatom,natoms,npts,ndim,n,n1,n2,ilaux,icaux,
     &                  radius,cent,alpha,coeffs,coeffp,coeffd,x,y,z,q,
     &                  extent2,extent3,dcc,xx,yy,zz,rr,sfac)

c  this subroutine written by alain st-amant of the department of
c  pharmaceutical chemistry, university of california, san francisco.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      dimension ilaux(*),icaux(*)

      dimension cent(3,*),alpha(*),coeffs(*),coeffp(*),coeffd(*)
      dimension x(*),y(*),z(*),q(ndim,*),xx(*),yy(*),zz(*),rr(*),sfac(*)
      dimension extent2(*),extent3(*),dcc(natoms,*)

      data zero,one,three,expcut/0.0d0,1.0d0,3.0d0,50.0d0/

      sqrt3=sqrt(three)

      do 1001 i=1,n
 1001 extent3(i)=zero

c  now for lone-s exponent xc fitting functions

      k=1

      do 1002 i=1,n1

      ilabel=ilaux(i)

      nc=icaux(ilabel)

      distmin=dcc(nc,iatom)-radius
      if(distmin.gt.extent2(i)) goto 1002

      do 1003 j=1,npts
      xx(j)=x(j)-cent(1,nc)
      yy(j)=y(j)-cent(2,nc)
 1003 zz(j)=z(j)-cent(3,nc)

      do 1004 j=1,npts
 1004 rr(j)=(xx(j)*xx(j)+yy(j)*yy(j)+zz(j)*zz(j))*alpha(ilabel)

      do 1005 j=1,npts
      if(rr(j).lt.expcut) then
                                q(j,k)=exp(-rr(j))*coeffs(ilabel)
                          else
                                q(j,k)=zero
                          endif
 1005 continue

      extent3(k)=sumsum(npts,q(1,k))/real64(npts)

 1002 k=k+1

      if(n2.eq.n1) return

c  now for shared exponent xc fitting functions

      do 1006 i=n1+1,n2

      ilabel=ilaux(i)

      nc=icaux(ilabel)

      distmin=dcc(nc,iatom)-radius
      if(distmin.gt.extent2(i)) goto 1006

      do 1007 j=1,npts
      xx(j)=x(j)-cent(1,nc)
      yy(j)=y(j)-cent(2,nc)
 1007 zz(j)=z(j)-cent(3,nc)

      do 1008 j=1,npts
 1008 rr(j)=(xx(j)*xx(j)+yy(j)*yy(j)+zz(j)*zz(j))*alpha(ilabel)

      do 1009 j=1,npts
      if(rr(j).lt.expcut) then

                   sfac(j)=exp(-rr(j))

                   q(j,k+0)=sfac(j)*coeffs(ilabel)

                   q(j,k+1)=xx(j)*sfac(j)*coeffp(ilabel)
                   q(j,k+2)=yy(j)*sfac(j)*coeffp(ilabel)
                   q(j,k+3)=zz(j)*sfac(j)*coeffp(ilabel)

                   q(j,k+4)=xx(j)*xx(j)*sfac(j)*coeffd(ilabel)/sqrt3
                   q(j,k+5)=xx(j)*yy(j)*sfac(j)*coeffd(ilabel)
                   q(j,k+6)=xx(j)*zz(j)*sfac(j)*coeffd(ilabel)
                   q(j,k+7)=yy(j)*yy(j)*sfac(j)*coeffd(ilabel)/sqrt3
                   q(j,k+8)=yy(j)*zz(j)*sfac(j)*coeffd(ilabel)
                   q(j,k+9)=zz(j)*zz(j)*sfac(j)*coeffd(ilabel)/sqrt3

                          else

                   q(j,k+0)=zero

                   q(j,k+1)=zero
                   q(j,k+2)=zero
                   q(j,k+3)=zero

                   q(j,k+4)=zero
                   q(j,k+5)=zero
                   q(j,k+6)=zero
                   q(j,k+7)=zero
                   q(j,k+8)=zero
                   q(j,k+9)=zero

                          endif
 1009 continue

      extent3(k+0)=sumsum(npts,q(1,k+0))/real64(npts)
      extent3(k+1)=sumsum(npts,q(1,k+1))/real64(npts)
      extent3(k+2)=sumsum(npts,q(1,k+2))/real64(npts)
      extent3(k+3)=sumsum(npts,q(1,k+3))/real64(npts)
      extent3(k+4)=sumsum(npts,q(1,k+4))/real64(npts)
      extent3(k+5)=sumsum(npts,q(1,k+5))/real64(npts)
      extent3(k+6)=sumsum(npts,q(1,k+6))/real64(npts)
      extent3(k+7)=sumsum(npts,q(1,k+7))/real64(npts)
      extent3(k+8)=sumsum(npts,q(1,k+8))/real64(npts)
      extent3(k+9)=sumsum(npts,q(1,k+9))/real64(npts)

 1006 k=k+10

      return
      end
