      subroutine numxcc(nreal,ninteger,nmaxatoms,nmaxpoints,natoms,
     &                  nconts,ncontp,ncontd,ncds,ncdspd,nxcs,nxcspd,
     &                  nbsystems,level,iwkvec,mtloca,ilfunc,ngauss,
     &                  icfunc,nptrs,nptrp,nptrd,natomtype,nptsatom,
     &                  icdlfunc,icdcfunc,ixclfunc,ixccfunc,isys,cutoff,
     &                  wkvec,coord,core,fock1,fock2,fock3,alpha,coeff,
     &                  alphacd,coefscd,coefpcd,coefdcd,alphaxc,coefsxc,
     &                  coefpxc,coefdxc,cdfitc,xcfit1,xcfit2,xcfit3)

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 iwkvec(*),mtloca(*),ilfunc(*),ngauss(*),icfunc(*),
     &          nptrs(*),nptrp(*),nptrd(*),natomtype(*),nptsatom(*),
     &          icdlfunc(*),icdcfunc(*),ixclfunc(*),ixccfunc(*),
     &          isys(nmaxatoms,*)

      dimension wkvec(*),core(*),fock1(*),fock2(*),fock3(*),coord(3,*),
     &          alpha(*),coeff(*),alphacd(*),coefscd(*),coefpcd(*),
     &          coefdcd(*),alphaxc(*),coefsxc(*),coefpxc(*),coefdxc(*),
     &          cdfitc(*),xcfit1(*),xcfit2(*),xcfit3(*)

      data zero,one/0.0d0,1.0d0/

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

      ncntrt=nconts+3*ncontp+6*ncontd

      ncntsq=ncntrt*(ncntrt+1)/2

      do 1001 i=1,ncntsq
      fock1(i)=zero
      fock2(i)=zero
 1001 fock3(i)=zero

      ncdfunc=ncds+10*ncdspd
      nxcfunc=nxcs+10*nxcspd

      call twoemm(nreal-3*ncntsq,ninteger,nconts,ncontp,ncontd,ncds,
     &            ncdspd,nxcs,nxcspd,nptrs,nptrp,nptrd,mtloca,ilfunc,
     &            ngauss,icfunc,icdlfunc,icdcfunc,ixclfunc,ixccfunc,
     &            iwkvec,cutoff,coord,alpha,coeff,alphacd,coefscd,
     &            coefpcd,coefdcd,alphaxc,coefsxc,coefpxc,coefdxc,
     &            cdfitc,xcfit1,xcfit2,xcfit3,wkvec(1),fock1,fock2,
     &            fock3,wkvec(2*ncntsq+1),wkvec(3*ncntsq+1))

      do 1002 i=1,ncntsq
      fock1(i)=fock1(i)+core(i)
      fock2(i)=fock2(i)+core(i)
 1002 fock3(i)=fock3(i)+core(i)

      nstart=min(512,32*numprocs)

      npoints=nstart+16

 1003 npoints=npoints-16

      if(npoints.le.0) stop 'numxcc real overflow'

      ja=1
      jb=ja+nmaxpoints
      jc=jb+nmaxpoints
      jd=jc+nmaxpoints
      je=jd+nmaxpoints
      jf=je+nmaxpoints
      jg=jf+nmaxpoints
      jh=jg+nmaxpoints
      ji=jh+npoints
      jj=ji+npoints
      jk=jj+npoints
      jl=jk+ncntrt*npoints
      jm=jl+ncntrt*npoints
      jn=jm+nconts+ncontp+ncontd
      jo=jn+natoms
      jp=jo+ncntrt
      jq=jp+natoms*natoms
      jr=jq+(ncntrt*(ncntrt+1))/2
      js=jr+(ncntrt*(ncntrt+1))/2
      jt=js+(ncntrt*(ncntrt+1))/2

      if((jt+12*npoints).gt.nreal) goto 1003

      call numxc(npoints,natoms,ncntrt,nconts,ncontp,ncontd,ilfunc,
     &           icfunc,ngauss,nptrs,nptrp,nptrd,natomtype,nptsatom,
     &           iwkvec,coord,alpha,coeff,fock1,fock2,fock3,wkvec(ja),
     &           wkvec(jb),wkvec(jc),wkvec(jd),wkvec(je),wkvec(jf),
     &           wkvec(jg),wkvec(jh),wkvec(ji),wkvec(jj),wkvec(jk),
     &           wkvec(jl),wkvec(jm),wkvec(jn),wkvec(jo),wkvec(jp),
     &           wkvec(jq),wkvec(jr),wkvec(js),wkvec(jt))

      return
      end
