      subroutine gridgrad(iatom,npoints,nmaxpts,ncenters,nconts,ncontp,
     &                    ncontd,ilfunc,icfunc,ngauss,nptrs,nptrp,nptrd,
     &                    natomtype,istore,iwkvec,cent,alpha,coeff,
     &                    wtbecke,extent,dcp,dcc,tmat,drvtv,x,y,z,
     &                    weight)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

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

      dimension ilfunc(*),icfunc(*),ngauss(*),nptrs(*),nptrp(*),nptrd(*)
      dimension natomtype(*),istore(*),iwkvec(*)

      dimension cent(3,*),alpha(*),coeff(*)

      dimension wtbecke(*),extent(*),dcp(*),dcc(ncenters,*)
      dimension tmat(ncenters,*),drvtv(nmaxpts,*)
      dimension x(*),y(*),z(*),weight(*)

      data zero,half,one,two,four,five/0.0d0,0.5d0,1.0d0,
     &                                 2.0d0,4.0d0,5.0d0/
      data sixteen,twentyone,thirtyfive,one_o_5/16.0d0,21.0d0,
     &                                          35.0d0,105.0d0/
      data tolerance,tol,bohr/0.0000000001d0,0.000001d0,0.529177d0/
      data afactor,pi/0.64d0,3.1415926535898d0/
      data alpha1,alpha2,alpha3/0.40d0,1.20d0,2.00d0/

      do 1001 i=1,ncenters
 1001 iwkvec(i)=0

      do 1002 j=1,3*ncenters
      do 1002 i=1,npoints
 1002 drvtv(i,j)=zero

      call diffuse1(ncenters,nconts,ncontp,ncontd,ilfunc,icfunc,ngauss,
     &              nptrs,nptrp,nptrd,tolerance,alpha,coeff,cent,extent)

      do 1003 katom=1,ncenters
      do 1003 jatom=1,ncenters
 1003 dcc(jatom,katom)=dist(cent(1,jatom),cent(1,katom))

      smallest=one/tolerance

      do 1004 jatom=1,ncenters
      if((jatom.eq.iatom).or.(natomtype(jatom).eq.0)) goto 1004
      if(dcc(iatom,jatom).le.smallest) smallest=dcc(iatom,jatom)
 1004 continue

      smallest=smallest*(one-afactor)*half

      do 1005 ipoint=1,npoints

      do 1006 icenter=1,ncenters
 1006 wtbecke(icenter)=zero

      dcp(iatom)=sqrt((cent(1,iatom)-x(ipoint))**2
     &                 +(cent(2,iatom)-y(ipoint))**2
     &                 +(cent(3,iatom)-z(ipoint))**2)

      if(dcp(iatom).lt.smallest) goto 1005

      ihits=0

      do 1007 icenter=1,ncenters
      rr=(cent(1,icenter)-x(ipoint))**2
     &  +(cent(2,icenter)-y(ipoint))**2
     &  +(cent(3,icenter)-z(ipoint))**2
      if((rr.le.extent(icenter)).and.(natomtype(icenter).ne.0)) then

                                 ihits=ihits+1
                                 istore(ihits+1)=icenter
                                 dcp(icenter)=sqrt(rr)

                                                                endif
 1007 continue

      if(ihits.eq.1) goto 1005

      istore(1)=ihits

      do 1008 icenter=1,ihits
      jcenter=istore(icenter+1)

      iwkvec(jcenter)=1
      iwkvec(iatom)=1

      partition=one

      do 1009 kcenter=1,ihits
      lcenter=istore(kcenter+1)
      if(lcenter.eq.jcenter) goto 1009

      u=(dcp(jcenter)-dcp(lcenter))/dcc(jcenter,lcenter)

      if(abs(u).lt.afactor) then
                                  zz=(thirtyfive*(u/afactor)
     &                               -thirtyfive*(u/afactor)**3
     &                               +twentyone*(u/afactor)**5
     &                               -five*(u/afactor)**7)/sixteen
                            else
                                  if(u.le.(-afactor)) zz=-one
                                  if(u.ge.(+afactor)) zz=+one
                            endif

      s=half*(one-zz)

      if(abs(u).gt.afactor) then

         tmat(jcenter,lcenter)=zero

                            else

         if(s.gt.tol) then

              tmat(jcenter,lcenter)=(half/s)
     &                             *(-one/sixteen)
     &                             *((thirtyfive/afactor)
     &                              -((one_o_5*(u**2))/(afactor**3))
     &                              +((one_o_5*(u**4))/(afactor**5))
     &                              -((thirtyfive*(u**6))/(afactor**7)))

                      else

              tmat(jcenter,lcenter)=zero

                      endif

                            endif

      partition=partition*s

 1009 continue

      wtbecke(jcenter)=partition

 1008 continue

      totalwtbecke=zero

      do 1010 icenter=1,ihits
      jcenter=istore(icenter+1)
 1010 totalwtbecke=totalwtbecke+wtbecke(jcenter)

      do 1011 icenter=1,ihits
      jcenter=istore(icenter+1)
      if(jcenter.eq.iatom) goto 1011

      jj=3*(jcenter-1)

      factor1=-(one/dcc(iatom,jcenter))/dcp(jcenter)
      factor2=(dcp(iatom)-dcp(jcenter))/(dcc(iatom,jcenter)**3)

      duabx=factor1*(cent(1,jcenter)-x(ipoint))
     &     +factor2*(cent(1,iatom)-cent(1,jcenter))
      duaby=factor1*(cent(2,jcenter)-y(ipoint))
     &     +factor2*(cent(2,iatom)-cent(2,jcenter))
      duabz=factor1*(cent(3,jcenter)-z(ipoint))
     &     +factor2*(cent(3,iatom)-cent(3,jcenter))

      factor=wtbecke(iatom)*tmat(iatom,jcenter)/totalwtbecke

      drvtv(ipoint,jj+1)=drvtv(ipoint,jj+1)+factor*duabx
      drvtv(ipoint,jj+2)=drvtv(ipoint,jj+2)+factor*duaby
      drvtv(ipoint,jj+3)=drvtv(ipoint,jj+3)+factor*duabz

      do 1012 kcenter=1,ihits
      lcenter=istore(kcenter+1)
      if(lcenter.eq.jcenter) goto 1012

      jj=3*(jcenter-1)

      factor1=-(one/dcc(lcenter,jcenter))/dcp(jcenter)
      factor2=(dcp(lcenter)-dcp(jcenter))/(dcc(lcenter,jcenter)**3)

      ducbx=factor1*(cent(1,jcenter)-x(ipoint))
     &     +factor2*(cent(1,lcenter)-cent(1,jcenter))
      ducby=factor1*(cent(2,jcenter)-y(ipoint))
     &     +factor2*(cent(2,lcenter)-cent(2,jcenter))
      ducbz=factor1*(cent(3,jcenter)-z(ipoint))
     &     +factor2*(cent(3,lcenter)-cent(3,jcenter))

      factor=wtbecke(iatom)*wtbecke(lcenter)
     &      *tmat(lcenter,jcenter)/(totalwtbecke**2)

      drvtv(ipoint,jj+1)=drvtv(ipoint,jj+1)-factor*ducbx
      drvtv(ipoint,jj+2)=drvtv(ipoint,jj+2)-factor*ducby
      drvtv(ipoint,jj+3)=drvtv(ipoint,jj+3)-factor*ducbz

 1012 continue

      do 1013 kcenter=1,ihits
      lcenter=istore(kcenter+1)
      if(lcenter.eq.jcenter) goto 1013

      jj=3*(jcenter-1)

      factor1=-(one/dcc(jcenter,lcenter))/dcp(lcenter)
      factor2=(dcp(jcenter)-dcp(lcenter))/(dcc(jcenter,lcenter)**3)

      dubcx=factor1*(cent(1,lcenter)-x(ipoint))
     &     +factor2*(cent(1,jcenter)-cent(1,lcenter))
      dubcy=factor1*(cent(2,lcenter)-y(ipoint))
     &     +factor2*(cent(2,jcenter)-cent(2,lcenter))
      dubcz=factor1*(cent(3,lcenter)-z(ipoint))
     &     +factor2*(cent(3,jcenter)-cent(3,lcenter))

      factor=wtbecke(iatom)*wtbecke(jcenter)
     &      *tmat(jcenter,lcenter)/(totalwtbecke**2)

      drvtv(ipoint,jj+1)=drvtv(ipoint,jj+1)+factor*dubcx
      drvtv(ipoint,jj+2)=drvtv(ipoint,jj+2)+factor*dubcy
      drvtv(ipoint,jj+3)=drvtv(ipoint,jj+3)+factor*dubcz

 1013 continue

 1011 continue

      do 1014 j=1,3
      jj=3*(iatom-1)+j
      do 1014 icenter=1,ihits
      jcenter=istore(icenter+1)
      if(jcenter.eq.iatom) goto 1014
      kk=3*(jcenter-1)+j
      drvtv(ipoint,jj)=drvtv(ipoint,jj)-drvtv(ipoint,kk)
 1014 continue

      wght_adjust=zero

      if(wtbecke(iatom).gt.tol) wght_adjust=weight(ipoint)*totalwtbecke
     &                                     /wtbecke(iatom)

      do 1015 icenter=1,ihits
      jcenter=istore(icenter+1)
      do 1015 j=1,3
      jj=3*(jcenter-1)+j
 1015 drvtv(ipoint,jj)=drvtv(ipoint,jj)*wght_adjust

 1005 continue

      return
      end
