Best Practice DFT Protocols for Basic Molecular Computational Chemistry

Markus Bursch,∗a Jan-Michael Mewes,∗b , Andreas Hansen∗c , and Stefan Grimme∗d

Abstract: electronic-structure level, which is possible only through chem-

Nowadays, many chemical investigations can be supported the- ical synthesis, spectroscopy, and quantum chemical calculation.
oretically by routine molecular structure calculations, conformer The latter of these three, computational chemistry and specifi-
ensembles, reaction energies, barrier heights, and predicted spec- cally Kohn-Sham density functional theory (DFT), has firmly con-
troscopic properties. Such standard computational chemistry ap- solidated its position as a third work-horse besides synthesis and
plications are most often conducted with density functional the- spectroscopy in recent decades. We argue that the general impor-
ory (DFT) and atom-centered atomic orbital basis sets imple- tance of computational chemistry and DFT, in particular, stems
mented in many standard quantum chemistry software packages. from its outstanding effort-to-insight and cost-to-accuracy ratios
This work aims to provide general guidance on the various tech- compared to related approaches, in other words, their efficiency
nical and methodological aspects of DFT calculations for molec- (vide infra).
ular systems, and how to achieve an optimal balance between From a more fundamental perspective, DFT is a formally ex-
accuracy, robustness, and computational efficiency through multi- act but practically empirical "first-principles" electronic-structure
level approaches. The main points discussed are the density func- approach to solve the fermionic many-electron problem that un-
tional, the atomic orbital basis sets, and the computational proto- derlies most of chemistry and large parts of biology and physics.
col to describe and predict experimental behavior properly. This is When applied together with a mixed quantum-classical treatment
done in three main parts: Firstly, in the form of a step-by-step de- for the nuclei using molecular dynamics (MD) or harmonic ap-
cision tree to guide the overall computational approach depend- proximations for the potential energy surface (PES), DFT can ad-
ing on the problem; secondly, using a recommendation matrix dress many problems in (bio)chemistry and physics with suffi-
that addresses the most critical aspects regarding the functional cient accuracy to derive meaningful insight. The subtle theoret-
and basis set depending on the computational task at hand (struc- ical and technical details of DFT are well understood and have
ture optimization, reaction energy calculations etc.); and thirdly, recently been discussed in an extensive open discussion type of
by applying all steps to some representative examples to illustrate review 1 . This open review addresses all more detailed and theo-
the recommended protocols and effect of methodological choices. retical questions about DFT which are raised here but whose an-
swer is beyond the scope of this practically-oriented work. For a
recent discussion of the grand challenges in theoretical and com-
putational chemistry, see Refs. 2,3.
DFT offers an excellent compromise between required com-
putation time and the quality of the results in comparison to
Introduction the alternatives, which are less accurate and robust but much
Chemistry and chemical synthesis are indispensable tools for hu- faster semi-empirical quantum mechanics 4–6 (often termed SQM)
mankind in addressing the most urgent current and future chal- on the one hand, and on the other more accurate and robust
lenges, such as efficient energy storage and conversion, sustain- but slower wavefunction theory (WFT) based approaches such
able food supply, and affordable medication and health care. In as coupled-cluster, (see Fig. 1). Further, human and power
all of these examples, a rational design of molecules and ma- resources are spared. Moreover, and just as important, DFT
terials, e.g., new catalysts, electrolytes for batteries, hosts and can be considered a robust theory in that a breakdown in the
emitters for organic electronics, and new drugs, takes a central form of entirely wrong results is scarce, even when applied
role. Here, it is crucial to understand matter at the atomic and to challenging molecules or exotic chemistry. This contrasts
semi-empirical quantum mechanics and other, even empirical
approaches, which require much more careful sanity-checking,
[a] Max-Planck-Institut für Kohlenforschung, Kaiser-Wilhelm-Platz 1, D-45470 Mül- which typically means a comparison to DFT. These properties
heim an der Ruhr, Germany. E-mail: bursch@kofo.mpg.de grant DFT the role of a black-box method that non-experts can
[b] E-mail: janmewes@janmewes.de apply to many chemical problems. Such applications typically
[c] E-mail: hansen@thch.uni-bonn.de
involve sanity checking suggested structures and reaction mech-
[d] E-mail: grimme@thch.uni-bonn.de
[b,c,d] Mulliken Center for Theoretical Chemistry, Institut für Physikalische und Theo- anisms by synthetic chemists or visualizing the frontier orbitals
retische Chemie, Universität Bonn, Beringstraße 4, 53115 Bonn, Germany. and HOMO and LUMO energies to interpret electrochemical or

1
optical experiments in materials science. interactions in the natural bonding orbital (NBO) 22 framework.
Nevertheless, the choice of a reasonable, efficient yet accurate The conclusions herein are based on more than 25 years of ex-
quantum chemistry treatment for a wide range of chemical prob- perience in the field of DFT and functional development, ther-
lems is still a challenging task, even for experienced computa- mochemical benchmarking, as well as hundreds of collabora-
tional chemists. This is also due to a vast number of available tive chemical applications in mechanistic or spectroscopic studies.
method combinations that have been developed and presented Our recommendations are mostly based on hard evidence, that is,
in recent years. Already for the fundamental choice of the den- they rely on large-scale comparisons of approximate DFT results
sity functional/atomic orbital basis set combination, hundreds or for a wide range of chemical properties with those from experi-
even thousands of combinations are possible in typical programs, ments or highly accurate and robust coupled-cluster theory, the
many of which are in common use. While this is uncritical for so-called gold-standard in computational chemistry (benchmark-
a small molecule where one can use a "sledgehammer to crack ing). However, some conclusions for the performance of a model,
a nut" (that is, use expensive double-hybrid density functionals basis set or particular density functional are numerically and sta-
and large atomic orbital basis sets), the treatment of systems with tistically challenging to quantify. Hence, our recommendations
50-100 atoms or many relevant low-energy conformers demands include a personal experience-based flavor. In this context, we
critical compromises in methodological choices in order to keep put more emphasis on the robustness of a method than its "peak-
the computational cost manageable. performance" reflected, for example, in the ranking in standard
Unfortunately, in some QM programs, the default methods thermochemical benchmark sets such as the GMTKN55 23 . This
are outdated, which may tempt inexperienced users to apply is because in our experience in predictive applications, robust-
no longer recommended methods to circumvent these compli- ness and reliability, that is, avoiding large and unexpected er-
cations. A prominent example is the popular B3LYP 7,8 /6-31G* rors, is more important than getting the numbers right to the last
functional/atomic orbital basis set combination that is still fre- kcal·mol−1 .
quently used even though it is known to perform poorly even Finally, we want to point out that in computational chemistry,
for simple cases. 9–11 The knowledge that B3LYP/6-31G* suffers the wording "thermochemistry" also includes the calculation of
from severe inherent errors, namely missing London dispersion reaction barrier heights, although they might be considered sepa-
effects ("over-repulsiveness") and strong basis set superposition rately under the common term chemical kinetics.
error (BSSE), seems to "diffuse" exceptionally slowly from the the-
oretical to the computational chemist community. In the last 10- 1 General considerations
20 years, the availability of better functionals, 12,13 standardized Decision making in computational chemistry
dispersion corrections, 14 and empirical corrections for BSSE 15,16
Defining a suitable set of theoretical methods is the key to accu-
made B3LYP/6-31G* computations obsolete. Today, much more
rately describing a chemical system. This includes not only the
accurate, robust, and sometimes even computationally cheaper
selection of a quantum chemistry model for the basic electronic
alternatives exist, e.g., in the form of composite methods like
structure but also the choice of an appropriate model system to
B3LYP-3c 17 , r2 SCAN-3c 18 , B3LYP-D3-DCP 19 , or B97M-V/def2-
describe the physico-chemical totality of the problem. Generally,
SVPD/DFT-C 16 to name just a few. Such methods use new de-
these choices comprise a complex set of fundamental decisions.
velopments to eliminate the systematic errors of B3LYP/6-31G*
An exemplary flowchart that illustrates this decision making and
without increasing the computational cost. Fig. 1 illustrates the
is applicable to a large part of typical quantum chemical applica-
relative computational demands of the discussed approaches, and
tions is shown in Fig. 2.
section 3.2 compares some of them to B3LYP/6-31G* for a non-
covalently bound complex. The main goal of this work is to in-
troduce these new methods and developments to interested non- Electronic Structure
experts and chemists with some background in theory. The possibly most fundamental aspect that decides if common
To this end, best practice recommendations are provided for DFT is applicable is whether the system under consideration is
the most common workflows encountered in typical applications, well-represented by a single-determinant wavefunction and thus
i.e., structure and ensemble determination, computation of re- has single-reference character, or if multiple determinants are re-
action energies, barriers, free energies and solvation effects are quired and the system has multi-reference character. Luckily, the
illustrated with a few typical examples or case studies. We will most common examples fall into the first category. These sys-
also briefly discuss the "embedding" of properly conducted DFT tems possess a single-reference electronic structure and are thus
calculations into multi-level workflows to determine solvated and readily describable by the common DFT methods discussed in this
thermally averaged conformer, protomer, or tautomer ensembles. work. This holds true in particular for diamagnetic closed-shell
However, due to quantum chemistry’s broad range and complex- (organic) molecules, the vast majority of which possess single-
ity, some interesting but less common topics like theoretical spec- reference character. The few exceptions, such as biradicals, often
troscopy, excited states, periodic systems, or heavy element chem- have low-lying triplet states, which can be checked with an unre-
istry are not considered here or touched only very briefly. Another stricted broken-symmetry 27,28 DFT calculation.
important area not covered here is the analysis and understand- Systems where multi-reference character should be expected
ing of atomic and fragment interactions by, e.g., energy decompo- are radicals, 3d transition metal complexes, low band-gap
sition analysis, 20,21 , or wavefunction composition, e.g., by orbital systems (see below), and transition states of open-shell dis-

2
 a discussion of multi-reference vs. related multi-determinantal
cases appearing in low-spin open-shell systems see Ref. 37.

< 0.5 CC Recommendations:

• Check for multi-reference character through simple indi-

accuracy / kcal·mol−1

double-hybrid cators (HOMO-LUMO gap, fractional-occupation-number-

1 "chemical accuracy"
weighted density).
hybrid
r2SCAN-3c
• Be more careful with open-shell systems (low-spin in partic-
(m)GGA
5 ular).

B3LYP/6-31G*
• Do not apply single-reference methods like DFT to multi-
SQM reference systems.
> 10
• ≥ TZ
• def2-mTZVPP Solvation
• minimal
The next fundamental question is that of the present state of ag-
< 0.01 1 10 100 gregation or in which form of a substance mixture the molecule
relative computation time to be examined is present. The neighboring molecules in a solid
or for a solute in solution can have drastic effects on the struc-
Figure 1 Accuracy in typical thermochemical applications vs. computa- ture and properties of the entire system. Accordingly, for con-
tion time (logarithmic scale) for common quantum chemistry methods. densed phase chemistry, a suitable solvent model should be ap-
Wavefunction based coupled-cluster (CC) methods like CCSD(T) (possi- plied in any case. The most common approach in a DFT con-
bly in combination with local approximations) provide benchmark quality
results of better than 1 kcal·mol−1 for common chemical energy changes
text is to use continuum solvation models that include interac-
of almost any well-behaved single-reference system (for details see text). tion of the molecule with the solvent implicitly via an effective
Lower rungs of the DFT hierarchy of methods from small basis set com- potential in the Hamiltonian. This means that no actual solvent
posite approaches 24–26 (e.g., r2 SCAN-3c 18 ), to meta-GGA or hybrid molecules are present in the calculation. Prominent representa-
functionals provide systematically improved results when coupled with
tives of this class include the conductor-like polarizable contin-
large atomic orbital basis sets. The most sophisticated double-hybrid
functionals often yield results close to a coupled-cluster reference level. uum model (CPCM) 38 , the solvation model based on the molecu-
See the "Choice of functional" section for a more detailed discussion of lar electron density (SMD) 39 , the conductor-like screening model
the density functional classes. SQM = semi-empirical quantum mechan- (COSMO) 40 , the conductor-like screening model for real solvents
ical method. (COSMO-RS) 41 , and the direct conductor-like screening model
for real solvents (DCOSMO-RS) 42 . The mentioned methods dif-
fer in various aspects. CPCM und COSMO are purely electro-
sociation processes. In all of these, multiple near-degenerate static models, lacking contributions from cavity creation that cost
electronic configurations (with different orbital occupations) can energy in the solvent, and attractive van-der-Waals interactions
be present, leading to complicated electronic states. Accordingly, with the solvent, which lead to substantial errors if the solvent-
for these cases, it should be checked in advance whether a accessible surface area is changing significantly. SMD, COSMO-
so-called multi-reference case is present. While this test is RS, and DCOSMO-RS include such contributions and are thus rec-
generally rather complicated, there are simple hints that allow an ommended (cf. Section 3.2). Nevertheless, it should be noted
initial yet often sufficient assessment of possible multi-reference that COSMO-RS cannot be used in geometry optimizations or fre-
character. For example, a very small gap between the highest quency calculations and has to be replaced by DCOSMO-RS for
occupied and lowest unoccupied molecular orbitals, the so-called this purpose. Further noteworthy implicit solvation methods that
HOMO-LUMO gap, of <0.5-1 eV in a test GGA calculation (see can also be used with semi-empirical quantum mechanical and
below) and exceptionally slow self-consistent field convergence force-field methods are the generalized Born model with solvent
already give a first indication of unusual electronic complexity. If accessible surface area (GBSA) 43,44 and the analytical linearized
this first crude hint emerges, more sophisticated measures should Poisson-Boltzmann model (ALPB) 45,46 .
be considered. Most available approaches are wavefunction Nevertheless, for specific cases, the inclusion of explicit solvent
theory based 29–33 and thus limited in their applicability to molecules may be necessary and implicit solvation models
small molecules. The alternative fractional-occupation-number- become insufficient. 47 In micro-solvation approaches, actual
weighted density (FOD) 34–36 represents an easily applicable solvent molecules are placed at important, most strongly bound,
DFT-based estimate. Here, a moderate artificial increase in the positions of a system. 48 However, explicit solvation also has
electronic temperature can be used to populate and visualize its caveats as it can be very difficult to converge properties
low-lying, possibly problematic electronic states. With any with the number of explicit solvent molecules. Moreover, the
indication of a significant multi-reference character, application potential-energy surface of explicitly solvated systems is often
of standard DFT methods is not recommended and experts for times flat and peppered with local minima of different solvent
sophisticated multi-reference theory should be consulted. For structures, making optimizations lengthy and tedious. Therefore,

