      subroutine cdovxx(ncds,ncdspd,ndim,ilfunc,icfunc,alpha,coeffs,
     &                  coeffp,coeffd,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(*),coeffs(*),coeffp(*),coeffd(*),cdsmat(ndim,*)
      dimension cent(3,*)

      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)

      dnorm1=sqrt(forty9/three)
      dnorm2=    (forty9/three)

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

      n=(ncdspd*(ncdspd+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,ncdspd
      do 1002 j=1,i

      ij=ij+1

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

      ilabel=ilfunc(ncds+i)
      jlabel=ilfunc(ncds+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)

      wxqx=(zetaa*cent(1,nc1)+zetab*cent(1,nc2))*zetaab-cent(1,nc2)
      wyqy=(zetaa*cent(2,nc1)+zetab*cent(2,nc2))*zetaab-cent(2,nc2)
      wzqz=(zetaa*cent(3,nc1)+zetab*cent(3,nc2))*zetaab-cent(3,nc2)

      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(5,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)
      sasa3=g(4)*factor*sqrt(zetaab)/(zetaa*zetab)
      sasa4=g(5)*factor*sqrt(zetaab)/(zetaa*zetab)

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

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

      pasa2=wxpx*sasa3
      pbsa2=wypy*sasa3
      pcsa2=wzpz*sasa3

      pasa3=wxpx*sasa4
      pbsa3=wypy*sasa4
      pcsa3=wzpz*sasa4

      sapa0=wxqx*sasa1
      sapb0=wyqy*sasa1
      sapc0=wzqz*sasa1

      sapa1=wxqx*sasa2
      sapb1=wyqy*sasa2
      sapc1=wzqz*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

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

      dasa1=wxpx*pasa2+corr
      dbsa1=wypy*pasa2
      dcsa1=wzpz*pasa2
      ddsa1=wypy*pbsa2+corr
      desa1=wzpz*pbsa2
      dfsa1=wzpz*pcsa2+corr

      corr=(sasa2-zetab*zetaab*sasa3)/(two*zetaa)

      dasa2=wxpx*pasa3+corr
      dbsa2=wypy*pasa3
      dcsa2=wzpz*pasa3
      ddsa2=wypy*pbsa3+corr
      desa2=wzpz*pbsa3
      dfsa2=wzpz*pcsa3+corr

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

      sada0=wxqx*sapa1+corr
      sadb0=wyqy*sapa1
      sadc0=wzqz*sapa1
      sadd0=wyqy*sapb1+corr
      sade0=wzqz*sapb1
      sadf0=wzqz*sapc1+corr

      corr=zetaab*sasa1/two

      papa0=wxqx*pasa1+corr
      papb0=wyqy*pasa1
      papc0=wzqz*pasa1
      pbpa0=wxqx*pbsa1
      pbpb0=wyqy*pbsa1+corr
      pbpc0=wzqz*pbsa1
      pcpa0=wxqx*pcsa1
      pcpb0=wyqy*pcsa1
      pcpc0=wzqz*pcsa1+corr

      corr=zetaab*sasa2/two

      papa1=wxqx*pasa2+corr
      papb1=wyqy*pasa2
      papc1=wzqz*pasa2
      pbpa1=wxqx*pbsa2
      pbpb1=wyqy*pbsa2+corr
      pbpc1=wzqz*pbsa2
      pcpa1=wxqx*pcsa2
      pcpb1=wyqy*pcsa2
      pcpc1=wzqz*pcsa2+corr

      corrx=zetaab*pasa1/two
      corry=zetaab*pbsa1/two
      corrz=zetaab*pcsa1/two

      dapa0=wxqx*dasa1+corrx*two
      dapb0=wyqy*dasa1
      dapc0=wzqz*dasa1
      dbpa0=wxqx*dbsa1+corry
      dbpb0=wyqy*dbsa1+corrx
      dbpc0=wzqz*dbsa1
      dcpa0=wxqx*dcsa1+corrz
      dcpb0=wyqy*dcsa1
      dcpc0=wzqz*dcsa1+corrx
      ddpa0=wxqx*ddsa1
      ddpb0=wyqy*ddsa1+corry*two
      ddpc0=wzqz*ddsa1
      depa0=wxqx*desa1
      depb0=wyqy*desa1+corrz
      depc0=wzqz*desa1+corry
      dfpa0=wxqx*dfsa1
      dfpb0=wyqy*dfsa1
      dfpc0=wzqz*dfsa1+corrz*two

      corrx=zetaab*pasa2/two
      corry=zetaab*pbsa2/two
      corrz=zetaab*pcsa2/two

      dapa1=wxqx*dasa2+corrx*two
      dapb1=wyqy*dasa2
      dapc1=wzqz*dasa2
      dbpa1=wxqx*dbsa2+corry
      dbpb1=wyqy*dbsa2+corrx
      dbpc1=wzqz*dbsa2
      dcpa1=wxqx*dcsa2+corrz
      dcpb1=wyqy*dcsa2
      dcpc1=wzqz*dcsa2+corrx
      ddpa1=wxqx*ddsa2
      ddpb1=wyqy*ddsa2+corry*two
      ddpc1=wzqz*ddsa2
      depa1=wxqx*desa2
      depb1=wyqy*desa2+corrz
      depc1=wzqz*desa2+corry
      dfpa1=wxqx*dfsa2
      dfpb1=wyqy*dfsa2
      dfpc1=wzqz*dfsa2+corrz*two

      xcorr=(pasa0-zetaa*zetaab*pasa1)/(two*zetab)
      ycorr=(pbsa0-zetaa*zetaab*pbsa1)/(two*zetab)
      zcorr=(pcsa0-zetaa*zetaab*pcsa1)/(two*zetab)

      corrx=zetaab*sapa1/two
      corry=zetaab*sapb1/two
      corrz=zetaab*sapc1/two

      pada0=wxqx*papa1+xcorr+corrx
      padb0=wyqy*papa1
      padc0=wzqz*papa1
      padd0=wyqy*papb1+xcorr
      pade0=wzqz*papb1
      padf0=wzqz*papc1+xcorr
      pbda0=wxqx*pbpa1+ycorr
      pbdb0=wyqy*pbpa1      +corrx
      pbdc0=wzqz*pbpa1
      pbdd0=wyqy*pbpb1+ycorr+corry
      pbde0=wzqz*pbpb1
      pbdf0=wzqz*pbpc1+ycorr
      pcda0=wxqx*pcpa1+zcorr
      pcdb0=wyqy*pcpa1
      pcdc0=wzqz*pcpa1      +corrx
      pcdd0=wyqy*pcpb1+zcorr
      pcde0=wzqz*pcpb1      +corry
      pcdf0=wzqz*pcpc1+zcorr+corrz

      xxcorr=(dasa0-zetaa*zetaab*dasa1)/(two*zetab)
      xycorr=(dbsa0-zetaa*zetaab*dbsa1)/(two*zetab)
      xzcorr=(dcsa0-zetaa*zetaab*dcsa1)/(two*zetab)
      yycorr=(ddsa0-zetaa*zetaab*ddsa1)/(two*zetab)
      yzcorr=(desa0-zetaa*zetaab*desa1)/(two*zetab)
      zzcorr=(dfsa0-zetaa*zetaab*dfsa1)/(two*zetab)

      xcorrx=zetaab*papa1/two
      xcorry=zetaab*papb1/two
      xcorrz=zetaab*papc1/two
      ycorrx=zetaab*pbpa1/two
      ycorry=zetaab*pbpb1/two
      ycorrz=zetaab*pbpc1/two
      zcorrx=zetaab*pcpa1/two
      zcorry=zetaab*pcpb1/two
      zcorrz=zetaab*pcpc1/two

      dada0=wxqx*dapa1+xxcorr+xcorrx*two
      dadb0=wyqy*dapa1
      dadc0=wzqz*dapa1
      dadd0=wyqy*dapb1+xxcorr
      dade0=wzqz*dapb1
      dadf0=wzqz*dapc1+xxcorr
      dbda0=wxqx*dbpa1+xycorr+ycorrx
      dbdb0=wyqy*dbpa1       +xcorrx
      dbdc0=wzqz*dbpa1
      dbdd0=wyqy*dbpb1+xycorr+xcorry
      dbde0=wzqz*dbpb1
      dbdf0=wzqz*dbpc1+xycorr
      dcda0=wxqx*dcpa1+xzcorr+zcorrx
      dcdb0=wyqy*dcpa1
      dcdc0=wzqz*dcpa1       +xcorrx
      dcdd0=wyqy*dcpb1+xzcorr
      dcde0=wzqz*dcpb1       +xcorry
      dcdf0=wzqz*dcpc1+xzcorr+xcorrz
      ddda0=wxqx*ddpa1+yycorr
      dddb0=wyqy*ddpa1       +ycorrx*two
      dddc0=wzqz*ddpa1
      dddd0=wyqy*ddpb1+yycorr+ycorry*two
      ddde0=wzqz*ddpb1
      dddf0=wzqz*ddpc1+yycorr
      deda0=wxqx*depa1+yzcorr
      dedb0=wyqy*depa1       +zcorrx
      dedc0=wzqz*depa1       +ycorrx
      dedd0=wyqy*depb1+yzcorr+zcorry
      dede0=wzqz*depb1       +ycorry
      dedf0=wzqz*depc1+yzcorr+ycorrz
      dfda0=wxqx*dfpa1+zzcorr
      dfdb0=wyqy*dfpa1
      dfdc0=wzqz*dfpa1       +zcorrx*two
      dfdd0=wyqy*dfpb1+zzcorr
      dfde0=wzqz*dfpb1       +zcorry*two
      dfdf0=wzqz*dfpc1+zzcorr+zcorrz*two

      m=ncds+10*(i-1)
      n=ncds+10*(j-1)

      cdsmat(m+01,n+01)=coeffs(ilabel)*coeffs(jlabel)*sasa0
      cdsmat(m+01,n+02)=coeffs(ilabel)*coeffp(jlabel)*sapa0
      cdsmat(m+01,n+03)=coeffs(ilabel)*coeffp(jlabel)*sapb0
      cdsmat(m+01,n+04)=coeffs(ilabel)*coeffp(jlabel)*sapc0
      cdsmat(m+01,n+05)=coeffs(ilabel)*coeffd(jlabel)*sada0/dnorm1
      cdsmat(m+01,n+06)=coeffs(ilabel)*coeffd(jlabel)*sadb0
      cdsmat(m+01,n+07)=coeffs(ilabel)*coeffd(jlabel)*sadc0
      cdsmat(m+01,n+08)=coeffs(ilabel)*coeffd(jlabel)*sadd0/dnorm1
      cdsmat(m+01,n+09)=coeffs(ilabel)*coeffd(jlabel)*sade0
      cdsmat(m+01,n+10)=coeffs(ilabel)*coeffd(jlabel)*sadf0/dnorm1

      cdsmat(m+02,n+01)=coeffp(ilabel)*coeffs(jlabel)*pasa0
      cdsmat(m+02,n+02)=coeffp(ilabel)*coeffp(jlabel)*papa0
      cdsmat(m+02,n+03)=coeffp(ilabel)*coeffp(jlabel)*papb0
      cdsmat(m+02,n+04)=coeffp(ilabel)*coeffp(jlabel)*papc0
      cdsmat(m+02,n+05)=coeffp(ilabel)*coeffd(jlabel)*pada0/dnorm1
      cdsmat(m+02,n+06)=coeffp(ilabel)*coeffd(jlabel)*padb0
      cdsmat(m+02,n+07)=coeffp(ilabel)*coeffd(jlabel)*padc0
      cdsmat(m+02,n+08)=coeffp(ilabel)*coeffd(jlabel)*padd0/dnorm1
      cdsmat(m+02,n+09)=coeffp(ilabel)*coeffd(jlabel)*pade0
      cdsmat(m+02,n+10)=coeffp(ilabel)*coeffd(jlabel)*padf0/dnorm1

      cdsmat(m+03,n+01)=coeffp(ilabel)*coeffs(jlabel)*pbsa0
      cdsmat(m+03,n+02)=coeffp(ilabel)*coeffp(jlabel)*pbpa0
      cdsmat(m+03,n+03)=coeffp(ilabel)*coeffp(jlabel)*pbpb0
      cdsmat(m+03,n+04)=coeffp(ilabel)*coeffp(jlabel)*pbpc0
      cdsmat(m+03,n+05)=coeffp(ilabel)*coeffd(jlabel)*pbda0/dnorm1
      cdsmat(m+03,n+06)=coeffp(ilabel)*coeffd(jlabel)*pbdb0
      cdsmat(m+03,n+07)=coeffp(ilabel)*coeffd(jlabel)*pbdc0
      cdsmat(m+03,n+08)=coeffp(ilabel)*coeffd(jlabel)*pbdd0/dnorm1
      cdsmat(m+03,n+09)=coeffp(ilabel)*coeffd(jlabel)*pbde0
      cdsmat(m+03,n+10)=coeffp(ilabel)*coeffd(jlabel)*pbdf0/dnorm1

      cdsmat(m+04,n+01)=coeffp(ilabel)*coeffs(jlabel)*pcsa0
      cdsmat(m+04,n+02)=coeffp(ilabel)*coeffp(jlabel)*pcpa0
      cdsmat(m+04,n+03)=coeffp(ilabel)*coeffp(jlabel)*pcpb0
      cdsmat(m+04,n+04)=coeffp(ilabel)*coeffp(jlabel)*pcpc0
      cdsmat(m+04,n+05)=coeffp(ilabel)*coeffd(jlabel)*pcda0/dnorm1
      cdsmat(m+04,n+06)=coeffp(ilabel)*coeffd(jlabel)*pcdb0
      cdsmat(m+04,n+07)=coeffp(ilabel)*coeffd(jlabel)*pcdc0
      cdsmat(m+04,n+08)=coeffp(ilabel)*coeffd(jlabel)*pcdd0/dnorm1
      cdsmat(m+04,n+09)=coeffp(ilabel)*coeffd(jlabel)*pcde0
      cdsmat(m+04,n+10)=coeffp(ilabel)*coeffd(jlabel)*pcdf0/dnorm1

      cdsmat(m+05,n+01)=coeffd(ilabel)*coeffs(jlabel)*dasa0/dnorm1
      cdsmat(m+05,n+02)=coeffd(ilabel)*coeffp(jlabel)*dapa0/dnorm1
      cdsmat(m+05,n+03)=coeffd(ilabel)*coeffp(jlabel)*dapb0/dnorm1
      cdsmat(m+05,n+04)=coeffd(ilabel)*coeffp(jlabel)*dapc0/dnorm1
      cdsmat(m+05,n+05)=coeffd(ilabel)*coeffd(jlabel)*dada0/dnorm2
      cdsmat(m+05,n+06)=coeffd(ilabel)*coeffd(jlabel)*dadb0/dnorm1
      cdsmat(m+05,n+07)=coeffd(ilabel)*coeffd(jlabel)*dadc0/dnorm1
      cdsmat(m+05,n+08)=coeffd(ilabel)*coeffd(jlabel)*dadd0/dnorm2
      cdsmat(m+05,n+09)=coeffd(ilabel)*coeffd(jlabel)*dade0/dnorm1
      cdsmat(m+05,n+10)=coeffd(ilabel)*coeffd(jlabel)*dadf0/dnorm2

      cdsmat(m+06,n+01)=coeffd(ilabel)*coeffs(jlabel)*dbsa0
      cdsmat(m+06,n+02)=coeffd(ilabel)*coeffp(jlabel)*dbpa0
      cdsmat(m+06,n+03)=coeffd(ilabel)*coeffp(jlabel)*dbpb0
      cdsmat(m+06,n+04)=coeffd(ilabel)*coeffp(jlabel)*dbpc0
      cdsmat(m+06,n+05)=coeffd(ilabel)*coeffd(jlabel)*dbda0/dnorm1
      cdsmat(m+06,n+06)=coeffd(ilabel)*coeffd(jlabel)*dbdb0
      cdsmat(m+06,n+07)=coeffd(ilabel)*coeffd(jlabel)*dbdc0
      cdsmat(m+06,n+08)=coeffd(ilabel)*coeffd(jlabel)*dbdd0/dnorm1
      cdsmat(m+06,n+09)=coeffd(ilabel)*coeffd(jlabel)*dbde0
      cdsmat(m+06,n+10)=coeffd(ilabel)*coeffd(jlabel)*dbdf0/dnorm1

      cdsmat(m+07,n+01)=coeffd(ilabel)*coeffs(jlabel)*dcsa0
      cdsmat(m+07,n+02)=coeffd(ilabel)*coeffp(jlabel)*dcpa0
      cdsmat(m+07,n+03)=coeffd(ilabel)*coeffp(jlabel)*dcpb0
      cdsmat(m+07,n+04)=coeffd(ilabel)*coeffp(jlabel)*dcpc0
      cdsmat(m+07,n+05)=coeffd(ilabel)*coeffd(jlabel)*dcda0/dnorm1
      cdsmat(m+07,n+06)=coeffd(ilabel)*coeffd(jlabel)*dcdb0
      cdsmat(m+07,n+07)=coeffd(ilabel)*coeffd(jlabel)*dcdc0
      cdsmat(m+07,n+08)=coeffd(ilabel)*coeffd(jlabel)*dcdd0/dnorm1
      cdsmat(m+07,n+09)=coeffd(ilabel)*coeffd(jlabel)*dcde0
      cdsmat(m+07,n+10)=coeffd(ilabel)*coeffd(jlabel)*dcdf0/dnorm1

      cdsmat(m+08,n+01)=coeffd(ilabel)*coeffs(jlabel)*ddsa0/dnorm1
      cdsmat(m+08,n+02)=coeffd(ilabel)*coeffp(jlabel)*ddpa0/dnorm1
      cdsmat(m+08,n+03)=coeffd(ilabel)*coeffp(jlabel)*ddpb0/dnorm1
      cdsmat(m+08,n+04)=coeffd(ilabel)*coeffp(jlabel)*ddpc0/dnorm1
      cdsmat(m+08,n+05)=coeffd(ilabel)*coeffd(jlabel)*ddda0/dnorm2
      cdsmat(m+08,n+06)=coeffd(ilabel)*coeffd(jlabel)*dddb0/dnorm1
      cdsmat(m+08,n+07)=coeffd(ilabel)*coeffd(jlabel)*dddc0/dnorm1
      cdsmat(m+08,n+08)=coeffd(ilabel)*coeffd(jlabel)*dddd0/dnorm2
      cdsmat(m+08,n+09)=coeffd(ilabel)*coeffd(jlabel)*ddde0/dnorm1
      cdsmat(m+08,n+10)=coeffd(ilabel)*coeffd(jlabel)*dddf0/dnorm2

      cdsmat(m+09,n+01)=coeffd(ilabel)*coeffs(jlabel)*desa0
      cdsmat(m+09,n+02)=coeffd(ilabel)*coeffp(jlabel)*depa0
      cdsmat(m+09,n+03)=coeffd(ilabel)*coeffp(jlabel)*depb0
      cdsmat(m+09,n+04)=coeffd(ilabel)*coeffp(jlabel)*depc0
      cdsmat(m+09,n+05)=coeffd(ilabel)*coeffd(jlabel)*deda0/dnorm1
      cdsmat(m+09,n+06)=coeffd(ilabel)*coeffd(jlabel)*dedb0
      cdsmat(m+09,n+07)=coeffd(ilabel)*coeffd(jlabel)*dedc0
      cdsmat(m+09,n+08)=coeffd(ilabel)*coeffd(jlabel)*dedd0/dnorm1
      cdsmat(m+09,n+09)=coeffd(ilabel)*coeffd(jlabel)*dede0
      cdsmat(m+09,n+10)=coeffd(ilabel)*coeffd(jlabel)*dedf0/dnorm1

      cdsmat(m+10,n+01)=coeffd(ilabel)*coeffs(jlabel)*dfsa0/dnorm1
      cdsmat(m+10,n+02)=coeffd(ilabel)*coeffp(jlabel)*dfpa0/dnorm1
      cdsmat(m+10,n+03)=coeffd(ilabel)*coeffp(jlabel)*dfpb0/dnorm1
      cdsmat(m+10,n+04)=coeffd(ilabel)*coeffp(jlabel)*dfpc0/dnorm1
      cdsmat(m+10,n+05)=coeffd(ilabel)*coeffd(jlabel)*dfda0/dnorm2
      cdsmat(m+10,n+06)=coeffd(ilabel)*coeffd(jlabel)*dfdb0/dnorm1
      cdsmat(m+10,n+07)=coeffd(ilabel)*coeffd(jlabel)*dfdc0/dnorm1
      cdsmat(m+10,n+08)=coeffd(ilabel)*coeffd(jlabel)*dfdd0/dnorm2
      cdsmat(m+10,n+09)=coeffd(ilabel)*coeffd(jlabel)*dfde0/dnorm1
      cdsmat(m+10,n+10)=coeffd(ilabel)*coeffd(jlabel)*dfdf0/dnorm2

 1002 continue

      return
      end
