      subroutine annihilate(nmax,ncontractions,nocc,istart,iend,dmaxpc,
     &                      tol,dmat,fock,coeff,fockmat,veci,vecj,
     &                      wkmat1,wkmat2,wkvec,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(3,*),iend(3,*)

      dimension dmat(*),fock(*),coeff(nmax,*),fockmat(ncontractions,*),
     &          veci(*),vecj(*),wkmat1(ncontractions,*),
     &          wkmat2(ncontractions,*),wkvec(*),s(*),tau(*)

      data loopmax/500/

      data small,tol_min,tol_max/0.0001d0,0.0001d0,0.0000000001d0/
      data zero,half,one,two,four/0.0d0,0.5d0,1.0d0,2.0d0,4.0d0/

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

      tolerance=tol
      if(tolerance.gt.tol_min) tolerance=tol_min
      if(tolerance.lt.tol_max) tolerance=tol_max

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

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

      isent=nocc
      do 1003 islave=0,numprocs-1
      ia=ncontractions-isent
      ib=numprocs-islave
      nloop=ia/ib
      if(mod(ia,ib).ne.0) nloop=nloop+1
      istart(3,islave+1)=isent+1
      iend(3,islave+1)=isent+nloop
 1003 isent=isent+nloop

      do 1004 j=1,ncontractions
      do 1004 i=1,ncontractions
      ij=((max(i,j)*(max(i,j)-1))/2)+min(i,j)
 1004 fockmat(i,j)=fock(ij)

      do 1005 j=1,ncontractions
      do 1005 i=istart(1,myid+1),iend(1,myid+1)
      wkmat1(i,j)=zero
      do 1005 k=1,ncontractions
 1005 wkmat1(i,j)=wkmat1(i,j)+coeff(k,j)*fockmat(k,i)

      do 1006 j=1,ncontractions
      do 1006 i=1,ncontractions
      wkmat2(i,j)=zero
      do 1006 k=istart(1,myid+1),iend(1,myid+1)
 1006 wkmat2(i,j)=wkmat2(i,j)+coeff(k,i)*wkmat1(k,j)

      call mpi_reduce(wkmat2,fockmat,ncontractions**2,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
      call mpi_bcast(fockmat,ncontractions**2,
     &               mpi_real8,0,mpi_comm_world,ierr)

      do 1007 loop=1,loopmax

      f_ij_max=zero

      do 1008 i=1,nocc

      do 1009 j=nocc+1,ncontractions
      wkvec(j)=zero
      do 1009 k=istart(1,myid+1),iend(1,myid+1)
 1009 wkvec(j)=wkvec(j)+coeff(k,i)*wkmat1(k,j)

      call mpi_reduce(wkvec(nocc+1),wkmat2(nocc+1,1),ncontractions-nocc,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
      call mpi_bcast(wkmat2(nocc+1,1),ncontractions-nocc,
     &               mpi_real8,0,mpi_comm_world,ierr)

      do 1010 j=nocc+1,ncontractions
 1010 fockmat(i,j)=wkmat2(j,1)

      do 1011 j=nocc+1,ncontractions
      wkvec(j)=zero
      do 1011 k=istart(1,myid+1),iend(1,myid+1)
 1011 wkvec(j)=wkvec(j)+coeff(k,j)*wkmat1(k,j)

      call mpi_reduce(wkvec(nocc+1),wkmat2(nocc+1,1),ncontractions-nocc,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
      call mpi_bcast(wkmat2(nocc+1,1),ncontractions-nocc,
     &               mpi_real8,0,mpi_comm_world,ierr)

      do 1012 j=nocc+1,ncontractions
 1012 fockmat(j,j)=wkmat2(j,1)

      f_ii=zero
      do 1013 k=istart(1,myid+1),iend(1,myid+1)
 1013 f_ii=f_ii+coeff(k,i)*wkmat1(k,i)

      call mpi_reduce(f_ii,fockmat(i,i),1,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
      call mpi_bcast(fockmat(i,i),1,mpi_real8,0,mpi_comm_world,ierr)

      do 1014 j=istart(3,myid+1),iend(3,myid+1)

      f_ij=fockmat(i,j)
      f_ii=fockmat(i,i)
      f_jj=fockmat(j,j)

      if(abs(f_ij).gt.f_ij_max) f_ij_max=abs(f_ij)

      d=f_jj-f_ii

      if(abs(d).le.tol_min) then
                                  s(j)=zero
                                  tau(j)=zero

                                  goto 1014
                            endif

      c=sqrt(four*f_ij*f_ij+d*d)

      diff=abs(c-d)

      if(diff.gt.small) then
                              s(j)=sqrt(half*diff/c)
                        else
                              s(j)=sqrt((f_ij*f_ij)/(d*c))
                        endif

      if(f_ij.ge.zero) s(j)=-s(j)

      c=sqrt(max(zero,one-s(j)*s(j)))

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

 1014 continue

      do 1015 n=1,numprocs
      nloop=iend(3,n)-istart(3,n)+1
      call mpi_bcast(s(istart(3,n)),
     &               nloop,mpi_real8,n-1,mpi_comm_world,ierr)
      call mpi_bcast(tau(istart(3,n)),
     &               nloop,mpi_real8,n-1,mpi_comm_world,ierr)
 1015 continue

      do 1008 j=nocc+1,ncontractions

      do 1016 k=istart(1,myid+1),iend(1,myid+1)
 1016 veci(k)=coeff(k,i)

      do 1017 k=istart(1,myid+1),iend(1,myid+1)
 1017 vecj(k)=coeff(k,j)

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

      do 1019 k=istart(1,myid+1),iend(1,myid+1)
 1019 coeff(k,j)=coeff(k,j)-s(j)*(veci(k)+tau(j)*vecj(k))

      do 1020 k=istart(1,myid+1),iend(1,myid+1)
 1020 veci(k)=wkmat1(k,i)

      do 1021 k=istart(1,myid+1),iend(1,myid+1)
 1021 vecj(k)=wkmat1(k,j)

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

      do 1023 k=istart(1,myid+1),iend(1,myid+1)
 1023 wkmat1(k,j)=wkmat1(k,j)-s(j)*(veci(k)+tau(j)*vecj(k))

 1008 continue

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

      if(myid.eq.0) then

              f_ij_max_all=zero

              do 1024 i=0,numprocs-1
              call mpi_recv(f_ij_max,1,mpi_real8,i,
     &                      mpi_any_tag,mpi_comm_world,istat,ierror)
 1024         if(f_ij_max.gt.f_ij_max_all) f_ij_max_all=f_ij_max

                    endif

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

      if(f_ij_max_all.lt.tolerance) goto 1025

      if((loop.eq.loopmax).and.(myid.eq.0)) write(6,1026) f_ij_max
 1026 format(/,' annihilation routine did not fully converge,',/,
     &         ' maximum value = ', f17.13,/,
     &         ' (do not worry unless this is the final iteration)')

 1007 continue

 1025 continue

      do 1027 n=1,numprocs
      do 1027 j=1,ncontractions
      nloop=iend(1,n)-istart(1,n)+1
 1027 call mpi_bcast(coeff(istart(1,n),j),nloop,
     &               mpi_real8,n-1,mpi_comm_world,ierr)

      do 1028 j=1,ncontractions
      do 1028 i=1,ncontractions
      ij=((max(i,j)*(max(i,j)-1))/2)+min(i,j)
 1028 fockmat(i,j)=fock(ij)

      do 1029 j=1,nocc
      do 1029 i=istart(1,myid+1),iend(1,myid+1)
      wkmat1(i,j)=zero
      do 1029 k=1,ncontractions
 1029 wkmat1(i,j)=wkmat1(i,j)+coeff(k,j)*fockmat(k,i)

      do 1030 j=1,nocc
      do 1030 i=1,nocc
      fockmat(i,j)=zero
      do 1030 k=istart(1,myid+1),iend(1,myid+1)
 1030 fockmat(i,j)=fockmat(i,j)+coeff(k,i)*wkmat1(k,j)

      do 1031 i=1,nocc
      call mpi_reduce(fockmat(1,i),wkmat2(1,i),
     &                nocc,mpi_real8,mpi_sum,0,mpi_comm_world,ierr)
 1031 call mpi_bcast(wkmat2(1,i),nocc,mpi_real8,0,mpi_comm_world,ierr)

      do 1032 i=1,nocc
 1032 if(myid.eq.0) call fastwr(51,wkmat2(1,i),nocc)

      ncntsq=(ncontractions*(ncontractions+1))/2

      do 1033 i=1,ncntsq
 1033 wkvec(i)=zero

      do 1034 i=1,ncontractions
      do 1034 j=1,i
      imax=max(i,j)
      imin=min(i,j)
      ij=((imax*(imax-1))/2)+imin
      temp_value=zero
      do 1035 k=istart(2,myid+1),iend(2,myid+1)
 1035 temp_value=temp_value+coeff(i,k)*coeff(j,k)
      if(i.ne.j) temp_value=two*temp_value
 1034 wkvec(ij)=temp_value

      call mpi_reduce(wkvec,fockmat,ncntsq,
     &                mpi_real8,mpi_sum,0,mpi_comm_world,ierr)

      if(myid.eq.0) then

              call stringcopy(ncntsq,wkvec,fockmat)

              do 1036 i=1,ncntsq
              difference=abs(wkvec(i)-dmat(i))
              if(difference.gt.dmaxpc) dmaxpc=difference
 1036         dmat(i)=wkvec(i)

                    endif

      call mpi_bcast(dmat,ncntsq,mpi_real8,0,mpi_comm_world,ierr)
      call mpi_bcast(dmaxpc,1,mpi_real8,0,mpi_comm_world,ierr)

      return
      end
