C
C $Id: diisub.F,v 1.5 1998/07/16 16:39:37 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE DIISUB(DIIS,GRAD1,EHEAT1,DIRECT1,DNORM,GNORM1,ALPHA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
      LOGICAL DIIS
      DIMENSION GRAD1(*),DIRECT1(*)
C
C     LOCAL:
C
      PARAMETER (M1=MXSTOR+1)
      PARAMETER (M1SQR=M1*M1)
C
C     MXSTOR IS SET IN divcon.dim.  RECOMMENDED VALUE IS 20 -- I.E., WILL
C     STORE UP TO 20 SETS OF COORDINATES AND GRADIENTS.
C
      DIMENSION B(M1SQR),C(M1),RHS(M1),WTEMP(M1),CRDSET(MAXPAR,MXSTOR),
     .          GRDSET(MAXPAR,MXSTOR),ESET(MXSTOR),HGSET(MAXPAR,MXSTOR),
     .          DTEST(MAXPAR),GTEST(MAXPAR)
      DATA NSETS,ISTORE /0,0/
      SAVE NSETS,ISTORE
      SAVE CRDSET,GRDSET,ESET,HGSET
C
C
C     ASSIGN LIMITS AND POINTERS:
C
C         MAXSET = THE MAXIMUM NUMBER OF SETS OF COORDINATES, GRADIENTS,
C                  AND ENERGIES THAT WILL EVER BE USED TO GENERATE DIIS
C                  STEP.  SHOULD BE LESS THAN NPAR TO PREVENT GETTING A
C                  SINGULAR B MATRIX.
C
C         NSETS  = THE NUMBER OF SETS WE ARE USING FOR THE CURRENT ITERATION.
C
C         ISTORE = THE POSITION IN CRDSET, GRDSET AND ESET TO STORE THE
C                  LATEST INFORMATION.  THE MAXSET POSITIONS ARE FIRST
C                  FILLED UP, THEN CYCLED BACK THROUGH IF NECESSARY.
C
      MAXSET = MIN(MXSTOR,NPAR/2)
      NSETS = NSETS + 1
      NSETS = MIN(NSETS,MAXSET)
      ISTORE = ISTORE + 1
      IF(ISTORE.GT.MAXSET) ISTORE = ISTORE - MAXSET
C
C     DIIS WILL BE .FALSE. IF A FULL LINE MINIMIZATION HAS ALREADY BEEN
C     DONE.  IN THIS CASE, JUST RETURN AFTER STORING COORDINATES, GRADIENT
C     AND ENERGY.  SAME THING HOLDS ON THE FIRST CALL TO DIISUB.
C
      CALL SAVSET(ISTORE,ZMAT,GRAD1,EHEAT1,CRDSET,GRDSET,ESET)
      IF(.NOT.DIIS.OR.NSETS.EQ.1) RETURN
C
C     ASSIGN HINV*GRAD1(KSET) FOR THE NSETS WORTH OF GRADIENT VECTORS.
C
      DO 200 KSET=1,NSETS
        IJ = 0
        II = 0
        DO 180 I=1,NPAR
          HGI = 0.0D0
          DO 140 J=1,I
            IJ = IJ + 1
            HGI = HGI + HINV(IJ)*GRDSET(J,KSET)
 140      CONTINUE
          IF(I.LT.NPAR)THEN
            II = II + I
            JI = II + I
            DO 160 J=I+1,NPAR
              HGI = HGI + HINV(JI)*GRDSET(J,KSET)
              JI = JI + J
 160        CONTINUE
          ENDIF
          HGSET(I,KSET) = HGI
 180    CONTINUE
 200  CONTINUE
C
C     ASSIGN THE B MATRIX:
C
C          _                                                        _
C         |                                                          |
C         |  HG(1,1)      HG(1,2)     . . .     HG(1,NSETS)      -1  |
C         |  HG(2,1)      HG(2,2)     . . .     HG(2,NSETS)      -1  |
C         |  .            .                     .                 .  |
C     B = |  .            .                     .                 .  |
C         |  .            .                     .                 .  |
C         |  HG(NSETS,1)  HG(NSETS,2) . . .     HG(NSETS,NSETS)  -1  |
C         | -1          -1          . . .      -1                 0  |
C         |_                                                        _|
C
C
C     WHERE HG(ISET,JSET) = (HINV*GRAD1(ISET))'*(HINV*GRAD1(JSET))
C
      NDIM = NSETS + 1
      JI0 = 0
      DO 300 ISET=1,NSETS
        IJ0 = 0
        DO 280 JSET=1,ISET
          IJ = IJ0 + ISET
          JI = JI0 + JSET
          BIJ = 0.0D0
          DO 260 K=1,NPAR
            BIJ = BIJ + HGSET(K,ISET)*HGSET(K,JSET)
 260      CONTINUE
          B(IJ) = BIJ
          B(JI) = BIJ
          IJ0 = IJ0 + NDIM
 280    CONTINUE
        JI0 = JI0 + NDIM
 300  CONTINUE
      IJ = NDIM
      JI = NSETS*NDIM+1
      DO 320 JSET=1,NSETS
        B(IJ) = -1.0D0
        B(JI) = -1.0D0
        RHS(JSET) = 0.0D0
        IJ = IJ + NDIM
        JI = JI + 1
 320  CONTINUE
      B(IJ+NDIM) = 0.0D0
      RHS(NDIM) = -1.0D0
C
C     GET THE VECTOR C BY SOLVING THE LINEAR SYSTEM B*C = RHS.
C
      THRESH = 1.0D-12
      IERROR = 0
      CALL LSOLVE(NDIM,B,RHS,WTEMP,THRESH,C,IERROR)
C
C     DON'T BOTHER WITH THE DIIS STEP IF B WAS SINGULAR.  NOT A FATAL
C     ERROR, THOUGH, BECAUSE WE CAN JUST KEEP THE RESULTS OF THE
C     ORIGINAL STEP.
C
      IF(IERROR.NE.0) RETURN
C
C     GET THE INTERPOLATED STEP (=DIRECT), GRADIENT, AND ENERGY.
C
      DNTEST = 0.0D0
      GNTEST = 0.0D0
      DO 400 I=1,NPAR
C
C       INTERPOLATED DIRECT IS JUST THE INTERPOLATED COORDINATES MINUS
C       THE ORIGINAL COORDINATES.
C
        DI = -CRDSET(I,ISTORE)
        GI = 0.0D0
        DO 380 KSET=1,NSETS
          DI = DI + C(KSET)*CRDSET(I,KSET)
          GI = GI + C(KSET)*GRDSET(I,KSET)
 380    CONTINUE
        DTEST(I) = DI
        GTEST(I) = GI
        DNTEST = DNTEST + DI*DI
        GNTEST = GNTEST + GI*GI
 400  CONTINUE
      DNTEST = DSQRT(DNTEST)
      GNTEST = DSQRT(GNTEST)
      ETEST = 0.0D0
      DO 420 KSET=1,NSETS
        ETEST = ETEST + C(KSET)*ESET(KSET)
 420  CONTINUE
C
C     USE INTERPOLATED DIRECT VECTOR WITH UNIT STEP TO GET INTERPOLATED
C     COORDINATES, BUT ONLY IF THE INTERPOLATED ENERGY AND GRADIENT
C     NORM ARE LOWER THAN THE ORIGINAL ONES.
C
      IF(ETEST.LT.EHEAT1.AND.GNTEST.LT.GNORM1)THEN
        DNORM = DNTEST
        GNORM1 = GNTEST
        DO 440 I=1,NPAR
          DIRECT1(I) = DTEST(I)
          GRAD1(I) = GTEST(I)
 440    CONTINUE
        ALPHA = 1.0D0
        CALL GETCRD(CRDSET(1,ISTORE),DIRECT1,ALPHA,IERROR)
      ENDIF
      RETURN
      END
C
C
C
      SUBROUTINE SAVSET(ISTORE,ZMAT1,GRAD1,EHEAT1,CRDSET,GRDSET,ESET)
C
C     STORES COORDINATES (ZMAT1), GRADIENT (GRAD1), AND ENERGY (EHEAT1)
C     IN POSTION ISTORE OF CRDSET, GRDSET, AND ESET, RESPECTIVELY.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
      DIMENSION ZMAT1(3,*),GRAD1(*),CRDSET(MAXPAR,*),GRDSET(MAXPAR,*),
     .          ESET(*)
C
      DO 20 I=1,NPAR
        ITYP = IPAR(1,I)
        IATM = IPAR(2,I)
        CRDSET(I,ISTORE) = ZMAT1(ITYP,IATM)
        GRDSET(I,ISTORE) = GRAD1(I)
 20   CONTINUE
      ESET(ISTORE) = EHEAT1
      RETURN
      END

CC---------------------------------------------------------------------CC