3
 methodological decision
single-reference electronic structure?
key question
yes no

gas-phase? apply sophisticated multi-reference methods

yes no

choose solvation model

describable by a single structure/conformer?

yes no

choose quantum chemistry code

apply conformational search algorithm or
molecular dynamics averaging and use ensemble

choose functional

heavy elements (Z > 36) involved?

yes no
consider relativistic treatment

choose basis set

optimize geometry or ensemble

compute final energies or desired property

enthalpy or Gibbs free energy desired?

yes no

apply zero-point vibrational energy and use computed energies

thermostatistical corrections

Figure 2 Conceptual flowchart of decision making in elementary steps in typical computational chemistry calculations.

explicit solvation should be used with care. Neglecting solvation • Consider explicit solvation if necessary.
effects, specifically for polar or charged molecules, can result
in large deviations in thermochemical calculations, and even
fundamentally wrong electronic structures, e.g., for zwitter-ions. Molecular Flexibility

Another important aspect is the structural flexibility of the

Recommendations: system. For highly flexible structures, a molecular property, such
as energy, nuclear magnetic resonance spectra or optical rotation
• Choose a model/state of aggregation close to the experi-
values, may not be sufficiently described by a single structure.
ment.
At finite temperatures, various conformers are populated and
• Apply implicit solvation models for a molecule in solution. the overall property must be described as thermal average over
Best use physically complete models, such as COSMO-RS or the unique property values of each conformer. Accordingly, it is
SMD. recommended to evaluate the flexibility and the accessibility of
relevant conformers for any given system. The flexibility of a
• Be careful with charged systems where continuum models molecule may be roughly categorized by the number of conform-
may be inaccurate (the higher the charge density the more ers in an energy window of 3 kcal·mol−1 , which corresponds to
inaccurate). five times the thermal at room temperature, 5 × RT , with respect

4
to the lowest energy conformer. Systems with few conformers
(≈ 1-3) in this window may be considered relatively rigid, those
"heaven" of chemical accuracy
with dozens as intermediate cases, and with hundreds as very
flexible. In any case, finding the overall lowest conformer (global
minimum) in the given chemical environment is important. 49
For example, a conformer found in the solid, e.g., determined
by X-ray crystallography, may not be the most favorable one in
5 virt. ϕ(r) double-hybrid PWPB95
B2PLYP

solution or gas-phase. 50 In general, even for only medium sized

computational cost / accuracy

4 (Fock-exc.)
occ. ϕ(r) hybrid ωB97X-D
molecules (30-50 atoms), a sophisticated conformational search (global, range-sep., local) B3LYP
is not trivial, computationally demanding, and usually prohibitive
at a pure DFT level. Hence, multi-level approaches (see section
5) involving efficient semi-empirical quantum mechanical or
force-field (FF) methods are necessary. The CREST 51 /CENSO 52
3 𝜏(r) meta-GGA
TPSS
r2SCAN

approach represents a valuable, easily applicable tool for semi-

automated conformation sampling and subsequent energetic
ranking of conformer-rotamer ensembles (CREs). Further, the
2 ∇ρ(r) GGA
PBE
B97-D

flexibility index given by the CREST program can be used as

additional indicator for the molecular flexibility. Besides CREST,
alternative, less general conformer generation procedures are 1 ρ(r) LSDA SVWN

described in the literature. 53–59

Hartree "hell"
Recommendations:

• Check for the role of structural flexibility/conformations.

Figure 3 Functional categorization according to Perdew’s “Jacob’s lad-
der”. ρ = electron density, τ = kinetic energy density, φ = molecular or-
• Apply automated conformation search algorithms, e.g., bital, Fock-exc. = Fock exchange.
CREST.

• Try to find and verify the lowest energy conformer. Even though these categories are based on fundamental theo-
retical aspects, the functionals within each rung can vary strongly
• Consider Boltzmann-averaged property calculations. in accuracy. Accordingly, comprehensive benchmark studies that
assess the performance of any functional with respect to the de-
Choice of Functional sired target property are indispensable. Nevertheless, the Jacob’s
A critical issue in DFT is the choice of the basic exchange- ladder categorization allows a crude estimation of some system-
correlation energy functional, often simplistically called func- atic functional errors. Moreover, it allows to categorize function-
tional. It aims to absorb all extremely complicated many-particle als by their numerical efficiency. In this regard, the perhaps most
correlation and fermionic (exchange) effects into a seemingly important distinction is made based on the inclusion of Fock ex-
simple but formally exact and theoretically existing mean-field change, also termed non-local or exact exchange. The (meta-
electronic energy and potential. During the last decades, hun- )GGA functionals of up to rung 3 do not include Fock exchange
dreds of different functionals were constructed which vary in their and are called local or semi-local functionals, whereas (double-
conception, target application, and their overall quality. 12 Ac- )hybrid functionals on rungs 4 and 5 include Fock-exchange.
cordingly, there are general-purpose as well as task-specific func- Since the calculation of Fock-exchange is a computational bot-
tionals, highly parameterized and fundamental first-principle- tleneck, semi-local functionals are generally more efficient than
based ones. A detailed discussion of those aspects is beyond hybrid functionals (see below).
the scope of this work as we focus on the suitability for appli- Two of the most critical and thus prominent errors in actual
cations in thermochemical calculations. The most prominent at- DFT approximations are the so called self-interaction error (SIE),
tempt to categorize density functionals based on their physical and missing long-range correlation effects that give rise to Lon-
ingredients represents Perdew’s “Jacob’s ladder” (Fig. 3). 60 Here, don dispersion. Although the lacking description of long-range
functionals are ranked according to their degree of approximation correlation is a fundamental shortcoming of DFT, it can nowa-
as measured by the included electron density descriptors for the days easily be fixed by including one of several available proven
exchange-correlation term (occupied orbitals via the density, first dispersion corrections (we recommend D4, D3, or VV10). 14 We
derivative of the density, second derivative, occupied orbitals via argue that nowadays, this is indispensable in any DFT treatment,
Fock exchange, virtual orbitals via MP2, ...) and thus the expected and do not see any context in which the dispersion correction
accuracy. The most relevant categories ordered by increasing ac- should be left out. See section 2.7 for a discussion and example
curacy include the generalized-gradient-approximations- (GGA), 3.2 and some illustrative numbers.
meta-GGA- (mGGA), hybrid-, and double-hybrid functionals. The SIE 61–63 results from an artificial interaction of an electron

5
with its own mean electric field by the approximate exchange- ficiently to Fock exchange, which thus dominates the calculation
correlation functional and is more difficult to repair. In con- time. However, there are semi-numerical integration techniques
trast, Hartree-Fock theory is SIE-free, because of an exact math- such as chain-of-spheres-exchange (COSX) 77 that also enable sig-
ematical relation between the integrals describing the Coulomb- nificant speedups in hybrid calculations for basis sets of TZ quality
and (exact) exchange interaction: The two exactly cancel each or larger (factor of 3-10) if combined with RI. Thus, we strongly
other when an electron formally interacts with itself. In prac- recommend to use both RI and COSX if available (e.g., in ORCA
tice, SIE in DFT typically results in an over-delocalization of the as RIJCOSX and TURBOMOLE as SENEX 78 ). The RI technique
electron density and artificial energy stabilization for delocalized also speeds up the MP2 calculation part in double-hybrid calcu-
electronic situations (with bond stretched H+ 2 as prime example). lations which then become only slightly more costly than hybrids
The SIE is present in all semi-local (m)GGA functionals (rungs 1- and hence this safe approximation is strongly recommended. In
3, e.g., PBE 64 ). Hybrid functionals attempt to reduce the SIE by any case, consistency checks of the chosen functional are recom-
replacing a fraction of approximate DFT exchange with SIE-free mended by comparing the results of a few functionals from each
Fock exchange (e.g., B3LYP 7,8 with 20% and PBE0 65 with 25%). class for a representative model system.
However, this only reduces the error but does not eliminate it, A note on Minnesota functionals – Beginning in 2005 a series
such that functionals employing small fractions of Fock exchange of functionals was developed by the Truhlar group in Minnesota,
are still prone to SIE. Nevertheless, simply employing 100% of which are typically termed MXX-Y where XX stands for the
Fock exchange also leads to a very poor performance. One way year and Y informs about the purpose (e.g., HF for 100% Fock
around this dilemma is provided by range-separated hybrid func- exchange, "2X" for twice the Fock exchange, "L" for local). 79–83
tionals, which do not admix a constant amount of Fock exchange, These functionals are widely used and perform very well on large
but make the admixture dependent on the inter-electronic dis- main-group benchmark sets like the GMTKN55, where M06-2X
tance. 66,67 In general, Fock exchange admixtures of 5-20% are and M05-2X are actually some of the best-performing global
typically considered small, 20-30% moderate, and >30% high. hybrids, and M06-L is one of the best meta-GGAs. These methods
The SIE can be viewed as a special form of the more general represent a significant step forward in functional development
delocalization error of semi-local functionals resulting from the compared to B3LYP, and solve some complicated electron corre-
incorrect description of regions in a molecule with (effectively) lation problems in DFT such as alkane branching. 84 However,
fractional charges. 63,68,69 Regarding the understanding of these we do not explicitly recommend them for several reasons: Firstly
problems, some progress has been achieved recently by distin- and perhaps most importantly, these functionals are not as
guishing density driven errors caused by the erroneous (m)GGA robust and black-box as others and thus require the user to pay
exchange-correlation potential and inherent, usually smaller er- attention to technical details. This is because they are often very
rors of the energy functional. 70 sensitive to the size of the integration grid and the basis set, 85
Double-hybrid functionals 71–73 represent the highest rung and which may lead to discontinuities in potential energy surfaces
additionally introduce a wavefunction-theory-based correction to and, in turn, problems with geometry optimizations 86 . Secondly,
the correlation energy. This is most commonly achieved by per- their performance strongly depends on the chemical system. This
turbation theory methods 74 such as second order Møller-Plesset is because their very good performance for typical main-group
theory (MP2) 75 or its DFT variant 74 . Their typically very high chemistry relies on a quite extensive parameterization, which
amounts of Fock exchange (>50%), make them particularly re- may lead to problems for less common systems reflected in the
silient towards SIE and negative influences of large Fock exchange "mindless" benchmark, 23 or for transition-metal chemistry. 36,87
admixtures are balanced by the explicit, virtual orbital dependent Thirdly, these functionals can be problematic for noncovalent
correlation treatment. Nevertheless, these methods have an in- interactions. Although they are designed to include dispersion
creased computational cost and further introduce some restric- effects at an electronic level, and do, in fact, work quite well for
tions in their general applicability. For example, the vulnerability weakly bound systems at their equilibrium distances (see exam-
of the MP2 part for small gap systems may also be problematic ple 3.2), they cannot recover the correct asymptotic behaviour
in MP2 based double-hybrid functionals and hence, treating such of London dispersion in the long intermolecular distance regime.
systems in this way requires some caution. To mitigate this, a dispersion correction needs to be added in
Opposing the expected higher accuracy for fourth-rung hybrid certain situations, which, however, can also result in overbinding
and fifth-rung double-hybrid functionals, their increased compu- in others. 88
tational cost and less favorable scaling are important criteria,
specifically for large systems. (m)GGA treatments formally scale Recommendations:
with the system size (N) as N 3 if the resolution of the identity
(RI) approximation, also known as density-fitting, 76 is applied. • Choose a functional with caution based on the chemical sys-
For local (m)GGA functionals, RI can yield speedups by a fac- tem under investigation and the task at hand and not based
tor of 5-30, depending on system and basis set. Hybrids already on popularity.
scale with N 4 and MP2 based double-hybrids with N 5 and accord- • Always include a dispersion correction.
ingly, some higher rung functionals become unfeasible for large
systems (>200-300 atoms) on common hardware. Hybrid calcu- • Check for reliable cost-benefit-combinations and consider
lations do not profit from RI as much since it cannot be applied ef- (m)GGAs. Hybrid functionals are more accurate but also

6
 much more expensive than (m)GGAs. covalent bonds, the so-called counter-poise correction can be ap-
