      SUBROUTINE EF(XPARAM, NVAR, FUNCT)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension XPARAM(MAXPAR)
**********************************************************************
*
*   BAKER IS A QUASI NEWTON RAPHSON OPTIMIZATION ROUTINE BASED ON
*         JON BAKERS ALGORITHM (J.COMP.CHEM. 7, 385).
*
*  ON ENTRY XPARAM = VALUES OF PARAMETERS TO BE OPTIMISED.
*           NVAR   = NUMBER OF PARAMETERS TO BE OPTIMISED.
*
*  ON EXIT  XPARAM = OPTIMISED PARAMETERS.
*           FUNCT  = HEAT OF FORMATION IN KCAL/MOL.
*
*                    IMPLEMENTED IN MOPAC AUGUST 1990 BY JJPS
*
**********************************************************************
C
C   THE BAKER ROUTINE IS CONTROLED BY THE FOLLOWING KEYWORDS :
C   EF,TS           : SEARCH FOR MINIMUM OR TRANSITION STRUCTURE
C   CYCLES          : MAXIMUM NUMBER OF ITERATIONS IN OPTIMIZATION
C   RECALC=N        : NUMBER OF CYCLES BETWEEN HESSIAN CALCULATION
C   DMAX=N.N        : MAXIMUM STEPSIZE IN ANG OR RAD
C   IUPD=0,1,2      : HESSIAN UPDATING ; NONE, POWELL, BFGS
C   HESS=0,1,2,3    : START HESSIAN ; DIAGONAL, CALCULATE,
C                     READ FROM DISK, READ FROM INPUT.
C                     NOTE THAT CALCULATION OF SELECTIVE ROWS (COLOUMNS)
C                     OF THE HESSIAN COMBINED WITH DIAGONAL GUESS IS
C                     CONTROLED BY SETTING OPTIMIZATION PARAMETERS = 2.
C                     IF ANY OPT PARAMETERS = 2, then ABOVE HESS OPTION
C                     IS OVERRIDDEN.
C   MODE=N          : FOLLOW THE N'TH LOWEST HESSIAN MODE TOWARD TS
C   NONR            : DO NOT USE NEWTON-RAPHSON STEP IF FORMD FAILS AND
C                     MAX GRADIENT IS LESS THAN RCUT.
C   EIGINV          : USE HESSIAN EIGENVALUE REVERSION RATHER THAN
C                     DIAGONAL SHIFT IF FORMD FAILS AND NONR OR GNORM
C                     IS LARGER THAN RCUT.
C   ORIDE           : USE WHATEVER LAMDA(S) FORMD PRODUCES, EVEN IF THEY
C                     ARE 'UNACCEPTABLE'.
C   GNORM=N.N       : EXIT WHEN LARGEST GRADIENT IS LESS THAN N.N
C   PRNT=0-5        : 0 : PRINT SUMMERY EVERY CYCLE
C                     1 : + GEOMETRY AND TS MODE
C                     2 : + GRADIENTS
C                     3 : + HESSIAN EIGENVALUES AND -VECTORS
C                     4 : + HESSIAN MATRIX
C                     5 : MAXIMUM PRINTING = DEBUG PRINTING
C
C   DEFAULTS :
C   PRNT     : 0
C   CYCLES   : 100
C   DMAX     : 0.20
C   RECALC   : 999999
C   IUPD     : 2 FOR MINIMUM SEARCH, 1 FOR TS SEARCH
C   HESS     : 0 FOR MINIMUM SEARCH, 1 FOR TS SEARCH
C   MODE     : 1, ONLY USED IN TS SEARCH
C   GNORM    : 0.4 OR 0.01 IF "PRECISE" IS SPECIFIED
C   NONR     : FALSE
C   EIGINV   : FALSE
C   ORIDE    : FALSE
C
C   HINTS AND COMMENTS :
C   NORMALLY THE DEFAULT OPTIONS ARE SUFFICIENT. NONR, EIGINV AND
C      ORIDE OPTIONS ARE SPECIAL AND SHOULD RARELY BE TOUCHED.
C   CONVERGENCE CRITERIA IS THE SAME AS FOR "SIGMA" OPTIMIZATION,
C      NAMELY THAT ALL GRADIENTS MUST BE BELOW 0.4 OR 0.01 IF
C      "PRECISE" IS SPECIFIED. RESULTS OBTAINED WITHOUT THE "PRECISE"
C      OPTION TURNED ON MAY BE UNRELIABLE DUE TO ACCUMULATED ERRORS
C      FROM VARIOUS PARTS OF THE PROGRAM (SEE E.G. J.COMP.CHEM. 9,X)
C      BUT CAN BE USED IN PRELIMINARY STAGES. ONLY RESULTS FROM
C      PRECISE RUNS SHOULD EVER BE PUBLISHED!
C      WITH THE PRECISE KEYWORD GEOMETRIES ARE USUALLY CONVERGED TO
C      0.0001 IN BOND LENGTH, 0.01 IN BOND ANGLES AND 0.1 IN TORSIONAL
C      ANGLES. ENERGIES ARE CONVERGED TO AT LEAST 0.001 KCAL/MOL.
C   CONVERGENCE CAN BE MONITORED IN TWO WAYS : THE GRADIENT NORM SHOULD
C      GO TOWARD ZERO, AND THE LAMDA VALUES SHOULD ALSO GO TO ZERO.
C   HESSIAN AND GEOMETRY WILL BE SAVED ON DISK (AND CAN BE RESTARTED) :
C      1) AFTER HESSIAN CALCULATION (PARTLY OR WHOLE)
C      2) IN AN OPTIMIZATION EVERY X MIN. CPU USED (MACHINE DEPENDENT)
C      3) BEFORE TERMINATION, EITHER DUE TO CONVERGENCE, TIME OR CYCLES
C   IF HESS=2 THEN THE GEOMETRY IS TAKEN FROM THE INPUT WHILE THE HESSIA
C      IS READ FROM DISK. UPON RESTART BOTH GEOMETRY AND HESSIAN IS
C      AUTOMATICALLY PICKED UP FROM DISK.
C   TS SEARCH REQUIRES HESS=1, 2 OR 3, IT CANNOT WORK WITH HESS=0
C   HESS=1 FOR BAKERMIN IS RECOMMENDED FOR DIFFICULT CASES (E.G. CYCLIC)
C   HESS=3 WILL READ THE HESSIAN FROM INPUT, MOSTLY USED FOR TESTING.
C      TO READ IT IN, PLACE THE WORD 'HESSIAN' SOMEWHERE IN THE INPUT
C      FILE IN format(A7), FOLLOWED BY THE HESSIAN ELEMENTS IN
C      format(8F10.2).
C   THE ROUTINE (ESPECIALLY TS SEARCH) WORKS BEST IN INTERNAL COORDINATE
C   MOST COMMON ERROR FROM BAKER IS "UNABLE TO DETERMINE LAMDA IN
C      FORMD". THIS WILL OFTEN HAPPENS IF :
C      1) SYMMETRY CONSTRAINED VARIABLES ARE MARKED FOR OPTIMIZATION
C      2) VARIABLE DETERMINED BY SYMMETRY IS MARKED FOR OPTIMIZATION
C      3) MORE VARIABLES THAN INTERNAL DEGREES OF FREEDOM SPECIFIED
C      4) TOO LARGE STEP SIZE HAS BEEN CHOSEN, UPDATING GOES HAYWIRE
C      5) BAD HESSIAN IS USED AS START HESSIAN
C   EVEN WHEN THESE REASONS HAVE BEEN ELIMINATED, THE FAIL OCCURS QUITE
C      FREQUENTLY IN THE ***** ENVIROMENT, AND USUALLY IF SOFT MODES
C      SUCH AS METHYL ROTATIONS ARE PRESENT. SOME PRELIMINARY PATCHES
C      FOR CIRCUMVENTING THESE FAILURES HAVE BEEN IMPLEMENTED IN THIS
C      VERSION. IF THE STRUCTURE IS CLOSE TO OPTIMIZED (MAX GRAD.LT.RCUT
C      = 10), THEN SIMPLY TAKE THE NEWTON-RAPHSON STEP, UNLESS NONR IS
C      TURNED ON. IF MAX GRAD.GT.RCUT OR NONR IS ON, THEN USE EITHER
C      A HESSIAN SHIFT OR A HESSIAN EIGENVALUE INVERSION.
C      THE DEFAULT, AND ONLY OPTION FOR MINIMUM SEARCHES, IS TO
C      INCREASE THE HESSIAN DIAGONAL WITH LOWEST HESSIAN EIGENVALUE
C      (SECOND LOWEST EIGENVALUE IF TS SEARCH) + 2*EIGMIN.
C      THIS SHOULD ALWAYS WORK FOR A MIN SEARCH, BUT MAY GO WRONG IN A
C      TS SEARCH. IF THIS HAPPENS ONE CAN TRY TO TURN ON THE EIGINV
C      OPTION. THIS WILL REVERSE THE SIGN OF ALL NEGATIVE EIGENVALUES
C      OTHER THAN THE TS MODE BEING FOLLOWED. AGAIN THERE IS NO GARANTI
C      THAT THIS IS WILL LEAD TO A 'BETTER' GEOMETRY.
C      FINALLY AN OVERRIDE SWITCH IS IMPLEMENTED. IF KEYWORD ORIDE IS
C      USED THEN WHATEVER LAMDA'S FORMD PRODUCES WILL BE USED,
C      EVEN IF THEY ARE 'UNACCEPTABLE', AS DISCUSSED BY SIMONS.
C   THESE PATCHES ARE AS MENTIONED STILL EXPERIMENTAL, BUT HAVE
C      BEEN FOUND TO WORK IN MOST CASES.
C   INCREASING DMAX CAN LEAD TO FASTER CONVERGENCE BUT CAN ALSO MAKE
C      THE OPTIMIZATION GO BAD VERY FAST. FURTHERMORE THE HESSIAN
C      UPDATING DETERIORATE WHEN USING LARGE STEPSIZES. REDUCING THE
C      STEPSIZE TO 0.10 OR 0.05 IS RECOMMENDED WHEN ENCOUNTERING
C      CONVERGENCE PROBLEMS.
C   RECALC=N CALCULATES THE HESSIAN EVERY N STEPS IN THE OPTIMIZATION.
C      FOR SMALL N THIS IS COSTLY BUT IS ALSO VERY EFFECTIVE IN
C      TERMS OF CONVERGENCE. RECALC=10 AND DMAX=0.10 CAN BE USEFUL
C      FOR DIFFICULT CASES. IN EXTREME CASES RECALC=1 AND DMAX=0.05
C      WILL ALWAYS FIND A STATIONARY POINT, IF IT EXIST.
C   IF READING IN GNORM=N.N NOTICE THAT IT IS MAX GRADIENT THRESHOLD
C      THAT IS BEING MODIFIED, NOT GRADIENT NORM. IGNORE THE MESSAGE
C      ABOUT "FLEPO EXIT WHEN...".
C   IUPD SHOULD VERY RARELY BE TOUCHED. IUPD=1 CAN BE USED IN MINIMUM
C      SEARCHES IF THE THE MESSAGE "HEREDITARY POSITIVE DEFINITENESS
C      ENDANGERED. UPDATE SKIPPED THIS CYCLE" OCCURS EVERY CYCLE FOR
C      10-20 ITERATIONS. NEVER USE IUPD=2 FOR TS SEARCH!
C
C   **** NOTICE THAT VERY LITTLE ERROR CHECKING IS DONE.
C   **** ONLY RUNS REQUESTING INCOMPATIBLE OPTIONS ARE TERMINATED.
C   **** A MINIMUM KNOWLEDGE OF THE THEORY BEHIND THE ROUTINE
C   **** IS RECOMMENDED TO GET THE OPTIMUM USAGE OUT OF IT.
C
      common /MESAGE/ IFLEPO,ISCF
      common /GEOVAR/ NDUM,LOC(2,MAXPAR), IDUMY, XARAM(MAXPAR)
      common /GEOM  / GEO(3,NUMATM)
      common /GEOSYM/ NDEP,LOCPAR(MAXPAR),IDEPFN(MAXPAR),LOCDEP(MAXPAR)
      common /ISTOPE/ AMS(107)
      common /LAST  / LAST
      common /KEYWRD/ KEYWRD
      common /TIME  / TIME0
      common /GRADNT/ GRAD(MAXPAR),GNFINA
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      common /NUMCAL/ NUMCAL
      common /TIMDMP/ TLEFT, TDUMP
      common /SIGMA2/ GNEXT1(MAXPAR), GMIN1(MAXPAR)
      common /NLLCOM/ HESS(MAXPAR,MAXPAR),BMAT(MAXPAR,MAXPAR),
     1                PMAT(MAXPAR**2)
      common /SCRACH/ PVEC
      common /SCFTYP/ EMIN, LIMSCF
      common /OPTEF/  OLDF(MAXPAR),D(MAXPAR),VMODE(MAXPAR),
     1                U(MAXPAR,MAXPAR),DD,MODE,NSTEP,NEGREQ,IPRNT
      dimension IPOW(9), EIGVAL(MAXPAR),TVEC(MAXPAR),SVEC(MAXPAR),
     1FX(MAXPAR),HESSC(MAXHES),UC(MAXPAR**2), PVEC(MAXPAR*MAXPAR)
      logical RESTRT,SCF1,FAIL,HSSHFT,EIGINV,NR,ORIDE
      logical AIGEO, LIMSCF, LOG
      character KEYWRD*241
      character HESWRD*7,WRD*7
      character *80 getnam
      equivalence(IPOW(1),IHESS)
      data HESWRD   /'HESSIAN'/
      data  ICALCN,ZZERO,ONE,TWO    /0,0.D0,1.D0,2.D0/
      data EIGMIN,RCUT /1.D-1,1.D+01/
