      subroutine diffuse1(ncenters,nconts,ncontp,ncontd,ilfunc,icfunc,
     &                    ngauss,nptrs,nptrp,nptrd,tolerance,alpha,
     &                    coeff,cent,extent)

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 ilfunc(*),icfunc(*),ngauss(*),nptrs(*),nptrp(*),nptrd(*)

      dimension alpha(*),coeff(*),cent(3,*),extent(*)

      data nextent/50/
      data zero,one,expcut/0.0d0,1.0d0,50.0d0/

      do 1001 i=1,ncenters
 1001 extent(i)=zero

      do 1002 iloop=1,nextent

      xx=real64(iloop)
      rr=xx*xx

c  loop over s functions

      do 1003 i=1,nconts

      s=zero

      ilabel=nptrs(i)

      iptr1=ilfunc(ilabel)
      iptr2=iptr1+ngauss(ilabel)-1

      nc=icfunc(ilabel)

      do 1004 j=iptr1,iptr2
      rralpha=rr*alpha(j)
      if(rralpha.gt.expcut) goto 1004
      s=s+coeff(j)*exp(-rralpha)
 1004 continue

      if(abs(s).gt.tolerance) then
                                    if(xx.gt.extent(nc)) extent(nc)=xx
                              endif

 1003 continue

c  loop over p functions

      do 1005 i=1,ncontp

      px=zero

      ilabel=nptrp(i)

      iptr1=ilfunc(ilabel)
      iptr2=iptr1+ngauss(ilabel)-1

      nc=icfunc(ilabel)

      do 1006 j=iptr1,iptr2
      rralpha=rr*alpha(j)
      if(rralpha.gt.expcut) goto 1006
      s=coeff(j)*exp(-rralpha)
      px=px+xx*s
 1006 continue

      if(abs(px).gt.tolerance) then
                                     if(xx.gt.extent(nc)) extent(nc)=xx
                               endif

 1005 continue

c  loop over d functions

      do 1007 i=1,ncontd

      dxx=zero

      ilabel=nptrd(i)

      iptr1=ilfunc(ilabel)
      iptr2=iptr1+ngauss(ilabel)-1

      nc=icfunc(ilabel)

      do 1008 j=iptr1,iptr2
      rralpha=rr*alpha(j)
      if(rralpha.gt.expcut) goto 1008
      s=coeff(j)*exp(-rralpha)
      dxx=dxx+xx*xx*s
 1008 continue

      if(abs(dxx).gt.tolerance) then
                                      if(xx.gt.extent(nc)) extent(nc)=xx
                                endif

 1007 continue

 1002 continue

      do 1009 i=1,ncenters
 1009 if(extent(i).ge.real64(nextent-1)) stop 'extent may be a problem'

      do 1010 i=1,ncenters
 1010 extent(i)=(extent(i)+one)**2

      return
      end
