      SUBROUTINE FLEPO (XPARAM,NVAR,FUNCT1)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension XPARAM(*)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /KEYWRD/ KEYWRD
      common /NUMSCF/ NSCF
      common /LAST  / LAST
      common /GRAVEC/ COSINE
      common /PATH  / LATOM,LPARAM,REACT(200)
      common /GRADNT/ GRAD(MAXPAR),GNORM
      common /MESAGE/ IFLEPO,ISCF
      common /TIME  / TIME0
      common /FMATRX/ HESINV(MAXHES)
      common /SCFTYP/ EMIN, LIMSCF
      common /TIMDMP/ TLEFT, TDUMP
      common /NUMCAL/ NUMCAL
      character*241 KEYWRD
C
C
C     *
C     THIS SUBROUTINE ATTEMPTS TO MINIMIZE A REAL-VALUED FUNCTION OF
C     THE N-COMPONENT REAL VECTOR XPARAM ACCORDING TO THE
C     BFGS FORMULA. RELEVANT REFERENCES ARE
C
C     BROYDEN, C.G., JOURNAL OF THE INSTITUTE FOR MATHEMATICS AND
C                     APPLICATIONS, VOL. 6 PP 222-231, 1970.
C     FLETCHER, R., COMPUTER JOURNAL, VOL. 13, PP 317-322, 1970.
C
C     GOLDFARB, D. MATHEMATICS OF COMPUTATION, VOL. 24, PP 23-26, 1970.
C
C     SHANNO, D.F. MATHEMATICS OF COMPUTATION, VOL. 24, PP 647-656
C                    1970.
C
C   SEE ALSO SUMMARY IN
C
C    HEAD, J.D.; AND ZERNER, M.C., CHEMICAL PHYSICS LETTERS, VOL. 122,
C          264 (1985).
C    SHANNO, D.F., J. OF OPTIMIZATION THEORY AND APPLICATIONS
C          VOL.46, NO 1 PP 87-94 1985.
C     *
C     THE FUNCTION CAN ALSO BE MINIMIZED USING THE
C     DAVIDON-FLETCHER-POWELL ALGORITHM (COMPUTER JOURNAL, VOL. 6,
C     P. 163).
C
C     THE USER MUST SUPPLY THE SUBROUTINE
C     COMPFG(XPARAM,.TRUE.,FUNCT,.TRUE.,GRAD,LGRAD)
C     WHICH COMPUTES FUNCTION VALUES  FUNCT AT GIVEN VALUES FOR THE
C     VARIABLES XPARAM, AND THE GRADIENT GRAD IF LGRAD=.TRUE.
C     THE MINIMIZATION PROCEEDS BY A SEQUENCE OF ONE-dimensionAL
C     MINIMIZATIONS.  THESE ARE CARRIED OUT WITHOUT GRADIENT COMPUTATION
C     BY THE SUBROUTINE LINMIN, WHICH SOLVES THE SUBPROBLEM OF
C     MINIMIZING THE FUNCTION FUNCT ALONG THE LINE XPARAM+ALPHA*PVECT,
C     WHERE XPARAM
C     IS THE VECTOR OF CURRENT VARIABLE VALUES,  ALPHA IS A SCALAR
C     VARIABLE, AND  PVECT  IS A SEARCH-DIRECTION VECTOR PROVIDED BY THE
C     BFGS OR DAVIDON-FLETCHER-POWELL ALGORITHM.  EACH ITERATION STEP CA
C     OUT BY FLEPO PROCEEDS BY LETTING LINMIN FIND A VALUE FOR ALPHA
C     WHICH MINIMIZES  FUNCT  ALONG  XPARAM+ALPHA*PVECT, BY
C     UPDATING THE VECTOR  XPARAM  BY THE AMOUNT ALPHA*PVECT, AND
C     FINALLY BY GENERATING A NEW VECTOR  PVECT.  UNDER
C     CERTAIN RESTRICTIONS (POWELL, J.INST.MATHS.APPLICS.(1971),
C     V.7,21-36)  A SEQUENCE OF FUNCT VALUES CONVERGING TO SOME
C     LOCAL MINIMUM VALUE AND A SEQUENCE OF
C     XPARAM VECTORS CONVERGING TO THE CORRESPONDING MINIMUM POINT
C     ARE PRODUCED.
C                          CONVERGENCE TESTS.
C
C     HERBERTS TEST: THE ESTIMATED DISTANCE FROM THE CURRENT POINT
C                    POINT TO THE MINIMUM IS LESS THAN TOLERA.
C
C                    "HERBERTS TEST SATISFIED - GEOMETRY OPTIMIZED"
C
C     GRADIENT TEST: THE GRADIENT NORM HAS BECOME LESS THAN TOLERG
C                    TIMES THE SQUARE ROOT OF THE NUMBER OF VARIABLES.
C
C                    "TEST ON GRADIENT SATISFIED".
C
C     XPARAM TEST:  THE RELATIVE CHANGE IN XPARAM, MEASURED BY ITS NORM,
C                   OVER ANY TWO SUCCESSIVE ITERATION STEPS DROPS BELOW
C                   TOLERX.
C
C                    "TEST ON XPARAM SATISFIED".
C
C     FUNCTION TEST: THE CALCULATED VALUE OF THE HEAT OF FORMATION
C                    BETWEEN ANY TWO CYCLES IS WITHIN TOLERF OF
C                    EACH OTHER.
C
C                    "HEAT OF FORMATION TEST SATISFIED"
C
C     FOR THE GRADIENT, FUNCTION, AND XPARAM TESTS A FURTHER CONDITION,
C     THAT NO INDIVIDUAL COMPONENT OF THE GRADIENT IS GREATER
C     THAN TOLERG, MUST BE SATISFIED, IN WHICH CASE THE
C     CALCULATION EXITS WITH THE MESSAGE
C
C                     "PETERS TEST SATISFIED"
C
C     WILL BE PRINTED, AND FUNCT AND XPARAM WILL CONTAIN THE LAST
C     FUNCTION VALUE CUM VARIABLE VALUES REACHED.
C
C
C     THE BROYDEN-FLETCHER-GOLDFARB-SHANNO AND DAVIDON-FLETCHER-POWELL
C     ALGORITHMS CHOOSE SEARCH DIRECTIONS
C     ON THE BASIS OF LOCAL PROPERTIES OF THE FUNCTION.  A MATRIX  H,
C     WHICH IN FLEPO IS PRESET WITH THE IDENTITY, IS MAINTAINED AND
C     UPDATED AT EACH ITERATION STEP.  THE MATRIX DESCRIBES A LOCAL
C     METRIC ON THE SURFACE OF FUNCTION VALUES ABOVE THE POINT XPARAM.
C     THE SEARCH-DIRECTION VECTOR  PVECT  IS SIMPLY A TRANSFORMATION
C     OF THE GRADIENT  GRAD  BY THE MATRIX H.
C
      dimension XVAR(MAXPAR), GVAR(MAXPAR), XD(MAXPAR), GD(MAXPAR),
     1GLAST(MAXPAR), XLAST(MAXPAR), GG(MAXPAR), PVECT(MAXPAR)
      dimension MDFP(9),XDFP(9), XTEMP(MAXPAR), GTEMP(MAXPAR)
      logical OKF, PRINT,  TIME, RESTRT, MINPRT, SADDLE, GEOOK, LOG
     1        ,RESET, RESFIL, LGRAD, DFP, LDIIS, THIEL, DIISOK, FRST,
     2         LIMSCF
      save  ICALCN
      save  RST, TDEL, SFACT, DELL, EINC, IGG1, DEL
      save  RESTRT, GEOOK, DFP, CONST
      save  SADDLE, MINPRT, ROOTV, PRINT, DELHOF, TOLERF, TOLERG
      save  TOLERX, DROP, FREPF, IHDIM, CNCADD, ABSMIN, ITRY1, ITRY2
      save  OKF, JCYC, LNSTOP, IREPET, ALPHA, PNORM, JNRST, CYCMX
      save  COS, NCOUNT, RESFIL, MDFP, TX1, TX2, TLAST
      save  TIME, TOTIME
      equivalence (MDFP(1),JCYC  ),(MDFP(2),JNRST),(MDFP(3),NCOUNT),
     1            (MDFP(4),LNSTOP),(XDFP(1),ALPHA),(XDFP(2),COS   ),
     2            (XDFP(3),PNORM ),(XDFP(4),DROP ),(XDFP(5),DEL   ),
     3            (XDFP(6),FREPF ),(XDFP(7),CYCMX),(XDFP(8),TOTIME)
      data ICALCN /0/
