C
C $Id: linmin.F,v 1.6 1998/07/16 16:40:03 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      SUBROUTINE LINMIN(DIIS,DIRECT1,GRAD1,STEP,EHEAT1,
     &     ALPHA,IERROR)
C
C     CARRIES OUT A LINE MINIMIZATION IN THE DIRECTION OF THE VECTOR
C     DIRECT.  THE ENERGY FROM A FORWARD HALF STEP IS COMBINED WITH
C     THE REFERENCE GRADIENT TO CONSTRUCT A PARABOLA IN THE COORDINATE
C     ALPHA*DIRECT.  THE PREDICTED MINIMUM COORDINATE IS RETURNED
C     IN ALPHA, AND THE COORESPONDING ENERGY IN EHEAT1.
C
C     IF THE DIIS FLAG IS TURNED ON, THEN WE WILL RETURN AFTER THE
C     INITIAL STEP IF THE ENERGY DROPS.  IF THE ENERGY INCREASES,
C     HOWEVER, THEN THE FULL LINE MINIMIZATION WILL PROCEED, AND DIIS
C     WILL BE SET TO .FALSE. TO TELL THE CALLING ROUTINE NOT TO USE
C     THE DIIS PROCEDURE THIS CYCLE.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
#include "divcon.dim"
#include "divcon.h"
      LOGICAL DIIS, newsub
      DIMENSION DIRECT1(*),GRAD1(*)
C
C     LOCAL:
C
      DIMENSION COORD0(MAXPAR)
C
      IERROR = 0
      IFAIL = 0
C
C     STORE THE ORIGINAL ENERGY AND GEOMETRY, THEN SHIFT THE GEOMETRY
C     FORWARD BY HALF A STEP.
C
      DMAX = 0.0D0
      EHEAT0 = EHEAT1
      IF(XYZSPC)THEN
        IJ = 0
        DO 20 I=1,NATOMS
          DO 10 J=1,3
            IJ = IJ + 1
            COORD0(IJ) = XYZ(J,I)
            ABSD = ABS(DIRECT1(IJ))
            IF(ABSD.GT.DMAX) DMAX = ABSD
 10       CONTINUE
 20     CONTINUE
      ELSE
        DO 40 I=1,NPAR
          ITYP = IPAR(1,I)
          IATM = IPAR(2,I)
          COORD0(I) = ZMAT(ITYP,IATM)
          ABSD = ABS(DIRECT1(I))
          IF(ABSD.GT.DMAX) DMAX = ABSD
 40     CONTINUE
      ENDIF
      IF(DIIS)THEN
        STEP1 = 0.5D0
      ELSE
        STEP1 = 0.5D0*STEP
      ENDIF
      CALL GETCRD(COORD0,DIRECT1,STEP1,IERROR)
      IF(IERROR.NE.0) RETURN
C
C     GET THE ENERGY OF THIS NEW GEOMETRY.
C
      newsub = .true.
C     do setup for pbc's if they have been requested
      if(PBC)then
         call setbox(ierror)
         if(ierror.ne.0) return
      endif
      if (newsub) then
         CALL GENSUB(IERROR)
         IF(IERROR.NE.0) RETURN
      endif
      CALL ENERGY(newsub,EHEAT2,IERROR)
      IF(IERROR.NE.0) RETURN
C
C     IF DIIS IS BEING USED, RETURN IF THE ENERGY DROPPED.  OTHERWISE,
C     PROCEED WITH THE FULL LINE MINIMIZATION.
C
      IF(DIIS.AND.EHEAT2.LT.EHEAT0)THEN
        EHEAT1 = EHEAT2
        ALPHA = STEP1
        RETURN
      ELSE
        DIIS = .FALSE.
      ENDIF
C
C     DETERMINE THE PARAMETERS FOR THE QUADRATIC: E = A + B*X + C*X**2.
C     IF C IS NOT GREATER THAN ZERO, THEN THE CURVE WILL NOT BE
C     U-SHAPED AND WE MUST REJECT THIS FORMULA.
C
      B = 0.0D0
      DO 60 I=1,NPAR
        B = B + GRAD1(I)*DIRECT1(I)
 60   CONTINUE
      C = (EHEAT2 -  EHEAT0 - B*STEP1)/STEP1**2
      IF(C.GT.0.0D0)THEN
