      subroutine ddiis(iunit,niteration,ncontractions,dmix,emax,focka,
     &                 fockb,amat,bmat,wkmat1,wkmat2,wkmat3,wkmat4,
     &                 overlap,ematrix,wkvec)

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"

      parameter(maxdiis=200)

      dimension istat(mpi_status_size)

      dimension bmatrix(maxdiis,maxdiis),bbmatrix(maxdiis,maxdiis)

      dimension focka(*),fockb(*),amat(*),bmat(*)

      dimension wkmat1(ncontractions,*),wkmat2(ncontractions,*),
     &          wkmat3(ncontractions,*),wkmat4(ncontractions,*),
     &          overlap(ncontractions,*),ematrix(ncontractions,*),
     &          wkvec(*)

      data zero,half,one,two,startdiis/0.0d0,0.5d0,1.0d0,2.0d0,0.2d0/

      save ndiiss,ndiisw,bmatrix

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

      if(niteration.gt.1) goto 1001

      ndiiss=0

      do 1002 i=1,maxdiis
      bmatrix(i,1)=-one
 1002 bmatrix(1,i)=-one

      bmatrix(1,1)=zero

      return

 1001 continue

      if(ndiiss.eq.0) then

                rewind 45

                do 1003 i=1,ncontractions
 1003           if(myid.eq.0) call fastrd(45,overlap(1,i),ncontractions)

                call mpi_bcast(overlap,ncontractions**2,
     &                         mpi_real8,0,mpi_comm_world,ierr)

                ij=0
                do 1004 i=1,ncontractions
                do 1004 j=1,i
                ij=ij+1
                wkmat1(i,j)=focka(ij)
 1004           wkmat1(j,i)=focka(ij)

                ij=0
                do 1005 i=1,ncontractions
                do 1005 j=1,i
                ij=ij+1
                wkmat2(i,j)=half*amat(ij)
 1005           wkmat2(j,i)=half*amat(ij)

                do 1006 i=1,ncontractions
 1006           wkmat2(i,i)=wkmat2(i,i)*two

                call mxmy_mpi(ncontractions,wkmat2,wkmat1,wkmat3)
                call mxmy_mpi(ncontractions,overlap,wkmat3,ematrix)
                call mxmy_mpi(ncontractions,wkmat2,overlap,wkmat3)
                call mxmy_mpi(ncontractions,wkmat1,wkmat3,wkmat4)

                do 1007 j=1,ncontractions
                do 1007 i=1,ncontractions
 1007           ematrix(i,j)=wkmat4(i,j)-ematrix(i,j)

                rewind iunit

                if(myid.eq.0) call fastwr(iunit,
     &                                    wkmat1,ncontractions**2)

                ij=0
                do 1008 i=1,ncontractions
                do 1008 j=1,i
                ij=ij+1
                wkmat1(i,j)=fockb(ij)
 1008           wkmat1(j,i)=fockb(ij)

                ij=0
                do 1009 i=1,ncontractions
                do 1009 j=1,i
                ij=ij+1
                wkmat2(i,j)=half*bmat(ij)
 1009           wkmat2(j,i)=half*bmat(ij)

                do 1010 i=1,ncontractions
 1010           wkmat2(i,i)=wkmat2(i,i)*two

                call mxmy_mpi(ncontractions,wkmat2,wkmat1,wkmat3)
                call mxmy_mpi(ncontractions,overlap,wkmat3,wkmat4)

                do 1011 j=1,ncontractions
                do 1011 i=1,ncontractions
 1011           ematrix(i,j)=ematrix(i,j)-wkmat4(i,j)

                call mxmy_mpi(ncontractions,wkmat2,overlap,wkmat3)
                call mxmy_mpi(ncontractions,wkmat1,wkmat3,wkmat4)

                do 1012 j=1,ncontractions
                do 1012 i=1,ncontractions
 1012           ematrix(i,j)=ematrix(i,j)+wkmat4(i,j)

                rewind iunit+2

                if(myid.eq.0) call fastwr(iunit+2,
     &                                    wkmat1,ncontractions**2)

                rewind iunit+1

                if(myid.eq.0) call fastwr(iunit+1,
     &                                    ematrix,ncontractions**2)

                emax=zero

                do 1013 i=1,ncontractions
                do 1013 j=1,ncontractions
 1013           if(abs(ematrix(i,j)).gt.emax) emax=abs(ematrix(i,j))

                if(emax.gt.startdiis) return

                dmix=one

                ndiiss=1

                ndiisw=niteration

                bmatrix(2,2)=zero

                do 1014 j=1,ncontractions
                do 1014 i=1,ncontractions
 1014           bmatrix(2,2)=bmatrix(2,2)+ematrix(i,j)**2

                return

                      endif

      rewind 45

      do 1015 i=1,ncontractions
 1015 if(myid.eq.0) call fastrd(45,overlap(1,i),ncontractions)

      call mpi_bcast(overlap,ncontractions**2,
     &               mpi_real8,0,mpi_comm_world,ierr)

      ij=0
      do 1016 i=1,ncontractions
      do 1016 j=1,i
      ij=ij+1
      wkmat1(i,j)=focka(ij)
 1016 wkmat1(j,i)=focka(ij)

      ij=0
      do 1017 i=1,ncontractions
      do 1017 j=1,i
      ij=ij+1
      wkmat2(i,j)=half*amat(ij)
 1017 wkmat2(j,i)=half*amat(ij)

      do 1018 i=1,ncontractions
 1018 wkmat2(i,i)=wkmat2(i,i)*two

      call mxmy_mpi(ncontractions,wkmat2,wkmat1,wkmat3)
      call mxmy_mpi(ncontractions,overlap,wkmat3,ematrix)
      call mxmy_mpi(ncontractions,wkmat2,overlap,wkmat3)
      call mxmy_mpi(ncontractions,wkmat1,wkmat3,wkmat4)

      do 1019 j=1,ncontractions
      do 1019 i=1,ncontractions
 1019 ematrix(i,j)=wkmat4(i,j)-ematrix(i,j)

      if(myid.eq.0) endfile iunit
      if(myid.eq.0) backspace iunit

      if(myid.eq.0) call fastwr(iunit,wkmat1,ncontractions**2)

      ij=0
      do 1020 i=1,ncontractions
      do 1020 j=1,i
      ij=ij+1
      wkmat1(i,j)=fockb(ij)
 1020 wkmat1(j,i)=fockb(ij)

      ij=0
      do 1021 i=1,ncontractions
      do 1021 j=1,i
      ij=ij+1
      wkmat2(i,j)=half*bmat(ij)
 1021 wkmat2(j,i)=half*bmat(ij)

      do 1022 i=1,ncontractions
 1022 wkmat2(i,i)=wkmat2(i,i)*two

      call mxmy_mpi(ncontractions,wkmat2,wkmat1,wkmat3)
      call mxmy_mpi(ncontractions,overlap,wkmat3,wkmat4)

      do 1023 j=1,ncontractions
      do 1023 i=1,ncontractions
 1023 ematrix(i,j)=ematrix(i,j)-wkmat4(i,j)

      call mxmy_mpi(ncontractions,wkmat2,overlap,wkmat3)
      call mxmy_mpi(ncontractions,wkmat1,wkmat3,wkmat4)

      do 1024 j=1,ncontractions
      do 1024 i=1,ncontractions
 1024 ematrix(i,j)=ematrix(i,j)+wkmat4(i,j)

      nn_start=ndiisw
      nn_end=niteration

      if(myid.eq.0) then
                          endfile iunit+2
                          backspace iunit+2

                          call fastwr(iunit+2,wkmat1,ncontractions**2)

                          endfile iunit+1
                          backspace iunit+1

                          call fastwr(iunit+1,ematrix,ncontractions**2)

                          rewind iunit+1
                    endif

      do 1025 k=nn_start,nn_end

      if(myid.eq.0) call fastrd(iunit+1,wkmat1,ncontractions**2)

      call mpi_bcast(wkmat1,ncontractions**2,
     &               mpi_real8,0,mpi_comm_world,ierr)

      factor=zero

      do 1026 i=1,ncontractions
      do 1026 j=1,ncontractions
 1026 factor=factor+wkmat1(j,i)*ematrix(j,i)

      bmatrix(k-nn_start+2,niteration-nn_start+2)=factor
      bmatrix(niteration-nn_start+2,k-nn_start+2)=factor

 1025 continue

      nn_stop=niteration-ndiisw+2

      do 1027 i=1,nn_stop
      do 1027 j=1,nn_stop
 1027 bbmatrix(j,i)=bmatrix(j,i)

      do 1028 i=1,nn_stop
 1028 wkvec(i)=zero

      wkvec(1)=-one

      emax=zero

      do 1029 i=1,ncontractions
      do 1029 j=1,ncontractions
 1029 if(abs(ematrix(j,i)).gt.emax) emax=abs(ematrix(j,i))

      ndiisd=nn_stop

      if(nn_stop.gt.02) ndiisd=2
      if(nn_stop.gt.04) ndiisd=3
      if(nn_stop.gt.06) ndiisd=4
      if(nn_stop.gt.08) ndiisd=min(10,nn_stop-4)

      if(nn_stop.gt.ndiisd) then
                                  nnn=nn_stop-ndiisd

                                  do 1030 i=2,ndiisd
                                  do 1030 j=2,ndiisd
 1030                             bbmatrix(i,j)=bmatrix(i+nnn,j+nnn)
                            endif

      ndim=min(ndiisd,nn_stop)

      call matinv(bbmatrix,wkvec,ndim,2,maxdiis)

      if(myid.eq.0) rewind iunit
      if(myid.eq.0) rewind iunit+2

      if(nn_stop.gt.ndiisd) then

                    nnn=nn_stop-ndiisd

                    do 1031 k=1,nnn
 1031               if(myid.eq.0) call fastrd(iunit,wkmat1,1)

                    do 1032 k=1,nnn
 1032               if(myid.eq.0) call fastrd(iunit+2,wkmat1(1,i),1)

                            endif

      do 1033 i=1,ncontractions
      do 1033 j=1,ncontractions
 1033 wkmat1(j,i)=zero

      do 1034 i=1,ncontractions
      do 1034 j=1,ncontractions
 1034 wkmat2(j,i)=zero

      do 1035 k=2,ndim

      if(myid.eq.0) call fastrd(iunit,wkmat3,ncontractions**2)

      call mpi_bcast(wkmat3,ncontractions**2,
     &               mpi_real8,0,mpi_comm_world,ierr)

      if(myid.eq.0) call fastrd(iunit+2,wkmat4,ncontractions**2)

      call mpi_bcast(wkmat4,ncontractions**2,
     &               mpi_real8,0,mpi_comm_world,ierr)

      do 1036 i=1,ncontractions
      do 1036 j=1,ncontractions
 1036 wkmat1(j,i)=wkmat1(j,i)+wkvec(k)*wkmat3(j,i)

      do 1037 i=1,ncontractions
      do 1037 j=1,ncontractions
 1037 wkmat2(j,i)=wkmat2(j,i)+wkvec(k)*wkmat4(j,i)

 1035 continue

      ij=0

      do 1038 i=1,ncontractions
      do 1038 j=1,i
      ij=ij+1
 1038 focka(ij)=wkmat1(j,i)

      ij=0

      do 1039 i=1,ncontractions
      do 1039 j=1,i
      ij=ij+1
 1039 fockb(ij)=wkmat2(j,i)

      return
      end
