```

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
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

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
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:

----      ----------   -----------
LINEAR        1             1
stabilized
LINEAR        3             2
QUAD          1             1  (+ reuse of historical
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.

```

created on 7/7/2017