Get a reaction: Intrinsic Reaction Coordinate

Lucas requested a post on computing an intrinsic reaction coordinate (IRC). An IRC is a special case of the minimum energy path (MEP) connecting reactant, transition state, and product (more about this at the end of the post).

There are many uses for the IRC, the most basic of which is to verify that the transition state you found actually connects the reactants and products you think it does (this is not always obvious from the normal mode associated with the imaginary frequency). It is also a good way to identify any reaction intermediates you haven't thought about. The IRC is also a way generate an animation of the reaction (see for example ChemTube3D).

An IRC is generated by following the gradient down-hill from the transition state to products and reactants. (It's a bit like a geometry optimization initiated at the transition state, but the algorithm used to follow the gradient is different). Computing a complete IRC thus requires a minimum of two separate runs, and you must find the TS before you can compute an IRC.

The screencast shows how to generate the IRC associated with the 3TSa transition state I found in a previous post. The first input file I generated contain the 3TSa coordinates and the following keywords:

$contrl runtyp=irc $end
$basis gbasis=pm3 $end
$irc saddle=.t. tsengy=.t. forwrd=.t. npoint=100 stride=0.2 opttol=0.0005 $end

saddle=.t. tells GAMESS that the coordinates in $data corresponds to a TS. Because the TS is a stationary point the gradient is zero (or very small), so there is no gradient to follow down-hill. Therefore, the IRC is initiated by displacing the geometry along the normal mode corresponding to the imaginary frequency. Therefore, saddle=.t. GAMESS expects a $HESS group, which contains the frequency information.

Whether the normal mode points towards the reactants or products is completely arbitrary, and the forwrd=.f. keyword is used to change the direction in which the normal mode points. forwrd=.t. is the default and I only include it in the first run for completeness.

tsengy=.t. tells GAMESS to re-compute the energy of the TS (re-compute because you already computed it when you computed the frequencies). If you want to make an energy plot in MacMolPlt as I show in the screen cast, you must set tsengy=.t.

The displaced geometry will have a non-zero gradient, the negative of which points downhill. The geometry is displaced in this direction, the gradient is recomputed, and the procedure is repeated npoint=100 times or until the (root-mean-square) of the gradient is below opttol=0.0005 atomic units, at which point you are at a minimum. (Notice that the default for npoint is 1, and must always be changed).

stride=0.2 tells GAMESS how far along the gradient it should displace at each step (in this case 0.2 atomic units). The default is 0.3, but when I tried this, the IRC calculation in the forward direction failed after 7 steps, because the coordinate deviated too much from the minimum energy path, so I had to decrease it. This means I will have to take more steps to reach the minimum.

As you saw in the screencast, GAMESS takes 100 steps in each direction, which means the IRC does not go all the way to the minimum on each side of the TS. GAMESS does print out restart information to continue the IRC, which I didn't do in this case.

Finally a note on MacMolPlt. As I show in the screencast, MacMolPlt can be used to combine the two IRC runs (forwrd=.t. and forwrd=.f.) by using Files > Add Frames from File... You should first load the file that goes to product and then add the file from the IRC that goes to reactants. Here you should tell MacMolPlt to reverse the order of the structures in the file, so the animation goes from reactions to transition state to product.

This is done by selecting "make these points negative". "Negative" refers to the x-coordinate in the energy plot. The x-coordinate is what separates an IRC from a minimum energy path. An IRC is a plot of the energy of each structure along the minimum energy path vs the root-mean-square change in the mass weighted Cartesian coordinates relative to the TS structure. The TS thus has a x-coordinate of zero, while structures leading to the reactants are defined by negative RMSD values.

If the energy is plotted against something else, such as a bond-length, then it is simply a minimum energy path. Because an IRC uses mass-weighted coordinates, using different isotopes will lead to different IRCs, and one use of IRC is to determine tunneling corrections to the reaction rate.