Transition+State+Modeling

=Transition state modeling in quantum chemistry=

Tips and tricks learned from my stint at Proctor & Gamble -- NOTE: No data contained herein is confidential and is taken from tutorials from Spartan or from previously published materials.
While a huge fan of opensource software such as NWChem, Psi4, GAMESS, etc. we must to Spartan/Qchem submit for their easy-to-use interface in modeling.

So, we have a reaction, we know reactants, we know products. How are these molecules linked? That's what we want to know! We should find one, two, three, maybe many more possible structures that conceivably link the potential energy surface of the reactants with that of the products. If we were to model the potential energy surface, the structure(s) we hope to find fall within this space at points where we achieve a local minimum in energy in the x and y axes, while a maximum in the z. These are simply curvatures in the second-derivative of the energy where two curvatures are negative (minima), one positive (maximum). It looks like a potato chip. However, a true potential energy surface is many-dimensional and we can have many-order saddle points. If we were to compute the vibrational spectrum for a molecule at a saddle point, we'd end up with negative frequencies. The number of negative frequencies we find is equal to the order of the saddle point. Anything beyond exactly 1 is not a transition state candidate in any way. If you compute exactly one negative (and in more appropriate terms it's a complex number, so imaginary) frequency **AND** the forces are well-converged, meaning very close to zero, then you've found a transition state of some kind. It is important to realize that not every transition state you find will necessarily connect your reactants to your expected products. It could be a transition state to other chemistry. To check this, Spartan has a wonderful tool called "MakeList" which will take the imaginary frequency vibration and make a set of coordinates for some number of steps along the vibration. You then simply run a geometry optimization on the point to the left (let's call this the reactant side) and to the right (let's call this the product side) of the saddle point. What you should recover is the reactant and product connected by this transition state. The differences in the electronic energies of reactants and the transition state, corrected for the zero point energy (ZPE), neglecting thermal effects, is the activation barrier for the reaction. The barrier to head backwards simply involves the energy of the products then. The ZPE gives you the energy at 0K, typically there is an additional thermal correction to go from 0K to say 298K, but it is usually smaller and often neglected. Using the activation barrier, you can use experimental data and the Arrhenius equation to compute the kinetics! Experimental data will help in determining an appropriate constant to fit to your model.

So at P&G, what did I learn about modeling transition states?

1.) Pen and paper chemistry will get you pretty far, believe it or not. Organic chemistry isn't just to weed out the slackers, that stuff's useful! 2.) Transition states love to form 6-membered rings and just finding a way to form a ring that incorporates the reacting bits of the molecule can lead to a quick set of candidates. 3.) If the reaction involves a typical coordination, like SN2 does for 'backside attack' you should probably run a calculation where the nucleophile is placed in an axial position, along with the leaving group. 4.) If you have a number of non-covalently interacting species present, we found molecular dynamics to an amazing source of potential new chemistries to explore. Keep in mind that classical mechanics may get it wrong or that the structure you see is due to larger scale effects you may not recover in gas phase quantum mechanics on the smaller subset of atoms. 5.) QST2 and QST3 from Gaussian are useful as well. QST3 is more so than QST2 because you have a TS guess already. Effectively all you do is pass a rope between the structures and wiggle it around until you find the right path to link them. In QST2 you give reactants and products. In QST3 you give reactants, products, and a TS guess (in that order!). We have found this to be most successful when the reactants and products are reasonably close to the TS structure as well. If you found your TS by performing an energy profile, grab a structure from either side of your TS guess that's somewhere a little ways up the energy hill towards the TS. You can also perform an IRC which will attempt to follow the energy downhill on both sides to reactants and products. The "MakeList" method I mentioned earlier (without optimizing) is the 'cheater' method to get an IRC.

Okay, I think I have a candidate, how should I optimize it?

This is tricky. I have found just throwing caution to the wind can result in running calculations for needlessly long times, only to never have it converge. So, here's what I do in Spartan (you can do something similar in Gaussian).

Spartan technique...

 * 1) Transition state optimization at the lowest level you think you can get away with (I usually do HF/3-21G* or HF/6-31G(d))
 * 2) Set OPTCYCLE=100 to let it roll around in the potential energy surface for awhile
 * 3) If you're lucky, you're done. However, chances are you ran out of steps. The geometry likely changed a bit, maybe some bonds are longer, some are shorter. Maybe the approximate Hessian matrix used by default in Spartan isn't terribly effective here. If you go out to say 999 steps and you notice the calculation gets hung up and just starts oscillating between 2-3 structures -- you're Hessian has run its course. Time for the legitimate one. So for step 3, setup an energy calculation at the same level of theory you want to get a TS at. Click the button to compute the IR spectrum, which requires the Hessian matrix.
 * 4) Re-run the same calculation you did in step 2, except now Spartan will read the Hessian matrix from step 3! If your calculation was stuck bouncing between a couple structures, it's likely fixed! If not, you'll need to consider bumping up the convergence criteria or using the DIIS algorithm. Typically, I also set my force and displacement criteria to what is equivalent to opt=tight in Gaussian. Keywords are: GRADIENTTOLERANCE=0.000001 DISTANCETOLERANCE=0.000004 SCFTOLERANCE=HIGH. Once major difference between Spartan/Qchem and Gaussian -- only the forces and one of either the energy or displacements must converge for the optimization to succeed. In Gaussian, all three must meet their criteria. There is no way I know of to change this, other than to re-run in Gaussian.

Gaussian technique...

 * 1) Transition state optimization at the lowest level you think you can get away with (I usually do HF/3-21G* or HF/6-31G(d))
 * 2) opt=(ts,tight,calcfc) --> this is the equivalent of running the energy calculation we discussed for Spartan, then using the resulting Hessian. The force constants are computed only at the initial step. I usually do 50-100 cycles, then re-run the calculation with the new coordinates (you can use geom=allcheck). Got a truly difficult case? Consider calcall instead of calcfc which recomputes force constants at every step. This is eye-bleedingly slow in explicitly correlated methods (i.e., MP2, CC, CI). Use as a last resort.
 * 3) Just repeat until it converges. May take several cycles. Save some time and leave the freq keyword off so it doesn't do the whole thermo-shebang for a poor structure. Add it once you get pretty close or once you get a good geometry.

MORE TO COME!

Author: TPP 05/03/2013