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.