      subroutine gradxc(ndim,ngridmax,ngridtype,iptsweight,ncntrt,
     &                  ncentr,nconts,ncontp,ncontd,natnumber,nptsatom,
     &                  nfuncatom,ilfunc,icfunc,ngauss,nptrs,nptrp,
     &                  nptrd,iwv,wkvec,wv,dmta,dmtb,cent,grad,alpha,
     &                  coeff)

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)

      include "mpif.h"

      dimension natnumber(*),nptsatom(*),nfuncatom(*)
      dimension ilfunc(*),icfunc(*),ngauss(*),nptrs(*),nptrp(*),nptrd(*)
      dimension iwv(*)

      dimension wkvec(*),wv(*),dmta(*),dmtb(*),cent(3,*),grad(3,*)
      dimension alpha(*),coeff(*)

      data zero/0.0d0/

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      ii=0

      do 1001 i=1,nconts
      nc=icfunc(nptrs(i))
      nfuncatom(ii+1)=nc
 1001 ii=ii+1

      do 1002 i=1,ncontp
      nc=icfunc(nptrp(i))
      do 1003 k=1,3
 1003 nfuncatom(ii+k)=nc
 1002 ii=ii+3

      do 1004 i=1,ncontd
      nc=icfunc(nptrd(i))
      do 1005 k=1,6
 1005 nfuncatom(ii+k)=nc
 1004 ii=ii+6

      do 1006 i=1,3*ncentr
 1006 wkvec(i)=zero

      nstart=min(512,32*numprocs)

      ngridbatch=nstart+16

 1007 ngridbatch=ngridbatch-16

      if(ngridbatch.le.0) stop 'gradxc going into xcgrad'

      ja=1
      jb=ja+ngridmax
      jc=jb+ngridmax
      jd=jc+ngridmax
      je=jd+ngridmax
      jf=je+ngridmax
      jg=jf+ngridmax
      jh=jg+ngridmax
      ji=jh+ngridbatch*ncntrt
      jj=ji+ngridbatch*ncntrt
      jk=jj+ngridbatch*ncntrt
      jl=jk+ngridbatch*ncntrt
      jm=jl+ncentr*ncentr
      jn=jm+ncntrt
      jo=jn+ncntrt
      jp=jo+ngridbatch*ncentr*3
      jq=jp+2*ncentr+2*ncentr*ncentr+ncentr*ngridmax

      ka=1
      kb=ka+ncentr+1

      if(jq.gt.ndim) goto 1007

      call xcgrad(ngridbatch,ngridmax,ngridtype,iptsweight,ncntrt,
     &            nconts,ncontp,ncontd,ncentr,ilfunc,icfunc,ngauss,
     &            nptrs,nptrp,nptrd,nptsatom,natnumber,nfuncatom,
     &            iwv(ka),iwv(kb),alpha,coeff,dmta,dmtb,cent,wkvec,
     &            wv(ja),wv(jb),wv(jc),wv(jd),wv(je),wv(jf),wv(jg),
     &            wv(jh),wv(ji),wv(jj),wv(jk),wv(jl),wv(jm),wv(jn),
     &            wv(jo),wv(jp))

      call mpi_reduce(wkvec,wkvec(3*ncentr+1),
     &                3*ncentr,mpi_real8,mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) then
                          do 1008 i=1,ncentr
                          do 1008 j=1,3
 1008                     grad(j,i)=grad(j,i)+wkvec(3*ncentr+3*(i-1)+j)
                    endif

      return
      end
