Is it okay to use semi-empirical methods to optimize structures if I perform high-level energy calculations?

Question

For non-covalent systems comprised of small-molecule monomers, folks are usually interested in understanding the interaction and binding energies of small (gas-phase) clusters ... things like $\ce{X(H2O)_N}$ where $N$ typically runs from 1-10.

For pure water or small molecules in water, I would say the current "benchmark structures" are usually optimized with something like RI-MP2 or DF-MP2. For larger datasets, I would say B3LYP-D3(BJ) is still quite popular, still with $N\approx$ 10-12.

Another approach that is popular, is to perform calculations on structures that were extracted from MD simulations. This makes sense when averaging over many clusters or when $N$ is large.

I'm wondering what you think about something "in between"... like interaction/binding energies using something "high-level" [e.g. DLPNO-CCSD(T)] over 1 set of structures $\ce{X(H2O)_N}$ where $N$ = 1 - 20, but were optimized using a semi-empirical method like DFT-3c [e.g. PBEh-3c]

In my opinion, this seems like an okay first approximation since the analysis could be refined if needed or complemented by selecting one of the two typical approaches. But on the other hand, since this is somehow in between, should I just pick one of the other approaches?

Part of my initial thought is that -3c methods are made to save us time while giving decent structures, so why not? Thoughts?

Answer 1

I generally try to use the most accurate methods that my resources will allow me to use, within a reasonable budget for the desired accuracy.

If you're interested in 20 water molecules or 60 atoms (before including your $\textrm{X}$ atom), then even the cost of DLPNO-CCSD(T) might be quite brutal, so it might as well be done on a decent geometry the first time, so that you don't have to redo the calculation when you later learn that your original geometry was badly flawed. However, I would certainly not want to use RI-MP2 or DF-MP2 for the geometry optimization for such a large system.

This would seem to limit us to the last two options that you mentioned for the geometry optimization: MD and semi-empirical methods. The MD would be quite fast for 60 atoms, so you can obtain the MD structures for $N$ = 1-20 anyway, with what I think would be negligible cost compared to the DLPNO-CCSD(T) calculations that you're about to do. You could then try optimizing the structures with a semi-empirical method (I would pick one and stick with it, rather than repeating potentially expensive calculations too many times) with N=1,2,3..., checking each time how different the results are compared to the MD calculations. If the results are mostly unaffected by the optimization method, then there's no need to do heavy calculations all the way up to N=20, since the cheaper calculations are probably just as good and are already available. If the results are drastically different, then you have another question to ask: how much more accurate are they? Personally I think that there's already so many more approximations going on than you have even mentioned:

• which basis set are you using?
• which basis set extrapolation scheme are you using?
• does the $\textrm{X}$ atom require a relativistic treatment or does it cause the need for considering multi-reference character that is not amenable to using CCSD(T)?
• how big is the CCSD(T) error from FCI?
• how big is the DLPNO-CCSD(T) error from CCSD(T)?
• how close is the optimized semi-empirical structure from the real structure or the DLPNO-CCSD(T) structure?

After considering all of those factors, if you still think a more expensive geometry optimization is worth the computer time and global warming that it causes, and that the delay in publishing the paper is worth it because of a meaningful gain in accuracy, then by all means do it! If your question is whether or not it's "okay" to do it, I see no reason why it wouldn't be okay! But considering that other errors such as the basis set incompleteness error and the correlation incompleteness error will still be there no matter what, I might just be inclined to stick with the more economical geometry optimization method (although you'll know better about the expected CPU time and value of it once you start doing the actual calculations: maybe I'm worried about nothing and the geometry optimizations will be lightning fast with both options, or maybe you'll see it being prohibitively slow at a small $N$ value and your decision will become immediately obvious).

You can also just pick a small basis set to test everything (all candidate methods that you're considering) and see how much things matter with the small basis set, where everything is lightning fast (comparatively). This will give you very valuable guidance on what to do for the "real" calculations with the bigger basis sets (if you care about DLPNO-CCSD(T) type accuracy, then you're probably planning to do at least a DZ,TZ->CBS extrapolation, but you can practice and explore things with things like 6-31G and 3-21G first).