      SUBROUTINE DIIS(XP, XPARAM, GP, GRAD, HP, HEAT, HS, NVAR, FRST)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension XP(NVAR), XPARAM(NVAR), GP(NVAR),
     1          GRAD(NVAR), HS(NVAR*NVAR)
      logical FRST
************************************************************************
*                                                                      *
*     DIIS PERFORMS DIRECT INVERSION IN THE ITERATIVE SUBSPACE         *
*                                                                      *
*     THIS INVOLVES SOLVING FOR C IN XPARAM(NEW) = XPARAM' - HG'       *
*                                                                      *
*  WHERE XPARAM' = SUM(C(I)XPARAM(I), THE C COEFFICIENTES COMING FROM  *
*                                                                      *
*                   | B   1 | . | C | = | 0 |                          *
*                   | 1   0 |   |-L |   | 1 |                          *
*                                                                      *
* WHERE B(I,J) =GRAD(I)H(T)HGRAD(J)  GRAD(I) = GRADIENT ON CYCLE I     *
*                              H    = INVERSE HESSIAN                  *
*                                                                      *
*                          REFERENCE                                   *
*                                                                      *
*  P. CSASZAR, P. PULAY, J. MOL. STRUCT. (THEOCHEM), 114, 31 (1984)    *
*                                                                      *
************************************************************************
************************************************************************
*                                                                      *
*     GEOMETRY OPTIMIZATION USING THE METHOD OF DIRECT INVERSION IN    *
*     THE ITERATIVE SUBSPACE (GDIIS), COMBINED WITH THE BFGS OPTIMIZER *
*     (A VARIABLE METRIC METHOD)                                       *
*                                                                      *
*     WRITTEN BY PETER L. CUMMINS, UNIVERSITY OF SYDNEY, AUSTRALIA     *
*                                                                      *
*                              REFERENCE                               *
*                                                                      *
*      "COMPUTATIONAL STRATEGIES FOR THE OPTIMIZATION OF EQUILIBRIUM   *
*     GEOMETRIES AND TRANSITION-STATE STRUCTURES AT THE SEMIEMPIRICAL  *
*     LEVEL", PETER L. CUMMINS, JILL E. GREADY, J. COMP. CHEM., 10,    *
*     939-950 (1989).                                                  *
*                                                                      *
*     MODIFIED BY JJPS TO CONFORM TO EXISTING MOPAC CONVENTIONS        *
*                                                                      *
************************************************************************
      common /KEYWRD/ KEYWRD
      PARAMETER (MRESET=15, M2=(MRESET+1)*(MRESET+1))
      dimension  XSET(MRESET*MAXPAR),GSET(MRESET*MAXPAR), ESET(MRESET)
      dimension DX(MAXPAR),GSAVE(MAXPAR),
     1 ERR(MRESET*MAXPAR),B(M2),BS(M2),BST(M2)
      logical DEBUG, PRINT
      character*241 KEYWRD
      DEBUG=.FALSE.
      PRINT=(INDEX(KEYWRD,' DIIS').NE.0)
      if (PRINT) DEBUG=(INDEX(KEYWRD,'DEBUG').NE.0)
      if (PRINT)  write(MIOB,'(/,''      ***** BEGIN GDIIS ***** '')')
C
C  SPACE SIMPLY LOADS THE CURRENT VALUES OF XPARAM AND GNORM INTO
C  THE ARRAYS XSET AND GSET
C
      call SPACE(MRESET,MSET,XPARAM, GRAD, HEAT, NVAR, XSET, GSET, ESET
     1, FRST)
C
C     INITIALIZE SOME VARIABLES AND CONSTANTS
C
      NDIIS = MSET
      MPLUS = MSET + 1
      MM = MPLUS * MPLUS
C
C     COMPUTE THE APPROXIMATE ERROR VECTORS
C
      INV=-NVAR
      do 30 I=1,MSET
         INV = INV + NVAR
         do 30 J=1,NVAR
            S = 0.D0
            KJ=(J*(J-1))/2
            do 10 K=1,J
               KJ = KJ+1
   10       S = S - HS(KJ) * GSET(INV+K)
            do 20 K=J+1,NVAR
               KJ = (K*(K-1))/2+J
   20       S = S - HS(KJ) * GSET(INV+K)
   30 ERR(INV+J) = S
C
C     CONSTRUCT THE GDIIS MATRIX
C
      do 40 I=1,MM
   40 B(I) = 1.D0
      JJ=0
      INV=-NVAR
      do 50 I=1,MSET
         INV=INV+NVAR
         JNV=-NVAR
         do 50 J=1,MSET
            JNV=JNV+NVAR
            JJ = JJ + 1
            B(JJ)=0.D0
            do 50 K=1,NVAR
   50 B(JJ) = B(JJ) + ERR(INV+K) * ERR(JNV+K)
C
      do 60 I=MSET-1,1,-1
         do 60 J=MSET,1,-1
   60 B(I*MSET+J+I) = B(I*MSET+J)
      do 70 I=1,MPLUS
         B(MPLUS*I) = 1.D0
   70 B(MPLUS*MSET+I) = 1.D0
      B(MM) = 0.D0
