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 14 MAR 05 - HL  - ADD CHARGE TRANSFER ARG TO EFDSPL, EFSAVE:DYN.MEMORY
C  7 MAR 05 - IA  - EFGEOM: DYN.POL. ARG FOR EFDSPL, EFSAVE: DYN.POL.
C 13 FEB 05 - JS  - ADD EFP MP2 CHANGES TO EFGEOM
C  5 FEB 05 - MWS - PAD COMMON MCINP
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 26 MAR 03 - JMR - PUT TINKER COORDS IN IRC FILE
C 14 JAN 03 - JI  - ADD HOOK TO ORMAS CI STEP
C  7 AUG 02 - HT  - USE DYNAMIC MEMORY FOR EFP PAULI REPULSION
C 17 APR 02 - MWS - EFRMSS: SUPPRESS SOME PRINTING
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C  1 AUG 01 - JI  - TWEAKS FOR GENERAL CI
C 25 JUN 01 - MWS - ALTER COMMON BLOCK SCFOPT
C 20 FEB 01 - MWS - PAD PAULMO COMMON
C 11 OCT 00 - MWS - UPDATE THE DETWFN COMMON
C 25 MAR 00 - CHC - FIX QMMM IRC RUN
C 21 DEC 99 - MWS - MAKE GUGWFN COMMON CONSISTENT
C 29 AUG 99 - CHC - PATCHES FOR QMMM
C 13 FEB 99 - MAF - IRCMOS: PRINT/PUNCH ONLY L0 ORBITALS AT MOST
C 12 NOV 98 - MWS - KEEP COMMON MCINP CONSISTENT
C 22 APR 98 - MWS - IRCMOS: COUNT MCSCF OCCUPIED ORBITALS CORRECTLY
C 13 APR 98 - MWS - IRCMOS: READ SYMMETRY LABELS FROM DAF
C  6 JAN 98 - JAB - ADD A VARIABLE TO ZMAT COMMON BLOCK
C 16 JUL 97 - GNM - INIIRC,IRCGS2,IRCX: CHANGES FOR FRAGONLY RUNS
C 18 DEC 96 - MWS - EFGEOM: UPDATE CALL TO EFDSPL
C 17 OCT 96 - SPW - IRCGS2,IRCSTP,STBIMK: ADD CALL TO CIGRAD FOR CI RUNS
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C 24 MAY 96 - WC  - REMOVE EFPEX IN COMMOM EFPPAR
C  1 APR 96 - MWS - IRCGS2: TRIPLE THE RTOL TEST COMPARED TO 29 JAN 96
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C  9 FEB 96 - WC  - FRAGMENT DECLARATIONS CORRECTED
C 29 JAN 96 - MWS - IRCGS2: MAKE RTOL TEST PROPORTIONAL TO NO. OF ATOMS
C  9 JAN 96 - WC  - ADD FRAGMENT OUTPUT TO FILE IRCDATA
C 15 NOV 95 - FRJ - PROJECT HESSIAN IN IRCGS2 FOR CARTESIAN,
C                   IMPROVE LAMBDA BRACKETING IN IRCGS2
C 28 SEP 95 - MWS - IRCGS2: PRINT INTERNUCLEAR DISTANCES FOR IRC POINTS
C 10 MAY 95 - JAB - IRCX: EXPAND FORMAT WIDTHS
C 21 APR 95 - MWS - INIIRC,IRCGS2,IRCSTP,IRCX: PRINT ALL+UNIQUE COORDS
C  1 FEB 95 - MWS - IRCGS2: RESTORE TERMINATION MESSAGE 9350 TO USE
C  5 JAN 95 - PND - EFSAVE: PAD COMMON EFPPAR
C 28 NOV 94 - WC  - IRCGS2: FIX RESTART INFO FOR FRAGMENTS
C 17 NOV 94 - WC  - IMPLEMENT IRC FOR EFFECTIVE FRAGMENTS
C 11 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 26 MAY 93 - JHJ - FIX BUG WITH PARALLEL EXEC OF IRC RESTART PUNCH
C 19 APR 93 - SJS - PUNCH CURRENT IRC POINT (NOT NEXT CIRCLE POINT)
C 16 MAR 93 - MWS - IMPROVE IRC RESTART INFO, MAKE GS2 THE DEFAULT
C 12 DEC 92 - MWS - IRCGS2: FIX PARALLEL GS2 BUG, SHORTEN LINES
C 24 NOV 92 - MWS - IRCSTP,SEABIS,STBIMK: PAD OUT IRCDAT COMMON
C 15 APR 92 - SJS - GONZALEZ-SCHLEGEL 2ND-ORDER ALGORITHM INTEGRATED
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C 10 JAN 92 - MWS,TLW - PARALLEL I/O CHANGES
C 13 NOV 91 - MWS - IRCX: PUNCH RESTART COORDINATES IN ANGSTROMS.
C 13 FEB 91 - TLW - SET STABLZ .FALSE. IF PACE IS NOT LINEAR
C 15 AUG 90 -  FJ - PILOT CODE FOR SCHLEGEL METHOD FOR IRC FOLLOWING.
C  9 DEC 89 - MWS - CORRECT NUCLEAR CHARGES FOR ECP CASES.
C 23 OCT 89 - MWS - DELETE /FUNCT/ SAVES
C 20 FEB 89 - STE - IRCSTP: SIMPLIFY STATEMENT LABELS A BIT
C  2 FEB 89 - MWS - PASS NUCLEAR MASSES VIA NEW COMMON /MASSES/
C 12 JAN 89 - MWS - USE FCMIN INSTEAD OF RDHES
C 02 DEC 88 - JAB - PRINT ANGLE BETWEEN SUCCESSIVE GRADIENTS FOR LINEAR
C 20 NOV 88 - MWS - TURN OFF STBSKP IN IRCSTP
C 25 AUG 88 - MWS - IMPROVE PRINTOUT WHEN ELBOW SKIPS STABILIZATION,
C                   ABORT JOB IF ENERGY RISES ON THE IRC.
C  8 AUG 88 - MWS - SKIP E,G AT PREDICTED RK4 POINT. FIX RESTART PUNCH
C 22 JUL 88 - KIM - ADDED CODE FOR HIGHER ORDER IRC CALCULATION.
C  4 AUG 87 - MWS - FIX ERROR IN INPUT OF CMODE
C 15 MAY 87 - MWS(STE) - DIAG SCR TO 8 VECTORS IN INIIRC,INIIRM
C 31 JAN 87 - MWS - ADD OPTION FOR ELBOW TEST, REORGANIZE LOGIC,
C                   IMPROVE RESTART PUNCHOUT, DELETE CUBIC STUFF.
C  4 NOV 86 - STE - USE PARAMETERS, DELETE 9020 IN IRCSTP
C 14 AUG 86 - MWS - ALPHABETIZE ROUTINES, ROUTINE OPNIRC CALLED
C  7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
C 17 JUN 86 - MWS - CHANGE CONVERSION FACTOR TO WAVENUMBERS,
C                   DEFAULT FOR STABLZ MADE .FALSE.
C  2 MAY 86 - MWS - ADD QUADRATIC PACE (AND CUBIC STUFF)
C 23 APR 86 - MWS - ALLOW PUNCH OF IRC INFO TO SEPARATE FILE
C 11 APR 86 - MWS - REACTION PATH STABILIZATION MADE OPTIONAL
C  3 DEC 85 - STE - USE GENERIC ABS,SQRT; CLEAN UP LOGIC
C 11 JUL 85 - MWS - FIX A FORMAT STATEMENT
C 26 FEB 85 - MWS - REMOVE LOGICAL*1 INPUT VARIABLES, STRIDE TO 0.15,
C                   SWITCH TO PZANDG
C 31 JAN 85 - MWS - NO. OF INTERNALS FIXED FOR LINEAR BEND CASE
C  6 DEC 84 - MWS - CHANGE MASS WEIGHTING OF STEP FROM QB
C 28 OCT 83 - MWS - INCORPORATE MICHEL'S IDEAS ON THE LINEAR SEARCH
C 29 AUG 83 - MWS - FIX ERROR IN CHOICE OF SEARCH DIRECTION
C 17 AUG 83 - MWS - FIX UP PRINTOUT SOME
C  3 JUL 83 - MWS - ALLOW INPUT OF THE IMAGINARY MODE
C 24 JUN 83 - MWS - IMPLEMENT IRC IN GA0ESS
C*MODULE RXNCRD  *DECK EFDST
      SUBROUTINE EFDST(NFRG,QBM,QCM,SAC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION QBM(*),QCM(*)
      DATA ZERO/0.0D+00/
C
      IF (NFRG.EQ.0) RETURN
C
      SAC=SAC**2
      JDX=0
      DO 10 I=1,NFRG
      DEF=ZERO
      DO 20 K=1,6
   20 DEF=DEF+(QBM(JDX+K)+QCM(JDX+K))**2
      SAC=SAC+DEF
      JDX=JDX+6
   10 CONTINUE
      SAC=SQRT(SAC)
C
      RETURN
      END
C*MODULE RXNCRD  *DECK EFGEOM
      SUBROUTINE EFGEOM(NC,NCF,A,B,RM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION A(*),B(*),RM(*)
C
      PARAMETER (MXFRG=1050, MXPT=2000, MXFGPT=12000,
     *           MXDFG=5, MXAO=8192, MXSH=5000, MXGTOT=20000,
     *           MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /EFPBAS/ EX(MXGTOT,MXDFG),CS(MXGTOT,MXDFG),
     *                CP(MXGTOT,MXDFG),CD(MXGTOT,MXDFG),
     *                CF(MXGTOT,MXDFG),CG(MXGTOT,MXDFG),
     *                PRNAME(MXFGPT),PRCORD(3,MXFGPT),EFZNUC(MXFGPT),
     *                KSTART(MXSH,MXDFG),KATOM(MXSH,MXDFG),
     *                KTYPE(MXSH,MXDFG),KNG(MXSH,MXDFG),
     *                KLOC(MXSH,MXDFG),KMIN(MXSH,MXDFG),
     *                KMAX(MXSH,MXDFG),NSHELL(MXDFG),NGAUSS(MXDFG),
     *                NATEF(MXDFG),NUMEF(MXDFG),NTPATM
      COMMON /EFPCMO/ NOCC(MXDFG),NVIR(MXDFG),NAO(MXDFG),NTAO,
     *                CANFOK(MXAO,MXDFG)
      COMMON /FMCOM / XX(1)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGTYP/ NDFRG,ISET(MXFRG),NAMIDX(MXDFG)
      COMMON /PAULMO/ CENTNM(MXFGPT),CENTCD(3,MXFGPT),NORB(MXFRG),
     *                NPBF(MXFRG),NTMO
      COMMON /POLTMP/ POLT2(9),
     *                ACHGHF(2,MXPT),ADIPHF(3,MXPT),AQADHF(6,MXPT),
     *                AMULTMP(MXPT),BMULTMP(MXPT),
     *                AMULTMP2(MXPT),BMULTMP2(MXPT),
     *                TCLPD(4*MXPT),TZLPD(4*MXPT),
     *                JNLPR(4*MXPT),JKFR(MXPT),JKLR(MXPT),
     *                NDPTS(MXFRG),
     *                IEFP2,IEFMON,IEFSCR,IEFSCR2,IEFDIP,IEFQUD,IEFREP,
     *                IEFPOL,IEFPOL2,NDTTPT
C
      IF (NFRG.EQ.0) RETURN
C
C     GET CART DISPLACEMENT FOR ROTATION
C
      CALL EFGRDM(NC,NCF,NFRG,A,RM,B,0)
C
      INDX = NC
      NMST = 0
      NPST = 0
      NDPST = 0
      NRST = 0
      NDST = 0
      NATST = 0
      NORST = 0
      NAOST = 0
C
      MXBF=0
      DO I = 1, NFRG
         MXBF=MAX(MXBF,NPBF(I))
      END DO
      CALL VALFM(LOADFM)
      LPROVEC=LOADFM+1
      LCANVEC=LPROVEC+MXBF*NTMO
      LAST=LCANVEC+MXBF*NTAO
      NEED=LAST-LOADFM-1
      CALL GETFM(NEED)
C
      LENPV=MXBF*NTMO
      LENCV=MXBF*NTAO
      IF (NFRG.GT.0.AND.NTMO.GT.0)
     *   CALL DAREAD(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
      IF (NFRG.GT.0.AND.NTAO.GT.0)
     *   CALL DAREAD(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
C
      DO 100 I=1,NFRG
        CALL EFDSPL(I,B(INDX+1),B(INDX+2),B(INDX+3),B(INDX+4),B(INDX+5),
     *              B(INDX+6),NMST,NPST,NDPST,NRST,NDST,NATST,NORST,
     *              XX(LPROVEC),NAOST,XX(LCANVEC),MXBF)
C
        INDX = INDX + 6
        NMST  = NMST  + NMPTS(I)
        NPST  = NPST  + NPPTS(I)
        NDPST = NDPST + NDPPTS(I)
        NRST  = NRST  + NRPTS(I)
        NDST  = NDST  + NDPTS(I)
        NATST = NATST + NATEF(ISET(I))
        NORST = NORST + NORB(I)
        NAOST = NAOST + NAO(ISET(I))
  100 CONTINUE
C
      IF (NFRG.GT.0.AND.NTMO.GT.0)
     *     CALL DAWRIT(IDAF,IODA,XX(LPROVEC),LENPV,267,0)
      IF (NFRG.GT.0.AND.NTAO.GT.0)
     *     CALL DAWRIT(IDAF,IODA,XX(LCANVEC),LENCV,269,0)
      CALL RETFM(NEED)
C
C     GET NEW FRGMI ETC.
C
      CALL EFRMSS(NFRG,RM(NC+1))
C
      RETURN
      END
C*MODULE RXNCRD  *DECK EFRMSS
      SUBROUTINE EFRMSS(NFRG,RM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,MASWRK,DSKWRK,SOME
      DIMENSION RM(*)
      DIMENSION TEMP(6),SCR(24),EIG(3),VEC(3,3),IWORK(3)
      CHARACTER*6 FRGNAM
      CHARACTER*8 PTNAM
      PARAMETER (MXFRG=1050, MXFGPT=12000)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, SMALL=1.0D-08)
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /FRAGMT/ XCRD(3,MXFRG),YCRD(3,MXFRG),ZCRD(3,MXFRG),
     *                PTNAM(3,MXFRG),FRGNAM(MXFRG)
C
      IF (NFRG.EQ.0) RETURN
      SOME = MASWRK  .AND.  NPRINT.NE.-5
C
C     THIS ROUTINE RETURNS THREE PIECES OF INFORMATION:
C     A) RECIPROCAL TOTAL MASSES OF EACH FRAGMENT
C        ELEMENTS 1 TO  3*NFRG
C     B) MOMENTS OF INERTIA IN FRAG'S ORIGINAL COORDINATES
C          3*NFRG+1 TO 12*NFRG
C     C) MOMENTS OF INERTIA IN FRAG'S LOCAL PRINCIPAL AXES
C         12*NFRG+1 TO 21*NFRG
C     THEREFORE, TOTAL STORAGE NEEDED FOR -RM- IS 21*NFRG.
C
      CALL EFCM
      IF(SOME) WRITE(IW,9230)
      I0=0
      DO 100 INF=1,NFRG
        IF(SOME) WRITE(IW,9240) INF,FRGNAM(INF),FRGMAS(INF)
        DMF=ONE/SQRT(FRGMAS(INF))
        RM(I0+1)=DMF
        RM(I0+2)=DMF
        RM(I0+3)=DMF
        I0=I0+3
  100 CONTINUE
C
      IF(SOME) WRITE(IW,9250)
      L=3*NFRG+1
      M=12*NFRG+1
      DO 300 INF=1,NFRG
        IF(SOME) WRITE(IW,9260) INF,FRGNAM(INF),
     *               FRGMI(1,INF),FRGMI(2,INF),FRGMI(3,INF),
     *               FRGMI(4,INF),FRGMI(5,INF),FRGMI(6,INF)
        CALL DCOPY(6,FRGMI(1,INF),1,TEMP,1)
        CALL GLDIAG(3,3,3,TEMP,SCR,EIG,VEC,IGERR,IWORK)
        DO 101 J=1,3
          IF (ABS(EIG(J)).LT.SMALL) THEN
            TEMP(J)=ZERO
            TEMP(J+3)=ZERO
          ELSE IF (EIG(J).GE.SMALL) THEN
            TEMP(J+3)=SQRT(EIG(J))
            TEMP(J)=ONE/TEMP(J+3)
          ELSE
            WRITE(IW,9270)
            WRITE(IW,9280) INF,J,EIG(J)
            CALL ABRT
          END IF
  101   CONTINUE
        DO 230 I=1,3
          DO 220 J=1,3
            RM(L) = ZERO
            RM(M) = ZERO
            DO 210 K=1,3
              RM(L) = RM(L) + VEC(I,K)*VEC(J,K)*TEMP(K)
              RM(M) = RM(M) + VEC(I,K)*VEC(J,K)*TEMP(K+3)
  210       CONTINUE
            L=L+1
            M=M+1
  220     CONTINUE
  230   CONTINUE
  300 CONTINUE
C
 9230 FORMAT(/1X,'FRAGMENT MASSES USED ARE (IN AMU)')
 9250 FORMAT(/1X,'FRAGMENT MOMENTS OF INERTIA USED (LOW TRIANGLE)',
     *       ' ARE (IN AMU)')
 9240 FORMAT(3X,I2,2X,A6,1X,F11.6)
 9260 FORMAT(3X,I2,2X,A6,1X,6F11.6)
 9270 FORMAT(' MOMENTS OF INERTIA ON PRINCIPLE AXES CAN NOT BE',
     *       ' NEGITIVE')
 9280 FORMAT(' IFRG =',I3,2X,'J =',F11.6)
C
      RETURN
      END
C*MODULE RXNCRD  *DECK EFSAVE
      SUBROUTINE EFSAVE(IOPT,CSAVE,DSAVE,QSAVE,OSAVE,PSAVE,SAVEP,
     *                  DPSAVE,SAVEDP,RSAVE,SAVEC,SAVEM)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
C
      PARAMETER (MXFRG=1050, MXFGPT=12000, MXIFRQ=12,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      CHARACTER*8 FRGNME,POLNAM,DPOLNAM,REPNAM
C
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /EFPPAR/ EFP(3,MXFGPT),EFPOL(9,MXFGPT),
     *                EFDP(3,MXIFRQ*MXFGPT),EFDPOL(9,MXIFRQ*MXFGPT),
     *                ENO,DIND(3,MXFGPT),DINDD(3,MXFGPT),
     *                POLNAM(MXFGPT),DPOLNAM(MXIFRQ*MXFGPT)
      COMMON /REPPAR/ CREP(3,MXFGPT),CLPR(4*MXFGPT),ZLPR(4*MXFGPT),
     *                NLPR(4*MXFGPT),KFR(MXFGPT),KLR(MXFGPT),
     *                REPNAM(MXFGPT)
      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
C
      IF (NFRG.EQ.0) RETURN
C
C         SAVE EFFECTIVE FRAGMENT POTENTIAL TO A COMMON BLOCK
C               WHAT A WASTE OF MEMORY THIS TURKEY IS!
C
      NFRG3M  = 3*NMTTPT
      NFRG3P  = 3*NPTTPT
      NFRG3DP = 3*NDPTTPT*MXIFRQ
      NFRG3R  = 3*NRTTPT
      NFRG6M  = 2*NFRG3M
      NFRG9P  = 3*NFRG3P
      NFRG9DP = 9*NDPTTPT*MXIFRQ
      NFRG10M = 10*NMTTPT
      NFRG3=NFRG*3
      NFRG6=NFRG*6
C
      IF (IOPT.EQ.0) THEN
         CALL DCOPY(NFRG3M,EFC,1,CSAVE,1)
         CALL DCOPY(NFRG3M,EFDIP,1,DSAVE,1)
         CALL DCOPY(NFRG6M,EFQAD,1,QSAVE,1)
         CALL DCOPY(NFRG10M,EFOCT,1,OSAVE,1)
         CALL DCOPY(NFRG3P,EFP,1,PSAVE,1)
         CALL DCOPY(NFRG9P,EFPOL,1,SAVEP,1)
         CALL DCOPY(NFRG3DP,EFDP,1,DPSAVE,1)
         CALL DCOPY(NFRG9DP,EFDPOL,1,SAVEDP,1)
         CALL DCOPY(NFRG3R,CREP,1,RSAVE,1)
         CALL DCOPY(NFRG3,EFCENT,1,SAVEC,1)
         CALL DCOPY(NFRG6,FRGMI,1,SAVEM,1)
      ELSE
         CALL DCOPY(NFRG3M,CSAVE,1,EFC,1)
         CALL DCOPY(NFRG3M,DSAVE,1,EFDIP,1)
         CALL DCOPY(NFRG6M,QSAVE,1,EFQAD,1)
         CALL DCOPY(NFRG10M,OSAVE,1,EFOCT,1)
         CALL DCOPY(NFRG3P,PSAVE,1,EFP,1)
         CALL DCOPY(NFRG9P,SAVEP,1,EFPOL,1)
         CALL DCOPY(NFRG3DP,DPSAVE,1,EFDP,1)
         CALL DCOPY(NFRG9DP,SAVEDP,1,EFDPOL,1)
         CALL DCOPY(NFRG3R,RSAVE,1,CREP,1)
         CALL DCOPY(NFRG3,SAVEC,1,EFCENT,1)
         CALL DCOPY(NFRG6,SAVEM,1,FRGMI,1)
      END IF
C
      RETURN
      END
C*MODULE RXNCRD  *DECK INIIRC
      SUBROUTINE INIIRC(EVIB,FORWRD,NCOORD,NCF,NC2,FREQ,CMODE,NHOPT,
     *                  HESS,A,IA,SCR,VEC,E)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION K0,K1,K2,K3
C
      LOGICAL FORWRD,STBSKP,GOTEH,MASWRK,GOPARR,DSKWRK,GSCHL2,FGONLY
C
      DIMENSION HESS(NCOORD+NCF,NCOORD+NCF),A(NC2),E(NCOORD+NCF),
     *          IA(NCOORD+NCF),SCR(NCOORD+NCF,8),
     *          VEC(NCOORD+NCF,NCOORD+NCF),CMODE(NCOORD+NCF)
C
      PARAMETER (MXATM=2000, MXA3=3*MXATM,
     *           MXFRG=1050, MXNCF=6*MXFRG,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      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 /IRCDAT/ EA,EB,EC,SAB,SBC,STOTAL,GANORM,GBNORM,GCNORM,
     *                GBIS,GCUT,RCUT,ACUT,STRIDE,OPTTOL,
     *                QA(MXA3+MXNCF),QB(MXA3+MXNCF),QC(MXA3+MXNCF),
     *                GA(MXA3+MXNCF),GB(MXA3+MXNCF),GC(MXA3+MXNCF),
     *                QCM(MXA3+MXNCF),QBM(MXA3+MXNCF),QAM(MXA3+MXNCF),
     *                GA2(MXA3+MXNCF),GA1(MXA3+MXNCF),GA0(MXA3+MXNCF),
     *                QCNM(MXA3+MXNCF),QCN(MXA3+MXNCF),TGP(MXA3+MXNCF),
     *                BIS(MXA3+MXNCF),EBIS(10),XBIS(10),
     *                RTRMS(MXA3+MXFRG*21),K0(MXA3+MXNCF),
     *                K1(MXA3+MXNCF),K2(MXA3+MXNCF),K3(MXA3+MXNCF),
     *                MXOPT,IHUPD,IPIRC,NPIRC,NPRT,NPUN,STBSKP,GSCHL2
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, SMALL=1.0D-05
     *          ,FACT=2.642461D+07, HALF=0.5D+00)
C
C     ----- THIS ROUTINE GENERATES AN INITIAL POINT WITH A NONZERO
C           GRADIENT VECTOR BY FOLLOWING THE IMAGINARY NORMAL MODE
C           A SHORT DISTANCE AWAY FROM THE SADDLE POINT -----
C
      IF(MASWRK) WRITE(IW,9000)
C
C     ----- STORE SADDLE POINT COORDS IN QA, WITH ZERO GRADIENT ----
C
      FGONLY = NAT.EQ.0  .AND.  NFRG.GT.0
      NC1 = NCOORD + NCF
      NC3 = NC1 * NC1
      NCP1 = NCOORD + 1
      STOTAL = ZERO
      CALL DCOPY(NCOORD,C,1,QA,1)
      CALL VCLR(GA,1,NC1)
      IF(NHOPT.EQ.0) GO TO 200
C
C     ----- OBTAIN IMAGINARY MODE FROM THE INPUT FORCE FIELD -----
C
C     ----- READ FORCE CONSTANT MATRIX, WHICH IS RETURNED
C           IN CARTESIAN COORDINATES -----
C
      CALL FCMIN(HESS,NC1,GOTEH)
      IF(.NOT.GOTEH) THEN
         IF(MASWRK) WRITE(IW,9010)
         CALL ABRT
      END IF
C
      CALL DAWRIT(IDAF,IODA,HESS,NC3,4,0)
C
C     ----- MASS WEIGHT THE HESSIAN, STORE IN SYMMETRIC FORM -----
C
      CALL HESMAS(NCOORD,NCF,NFRG,HESS,RTRMS,A,0)
C
C     ----- OBTAIN NORMAL MODES AND FREQUENCIES -----
C
      IERR = 0
      CALL GLDIAG(NC1,NC1,NC1,A,SCR,E,VEC,IERR,IA)
      IF (IERR .NE. 0) CALL ABRT
      EONE = E(1)
      IF(EONE.GE.ZERO) THEN
         IF(MASWRK) WRITE(IW,9020) EONE
         CALL ABRT
      END IF
      FREQ = SQRT(-FACT*EONE)
C
C     ----- CONVERT THE IMAGINARY MODE FROM MASS WEIGHTED
C           TO ORDINARY CARTESIAN COORDINATES -----
C
      DO 150 I=1,NCOORD
         CMODE(I) = RTRMS(I) * VEC(I,1)
  150 CONTINUE
C     NOTE CMODE FOR FRAG IN M-W.
      CALL DCOPY(NCF,VEC(NCP1,1),1,CMODE(NCP1),1)
      GO TO 300
C
C     ----- NORMALIZE THE INPUT NORMAL MODE SO ITS LENGTH
C           IN MASS WEIGHTED CARTESIAN SPACE IS UNITY -----
C
  200 CONTINUE
      EONE = -FREQ*FREQ/FACT
      IF(EONE.GE.ZERO) CALL ABRT
      VNORM = ZERO
      DO 220 I=1,NCOORD
         V = CMODE(I)
         V = V/RTRMS(I)
         VNORM = VNORM + V * V
  220 CONTINUE
      VNORM = SQRT(VNORM)
      IF(VNORM.LT.SMALL) CALL ABRT
      VNORM = ONE/VNORM
      IF(ABS(ONE-VNORM).GT.SMALL .AND. MASWRK) WRITE(IW,9100) VNORM
      DO 240 I=1,NCOORD
         CMODE(I) = VNORM * CMODE(I)
  240 CONTINUE
C
C     ----- AT THIS POINT, CMODE IS IN CARTESIAN COORDINATES,
C           WITH UNIT LENGTH IN MASS-WEIGHTED CARTESIAN SPACE.
C           EONE IS THE CORRESPONDING NEGATIVE CURVATURE. -----
C
C     ----- PICK DIRECTION ALONG THE IMAGINARY MODE,
C           AND THE MAGNITUDE OF THE HARMONIC STEP -----
C
  300 CONTINUE
      IF(GSCHL2) THEN
         STEP = HALF*STRIDE
      ELSE
         STEP = SQRT(-TWO * EVIB/EONE)
      END IF
C
C     FOLLOW MAX MODE IN AB INITIO PART.
C
      IF (.NOT. FGONLY) THEN
         IMAX = IDAMAX(NCOORD,CMODE,1)
         IF(     FORWRD  .AND.  CMODE(IMAX).LT.ZERO) STEP = -STEP
         IF(.NOT.FORWRD  .AND.  CMODE(IMAX).GT.ZERO) STEP = -STEP
      ELSE
         IF (.NOT. FORWRD) STEP = -STEP
      END IF
C
C     ----- FORM A SYMMETRIZED STEP ALONG THE NORMAL MODE -----
C
      CALL DSCAL(NC1,STEP,CMODE,1)
      IF (NFRG.EQ.0) CALL SYMDR(CMODE)
      IF(MASWRK) THEN
         WRITE(IW,9040) STEP,FREQ
         WRITE(IW,9060) (CMODE(I),I=1,NCOORD)
         IF (NFRG.GT.0) WRITE(IW,9070) (CMODE(I),I=NCP1,NC1)
      END IF
C
C     ----- GENERATE CARTESIAN COORDINATES OF THE POINT THAT
C           IS DISPLACED FROM THE SADDLE POINT -----
C
      CALL VADD(CMODE,1,QA,1,QB,1,NCOORD)
      IF(GSCHL2) CALL VADD(CMODE,1,QB,1,QC,1,NCOORD)
C     NOTE THAT START AT CMODE(NC+1), THEY ARE IN MASS-WEIGHTED
      CALL DCOPY(NCF,CMODE(NCP1),1,QBM(NCP1),1)
      IF(GSCHL2) CALL DCOPY(NCF,CMODE(NCP1),1,QCM(NCP1),1)
      DIST = ZERO
      IF(GSCHL2) RETURN
C
      CALL IRCDST(QA,QB,RTRMS,NCOORD,DIST)
      IF(MASWRK) WRITE(IW,9080) DIST
      SAB = DIST
      STOTAL = STOTAL + DIST
      RETURN
C
 9000 FORMAT(//1X,'JUMPING OFF SADDLE POINT ALONG THE IMAGINARY',
     *          ' NORMAL MODE...')
 9010 FORMAT(1X,'**** ERROR, YOU MUST INPUT EITHER 1) A $HESS GROUP'/
     *       1X,'OR 2) BOTH FREQ AND CMODE IN $IRC.')
 9020 FORMAT(1X,'**** ERROR, THE INPUT HESSIAN DOES NOT',
     *          ' HAVE AN IMAGINARY MODE, E(1)=',F16.6)
 9040 FORMAT(1X,'INITIAL DISPLACEMENT=',F6.2,' ALONG MODE WITH ',
     *           'FREQUENCY=',F12.3,' I CM**-1'/
     *       1X,'THE COMPONENTS OF THE DISPLACEMENT ARE (BOHR) -')
 9060 FORMAT(1X,3F11.6,4X,3F11.6,4X,3F11.6)
 9070 FORMAT(1X,'THE COMPONENTS OF THE MASS_WEIGHTED DISPLACEMENT',
     *       '  FOR EFFECTIVE FRAGMENT(S) ARE (BORE*AMU**0.5) -',/,
     *       1X,3F11.6,4X,3F11.6,4X,3F11.6)
 9080 FORMAT(1X,'THE DISTANCE FROM THE SADDLE POINT TO',
     *          ' THE FIRST IRC POINT IS',F10.5,' SQRT(AMU)-BOHR)')
 9100 FORMAT(1X,'*** INPUT NORMAL MODE NORMALIZED BY FACTOR OF',F10.4)
      END
C*MODULE RXNCRD  *DECK INIIRM
      SUBROUTINE INIIRM(EVIB,FORWRD,NCOORD,NCF,FREQ,CMODE,NHOPT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL FORWRD
C
      DIMENSION CMODE(*)
C
      COMMON /FMCOM / XX(1)
C
C     ----- GROW MEMORY FOR INIIRC -----
C
      NC1 = NCOORD + NCF
      NC2 = (NC1*NC1 + NC1)/2
      NC3 = NC1 * NC1
      LOADFM = 0
      CALL VALFM(LOADFM)
      I10 =   1 + LOADFM
      I20 = I10 + NC3
      I30 = I20 + NC2
      I40 = I30 + NC1
      I50 = I40 + NC1 * 8
      I60 = I50 + NC3
      LAST= I60 + NC1
      NEED = LAST - LOADFM
      CALL GETFM(NEED)
      CALL INIIRC(EVIB,FORWRD,NCOORD,NCF,NC2,FREQ,CMODE,NHOPT,
     *            XX(I10),XX(I20),XX(I30),XX(I40),XX(I50),XX(I60))
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE RXNCRD  *DECK IRCDST
      SUBROUTINE IRCDST(A,B,RTRMS,NCOORD,DIST)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(NCOORD),B(NCOORD),RTRMS(NCOORD)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- COMPUTE MASS-WEIGHTED DISTANCE FROM POINT A TO B -----
C           RTRMS CONTAINS INVERSE SQUARE ROOT OF MASSES
C
      DIST = ZERO
      DO 100 I=1,NCOORD
         DIFF = (A(I) - B(I))/RTRMS(I)
         DIST = DIST + DIFF * DIFF
  100 CONTINUE
      DIST = SQRT(DIST)
      RETURN
      END
C*MODULE RXNCRD  *DECK IRCGS2
      SUBROUTINE IRCGS2(FC,A,VEC,EIG,SCR,ISCR,BUF1,NC,NCF,FIRST,LAST,
     *                  SADDLE,CVGED)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION K0,K1,K2,K3,LAMDA
C
      LOGICAL FIRST,LAST,SADDLE,MASWRK,GOPARR,DSKWRK,
     *        STBSKP,GSCHL2,GOTEH,CVGED,EXTPRT,EXTPUN,ALLPUN,
     *        FRODO1,FRODO2,FGONLY,MMONLY,QMMM
C
      DIMENSION FC(*),A(NC+NCF,NC+NCF),VEC(NC+NCF,NC+NCF),EIG(NC+NCF),
     *          SCR(NC+NCF,8),ISCR(NC+NCF),BUF1(NC+NCF,NC+NCF)
      DIMENSION CMASS(3)
C
C        NOTE THAT -FC- AND -A- ARE EQUIVALENT BY CALL
C
      PARAMETER (MXATM=2000, MXSH=5000, MXAO=8192,
     *           MXPT=2000, MXFRG=1050, MXFGPT=12000, MXIFRQ=12)
      PARAMETER (MXNCF=6*MXFRG, MXA3=3*MXATM,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      CHARACTER*8 FRGNME
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3*MXFRG),TORQ(3*MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FMCOM / XX(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      COMMON /FUNCT / E,EG(MXA3)
      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 /IRCDAT/ EA,EB,EC,SAB,SAC,STOTAL,GANORM,GBNORM,GCNORM,
     *                GBIS,GCUT,RCUT,ACUT,STRIDE,OPTTOL,
     *                QA(MXA3+MXNCF),QB(MXA3+MXNCF),QC(MXA3+MXNCF),
     *                GA(MXA3+MXNCF),GB(MXA3+MXNCF),GC(MXA3+MXNCF),
     *                QCM(MXA3+MXNCF),QBM(MXA3+MXNCF),QAM(MXA3+MXNCF),
     *                GA2(MXA3+MXNCF),GA1(MXA3+MXNCF),GA0(MXA3+MXNCF),
     *                QCNM(MXA3+MXNCF),QCN(MXA3+MXNCF),TGP(MXA3+MXNCF),
     *                BIS(MXA3+MXNCF),EBIS(10),XBIS(10),
     *                RTRMS(MXA3+MXFRG*21),K0(MXA3+MXNCF),
     *                K1(MXA3+MXNCF),K2(MXA3+MXNCF),K3(MXA3+MXNCF),
     *                MXOPT,IHUPD,IPIRC,NPIRC,NPRT,NPUN,STBSKP,GSCHL2
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      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 /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TINOPT/ MMONLY,QMMM
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /MASSES/ ZMASS(MXATM)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00,
     *           THREE=3.0D+00, HALF=0.5D+00,
     *           DTOL=1.0D-08, RTOL=1.0D-04,
     *           BIG=1.0D+03, DEVI=5.0D-02)
C
      PARAMETER (TMNIN=1.0D-09, EPS=1.0D-16, TEN=1.0D+01)
C
      DATA CHECK/8HCHECK   /
      DATA GUGA   /8HGUGA    /
C
C     THIS REPRESENTS AN IMPLEMENTATION OF GONZALEZ-SCHLEGEL'S 2ND-ORDER
C     ALGORITHM (J.CHEM.PHYS. 90,2154 + 95,5853. J.PHYS.CHEM. 94,5523).
C
      EXTPRT = NPRT.EQ.1  .OR.  (NPRT.EQ. 0  .AND. (FIRST.OR.LAST))
      EXTPUN = NPUN.GE.0  .OR.  (NPUN.EQ.-1  .AND. (FIRST.OR.LAST))
      ALLPUN = FIRST  .OR.  LAST  .OR.  (NPUN.EQ.1)
C
C     QC IS CURRENT ESTIMATE OF THE CONSTRAINED OPTIMIZED GEOMETRY.
C     INITIAL VALUE IS 1/2S ALONG THE GRADIENT GA FROM QB(PIVOT POINT).
C
      FGONLY = NAT.EQ.0  .AND.  NFRG.GT.0
      NCP1 = NC + 1
      NC1 = NC + NCF
      NC3 = NC1*NC1
      TODEG = 180.0D+00/ACOS(-ONE)
      HAFS  = HALF*STRIDE
C
C     ----- IF NOT STARTING FROM SADDLE POINT, PREVIOUS IRC GEOMETRY,
C           GRADIENT AND HESSIAN ARE READ IN OR MUST BE GENERATED -----
C     WHEN STARTING FROM THE SADDLE, THE INITIAL CIRCLE CENTER WAS
C     ALREADY GENERATED BY FOLLOWING IMAGINARY MODE BY HALF OF STRIDE.
C
      IF(FIRST .AND. .NOT.SADDLE) THEN
         GRAD = DDOT(NC1,GA,1,GA,1)
         IF(GRAD.EQ.ZERO) THEN
            IF (NPIRC.EQ.0) NPIRC=1
            IF(MASWRK) THEN
               WRITE(IW,9200)
               WRITE(IW,9000) NPIRC-1,STOTAL,ZERO
               CALL PRATM(IW,2)
               CALL INTR
            END IF
            NPRINT=-5
            IF(NPRT.GT.0)                         NPRINT= 7
            IF(NPRT.EQ.0  .AND.  (FIRST.OR.LAST)) NPRINT= 7
            IF(NPRT.LT.0)                         NPRINT=-5
            NPUNCH=0
            IF(NPUN.GT.0)                         NPUNCH= 2
            IF(NPUN.EQ.0  .AND.  (FIRST.OR.LAST)) NPUNCH= 2
            IF(NPUN.LT.0)                         NPUNCH= 0
            E = ZERO
            CALL ENERGX
            IF(E .EQ. ZERO  .AND.  EXETYP .NE. CHECK) THEN
               IF(MASWRK) WRITE(IW,9080)
               CALL ABRT
            END IF
            EA = E
            CALL VCLR(EG,1,NC)
            IF (CITYP.EQ.GUGA) THEN
               CALL CIGRAD
            ELSE
               CALL HFGRAD
            END IF
            IF (QMMM) CALL TOYS
            CALL EGMOVE(GA,EG,DEFT,TORQ,NC,NFRG)
         END IF
C
         IF(EXETYP.EQ.CHECK) RETURN
C                            ******
C
C              GENERATE INITIAL CIRCLE CENTER
C
         CALL DCOPY(NC,C,1,QA,1)
         DO 10 I = 1, NC
            GA0(I) = GA(I) * RTRMS(I)
   10    CONTINUE
         CALL EFGRDM(NC,NCF,NFRG,GA,RTRMS,GA0,0)
         DO 15 I = 1, NC
            QAM(I) = QA(I) / RTRMS(I)
   15    CONTINUE
         GALEN = DDOT(NC1,GA0,1,GA0,1)
         GALEN = SQRT(GALEN)
         FACT  = HAFS/GALEN
         CALL DCOPY(NC1,GA0,1,K0,1)
         CALL DSCAL(NC1,-FACT,K0,1)
         DO 20 I = 1,NC
            K1(I) = K0(I)*RTRMS(I)
   20    CONTINUE
         CALL VADD(K0,1,QAM,1,QBM,1,NC)
         IF(NCF.GT.0) CALL DCOPY(NCF,K0(NCP1),1,QBM(NCP1),1)
         CALL VADD(K1,1,QA,1,QB,1,NC)
         CALL EFGEOM(NC,NCF,QBM,QB,RTRMS)
         CALL VADD(K0,1,QBM,1,QCM,1,NC)
         IF(NCF.GT.0) CALL DCOPY(NCF,K0(NCP1),1,QCM(NCP1),1)
         CALL VADD(K1,1,QB,1,QC,1,NC)
         GOTEH = .FALSE.
         CALL FCMIN(A,NC1,GOTEH)
         IF(GOTEH) THEN
C
            CALL DAWRIT(IDAF,IODA,A,NC3,4,0)
         ELSE
            IF(MASWRK) WRITE(IW,9210)
            CALL ABRT
         END IF
      GO TO 40
      END IF
C
      NPRINT=-5
      NPUNCH= 0
C
      DO 30 I = 1,NC
         RTRMSI = RTRMS(I)
         QAM(I) = QA(I) / RTRMSI
         QBM(I) = QB(I) / RTRMSI
         QCM(I) = QC(I) / RTRMSI
   30 CONTINUE
      CALL EFGEOM(NC,NCF,QBM,QB,RTRMS)
C
   40 CONTINUE
      FIRST = .FALSE.
      CALL DCOPY(NC1,QCM,1,BIS,1)
C
      CALL VALFM(LOADFM)
      LCSV   = LOADFM + 1
      LDSV   = LCSV   + 3*MXFGPT
      LQSV   = LDSV   + 3*MXFGPT
      LOSV   = LQSV   + 6*MXFGPT
      LPSV   = LOSV   +10*MXFGPT
      LSVP   = LPSV   + 3*MXFGPT
      LRSV   = LSVP   + 9*MXFGPT
      LDPSV  = LSVP   + 9*MXFGPT
      LSVDP  = LDPSV  + 3*MXFGPT*MXIFRQ
      LDPSV  = LSVDP  + 9*MXFGPT*MXIFRQ
      LSVC   = LDPSV  + 3*MXFGPT
      LSVM   = LSVC   + 3*MXFRG
      LASTSV = LSVM   + 6*MXFRG
      NEEDSV = LASTSV - LOADFM -1
      IF(NFRG.GT.0) CALL GETFM(NEEDSV)
C
      CALL EFSAVE(0,XX(LCSV),XX(LDSV),XX(LQSV),XX(LOSV),
     *            XX(LPSV),XX(LSVP),XX(LDPSV),XX(LSVDP),
     *            XX(LRSV),XX(LSVC),XX(LSVM))
      CALL EFGEOM(NC,NCF,QCM,QC,RTRMS)
C
C     ----- START OF CONSTRAINED OPTIMIZATION LOOP -----
C     QA:   PREVIOUS POINT ON THE PATH, NO DEFINITION FOR EF
C                                       USED AS GA2 FOR EF
C     QAM:  PREVIOUS POINT ON THE PATH (M-W); NO DEFINITION FOR EF,
C                                       USED AS STORGE OF QCM FOR EF
C     QB:   PIVOT COORDINATES; FOR EF AS QB-QA
C     QBM:  PIVOT COORDINATES (M-W); FOR EF AS QBM-QAM
C     QC:   CURRENT COORDINATES; FOR EF AS QC-QB
C     QCM:  CURRENT COORDINATES (M-W); FOR EF AS QCM-QBM
C     QCN:  PREDICT COORDINATES; FOR EF AS QCN-QB
C     QCNM: PREDICT COORDINATES (M-W); FOR EF AS QCNM-QBM
C     GA1:  PREVIOUS COORDINATE; FOR EF AS QCM'-QCM
C     GC:  CURRENT GRADIENT
C     GB:  PREVIOUS GRADIENT
C     GA2: CURRENT GRADIENT (M-W)
C     GA0: PREVIOUS TANGENT GRADIENT (M-W)
C     TGP: CURRENT TANGENT GRADIENT (M-W)
C     K0:  QCM-QBM
C     K1:  SCRETCH
C     K2:  SCRETCH
C     K3:  SCRETCH
C
      IF(MASWRK) WRITE(IW,9220)
      DO 500 II = 1,MXOPT
         NCYC=II
         IF(MASWRK) WRITE(IW,8010) II,NPIRC
         CALL DCOPY(NC,QC,1,C,1)
C
         CALL PRATM(IW,2)
C
         E = ZERO
         CALL ENERGX
         IF(E.EQ.ZERO  .AND.  EXETYP.NE.CHECK) THEN
            IF(MASWRK) WRITE(IW,9080)
            CALL ABRT
         END IF
         EC = E
         IF(II.EQ.1) ETMP = EC
C
         CALL VCLR(EG,1,NC)
         IF (CITYP.EQ.GUGA) THEN
            CALL CIGRAD
         ELSE
            IF (FGONLY) CALL EFGRAD
            IF (.NOT. FGONLY) CALL HFGRAD
         END IF
         IF (QMMM) THEN
            CALL TOYS
            EC = E
            IF(II.EQ.1) ETMP = EC
         END IF
C
         IF(EXETYP.EQ.CHECK) THEN
             IF(NFRG.GT.0) CALL RETFM(NEEDSV)
             RETURN
         END IF
C
         CALL EGMOVE(GC,EG,DEFT,TORQ,NC,NFRG)
         IMAX = IDAMAX(NC1,GC,1)
         GMAX = ABS(GC(IMAX))
         GRMS = DDOT(NC1,GC,1,GC,1)
         GRMS = SQRT(GRMS/NC1)
         IF(GMAX.LT.OPTTOL  .AND.
     *      GRMS.LT.OPTTOL/THREE  .AND.
     *      STOTAL.GT.ONE) THEN
            IF(MASWRK) WRITE(IW,9350) GMAX,GRMS
            CVGED = .TRUE.
         ELSE
            CVGED = .FALSE.
         END IF
C
         DO 50 I=1,NC
            GA2(I) = GC(I) * RTRMS(I)
   50    CONTINUE
         CALL EFGRDM(NC,NCF,NFRG,GC,RTRMS,GA2,0)
         GCNORM = DDOT(NC1,GA2,1,GA2,1)
         GCNORM = SQRT(GCNORM)
C
C     ----- CONSTRUCT A POINT P' ON CIRCLE AND TANGENT GRADIENT, AND
C     CHECK RNORM. IF RNORM DIFFER SIGNIFICANTLY FROM 1/2S YOU ARE
C     HAVING NUMERICAL PROBLEMS. THIS HAPPENS WHEN LAMDA IS CLOSE TO
C     EONE SINCE STEP IS SENSITIVE TO ROUND OFF ERRORS. THIS USUALLY
C     INDICATE THAT THE NEXT POINT IS VERY CLOSE TO A MINIMUM. -----
C
         CALL VSUB(QBM,1,QCM,1,K0,1,NC)
         IF(NCF.GT.0) CALL DCOPY(NCF,QCM(NCP1),1,K0(NCP1),1)
         CALL DCOPY(NC1,K0,1,K1,1)
         RNORM = DDOT(NC1,K1,1,K1,1)
         RNORM = SQRT(RNORM)
         IF (FGONLY) TMPTOL = RTOL
         IF (.NOT. FGONLY) TMPTOL = NAT*RTOL
         IF (ABS(RNORM-HAFS) .GT. TMPTOL) THEN
            IF(MASWRK) WRITE(IW,9230) RNORM,HAFS
            CALL ABRT
         END IF
         RNM1  = ONE/RNORM
         CALL DSCAL(NC1,RNM1,K1,1)
         GALONG = DDOT(NC1,K1,1,GA2,1)
         DO 60 I = 1,NC1
            TGP(I) = GA2(I) - GALONG*K1(I)
   60    CONTINUE
         GTNORM = DDOT(NC1,TGP,1,TGP,1)
         GTNORM = SQRT(GTNORM)
C
C     ----- AFTER THE FIRST OPTIMIZATION POINT UPDATE HESSIAN BY EITHER
C           POWELL OR BFGS (DEFAULT) -----
C
         IF(II.GT.1) THEN
            CALL VSUB(GB,1,GC,1,K3,1,NC1)
            CALL VSUB(GA1,1,QC,1,K1,1,NC)
            IF(NCF.GT.0) CALL DCOPY(NCF,GA1(NCP1),1,K1(NCP1),1)
            CALL DAREAD(IDAF,IODA,A,NC3,4,0)
            CALL MRARBR(A,NC1,NC1,NC1,K1,NC1,1,K2,NC1)
C
C     ----- POWELL UPDATE -----
C
            IF(IHUPD.EQ.1) THEN
               CALL VSUB(K2,1,K3,1,K2,1,NC1)
               DDS = DDOT(NC1,K1,1,K1,1)
               DDTD = DDOT(NC1,K2,1,K1,1)
               DDTD = DDTD/DDS
               DO 110 I = 1,NC1
                  DO 100 J = 1,I
                     TMP = K2(I)*K1(J) + K1(I)*K2(J) - K1(I)*DDTD*K1(J)
                     A(I,J) = A(I,J) + TMP/DDS
                     A(J,I) = A(I,J)
  100             CONTINUE
  110          CONTINUE
            END IF
C
C     ----- BFGS UPDATE -----
C
            IF(IHUPD.EQ.2) THEN
               DDS = DDOT(NC1,K3,1,K1,1)
               IF (DDS .LT. ZERO .AND. MASWRK) WRITE(IW,9240)
               DDTD = DDOT(NC1,K2,1,K1,1)
               DO 130 I = 1,NC1
                  DO 120 J = 1,I
                     TMP = (K3(I)*K3(J))/DDS - (K2(I)*K2(J))/DDTD
                     A(I,J) = A(I,J) + TMP
                     A(J,I) = A(I,J)
  120             CONTINUE
  130          CONTINUE
            END IF
            CALL DAWRIT(IDAF,IODA,A,NC3,4,0)
            IF (CVGED) GO TO 600
C
C     ----- LINEAR INTERPOLATION ALONG PREVIOUS GEOMETRY STEP -----
C
            DO 140 I = 1,NC
               K1(I) = GA1(I)/RTRMS(I)
  140       CONTINUE
            CALL VSUB(QBM,1,K1,1,K1,1,NC)
            IF(NCF.GT.0) CALL DCOPY(NCF,QAM(NCP1),1,K1(NCP1),1)
            RPNORM = DDOT(NC1,K1,1,K1,1)
            RPNORM = SQRT(RPNORM)
            PPDOT  = DDOT(NC1,K0,1,K1,1)
            PANG   = PPDOT/(RNORM*RPNORM)
            IF (ABS(PANG) .GT. ONE) PANG = PANG/ABS(PANG)
            PANG   = ACOS(PANG)
C NOTE
            CALL VSUB(K1,1,K0,1,K2,1,NC1)
C           G2PRIM = DDOT(NC1,GA0,1,GA0,1)
C           G2PRIM = SQRT(G2PRIM)
            G2PRIM = DDOT(NC1,K2,1,GA0,1)
C           G1PRIM = GTNORM
            G1PRIM = DDOT(NC1,K2,1,TGP,1)
            THETA  = PANG*G2PRIM/(G2PRIM-G1PRIM)
C
            IF(MASWRK) WRITE(IW,9250) PANG*TODEG,THETA*TODEG
            TCOM = PANG-THETA
            IF((ABS(THETA).GT.ACUT) .AND. (ABS(TCOM).GT.ACUT)) THEN
               IF(MASWRK) WRITE(IW,9260) ACUT*TODEG
               GO TO 165
            END IF
C
            IF((THETA.LT.ZERO) .AND. (PANG.LT.ACUT)) THEN
              IF(MASWRK) WRITE(IW,9265)
              GO TO 165
            END IF
C           IF((THETA.GT.PANG) .AND. (((THETA-PANG)/PANG).GT.0.1)) THEN
C             IF(MASWRK) WRITE(IW,9266)
C             THETA  = 1.1*PANG
C             TCOM = PANG-THETA
C             IF(MASWRK) WRITE(IW,9250) PANG*TODEG,THETA*TODEG
C           END IF
C
            TMP1 = COS(THETA) - SIN(THETA)*COS(PANG)/SIN(PANG)
            TMP2 = SIN(THETA)/SIN(PANG)
            TMP3 = THETA/PANG
            DO 150 I = 1,NC
               RTRMSI = RTRMS(I)
               K2(I) = K1(I)*TMP1 + K0(I)*TMP2
               GA2(I) = GB(I)*RTRMSI*(1-TMP3) + GA2(I)*TMP3
               GC(I) = GA2(I)/RTRMSI
C              GC(I) = GB(I)*(1-TMP3) + GC(I)*TMP3
C              GA2(I) = GC(I)*RTRMSI
  150       CONTINUE
            DO 155 I = NCP1,NC1
               K2(I) = K1(I)*TMP1 + K0(I)*TMP2
               GC(I) = GB(I)*(1-TMP3) + GC(I)*TMP3
               GA2(I) = QA(I)*(1-TMP3) + GA2(I)*TMP3
  155       CONTINUE
            CALL VADD(K2,1,QBM,1,QCM,1,NC)
            IF(NCF.GT.0) CALL DCOPY(NCF,K2(NCP1),1,QCM(NCP1),1)
            DO 160 I = 1,NC
               QC(I) = QCM(I)*RTRMS(I)
  160       CONTINUE
            CALL EFSAVE(1,XX(LCSV),XX(LDSV),XX(LQSV),XX(LOSV),
     *                  XX(LPSV),XX(LSVP),XX(LDPSV),XX(LSVDP),
     *                  XX(LRSV),XX(LSVC),XX(LSVM))
            CALL EFGEOM(NC,NCF,QCM,QC,RTRMS)
            CALL DCOPY(NC1,K2,1,K0,1)
         END IF
C
C     ----- READ IN HESSIAN, DIAGONALIZE AND PROJECT K0 AND GA2 ALONG
C           HESSIAN EIGENVECTORS -----
C
  165    CONTINUE
         CALL DAREAD(IDAF,IODA,VEC,NC3,4,0)
         CALL HESMAS(NC,NCF,NFRG,VEC,RTRMS,FC,0)
C
C        PROJECT HESSIAN
C
         CALL CENMAS(NAT,NFRG,C,EIG,ZMASST,CMASS,ZMASS)
         CALL CPYTSQ(FC,VEC,NC1,1)
         CALL PRJFC(.FALSE.,.TRUE.,ZMASST,VEC,SCR,EIG,RTRMS,FRGMI,
     *              BUF1,FC,NAT,NC,NC1)
         CALL CPYSQT(VEC,FC,NC1,1)
C
         IERR = 0
         CALL GLDIAG(NC1,NC1,NC1,FC,SCR,EIG,VEC,IERR,ISCR)
         IF (IERR .NE. 0) CALL ABRT
C
         DO 190 I = 1,NC1
            K1(I) = DDOT(NC1,K0,1,VEC(1,I),1)
            K2(I) = DDOT(NC1,GA2,1,VEC(1,I),1)
            IF (ABS(EIG(I)).LT.TMNIN) THEN
               EIG(I)=ZERO
               K1(I)=ZERO
               K2(I)=ZERO
            END IF
  190    CONTINUE
C
C     ----- ITERATION FOR LAMBDA -----
C
         DO 554 I=1,NC1
            IF (EIG(I).NE.ZERO) EONE=EIG(I)
            IF (EIG(I).NE.ZERO) GO TO 555
  554    CONTINUE
  555    CONTINUE
         SSMIN=MAX(ABS(EONE)*EPS,(TEN*EPS))
         SSMAX=MAX(BIG,ABS(EONE))
         SSMAX=SSMAX*BIG
         BULI=EONE
         FRODO1=.FALSE.
         FRODO2=.FALSE.
         LAMDA = EONE - DEVI
         BL = LAMDA - DEVI
         BU = LAMDA + DEVI*HALF
         HAFS2 = HAFS*HAFS
  195    FL = ZERO
         FU = ZERO
         DO 200 I = 1,NC1
            TMP1 = EIG(I)*K1(I) - K2(I)
            TMP2 = EIG(I) - BL
            FL   = FL + (TMP1*TMP1)/(TMP2*TMP2)
            TMP2 = EIG(I) - BU
            FU   = FU + (TMP1*TMP1)/(TMP2*TMP2)
  200    CONTINUE
         FL = FL - HAFS2
         FU = FU - HAFS2
         IF (FL*FU .LT. ZERO) GO TO 205
         BL = BL - (EONE-BL)
         BU = BU + HALF*(EONE-BU)
         IF (BL.LE.-SSMAX) THEN
            BL = -SSMAX
            FRODO1=.TRUE.
         END IF
         IF (ABS(BULI-BU).LE.SSMIN) THEN
            BU = BULI-SSMIN
            FRODO2=.TRUE.
         END IF
         IF (FRODO1.AND.FRODO2) THEN
            IF(MASWRK) WRITE(IW,9270) EONE,BL,BU,FL,FU
            CALL ABRT
         END IF
         GO TO 195
C
  205    CONTINUE
         NCNT = 0
         XLAMDA = ZERO
  215    CONTINUE
         FL = ZERO
         FU = ZERO
         FM = ZERO
         LAMDA = HALF*(BU+BL)
         DO 220 I = 1,NC1
            EIGN = EIG(I)
            TMP1 = EIGN*K1(I) - K2(I)
            TMP2 = EIGN - BL
            FL   = FL + (TMP1*TMP1)/(TMP2*TMP2)
            TMP2 = EIGN - BU
            FU   = FU + (TMP1*TMP1)/(TMP2*TMP2)
            TMP2 = EIGN - LAMDA
            FM   = FM + (TMP1*TMP1)/(TMP2*TMP2)
  220    CONTINUE
         FL = FL - HAFS2
         FU = FU - HAFS2
         FM = FM - HAFS2
         IF (ABS(XLAMDA-LAMDA) .LT. DTOL) GO TO 225
         NCNT = NCNT + 1
         IF(NCNT.GT.1000) THEN
            IF(MASWRK) WRITE(IW,9280) BL,BU,LAMDA,FM
            CALL ABRT
         END IF
         XLAMDA = LAMDA
         IF (FM*FU .LT. ZERO) BL = LAMDA
         IF (FM*FL .LT. ZERO) BU = LAMDA
         GO TO 215
C
  225    CONTINUE
         IF (LAMDA .GT. EONE .AND. MASWRK) WRITE(IW,9290) LAMDA,EONE
         IF (ABS(EONE-LAMDA) .LT. DTOL) THEN
            IF(MASWRK) WRITE(IW,9300) LAMDA,EONE
            CALL ABRT
         END IF
C
C     ----- FORM NR STEP -----
C
         CALL VCLR(K3,1,NC1)
         DO 240 I = 1,NC1
            IF (EIG(I).EQ.ZERO) THEN
               TMP=ZERO
            ELSE
               TMP = (LAMDA*K1(I)-K2(I))/(EIG(I)-LAMDA)
            END IF
            DO 230 J = 1,NC1
               K3(J) = K3(J) + TMP*VEC(J,I)
  230       CONTINUE
  240    CONTINUE
         IF (NFRG.EQ.0) CALL SYMDR(K3)
C
         CALL VADD(K3,1,QCM,1,QCNM,1,NC1)
         DO 250 I = 1,NC
            K3(I) = K3(I)*RTRMS(I)
  250    CONTINUE
C   GA1(NCP1) IS CHANGE OF P' -> P
         CALL EFGRDM(NC,NCF,NFRG,K3,RTRMS,GA1,0)
C
         IF(NCF.GT.0) CALL DCOPY(NCF,GA1(NCP1),1,K3(NCP1),1)
         CALL VADD(K3,1,QC,1,QCN,1,NC)
         DISP = DDOT(NC1,K3,1,K3,1)
         DISP = SQRT(DISP)
C
C     ----- CONVERGENCY CHECK ON TANGENT VECTOR AND DISPLACEMENT -----
C
         CALL VSUB(QBM,1,QCNM,1,K1,1,NC)
         IF(NCF.GT.0) CALL DCOPY(NCF,QCNM(NCP1),1,K1(NCP1),1)
C
         X1 = DDOT(NC1,K0,1,K0,1)
         X2 = DDOT(NC1,K1,1,K1,1)
         X3 = DDOT(NC1,K0,1,K1,1)
         X4 = X3/SQRT(X1*X2)
         AX = ABS(X4)
         IF (AX .GT. ONE) X4 = X4/AX
         X4 = TODEG*ACOS(X4)
         IF(MASWRK) WRITE(IW,9310) EC,GTNORM,X4,DISP
         IF(GTNORM.LT.GCUT  .AND.  DISP.LT.RCUT) GO TO 600
C
C             ARE WE IN TROUBLE WITH ENERGY GOING UP?
C
         IF(EC.GT.EB   .AND.  MASWRK) WRITE(IW,9330)
C
C     ----- CIRCLE OPTIMIZATION NOT CONVERGED -----
C     MOVE TO NEXT CIRCLE POINT
C
         EB = EC
C
         CALL EFSAVE(1,XX(LCSV),XX(LDSV),XX(LQSV),XX(LOSV),
     *               XX(LPSV),XX(LSVP),XX(LDPSV),XX(LSVDP),
     *               XX(LRSV),XX(LSVC),XX(LSVM))
         CALL EFGEOM(NC,NCF,QCNM,QCN,RTRMS)
C
C QAM FOR M-W QCM (FRAG)
C QA FOR M-W GA2 (FRAG)
C
         IF(NCF.GT.0) CALL DCOPY(NCF,QCM(NCP1),1,QAM(NCP1),1)
         IF(NCF.GT.0) CALL DCOPY(NCF,GA2(NCP1),1,QA(NCP1),1)
         CALL DCOPY(NC1,QCNM,1,QCM,1)
         CALL DCOPY(NC1,GC,1,GB,1)
         CALL DCOPY(NC1,TGP,1,GA0,1)
         CALL DCOPY(NC,QC,1,GA1,1)
         CALL DCOPY(NC1,QCN,1,QC,1)
  500 CONTINUE
C
C     ----- END OF CONSTRAINED OPTIMIZATION LOOP -----
C     TOO MANY ITERATIONS IS AN ABEND
C
      IF(MASWRK) WRITE(IW,9340) MXOPT
      CALL ABRT
      STOP
C
  600 CONTINUE
      IF(NFRG.GT.0) CALL RETFM(NEEDSV)
C
      CALL VSUB(QBM,1,QCM,1,K0,1,NC)
      IF(NCF.GT.0) CALL DCOPY(NCF,QCM(NCP1),1,K0(NCP1),1)
      CALL VSUB(QBM,1,BIS,1,K1,1,NC)
      IF(NCF.GT.0) CALL DCOPY(NCF,BIS(NCP1),1,K1(NCP1),1)
      X1 = DDOT(NC1,K0,1,K0,1)
      X2 = DDOT(NC1,K1,1,K1,1)
      X3 = DDOT(NC1,K0,1,K1,1)
      X4 = X3/SQRT(X1*X2)
      AX = ABS(X4)
      IF (AX .GT. ONE) X4 = X4/AX
      TANG = TODEG*ACOS(X4)
      IF(MASWRK) WRITE(IW,9320) TANG,ETMP-EC
C
      IF(EC.GT.EA) THEN
         IF(MASWRK) WRITE(IW,9090) EC,EA
         CALL ABRT
      END IF
C
C     ----- WRITE OUT GEOMETRY AND GRADIENT -----
C
      CALL IRCDST(QA,QC,RTRMS,NC,SAC)
      CALL EFDST(NFRG,QBM(NCP1),QCM(NCP1),SAC)
      STOTAL = STOTAL + SAC
      IF(MASWRK) THEN
         WRITE(IW,9180)
         WRITE(IW,9160) SAC
         WRITE(IW,9000) NPIRC,STOTAL,EC
         WRITE(IW,9010) NCYC
         WRITE(IW,*) ' '
         CALL PRATM(IW,1)
         WRITE(IW,*) ' '
         IF(NT.GT.1) CALL PRATM(IW,2)
         IF ((NPRT.GT.-2 .OR. FIRST .OR. LAST) .AND. .NOT.FGONLY)
     *      CALL INTR
         IF (.NOT.FGONLY) WRITE(IW,9100)
         I0=0
         DO 670 I=1,NAT
            WRITE(IW,9060) ANAM(I),BNAM(I),(GC(I0+J),J=1,3)
            I0=I0+3
  670    CONTINUE
C
         WRITE(IP,8005) NPIRC,STOTAL
      END IF
C
C        OPTIONAL PRINT/PUNCH OF ORBITALS
C
      IF ((EXTPRT .OR. EXTPUN) .AND. .NOT.FGONLY)
     *   CALL IRCMOS(EXTPRT,EXTPUN,ALLPUN,NPIRC,STOTAL,EC)
C
C     ----- PUNCH DATA FOR POSSIBLE VTST USE -----
C
      FACT = ONE/GCNORM
      CALL DSCAL(NC1,FACT,GA2,1)
      IF(MASWRK) THEN
         WRITE(IPIRC,8000) NPIRC,STOTAL,EC
C            DUMP COORDINATES
         I0=0
         DO 680 I=1,NAT
            ZNUC = ZAN(I) + IZCORE(I)
            WRITE(IPIRC,8020) ANAM(I),BNAM(I),ZNUC,(C(I0+J),J=1,3)
            I0=I0+3
  680    CONTINUE
         IF (NFRG.EQ.0) GO TO 682
         DO 681 I=1,NMTTPT
            IF (FGNUC(I).EQ.ZERO) GO TO 681
            WRITE(IPIRC,8021) FRGNME(I),FGNUC(I),(EFC(J,I),J=1,3)
  681    CONTINUE
  682    CONTINUE
         IF (QMMM) CALL PRTXYZ(IPIRC)
C            DUMP GRADIENTS
         WRITE(IPIRC,8040) GCNORM
         I0=0
         DO 690 I=1,NAT
            ZNUC = ZAN(I) + IZCORE(I)
            WRITE(IPIRC,8020) ANAM(I),BNAM(I),ZNUC,(GA2(I0+J),J=1,3)
            I0=I0+3
  690    CONTINUE
         CALL FLSHBF(IPIRC)
      END IF
C
C     ----- MOVE DATA AND RETURN TO SEE IF WE NEED MORE POINTS -----
C
      EA = EC
      CALL DCOPY(NC1,GC,1,GA,1)
      CALL DCOPY(NC,QC,1,QA,1)
C
C     ----- CONSTRUCT NEW PIVOT POINT AND NEW CIRCLE GUESS -----
C
      GALEN = DDOT(NC1,GA2,1,GA2,1)
      GALEN = SQRT(GALEN)
      FACT = HAFS/GALEN
      DO 710 I = 1,NC
         K1(I) = -GA2(I)*FACT*RTRMS(I)
  710 CONTINUE
C MASS-WEIGHTED
      DO 720 I = NCP1,NC1
         K1(I) = -GA2(I)*FACT
  720 CONTINUE
      CALL VADD(K1,1,QA,1,QB,1,NC)
      CALL VADD(K1,1,QB,1,QC,1,NC)
      IF(NCF.GT.0) CALL DCOPY(NCF,K1(NCP1),1,QBM(NCP1),1)
      IF(NCF.GT.0) CALL DCOPY(NCF,K1(NCP1),1,QCM(NCP1),1)
      NPIRC = NPIRC+1
      RETURN
C
 8000 FORMAT('***** BEGIN IRC INFORMATION PACKET *****'/
     *    'POINT=',I4,' STOTAL=',F10.5,'  E=',F20.10/
     *       'CARTESIAN COORDINATES (BOHR)')
 8005 FORMAT('----- DATA FOR IRC POINT',I6,' AT STOTAL=',F10.5,' -----')
 8010 FORMAT(' CONSTRAINED OPTIMIZATION POINT',I4,' FOR PATH POINT',I4)
 8020 FORMAT(A8,A2,F5.1,3F20.10)
 8021 FORMAT(A8,2X,F5.1,3F20.10)
 8040 FORMAT('MASS-WEIGHTED GRADIENT - ORIG.NORM=',1P,E13.6,
     *       ' (HARTREE/BOHR-SQRT(AMU))')
C
 9000 FORMAT(/'1',9X,'POINT',I4,' ON THE REACTION PATH'/10X,30(1H*)/
     *       10X,'AT PATH DISTANCE STOTAL =',F10.5,4X,' SQRT(AMU)-BOHR'/
     *       10X,'TOTAL ENERGY            =',F14.7,' HARTREE')
 9010 FORMAT(10X,I3,' CONSTRAINED OPTIMIZATION STEPS WERE PERFORMED')
 9060 FORMAT(10X,A8,A2,4X,3F15.10)
 9080 FORMAT(1X,'**** ERROR, SCF FAILED TO CONVERGE')
 9090 FORMAT(1X,'CURRENT ENERGY=',F20.10,', IS HIGHER THAN THE'/
     *       1X,'   LAST ENERGY=',F20.10/
     *       1X,'YOU ARE NO LONGER ON A STEEPEST DESCENT PATH!'/
     *       1X,'YOU MAY JUST HAVE PASSED THROUGH A MINIMUM')
 9100 FORMAT(/10X,'CURRENT CARTESIAN GRADIENT'/
     *        10X,'ATOM',15X,'DE/DX',10X,'DE/DY',10X,'DE/DZ'/
     *        10X,59(1H-))
 9160 FORMAT(1X,'STEP DISTANCE IS',F10.5,' SQRT(AMU)-BOHR)')
 9180 FORMAT(1X,'***** NEXT POINT ON IRC FOUND *****')
 9200 FORMAT(1X,'WARNING! GRADIENTS FROM PREVIOUS IRC POINT NOT FOUND',
     *     ' IN $IRC GROUP. GENERATING NEW GRADIENTS AT CURRENT POINT')
 9210 FORMAT(1X,'GONZALEZ-SCHLEGEL IRC PACE REQUIRES HESSIAN INPUT'/)
 9220 FORMAT(/1X,45(1H*)/
     *        1X,'***** STARTING CONSTRAINED OPTIMIZATION *****'/
     *        1X,'*** TO FIND THE NEXT REACTION PATH POINT ****'/
     *        1X,45(1H*))
 9230 FORMAT(/1X,'RADIUS IN CIRCLE OPTIMIZATION',F10.7,' DEVIATES',
     *       ' SIGNIFICANTLY FROM CONSTRAIN CONDITION',F10.7/
     *       ' IT IS POSSIBLE THAT THE NEXT IRC POINT IS CLOSE TO',
     *       ' A MINIMUM')
 9240 FORMAT(/1X,'WARNING! POSITIVE DEFINITENESS ENDANGERED IN BFGS',
     *       ' UPDATE OF HESSIAN. THIS MAY CAUSE TROUBLE LATER....')
 9250 FORMAT(/1X,'ANGLE TO PREVIOUS OPTIMIZATION POINT IS',F8.3/
     *        1X,'LINEAR INTERPOLATION GIVES ANGLE       ',F8.3)
 9260 FORMAT(/1X,'INTERPOLATION TOO FAR FROM END POINTS. ACUT =',F8.3/
     *        1X,'SKIPPING INTERPOLATION THIS STEP')
 9265 FORMAT(/1X,'INTERPOLATION BECOMES EXTERPOLATION',
     *        1X,'SKIPPING INTERPOLATION THIS STEP')
 9270 FORMAT(/1X,'!!!!! NUMERICAL PROBLEMS IN BRACKETING LAMDA'/
     *        1X,'LOWEST HESSIAN EIGENVALUE, LOWER AND UPPER SEARCH',
     *           ' BOUNDS =',3G20.15)
 9280 FORMAT(/1X,'!!!!! TOO MANY ITERATIONS IN LAMDA SEARCH',/,4F18.10)
 9290 FORMAT(/1X,'WARNING! LAMDA =',F15.10,' LARGER THAN LOWEST',
     *       ' HESSIAN EIGENVALUE =',F15.10)
 9300 FORMAT(/1X,'!!!!! ERROR! LAMDA =',F15.10,' IS SO CLOSE TO LOWEST',
     *       ' HESSIAN EIGENVALUE =',F15.10,' THAT STEP CANNOT BE',
     *       ' CALCULATED ACCURATELY')
 9310 FORMAT(/1X,'CURRENT OPTIMIZATION POINT :'/
     *        1X,'ENERGY',33X,'=',F14.7,' HARTREE'/
     *        1X,'MASS WEIGHTED TANGENT GRADIENT NORM    =',F14.7,
     *        ' HARTREE/BOHR-SQRT(AMU)'/
     *        1X,'NEXT GEOMETRY STEP HAS PREDICTED ANGLE =',6X,F8.3,
     *        ' DEGREES',/
     *        1X,'AND CARTESIAN RMS DISPLACEMENT         =',F14.7,
     *        ' BOHR'/)
 9320 FORMAT(1X,43(1H*)/
     *       1X,'***** END OF CONSTRAINED OPTIMIZATION *****'/
     *       1X,43(1H*)//
     *       1X,'ANGLE BETWEEN GRADIENTS AT CURRENT AND PREVIOUS IRC',
     *        ' POINT =',F8.3,/
     *        1X,'ENERGY WAS LOWERED BY',F18.10,' HARTREE DURING',
     *        ' OPTIMIZATION'/)
 9330 FORMAT(1X,'......ENERGY IS GOING UP......'/)
 9340 FORMAT(/1X,'!!!!! TOO MANY STEPS TAKEN IN CONSTRAINED',
     *        ' OPTIMIZATION',/,1X,'MAX GIVEN (MXOPT) =',I6)
 9350 FORMAT(/1X,'STOPPING BECAUSE GRADIENT IS BELOW -OPTTOL-'/
     *       1X,'YOUR RUN IS PROBABLY VERY NEAR A PES MINIMUM.'/
     *       1X,'MAX GRADIENT=',F9.6,' RMS GRADIENT=',F9.6//)
      END
C*MODULE RXNCRD  *DECK IRCGSM
      SUBROUTINE IRCGSM(NC,NCF,FIRST,LAST,SADDLE,CVGED)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL FIRST,LAST,SADDLE,CVGED
C
      COMMON /FMCOM / XX(1)
C
C     ----- GROW MEMORY FOR IRCGS2 ----
C       -FC-      AT I10
C       -A-       AT I10
C       -VEC-     AT I20
C       -EIGEN-   AT I30
C       -SCR-     AT I40
C       -IA-      AT I50
C
      NC1=NC+NCF
      NC3 = NC1*NC1
      LOADFM = 0
      CALL VALFM(LOADFM)
      I10 = 1 + LOADFM
      I20 = I10 + NC3
      I30 = I20 + NC3
      I40 = I30 + NC1
      I50 = I40 + 8*NC1
      I60 = I50 + NC1
      LASTM= I60 + NC3
      NEED = LASTM - LOADFM - 1
      CALL GETFM(NEED)
C
C     INFORMATION FOR EFFECTIVE FRAGMENTS WILL BE CARRIED BY
C     COMMON BLOCKS -FRGINF- AND -FGRAD-
C
      CALL IRCGS2(XX(I10),XX(I10),XX(I20),XX(I30),XX(I40),
     *            XX(I50),XX(I60),NC,NCF,FIRST,LAST,SADDLE,CVGED)
C
C     ----- RETURN FAST MEMORY -----
C
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE RXNCRD  *DECK IRCSTP
      SUBROUTINE IRCSTP(STRIDE,NCOORD,PACE,SADDLE,STABLZ,DELTA,ELBOW,
     *                  NAMPC,FIRST,LAST,PENULT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION K0,K1,K2,K3
C
      LOGICAL STABLZ,SADDLE,STBSKP,LINMOL,MASWRK,GOPARR,DSKWRK,
     *        GSCHL2,FIRST,LAST,PENULT,MMONLY,QMMM
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00,
     *           ONE80=180.0D+00, GTOL=1.0D-08, TOL5=5.0D-06)
C
      PARAMETER (MXATM=2000, MXA3=3*MXATM, MXSH=5000, MXAO=8192)
      PARAMETER (MXFRG=1050, MXNCF=6*MXFRG)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /FUNCT / E,EG(MXA3)
      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 /IRCDAT/ EA,EB,EC,SAB,SAC,STOTAL,GANORM,GBNORM,GCNORM,
     *                GBIS,GCUT,RCUT,ACUT,DUMMY1,OPTTOL,
     *                QA(MXA3+MXNCF),QB(MXA3+MXNCF),QC(MXA3+MXNCF),
     *                GA(MXA3+MXNCF),GB(MXA3+MXNCF),GC(MXA3+MXNCF),
     *                QCM(MXA3+MXNCF),QBM(MXA3+MXNCF),QAM(MXA3+MXNCF),
     *                GA2(MXA3+MXNCF),GA1(MXA3+MXNCF),GA0(MXA3+MXNCF),
     *                QCNM(MXA3+MXNCF),QCN(MXA3+MXNCF),TGP(MXA3+MXNCF),
     *                BIS(MXA3+MXNCF),EBIS(10),XBIS(10),
     *                RTRMS(MXA3+MXFRG*21),K0(MXA3+MXNCF),
     *                K1(MXA3+MXNCF),K2(MXA3+MXNCF),K3(MXA3+MXNCF),
     *                MXOPT,IHUPD,IPIRC,NPIRC,NPRT,NPUN,STBSKP,GSCHL2
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      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 /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TINOPT/ MMONLY,QMMM
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINMOL
C
      DATA ALINR/8HLINEAR  /, QUAD/8HQUAD    /, CHECK/8HCHECK   /
      DATA AMPC4/8HAMPC4   /, RK4/8HRK4      /
      DATA GUGA   /8HGUGA    /
C
C     ----- JUMP FROM POINT QB ON THE IRC TO A POINT QC BY
C           A LINEAR (POSSIBLY STABILIZED) OR QUADRATIC PACE -----
C
C     ----- IF A LINEAR PACE WITH STABILIZATION HAS BEEN CHOSEN,
C           AND THE PATH IS NOT CURVED VERY MUCH, WE MAY HAVE
C           SKIPPED THE STABILIZATION, AND ALREADY KNOW EB, GB -----
C
      IMARK = 0
      IF(STBSKP .AND. MASWRK) WRITE(IW,9010) NPIRC,STOTAL
      IF(STBSKP) GO TO 300
C
C     ----- BEGINNING POINT QB ON THE REACTION PATH -----
C
      NPRINT = -5
      IF(NPRT.EQ.1) NPRINT=7
      IF(NPRT.EQ.0  .AND.  (FIRST.OR.LAST)) NPRINT=7
C
      NPUNCH = 0
      IF(NPUN.EQ.1) NPUNCH = 2
      IF(NPUN.EQ.0) NPUNCH = 1
      IF(NPUN.EQ. 0  .AND.  (FIRST.OR.LAST)) NPUNCH=2
      IF(NPUN.EQ.-1  .AND.  (FIRST.OR.LAST)) NPUNCH=2
C
      IF(MASWRK) WRITE(IW,9000) NPIRC,STOTAL
      IF(MASWRK) WRITE(IP,9020) NPIRC
      DO 130 I=1,NCOORD
         C(I) = QB(I)
  130 CONTINUE
      CALL PRATM(IW,1)
      IF(NT.GT.1) CALL PRATM(IW,2)
      IF(NPRT.GT.-2  .OR.  FIRST  .OR.  LAST) CALL INTR
C
C     ----- ENERGY AT POINT QB -----
C
      E = ZERO
      CALL ENERGX
      IF(E .EQ. ZERO  .AND.  EXETYP .NE. CHECK) THEN
         IF(MASWRK) WRITE(IW,9080)
         CALL ABRT
      END IF
      EB = E
C
      IF(EB.GT.EA) THEN
         IF(MASWRK) WRITE(IW,9090) EB,EA
         CALL ABRT
      END IF
C
C     ----- GRADIENT AT POINT QB -----
C
      DO 200 I=1,NCOORD
         EG(I) = ZERO
  200 CONTINUE
      IF (CITYP.EQ.GUGA) THEN
         CALL CIGRAD
      ELSE
         CALL HFGRAD
      END IF
      IF (QMMM) CALL TOYS
      IF (EXETYP .EQ. CHECK) RETURN
      DO 220 I=1,NCOORD
         GB(I) = EG(I)
  220 CONTINUE
C
C     ----- PRINT INTERNAL COORDINATES, WITH GRADIENT -----
C
      IF(NZMAT.GT.0) THEN
         CALL BANDBI
         CALL TRANG(EG,NVAR,NCOORD)
         CALL PZANDG(EG,1)
      END IF
C
C     ----- NORMALIZED MASS WEIGHTED GRADIENT AT POINT QB -----
C
      GBNORM = ZERO
      DO 240 I=1,NCOORD
         G = RTRMS(I) * GB(I)
         IF(ABS(G).LT.GTOL) G=ZERO
         GB(I) = G
         GBNORM = GBNORM + G * G
  240 CONTINUE
      GBNORM = SQRT(GBNORM)
      IF(MASWRK) WRITE(IW,9100) GBNORM
      FACT = ONE/GBNORM
      DO 250 I=1,NCOORD
         GB(I) = FACT * GB(I)
  250 CONTINUE
      IF(MASWRK) THEN
      WRITE(IW,9120)
      I0 = 0
      DO 260 I=1,NAT
         WRITE(IW,9140) ANAM(I),BNAM(I),(GB(I0+J),J=1,3)
         I0 = I0 + 3
  260 CONTINUE
      END IF
C
C     ----- PUNCH DATA FROM POINT QB ON IRC FOR POSSIBLE VTST USE -----
C
  300 CONTINUE
      STBSKP=.FALSE.
      IF(MASWRK) THEN
         WRITE(IPIRC,8000) NPIRC,STOTAL,EB
C            DUMP COORDINATES
         I0=0
         DO 310 I=1,NAT
            ZNUC = ZAN(I) + IZCORE(I)
            WRITE(IPIRC,8020) ANAM(I),BNAM(I),ZNUC,(C(I0+J),J=1,3)
            I0=I0+3
  310    CONTINUE
         IF(QMMM) CALL PRTXYZ(IPIRC)
C            DUMP GRADIENT
         WRITE(IPIRC,8040) GBNORM
         I0=0
         DO 320 I=1,NAT
            ZNUC = ZAN(I) + IZCORE(I)
            WRITE(IPIRC,8020) ANAM(I),BNAM(I),ZNUC,(GB(I0+J),J=1,3)
            I0=I0+3
  320    CONTINUE
         CALL FLSHBF(IPIRC)
      END IF
C
C     ----- STEP TO POINT QC BY THE APPROPRIATE PACE -----
C
      NPUNCH = 0
      NPRINT = -5
      IF(MASWRK) WRITE(IW,9150) PACE
      IF(PACE.EQ.RK4)   GO TO 660
      IF(PACE.EQ.AMPC4)  GO TO 550
      IF(SADDLE .AND. NPIRC.EQ.1) GO TO 400
      IF(PACE.EQ.ALINR) GO TO 400
      IF(PACE.EQ.QUAD)  GO TO 500
      CALL ABRT
C
C     ----- LINEAR GRADIENT FOLLOWING -----
C
  400 CONTINUE
      A = -STRIDE
      DO 420 I=1,NCOORD
         QC(I) = A * GB(I)
  420 CONTINUE
      CALL SYMDR(QC)
      DO 440 I=1,NCOORD
         QC(I) = QB(I) + QC(I)*RTRMS(I)
  440 CONTINUE
      IF(STABLZ) THEN
         CALL STBIMK(DELTA,ELBOW,NCOORD,FIRST,PENULT)
         IF (EXETYP .EQ. CHECK) RETURN
      ELSE
C
C     ----- IF NOT JUMPING FROM SADDLE POINT, CALCULATE ANGLE BETWEEN
C     ----- SUCCESSIVE GRADIENTS.
C
         IF (NPIRC .GT. 1) THEN
            DOTPRD = ZERO
            DO 450 I=1,NCOORD
               DOTPRD = DOTPRD + GA(I)*GB(I)
  450       CONTINUE
            THETA = ACOS(DOTPRD)
            PI = ACOS(-ONE)
            THETA = ONE80*THETA/PI
            THETA = ONE80 - THETA
            IF(MASWRK) WRITE(IW,9200) THETA
         END IF
      END IF
      GO TO 700
C
C     ----- QUADRATIC GRADIENT FOLLOWING -----
C
  500 CONTINUE
      A = -STRIDE
      B = (STRIDE*STRIDE)/(TWO*SAB)
      DO 510 I=1,NCOORD
         QC(I) = (A-B)*GB(I) - B*GA(I)
  510 CONTINUE
      CALL SYMDR(QC)
      DO 520 I=1,NCOORD
         QC(I) = QB(I) + QC(I)*RTRMS(I)
  520 CONTINUE
      GO TO 700
C
C     ----- ADAMS-MOULTON PREDICTOR-CORRECTOR ORDER FOUR
C           THE INFORMATION NEEDED TO START THIS METHOD
C           IS OBTAINED BY DOING 1 STEP USING RK4.
C
  550 CONTINUE
      H = -STRIDE
         GO TO (660,564) NAMPC
C
  564    CONTINUE
         DO 578 I=1,NCOORD
            QC(I) = (H/24.0D+00) *
     *              (55.0D+00*GB(I) - 59.0D+00*GA0(I) +
     *               37.0D+00*GA1(I) - 9.0D+00*GA2(I))
  578    CONTINUE
C
      CALL SYMDR(QC)
      DO 580 I=1,NCOORD
           QCNM(I) = QB(I) + QC(I)*RTRMS(I)
  580 CONTINUE
C
C     ----- ENERGY AT PREDICTED POINT
C
      DO 590 I=1,NCOORD
         C(I) = QCNM(I)
  590 CONTINUE
      E = ZERO
      CALL ENERGX
      IF(E .EQ. ZERO  .AND.  EXETYP .NE. CHECK) THEN
         IF(MASWRK) WRITE(IW,9080)
         CALL ABRT
      END IF
C
C     ----- GRADIENT AT PREDICTED POINT
C
      DO 600 I=1,NCOORD
         EG(I) =  ZERO
  600 CONTINUE
      IF (CITYP.EQ.GUGA) THEN
         CALL CIGRAD
      ELSE
         CALL HFGRAD
      END IF
      IF (QMMM) CALL TOYS
      IF(EXETYP .EQ. CHECK) RETURN
C
C     ----- NORMALIZED MASS WEIGHTED GRADIENT AT PREDICTED POINT -----
C
      GBNORM = ZERO
      DO 610 I=1,NCOORD
         G = RTRMS(I) * EG(I)
         IF(ABS(G).LT.GTOL) G=ZERO
         EG(I) = G
         GBNORM = GBNORM + G * G
  610 CONTINUE
      GBNORM = SQRT(GBNORM)
      FACT = ONE/GBNORM
      DO 620 I=1,NCOORD
         TGP(I) = FACT * EG(I)
  620 CONTINUE
C
C     ----- CORRECTOR
C
      DO 630 I=1,NCOORD
         QC(I) = (H/24.0D+00) *
     *           (9.0D+00*TGP(I)+19.0D+00*GB(I)-5.0D+00*GA0(I)+GA1(I))
  630 CONTINUE
      CALL SYMDR(QC)
      DO 640 I=1,NCOORD
         QCN(I) = QB(I) + QC(I)*RTRMS(I)
  640 CONTINUE
C
C     ----- VARIABLE STEPPER
C
      DIFF = ZERO
      SIGMA = ZERO
      DO 650 I=1,NCOORD
         DIFF = (QCNM(I) - QCN(I))/RTRMS(I)
         SIGMA = SIGMA + DIFF*DIFF
  650 CONTINUE
      SIGMA = SQRT(SIGMA)/14.0D+00
      POWER = 1.0D+00/6.0D+00
      HCORR = ((TOL5/3.0D+00)/SIGMA)**POWER
      IF((SIGMA.LT.TOL5).OR.(IMARK.EQ.1)) GO TO 700
      IMARK = 1
      H = HCORR*H
      GO TO 564
C
C     ----- RUNGE-KUTTA ORDER FOUR
C
  660 CONTINUE
      NRK4 = 1
      H = -STRIDE
  665 CONTINUE
         GO TO (672,674,676,678) NRK4
C
  672    CONTINUE
         DO 682 I=1,NCOORD
            K0(I) = H*GB(I)
            QC(I) = 0.5D+00*K0(I)
  682    CONTINUE
         GO TO 670
C
  674    CONTINUE
         DO 684 I=1,NCOORD
            IF(PACE.EQ.AMPC4) GA2(I) = EG(I)
            K1(I) = H*EG(I)
            QC(I) = 0.5D+00*K1(I)
  684    CONTINUE
         GO TO 670
C
  676    CONTINUE
         DO 686 I=1,NCOORD
            IF(PACE.EQ.AMPC4) GA1(I) = EG(I)
            K2(I) = H*EG(I)
            QC(I) = K2(I)
  686    CONTINUE
         GO TO 670
C
  678    CONTINUE
         DO 688 I=1,NCOORD
            IF(PACE.EQ.AMPC4) GA0(I) = EG(I)
            K3(I) = H*EG(I)
            QC(I) = (1.0D+00/6.0D+00) *
     *              (K0(I)+2.0D+00*K1(I)+2.0D+00*K2(I)+K3(I))
  688    CONTINUE
C
  670 CONTINUE
C
      NRK4 = NRK4 + 1
      CALL SYMDR(QC)
      DO 690 I=1,NCOORD
         QC(I) = QB(I) + QC(I)*RTRMS(I)
  690 CONTINUE
      IF(NRK4.GT.4) GO TO 700
C
C    ----- ENERGY AT INTERMEDIATE R-K POINT -----
C
      DO 750 I=1,NCOORD
         C(I) = QC(I)
  750 CONTINUE
      E = ZERO
      CALL ENERGX
      IF(E .EQ. ZERO  .AND.  EXETYP .NE. CHECK) THEN
         IF(MASWRK) WRITE(IW,9080)
         CALL ABRT
      END IF
C
C     ----- GRADIENT AT INTERMEDIATE R-K POINT
C
      DO 692 I=1,NCOORD
        EG(I) = ZERO
  692 CONTINUE
      IF (CITYP.EQ.GUGA) THEN
         CALL CIGRAD
      ELSE
         CALL HFGRAD
      END IF
      IF (QMMM) CALL TOYS
      IF(EXETYP .EQ. CHECK)  RETURN
C
C     ----- NORMALIZED MASS WEIGHTED GRADIENT AT PREDICTED POINT -----
C
      GBNORM = ZERO
      DO 694 I=1,NCOORD
         G = RTRMS(I) * EG(I)
         IF(ABS(G).LT.GTOL) G=ZERO
         EG(I) = G
         GBNORM = GBNORM + G * G
  694 CONTINUE
      GBNORM = SQRT(GBNORM)
      FACT = ONE/GBNORM
      DO 696 I=1,NCOORD
         IF(NRK4.LE.4) THEN
            EG(I) = FACT * EG(I)
         ELSE
            GB(I) = FACT * EG(I)
         END IF
  696 CONTINUE
      IF(NRK4.LE.4) THEN
         GO TO 665
      ELSE
         NRK4 = 0
      END IF
C
C     ----- QC IS THE PREDICTED NEXT POINT ON THE IRC -----
C
  700 CONTINUE
      IF((PACE.EQ.AMPC4).AND.(NAMPC.GT.1)) THEN
         CALL IRCDST(QB,QCN,RTRMS,NCOORD,SBC)
      ELSE
         CALL IRCDST(QB,QC,RTRMS,NCOORD,SBC)
      END IF
      STOTAL = STOTAL + SBC
      IF(MASWRK) THEN
         WRITE(IW,9160) SBC,STOTAL
         WRITE(IW,9180)
      END IF
C
C     ----- CURRENT POINT QB BECOMES OLD POINT QA,
C           NEXT POINT QC BECOMES CURRENT POINT QB. -----
C
      NPIRC = NPIRC+1
      EA = EB
      EB = EC
      SAB = SBC
      GANORM = GBNORM
      GBNORM = GCNORM
      DO 740 I=1,NCOORD
         IF(PACE.EQ.AMPC4) THEN
           IF(NAMPC.EQ.1) THEN
              QB(I)  = QC(I)
              GB(I)  = GC(I)
           ELSE
              GA2(I) = GA1(I)
              QCM(I) = QBM(I)
              GA1(I) = GA0(I)
              QBM(I) = QAM(I)
              GA0(I) = GB(I)
              QAM(I) = QB(I)
              QB(I) = QCN(I)
           END IF
         ELSE
            QA(I) = QB(I)
            QB(I) = QC(I)
            GA(I) = GB(I)
            GB(I) = GC(I)
         END IF
  740 CONTINUE
      IF(NAMPC.EQ.1) NAMPC=NAMPC+1
      RETURN
C
 8000 FORMAT('***** BEGIN IRC INFORMATION PACKET *****'/
     *    'POINT=',I4,' STOTAL=',F10.5,'  E=',F20.10/
     *       'CARTESIAN COORDINATES (BOHR)')
 8020 FORMAT(A8,A2,F5.1,3F20.10)
 8040 FORMAT('MASS-WEIGHTED GRADIENT - ORIG.NORM=',1P,E13.6,
     *       ' (HARTREE/BOHR-SQRT(AMU))')
C
 9000 FORMAT('1',9X,'POINT',I4,' ON THE REACTION PATH'/10X,30(1H*)/
     *       10X,'AT PATH DISTANCE STOTAL=',F10.5//
     *       10X,'THE COORDINATES OF POINT QB ARE')
 9010 FORMAT(1X,'BECAUSE OF THE ELBOW CRITERIA, THE PREVIOUS GEOMETRY',
     *          ' IS IRC POINT',I5/1X,'AT PATH DISTANCE STOTAL=',F10.5)
 9020 FORMAT('----- DATA FOR IRC POINT',I5,' FOLLOWS -----')
 9080 FORMAT(1X,'**** ERROR, SCF FAILED TO CONVERGE')
 9090 FORMAT(1X,'CURRENT ENERGY=',F20.10,', IS HIGHER THAN THE'/
     *       1X,'   LAST ENERGY=',F20.10/
     *       1X,'YOU ARE NO LONGER ON A STEEPEST DESCENT PATH!'/
     *       1X,'THIS IRC SHOULD BE RUN WITH SMALLER STEPS AND/OR',
     *          ' A BETTER INTEGRATION METHOD.')
 9100 FORMAT(/5X,'NORMALIZED MASS WEIGHTED GRADIENT AT POINT QB'/
     *   5X,'(ORIGINAL GBNORM=',F10.6,' HARTREE/BOHR-SQRT(AMU))')
 9120 FORMAT(/10X,'ATOM',12X,'DE/DX',7X,'DE/DY',7X,'DE/DZ'/
     *        10X,46(1H-))
 9140 FORMAT(10X,A8,A2,3(2X,F10.6))
 9150 FORMAT(//1X,'*** JUMPING TO POINT QC VIA A ',A8,'PACE ***')
 9160 FORMAT(1X,'THE DISTANCE FROM QB TO QC (OR QC'') IS',
     *          F10.5,' SQRT(AMU)-BOHR)'/
     *       1X,' TOTAL PATH DISTANCE TO QC (OR QC'') IS',
     *          F10.5,' SQRT(AMU)-BOHR)')
 9180 FORMAT(1X,'......NEXT POINT ON IRC FOUND......',
     *          'RENAMING QC(QC'')-->QB-->QA.')
 9200 FORMAT(/5X,'THE ANGLE BETWEEN GRADIENTS GA AND GB IS',F8.2,
     *         ' DEGREES')
      END
C
C*MODULE RXNCRD  *DECK IRCX
      SUBROUTINE IRCX
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION K0,K1,K2,K3
C
      LOGICAL FORWRD,SADDLE,STABLZ,TSENGY,STBSKP,LINMOL,MASWRK,GOPARR,
     *        DSKWRK,GSCHL2,FIRST,LAST,PENULT,CVGED
      LOGICAL FGONLY
C
      PARAMETER (NNAM=31)
      PARAMETER (MXATM=2000, MXA3=3*MXATM, MXAO=8192,
     *           MXFRG=1050, MXNCF=6*MXFRG,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
C
      COMMON /FMCOM / X(1)
      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 / E,EG(MXA3)
      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 /IRCDAT/ EA,EB,EC,SAB,SAC,STOTAL,GANORM,GBNORM,GCNORM,
     *                GBIS,GCUT,RCUT,ACUT,STRIDE,OPTTOL,
     *                QA(MXA3+MXNCF),QB(MXA3+MXNCF),QC(MXA3+MXNCF),
     *                GA(MXA3+MXNCF),GB(MXA3+MXNCF),GC(MXA3+MXNCF),
     *                QCM(MXA3+MXNCF),QBM(MXA3+MXNCF),QAM(MXA3+MXNCF),
     *                GA2(MXA3+MXNCF),GA1(MXA3+MXNCF),GA0(MXA3+MXNCF),
     *                QCNM(MXA3+MXNCF),QCN(MXA3+MXNCF),TGP(MXA3+MXNCF),
     *                BIS(MXA3+MXNCF),EBIS(10),XBIS(10),
     *                RTRMS(MXA3+MXFRG*21),K0(MXA3+MXNCF),
     *                K1(MXA3+MXNCF),K2(MXA3+MXNCF),K3(MXA3+MXNCF),
     *                MXOPT,IHUPD,IPIRC,NPIRC,NPRT,NPUN,STBSKP,GSCHL2
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINMOL
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      DATA SIRC/8HIRC     /
      DATA QNAM/8HSADDLE  ,8HFORWRD  ,8HSTABLZ  ,8HEVIB    ,
     *          8HSTRIDE  ,8HDELTA   ,8HSAB     ,8HSTOTAL  ,
     *          8HFREQ    ,8HCMODE   ,8HNPOINT  ,8HNEXTPT  ,
     *          8HPACE    ,8HTSENGY  ,8HELBOW   ,8HGA      ,
     *          8HREADQB  ,8HEB      ,8HGBNORM  ,8HGB      ,
     *          8HGA0     ,8HGA1     ,8HGA2     ,8HGCUT    ,
     *          8HRCUT    ,8HACUT    ,8HMXOPT   ,8HIHUPD   ,
     *          8HOPTTOL  ,8HNPRT    ,8HNPUN    /
      DATA KQNAM/0,0,0,3,3,3,3,3,3,-3,1,1,5,0,3,-3,0,3,3,-3,
     *           -3,-3,-3,3,3,3,1,1,3,1,1/
C
      DATA CHECK/8HCHECK   /
      DATA ALINR/8HLINEAR  /, QUAD/8HQUAD    /
      DATA AMPC4/8HAMPC4   /,  RK4/8HRK4     /,  GS2/8HGS2     /
C
C     ----- INTRINSIC REACTION COORDINATE PATHFINDER -----
C
C                    WRITTEN BY MIKE SCHMIDT
C          NORTH DAKOTA STATE UNIVERSITY --- JUNE 8,1983
C
C     GIVEN A POINT ON THE INTRINSIC REACTION COORDINATE, THIS
C     CODE FINDS THE NEXT POINT ON THE IRC, BY THE ALGORITHM OF
C           K.ISHIDA, K.MOROKUMA, A.KOMORNICKI
C               J.CHEM. PHYS. 66,2153-2156(1977)
C     WITH MODIFICATIONS AS DESCRIBED BY
C           M.W.SCHMIDT, M.S.GORDON, AND M.DUPUIS,
C               J. AM. CHEM. SOC., 107,2585-2589(1985)
C     OR BY THE AMPC4 OR RK4 ALGORITHMS DESCRIBED BY
C           B.C.GARRETT, M.J.REDMON, R.STECKLER, D.G.TRUHLAR,
C           K.K.BALDRIDGE, D.BARTOL, M.W.SCHMIDT, M.S.GORDON
C               J.PHYS.CHEM. 92, 1476-1488(1988)
C           K.K.BALDRIDGE, M.S.GORDON, R.STECKLER, D.G.TRUHLAR,
C               J.PHYS.CHEM. 93, 5107-5119(1989)
C     WHICH WAS IMPLEMENTED BY KIM BALDRIDGE,
C     OR BY THE SECOND-ORDER ALGORITHM DESCRIBED BY
C           C.GONZALEZ AND H.B.SCHLEGEL
C               J.CHEM.PHYS.  90, 2154-2161(1989)
C               J.PHYS.CHEM.  94, 5523-5527(1990)
C               J.CHEM.PHYS.  95, 5853-5860(1991)
C     WHICH WAS IMPLEMENTED BY FRANK JENSEN AND SHUJUN SU.
C
      FGONLY = NUM.EQ.0  .AND.  NAT.EQ.0
C
      DO 10 I=1,NNAM
         IF(KQNAM(I).EQ.-3) KQNAM(I) = 10*(MXA3+MXNCF) + 3
   10 CONTINUE
      NCOORD = 3 * NAT
      NCF=6*NFRG
      NC1=NCOORD + NCF
      CALL VCLR(EA,1,21*(MXA3+MXNCF)+15*MXFRG+35)
C
C           READ NAMELIST $IRC
C
      PACE   = GS2
      STRIDE = 0.30D+00
      NPOINT = 1
      NEXTPT = 0
C
      SADDLE = .FALSE.
      FORWRD = .TRUE.
      EVIB = 0.0005D+00
      TSENGY = .FALSE.
      FREQ = ZERO
C
      STABLZ = .TRUE.
      STBSKP = .FALSE.
      DELTA = 0.025D+00
      ELBOW = 175.0D+00
C
      GCUT = ZERO
      RCUT = ZERO
      ACUT = 5.0D+00
      MXOPT = 20
      IHUPD = 2
      OPTTOL=0.0001D+00
C
      NPRT=0
      NPUN=0
C                      BIS IS USED AS TEMP STORAGE OF CMODE
      JRET = 0
      CALL NAMEIO(IR,JRET,SIRC,NNAM,QNAM,KQNAM,
     *            SADDLE,FORWRD,STABLZ,EVIB,STRIDE,DELTA,
     *            SAB,STOTAL,FREQ,BIS,NPOINT,NEXTPT,PACE,
     *            TSENGY,ELBOW,GA,STBSKP,EB,GBNORM,GB,
     *            GA0,GA1,GA2,GCUT,RCUT,ACUT,MXOPT,IHUPD,
     *            OPTTOL,NPRT,NPUN,
     *            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)
      EVIB = ABS(EVIB)
      STRIDE = ABS(STRIDE)
      GSCHL2 = PACE.EQ.GS2
      IF(PACE.EQ.GS2) THEN
         IF(GCUT.EQ.ZERO) THEN
            GCUT = 10.0D-05 * (STRIDE/0.3D+00)
            IF(GCUT.GT.20.0D-05) GCUT=20.0D-05
            IF(GCUT.LT. 5.0D-05) GCUT= 5.0D-05
         END IF
         IF(RCUT.EQ.ZERO) THEN
            RCUT = 10.0D-04 * (STRIDE/0.3D+00)
            IF(RCUT.GT.20.0D-04) RCUT=20.0D-04
            IF(RCUT.LT. 5.0D-04) RCUT= 5.0D-04
         END IF
      END IF
      IF(PACE.NE.ALINR) STABLZ = .FALSE.
      IF(NPRT.GT. 1) NPRT= 1
      IF(NPRT.LT.-2) NPRT=-2
      IF(NPUN.GT. 1) NPUN= 1
      IF(NPUN.LT.-2) NPUN=-2
C
      IF(MASWRK) THEN
         WRITE(IW,9000) PACE,STRIDE,NPOINT,NEXTPT,STOTAL,SAB,
     *                  NPRT,NPUN
         WRITE(IW,9020) SADDLE,TSENGY,FORWRD,EVIB,FREQ
         WRITE(IW,9040) STABLZ,ELBOW,DELTA,STBSKP,EB,GBNORM
         WRITE(IW,9060) GCUT,RCUT,ACUT,MXOPT,IHUPD
      END IF
      IF(PACE.NE.ALINR  .AND.  PACE.NE.QUAD  .AND.
     *   PACE.NE.AMPC4  .AND.  PACE.NE.RK4  .AND.  PACE.NE.GS2)  THEN
            IF(MASWRK) WRITE(IW,9100)
            CALL ABRT
      END IF
      ACUT = ACUT*ACOS(-ONE)/180.0D+00
      IF(PACE.EQ.AMPC4) NAMPC = 1
      IF(JRET.EQ.2) CALL ABRT
C
C     ----- OPEN SUPPLEMENTAL PUNCH FILE FOR IRC RESULTS -----
C
      IPIRC=4
      CALL SEQOPN(IPIRC,'IRCDATA','NEW',.FALSE.,'FORMATTED')
C
C     ----- PRINT ATOMIC MASSES USED, THEN GET TRIPLES
C           OF INVERT SQUARE ROOTS OF THE MASSES -----
C
      IF(MASWRK) WRITE(IW,9200)
      I0=0
      DO 100 I=1,NAT
         IF(MASWRK) WRITE(IW,9220) ANAM(I),BNAM(I),ZMASS(I)
         TEMP = ONE/SQRT(ZMASS(I))
         RTRMS(I0+1)=TEMP
         RTRMS(I0+2)=TEMP
         RTRMS(I0+3)=TEMP
         I0=I0+3
  100 CONTINUE
C
C     CALCULATE EFFECT FRAGMENTS MASS AND MOMENTS OF INERTIA
C
      CALL EFRMSS(NFRG,RTRMS(NCOORD+1))
C
C     ----- MOVE COORDS OF INPUT $DATA GROUP TO CURRENT POINT QB -----
C
      CALL DCOPY(NCOORD,C,1,QB,1)
C
C      ----- GENERATE A FIRST POINT ON THE IRC JUST OFF
C            THE SADDLE POINT (BIS HOLDS CMODE)-----
C
      IF(SADDLE) THEN
         CMODLN = SQRT(DDOT(NC1,BIS,1,BIS,1))
         NHOPT=1
         IF(CMODLN .GT. 0.001D+00) NHOPT=0
         IF(MASWRK) THEN
            IF(NHOPT.EQ.0) THEN
               WRITE(IW,9300) FREQ
               WRITE(IW,9320) (BIS(I),I=1,NC1)
            ELSE
               WRITE(IW,9340)
            END IF
         END IF
         IF(TSENGY) THEN
            E=ZERO
            IF(MASWRK) WRITE(IW,9360)
            CALL ENERGX
            IF (EXETYP .EQ. CHECK) GO TO 500
            IF(E.EQ.ZERO) CALL ABRT
            EA = E
         END IF
         CALL INIIRM(EVIB,FORWRD,NCOORD,NCF,FREQ,BIS,NHOPT)
         NEXTPT=1
      END IF
C
C     ----- GENERATE ADDITIONAL POINTS ON THE IRC -----
C     GENERATE OUTPUT AT NEXT-TO-LAST POINT IF THE LAST POINT
C     IS GOING TO BE PREDICTED, NOT COMPUTED.
C
      CVGED = .FALSE.
      NPIRC = NEXTPT
  200 CONTINUE
      IF(PACE.EQ.GS2) THEN
         FIRST = NPIRC.EQ.NEXTPT
         LAST  = NPIRC-NEXTPT+1 .EQ. NPOINT
         CALL IRCGSM(NCOORD,NCF,FIRST,LAST,SADDLE,CVGED)
      ELSE
         FIRST = NPIRC.EQ.NEXTPT
         PENULT= NPIRC-NEXTPT+2 .EQ. NPOINT
         LAST  = NPIRC-NEXTPT+1 .EQ. NPOINT
         IF(PACE.EQ.RK4) LAST = PENULT
         CALL IRCSTP(STRIDE,NCOORD,PACE,SADDLE,STABLZ,DELTA,ELBOW,
     *               NAMPC,FIRST,LAST,PENULT)
         CALL DCOPY(NCOORD,QB,1,C,1)
      END IF
C
C     ----- PRINT/PUNCH PREDICTED COORDINATES OF NEXT IRC POINT ---
C     THE PREDICTED NEXT POINT IS NOW IN C.
C     EXCEPT FOR GS2 METHOD, QA HOLDS THE LAST IRC POINT WHERE
C     ENERGY/GRADIENT WAS ACTUALLY COMPUTED.
C
      IF(EXETYP.EQ.CHECK) GO TO 600
C
      EN = ENUC(NAT,ZAN,C)
      IF(NZMAT.GT.0) THEN
         CALL BANDBI
         CALL PZANDG(0,0)
      END IF
C
      IF (MASWRK) THEN
         IF (.NOT.FGONLY) WRITE(IW,9420) EN
         WRITE(IW,9400)
         WRITE(IP,9400)
         CALL PRATM(IW,2)
         CALL PRATM(IP,2)
      END IF
C
C                $IRC FOR GS2, AND A HESSIAN
C
      IF(GSCHL2) THEN
         IF (MASWRK) THEN
            WRITE(IW,9430) PACE,STRIDE,STOTAL,NPIRC,NPRT,NPUN
            WRITE(IP,9430) PACE,STRIDE,STOTAL,NPIRC,NPRT,NPUN
            I0 = 0
            DO 330 I=1,NAT+NCF/3
               I1 = I0 + 1
               WRITE(IW,9460) I1,(GA(I0+J),J=1,3)
               WRITE(IP,9460) I1,(GA(I0+J),J=1,3)
               I0 = I0 + 3
  330       CONTINUE
            WRITE(IW,9490)
            WRITE(IP,9490)
C
            WRITE(IW,9410)
         END IF
         CALL VALFM(LOADFM)
         LHSS = LOADFM + 1
         LASTM = LHSS + NC1**2
         NEED = LASTM - LOADFM - 1
         CALL GETFM(NEED)
         CALL DAREAD(IDAF,IODA,X(LHSS),NC1**2,4,0)
         CALL FCMPUN(X(LHSS),NC1)
         CALL RETFM(NEED)
         IF(CVGED) GO TO 600
         GO TO 500
      END IF
C
      IF (.NOT. MASWRK) GO TO 500
C
C            $IRC FOR ALL OTHER IRC STEPS
C
      WRITE(IW,9430) PACE,STRIDE,STOTAL,NPIRC,NPRT,NPUN
      WRITE(IP,9430) PACE,STRIDE,STOTAL,NPIRC,NPRT,NPUN
      IF(STABLZ) THEN
         WRITE(IW,9440) DELTA,ELBOW
         WRITE(IP,9440) DELTA,ELBOW
      END IF
      IF(PACE.EQ.QUAD) THEN
         WRITE(IW,9450) SAB
         WRITE(IP,9450) SAB
         I0 = 0
         DO 440 I=1,NAT
            I1 = I0 + 1
            WRITE(IW,9460) I1,(GA(I0+J),J=1,3)
            WRITE(IP,9460) I1,(GA(I0+J),J=1,3)
            I0 = I0 + 3
  440    CONTINUE
      END IF
      IF(PACE.EQ.AMPC4) THEN
         WRITE(IW,9450) SAB
         WRITE(IP,9450) SAB
         I0 = 0
         DO 442 I=1,NAT
            I1 = I0 + 1
            WRITE(IW,9462) I1,(GA0(I0+J),J=1,3)
            WRITE(IP,9462) I1,(GA0(I0+J),J=1,3)
            I0 = I0 + 3
  442    CONTINUE
         I0 = 0
         DO 444 I=1,NAT
            I1 = I0 + 1
            WRITE(IW,9464) I1,(GA1(I0+J),J=1,3)
            WRITE(IP,9464) I1,(GA1(I0+J),J=1,3)
            I0 = I0 + 3
  444    CONTINUE
      I0 = 0
      DO 446 I=1,NAT
         I1 = I0 + 1
         WRITE(IW,9466) I1,(GA2(I0+J),J=1,3)
         WRITE(IP,9466) I1,(GA2(I0+J),J=1,3)
         I0 = I0 + 3
  446 CONTINUE
      I0 = 0
      DO 448 I=1,NAT
         I1 = I0 + 1
         WRITE(IW,9468) I1,(GB(I0+J),J=1,3)
         WRITE(IP,9468) I1,(GB(I0+J),J=1,3)
         I0 = I0 + 3
  448 CONTINUE
      END IF
      IF(STBSKP) THEN
         WRITE(IW,9470) EB,GBNORM
         WRITE(IP,9470) EB,GBNORM
         I0 = 0
         DO 450 I=1,NAT
            I1 = I0 + 1
            WRITE(IW,9480) I1,(GB(I0+J),J=1,3)
            WRITE(IP,9480) I1,(GB(I0+J),J=1,3)
            I0 = I0 + 3
  450    CONTINUE
      END IF
      WRITE(IW,9490)
      WRITE(IP,9490)
C
  500 CONTINUE
      IF(PACE.EQ.GS2) THEN
         IF(NPIRC-NEXTPT  .LT.NPOINT) GO TO 200
      ELSE
         IF(NPIRC-NEXTPT+1.LT.NPOINT) GO TO 200
      END IF
C
  600 CONTINUE
      IF(MASWRK) WRITE(IW,9500)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/10X,29(1H-)/10X,'INTRINSIC REACTION COORDINATE'/
     *        10X,29(1H-)/
     *        10X,'PACE   =',A8/
     *        10X,'STRIDE =',F8.5,' BOHR*SQRT(AMU)'/
     *        10X,'NPOINT =',I8/
     *        10X,'NEXTPT =',I8/
     *        10X,'STOTAL =',F8.4,' BOHR*SQRT(AMU)'/
     *        10X,'SAB    =',F8.4,' BOHR*SQRT(AMU)'/
     *        10X,'NPRT   =',I8/
     *        10X,'NPUN   =',I8)
 9020 FORMAT(/10X,'SADDLE =',L8/
     *        10X,'TSENGY =',L8/
     *        10X,'FORWRD =',L8/
     *        10X,'EVIB   =',F8.5,' HARTREE'/
     *        10X,'FREQ   =',F8.1,' CM**-1')
 9040 FORMAT(/10X,'STABLZ =',L8/
     *        10X,'ELBOW  =',F8.2,' DEGREES'/
     *        10X,'DELTA  =',F8.5,' BOHR'//
     *        10X,'READQB =',L8/
     *        10X,'EB     =',F20.10,' HARTREES'/
     *        10X,'GBNORM =',F20.10,' HARTREES/SQRT(AMU)-BOHR')
 9060 FORMAT(/10X,'GCUT   =',F10.7,' HARTREES/SQRT(AMU)-BOHR'/
     *        10X,'RCUT   =',F10.7,' BOHR'/
     *        10X,'ACUT   =',F8.2,' DEGREES'/
     *        10X,'MXOPT  =',I8/
     *        10X,'IHUPD  =',I8)
C
 9100 FORMAT(1X,'*** ILLEGAL PACE WAS REQUESTED')
C
 9200 FORMAT(/1X,'ATOMIC MASSES USED (IN AMU) ARE')
 9220 FORMAT(5X,A8,A2,1X,F10.4)
C
 9300 FORMAT(/1X,'JUMP AWAY FROM SADDLE POINT ALONG MODE WITH',
     *       ' FREQUENCY',F11.5,' INPUT CMODE=')
 9320 FORMAT(1X,3F15.8,5X,3F15.8)
 9340 FORMAT(/1X,'JUMP AWAY FROM SADDLE POINT WILL BE FOUND FROM',
     *       ' THE INPUT HESSIAN MATRIX')
 9360 FORMAT(/1X,'COMPUTING WAVEFUNCTION AT TRANSITION STATE...')
C
 9400 FORMAT(/1X,'* * * IRC RESTART INFORMATION * * *'/
     *       1X,'COORDINATES FOR $DATA AND ENTIRE $IRC GROUP FOLLOW'/)
 9410 FORMAT(1X,'A GS2 RESTART ALSO REQUIRES THE UPDATED',
     *          ' $HESS GROUP FROM THE PUNCH FILE')
 9420 FORMAT(/1X,'THE NUCLEAR REPULSION AT THE IRC POINT IS',F20.10)
 9430 FORMAT(' $IRC   PACE=',A8,'   STRIDE=',F8.5,'   NPOINT=??'/
     *       '        STOTAL=',F10.6,'   NEXTPT=',I4/
     *       '        NPRT=',I2,'   NPUN=',I2,'   SADDLE=.FALSE.')
 9440 FORMAT('        STABLZ=.T.   DELTA=',F8.5,'   ELBOW=',F8.3)
 9450 FORMAT('        SAB=',F10.6)
 9460 FORMAT('        GA(',I3,')=',1P,E16.9,',',E16.9,',',E16.9)
 9462 FORMAT('        GA0(',I3,')=',1P,E16.9,',',E16.9,',',E16.9)
 9464 FORMAT('        GA1(',I3,')=',1P,E16.9,',',E16.9,',',E16.9)
 9466 FORMAT('        GA2(',I3,')=',1P,E16.9,',',E16.9,',',E16.9)
 9468 FORMAT('        GB(',I3,')=',1P,E16.9,',',E16.9,',',E16.9)
 9470 FORMAT('        READQB=.T.    EB=',F20.10,'   GBNORM=',1P,E16.9)
 9480 FORMAT('        GB(',I3,')=',1P,E16.9,',',E16.9,',',E16.9)
 9490 FORMAT(' $END   ')
 9500 FORMAT(1X,'......DONE WITH THIS PORTION OF THE IRC......')
      END
C*MODULE RXNCRD  *DECK QFIT1V
      SUBROUTINE QFIT1V(E1,G1,E2,C0,C1,C2,STEP,IW,FIT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL FIT,MASWRK,GOPARR,DSKWRK
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- FIT PARABOLA TO ENERGY AND GRADIENT AT QC,
C           AND ENERGY AT STEP ALONG VECTOR BIS -----
C
      C0 = E1
      C1 = G1
      C2 = ((E2-E1) - STEP*G1)/(STEP*STEP)
      FIT = .TRUE.
C
C     ----- TEST FIT FOR RELIABILITY -----
C
      XMIN = -C1/(C2+C2)
      IF(XMIN.GT.(STEP+STEP)) FIT = .FALSE.
      IF(XMIN.LT.-STEP)       FIT = .FALSE.
      IF(FIT) RETURN
C
      EMIN = C2 * XMIN*XMIN  + C1 * XMIN  + C0
      IF(MASWRK) WRITE(IW,9000) XMIN,EMIN,C0,C1,C2
      RETURN
C
 9000 FORMAT(10X,'PARABOLIC FIT TO 2 ENERGIES AND 1 GRADIENT',
     *      ' REJECTED DUE TO LARGE XMIN.'/
     *    10X,'XMIN,EMIN=',F12.6,F16.9,2X,
     *     'C0,C1,C2=',3F16.10/
     *   10X,'EXPLICIT LINEAR SEARCH WILL BE PERFORMED')
      END
C*MODULE RXNCRD  *DECK SEABIS
      SUBROUTINE SEABIS(DELTA,NCOORD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION K0,K1,K2,K3
C
      LOGICAL FIT,STBSKP,MASWRK,GOPARR,DSKWRK,GSCHL2
C
      PARAMETER (MXATM=2000, MXA3=3*MXATM, MXAO=8192,
     *           MXFRG=1050, MXNCF=6*MXFRG)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /FUNCT / E,EG(MXA3)
      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 /IRCDAT/ EA,EB,EC,SAB,SAC,STOTAL,GANORM,GBNORM,GCNORM,
     *                GBIS,GCUT,RCUT,ACUT,STRIDE,OPTTOL,
     *                QA(MXA3+MXNCF),QB(MXA3+MXNCF),QC(MXA3+MXNCF),
     *                GA(MXA3+MXNCF),GB(MXA3+MXNCF),GC(MXA3+MXNCF),
     *                QCM(MXA3+MXNCF),QBM(MXA3+MXNCF),QAM(MXA3+MXNCF),
     *                GA2(MXA3+MXNCF),GA1(MXA3+MXNCF),GA0(MXA3+MXNCF),
     *                QCNM(MXA3+MXNCF),QCN(MXA3+MXNCF),TGP(MXA3+MXNCF),
     *                BIS(MXA3+MXNCF),EBIS(10),XBIS(10),
     *                RTRMS(MXA3+MXFRG*21),K0(MXA3+MXNCF),
     *                K1(MXA3+MXNCF),K2(MXA3+MXNCF),K3(MXA3+MXNCF),
     *                MXOPT,IHUPD,IPIRC,NPIRC,NPRT,NPUN,STBSKP,GSCHL2
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C
C     ----- SEARCH ALONG THE LINE BIS FOR A MINIMUM ENERGY.
C           THIS MINIMUM WILL BE POINT QC' ON THE IRC -----
C
      NP = 1
      XBIS(NP) = ZERO
      EBIS(NP) = EC
      STEP = DELTA
C
C     ----- BEGIN LINEAR SEARCH -----
C
  100 CONTINUE
      NP = NP + 1
      IF(NP.GT.10 .AND. MASWRK) WRITE(IW,9000)
      IF(NP.GT.10) CALL ABRT
      DO 120 I=1,NCOORD
         C(I) = QC(I) + STEP * BIS(I)
  120 CONTINUE
      IF(MASWRK) THEN
         WRITE(IW,9020) NP,STEP
         WRITE(IW,9040)
         I0 = 0
         DO 140 I=1,NAT
            WRITE(IW,9060) ANAM(I),BNAM(I),(C(I0+J),J=1,3)
            I0 = I0 + 3
  140    CONTINUE
      END IF
C
C     ----- ENERGY EVALUATION -----
C
      E = ZERO
      CALL ENERGX
      IF(E.EQ.ZERO .AND. MASWRK) WRITE(IW,9080)
      IF(E.EQ.ZERO) CALL ABRT
      EBIS(NP) = E
      XBIS(NP) = STEP
      FIT = .FALSE.
      C0 = ZERO
      C1 = ZERO
      C2 = ZERO
      IF(NP.EQ.2) CALL QFIT1V(EC,GBIS,E,C0,C1,C2,STEP,IW,FIT)
      IF(FIT) GO TO 320
      IF(EBIS(NP).GT.EBIS(NP-1)) GO TO 160
         STEP = STEP + STEP
         GO TO 100
C
  160 CONTINUE
      IF(NP.GE.3) GO TO 300
         STEP = -STEP
         E = EBIS(1)
         X = XBIS(1)
         EBIS(1) = EBIS(2)
         XBIS(1) = XBIS(2)
         EBIS(2) = E
         XBIS(2) = X
         GO TO 100
C
C     ----- MINIMUM ENERGY ON THE LINE IS BRACKETED.
C           FIT PARABOLA TO LAST THREE POINTS -----
C
  300 CONTINUE
      X1 = XBIS(NP-2)
      F1 = EBIS(NP-2)
      X2 = XBIS(NP-1)
      F2 = EBIS(NP-1)
      X3 = XBIS(NP  )
      F3 = EBIS(NP  )
C
      C2 = ((F1-F2)*(X2-X3) - (F2-F3)*(X1-X2))
      C2 = C2/((X1*X1-X2*X2)*(X2-X3) - (X2*X2-X3*X3)*(X1-X2))
      C1 = (F1 - F2 + C2 * (X2*X2-X1*X1))/(X1-X2)
      C0 = F1 - C2*X1*X1 - C1*X1
C
C     ----- PREDICT EXTREMUM OF THE PARABOLA -----
C
  320 CONTINUE
      IF(MASWRK) WRITE(IW,9100) C0,C1,C2
      XMIN = -C1/(C2+C2)
      EMIN = C2 * XMIN*XMIN + C1 * XMIN + C0
      IF(MASWRK) WRITE(IW,9120) XMIN,EMIN
C
C     ----- GENERATE COORDINATES OF POINT QC' (IN QC) AT THE
C           PREDICTED MINIMUM OF THE LINEAR SEARCH -----
C
      DO 400 I=1,NCOORD
         QC(I) = QC(I) + XMIN * BIS(I)
  400 CONTINUE
      RETURN
C
 9000 FORMAT(1X,'**** ERROR, LINEAR SEARCH TAKES TOO',
     *          ' MANY POINTS TO FIND MINIMUM')
 9020 FORMAT(///1X,9X,'POINT',I3,' ON LINEAR SEARCH...',
     *          'STEP=',F9.5,' BOHR')
 9040 FORMAT(/10X,'ATOM',16X,'X',14X,'Y',14X,'Z'/
     *        10X,55(1H-)/)
 9060 FORMAT(10X,A8,A2,3F15.6)
 9080 FORMAT(1X,'***** ERROR, SCF FAILED TO CONVERGE')
 9100 FORMAT(//10X,'THE PARABOLIC FIT IS -'/
     *       10X,F16.10,' + ',F16.10,' * STEP +',F16.10,
     *           ' * STEP**2 = E(STEP)')
 9120 FORMAT(/10X,'THE MINIMUM OF THE FIT IS AT STEP=',
     *        F12.6,' WITH A PREDICTED ENERGY OF',F16.9/
     *        10X,'STABILIZED POINT QC'' ON THE IRC HAS BEEN FOUND')
      END
C*MODULE RXNCRD  *DECK STBIMK
      SUBROUTINE STBIMK(DELTA,ELBOW,NCOORD,FIRST,PENULT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION K0,K1,K2,K3
C
      LOGICAL STBSKP,LINMOL,MASWRK,GOPARR,DSKWRK,GSCHL2,
     *        FIRST,PENULT,EXTPRT,EXTPUN,ALLPUN,MMONLY,QMMM
C
      PARAMETER (MXATM=2000, MXA3=3*MXATM, MXAO=8192,
     *           MXFRG=1050, MXNCF=6*MXFRG)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, ONE80=1.8D+02, GTOL=1.0D-08)
C
      COMMON /FUNCT / E,EG(MXA3)
      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 /IRCDAT/ EA,EB,EC,SAB,SAC,STOTAL,GANORM,GBNORM,GCNORM,
     *                GBIS,GCUT,RCUT,ACUT,STRIDE,OPTTOL,
     *                QA(MXA3+MXNCF),QB(MXA3+MXNCF),QC(MXA3+MXNCF),
     *                GA(MXA3+MXNCF),GB(MXA3+MXNCF),GC(MXA3+MXNCF),
     *                QCM(MXA3+MXNCF),QBM(MXA3+MXNCF),QAM(MXA3+MXNCF),
     *                GA2(MXA3+MXNCF),GA1(MXA3+MXNCF),GA0(MXA3+MXNCF),
     *                QCNM(MXA3+MXNCF),QCN(MXA3+MXNCF),TGP(MXA3+MXNCF),
     *                BIS(MXA3+MXNCF),EBIS(10),XBIS(10),
     *                RTRMS(MXA3+MXFRG*21),K0(MXA3+MXNCF),
     *                K1(MXA3+MXNCF),K2(MXA3+MXNCF),K3(MXA3+MXNCF),
     *                MXOPT,IHUPD,IPIRC,NPIRC,NPRT,NPUN,STBSKP,GSCHL2
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /TINOPT/ MMONLY,QMMM
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINMOL
C
      DATA CHECK/8HCHECK   /
      DATA GUGA   /8HGUGA    /
C
C     ----- STABILIZE THE REACTION PATH BY SETTING UP THE
C           SEARCH ALONG THE BISECTOR AS DESCRIBED BY
C           ISHIDA, MOROKUMA, AND KOMORNICKI -----
C
      EXTPRT = NPRT.EQ.1  .OR.  (NPRT.EQ. 0  .AND. (FIRST.OR.PENULT))
      EXTPUN = NPUN.GE.0  .OR.  (NPUN.EQ.-1  .AND. (FIRST.OR.PENULT))
C
      IF(MASWRK) THEN
         WRITE(IW,9000)
         I0 = 0
         DO 100 I=1,NAT
            WRITE(IW,9020) ANAM(I),BNAM(I),(QC(I0+J),J=1,3)
            I0 = I0 + 3
  100    CONTINUE
      END IF
C
C     ----- ENERGY AT POINT QC -----
C     QC IS THE POINT NOT QUITE ON THE PATH, WHICH THE BISECTOR
C     WILL PASS THROUGH.
C
      DO 120 I=1,NCOORD
         C(I) = QC(I)
  120 CONTINUE
      E = ZERO
      CALL ENERGX
      IF (E .EQ. ZERO) THEN
         IF(MASWRK) WRITE(IW,9040)
         CALL ABRT
      END IF
      EC = E
C
C     ----- GRADIENT AT POINT QC -----
C
      DO 400 I=1,NCOORD
         EG(I) = ZERO
  400 CONTINUE
      IF (CITYP.EQ.GUGA) THEN
         CALL CIGRAD
      ELSE
         CALL HFGRAD
      END IF
      IF (QMMM) CALL TOYS
      IF (EXETYP .EQ. CHECK) RETURN
C
C     ----- NORMALIZED, MASS WEIGHTED GRADIENT AT POINT QC -----
C
      GCNORM = ZERO
      DO 420 I=1,NCOORD
         G = RTRMS(I) * EG(I)
         IF(ABS(G).LT.GTOL) G = ZERO
         GCNORM = GCNORM + G * G
         GC(I) = G
  420 CONTINUE
      GCNORM = SQRT(GCNORM)
      IF(MASWRK) WRITE(IW,9060) GCNORM
      FACT = ONE/GCNORM
      DO 440 I=1,NCOORD
         GC(I) = FACT * GC(I)
  440 CONTINUE
      IF(MASWRK) THEN
         WRITE(IW,9080)
         I0 = 0
         DO 460 I=1,NAT
            WRITE(IW,9100) ANAM(I),BNAM(I),(GC(I0+J),J=1,3)
            I0 = I0 + 3
  460    CONTINUE
      END IF
C
C     ----- TEST GRADIENT AT POINT QC AND QB FOR COLLINEARITY -----
C     IF THESE ARE SUFFICIENTLY PARALLEL, THE STABILIZATION IS SKIPPED.
C     NOTE THAT GB AND GC ARE ALREADY NORMALIZED.
C
      DOTPRD = ZERO
      DO 500 I=1,NCOORD
         DOTPRD = DOTPRD + GC(I)*GB(I)
  500 CONTINUE
      THETA = ACOS(DOTPRD)
      PI = ACOS(-ONE)
      THETA = ONE80*THETA/PI
      THETA = ONE80 - THETA
      IF(MASWRK) WRITE(IW,9120) THETA
      IF(THETA.LT.ELBOW) GO TO 600
C
C     ----- STABILIZATION IS NOT NECESSARY, AS WE NOW HAVE
C           FOUND THE NEXT IRC POINT.  PRINT AND EXIT -----
C
      NEXT = NPIRC + 1
      IF(MASWRK) WRITE(IW,9140) ELBOW,NEXT
      STBSKP = .TRUE.
      IF(NZMAT.GT.0) THEN
         IF(MASWRK) WRITE(IW,9160)
         CALL BANDBI
         CALL TRANG(EG,NVAR,NCOORD)
         CALL PZANDG(EG,1)
      END IF
      IF(EXTPRT .OR. EXTPUN) THEN
         ALLPUN = FIRST  .OR.  PENULT  .OR.  (NPUN.EQ.1)
         CALL IRCMOS(EXTPRT,EXTPUN,ALLPUN,NPIRC+1,STOTAL+STRIDE,EC)
      END IF
      RETURN
C     ******
C
C     ----- DETERMINE UNIT VECTOR ALONG BISECTOR
C           OF THE UNIT VECTORS ALONG -GB AND -GC -----
C
  600 CONTINUE
      IF(MASWRK) WRITE(IW,9180)
      STBSKP = .FALSE.
      DO 610 I=1,NCOORD
         BIS(I) = GB(I) - GC(I)
  610 CONTINUE
      CALL SYMDR(BIS)
      BISNRM = ZERO
      DO 620 I=1,NCOORD
         BISNRM = BISNRM + BIS(I) * BIS(I)
  620 CONTINUE
      BISNRM = SQRT(BISNRM)
      IF(BISNRM.LT.GTOL) CALL ABRT
      FACT = ONE/BISNRM
      DO 640 I=1,NCOORD
         BIS(I) = FACT * BIS(I)
  640 CONTINUE
C
C     ----- GET COMPONENT OF GRADIENT ALONG BISECTOR DIRECTION -----
C     EG STILL HOLDS UN-MASS-WEIGHTED, UN-NORMALIZED GRADIENT AT QC.
C
      GBIS = ZERO
      DO 700 I=1,NCOORD
         GBIS = GBIS + EG(I) * BIS(I)
  700 CONTINUE
      IF(MASWRK) THEN
         WRITE(IW,9200) GBIS
         WRITE(IW,9220)
         I0 = 0
         DO 720 I=1,NAT
            WRITE(IW,9100) ANAM(I),BNAM(I),(BIS(I0+J),J=1,3)
            I0 = I0 + 3
  720    CONTINUE
      END IF
C
C     ----- SEARCH ALONG THE BISECTOR FOR THE ENERGY MINIMUM QC' -----
C
      CALL SEABIS(DELTA,NCOORD)
      RETURN
C
 9000 FORMAT(/10X,'THE COORDINATES OF POINT QC (BOHRS) ARE'/
     *        10X,'ATOM',16X,'X',14X,'Y',14X,'Z'/
     *        10X,59(1H-)/)
 9020 FORMAT(10X,A8,A2,4X,3F15.10)
 9040 FORMAT(1X,'**** ERROR, SCF FAILED TO CONVERGE')
 9060 FORMAT(/10X,'THE NORMALIZED MASS WEIGHTED GRADIENT AT POINT QC'/
     *        10X,'(ORIGINAL GCNORM=',F10.6,' HARTREE/BOHR-SQRT(AMU))')
 9080 FORMAT(/10X,'ATOM',12X,'DE/DX',7X,'DE/DY',7X,'DE/DZ'/
     *        10X,46(1H-))
 9100 FORMAT(10X,A8,A2,3(2X,F10.6))
 9120 FORMAT(/5X,'THE ANGLE BETWEEN GRADIENTS GB AND GC IS',F8.2,
     *         ' DEGREES')
 9140 FORMAT(/5X,'* * * STABILIZATION IS BEING SKIPPED * * *'/
     *        5X,'      (POINT QC IS TAKEN TO BE QC''.)'/
     *        5X,'BECAUSE THE ANGLE BETWEEN GRADIENTS EXCEEDS',
     *           ' ELBOW ANGLE',F8.2/
     *        5X,'THE PRECEEDING POINT IS THEREFORE POINT',I4,
     *           ' ON THE REACTION PATH.'/
     *        5X,'WHEN THIS POINT IS RELABELED FROM QC'' TO QB,'/
     *        5X,'ITS ENERGY AND GRADIENT WILL BE ALREADY KNOWN,'/
     *        5X,'SO THAT WE CAN PROCEED IMMEDIATELY ',
     *           'TO LOOK FOR THE NEXT IRC POINT.')
 9160 FORMAT(5X,'THE INTERNAL COORDINATES OF THE PRECEEDING POINT ARE')
 9180 FORMAT(///1X,'*** ISHIDA, MOROKUMA, KOMORNICKI STABILIZATION',
     *              ' OF THE REACTION PATH ***'/
     *       5X,'THE POINT QC ABOVE IS AT STEP=0.0 BOHR'/
     *       5X,'LOOKING FOR MINIMUM ENERGY POINT QC'' ALONG THE',
     *           ' GRADIENT BISECTOR')
 9200 FORMAT(5X,'THE COMPONENT OF THE GRADIENT AT QC ALONG',
     *     ' THE BISECTOR DIRECTION IS',F11.6,' HARTREE/BOHR')
 9220 FORMAT(/10X,'THE UNIT VECTOR ALONG THE BISECTOR',
     *            ' DIRECTION IS (DIMENSIONLESS)'//
     *            10X,'ATOM',12X,'X',11X,'Y',11X,'Z'/
     *            10X,50(1H-))
      END
C*MODULE RXNCRD  *DECK IRCMOS
      SUBROUTINE IRCMOS(EXTPRT,EXTPUN,ALLPUN,NPIRC,STOTAL,ENERGY)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      LOGICAL EXTPRT,EXTPUN,ALLPUN,GOPARR,DSKWRK,MASWRK,
     *        CANONC,FCORE,FORS,NOCI,EKT,LINSER
C
      PARAMETER (MXATM=2000, MXNORO=250, MXRT=100)
C
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRITDT,PRTTOL,SDET,SZDET,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCORDT,NACTDT,
     *                NORBDT,NADET,NBDET,KDET,KSTDET,IROOT,IPURES,
     *                MAXW1,NITDET,MAXPDT,NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /FMCOM / X(1)
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      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 /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,NORB,
     *                NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DATA VECWRD,ENDWRD/8H $VEC   ,8H $END   /, UHF/8HUHF     /
      DATA RMCSCF/8HMCSCF   /
      DATA GUGA,ALDET,GENCI,ORMAS
     *       /8HGUGA    ,8HALDET   ,8HGENCI   ,8HORMAS   /
C
C        PRINT/PUNCH MOLECULAR ORBITALS AT AN IRC POINT
C        WHOSE WAVEFUNCTION HAS ALREADY BEEN EVALUATED WITH
C        PRINT AND PUNCH OPTIONS TURNED OFF.  THE ORBITALS
C        SHOULD ALREADY NOW BE ON THE DAF FILE.
C
      L0 = NQMT
      L1 = NUM
      L3 = L1*L1
C
      CALL VALFM(LOADFM)
      LVEC  = LOADFM + 1
      LEIG  = LVEC   + L3
      LSYMB = LEIG   + L1
      LAST  = LSYMB  + L1
      NEED = LAST - LOADFM -1
      CALL GETFM(NEED)
C
      NMOPRT = L0
      NMOPUN = NA
      IF(SCFTYP.EQ.RMCSCF) THEN
         IF(CISTEP.EQ.GUGA)  NMOPUN=NMCC+NDOC+NALP+NAOS+NBOS+NVAL
         IF(CISTEP.EQ.ALDET) NMOPUN=NCORSV+NACTDT
         IF(CISTEP.EQ.ORMAS) NMOPUN=NCORSV+NACTDT
         IF(CISTEP.EQ.GENCI) NMOPUN=NCORSV+NACTDT
      END IF
      IF(CITYP.EQ.GUGA  .AND.  ICICI.EQ.1) NMOPUN=L0
      IF(ALLPUN) NMOPUN=L0
C
      IF(MASWRK) THEN
         IF(EXTPRT) THEN
            WRITE(IW,*) ' '
            WRITE(IW,*) 'CONVERGED ORBITALS AT IRC POINT',NPIRC
            IF(SCFTYP.EQ.UHF) WRITE(IW,*) ' '
            IF(SCFTYP.EQ.UHF) WRITE(IW,*) 'ALPHA ORBITALS'
         END IF
         IF(EXTPUN) THEN
            WRITE(IP,8500) NPIRC,STOTAL,ENERGY,SCFTYP
            WRITE(IP,8510) VECWRD
         END IF
      END IF
C
      NDAF = 15
      IPASS= 0
  675 CONTINUE
      IPASS=IPASS+1
      CALL DAREAD(IDAF,IODA,X(LVEC),L3,NDAF  ,0)
      IF(EXTPRT) THEN
         IF(IPASS.EQ.1) CALL DAREAD(IDAF,IODA,X(LSYMB),L1,255,1)
         IF(IPASS.EQ.2) CALL DAREAD(IDAF,IODA,X(LSYMB),L1,256,1)
         CALL DAREAD(IDAF,IODA,X(LEIG),L1,NDAF+2,0)
         CALL PREVS(X(LVEC),X(LEIG),X(LSYMB),NMOPRT,L1,L1,1)
      END IF
      IF(EXTPUN) CALL PUSQL(X(LVEC),NMOPUN,L1,L1)
      IF(SCFTYP.EQ.UHF) THEN
         NMOPUN=NB
         IF(ALLPUN) NMOPUN=L1
         NDAF = 19
         IF(MASWRK) WRITE(IW,*) ' '
         IF(MASWRK) WRITE(IW,*) 'BETA ORBITALS'
         IF(IPASS.LT.2) GO TO 675
      END IF
      IF(MASWRK  .AND.  EXTPUN) WRITE(IP,8510) ENDWRD
      CALL RETFM(NEED)
      RETURN
C
 8500 FORMAT('ORBITALS FROM IRC POINT',I5,', STOTAL=',F10.5/
     *       'ENERGY=',F20.10,', SCFTYP=',A8)
 8510 FORMAT(A8)
      END
