      subroutine setupgdp(m,n,ncontd,ncontp,nsheld,nshelp,mtloca,
     &                    iprimvc,nprimvc,icentvc,icontr,inddyy,icentr,
     &                    coord,alpha,coeff,zeta,zetaa,zetab,twozetaa,
     &                    twozetab,aax,aay,aaz,eta,eta2,etarr,base,px,
     &                    py,pz,pxax,pyay,pzaz,pxbx,pyby,pzbz,axbx,ayby,
     &                    azbz,buffer)

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

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

      data zero,two,expcut/0.0d0,2.0d0,100.0d0/

      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,ncontp

      icnptr=ncontp*(i-1)+j

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

      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

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

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

      n=n+1

      icontr(n)=icnptr

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

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

      twozetaa(n)=two*zetaa(n)
      twozetab(n)=two*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)

      eta2(n)=two*eta(n)

      etarr(n)=eta(n)*rr

      base(n)=zero

      if(etarr(n).lt.expcut) then
                                   base(n)=factor*exp(-etarr(n))
                             else
                                   n=n-1
                             endif

 1003 continue
 1002 continue

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

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

      inddyy(ncontp*(i-1)+j,01)=i4
      inddyy(ncontp*(i-1)+j,02)=i4+1
      inddyy(ncontp*(i-1)+j,03)=i4+2
      inddyy(ncontp*(i-1)+j,04)=i4+i1
      inddyy(ncontp*(i-1)+j,05)=i4+i1+1
      inddyy(ncontp*(i-1)+j,06)=i4+i1+2
      inddyy(ncontp*(i-1)+j,07)=i4+i1*2+1
      inddyy(ncontp*(i-1)+j,08)=i4+i1*2+2
      inddyy(ncontp*(i-1)+j,09)=i4+i1*2+3
      inddyy(ncontp*(i-1)+j,10)=i4+i1*3+3
      inddyy(ncontp*(i-1)+j,11)=i4+i1*3+4
      inddyy(ncontp*(i-1)+j,12)=i4+i1*3+5
      inddyy(ncontp*(i-1)+j,13)=i4+i1*4+6
      inddyy(ncontp*(i-1)+j,14)=i4+i1*4+7
      inddyy(ncontp*(i-1)+j,15)=i4+i1*4+8
      inddyy(ncontp*(i-1)+j,16)=i4+i1*5+10
      inddyy(ncontp*(i-1)+j,17)=i4+i1*5+11
      inddyy(ncontp*(i-1)+j,18)=i4+i1*5+12

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

      icnptr=ncontp*(i-1)+j

      icentr(icnptr,1)=icenta
      icentr(icnptr,2)=icentb

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

 1004 continue

      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),twozetaa)
      call stringcopy(n,buffer(01*n+1),twozetab)
      call stringcopy(n,buffer(02*n+1),etarr)
      call stringcopy(n,buffer(03*n+1),zeta)
      call stringcopy(n,buffer(04*n+1),base)
      call stringcopy(n,buffer(05*n+1),eta2)
      call stringcopy(n,buffer(06*n+1),pxax)
      call stringcopy(n,buffer(07*n+1),pyay)
      call stringcopy(n,buffer(08*n+1),pzaz)
      call stringcopy(n,buffer(09*n+1),pxbx)
      call stringcopy(n,buffer(10*n+1),pyby)
      call stringcopy(n,buffer(11*n+1),pzbz)
      call stringcopy(n,buffer(12*n+1),eta)
      call stringcopy(n,buffer(13*n+1),px)
      call stringcopy(n,buffer(14*n+1),py)
      call stringcopy(n,buffer(15*n+1),pz)

      call mpi_bsend(buffer,16*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,16*nsub,mpi_real8,
     &                         i,03,mpi_comm_world,istat,ierr)

                 call stringcopy(nsub,twozetaa(n+1),buffer(00*nsub+1))
                 call stringcopy(nsub,twozetab(n+1),buffer(01*nsub+1))
                 call stringcopy(nsub,etarr(n+1),buffer(02*nsub+1))
                 call stringcopy(nsub,zeta(n+1),buffer(03*nsub+1))
                 call stringcopy(nsub,base(n+1),buffer(04*nsub+1))
                 call stringcopy(nsub,eta2(n+1),buffer(05*nsub+1))
                 call stringcopy(nsub,pxax(n+1),buffer(06*nsub+1))
                 call stringcopy(nsub,pyay(n+1),buffer(07*nsub+1))
                 call stringcopy(nsub,pzaz(n+1),buffer(08*nsub+1))
                 call stringcopy(nsub,pxbx(n+1),buffer(09*nsub+1))
                 call stringcopy(nsub,pyby(n+1),buffer(10*nsub+1))
                 call stringcopy(nsub,pzbz(n+1),buffer(11*nsub+1))
                 call stringcopy(nsub,eta(n+1),buffer(12*nsub+1))
                 call stringcopy(nsub,px(n+1),buffer(13*nsub+1))
                 call stringcopy(nsub,py(n+1),buffer(14*nsub+1))
                 call stringcopy(nsub,pz(n+1),buffer(15*nsub+1))

                n=n+nsub

 1005           continue

                    endif

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

      return
      end
