      subroutine cdovxs(ncds,ncdspd,ndim,ilfunc,icfunc,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 g(13)

      dimension ilfunc(*),icfunc(*)

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

      data one,two,three/1.0d0,2.0d0,3.0d0/
      data four,five,forty9/4.0d0,5.0d0,49.0d0/
      data pi/3.1415926535898d0/

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

      sqrt493=sqrt(forty9/three)

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

      n=ncds*ncdspd

      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,ncdspd
      do 1002 j=1,ncds

      ij=ij+1

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

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

      nc1=icfunc(ilabel)
      nc2=icfunc(jlabel)

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

      zetaab=one/(zetaa+zetab)

      wxpx=(zetaa*cent(1,nc1)+zetab*cent(1,nc2))*zetaab-cent(1,nc1)
      wypy=(zetaa*cent(2,nc1)+zetab*cent(2,nc2))*zetaab-cent(2,nc1)
      wzpz=(zetaa*cent(3,nc1)+zetab*cent(3,nc2))*zetaab-cent(3,nc1)

      axbx=cent(1,nc1)-cent(1,nc2)
      ayby=cent(2,nc1)-cent(2,nc2)
      azbz=cent(3,nc1)-cent(3,nc2)

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

      call augg(3,u,g)

      sasa0=g(1)*factor*sqrt(zetaab)/(zetaa*zetab)
      sasa1=g(2)*factor*sqrt(zetaab)/(zetaa*zetab)
      sasa2=g(3)*factor*sqrt(zetaab)/(zetaa*zetab)

      pasa0=wxpx*sasa1
      pbsa0=wypy*sasa1
      pcsa0=wzpz*sasa1

      pasa1=wxpx*sasa2
      pbsa1=wypy*sasa2
      pcsa1=wzpz*sasa2

      corr=(sasa0-zetab*zetaab*sasa1)/(two*zetaa)

      dasa0=wxpx*pasa1+corr
      dbsa0=wypy*pasa1
      dcsa0=wzpz*pasa1
      ddsa0=wypy*pbsa1+corr
      desa0=wzpz*pbsa1
      dfsa0=wzpz*pcsa1+corr

      n=ncds+10*(i-1)

      cdsmat(n+01,j)=coeffs(ilabel)*coeffs(jlabel)*sasa0
      cdsmat(n+02,j)=coeffp(ilabel)*coeffs(jlabel)*pasa0
      cdsmat(n+03,j)=coeffp(ilabel)*coeffs(jlabel)*pbsa0
      cdsmat(n+04,j)=coeffp(ilabel)*coeffs(jlabel)*pcsa0
      cdsmat(n+05,j)=coeffd(ilabel)*coeffs(jlabel)*dasa0/sqrt493
      cdsmat(n+06,j)=coeffd(ilabel)*coeffs(jlabel)*dbsa0
      cdsmat(n+07,j)=coeffd(ilabel)*coeffs(jlabel)*dcsa0
      cdsmat(n+08,j)=coeffd(ilabel)*coeffs(jlabel)*ddsa0/sqrt493
      cdsmat(n+09,j)=coeffd(ilabel)*coeffs(jlabel)*desa0
      cdsmat(n+10,j)=coeffd(ilabel)*coeffs(jlabel)*dfsa0/sqrt493

 1002 continue

      return
      end