C     DGHSX IS HESSIAN DIAGONAL FOR HESS=0. FOR STRETCHING, ANGLE,
C     DIHEDRAL.  THE VALUES GIVEN BELOW ARE 'OPTIMUM' FOR CYCLOHEXANONE
      DGHSS=2500.D0
      DGHSA= 900.D0
      DGHSD=1000.D0
C     GET ALL INITIALIZATION DATA
      if (ICALCN.NE.NUMCAL) then
         NVAR=ABS(NVAR)
         LDUMP=0
         ICALCN=NUMCAL
         AIGEO=.FALSE.
         RESTRT=(INDEX(KEYWRD,'RESTART') .NE. 0)
         LOG    = INDEX(KEYWRD,'NOLOG').EQ.0
         SCF1=(INDEX(KEYWRD,'1SCF') .NE. 0)
         ORIDE=(INDEX(KEYWRD,'ORIDE') .NE. 0)
         NSTEP=0
         IHESS=0
         MODE =0
         LAST=0
         NTIME=0
         ILOOP=1
         XINC=0.00529167D0
         IMIN=INDEX(KEYWRD,' EF')
         if (IMIN.NE.0) then
            IGTHES=0
            IUPD  =2
            NEGREQ=0
         end if
         LIMSCF=.TRUE.
         ITS=INDEX(KEYWRD,' TS')
         if (ITS.NE.0) then
            IGTHES=1
            IUPD  =1
            NEGREQ=1
         end if
         I=INDEX(KEYWRD,' DEBUG')*INDEX(KEYWRD,' BAKER')
         IPRNT=0
         if (I.NE.0) IPRNT=2
         IRECLC=999999
         I=INDEX(KEYWRD,' RECALC=')
         if (I.NE.0) IRECLC = READA(KEYWRD,I)
         I=INDEX(KEYWRD,' IUPD=')
         if (I.NE.0) IUPD = READA(KEYWRD,I)
         I=INDEX(KEYWRD,' MODE=')
         if (I.NE.0) MODE = READA(KEYWRD,I)
         DMAX=0.2
         I=INDEX(KEYWRD,' DMAX=')
         if (I.NE.0) DMAX = READA(KEYWRD,I)
         TOL2=4.D-1
         if (INDEX(KEYWRD,' PREC') .NE. 0) TOL2=1.D-2
         I=INDEX(KEYWRD,' GNORM=')
         if (I.NE.0) TOL2 = READA(KEYWRD,I)
         if (INDEX(KEYWRD,' LET').EQ.0.AND.TOL2.LT.0.01D0) then
            write(MIOB,'(/,A)')
     1                  '  GNORM HAS BEEN SET TOO LOW, RESET TO 0.01'
            TOL2=0.01D0
         end if
         HSSHFT=.TRUE.
         NR=.TRUE.
         if (INDEX(KEYWRD,'NONR') .NE. 0) NR=.FALSE.
         EIGINV=.FALSE.
         EIGINV=(INDEX(KEYWRD,'EIGINV') .NE. 0)
         if (EIGINV) then
            HSSHFT=.FALSE.
            if (IPRNT.EQ.2) then
               write(MIOB,*) ' '
               write(MIOB,*)
     %                '    HESSIAN EIGENVALUE INVERSION TURNED ON'
            end if
         end if
         I=INDEX(KEYWRD,' HESS=')
         if (I.NE.0) IGTHES = READA(KEYWRD,I)
         TIME1=TIME0
         TIME2=TIME1
