Correct jobs sequence for transition-state (TS) search

Question

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**)

?

Answer 1

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.

Answer 2

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.

Answer 3

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.