# Reaction Pathways: What to do when IRC geometries are relaxed (optimized) and the energy changes significantly

The Story:

Reaction pathways can be difficult things to deal with using electronic structure methods. Intrinsic reaction coordinates (IRCs) can be determined a number of ways (a discussion of these methods lie way outside the scope of this post). Suffice to say, the reaction pathway is an 'energy' pathway of reaction, connecting reactants and products, usually through some intermediate and/or transition state (TS).

For the purposes of this question, imagine implementing some sort of IRC technique which begins with a previously resolved transition state (Figure 1). The computation follows the reaction coordinate in both the forward and reverse directions from the TS energetically downhill until it reaches a (for the sake of simplicity) a minimum energy point with features which satisfy the method used. In general, the structure and energy of this final IRC geometry should be very representative of the structure which is connected to the TS via the pathway determined. If you were to optimize this final geometry, you should get good agreement between the optimized structure and the IRC configuration.

Now consider this scenario. An optimization was performed on the IRC resolved structure B and the energy relaxation was significant and the structure significantly changed, giving rise to a new (and lower) energy. One could no longer (within the formalism of IRC methods) say that the optimized structure B lies on the reaction pathway, and therefore is not a representative configuration at point B on the reaction pathway associated with the TS.

The Question:

Does anyone know of an example in the literature where this type of behavior occurs? How does one 'deal' with a situation like this?

Other Information:

Let us assume the following:

1. That this scenario occurs with a few IRC procedures (not just one).

2. That the IRC computations were given a ridiculously large amount of cycles and steps

3. That various stepsizes were sampled in parallel IRC jobs with no significant difference in the results

(I will add more as I think of them)

Imagine two somewhat complex molecules coming together to form a single covalent bond. You run an IRC from the TS (which describes the bond formation) in the reverse direction and the intermolecular separation begins to increase. This separation continues to increase until the IRC job terminates properly. Now imagine running an optimization on these two (slightly complex) molecules which are separated by maybe 4.0 Å. Obviously these molecules are going to interact because they are not infinitely separated. You observe the (slightly complex) molecules moving toward each other during the optimization procedure (rotating, twisting, gyrating along the way) until you now have a minimum energy structure which is much lower in energy than the resulting IRC computation. These two molecules interact via intermolecular interactions (no covalent bonds are formed between the molecules). This optimized structure is analogous to a reactant complex (as opposed to separated reactants) and has an energy a few kcal mol$^{-1}$ lower than the corresponding IRC geometry. How on earth does this fit in to the reaction pathway?

• I am not entirely familiar with the techniques that you are using but my hunch is that this is a matter of kinetics =/= thermodynamics. When running this IRC simulation, are you forcing the reaction to take place? In your last paragraph it sounds like you're taking a TS, separating the two molecules and then doing an energy minimization. The resulting complex now has a lower energy state than the covalently bonded version of the two molecules. Is this correct? – CTKlein Apr 7 '14 at 19:34
• The cartoon I've included is not quantitative (my apologies). An IRC simply follows the imaginary mode of vibration in the TS downhill along the corresponding coordinate (its not forced). The resulting reactant complex may or may not be lower in energy than the covalently bound molecule (the products). – LordStryker Apr 7 '14 at 19:52
• I guess this is a known problem with Angstroms (tex.stackexchange.com/questions/24255/angstrom-not-working), so I pasted it in. Stellar question! – jonsca Apr 7 '14 at 21:09
• @LordStryker but the chart does qualitatively show that the associated complex has a lower free energy as compared to the covalently bound complex? I could certainly see that being a plausible result... – CTKlein Apr 7 '14 at 23:17
• I have not had the time to read it, but it might help you anyway: mathematik.uni-leipzig.de/~quapp/JTCCpreprint.pdf DOI: 10.1142/S0219633603000604 – Martin - マーチン Apr 22 '14 at 3:25

IRC is a gradient following algorithm that tries to follow the imaginary mode. Check the algorithm (several one implemented, so really depends on which you are talking about) and you will see that unlike in optimisation algorithms, there is a big emphasis on following the real gradient (and not just jump to the endpoint as fast as possible).

It works well around the TS where the gradient is large and has a well define direction. However getting closer to the optimised geometry, the requirements of following the exact gradient and having an effective optimization algorithm are contradicting. In analytical sense it may be not obvious, but when you have numeric gradients, small second derivatives (these are almost always very soft modes), the direction of the gradient may have big errors and you should keep a small stepsize to check it thoroughly.

In summary, the situation you wrote down is more or less common at some extent in almost all the IRC calculations due to the numerical nature of the problem. Harder imaginary frequencies may help.

The IRC algorithm usually follows the gradient from TS to product state B and by definition has to go always downhill in energy. In best case, it treats the Potential Energy Surface (PES) around the given point as quadratic and just integrates Newtonian equation of motion and does it until it reaches local minimum.

This procedure is different from geometry optimization, which can use powerful methods to find energy minimum.

Static picture, T = 0K shows that the reaction will end in the "upper" B state. If you then ask for energy minimum, reorientation of any kind (intermolecular or conformational) takes place and you find "lower" B. Do not be decepted by your plot, the optimized B lies on different point of reaction coordinate (however you define it) than the IRC-type B.

Dynamic picture, T > 0K shows that IRC is not the only way to proceed from TS to product, although it is still important representative. Upon running molecular dynamics (MD) from TS, you will end up in the equilibrium state, probably well represented by "lower" B.

It would be quite difficult to find exactly this type of behavior in literature, because such technical topics are quite often swept under the carpet and only indirect hints remain. Similar topics are discussed in Levine: Molecular Reaction Dynamics. As a simple example take the E2 elimination from cyclohexane, where the leaving groups must be axial, so the reactive species is different from energy minimum.

And how to deal with it? It is up to you, where you set zero to your relative energy scale. It will depend on which energies you want to compare, if you propose some series of reactions.

• It would be quite difficult to find exactly this type of behavior in literature, because such technical topics are quite often swept under the carpet and only indirect hints remain. This is what I was afraid of. Wordsmithing is as much an art as it is a skill when it comes to publishing. – LordStryker Apr 8 '14 at 11:58
• On the other hand, it really IS a just a technical question, where to set zero and how to treat the pre-reaction complex. As an illustration how flexibly is the IRC used you can look at the X-axis labeling in the papers. This ranges from no labels to weird labels, while only occasionally you can see some units and values. – ssavec Apr 8 '14 at 16:15
• Well if someone is proposing a reaction mechanism and an optimized geometry is drastically different from the one on the IRC, how does one say that 'this structure' is representative of the one on the actual pathway? Doesn't this sort of 'sink' the whole "We've mapped out the reaction pathway..." conclusion? – LordStryker Apr 8 '14 at 16:30
• It is possible to map the PES around the reactant state and find transition states for the conformational (or similar) change. You often find, that the energies involved in this process are order of magnitude smaller than for the main reaction and treat it using steady-state approximation. – ssavec Apr 8 '14 at 16:42
• @LordStryker "Although the reaction coordinate' does not correspond to any actual trajectory on a given potential energy surface, it gives a path about which many actual trajectories, especially those of low energy, are likely to travel." IRC should be used to map/ approximate the immediate surrounding of a TS, therefore ensuring that your geometry optimization will end up in the most likely/ closest local minimum. And that is also enough, because all of the local minima/ TS are on the same PES. – Martin - マーチン Apr 9 '14 at 3:17