C 14 NOV 05 - DGF - PAD COMMON BLOCK SCFOPT
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - DGF - GRXINI: SANITIZE A COUPLE OF CONSTANTS
C 30 APR 05 - DGF - SYNCHRONISE OPTGRD COMMON BLOCK
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK FRGINF
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 17 APR 02 - MWS - SYNCH UP FRGINF COMMON
C 26 MAR 02 - KRG - END JOB PROPERLY (EXCEPT IT SHOULD RETURN TO MAIN!)
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C 25 JUN 01 - MWS - ALTER COMMON BLOCK SCFOPT
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C  6 JAN 98 - JAB - ADD A VARIABLE TO ZMAT COMMON BLOCK
C 18 DEC 96 - FRJ - GRXSR: TWEAK NUMERICAL CONTROLS BASED ON EXPERIENCE
C 17 OCT 96 - SPW - GRXSR: ADD CALL TO CIGRAD FOR CI RUNS
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C 20 SEP 95 - FRJ - CODED UP GRADIENT EXTREMAL ROUTINES
C
C*MODULE GRADEX  *DECK GRADEX
      SUBROUTINE GRADEX
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DOUBLE PRECISION METHOD
C
      LOGICAL OUT,GOPARR,DSKWRK,MASWRK,STPT
      LOGICAL GOFRST,HSDFDB
C
      PARAMETER (MXATM=2000)
C
      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      DATA CHECK /8HCHECK   /, DEBUG/8HDEBUG   /, DBUGME/8HGRADEX  /
      DATA GEJJH /8HJJH     /, GESR /8HSR      /
C
C     ----- GRADIENT EXTREMUM FOLLOWING -----
C
C     CURRENT VERSION IMPLEMENTED IN GAMESS
C     SEPTEMBER 1993 BY F.JENSEN, ODENSE, DK
C
C     ----- INITIALIZE THE STATIONARY POINT OPTIMIZATION -----
C
      OUT = NPRINT.EQ.-6 .AND. MASWRK
      IF(MASWRK  .AND.  EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)
     *        OUT=.TRUE.
C
      NSERCH = 0
      NPUNCH=1
C
      CALL GRXSET(NCVAL,NCOORD,METHOD,OUT,
     *                  SNUMH,STPT,STSTEP,GOFRST,DELCOR,SNRMAX,
     *                  MYSTEP,HSDFDB)
C
      IF(EXETYP.EQ.CHECK) RETURN
      NPRINT=-5
C
      IF(STPT) CALL GRXINM(NCOORD,NCVAL,NCVAL,STSTEP,OUT)
C
      IF(METHOD.EQ.GEJJH) CALL GRXJJM(NCVAL,NCOORD,STPT,GOFRST,OUT,
     *                     SNRMAX)
      IF(METHOD.EQ.GESR)  CALL GRXSRM(NCOORD,STPT,GOFRST,OUT,
     *                     SNUMH,DELCOR,STSTEP,SNRMAX,MYSTEP,HSDFDB)
C
      RETURN
C
      END