plied to correct for BSSE. 91 An efficient alternative to this com-
• Check the consistency between different functional classes
putationally demanding correction is provided by approximate,
(e.g., compare a hybrid and a mGGA).
empirical correction schemes that are based on the molecular
• In critical cases, test hybrids including different amounts of structure, such as the geometric counter-poise correction (gCP),
Fock exchange. or employ specially adapted effective core potentials 92 In con-
trast to the full counterpoise corrections, these are always appli-
• Consider proven multi-level approaches and composite cable and computationally cheap, and thus can also employed to
methods for larger systems (see Section 2). correct for the intramolecular BSSE. Such approximate counter-
poise corrections can repair the most drastic effects of BSSE, e.g.,
Choice of Basis Set in geometry optimizations with small basis sets. Their applica-
Another important aspect regarding the computational tion is similarly straightforward as that of dispersion corrections.
speed/accuracy compromise is the applied atomic orbital Thus, we recommend the gCP approach of Kruse and Grimme 15
basis set. From a fundamental point of view, this is merely a tech- for geometry optimizations with small basis sets that supports HF
nical aspect, because DFT calculations can at least in principle and DFT as well as many basis sets. The related DFT-C approach
be numerically converged to the complete basis set (CBS) limit of Witte and Head-Gordon is further recommended for accurate
where this influence is eliminated. In practice, however, this is noncovalent interaction energy calculations (adapted specifically
rarely done and finite basis sets are applied, introducing some for DFT/def2-SVPD calculations, see example 3.2), 16 . Other no-
errors. However, basis set related errors in DFT are typically table concepts are the proximity-function of Faver and Merz for
much smaller than for correlated wavefunction theory-based large biomolecules, 93 as well as the ACP-n approach of Jensen. 94
methods. This weaker basis set dependency compared with The most commonly used Gaussian type contracted basis
wavefunction theory is a strong point of DFT. Nevertheless, even sets belong to the Pople (e.g., 6-31G), 95 Dunning (cc-pVXZ) 96 ,
the best functional will yield bad results if evaluated in a small Jensen (pc(seg)-X) 97,98 , and Ahlrichs (def2-XVP) 99,100 basis set
and insufficient basis set. families. As a somewhat technical side note we want to mention
The most important characteristic of basis sets is their com- that Pople-type basis sets such as 6-31G* or 6-311G** as well
pleteness, often referred to as basis set size reflecting the number as the Dunning-type sets cc-pVXZ (X=D, T , ...) are not recom-
of functions to represent a given electron. An important error for mended here for standard DFT treatments, mainly because more
too small basis sets is the so-called basis set incompleteness error efficient and consistent alternatives by Ahlrichs and co-workers
(BSIE) 89 , which results from an insufficient function space in the are available. 101 All basis sets may be augmented by additional
linear combination of atomic orbitals expansion. In short, BSIE polarization that have a higher angular momentum or diffuse
arises when the employed basis set is not flexible enough to de- functions with small exponents to introduce more flexibility if
scribe the fine details of the electron density. In most cases, the necessary, e.g., for anions, dipole moments or electric polariz-
description of the valence electrons is most crucial and thus ba- abilities. For the recommended Ahlrichs basis sets, examples are
sis sets are usually categorized according to the so-called cardinal the def2-TZVPP basis set with added polarization functions and
number that indicates the number of independent basis functions def2-TZVPD 102 or ma-TZVP 103 (In some contexts also denoted
per occupied valence orbital. The corresponding size is usually as ma-def2-XVP) that are augmented with diffuse basis functions.
referred to as double- (DZ), triple- (TZ), quadruple- (QZ),...,-zeta To test if the results of a calculation are converged with respect to
where the term zeta refers to the number of filled electron shells the basis set size, it should be investigated how they change when
in the neutral atom (i.e., 1s, 2s, 2p, 3s, 3p, 3d, ...). Another ba- the cardinal number is increased (def2-TZVP to def2-QZVP), and
sis set related error arises if the basis set is too small and hence when polarization function are added or removed (def2-TZVP to
rather incomplete: Spatially close atoms and fragments start to def2-TZVPP).
“borrow” basis functions from each other, resulting in an artificial If heavy elements are involved, def2-XVP basis sets use the
energy lowering for more compact structures, which is known as matching Stuttgart-Cologne effective core potentials (def2-ECPs,
the basis set superposition error (BSSE). 90 BSSE is the practically by default for Z > 36) 104 to replace the inner core electrons.
most relevant error and commonly associated only with weak in- This not only reduces the computation time for very heavy
teractions and noncovalently bound inter-molecular complexes, elements, but also increases the accuracy due to the implicit
for which it can become relatively large being in the magnitude inclusion of (scalar)-relativistic effects, which mainly affect the
of the interaction energy. However, it is important to recognize core electrons. The application of robust small-core effective
that BSSE is always present and also affects, e.g., conformational core potentials is typically sufficient for most thermochemical
energies and even molecular structures if too small basis sets are property calculations. Nevertheless, for certain properties of very
used. Accordingly, knowing which computations are specifically heavy nuclei and properties involving core electrons such as NMR
sensitive to basis set size and where it is worth to invest the in- shieldings, explicitly relativistic all-electron calculations with
creased computational demand resulting from a larger basis set special Hamiltonians, such as X2C 105,106 , ZORA 107 , or DKH 108
to significantly improve the result is fundamental. This aspect is may be necessary, but are beyond the scope of this work.
the topic of the following Section 2 and Fig. 5. For BSSE prone
systems with clearly separated fragments with no inter-fragment Recommendations:

7
 • Be aware of BSIE and BSSE. Quantum Chemistry Program Packages
In the last decades, various more or less specialized program
• Try to approach a reasonable basis-set size (≥TZ) for energy- packages were developed. All of them provide the same basic
related properties. DFT functionality and relevant differences concern their ease of
use, their efficiency, and their availability (free or commercial).
• Consider adding polarization functions for flexibility (def2- In the following, a selection of common program packages with
XVPP). versatile functionality in molecular applications is given and
individual strengths and specialties are highlighted briefly. These
• Consider adding diffuse functions for anions, dipole mo- programs proved to be reliable for a wide range of quantum
ments, and polarizabilities (ma-XVP, def2-XVPD). chemistry applications in our group.

• Check for basis set convergence (increase/decrease cardinal TURBOMOLE 116

number by one).
• Very fast and robust
• If heavy atoms are present (Z > 36), apply effective core • Technically advanced and state-of-the-art algorithms
potentials.
• Molecular symmetry handling
Comparing Apples with Apples • Fast hybrid functional implementation by semi-numerical
Finally, a remaining critical point is the consideration of finite Fock-exchange approximation (SENEX)
temperature effects. Most computational models per construc-
ORCA 117–119
tion yield results valid for the absolute zero temperature (T =
0 K,−273.15°C), nuclear equilibrium scenario. Therefore, a sys- • Efficient implementation of DFT and wavefunction theory
tematic deviation compared to experimental data at finite temper-
• Fast hybrid functional implementation by semi-numerical
atures is expected as standard equilibrium structure treatments,
Fock-exchange approximation (RIJCOSX)
such as geometry optimizations, do not include effects such as
nuclear zero-point vibrational energy (ZPVE), vibrationally (ther- • Large toolkit for molecular spectroscopy
mally) elongated bonds, or molecular entropy. Accordingly, a per-
fect agreement between experimental structures obtained at finite • Easy and intuitive input structure
temperature and those calculated at T=0 K is not necessarily de-
• Free of charge for academic use
sirable since – sometimes substantial 109–111 – finite-temperature
effects can cause a significant bias. However, for common cova- Q-Chem 120
lent bonds between typical atoms, these effects are too small to
• Large number of implemented density functionals
have a significant influence.
More importantly, bare electronic energies calculated in the gas • Large variety of specialized DFT treatments (time-
phase (E) cannot be compared enthalpies (H) or free/Gibb’s ener- dependent-DFT, constrained-DFT) and analysis tools (en-
gies (G), which are typically measured in solution. This is because ergy decomposition analysis, EDA)
the zero-point vibrational energy that is contained in the internal
• Special DFT methods for NCI (SAPT, DFT-C)
energy U, thermostatistical corrections contained in H, entropic
corrections and free energy of solvation contained in G are of- Psi4 121
ten substantial for typical reactions and can take on the order of
• Modular code with great interfacing/scripting capabilities
several kcal·mol−1 (see examples in Section 3).
Recommendations: • Special DFT methods for NCI (SAPT)

• Try to model the experiment as closely as possible. • Free of charge

Molpro 122
• Apply zero-point vibrational energy and thermostatistical
corrections to enthalpy or free energy if necessary. • Advanced DFT (RPA, ACFDT) and embedding (WFT-in-DFT)
techniques
• For large systems with many small vibrational frequencies
• Special DFT methods for noncovalent interactions (SAPT)
<50-100 cm−1 , consider the robust mRRHO model for the
entropy part of the free energy. • Free of charge

As many quantum chemical applications typically involve man- AMS 123

ifold methods based on various physical concepts, further reading
• Use of Slater-type orbital (STO) basis sets
on general computational chemistry as well as specific aspects is
recommended. 5,26,73,89,112–115 • Good relativistic treatments

8
 • Huge quantum chemical toolbox reference (W1−F12)
−10
mGGA (r2SCAN−D4)
• Great graphical user interface (GUI) hybrid (r2SCAN0−D4)
−12
The widely used Gaussian program 124 package is a promi-
nent representative that is not specially recommended here as it −14

Eint. / kcal⋅mol−1
seems to be computationally not very efficient, lacking important
−16
recent technical advances, and better, even freely available,
alternatives with a similar or even higher functionality exist in 0.14e− 0.18e−
−18
our opinion.
−20
Specifically for large systems or very demanding computational
tasks, DFT methods can be complemented by semi-empirical −22 SIE
quantum mechanical methods. Even though some of these are
−24 Re
implemented in the mentioned quantum chemistry packages, sev-
eral features are only available in the corresponding original 2.6 2.8 3 3.2 3.4
codes. The most frequently used semi-empirical quantum me- RCMA / Å
chanics programs are xtb 125 , DFTB+ 126 , and MOPAC 127 . For
example, the recent single-point hessian (SPH) 128 approach can Figure 4 Potential energy surface along the cyclopropenyl–anthracene
center-of-mass distance RCMA for the mGGA r2 SCAN-D4, its hybrid vari-
be utilized with semi-empirical methods by using the native im- ant with 25% Fock exchange r2 SCAN0-D4 130 , and the W1-F12 refer-
plementation in the program. Further, the xtb program can be ence. W1-F12 denotes a highly accurate wavefunction theory-based ref-
used as a driver for other codes (using the xtb functionality with erence level. All DFT data calculated with the def2-QZVPP basis set.
parts of other codes to employ, e.g., DFT in the SPH approach) Re = equilibrium distance. Colored arrows indicate the charge transfer
from anthracene to the cyclopropenyl cation for the respective theoretical
such as ORCA. A detailed documentation can be found in Ref. 129 .
level.
2 The right tool for the task
A central task in computational chemistry is to balance compu- proach. One example for this is the anthracene-cyclopropenyl
tational demands against methodological accuracy and robust- cation potential-energy surface from our recent article introduc-
ness. To this end, it is important to be aware of strengths and ing r2 SCAN-3c 18 that is shown in Fig. 4. Although the mGGA-
weaknesses of specific methods, to consider the system size, and based composite method overestimates the interaction energy
the required accuracy of the target quantities such as structures, due to SIE, the equilibrium intermolecular distance is in very good
reaction energies, or conformational energies. Practically rel- agreement with the high-level coupled-cluster reference.
evant systematic errors in this respect are the aforementioned Recognizing that structure optimizations are generally much
self-interaction error (SIE), basis-set superposition (BSSE) and in- less sensitive to the level of theory than energy and many other
completeness errors (BSIE), as well as the lacking description of property calculations, multi-level approaches enables large com-
London dispersion by most functionals. The relevance of these putational savings without any significant loss of accuracy. To
aspects depends not only on the system under investigation but generalize from this example, we summarized the most suitable
even more so on the task at hand. Thus, it is instructive to discuss functional/basis-set combinations for typical steps of a computa-
available methodological choices in the framework of the most tional study in Fig. 5. To guide the choices visually, we marked
typical steps of a computational investigation. The aim of this the recommended level of theory with the best balance between
section is thus to guide the choice of the methodological tool set, computational effort and robustness in blue, accurate and robust
to provide the means to adapt it to the task at hand. Simply put: but not necessarily efficient choices in green, and methods that
We want to explain how one can cut corners where it does not should be avoided in red. Less clear-cut cases (in between red
hurt through clever methodological choices and the use of multi- and green) are marked in yellow. These methods can under cir-
level approaches. cumstances be a good choice but are not as robust as green ones.
An illustrative example for multi-level approaches is the use Hence, results obtained at this level should be checked for system-
of efficient semi-local (GGA and mGGA, i.e., no Fock exchange) atic errors as indicated in the field, and often better alternatives
functionals or composite methods for structure optimizations and (blue fields) are recommended instead. In all cases, the most se-
vibrational frequency calculations, combined with single-point vere systematic errors and drawbacks we expect based on our ex-
energy calculations at a hybrid/QZ level. The underlying idea perience are shown as text. In the following, we will discuss this
is that although energies calculated with (m)GGA functionals are table column (task)-wise to motivate and explain our choices, and
susceptible to SIE, they still provide very reasonable structures at provide additional details. To streamline this discussion, we fully
a small fraction of the computational cost of the more advanced separate the theoretical tasks in the computational context. This
non-local (hybrid and double-hybrid) functionals. This holds true means that, e.g., the discussion of methods for conformational
even for SIE-prone systems, such that single-point calculations energies only refers to calculation of electronic energies and ex-
with hybrid functionals on (m)GGA structures often provide en- cludes vibrational and entropic contributions, which we discuss
ergies and properties as accurate as a fully optimized hybrid ap- separately.

9
 functional class / basis reaction
structures frequencies conformers barriers NCI
recommendation size energies

DZ structureb BSIE SIE, BSIE SIE, BSIE SIE, BSIE

(meta-)GGA
TZ SIE, BSSE SIE, BSSE SIE, BSSE
B97M-V, r2SCAN a, TPSS
QZ excessive excessive SIE SIE SIE

DZ structureb BSIE BSIE SIEc BSIE

hybrid
TZ BSSE SIEc BSSE
PW6B95, PBE0, B3LYP
QZ excessive excessive SIEc

DZ structureb BSIE BSIE BSIE BSIE

range-separated hybrid
TZ BSSE BSIE, BSSE BSSE
ωB97M/X-V, ωB97M/X-D4
QZ excessive excessive
structureb,
DZ BSIE BSIE BSIE BSIE
double-hybrid excessive
TZ excessive excessive BSSE BSSE BSSE BSSE
PWPB95, ωB97M(2),
revDSD-PBEP86-D4
QZ excessive excessive

composite methods

B97-3c TZ SIE SIE SIE

r2SCAN-3c a TZ SIE SIE SIE

PBEh-3c d DZ BSIE BSIE BSIE

commonly used
combinations
structureb,
B3LYP/6-31G* DZ BSSE, no-D BSSE, no-D BSSE, no-D BSSE, no-D BSSE, no-D
no-D
BP86-D3 TZ SIE, BSSE SIE, BSSE SIE, BSSE

M06-2X/6-311G** a,d TZ
a
High dependency on the size of the numerical integration grid.
b
If the structure is reasonable (DZ ≈ TZ), GGA/DZ or even semi-empirical frequencies are sufficient. If not, use SPH.
c
Self-interaction error can limit the accuracy for hybrid functionals with small amounts of Fock exchange (< 20%).
d
The large amount of Fock exchange can cause pronounced errors for transition metal complexes.

not recommended intermediate accurate recommended / good cost-benefit-ratio

Figure 5 Decision matrix to guide the choice for a method combination (functional class/basis set) for common computational tasks. All considerations
imply the use of a dispersion-correction. Red, yellow, and green color code indicate the accuracy and reliability of a method for the given task, while
blue marks our recommendation based on a good cost/accuracy ratio. Text in the fields points out the most relevant, systematic errors we expect at
this level (SIE, BSSE, BSIE) or if a theory level is unnecessarily demanding (excessive). Excessive method combinations probably do not yield any
significant increase in accuracy justifying the much increased computational cost in the respective application. Selected recommended methods are
given with the respective method class. DZ = double-zeta; TZ = triple-zeta; QZ = quadruple-zeta basis set. BSIE = basis set incompleteness error;
BSSE = basis set superposition error; no-D = missing London dispersion ; SIE = self-interaction error; Fock exchange = non-local exchange (also
termed exact exchange) from wavefunction theory.

