      SUBROUTINE REACT1(ESCF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL PASTOR, PBSTOR
      INCLUDE 'SIZES'
      COMMON /GEOM  / GEO(3,NUMATM), XCOORD(3,NUMATM)
      DIMENSION GEOA(3,NUMATM), GEOVEC(3,NUMATM),
     1          PASTOR(MPACK),
     2          PBSTOR(MPACK), XOLD(MAXPAR), GROLD(MAXPAR)
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1                NA(NUMATM), NB(NUMATM), NC(NUMATM)
      COMMON /DENSTY/ P(MPACK),PA(MPACK),PB(MPACK)
      COMMON /GEOSYM/ NDEP,LOCPAR(MAXPAR),IDEPFN(MAXPAR),LOCDEP(MAXPAR)
      COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
      COMMON /GRADNT/ GRAD(MAXPAR),GNORM
      COMMON /ISTOPE/ AMS(107)
      COMMON /GRAVEC/ COSINE
      COMMON /KEYWRD/ KEYWRD
      COMMON /MESAGE/ IFLEPO,ISCF
     1       /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     2                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     3                NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /REACTN/ STEP, GEOA, GEOVEC,CALCST
      LOGICAL GRADNT, FINISH, XYZ, INT, GOK(2)
      SAVE GRADNT, FINISH, XYZ,INT, GOK
************************************************************************
*
*  REACT1 DETERMINES THE TRANSITION STATE OF A CHEMICAL REACTION.
*
*   REACT WORKS BY USING TWO SYSTEMS SIMULTANEOUSLY, THE HEATS OF
*   FORMATION OF BOTH ARE CALCULATED, THEN THE MORE STABLE ONE
*   IS MOVED IN THE DIRECTION OF THE OTHER. AFTER A STEP THE
*   ENERGIES ARE COMPARED, AND THE NOW LOWER-ENERGY FORM IS MOVED
*   IN THE DIRECTION OF THE HIGHER-ENERGY FORM. THIS IS REPEATED
*   UNTIL THE SADDLE POINT IS REACHED.
*
*   IF ONE FORM IS MOVED 3 TIMES IN SUCCESSION, THEN THE HIGHER ENERGY
*   FORM IS RE-OPTIMIZED WITHOUT SHORTENING THE DISTANCE BETWEEN THE TWO
*   FORMS. THIS REDUCES THE CHANCE OF BEING CAUGHT ON THE SIDE OF A
*   TRANSITION STATE.
*
************************************************************************
      DIMENSION IDUM1(NUMATM), IDUM2(NUMATM), XSTORE(MAXPAR),
     1IDUM3(NUMATM), COORD(3,NUMATM), IROT(2,3)
      DIMENSION IDUMMY(3*NUMATM)
      SAVE IROT
      CHARACTER*241 KEYWRD
      EQUIVALENCE (IDUMMY,COORD)
      DATA IROT/1,2,1,3,2,3/
      GOLD=0.D0
      LINEAR=0
      IFLAG=1
      GOK(1)=.FALSE.
      GOK(2)=.FALSE.
      XYZ=(INDEX(KEYWRD,' XYZ') .NE. 0)
      GRADNT=(INDEX(KEYWRD,'GRAD') .NE. 0)
      I=(INDEX(KEYWRD,' BAR'))
      STEPMX=0.15D0
      IF(I.NE.0) STEPMX=READA(KEYWRD,I)
      MAXSTP=1000
C
C    READ IN THE SECOND GEOMETRY.
C
      IF(XYZ) THEN
         CALL GETGEO(5,LABELS,GEOA,LOC,NA,NB,NC,AMS,NATOMS,INT)
      ELSE
         CALL GETGEO(5,IDUM1,GEOA,IDUMMY,
     1         IDUM1,IDUM2,IDUM3,AMS,NATOMS,INT)
C
C  IF INTERNAL COORDINATES ARE TO BE USED, CHECK THE CONNECTIVITY
C
         L=0
         DO 10 I=1,NATOMS
            IF(IDUM1(I).NE.NA(I))THEN
               L=L+1
               IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY
     1'')')
               WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE BOND LAB
     1ELS ARE DIFFERENT:      ',IDUM1(I),' AND',NA(I)
            ENDIF
            IF(IDUM2(I).NE.NB(I))THEN
               L=L+1
               IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY
     1'')')
               WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE BOND ANG
     1LE LABELS ARE DIFFERENT:',IDUM2(I),' AND',NB(I)
            ENDIF
            IF(IDUM3(I).NE.NC(I))THEN
               L=L+1
               IF(L.EQ.1)WRITE(6,'(10X,''ERRORS DETECTED IN CONNECTIVITY
     1'')')
               WRITE(6,'(A,I3,A,I3,A,I3,A)')' FOR ATOM',I,' THE DIHEDRAL
     1 LABELS ARE DIFFERENT:  ',IDUM3(I),' AND',NC(I)
            ENDIF
   10    CONTINUE
         IF(L.NE.0)WRITE(6,'(10X,A)')' CORRECT BEFORE RESUBMISSION'
         IF(L.NE.0)STOP
      ENDIF
      TIME0= SECOND()
