# Extending density functional theory with near chemical accuracy beyond pure
water

Suhwan Songa,b, Stefan Vuckovicc,d, Youngsam Kima, Hayoung Yua,
Eunji Sima,, and Kieron Burkeb,e
aDepartment of Chemistry, Yonsei University, 50 Yonsei-ro Seodaemun-gu, Seoul
03722, Korea
bDepartment of Chemistry, University of California, Irvine, CA 92697, USA
cInstitute for Microelectronics and Microsystems (CNR-IMM), Via
Monteroni,Campus Unisalento, 73100 Lecce, Italy
dDepartment of Chemistry&Pharmaceutical Sciences and Amsterdam Institute of
Molecular and Life Sciences (AIMMS),
Faculty of Science, Vrije Universiteit, De Boelelaan 1083, 1081HV Amsterdam,
The Netherlands
eDepartments of Physics & Astronomy, University of California, Irvine, CA
92697, USA
<EMAIL_ADDRESS>
## I. Introduction

The overwhelming importance of simulating water: The properties of water, such
as the uniqueness of its phase diagram, never stop surprising scientific
communities.[1] Given the vital importance of water in fields that vary from
material science to biology, there has been a recent surge in the development
and competition of different electronic structure methods for simulating
water.[2, 3, 4, 5, 6, 7, 8, 9] As ab initio quantum-chemical methods are too
expensive for large systems, Kohn-Sham density functional theory (KS-DFT) has
become a workhorse of electronic structure methods for running water
calculations.[10, 11, 12, 13, 14] But, despite an excellent accuracy to cost
ratio, historically KS-DFT has been unable to deliver sufficiently high
accuracy in water simulations to reproduce experimental data.[15, 16, 17, 18]

A recent breakthrough in this direction by Dasgupta et al. showed that the
strongly constrained and appropriately normed (SCAN) functional, when used in
tandem with ‘density-corrected DFT’ (DC-DFT), is a game changer for water
simulation, because it brings KS-DFT close to chemical accuracy.[3, 4] The
role of water in a chemical or biochemical reaction goes beyond providing an
environment to help a reaction in an aqueous solvation and is often explicitly
involved in the mechanism. For this reason, a complete understanding of the
reaction is possible only when the interaction between water and other
molecules is accurately described. Figure 1 shows how an integratively
designed DC-DFT procedure, HF-r2SCAN-DC4, describes not only the interactions
between water-water, water-organic molecules, and water-biochemical molecules
in various situations, but also the interactions of noncovalent complexes at
chemical accuracy or better.

The importance of the density: DC-DFT is a general framework that separates
errors of any DFT calculations into a contribution coming from the approximate
"D" (density) and the ’true’ error coming from the approximate "F"
(functional).[19, 20, 21, 22] In addition to being a rigorous exact theory,
DC-DFT gives practical guidance on when and how it can be used to reduce
errors in DFT simulation.[23, 24, 25, 26] Standard DFT calculations are
performed self-consistently (SC). The simplest form of practical DC-DFT is HF-
DFT, where density functionals are evaluated instead on Hartree-Fock (HF)
densities and orbitals.[27, 28, 29, 30, 31] While in most cases, SC-DFT gives
the best answer, in some errors in specific cases SC-DFT suffers from large
energetic errors due to the approximate density (density-driven errors).[19,
22] In such cases, HF-DFT typically yields significant improvements over SC-
DFT, and these include a number of chemical domains (barrier heights, some
torsional barriers, halogen bonds, anions, etc.).[26]

The importance of the functional: SCAN is a non-empirical meta-GGA functional
designed to satisfy 17 exact physical constraints, and to recover several
nonbonded ‘norms’.[32] Meta-GGA’s use the KS kinetic energy density as an
ingredient, but are not hybrid functionals like B3LYP[33, 34, 35, 36], which
include some fraction of exact exchange from a HF calculation.[35] In terms of
accuracy, SCAN is often on par with highly empirical more expensive density
functionals designed for molecules. At the same time, it enjoys great
successes for simulations of extended systems, making it one of the most-used
general-purpose functionals developed over the last 10 years.[37, 38, 39, 40,
41] Earlier works have shown that

Figure 1:  Performance of HF-r2SCAN-DC4 relative to HF-SCAN for various
chemical reactions: (a) the interaction energy of various configurations of
the stacked cytosine dimer, where HF-SCAN underbinds by 2-3 kcal/mol; (b)
energies of water hexamer relative to the lowest-lying prism isomer, with HF-
SCAN underestimating by almost 1 kcal/mol; (c) errors in binding energy of
WATER27 complexes as a function of density sensitivity (how much a DFT energy
changes when the density is changed), showing how large errors can be without
using the HF density. One cluster, H3O+(H2O)6OH- (at $x$ close to 4 kcal/mol)
is an outlier argued to exhibit a significant multiconfigurational
character[4]; (d) relative energies of water 20-mer isomers (not density
sensitive) from WATER27, where self-consistent SC-r2SCAN-D4 performs best, but
using the HF density introduces little error; (e) errors in interaction
energies in the water$\cdots$aspirin dimer structures from an MD simulation at
T=298.15 K; (f) mean-absolute-errors (MAEs) for intra- and inter-molecular
noncovalent interactions datasets from the GMTKN55 database. For more details,
see the main text and supporting information.

standard (SC) DFT calculations of water clusters suffer badly from density-
driven errors, which explains why HF-SCAN is much more accurate than its SC
counterpart for simulations of water.[3, 25] In addition to water clusters,
Dasgupta et al. used HF-SCAN in tandem with many-body potential energy
function related to the highly popular MB-pol[10, 11, 12] to run molecular
dynamics (MD) simulations of liquid water and obtained results in excellent
agreement with the experimental data. These were the first successful DFT-
based simulations able to correctly describe the condensation of water.

Nevertheless, the convergence of the SCAN functional can be painfully slow
with respect to the size of molecular grids, due to either the size of a
system or because it would require grids larger than those available in most
of the standard-quantum chemical codes.[31] Larger grids also lead to longer
computational times. To address these issues of SCAN, Perdew and co-workers
developed the regularized-restored SCAN functional (r2SCAN), which regularizes
SCAN but restores SCAN’s adherence to exact constraints.[42] But, as we show
below, a standalone version of HF-r2SCAN is much less accurate for water
simulations than HF-SCAN.

The vital importance of dispersion: Despite enormous success in modelling
water, HF-SCAN is not a panacea. In their water simulations, Dasgupta et al.
used HF-SCAN without dispersion correction, as they found that the standard
dispersion corrections, such as those of Grimme[43], worsen the original
results of HF-SCAN for water. But such dispersion corrections have long been
known to be necessary for noncovalent interactions (NCIs).[44, 45, 46, 47, 48,
49, 50, 51, 52, 53] So, despite delivering a high accuracy for pure water
simulations, HF-SCAN without a dispersion correction cannot describe
accurately long-range dispersion interactions. For this reason, the errors of
HF-SCAN are several times larger than those of DFT enhanced by a dispersion
correction for the standard noncovalent datasets.[40] The challenge is then to
construct an efficient density functional that correctly describes NCIs of
different nature, while recovering or even improving the accuracy of HF-SCAN
for water simulations.

HF-r2SCAN-DC4, an integratively designed DC-DFT procedure: In the present
paper, we resolve these issues by using the principles of DC-DFT to carefully
parameterize a dispersion correction for HF-r2SCAN. This yields HF-r2SCAN-DC4,
which produces the following key results: (i) HF-r2SCAN-DC4 improves upon HF-
SCAN for pure water simulations, by up to 0.7 kcal/mol for relative energies
of water hexamers, and up to 2.4 kcal/mol for those of water 20-mers; (ii)
HF-r2SCAN-DC4 is far more accurate than HF-SCAN for interactions of water with
other molecules and for NCIs in general, because of the inclusion of explicit
dispersion corrections; (iii) HF-r2SCAN-DC4 can be routinely and efficiently
used in calculations because, unlike HF-SCAN[31], HF-r2SCAN-DC4 has no grid
convergence issues. In our HF-r2SCAN-DC4, each of the three ingredients is
vitally important: The "HF" part reduces density-driven errors, while r2SCAN
fixes the grid issues of SCAN. But most importantly, the way in which we
parametrize the D4 corrections by using the DC-DFT principles is vital, as an
unwitting fitting of D4 ruins the accuracy for water simulations. If we drop
any of those elements of HF-r2SCAN-DC4, at least one of its three appealing
results will be lost.

Figure 2:  Water dimer interaction energies for (a) Smith stationary
points[54] and (b) MD simulated water dimers with the oxygen-oxygen distance.
For (a), MAEs of each functional are (following the order in the legend) 0.25,
0.11, 0.09, 0.17, and 0.08 kcal/mol. DLPNO-CCSD(T)-F12 has been used as a
reference. For (b), MAE of each functionals are 0.25, 0.08, 0.20, 0.30 and
0.08 kcal/mol. Figure S1 shows the corresponding density sensitivities and
Figure S2 shows the errors of approximations for the Smith dimers and
interaction energies for MD dimers.

To illustrate all these points, and how they work together, we created Figure
1. We show how HF-r2SCAN-DC4 is better than HF-SCAN for interactions of
nucleobases [panel (a)], water molecules with one another [panels (b), (c),
(d)], water with other molecules [panel (e)], and NCIs in general [panel (f)].