10
2.1 Structure tion (D3/D4 115,135,136 , VV10 137 ). Note that BLYP 138,139 is not
recommended because the GGA-typical slight systematic overes-
Structure optimizations are the first step of most quantum-
timation of covalent bond-lengths is more pronounced than with
chemical investigations, either starting from an experimental ref-
the recommended GGAs.
erence or a low-level guess that may employ semi-empirical quan-
tum mechanics or force-field methods. At this point, the most typ-
ical user errors concern the starting structure(s) and mistakes in 2.2 Frequencies
the input resulting in wrong charge, spin- or protonation states, The calculation of vibrational frequencies is indispensable to ob-
which drastically alter the results, whereas wrong choices of the tain zero-point vibration energies, thermostatistical corrections to
theoretical model are less common and often less severe. Hence, enthalpy and free energy, and also for the prediction of IR/Raman
we recommend to carefully check the input before setting the spectra, which are not in the focus here. Since the underlying cal-
computational machinery in motion. Self-consistent field conver- culation of second energy derivatives quickly becomes very de-
gence issues are often a strong hint towards problems with the manding for larger systems, it is desirable to conduct such calcu-
input. As already discussed above, structures are much less sensi- lations at the lowest level that still provides reasonable results.
tive to the functional and basis set than energies and properties. Therefore, we recommend to use efficient composite methods or,
Therefore, a well-balanced TZ basis set (def2-TZVP) is sufficient, alternatively, mGGA or hybrid functionals with a DZ basis set in
and going higher to QZ level is generally a waste of computational combination with the gCP correction. Similar to structures, error-
resources. Even well-balanced DZ basis-sets (def2-SVP) can pro- cancellation is very stable and, since the thermostatistical cor-
vide useful results when combined with the empirical geometric rection to a reaction energy is typically much smaller than the
counterpoise (gCP) 15 correction to mitigate the sometimes sig- electronic part from DFT, the impact of a lower level is naturally
nificant structural impact of BSSE. In fact, combining small but limited. In many cases, it is even sufficient to obtain thermosta-
well-balanced basis sets with the gCP correction is the basis of the tistical corrections at a semi-empirical quantum mechanics level,
PBEh-3c 25 (DZ, def2-mSVP) and r2 SCAN-3c (TZ, mTZVPP) com- e.g., with GFN2-xTB 5,140 , as shown in example 3.2. We note that
posite methods, which are tailor-made for the task of structure application of low-level semi-empirical methods should be care-
optimization and thus strongly recommended. In our opinion, it fully checked if exotic bonds, transition metals, or heavy main-
is rarely required to go beyond this level for structure optimiza- group elements are present.
tions in standard thermochemical studies. The only complication that arises concerning the use of lower-
Concerning the functional, (m)GGAs are typically sufficient level methods for frequencies is that the structure for which the
(e.g., r2 SCAN-D4 131,132 , TPSS 133 -D4, or even PBE-D4) as al- vibrational frequencies are calculated has to be an energy mini-
ready discussed in the beginning of this section. This provides mum at this very level (fully optimized, vanishing atomic forces).
the additional advantage to fully exploit the usually very accu- Otherwise, artificial imaginary frequencies will severely impact
rate resolution-of-the-identity (RI) approximation 76 , also called the accuracy of the calculated thermostatistical corrections. For
density-fitting, which makes (m)GGA calculations much more this reason, the same theory level has to be used for structure
affordable than using hybrids. As a result, mGGA function- optimization and vibrational frequency calculations or a second
als can be employed with larger basis sets at essentially the set of optimized structures obtained at this lower level has to be
same or even lower cost than hybrids with small basis sets like used ). This limitation can be overcome with the recently pre-
B3LYP/6-31G*. In turn, we argue that hybrid functionals should sented single-point Hessian (SPH) approach, which allows to ob-
only be used for structural optimizations if there is a good rea- tain frequencies and reasonable thermostatistical corrections for
son, such as strong SIE, application to transition state searches, any non-equilibrium structure through the application of a spe-
weakly bound electrons in anions or presence of heavy main- cific biasing potential. 128 While originally developed for the semi-
group elements. To investigate if hybrid functionals have a sig- empirical method GFNn-xTB, the SPH approach can also be used
nificant influence, we suggest to begin with comparing results for DFT if the xtb program is used as a driver for a quantum-
obtained with the r2 SCAN-3c and PBEh-3c composite methods. chemistry program like ORCA (see xtb documentation). 129
A very robust but also significantly more expensive choice is Note that in general, a presence of a few low-energy imaginary
PBE0-D4/TZ, which can already be regarded as a benchmark frequencies <50-100 cm-1 is quite typical for large systems and
level for structural optimizations. B3LYP-D4/TZ achieves a sim- often technically/numerically related to the DFT grid. This is per
ilar level of accuracy but is – in our experience – not quite as ro- se not a problem if properly dealt with, that is, the frequencies
bust as PBE0-D4/TZ, in particular for transition-metal containing should be inverted (multiplied by −i) and treated as normal real
systems. Also range-separated hybrids and double-hybrid func- frequencies, like in the xtb program. 125 However, if these small
tionals have been shown to provide very accurate structures. 134 imaginary modes are not inverted or removed, or too many are
However, application of double-hybrids is typically excessive and present, they can significantly bias the thermostatistical entropy
not necessary for standard applications due to beneficial error- corrections. In contrast, the appearance of high-energy imaginary
cancellation effects for structures. Even for very large systems, modes (> 100 /cm) in frequency calculations for fully optimized
structural optimizations should not go below polarized DZ ba- structures indicates that the given structure is no minimum on
sis set level. The functional should be at least of good (m)GGA the respective potential energy surface and requires further re-
quality (PBE or TPSS) and always include a dispersion correc- finement, e.g., by manual distortion and re-optimization.

11
 Only if the individual vibrational frequencies and IR/Raman in cases where properties are weighted for several conformers,
