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. 

Is your reaction unimolecular?
$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.



Is your reaction bimolecular?
$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. 

Sunday, January 8, 2012

The ideal app for an ideal gas


I recently came across a very cool app called Atoms in Motion ($2.99).  With this app you can perform molecular dynamics simulations of He, Ne, Ar, and Kr (and mixtures thereof) at different temperatures, volumes, and number of atoms.


It is effortless to set up a simulation and change the conditions: swipe to change T (or shake the iPad!) and pinch to change V.  You can also flick an atom and, frankly, the most fun is to build a few clusters at 1 K and then move an atom carefully in position and flick it at a cluster! Something I have dubbed "atomic billiards".


Teaching
Since not every student has an iPad (yet), the main teaching-use will probably be in lecture, projection the simulation on using a VGA adapter (hint: simulations make for great peer instruction questions).  Because it is based on real simulations, and not pre-recorded movies, the app is very versatile and can be used to demonstrate very simple things like solids, liquids, and gases or an atomic view of temperature or relatively complex things like defects in solids and deviations from ideal gas behavior (using the plotting tool). 

However, the ideal use of this app would really be self-guided discovery.  Most features are very intuitive (the plotting tool takes some practice) and I'd love to see what a class of 10-year olds iPad users would get out of this. Side note: it would really have been great if the iPad could vibrate and the pressure could be linked to the vibration.


Documentation
Another strength of this app is the quality of documentation and extra information available in the app and on the top bar of the website.  I'm particularly impressed by the animated figure on "how big is an atom", which reminded me of the orders-of-ten approach.  I would be lovely to see a similar explanation of the nano-second.

Disclaimer: I was made aware of this app by the developer, but I bought my own copy, and was not asked to write a review on the blog.