C
C       FIT THE ENERGIES TO THE CURVE:  E = A + B*X + C*X*X.  THE
C       MINIMUM OCCURS WHEN X = ALPHA = -B/(2*C)
C
        ALPHA = -0.5D0*B/C
C
C       PUT A BOUND ON ALPHA.
C
        ALPHMX = 0.05D0/DMAX
        IF(ABS(ALPHA).GT.ALPHMX) ALPHA = SIGN(ALPHMX,ALPHA)
      ELSE
C
C       THE ENERGY CURVE IS NOT STRICTLY U-SHAPED.  THIS MAY BE CAUSED
C       BY BEING AT THE ASYMPTOTIC END OF A MORSE-TYPE CURVE WHERE
C       THE CURVATURE BECOMES NEGATIVE, I.E. WHERE THE * IS BELOW.
C
C                        |
C                        |           <--DIRECT
C                        |.                  .
C                      E |.             *
C                        | .         .
C                        |  .      .
C                        |    . .
C                        |
C                        |____________________
C                                  X
C
C
C     A LARGER STEP IN THE ENERGY LOWERING DIRECTION MAY BE NECESSARY.  IF
C     THE FIRST STEP CAUSED THE ENERGY TO GO UP, HOWEVER, THEN SET ALPHA
C     EQUAL TO HALF OF THE FIRST STEP.
C
        IF(EHEAT2.LT.EHEAT0)THEN
          ALPHA = 5.0D0*STEP
        ELSE
          ALPHA = 0.5D0*STEP1
        ENDIF
      ENDIF
C
C     GET COORDINATES AND ENERGY CORRESPONDING TO THE NEW VALUE OF ALPHA.
C
      CALL GETCRD(COORD0,DIRECT1,ALPHA,IERROR)
      IF(IERROR.NE.0)THEN
C
C       ALPHA WAS SO BIG THAT IT BECAME IMPOSSIBLE TO COMPUTE CARTESIAN
C       COORDINATES FROM INTERNAL COORDINATES.  WILL STICK WITH THE FIRST
C       STEP IF IT LOWERED THE ENERGY.
C
        IFAIL = 1
        GO TO 1000
      ENDIF
C     do setup for pbc's if they have been requested
      if(PBC)then
         call setbox(ierror)
         if(ierror.ne.0) return
      endif
      if (newsub) then
         CALL GENSUB(IERROR)
         IF(IERROR.NE.0) RETURN
      endif
      CALL ENERGY(newsub,EHTEST,IERROR)
      IF(IERROR.NE.0)THEN
C
C       ALPHA WAS SO BIG THAT THE RESULTING GEOMETRY DID NOT YIELD AN SCF.
C       STICK WITH THE FIRST STEP IF IT LOWERED THE ENERGY.
C
        IFAIL = 1
        GO TO 1000
      ENDIF
      IF(EHTEST.GT.EHEAT0)THEN
C
C       THE COMPUTED ALPHA CAUSED THE ENERGY TO GO UP.  WILL HAVE TO
C       CHOOSE ALPHA EQUAL TO STEP1 OR SOMETHING SMALLER.
C
        IFAIL = 1
      ELSE
C
C       ACCEPT THIS ALPHA AND THE ENERGY.
C
        EHEAT1 = EHTEST
      ENDIF
 1000 IF(IFAIL.EQ.1)THEN
C
C       EITHER THE CURVE IS NOT U-SHAPED OR THE COMPUTED ALPHA CANNOT
C       BE USED.  CHOOSE ALPHA=STEP1 IF IT LOWERED THE ENERGY.  OTHERWISE,
C       MAKE ALPHA VERY SMALL.
C
        IERROR = 0
        IF(EHEAT2.LT.EHEAT0)THEN
          ALPHA = STEP1
        ELSE
          ALPHA = 0.1D0*STEP1
        ENDIF
C
C       GET ENERGY USING THE BEST STEPSIZE.
C
        CALL GETCRD(COORD0,DIRECT1,ALPHA,IERROR)
C     . do setup for pbc's if they have been requested
        if(PBC)then
           call setbox(ierror)
           if(ierror.ne.0) return
        endif
        if (newsub) then
           CALL GENSUB(IERROR)
           IF(IERROR.NE.0) RETURN
        endif
        CALL ENERGY(newsub,EHEAT1,IERROR)
      ENDIF
      RETURN
      END