# Correct jobs sequence for transition-state (TS) search

I've seen many papers analyzing reaction pathway using a small basis set like 6-31G* for all optimization steps and later single point calculations in bigger one (like 6-311++G**). What is the correct sequence of jobs for TS search in such a case? Should it be:

1. Frequency (6-31G*)
2. TS_search (6-31G*)
3. Single Point Energy (6-311++G**)
4. Frequency (6-311++G**)

or

1. Frequency (6-31G*)
2. TS_search (6-31G*)
3. Frequency (6-31G*)
4. Single point Energy (6-311++G**)
5. Frequency (6-311++G**)

?

J1 Freq (6-31G*); J2 TS_search (6-31G*); J3 Freq (6-31G*);

is good because it let's you optimize the saddle point and then get the thermodynamic data

For the single-point correction don't recalculate the frequencies because the thermodynamic corrections won't change much anyway. You can also include implicit solvent in J4.

A common and not quite unreasonable approach is to do molecular structure optimisations at a low level, establish the path(s) at this level, and obtain better energies at a higher level. It really depends on what you need and what you can afford. However, there are plenty of things to consider, and doing this approach without checking may lead to wrong conclusions.

1. Freq (6-31G*)
2. TS_search (6-31G*)
3. SP (6-311++G**)
4. Freq (6-311++G**)

$$\strut{}$$

1. Freq (6-31G*)
2. TS_search (6-31G*)
3. Freq (6-31G*)
4. SP (6-311++G**)
5. Freq (6-311++G**)?

Both approaches are basically the same, at least for the transition state, because 2 will include frequencies, otherwise you will not find the TS.

Therefore 1, 2, 2 + Freq. is fine; 4 (top) or 5 (bottom) is wrong and will almost certainly lead to less accurate energies.
I also generally consider Pople basis sets outdated and underperforming, and especially inconsistent for any interpolation. Use something more reasonable, there are plenty of options. For DFT I recommend the def2-XVZP family as it is available in most QC packages. See the incredible https://www.basissetexchange.org/ for an extensive list and file formats.

Back to the job 5 problem. With a frequency calculation you are almost always employing the harmonic approximation. Therefore the obtained modes are only meaningful, if they are computed as the same level of theory as the optimised structure. A valid approach, however, is to recompute energies at a higher level and add thermal corrections from the optimisation at a lower level.
For example, in generic terms Density Functional Approximation (DFA), Split Valence Basis Set (SVBS), X-tuple Zeta Basis Set (X > 2, XZBS):