intensities are desired and not just thermostatistical corrections, which can vary widely (see example 3.3).
it might be useful to move to hybrid level and larger basis sets as Thus, conformational energies should at least be obtained with
this improves the calculated spectral intensities. One functional TZ basis sets, and for this reason, the PBEh-3c composite method
with a proven track record for the computation of IR spectra is (based on a DZ basis) is not recommended for the task. Con-
B3LYP, or its dispersion corrected low-cost variant B3LYP-3c, see formational energies are particularly sensitive to mid- and long-
Refs. 7,17. For a comprehensive analysis of the performance of range electron-correlation effects, such that the dispersion cor-
various functionals, basis sets, and the influence of Fock exchange rection takes an important role. In particular, for metal-organic
on vibrational frequencies, we refer to the work of Radom and systems, density or charge-dependent corrections like VV10 or D4
coworkers. 141 should be preferred over the charge-independent D3 scheme. Ac-
cording to the conformational subsets of the GMTKN55 bench-
2.3 General Remarks on Energy Calculations mark, r2 SCAN-D4 and the r2 SCAN-3c and B97-3c composite
For the following four categories, conformers, reaction ener- methods are particularly well-suited for conformation energies.
gies, barriers, and noncovalent interactions, the calculated If the systems are small and higher-level calculations are afford-
electronic energies become the central quantity, which requires able; hybrid functionals can be employed with large TZ basis sets.
larger basis sets for converged results than structures or frequen- Here, we recommend the ωB97X-V and ωB97M-V approaches of
cies. Therefore, DZ basis sets (like 6-31G** or def2-SVP) are no Mardirossian and Head-Gordon or the respective D4 or D3 ana-
longer sufficient, and we strongly advise against using them ex- logues. 145,146 Also, a common dispersion-corrected B3LYP/TZ or
cept if they are part of composite schemes made for this pur- QZ approach can provide accurate conformational energies, yet
pose. Even in combination with full counter-poise corrections this method has been outperformed in a recent benchmark by the
or gCP, the residual BSSE and BSIE of DZ basis sets is substan- r2 SCAN-3c composite method at a small fraction of the computa-
tial. One possible exemption is the calculation of NCIs with the tional cost (see also example 3.3). Double hybrid functionals can
def2-SVPD basis set and tailor-made DFT-C correction for BSSE be used for small systems and maximum robustness and accuracy.
of Head-Gordon and Witte, 16 which we showcase for the NCI ex- The most accurate functionals in the GMTKN55 benchmark are
ample in Section 3.2. Common TZ basis sets often yield results ωB97M(2) 147 of Mardirossian and Head-Gordon and revDSD-
already reasonably close to the basis set limit, but convergence PBEP86-D4 148 of Martin and coworkers. A particularly robust
should be checked at a QZ level in critical cases. Before reduc- and widely available double hybrid we also want to recommend
ing the basis set size to the bare minimum, we suggest moving to is PWPB95-D4 149 . On a side note, we want to mention that
more efficient (m)GGA functionals instead of hybrids or even TZ- semi-empirical quantum mechanics and force-field-based meth-
based composite methods, which are purpose-made to perform ods yield much less reliable conformational energies than DFT
well with smaller basis sets. and hence can only be applied in initial steps of multi-level work-
For double-hybrid functionals, the limitations of TZ basis are flows and by using a very conservative (large) energy selection
more severe since the MP2 component of the calculation typically window, and are best coupled with DFT-based energy re-ranking
has a stronger basis-set dependency than the (m)GGA or hybrid (see workflow given in Fig. 7). 52
DFT part. Here, complete-basis-set-extrapolation (CBS) from TZ
to QZ should be considered 142,143 , but the details are beyond 2.5 Reaction Energies
the scope of this work. Lastly, we suggest deciding on one ap- Reaction energies refer to the difference in the total electronic
proach for all energy-related properties to retain a certain degree energies between reactants, products, and possible intermediates
of consistency and comparability, meaning that reaction energies, that constitute minima on the potential energy surface. Since the
barrier heights, and association energies of reactants in one reac- geometric and electronic structure of the molecules differs more
tion or reaction network should be obtained consistently with one widely, error compensation is weaker than for conformational en-
method combination, e.g., ωB97M-V/QZ 144 single-point energies ergies. However, the practically acceptable error for reaction en-
on r2 SCAN-3c structures (ωB97M-V/QZ//r2 SCAN-3c) Therefore, ergies is typically also larger than for conformational energies,
a compromise needs to be made considering all the requirements and an accuracy of about 1-2 kcal·mol−1 , which is difficult to ob-
stated in the following. tain experimentally, is often considered to be sufficient. While
all of this strongly depends on the reaction under investigation, a
2.4 Conformers number of general findings hold true for the vast majority of ex-
Conformational energies refer to the differences in the electronic amples: Regarding the basis set, reaction energies require large
energy of different conformers of a given molecule with a fixed TZ or QZ basis sets for converged results (see examples). The
covalent bond topology. Conformational energies are more for- basis-set dependency should be carefully investigated by compar-
giving in the energy-related categories due to the typically similar ing single-point energies obtained with the next-smaller basis set
structures as they strongly profit from error compensation. At (e.g., def2-QZVP to def2-TZVP). Even if computational resources
the same time, however, they need to be accurate within a range are at the limit, the basis set size should not be reduced below TZ
of about 0.1-0.2 kcal·mol−1 to predict Boltzmann populations at quality.
room temperature reasonably well. This is particularly important Instead of reducing the basis set size, it should always be

12
considered to switch to semi-local (m)GGA functionals if SIE is heavy-atom peri-cyclic reactions, and almost absent for covalent,
not a major concern. Surprisingly, one of the most efficient ap- bond-conserving inversion or conformational processes. Hence,
proaches for accurate reaction energies is the r2 SCAN-3c compos- cheap (m)GGA functionals may be used in special cases such as
ite method, which outperforms the recommended B3LYP-D4/QZ conformational or inversion barriers after careful testing but, in
and PW6B95/QZ hybrid methods on the reaction energy subset general, barrier height calculations are the only category in which
of the GMTKN55 benchmark (see Fig. 4 in Ref. 18). we advise against using semi-local (m)GGA functionals. Since
In general, however, the accuracy and, in particular, the ro- the calculation of barrier heights is perhaps the most challenging
bustness of predicted reaction energies of non-metallic systems task, we provide an extended discussion of the aspects mentioned
profits from an admixture of Fock exchange. The optimum value above in the example shown in Section 3.4.
for reaction energies with global hybrids is at ≈ 25% Fock ex- In general, to mitigate the errors related to SIE, range-
change. In contrast, barrier-heights often profit from even higher separated hybrids are strongly recommended for barrier heights.
amounts of ≈ 50% depending on the reaction type (see details Specifically, ωB97M-V (also ωB97X-V) are performing very well
in Section 2.6). Therefore, if both reaction energies and barrier on the barrier height subsets of the GMTKN55 and other bench-
heights are to be calculated, preferably with the same functional mark sets. If global hybrids need to be employed for some rea-
as discussed above, a compromise must be made with global hy- son, those with an increased Fock exchange (>30%) should be
brid functionals. In this regard, we recommend Truhlar’s PW6B95 preferred, and the results should be compared against range-
(28% Fock-exchange) 150 hybrid with D3 or D4 dispersion correc- separated hybrids. Global hybrid functionals that have been
tion, which provides accurate and robust thermochemistry, and specifically designed for the prediction of barrier heights, like
has been our default hybrid functional for mechanistic studies for the BMK 151 or MPWB1K 152 functionals, typically use 40-50%
many years. Fock exchange (BMK 42%, MPWB1K 44%). Note that such high
An approach to solve the problematic of high or low Fock- amounts of Fock exchange typically deteriorate the performance
exchange is provided by range-separated functionals such as of reaction energies. Also, double hybrid functionals are well
ωB97X-V and ωB97M-V. The variable admixture of Fock exchange suited for this task since they typically use a much larger fraction
in these functionals makes it possible to get the best of both of Fock exchange (>50%) than global hybrids and nevertheless
worlds: good reaction energies and barrier heights. Accordingly, provide very accurate reaction energies. Due to their challenging
the best-performing hybrid functionals on the GMTKN55 bench- electronic structure, basis set convergence may also be slower for
mark set are the above-mentioned range-separated ones being the transition states/barrier heights than for other properties, and QZ
only hybrids with a WTMAD2 below 4 kcal·mol−1 . Double hybrid basis sets should be considered.
functionals are even more accurate and robust choices but also The London dispersion energy contribution to typical chemical
most computationally demanding because they require larger ba- reaction energies or barriers computed with standard function-
sis sets. als can be large (or even decisive, see Ref. 153 ), especially for
Lastly, to complete this discussion of the optimal amount of molecules with >20-30 atoms, and hence its explicit considera-
Fock exchange, we mention that this value also depends on the tion is very strongly recommended.
type of the studied reaction. While the arguments and val-
ues presented above apply to main-group chemistry, transition- 2.7 Noncovalent Interactions
metal compounds typically require lower amounts (roughly about
Noncovalent Interactions refer to the difference in electronic en-
half of the values given above). Accordingly, functionals with
ergies between a non-covalently interacting complex and its iso-
very high amounts of Fock exchange and range-separated hybrids
lated molecular fragments. Since, by definition, no covalent
should be applied with care. An example is the prominent M06-
bonds are changed upon association of the fragments, noncova-
2X (54%) functional and the PBEh-3c composite method, which
lent interaction (NCI) energies strongly profit from error com-
showed larger errors in recent transition-metal reaction bench-
pensation, even more so than conformational energies. However,
marks. 36,87
at the same time, NCI energies are often relatively small on an
atom pair-wise basis, and therefore the required accuracy is of-
2.6 Barrier Heights ten higher than for reaction energies. Moreover, small NCI en-
Barrier heights refer to the electronic energy difference between ergies can cause BSSE and SIE to become relatively large, and
the transition state and the corresponding reactants/products. their influence should be carefully investigated. Since London
Many transition state structures typically involve at least one dispersion is usually a dominant contribution to the binding in
stretched bond and, in turn, near-degenerate orbitals with weakly NCI complexes, the dispersion correction is particularly impor-
bound electrons, which gives rise to a particularly challenging tant. While VV10 can have a slight edge over D4 in exotic and
electronic structure. As a result, transition states are typically charged systems due to its dependence on the density, the inher-
prone to SIE, which often leads to a systematic underestimation ently better C6 coefficients 14,136,137 as well as the inclusion of
of their electronic energy and, in turn, barrier heights by (m)GGA three-body-terms 154,155 in D4 (and D3-ATM in PBEh-3c), makes
functionals. This error strongly depends on the one-electron char- D3-ATM and D4 more accurate in highly polarizable, small-gap
acter of the breaking bond(s) in the transition state: It is largest systems (e.g., bucky-balls or graphene sheets, cf. the L7 and S30L
for hydrogen-transfer or dissociation reactions, smaller, e.g., for benchmarks 156,157 ).

13
 All "3c" composite methods can be recommended to efficiently tween the close aromatic rings. Accordingly, a large deviation
calculate noncovalent interaction energies, as they have been de- is obtained at the dispersion-uncorrected B3LYP/QZ level, which
signed with this purpose in mind. For most systems, r2 SCAN-3c strongly underestimates the reaction enthalpy (∆H B3LYP/QZ =
is the best choice, while for SIE prone, typically highly polar sys- −25.2 kcal·mol−1 ). Including the D4 dispersion correction results
tems, PBEh-3c can be superior. However, the small basis of PBEh- in a much improved value of ∆H B3LYP-D4/QZ = −43.4 kcal·mol−1 .
3c can be problematic in such cases due to BSSE and BSIE. An- Nevertheless, the deviation from the W1-F12 reference value
other composite approach has been developed by Head-Gordon still amounts to 5.1 kcal·mol−1 . An even better agreement with
and coworkers that combines the accurate B97M-V mGGA func- the reference value is obtained employing the PWPB95-D4/def2-
tional with the def2-SVPD basis set and a tailor-made gCP-derived QZVP double-hybrid. The respective reaction enthalpy is cal-
correction termed DFT-C. 16 This approach provided NCI energies culated to ∆H PWPB95-D4/QZ = −48.5 kcal·mol−1 in near-perfect
with an accuracy comparable to B97M-V/QZ results in their tests. agreement with the reference.
We demonstrate this approach in the NCI example in Section 3.2. The strain-reducing isomerization is less prone to intrinsic func-
Finally, we want to mention that due to the good error tional errors, mainly due to beneficial error compensation re-
compensation for noncovalent interaction energies, even semi- sulting from the chemical similarity of [2.2]paracyclophane and
empirical quantum mechanical methods like GFN2-xTB or PM6- [2.2]metacyclophane. Accordingly, all tested DFT methods are in
D3H4 158,159 can provide reasonable results for interactions in- reasonable agreement and again, the PWPB95-D4/QZ result is in
cluding complex geometries at a fraction of the computation cost perfect agreement with the reference value. In both cases, the
of a DFT calculation. Combining these very "low-cost" structures energy to enthalpy corrections are small compared to the relative
and frequencies with single-point calculations at a composite-DFT electronic energies, only contributing by 4.2 and 0.5 kcal·mol−1 ,
level is perhaps the best way to estimate noncovalent interaction respectively.
energies for very large systems with hundreds or thousands of
atoms. 3.2 Noncovalent interactions
On the high-end of the methodological spectrum, ωB97M-V Noncovalent interactions (NCIs) play an important role in chem-
and ωB97X-V provide exceedingly accurate noncovalent interac- istry, particularly in bio- and supramolecular systems. 165,166 The
tion energies if combined with a large TZ/QZ basis. Also, double theoretical description of supramolecular complexes is therefore
hybrid functionals are very accurate and robust for noncovalent of considerable importance, and particularly the prediction of
interaction energies, but also more basis set dependent and de- binding free energies (∆G). Since chemically relevant systems
manding. Noncovalent interactions in or with very small gap are often quite large and also flexible, this task is challenging
systems (metals) are generally not well understood and require for computational chemistry. To compare calculated values and
special treatment. 160 experimental data such as binding free energies measured in so-
lution, the thermostatistical corrections must also include the en-
3 Examples tropy terms. This is specifically the case for bimolecular reactions.
3.1 Formation and Isomerization of [2.2]Paracyclophane The common approach to calculating binding free energies for the
formation of a NCI complex at a given temperature is shown in
The dimerization of 3,6-dimethylidenecyclohexa-1,4-diene to
Eq. (1) 167
[2.2]paracyclophane and its subsequent isomerization to
∆G = ∆E + ∆δ Gsolv + ∆GmRRHO , (1)
[2.2]metacyclophane shown in Fig. 6 represents a funda-
mental chemical transformation in organic chemistry. For where ∆E refers to the difference of the total electronic gas-phase
both reactions, experimentally determined standard heats of energies, ∆δ Gsolv to the difference in solvation free energies, and
formation can be used to derive the reaction enthalpies ∆GmRRHO corresponds to the difference in the thermostatistical
(∆H exptl. ). 161–163 Nevertheless, due to high uncertainties for the contributions. Depending on details of the complex in question
3,6-dimethylidenecyclohexa-1,4-diene monomer, the reaction en- (e.g., charged vs. neutral, H-bonds vs. π-π stacking, solvent),
thalpy of the dimerization amounts to −41.9±10.6 kcal·mol−1 . these individual contributions may vary significantly in size and
This error estimate is much smaller for the isomerization for can be of different sign. Regardless, all three contributions to ∆G
which an reaction enthalpy of −17.9±2.9 kcal·mol−1 is derived. must be described fairly precisely, as these are typically large and
To obtain more reliable reference values, highly accurate W1- only partially cancel each other to give the typically rather small
F12 164 electronic energies were calculated and combined with experimental ∆Gs (−1 to −15 kcal·mol−1 ). 157
PBE0-D4/def2-TZVP zero-point vibrational energy and enthalpy Recently, we have proposed a workflow 49 to compute the con-
corrections. Geometries calculated on the same DFT level were tributions to Equation (1), which is summarized in Fig. 7 and
used throughout (Fig. 6). The calculated reference reaction en- showcased for this example. The example complex is a +4
thalpies amount to −48.7 and −19.0 kcal·mol−1 , and thus lie charged macrocyclic host (termed CBPQT4+) with bromoben-
within the large error bars of the experimental values. Specifi- zene as guest, whose ∆G was experimentally determined to be
cally, the dimerization reaction represents a challenging task for −5 kcal·mol−1 in aqueous solution. 168 The first challenge is to
computational methods as two new strained single bonds are model the actual molecular structure as realistically as possible,
formed and the product includes pronounced intramolecular Lon- ideally in solution at finite temperature. To this end, a con-
don dispersion as well as exchange-correlation interactions be- former search using the automatic CREST and CENSO approaches

14
 dimerization isomerization
ΔHexptl. = −42 ± 11 ΔHexptl. = −18 ± 3
2
ΔHW1-F12 = −48.7 ± 0.5 ΔHW1-F12 = −19.0 ± 0.5
ΔHB3LYP/QZ = −25.2 ΔHB3LYP/QZ = −20.7
ΔHB3LYP-D4/QZ = −43.4 ΔHB3LYP-D4/QZ = −17.4
ΔHPWPB95-D4/QZ = −48.5 ΔHPWPB95-D4/QZ = −19.1

[2.2]paracyclophane [2.2]metacyclophane

Figure 6 Formation of [2.2]paracyclophane 161 from 3,6-dimethylidenecyclohexa-1,4-diene 162 and subsequent isomerization to [2.2]metacyclo-
phane 163 . All DFT reaction enthalpies employ PBE0-D4/def2-TZVP geometries, zero-point vibrational energy and enthalpy corrections at T = 298.15 K.
All values in kcal·mol−1 . W1-F12 denotes a highly accurate wavefunction-theory-based reference level. QZ = def2-QZVP.

the conformational entropy should be explicitly considered, 52,169

4+ automatic CREST which significantly increases the cost of the calculations.
(GFN2-xTB(ALPB))
For the final geometry optimization (and conveniently also the
final energy calculation, see below) of the energetically lowest
N
standard CENSO workflow conformer, we recommend the composite DFT method r2 SCAN-
Br 3c with DCOSMO-RS as implicit solvation model. Alternatively,
part 0: SMD, CPCM, and COSMO with descending preference from first
exclude conformers > 4 kcal·mol−1
B97-D3/def2-SV(P) // GFN2-xTB to last can be used. Note that COSMO-RS cannot be used in ge-
ometry optimizations. For the example discussed here, all four
part 1: mentioned implicit solvation models are suitable for the geome-
exclude conformers > 3.5 kcal·mol−1 try optimization as evident from Fig. 8. Moreover, we note that
r2SCAN-3c + COSMO-RS
+ GmRRHO(GFN2[ALPB]-SPH) also our previous group-default for geometries, a TPSS-D3/def2-
TZVP optimization in gas phase, gives accurate results for this
part 2: example and also in general. This is because results are often not
final geometry optimization
very sensitive to the employed structures, as already discussed in
r2SCAN-3c[DCOSMO-RS]
bromobenzene@CBPQT(4+) Section 2.1.
final ΔG ranking The thermostatistical corrections ∆GmRRHO typically provide a
ΔGcalc. = −5.4 kcal·mol−1 E(r2SCAN-3c) + Gsolv(COSMO-RS) significant positive, repulsive contribution to the free binding en-
+ GmRRHO(GFN2-xTB[ALPB]-SPH)
ΔGexptl. = −5.0 ± 0.5 kcal·mol−1 ergy (blue bars in Fig. 8). We calculate them in the modified
rigid-rotor harmonic oscillator (mRRHO) approximation, which
includes a special treatment for low-frequency modes, zero-point
Figure 7 Suggested computational workflow to calculate the binding free
harmonic vibrational energy (ZPVE) and heat/volume work cor-
energy (∆G) of bromobenzene to CBPQT(4+) in water. For details on the
CREST/CENSO workflow see Refs. 51,52, and 49 rections, but neglects the conformational entropy. This can be a
good approximation if the involved molecules are mostly rigid as
discussed above. Our default is to compute vibrational frequen-
is important, particularly for the complex but also for hosts and cies with fast semi-empirical methods such as GFN2-xTB, which
guests if they are flexible molecules. In Ref. 46, we demon- is typically in good agreement with DFT reference values (devi-
strated that the efficient semi-empirical extended tight-binding ation ≤1-2 kcal·mol−1 ). 170 This requires a suitable implicit sol-
method GFN2-xTB with a specially adapted implicit solvation vation model like ALPB as well as the single point Hessian (SPH)
model (ALPB) is well-suited for this purpose. For the complex ansatz, 128 which uses a biasing potential approach to create an
discussed here, twelve conformers within a relative energy range artificial minimum at the DFT structure. The comparison of the
of 0–2.5 kcal·mol−1 are found after the final ∆G ranking even default GFN2-xTB[ALPB]-SPH ∆GmRRHO contributions with cor-
though for rigid molecules like this one, it is often sufficient to responding DFT values obtained in gas-phase at the TPSS-D3/TZ
consider only the energetically lowest conformer. However, for level shown in Fig. 8 confirms that the semi-empirical approach
very flexible host or guest structures, the consideration of further yields practically identical results but in a few minutes vs. several
conformers by means of Boltzmann weighting is strongly recom- hours of computation time for the DFT calculation. The alter-
mended. 49 The reason is that if the configuration space is signif- native to SPH calculations, that is, a full re-optimization of the
icantly narrowed in the complex compared to a free host/guest, geometry with GFN2-xTB (xtb-keyword: -OHESS), yields slightly
e.g., because an alkyl chain is fixed in one position in the complex, larger deviations from the gas-phase DFT reference (see Fig. 8),
this gives rise to a large entropy penalty and, in turn, temperature but the SPH calculations benefit also from error cancellation due
dependency of the binding due to the T · S term. In such cases, to the use of an implicit solvation model. Hence, we generally

15
recommend to calculate the frequencies with an implicit solva- to within 1-3 kcal·mol−1 from experimentally determined val-
tion model and the SPH algorithm. ues, 157 which is much less than the sum of the maximum errors
To calculate the binding energy (∆E) contribution to ∆G, our of the individual contributions would suggest. Also for this ex-
default protocol uses the r2 SCAN-3c single-point energy. In gen- ample, the calculated ∆G value of −5.4 kcal·mol−1 agrees very
eral, composite DFT methods such as r2 SCAN-3c or B97M-V/def2- well with the experimentally determined one for a wide range of
SVPD/DFT-C are efficient alternatives to numerically converged method combinations. Moreover, the presented workflow has the
QZ basis set DFT calculations and approximately 50 times faster advantage to be computationally quite fast, requiring only about
for this system size. As evident from Fig. 8, r2 SCAN-3c (default, 30 hours in total on a common eight-core CPU for the shown
leftmost bar), B3LYP-3c, and B97M-V/DFT-C are all in good agree- example. This efficiency largely results from the use of fast semi-
ment with the reference. However, it should be noted that when empirical methods for the frequencies and efficient mGGA-based
the systems are highly charged or feature exotic chemical interac- composite DFT methods for geometries and single point energies.
tions, it is generally advisable to compare ∆E obtained with com- With such an efficient multi-level scheme, reliable affinity predic-
posite methods to a more robust DFT/QZ calculation with a well- tions are possible for much larger complexes with up to 200-300
performing hybrid (e.g., ωB97M-V) to be on the safe side. The re- atoms in practical computation times. Finally, we want to men-
maining scatter of typically ±2 kcal·mol−1 for ∆E can usually be tion that the protocol shown in Fig. 7 is fully automated via the
attributed in about half to the errors of the density functional and freely available CREST and CENSO programs, and can thus be
in the other half to the errors of the dispersion correction, which is invoked via two simple UNIX commands as described in the doc-
particularly important in NCI examples (see discussion in Section umentation. 172,173
2.7). For this example, VV10/NL is actually somewhat more accu-
rate in combination with B3LYP than D4 (see Fig. 8). Concerning
3.3 Optical rotation of α-/β -D-glucopyranose
the basis set for the DFT calculation, the residual basis set errors
are smaller than the intrinsic functional error with QZ basis sets The calculation of relative electronic energies is fundamental
and often already with large TZ basis sets. However, this changes and can indirectly have a crucial impact on property calculations
when the basis set size is further reduced to DZ, as evident from when ensemble-averaged treatments are used. This is for exam-
the hierarchy of B3LYP-based results summarized in Fig. 8: With ple the case for computed specific optical rotation values of a
the TZ basis the B3LYP-D4 results change only slightly compared mixture of α- and β −D-glucopyranose depicted in Figure 9. 174
to QZ (and also M06-2X provides very good agreement with the The optical rotation strongly depends on details of the molecular
TZ basis), whereas the DZ basis 6-31G* with B3LYP-D3 exhibits structure, and for flexible systems, on a reliable calculation of the
strong BSSE-induced overbinding. In contrast, the still widely relevant populated conformers. Further, the values of α- and β -D-
used B3LYP/6-31G* approach drastically underestimates the ∆E glucopyranose isomers differ strongly from each other, and thus
contribution due to the lack of a dispersion correction, and thus the accurate description of the thermodynamic equilibrium be-
the resulting ∆G is off by more than 20 kcal·mol−1 . This demon- tween both forms is crucial. Here, the impact of the choice of the
strates that the often assumed error compensation between lack- theoretical level for the free energy calculations on the final op-
ing London dispersion and BSSE in B3LYP/6-31G* is not to be tical rotation values calculated at the PBE/def2-SVPD[COSMO]
trusted. To show how the shortcomings of B3LYP/6-31G* can be level is demonstrated. For both anomers, conformer-rotamer
fixed at essentially no extra cost, we included the value obtained ensembles (CRE) were generated employing the CREST/CENSO
with the B3LYP-3c approach, which combines B3LYP-D3 with a approach, and the finally obtained conformers were originally
DZ basis (def2-SVP) and the default gCP correction for BSSE (see ranked regarding their conformational free energies on the
Fig. 8). Evidently, this physically sound approach gives results r2 SCAN-3c[COSMO-RS]//r2 SCAN-3c[DCOSMO-RS] level of the-
very close to B3LYP-D4/QZ at the same cost as B3LYP/6-31G*, ory. The same CREs are here re-ranked employing B3LYP method
demonstrating once again that there is no reason for using this combinations (B3LYP/QZ, B3LYP-D4/QZ, B3LYP/6-31G*, QZ =
outdated but still popular method. 9 def2-QZVP), as well as BP86 175,176 -D4/QZ for the electronic en-
The most challenging contribution in the entire workflow is the ergy contribution. All others (zero-point vibrational energy, ther-
solvation free energy contribution ∆Gsolv , especially for higher or mostatistical, and solvation corrections) were taken from the
negatively charged systems in polar, H-bonding solvents. Here, original r2 SCAN-3c calculation. In agreement with results from
an error range of 2-3 kcal·mol−1 is realistic due to a large ∆δ Gsolv benchmark studies, the r2 SCAN-3c composite method yields very
value, which has an estimated intrinsic error of 10-20%, even good energetic rankings, and the Boltzmann-weighted optical ro-
with the best implicit solvation models available. Among them, tation value for the α- and β -D-glucopyranose equilibrium is in
according to our experience, COSMO-RS yields the most reliable excellent agreement with the experimental value of α exptl. = 52.7
results, while SMD can serve as a good alternative. In contrast, (all values in the following in the usual degree dm(g/cm3 )−1
purely electrostatic models like COSMO and CPCM often perform units). The same holds for the value that is based on a B3LYP-
worse because they neglect all non-electrostatic terms, which are D4/QZ ranking amounting to α B3LYP-D4/QZ = 54.6. For the plain
particularly important if the solvent accessible molecular surfaces B3LYP/QZ level, a worse result of α B3LYP/QZ = 41.2 is obtained,
change during the reaction (as is the case here, see Fig. 8). underlining the indispensability of a London dispersion correc-
Due to fortuitous error compensation effects, the free binding tion even for relatively small molecules. The frequently used
energy produced by our standard workflow is typically accurate B3LYP/6-31G* approach relies on a difficult-to-control error com-

16
 method for
solvation model implicit solvation
change from default: frequency DFT method for ΔE calculation
in geo.-opt. for ΔδGsolv.
calculation in ΔGmRRHO

20 ΔE ΔδGsolv. ΔGmRRHO ΔGcalc.

15 ✗

10
size of contribution / kcal·mol−1

0
✗
ΔGexptl.
-5 ✓ ✓ ✓ ✓ ✓ ✓ ✓ ✓
✓ ✓
✓ ✗ ✓
✗
-10

-15 ✗

-20
−31.3

-25
t

M
O

3c

TZ
S)

