9.5 Intrinsic Reaction Coordinates
The concept of a reaction path, although seemingly well-defined chemically (i.e., how the atoms in the system move to get from reactants to products), is somewhat ambiguous mathematically because, using the usual definitions, it depends on the coordinate system. Stationary points on a potential energy surface are independent of coordinates, but the path connecting them is not, and so different coordinate systems will produce different reaction paths. There are even different definitions of what constitutes a “reaction path”; the one used in Q-Chem is based on the intrinsic reaction coordinate (IRC), first defined in this context by Fukui [423]. This is essentially a series of steepest descent paths going downhill from the transition state.
The reaction path is most unlikely to be a straight line and so by taking a finite step length along the direction of the gradient you will leave the “true” path. A series of small steepest descent steps will zig-zag along the actual reaction path (this is known as “stitching”). Ishida et al. [424] developed a predictor-corrector algorithm, involving a second gradient calculation after the initial steepest descent step, followed by a line search along the gradient bisector to get back on the path; this was subsequently improved by Schmidt et al. [425], and is the method we have adopted. For the first step downhill from the transition state this approach cannot be used (as the gradient is zero); instead a step is taken along the Hessian mode corresponding to the imaginary frequency.
The reaction path can be defined and followed in Z-matrix coordinates, Cartesian coordinates or mass-weighted Cartesians. The latter represents the “true” IRC as defined by Fukui [423]. However, if the main reason for following the reaction path is simply to determine which minima a given transition state connects (perhaps the major use), then it doesn’t matter which coordinates are used. In order to use the IRC code the transition state geometry and the exact Hessian must be available. These must be computed via transition state (JOBTYPE = TS) and frequency calculation (JOBTYPE = FREQ) respectively.
9.5.1 Job Control
An IRC calculation is invoked by setting the JOBTYPE $rem to RPATH.
RPATH_COORDS
Determines which coordinate system to use in the IRC search.
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
Use mass-weighted coordinates.
1
Use Cartesian coordinates.
2
Use Z-matrix coordinates.
RECOMMENDATION:
Use default.
RPATH_DIRECTION
Determines the direction of the eigen mode to follow. This will not usually be known prior to the Hessian diagonalization.
TYPE:
INTEGER
DEFAULT:
1
OPTIONS:
1
Descend in the positive direction of the eigen mode.
-1
Descend in the negative direction of the eigen mode.
RECOMMENDATION:
It is usually not possible to determine in which direction to go a priori, and therefore both directions will need to be considered.
RPATH_MAX_CYCLES
Specifies the maximum number of points to find on the reaction path.
TYPE:
INTEGER
DEFAULT:
20
OPTIONS:
User-defined number of cycles.
RECOMMENDATION:
Use more points if the minimum is desired, but not reached using the default.
RPATH_MAX_STEPSIZE
Specifies the maximum step size to be taken (in thousandths of a.u.).
TYPE:
INTEGER
DEFAULT:
150
corresponding to a step size of 0.15 a.u..
OPTIONS:
Step size =
RECOMMENDATION:
None.
RPATH_TOL_DISPLACEMENT
Specifies the convergence threshold for the step. If a step size is chosen by the algorithm that is smaller than this, the path is deemed to have reached the minimum.
TYPE:
INTEGER
DEFAULT:
5000
Corresponding to 0.005 a.u.
OPTIONS:
User-defined. Tolerance =
RECOMMENDATION:
Use default. Note that this option only controls the threshold for ending the RPATH job and does nothing to the intermediate steps of the calculation. A smaller value will provide reaction paths that end closer to the true minimum. Use of smaller values without adjusting RPATH_MAX_STEPSIZE, however, can lead to oscillations about the minimum.
RPATH_PRINT
Specifies the print output level.
TYPE:
INTEGER
DEFAULT:
2
OPTIONS:
RECOMMENDATION:
Use default, little additional information is printed at higher levels. Most of the output arises from the multiple single point calculations that are performed along the reaction pathway.
9.5.2 Example
Example 9.197
$molecule
0 1
C
H 1 1.20191
N 1 1.22178 2 72.76337
$end
$rem
JOBTYPE freq
BASIS sto-3g
METHOD hf
$end
@@@
$molecule
read
$end
$rem
JOBTYPE rpath
BASIS sto-3g
METHOD hf
SCF_GUESS read
RPATH_MAX_CYCLES 30
$end