C   DONE WITH ALL INITIALIZING STUFF.
C   CHECK THAT OPTIONS REQUESTED ARE RESONABLE
         if (NVAR.GT.(3*NUMATM-6)) write(MIOB,10)
   10    format(/,'*** WARNING! MORE VARIABLES THAN DEGREES OF FREEDOM',
     1/)
         if ((ITS.NE.0).AND.(IUPD.EQ.2)) then
            write(MIOB,*) ' TS SEARCH AND BFGS UPDATE WILL NOT WORK'
            call filclo ! J/M
            stop ' in EF - TS search and BFGS update won''t work'
         end if
         if ((ITS.NE.0).AND.(IGTHES.EQ.0)) then
            write(MIOB,*)
     %         ' TS SEARCH REQUIRE BETTER THAN DIAGONAL HESSIAN'
            call filclo ! J/M
            stop ' in EF - TS search requires more than diag. Hessian'
         end if
         if ((IMIN.NE.0).AND.(EIGINV)) then
            write(MIOB,*) ' MIN SEARCH AND EIGINV WILL NOT WORK'
            call filclo ! J/M
            stop ' in EF - MIN SEARCH AND EIGINV WILL NOT WORK'
         end if
         if ((IGTHES.LT.-1).OR.(IGTHES.GT.3)) then
            write(MIOB,*) ' UNRECOGNIZED HESS OPTION',IGTHES
            call filclo ! J/M
            stop ' in EF - UNRECOGNIZED HESS OPTION'
         end if
         if (IGTHES.EQ.2) go to 90
         if (RESTRT) then
