Intrinsic Reaction Coordinate Methods
 
    The Intrinsic Reaction Coordinate (IRC) is defined as 
the minimum energy path connecting the reactants to 
products via the transition state.  In practice, the IRC is 
found by first locating the transition state for the 
reaction.  The IRC is then found in halves, going forward 
and backwards from the saddle point, down the steepest 
descent path in mass weighted Cartesian coordinates.  This 
is accomplished by numerical integration of the IRC 
equations, by a variety of methods to be described below.
 
    The IRC is becoming an important part of polyatomic 
dynamics research, as it is hoped that only knowledge of 
the PES in the vicinity of the IRC is needed for prediction 
of reaction rates, at least at threshhold energies.  The 
IRC has a number of uses for electronic structure purposes 
as well.  These include the proof that a certain transition 
structure does indeed connect a particular set of reactants 
and products, as the structure and imaginary frequency 
normal mode at the saddle point do not always unambiguously 
identify the reactants and products.  The study of the 
electronic and geometric structure along the IRC is also of 
interest.  For example, one can obtain localized orbitals 
along the path to determine when bonds break or form.
 
    The accuracy to which the IRC is determined is dictated 
by the use one intends for it.  Dynamical calculations 
require a very accurate determination of the path, as 
derivative information (second derivatives of the PES at 
various IRC points, and path curvature) is required later. 
Thus, a sophisticated integration method (such as GS2), and 
small step sizes (STRIDE=0.05, 0.01, or even smaller) may 
be needed.  In addition to this, care should be taken to 
locate the transition state carefully (perhaps decreasing 
OPTTOL by a factor of 10, to OPTTOL=1D-5), and in the 
initiation of the IRC run.  The latter might require a 
hessian matrix obtained by double differencing, certainly 
the hessian should be PROJCT'd or PURIFY'd.  Note also that 
EVIB must be chosen carefully, as decribed below.
 
    On the other hand, identification of reactants and 
products allows for much larger step sizes, and cruder 
integration methods.  In this type of IRC one might want to 
be careful in leaving the saddle point (perhaps STRIDE 
should be reduced to 0.10 or 0.05 for the first few steps 
away from the transition state), but once a few points have 
been taken, larger step sizes can be employed.  In general, 
the defaults in the $IRC group are set up for this latter, 
cruder quality IRC.  The STRIDE value for the GS2 method 
can usually be safely larger than for other methods, no 
matter what your interest in accuracy is. 

     The next few paragraphs describe the various 
integrators, but note that GS2 is superior to the others.

     The simplest method of determining an IRC is linear 
gradient following, PACE=LINEAR.  This method is also known 
as Euler's method.  If you are employing PACE=LINEAR, you 
can select "stabilization" of the reaction path by the 
Ishida, Morokuma, Komornicki method.  This type of 
corrector has no apparent mathematical basis, but works 
rather well since the bisector usually intersects the 
reaction path at right angles (for small step sizes).  The 
ELBOW variable allows for a method intermediate to LINEAR 
and stabilized LINEAR, in that the stabilization will be 
skipped if the gradients at the original IRC point, and at 
the result of a linear prediction step form an angle 
greater than ELBOW.  Set ELBOW=180 to always perform the 
stabilization.
 
     A closely related method is PACE=QUAD, which fits a 
quadratic polynomial to the gradient at the current and 
immediately previous IRC point to predict the next point. 
This pace has the same computational requirement as LINEAR, 
and is slightly more accurate due to the reuse of the old 
gradient.  However, stabilization is not possible for this 
pace, thus a stabilized LINEAR path is usually more 
accurate than QUAD.
 
    Two rather more sophisticated methods for integrating 
the IRC equations are the fourth order Adams-Moulton 
predictor-corrector (PACE=AMPC4) and fourth order Runge- 
Kutta (PACE=RK4).  AMPC4 takes a step towards the next IRC 
point (prediction), and based on the gradient found at this 
point (in the near vincinity of the next IRC point) obtains 
a modified step to the desired IRC point (correction). 
AMPC4 uses variable step sizes, based on the input STRIDE. 
RK4 takes several steps part way toward the next IRC point, 
and uses the gradient at these points to predict the next 
IRC point.  RK4 is one of the most accurate integration 
method implemented in GAMESS, and is also the most time 
consuming.

    The Gonzalez-Schlegel 2nd order method (PACE=GS2) finds 
