      subroutine xcgrad(nptsbatch,nmaxpts,ngridtype,iptsweight,ncntrt,
     &                  nconts,ncontp,ncontd,ncentr,ilfunc,icfunc,
     &                  ngauss,nptrs,nptrp,nptrd,nptsatom,natnumber,
     &                  nfuncatom,istore,iwkvec,alpha,coeff,dmta,dmtb,
     &                  cent,grad,x,y,z,weight,xcalph,xcbeta,xcenrg,qf,
     &                  qx,qy,qz,dcc,extent1,extent3,drvtv,wv)

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 istat(mpi_status_size)

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

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

      dimension x(*),y(*),z(*),weight(*),xcalph(*),xcbeta(*),xcenrg(*)
      dimension qf(nptsbatch,*)
      dimension qx(nptsbatch,*),qy(nptsbatch,*),qz(nptsbatch,*)
      dimension dcc(ncentr,*),extent1(*),extent3(*),drvtv(nptsbatch,*)
      dimension wv(*)

      data cutoff,tolerance,two/1.0d-16,1.0d-18,2.0d0/

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

      if(myid.eq.0) then
                          rewind 11
                          rewind 12
                    endif

      do 1001 i=1,ncentr
      do 1001 j=1,ncentr
 1001 dcc(j,i)=dist(cent(1,j),cent(1,i))

      call diffuse2(nconts,ncontp,ncontd,ilfunc,icfunc,ngauss,nptrs,
     &              nptrp,nptrd,cutoff,alpha,coeff,cent,extent1)

      do 1002 iatom=1,ncentr

      if(natnumber(iatom).eq.0) goto 1002

      if(myid.eq.0) then
                          call fastrd(11,x,nptsatom(iatom))
                          call fastrd(11,y,nptsatom(iatom))
                          call fastrd(11,z,nptsatom(iatom))
                          call fastrd(11,weight,nptsatom(iatom))
                          call fastrd(12,xcalph,nptsatom(iatom))
                          call fastrd(12,xcbeta,nptsatom(iatom))
                          call fastrd(12,xcenrg,nptsatom(iatom))
                    endif

      do 1003 ia=1,nptsatom(iatom),nptsbatch

      radius=sqrt((cent(1,iatom)-x(ia+nptsbatch-1))**2
     &           +(cent(2,iatom)-y(ia+nptsbatch-1))**2
     &           +(cent(3,iatom)-z(ia+nptsbatch-1))**2)

      call mpi_bcast(radius,1,mpi_real8,0,mpi_comm_world,ierr)

      ib=min(ia+nptsbatch-1,nptsatom(iatom))

      npts=ib-ia+1

      nmaxloop=(npts/numprocs)+1
      if((nmaxloop*numprocs).lt.npts) nmaxloop=nmaxloop+1

      if(myid.eq.0) then

              ja=ia

              islave=0

              do 1004 i=1,npts,nmaxloop

              nloop=min(nmaxloop,npts-i+1)

              call mpi_bsend(nloop,1,mpi_integer,islave,
     &                       islave,mpi_comm_world,ierr)
              call mpi_bsend(ja,1,mpi_integer,islave,
     &                       islave,mpi_comm_world,ierr)

              call mpi_bsend(x(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(y(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(z(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(weight(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(xcalph(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(xcbeta(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(xcenrg(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)

              ja=ja+nloop

              last_slave=islave
                           
 1004         islave=islave+1

              do 1005 i=islave,numprocs-1
              ii=0
 1005         call mpi_bsend(ii,1,mpi_integer,i,i,mpi_comm_world,ierr)

                    endif

      call mpi_recv(npts,1,
     &              mpi_integer,0,mpi_any_tag,mpi_comm_world,istat,ierr)

      if(npts.eq.0) goto 1003

      call mpi_recv(ja,1,
     &              mpi_integer,0,mpi_any_tag,mpi_comm_world,istat,ierr)

      call mpi_recv(x(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(y(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(z(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(weight(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(xcalph(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(xcbeta(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(xcenrg(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)

      ka=1
      kb=ka+ncentr
      kc=kb+ncentr
      kd=kc+ncentr
      ke=kd+ncentr*ncentr

      if(iptsweight.eq.1) then

                   call gridgrad(iatom,npts,nptsbatch,ncentr,nconts,
     &                           ncontp,ncontd,ilfunc,icfunc,ngauss,
     &                           nptrs,nptrp,nptrd,natnumber,istore,
     &                           iwkvec,cent,alpha,coeff,wv(ka),wv(kb),
     &                           wv(kc),wv(kd),wv(ke),drvtv,x(ja),y(ja),
     &                           z(ja),weight(ja))

                          endif

      do 1006 i=ja,ja+npts-1
      xcalph(i)=xcalph(i)*weight(i)
 1006 xcbeta(i)=xcbeta(i)*weight(i)

      call compg1(iatom,ncntrt,npts,nptsbatch,ncentr,nconts,ncontp,
     &            ncontd,ilfunc,icfunc,ngauss,nptrs,nptrp,nptrd,radius,
     &            cent,alpha,coeff,x(ja),y(ja),z(ja),qf,qx,qy,qz,
     &            extent1,extent3,dcc,wv(00*npts+1),wv(01*npts+1),
     &            wv(02*npts+1),wv(03*npts+1),wv(04*npts+1),
     &            wv(05*npts+1),wv(06*npts+1),wv(07*npts+1),
     &            wv(08*npts+1),wv(09*npts+1),wv(10*npts+1),
     &            wv(11*npts+1),wv(12*npts+1),wv(13*npts+1),
     &            wv(14*npts+1),wv(15*npts+1),wv(16*npts+1),
     &            wv(17*npts+1),wv(18*npts+1),wv(19*npts+1),
     &            wv(20*npts+1),wv(21*npts+1),wv(22*npts+1),
     &            wv(23*npts+1),wv(24*npts+1))

      xcmax=dmaxvec(xcenrg(ja),npts)

      ihits=0

      do 1007 i=1,ncntrt
      if(extent3(i).gt.tolerance) then
                                        ihits=ihits+1
                                        iwkvec(ncentr+ihits)=i
                                  endif
 1007 continue

      do 1008 i=1,ihits
      ii=iwkvec(ncentr+i)

      icenter=nfuncatom(ii)

      do 1008 j=1,ihits
      jj=iwkvec(ncentr+j)

      ijmax=max(ii,jj)
      ijmin=min(ii,jj)
      ij=(ijmax*(ijmax-1)/2)+ijmin

      if(ii.eq.jj) then
                         dnstya=dmta(ij)*two
                         dnstyb=dmtb(ij)*two
                   else
                         dnstya=dmta(ij)
                         dnstyb=dmtb(ij)
                   endif

      do 1009 k=1,npts
 1009 wv(k)=(dnstya*xcalph(ja+k-1)+dnstyb*xcbeta(ja+k-1))*qf(k,jj)

      do 1010 k=1,npts
      grad(1,icenter)=grad(1,icenter)-qx(k,ii)*wv(k)
      grad(2,icenter)=grad(2,icenter)-qy(k,ii)*wv(k)
 1010 grad(3,icenter)=grad(3,icenter)-qz(k,ii)*wv(k)

      if((jj.gt.ii).or.(iptsweight.eq.0)) goto 1008

      dnstya=dmta(ij)
      dnstyb=dmtb(ij)

      check=(dnstya+dnstyb)*extent3(ii)*extent3(jj)*xcmax

      if(abs(check).lt.tolerance) goto 1008

      if(iwkvec(iatom).eq.1) then
                                   do 1011 k=1,npts
 1011                              wv(k)=(dnstya+dnstyb)*xcenrg(ja+k-1)
     &                                  *(qf(k,ii)*qf(k,jj))
                             endif

      do 1012 l=1,ncentr
      if(iwkvec(l).eq.1) then
                               do 1013 m=1,3
                               n=3*(l-1)+m
                               do 1013 k=1,npts
 1013                          grad(m,l)=grad(m,l)+wv(k)*drvtv(k,n)
                         endif
 1012 continue

 1008 continue
 1003 continue
 1002 continue

      return
      end