C
C   RESTORE DATA. I INDICATES (HESSIAN RESTART OR OPTIMIZATION
C   RESTART). IF I .GT. 0 then HESSIAN RESTART AND I IS LAST
C   STEP CALCULATED IN THE HESSIAN. IF I .LE. 0 then J (NSTEP)
C   IN AN OPTIMIZATION HAS BEEN DONE.
C
            IPOW(9)=0
            call BKRSAV(TT0,HESS,FUNCT,GRAD,XPARAM,PMAT,I,J,BMAT,IPOW)
            K=TT0/1000000.D0
            TIME0=TIME0-TT0+K*1000000.D0
            if (I .GT. 0) then
               ILOOP=I
               IGTHES=4
               NSTEP=J
               write(MIOB,
     %              '(10X,''RESTARTING HESSIAN AT POINT'',I4)') ILOOP
               if (NSTEP.NE.0)
     %             write(MIOB,'(10X,''IN OPTIMIZATION STEP'',I4)')
     1             NSTEP
            else
               NSTEP=J
               write(MIOB,
     %          '(//10X,''RESTARTING OPTIMIZATION AT STEP'',I4)') NSTEP
               do 20 I=1,NVAR
   20          GRAD(I)=ZZERO
               call COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
               go to 230
            end if
         else
C   NOT A RESTART, WE NEED TO GET THE GRADIENTS
            do 30 I=1,NVAR
   30       GRAD(I)=ZZERO
            call COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
         end if
      end if
      if (SCF1) then
         GNFINA=SQRT(DOT(GRAD,GRAD,NVAR))
         IFLEPO=1
         return
      end if
      if (ITS.EQ.0) then
C     CHECK THAT GEOMETRY IS NOT ALREADY OPTIMIZED
         RMX=ZZERO
         do 40 I=1,NVAR
            if (ABS(GRAD(I)).GT.RMX) RMX=ABS(GRAD(I))
   40    continue
         if (RMX.LT.TOL2) then
            IFLEPO=2
            LAST=1
            return
         end if
      end if