TZ

TZ

*
ul

-V
G

G
Q

Q
SM

SM
SM
PC

PC

ES
fa

P-
31

31

X/
3/

4/

7M
4/

L/

4/
de

-D

LY

-2
C

-D

D
O

6-

6-
H

B9
P-
P-

P-
C

06
O

SS

3/

P/

B3
95

LY
B(

LY

LY

LY

M
PB
TP

P-
xT

B3
B3

B3

B3
PW

LY
2-
FN

B3
G

Figure 8 Breakdown of the contributions to the overall calculated ∆Gcalc. with the final result represented by a black line, and experimental reference
given as grey dashed line with the shaded areas marking the estimated error range of ±1.0 kcal·mol−1 and a more conservative error estimate of
±2.5 kcal·mol−1 . Methods that do not reach the ±1.0 kcal·mol−1 window but are within ±2.5 kcal·mol−1 of the experimental reference value are
marked by yellow ticks. The leftmost bar represents the default approach (∆E(r2 SCAN-3c) + ∆δ Gsolv (COSMO-RS) + ∆GmRRHO (GFN2-xTB[ALPB]-
SPH) at r2 SCAN-3c[DCOSMO-RS] geometry), while the other bars illustrate the effect of selected method variations as indicated at the top (solvation,
frequencies, electronic energy). Contributions that are not affects by these variations are depicted in brighter colors. B97M-V* = B97M-V/def2-
SVPD/DFT-C. B3LYP-NL utilizes VV10 with refitted parameters. 171 OHESS = GFN2-xTB[ALPB] frequencies with a re-optmized structure instead of
SPH with the DFT structure.

pensation of neglected London dispersion and BSSE. Neverthe- types in organic chemistry. Even though, these basic reaction
less, the contribution of α-D-glucopyranose is overestimated, re- types seem comparably simple, their theoretical description still
sulting in a relatively bad value of α B3LYP/6-31G* = 61.7. The requires a profound choice of the applied quantum chemical
pure GGA functional BP86-D4/QZ, which in general does not method. For the [2+4] cycloaddition of ethylene to cyclopen-
perform well for conformational energies, yields an even worse tadiene (Fig. 10a) and the nucleophilic attack of OH− to fluo-
agreement with the experiment (α BP86-D4/QZ = 72.2). This ex- romethane (Fig. 10b), accurate reference electronic activation en-
ample clearly demonstrates that the best results are obtained if all ergies (∆E ‡ ) are available. 23 For both examples, the BP86-D3/QZ
physically relevant effects, such as London dispersion and basis- (QZ = def2-QZVP) GGA strongly underestimates the reaction bar-
set completeness, are properly taken into account. Furthermore, rier. Specifically, the reaction barrier of the SN 2 reaction is under-
conformation-sensitive properties like optical rotation can indi- estimated by 10.2 kcal·mol−1 . Here, the pronounced charge de-
rectly be used to assess the quality of theoretical approximations. localization in the transition state causes SIE related issues with
the GGA method. The mGGA M06-L/QZ yields much improved
results, which may be attributed to its highly empirical charac-
3.4 Reaction barriers
ter with the functional being also trained on reproducing barrier
3.4.1 SN 2 and Diels–Alder reaction heights in similar systems. Nevertheless, hybrid functionals such
as PBE0 are expected to yield improved results over the (m)GGA
The reliable calculation of reaction barriers is crucial for the inves-
methods by an enhanced physical description and for both reac-
tigation of complex reaction mechanisms. They allow for a deeper
tions, PBE0-D3/QZ yields a significant improvement over BP86-
understanding of key reactions and thus a computer-aided design
D3/QZ. Finally, the range-separated hybrid ωB97X-V/QZ further
of novel catalysts and a targeted tuning of chemical reactions. 177
improves the results compared with the global hybrid and the
SN 2 and Diels–Alder reactions represent well-known reaction

17
 𝛼-D-glucopyranose 𝛽-D-glucopyranose
a
‡
H OH H OH

H O H O
HO
H
HO
OH
+
HO HO
H OH H OH
80
specific optical rotation / °[dm(g/cm3)]−1

H OH H H

✗
70 20
✓ ✓

ΔE‡ / kcal·mol−1
✗ ΔE‡ref. = 18.0 kcal·mol−1
60 15
✓ ✓ ✗
50 𝛼exptl. = 52.7
10 ✗
40 ✗

5
30

20 0
b
10 ‡−
OH− F−
O F
0 + +
CH3F CH3OH
c
Z

-3
G
Q

/Q

20
31

AN
4/

4/
P
-D

D
6-
LY

C
-
P/
P

86
B3

r S2

ΔE‡ / kcal·mol−1
LY

LY

ΔE‡ref. = 17.6 kcal·mol−1

✓
BP
B3

15
B3

✗
✗
Figure 9 Calculated, Boltzmann-weighted optical rotation for the equilib-
10
rium of α- and β -D-glucopyranose in water at 20°C. For each anomer,
15-18 conformers are considered in the example. QZ = def2-QZVP. ✗
5

GGA, underestimating the SN 2 reference reaction barrier only 0

slightly by −1.8 kcal·mol−1 . It also yields comparably good agree-
3

3
-L

-V
-D

-D
06

7X
86

E0
ment for the Diels–Alder activation reaction overestimating the
M

B9
BP

PB

ω
reference value by 1.8 kcal·mol−1 .

3.4.2 Hydrogen atom transfer reaction Figure 10 Calculated gas-phase electronic reaction barriers (∆E ‡ ) for se-
lected functionals using the def2-QZVP basis set 23 for a) the Diels–Alder
Another challenging example reaction is depicted in Fig. 11. cycloaddition of cyclopentadiene and ethylene (reference level: accurate
Here, a hydrogen atom is transferred from a molybdenum hy- W1-F12 wavefunction theory), and b) the SN 2 reaction of fluoromethane
dride complex TpMo(CO)3 H (Tp = tris(pyrazolyl)borate) to with a hydroxide anion (reference level: very accurate W2-F12 wavefunc-
a Gomberg-type (tBu-4-C6 H4 )3 C• radical. 178 The experimental tion theory). The light-grey area indicates the expected error-range of the
reference method.
free energy activation barrier was determined to ∆G‡ = 19.2
kcal·mol−1 . As discussed in the previous sections and the elec-
tronic activation energy examples, transition states with stretched the transition state does not represent a stationary point at the
covalent bonds are prone to SIE, and thus a significant underes- corresponding potential energy surface. B3LYP-D4/QZ yields a
timation of the reaction barrier by (m)GGA methods is expected. small, yet reasonable ∆E ‡ value of 3.0 kcal·mol−1 . The small
This holds especially for bonds involving hydrogen where the one- electronic activation energy further underlines the importance of
electron character is usually large. Compared with the previous solvation, enthalpy and free energy corrections for a direct com-
example, the desired property is an activation free energy in solu- parison to the experiment. The free activation barrier is only well-
tion, involving additionally ZPVE, thermal, solvation and entropy described upon inclusion of all relevant contributions (corrections
effects. For this reaction, the BP86-D4/QZ functional used for computed at the geometry optimization level B97-3c, solvation
the electronic energy underestimates ∆G‡ by almost 4 kcal·mol−1 . corrections by COSMO-RS) leading to an increase of the barrier
This is even more pronounced upon decreasing the basis set size from 3.0 (∆E ‡ ) to 7.2 (∆H ‡ ) and finally 20.6 kcal·mol−1 (∆G‡ ).
from QZ to TZ to DZ with a severe underestimation by over
7 kcal·mol−1 for BP86-D4/DZ. With the B3LYP-D4 hybrid func- 4 Perspective
tional, ∆G‡ is slightly overestimated by 1.4 kcal·mol−1 yet close The development of quantum chemistry over the last 20-30 years,
to the chemical accuracy window of 1 kcal·mol−1 . Nevertheless, and foremost DFT, is a great success story. The fact that nowa-
also for hybrid functionals small basis sets result in an underes- days, non-experts can do reasonable quantum chemistry calcula-
timation of the reaction barrier which seems to be rather sensi- tions for large chemically relevant systems on desktop computers
tive here. The GGA functional even yields a negative electronic is fantastic. The influence of having widely available computa-
activation energy of ∆E ‡ = −2.2 kcal·mol−1 (BP86-D4/QZ) and tional tools on chemical research can not be understated, and