C
C   START OF ONCE-ONLY SECTION
C
      EMIN=0.D0
      if (ICALCN.NE.NUMCAL) then
C
C   THE FOLLOWING CONSTANTS SHOULD BE SET BY THE USER.
C
         RST   = 0.05D0
         IPRT  = MIOB
         TDEL  = 0.06D0
         NRST  = 30
         SFACT = 1.5
         DELL  = 0.01D0
         EINC  = 0.3D0
         IGG1  = 3
         DEL=DELL
C
C    THESE CONSTANTS SHOULD BE SET BY THE PROGRAM.
C
         RESTRT = INDEX(KEYWRD,'RESTAR').NE.0
         THIEL  = INDEX(KEYWRD,'NOTHIE').EQ.0
         GEOOK  = INDEX(KEYWRD,'GEO-OK').NE.0
         LOG    = INDEX(KEYWRD,'NOLOG').EQ.0
         LDIIS  = INDEX(KEYWRD,'NODIIS').EQ.0
         SADDLE = INDEX(KEYWRD,'SADDLE').NE.0
         MINPRT = .NOT.SADDLE
         TIME   = INDEX(KEYWRD,'TIME').NE.0
         CONST=1.D0
C
C      THE DAVIDON-FLETCHER-POWELL METHOD IS NOT RECOMMendED
C      BUT CAN BE INVOKED BY USING THE KEYWORD 'DFP'
C
         DFP=INDEX(KEYWRD,'DFP').NE.0
