      subroutine cdovss(ncds,ndim,ilfunc,icfunc,alpha,coeff,cdsmat,cent)

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 g(13)

      dimension ilfunc(*),icfunc(*)

      dimension alpha(*),coeff(*),cdsmat(ndim,*),cent(3,*)

      data one,two,pi/1.0d0,2.0d0,3.1415926535898d0/
      data four,five/4.0d0,5.0d0/

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

      factor=(pi**(five/four)*sqrt(two))**2

      n=(ncds*(ncds+1))/2

      if(myid.eq.0) then
                          isent=0

                          do 1001 islave=0,numprocs-1

                          ia=n-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

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

 1001                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(isent,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)

      istart=isent+1
      iend=isent+nloop

      ij=0

      do 1002 i=1,ncds
      do 1002 j=1,i

      ij=ij+1

      if((ij.lt.istart).or.(ij.gt.iend)) goto 1002

      ilabel=ilfunc(i)
      jlabel=ilfunc(j)

      zetaa=alpha(ilabel)
      zetab=alpha(jlabel)

      zetaab=one/(zetaa+zetab)

      axbx=cent(1,icfunc(ilabel))-cent(1,icfunc(jlabel))
      ayby=cent(2,icfunc(ilabel))-cent(2,icfunc(jlabel))
      azbz=cent(3,icfunc(ilabel))-cent(3,icfunc(jlabel))

      u=zetaa*zetab*zetaab*(axbx*axbx+ayby*ayby+azbz*azbz)

      call augg(1,u,g)

      cdsmat(i,j)=coeff(ilabel)*coeff(jlabel)*g(1)*factor*sqrt(zetaab)
     &           /(zetaa*zetab)

 1002 continue

      return
      end