C     WE NOW HAVE THE GRADIENTS FOR THE CURRENT GEOMETRY,
C     GET THE HESSIAN. DEPENDING ON IGTHES WE GET IT FROM :
C
C     -1 : DIAGONAL MATRIX, DGHSX*I, BUT ROWS AND COLOUMNS OF
C          HESSIAN CALCULATED ACCORDING TO VARIABLES MARKED '2'
C      0 : DIAGONAL MATRIX, DGHSX*I (DEFAULT FOR MIN-SEARCH)
C      1 : CALCULATE IT NUMERIcallY (DEFAULT FOR TS-SEARCH)
C      2 : READ IN FROM FTN009
C      3 : READ IT IN FROM INPUT FILE (FTN005)
C      4 : READ IN FROM FTN009 (DURING RESTART, PARTLY OR WHOLE,
C          ALREADY DONE AT THIS POINT)
C
      if (ILOOP .LE. 0) go to 230
C     if ILOOP IS .LE.0 THIS IS AN OPTIMIZATION RESTART
C     AND HESSIAN SHOULD ALREADY BE AVAILABLE
      if ((IGTHES.EQ.0).OR.(IGTHES.EQ.-1)) then
         if (IPRNT.EQ.2) then
            if (IGTHES.EQ.-1) write(MIOB,50)
            if (IGTHES.EQ. 0) write(MIOB,60)
         end if
   50    format(/,10X,'HESSIAN PARTLY CALCULATED',/)
   60    format(/,10X,'DIAGONAL MATRIX USED AS START HESSIAN',/)
         do 70 I=1,NVAR
            do 70 J=1,NVAR
               HESS(I,J)=ZZERO
   70    continue
         IJ=1
         do 80 J=1,NUMATM
            do 80 I=1,3
               if (LOC(2,IJ).EQ.I.AND.LOC(1,IJ).EQ.J) then
                  if (I.EQ.1) HESS(IJ,IJ)=DGHSS
                  if (I.EQ.2) HESS(IJ,IJ)=DGHSA
                  if (I.EQ.3) HESS(IJ,IJ)=DGHSD
                  IJ=IJ+1
               end if
   80    continue
         IJ=IJ-1
         if (IJ.NE.NVAR)
     1      WRITE(*,*) 'ERROR IN IGTHES=-1,IJ,NVAR',IJ,NVAR
      end if
C
   90 if (IGTHES.EQ.2) then
         write(MIOB,100)
  100    format(/,10X,'HESSIAN READ FROM DISK',/)
         IPOW(9)=0
C        USE DUMMY ARRAY FOR call EXCEPT FOR HESSIAN
         call BKRSAV(TDM,HESS,FDMY,GNEXT1,GMIN1,PMAT,IIDUM,J,BMAT,IPOW)
      end if
      if (IGTHES .EQ.3) then
         write(MIOB,110)
  110    format(/,10X,'HESSIAN READ FROM INPUT',/)
c
      open(unit=MION,file=GETNAM('FOR008'),status='OLD',
     1      form='FORMATTED')
c
         rewind MION
  120    read(MION,150,end=630,ERR=630) WRD
         if (WRD .EQ. HESWRD) go to 140
         go to 120
  130    format(1X,A7)
  140    do 160 I=1,NVAR
  150       format(A7)
            read(MION,170,end=160,ERR=630) (HESS(I,J),J=1,NVAR)
  160    continue
  170    format(8F10.2)
      end if
  180 if ((IGTHES.EQ.1).OR.(IGTHES.EQ.4)) then
C       IF IGTHES IS .EQ. 4, THEN THIS IS A HESSIAN RESTART.
C       USE GNEXT1 AND DUMMY FOR callS TO COMPFG DURING HESSIAN
C       CALCULATION
         if (IGTHES.EQ.1) write(MIOB,190)
  190    format(/,10X,'HESSIAN CALCULATED NUMERICALLY',/)
         TSTORE=TIME1
         do 210 I=ILOOP,NVAR
            XPARAM(I)=XPARAM(I) + XINC
            call COMPFG(XPARAM, .TRUE., DUMMY, .TRUE., GNEXT1, .TRUE.)
            if (IPRNT.GE.5) write(MIOB,'(I3,12(8F9.4,/3X))')
     1                                        I,(GNEXT1(IF),IF=1,NVAR)
            XPARAM(I)=XPARAM(I) - XINC
            do 200 J=1,NVAR
  200       HESS(I,J)= (GNEXT1(J)-GRAD(J))/XINC
            TIME2=SECOND()
            TSTEP=TIME2-TIME1
            TLEFT=TLEFT-TSTEP
            TIME1=TIME2
            if ( TLEFT .LT. TSTEP*TWO) then
C
C  STORE PARTIAL HESSIAN PATRIX
C  STORE GRADIENTS FOR GEOMETRY AND ILOOP AS POSITIVE
               write(MIOB,'(A)') ' NOT ENOUGH TIME TO COMPLETE HESSIAN'
               write(MIOB,'(A,I4)')
     %                          ' stopPING IN HESSIAN AT COORDINATE:',I
               IPOW(9)=1
               TT0=SECOND()-TIME0
               call BKRSAV(TT0,HESS,FUNCT,GRAD,XPARAM,PMAT,I,NSTEP,BMAT,
     1IPOW)
               call filclo  ! J/M
               stop ' in EF -  NOT ENOUGH TIME TO COMPLETE HESSIAN'
            end if
  210    continue
         TIME1=TSTORE
      end if