C
C  ORDER OF PRECISION:   'GNORM' TAKES PRECEDENCE OVER 'FORCE', WHICH
C                        TAKES PRECEDENCE OVER 'PRECISE'.
         TOLERG=1.0D0
         if (INDEX(KEYWRD,'PREC') .NE. 0) TOLERG=0.2D0
         if (INDEX(KEYWRD,'FORCE') .NE. 0) TOLERG = 0.1D0
C
C      READ IN THE GRADIENT-NORM LIMIT, if SPECIFIED
C
         if (INDEX(KEYWRD,'GNORM=').NE.0) then
            ROOTV=1.D0
            CONST=1.D-20
            TOLERG = READA(KEYWRD,INDEX(KEYWRD,'GNORM='))
            if (INDEX(KEYWRD,' LET').EQ.0.AND.TOLERG.LT.1.D-2) then
               write(MIOB,'(/,A)')
     1                    '  GNORM HAS BEEN SET TOO LOW, RESET TO 0.01'
               TOLERG=1.D-2
            end if
         else
            ROOTV=SQRT(NVAR+1.D-5)
         end if
         TOLERX = 0.0001D0*CONST
         DELHOF = 0.0010D0*CONST
         TOLERF = 0.002D0*CONST
         TOLRG  = TOLERG
C
C  MINOR BOOK-KEEPING
C
         TLAST=TLEFT
         TX2=SECOND()
         TLEFT=TLEFT-TX2+TIME0
         PRINT  = (INDEX(KEYWRD,'FLEPO').NE.0)
C
C   THE FOLLOWING CONSTANTS SHOULD BE SET TO SOME ARBITARY LARGE VALUE.
C
         DROP  = 1.D15
         FREPF = 1.D15
C
C     AND FINALLY, THE FOLLOWING CONSTANTS ARE CALCULATED.
C
         IHDIM=(NVAR*(NVAR+1))/2
         CNCADD=1.0D00/ROOTV
         if (CNCADD.GT.0.15D00) CNCADD=0.15D00
         ICALCN=NUMCAL
         if (RESTRT) then
            JNRST=1
            MDFP(9)=0
            call DFPSAV(TOTIME,XPARAM,GD,XLAST,FUNCT1,MDFP,XDFP)
            I=TOTIME/1000000.D0
            TOTIME=TOTIME-I*1000000.D0
            TIME0=TIME0-TOTIME
            NSCF=MDFP(5)
            write(IPRT,'(//10X,''TOTAL TIME USED SO FAR:'',
     1    F13.2,'' SECONDS'')') TOTIME
            if (INDEX(KEYWRD,'1SCF') .NE. 0) then
               LAST=1
               LGRAD= INDEX(KEYWRD,'GRAD').NE.0
               call COMPFG (XPARAM,.TRUE.,FUNCT1,.TRUE.,GRAD,LGRAD)
               IFLEPO=13
               EMIN=0.D0
               return
            end if
         end if
C
C   end OF ONCE-ONLY SETUP
C
      end if
C
C     FIRST, WE INITIALIZE THE VARIABLES.
C
      DIISOK=.FALSE.
      IRESET=0
      ABSMIN=1.D6
      FRST=.TRUE.
      ITRY1=0
      ITRY2=0
      JCYC=0
      LNSTOP=1
      IREPET=1
      LIMSCF=.TRUE.
      ALPHA = 1.0D00
      PNORM=1.0D00
      JNRST=0
      CYCMX=0.D0
      COS=0.0D00
      TOTIME=0.D0
      NCOUNT=1
      if ( SADDLE) then
*
*   WE DON'T NEED HIGH PRECISION DURING A SADDLE-POINT CALCULATION.
*
         if (NVAR.GT.0)      GNORM=SQRT(DOT(GRAD,GRAD,NVAR))-3.D0
         if (GNORM.GT.10.D0) GNORM =10.D0
         if (GNORM.GT.1.D0)  TOLERG=TOLRG*GNORM
         write(IPRT,'('' GRADIENT CRITERION IN FLEPO ='',F12.5)')
     1         TOLERG
      end if
      if (NVAR.EQ.1) then
         PVECT(1)=0.01D0
         ALPHA=1.D0
         go to 300
      end if
      TOTIME=0.D0
