      subroutine s_annihilate(nmax,n,istart,iend,overlap,coeffa,coeffb,
     &                        smat,wkmat1,wkmat2,veci,vecj,s,tau)

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)

      include "mpif.h"

      dimension istat(mpi_status_size)

      dimension istart(*),iend(*)

      dimension overlap(nmax,*),coeffa(nmax,*),coeffb(nmax,*),
     &          smat(n,*),wkmat1(n,*),wkmat2(n,*),veci(*),vecj(*),s(*),
     &          tau(*)

      data loopmax/100/

      data big,tol/10000.0d0,0.0000000001d0/
      data zero,half,one,two,hundred/0.0d0,0.5d0,1.0d0,2.0d0,100.0d0/

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      isent=0
      do 1001 islave=0,numprocs-1
      ia=n-isent
      ib=numprocs-islave
      nloop=ia/ib
      if(mod(ia,ib).ne.0) nloop=nloop+1
      istart(islave+1)=isent+1
      iend(islave+1)=isent+nloop
 1001 isent=isent+nloop

      do 1002 j=1,n
      do 1003 i=istart(myid+1),iend(myid+1)
 1003 coeffa(i,j)=zero
 1002 coeffa(j,j)=one

      do 1004 j=1,n
      do 1004 i=istart(myid+1),iend(myid+1)
 1004 smat(i,j)=overlap(i,j)

      do 1005 j=1,n
      do 1005 i=1,n
      wkmat1(i,j)=zero
      do 1005 k=istart(myid+1),iend(myid+1)
 1005 wkmat1(i,j)=wkmat1(i,j)+coeffa(k,j)*smat(k,i)

      call mpi_reduce(wkmat1,wkmat2,
     &                n*n,mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
      call mpi_bcast(wkmat2,n*n,mpi_real8,0,mpi_comm_world,ierr)

      prev_s_ij_max=two

      do 1006 loop=1,loopmax

      s_ij_max=zero

      do 1007 j=1,n

      do 1008 i=1,n
      wkmat1(i,j)=zero
      do 1008 k=istart(myid+1),iend(myid+1)
 1008 wkmat1(i,j)=wkmat1(i,j)+coeffa(k,j)*smat(k,i)

      call mpi_reduce(wkmat1(1,j),wkmat2(1,j),
     &                n,mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
      call mpi_bcast(wkmat2(1,j),n,mpi_real8,0,mpi_comm_world,ierr)

      do 1009 i=istart(myid+1),iend(myid+1)

      if(i.eq.j) then
                       s(i)=zero
                       tau(i)=zero

                       goto 1009
                 endif

      s_ij=wkmat2(i,j)

      if(abs(s_ij).gt.s_ij_max) s_ij_max=abs(s_ij)

      t01=tol/hundred

      if(abs(s_ij).gt.t01) then
                                 s_ii=wkmat2(i,i)
                                 s_jj=wkmat2(j,j)

                                 theta=(s_jj-s_ii)/(two*s_ij)
                           else
                                 s(i)=zero
                                 tau(i)=zero

                                 goto 1009
                           endif

      if(theta.lt.zero) then
                              tan_phi=theta+sqrt(theta*theta+one)
                        else
                              tan_phi=theta-sqrt(theta*theta+one)
                        endif

      if(abs(theta).gt.big) tan_phi=-(half/theta)

      phi=atan(tan_phi)

      if(prev_s_ij_max.gt.(one/hundred)) then
                                               c=cos(half*phi)
                                               s(i)=sin(half*phi)
                                         else
                                               c=cos(phi)
                                               s(i)=sin(phi)
                                         endif

      tau(i)=s(i)/(one+c)

 1009 continue

      do 1010 l=1,numprocs
      m=iend(l)-istart(l)+1
      call mpi_bcast(s(istart(l)),m,mpi_real8,l-1,mpi_comm_world,ierr)
 1010 call mpi_bcast(tau(istart(l)),m,mpi_real8,l-1,mpi_comm_world,ierr)

      cutoff=prev_s_ij_max/big

      do 1007 i=1,n

      if(abs(s(i)).lt.cutoff) goto 1007

      do 1011 k=istart(myid+1),iend(myid+1)
 1011 veci(k)=coeffa(k,i)

      do 1012 k=istart(myid+1),iend(myid+1)
 1012 vecj(k)=coeffa(k,j)

      do 1013 k=istart(myid+1),iend(myid+1)
 1013 coeffa(k,i)=coeffa(k,i)+s(i)*(vecj(k)-tau(i)*veci(k))

      do 1014 k=istart(myid+1),iend(myid+1)
 1014 coeffa(k,j)=coeffa(k,j)-s(i)*(veci(k)+tau(i)*vecj(k))

      do 1015 k=istart(myid+1),iend(myid+1)
 1015 veci(k)=smat(k,i)

      do 1016 k=istart(myid+1),iend(myid+1)
 1016 vecj(k)=smat(k,j)

      do 1017 k=istart(myid+1),iend(myid+1)
 1017 smat(k,i)=smat(k,i)+s(i)*(vecj(k)-tau(i)*veci(k))

      do 1018 k=istart(myid+1),iend(myid+1)
 1018 smat(k,j)=smat(k,j)-s(i)*(veci(k)+tau(i)*vecj(k))

 1007 continue

      call mpi_bsend(s_ij_max,1,mpi_real8,0,0,mpi_comm_world,ierr)

      if(myid.eq.0) then

              s_ij_max_all=zero

              do 1019 i=0,numprocs-1
              call mpi_recv(s_ij_max,1,mpi_real8,i,
     &                      mpi_any_tag,mpi_comm_world,istat,ierror)
 1019         if(s_ij_max.gt.s_ij_max_all) s_ij_max_all=s_ij_max

                    endif

      call mpi_bcast(s_ij_max_all,1,mpi_real8,0,mpi_comm_world,ierr)

      if(s_ij_max_all.lt.tol) goto 1020

      if((loop.eq.loopmax).and.(myid.eq.0)) write(6,1021)
 1021 format(/,' s_annihilation routine did not fully converge')

      if(loop.eq.loopmax) stop

 1006 prev_s_ij_max=s_ij_max_all

 1020 continue

      do 1022 j=1,n
      do 1022 i=1,n
      wkmat1(i,j)=zero
      do 1022 k=istart(myid+1),iend(myid+1)
 1022 wkmat1(i,j)=wkmat1(i,j)+coeffa(k,j)*smat(k,i)

      call mpi_reduce(wkmat1,wkmat2,
     &                n*n,mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
      call mpi_bcast(wkmat2,n*n,mpi_real8,0,mpi_comm_world,ierr)

      do 1023 i=1,n

      factor=sqrt(one/wkmat2(i,i))

      do 1024 j=istart(myid+1),iend(myid+1)
 1024 coeffa(j,i)=coeffa(j,i)*factor

      do 1025 j=istart(myid+1),iend(myid+1)
 1025 coeffb(j,i)=coeffa(j,i)

 1023 continue

      do 1026 k=1,numprocs
      do 1026 j=1,n
      nloop=iend(k)-istart(k)+1
      call mpi_bcast(coeffa(istart(k),j),iend(k)-istart(k)+1,
     &               mpi_real8,k-1,mpi_comm_world,ierr)
 1026 call mpi_bcast(coeffb(istart(k),j),iend(k)-istart(k)+1,
     &               mpi_real8,k-1,mpi_comm_world,ierr)

      rewind 44

      do 1027 i=1,n
 1027 if(myid.eq.0) call fastwr(44,coeffa(1,i),n)

      return
      end