\small \begin{array}{llcc} \text{Optimisation:} & \text{DFA/SVBS} & \leadsto & E_\mathrm{el, SV}\\ \text{Anal. 2nd deriv.:} & \text{DFA/SVBS} & \leadsto & E_\mathrm{ZPE, SV}, \dots, H_{\mathrm{corr., SV}, T}, G_{\mathrm{corr., SV}, T}\\\hline \mathrm{Sum_{SV}} & & & E_\mathrm{o, SV}, \dots, H_{\mathrm{SV}, T}, G_{\mathrm{SV}, T}\\[2ex] \text{Single Point} & \text{DFA/3ZBS//SVBS} & \leadsto & E_\mathrm{el, 3Z}\\\hline \text{Final Energies} & \text{as DFA/3ZBS//SVBS} & & \begin{align} E_\mathrm{o, 3Z/SV} &= E_\mathrm{el, 3Z} + E_\mathrm{ZPE, SV}\\ \vdots\\ H_{\mathrm{3Z/SV}, T} &= E_\mathrm{el, 3Z} + H_{\mathrm{corr., SV}, T}\\ G_{\mathrm{3Z/SV}, T} &= E_\mathrm{el, 3Z} + G_{\mathrm{corr., SV}, T} \end{align}\\\hline\hline \end{array}

Some more practical tipps:

There are also a couple of modern tricks that may help you avoid blind guessing, and also avoid expensive calculations. The escalation scheme of above is probably already more than 20 years old. Thankfully not only our computers, but also our methodology has greatly improved. State of the art DFA from the last millenia can now run in a tiny fraction of the time. A consequence of that is treating larger systems is becoming increasingly popular. Unfortunately, the art of doing calculations efficiently has suffered from this. Here are a few tips, I hope you will consider for the future:

• Approximate wherever you can to generate starting structures

• Minimise conformational space, replace large moieties with methyl groups when possible. Especially isopropyl groups can be troublesome, or cyclohexyl groups and derivatives.

• Use cheap methodology. Semi-empirics has made a huge comeback, and it is fast and quite reliable nowadays. I cannot recommend xTB enough, see GitHub.

• Check the conformational space at that level. You don't want to be stuck with a guess, which is a high energy structure and not realise it.

• Use cheap pure DFA as workhorses (workhorse functional = WhF). Popular examples are BP86, PBE, M06L (or newer), TPSS. If your QC package can handle it, use density fitting (DF) or similarly the Resolution of the Identity approximation (RI).

For DFA, always use a large grid; it will cost more, but lead to faster convergence. Similarly, you always want to emplay tight criteria for the density. However, you can play with loose optimisation criteria. A loosely converged structure is always a better guess than a constructed structure.

Use these structures and validate/calibrate them against other functionals; again some popular examples are PBE0, M06 (or newer), TPSSh, B3LYP (and other derivatives). You don't want to have done all the work just to realise the one functional you have picked as a workhorse, produces significantly other results than the others.

• Optimise key intermediates first. This will help you find the transition states, as you have reasonable approximations for bond lengths etc.
Gaussian has a method called QST2/3, which will try to automatically compute a TS between two minima. Sure, the calculations are a bit more extensive, but you might find the TS at the first try, and not have blindly guess. You can or even should do that with semi-empirics, too.
In many other packages (obviously not Gaussian) you have nudged elastic band methods. Use them with semi-empirics. For example Orca (tutorial from the input library) can be coupled with xtb and provide a very reasonable guess for a TS within the hour.

• At this stage you'll probably have $$\ce{A -> [AB]^\dagger -> B}$$, etc. optimised at your WhF/SVBS level of theory. Validate! Then run higher level calculations with larger basis sets. Make sure your paths are well connected and make sense; run intrinsic reaction coordinate analyses, or other displacement algorithms.

• Validate, again. Yes I know this sounds like a broken record, but you want to justify your results. Obviously, your mileage may vary, so tailor this to your applications. You may want to reoptimise the key steps of your path, check for deviations in structure and whether the energies hold. If you calculate catalytic cycles, I recommend reading about the energy span model.[1]

• If you can afford it (and want to make double-dead sure), also validate against wave-function based methods, e.g. MP2, DLPNO-CC, CCSD-F12, etc..

1. Kozuch, S.; Shaik, S. How to Conceptualize Catalytic Cycles? The Energetic Span Model. Acc. Chem. Res. 2011, 44 (2), 101–110. DOI: 10.1021/ar1000956.
• A follow up question: How in this case (especially when compare several levels of theory) treat imaginary frequencies around i10 - i80 that sometimes appear? They do not affect further IRC calculation, but still exist. @CodyAldaz, can you also comment? Apr 30, 2020 at 20:20
• @mit.eremin If you have these small modes, your grid might be too small. Or numerical noise, try disturbing the structure to get them out. I'd generally treat them as local minima. It's probably case by case... (Also, note that you can only notify one person in a comment, and that person has to be part of the conversation. So you did not reach Cody.) Apr 30, 2020 at 22:17
• Very good practical pieces of advice. Most people misses a lot out not using “workhorse DFAs”, proper grid etc
– Greg
Feb 14 at 3:55

If you're using Gaussian you can also use the QST2 and QST3 methods: the first uses reagents and products as input and try to think of a reaction coordinate to find the MEP; the former uses reagent, product and a guess of the TS as inputs instead.

If you want to change computational software, on the other hand, there is ORCA which allows you to compute an NEB simulation that takes as input the reagent and the product and some images (that will be computed by the software) and tries to fit this on the PES of the system finding the maximum in energy in quite a few iterations.

• This is a great suggestion and a good comment. Unfortunately I think it doesn't really answer the question. To stand alone, there should be something about the level of theory in here. Doing NEB with semi empirics is quite fast and produces good starting structures; using something else could take forever and yield nothing "better". Feb 13 at 13:32
• @Andrea I think you meant QST2 and QST3 instead of QT2 and QT3. Feb 13 at 15:11
• @SRMaiti yes, thanks to have it pointed out Feb 13 at 17:57
• @Martin-マーチン you can do a NEB with every theory level you want. It's a bit more computational expencive Feb 13 at 17:58
• Yes, obviously, but it really depends on the focus whether this is necessary or not. It you're interested in a transition state, a fully converged NEB is most likely overkill. By the way, I thought of including this advice in my answer, but I realised: it is already there. I retain my initial statement that this is not a complete answer and it doesn't even add anything new. Feb 14 at 19:17