C
C CALCULATE THE VALUE OF THE FUNCTION -> FUNCT1, AND GRADIENTS -> GRAD.
C NORMAL SET-UP OF FUNCT1 AND GRAD, DONE ONCE ONLY.
C
      call COMPFG (XPARAM,.TRUE.,FUNCT1,.TRUE.,GRAD,.TRUE.)
      call SCOPY (NVAR,GRAD,1,GD,1)
      if (NVAR.NE.0) then
         GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
         GNORMR=GNORM
         if ( LNSTOP.NE.1 .AND. COS.GT.RST .AND.
     1       (JNRST.LT.NRST.OR..NOT.DFP)   .AND.
     2        RESTRT) then
            call SCOPY (NVAR,GD,1,GLAST,1)
         else
            call SCOPY (NVAR,GRAD,1,GLAST,1)
         end if
         GLNORM=GNORM
      end if
      if (GNORM.LT.TOLERG.OR.NVAR.EQ.0) then
         ifLEPO=2
         if (RESTRT) then
            call COMPFG (XPARAM,.TRUE.,FUNCT1,.TRUE.,GRAD,.TRUE.)
         else
            call COMPFG (XPARAM,.TRUE.,FUNCT1,.TRUE.,GRAD,.FALSE.)
         end if
         EMIN=0.D0
         return
      end if
      TX1 =  SECOND()
      TLEFT=TLEFT-TX1+TX2
C     *
C     START OF EACH ITERATION CYCLE ...
C     *
C
      RESET=.FALSE.
      GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
      if (GNORMR.LT.1.D-10)GNORMR=GNORM
      go to 30
   10 continue
      if (COS .LT. RST) then
         do 20 I=1,NVAR
   20    GD(I)=0.5D0
      end if
   30 continue
      JCYC=JCYC+1
      JNRST=JNRST+1
      I80=0
   40 continue
      if (I80.EQ.1    .OR.
     1    LNSTOP.EQ.1 .OR.
     2     COS.LE.RST .OR.(JNRST.GE.NRST.AND.DFP)) then
C
C     *
C     RESTART SECTION
C     *
C
         RESET=.TRUE.
         do 50 I=1,NVAR
C
C  MAKE THE FIRST STEP A WEAK FUNCTION OF THE GRADIENT
C
            STEP=ABS(GRAD(I))*0.0002D0
            STEP=MAX(0.01D0,MIN(0.04D0,STEP))
C#         XD(I)=XPARAM(I)-SIGN(STEP,GRAD(I))
            XD(I)=XPARAM(I)-SIGN(DEL,GRAD(I))
   50    continue
C#      write(MIOB,'(10F8.3)') (XD(I)-XPARAM(I),I=1,NVAR)
C
C THIS call OF COMPFG IS USED TO CALCULATE THE SECOND-ORDER MATRIX IN H
C IF THE NEW POINT HAPPENS TO IMPROVE THE RESULT, THEN IT IS KEPT.
C OTHERWISE IT IS SCRAPPED, BUT STILL THE SECOND-ORDER MATRIX IS O.K.
C
C#      write(MIOB,*) ' RESET HESSIAN'
         call COMPFG (XD,.TRUE.,FUNCT2,.TRUE.,GD,.TRUE.)
         if (.NOT. GEOOK .AND. SQRT(DOT(GD,GD,NVAR))/GNORM.GT.10.
     1        AND.GNORM.GT.20.AND.JCYC.GT.2) then
C
C  THE GEOMETRY IS BADLY SPECIFIED IN THAT MINOR CHANGES IN INTERNAL
C  COORDINATES LEAD TO LARGE CHANGES IN CARTESIAN COORDINATES, AND THESE
C  LARGE CHANGES ARE BETWEEN PAIRS OF ATOMS THAT ARE CHEMIcallY BONDED
C  TOGETHER.
            write(IPRT,'('' GRADIENTS OF OLD GEOMETRY, GNORM='',F13.6)')
     1              GNORM
            write(IPRT,'(6F12.6)') (GRAD(I),I=1,NVAR)
            GDNORM=SQRT(DOT(GD,GD,NVAR))
            write(IPRT,'('' GRADIENTS OF NEW GEOMETRY, GNORM='',F13.6)')
     1              GDNORM
            write(IPRT,'(6F12.6)') (GD(I),I=1,NVAR)
            write(IPRT,'(///20X,''CALCULATION ABANDONED AT THIS POINT!''
     1)')
            write(IPRT,'(//10X,'' SMALL CHANGES IN INTERNAL COORDINATES
     1ARE   '',/10X,'' CAUSING A LARGE CHANGE IN THE DISTANCE BETWEEN'',
     2/   10X,'' CHEMIcallY-BOUND ATOMS. THE GEOMETRY OPTIMIZATION'',/
     3   10X,'' PROCEDURE WOULD LIKELY PRODUCE INCORRECT RESULTS'')')
            call GEOUT(1)
            call filclo  !  J/M
            stop ' FLEPO - SM. changes vs. large changes'
         end if
         NCOUNT=NCOUNT+1
         do 60 I=1,IHDIM
   60    HESINV(I)=0.0D00
         II=0
         do 90 I=1,NVAR
            II=II+I
            DELTAG=GRAD(I)-GD(I)
            DELTAX=XPARAM(I)-XD(I)
            if (ABS(DELTAG).LT.1.D-12) go to 70
            GGD=ABS(GRAD(I))
            if (FUNCT2.LT.FUNCT1) GGD=ABS(GD(I))
            HESINV(II)=DELTAX/DELTAG
            if (HESINV(II).LT.0.0D00.AND.GGD.LT.1.D-12) go to 70
            if (HESINV(II).LT.0.0D00) HESINV(II)=TDEL/GGD
            go to 80
   70       HESINV(II)=0.01D00
   80       continue
            if (GGD.LT.1.D-12) GGD=1.D-12
            PMSTEP=ABS(0.1D0/GGD)
            if (HESINV(II).GT.PMSTEP) HESINV(II)=PMSTEP
   90    continue
         JNRST=0
         if (JCYC.LT.2) COSINE=1.D0
         if (FUNCT2 .GE. FUNCT1) then
            if (PRINT) write (IPRT,100) FUNCT1,FUNCT2
  100       format(' FUNCTION VALUE=',F13.7,
     1           '  WILL NOT BE REPLACED BY VALUE=',F13.7,/10X,
     2           'CALCULATED BY RESTART PROCEDURE',/)
            COSINE=1.D0
         else
            if ( PRINT ) write (IPRT,110) FUNCT1,FUNCT2
  110       format(' FUNCTION VALUE=',F13.7,
     1           ' IS BEING REPLACED BY VALUE=',F13.7/10X,
     2           ' FOUND IN RESTART PROCEDURE'/6X,'THE CORRESPONDING'
     3           ' X VALUES AND GRADIENTS ARE ALSO BEING REPLACED'/)
            FUNCT1=FUNCT2
            call SCOPY (NVAR,XD,1,XPARAM,1)
            call SCOPY (NVAR,GD,1,GRAD  ,1)
            GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
            if (GNORMR.LT.1.D-10) GNORMR=GNORM
         end if
      else