C
C  SWAP FIRST AND SECOND GEOMETRIES AROUND
C  SO THAT GEOUT CAN OUTPUT DATA ON SECOND GEOMETRY.
C
      NUMAT2=0
      DO 20 I=1,NATOMS
         IF(LABELS(I).NE.99) NUMAT2=NUMAT2+1
         X=GEOA(1,I)
         GEOA(1,I)=GEO(1,I)
         GEO(1,I)=X
         X=GEOA(2,I)*0.0174532925D0
         GEOA(2,I)=GEO(2,I)
         GEO(2,I)=X
         X=GEOA(3,I)*0.0174532925D0
         GEOA(3,I)=GEO(3,I)
         GEO(3,I)=X
   20 CONTINUE
      IF(NUMAT2.NE.NUMAT) THEN
         WRITE(6,'(//10X,'' NUMBER OF ATOMS IN SECOND SYSTEM IS '',
     1''INCORRECT'',/)')
         WRITE(6,'('' NUMBER OF ATOMS IN FIRST  SYSTEM ='',I4)')NUMAT
         WRITE(6,'('' NUMBER OF ATOMS IN SECOND SYSTEM ='',I4)')NUMAT2
         GOTO 280
      ENDIF
      WRITE(6,'(//10X,'' GEOMETRY OF SECOND SYSTEM'',/)')
      IF(NDEP.NE.0) CALL SYMTRY
      CALL GEOUT(1)
C
C     CONVERT TO CARTESIAN, IF NECESSARY
C
      IF(   XYZ   )THEN
         CALL GMETRY(GEO,COORD)
         SUMX=0.D0
         SUMY=0.D0
         SUMZ=0.D0
         DO 30 J=1,NUMAT
            SUMX=SUMX+COORD(1,J)
            SUMY=SUMY+COORD(2,J)
   30    SUMZ=SUMZ+COORD(3,J)
         SUMX=SUMX/NUMAT
         SUMY=SUMY/NUMAT
         SUMZ=SUMZ/NUMAT
         DO 40 J=1,NUMAT
            GEO(1,J)=COORD(1,J)-SUMX
            GEO(2,J)=COORD(2,J)-SUMY
   40    GEO(3,J)=COORD(3,J)-SUMZ
         WRITE(6,'(//,''  CARTESIAN GEOMETRY OF FIRST SYSTEM'',//)')
         WRITE(6,'(3F14.5)')((GEO(J,I),J=1,3),I=1,NUMAT)
         SUMX=0.D0
         SUMY=0.D0
         SUMZ=0.D0
         DO 50 J=1,NUMAT
            SUMX=SUMX+GEOA(1,J)
            SUMY=SUMY+GEOA(2,J)
   50    SUMZ=SUMZ+GEOA(3,J)
         SUM=0.D0
         SUMX=SUMX/NUMAT
         SUMY=SUMY/NUMAT
         SUMZ=SUMZ/NUMAT
         DO 60 J=1,NUMAT
            GEOA(1,J)=GEOA(1,J)-SUMX
            GEOA(2,J)=GEOA(2,J)-SUMY
            GEOA(3,J)=GEOA(3,J)-SUMZ
            SUM=SUM+(GEO(1,J)-GEOA(1,J))**2
     1           +(GEO(2,J)-GEOA(2,J))**2
     2           +(GEO(3,J)-GEOA(3,J))**2
   60    CONTINUE
         DO 110 L=3,1,-1