Stacking interactions in nucleobases are of vital importance in biology as
their energetics is essential to describe the formation and stability of DNA
and RNA.[55, 56] In Figure 1(a), we compare the accuracy of HF-SCAN and
HF-r2SCAN-DC4 for interaction energies of stacked cytosine dimers at different
configurations. As we can see from Figure 1(a), our HF-r2SCAN-DC4 essentially
greatly reduces the errors of HF-SCAN that systematically underbinds these
stacked complexes by about 2.5 kcal/mol. This demonstrates that despite its
success for modeling water, HF-SCAN misses most of dispersion and thus cannot
compete with our HF-r2SCAN-DC4 in modelling NCIs. This is especially the case
for NCIs dominated by dispersion interactions as those present in stacked
nucleobases. (See Figure S3 in the supporting information for the errors in
interaction energies.) We note that the mean absolute error (MAE) of
HF-r2SCAN-DC4 (0.4 kcal/mol) is very good relative to HF-SCAN, but not very
impressive relative to B3LYP-D3(BJ) (less than 0.2 kcal/mol).[56] But such
functionals include only a fraction of HF exchange, and so still suffer from
large density-driven errors in water, and so have larger errors for pure water
(as shown below).

Figure 3:  K-body interaction energy errors with (a) $K=2$, (b) $K=3$, and (c)
total, and (d) the interaction energy for 8 water hexamers. (For higher order
$K$-body interaction energies, see Figure S5 and S6.) Geometries and
CCSD(T)/CBS reference interaction energies are from Ref. [13]. The MAEs of
HF-r2SCAN-DC4 and HF-SCAN are 0.19 kcal/mol and 0.22 kcal/mol, respectively.

Water hexamers, "the smallest drops of water"[57, 58], are important, as they
represent the transition from two-dimensional to three-dimensional hydrogen-
bonding networks.[59, 60, 61] The energy differences between two adjacent
isomers of water hexamers are tiny, making even the ordering of isomers a very
challenging test for quantum-chemical methods.[59, 62] In Figure 1(b), we
compare the energies of water hexamer isomers relative to the energy of the
prism, as the lowest-lying isomer.[59, 63, 64] Despite being more accurate for
water hexamers than most DFT methods available on the market, HF-SCAN mistakes
the ordering of the isomers, as it predicts too low energies of the chair
isomer. Our HF-r2SCAN-DC4 is also here superior to HF-SCAN, as it not only
gives the right ordering of isomers, but essentially reproduces the reference
values for the relative energies of isomers. If D4 is fitted by not accounting
for the DC-DFT principles (see below), the accuracy of HF-r2SCAN-DC4 for the
water simulation is lost. This happened in Ref. [65] and will be discussed in
the method section.

We use the WATER27 dataset to illustrate the importance (and subtlety) of DC-
DFT for water simulations. WATER27 is a standard dataset for binding energies
of water clusters. Density sensitivity, $\tilde{S}$, is a measure for how
sensitive a given DFT simulation is to errors in densities (see Section S2 in
the supporting information for further details and specific definitions).[23]
Typically, the errors of SC-DFT calculations grow with $\tilde{S}$, indicating
the presence of large density-driven errors.[26, 31, 66] DC-DFT reduces these
large density-driven errors of SC-DFT and thus the errors of DC-DFT do not
grow with $\tilde{S}$. In Figure 1(c) we plot WATER27 errors as a function of
density sensitivity. As the errors of SC-r2SCAN-D4 grow with $\tilde{S}$, so
also does the energetic improvement of HF-r2SCAN-DC4 over SC-r2SCAN-D4.
Furthermore, sometimes dispersion corrections worsen SC-DFT for cases with
large density-driven errors.[25, 26] This is also the case here, as
SC-r2SCAN-D4 significantly deteriorates the accuracy of SC-r2SCAN (see Figure
S4). The errors of HF-SCAN are also substantially lower than those of
SC-r2SCAN-D4, and for most of the binding energies of the WATER27 clusters,
HF-SCAN is comparable to HF-r2SCAN-DC4. But, for the four clusters with the
largest sensitivities, HF-r2SCAN-DC4 outperforms HF-SCAN by $\sim$4 kcal/mol.

WATER27 is a part of the GMTKN55[40], a database that we use to train the D4
parameters in HF-r2SCAN-DC4 (see methods). But, according to the principles of
DC-DFT, we exclude those WATER27 clusters that are density-sensitive, as their
energetic errors are dominated by the errors in their densities.[26] Thus none
of the clusters that are to the right of the vertical dashed line placed at
$\tilde{S}$ =2 kcal/mol (see Method Section for the details on this reasoning)
are used in the fitting, which means HF-r2SCAN-DC4 makes genuinely accurate
predications for a vast majority of these water clusters. Not only does it
recover HF-SCAN for binding energies of the water clusters, but also provides
substantial improvements for the most challenging clusters.

An important question is whether or not one should always correct the density.
The general principles of DC-DFT say that one should only correct the density
in cases of substantial density-driven errors. In density insensitive cases,
the effect of correcting the density should be small, and may actually worsen
energetics. Figure 1(d) shows energies of water 20-mers relative to the energy
of the lowest of the four 20-mers. Here SC-r2SCAN-D4 beats its DC counterpart
every time. In contrast to large $\tilde{S}$ for binding energies of the four
20-mers (the last four datapoints in Figure 1(c)), the sensitivities
corresponding to their relative isomer energies are about twenty times smaller
(see Figure S7). Thus the higher accuracy of SC-r2SCAN-D4 over HF-r2SCAN-DC4
does not come as a surprise. But the crucial point is that, even in this low-
sensitivity scenario, the errors introduced by the HF density are far smaller
than those of HF-SCAN, and remain tiny on a per molecule basis.

A crucial figure of merit is how accurate energetics are for water molecules
in the vicinity of an organic molecule, especially if it is polar. In Figure
1(e), we show errors in the interaction energies between water and aspirin
from structures that we extracted from an MD simulation at T=298.15 K (see
Section S6 for further details on the MD simulation). The structures are
sorted by the distance between the oxygen atom in water and the specified
oxygen atom in the carboxyl group of aspirin. The errors of HF-r2SCAN-DC4 are
much smaller than those of HF-SCAN. They are also substantially smaller than
those of SC-r2SCAN-D4 (see Figure S8), demonstrating again the importance of
both the D4 and DC components in our method.

Getting NCI right across a broad range of molecules is important, even in the
absence of water. The GMTKN55 collection of 55 databases has become a standard
benchmark[40] and includes many databases for NCIs. In Figure 1(f), we compare
the MAEs of HF-SCAN and HF-r2SCAN-DC4 for the standard datasets with intra-
and intermolecular NCIs[40]. Despite its high accuracy for water clusters, HF-
SCAN does not capture long-ranged dispersion interactions. This is why it is
far less accurate than HF-r2SCAN-DC4 for noncovalent datasets. We can see that
HF-r2SCAN-DC4 is highly accurate here, and on average it beats SC-r2SCAN-D4
for both inter- and intramolecular NCIs (see Table S1 in the supporting
information comparing the metrics for overall performance).

## II. Results

### i. Interaction energies for water dimers

As discussed already, HF-SCAN performs incredibly well for interactions in
pure water. In this section, we look at select water dimers that are relevant
to water simulations, and show how HF-r2SCAN-DC4 reproduces (or even exceeds)
this accuracy. More importantly, we show how each aspect of its construction
(density correction, regularization of SCAN, and dispersion correction) are
vital to its accuracy for water. Later we will show that no other
approximation at this level of cost comes close to this performance for water.

Figure 2 shows the interaction energies for many water dimers (the difference
in the energies of a dimer and two monomers). (a) shows the interaction
energies at Smith stationary points, some of which resemble geometries from
dense ice structures.[18] (b) shows the errors of approximations in
interaction energies for water dimers as a function of the distance between
the two oxygen atoms. The underlying structures were extracted from an MD
simulation at T=298.15K (see Section S6 for further details on the
simulation). For the interaction energies of these water dimers, HF-SCAN
without a dispersion correction already provides a very high accuracy (with
MAEs of less than 0.1 kcal/mol). Our HF-r2SCAN-DC4 essentially recovers this
high accuracy of HF-SCAN. Similar patterns observed for binding energies of
water clusters are also seen here.

By studying the various plots, one can assess the importance of the relative
contributions to HF-r2SCAN-DC4. First, the purple points give HF-r2SCAN, to be
contrasted with HF-SCAN. We see that HF-r2SCAN significantly (on this scale)
underestimates the interaction energy. Even though r2SCAN was designed to
reproduce the results of SCAN, these differences are so small as to be
negligble for most purposes. However, they are clearly significant here,
showing HF-r2SCAN is noticeably less accurate for these dimers. The addition
of the D4 correction, however, makes their errors comparable.

On the other hand, we may also consider the importance of density correction.
We see that SC-r2SCAN-D4 considerably overestimates interaction energies. In
fact, SC-r2SCAN does rather well, as the errors due to poor density and
missing dispersion cancel.

