      SUBROUTINE REACT1(ESCF)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      common /GEOM  / GEO(3,NUMATM)
      dimension GEOA(3,NUMATM), GEOVEC(3,NUMATM),
     1          P1STOR(MPACK), P2STOR(MPACK),
     2          P3STOR(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),
     1          IDUM3(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(MIOA,LABELS,GEOA,LOC,NA,NB,NC,AMS,NATOMS,INT)
      else
         call GETGEO(MIOA,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(MIOB,'(10X,''ERRORS DETECTED IN CONNECTIVITY'')')
               write(MIOB,'(A,I3,A,I3,A,I3,A)')
     %           ' FOR ATOM',I,
     %           ' THE BOND LABELS ARE DIFFERENT:      ',
     %            IDUM1(I),' AND',NA(I)
            end if
            if (IDUM2(I).NE.NB(I)) then
               L=L+1
               if (L.EQ.1)
     %           write(MIOB,'(10X,''ERRORS DETECTED IN CONNECTIVITY'')')
               write(MIOB,'(A,I3,A,I3,A,I3,A)')
     %           ' FOR ATOM',I,
     %           ' THE BOND ANGLE LABELS ARE DIFFERENT:',
     %            IDUM2(I),' AND',NB(I)
            end if
            if (IDUM3(I).NE.NC(I)) then
               L=L+1
               if (L.EQ.1)
     %           write(MIOB,'(10X,''ERRORS DETECTED IN CONNECTIVITY'')')
               write(MIOB,'(A,I3,A,I3,A,I3,A)')
     %           ' FOR ATOM',I,
     %           ' THE DIHEDRAL LABELS ARE DIFFERENT:  ',
     %            IDUM3(I),' AND',NC(I)
            end if
   10    continue
         if (L.NE.0) then
           write(MIOB,'(10X,A)')' CORRECT BEFORE RESUBMITTING'
           call filclo  ! J/M
           stop ' in REACT1 - connectivity problems'
         end if
      end if
      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(MIOB,'(//10X,'' NUMBER OF ATOMS IN SECOND SYSTEM IS '',
     %                      ''INCORRECT'',/)')
         write(MIOB,'('' NUMBER OF ATOMS IN FIRST  SYSTEM ='',I4)')NUMAT
         write(MIOB,'('' NUMBER OF ATOMS IN SECOND SYSTEM ='',I4)')
     %                                                            NUMAT2
         go to 280
      end if
      write(MIOB,'(//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(MIOB,'(//,''  CARTESIAN GEOMETRY OF FIRST SYSTEM'',//)')
         write(MIOB,'(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
                        go to 100
                     end if
                     SA=-SA
                  end if
   90          SUM=SUMM
  100       continue
  110    continue
         write(MIOB,'(//,''  CARTESIAN GEOMETRY OF SECOND SYSTEM'',//)')
         write(MIOB,'(3F14.5)') ((GEOA(J,I),J=1,3),I=1,NUMAT)
         write(MIOB,'(//,''   "DISTANCE":'',F13.6)') SUM
         write(MIOB,'(//,''  REACTION COORDINATE VECTOR'',//)')
         write(MIOB,'(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)
            end if
  130    continue
         NATOMS=NUMAT
      end if
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(MIOB,'(///10X,''THERE ARE NO VARIABLES IN THE SADDLE'',
     %                       '' CALCULATION!'')')
         call filclo  !  J/M
         stop ' in REACT1 - no variables in the saddle calculation'
      end if
      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(MIOB,'(//,3(5X,A/))')  ' BOTH GEOMETRIES ARE IDENTICAL',
     %       ' A SADDLE CALCULATION INVOLVES A REACTANT AND A PRODUCT',
     2       ' THESE MUST BE DIFFERENT GEOMETRIES'
         call filclo  ! J/M
         stop ' in REACT1 - identical substrate & product geometries'
      end if
      ONE=1.D0
      DELL=0.1D0
      EOLD=-2000.D0
      TIME1=SECOND()
      SWAP=0
      do 240 ILOOP=1,MAXSTP
         write(MIOB,'('' '',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.1
         write(MIOB,'(''  BAR SHORTENED BY'',F12.7,'' PERCENT'')')
     %                                               STEP/STEP0*100.D0
         STEP0=STEP0-STEP
         if (STEP0.LT.0.01D0) go to 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
               P1STOR(I)=P(I)
               P2STOR(I)=PA(I)
  160       P3STOR(I)=PB(I)
         end if
         do 170 I=1,NVAR
  170    XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
         if (IFLAG.EQ.1) then
            write(MIOB,
     %        '(//10X,''FOR POINT'',I3,'' SECOND STRUCTURE'')') ILOOP
         else
            write(MIOB,
     %        '(//10X,''FOR POINT'',I3,'' FIRST  STRUCTURE'')') ILOOP
         end if
         write(MIOB,'('' DISTANCE A - B  '',F12.6)')            STEP
         write(*,'('' 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(MIOB,'(''  ACTUAL GRADIENTS OF THIS POINT'')')
            write(MIOB,'(8F10.4)')  (GRAD(I),I=1,NVAR)
         end if
         write(MIOB,'('' HEAT            '',F12.6)')  FUNCT1
         GNORM=SQRT(DOT(GRAD,GRAD,NVAR))
         write(MIOB,'('' GRADIENT NORM   '',F12.6)')  GNORM
         COSINE=COSINE*ONE
         write(MIOB,'('' DIRECTION COSINE'',F12.6)')  COSINE
         call GEOUT(MIOB)
         if (   SWAP.GT.2.9D0 .OR. ILOOP .GT. 3 .AND.
     1      COSINE .LT. 0.D0 .OR. ESCF .GT. EOLD)      then
            if (SWAP.GT.2.9D0) then
               SWAP=0.D0
            else
               SWAP=0.5D0
            end if
C
C   SWAP REACTANT AND PRODUCT AROUND
C
            FINISH=(GOK(1).AND.GOK(2) .AND. COSINE .LT. 0.D0)
            if (FINISH) then
               write(MIOB,
     %           '(//10X,'' BOTH SYSTEMS ARE ON THE SAME SIDE OF THE '',
     %                   ''TRANSITION STATE -''/10X,
     %                   '' GEOMETRIES OF THE SYSTEMS'',
     %                   '' 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(MIOB,
     %            '(//10X,'' GEOMETRY ON ONE SIDE OF THE TRANSITION'',
     %                    '' STATE'')')
               call WRITMO(TIME0,FUNCT1)
            end if
            TIME2=SECOND()
            write(MIOB,'('' TIME='',F9.2)')  TIME2-TIME1
            TIME1=TIME2
            write(MIOB,'(''  REACTANTS AND PRODUCTS SWAPPED AROUND'')')
            IFLAG=1-IFLAG
            ONE=-1.D0
            EOLD=ESCF
            SUM=GOLD
            GOLD=GNORM
            I=1.7+ONE*0.5
            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=P1STOR(I)
               P1STOR(I)=P(I)
               P(I)=X
               X=P2STOR(I)
               P2STOR(I)=PA(I)
               PA(I)=X
               X=P3STOR(I)
               P3STOR(I)=PB(I)
               PB(I)=X
  230       continue
            if (FINISH) go to 250
         else
            ONE=1.D0
         end if
  240 continue
  250 continue
      write(MIOB,'('' 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(MIOB,
     %  '('' BEST ESTIMATE GEOMETRY OF THE TRANSITION STATE'')')
      write(MIOB,'(//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
