      subroutine numxc(npoints,natoms,ncntrt,nconts,ncontp,ncontd,
     &                 ilfunc,icfunc,ngauss,nptrs,nptrp,nptrd,natomtype,
     &                 nptsatom,iwkvec,cent,alpha,coeff,fock1,fock2,
     &                 fock3,xpoints,ypoints,zpoints,weights,xcpota,
     &                 xcpotb,xcnrgy,xcpot1,xcpot2,xcenergy,qf,qftrans,
     &                 extent1,extent2,extent3,dcc,dump1,dump2,dump3,
     &                 wkvec)

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 natomtype(*),nptsatom(*),iwkvec(*)

      dimension cent(3,*),alpha(*),coeff(*),fock1(*),fock2(*),fock3(*)

      dimension xpoints(*),ypoints(*),zpoints(*),weights(*)
      dimension xcpota(*),xcpotb(*),xcnrgy(*)
      dimension xcpot1(*),xcpot2(*),xcenergy(*)
      dimension qf(npoints,*),qftrans(ncntrt,*)
      dimension extent1(*),extent2(*),extent3(*),dcc(natoms,*)
      dimension dump1(*),dump2(*),dump3(*),wkvec(*)

      data zero,cutoff,tolerance/0.0d0,1.0d-18,1.0d-16/

      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

      ncntsq=(ncntrt*(ncntrt+1))/2

      do 1001 i=1,ncntsq
      dump1(i)=zero
      dump2(i)=zero
 1001 dump3(i)=zero

      do 1002 i=1,natoms
      do 1002 j=1,natoms
 1002 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 1003 iatom=1,natoms

      if(natomtype(iatom).eq.0) goto 1003

      if(myid.eq.0) then
                          call fastrd(11,xpoints,nptsatom(iatom))
                          call fastrd(11,ypoints,nptsatom(iatom))
                          call fastrd(11,zpoints,nptsatom(iatom))

                          call fastrd(11,weights,nptsatom(iatom))

                          call fastrd(12,xcpota,nptsatom(iatom))
                          call fastrd(12,xcpotb,nptsatom(iatom))

                          call fastrd(12,xcnrgy,nptsatom(iatom))
                    endif

      do 1004 istart=1,nptsatom(iatom),npoints

      iend=min(istart+npoints-1,nptsatom(iatom))

      radius=sqrt((cent(1,iatom)-xpoints(iend))**2
     &           +(cent(2,iatom)-ypoints(iend))**2
     &           +(cent(3,iatom)-zpoints(iend))**2)

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

      npts=iend-istart+1

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

      if(myid.eq.0) then

              ja=istart

              islave=0

              do 1005 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(xpoints(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(ypoints(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(zpoints(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(weights(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(xcpota(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(xcpotb(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)
              call mpi_bsend(xcnrgy(ja),nloop,mpi_real8,
     &                       islave,islave,mpi_comm_world,ierr)

              ja=ja+nloop

              last_slave=islave

 1005         islave=islave+1

              do 1006 i=islave,numprocs-1
              ii=0
 1006         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 1004

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

      jb=ja+npts-1

      call mpi_recv(xpoints(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(ypoints(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(zpoints(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(weights(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(xcpota(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(xcpotb(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(xcnrgy(ja),npts,
     &              mpi_real8,0,mpi_any_tag,mpi_comm_world,istat,ierr)

      do 1007 i=ja,jb
      xcpot1(i-ja+1)=xcpota(i)*weights(i)
      xcpot2(i-ja+1)=xcpotb(i)*weights(i)
 1007 xcenergy(i-ja+1)=xcnrgy(i)*weights(i)

      call compg0(iatom,ncntrt,npts,npoints,natoms,nconts,ncontp,ncontd,
     &            ilfunc,icfunc,ngauss,nptrs,nptrp,nptrd,radius,cent,
     &            alpha,coeff,xpoints(ja),ypoints(ja),zpoints(ja),qf,
     &            extent1,extent3,dcc,wkvec(00*npoints+1),
     &            wkvec(01*npoints+1),wkvec(02*npoints+1),
     &            wkvec(03*npoints+1),wkvec(04*npoints+1),
     &            wkvec(05*npoints+1),wkvec(06*npoints+1),
     &            wkvec(07*npoints+1),wkvec(08*npoints+1),
     &            wkvec(09*npoints+1),wkvec(10*npoints+1),
     &            wkvec(11*npoints+1))

      ihits=0

      do 1008 i=1,ncntrt
      if(extent3(i).gt.tolerance) then

                       ihits=ihits+1
                       iwkvec(ihits)=i

                       ii=i*(i-1)/2

                       ka=0*npts
                       kb=1*npts
                       kc=2*npts

                       do 1009 k=1,npts
                       wkvec(ka+k)=qf(k,i)*xcpot1(k)
                       wkvec(kb+k)=qf(k,i)*xcpot2(k)
 1009                  wkvec(kc+k)=qf(k,i)*xcenergy(k)

                       do 1010 j=1,ihits
                       jj=iwkvec(j)
                       if(extent3(i)*extent3(jj).lt.tolerance) goto 1010
                       ij=ii+jj
                       do 1011 k=1,npts
 1011                  dump1(ij)=dump1(ij)+qf(k,jj)*wkvec(ka+k)
                       do 1012 k=1,npts
 1012                  dump2(ij)=dump2(ij)+qf(k,jj)*wkvec(kb+k)
                       do 1013 k=1,npts
 1013                  dump3(ij)=dump3(ij)+qf(k,jj)*wkvec(kc+k)
 1010                  continue

                                  endif
 1008 continue

 1004 continue
 1003 continue

      call mpi_reduce(dump1,wkvec,ncntsq,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) then
                          do 1014 i=1,ncntsq
 1014                     fock1(i)=fock1(i)+wkvec(i)
                    endif

      call mpi_bcast(fock1,ncntsq,mpi_real8,0,mpi_comm_world,ierr)

      call mpi_reduce(dump2,wkvec,ncntsq,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) then
                          do 1015 i=1,ncntsq
 1015                     fock2(i)=fock2(i)+wkvec(i)
                    endif

      call mpi_bcast(fock2,ncntsq,mpi_real8,0,mpi_comm_world,ierr)

      call mpi_reduce(dump3,wkvec,ncntsq,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) then
                          do 1016 i=1,ncntsq
 1016                     fock3(i)=fock3(i)+wkvec(i)
                    endif

      call mpi_bcast(fock3,ncntsq,mpi_real8,0,mpi_comm_world,ierr)

      return
      end