We can also observe from Figure 2(b) that the improvement of HF-r2SCAN-DC4
over SC-r2SCAN-D4 decreases with the distance between the two oxygen atoms in
water dimers. This can be understood in terms of underlying density
sensitivity which also decreases with the O-O distance (see Figure S1).

### ii. Many-body interactions in larger water clusters

In Figure 3 we compare errors of HF-r2SCAN-DC4 and HF-SCAN for the interaction
energies of the eight standard water hexamers.[59, 60] In addition to total
interaction energies, we also use the many-body expansion (MBE) to show the
$K$-body contributions to these energies (with $K$ in between 2 and 6). This
is a standard methodology for understanding the origins of errors in water
models.[3, 13, 67] The energetic importance of the $K$-body contributions
decreases rapidly with $K$ (Figure S6), making the 2-body contributions by far
the most important, and these are where significant differences emerge when
the density is corrected. But in order to reach chemical accuracy, a proper
description of the higher-order contributions also matters. The 2-body plot
shows that HF-SCAN has a rather systemative overestimate of about 0.5
kcal/mol, whereas HF-r2SCAN-DC4 is substantially less for about half the
clusters. The 3-body plot shows them being almost identical. But in the total
error, we see that HF-r2SCAN-DC4 is far more systematic, as HF-SCAN makes
errors of opposite sign, while HF-r2SCAN-DC4 is always an overestimate of
about 0.2 kcal/mol.