C
C     *
C     UPDATE VARIABLE-METRIC MATRIX
C     *
C
         do 120 I=1,NVAR
            XVAR(I)=XPARAM(I)-XLAST(I)
  120    GVAR(I)=GRAD(I)-GLAST(I)
         call SUPDOT(GG,HESINV,GVAR,NVAR,1)
         YHY=DOT(GG,GVAR,NVAR)
         SY =DOT(XVAR,GVAR,NVAR)
         K=0
C
C    UPDATE ACCORDING TO DAVIDON-FLETCHER-POWELL
C
         if (DFP) then
            do 130 I=1,NVAR
               XVARI=XVAR(I)/SY
               GGI=GG(I)/YHY
               do 130 J=1,I
                  K=K+1
  130       HESINV(K)=HESINV(K)+XVAR(J)*XVARI-GG(J)*GGI
C
C     UPDATE USING THE BFGS FORMALISM
C
         else
            YHY=1.0D0 + YHY/SY
            do 140 I=1,NVAR
               XVARI=XVAR(I)/SY
               GGI=GG(I)/SY
               do 140 J=1,I
                  K=K+1
  140       HESINV(K)=HESINV(K)-GG(J)*XVARI-XVAR(J)*GGI + YHY*XVAR(J)*XV
     1ARI
         end if
      end if
C#      do 191 I=1,IHDIM
C#  191 HTEMP(I)=HESINV(I)
C#      call HQRII(HTEMP, NVAR, NVAR, XTEMP, VECTS)
C#      J=0
C#      do 193 I=1,NVAR
C#      if (XTEMP(I).LT.0.0D0) then
C#      J=J+1
C#      XTEMP(I)=0.00002D0
C#      end if
C#  193 continue
C#      if (J.NE.0) then
C#      do 194 I=1,IHDIM
C#  194 HTEMP(I)=HESINV(I)
C#      call HREFM(NVAR,VECTS,XTEMP,HESINV)
C#      write(MIOB,*) ' ORIGINAL HESSIAN'
C#      call VECPRT(HTEMP,NVAR)
C#      write(MIOB,*) ' REFORMED HESSIAN'
C#      call VECPRT(HESINV,NVAR)
C#      end if
C#      write(MIOB,*) ' EIGENVALUES OF HESSIAN MATRIX'
C#      write(MIOB,'(1X,5G12.6)') (6.951D-3/XTEMP(I),I=1,NVAR)
C
C     *
C     ESTABLISH NEW SEARCH DIRECTION
C     *
      PNLAST=PNORM
C#      call vecprt(hesinv,nvar)
      call SUPDOT(PVECT,HESINV,GRAD,NVAR,1)
      PNORM=SQRT(DOT(PVECT,PVECT,NVAR))
      if (PNORM.GT.1.5D0*PNLAST) then
C
C  TRIM PVECT BACK
C
         do 150 I=1,NVAR
  150    PVECT(I)=PVECT(I)*1.5D0*PNLAST/PNORM
         PNORM=1.5D0*PNLAST
      end if
      DOTT=-DOT(PVECT,GRAD,NVAR)
      do 160 I=1,NVAR
  160 PVECT(I)=-PVECT(I)
      COS=-DOTT/(PNORM*GNORM)
      if (JNRST.EQ.0) go to 190
      if (COS.LE.CNCADD.AND.DROP.GT.1.0D00) go to 170
      if (COS.LE.RST) go to 170
      go to 190
  170 continue
