## Sunday, January 22, 2012

### Negative activation energies in molecular modeling: diagnosis and cures

This post is inspired by a recent discussion with Paolo in the comments section of this post.

So, you've found that the energy of your transition state (TS) is lower than your reactants (i.e. you have a negative activation energy). Don't panic, instead pretend you're Dr House and go through it methodically:

Did you actually find the right TS?
A TS is a completely optimized geometry (no constraints, "zero" gradient) with one has one and only one imaginary frequency. Yes?  OK, how big is the imaginary frequency and does the normal mode look like what you would expect? Sometimes minima can have small (< ca 100 cm-1) imaginary frequencies due to numerical "noise", or the TS search algorithm will find a TS for, say, methyl rotation. No, looks OK?   (Consider computing an IRC do be really sure, just in case.)  Let's move on.

How much lower?
Some processes have very low (< 1 kcal/mol) barriers.  Numerical "noise" can cause these barriers to come out slightly negative (-0.1 to -1.0 kcal/mol) instead. If so, consider your reaction barrier-less for all intends and purposes. No, much more negative?  Let's move on.

What energy are we talking about here?
If the activation energy is calculated using single point energies (for example B3LYP/6-31G(d)//RHF/3-21G) then the B3LYP/6-31G(d) PES may not have a barrier or the B3LYP/6-31G(d) TS geometry looks very different from the RHF/3-21G TS geometry.

The easiest way to check for this particular case is to geometry optimize your reactant at the B3LYP/6-31G(d) level of theory.  If the optimization results in the product geometry, then there is (very likely) no barrier to the reaction.  If the optimization gives you a reactant geometry, then one explanation is that 3-21G is not a reliable method for finding the TS and the solution is to use B3LYP/6-31G(d) in your TS search. There are other explanations, so read on.

If the activation energy is calculated using zero point energy or free energy corrections and is negative (this is not uncommon for reactions whose PES barriers are less than 2-3 kcal/mol) then you can consider the reaction barrier-less.  If you want to be really sure, you can add these corrections along the entire IRC.  However, there are also other explanations so read on.

$R \rightarrow TS \rightarrow P$
OK, so you have found a reasonable-looking TS using some method and the electronic energy barrier computed using that method is quite negative.  If your reaction is unimolecular then the most likely explanation is that you have not found the lowest energy conformation of your reactant.  If can't find a lower energy structure for your reaction, compute an IRC.  An IRC follows the minimum energy path downhill to your reactant and product, so you will find structures with a lower energy than your TS.

$R_1 + R_2 \rightarrow R_1/R_2 \rightarrow TS \rightarrow P$

OK, so you have found a reasonable-looking TS using some method and the electronic energy barrier computed using that method is quite negative.  The first question is how did you compute the barrier?

If you computed the barrier as the energy difference between the TS and the bound complex of your two reactants (R1/R2) then the most likely explanation is that you have not found the lowest energy conformation of this complex (see the section on unimolecular reactions).

If you computed the barrier as the energy difference between the TS and the sum of the energies and your two isolated reactants [E(R1)+E(R2)] then there is not necessarily anything wrong with your calculations.  Read on.

How should I compute the activation energy of a bimolecular reaction?
(This sub-section has been updated 2012.05.05)
The activation energy for a biomolecular reaction should be computed as

$E_a = E(TS) - [E(R_1) + E(R_2)]$

If this energy is negative then you would say that this is a barrier-less reaction and the reaction rate is proportional to the collision frequency.