18
 ΔE‡ ΔG‡mRRHO, solv. ΔH‡mRRHO ∑ However, also on the electronic structure side of the problem,
there are challenges ahead. Although many very relevant chem-
‡•
ical properties and problems nowadays can be solved by stan-
dard DFT treatments as described here, there are still problem-
TpMo(CO)3H Mo TpMo(CO)3•
H
+ + atic systems (e.g., open-shell transition metal complexes), open
C(tBuPh)3• HC(tBuPh)3 questions (e.g., how to treat strongly solvated, highly charged
systems?) and problems which are fundamentally difficult (en-
ΔE‡ / ΔH‡ / ΔG‡ / kcal·mol−1

25 tropy). For those aspects, non-standard treatments and expert

knowledge are often required and should be involved.
20 ✓ ✓ ΔG‡exptl. = 19.2±1.5 kcal·mol−1
In our opinion, the positive development of the basic den-
✗ sity functionals over the last decade has slowed down and en-
15 ✗ ✗
tered a kind of saturation regime in terms of the "peak" accu-
✗ racy achieved. Nevertheless, we note some promising new de-
10 velopments which are already applicable such as local hybrid
functionals 179–181 or non-self-consistent field treatments with
5 cheap (m)GGAs to avoid exchange-correlation potential driven
errors 182,183 . Currently, machine-learned functionals 184–186 are
0 still in their infancy, but the potential of such extremely empiri-
cal (and in no way "cheap") methods could be high, at least for
-5 organic and main group compounds. In turn, DFT is the method
of choice for generating extremely large data bases with millions
of compound entries for machine learning algorithms 187–189 . In
Z

Z
TZ

TZ

Z
D

D
Q

Q
4/

4/
4/

4/
4/

4/
D

-D
D

-D

this context, we still see a great potential to make the best per-
D

-D
P-
P-

P-

86
86

86
LY
LY

LY

BP
BP

forming functionals (which are sufficiently accurate for about

BP
B3
B3

B3

>90% of all chemical applications) significantly faster but with-

Figure 11 Calculated contributions to the activation free energy ∆G‡ for out losing robustness or numerical accuracy. Central here are
a metal-centered hydrogen-transfer reaction based on hybrid (B3LYP- accurate yet efficient approximations of the non-local (Fock) ex-
D4) and GGA (BP86-D4) electronic energies (grey bars), B97-3c ther-
mostatistical contributions (yellow bars), and COSMO-RS solvation and
change part of hybrid or double hybrid density functionals and,
entropic corrections (blue bars). Black lines represent the sum of all con- for lower-rungs (mGGA), the numerical integration of the semi-
tributions yielding ∆G‡ . DZ = def2-SVP; TZ = def2-TZVP; QZ = def2- local exchange-correlation energy. Because DFT is a rather gen-
QZVPP basis set. The light-grey area indicates the expected error-range eral "first-principles" approach, we believe that most of our con-
of the reference method.
clusions also hold true for computations of (molecular) solids and
liquids under periodic boundary conditions, although we do not
explicitly considered these here. For more special cases like low-
will presumably grow even larger in the future. We hope that the bandgap systems, such as metallic solids, however, things may
guidelines and recommendations given in this work help to in- change more drastically, such that a consideration of special tech-
crease the reliability of DFT-based quantum chemistry predictions nical settings may be required to avoid a fundamental breakdown
in the daily work of many chemists. of approximations.
We strongly emphasized the aspect of finding the right method-
ological compromise between computational effort (speed) and
Conflicts of interest
desired accuracy while still keeping as close as possible "the right There are no conflicts to declare.
answer for the right reason". An often overlooked aspect in this
balance is that faster theoretical methods enable a more exten-
Acknowledgements
sive – and hence more reliable – study of the system under in- The German Science Foundation (DFG) is gratefully acknowl-
vestigation, for example regarding its conformational behavior, edged for financial support through a Gottfried Wilhelm Leib-
molecular dynamics, or explicit solvation issues. This is of par- niz prize to S.G. and the SPP 2363 ”Utilization and Development
ticular importance since, in our experience, errors and deviations of Machine Learning for Molecular Applications – Molecular Ma-
due to the neglect of important low-lying conformers (ensemble chine Learning”. Further, S.G. and M.B. gratefully acknowledge fi-
properties vs. individual molecule property) can be even larger nancial support of the Max Planck Society through the Max Planck
than the errors in the electronic energy by the functional or ba- fellow program. The authors further thank the ORCA, TURBO-
sis set approximations. We thus want to motivate the reader to MOLE (especially Uwe Huniar), and AMS developer teams for
conduct systematic conformational searches, explore the dynam- steady support as well as Johannes Gorges for his help regard-
ical behavior by means of MD simulations, consider explicit sol- ing some calculations for this paper. We also thank Uwe Hohm,
vation treatments more routinely with currently developed meth- Tunga Salthammer, Wolf Palm, Sigrid Peyerimhoff, Marcel Bam-
ods, and to employ efficient multi-level approaches for chemically berg, Matthias Wagner, Demyan Prokopchuk, and Julia Bursch-
realistic models. Pöppinghaus for comments and fruitful discussions. Further, the

19
authors thank the whole Grimme group for providing a cornu- J. Chem. Phys. 2015, 143, 054107.
copia of diverse experiences and for proofreading. 26 Caldeweyher, E.; Brandenburg, J. G. J. Phys. Condens. Mat-
ter 2018, 30, 213001.
Notes and references
27 Harju, A.; Räsänen, E.; Saarikoski, H.; Puska, M. J.; Niemi-
1 Adamo, C. et al. “DFT Exchange: Sharing Perspectives on the nen, R. M.; Niemelä, K. Phys. Rev. B - Condens. Matter Mater.
Workhorse of Quantum Chemistry and Materials Science”, Phys. 2004, 69, 153101.
2022 , to be submitted to ChemRxiv. 28 Noodleman, L. J. Chem. Phys. 1998, 74, 5737.
2 Houk, K. N.; Liu, F. Acc. Chem. Res. 2017, 50, 539–543. 29 Janssen, C. L.; Nielsen, I. M. Chem. Phys. Lett. 1998, 290,
3 Grimme, S.; Schreiner, P. R. Angew. Chem. Int. Ed. 2018, 57, 423–430.
4170–4176. 30 Lee, T. J.; Taylor, P. R. Int. J. Quantum Chem. 1989, 36,
4 Husch, T.; Vaucher, A. C.; Reiher, M. Int. J. Quantum Chem. 199–207.
2018, 118, e25799. 31 Tishchenko, O.; Zheng, J.; Truhlar, D. G. J. Chem. Theory
5 Bannwarth, C.; Caldeweyher, E.; Ehlert, S.; Hansen, A.; Comput. 2008, 4, 1208–1219.
Pracht, P.; Seibert, J.; Spicher, S.; Grimme, S. Wiley Inter- 32 Karton, A.; Daon, S.; Martin, J. M. Chem. Phys. Lett. 2011,
discip. Rev. Comput. Mol. Sci. 2021, 11, e1493. 510, 165–178.
6 Elstner, M.; Seifert, G. Philos. Trans. R. Soc. A Math. Phys. 33 Langhoff, S. R.; Davidson, E. R. Int. J. Quantum Chem.
Eng. Sci. 2014, 372,. 1974, 8, 61–72.
7 Axel D. Becke, J. Chem. Phys. 1993, 98, 5648–5656. 34 Bauer, C. A.; Hansen, A.; Grimme, S. Chem. Eur. J. 2017,
8 Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; 23, 6150–6164.
Frisch, M. J. J. Phys. Chem. 1994, 98, 11623–11627. 35 Grimme, S.; Hansen, A. Angew. Chem. Int. Ed. 2015, 54,
9 Castro-Alvarez, A.; Carneros, H.; Sánchez, D.; Vilarrasa, J. 12308–12313.
J. Org. Chem. 2015, 80, 11977–11985. 36 Maurer, L. R.; Bursch, M.; Grimme, S.; Hansen, A. J. Chem.
10 Schreiner, P. R.; Fokin, A. A.; Pascal, R. A.; De Meijere, A. Theory Comput. 2021, 17, 6134–6151.
Org. Lett. 2006, 8, 3635–3638. 37 Neese, F. J. Biol. Inorg. Chem. 2006, 11, 702–711.
11 Wodrich, M. D.; Corminboeuf, C.; Schleyer, P. V. R. Org. 38 Barone, V.; Cossi, M. J. Phys. Chem. A 1998, 102, 1995–
Lett. 2006, 8, 3631–3634. 2001.
12 Mardirossian, N.; Head-Gordon, M. Mol. Phys. 2017, 115, 39 Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem.
2315–2372. B 2009, 113, 6378–6396.
13 Verma, P.; Truhlar, D. G. Trends Chem. 2020, 2, 302–318. 40 Klamt, A.; Schüürmann, G. J. Chem. Soc., Perkin Trans. 2
14 Grimme, S.; Hansen, A.; Brandenburg, J. G.; Bannwarth, C. 1993, 0, 799–805.
Chem. Rev. 2016, 116, 5105–5154. 41 Klamt, A. J. Phys. Chem. 1995, 99, 2224–2235.
15 Kruse, H.; Grimme, S. J. Chem. Phys. 2012, 136, 154101. 42 Klamt, A.; Diedenhofen, M. J. Phys. Chem. A 2015, 119,
16 Witte, J.; Neaton, J. B.; Head-Gordon, M. J. Chem. Phys. 5439–5445.
2017, 146, 234105. 43 Clark Still, W.; Tempczyk, A.; Hawley, R. C.; Hendrick-
17 Pracht, P.; Grant, D. F.; Grimme, S. J. Chem. Theory Comput. son, T. J. Am. Chem. Soc. 1990, 112, 6127–6129.
2020, 16, 7044–7060. 44 Lange, A. W.; Herbert, J. M. J. Chem. Theory Comput. 2012,
18 Grimme, S.; Hansen, A.; Ehlert, S.; Mewes, J. M. J. Chem. 8, 1999–2011.
Phys. 2021, 154, 064103. 45 Sigalov, G.; Fenley, A.; Onufriev, A. J. Chem. Phys. 2006,
19 Van Santen, J. A.; Dilabio, G. A. J. Phys. Chem. A 2015, 119, 124, 124902.
6703–6713. 46 Ehlert, S.; Stahn, M.; Spicher, S.; Grimme, S. J. Chem.
20 Mao, Y.; Loipersberger, M.; Horn, P. R.; Das, A.; De- Theory Comput. 2021, 17, 4250–4261.
merdash, O.; Levine, D. S.; Prasad Veccham, S.; Head- 47 Sicinska, D.; Paneth, P.; Truhlar, D. G. J. Phys. Chem. B
Gordon, T.; Head-Gordon, M. Ann. Rev. Phys. Chem. 2021, 2002, 106, 2708–2713.
72, 641-666 PMID: 33636998.
48 Sure, R.; el Mahdali, M.; Plajer, A.; Deglmann, P. J. Comput.
21 von Hopffgarten, M.; Frenking, G. Wiley Interdiscip. Rev. Aided. Mol. Des. 2021, 35, 473–492.
Comput. Mol. Sci. 2012, 2, 43–62.
49 Grimme, S.; Bohle, F.; Hansen, A.; Pracht, P.; Spicher, S.;
22 Glendening, E. D.; Landis, C. R.; Weinhold, F. Wiley Inter- Stahn, M. J. Phys. Chem. A 2021, 125, 4039–4054.
discip. Rev. Comput. Mol. Sci. 2012, 2, 1–42.
50 Bursch, M.; Hansen, A.; Pracht, P.; Kohn, J. T.; Grimme, S.
23 Goerigk, L.; Hansen, A.; Bauer, C.; Ehrlich, S.; Najibi, A.; Phys. Chem. Chem. Phys. 2021, 23, 287–299.
Grimme, S. Phys. Chem. Chem. Phys. 2017, 19, 32184–
51 Grimme, S. J. Chem. Theory Comput. 2019, 15, 2847–2862.
32215.
52 Grimme, S.; Bohle, F.; Hansen, A.; Pracht, P.; Spicher, S.;
24 Brandenburg, J. G.; Bannwarth, C.; Hansen, A.; Grimme, S.
Stahn, M. J. Phys. Chem. A 2021, 125, 4039–4054.
J. Chem. Phys. 2018, 148, 064104.
53 Sobez, J. G.; Reiher, M. J. Chem. Inf. Model. 2020, 60, 3884–
25 Grimme, S.; Brandenburg, J. G.; Bannwarth, C.; Hansen, A.

20
 3900. 2810–2817.
