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

Question

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)

Potentially Helpful Imagery

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?

Answer 1

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.

Answer 2

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.