This consistency is important on the plot (d), showing the interaction energy
of the 8 hexamers. Because HF-r2SCAN-DC4 is so consistent, it gets the
ordering in interaction energies of all clusters correct, whereas HF-SCAN
incorrectly predicts that the interaction energy in the bag is higher than
that of the chair. The MAE of HF-r2SCAN-DC4 is 0.19 kcal/mol, lower than 0.22
kcal/mol for HF-SCAN. On average, HF-r2SCAN-DC4 also improves individual
$K$-body contributions to the interaction energies, except for $K=4$, where
both are marginally small (Figure S5). This MBE test shows us that the
improvement of HF-r2SCAN-DC4 over HF-SCAN for the water hexamer interaction
energies (seen also for the relative isomer energies (Figure 1(b)) is
systematic and does not result from the error cancellations between different
$K$-body contributions (for the detailed information of water hexamer
isomerization energy in Figure 1, see Figure S9).

### iii. Water$\cdots$cytosine interaction energies

Figure 4:  Errors in interaction energies of water$\cdots$cytosine complexes
sorted by the distance between the oxygen atom in cytosine and the oxygen atom
in water. Reference interaction energies have been computed at the DLPNO-
CCSD(T)-F12/AVQZ level of theory.  Figure 5:  (a) The mean absolute error
(MAE) for the water-based reactions appear in this work (hexamer isomer
energies, water 20-mers isomer energies, WATER27 binding energies, water-small
organic molecule interaction energies, and water dimer interaction energies)
versus the weighted-mean-absolute-deviation-2 (WTMAD-2) for the GMTKN55
database for selected functionals. For a further description of the reactions
used in the $y$-axis, see Section S8 in the supporting information. HF-SCAN-D4
functional used here is from Ref. [68]. (b) The hexagon plot with MAEs for
selected water-based datasets and WTMAD-2 values for the whole GMTKN55
databases (for WTMAD-2 values for other GMTKN55 database, see Figure S10).
Abbreviations of isomerization (I), binding (B), and interaction (T) energy
are noted in the vertex caption. MAEs of HF-r2SCAN-DC4 for individual GMTKN55
datasets are shown in Table S1. In Figure S11, we give further details about
the interaction energies used in the "water-small organic molecule" dataset.

IIn Figure 4, we study the performance of different variations for
microhydration of cytosine, by specifically focusing on the interaction
energies in water$\cdots$cytosine complexes. We generate these complexes as
described in Section S7, and in all of them, water interacts with cytosine
through the hydrogen bond formed between the hydrogen atom in water and the
oxygen atom in cytosine. For each complex, the errors of HF-r2SCAN-DC4 are
small, and with the MAE of 0.09 kcal/mol, it is the best performer in Figure
4.

The errors of HF-SCAN are much smaller here than for cytosine dimers (Figure
1(a)), in which the role of dispersion is more important. Nevertheless,
HF-r2SCAN-DC4 provides here a significant improvement over HF-SCAN. It is also
interesting to observe what happens after we add the dispersion correction to
HF-r2SCAN and its SC counterpart. In the case of HF-r2SCAN, the errors in the
interaction energies are greatly reduced (roughly by a factor of 6 on
average). In stark contrast, adding D4 to SC-r2SCAN significantly deteriorates
its accuracy, as SC-r2SCAN already overbinds water$\cdots$cytosine complexes
and D4 makes the overbinding stronger.

### iv. Wide applicability of HF-r2SCAN-DC4

A functional that works extremely well for pure water but nothing else is not
widely applicable. Recently, GMTKN55 of 55 databases has become a popular
benchmark for testing the accuracy of density functionals for main-group
chemistry. Figure 5 has been designed to illustrate performance of functionals
for both pure water and on the GMTKN55 database simultaneously. The water
metric ($y$-axis on the left) combines most of the reactions with water used
in this paper, and is carefully defined in the supporting information Section
S8.

Figure 5(a) shows errors on GMTKN55 on the $x$-axis and errors on the water
metric on the $y$-axis, each in kcal/mol. The $x$-axis ranges from about 3-10
kcal/mol, spanning the performance of modern approximations for main group
chemistry, such as atomization energies. The $y$-axis range is much smaller,
running less than 4.0 kcal/mol, reflecting the much smaller magnitude of NCIs
in water, and how high accuracy needs to be in order to have an accurate model
for water. Here, HF-SCAN sets a high standard, with a water error near 1.0
kcal/mol (the chemical accuracy claimed in Ref. [3]), while most standard-use
functionals cannot compete. On the other hand, SCAN is designed mainly to
improve materials calculations without the cost of a hybrid functional, and
HF-SCAN has a high error on GMTKN55 (about 9 kcal/mol). Popular functionals
have much smaller GMTKN55 errors, but perform worse on water. We also show the
many combinations of HF-r2SCAN-DC4 that do not include all the right
ingredients, showing they all perform less well on water than HF-SCAN. We
finally include $\omega$B97M-V functional [69], which might be considered the
DFT gold-standard here, with the smallest errors for both water and main-group
chemistry. But this range-separated functional with nonlocal correlation
functional is far more expensive to compute than most functionals ,including
its own D4 variant,[70] and is less practical for DFT-MD simulations than
e.g., SCAN. We have included it here only to show what is possible in
principle with DFT.

Figure 6:  Comparison between HF-r2SCAN-DC4 and SM21 for (a) water hexamer
isomerization energy, (b) water hexamer interaction energy, (c) water 20-mer
isomerization energy, and (d) hexgonal plot same as Figure 5(b). SM21 (green)
is HF-r2SCAN but with different D4 parameters obtained in Ref. [65].

But the performance of HF-r2SCAN-DC4 is remarkable. Its errors on both water
and the GMTKN55 dataset are almost half of those of HF-SCAN. No other
functional in our collection comes close for water. Clearly, all the
chemically-inclined approximations which are comparable for main-group
chemistry do much worse.

In Figure 5(b), we show the hexagon plots comparing the MAEs of several
density functionals, where the position of five vertices denote the MAEs for
individual water-based datasets, while the sixth vertex denotes the overall
performance of the functionals for the whole GMTKN55 databases, as measured by
the weighted-mean-absolute-deviation-2 (WTMAD-2). It is the MAE for all the
reactions from these five water-based datasets that we use as the quantity on
the $y$-axis in Figure 5(a). The size of the hexagon of HF-r2SCAN-DC4 is the
closest to that of more costly $\omega$B97M-V. We can also see that the
performance of HF-r2SCAN-DC4 is far superior to that of HF-SCAN. M062X-D3(0),
a meta-hybrid that is very accurate for small organic molecules,[40] and
yields WTMAD-2 which is slightly lower than that of HF-r2SCAN-DC4. But, for
water simulations, M062X-D3(0) is nowhere close to HF-r2SCAN-DC4, as can be
seen from the position of the remaining five vertices.

## III. Methods

The basic principles of DC-DFT are covered elsewhere in the literature[19,
22], and reviewed in the supporting information. In most KS-DFT calculations,
the error in the density has a negligible effect on the energy errors. But
sometimes the error in a SC density leads to a noticeable contribution, which
can be reduced if a more accurate density is used instead. For many semilocal
exchange-correlation approximations in molecular calculations, when a
calculation is density sensitive, often the HF density then yields
significantly smaller energy errors. These principles have led to improved
energetics in reaction barrier heights, electron affinities, and also for the
ground state geometries of noncovalent interaction systems, etc.[20, 71, 72,
73]

Application of the principles of DC-DFT is subtle in the case of r2SCAN-D4,
because of the need to separate out the error due to density correction from
the fitting of the D4 corrections. For example, for halogen bonds, the
density-driven errors are far larger than dispersion corrections, so all
fitting must be done on density-corrected energetics. Moreover, when empirical
functionals contain parameters, such parameters should be fit only on density-
insensitive calculations, so that the parameters optimize the ‘true’
functional error.

With these principles in mind, we find the parameters for HF-r2SCAN-DC4 using
the density-insensitive calculations in the GMTKN55 dataset as a training set
while using water$\cdots$water pair interaction energy as a validation set. We
find their optimum values by minimizing MAE values over all such cases. This
is detailed in the supporting information section S4. This is why we use the
acronym DC4 instead of D4, meaning that we use the principles of DC-DFT to
find the underlying D4 parameters.

In Refs. [31] and [66], we proposed DC(HF)-DFT, a DC-DFT procedure that
discriminately uses HF densities based on the density sensitivity criterion.
The main idea of DC(HF)-DFT is to use HF-DFT for density-sensitive (DS)
reactions and SC-DFT for density-insensitive (DI) reactions (possible spin-
contaminations of the HF results are also taken into account as detailed in
Section S3). While we consider DC(HF)-DFT a state-of-the-art DC-DFT procedure,
for our HF-r2SCAN-DC4 we use HF-DFT, meaning that the functional is always
evaluated on the HF density regardless of the sensitivity criterion. To use
DC(HF)-DFT, we need to compute density sensitivity for each reaction of
interest and possibly make adjustments to its cutoff value which is used to
declare whether a given reaction is DS or DI.[31, 66] This would also require
having two sets of D4 parameters, one for DS and the other for DI reactions.
All these efforts would undermine the ease of use of r2SCAN, which is a
general-purpose functional. For this reason and encouraged by the very good
performance of HF-DFT with SCAN-like functionals[32, 42], we employ HF-DFT[74]
as a DC-DFT procedure for HF-r2SCAN-DC4. While our HF-r2SCAN-DC4 can be
routinely used by applying it to HF orbitals without ever needing to calculate
density sensitivity of a given reaction, the use of DC-DFT principles and
density sensitivity is vital for our training of HF-r2SCAN-DC4 as explained
above.

To illustrate what can happen when these principles are not applied, we show
results from Ref. [65]. This is a version of HF-r2SCAN-D4, but where all
reactions in GMTKN55 were used, and the WTMAD-2[40] was used as the cost
function instead. Figure 6 illustrates the results for the larger water
clusters. In every case, they are noticeably worse than ours. Moreover, (d)
shows that, apart from matching on WTMAD-2 measure, HF-r2SCAN-DC4 yields more
accurate results in every other case.

## IV. Conclusions

The work of Ref. [3] was a breakthrough in models for water, showing that, by
using the principles of DC-DFT, a moderate-cost density functional
approximation approached chemical accuracy for many relevant properties of
small water clusters. However that functional is lacking in dispersion
corrections, yielding large errors for energetics between organic and
biological molecules. It also inherits some of the numerical issues of the
original SCAN functional, which have been eliminated by using r2SCAN instead
in most other applications. However, the small differences between these two
wreak havoc on the much smaller scale of subtle energy differences of water
clusters.

The present work shows that, by a very careful application of the principles
of DC-DFT, all these difficulties can be overcome, and even greater accuracy
achieved for pure water, while still including dispersion for other molecules
where it can be vital. Finding the correct parameters depends crucially on
training on only density-insensitive chemical reactions, as inclusion of
density-sensitive reactions yields suboptimal values for the parameters.

Even if HF-r2SCAN-DC4 could be run at close to meta-GGA cost, KS-DFT MD
simulations are typically far more costly than MD with machine learning (ML)
interatomic potentials. But accurate force-field generation requires highly
accurate reference energetics data as a training set, and CCSD(T) or Quantum
Monte Carlo (QMC) methods are frequently used as reference methods these
days.[12] Due to the large computational cost for such ab initio calculation,
a more practical yet accurate method is in demand, and HF-r2SCAN-DC4 can
replace them for calculating moderately large biomolecular systems. We suggest
HF-r2SCAN-DC4 be tested and applied in solution wherever practical.

Acknowledgement

ES, SS, YK, and HY are grateful for support from the National Research
Foundation of Korea (NRF-2020R1A2C2007468 and NRF-2020R1A4A1017737) and
Samsung Electronics (IO211126-09176-01). KB acknowledges funding from NSF
(CHE-2154371). SV acknowledges funding from the Marie Skłodowska-Curie grant
101033630 (EU’s Horizon 2020 programme). We thank John Perdew and Francesco
Paesani and his group for many useful discussions.

## References

  * [1] Linfeng Zhang, Han Wang, Roberto Car, and E Weinan.  Phase diagram of a deep potential water model.  Physical Review Letters, 126(23):236001, 2021.
  * [2] Kamal Sharkas, Kamal Wagle, Biswajit Santra, Sharmin Akter, Rajendra R Zope, Tunna Baruah, Koblar A Jackson, John P Perdew, and Juan E Peralta.  Self-interaction error overbinds water clusters but cancels in structural energy differences.  Proceedings of the National Academy of Sciences, 117(21):11283–11288, 2020.
  * [3] Saswata Dasgupta, Eleftherios Lambros, John P Perdew, and Francesco Paesani.  Elevating density functional theory to chemical accuracy for water simulations through a density-corrected many-body formalism.  Nature Communications, 12(1):1–12, 2021.
  * [4] Saswata Dasgupta, Chandra Shahi, Pradeep Bhetwal, John P. Perdew, and Francesco Paesani.  How good is the density-corrected scan functional for neutral and ionic aqueous systems, and what is so right about the Hartree-Fock density?  Journal of Chemical Theory and Computation, 18(8):4745–4761, 2022\.
  * [5] Chunyi Zhang, Fujie Tang, Mohan Chen, Jianhang Xu, Linfeng Zhang, Diana Y Qiu, John P Perdew, Michael L Klein, and Xifan Wu.  Modeling liquid water by climbing up jacob’s ladder in density functional theory facilitated by using deep neural network potentials.  Journal of Physical Chemistry B, 125(41):11444–11456, 2021.
  * [6] Alberto Torres, Luana S Pedroza, Marivi Fernandez-Serra, and Alexandre R Rocha.  Using neural network force fields to ascertain the quality of ab initio simulations of liquid water.  Journal of Physical Chemistry B, 125(38):10772–10778, 2021.
  * [7] Eleftherios Lambros, Jie Hu, and Francesco Paesani.  Assessing the accuracy of the scan functional for water through a many-body analysis of the adiabatic connection formula.  Journal of Chemical Theory and Computation, 17(6):3739–3749, 2021\.
  * [8] Etienne Palos, Eleftherios Lambros, Steven Swee, Jie Hu, Saswata Dasgupta, and Francesco Paesani.  Assessing the interplay between functional-driven and density-driven errors in dft models of water.  Journal of Chemical Theory and Computation, 18(6):3410–3426, 2022\.
  * [9] Etienne Palos, Eleftherios Lambros, Saswata Dasgupta, and Francesco Paesani.  Density functional theory of water with the machine-learned dm21 functional.  Journal of Chemical Physics, 156(16):161103, 2022.
  * [10] Volodymyr Babin, Claude Leforestier, and Francesco Paesani.  Development of a “first principles” water potential with flexible monomers: Dimer potential energy surface, vrt spectrum, and second virial coefficient.  Journal of Chemical Theory and Computation, 9(12):5395–5403, 2013\.
  * [11] Volodymyr Babin, Gregory R Medders, and Francesco Paesani.  Development of a “first principles” water potential with flexible monomers. ii: Trimer potential energy surface, third virial coefficient, and small clusters.  Journal of Chemical Theory and Computation, 10(4):1599–1607, 2014\.
  * [12] Gregory R Medders, Volodymyr Babin, and Francesco Paesani.  Development of a “first principles” water potential with flexible monomers. iii. liquid phase properties.  Journal of Chemical Theory and Computation, 10(8):2906–2910, 2014\.
  * [13] Sandeep K Reddy, Shelby C Straight, Pushp Bajaj, C. Huy Pham, Marc Riera, Daniel R Moberg, Miguel A Morales, Chris Knight, Andreas W Götz, and Francesco Paesani.  On the accuracy of the mb-pol many-body potential for water: Interaction energies, vibrational frequencies, and classical thermodynamic and dynamical properties from clusters to liquid water and ice.  Journal of Chemical Physics, 145(19):194504, 2016.
  * [14] Mohan Chen, Hsin-Yu Ko, Richard C Remsing, Marcos F Calegari Andrade, Biswajit Santra, Zhaoru Sun, Annabella Selloni, Roberto Car, Michael L Klein, John P Perdew, et al.  Ab initio theory and modeling of water.  Proceedings of the National Academy of Sciences, 114(41):10846–10851, 2017.
  * [15] I-Feng W Kuo, Christopher J Mundy, Matthew J McGrath, J Ilja Siepmann, Joost VandeVondele, Michiel Sprik, Jürg Hutter, Bin Chen, Michael L Klein, Fawzi Mohamed, et al.  Liquid water from first principles: Investigation of different sampling approaches.  Journal of Physical Chemistry B, 108(34):12990–12998, 2004.
  * [16] Thomas D K$\ddot{u}$hne, Matthias Krack, and Michele Parrinello.  Static and dynamical properties of liquid water from first principles by a novel car- parrinello-like approach.  Journal of Chemical Theory and Computation, 5(2):235–241, 2009\.
  * [17] Biswajit Santra, Jiří Klimeš, Dario Alfe, Alexandre Tkatchenko, Ben Slater, Angelos Michaelides, Roberto Car, and Matthias Scheffler.  Hydrogen bonds and van der waals forces in ice at ambient and high pressures.  Physical Review Letters, 107(18):185701, 2011.
  * [18] Michael J Gillan, Dario Alfe, and Angelos Michaelides.  Perspective: How good is dft for water?  Journal of Chemical Physics, 144(13):130901, 2016.
  * [19] Min-Cheol Kim, Eunji Sim, and Kieron Burke.  Understanding and reducing errors in density functional calculations.  Physical Review Letters, 111(7):073003, 2013.
  * [20] Min-Cheol Kim, Eunji Sim, and Kieron Burke.  Ions in solution: Density corrected density functional theory (dc-dft).  Journal of Chemical Physics, 140(18):18A528, 2014.
  * [21] Adam Wasserman, Jonathan Nafziger, Kaili Jiang, Min-Cheol Kim, Eunji Sim, and Kieron Burke.  The importance of being self-consistent.  Annual Review of Physical Chemistry, 68(1):555–581, 2017.
  * [22] Stefan Vuckovic, Suhwan Song, John Kozlowski, Eunji Sim, and Kieron Burke.  Density functional analysis: The theory of density-corrected dft.  Journal of Chemical Theory and Computation, 15(12):6636–6646, 2019\.
  * [23] Eunji Sim, Suhwan Song, and Kieron Burke.  Quantifying density errors in dft.  Journal of Physical Chemistry Letters, 9(22):6385–6392, 2018.
  * [24] Seungsoo Nam, Suhwan Song, Eunji Sim, and Kieron Burke.  Measuring density-driven errors using kohn–sham inversion.  Journal of Chemical Theory and Computation, 16(8):5014–5023, July 2020.
  * [25] Subrata Jana, Abhilash Patra, Szymon Śmiga, Lucian A Constantin, and Prasanjit Samal.  Insights from the density functional performance of water and water–solid interactions: Scan in relation to other meta-ggas.  Journal of Chemical Physics, 153(21):214116, 2020.
  * [26] Suhwan Song, Stefan Vuckovic, Eunji Sim, and Kieron Burke.  Density sensitivity of empirical functionals.  Journal of Physical Chemistry Letters, 12(2):800–807, 2021.
  * [27] Min-Cheol Kim, Eunji Sim, and Kieron Burke.  Communication: Avoiding unbound anions in density functional calculations.  Journal of Chemical Physics, 134(17):171103, 2011.
  * [28] Min-Cheol Kim, Hansol Park, Suyeon Son, Eunji Sim, and Kieron Burke.  Improved dft potential energy surfaces via improved densities.  Journal of Physical Chemistry Letters, 6(19):3802–3807, 2015.
  * [29] Suhwan Song, Min-Cheol Kim, Eunji Sim, Anouar Benali, Olle Heinonen, and Kieron Burke.  Benchmarks and reliable dft results for spin gaps of small ligand fe (ii) complexes.  Journal of Chemical Theory and Computation, 14(5):2304–2311, 2018\.
  * [30] Yeil Kim, Suhwan Song, Eunji Sim, and Kieron Burke.  Halogen and chalcogen binding dominated by density-driven errors.  Journal of Physical Chemistry Letters, 10(2):295–301, 2018.
  * [31] Suhwan Song, Stefan Vuckovic, Eunji Sim, and Kieron Burke.  Density-corrected dft explained: Questions and answers.  Journal of Chemical Theory and Computation, 18(2):817–827, 2022\.
  * [32] Jianwei Sun, Adrienn Ruzsinszky, and John P Perdew.  Strongly constrained and appropriately normed semilocal density functional.  Physical Review Letters, 115(3):036402, 2015.
  * [33] Axel D Becke.  Density-functional exchange-energy approximation with correct asymptotic behavior.  Physical Review A, 38(6):3098–3100, 1988.
  * [34] Chengteh Lee, Weitao Yang, and Robert G Parr.  Development of the colle-salvetti correlation-energy formula into a functional of the electron density.  Physical Review B, 37(2):785–789, 1988.
  * [35] Axel D Becke.  Density-functional thermochemistry. iii. the role of exact exchange.  Journal of Chemical Physics, 98(7):5648–5652, 1993.
  * [36] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch.  Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields.  Journal of Physical Chemistry, 98(45):11623–11627, 1994.
  * [37] Jianwei Sun, Richard C Remsing, Yubo Zhang, Zhaoru Sun, Adrienn Ruzsinszky, Haowei Peng, Zenghui Yang, Arpita Paul, Umesh Waghmare, Xifan Wu, et al.  Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional.  Nature Chemistry, 8(9):831–836, 2016.
  * [38] Fabien Tran, Julia Stelzl, and Peter Blaha.  Rungs 1 to 4 of dft jacob’s ladder: Extensive test on the lattice constant, bulk modulus, and cohesive energy of solids.  Journal of Chemical Physics, 144(20):204120, 2016.
  * [39] Narbe Mardirossian and Martin Head-Gordon.  Thirty years of density functional theory in computational chemistry: an overview and extensive assessment of 200 density functionals.  Molecular Physics, 115(19):2315–2372, 2017.
  * [40] Lars Goerigk, Andreas Hansen, Christoph Bauer, Stephan Ehrlich, Asim Najibi, and Stefan Grimme.  A look at the density functional theory zoo with the advanced gmtkn55 database for general main group thermochemistry, kinetics and noncovalent interactions.  Physical Chemistry Chemical Physics, 19(48):32184–32215, 2017.
  * [41] Yubo Zhang, Jianwei Sun, John P Perdew, and Xifan Wu.  Comparative first-principles studies of prototypical ferroelectric materials by lda, gga, and scan meta-gga.  Physical Review B, 96(3):035143, 2017.
  * [42] James W Furness, Aaron D Kaplan, Jinliang Ning, John P Perdew, and Jianwei Sun.  Accurate and numerically efficient r2scan meta-generalized gradient approximation.  Journal of Physical Chemistry Letters, 11(19):8208–8215, 2020.
  * [43] Stefan Grimme, Jens Antony, Stephan Ehrlich, and Helge Krieg.  A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu.  Journal of Chemical Physics, 132(15):154104, 2010.
  * [44] Xueyuan Wu, MC Vargas, S Nayak, V Lotrich, and G Scoles.  Towards extending the applicability of density functional theory to weakly bound systems.  Journal of Chemical Physics, 115(19):8748–8757, 2001.
  * [45] Evert Jan Meijer and Michiel Sprik.  A density-functional study of the intermolecular interactions of benzene.  Journal of Chemical Physics, 105(19):8684–8689, 1996.
  * [46] Stefan Grimme.  Accurate description of van der waals complexes by density functional theory including empirical corrections.  Journal of Computational Chemistry, 25(12):1463–1473, 2004.
  * [47] Axel D Becke and Erin R Johnson.  Exchange-hole dipole moment and the dispersion interaction.  Journal of Chemical Physics, 122(15):154104, 2005.
  * [48] Stefan Grimme.  Semiempirical gga-type density functional constructed with a long-range dispersion correction.  Journal of Computational Chemistry, 27(15):1787–1799, 2006.
  * [49] Axel D Becke and Erin R Johnson.  Exchange-hole dipole moment and the dispersion interaction revisited.  Journal of Chemical Physics, 127(15):154108, 2007.
  * [50] Alexandre Tkatchenko and Matthias Scheffler.  Accurate molecular van der waals interactions from ground-state electron density and free-atom reference data.  Physical Review Letters, 102(7):073005, 2009.
  * [51] Oleg A Vydrov and Troy Van Voorhis.  Nonlocal van der waals density functional: The simpler the better.  Journal of Chemical Physics, 133(24):244103, 2010.
  * [52] Jan Rezac and Pavel Hobza.  Benchmark calculations of interaction energies in noncovalent complexes and their applications.  Chemical Reviews, 116(9):5038–5071, 2016.
  * [53] Alastair JA Price, Kyle R Bryenton, and Erin R Johnson.  Requirements for an accurate dispersion-corrected density functional.  Journal of Chemical Physics, 154(23):230902, 2021.
  * [54] Brian J Smith, David J Swanton, John A Pople, Henry F Schaefer III, and Leo Radom.  Transition structures for the interchange of hydrogen atoms within the water dimer.  Journal of Chemical Physics, 92(2):1240–1247, 1990.
  * [55] Ilyas Yildirim and Douglas H Turner.  Rna challenges for computational chemists.  Biochemistry, 44(40):13225–13234, 2005.
  * [56] Holger Kruse and Jiri Sponer.  Revisiting the potential energy surface of the stacked cytosine dimer: Fno-ccsd (t) interaction energies, sapt decompositions, and benchmarking.  Journal of Physical Chemistry A, 123(42):9209–9222, 2019.
  * [57] K Nauta and RE Miller.  Formation of cyclic water hexamer in liquid helium: The smallest piece of ice.  Science, 287(5451):293–295, 2000.
  * [58] Yimin Wang, Volodymyr Babin, Joel M Bowman, and Francesco Paesani.  The water hexamer: cage, prism, or both. full dimensional quantum simulations say both.  Journal of the American Chemical Society, 134(27):11116–11119, 2012\.
  * [59] Desiree M Bates and Gregory S Tschumper.  Ccsd (t) complete basis set limit relative energies for low-lying water hexamer structures.  Journal of Physical Chemistry A, 113(15):3555–3559, 2009.
  * [60] Yiming Chen and Hui Li.  Intermolecular interaction in water hexamer.  Journal of Physical Chemistry A, 114(43):11719–11724, 2010.
  * [61] A Otero-De-La-Roza and Erin R Johnson.  Non-covalent interactions and thermochemistry using xdm-corrected hybrid and range-separated hybrid density functionals.  Journal of Chemical Physics, 138(20):204109, 2013.
  * [62] Biswajit Santra, Angelos Michaelides, Martin Fuchs, Alexandre Tkatchenko, Claudia Filippi, and Matthias Scheffler.  On the accuracy of density-functional theory exchange-correlation functionals for h bonds in small water clusters. ii. the water hexamer and van der waals interactions.  Journal of Chemical Physics, 129(19):194111, 2008.
  * [63] Ryan M Olson, Jonathan L Bentz, Ricky A Kendall, Michael W Schmidt, and Mark S Gordon.  A novel approach to parallel coupled cluster calculations: Combining distributed and shared memory techniques for modern cluster based systems.  Journal of Chemical Theory and Computation, 3(4):1312–1328, 2007\.
  * [64] MJ Gillan, FR Manby, MD Towler, and D Alfè.  Assessing the accuracy of quantum monte carlo and density functional theory for energetics of small water clusters.  Journal of Chemical Physics, 136(24):244105, 2012.
  * [65] Golokesh Santra and Jan ML Martin.  Pure and hybrid scan, rscan, and r2scan: Which one is preferred in ks-and hf-dft calculations, and how does d4 dispersion correction affect this ranking?  Molecules, 27(1):141, 2021.
  * [66] Eunji Sim, Suhwan Song, Stefan Vuckovic, and Kieron Burke.  Improving results by improving densities: Density-corrected density functional theory.  Journal of the American Chemical Society, 144(15):6625–6639, 2022\.
  * [67] Urszula Góra, Rafał Podeszwa, Wojciech Cencek, and Krzysztof Szalewicz.  Interaction energies of large clusters from many-body expansion.  Journal of Chemical Physics, 135(22):224102, 2011.
  * [68] Golokesh Santra and Jan ML Martin.  What types of chemical problems benefit from density-corrected dft? a probe using an extensive and chemically diverse test suite.  Journal of Chemical Theory and Computation, 17(3):1368–1379, 2021\.
  * [69] Narbe Mardirossian and Martin Head-Gordon.  $\omega$b97m-v: A combinatorially optimized, range-separated hybrid, meta-gga density functional with vv10 nonlocal correlation.  Journal of Chemical Physics, 144(21):214110, 2016.
  * [70] Narbe Mardirossian and Martin Head-Gordon.  Survival of the most transferable at the top of jacob’s ladder: Defining and testing the $\omega$ b97m (2) double hybrid density functional.  Journal of Chemical Physics, 148(24):241736, 2018.
  * [71] Benjamin G Janesko and Gustavo E Scuseria.  Hartree–fock orbitals significantly improve the reaction barrier heights predicted by semilocal density functionals.  Journal of Chemical Physics, 128(24):244112, 2008.
  * [72] Donghyung Lee, Filipp Furche, and Kieron Burke.  Accuracy of electron affinities of atoms in approximate density functional theory.  Journal of Physical Chemistry Letters, 1(14):2124–2129, 2010.
  * [73] Donghyung Lee and Kieron Burke.  Finding electron affinities with approximate density functionals.  Molecular Physics, 108(19-20):2687–2701, 2010.
  * [74] Yonsei University Theoretical and Computational Chemistry Laboratory.  Density Corrected-Density Functional Theory.  http://tccl.yonsei.ac.kr/mediawiki/index.php/DC-DFT, 2022.  (accessed 2022-07-07).
  * [75] Walter Kohn and Lu Jeu Sham.  Self-consistent equations including exchange and correlation effects.  Physical Review, 140(4A):A1133, 1965.
  * [76] John P Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple.  Physical Review Letters, 77(18):3865, 1996.
  * [77] Enrico Clementi and Subhas J Chakravorty.  A comparative study of density functional models to estimate molecular atomization energies.  Journal of Chemical Physics, 93(4):2591–2602, 1990.
  * [78] Peter MW Gill, Benny G Johnson, John A Pople, and Michael J Frisch.  An investigation of the performance of a hybrid of hartree-fock and density functional theory.  International Journal of Quantum Chemistry, 44(S26):319–331, 1992\.
  * [79] Peter MW Gill, Benny G Johnson, John A Pople, and Michael J Frisch.  The performance of the Becke-Lee-Yang-Parr (B-LYP) density functional theory with various basis sets.  Chemical Physics Letters, 197(4-5):499–505, 1992.
  * [80] Gustavo E Scuseria.  Comparison of coupled-cluster results with a hybrid of hartree–fock and density functional theory.  Journal of chemical physics, 97(10):7528–7530, 1992.
  * [81] Nevin Oliphant and Rodney J Bartlett.  A systematic comparison of molecular properties obtained using hartree–fock, a hybrid hartree–fock density-functional-theory, and coupled-cluster methods.  Journal of Chemical Physics, 100(9):6550–6561, 1994.
  * [82] Adam Rettig, Diptarka Hait, Luke W Bertels, and Martin Head-Gordon.  Third-order møller–plesset theory made more useful? the role of density functional theory orbitals.  Journal of Chemical Theory and Computation, 16(12):7473–7489, 2020\.
  * [83] Eike Caldeweyher, Sebastian Ehlert, Andreas Hansen, Hagen Neugebauer, Sebastian Spicher, Christoph Bannwarth, and Stefan Grimme.  A generally applicable atomic-charge dependent london dispersion correction.  Journal of Chemical Physics, 150(15):154122, 2019.
  * [84] Christoph Bannwarth, Sebastian Ehlert, and Stefan Grimme.  Gfn2-xtb-an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions.  Journal of Chemical Theory and Computation, 15(3):1652–1671, 2019\.
  * [85] Sebastian Spicher and Stefan Grimme.  Robust atomistic modeling of materials, organometallic, and biochemical systems.  Angewandte Chemie International Edition, 59(36):15665–15673, 2020\.
  * [86] Frank Neese, Frank Wennmohs, Ute Becker, and Christoph Riplinger.  The orca quantum chemistry program package.  Journal of Chemical Physics, 152(22):224108, 2020.
  * [87] Sebastian Ehlert, Uwe Huniar, Jinliang Ning, James W Furness, Jianwei Sun, Aaron D Kaplan, John P Perdew, and Jan Gerit Brandenburg.  r2scan-d4: Dispersion corrected meta-generalized gradient approximation for general chemical applications.  Journal of Chemical Physics, 154(6):061101, 2021.
  * [88] Jan Rezác.  Non-covalent interactions atlas benchmark data sets: Hydrogen bonding.  Journal of Chemical Theory and Computation, 16(4):2355–2368, 2020\.

Supporting Information

## S1. DFT and DC-DFT

In KS-DFT [75], the ground state energy is obtained by minimizing the
following functional over densities,

$E_{v}[n]=T_{{}_{\scriptscriptstyle\rm S}}[n]+U_{\scriptscriptstyle\rm
H}[n]+E_{\scriptscriptstyle\rm XC}[n]+\int d^{3}r\,n({\bf r})v({\bf r}),$ (1)

where the minimizing $n_{v}({\bf r})$ is the ground-state density,
$T_{{}_{\scriptscriptstyle\rm S}}[n]$ is the KS noninteracting kinetic energy
functional, $U_{\scriptscriptstyle\rm H}[n]$ the Hartree energy, and
$E_{{}_{\scriptscriptstyle\rm XC}}[n]$ is the exchange-correlation (XC)
functional. In practical KS calculations, the only term that of Eq. 1 that is
approximated is $E_{{}_{\scriptscriptstyle\rm XC}}[n]$. Whatever approximation
we use for XC, it yields an error in both energies and the densities, as the
minimizing density in Eq. 1 is different from the exact when an approximate XC
functional is used in place of its exact counterpart. We can write the error
of any KS calculation as: $\Delta E=\tilde{E}[\tilde{n}]-E[n]$ as, where
$E[n]$ is the exact functional (Eq. 1) and $n$ is the exact density (we drop
the $v$ subscript for brevity), while tildes denote approximate quantities.
The main idea of DC-DFT is to separate the errors in $\Delta E$ into a
functional error, which is present even with the exact density [19]:

$\Delta E_{\scriptscriptstyle\rm
F}=\tilde{E}[n]-E[n]=\tilde{E}_{\scriptscriptstyle\rm
XC}[n]-E_{\scriptscriptstyle\rm XC}[n],$ (2)

and the density-driven error is the remainder of $\Delta E$:

$\Delta E_{\scriptscriptstyle\rm D}=\tilde{E}[\tilde{n}]-\tilde{E}[n].$ (3)

In most KS-DFT calculations, $\Delta E_{\scriptscriptstyle\rm F}$ strongly
dominates $\Delta E$, implying that approximate KS densities are very good as
measured by the impact on energies. A very simply example for such case is the
total energy of the helium atom, whose energy computed from the PBE
functional[76] barely changes after the PBE density is replaced with the exact
one.[19] However, for a significant number of chemical domains (e.g., anions
and barrier heights), $\Delta E_{\scriptscriptstyle\rm D}$ can be much larger
than $\Delta E_{\scriptscriptstyle\rm F}$[22, 24]. A good example for such a
case is H- with the PBE functional, which gives excellent energies when
evaluated on the exact densities, whereas the self-consistent PBE density
cannot even bind two electrons for this anion.

## S2. Spotting and curing large density-driven errors

The following two questions are of key importance in DC-DFT: (i) How do we
spot cases where $\Delta E_{\scriptscriptstyle\rm D}$ is large?; (ii) What do
we do in such cases to reduce large density-driven errors? With access to
exact densities we can easily answer (i), as with those we can easily measure
$\Delta E_{\scriptscriptstyle\rm D}$ by Eq. 3. At the same time, we answered
(ii), as $\Delta E_{\scriptscriptstyle\rm D}$ entirely vanishes with exact
densities. However, exact densities are available only for small systems [19,
23], and if we always had access to exact energies and densities from highly
accurate wavefunction theories, we would not even bother with KS-DFT. Thus one
needs to find more practical ways to answer (i) and (ii). In relation to (i),
the following quantity has been introduced:

$\tilde{S}=\left|\tilde{E}[n^{\rm LDA}]-\tilde{E}[n^{\rm HF}]\right|,$ (4)

and is called density sensitivity. $\tilde{S}$ requires two nonempirical
densities: the HF densities which are typically overlocalized and the local
density approximation (LDA) densities which are typically delocalized.
$\tilde{S}$ serves as a practical measure of density sensitivity of a given
reaction and approximate functional. For small molecules, $\tilde{S}$ greater
than the heuristic cutoff of 2 kcal/mol implies density sensitivity,
indicating that the calculation may suffer from a large $\Delta
E_{\scriptscriptstyle\rm D}$. What should we do then to reduce large $\Delta
E_{\scriptscriptstyle\rm D}$? In these cases, evaluating an approximate
functional on the HF in place of self-consistent densities will likely reduce
$\Delta E_{\scriptscriptstyle\rm D}$ and likely improve the functional’s
performance. This procedure, called HF-DFT, is the practice of evaluating an
XC approximation on the HF density and orbitals. It had been used a long
before DC-DFT was proposed [77, 78, 79, 80, 81] but only because HF densities
were more convenient than self-consistent densities. KS-DFT does not always
benefit from HF densities (e.g, cases where $\Delta E_{\scriptscriptstyle\rm
D}$ is small) and in Refs. [26, 31, 66] we discuss in more details formal and
practical (dis)advantages of HF-DFT over SC-DFT.

## S3. HF-DFT and related DC-DFT procedures

Following the idea that in some cases HF-DFT works better and SC-DFT in
others, DC(HF)-DFT has been proposed [31, 66]. It is a procedure that
discriminately uses HF densities, as DC(HF)-DFT becomes HF-DFT for cases that
are both density-sensitive ($\tilde{S}$ above a given cut-off value, 2
kcal/mol as discussed in Ref. [23].) and whose HF solution is not severely
spin-contaminated. Otherwise, DC(HF)-DFT reverts to SC-DFT. Since HF-DFT uses
the HF density as a proxy for the exact density, we only use it when there is
little or no spin contamination. We calculate the expectation values of the
spin-squared operator, $S^{2}$, and only use the HF density if the <$S^{2}$>
from the HF calculation deviates less than 10% from the exact <$S^{2}$> as
discussed in Refs. [31] and [82]. Otherwise, we use the self-consistent
density. As it combines the best of both of them, DC(HF)-DFT comes with a
range of advantages over both HF-DFT and SC-DFT as further detailed in Ref.
[31]. These advantages come at a small extra cost, as for DC(HF)-DFT we need
to run up to three distinct self-consistent cycles to obtain the three
densities (the SC density for a given functional and those from HF and LDA
needed to calculate $\tilde{S}$). While we consider DC(HF)-DFT the state-of-
the-art DC-DFT-based procedure, in the present work we want a simple DC
framework that can be applied easily and routinely. r2SCAN and SCAN are
general-purpose functionals and the ease of their use would be undermined if
tandem with DC(HF)-DFT, which would require always calculating $\tilde{S}$ and
possibly making adjustments to its cut-off value. For this reason and
encouraged by the very good performance of HF-DFT with SCAN-like
functionals[32, 42], we employ HF-DFT as a DC-DFT procedure throughout this
work. As said, constructing the robust and accurate HF-r2SCAN-DC4 is the
central objective of this work. While the resulting HF-r2SCAN-DC4 can be
routinely used by applying it to HF orbitals without ever needing to calculate
$\tilde{S}$ of a given reaction, the use of $\tilde{S}$ is vital for our
training of HF-r2SCAN-DC4. Specifically, we use density-sensitivities of the
training reactions to fit the D4 part of HF-r2SCAN. Further technical details
of this fitting procedure will be given in the next section.

## S4. Optimizing dispersion parameters

D4 stands for the generally applicable atomic-charge dependent London
dispersion correction term developed by Grimme and co-workers.[83]. It has 4
functional-dependent parameters $s_{6}$, $s_{8}$, $a_{1}$, and $a_{2}$.
Following Refs. [43] and [83], we set $s_{6}$ to unity as is common for
functionals that do not capture long-range dispersion interactions. We
optimized the $s_{8}$, $a_{1}$, and $a_{2}$ parameters by minimizing the mean
absolute error (MAE) for the density-insensitive GMTKN55 reactions by
following the DC-DFT ideas of Ref. [26]. However, the density-insensitive
reactions in GMTKN55 largely fall into two distinct parameter groups for HF-
r2SCAN: $s_{8}$ has a negative value for noncovalent interactions, but is
positive for the rest. The difference in MAE of density-insensitive cases
between those two groups is miniscule (below 0.01 kcal/mol). For example,
($s_{8}$,$a_{1}$,$a_{2}$)=(-0.20,0.07,6.50) gives 1.209 kcal/mol for the
density-insensitive MAE while (0.39,0.09,7.02) gives 1.210 kcal/mol. Such a
difference is not meaningful. Small changes in computational details such as
DFT grid information, two-electron operator fitting scheme, etc. changes the
values of the parameters, since reaction energy errors and density sensitivity
values can be changed by 0.01 kcal/mol with those changes. To eliminate this
ambiguity while ensuring accuracy in water interactions, we include the
density-insensitive water$\cdots$water pair interaction energy as a validation
set. The two most stable water hexamers, the prism and the cage, are used to
calculate the water$\cdots$water 2-body interaction energy error per dimer,
relative to CCSD(T)/CBS in Ref. [13]. We multiply its weight by 7 in our loss
function to produce a better defined minimum and regularize the result (if we
used 1, it has no effect; if we used 1000, we simply fit to this data). We can
rationalize this value by noting that the mean density sensitivity of these
pairs is 0.27 kcal/mol, which is about 1/7th of our density sensitivity
cutoff. The resulting values for the three parameters are: -0.36, 0.23, 5.23
for $s_{8}$, $a_{1}$, and $a_{2}$ each.

## S5. Additional results for the GMTKN55 database

In Table S1, we list MAE (kcal/mol) of SC-r2SCAN and HF-r2SCAN functional with
and without the dispersion correction for the chemically diverse GMTKN55
database [40]. def2-QZVPPD basis set is used.

| HF-r2SCAN | HF-r2SCAN-DC4 | SC-r2SCAN | SC-r2SCAN-D4
---|---|---|---|---
MAE (all) | 2.82 | 2.42 | 2.98 | 2.91
MoM (all) | 2.77 | 2.32 | 3.24 | 3.10
WTMAD-2 (all)† | 8.54 | 5.18 | 8.65 | 7.10
basic† | 4.19 | 4.09 | 4.92 | 4.89
react.† | 8.10 | 5.85 | 9.48 | 8.11
barriers† | 7.18 | 8.04 | 13.11 | 13.56
inter. NCI† | 13.27 | 4.88 | 10.96 | 6.77
intra. NCI† | 11.96 | 4.80 | 8.67 | 5.87
basic
W4 | 6.92 | 6.41 | 3.81 | 3.86
G21EA | 4.59 | 4.55 | 3.88 | 3.86
G21IP | 4.56 | 4.54 | 4.66 | 4.64
DIPCS10 | 4.42 | 4.34 | 5.18 | 5.14
PA26 | 1.85 | 1.80 | 2.44 | 2.40
SIE4x4 | 12.14 | 12.24 | 17.93 | 17.98
ALKBDE10 | 4.71 | 4.66 | 5.01 | 5.00
YBDE18 | 3.89 | 3.67 | 3.89 | 3.36
AL2X6 | 1.02 | 1.23 | 0.93 | 1.58
HEAVYSB11 | 5.60 | 4.32 | 3.91 | 3.15
NBPRC | 1.24 | 1.24 | 1.60 | 1.52
ALK8 | 1.73 | 1.77 | 2.73 | 2.88
RC21 | 2.82 | 2.39 | 4.57 | 4.95
G2RC | 4.24 | 4.56 | 5.38 | 5.55
BH76RC | 2.57 | 2.56 | 2.97 | 2.97
FH51 | 1.71 | 1.72 | 2.19 | 2.16
TAUT15 | 1.11 | 1.10 | 1.58 | 1.57
DC13 | 8.96 | 7.97 | 8.63 | 7.72
react.
MB16-43 | 10.48 | 10.79 | 12.59 | 14.08
DARC | 3.82 | 2.00 | 3.71 | 2.70
RSE43 | 0.96 | 0.96 | 1.55 | 1.51
BSR36 | 3.20 | 0.16 | 2.32 | 0.48
CDIE20 | 1.21 | 1.13 | 1.63 | 1.61
ISO34 | 1.52 | 1.36 | 1.36 | 1.29
ISOL24 | 4.34 | 3.01 | 4.96 | 4.10
C60ISO | 3.52 | 3.88 | 5.35 | 5.57
PArel | 1.14 | 1.15 | 1.55 | 1.54
barriers
BH76 | 2.85 | 2.90 | 6.87 | 6.98
BHPERI | 4.14 | 5.95 | 3.86 | 4.65
BHDIV10 | 3.49 | 3.89 | 5.98 | 6.11
INV24 | 1.35 | 1.34 | 1.22 | 1.14
BHROT27 | 0.62 | 0.64 | 0.76 | 0.76
PX13 | 4.74 | 5.03 | 8.75 | 8.83
WCPT18 | 2.73 | 3.18 | 5.81 | 5.99
inter. NCI
RG18 | 0.26 | 0.10 | 0.23 | 0.16
ADIM6 | 2.65 | 0.45 | 1.98 | 0.34
s22 | 1.55 | 0.42 | 1.18 | 0.24
S66 | 1.42 | 0.30 | 1.02 | 0.26
HEAVY28 | 0.71 | 0.37 | 0.52 | 0.30
WATER27 | 3.94 | 1.01 | 4.24 | 6.30
CARBHB12 | 0.63 | 0.60 | 0.88 | 1.06
PNICO23 | 0.71 | 0.29 | 0.64 | 0.76
HAL59 | 1.01 | 0.41 | 0.99 | 0.80
AHB21 | 0.63 | 0.68 | 1.15 | 1.35
CHB6 | 0.46 | 0.58 | 0.49 | 0.52
IL16 | 1.88 | 0.43 | 0.33 | 0.64
intra. NCI
IDISP | 6.83 | 1.58 | 10.67 | 7.21
ICONF | 0.30 | 0.23 | 0.32 | 0.29
ACONF | 0.54 | 0.16 | 0.38 | 0.18
Amino20x4 | 0.40 | 0.27 | 0.26 | 0.19
PCONF | 1.22 | 0.44 | 1.05 | 0.41
MCONF | 0.94 | 0.21 | 0.63 | 0.45
SCONF | 0.57 | 0.19 | 0.37 | 0.51
UPU23 | 1.18 | 0.38 | 0.95 | 0.41
BUT14DIOL | 0.40 | 0.14 | 0.14 | 0.23

Table S1:  Mean-absolute-error (MAE) of GMTKN55 for selected functionals for
the individual datasets in the GMTKN55 database. WTMAD-2 value from Ref. [40]
and mean of means (MoM) values are also noted. WTMAD-2 values for individual
GMTKN55 categories †(basic properties and reaction energies for small systems
(basic), reaction energies for large systems and isomerisation reactions
(react.), reaction barrier heights (barriers), intermolecular noncovalent
interactions (inter. NCI), and intramolecular noncovalent interactions (intra.
NCI)) are also added. All units are kcal/mol.

## S6. Additional results for water clusters

### i. MD generated dimer structures

The structures used in Figures 1(e) and 2(b), have been obtained from the
molecular dynamics (MD) simulation. The simulation is performed within the XTB
package[84] and the GFN-FF force field[85], enabling us to generate the
various dimer configurations. The total simulation time is 50 ps, while the
integration time step is 4.0 fs using a Berendsen thermostat at 298K in the
NVT ensemble. We use the SHAKE algorithm to constrain bonds, for all bonds
with 4 amu for the hydrogen atom mass. Then, we randomly selected 110
different configurations for the water$\cdots$water dimer and 80 for the
water$\cdots$Aspirin dimers. The reference interaction energies are then
calculated with DLPNO-CCSD(T)-F12/TightPNO method with the aug-cc-pvqz basis
set for water$\cdots$water dimers and aug-cc-pvtz basis set for
water$\cdots$Aspirin dimers.

### ii. Many-body expansion of the interaction energy

The interaction energy can be decomposed into 2-body, 3-body, etc. by using
the many-body expansion.[67] For example, the interaction energy of the water
hexamer can be divided into $K$-body contributions,

$E_{int}=E_{int}^{2-body}+E_{int}^{3-body}+\cdots+E_{int}^{6-body}$ (5)

where $E_{int}^{K-body}$ is the K-body interaction energy which can be
calculated from the total energy of the subcluster of the $N$-mer cluster:
[67]

$E_{int}^{K-body}=\sum_{i=1}^{K}(-1)^{K-i}\begin{pmatrix}N-i\\\
K-i\end{pmatrix}S_{tot}^{i}$ (6)

where $S_{tot}^{i}$ stands for the total energy summation of the $i$-th
monomer subcluster.

## S7. Additional results for the complexes with cytosine

In Figure 4, we plot the cytosine$\cdots$water and cytosine$\cdots$cytosine
interaction energy error plots, respectively. The 14 different
cytosine$\cdots$cytosine configurations and reference interaction energies are
from Ref. [56]. For the cytosine$\cdots$water interaction, we place two water
molecules around 14 different cytosine$\cdots$cytosine structures and
optimized the water molecular coordinates while fixing the
cytosine$\cdots$cytosine coordinates. B3LYP functional is used for the
geometry optimization. DLPNO-CCSD(T)-F12/aug-cc-pvqz with TightPNO is used as
a reference cytosine$\cdots$water interaction energy with the ORCA
package.[86]

## S8. Interactions including water

For the calculations shown in Figure 5(a), we combined energies of 145
reactions involving water. They include: the water hexamer isomerization in
Figure 1(b), the water binding energy of WATER27 dataset in Figure 1(c), the
water 20-mer isomerization in Figure 1(d), the water dimer stationary point
geometry interaction energy in Figure 2(a), and water$\cdots$small organic
molecule interaction energy in Figure S11.

Figure S1:  Density sensitivity plot of (left) Smith dimer configuration and
(right) MD generated dimer configuration in Figure 2.

Figure S2:  (left) Smith dimer interaction energy errors and (right) reference
dimer interaction energies of MD generated dimer corresponding to Figure 2.
Figure S3:  Interaction energy error plot of cytosine$\cdots$cytosine
compounds for HF-r2SCAN-DC4/aug-cc-pvqz and HF-SCAN. Reference
HF/7Z-CP+MP2/CBS(6,7)-CP+dCC(cc-pVTZ-F12) energies are from Ref. [56].  Figure
S4:  Binding energy error of WATER27 dataset for SC-r2SCAN and SC-r2SCAN-D4.
The x-axis indicates the reaction number in WATER27 and the detailed
information of reaction including geometries can be found in the GMTKN55
database.[40] D4 parameters of SC-r2SCAN-D4 are from Ref. [87].  Figure S5:
The K-body energy plot corresponding to Figure 3. For comparison, SC-SCAN and
HF-SCAN-DC4 is additionally plotted in the total subplot. def2qzvppd basis set
is used. Figure S6:  $K$-body percentage contributions to reference
interaction energies of water hexamers (Figure 3).  Figure S7:  Density
sensitivity $\tilde{S}$ value for isomer energies for the water 20-mers
corresponding to Figure 1(d).  Figure S8:  Aspirin$\cdots$water interaction
including SC-r2SCAN-D4. The x-axis is the oxygen$\cdots$oxygen distance
between oxygen in the water and the specified oxygen in Aspirin (see inset of
Fig .1(e)). aug-cc-pvtz basis set is used.  Figure S9:  (left) Water hexamer
isomer relative energy compared to its global minimum geometry prism
structure; (right) their $\tilde{S}$ value calculated from the r2SCAN
functional. Geometries are from Ref. [13]. MAEs of each functional are 0.61,
0.04, 0.13, 0.36, and 0.43 kcal/mol compared to the CCSD(T)/CBS from Ref.
[61]. The ordering is the same as the legend ordering. For HF-r2SCAN-DC4, we
used the dispersion parameters from Section S4, whilst for SC-r2SCAN-D4, we
took the D4 parameters from Ref. [87]. The aug-cc-pvqz basis set was used for
the calculations.  Figure S10:  The hexagon plot with WTMAD-2 (kcal/mol) for
all GMTKN55 and its categories for selected functionals. (See Ref. [40] for
the detailed description of the categories).  Figure S11:  Interaction energy
curve for two non-hydrogen bonded systems. Geometries are from the HB375
dataset described in Ref. [88]. The x-axis show the scaled distance factor,
$f$, which is used to make a translated vector $t$,
$t=\frac{v}{|v|}(f-1)r_{ref}$ where $v$ is the bond direction vector and
$r_{ref}$ is the distance between the hydrogen and the electron-donor atom.
Detailed information about the x-axis can be found in Ref. [88].