C
C     DOCKING IS DONE IN STEPS OF 16, 4, AND 1 DEGREES AT A TIME.
C
            CA=COS(4.D0**(L-1)*0.01745329D0)
            SA=SQRT(ABS(1.D0-CA**2))
            DO 100 J=1,3
               IR=IROT(1,J)
               JR=IROT(2,J)
               DO 90 I=1,10
                  SUMM=0.D0
                  DO 70 K=1,NUMAT
                     X         = CA*GEOA(IR,K)+SA*GEOA(JR,K)
                     GEOA(JR,K)=-SA*GEOA(IR,K)+CA*GEOA(JR,K)
                     GEOA(IR,K)=X
                     SUMM=SUMM+(GEO(1,K)-GEOA(1,K))**2
     1                         +(GEO(2,K)-GEOA(2,K))**2
     2                         +(GEO(3,K)-GEOA(3,K))**2
   70             CONTINUE
                  IF(SUMM.GT.SUM) THEN
                     IF(I.GT.1)THEN
                        SA=-SA
                        DO 80 K=1,NUMAT
                           X         = CA*GEOA(IR,K)+SA*GEOA(JR,K)
                           GEOA(JR,K)=-SA*GEOA(IR,K)+CA*GEOA(JR,K)
                           GEOA(IR,K)=X
   80                   CONTINUE
                        GOTO 100
                     ENDIF
                     SA=-SA
                  ENDIF
   90          SUM=SUMM
  100       CONTINUE
  110    CONTINUE
         WRITE(6,'(//,''  CARTESIAN GEOMETRY OF SECOND SYSTEM'',//)')
         WRITE(6,'(3F14.5)')((GEOA(J,I),J=1,3),I=1,NUMAT)
         WRITE(6,'(//,''   "DISTANCE":'',F13.6)')SUM
         WRITE(6,'(//,''  REACTION COORDINATE VECTOR'',//)')
         WRITE(6,'(3F14.5)')((GEOA(J,I)-GEO(J,I),J=1,3),I=1,NUMAT)
         NA(1)=99
         J=0
         NVAR=0
         DO 130 I=1,NATOMS
            IF(LABELS(I).NE.99)THEN
               J=J+1
               DO 120 K=1,3
                  NVAR=NVAR+1
                  LOC(2,NVAR)=K
  120          LOC(1,NVAR)=J
               LABELS(J)=LABELS(I)
            ENDIF
  130    CONTINUE
         NATOMS=NUMAT
      ENDIF
C
C   XPARAM HOLDS THE VARIABLE PARAMETERS FOR GEOMETRY IN GEO
C   XOLD   HOLDS THE VARIABLE PARAMETERS FOR GEOMETRY IN GEOA
C
      IF(NVAR.EQ.0)THEN
         WRITE(6,'(///10X,''THERE ARE NO VARIABLES IN THE SADDLE'',
     1'' CALCULATION!'')')
         STOP
      ENDIF
      SUM=0.D0
      DO 140 I=1,NVAR
         GROLD(I)=1.D0
         XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
         XOLD(I)=GEOA(LOC(2,I),LOC(1,I))
  140 SUM=SUM+(XPARAM(I)-XOLD(I))**2
      STEP0=SQRT(SUM)
      IF(STEP0.LT.1.D-5)THEN
         WRITE(6,'(//,3(5X,A,/))')' BOTH GEOMETRIES ARE IDENTICAL',
     1' A SADDLE CALCULATION INVOLVES A REACTANT AND A PRODUCT',
     2' THESE MUST BE DIFFERENT GEOMETRIES'
         STOP
      ENDIF
      ONE=1.D0
      DELL=0.1D0
      EOLD=-2000.D0
      TIME1=SECOND()
      SWAP=0
      DO 240 ILOOP=1,MAXSTP
         WRITE(6,'('' '',40(''*+''))')
C
C   THIS METHOD OF CALCULATING 'STEP' IS QUITE ARBITARY, AND NEEDS
C   TO BE IMPROVED BY INTELLIGENT GUESSWORK!
C
         IF (GNORM.LT.1.D-3)GNORM=1.D-3
         STEP=MIN(SWAP,0.5D0, 6.D0/GNORM, DELL,STEPMX*STEP0+0.005D0)
         STEP=MIN(0.2D0,STEP/STEP0)*STEP0
         SWAP=SWAP+1.D0
         DELL=DELL+0.1D0
         WRITE(6,'(''  BAR SHORTENED BY'',F12.7,'' PERCENT'')')
     1STEP/STEP0*100.D0
         STEP0=STEP0-STEP
         IF(STEP0.LT.0.01D0) GOTO 250
         STEP=STEP0
         DO 150 I=1,NVAR
  150    XSTORE(I)=XPARAM(I)
         CALL FLEPO(XPARAM, NVAR, ESCF)
         IF(LINEAR.EQ.0)THEN
            LINEAR=(NORBS*(NORBS+1))/2
            DO 160 I=1,LINEAR
               PASTOR(I)=PA(I)
  160       PBSTOR(I)=PB(I)
         ENDIF
         DO 170 I=1,NVAR
  170    XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
         IF(IFLAG.EQ.1)THEN
            WRITE(6,'(//10X,''FOR POINT'',I3,'' SECOND STRUCTURE'')')ILO
     1OP
         ELSE
            WRITE(6,'(//10X,''FOR POINT'',I3,'' FIRST  STRUCTURE'')')ILO
     1OP
         ENDIF
         WRITE(6,'('' DISTANCE A - B  '',F12.6)')STEP
