      subroutine setuppp(m,n,ncontp,nshelp,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)

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)

      dimension nshelp(*),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(*)

      data zero,expcut/0.0d0,50.0d0/

      icnptr=0

      n=0

      do 1001 i=1,ncontp
      do 1001 j=1,i

      icnptr=icnptr+1

      iminj(icnptr)=i-j

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

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

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

      inddyy(icnptr,1)=i4
      inddyy(icnptr,2)=i4+1
      inddyy(icnptr,3)=i4+2
      inddyy(icnptr,4)=i4+i1
      inddyy(icnptr,5)=i4+i1+1
      inddyy(icnptr,6)=i4+i1+2
      inddyy(icnptr,7)=i4+i1+i1+1
      inddyy(icnptr,8)=i4+i1+i1+2
      inddyy(icnptr,9)=i4+i1+i1+3

      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 1001

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

      test=alpha(ii)*alpha(jj)*rr/(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
 1002 continue
 1001 continue

      return
      end