C
C     SYMMETRIZE HESSIAN
      do 220 I=1,NVAR
         do 220 J=1,I-1
            HESS(I,J)=(HESS(I,J)+HESS(J,I))/TWO
            HESS(J,I)=HESS(I,J)
  220 continue
  230 continue
C     START OF MAIN LOOP
C     WE NOW HAVE GRADIENTS AND A HESSIAN. IF THIS IS THE FIRST
C     TIME THROUGH DON'T UPDATE THE HESSIAN. FOR LATER LOOPS ALSO
C     CHECK IF WE NEED TO RECALCULATE THE HESSIAN
      TIME2=SECOND()
      TSTEP=TIME2-TIME1
      TLEFT=TLEFT-TSTEP
      TIME1=TIME2
      IFLEPO=0
  240 continue
      if (IHESS.GE.IRECLC.AND.IFLEPO.NE.15) then
         ILOOP=1
         IHESS=0
         IGTHES=1
         go to 180
      else
         if (IFLEPO.EQ.15) go to 270
         if (IPRNT.GE.2) then
            write(MIOB,'('' XPARAM '')')
            write(MIOB,'(5(2I3,F10.4))')
     %                          (LOC(1,I),LOC(2,I),XPARAM(I),I=1,NVAR)
            write(MIOB,'('' GRADIENTS'')')
            write(MIOB,'(3X,8F9.3)')                (GRAD(I),I=1,NVAR)
         end if
C
C        PRINT RESULTS IN CYCLE
         GNFINA=SQRT(DOT(GRAD,GRAD,NVAR))
         if (LDUMP.EQ.0) then
            write(MIOB,250)            NSTEP+1,MIN(TSTEP,9999.99D0),
     %           MIN(TLEFT,9999999.9D0),MIN(GNFINA,999999.999D0),FUNCT
            write(*,250)            NSTEP+1,MIN(TSTEP,9999.99D0),
     %           MIN(TLEFT,9999999.9D0),MIN(GNFINA,999999.999D0),FUNCT
            if (LOG) write(MIOE,250)   NSTEP+1,MIN(TSTEP,9999.99D0),
     %           MIN(TLEFT,9999999.9D0),MIN(GNFINA,999999.999D0),FUNCT
  250       format(' CYCLE:',I4,' TIME:',F7.2,' TIME LEFT:',F9.1,
     1             ' GRAD.:',F10.3,' HEAT:',G13.7)
         else
            write(MIOB,260)          MIN(TLEFT,9999999.9D0),
     %                       MIN(GNFINA,999999.999D0),FUNCT
            if (LOG) write(MIOE,260) MIN(TLEFT,9999999.9D0),
     %                       MIN(GNFINA,999999.999D0),FUNCT
  260       format(' RESTART FILE WRITTEN,   TIME LEFT:',F9.1,
     1             ' GRAD.:',F10.3,' HEAT:',G13.7)
         end if
C
C        TEST FOR CONVERGENCE
C
         RMX=SQRT(DOT(GRAD,GRAD,NVAR))
         if (RMX.LT.TOL2) go to 580
  270    continue
         IHESS=IHESS+1
         NSTEP=NSTEP+1
         if (IHESS.GT.1) call UPDHES(SVEC,TVEC,NVAR,IUPD)
      end if
C
      NHSINC=0
  280 IJ=0
      do 290 I=1,NVAR
         do 290 J=1,I
            IJ=IJ+1
            HESSC(IJ)=HESS(J,I)
  290 continue
      call HQRII(HESSC,NVAR,NVAR,EIGVAL,UC)
      IJ=0
      do 300 I=1,NVAR
         do 300 J=1,NVAR
            IJ=IJ+1
            U(J,I)=UC(IJ)
  300 continue
      if (IPRNT.GE.3) then
         if (IPRNT.GE.4) then
            write(MIOB,*) ' '
            write(MIOB,*) '              HESSIAN MATRIX'
            LOW=1
            NUP=8
  310       NUP=MIN(NUP,NVAR)
            write(MIOB,330) (I,I=LOW,NUP)
            do 320 I=1,NVAR
               write(MIOB,340) I,(HESS(I,J),J=LOW,NUP)
  320       continue
            NUP=NUP+8
            LOW=LOW+8
            if (LOW.LT.NVAR) go to 310
  330       format(/,3X,8I9)
  340       format(1X,I3,8F9.1)
  350       format(/,4X,8F9.1,/)
  360       format(1X,I3,8F9.4)
         end if
         write(MIOB,*) ' '
         write(MIOB,*) '              HESSIAN EIGENVALUES AND -VECTORS'
         LOW=1
         NUP=8
  370    NUP=MIN(NUP,NVAR)
         write(MIOB,330) (I,I=LOW,NUP)
         write(MIOB,350) (EIGVAL(I),I=LOW,NUP)
         do 380 I=1,NVAR
            write(MIOB,360) I,(U(I,J),J=LOW,NUP)
  380    continue
         NUP=NUP+8
         LOW=LOW+8
         if (LOW.LT.NVAR) go to 370
      end if
      do 390 I=1,NVAR
         if (ABS(EIGVAL(I)).LT.EIGMIN) then
            if (EIGVAL(I).LT.ZZERO) EIGVAL(I)=-EIGMIN
            if (EIGVAL(I).GT.ZZERO) EIGVAL(I)=EIGMIN
            write(MIOB,400) I,EIGVAL(I)
         end if
  390 continue
  400 format(5X,'WARNING! EIGENVALUE',I3,' TOO SMALL, REPLACED BY',
     1F7.3)
      NEG=0
      do 410 I=1,NVAR
         if (EIGVAL(I) .LT. ZZERO) NEG=NEG+1
  410 continue
      if (RMX.LT.TOL2) go to 580
      if (IHESS.EQ.1.AND.IPRNT.EQ.2) write(MIOB,420) NEG
  420 format(/,10X,'HESSIAN HAS',I3,' NEGATIVE EIGENVALUE(S)',/)
      do 430 I=1,NVAR
         FX(I)=DOT(U(1,I),GRAD,NVAR)
  430 continue