C
C   NOW TO CALCULATE THE "CORRECT" GRADIENTS, SWITCH OFF 'STEP'.
C
         STEP=0.D0
         DO 180 I=1,NVAR
  180    GRAD(I)=GROLD(I)
         CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.)
         DO 190 I=1,NVAR
  190    GROLD(I)=GRAD(I)
         IF (GRADNT) THEN
            WRITE(6,'(''  ACTUAL GRADIENTS OF THIS POINT'')')
            WRITE(6,'(8F10.4)')(GRAD(I),I=1,NVAR)
         ENDIF
         WRITE(6,'('' HEAT            '',F12.6)')FUNCT1
         GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
         WRITE(6,'('' GRADIENT NORM   '',F12.6)')GNORM
         COSINE=COSINE*ONE
         WRITE(6,'('' DIRECTION COSINE'',F12.6)')COSINE
         CALL GEOUT(6)
         IF(SWAP.GT.2.9D0 .OR. ILOOP .GT. 3 .AND. COSINE .LT. 0.D0
     1  .OR. ESCF .GT. EOLD)
     2  THEN
            IF(SWAP.GT.2.9D0)THEN
               SWAP=0.D0
            ELSE
               SWAP=0.5D0
            ENDIF
C
C   SWAP REACTANT AND PRODUCT AROUND
C
            FINISH=(GOK(1).AND.GOK(2) .AND. COSINE .LT. 0.D0)
            IF(FINISH) THEN
               WRITE(6,'(//10X,'' BOTH SYSTEMS ARE ON THE SAME SIDE OF T
     1HE '',''TRANSITION STATE -'',/10X,'' GEOMETRIES OF THE SYSTEMS'',
     2'' ON EACH SIDE OF THE T.S. ARE AS FOLLOWS'')')
               DO 200 I=1,NVAR
  200          XPARAM(I)=XSTORE(I)
               CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.)
               WRITE(6,'(//10X,'' GEOMETRY ON ONE SIDE OF THE TRANSITION
     1'','' STATE'')')
               CALL WRITMO(TIME0,FUNCT1)
            ENDIF
            TIME2=SECOND()
            WRITE(6,'('' TIME='',F9.2)')TIME2-TIME1
            TIME1=TIME2
            WRITE(6,'(''  REACTANTS AND PRODUCTS SWAPPED AROUND'')')
            IFLAG=1-IFLAG
            ONE=-1.D0
            EOLD=ESCF
            SUM=GOLD
            GOLD=GNORM
            I=1.7D0+ONE*0.5D0
            IF(GNORM.GT.10.D0)GOK(I)=.TRUE.
            GNORM=SUM
            DO 210 I=1,NATOMS
               DO 210 J=1,3
                  X=GEO(J,I)
                  GEO(J,I)=GEOA(J,I)
  210       GEOA(J,I)=X
            DO 220 I=1,NVAR
               X=XOLD(I)
               XOLD(I)=XPARAM(I)
  220       XPARAM(I)=X
C
C
C    SWAP AROUND THE DENSITY MATRICES.
C
            DO 230 I=1,LINEAR
               X=PASTOR(I)
               PASTOR(I)=PA(I)
               PA(I)=X
               X=PBSTOR(I)
               PBSTOR(I)=PB(I)
               PB(I)=X
               P(I)=PA(I)+PB(I)
  230       CONTINUE
            IF(FINISH) GOTO 250
         ELSE
            ONE=1.D0
         ENDIF
  240 CONTINUE
  250 CONTINUE
      WRITE(6,'('' AT END OF REACTION'')')
      GOLD=SQRT(DOT(GRAD,GRAD,NVAR))
      CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.)
      GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
      DO 260 I=1,NVAR
  260 GROLD(I)=XPARAM(I)
      CALL WRITMO(TIME0,FUNCT1)
*
* THE GEOMETRIES HAVE (A) BEEN OPTIMIZED CORRECTLY, OR
*                     (B) BOTH ENDED UP ON THE SAME SIDE OF THE T.S.
*
*  TRANSITION STATE LIES BETWEEN THE TWO GEOMETRIES
*
      C1=GOLD/(GOLD+GNORM)
      C2=1.D0-C1
      WRITE(6,'('' BEST ESTIMATE GEOMETRY OF THE TRANSITION STATE'')')
      WRITE(6,'(//10X,'' C1='',F8.3,''C2='',F8.3)')C1,C2
      DO 270 I=1,NVAR
  270 XPARAM(I)=C1*GROLD(I)+C2*XOLD(I)
      STEP=0.D0
      CALL COMPFG (XPARAM, .TRUE., FUNCT1,.TRUE.,GRAD,.TRUE.)
      CALL WRITMO(TIME0,FUNCT1)
  280 RETURN
      END