C*MODULE GRADEX  *DECK GRXINM
      SUBROUTINE GRXINM(NCOORD,NCVAL,NCVF,STSTEP,OUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL OUT
C
      COMMON /FMCOM / XX(1)
C
C     ----- GROW MEMORY FOR GRXINI -----
C       -FC-      AT I10
C       -A-       AT I10
C       -VEC-     AT I20
C       -SCR-     AT I30
C       -EIGEN-   AT I40
C       -IA-      AT I50
C       -BUF1-    AT I60
C       -COMC-    AT I70
C       -RM-      AT I80
C
      NC1 = NCVF
      NC3 = NC1*NC1
      LOADFM = 0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I20 = I10 + NC3
      I30 = I20 + NC3
      I40 = I30 + 8*NC1
      I50 = I40 + NC1
      I60 = I50 + NC1
      I70 = I60 + NC3
      I80 = I70 + NC1
      I90 = I80 + NC1
      LAST = I90
      NEED = LAST - I10
      CALL GETFM(NEED)
C
      CALL GRXINI(NCOORD,NCVAL,NCVF,STSTEP,OUT,XX(I10),XX(I10),
     *           XX(I20),XX(I30),XX(I40),XX(I50),XX(I60),
     *           XX(I70),XX(I80))
C
C     ----- RETURN FAST MEMORY -----
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE GRADEX  *DECK GRXINI
      SUBROUTINE GRXINI(NCOORD,NCVAL,NCVF,STSTEP,OUT,
     *                  FC,A,VEC,SCR,EIGEN,IA,BUF1,COMC,RM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION VEC(NCVF,NCVF),A(NCVF,NCVF),FC(1),
     *          SCR(8,NCVF),EIGEN(NCVF),IA(NCVF),BUF1(NCVF,NCVF),
     *          COMC(NCVF),RM(NCVF)
      DIMENSION CMASS(3)
C
      LOGICAL OUT,LINEAR
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
      PARAMETER (MXATM=2000, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FUNCT / ENERGY,F(3*MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTEF / D(3*MXATM+6*MXFRG),OLDF(3*MXATM+6*MXFRG),
     *                VMODE(3*MXATM+6*MXFRG),RADIUS
      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGF
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
C
C     DIAGONALIZE INITIAL HESSIAN AND SELECT INITIAL STEP DIRECTION
C     AS IFOLOW EIGENVECTOR. STUFF THIS INTO VMODE AND CALCULATE
C     INITIAL POINT AT STSTEP ALONG THIS DIRECTION
C
C     ----- PROJECT OUT T,R MODES -----
C     USE UNIT MASSES TO PROJECT OUT DIRECTLY IN CARTESIANS
C
      LRM=3*NAT+21*NFRG
      LFRG=6*NFRG
      CALL VCLR(RM,1,LRM)
      CALL VCLR(FRGMI,1,LFRG)
      I0=1
      DO 10 I=1,3*(NAT+NFRG)
         RM(I0  ) = ONE
         I0=I0+1
   10 CONTINUE
      I0=3*(NAT+NFRG)+1
      I1=3*(NAT+4*NFRG)+1
      DO 20 I=1,NFRG
         RM(I0  ) = ONE
         RM(I0+4) = ONE
         RM(I0+8) = ONE
         RM(I1  ) = ONE
         RM(I1+4) = ONE
         RM(I1+8) = ONE
         FRGMI(1,I) = ONE
         FRGMI(3,I) = ONE
         FRGMI(6,I) = ONE
         I0=I0+9
         I1=I1+9
   20 CONTINUE
      ZMASST=NAT+NFRG
C
      CALL VCLR(CMASS,1,3)
      DO 40 I=1,NAT
         DO 30 J=1,3
            CMASS(J)=CMASS(J)+C(J,I)
   30    CONTINUE
   40 CONTINUE
      DO 45 I=1,NFRG
         DO 44 J=1,3
            CMASS(J)=CMASS(J)+EFCENT(J,I)
   44    CONTINUE
   45 CONTINUE
      DO 50 I=1,3
         CMASS(I)=CMASS(I)/ZMASST
   50 CONTINUE
C
      KK=0
      DO 70 I=1,NAT
         DO 60 J=1,3
            KK=KK+1
            COMC(KK) = C(J,I)-CMASS(J)
   60    CONTINUE
   70 CONTINUE
      DO 90 I=1,NFRG
         DO 80 J=1,3
            KK=KK+1
            COMC(KK) = EFCENT(J,I)-CMASS(J)
   80    CONTINUE
   90 CONTINUE
C
      CALL DAREAD(IDAF,IODA,FC,NCVF*NCVF,4,0)
C     IF(.NOT.OUT)MASWRK=.FALSE.
      CALL PRJFC(.FALSE.,.TRUE.,ZMASST,FC,D,COMC,RM,FRGMI,
     *              BUF1,VEC,NAT,NCVAL,NCVF)
      IJ = 0
      DO 150 I = 1,NCOORD
         DO 140 J = 1,I
            IJ = IJ + 1
            FC(IJ) = A(J,I)
  140    CONTINUE
  150 CONTINUE
C
C     ----- OBTAIN CURRENT NORMAL MODE DIRECTIONS -----
C
      IF(OUT) WRITE(IW,9000)
      IF(OUT) CALL PRTRI(FC,NCOORD)
      IERR = 0
      CALL GLDIAG(NCVF,NCVF,NCVF,FC,SCR,EIGEN,VEC,IERR,IA)
C     SHIFT T, R MODES TO BE THE LAST 6(5)
      NTR=0
      DO 151 I=NCVF,1,-1
         IF(ABS(EIGEN(I)).LT.1.0D-09) NTR=NTR+1
         IF(ABS(EIGEN(I)).LT.1.0D-09) ISHIFT=I
  151 CONTINUE
      IF (NTR.NE.5 .AND. NTR.NE.6) THEN
         WRITE(IW,*) 'ERROR! NOT 5 OR 6 TR MODES',(EIGEN(I),I=1,NCVF)
         CALL ABRT
      END IF
      JSHIFT=ISHIFT-1
      DO 153 J=1,NTR
         JSHIFT=JSHIFT+1
         SCR(8,J)=EIGEN(JSHIFT)
         DO 152 I=1,NCVF
            SCR(J,I)=VEC(I,JSHIFT)
  152    CONTINUE
  153 CONTINUE
      DO 155 J=ISHIFT,NCVF-NTR
         EIGEN(J)=EIGEN(J+NTR)
         DO 154 I=1,NCVF
            VEC(I,J)=VEC(I,J+NTR)
  154    CONTINUE
  155 CONTINUE
      JTMP=0
      DO 157 J=NCVF-NTR+1,NCVF
         JTMP=JTMP+1
         EIGEN(J)=SCR(8,JTMP)
         DO 156 I=1,NCVF
            VEC(I,J)=SCR(JTMP,I)
  156    CONTINUE
  157 CONTINUE
      IF(OUT) WRITE(IW,9010)
      IF(OUT) CALL PREVNL(VEC,EIGEN,NCVF,NCVF,NCVF)
C
C     ----- PRINT (UP TO 10) INITIAL MODES  -----
C
      WRITE(IW,9012)
      KMAX = MIN(2,1 + (NCVF-1)/5)
      DO 170 K=1,KMAX
         MINJ = 5*(K-1) + 1
         MAXJ = MIN(K*5,NCVF)
         WRITE(IW,9013)
         WRITE(IW,9014) (I,I=MINJ,MAXJ)
         WRITE(IW,9015) (EIGEN(I),I=MINJ,MAXJ)
         WRITE(IW,9014)
         DO 160 I=1,NCVF
            WRITE(IW,9016) I,(VEC(I,J),J=MINJ,MAXJ)
  160    CONTINUE
  170 CONTINUE
C
      GNORM=SQRT(DDOT(NCVF,F,1,F,1))
      IF (GNORM.GT.CONVF) WRITE(IW,9019)GNORM,CONVF
      ITMP=ABS(IFOLOW)
      WRITE(IW,9020)STSTEP,IFOLOW,EIGEN(ITMP)
      TEST=ZERO
      DO 175 I=1,NCVF
         VMODE(I)=VEC(I,ITMP)
         IF (ABS(VMODE(I)).GT.(ABS(TEST)+1.0D-06))TEST=VMODE(I)
  175 CONTINUE
      IF(TEST.LT.ZERO) CALL DSCAL(NCVF,-ONE,VMODE,1)
      IF(IFOLOW.LT.ZERO) CALL DSCAL(NCVF,-ONE,VMODE,1)
      WRITE(IW,9021) (VMODE(I),I=1,NCVF)
      IF (ABS(EIGEN(ITMP)).LT.1.0D-09) THEN
         WRITE(IW,9022)
         CALL ABRT
      END IF
      CALL DCOPY(NCVF,VMODE,1,D,1)
      CALL DSCAL(NCVF,STSTEP,D,1)
      RADIUS=STSTEP
      IFOLOW=ABS(IFOLOW)
C
      IF(NZMAT.GT.0) THEN
         CALL DCOPY(NCVF,VMODE,1,EIGEN,1)
         CALL TRANG(EIGEN,NZVAR,NCOORD)
         WRITE(IW,9023)
         WRITE(IW,9021) (EIGEN(I),I=1,NZVAR)
      END IF
C
      CALL DCOPY(NCVF,C,1,X,1)
      CALL SYMDR(D)
      CALL UPCOOR(D,NCOORD,OUT)
C
 9000 FORMAT(16X,'FORCE CONSTANT MATRIX')
 9010 FORMAT(/10X,39HFORCE CONSTANT EIGENVALUES AND VECTORS ,/)
 9012 FORMAT(10X,'THE FIRST FEW LOCAL MODES AND EIGENVALUES ARE')
 9013 FORMAT(1X)
 9014 FORMAT(1X,' MODE',5I12)
 9015 FORMAT(6X,5F12.6)
 9016 FORMAT(1X,I3,2X,5F12.6)
 9019 FORMAT(/,10X,'WARNING! INITIAL GRADIENT NORM =',F10.6,/
     *       30X,'BUT OPTTOL =',F10.6,/,
     *       10X,'ARE YOU SURE YOUR START GEOMETRY IS A ',
     *           'STATIONARY POINT?')
 9020 FORMAT(/,10X,'INPUT GEOMETRY TAKEN AS STATIONARY',/,
     *       10X,'GENERATING INITIAL GEOMETRY BY TAKING A STEP OF',
     *       F8.5,1X,'ALONG HESSIAN MODE',I4,/,
     *       10X,'WITH EIGENVALUE',F10.6,1X,'AND CARTESIAN COMPONENTS')
 9021 FORMAT(9F10.5)
 9022 FORMAT(/,10X,'TRYING TO FOLLOW MODE WITH EIGENVALUE = 0',/,
     *       10X,'MAYBE THIS IS A TRANSLATIONAL OR ROTATIONAL MODE?')
 9023 FORMAT(/,10X,'INTERNAL COMPONENTS')
      RETURN
      END
C*MODULE GRADEX  *DECK GRXJJM
      SUBROUTINE GRXJJM(NCVAL,NCOORD,STPT,GOFRST,OUT,SNRMAX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL STPT,GOFRST,OUT
C
      COMMON /FMCOM / X(1)
C
C     ----- GROW MEMORY FOR GRXJJH -----
C       -RM-      AT I10
C       -EIGEN-   AT I20
C       -TEMP-    AT I30
C       -WORK-    AT I35
C       -IA-      AT I40
C       -SCR-     AT I50
C       -HESSL-   AT I55
C       -HESS-    AT I60
C       -UMAT-    AT I70
C       -XMAT-    AT I80
C       -YMAT-    AT I90
C
      NC1 = NCOORD
      NC3 = NCOORD*NCOORD
      LOADFM = 0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I20 = I10 + NC1
      I30 = I20 + NC1
      I35 = I30 + NC1
      I40 = I35 + NC1
      I50 = I40 + NC1
      I55 = I50 + 8*NC1
      I60 = I55 + NC3
      I70 = I60 + NC3
      I80 = I70 + NC3
      I90 = I80 + NC3
      I100 = I90 + NC3
      LAST = I100
      NEED = LAST - I10
      CALL GETFM(NEED)
C
      CALL GRXJJH(NCVAL,NC1,STPT,GOFRST,OUT,SNRMAX,
     *           X(I10),X(I20),X(I30),X(I35),
     *           X(I40),X(I50),X(I55),X(I60),X(I70),X(I80),X(I90))
C
C     ----- RETURN FAST MEMORY -----
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE GRADEX  *DECK GRXJJH
      SUBROUTINE GRXJJH(NCVAL,NC1,STPT,GOFRST,OUT,SNRMAX,
     *                 RM,EIGEN,TEMP,WORK,
     *                 IA,SCR,HESSL,HESS,UMAT,XMAT,YMAT)
C
C     THIS REPRESENT AN IMPLEMENTATION OF THE JORGENSEN,
C     JENSEN AND  HELGAKER ALGORITHM FOR FOLLOWING GRADIENT EXTREMALS.
C     P. JORGENSEN, H. J. AA. JENSEN, T. HELGAKER
C     THEO. CHIM. ACTA 73, 55 (1988).
C
C     IMPLEMENTED IN GAMESS SEPT. 1995
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION TIMSTR(3)
C
      LOGICAL OUT,GOPARR,DSKWRK,MASWRK,LINEAR
      LOGICAL STPT,GOFRST
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, PT5=0.5D+00, TWO=2.0D+00,
     *           THREE=3.0D+00, PT75=0.75D+00, FOUR=4.0D+00,
     *           TMONE=1.0D-01, TMSIX=1.0D-06)
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION CMASS(3)
      DIMENSION RM(NC1),EIGEN(NC1),TEMP(NC1),IA(NC1),WORK(NC1)
      DIMENSION HESSL(NC1*NC1),HESS(NC1,NC1),SCR(8,NC1)
      DIMENSION UMAT(NC1,NC1),XMAT(NC1,NC1),YMAT(NC1,NC1)
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3*MXFRG),TORQ(3*MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FUNCT / ENERGY,EG(3*MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /OPTEF / D(3*MXATM+6*MXFRG),OLDF(3*MXATM+6*MXFRG),
     *                VMODE(3*MXATM+6*MXFRG),RADIUS
      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGF
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
      SAVE DEPRE
C
      DATA DEPRE/0.0D+00/
C
C     ----- GRADIENT EXTREMUM FOLLOWING USING THE ALGORITHM
C     OF JORGENSEN, JENSEN AND HELGAKAR (THEO.CHIM.ACTA 73, 55) -----
C
C     CURRENT VERSION IMPLEMENTED IN GAMESS
C     NOVEMBER 1993 BY F.JENSEN, ODENSE, DK
C
C     NC1 IS = NCOORD
      TODEG=180.0D+00/ACOS(-ONE)
C
      DO 10 I=1,NC1
         ZMASS(I)=ONE
         RM(I)=ONE
   10 CONTINUE
      OLDE=ZERO
C
 777  CONTINUE
      NSERCH=NSERCH+1
      WRITE(IW,9100)
      CALL PRATM(IW,2)
      IF(NZMAT.GT.0) CALL PZANDG(EG,0)
C     CALCULATE GRADIENT AND HESSIAN
      WRITE(IW,9110)
C     CALL DCOPY(NC1,X,1,C,1)
      CALL HESSX(.FALSE.,.FALSE.)
      WRITE(IW,9120)
      IF(NSERCH.GT.1) THEN
         DEACT = ENERGY-OLDE
         RATIO = DEACT/DEPRE
         IF(MASWRK) WRITE(IW,9060) DEACT,DEPRE,RATIO
C        CURRENTLY TRUST UPDATE IS INACTIVE
         IF(ITRUPD.EQ.1) THEN
C
C           IF THIS WAS A GOOD STEP, INCREASE THE ALLOWED STEP LENGTH
C           REDUCE IT IF IT WAS A POOR STEP
C
            IF (RATIO.LE.PT5 .OR. RATIO.GE.TWO)
     *          DXMAXT=DXMAXT/TWO
            IF (RATIO.GE.PT75 .AND. RATIO.LE.(FOUR/THREE))
     *          DXMAXT=DXMAXT*SQRT(TWO)
C
C           BE BRAVE IF  0.90 < RATIO < 1.10
C
            IF (ABS(RATIO-ONE).LT.TMONE) DXMAXT=DXMAXT*SQRT(TWO)
            DXMAXT=MAX(DXMAXT,TRMIN)
            DXMAXT=MIN(DXMAXT,TRMAX)
         END IF
      END IF
      OLDE=ENERGY
 9060 FORMAT(/,10X,'   ACTUAL ENERGY CHANGE WAS',F15.10/
     *       10X,'PREDICTED ENERGY CHANGE WAS',F15.10,' RATIO=',F7.3)
      WRITE (IW,9001) NSERCH,ENERGY
      WRITE (IW,9011)
      DO 105 IAT = 1,NAT
         N = 3*(IAT-1)
         ZNUC = ZAN(IAT) + IZCORE(IAT)
         WRITE (IW,9012) IAT,ANAM(IAT),BNAM(IAT),ZNUC,
     *                   (C(J,IAT),J=1,3), (EG(N+I),I=1,3)
  105 CONTINUE
      IF(NZMAT.GT.0) THEN
         CALL DCOPY(NC1,EG,1,TEMP,1)
         CALL TRANG(TEMP,NZVAR,NC1)
         CALL PZANDG(TEMP,1)
      END IF
      GNORM=SQRT(DDOT(NC1,EG,1,EG,1))
      IF (GNORM.LT.CONVF) THEN
         WRITE(IW,9021) NSERCH,ENERGY,GNORM
         WRITE(IW,9900) GNORM,CONVF
         CALL SIGEND(.TRUE.,NSERCH,NCVAL,NC1,0,0)
         RETURN
      END IF
C
 9900 FORMAT(1X,'GRADIENT NORM =',F15.10,' LESS THAN CUTOFF =',
     *          F15.10,/,1X,'GEOMETRY CONVERGED')
C
C     DIAGONALIZE HESSIAN AND CHECK FOR PROPER NUMBER OF NTR
C     AND NEGATIVE EIGENVALUES
C
      CALL DAREAD(IDAF,IODA,HESS,NC1*NC1,4,0)
      IF(OUT) WRITE(IW,*) 'HESS'
      IF(OUT) CALL PRSQ(HESS,NC1,NC1,NC1)
      CALL CENMAS(NAT,NFRG,C,TEMP,ZMASST,CMASS,ZMASS)
      CALL PRJFC(.FALSE.,.TRUE.,ZMASST,HESS,WORK,TEMP,RM,FRGMI,
     *              XMAT,YMAT,NAT,NC1,NC1)
      CALL CPYSQT(HESS,HESSL,NC1,1)
      CALL GLDIAG(NC1,NC1,NC1,HESSL,SCR,EIGEN,UMAT,IERR,IA)
      IF(OUT) WRITE(IW,*) 'HESSIAN EIGENVALUES',(EIGEN(I),I=1,NC1)
      NTR=0
      DO 270 I=NC1,1,-1
         IF (ABS(EIGEN(I)).LT.TMSIX) EIGEN(I)=ZERO
         IF (EIGEN(I).EQ.ZERO) NTR=NTR+1
         IF (EIGEN(I).EQ.ZERO) ISHIFT=I
  270 CONTINUE
      IF (NTR.NE.5 .AND. NTR.NE.6) THEN
         WRITE(IW,*) 'WARNING! NOT 5 OR 6 TR MODES',(EIGEN(I),I=1,NC1)
      END IF
C
C     SHIFT T,R MODES TO BE THE LAST NTR
C
      JSHIFT=ISHIFT-1
      DO 221 J=1,NTR
         JSHIFT=JSHIFT+1
         XMAT(9,J)=EIGEN(JSHIFT)
         DO 211 I=1,NC1
            XMAT(J,I)=UMAT(I,JSHIFT)
  211    CONTINUE
  221 CONTINUE
      DO 241 J=ISHIFT,NC1-NTR
         EIGEN(J)=EIGEN(J+NTR)
         DO 231 I=1,NC1
            UMAT(I,J)=UMAT(I,J+NTR)
  231    CONTINUE
  241 CONTINUE
      JTMP=0
      DO 261 J=NC1-NTR+1,NC1
         JTMP=JTMP+1
         EIGEN(J)=XMAT(9,JTMP)
         DO 251 I=1,NC1
            UMAT(I,J)=XMAT(JTMP,I)
  251    CONTINUE
  261 CONTINUE
C
      IF(.NOT.STPT) THEN
         WRITE(IW,*) 'START GEOMETRY IS NOT A STATIONARY POINT,'
         WRITE(IW,*) 'SELECTING CURRENT HESSIAN EIGENVECTOR AS'
         WRITE(IW,*) 'THE GRADIENT EXTREMAL DIRECTION.'
         KK=IFOLOW
         IF (ABS(EIGEN(ABS(KK))).EQ.ZERO) THEN
            WRITE(IW,*) 'INITIAL DIRECTION HAS ZERO EIGENVALUE'
            WRITE(IW,*) 'MAYBE THIS IS A T OR R MODE ???',KK
            CALL ABRT
         END IF
         IF(IFOLOW.LT.0)THEN
            IFOLOW=-IFOLOW
            DIR=-ONE
         ELSE
            DIR=ONE
         END IF
         DO 271 I=1,NC1
            VMODE(I)=DIR*UMAT(I,IFOLOW)
  271    CONTINUE
         STPT=.TRUE.
      END IF
      IF(OUT) WRITE(IW,*) 'HESSIAN EIGENVALUES',(EIGEN(I),I=1,NC1)
      IF(OUT) WRITE(IW,*) 'HESSIAN EIGENVECTORS'
      IF(OUT) CALL PRSQ(UMAT,NC1,NC1,NC1)
      TMPM=ZERO
      TMPP=ZERO
      IMODE=0
      IPRO =0
      WRITE(IW,9000)
      DO 292 I=1,NC1-NTR
         TMP1=ABS(DDOT(NC1,EG,1,UMAT(1,I),1))
         TMP1=TMP1/GNORM
         IF (ABS(TMP1).GT.ABS(TMPP)) IPRO=I
         IF (ABS(TMP1).GT.ABS(TMPP)) TMPP=TMP1
         TMP2=(DDOT(NC1,VMODE,1,UMAT(1,I),1))
         IF (ABS(TMP2).GT.ABS(TMPM)) IMODE=I
         IF (ABS(TMP2).GT.ABS(TMPM)) TMPM=TMP2
         WRITE(IW,9010)I,EIGEN(I),TMP1,ABS(TMP2)
  292 CONTINUE
 9000 FORMAT(/,1X,'PROJECTION OF GRADIENT ALONG HESSIAN EIGENVECTORS',
     *          ' AND MODE OVERLAPS',/,
     *      15X,'EIGENVALUE',5X,'PROJECTION',7X,'OVERLAP',/)
 9010 FORMAT(5X,I5,3F15.8)
      DIR=ONE
      IF (ABS(TMPM).GT.ONE)TMPM=TMPM/ABS(TMPM)
      IF (TMPM.LT.ZERO) DIR=-ONE
      IF (TMPM.LT.ZERO) TMPM=-TMPM
      ANGLE=TODEG*ACOS(TMPM)
      DO 291 I=1,NC1
         VMODE(I)=DIR*UMAT(I,IMODE)
         UMAT(I,IMODE)=VMODE(I)
  291 CONTINUE
      WRITE(IW,9015)IMODE,TMPM
      IF(OUT) WRITE(IW,9025) (VMODE(I),I=1,NC1)
      IF(NZMAT.GT.0 .AND. OUT) THEN
         CALL DCOPY(NC1,VMODE,1,TEMP,1)
         CALL TRANG(TEMP,NZVAR,NC1)
         WRITE(IW,9026) (TEMP(J),J=1,NZVAR)
      END IF
      IF(IMODE.NE.IFOLOW) WRITE(IW,9016)IFOLOW,IMODE
      IFOLOW=IMODE
      IF (IPRO.NE.IFOLOW) THEN
         WRITE(IW,9019)
C        CALL ABRT
      END IF
 9025 FORMAT(1X,'CARTESIAN COMPONENTS',5(/,8F10.5),/)
 9026 FORMAT(1X,'INTERNAL  COMPONENTS',5(/,8F10.5),/)
 9015 FORMAT(/,1X,'LARGEST OVERLAP WITH NEW HESSIAN IS MODE',I5,3X,
     *          'WITH OVERLAP',F8.5,/)
 9016 FORMAT(1X,'WARNING MODE SWITCHING, FROM',I5,3X,'TO',I5)
 9018 FORMAT(1X,'INSIDE NR REGION, STEPLENGTH =',F10.5,' TAKING',
     *          ' PURE NR STEP')
 9019 FORMAT(1X,'WARNING!!! ALGORITHM SEEMS TO BE OFF TRACK',/,
     *1X,'LARGEST GRADIENT COMPONENT IS NOT ALONG CORRECT EIGENVECTOR')
C
C     FORM PURE NR STEP PERPENDICULAR TO VMODE DIRECTION
C     IF AN EIGENVALUE HAS BEEN ZERO OUT THEN ZERO CORRESPONDING TEMP TO
C     ALLOW FORMATION OF STEP WITHOUT THESE CONTRIBUTIONS.
      DO 220 I=1,NC1
         TEMP(I)=DDOT(NC1,EG,1,UMAT(1,I),1)
         IF (EIGEN(I).EQ.ZERO) TEMP(I)=ZERO
  220 CONTINUE
      CALL VCLR(D,1,NC1)
      CALL VCLR(HESSL,1,NC1)
      DEPRE=ZERO
      DO 230 I=1,NC1
         IF (EIGEN(I).EQ. ZERO) THEN
            TMP=ZERO
         ELSE
            TMP=-TEMP(I)/EIGEN(I)
         END IF
         IF(OUT) WRITE(IW,444)I,TEMP(I),EIGEN(I),TMP
444      FORMAT(1X,'JJH, NR STEP',I5,3F15.10)
         TTMP = TMP*TEMP(I) + PT5*TMP*TMP*EIGEN(I)
         DEPRE = DEPRE+TTMP
         IF(I.EQ.IFOLOW)DEIFOL=TTMP
         DO 222 J=1,NC1
            D(J)=D(J)+TMP*UMAT(J,I)
            IF(I.EQ.IFOLOW)HESSL(J)=HESSL(J)+TMP*UMAT(J,I)
  222    CONTINUE
  230 CONTINUE
      DD = SQRT(DDOT(NC1,D,1,D,1))
      WRITE(IW,9090)DD
C
      IF (GOFRST .AND. DD.LE.SNRMAX .AND.
     *     ( (STPT .AND. NSERCH.GT.3) .OR. .NOT.STPT) ) THEN
          WRITE(IW,9018)DD
          CALL DCOPY(NC1,D,1,TEMP,1)
          GO TO 461
      ELSE
          DO 223 J=1,NC1
             D(J)=D(J)-HESSL(J)
  223    CONTINUE
         DEPRE=DEPRE-DEIFOL
         DD = SQRT(DDOT(NC1,D,1,D,1))
C        WRITE(IW,*) 'PROJECTED STEP LENGTH',DD
      END IF
      IF (DD.GE.DXMAXT) THEN
         ALPHA=ZERO
         SKAL=DXMAXT/DD
         WRITE(IW,*) ' SCALING PROJECTED NR STEP BY',SKAL
         DD=DXMAXT
      ELSE
C        ALPHA IS DETERMINED FROM PYTHAGORAS. IT IS POSITIVE AS
C        THE DIRECTION OF VMODE IS "CORRECT"
         ALPHA=SQRT(DXMAXT*DXMAXT-DD*DD)
         SKAL=ONE
      END IF
      WRITE(IW,9030)ALPHA,DD
      TTMP = ALPHA*TEMP(IFOLOW) + PT5*ALPHA*ALPHA*EIGEN(IFOLOW)
      DEPRE = DEPRE+TTMP
      DO 460 I=1,NC1
         TEMP(I)=SKAL*D(I)+ALPHA*VMODE(I)
  460 CONTINUE
  461 CONTINUE
      IF(OUT) WRITE(IW,9025) (TEMP(I),I=1,NC1)
      IF(OUT .AND. NZMAT.GT.0) THEN
         CALL DCOPY(NC1,TEMP,1,WORK,1)
         CALL TRANG(WORK,NZVAR,NC1)
         WRITE(IW,9026) (WORK(J),J=1,NZVAR)
      END IF
      CALL SYMDR(TEMP)
      CALL UPCOOR(TEMP,NC1,OUT)
      WRITE(IW,9020)NSERCH,ENERGY,GNORM,ANGLE
C     TOTX=SQRT(DDOT(NC1,TEMP,1,TEMP,1))
C     WRITE(IW,9040)NSERCH,DXMAXT,TOTX
C9040 FORMAT(5X,'CYCLE',I5,3X,'TRUST RADIUS',F12.8,3X,'STEP LENGTH',
C    *          F12.8)
      IF (NSERCH.GE.NSTEP) THEN
         WRITE(IW,9070)
         GO TO 999
      END IF
      GO TO 777
C
  999 CONTINUE
      CALL BIGFM(MAXMEM)
      CALL TMDATE(TIMSTR)
      IF (MASWRK) WRITE(IW,9002) TIMSTR
      CALL ENDING
      STOP
C
 9002 FORMAT(' EXECUTION OF GAMESS TERMINATED NORMALLY ',3A8)
 9001 FORMAT(/10X,'NSERCH=',I3,5X,'ENERGY=',F16.7)
 9011 FORMAT(/30X,18(1H-),25X,23(1H-),/30X,18HCOORDINATES (BOHR), 25X,
     *  23HGRADIENT (HARTREE/BOHR)/30X,18(1H-),25X,23(1H-)// ,8X,
     *  13HATOM     ZNUC,7X,1HX,13X,1HY,13X,1HZ,20X,1HX,13X,1HY, 13X,
     *  1HZ/1X,112(1H-)/)
 9012 FORMAT(I3,2X,A8,A2,F6.1,3F14.7,8X,3F14.7)
 9020 FORMAT(//,5X,'CYCLE',I5,3X,'ENERGY',F15.8,3X,'GNORM',F15.8,3X,
     $          'ANGLE',F15.8)
 9021 FORMAT(//,5X,'CYCLE',I5,3X,'ENERGY',F15.8,3X,'GNORM',F15.8,3X)
 9030 FORMAT(/10X,'STEPLENGTH ALONG            UPHILL MODE',F12.8,
     *       /10X,'STEPLENGTH PERPENDICULAR TO UPHILL MODE',F12.8)
 9100 FORMAT(/,20X,'CURRENT GEOMETRY',/)
 9110 FORMAT(//,5X,34(1H-),/,5X,
     $       '***     CALCULATING HESSIAN    ***',
     $       /,5X,34(1H-),//)
 9120 FORMAT(//,5X,34(1H-),/,5X,
     $       '*** END OF HESSIAN CALCULATION ***',
     $       /,5X,34(1H-),//)
 9070 FORMAT(5X,'EXCESS NUMBER OF STEPS')
 9090 FORMAT(1X,'PURE NR STEP HAS LENGTH',F12.6)
      END
C*MODULE GRADEX  *DECK GRXSET
      SUBROUTINE GRXSET(NCVAL,NCOORD,METHOD,OUT,
     *                  SNUMH,STPT,STSTEP,GOFRST,DELCOR,SNRMAX,
     *                  MYSTEP,HSDFDB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DOUBLE PRECISION METHOD
C
      LOGICAL RESTAR,OUT,HSSEND,GOPARR,DSKWRK,MASWRK
      LOGICAL GOTEG,PURIFY,STPT,GOFRST,HSDFDB
C
      PARAMETER (ZERO=0.0D+00)
      PARAMETER (NNAM=15)
      PARAMETER (MXATM=2000)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM),GRDXTR(1)
C
      COMMON /HSSCTL/ IHESSM,IHREP,HSSEND
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGF
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C
      DATA GRDXTR /8HGRADEX  /
      DATA QNAM /8HMETHOD  ,8HNSTEP   ,8HIFOLOW  ,8HDPRED   ,
     *           8HRESTAR  ,8HOPTTOL  ,8HSNUMH   ,8HHESS    ,
     *           8HSTPT    ,8HSTSTEP  ,8HSNRMAX  ,8HMYSTEP  ,
     *           8HGOFRST  ,8HDELCOR  ,8HHSDFDB  /
      DATA KQNAM/5,1,1,3,0,3,3,5,0,3,3,1,0,3,0/
C
      DATA GEJJH /8HJJH     /, GESR /8HSR      /
      DATA CARDS /8HREAD    /, CALC /8HCALC    /
C
C     ----- INITIALIZE THE GEOMETRY SEARCH -----
C     NCOORD IS THE FULL 3*NAT CARTESIAN SPACE, WHILE NCVAL IS
C     THE DIMENSIONALITY OF THE GEOMETRY SEARCH. CURRENTLY
C     THE ONLY OPTION IS NCVAL=NCOORD
C
      IF (MASWRK) WRITE(IW,9000)
      NCOORD = 3*NAT
      NCVAL = NCOORD
C     IF(NZMAT .GT. 0) NCVAL = NVAR
      GOTEG=.FALSE.
      PURIFY=.FALSE.
C
C     ----- READ NAMELIST -$GRADEX- -----
C
      METHOD = GESR
      NSTEP  = 50
      IHREP  = 1
      IFOLOW = 1
      DPRED  = 1.0D-01
      SNUMH  = 1.0D-04
      DELCOR = 1.0D-03
      SNRMAX = 0.0D+00
      STSTEP = 1.0D-02
      STPT   = .TRUE.
      GOFRST = .TRUE.
      RESTAR = .FALSE.
      HSDFDB = .TRUE.
      HESS   = CALC
      OPTTOL = 1.0D-04
      MYSTEP = 20
C
C
      JRET = 0
      CALL NAMEIO(IR,JRET,GRDXTR,NNAM,QNAM,KQNAM,
     *            METHOD,NSTEP,IFOLOW,DPRED,RESTAR,OPTTOL,SNUMH,
     *            HESS,STPT,STSTEP,SNRMAX,MYSTEP,GOFRST,DELCOR,HSDFDB,
     *            0,0,0,0,   0,0,0,0,0,
     *     0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *     0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
C
C    PRINT INFO MOVED ABOVE HESS CALL,
C    THIS MAY GIVE INCORRECT PRINT IF RESTART
C           ----- PRINT OPTIMIZATION INFO -----
C
      DXMAXT = DPRED
      IF (SNRMAX.EQ.ZERO) SNRMAX=MIN(0.10D+00,DXMAXT)
      IF (MASWRK) THEN
         WRITE(IW,9010) METHOD,NSTEP,IFOLOW,RESTAR,STPT,GOFRST
         WRITE(IW,9020) DPRED,DELCOR,SNUMH,OPTTOL,SNRMAX,MYSTEP,
     *                  STSTEP,HSDFDB
      END IF
 9010 FORMAT(/10X,'PARAMETERS CONTROLLING GEOMETRY SEARCH ARE'/
     *        10X,'METHOD =',A8,12X,'NSTEP  =',I10/
     *        10X,'IFOLOW =',I10,10X,'RESTAR =',L10/
     *        10X,'STPT   =',L10,10X,'GOFRST =',L10)
 9020 FORMAT(10X,'DPRED  =',1P,E10.3,10X,'DELCOR =',1P,E10.3/
     *       10X,'SNUMH  =',1P,E10.3,10X,'OPTTOL =',1P,E10.3/
     *       10X,'SNRMAX =',1P,E10.3,10X,'MYSTEP =',I10/
     *       10X,'STSTEP =',1P,E10.3,10X,'HSDFDB =',L10,/)
C
      NERR=0
      IF(JRET .GE. 2) NERR=NERR+1
      IF(METHOD.NE.GEJJH  .AND.  METHOD.NE.GESR) NERR=NERR+1
      IF (HESS.NE.CARDS .AND. HESS.NE.CALC) NERR=NERR+1
      IF(IFOLOW.EQ.0 .OR. ABS(IFOLOW).GT.NCOORD) NERR=NERR+1
      IF(DPRED.EQ.ZERO) NERR=NERR+1
C
C              RESTORE /OPTGRD/ FROM DAF (OLD $GRADEX VALUES)
C
      IF (.NOT.STPT) GO TO 100
      IF(RESTAR) THEN
         IF (MASWRK) WRITE(IW,9060)
         LENOPR = 3*MXATM + 34
         LENOPI = (31-1)/NWDVAR+1
         CALL DAREAD(IDAF,IODA,X ,LENOPR,37,0)
         CALL DAREAD(IDAF,IODA,IC,LENOPI,38,1)
         CALL DCOPY(NCOORD,X,1,C,1)
         IEXIT = 0
C
C              RESTORE NEW $GRADEX VALUES
C
      CALL NAMEIO(IR,JRET,GRDXTR,NNAM,QNAM,KQNAM,
     *            METHOD,NSTEP,IFOLOW,DPRED,RESTAR,OPTTOL,SNUMH,
     *            HESS,STPT,STSTEP,SNRMAX,MYSTEP,GOFRST,DELCOR,HSDFDB,
     *            0,0,0,0,  0,0,0,0,0,
     *    0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *    0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
         IF(JRET .GE. 2) NERR=NERR+1
         DXMAXT = DPRED
         IF (SNRMAX.EQ.ZERO) SNRMAX=MIN(0.10D+00,DXMAXT)
C
C           ----- OBTAIN INITIAL HESSIAN MATRIX -----
C
      ELSE
         IF (MASWRK) WRITE(IW,9005) HESS
         CALL HESSET(HESS,NCVAL,NCOORD,OUT,GOTEG,PURIFY)
C
C           ----- MOVE NUCLEAR COORDINATES INTO X -----
C
C        CALL DCOPY(NCOORD,C,1,X,1)
      END IF
  100 CONTINUE
      CALL DCOPY(NCOORD,C,1,X,1)
      CONVF = OPTTOL
C
      IF(NERR.EQ.0) RETURN
C
      IF (MASWRK) WRITE(IW,9055)
      CALL ABRT
C
 9000 FORMAT(//10X,21(1H-)/
     *       10X,'GRADIENT EXTREMAL RUN'/
     *       10X,21(1H-)/)
 9005 FORMAT(1X,'OBTAINING INITIAL HESSIAN, HESS=',A8)
 9055 FORMAT(1X,'***** ERROR IN $GRADEX INPUT *****')
 9060 FORMAT(1X,'***   RESTARTING GRADIENT EXTREMAL USING INFORMATION ',
     *       'FROM PREVIOUS RUN   ***')
      END
C*MODULE GRADEX  *DECK GRXSRM
      SUBROUTINE GRXSRM(NCOORD,STPT,GOFRST,OUT,
     *                  SNUMH,DELCOR,STSTEP,SNRMAX,MYSTEP,HSDFDB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL STPT,OUT,GOFRST,HSDFDB
C
      COMMON /FMCOM / X(1)
C
C     ----- GROW MEMORY FOR GRXSR -----
C     CURRENTLY IT IS QUITE 'PIGGY' WITH MEMORY,
C     MANY OF THE ARRAYS COULD BE OVERLAYED
C       -RM-      AT I10
C       -EIGEN-   AT I20
C       -TVEC-    AT I30
C       -EVEC-    AT I40
C       -ENEW-    AT I50
C       -ENEWU-   AT I60
C       -XVEC-    AT I70
C       -OLDX-    AT I80
C       -VOMODE-  AT I90
C       -DELX-    AT I100
C       -DELXU-   AT I110
C       -ALVEC-   AT I120
C       -APVEC-   AT I130
C       -AQVEC-   AT I140
C       -BLVEC-   AT I150
C       -BPVEC-   AT I160
C       -EVEC0-   AT I170
C       -EVEC1-   AT I180
C       -EINT-    AT I190
C       -AVEC-    AT I200
C       -BVEC-    AT I210
C       -JVEC-    AT I220
C       -KVEC-    AT I230
C       -WORK-    AT I240
C       -TEMP-    AT I250
C       -IA-      AT I260
C
C       -SCR-     AT I270
C
C       -HESSL-   AT I280
C       -HESS-    AT I290
C       -HSSP-    AT I300
C       -HSSM-    AT I310
C       -AMAT-    AT I320
C       -FMAT-    AT I330
C       -UMAT-    AT I340
C       -UPMAT-   AT I350
C       -UHUMAT-  AT I360
C       -UNRMAT-  AT I370
C       -TMAT-    AT I380
C       -TAMAT-   AT I390
C       -TPMAT-   AT I400
C       -XMAT-    AT I410
C       -YMAT-    AT I420
C       -ZMAT-    AT I430
C
      NC1 = NCOORD
      NC3 = NCOORD*NCOORD
      LOADFM = 0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I20 = I10 + NC1
      I30 = I20 + NC1
      I40 = I30 + NC1
      I50 = I40 + NC1
      I60 = I50 + NC1
      I70 = I60 + NC1
      I80 = I70 + NC1
      I90 = I80 + NC1
      I100 = I90 + NC1
      I110 = I100 + NC1
      I120 = I110 + NC1
      I130 = I120 + NC1
      I140 = I130 + NC1
      I150 = I140 + NC1
      I160 = I150 + NC1
      I170 = I160 + NC1
      I180 = I170 + NC1
      I190 = I180 + NC1
      I200 = I190 + NC1
      I210 = I200 + NC1
      I220 = I210 + NC1
      I230 = I220 + NC1
      I240 = I230 + NC1
      I250 = I240 + NC1
      I260 = I250 + NC1
      I270 = I260 + NC1
      I280 = I270 + 8*NC1
      I290 = I280 + NC3
      I300 = I290 + NC3
      I310 = I300 + NC3
      I320 = I310 + NC3
      I330 = I320 + NC3
      I340 = I330 + NC3
      I350 = I340 + NC3
      I360 = I350 + NC3
      I370 = I360 + NC3
      I380 = I370 + NC3
      I390 = I380 + NC3
      I400 = I390 + NC3
      I410 = I400 + NC3
      I420 = I410 + NC3
      I430 = I420 + NC3
      I440 = I430 + NC3
      LAST = I440
      NEED = LAST - I10
      CALL GETFM(NEED)
      CALL GRXSR(NC1,OUT,STPT,GOFRST,SNUMH,DELCOR,
     *           STSTEP,SNRMAX,MYSTEP,HSDFDB,
     *           X(I10),X(I20),X(I30),X(I40),X(I50),X(I60),
     *           X(I70),X(I80),X(I90),
     *           X(I100),X(I110),X(I120),
     *           X(I130),X(I140),X(I150),X(I160),X(I170),
     *           X(I180),X(I190),X(I200),X(I210),X(I220),
     *           X(I230),X(I240),X(I250),X(I260),X(I270),
     *           X(I280),X(I290),X(I300),X(I310),X(I320),
     *           X(I330),X(I340),X(I350),X(I360),X(I370),
     *           X(I380),X(I390),X(I400),X(I410),X(I420),
     *           X(I430))
C
C     ----- RETURN FAST MEMORY -----
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE GRADEX  *DECK GRXSR
      SUBROUTINE GRXSR(NC1,OUT,STPT,GOFRST,SNUMH,DELCOR,
     *                 STSTEP,SNRMAX,MYSTEP,HSDFDB,
     *                 RM,EIGEN,TVEC,EVEC,ENEW,ENEWU,
     *                 XVEC,OLDX,VOMODE,DELX,DELXU,
     *                 ALVEC,APVEC,AQVEC,BLVEC,BPVEC,
     *                 EVEC0,EVEC1,EINT,AVEC,BVEC,JVEC,KVEC,
     *                 WORK,TEMP,IA,SCR,
     *                 HESSL,HESS,HSSP,HSSM,
     *                 AMAT,FMAT,UMAT,UPMAT,UHUMAT,
     *                 UNRMAT,TMAT,TAMAT,TPMAT,
     *                 XMAT,YMAT,ZMAT)
C
C     THIS REPRESENT AN IMPLEMENTATION OF THE SUN AND RUEDENBERG
C     ALGORITHM FOR FOLLOWING GRADIENT EXTREMALS.
C     J. SUN, K. RUEDENBERG, J.CHEM.PHYS. 98, 9707 (1993)
C     WITH ADDITIONAL CAPABILITIES OF LOCATING BIFURCATION
C     AND TURNING POINTS
C     K. BONDENSGAARD, F. JENSEN, J.CHEM.PHYS. 104, 8025
C
C     IMPLEMENTED IN GAMESS SEPT. 1995
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL OUT,HSSEND,LINEAR
      LOGICAL STPT,GOFRST,LBIF,HSDFDB
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, PT5=0.5D+00, TWO=2.0D+00)
      PARAMETER (TM4=1.0D-04, TMSIX=1.0D-06, TMEIG=1.0D-08)
      PARAMETER (TMFIV=1.0D-05)
      PARAMETER (TMNIN=1.0D-09, TPSIX=1.0D+06, EPS=1.0D-14)
      PARAMETER (PT9=0.9D+00, TMT=1.0D-02)
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION CMASS(3)
      DIMENSION RM(NC1),EIGEN(NC1),TVEC(NC1),EVEC(NC1),ENEW(NC1),
     *          ENEWU(NC1),XVEC(NC1),OLDX(NC1),
     *          VOMODE(NC1),DELX(NC1),DELXU(NC1)
      DIMENSION ALVEC(NC1),APVEC(NC1),AQVEC(NC1),BLVEC(NC1),BPVEC(NC1),
     *          EVEC0(NC1),EVEC1(NC1),EINT(NC1),
     *          AVEC(NC1),BVEC(NC1),JVEC(NC1),KVEC(NC1)
      DIMENSION WORK(NC1),TEMP(NC1),IA(NC1),SCR(8,NC1)
      DIMENSION HESSL(NC1*NC1)
      DIMENSION HESS(NC1,NC1),HSSP(NC1,NC1),HSSM(NC1,NC1),
     *          AMAT(NC1,NC1),FMAT(NC1,NC1),
     *          UMAT(NC1,NC1),UPMAT(NC1,NC1),
     *          UHUMAT(NC1,NC1),UNRMAT(NC1,NC1),
     *          TMAT(NC1,NC1),TAMAT(NC1,NC1),TPMAT(NC1,NC1)
      DIMENSION XMAT(NC1,NC1),YMAT(NC1,NC1),ZMAT(NC1,NC1)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3*MXFRG),TORQ(3*MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FUNCT / ENERGY,EG(3*MXATM)
      COMMON /HSSCTL/ IHESSM,IHREP,HSSEND
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /OPTEF / D(3*MXATM+6*MXFRG),OLDF(3*MXATM+6*MXFRG),
     *                VMODE(3*MXATM+6*MXFRG),RADIUS
      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGF
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
C
      DATA SADPT  /8HSADPOINT/
      DATA GUGA   /8HGUGA    /
C
C     ----- GRADIENT EXTREMUM FOLLOWING USING THE ALGORITHM
C     OF SUN AND RUDENBERG (J.CHEM.PHYS. 98, 9707) ----
C     CURRENT VERSION RUNS IN CARTESIAN COORDINATES BUT ALLOWS PRINTING
C     OF THE GEOMETRY IN INTERNALS
C
C     CURRENT VERSION IMPLEMENTED IN GAMESS
C     SEPTEMBER 1995 BY F.JENSEN, ODENSE, DK
C
C     TVEC IS NORMALIZED GRADIENT
C     RM IS REDUCED MASSES, CURRENT = 1.0
C     EVEC, ENEW ARE CURRENT AND NEW GE TANGENTS IN CARTESIAN COORD.,
C         ENEWU IS ENEW IN LOCAL COORD
C     VOMODE IS PREVIOUS HESSIAN MODE BEING FOLLOWED
C     OLDX IS PREVIOUS GEOMETRY ON THE GE
C     DELX AND DELXU ARE CORRECTION STEP IN CARTESIAN AND LOCAL COORD.
C     ALVEC, APVEC, AQVEC ARE CURRENT, PREVIOUS AND PRE-PREVIOUS
C         A-MATRIX EIGENVALUES
C     BLVEC, BPVEC ARE CURRENT AND PREVIOUS B-VECTORS
C     JVEC AND KVEC KEEP TRACK OF ORDERING OF A-MATRIX EIGENVALUES
C     EVEC0 AND EVEC1 ARE GE TANGENTS AT BIFURCATION/TURNING POINTS
C     EINT IS THE INTERPOLATED GE TANGENT AT BIFURCATION/TURNING POINTS
C     AVEC AND BVEC ARE TEMP. STORAGE DURING BIRFUCATION/TURNING POINT
C         LOCATION
C     XVEC IS THE INTERPOLATION VECTOR DURING BIRFUCATION/TURNING POINT
C         LOCATION
C     HESS, HSSP, HSSM ARE HESSIANS AT CURRENT GEOMETRY AND THE TWO
C         DISPLACED GEOMETRIES ALONG THE GRADIENT
C     FMAT IS THE PARTIAL THIRD DERIVATIVE
C     AMAT IS THE 'A-MATRIX'
C     UMAT IS THE UNITARY TRANFORMATION TO THE LOCAL COORD. SYSTEM
C     UNRMAT IS THE UNITARY TRANSFORMATION WHICH DIAGONALIZE THE HESSIAN
C        UPMAT IS THE PREVIOUS UNRMAT
C     UHUMAT IS THE HESSIAN TRANSFORMED WITH UMAT
C     TMAT AND TAMAT ARE CURRENT A-MATRIX EIGENVECTORS IN LOCAL AND
C        CARTESIAN COORDINATES. TPMAT IS PREVIOUS TAMAT.
C     WORK,TEMP,SCR,IA,XMAT,YMAT,ZMAT ARE SCRATCH ARRAYS
C
C     NC1 IS = NCOORD
C
      TODEG=180.0D+00/ACOS(-ONE)
      LBIF=.FALSE.
      LCOUNT=0
      MICRO=0
      JAMODE = 0
      JNDEXA = 0
C
      DO 10 I=1,NC1
         ZMASS(I)=ONE
         RM(I)=ONE
   10 CONTINUE
      IF (STPT) THEN
         DO 20 I=1,NC1
            OLDX(I)=X(I)-STSTEP*VMODE(I)
            EVEC(I)=VMODE(I)
            VOMODE(I)=VMODE(I)
   20    CONTINUE
      ELSE
         CALL VCLR(OLDX,1,NC1)
         CALL VCLR(EVEC,1,NC1)
         CALL VCLR(VOMODE,1,NC1)
      END IF
C
C     THE SCF SHOULD BE CONVERGED TO AT LEAST ~10(-2) TIMES SNUMH
C     TO ENSURE STABILITY IN THE NUMERICAL HESSIAN DIFFERENTIATION
C
      IF(CONVHF.GT.0.01D+00*SNUMH) THEN
         CONVHF = 0.01D+00 * SNUMH
         WRITE(IW,9000) CONVHF
      END IF
      IF (.NOT.STPT) WRITE(IW,9010)
C
C     START OF MACRO LOOP
C
 100  CONTINUE
      IF(.NOT.LBIF)MICRO=-1
      IF(.NOT.LBIF)NSERCH=NSERCH+1
C
C     START OF MICRO LOOP
C
 110  CONTINUE
      IF(.NOT.LBIF)MICRO=MICRO+1
      IF (MICRO.EQ.0) THEN
         WRITE(IW,9020)
         CALL PRATM(IW,2)
         IF(NZMAT.GT.0) CALL PZANDG(EG,0)
      END IF
C     CALCULATE GRADIENT
      IF (HSDFDB) THEN
         WRITE(IW,9065)
         CALL ENERGX
         IF (ENERGY.EQ.ZERO) THEN
            WRITE(IW,9030)
            CALL ABRT
         END IF
         CALL VCLR(EG,1,NC1)
         IF (CITYP.EQ.GUGA) THEN
            CALL CIGRAD
         ELSE
            CALL HFGRAD
         END IF
         WRITE(IW,9066)
      ELSE
         WRITE(IW,9067)
         CALL HESSX(.FALSE.,.FALSE.)
         CALL DAREAD(IDAF,IODA,HESS,NC1*NC1,4,0)
         CALL DCOPY(NC1*NC1,HESS,1,HSSM,1)
         WRITE(IW,9068)
      END IF
      WRITE (IW,9040) NSERCH,ENERGY
      WRITE (IW,9050)
      DO 120 IAT = 1,NAT
         N = 3*(IAT-1)
         ZNUC = ZAN(IAT) + IZCORE(IAT)
         WRITE (IW,9060) IAT,ANAM(IAT),BNAM(IAT),ZNUC,
     *                   (C(J,IAT),J=1,3), (EG(N+I),I=1,3)
  120 CONTINUE
      IF(NZMAT.GT.0) THEN
         CALL DCOPY(NC1,EG,1,TEMP,1)
         CALL TRANG(TEMP,NZVAR,NC1)
         CALL PZANDG(TEMP,1)
      END IF
C     TVEC IS NORMALIZED GRADIENT
      GNORM=SQRT(DDOT(NC1,EG,1,EG,1))
      CALL DCOPY(NC1,EG,1,TVEC,1)
      TMP=ONE/GNORM
      CALL DSCAL(NC1,TMP,TVEC,1)
      IF(OUT) WRITE(IW,*) 'TVEC'
      IF(OUT) WRITE(IW,*) (TVEC(I),I=1,NC1)
C
C     CONSTRUCT X' POINTS FOR NUMERICAL DIFFERENTIATION OF H
C
      WRITE(IW,9070)
      CALL DAXPY(NC1,SNUMH,TVEC,1,X,1)
      CALL DCOPY(NC1,X,1,C,1)
      CALL HESSX(.FALSE.,.FALSE.)
      CALL DAREAD(IDAF,IODA,HSSP,NC1*NC1,4,0)
C
      IF (HSDFDB) THEN
         CALL DAXPY(NC1,-TWO*SNUMH,TVEC,1,X,1)
         CALL DCOPY(NC1,X,1,C,1)
         CALL HESSX(.FALSE.,.FALSE.)
         CALL DAREAD(IDAF,IODA,HSSM,NC1*NC1,4,0)
         CALL DAXPY(NC1,SNUMH,TVEC,1,X,1)
         CALL DCOPY(NC1,X,1,C,1)
         CALL VADD(HSSP,1,HSSM,1,HESS,1,NC1*NC1)
         CALL DSCAL(NC1*NC1,PT5,HESS,1)
      ELSE
         CALL DAXPY(NC1,-SNUMH,TVEC,1,X,1)
         CALL DCOPY(NC1,X,1,C,1)
      END IF
      IF(OUT) WRITE(IW,*) 'HESS'
      IF(OUT) CALL PRSQ(HESS,NC1,NC1,NC1)
      WRITE(IW,9080)
C
C     PROJECT HESSIAN. USE FMAT FOR STORING PROJECTED HESS
      CALL DCOPY(NC1*NC1,HESS,1,FMAT,1)
      CALL CENMAS(NAT,NFRG,C,TEMP,ZMASST,CMASS,ZMASS)
      CALL PRJFC(.FALSE.,.TRUE.,ZMASST,FMAT,WORK,TEMP,RM,FRGMI,
     *              XMAT,YMAT,NAT,NC1,NC1)
C     DIAGONALIZE HESS TO DETERMINE NUMBER OF T, R MODE
C     AND CALCULATE PURE NR STEPLENGTH
C     POSSIBLY SWITCH TO STANDARD OPTIMIZATION
C
      CALL CPYSQT(FMAT,HESSL,NC1,1)
      CALL GLDIAG(NC1,NC1,NC1,HESSL,SCR,EIGEN,UNRMAT,IERR,IA)
      NTR=0
      DO 125 I=NC1,1,-1
         IF (ABS(EIGEN(I)).LT.TMSIX) EIGEN(I)=ZERO
         IF (EIGEN(I).EQ.ZERO) NTR=NTR+1
         IF (EIGEN(I).EQ.ZERO) ISHIFT=I
  125 CONTINUE
C
C     SHIFT T,R MODES TO BE THE LAST NTR
C
      JSHIFT=ISHIFT-1
      DO 140 J=1,NTR
         JSHIFT=JSHIFT+1
         TEMP(J)=EIGEN(JSHIFT)
         DO 130 I=1,NC1
            XMAT(J,I)=UNRMAT(I,JSHIFT)
  130    CONTINUE
  140 CONTINUE
      DO 160 J=ISHIFT,NC1-NTR
         EIGEN(J)=EIGEN(J+NTR)
         DO 150 I=1,NC1
            UNRMAT(I,J)=UNRMAT(I,J+NTR)
  150    CONTINUE
  160 CONTINUE
      JTMP=0
      DO 180 J=NC1-NTR+1,NC1
         JTMP=JTMP+1
         EIGEN(J)=TEMP(JTMP)
         DO 170 I=1,NC1
            UNRMAT(I,J)=XMAT(JTMP,I)
  170    CONTINUE
  180 CONTINUE
C
      IF(OUT) WRITE(IW,*) 'HESSIAN EIGENVALUES',(EIGEN(I),I=1,NC1)
      RADIUS=ZERO
      DO 190 I=1,NC1-NTR
         TMP=DDOT(NC1,EG,1,UNRMAT(1,I),1)
         TMP=TMP/EIGEN(I)
         RADIUS=RADIUS+TMP*TMP
  190 CONTINUE
      RADIUS=SQRT(RADIUS)
      WRITE(IW,9090)RADIUS
      IF (GOFRST .AND. RADIUS.LE.SNRMAX .AND.
     *     ( (STPT .AND. NSERCH.GT.3) .OR. .NOT.STPT) ) THEN
          WRITE(IW,9100)
          IHREP = 1
          HSSEND = .FALSE.
          CALL DCOPY(NC1,VOMODE,1,VMODE,1)
          RUNTYP = SADPT
          FMAXT=1.0D+01
          NSTEP=10
          NZVAR=0
          NZMAT=0
C         CURRENTLY WE HAVE EITHER HSSP OR HSSM ON DISK
          CALL DAWRIT(IDAF,IODA,HESS,NC1*NC1,4,0)
          CALL SIGX(.TRUE.)
          RETURN
      END IF
C
      IF (NTR.NE.5 .AND. NTR.NE.6) THEN
         WRITE(IW,9110)
         WRITE(IW,*) (EIGEN(I),I=1,NC1)
         CALL ABRT
      END IF
C
C     DETERMINE THE UNITARY TRANSFORMATION WHICH TAKES THE CARTESIAN
C     COORDINATE SYSTEM TO ONE WHERE TVEC IS ALONG X1, THE T, R MODES
C     ALONG THE LAST 6(5), AND MAKES THE REST OF THE HESSIAN DIAGONAL
C
C     PROJECT OUT TVEC FROM HESSIAN
C
      DO 210 I=1,NC1
         DO 200 J=1,NC1
            ZMAT(I,J)=-TVEC(I)*TVEC(J)
  200    CONTINUE
  210 CONTINUE
      DO 220 I=1,NC1
         ZMAT(I,I)=ONE+ZMAT(I,I)
  220 CONTINUE
      CALL TFSQU(XMAT,FMAT,ZMAT,WORK,NC1,NC1)
      CALL CPYSQT(XMAT,HESSL,NC1,1)
      CALL GLDIAG(NC1,NC1,NC1,HESSL,SCR,WORK,UMAT,IERR,IA)
      NZ=0
      DO 230 I=NC1,1,-1
         IF(ABS(WORK(I)).LT.TMSIX)WORK(I)=ZERO
         IF(WORK(I).EQ.ZERO)NZ=NZ+1
         IF(WORK(I).EQ.ZERO)ISHIFT=I
  230 CONTINUE
      IF (NZ.NE.NTR+1) THEN
         WRITE(IW,*) ' ERROR! NOT NTR+1 EIGENVALUES OF TPHPT=0'
         WRITE(IW,*) (WORK(I),I=1,NC1)
         CALL ABRT
      END IF
C
C     FIRST MOVE ALL ZERO EIGENVALUES TO THE LAST NTR+1
C
      JSHIFT=ISHIFT-1
      DO 250 J=1,NZ
         JSHIFT=JSHIFT+1
         DO 240 I=1,NC1
            XMAT(J,I)=UMAT(I,JSHIFT)
  240    CONTINUE
  250 CONTINUE
      DO 270 J=ISHIFT,NC1-NZ
         DO 260 I=1,NC1
            UMAT(I,J)=UMAT(I,J+NZ)
  260    CONTINUE
  270 CONTINUE
      JTMP=0
      DO 290 J=NC1-NZ+1,NC1
         JTMP=JTMP+1
         DO 280 I=1,NC1
            UMAT(I,J)=XMAT(JTMP,I)
  280    CONTINUE
  290 CONTINUE
C
C     THEN MOVE ALL ONE UP, INSET TVEC AS THE FIRST COLUMN
C     AND GET LAST NTR TR MODES FROM UNRMAT
C
      DO 310 J=NC1-NZ,1,-1
         DO 300 I=1,NC1
            UMAT(I,J+1)=UMAT(I,J)
  300    CONTINUE
  310 CONTINUE
      DO 320 I=1,NC1
         UMAT(I,1)=TVEC(I)
  320 CONTINUE
      DO 340 J=NC1-NTR+1,NC1
         DO 330 I=1,NC1
            UMAT(I,J)=UNRMAT(I,J)
  330    CONTINUE
  340 CONTINUE
      IF(OUT) WRITE(IW,*) 'U-MATRIX'
      IF(OUT) CALL PRSQ(UMAT,NC1,NC1,NC1)
      CALL TFSQU(UHUMAT,FMAT,UMAT,WORK,NC1,NC1)
      IF(OUT) WRITE(IW,*) 'UHU-MATRIX'
      IF(OUT) CALL PRSQ(UHUMAT,NC1-NTR,NC1-NTR,NC1)
C
C     F = DH/DX1
C
      STEP=ONE/SNUMH
      IF (HSDFDB) STEP=STEP/TWO
      CALL VSUB(HSSM,1,HSSP,1,FMAT,1,NC1*NC1)
      CALL DSCAL(NC1*NC1,STEP,FMAT,1)
      IF(OUT) THEN
         CALL TFSQU(XMAT,FMAT,UMAT,WORK,NC1,NC1)
         WRITE(IW,*) 'F-MATRIX'
         CALL PRSQ(XMAT,NC1-NTR,NC1-NTR,NC1)
      END IF
C
C     CONSTRUCT A' MATRIX = DH/DX1*|G| + H**2 + H11*H
C
      DO 390 I=1,NC1
         DO 380 J=1,NC1
            TMP=ZERO
            DO 370 K=1,NC1
               TMP=TMP+HESS(I,K)*HESS(K,J)
  370       CONTINUE
            XMAT(I,J)=TMP
  380    CONTINUE
  390 CONTINUE
      TMP=UHUMAT(1,1)
      DO 410 I=1,NC1
         DO 400 J=1,NC1
            YMAT(I,J)=GNORM*FMAT(I,J)+XMAT(I,J)-TMP*HESS(I,J)
  400    CONTINUE
  410 CONTINUE
C
C     TRANSFORM A' TO LOCAL, AND REMOVE FIRST COMPONENT
C     A' IS NOW REGTANGULAR: NC1-NTR-1, NC1-NTR
C
      CALL TFSQU(XMAT,YMAT,UMAT,WORK,NC1,NC1)
      DO 412 I=2,NC1-NTR
         DO 411 J=1,NC1-NTR
            AMAT(I-1,J)=XMAT(I,J)
  411    CONTINUE
  412 CONTINUE
      IF(OUT) WRITE(IW,*) 'A-PRIME MATRIX'
      IF(OUT) CALL PRSQ(AMAT,NC1-NTR,NC1-NTR-1,NC1)
C
C     CALCULATE CORRECTION STEP
C     IN LOCAL COORDINATES THE FVEC IS SIMPLY=(0,H21,H31,...)
C
      DO 420 I=2,NC1-NTR
         TEMP(I-1)=GNORM*UHUMAT(I,1)
  420 CONTINUE
      DO 430 I=NC1-NTR,NC1
         TEMP(I)=ZERO
  430 CONTINUE
      IF(OUT) TMP=SQRT(DDOT(NC1,TEMP,1,TEMP,1))
C     WRITE(IW,*) 'THT NORM',TMP/GNORM
      IF(OUT) WRITE(IW,*) 'FVEC NORM',TMP
      IF(OUT) WRITE(IW,*) 'FVEC',(TEMP(I),I=1,NC1-NTR-1)
C
C     FORM (A'*A'T)**-1
C
      DO 460 I=1,NC1-NTR-1
         DO 450 J=1,NC1-NTR-1
            TMP=ZERO
            DO 440 K=1,NC1-NTR
               TMP=TMP+AMAT(I,K)*AMAT(J,K)
  440       CONTINUE
            XMAT(I,J)=TMP
  450    CONTINUE
  460 CONTINUE
      IF(OUT) WRITE(IW,*) 'AA-DAGGER MATRIX'
      IF(OUT) CALL PRSQ(XMAT,NC1-NTR-1,NC1-NTR-1,NC1)
C
C     MOVE TO SYSTEM WHERE AA-DAGGER IS DIAGONAL IN ORDER
C     TO IDENTIFY SINGULARITIES WHERE CORRECTOR STEP IS UNSTABLE
C
      IJ=0
      DO 480 I=1,NC1-NTR-1
         DO 470 J=1,I
            IJ=IJ+1
            HESSL(IJ)=XMAT(I,J)
  470    CONTINUE
  480 CONTINUE
      CALL GLDIAG(NC1,NC1-NTR-1,NC1-NTR-1,
     *            HESSL,SCR,WORK,ZMAT,IERR,IA)
      IF(OUT) WRITE(IW,*) 'AA-DAGGER EIGENVALUES AND -VECTORS'
      IF(OUT) WRITE(IW,*) (WORK(I),I=1,NC1-NTR-1)
      IF(OUT) CALL PRSQ(ZMAT,NC1-NTR-1,NC1-NTR-1,NC1)
C     USE XMAT FOR TEMPORARY MATRIX (A-DAGGER*ZMAT)
C     USE SCR(1,I) FOR TEMPORARY VECTOR (AA-DAGGER)-1 * FVEC
      CALL VCLR(SCR,1,8*NC1)
      DO 510 I=1,NC1-NTR
         DO 500 J=1,NC1-NTR-1
            SUM=ZERO
            DO 490 K=1,NC1-NTR-1
               SUM=SUM+AMAT(K,I)*ZMAT(K,J)
  490       CONTINUE
         XMAT(I,J)=SUM
  500    CONTINUE
  510 CONTINUE
C     REMOVE COMPONENTS CORRSPONDING TO NEAR-SINGULARITIES
      IF (ABS(WORK(1)).LT.TMEIG) WRITE(IW,9120)WORK(1)
      DO 540 I=1,NC1-NTR-1
         TMP=DDOT(NC1-NTR-1,TEMP,1,ZMAT(1,I),1)
         TMPE=WORK(I)
         IF (TMPE.NE.ZERO) THEN
            FSTEP=TMP/TMPE
         ELSE
            FSTEP=999.999D+00
         END IF
         SUM=ZERO
         DO 520 K=1,NC1-NTR
            SUM=SUM+(XMAT(K,I)*FSTEP)**2
  520    CONTINUE
         SUM=SQRT(SUM)
         IF(OUT) WRITE(IW,*) 'CORRECTOR STEP COMPONENT',I,SUM
         IF (ABS(TMPE).LT.TMEIG) THEN
            WRITE(IW,9130)I,TMPE,TMP,SUM
            FSTEP=ZERO
         END IF
         DO 530 J=1,NC1-NTR-1
            SCR(1,J)=SCR(1,J)+FSTEP*ZMAT(J,I)
  530    CONTINUE
  540 CONTINUE
      DO 560 I=1,NC1-NTR
         TMP=ZERO
         DO 550 J=1,NC1-NTR-1
            TMP=TMP+AMAT(J,I)*SCR(1,J)
  550    CONTINUE
         DELXU(I)=TMP
  560 CONTINUE
      IF(OUT) WRITE(IW,*) 'DELX IN LOCAL',(DELXU(I),I=1,NC1-NTR)
C
C     CALCULATE NEW GE TANGENT
C     FIRST COLOUMN OF A' MATRIX IS THE B VECTOR
C     REST IS THE A MATRIX TO BE INVERTED
C     NOTE THAT THERE IS A MISSING MINUS SIGN ON THE B VECTOR
C     IN THE PAPER
C
      DO 580 I=1,NC1-NTR-1
         TEMP(I)=-AMAT(I,1)
         DO 570 J=2,NC1-NTR
            XMAT(I,J-1)=AMAT(I,J)
  570    CONTINUE
  580 CONTINUE
      IF(OUT) WRITE(IW,*) 'A MATRIX'
      IF(OUT) CALL PRSQ(XMAT,NC1-NTR-1,NC1-NTR-1,NC1)
C
C     MOVE TO COORDINATE SYSTEM WHERE A MATRIX IS DIAGONAL
C
      IJ=0
      DO 600 I=1,NC1-NTR-1
         DO 590 J=1,I
            IJ=IJ+1
            HESSL(IJ)=XMAT(I,J)
  590    CONTINUE
  600 CONTINUE
      CALL GLDIAG(NC1,NC1-NTR-1,NC1-NTR-1,
     *            HESSL,SCR,ALVEC,TMAT,IERR,IA)
      WRITE(IW,9140)
      INDEXA=0
      DO 610 I=1,NC1-NTR-1
         BLVEC(I)=DDOT(NC1-NTR-1,TEMP,1,TMAT(1,I),1)
         WRITE(IW,9150)I,ALVEC(I),ABS(BLVEC(I))
         IF(ALVEC(I).LT.ZERO)INDEXA=INDEXA+1
  610 CONTINUE
      WRITE(IW,9160)INDEXA
      IF(OUT) WRITE(IW,*) 'A-MATRIX EIGENVECTORS'
      IF(OUT) WRITE(IW,*) (ALVEC(K),K=1,NC1-NTR-1)
      IF(OUT) CALL PRSQ(TMAT,NC1-NTR-1,NC1-NTR-1,NC1)
C
C     BACKTRANSFORM A-MAT EIGENVECTORS TO CARTESIAN
C     AND STORE FOR ESTABLISHING CONNECTION
C
      DO 640 K=1,NC1-NTR-1
         DO 630 I=1,NC1
            TMP=ZERO
            DO 620 J=2,NC1-NTR
               TMP=TMP+UMAT(I,J)*TMAT(J-1,K)
  620       CONTINUE
            TAMAT(I,K)=TMP
  630    CONTINUE
  640 CONTINUE
      CALL VCLR(ENEWU,1,NC1)
      IF(OUT) WRITE(IW,*) 'EVEC IN INTERNAL INTERNAL'
      DO 660 I=1,NC1-NTR-1
         IF (ABS(ALVEC(I)).GT.TMNIN) THEN
            TMP=BLVEC(I)/ALVEC(I)
         ELSE
            TMP=ZERO
         END IF
         DO 650 J=1,NC1-NTR-1
            ENEWU(J+1)=ENEWU(J+1)+TMP*TMAT(J,I)
  650    CONTINUE
  660 CONTINUE
      ENEWU(1)=ONE
      IF(OUT) WRITE(IW,*) 'ENEWU BEFORE NORMALIZATION'
      IF(OUT) WRITE(IW,9175) (ENEWU(I),I=1,NC1-NTR)
      ENTMP=ONE/SQRT(DDOT(NC1-NTR,ENEWU,1,ENEWU,1))
      CALL DSCAL(NC1-NTR,ENTMP,ENEWU,1)
      IF(OUT) WRITE(IW,*) 'GE TANGENT IN LOCAL COORDINATES'
      IF(OUT) WRITE(IW,9175) (ENEWU(I),I=1,NC1-NTR)
C
C     BACKTRANSFORM ENEWU AND DELXU VECTORS TO CARTESIAN COORDINATES
C
      DO 680 I=1,NC1
         TMP1=ZERO
         TMP2=ZERO
         DO 670 J=1,NC1-NTR
            TMP1=TMP1+UMAT(I,J)*ENEWU(J)
            TMP2=TMP2+UMAT(I,J)*DELXU(J)
  670    CONTINUE
         ENEW(I)=TMP1
         DELX(I)=-TMP2
  680 CONTINUE
      DNORM=SQRT(DDOT(NC1,DELX,1,DELX,1))
      IF(OUT) WRITE(IW,*) 'DNORM',DNORM
      IF(OUT) WRITE(IW,*) 'DELX',(DELX(I),I=1,NC1)
      IF(OUT) WRITE(IW,*) 'NEW EVEC',(ENEW(I),I=1,NC1)
      ANG=DDOT(NC1,ENEW,1,EVEC,1)
      IF (ABS(ANG).GT.ONE)ANG=ANG/ABS(ANG)
      IF (ANG.LT.ZERO) THEN
C     THE GE TANGENT SHOULD ALWAYS POINT IN ROUGHLY THE SAME DIRECTION
         IF(OUT) WRITE(IW,*) 'REVERSING E-VECTOR DIRECTION',ANG
         CALL DSCAL(NC1,-ONE,ENEW,1)
         ANG=-ANG
      END IF
      CALL VSUB(EVEC,1,ENEW,1,TEMP,1,NC1)
      GECURV=SQRT(DDOT(NC1,TEMP,1,TEMP,1))
      TANG=TODEG*ACOS(ANG)
      CALL VSUB(X,1,OLDX,1,TEMP,1,NC1)
      TOTX=SQRT(DDOT(NC1,TEMP,1,TEMP,1))
      GECURV=GECURV/TOTX
      IF(NZMAT.GT.0 .AND. OUT) THEN
         CALL DCOPY(NC1,ENEW,1,WORK,1)
         CALL TRANG(WORK,NZVAR,NC1)
         WRITE(IW,9170)
         WRITE(IW,9175) (WORK(I),I=1,NZVAR)
      END IF
C
      IF (.NOT.STPT .AND. NSERCH.EQ.1) THEN
         IF (DNORM.LT.DELCOR) THEN
            WRITE(IW,9180)IFOLOW
            TMP=DDOT(NC1,TVEC,1,ENEW,1)
C
C       IF TMP IS POSITIVE,
C       THE GE TANGENT POINTS IN THE UPHILL DIRECTION
C
            IF (TMP.LT.ZERO) CALL DSCAL(NC1,-ONE,ENEW,1)
            IF (IFOLOW.LT.0) CALL DSCAL(NC1,-ONE,ENEW,1)
            TMP=DDOT(NC1,TVEC,1,ENEW,1)
            IFOLOW=ABS(IFOLOW)
            TANG=99.99D+00
            TOTX=9.999999D+00
            GECURV=99.999D+00
            DO 690 I=1,NC1
               EVEC(I)=ENEW(I)
               VOMODE(I)=UNRMAT(I,IFOLOW)
  690       CONTINUE
         ELSE
            TANG=99.99D+00
         END IF
      END IF
C
C     CALCULATE CURVATURE OF EACH HESSIAN EIGENVECTOR
C
      CALL VCLR(WORK,1,NC1)
      IF (DNORM.LT.DELCOR .AND. NSERCH.GT.1) THEN
C        ESTABLISH CONNECTION AND CALCULATE CURVATURE
         DO 720 I=1,NC1-NTR
            TEST=ZERO
            DO 700 J=1,NC1
               TMP=DDOT(NC1,UNRMAT(1,I),1,UPMAT(1,J),1)
               IF (ABS(TMP).GT.ABS(TEST)) IJ=J
               IF (ABS(TMP).GT.ABS(TEST)) TEST=TMP
  700       CONTINUE
            TDIR=ONE
            IF(TEST.LT.ZERO)TDIR=-ONE
            TMP=ZERO
            DO 710 J=1,NC1
               TMP=TMP+(UNRMAT(J,I)-TDIR*UPMAT(J,IJ))**2
  710       CONTINUE
            WORK(I)=SQRT(TMP)/TOTX
  720    CONTINUE
      END IF
C
C     WRITE OUT PROJECTION OF GRADIENT, NEW EVEC, PREVIOUS TS-MODE
C     ALONG HESSIAN EIGENVECTORS, AND EIGENVECTOR CURVATURES
C
      TMPM=ZERO
      TMPP=ZERO
      IMODE=0
      IPRO =0
      WRITE(IW,9200)
      DO 730 I=1,NC1-NTR
         TMP1=DDOT(NC1,TVEC,1,UNRMAT(1,I),1)
         TMP2=DDOT(NC1,ENEW,1,UNRMAT(1,I),1)
         TMP3=DDOT(NC1,VOMODE,1,UNRMAT(1,I),1)
         WRITE(IW,9210)I,EIGEN(I),TMP1,TMP2,TMP3,WORK(I)
         IF (ABS(TMP1).GT.TMPP) IPRO=I
         IF (ABS(TMP1).GT.TMPP) TMPP=ABS(TMP1)
         IF (ABS(TMP3).GT.TMPM) IMODE=I
         IF (ABS(TMP3).GT.TMPM) TMPM=ABS(TMP3)
  730 CONTINUE
      IF(IMODE.NE.IFOLOW .AND. IMODE.NE.0) WRITE(IW,9220)IFOLOW,IMODE
C
C     ADD CORRECTION STEP
C
      IF (.NOT.LBIF) THEN
         CALL SYMDR(DELX)
         CALL UPCOOR(DELX,NC1,OUT)
         WRITE(IW,9230)
         CALL PRATM(IW,2)
         IF(NZMAT.GT.0) CALL PZANDG(EG,0)
      END IF
      WRITE(IW,9240)NSERCH,MICRO,DNORM,TANG,ENERGY
CTMP
      IF (DNORM.GT.DELCOR) THEN
         IF (MICRO.GE.MYSTEP) THEN
            WRITE(IW,9250)
            GO TO 900
         END IF
         WRITE(IW,9260)DNORM,DELCOR
         GO TO 110
      END IF
      IF(.NOT.LBIF) WRITE(IW,9261)DNORM,DELCOR
      IF (IPRO.NE.IMODE) THEN
         WRITE(IW,9270)
C        CALL ABRT
      END IF
      WRITE(IW,9280)NSERCH,ENERGY,GNORM,GECURV,
     *              NSERCH,TANG,
     *              NSERCH,DXMAXT,DNORM,TOTX,
     *              TODEG*ACOS(ENEWU(1))
      IFOLOW=IMODE
C
C     LOOK FOR GE-BIFURCATIONS
C     IF THE NUMBER OF NEGATIVE EIGENVALUES IN THE A MATRIX CHANGE,
C     FIND THE POINT WHERE ONE OF THE EIGENVALUES ARE ZERO AND
C     CHECK THE MAGNITUDE OF THE CORRESPONDING B-ELEMENT. IF THIS
C     ALSO IS "ZERO", WE HAVE A GE-BIFURCATION
C
      IAMODE=0
      NAMODE=0
      ALPH0=ZERO
      BETA0=ZERO
      IF (NSERCH.GT.1) THEN
C
C     ESTABLISH CONNECTION AND CHECK IF ALPHA VALUES ARE CHANGING SIGN
C     OR MAY HAVE TOUCHED ZERO
C
         DO 750 I=1,NC1-NTR-1
            TEST=ZERO
            DO 740 J=1,NC1-NTR-1
               TMP=DDOT(NC1,TAMAT(1,I),1,TPMAT(1,J),1)
               TMP=ABS(TMP)
               IF (TMP.GT.TEST)JTMP=J
               IF (TMP.GT.TEST)TEST=TMP
  740       CONTINUE
            WRITE(IW,9290)I,JTMP,TEST
            IF(.NOT.LBIF)KVEC(I)=JVEC(I)
            IF(.NOT.LBIF)JVEC(I)=JTMP
            IF (LBIF .AND. JTMP.EQ.JAMODE) THEN
               ALPH0=ALVEC(I)
               BETA0=ABS(BLVEC(I))
               IAMODE=I
            END IF
            IF (.NOT.LBIF .AND. ALVEC(I)*APVEC(JTMP).LT.ZERO) THEN
               IF(NAMODE.EQ.0)JAMODE=JTMP
               IF(NAMODE.EQ.0)IAMODE=I
               NAMODE=NAMODE+1
               CALL DCOPY(NC1,ALVEC,1,AVEC,1)
               CALL DCOPY(NC1,BLVEC,1,BVEC,1)
            END IF
  750    CONTINUE
      END IF
C     HAS AN ALPHA MODE MAYBE TOUCHED ZERO?
      IF (NSERCH.GT.2 .AND. .NOT.LBIF) THEN
         DO 760 I=1,NC1-NTR-1
            TMP0=ABS(AQVEC(KVEC(JVEC(I))))
            TMP1=ABS(APVEC(JVEC(I)))
            TMP2=ABS(ALVEC(I))
            IF ( (TMP1.LT.TMP0 .AND. TMP1.LT.TMP2
     *                         .AND. TMP1.LT.TMT) )
     *         WRITE(IW,9300)I,AQVEC(KVEC(JVEC(I))),APVEC(JVEC(I)),
     *                         ALVEC(I),ABS(BPVEC(JVEC(I)))
  760    CONTINUE
      END IF
C
C     ALPHA INDEX HAS CHANGED
C
      IF (NAMODE.NE.0) THEN
         WRITE(IW,9310)NAMODE,JNDEXA,INDEXA
         IF (NAMODE.GT.1) THEN
            WRITE(IW,9320)
C           CALL ABORT
         END IF
         WRITE(IW,9330)IAMODE
         LBIF=.TRUE.
         CALL VSUB(OLDX,1,X,1,XVEC,1,NC1)
         XNORM=SQRT(DDOT(NC1,XVEC,1,XVEC,1))
         CALL DCOPY(NC1,EVEC,1,EVEC0,1)
         CALL DCOPY(NC1,ENEW,1,EVEC1,1)
         IASAVE=INDEXA
         SINT0=ZERO
         SINT1=ONE
         ALP0=APVEC(JAMODE)
         ALP1=ALVEC(IAMODE)
         BET0=ABS(BPVEC(JAMODE))
         BET1=ABS(BLVEC(IAMODE))
         WRITE(IW,9340)
         WRITE(IW,9350)SINT0,ALP0,BET0
         WRITE(IW,9350)SINT1,ALP1,BET1
         SINTER=ALP0/(ALP0-ALP1)
         CALL DCOPY(NC1,OLDX,1,X,1)
         CALL DCOPY(NC1,XVEC,1,TEMP,1)
         CALL DSCAL(NC1,SINTER,TEMP,1)
         CALL UPCOOR(TEMP,NC1,OUT)
         DELCCR=DELCOR
         DELCOR=TPSIX
         GO TO 100
      END IF
C
C     USE FALSE POSITION TO HUNT DOWN THE ROOT.
C
      IF (LBIF) THEN
         LCOUNT=LCOUNT+1
         CALL VSUB(OLDX,1,X,1,TEMP,1,NC1)
         TMP=SQRT(DDOT(NC1,TEMP,1,TEMP,1))
         SINTER=TMP/XNORM
         WRITE(IW,9350)SINTER,ALPH0,BETA0
C
C        ROOT FINDING TERMINATING
C        WITH THE DEFAULT SETTINGS, THE ROOT FINDING QUITE OFTEN
C        HAVE TROUBLE LOCATING THE ROOT. THE PROBLEM IS THAT THE
C        INTERPOLATION IS DONE FROM GEOMETRIES WHICH HAVE THE
C        FINAL CORRECTION STEP ADDED, BUT THE GE INDEX IS FROM
C        THE GEOMETRY WITHOUT THE CORRECTION. THE SINTER = 0 AND 1
C        END POINTS ARE THUS NOT THE POINTS AT WHICH THE GE INDEX
C        CHANGE IS INFERRED. A REDUCED VALUE OF DELCOR SHOULD
C        FIX THIS PROBLEM. ALTERNATIVELY, THE FINAL CORRECTION STEP
C        COULD BE REMOVED. THE FIRST ALTERNATIVE IS EXPENSIVE, THE
C        LATTER IN GENERAL GIVES A LESS ACCURATE FOLLOWING. WE THUS
C        CHOOSE A LESS EXPENSIVE AND MORE ACCURATE GE FOLLOWING,
C        AT THE EXPENSE OF A LESS STABLE BF/TP LOCATION.
C
         IF (ABS(ALPH0).LT.TMFIV .OR. LCOUNT.GE.10) THEN
            IF(LCOUNT.LT.10) WRITE(IW,9370)IAMODE
            IF(LCOUNT.GE.10) WRITE(IW,9380)IAMODE
            TMP=ONE-SINTER
            CALL DSCAL(NC1,TMP,EVEC0,1)
            CALL DSCAL(NC1,SINTER,EVEC1,1)
            CALL VADD(EVEC0,1,EVEC1,1,EVEC1,1,NC1)
            TMP=DDOT(NC1,EVEC1,1,EVEC1,1)
            TMP=ONE/SQRT(TMP)
            CALL DSCAL(NC1,TMP,EVEC1,1)
            CALL DCOPY(NC1,EVEC1,1,EINT,1)
            ATMP=ABS(ALVEC(IAMODE))
            BTMP=ABS(BLVEC(IAMODE))
            IF (ATMP.LT.TMFIV) THEN
               IF (BTMP.LT.TMSIX) WRITE(IW,9400)
               IF (BTMP.GT.TM4) WRITE(IW,9410)
               IF (BTMP.LT.TM4 .AND. BTMP.GT.TMSIX) THEN
                  WRITE(IW,9420)
                  WRITE(IW,9430)
               END IF
            ELSE
               WRITE(IW,9430)
            END IF
C
C           IF THIS IS A BIFURCATION, THE WHOLE PLANE SPANNED
C           BY THE TWO FOLLOWING VECTORS IS A SOLUTION.
C           XI(IAMODE) = 0 : EVEC0
C           XI(IAMODE) = 1 : EVEC1
C           IF WE ARE FOLLOWING A SYMMETRIC GE AND HAVE A BIFURCATION,
C           EVEC0 IS THE SYMMETRIC SOLUTION AND
C           EVEC1 IS THE ASYMMETRIC SOLUTION (CROSSING GE)
C           IF WE HAVE A TURNING POINT, EVEC1 IS THE ACTUAL GE TANGENT
C           AND EVEC0 HAS NO MEANING.
C           THE CURRENT IMPLEMENTATION CANNOT SAFELY DISTINGUISH
C           BETWEEN A BIFURCATION OR TURNING POINT, THUS WE USE
C           THE INTERPOLATED GE TANGENT FOR THE NEXT PREDICTOR STEP.
C           SOME DAY WE MAY NEED BOTH EVEC0 AND EVEC1 ...
C           CURRENTLY WE CAN AT LEAST USE THEM FOR ANOTHER TEST
C           FOR WHETHER THIS IS A TURNING OR BIFURCATION POINT.
C           CONSTRUCT EVEC0
            CALL VCLR(TEMP,1,NC1)
            DO 780 I=1,NC1-NTR-1
               IF (ABS(ALVEC(I)).GT.EPS) THEN
                  TMP=BLVEC(I)/ALVEC(I)
               ELSE
                  TMP=ZERO
               END IF
               IF (I.EQ.IAMODE) TMP=ZERO
               DO 770 J=1,NC1-NTR-1
                  TEMP(J+1)=TEMP(J+1)+TMP*TMAT(J,I)
  770          CONTINUE
  780       CONTINUE
            TEMP(1)=ONE
            TMP=ONE/SQRT(DDOT(NC1-NTR,TEMP,1,TEMP,1))
            CALL DSCAL(NC1-NTR,TMP,TEMP,1)
            DO 800 I=1,NC1
               TMP=ZERO
               DO 790 J=1,NC1-NTR
                  TMP=TMP+UMAT(I,J)*TEMP(J)
  790          CONTINUE
               EVEC0(I)=TMP
  800       CONTINUE
C           EVEC1 IS GIVEN BY A-MAT EIGENVECTOR
            DO 810 I=1,NC1
               EVEC1(I)=TAMAT(I,IAMODE)
  810       CONTINUE
C
            TMP0=ABS(DDOT(NC1,EINT,1,EVEC0,1))
            TMP1=ABS(DDOT(NC1,EINT,1,EVEC1,1))
            WRITE(IW,9435)TMP0,TMP1
            IF (TMP0.GT.PT9 .OR. TMP1.GT.PT9) THEN
               IF (TMP0.GT.TMP1) THEN
                  WRITE(IW,9440)
               ELSE
                  WRITE(IW,9450)
               END IF
            ELSE
               WRITE(IW,9455)
            END IF
C
C           WRITE OUT ESTIMATE OF BIFURCATION/TURNING POINT GEOMETRY
C
            WRITE(IW,9360)
            CALL UPCOOR(DELX,NC1,OUT)
            CALL PRATM(IW,1)
            IF(NZMAT.GT.0) CALL PZANDG(EG,0)
C           GENERATE GEOMETRIES DISPLACED ALONG THE EVEC1 DIRECTION.
            CALL DCOPY(NC1,EVEC1,1,TEMP,1)
            CALL DSCAL(NC1,DXMAXT,TEMP,1)
            CALL UPCOOR(TEMP,NC1,OUT)
            WRITE(IW,9480)IAMODE,IPRO,DXMAXT
            CALL PRATM(IW,1)
            IF(NZMAT.GT.0) CALL PZANDG(EG,0)
            CALL DSCAL(NC1,-TWO,TEMP,1)
            CALL UPCOOR(TEMP,NC1,OUT)
            WRITE(IW,9480)-IAMODE,-IPRO,DXMAXT
            CALL PRATM(IW,1)
            IF(NZMAT.GT.0) CALL PZANDG(EG,0)
            CALL DSCAL(NC1,-PT5,TEMP,1)
            CALL UPCOOR(TEMP,NC1,OUT)
C
C           MAKE SELECTION OF WHICH GE TO FOLLOW
C           IF CROSSING GE, SUBSTITUTE TEMP FOR ENEW AND JUMP TO 870
C           CURRENTLY WE CONTINUE ALONG THE PREVIOUS DIRECTION
            WRITE(IW,9490)
            JNDEXA=IASAVE
            CALL DCOPY(NC1,AVEC,1,ALVEC,1)
            CALL DCOPY(NC1,BVEC,1,BLVEC,1)
            CALL DCOPY(NC1,EINT,1,ENEW,1)
            DELCOR=DELCCR
            LBIF=.FALSE.
            LCOUNT=0
            GO TO 870
         END IF
C
C     ROOT FINDING HAS NOT CONVERGED, CONSTRUCT NEW INTERPOLATION POINT
C
         IF (ALPH0*ALP0 .LT. ZERO) THEN
            ALP1=ALPH0
            SINT1=SINTER
         ELSE
            ALP0=ALPH0
            SINT0=SINTER
            DO 860 I=1,NC1
               DO 850 J=1,NC1
                  TPMAT(I,J)=TAMAT(I,J)
  850          CONTINUE
  860       CONTINUE
            JAMODE=IAMODE
         END IF
C
C        CONSTRUCT THE TWO CURRENT END POINT GEOMETRIES
         CALL DCOPY(NC1,XVEC,1,TEMP,1)
         CALL DSCAL(NC1,SINT0,TEMP,1)
         CALL VADD(OLDX,1,TEMP,1,TEMP,1,NC1)
         CALL DCOPY(NC1,TEMP,1,X,1)
         CALL DCOPY(NC1,XVEC,1,WORK,1)
         CALL DSCAL(NC1,SINT1,WORK,1)
         CALL VADD(OLDX,1,WORK,1,WORK,1,NC1)
         CALL VSUB(TEMP,1,WORK,1,TEMP,1,NC1)
         TMP=ALP0/(ALP0-ALP1)
         CALL DSCAL(NC1,TMP,TEMP,1)
         CALL UPCOOR(TEMP,NC1,OUT)
         GO TO 100
      END IF
C
C     END OF BIFURCATION/TURNING POINT LOCATION
C
C     MICRO ITERATION HAS CONVERGED
C     MOVE VARIOUS ARRAYS AND TAKE NEW PREDICTOR STEP
C
      JNDEXA=INDEXA
  870 CONTINUE
      DO 890 I=1,NC1
         AQVEC(I)=APVEC(I)
         APVEC(I)=ALVEC(I)
         BPVEC(I)=BLVEC(I)
         VOMODE(I)=UNRMAT(I,IFOLOW)
         EVEC(I)=ENEW(I)
         TEMP(I)=ENEW(I)
         OLDX(I)=X(I)
         DO 880 J=1,NC1
            UPMAT(I,J)=UNRMAT(I,J)
            TPMAT(I,J)=TAMAT(I,J)
  880    CONTINUE
  890 CONTINUE
      CALL DSCAL(NC1,DXMAXT,TEMP,1)
      CALL SYMDR(TEMP)
      CALL UPCOOR(TEMP,NC1,OUT)
      IF (NSERCH.GE.NSTEP) THEN
         WRITE(IW,9500)
         GO TO 900
      END IF
      GO TO 100
C
  900 CONTINUE
      RETURN
C
 9000 FORMAT(/10X,'WARNING! SCF CONVERGENCE IS TOO LOW FOR CHOSEN',
     *       ' SNUMH',/,10X,'CONVHF IS CHANGED TO',1P,E9.2)
 9010 FORMAT(1X,'STARTING GEOMETRY IS NOT A STATIONARY POINT, ',
     *       'CONSTRUCTING NEW GE TANGENT AND CORRECTION')
 9020 FORMAT(/,20X,'PREDICTOR GEOMETRY',/)
 9030 FORMAT(/,20X,'SCF IS UNCONVERGED, ABORTING RUN',/)
 9040 FORMAT(/10X,'NSERCH=',I3,5X,'ENERGY=',F16.7)
 9050 FORMAT(/30X,18(1H-),25X,23(1H-),/30X,18HCOORDINATES (BOHR), 25X,
     *  23HGRADIENT (HARTREE/BOHR)/30X,18(1H-),25X,23(1H-)// ,8X,
     *  13HATOM     ZNUC,7X,1HX,13X,1HY,13X,1HZ,20X,1HX,13X,1HY, 13X,
     *  1HZ/1X,112(1H-)/)
 9060 FORMAT(I3,2X,A8,A2,F6.1,3F14.7,8X,3F14.7)
 9065 FORMAT(//,5X,28(1H-),/,5X,
     $       '*** CALCULATING GRADIENT ***',
     $       /,5X,28(1H-),//)
 9066 FORMAT(//,5X,28(1H-),/,5X,
     $       '*** END OF      GRADIENT ***',
     $       /,5X,28(1H-),//)
 9067 FORMAT(//,5X,27(1H-),/,5X,
     $       '*** CALCULATING HESSIAN ***',
     $       /,5X,27(1H-),//)
 9068 FORMAT(//,5X,27(1H-),/,5X,
     $       '*** END OF      HESSIAN ***',
     $       /,5X,27(1H-),//)
 9070 FORMAT(//,5X,44(1H-),/,5X,
     $       '*** CALCULATING PARTIAL THIRD DERIVATIVE ***',
     $       /,5X,44(1H-),//)
 9080 FORMAT(//,5X,44(1H-),/,5X,
     $       '*** END OF      PARTIAL THIRD DERIVATIVE ***',
     $       /,5X,44(1H-),//)
 9090 FORMAT(1X,'PURE NR STEP HAS LENGTH',F12.6,/)
 9100 FORMAT(1X,'INSIDE NR REGION, SWITCHING TO PURE NR SEARCH')
 9110 FORMAT(1X,'ERROR! NOT 5 OR 6 T+R MODES',/,
     *       1X,'THE GEOMETRY MAY BE NEARLY LINEAR OR DISSOCIATING',
     *     /,1X,'THE CURRENT IMPLEMENTATION CANNOT HANDLE THIS!')
 9120 FORMAT(1X,'WARNING! AA-DAGGER MATRIX NEARLY SINGULAR,',
     *          ' LOWEST EIGENVALUE =',D12.3)
 9130 FORMAT(1X,'REMOVING CORRECTOR STEP COMPONENT',I4,/,
     *       1X,'EIGENVALUE =',D12.3,1X,'F-VEC =',D12.3,1X,
     *          'RESULTING CORRECTOR COMPONENT =',D12.3,/)
 9140 FORMAT(1X,'A-MATRIX EIGENVALUES AND B-VECTOR')
 9150 FORMAT(I5,2F15.6)
 9160 FORMAT(1X,'GE INDEX =',I4)
 9170 FORMAT(1X,'GE TANGENT IN INTERNAL COORDINATES')
 9175 FORMAT(20F10.5)
 9180 FORMAT(1X,'INITIALIZING GE RUN, IFOLOW =',I5)
 9200 FORMAT(/1X,'PROJECTIONS ALONG CURRENT HESSIAN EIGENVECTORS',
     *' AND THEIR CURVATURE'/
     *14X,'EIGENVALUE',4X,'GRADIENT',4X,'GE-TANGENT',4X,'E-VECTOR',
     * 4X,'E-CURVATURE')
 9210 FORMAT(5X,I5,6F13.6)
 9220 FORMAT(/,1X,'INFORMATION: ENERGY MODE SWITCHING, FROM',
     *       I5,3X,'TO',I5)
 9230 FORMAT(/,20X,'CORRECTOR GEOMETRY',/)
 9240 FORMAT(/,5X,'MACRO',I5,3X,'MICRO',I5,3X,
     *          'DNORM',F10.6,3X,'GE TAN ANGLE',F6.2,3X,'ENERGY',F15.8)
 9250 FORMAT(1X,'TO MANY MICRO ITERATIONS')
 9260 FORMAT(5X,'CORRECTION STEP',F10.6,' LARGER THAN LIMIT',F10.6,
     *          ' CONSTRUCTING IMPROVED GEOMETRY',/)
 9261 FORMAT(5X,'CORRECTION STEP',F10.6,' BELOW LIMIT',F10.6,
     *          ' TAKING NEW PREDICTOR STEP',/)
 9270 FORMAT(5X,'WARNING!! ALGORITHM SEEMS TO BE OFF TRACK',/,
     *5X,'LARGEST GRADIENT COMPONENT IS NOT ALONG THE',
     *' EIGENVECTOR WHICH HAS THE LARGEST OVERLAP')
 9280 FORMAT(//,5X,'CYCLE',I5,3X,'ENERGY',F15.8,3X,'GNORM',F10.6,3X,
     $          'GE CURVATURE',F8.3,/,
     $          5X,'CYCLE',I5,3X,'ANGLE BETWEEN GE TANGENTS',F6.2,/,
     $          5X,'CYCLE',I5,3X,'PREDICTOR',F10.6,3X,'CORRECTOR',
     $          F10.6,3X,'TOTAL',F10.6,/,
     $          18X,'ANGLE BETWEEN GE TANGENT AND GRADIENT',F6.2,//)
 9290 FORMAT(1X,'ALPHA MODE',I3,2X,'IS OLD MODE',I3,2X,
     *             'WITH OVERLAP',F8.5)
 9300 FORMAT(/,1X,'INFORMATION: ALPHA NORM HAS PASSED A MINIMUM',
     *       /,1X,'ALPHA MODE',I3,1X,'LAST THREE VALUES',3F12.6,1X,
     *       'BETA',F12.6,/)
 9310 FORMAT(1X,'GE IS CHANGING INDEX WITH',I3,2X,'FROM',I3,
     *       2X,'TO',I3)
 9320 FORMAT(1X,'CANNOT HANDLE MULTIPLE INDEX CHANGES AT PRESENT',
     *     /,1X,'TRY TO RESTART FROM THE PREVIOUS POINT WITH A',
     *      ' SMALLER STEPSIZE',//,1X,
     *       '*** INTERPOLATION WILL TRY TO LOCATE ONE ROOT ***',/)
 9330 FORMAT(1X,'SEARCHING FOR INDEX CHANGE FOR CURRENT MODE',
     *          ' NUMBER',I3,
     *          ' ALONG LINE CONNECTING CURRENT POINT WITH PREVIOUS')
 9340 FORMAT(1X,'INTERPOLATION END POINTS:')
 9350 FORMAT(/,1X,'SINTER = ',F8.5,3X,'ALPHA =',F10.6,5X,
     *             'BETA =',F10.6)
 9360 FORMAT(/,1X,'ADDING CORRECTION STEP TO PROVIDE BEST ESTIMATE',
     *           ' OF BIFURCATION/TURNING POINT GEOMETRY',/)
 9370 FORMAT(1X,'ZERO OF ALPHA FOUND FOR ALPHA MODE =',I5)
 9380 FORMAT(1X,'ROOT FINDING FAILED, ABORTING, ALPHA MODE =',I5)
 9400 FORMAT(/,1X,'VALUE OF BETA INDICATE A GE BIFURCATION POINT',/)
 9410 FORMAT(/,1X,'VALUE OF BETA INDICATE A GE TURNING POINT',/)
 9420 FORMAT(/,1X,'VALUE OF BETA IS INCONCLUSIVE',/)
 9430 FORMAT(1X,'PLEASE MAKE A MANUAL DECISSION WHETHER THIS IS',
     *          ' A BIFURCATION OR TURNING POINT',/,
     *       1X,'OR TRY TO REFINE THE INTERPOLATION WITH A SMALLER',
     *          ' STEPSIZE AND/OR TIGHTER FOLLOWING')
 9435 FORMAT(1X,'OVERLAPS OF INTERPOLATED GE TANGENT WITH LIMITING',
     *           ' GE TANGENT VECTORS ARE ',2F7.3)
 9440 FORMAT(1X,'OVERLAP INDICATES A GE BIFURCATION POINT',/)
 9450 FORMAT(1X,'OVERLAP INDICATES A GE TURNING POINT',/)
 9455 FORMAT(1X,'OVERLAP IS INCONCLUSIVE',
     *           ' WHETHER THIS IS A BIFURCATION OR TURNING POINT',/)
 9480 FORMAT(/,1X,'POSSIBLE CROSSING GE, ALPHA MODE =',I3,2X,
     *          'IFOLOW =',I3,2X,'DISPLACEMENT =',F8.5,/)
 9490 FORMAT(/,1X,'RETURNING TO ORIGINAL GE PATH',/,
     *       1X,'TAKING PREDICTOR STEP FROM BIFURCATION/TURNING',
     *          ' POINT WITH INTERPOLATED GE TANGENT')
 9500 FORMAT(5X,'EXCESS NUMBER OF STEPS')
      END