C  TAKE THE P-RFO STEP FOR A TS SEARCH AND
C  THE SIMPLE RFO STEP FOR A MINIMUM SEARCH
C
      if (IPRNT.GE.5.AND.NEG.EQ.NEGREQ.AND.NEGREQ.EQ.1) write(MIOB,690)
      if (IPRNT.GE.5.AND.NEG.EQ.NEGREQ.AND.NEGREQ.EQ.0) write(MIOB,700)
      if (IPRNT.GE.5.AND.NEG.NE.NEGREQ.AND.NEGREQ.EQ.1) write(MIOB,710)
      if (IPRNT.GE.5.AND.NEG.NE.NEGREQ.AND.NEGREQ.EQ.0) write(MIOB,720)
  440 call FORMD(EIGVAL,FX,NVAR,FAIL,ORIDE)
c
C     IF FAIL IS TRUE THEN NO ACCEPTABLE LAMDA VALUE COULD BE FOUND.
C     IF GNORM IS .LT. RCUT AND NR .EQ. TRUE, THEN TAKE SIMPLE NR-STEP.
C     IF GNORM IS .GT. RCUT OR NR .EQ. FALSE, THEN MODIFY HESSIAN.
C     IF HSSHFT IS ON, THEN INCREASE HESSIAN DIAGONAL ELEMENTS WITH
C     EIGVAL(1) + 2*EIGMIN FOR MIN SEARCH, OR EIGVAL(1+MODE) + 2*EIGMIN
C     FOR TS SEARCH. THEN CALL FORMD AGAIN.
C     IF EIGINV IS ON, THEN REVERSE THE SIGN OF ALL NEGATIVE EIGENVALUES
C     IN THE HESSIAN OTHER THAN THE TS MODE BEING FOLLOWED, AND call
C     FORMD AGAIN. IF NO ACCEPTABLE LAMDA CAN BE FOUND AFTER HESSIAN
C     HAS BEEN MODIFIED, WE ARE OUT OF LUCK...
C
      if (FAIL) then
         if (IFLEPO.EQ.7) return
         if (NR.AND.(RMX.LT.RCUT)) then
            write(MIOB,450) RCUT
  450       format(5X,'MAX GRADIENT IS LESS THAN RCUT =',F5.1,
     1             ' TAKING SIMPLE NR-STEP')
            do 460 I=1,NVAR
               D(I)=ZZERO
  460       continue
            do 480 I=1,NVAR
               TEMP=-FX(I)/EIGVAL(I)
               do 470 J=1,NVAR
                  D(J)=D(J)+TEMP*U(J,I)
  470          continue
  480       continue
            FAIL=.FALSE.
            go to 530
         end if
         if (HSSHFT) then
            NHSINC=NHSINC+1
            if (NHSINC.GT.1) then
               write(MIOB,640)
               write(MIOB,660)
               go to 620
            end if
            if (IMIN.NE.0) TEMP=TWO*EIGMIN+ABS(EIGVAL(1))
            if (ITS .NE.0) then
               if (MODE.EQ.0) TEMP=TWO*EIGMIN+ABS(EIGVAL(2))
               if (MODE.NE.0) TEMP=TWO*EIGMIN+ABS(EIGVAL(MODE+1))
            end if
            do 490 I=1,NVAR
               HESS(I,I)=HESS(I,I)+TEMP
               EIGVAL(I)=EIGVAL(I)+TEMP
  490       continue
            write(MIOB,500) TEMP
  500       format(5X,'HESSIAN DIAGONAL INCREASED BY',F8.1)
            FAIL=.FALSE.
            go to 440
         end if
         if (EIGINV) then
            NNEG=0
            do 520 I=1,NVAR
               if ((MODE.EQ.0).AND.(I.EQ.1)) go to 520
               if ((MODE.NE.0).AND.(I.EQ.MODE)) go to 520
               if (EIGVAL(I).LT.ZZERO) then
                  NNEG=NNEG+1
                  EIGVAL(I)=-EIGVAL(I)
                  write(MIOB,510) EIGVAL(I)
  510             format(5X,'HESSIAN EIGENVALUE REVERSED, NEW VALUE =',F
     18.1)
               end if
  520       continue
            if (NNEG.GT.0) then
               FAIL=.FALSE.
               go to 440
            else
               write(MIOB,650)
               write(MIOB,660)
               go to 620
            end if
         end if
      end if
C
C  WE NOW HAVE A NEW STEP IN D
C  CHECK THAT THE STEPSIZE DOES NOT EXCEED DMAX
C  IF SO, SCALE
C
  530 SKAL=ONE
      DD=DOT(D,D,NVAR)
      DD=SQRT(DD)