54 Watts, K. S.; Dalal, P.; Murphy, R. B.; Sherman, W.; Fries- 82 Peverati, R.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2012,
ner, R. A.; Shelley, J. C. J. Chem. Inf. Model. 2010, 50, 534– 14, 13171–13174.
546. 83 Yu, H. S.; He, X.; Li, S. L.; Truhlar, D. G. Chem. Sci. 2016,
55 Sauton, N.; Lagorce, D.; Villoutreix, B. O.; Miteva, M. A. 7, 5032–5051.
BMC Bioinformatics 2008, 9, 1–12. 84 Ess, D. H.; Liu, S.; De Proft, F. J. Phys. Chem. A 2010, 114,
56 Miteva, M. A.; Guyon, F.; Tufféry, P. Nucleic Acids Res. 2010, 12952–12957.
38, W622–W627. 85 Mardirossian, N.; Head-Gordon, M. J. Chem. Theory Com-
57 Riniker, S.; Landrum, G. A. J. Chem. Inf. Model. 2015, 55, put. 2013, 9, 4453–4461.
2562–2574. 86 Goerigk, L. J. Phys. Chem. Lett. 2015, 6, 3891–3896.
58 Vainio, M. J.; Johnson, M. S. J. Chem. Inf. Model. 2007, 47, 87 Dohm, S.; Hansen, A.; Steinmetz, M.; Grimme, S.; Checin-
2462–2474. ski, M. P. J. Chem. Theory Comput. 2018, 14, 2596–2608.
59 Hawkins, P. C.; Skillman, A. G.; Warren, G. L.; Elling- 88 Goerigk, L.; Kruse, H.; Grimme, S. ChemPhysChem 2011,
son, B. A.; Stahl, M. T. J. Chem. Inf. Model. 2010, 50, 572– 12, 3421–3433.
584. 89 Jensen, F. Introduction to Computational Chemistry; Wiley,
60 Perdew, J. P.; Schmidt, K. AIP Conf. Proc. 2001, 577, 1–20. 2007, Vol. 2.
61 Perdew, J. P.; Zunger, A. Phys. Rev. B 1981, 23, 5048–5079. 90 Balabin, R. M. J. Chem. Phys. 2008, 129, 164101.
62 Mori-Sánchez, P.; Cohen, A. J.; Yang, W. J. Chem. Phys. 91 van Duijneveldt, F. B.; van de Rijdt, J. G. D.; van
2006, 125, 201102. Lenthe, J. H. Chem. Rev. 1994, 94, 1873–1885.
63 Bao, J. L.; Gagliardi, L.; Truhlar, D. G. J. Phys. Chem. Lett. 92 Prasad, V. K.; Otero-de-la Roza, A.; DiLabio, G. A. J. Chem.
2018, 9, 2353–2358. Theory Comput. 2022, .
64 Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 93 Faver, J. C.; Zheng, Z.; Merz, K. M. Phys. Chem. Chem. Phys.
77, 3865–3868. 2012, 14, 7795–7799.
65 Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158– 94 Jensen, F. J. Chem. Theory Comput. 2010, 6, 100–106.
6170. 95 Ditchfield, R.; Hehre, W. J.; Pople, J. A. J. Chem. Phys.
66 Henderson, T. M.; Janesko, B. G.; Scuseria, G. E. J. Phys. 1971, 54, 720–723.
Chem. A 2008, 112, 12530–12542. 96 Dunning, T. H. J. Chem. Phys. 1989, 90, 1007–1023.
67 Brémond, É.; Pérez-Jiménez, Á. J.; Sancho-García, J. C.; 97 Jensen, F. J. Chem. Phys. 2001, 115, 9113-9125.
Adamo, C. J. Chem. Phys. 2019, 150, 201102. 98 Jensen, F. J. Chem. Phys. 2002, 116, 7372-7379.
68 Mori-Sánchez, P.; Cohen, A. J.; Yang, W. Phys. Rev. Lett. 99 Ahlrichs, R.; Weigend, F.; Häser, M.; Patzelt, H. Chem. Phys.
2008, 100, 146401. Lett 1998, 294, 143–152.
69 Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Science 2008, 321, 100 Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7,
792–794. 3297–3305.
70 Song, S.; Vuckovic, S.; Sim, E.; Burke, K. J. Chem. Theory 101 Weigend, F.; Furche, F.; Ahlrichs, R. J. Chem. Phys. 2003,
Comput. 2022, 18, 817–827. 119, 12753–12762.
71 Grimme, S. J. Chem. Phys. 2006, 124, 34108. 102 Rappoport, D.; Furche, F. J. Chem. Phys. 2010, 133, 134105.
72 Goerigk, L.; Grimme, S. Wiley Interdiscip. Rev. Comput. Mol. 103 Zheng, J.; Xu, X.; Truhlar, D. G. Theor. Chem. Acc. 2011,
Sci. 2014, 4, 576–600. 128, 295–305.
73 Sancho-García, J. C.; Adamo, C. Phys. Chem. Chem. Phys. 104 Andrae, D.; Häußermann, U.; Dolg, M.; Stoll, H.; Preuß, H.
2013, 15, 14581–14594. Theor. Chim. Acta 1990, 77, 123-141.
74 Görling, A.; Levy, M. Phys. Rev. A 1994, 50, 196-204. 105 Kutzelnigg, W.; Liu, W. J. Chem. Phys. 2005, 123, 241102.
75 Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618–622. 106 Saue, T. ChemPhysChem 2011, 12, 3077–3094.
76 Vahtras, O.; Almlöf, J.; Feyereisen, M. W. Chem. Phys. Lett. 107 van Lenthe, E.; Baerends, E. J.; Snijders, J. G. J. Chem. Phys.
1993, 213, 514–518. 1993, 99, 4597–4610.
77 Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Chem. 108 Reiher, M. Theor. Chem. Acc. 2006, 116, 241–252.
Phys. 2009, 356, 98–109. 109 Mewes, J.-M.; Hansen, A.; Grimme, S. Angew. Chem. Int.
78 Hellweg, A.; Huniar, U. “Accelerating Hartree–Fock ex- Ed. 2021, 60, 13144–13149.
change calculation using the TURBOMOLE program sys- 110 Ndongmouo, U. F.; Lee, M. S.; Rousseau, R.; Baletto, F.;
tem: different techniques for different purposes, DOI: Scandolo, S. J. Phys. Chem. A 2007, 111, 12810–12815.
10.48550/ARXIV.1610.07779”, 2016.
111 Mercurio, G.; Maurer, R. J.; Liu, W.; Hagen, S.;
79 Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215– Leyssner, F.; Tegeder, P.; Meyer, J.; Tkatchenko, A.;
241. Soubatch, S.; Reuter, K.; Tautz, F. S. Phys. Rev. B - Condens.
80 Zhao, Y.; Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101. Matter Mater. Phys. 2013, 88, 035421.
81 Peverati, R.; Truhlar, D. G. J. Phys. Chem. Lett. 2011, 2,

21
112 Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density 137 Vydrov, O. A.; Van Voorhis, T. J. Chem. Phys. 2010, 133,
Functional Theory; Wiley-VCH, 2001. 244103.
113 Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry; Dover 138 Becke, A. D. Phys. Rev. A 1988, 38, 3098–3100.
Publications, 1996. 139 Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785–
114 Helgaker, T.; Jørgensen, P.; Olsen, J. Molecular Electronic- 789.
Structure Theory; J. Wiley: New York, 2000. 140 Bannwarth, C.; Ehlert, S.; Grimme, S. J. Chem. Theory Com-
115 Grimme, S.; Hansen, A.; Brandenburg, J. G.; Bannwarth, C. put. 2019, 15, 1652–1671.
Chem. Rev. 2016, 116, 5105–5154. 141 Merrick, J. P.; Moran, D.; Radom, L. J. Phys. Chem. A 2007,
116 Furche, F.; Ahlrichs, R.; Hättig, C.; Klopper, W.; Sierka, M.; 111, 11683–11700.
Weigend, F. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 142 Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.;
91–100. Olsen, J. Chem. Phys. Lett. 1999, 302, 437–446.
117 Neese, F. Wiley Interdiscip. Rev.-Comput. Mol. Sci. 2012, 2, 143 Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.;
73–78. Koch, H.; Olsen, J.; Wilson, A. K. Chem. Phys. Lett. 1998,
118 Neese, F. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2018, 8, 286, 243–252.
e1327. 144 Mardirossian, N.; Head-Gordon, M. J. Chem. Phys. 2016,
119 Neese, F. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2022, 144, 214110.
e1606. 145 Lin, Y. S.; Li, G. D.; Mao, S. P.; Chai, J. D. J. Chem. Theory
120 Epifanovsky, E. et al. J. Chem. Phys. 2021, 155, 084801. Comput. 2013, 9, 263–272.
121 Smith, D. G. et al. J. Chem. Phys. 2020, 152, 184108. 146 Mardirossian, N.; Head-Gordon, M. Phys. Chem. Chem. Phys.
122 Werner, H. J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; 2014, 16, 9904–9924.
Schütz, M. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 147 Chai, J. D.; Head-Gordon, M. J. Chem. Phys. 2009, 131,
242–253. 174105.
123 te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca 148 Santra, G.; Sylvetsky, N.; Martin, J. M. J. Phys. Chem. A
Guerra, C.; van Gisbergen, S. J.; Snijders, J. G.; Ziegler, T. 2019, .
J. Comput. Chem. 2001, 22, 931–967. 149 Goerigk, L.; Grimme, S. J. Chem. Theory Comput. 2011, 7,
124 Frisch, M. J. et al. “Gaussian~16 Revision C.01”, 2016 Gaus- 291–309.
sian Inc. Wallingford CT. 150 Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2005, 109, 5656–
125 Semiempirical Extended Tight-Binding Program Package xtb, 5667.
Version 6.4.1., 2021, https://github.com/grimme-lab/xtb. 151 Boese, A. D.; Martin, J. M. J. Chem. Phys. 2004, 121, 3405–
126 Hourahine, B. et al. J. Chem. Phys. 2020, 152, 124101. 3416.
127 MOPAC2016, James J. P. Stewart, Stewart Com- 152 Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2004, 108, 6908–
putational Chemistry, Colorado Springs, CO, USA, 6918.
HTTP://OpenMOPAC.net (2016). 153 Grimme, S.; Schreiner, P. R. Angew. Chem. Int. Ed. 2011, 50,
128 Spicher, S.; Grimme, S. J. Chem. Theory Comput. 2021, 17, 12639–12642.
1701–1714. 154 Axilrod, B. M.; Teller, E. J. Chem. Phys. 1943, 11, 299–300.
129 User Guide to Semiempirical Tight Binding, https://xtb- 155 Muto, Y. 1943, 17, 629–631.
docs.readthedocs.io/en/latest/contents.html, access-date: 156 Sedlak, R.; Janowski, T.; Pitoňák, M.; Řezáč, J.; Pulay, P.;
2022-04-05. Hobza, P. J. Chem. Theory Comput. 2013, 9, 3364–3374.
130 Bursch, M.; Neugebauer, H.; Ehlert, S.; Grimme, S. J. Chem. 157 Sure, R.; Grimme, S. J. Chem. Theory Comput. 2015, 11,
Phys. 2022, 156, 10–12. 3785–3801.
131 Furness, J. W.; Kaplan, A. D.; Ning, J.; Perdew, J. P.; Sun, J. 158 Stewart, J. J. P. J. Mol. Model. 2007, 13, 1173.
J. Phys. Chem. Lett. 2020, 11, 8208–8215. 159 Řezáč, J.; Hobza, P. J. Chem. Theory Comput. 2012, 8, 141–
132 Ehlert, S.; Huniar, U.; Ning, J.; Furness, J. W.; Sun, J.; 151.
Kaplan, A. D.; Perdew, J. P.; Brandenburg, J. G. J. Chem. 160 Dobson, J. F. Int. J. Quantum Chem. 2014, 114, 1157–1161.
Phys. 2021, 154, 061101. 161 Nishiyama, K.; Sakiyama, N.; Seki, S.; Horita, H.; Ot-
133 Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. subo, T.; Misumi, S. Bull. Chem. Soc. Jpn. 1980, 53, 869–
Rev. Lett. 2003, 91, 146401. 877.
134 Risthaus, T.; Steinmetz, M.; Grimme, S. J. Comput. Chem. 162 Pollack, S. K.; Raine, B. C.; Hehre, W. J. J. Am. Chem. Soc.
2014, 35, 1509–1516. 1981, 103, 6308–6313.
135 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 163 Shieh, C. F.; McNally, D.; Boyd, R. H. Tetrahedron 1969, 25,
2010, 132, 154104. 3653–3665.
136 Caldeweyher, E.; Ehlert, S.; Hansen, A.; Neugebauer, H.; 164 Karton, A.; Martin, J. M. J. Chem. Phys. 2012, 136, 124114.
Spicher, S.; Bannwarth, C.; Grimme, S. J. Chem. Phys. 2019, 165 Kolesnichenko, I. V.; Anslyn, E. V. Chem. Soc. Rev. 2017, 46,
150, 154122.

22
 2385–2390. 178 Vibbert, H. B.; Neugebauer, H.; Norton, J. R.; Hansen, A.;
166 Riley, K. E.; Hobza, P. Wiley Interdiscip. Rev. Comput. Mol. Bursch, M.; Grimme, S. Can. J. Chem. 2021, 99, 216–220.
Sci. 2011, 1, 3–17. 179 Maier, T. M.; Arbuznikov, A. V.; Kaupp, M. Wiley Interdiscip.
167 Grimme, S. Chem. Eur. J. 2012, 18, 9955–9964. Rev. Comput. Mol. Sci. 2019, 9, e1378.
168 Dron, P. I.; Fourmentin, S.; Cazier, F.; Landy, D.; Sur- 180 Bahmann, H.; Kaupp, M. J. Chem. Theory Comput. 2015,
pateanu, G. Supramol. Chem. 2008, 20, 473–477. 11, 1540–1548.
169 Pracht, P.; Grimme, S. Chem. Sci. 2021, 12, 6551–6568. 181 Haasler, M.; Maier, T. M.; Grotjahn, R.; Gückel, S.; Ar-
170 Spicher, S.; Grimme, S. J. Phys. Chem. Lett. 2020, 11, 6606– buznikov, A. V.; Kaupp, M. J. Chem. Theory Comput. 2020,
6611. 16, 5645–5657.
171 Hujo, W.; Grimme, S. J. Chem. Theory Comput. 2011, 7, 182 Sim, E.; Song, S.; Burke, K. J. Phys. Chem. Lett. 2018, 9,
3866–3871. 6385–6392.
172 Introduction to CENSO – xtb documentation, https://xtb- 183 Kim, M. C.; Sim, E.; Burke, K. Phys. Rev. Lett. 2013, 111,
docs.readthedocs.io/en/latest/CENSO_docs/censo_usage.html 073003.
#calculate-free-energies-on-populated-dft-optimized- 184 Kirkpatrick, J. et al. Science 2021, 374, 1385–1389.
conformers, access-date: 2022-04-05. 185 Kalita, B.; Li, L.; McCarty, R. J.; Burke, K. Acc. Chem. Res.
173 Introduction to CREST – xtb documentation, https://xtb- 2021, 54, 818–826.
docs.readthedocs.io/en/latest/crestxmpl.html#imtd-gc, 186 Nagai, R.; Akashi, R.; Sugino, O. npj Comput. Mater. 2020,
access-date: 2022-04-05. 6, 1–8.
174 Bohle, F.; Seibert, J.; Grimme, S. J. Org. Chem. 2021, 86, 187 Schleder, G. R.; Padilha, A. C.; Acosta, C. M.; Costa, M.;
15522–15531. Fazzio, A. JPhys Mater. 2019, 2, 032001.
175 Becke, A. D. Phys. Rev. A 1988, 38, 3098–3100. 188 Belle, C. E.; Aksakalli, V.; Russo, S. P. J. Cheminform. 2021,
176 Perdew, J. P. Phys. Rev. B 1986, 33, 8822–8824. 13, 1–23.
177 Poree, C.; Schoenebeck, F. Acc. Chem. Res. 2017, 50, 605– 189 Blum, L. C.; Reymond, J. L. J. Am. Chem. Soc. 2009, 131,
608. 8732–8733.

23