the next IRC point by a constrained optimization on the 
surface of a hypersphere, centered at a point 1/2 STRIDE 
along the gradient vector leading from the previous IRC 
point.  By construction, the reaction path between two 
successive IRC points is a circle tangent to the two 
gradient vectors.  The algorithm is much more robust for 
large steps than the other methods, so it has been chosen 
as the default method.  Thus, the default for STRIDE is too 
large for the other methods.  The number of energy and 
gradients need to find the next point varies with the 
difficulty of the constrained  optimization, but is 
normally not very many points.  Taking more than 2-3 steps 
in this constrained optimization is indicative of reaction 
path curvature, and thus it may help to reduce the step 
size.  Use a small GCUT (same value as OPTTOL) when trying 
to integrate an IRC very accurately, to be sure the 
hypersphere optimizations are well converged.  Be sure to 
provide the updated hessian from the previous run when 
restarting PACE=GS2.

    The number of wavefunction evaluations, and energy 
gradients needed to jump from one point on the IRC to the 
next point are summarized in the following table:
 
     PACE      # energies   # gradients
     ----      ----------   -----------
    LINEAR        1             1
stabilized
    LINEAR        3             2
    QUAD          1             1  (+ reuse of historical
                                            gradient)
    AMPC4         2             2  (see note)
    RK4           4             4
    GS2          2-4           2-4 (equal numbers)
 
Note that the AMPC4 method sometimes does more than one 
correction step, with each such correction adding one more 
energy and gradient to the calculation.  You get what you 
pay for in IRC calculations:  the more energies and 
gradients which are used, the more accurate the path found.

    A description of these methods, as well as some others 
that were found to be not as good is geven by Kim Baldridge 
and Lisa Pederson, Pi Mu Epsilon J., 9, 513-521 (1993).
 
                         * * * 

    All methods are initiated by jumping from the saddle 
point, parallel to the normal mode (CMODE) which has an 
imaginary frequency.  The jump taken is designed to lower 
the energy by an amount EVIB.  The actual distance taken is 
thus a function of the imaginary frequency, as a smaller 
FREQ will produce a larger initial jump.  You can simply 
provide a $HESS group instead of CMODE and FREQ, which 
involves less typing.  To find out the actual step taken 
for a given EVIB, use EXETYP=CHECK.  The direction of the 
jump (towards reactants or products) is governed by FORWRD.  
Note that if you have decided to use small step sizes, you 
must employ a smaller EVIB to ensure a small first step.  
The GS2 method begins by following the normal mode by one 
half of STRIDE, and then performing a hypersphere 
minimization about that point, so EVIB is irrelevant to 
this PACE.
 
    The only method which proves that a properly converged 
IRC has been obtained is to regenerate the IRC with a 
smaller step size, and check that the IRC is unchanged. 
Again, note that the care with which an IRC must be 
obtained is highly dependent on what use it is intended 
for.
 
    Some key IRC references are:
K.Ishida, K.Morokuma, A.Komornicki
      J.Chem.Phys.  66, 2153-2156 (1977)
K.Muller
      Angew.Chem., Int.Ed.Engl. 19, 1-13 (1980)
M.W.Schmidt, M.S.Gordon, M.Dupuis
      J.Am.Chem.Soc.  107, 2585-2589 (1985)
B.C.Garrett, M.J.Redmon, R.Steckler, D.G.Truhlar,
K.K.Baldridge, D.Bartol, M.W.Schmidt, M.S.Gordon
      J.Phys.Chem.  92, 1476-1488(1988)
K.K.Baldridge, M.S.Gordon, R.Steckler, D.G.Truhlar
      J.Phys.Chem.  93, 5107-5119(1989)
C.Gonzalez, H.B.Schlegel
      J.Chem.Phys.  90, 2154-2161(1989)

    The IRC discussion closes with some practical tips:
 
    The $IRC group has a confusing array of variables, but 
fortunately very little thought need be given to most of 
them.  An IRC run is restarted by moving the coordinates of 
the next predicted IRC point into $DATA, and inserting the 
new $IRC group into your input file.  You must select the 
desired value for NPOINT.  Thus, only the first job which 
initiates the IRC requires much thought about $IRC.
 
    The symmetry specified in the $DATA deck should be the 
symmetry of the reaction path.  If a saddle point happens 
to have higher symmetry, use only the lower symmetry in the 
$DATA deck when initiating the IRC.  The reaction path will 
have a lower symmetry than the saddle point whenever the 
normal mode with imaginary frequency is not totally 
symmetric.  Be careful that the order and orientation of 
the atoms corresponds to that used in the run which 
generated the hessian matrix.
 
    If you wish to follow an IRC for a different isotope, 
use the $MASS group.  If you wish to follow the IRC in 
regular Cartesian coordinates, just enter unit masses for 
each atom.  Note that CMODE and FREQ are a function of the 
atomic masses, so either regenerate FREQ and CMODE, or more 
simply, provide the correct $HESS group.