      subroutine orbitals(n,ndim,nmos,nconts,ncontp,ncontd,ilfunc,
     &                    icfunc,ngauss,nptrs,nptrp,nptrd,nfuncatom,
     &                    coord,alpha,coeff,overlap,coeffmo,smat,
     &                    centroid,extent,exvec)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
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(*),nfuncatom(*)

      dimension coord(3,*),alpha(*),coeff(*),overlap(*),coeffmo(ndim,*),
     &          smat(n,*),centroid(3,*),extent(*),exvec(*)

      data zero,one,bohr/0.0d0,1.0d0,0.529177d0/
      data cutoff,tolerance/9999.0d0,0.0000000000000001d0/

      call diffuse2(nconts,ncontp,ncontd,ilfunc,icfunc,ngauss,nptrs,
     &              nptrp,nptrd,tolerance,alpha,coeff,coord,extent)

      do 1001 i=1,nconts
 1001 exvec(i)=extent(i)

      do 1002 i=1,ncontp
      exvec(nconts+3*(i-1)+1)=extent(nconts+i)
      exvec(nconts+3*(i-1)+2)=extent(nconts+i)
 1002 exvec(nconts+3*(i-1)+3)=extent(nconts+i)

      do 1003 i=1,ncontd
      exvec(nconts+3*ncontp+6*(i-1)+1)=extent(nconts+ncontp+i)
      exvec(nconts+3*ncontp+6*(i-1)+2)=extent(nconts+ncontp+i)
      exvec(nconts+3*ncontp+6*(i-1)+3)=extent(nconts+ncontp+i)
      exvec(nconts+3*ncontp+6*(i-1)+4)=extent(nconts+ncontp+i)
      exvec(nconts+3*ncontp+6*(i-1)+5)=extent(nconts+ncontp+i)
 1003 exvec(nconts+3*ncontp+6*(i-1)+6)=extent(nconts+ncontp+i)

      ij=0

      do 1004 i=1,n
      do 1004 j=1,i
      ij=ij+1
      smat(i,j)=overlap(ij)
 1004 smat(j,i)=overlap(ij)

      do 1005 j=1,nmos
      do 1005 i=1,n
      distance=dist(coord(1,nfuncatom(i)),centroid(1,j))
 1005 if((distance.gt.exvec(i)).and.
     &   (distance.gt.(cutoff/bohr))) coeffmo(i,j)=zero

      do 1006 k=1,nmos

      svalue=zero

      do 1007 j=1,n
      do 1007 i=1,n
 1007 svalue=svalue+coeffmo(i,k)*coeffmo(j,k)*smat(i,j)

      factor=one/sqrt(svalue)

      do 1008 i=1,n
 1008 coeffmo(i,k)=coeffmo(i,k)*factor

 1006 continue

      return
      end