C#      K=0
C#      do 222 I=1,NVAR
C#      do 223 J=1,I-1
C#      K=K+1
C#  223 HESINV(K)=HESINV(K)*0.75D0
C#      K=K+1
C#  222 HESINV(K)=HESINV(K)+0.005D0
C#      go to 241
      PNORM=PNLAST
      if ( PRINT ) write (IPRT,180) COS
  180 format(//,5X, 'SINCE COS=',F9.3,5X,'THE PROGRAM WILL go to RE',
     1'START SECTION',/)
      I80=1
      go to 40
  190 continue
      if ( PRINT ) write (IPRT,200) JCYC,FUNCT1
  200 format(1H , 'AT THE BEGINNING OF CYCLE',I5, '  THE FUNCTION VA
     1LUE IS ',F13.6/, '  THE CURRENT POINT IS ...')
      if (PRINT) write (IPRT,210) GNORM,COS
  210 format( '  GRADIENT NORM = ',F10.4/,'  ANGLE COSINE =',F10.4)
      if ( PRINT ) then
         write(MIOB,220)
  220    format('  THE CURRENT POINT IS ...')
         NTO6=(NVAR-1)/6+1
         IINC1=-5
         do 270 I=1,NTO6
            write(MIOB,'(/)')
            IINC1=IINC1+6
            IINC2=MIN(IINC1+5,NVAR)
            write(MIOB,230) (J,J=IINC1,IINC2)
            write(MIOB,240) (XPARAM(J),J=IINC1,IINC2)
            write(MIOB,250) (GRAD(J),J=IINC1,IINC2)
            write(MIOB,260) (PVECT(J),J=IINC1,IINC2)
  230       format(1H ,3X,  1HI,9X,I3,9(8X,I3))
  240       format(1H ,1X, 'XPARAM(I)',1X,F9.4,2X,9(F9.4,2X))
  250       format(1H ,1X, 'GRAD  (I)',F10.4,1X,9(F10.4,1X))
  260       format(1H ,1X, 'PVECT (I)',2X,F10.6,1X,9(F10.6,1X))
  270    continue
      end if
      LNSTOP=0
      ALPHA=ALPHA*PNLAST/PNORM
      call SCOPY (NVAR,GRAD,  1,GLAST,1)
      call SCOPY (NVAR,XPARAM,1,XLAST,1)
      GLNORM=GNORM
      if (JNRST.EQ.0) ALPHA=1.0D00
      DROP=ABS(ALPHA*DOTT)
      if (PRINT)write (IPRT,280) DROP
  280 format(1H , 13H -ALPHA.P.G =,F18.6,/)
      if (JNRST.NE.0.AND.DROP.LT.DELHOF) then
C
C   HERBERT'S TEST: THE PREDICTED DROP IN ENERGY IS LESS THAN DELHOF
C   IF PASSED, call COMPFG TO GET A GOOD SET OF EIGENVECTORS, THEN EXIT
C
         if (MINPRT) write (IPRT,290)
  290    format(//,10X,'HERBERTS TEST SATISFIED - GEOMETRY OPTIMIZED')
C
C   FLEPO IS endING PROPERLY. THIS IS IMMEDIATELY BEFORE THE return.
C
         LAST=1
         call COMPFG (XPARAM,.TRUE.,FUNCT,.TRUE.,GRAD,.FALSE.)
         IFLEPO=3
         TIME0=TIME0-TOTIME
         EMIN=0.D0
         return
      end if
      BETA =ALPHA
      SMVAL=FUNCT1
      DROPN=-ABS(DROP/ALPHA)
C
C    UPDATE GEOMETRY USING THE G-DIIS PROCEDURE
C
      if (DIISOK) then
         OKF=.TRUE.
         IC=1
      else
         OKF=.FALSE.
         IC=2
      end if
  300 call LINMIN(XPARAM,ALPHA,PVECT,NVAR,FUNCT1,OKF,IC,DROPN)
      if (NVAR.EQ.1) then
         write(MIOB,
     %      '('' ONLY ONE VARIABLE, THEREFORE ENERGY A MINIMUM'')')
         LAST=1
         LGRAD=(INDEX(KEYWRD,'GRAD').NE.0)
         call COMPFG (XPARAM,.FALSE.,FUNCT,.FALSE.,GRAD,LGRAD)
         IFLEPO=14
         EMIN=0.D0
         return
      end if
C   WE WANT ACCURATE DERIVATIVES AT THIS POINT
C
C   LINMIN DOES NOT GENERATE ANY DERIVATIVES, THEREFORE COMPFG MUST BE
C   CALLED TO end THE LINE SEARCH
C
C  IF THE DERIVATIVES ARE TO BE CALCULATED USING FULL SCF'S, THEN CHECK
C  WHETHER TO do FULL SCF'S (CRITERION FROM FLEPO: GRAD IS NULL).
C
      if (IRESET.GT.10    .OR.
     %    GNORM.LT.40.D0 .AND.
     %    GNORM/GNORMR.LT.0.33D0) then
         IRESET=0
         GNORMR=0.D0
         do 310 I=1,NVAR
  310    GRAD(I)=0.D0
      end if
      IRESET=IRESET+1
C
C
C     RESTORE TO STANDARD VALUE BEFORE COMPUTING THE GRADIENT
      RESET=.FALSE.
      if (THIEL) then
         call COMPFG (XPARAM, IC.NE.1,  SCRAP, .TRUE. ,GRAD,.TRUE.)
      else
         call COMPFG (XPARAM, .TRUE.,  FUNCT1, .TRUE. ,GRAD,.TRUE.)
      end if
      if (LDIIS) then
C
C  UPDATE GEOMETRY AND GRADIENT AFTER MAKING A STEP USING LINMIN
C
         do 320 I=1,NVAR
            XTEMP(I)=XPARAM(I)
  320    GTEMP(I)=GRAD(I)
         call DIIS(XTEMP, XPARAM, GTEMP, GRAD, HDIIS, FUNCT1,HESINV, NVA
     1R, FRST)
         if (HDIIS.LT.FUNCT1.AND.
     1SQRT(DOT(GTEMP,GTEMP,NVAR)) .LT. SQRT(DOT(GRAD,GRAD,NVAR))) then
            do 330 I=1,NVAR
               XPARAM(I)=XTEMP(I)
  330       GRAD(I)=GTEMP(I)
            DIISOK=.TRUE.
         else
            DIISOK=.FALSE.
         end if
      end if
      GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
      if (GNORMR.LT.1.D-10) GNORMR=GNORM
      NCOUNT=NCOUNT+1
      if ( .NOT. OKF) then
         LNSTOP = 1
         if (MINPRT)write (IPRT,'(/,20X, ''NO POINT LOWER IN ENERGY '',
     1    ''THAN THE STARTING POINT '',/,20X,''COULD BE FOUND '',
     2    ''IN THE LINE MINIMIZATION'')')
         FUNCT1=SMVAL
         ALPHA=BETA
         call SCOPY (NVAR,GLAST,1,GRAD  ,1)
         call SCOPY (NVAR,XLAST,1,XPARAM,1)
         if (JNRST.EQ.0) then
            write (IPRT,340)
  340       format(1H ,//,20X, 'SINCE COS WAS JUST RESET,THE SEARCH',
     1        ' IS BEING endED')
C
C           FLEPO IS endING BADLY. THIS IS IMMEDIATELY BEFORE THE return
C
            LAST=1
            call COMPFG (XPARAM, .TRUE., FUNCT, .TRUE. ,GRAD,.TRUE.)
            IFLEPO=4
            TIME0=TIME0-TOTIME
            EMIN=0.D0
            return
         end if
         if (PRINT) write (IPRT,350)
  350    format(1H ,20X, 'COS WILL BE RESET AND ANOTHER '
     1    ,'ATTEMPT MADE')
         COS=0.0D00
         go to 470
      end if
      XN=SQRT(DOT(XPARAM,XPARAM,NVAR))
      TX=ABS(ALPHA*PNORM)
      if (XN.NE.0.0D00) TX=TX/XN
      TF=ABS(SMVAL-FUNCT1)
      if (ABSMIN-SMVAL.LT.1.D-7) then
         ITRY1=ITRY1+1
         if (ITRY1.GT.10) then
            write(MIOB,
     %         '(//,'' HEAT OF FORMATION IS ESSENTIALLY STATIONARY'')')
            go to 460
         end if
      else
         ITRY1=0
         ABSMIN=SMVAL
      end if
      IOUT=0
      if (PRINT) write(MIOB,360) NCOUNT,COS,TX*XN,ALPHA,-DROP,-TF,GNORM
  360 format(/,'           NUMBER OF COUNTS =',I6,
     1'         COS    =',F11.4,/,
     2        '  ABSOLUTE  CHANGE IN X     =',F13.6,
     3'  ALPHA  =',F11.4,/,
     4        '  PREDICTED CHANGE IN F     =  ',G11.4,
     5'  ACTUAL =  ',G11.4,/,
     6        '  GRADIENT NORM             =  ',G11.4,//)
      if (TX.LE.TOLERX) then
         if (MINPRT) write (IPRT,370)
  370    format(' TEST ON X SATISFIED')
         go to 400
      end if
      if (TF.LE.TOLERF) then
C#         write(MIOB,*) TF,TOLERF
         if (MINPRT) write (IPRT,380)
  380    format(' HEAT OF FORMATION TEST SATISFIED')
         go to 400
      end if
      if (GNORM.LE.TOLERG*ROOTV) then
         if (MINPRT) write (IPRT,390)
  390    format(' TEST ON GRADIENT SATISFIED')
         go to 400
      end if
      go to 470
  400 do 440 I=1,NVAR
         if (ABS(GRAD(I)).GT.TOLERG) then
            IREPET=IREPET+1
            if (IREPET.GT.1) go to 410
            FREPF=FUNCT1
            COS=0.0D00
  410       if (MINPRT) write (IPRT,420) TOLERG
  420       format(20X,'HOWEVER, A COMPONENT OF GRADIENT IS ',
     1     'LARGER THAN',F6.2 ,/)
            if (ABS(FUNCT1-FREPF).GT.EINC) IREPET=0
            if (IREPET.GT.IGG1) then
               write (IPRT,430) IGG1,EINC
  430          format(10X,' THERE HAVE BEEN',I2,' ATTEMPTS TO REDUCE TH
     1E ',' GRADIENT.',/10X,' DURING THESE ATTEMPTS THE ENERGY DROPPED',
     2' BY LESS THAN',F4.1,' KCAL/MOLE',/
     310X,' FURTHER CALCULATION IS NOT JUSTIFIED AT THIS TIME.',/
     410X,' TO CONTINUE, START AGAIN WITH THE WORD "PRECISE"' )
               LAST=1
C!!
               call COMPFG (XPARAM,.TRUE.,FUNCT,.TRUE.,GRAD,.FALSE.)
               IFLEPO=8
               TIME0=TIME0-TOTIME
               EMIN=0.D0
               return
            else
               go to 470
            end if
         end if
  440 continue
      if (MINPRT) write (IPRT,450)
  450 format( 23H PETERS TEST SATISFIED )
  460 LAST=1
      call COMPFG (XPARAM,.TRUE.,FUNCT,.TRUE.,GRAD,.FALSE.)
      IFLEPO=6
      TIME0=TIME0-TOTIME
      EMIN=0.D0
      return
C
C   ALL TESTS HAVE FAILED, WE NEED TO do ANOTHER CYCLE.
C
  470 continue
      BSMVF=ABS(SMVAL-FUNCT1)
      if (BSMVF.GT.10.D00) COS = 0.0D00
      DEL=0.002D00
      if (BSMVF.GT.1.0D00) DEL=DELL/2.0D00
      if (BSMVF.GT.5.0D00) DEL=DELL
      TX2 = SECOND ()
      TCYCLE=TX2-TX1
      TX1=TX2
C
C end OF ITERATION LOOP, EVERYTHING IS STILL O.K. SO go to
C NEXT ITERATION, if THERE IS ENOUGH TIME LEFT.
C
      if (TCYCLE.LT.100000.D0) CYCMX=MAX(CYCMX,TCYCLE)
      TLEFT=TLEFT-TCYCLE
      if (TLEFT.LT.0)     TLEFT=-0.1D0
      if (TCYCLE.GT.1.D5) TCYCLE=0.D0
      if (TLAST-TLEFT.GT.TDUMP) then
         TOTIM=TOTIME   +   SECOND()-TIME0
         TLAST=TLEFT
         MDFP(9)=2
         RESFIL=.TRUE.
         MDFP(5)=NSCF
         call DFPSAV(TOTIM,XPARAM,GD,XLAST,FUNCT1,MDFP,XDFP)
      end if
      if (RESFIL) then
         if (MINPRT) write(MIOB,480) MIN(TLEFT,9999999.9D0),
     1                              MIN(GNORM,999999.999D0),FUNCT1
  480    format(' RESTART FILE WRITTEN,   TIME LEFT:',F9.1,
     1' GRAD.:',F10.3,' HEAT:',G13.7)
         RESFIL=.FALSE.
      else
         if (MINPRT) then
            write(MIOB,490) JCYC,MIN(TCYCLE,9999.99D0),
     1      MIN(TLEFT,9999999.9D0),MIN(GNORM,999999.999D0),FUNCT1
            write(*,490) JCYC,MIN(TCYCLE,9999.99D0),
     1      MIN(TLEFT,9999999.9D0),MIN(GNORM,999999.999D0),FUNCT1
         endif
         if (LOG) write(MIOE,490)    JCYC,MIN(TCYCLE,9999.99D0),
     1      MIN(TLEFT,9999999.9D0),MIN(GNORM,999999.999D0),FUNCT1
  490    format(' CYCLE:',I4,' TIME:',F7.2,' TIME LEFT:',F9.1,
     1' GRAD.:',F10.3,' HEAT:',G13.7)
      end if
      if (TLEFT.GT.SFACT*CYCMX) go to 10
      write(IPRT,500)
  500 format(20X, 42HTHERE IS NOT ENOUGH TIME FOR ANOTHER CYCLE,/,30X,
     118HNOW GOING TO FINAL)
      TOTIM=TOTIME   +   SECOND()-TIME0
      MDFP(9)=1
      MDFP(5)=NSCF
      call DFPSAV(TOTIM,XPARAM,GD,XLAST,FUNCT1,MDFP,XDFP)
      IFLEPO=-1
      return
C
C
      end
