      subroutine boys_loc(n,ndim,nocc,coeff,x_moment,y_moment,z_moment,
     &                    veci,vecj,wkmat,centroid)

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 coeff(ndim,*),x_moment(n,*),y_moment(n,*),z_moment(n,*),
     &          veci(*),vecj(*),wkmat(n,*),centroid(3,*)

      data nmaxloop/200/
      data zero,quarter,one/0.0d0,0.25d0,1.0d0/
      data pi,tol/3.1415926535898d0,0.00000001d0/

      dvalue=pi

      do 1001 j=1,n
      do 1001 i=1,n
      wkmat(i,j)=zero
      do 1001 k=1,n
 1001 wkmat(i,j)=wkmat(i,j)+coeff(k,j)*x_moment(k,i)

      do 1002 j=1,n
      do 1002 i=1,n
 1002 x_moment(i,j)=wkmat(i,j)

      do 1003 j=1,n
      do 1003 i=1,n
      wkmat(i,j)=zero
      do 1003 k=1,n
 1003 wkmat(i,j)=wkmat(i,j)+coeff(k,j)*y_moment(k,i)

      do 1004 j=1,n
      do 1004 i=1,n
 1004 y_moment(i,j)=wkmat(i,j)

      do 1005 j=1,n
      do 1005 i=1,n
      wkmat(i,j)=zero
      do 1005 k=1,n
 1005 wkmat(i,j)=wkmat(i,j)+coeff(k,j)*z_moment(k,i)

      do 1006 j=1,n
      do 1006 i=1,n
 1006 z_moment(i,j)=wkmat(i,j)

      do 1007 iloop=1,nmaxloop

      dvalue_save=dvalue

      dvalue=zero

      do 1008 i=1,nocc

      x_rii=zero
      y_rii=zero
      z_rii=zero

      do 1009 jj=1,n
      x_rii=x_rii+coeff(jj,i)*x_moment(jj,i)
      y_rii=y_rii+coeff(jj,i)*y_moment(jj,i)
 1009 z_rii=z_rii+coeff(jj,i)*z_moment(jj,i)

 1008 dvalue=dvalue+(x_rii*x_rii+y_rii*y_rii+z_rii*z_rii)

      if(abs(dvalue-dvalue_save).lt.tol) goto 1030

      do 1010 i=1,nocc
      do 1010 j=1,i

      if(i.eq.j) goto 1010

      x_rii=zero
      y_rii=zero
      z_rii=zero

      do 1011 jj=1,n
      x_rii=x_rii+coeff(jj,i)*x_moment(jj,i)
      y_rii=y_rii+coeff(jj,i)*y_moment(jj,i)
 1011 z_rii=z_rii+coeff(jj,i)*z_moment(jj,i)

      x_rjj=zero
      y_rjj=zero
      z_rjj=zero

      do 1012 jj=1,n
      x_rjj=x_rjj+coeff(jj,j)*x_moment(jj,j)
      y_rjj=y_rjj+coeff(jj,j)*y_moment(jj,j)
 1012 z_rjj=z_rjj+coeff(jj,j)*z_moment(jj,j)

      x_rij=zero
      y_rij=zero
      z_rij=zero

      do 1013 jj=1,n
      x_rij=x_rij+coeff(jj,j)*x_moment(jj,i)
      y_rij=y_rij+coeff(jj,j)*y_moment(jj,i)
 1013 z_rij=z_rij+coeff(jj,j)*z_moment(jj,i)

      aij=x_rij*x_rij
     &   +y_rij*y_rij
     &   +z_rij*z_rij
     &   -quarter*((x_rii-x_rjj)*(x_rii-x_rjj))
     &   -quarter*((y_rii-y_rjj)*(y_rii-y_rjj))
     &   -quarter*((z_rii-z_rjj)*(z_rii-z_rjj))

      bij=(x_rii-x_rjj)*x_rij
     &   +(y_rii-y_rjj)*y_rij
     &   +(z_rii-z_rjj)*z_rij

      factor=aij/(sqrt(aij**2+bij**2))

      if(abs(factor).gt.one) factor=one*factor/abs(factor)

      gamma=quarter*acos(-factor)

      if(bij.lt.zero) gamma=-gamma

      if(abs(gamma).gt.abs(quarter*pi)) stop 'boys localization'

      a=cos(gamma)
      b=sin(gamma)

      do 1014 ii=1,n
 1014 veci(ii)=coeff(ii,i)

      do 1015 ii=1,n
 1015 vecj(ii)=coeff(ii,j)

      do 1016 ii=1,n
 1016 coeff(ii,i)=a*veci(ii)+b*vecj(ii)

      do 1017 ii=1,n
 1017 coeff(ii,j)=a*vecj(ii)-b*veci(ii)

      do 1018 ii=1,n
 1018 veci(ii)=x_moment(ii,i)

      do 1019 ii=1,n
 1019 vecj(ii)=x_moment(ii,j)

      do 1020 ii=1,n
 1020 x_moment(ii,i)=a*veci(ii)+b*vecj(ii)

      do 1021 ii=1,n
 1021 x_moment(ii,j)=a*vecj(ii)-b*veci(ii)

      do 1022 ii=1,n
 1022 veci(ii)=y_moment(ii,i)

      do 1023 ii=1,n
 1023 vecj(ii)=y_moment(ii,j)

      do 1024 ii=1,n
 1024 y_moment(ii,i)=a*veci(ii)+b*vecj(ii)

      do 1025 ii=1,n
 1025 y_moment(ii,j)=a*vecj(ii)-b*veci(ii)

      do 1026 ii=1,n
 1026 veci(ii)=z_moment(ii,i)

      do 1027 ii=1,n
 1027 vecj(ii)=z_moment(ii,j)

      do 1028 ii=1,n
 1028 z_moment(ii,i)=a*veci(ii)+b*vecj(ii)

      do 1029 ii=1,n
 1029 z_moment(ii,j)=a*vecj(ii)-b*veci(ii)

 1010 continue

 1007 if(iloop.eq.nmaxloop) stop 'problem localizing orbitals'

 1030 continue

      write(6,1031)
 1031 format(//)

      do 1032 i=1,nocc

      x_rii=zero
      y_rii=zero
      z_rii=zero

      do 1033 jj=1,n
      x_rii=x_rii+coeff(jj,i)*x_moment(jj,i)
      y_rii=y_rii+coeff(jj,i)*y_moment(jj,i)
 1033 z_rii=z_rii+coeff(jj,i)*z_moment(jj,i)

      centroid(1,i)=x_rii
      centroid(2,i)=y_rii
      centroid(3,i)=z_rii

 1032 continue

      return
      end
