      subroutine cdcalls(nreal,ncds,ncdspd,ndim,nmaxatoms,natoms,
     &                   nbsystems,isys,icdlsave,ilfunc,icfunc,iwkvec,
     &                   cutoff,alpha,coeffs,coeffp,coeffd,cdsmat,cent,
     &                   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 isys(nmaxatoms,*),icdlsave(*),ilfunc(*),icfunc(*),
     &          iwkvec(*)

      dimension alpha(*),coeffs(*),coeffp(*),coeffd(*),cdsmat(ndim,*)
      dimension cent(3,*),wkvec(*)

      data zero/0.0d0/

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

      ntotal=ncds+10*ncdspd

      do 1001 i=1,ndim
      do 1001 j=1,i
 1001 cdsmat(i,j)=zero

      call cdovss(ncds,ntotal,ilfunc,icfunc,alpha,coeffs,cdsmat,cent)

      call cdovxs(ncds,ncdspd,ntotal,ilfunc,icfunc,alpha,coeffs,coeffp,
     &            coeffd,cdsmat,cent)

      call cdovxx(ncds,ncdspd,ntotal,ilfunc,icfunc,alpha,coeffs,coeffp,
     &            coeffd,cdsmat,cent)

      call mpi_reduce(cdsmat,wkvec,ntotal*ntotal,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) then
                          do 1002 i=1,ntotal
                          do 1002 j=1,i
                          ij=(j-1)*ntotal+i
                          cdsmat(i,j)=wkvec(ij)
 1002                     cdsmat(j,i)=wkvec(ij)

                          rewind 17

                          do 1003 i=1,ntotal
 1003                     call fastwr(17,cdsmat(1,i),ntotal)
                    endif

      ncds_save=ncds
      ncdspd_save=ncdspd

      do 1004 i=1,ndim
 1004 icdlsave(i)=ilfunc(i)

      rewind 18

      do 1005 i=1,nbsystems

      do 1006 j=1,natoms
 1006 iwkvec(j)=0

      do 1007 j=1,isys(1,i)
 1007 iwkvec(isys(j+1,i))=1

      do 1008 j=1,natoms
      do 1008 k=1,isys(1,i)
 1008 if(dist(cent(1,j),cent(1,isys(k+1,i))).le.cutoff) iwkvec(j)=1

      ncds=0

      do 1009 j=1,ncds_save
      if(iwkvec(icfunc(icdlsave(j))).eq.1) then

                               ncds=ncds+1
                               ilfunc(ncds)=icdlsave(j)

                                           endif
 1009 continue

      ncdspd=0

      do 1010 j=1,ncdspd_save
      if(iwkvec(icfunc(icdlsave(ncds_save+j))).eq.1) then

                               ncdspd=ncdspd+1
                               ilfunc(ncds+ncdspd)=icdlsave(ncds_save+j)

                                                     endif
 1010 continue

      ntotal=ncds+10*ncdspd

      do 1011 j=1,ndim
      do 1011 k=1,j
 1011 cdsmat(j,k)=zero

      call cdovss(ncds,ntotal,ilfunc,icfunc,alpha,coeffs,cdsmat,cent)

      call cdovxs(ncds,ncdspd,ntotal,ilfunc,icfunc,alpha,coeffs,coeffp,
     &            coeffd,cdsmat,cent)

      call cdovxx(ncds,ncdspd,ntotal,ilfunc,icfunc,alpha,coeffs,coeffp,
     &            coeffd,cdsmat,cent)

      call mpi_reduce(cdsmat,wkvec,ntotal*ntotal,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) then
                          do 1012 j=1,ntotal
                          do 1012 k=1,j
                          jk=(k-1)*ntotal+j
                          cdsmat(j,k)=wkvec(jk)
 1012                     cdsmat(k,j)=wkvec(jk)

                          do 1013 k=1,ntotal
                          do 1013 j=1,ntotal
 1013                     wkvec(ntotal*(k-1)+j)=cdsmat(j,k)
                    endif

      call mpi_bcast(wkvec,ntotal**2,mpi_real8,0,mpi_comm_world,ierr)

      ja=1
      jb=ja+ntotal*ntotal
      jc=jb+ntotal*ntotal
      jd=jc+ntotal*ntotal
      je=jd+ntotal*ntotal
      jf=je+ntotal
      jg=jf+ntotal
      jh=jg+ntotal
      ji=jh+ntotal

      if(ji.gt.nreal) stop 'real overflow in cdcalls'

      call x_annihilate(18,ntotal,iwkvec,iwkvec(numprocs+1),wkvec(ja),
     &                  cdsmat,wkvec(jb),wkvec(jc),wkvec(jd),wkvec(je),
     &                  wkvec(jf),wkvec(jg),wkvec(jh))

      ncds=ncds_save
      ncdspd=ncdspd_save

      do 1014 j=1,ndim
 1014 ilfunc(j)=icdlsave(j)

 1005 continue

      return
      end
