      subroutine setupcddd(m,n,ncontd,nsheld,mtloca,iprimvc,nprimvc,
     &                     icentvc,iminj,icontr,inddyy,cutoff,coord,
     &                     alpha,coeff,zeta,zetaa,zetab,aax,aay,aaz,eta,
     &                     etarr,base,px,py,pz,pxax,pyay,pzaz,pxbx,pyby,
     &                     pzbz,axbx,ayby,azbz,buffer,density)

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 nsheld(*),mtloca(*),iprimvc(*),nprimvc(*),icentvc(*),
     &          iminj(*),icontr(*),inddyy(m,*)

      dimension coord(3,*),alpha(*),coeff(*),zeta(*),zetaa(*),zetab(*),
     &          aax(*),aay(*),aaz(*),eta(*),etarr(*),base(*),px(*),
     &          py(*),pz(*),pxax(*),pyay(*),pzaz(*),pxbx(*),pyby(*),
     &          pzbz(*),axbx(*),ayby(*),azbz(*),buffer(*),density(*)

      data zero,expcut,tolerance/0.0d0,60.0d0,0.00000000000000000001d0/

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

      n=0

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

                          do 1001 islave=0,numprocs-1

                          ia=ncontd-isent
                          ib=numprocs-islave

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

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

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

 1001                     isent=isent+nloop
                    endif

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

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

      do 1002 i=iloop,iloop+nloop-1
      do 1002 j=1,i

      icnptr=(i*(i-1)/2)+j

      iminj(icnptr)=i-j

      nconta=nsheld(i)
      ncontb=nsheld(j)

      matlca=mtloca(nconta)
      matlcb=mtloca(ncontb)

      i1=matlca
      i2=(i1*(i1-1))/2
      i3=matlcb
      i4=i2+i3

      inddyy(icnptr,01)=i4
      inddyy(icnptr,02)=i4+1
      inddyy(icnptr,03)=i4+2
      inddyy(icnptr,04)=i4+3
      inddyy(icnptr,05)=i4+4
      inddyy(icnptr,06)=i4+5
      inddyy(icnptr,07)=i4+i1
      inddyy(icnptr,08)=i4+i1+1
      inddyy(icnptr,09)=i4+i1+2
      inddyy(icnptr,10)=i4+i1+3
      inddyy(icnptr,11)=i4+i1+4
      inddyy(icnptr,12)=i4+i1+5
      inddyy(icnptr,13)=i4+i1*2+1
      inddyy(icnptr,14)=i4+i1*2+2
      inddyy(icnptr,15)=i4+i1*2+3
      inddyy(icnptr,16)=i4+i1*2+4
      inddyy(icnptr,17)=i4+i1*2+5
      inddyy(icnptr,18)=i4+i1*2+6
      inddyy(icnptr,19)=i4+i1*3+3
      inddyy(icnptr,20)=i4+i1*3+4
      inddyy(icnptr,21)=i4+i1*3+5
      inddyy(icnptr,22)=i4+i1*3+6
      inddyy(icnptr,23)=i4+i1*3+7
      inddyy(icnptr,24)=i4+i1*3+8
      inddyy(icnptr,25)=i4+i1*4+6
      inddyy(icnptr,26)=i4+i1*4+7
      inddyy(icnptr,27)=i4+i1*4+8
      inddyy(icnptr,28)=i4+i1*4+9
      inddyy(icnptr,29)=i4+i1*4+10
      inddyy(icnptr,30)=i4+i1*4+11
      inddyy(icnptr,31)=i4+i1*5+10
      inddyy(icnptr,32)=i4+i1*5+11
      inddyy(icnptr,33)=i4+i1*5+12
      inddyy(icnptr,34)=i4+i1*5+13
      inddyy(icnptr,35)=i4+i1*5+14
      inddyy(icnptr,36)=i4+i1*5+15

      if(i.eq.j) then
                       dmax=max(abs(density(inddyy(icnptr,01))),
     &                          abs(density(inddyy(icnptr,07))),
     &                          abs(density(inddyy(icnptr,08))),
     &                          abs(density(inddyy(icnptr,13))),
     &                          abs(density(inddyy(icnptr,14))),
     &                          abs(density(inddyy(icnptr,15))),
     &                          abs(density(inddyy(icnptr,19))),
     &                          abs(density(inddyy(icnptr,20))),
     &                          abs(density(inddyy(icnptr,21))),
     &                          abs(density(inddyy(icnptr,22))),
     &                          abs(density(inddyy(icnptr,25))),
     &                          abs(density(inddyy(icnptr,26))),
     &                          abs(density(inddyy(icnptr,27))),
     &                          abs(density(inddyy(icnptr,28))),
     &                          abs(density(inddyy(icnptr,29))),
     &                          abs(density(inddyy(icnptr,31))),
     &                          abs(density(inddyy(icnptr,32))),
     &                          abs(density(inddyy(icnptr,33))),
     &                          abs(density(inddyy(icnptr,34))),
     &                          abs(density(inddyy(icnptr,35))),
     &                          abs(density(inddyy(icnptr,36))))

                       if(dmax.lt.tolerance) goto 1002
                 else
                       dmax=max(abs(density(inddyy(icnptr,01))),
     &                          abs(density(inddyy(icnptr,02))),
     &                          abs(density(inddyy(icnptr,03))),
     &                          abs(density(inddyy(icnptr,04))),
     &                          abs(density(inddyy(icnptr,05))),
     &                          abs(density(inddyy(icnptr,06))),
     &                          abs(density(inddyy(icnptr,07))),
     &                          abs(density(inddyy(icnptr,08))),
     &                          abs(density(inddyy(icnptr,09))),
     &                          abs(density(inddyy(icnptr,10))),
     &                          abs(density(inddyy(icnptr,11))),
     &                          abs(density(inddyy(icnptr,12))),
     &                          abs(density(inddyy(icnptr,13))),
     &                          abs(density(inddyy(icnptr,14))),
     &                          abs(density(inddyy(icnptr,15))),
     &                          abs(density(inddyy(icnptr,16))),
     &                          abs(density(inddyy(icnptr,17))),
     &                          abs(density(inddyy(icnptr,18))),
     &                          abs(density(inddyy(icnptr,19))),
     &                          abs(density(inddyy(icnptr,20))),
     &                          abs(density(inddyy(icnptr,21))),
     &                          abs(density(inddyy(icnptr,22))),
     &                          abs(density(inddyy(icnptr,23))),
     &                          abs(density(inddyy(icnptr,24))),
     &                          abs(density(inddyy(icnptr,25))),
     &                          abs(density(inddyy(icnptr,26))),
     &                          abs(density(inddyy(icnptr,27))),
     &                          abs(density(inddyy(icnptr,28))),
     &                          abs(density(inddyy(icnptr,29))),
     &                          abs(density(inddyy(icnptr,30))),
     &                          abs(density(inddyy(icnptr,31))),
     &                          abs(density(inddyy(icnptr,32))),
     &                          abs(density(inddyy(icnptr,33))),
     &                          abs(density(inddyy(icnptr,34))),
     &                          abs(density(inddyy(icnptr,35))),
     &                          abs(density(inddyy(icnptr,36))))

                       if(dmax.lt.tolerance) goto 1002
                 endif

      iprima=iprimvc(nconta)
      iprimb=iprimvc(ncontb)

      nprima=iprima+nprimvc(nconta)-1
      nprimb=iprimb+nprimvc(ncontb)-1

      icenta=icentvc(nconta)
      icentb=icentvc(ncontb)

      ax=coord(1,icenta)
      ay=coord(2,icenta)
      az=coord(3,icenta)

      bx=coord(1,icentb)
      by=coord(2,icentb)
      bz=coord(3,icentb)

      caxbx=ax-bx
      cayby=ay-by
      cazbz=az-bz

      axbx(icnptr)=caxbx
      ayby(icnptr)=cayby
      azbz(icnptr)=cazbz

      rr=caxbx*caxbx+cayby*cayby+cazbz*cazbz

      if(rr.gt.(cutoff**2)) goto 1002

      do 1003 ii=iprima,nprima
      do 1003 jj=iprimb,nprimb

      test=rr*alpha(ii)*alpha(jj)/(alpha(ii)+alpha(jj))

      if(test.lt.expcut) then
                               n=n+1

                               icontr(n)=icnptr

                               zetaa(n)=alpha(ii)
                               zetab(n)=alpha(jj)

                               zeta(n)=zetaa(n)+zetab(n)

                               factor=coeff(ii)*coeff(jj)

                               aax(n)=ax
                               aay(n)=ay
                               aaz(n)=az

                               px(n)=(zetaa(n)*ax+zetab(n)*bx)/zeta(n)
                               py(n)=(zetaa(n)*ay+zetab(n)*by)/zeta(n)
                               pz(n)=(zetaa(n)*az+zetab(n)*bz)/zeta(n)

                               pxax(n)=px(n)-ax
                               pyay(n)=py(n)-ay
                               pzaz(n)=pz(n)-az

                               pxbx(n)=px(n)-bx
                               pyby(n)=py(n)-by
                               pzbz(n)=pz(n)-bz

                               eta(n)=zetaa(n)*zetab(n)/zeta(n)

                               etarr(n)=eta(n)*rr

                               base(n)=factor*exp(-etarr(n))
                         endif
 1003 continue
 1002 continue

      do 1004 i=1,ncontd
      do 1004 j=1,i

      iminj((i*(i-1)/2)+j)=i-j

      i1=mtloca(nsheld(i))
      i2=(i1*(i1-1))/2
      i3=mtloca(nsheld(j))
      i4=i2+i3

      inddyy((i*(i-1)/2)+j,01)=i4
      inddyy((i*(i-1)/2)+j,02)=i4+1
      inddyy((i*(i-1)/2)+j,03)=i4+2
      inddyy((i*(i-1)/2)+j,04)=i4+3
      inddyy((i*(i-1)/2)+j,05)=i4+4
      inddyy((i*(i-1)/2)+j,06)=i4+5
      inddyy((i*(i-1)/2)+j,07)=i4+i1
      inddyy((i*(i-1)/2)+j,08)=i4+i1+1
      inddyy((i*(i-1)/2)+j,09)=i4+i1+2
      inddyy((i*(i-1)/2)+j,10)=i4+i1+3
      inddyy((i*(i-1)/2)+j,11)=i4+i1+4
      inddyy((i*(i-1)/2)+j,12)=i4+i1+5
      inddyy((i*(i-1)/2)+j,13)=i4+i1*2+1
      inddyy((i*(i-1)/2)+j,14)=i4+i1*2+2
      inddyy((i*(i-1)/2)+j,15)=i4+i1*2+3
      inddyy((i*(i-1)/2)+j,16)=i4+i1*2+4
      inddyy((i*(i-1)/2)+j,17)=i4+i1*2+5
      inddyy((i*(i-1)/2)+j,18)=i4+i1*2+6
      inddyy((i*(i-1)/2)+j,19)=i4+i1*3+3
      inddyy((i*(i-1)/2)+j,20)=i4+i1*3+4
      inddyy((i*(i-1)/2)+j,21)=i4+i1*3+5
      inddyy((i*(i-1)/2)+j,22)=i4+i1*3+6
      inddyy((i*(i-1)/2)+j,23)=i4+i1*3+7
      inddyy((i*(i-1)/2)+j,24)=i4+i1*3+8
      inddyy((i*(i-1)/2)+j,25)=i4+i1*4+6
      inddyy((i*(i-1)/2)+j,26)=i4+i1*4+7
      inddyy((i*(i-1)/2)+j,27)=i4+i1*4+8
      inddyy((i*(i-1)/2)+j,28)=i4+i1*4+9
      inddyy((i*(i-1)/2)+j,29)=i4+i1*4+10
      inddyy((i*(i-1)/2)+j,30)=i4+i1*4+11
      inddyy((i*(i-1)/2)+j,31)=i4+i1*5+10
      inddyy((i*(i-1)/2)+j,32)=i4+i1*5+11
      inddyy((i*(i-1)/2)+j,33)=i4+i1*5+12
      inddyy((i*(i-1)/2)+j,34)=i4+i1*5+13
      inddyy((i*(i-1)/2)+j,35)=i4+i1*5+14
      inddyy((i*(i-1)/2)+j,36)=i4+i1*5+15

      icenta=icentvc(nsheld(i))
      icentb=icentvc(nsheld(j))

      icnptr=(i*(i-1)/2)+j

      axbx(icnptr)=coord(1,icenta)-coord(1,icentb)
      ayby(icnptr)=coord(2,icenta)-coord(2,icentb)
 1004 azbz(icnptr)=coord(3,icenta)-coord(3,icentb)

      call mpi_bsend(n,1,mpi_integer,0,01,mpi_comm_world,ierr)

      call mpi_bsend(icontr,n,mpi_integer,0,02,mpi_comm_world,ierr)

      call stringcopy(n,buffer(00*n+1),zeta)
      call stringcopy(n,buffer(01*n+1),base)
      call stringcopy(n,buffer(02*n+1),aax)
      call stringcopy(n,buffer(03*n+1),aay)
      call stringcopy(n,buffer(04*n+1),aaz)
      call stringcopy(n,buffer(05*n+1),px)
      call stringcopy(n,buffer(06*n+1),py)
      call stringcopy(n,buffer(07*n+1),pz)
      call stringcopy(n,buffer(08*n+1),pxax)
      call stringcopy(n,buffer(09*n+1),pyay)
      call stringcopy(n,buffer(10*n+1),pzaz)

      call mpi_bsend(buffer,11*n,mpi_real8,0,03,mpi_comm_world,ierr)

      if(myid.eq.0) then

                    n=0

                    do 1005 i=0,numprocs-1

                    call mpi_recv(nsub,1,mpi_integer,
     &                            i,01,mpi_comm_world,istat,ierr)

                    call mpi_recv(icontr(n+1),nsub,mpi_integer,
     &                            i,02,mpi_comm_world,istat,ierr)

                    call mpi_recv(buffer,11*nsub,mpi_real8,
     &                            i,03,mpi_comm_world,istat,ierr)

                    call stringcopy(nsub,zeta(n+1),buffer(0*nsub+1))
                    call stringcopy(nsub,base(n+1),buffer(1*nsub+1))
                    call stringcopy(nsub,aax(n+1),buffer(2*nsub+1))
                    call stringcopy(nsub,aay(n+1),buffer(3*nsub+1))
                    call stringcopy(nsub,aaz(n+1),buffer(4*nsub+1))
                    call stringcopy(nsub,px(n+1),buffer(5*nsub+1))
                    call stringcopy(nsub,py(n+1),buffer(6*nsub+1))
                    call stringcopy(nsub,pz(n+1),buffer(7*nsub+1))
                    call stringcopy(nsub,pxax(n+1),buffer(08*nsub+1))
                    call stringcopy(nsub,pyay(n+1),buffer(09*nsub+1))
                    call stringcopy(nsub,pzaz(n+1),buffer(10*nsub+1))

                    n=n+nsub

 1005               continue

                    endif

      call mpi_bcast(n,1,mpi_integer,0,mpi_comm_world,ierr)

      return
      end