C
      if (DD.GT.DMAX) then
         SKAL=DMAX/DD
         do 540 I=1,NVAR
            D(I)=D(I)*SKAL
  540    continue
         DD=DMAX
         if (IPRNT.GE.2) write(MIOB,550) SKAL
  550    format(5X,'CALCULATED STEP SIZE TOO LARGE, SCALED WITH',F9.5)
      end if
C
      if (IPRNT.GE.2) write(MIOB,560) DD
  560 format(5X,'STEPSIZE USED IS',F9.5)
      if (IPRNT.GE.2) then
         write(MIOB,'('' CALCULATED STEP'')')
         write(MIOB,'(3X,8F9.5)') (D(I),I=1,NVAR)
      end if
C
C  FORM NEW XPARAM AND GEO, SAVE CURRENT GRADIENTS
C
      do 570 I=1,NVAR
         XPARAM(I)=XPARAM(I)+D(I)
         OLDF(I)=GRAD(I)
         K=LOC(1,I)
         L=LOC(2,I)
         GEO(L,K)=XPARAM(I)
  570 continue
      if (NDEP.NE.0) call SYMTRY
C     CHECK STEPS AND ENOUGH TIME FOR ANOTHER PASS
      TIME2=SECOND()
      TSTEP=TIME2-TIME1
      TLEFT=TLEFT-TSTEP
      TIME1=TIME2
      if (TSTEP.LT.ZZERO) TSTEP=ZZERO
      if (TLEFT .LT. TSTEP*TWO) go to 620
C     IN USER UNFRIendLY ENVIROMENT, SAVE RESULTS EVERY 1 CPU HRS
      ITTEST=AINT((TIME2-TIME0)/TDUMP)
      if (ITTEST.GT.NTIME) then
         LDUMP=1
         NTIME=MAX(ITTEST,(NTIME+1))
         IPOW(9)=2
         TT0=SECOND()-TIME0
         call BKRSAV(TT0,HESS,FUNCT,GRAD,XPARAM,PMAT,-NSTEP,NSTEP,BMAT,I
     1POW)
      else
         LDUMP=0
      end if
C     GET GRADIENT FOR NEW GEOMETRY AND return FOR ANOTHER CYCLE
      call COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .TRUE.)
      go to 240
  580 continue
      ifLEPO=15
C
C  CHECK: HAS EIGVAL BEEN SET, if NOT then GO BACK AND SET IT.
C
      if (DOT(EIGVAL,EIGVAL,NVAR).LT.1.D-15) go to 240
      if (NEG.NE.NEGREQ) then
         IFLEPO=16
         write(MIOB,590)
         write(MIOB,600) (EIGVAL(I),I=1,5)
  590    format(/,5X,'WARNING! HESSIAN DOES NOT HAVE THE REQUIRED ',
     1'STRUCTURE',/)
  600    format(5X,'LOWEST 5 EIGENVALUES OF HESSIAN',3X,5F8.2)
      end if
      LAST=1
      if (IPRNT.GE.2) write(MIOB,610) RMX,TOL2
  610 format(/,5X,'MAXIMUM GRADIENT =',F9.5,'  IS LESS THAN CUTOFF =',
     1F9.5,//)
C     call COMPFG TO CALCULATE ENERGY FOR FIXING MO-VECTOR BUG
      call COMPFG(XPARAM, .TRUE., FUNCT, .TRUE., GRAD, .FALSE.)
      return
  620 continue
C     WE RAN OUT OF TIME, OR FORMD FAIL. DUMP RESULTS
      if (TLEFT .LT. TSTEP*TWO) then
         write(MIOB,680)
      end if
      IPOW(9)=1
      TT0=SECOND()-TIME0
      call BKRSAV(TT0,HESS,FUNCT,GRAD,XPARAM,PMAT,-NSTEP,NSTEP,BMAT,I
     1POW)
      call filclo  !  J/M
      stop ' in EF - not enough time !'
  630 continue
      write(MIOB,*)'     ERROR DURING READ OF HESSIAN FROM INPUT'
      call filclo  !  J/M
      stop ' in EF -  ERROR DURING READ OF HESSIAN FROM INPUT'
  640 format(5X,'HESSIAN DIAGONAL INCREASED TWICE WITHOUT STEP',/,5X,
     1'STILL NO ACCEPTABLE LAMDA VALUE. TERMINATING')
  650 format(5X,'ALL NON-TS MODES HAVE POSITIVE EIGENVALUES',/,5X,
     1'YET FORMD CANNOT FIND ACCEPTABLE LAMDA. TERMINATING')
  660 format(5X,'EF OPTIMIZATION FAILED.... SORRY....')
  670 format(/,5X,'EXCESS NUMBER OF OPTIMIZATION CYCLES')
  680 format(/,5X,'NOT ENOUGH TIME FOR ANOTHER CYCLE')
  690 format(5X,'TS SEARCH. TAKING P-RFO STEP')
  700 format(5X,'MINIMUM SEARCH. TAKING SIMPLE RFO STEP')
  710 format(5X,'HESSIAN DOES NOT HAVE THE DESIRED LOCAL STRUCTURE'/
     1       5X,'TAKING P-RFO STEP')
  720 format(5X,'HESSIAN DOES NOT HAVE THE DESIRED LOCAL STRUCTURE'/
     1       5X,'TAKING SIMPLE RFO STEP')
      end