C
C     ELIMINATE ERROR VECTORS WITH THE LARGEST NORM
C
   80 continue
      do 90 I=1,MM
   90 BS(I) = B(I)
      if (NDIIS .EQ. MSET) go to 140
      do 130 II=1,MSET-NDIIS
         XMAX = -1.D10
         ITERA = 0
         do 110 I=1,MSET
            XNORM = 0.D0
            INV = (I-1) * MPLUS
            do 100 J=1,MSET
  100       XNORM = XNORM + ABS(B(INV + J))
            if (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) then
               XMAX = XNORM
               ITERA = I
               IONE = INV + I
            end if
  110    continue
         do 120 I=1,MPLUS
            INV = (I-1) * MPLUS
            do 120 J=1,MPLUS
               JNV = (J-1) * MPLUS
               if (J.EQ.ITERA) B(INV + J) = 0.D0
               B(JNV + I) = B(INV + J)
  120    continue
         B(IONE) = 1.0D0
  130 continue
  140 continue
C
      if (DEBUG) then
C
C     OUTPUT THE GDIIS MATRIX
C
         WRITE(*,'(/5X,'' GDIIS MATRIX'')')
         IJ = 0
         do 150 I=1,MPLUS
            do 150 J=1,I
               IJ = IJ + 1
  150    BST(IJ) = B( MPLUS * (J-1) + I)
         call VECPRT(BST,MPLUS)
      end if
C
C     SCALE DIIS MATRIX BEFORE INVERSION
C
      do 160 I=1,MPLUS
         II = MPLUS * (I-1) + I
  160 GSAVE(I) = 1.D0 / DSQRT(1.D-20+DABS(B(II)))
      GSAVE(MPLUS) = 1.D0
      do 170 I=1,MPLUS
         do 170 J=1,MPLUS
            IJ = MPLUS * (I-1) + J
  170 B(IJ) = B(IJ) * GSAVE(I) * GSAVE(J)
C
      if (DEBUG) then
C
C     OUTPUT SCALED GDIIS MATRIX
C
         WRITE(*,'(/5X,'' GDIIS MATRIX (SCALED)'')')
         IJ = 0
         do 180 I=1,MPLUS
            do 180 J=1,I
               IJ = IJ + 1
  180    BST(IJ) = B( MPLUS * (J-1) + I)
         call VECPRT(BST,MPLUS)
      end if
C
C     INVERT THE GDIIS MATRIX
C
      call MINV(B,MPLUS,DET)
C
      do 190 I=1,MPLUS
         do 190 J=1,MPLUS
            IJ = MPLUS * (I-1) + J
  190 B(IJ) = B(IJ) * GSAVE(I) * GSAVE(J)
C
C     COMPUTE THE INTERMEDIATE INTERPOLATED PARAMETER AND GRADIENT
C     VECTORS
C
      do 200 K=1,NVAR
         XP(K) = 0.D0
         GP(K) = 0.D0
         do 200 I=1,MSET
            INK = (I-1) * NVAR + K
            XP(K) = XP(K) + B(MPLUS*MSET+I) * XSET(INK)
  200 GP(K) = GP(K) + B(MPLUS*MSET+I) * GSET(INK)
      HP=0.D0
      do 210 I=1,MSET
  210 HP=HP+B(MPLUS*MSET+I)*ESET(I)
C
      do 220 K=1,NVAR
  220 DX(K) = XPARAM(K) - XP(K)
      XNORM = SQRT(DOT(DX,DX,NVAR))
      if (PRINT) then
         write(MIOB,'(/10X,''DEVIATION IN X '',F7.4,8X,''DETERMINANT '',
     %                G9.3)') XNORM,DET
         write(MIOB,'(10X,''GDIIS COEFFICIENTS'')')
         write(MIOB,'(10X,5F12.5)') (B(MPLUS*MSET+I),I=1,MSET)
      end if
C
C     THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY!
C
      THRES = MAX(10.D0**(-NVAR), 1.D-25)
      if (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) then
         if (PRINT) write(MIOB,'(10X,
     1      ''THE DIIS MATRIX IS ILL CONDITIONED''
     2 /10X,'' - PROBABLY, VECTORS ARE LINEARLY DEPENDENT - '',
     3 /10X,''THE DIIS STEP WILL BE REPEATED WITH A SMALLER SPACE'')')
         do 230 K=1,MM
  230    B(K) = BS(K)
         NDIIS = NDIIS - 1
         if (NDIIS .GT. 0) go to 80
         if (PRINT) WRITE(*,'(10X,''NEWTON-RAPHSON STEP TAKEN'')')
         do 240 K=1,NVAR
            XP(K) = XPARAM(K)
  240    GP(K) = GRAD(K)
C
      end if
      if (PRINT)  write(MIOB,'(/,''      ***** end GDIIS ***** '',/)')
C
      return
      end
