# Supporting information: A classical density functional theory for solvation
across length scales

Anna T. Bui  Yusuf Hamied Department of Chemistry, University of Cambridge,
Lensfield Road, Cambridge, CB2 1EW, United Kingdom  Stephen J. Cox
<EMAIL_ADDRESS>Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge, CB2 1EW, United Kingdom

(August 28, 2024)

###### Contents

  1. S1 Derivation of the theory
    1. S1.1 Separation of length scales
    2. S1.2 Functional for the slowly-varying density field
    3. S1.3 Functional for the rapidly-varying density field
    4. S1.4 Coarse-graining procedure
  2. S2 Square-gradient functional for the slowly-varying density field
    1. S2.1 Parameterization
    2. S2.2 Free liquid–vapor interface
    3. S2.3 Bulk compressibility
    4. S2.4 Adapting the reference slowly-varying density in the crossover regime
  3. S3 Simulation details for correlation functions
    1. S3.1 Simulation details
    2. S3.2 Interatomic potentials
    3. S3.3 Workflow for extracting the direct correlation function
    4. S3.4 Extracting solvation free energies of hard spheres
  4. S4 Coarse-graining procedure
    1. S4.1 Hard sphere splitting
    2. S4.2 Coarse-graining parameters
  5. S5 Numerical procedure for LCW-style cDFT minimisation
    1. S5.1 Self-consistent cycle
    2. S5.2 Solving for the slowly-varying density field
    3. S5.3 Solving for the rapidly-varying density field
    4. S5.4 External potentials
  6. S6 Comparison to mDFT
    1. S6.1 Theory
    2. S6.2 Numerical results
  7. S7 Measure of hydrophobicity and critical drying
    1. S7.1 Local compressibility
    2. S7.2 Water in contact with a planar hard wall
    3. S7.3 Water around hydrophobic spherical solutes
  8. S8 Temperature dependence of hydrophobic solvation
    1. S8.1 Parameterization
    2. S8.2 Entropic crossover

##  S1 Derivation of the theory

In classical density functional theory (cDFT), the grand potential of a fluid
at a fixed chemical potential $\mu$ in the presence of an external potential
$\phi(\bm{r})$ can be written as a functional of an average microscopic
density $\rho(\bm{r})$:

$\varOmega_{\phi}[\rho]=\mathcal{F}_{\mathrm{intr}}[\rho]-\mu\int\\!\\!\mathrm{d}\bm{r}\,\rho(\bm{r})+\int\\!\\!\mathrm{d}\bm{r}\,\rho(\bm{r})\phi(\bm{r}).$
(S1)

The intrinsic Helmholtz free energy functional $\mathcal{F}_{\mathrm{intr}}$
is independent of $\phi$ and can be split further into ideal and excess
contributions,

$\mathcal{F}_{\mathrm{intr}}[\rho]=\mathcal{F}^{(\mathrm{id})}_{\mathrm{intr}}[\rho]+\mathcal{F}^{(\mathrm{ex})}_{\mathrm{intr}}[\rho].$
(S2)

The ideal contribution $\mathcal{F}_{\mathrm{intr}}^{\mathrm{(id)}}$ is given
exactly as the intrinsic free energy of a non-interacting system of density
$\rho(\bm{r})$,

$\mathcal{F}_{\mathrm{intr}}^{\mathrm{(id)}}[\rho]=k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\rho(\bm{r})[\ln\big{(}\Lambda^{3}\rho(\bm{r})\big{)}-1],$
(S3)

where $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ is the temperature and
$\Lambda$ is the thermal de Broglie wavelength. The excess contribution
$\mathcal{F}_{\mathrm{intr}}^{\mathrm{(ex)}}$ accounts for correlations in the
fluid and is generally unknown. It acts as a generating functional for a
hierarchy of direct correlation functions. The one- and two-body direct
correlation functions of the fluid are defined as

$c^{(1)}(\bm{r})=-\beta\frac{\delta\mathcal{F}_{\mathrm{intr}}^{\mathrm{(ex)}}}{\delta\rho(\bm{r})},$
(S4)

and

$c^{(2)}(\bm{r},\bm{r}^{\prime})=-\beta\frac{\delta^{2}\mathcal{F}_{\mathrm{intr}}^{\mathrm{(ex)}}}{\delta\rho(\bm{r})\delta\rho(\bm{r}^{\prime})},$
(S5)

where $\beta=1/k_{\mathrm{B}}T$.

###  S1.1 Separation of length scales

We first suppose that an inhomogeneous slowly-varying density field
$\rho_{\mathrm{s}}(\bm{r})$ exists that acts as an appropriate reference for
the fluid of interest. Using thermodynamic integration along the linear
integration path

$\rho_{\alpha}(\bm{r})=\rho_{\mathrm{s}}(\bm{r})+\alpha\delta_{\mathrm{s}}\rho(\bm{r}),$
(S6)

where $\delta_{\mathrm{s}}\rho=\rho-\rho_{\mathrm{s}}$, the intrinsic excess
Helmholtz free energy of the fluid of interest can be expressed exactly as

$\mathcal{F}^{(\rm ex)}_{\rm intr}[\rho]=\mathcal{F}^{(\rm ex)}_{\rm
s,intr}[\rho_{\rm s}]-k_{\rm B}T\int\\!\\!\mathrm{d}\bm{r}\,c_{\rm
s}^{(1)}(\bm{r})\delta_{\rm s}\rho(\bm{r})-k_{\rm
B}T\int^{1}_{0}\\!\\!\mathrm{d}\alpha\,(1-\alpha)\int\\!\\!\mathrm{d}\bm{r}\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\rm
s}\rho(\bm{r})c_{\alpha}^{(2)}(\bm{r},\bm{r}^{\prime})\delta_{\rm
s}\rho(\bm{r}^{\prime}),$ (S7)

where $\mathcal{F}^{(\rm ex)}_{\rm s,intr}$ is the intrinsic excess Helmholtz
free energy, $c_{\rm s}^{(1)}$ is the one-body direct correlation function of
the slowly-varying density $\rho_{\mathrm{s}}$ and $c_{\alpha}^{(2)}$ is the
two-body direct correlation function of the fluid with the intermediate
density $\rho_{\alpha}$. If we approximate that $c_{\alpha}^{(2)}$ is equal to
the direct correlation function of the reference slowly-varying density
$c_{\mathrm{s}}^{(2)}$ at all values of $\alpha$, we arrive at an expansion
for $\mathcal{F}^{(\rm ex)}_{\rm intr}$ correct to second order in
$\delta_{\mathrm{s}}\rho$

$\mathcal{F}^{(\rm ex)}_{\rm intr}[\rho]=\mathcal{F}^{(\rm ex)}_{\rm
s,intr}[\rho_{\rm s}]-k_{\rm B}T\int\\!\\!\mathrm{d}\bm{r}\,c_{\rm
s}^{(1)}(\bm{r})\delta_{\rm s}\rho(\bm{r})-\frac{k_{\rm
B}T}{2}\int\\!\\!\mathrm{d}\bm{r}\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\rm
s}\rho(\bm{r})c_{\mathrm{s}}^{(2)}(\bm{r},\bm{r}^{\prime})\delta_{\rm
s}\rho(\bm{r}^{\prime})$ (S8)

Analogously, by introducing an auxiliary reference density
$\rho_{\mathrm{r}}$, an expansion for $\mathcal{F}^{(\rm ex)}_{\rm s,intr}$,
correct to second order in
$\delta_{\mathrm{r}}\rho_{\mathrm{s}}=\rho_{\mathrm{s}}-\rho_{\mathrm{r}}$, is
given as

$\mathcal{F}^{(\rm ex)}_{\rm s,intr}[\rho_{\rm s}]=\mathcal{F}^{(\rm ex)}_{\rm
r,intr}[\rho_{\rm r}]-k_{\rm B}T\int\\!\\!\mathrm{d}\bm{r}\,c_{\rm
r}^{(1)}(\bm{r})\delta_{\rm r}\rho_{\mathrm{s}}(\bm{r})-\frac{k_{\rm
B}T}{2}\int\\!\\!\mathrm{d}\bm{r}\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\rm
r}\rho_{\mathrm{s}}(\bm{r})c_{\mathrm{r}}^{(2)}(\bm{r},\bm{r}^{\prime})\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime})$
(S9)

The functional derivative of $\mathcal{F}^{(\rm ex)}_{\rm s,intr}$ with
respect to $\rho_{\rm s}$ gives

$c^{(1)}_{\mathrm{s}}(\bm{r})=c^{(1)}_{\mathrm{r}}(\bm{r})+\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}c^{(2)}_{\mathrm{r}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime}).$
(S10)

Upon substitution of Eq. S10 into Eq. S8, we obtain

$\begin{split}\mathcal{F}^{(\rm ex)}_{\rm intr}[\rho]=&\mathcal{F}^{(\rm
ex)}_{\rm s,intr}[\rho_{\rm s}]-k_{\rm B}T\int\\!\mathrm{d}\bm{r}\,c_{\rm
r}^{(1)}(\bm{r})\delta_{\rm s}\rho(\bm{r})\\\ &-\frac{k_{\rm
B}T}{2}\int\\!\mathrm{d}\bm{r}\\!\int\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\rm
s}\rho(\bm{r})c_{\rm s}^{(2)}(\bm{r},\bm{r}^{\prime})\delta_{\rm
s}\rho(\bm{r}^{\prime})-k_{\rm
B}T\int\\!\mathrm{d}\bm{r}\\!\int\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\rm
s}\rho(\bm{r})c_{\rm r}^{(2)}(\bm{r},\bm{r}^{\prime})\delta_{\rm r}\rho_{\rm
s}(\bm{r}^{\prime}).\end{split}$ (S11)

The final term in Eq. S11 acts to couple fluctuations in the full density with
those of the slowly varying density. The auxiliary reference fluid is at the
same chemical potential $\mu$ as the fluid of interest. Therefore, the excess
part of the intrinsic chemical potential for the auxiliary reference fluid is

$-k_{\mathrm{B}}Tc^{(1)}_{\mathrm{r}}(\bm{r})=\mu^{\mathrm{ex}}=\mu-
k_{\mathrm{B}}T\ln\big{(}\Lambda^{3}\rho_{\mathrm{r}}(\bm{r})\big{)}.$ (S12)

The expression for $\mathcal{F}^{(\rm ex)}_{\rm intr}$ becomes

$\begin{split}\mathcal{F}^{(\rm ex)}_{\rm intr}[\rho]=&\mathcal{F}^{(\rm
ex)}_{\rm s,intr}[\rho_{\rm
s}]+\mu\\!\int\\!\\!\mathrm{d}\bm{r}\,\delta_{\mathrm{s}}\rho(\bm{r})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\delta_{\mathrm{s}}\rho(\bm{r})\ln(\Lambda^{3}\rho_{\mathrm{r}}(\bm{r}))\\\
&-\frac{k_{\mathrm{B}}T}{2}\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{s}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r})\,\delta_{\mathrm{s}}\rho(\bm{r}^{\prime})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{r}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime}).\end{split}$
(S13)

Combining the ideal and excess contributions from Eq. S3 and Eq. S13 into Eq.
S2, the intrinsic Helmholtz free energy functional becomes

$\begin{split}\mathcal{F}_{\mathrm{intr}}[\rho]=&\,k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\rho(\bm{r})[\ln\big{(}\Lambda^{3}\rho(\bm{r})\big{)}-1]\\\
&+\mathcal{F}_{\mathrm{s,intr}}^{\mathrm{(ex)}}[\rho_{\mathrm{s}}]+\mu\\!\int\\!\\!\mathrm{d}\bm{r}\,\delta_{\mathrm{s}}\rho(\bm{r})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\delta_{\mathrm{s}}\rho(\bm{r})\ln\big{(}\Lambda^{3}\rho_{\mathrm{r}}(\bm{r})\big{)}\\\
&-\frac{k_{\mathrm{B}}T}{2}\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{s}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r})\,\delta_{\mathrm{s}}\rho(\bm{r}^{\prime})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{r}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime})\\\
&+\mathcal{F}_{\mathrm{s,intr}}^{\mathrm{(id)}}[\rho_{\mathrm{s}}]-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\rho_{\mathrm{s}}(\bm{r})[\ln\big{(}\Lambda^{3}\rho_{\mathrm{s}}(\bm{r})\big{)}-1],\end{split}$
(S14)

where in the last line, we have added and subtracted the intrinsic Helmholtz
free energy for an ideal gas of density $\rho_{\mathrm{s}}(\bm{r})$.
Rearranging, we obtain

$\begin{split}\mathcal{F}_{\mathrm{intr}}[\rho]=&\mathcal{F}_{\mathrm{s,intr}}[\rho_{\mathrm{s}}]+(\mu-
k_{\mathrm{B}}T)\\!\int\\!\\!\mathrm{d}\bm{r}\,\delta_{\mathrm{s}}\rho(\bm{r})\\\
&+k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\rho(\bm{r})\ln\left(\frac{\rho(\bm{r})}{\rho_{\mathrm{r}}(\bm{r})}\right)-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\rho_{\mathrm{s}}(\bm{r})\ln\left(\frac{\rho_{\mathrm{s}}(\bm{r})}{\rho_{\mathrm{r}}(\bm{r})}\right)\\\
&-\frac{k_{\mathrm{B}}T}{2}\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{s}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r})\,\delta_{\mathrm{s}}\rho(\bm{r}^{\prime})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{r}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime}).\end{split}$
(S15)

The grand potential functional defined by Eq. S1 is then given by

$\begin{split}\varOmega_{\phi}[\rho]=&\mathcal{F}_{\mathrm{s,intr}}[\rho_{\mathrm{s}}]-\mu\\!\int\\!\\!\mathrm{d}\bm{r}\rho_{\mathrm{s}}(\bm{r})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{r}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime})\\\
&+k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\left[\rho(\bm{r})\ln\left(\frac{\rho(\bm{r})}{\rho_{\mathrm{r}}(\bm{r})}\right)-\rho(\bm{r})\right]-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\left[\rho_{\mathrm{s}}(\bm{r})\ln\left(\frac{\rho_{\mathrm{s}}(\bm{r})}{\rho_{\mathrm{r}}(\bm{r})}\right)-\rho_{\mathrm{s}}(\bm{r})\right]\\\
&-\frac{k_{\mathrm{B}}T}{2}\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{s}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r})\,\delta_{\mathrm{s}}\rho(\bm{r}^{\prime})+\\!\int\\!\\!\mathrm{d}\bm{r}\,\phi(\bm{r})\rho(\bm{r}).\end{split}$
(S16)

To further separate the grand potential functional into contributions from the
slowly- and rapidly-varying components of the density field, we assume that
$c_{\rm r}^{(2)}$ can be split in a similar fashion, i.e.,

$c_{\rm r}^{(2)}(\bm{r},\bm{r}^{\prime})=c_{{\rm
r},0}^{(2)}(\bm{r},\bm{r}^{\prime})+c_{{\rm
r},1}^{(2)}(\bm{r},\bm{r}^{\prime}),$ (S17)

where $c_{{\rm r},0}^{(2)}(\bm{r},\bm{r}^{\prime})\simeq 0$ for
$|\bm{r}-\bm{r}^{\prime}|>\ell$, where $\ell$ is a molecular length scale.
Conversely, $c^{(1)}_{{\rm r},1}$ is slowly varying over length scales
comparable to $\ell$. By substituting Eq. S17 into Eq. S16 and separating out
the slowly-varying components of the density field into another functional
$\varOmega_{\mathrm{s}}$, we arrive at

$\begin{split}\varOmega_{\phi}[\rho]=&\varOmega_{\mathrm{s}}[\rho_{\mathrm{s}}]+\\!\int\\!\\!\mathrm{d}\bm{r}\,\phi(\bm{r})\rho(\bm{r})\\\
&+k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\left[\rho(\bm{r})\ln\left(\frac{\rho(\bm{r})}{\rho_{\mathrm{r}}(\bm{r})}\right)-\rho(\bm{r})\right]-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\left[\rho_{\mathrm{s}}(\bm{r})\ln\left(\frac{\rho_{\mathrm{s}}(\bm{r})}{\rho_{\mathrm{r}}(\bm{r})}\right)-\rho_{\mathrm{s}}(\bm{r})\right]\\\
&-\frac{k_{\mathrm{B}}T}{2}\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\mathrm{s}}\rho(\bm{r})\,c^{(2)}_{\mathrm{s}}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{s}}\rho(\bm{r}^{\prime})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\mathrm{s}}\rho(\bm{r})\,c^{(2)}_{\mathrm{r},0}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime}),\\\
\end{split}$ (S18)

where

$\varOmega_{\mathrm{s}}[\rho_{\mathrm{s}}]=\mathcal{F}_{\mathrm{intr,s}}[\rho_{\mathrm{s}}]-\mu\\!\int\\!\\!\mathrm{d}\bm{r}\rho_{\mathrm{s}}(\bm{r})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\mathrm{s}}\rho(\bm{r})\,c^{(2)}_{\mathrm{r},1}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime}).$
(S19)

###  S1.2 Functional for the slowly-varying density field

Since $\rho_{\mathrm{s}}$ is slowly varying over the length scale comparable
to $\ell$, we can take the van der Waals (square–gradient) functional for the
intrinsic Helmholtz free energy:

$\mathcal{F}_{\rm s,intr}[\rho_{\rm
s}]=\int\\!\mathrm{d}\bm{r}\,\left[f_{0}+f_{2}|\nabla\rho_{\rm
s}(\bm{r})|^{2}\right].$ (S20)

In general, both $f_{0}$ and $f_{2}$ are functions of
$\rho_{\mathrm{s}}(\bm{r})$. $f_{0}$ is the local Helmholtz free energy
density. Taking the approximation $f_{2}\approx\frac{1}{2}m$, where $m$ is a
positive constant independent of density, we introduce the local grand
potential density

$\displaystyle\omega\big{(}\rho_{\mathrm{s}}(\bm{r});\mu\big{)}$
$\displaystyle=f_{0}\big{(}\rho_{\rm s}(\bm{r})\big{)}-\mu\rho_{\rm s},$ (S21)
$\displaystyle=\omega_{\mathrm{coex}}(\rho_{\mathrm{s}})-\rho_{\mathrm{s}}\delta\mu,$
(S22)

where $\delta\mu=\mu-\mu_{\mathrm{coex}}$ describes the difference in chemical
potential from coexistence. In principle, $\omega_{\rm coex}$ can be obtained
from the equation of state of the bulk fluid; to our knowledge, this is not
known for the different intermolecular potentials for water that we explore in
this work. We therefore adopt a simple double-well form that broadly describes
the physics of liquid-vapor coexistence:

$\omega_{\rm coex}(\rho_{\rm
s})=\frac{C}{2}\left(\rho_{\mathrm{s}}-\rho_{\mathrm{l}}\right)^{2}\left(\rho_{\mathrm{s}}-\rho_{\mathrm{v}}\right)^{2}+\frac{D}{4}\left(\rho_{\mathrm{s}}-\rho_{\mathrm{l}}\right)^{4}\left(\rho_{\mathrm{s}}-\rho_{\mathrm{v}}\right)^{4}h(\rho_{\rm
s}-\rho_{\rm v})h(\rho_{\rm l}-\rho_{\rm s}),$ (S23)

where $\rho_{\mathrm{l}}$ and $\rho_{\mathrm{v}}$ are the liquid and vapor
densities at coexistence, locating the minima of
$\omega_{\mathrm{coex}}(\rho_{\mathrm{s}})$, and $h(\rho_{\rm s})$ is the
Heaviside step function. We will discuss the parameterization of $C$, $D$ and
$m$ in Section S2.

Equation S19 now reads

$\begin{split}\varOmega_{\rm s}[\rho_{\rm
s}]=\int\\!\mathrm{d}\bm{r}\,\left[\omega\big{(}\rho_{\rm
s}(\bm{r})\big{)}+\frac{m}{2}|\nabla\rho_{\rm
s}(\bm{r})|^{2}\right]-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\mathrm{s}}\rho(\bm{r})\,c^{(2)}_{\mathrm{r},1}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime}).\end{split}$
(S24)

Minimization of $\varOmega_{\rm s}$ with respect to $\rho_{\rm s}$ for a given
$\rho$ leads to

$\omega^{\prime}\big{(}\rho_{\rm s}(\bm{r})\big{)}-m\nabla^{2}\rho_{\rm
s}(\bm{r})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{r},1}(\bm{r},\bm{r}^{\prime})\,[\rho(\bm{r}^{\prime})-\rho_{\mathrm{s}}(\bm{r}^{\prime})]+k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{r},1}(\bm{r},\bm{r}^{\prime})\,[\rho_{\mathrm{s}}(\bm{r}^{\prime})-\rho_{\mathrm{r}}(\bm{r}^{\prime})]=0,$
(S25)

where $\omega^{\prime}$ indicates a partial differentiation with respect to
$\rho_{\mathrm{s}}$. Setting $\rho_{\rm r}=\rho_{\rm s}$, we recover Eq. 13 in
the main article

$\omega^{\prime}\big{(}\rho_{\rm s}(\bm{r})\big{)}-m\nabla^{2}\rho_{\rm
s}(\bm{r})+\psi(\bm{r})=0,$ (S26)

where we have defined

$\psi(\bm{r})=-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{r},1}(\bm{r},\bm{r}^{\prime})\,[\rho(\bm{r}^{\prime})-\rho_{\mathrm{s}}(\bm{r}^{\prime})].$
(S27)

This is formally similar to the “unbalancing potential” introduced by Weeks
and co-workers [1, 2].

###  S1.3 Functional for the rapidly-varying density field

We now minimize $\varOmega_{\phi}$ with respect to $\rho$ for fixed $\rho_{\rm
s}$ and in this case, we find

$\rho(\bm{r})=\rho_{\rm
r}(\bm{r})\exp\left[-\beta\phi(\bm{r})+\int\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\rm
s}(\bm{r},\bm{r}^{\prime})\,\delta_{\rm
s}\rho(\bm{r}^{\prime})+\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{r},0}(\bm{r},\bm{r}^{\prime})\,\delta_{\mathrm{r}}\rho_{\mathrm{s}}(\bm{r}^{\prime})\right].$
(S28)

Again setting $\rho_{\rm r}=\rho_{\rm s}$, we obtain Eq. 11 in the main
article:

$\rho(\bm{r})=\rho_{\rm
s}(\bm{r})\exp\left[-\beta\phi(\bm{r})+\int\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\rm
s}(\bm{r},\bm{r}^{\prime})\,\delta_{\rm s}\rho(\bm{r}^{\prime})\right].$ (S29)

For $c^{(2)}_{\rm s}$, we make an approximation

$c^{(2)}_{\mathrm{s}}(\bm{r},\bm{r}^{\prime})\approx
c^{(2)}_{\mathrm{u}}(|\bm{r}-\bm{r}^{\prime}|)\rho_{\mathrm{s}}(\bm{r})\rho_{\mathrm{s}}(\bm{r}^{\prime})\rho_{\mathrm{u}}^{-2},$
(S30)

which is exact in the limits of homogeneity and low density, and interpolates
smoothly in between.

###  S1.4 Coarse-graining procedure

The unbalancing potential as written in Eq. S27 is an integral of the rapidly-
varying component of the density field $\delta_{\rm s}\rho$ over a slowly-
varying function $c^{(2)}_{\mathrm{r,1}}$. This can be interpreted as a
coarse-graining of the density field. More concretely, we define a coarse-
graining procedure

$\overline{\rho}(\bm{r})=\frac{\int\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{{\rm
r},1}(\bm{r},\bm{r}^{\prime})\rho(\bm{r}^{\prime})}{\int\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{{\rm
r},1}(\bm{r},\bm{r}^{\prime})}.$ (S31)

Equation S26 then becomes

$\omega^{\prime}(\rho_{\rm s}(\bm{r}))=m\nabla^{2}\rho_{\rm
s}(\bm{r})+2a(\bm{r})[\overline{\rho}(\bm{r})-\overline{\rho}_{\rm
s}(\bm{r})],$ (S32)

with

$2a(\bm{r})=k_{\rm B}T\int\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{{\rm
r},1}(\bm{r},\bm{r}^{\prime}).$ (S33)

We can approximate the reference two-body direct correlation function
$c^{(2)}_{\mathrm{r}}$ with its translationally invariant counterpart of the
uniform fluid $c^{(2)}_{\mathrm{u}}$, such that

$\overline{\rho}(\bm{r})\approx\frac{\int\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{{\rm
u},1}(|\bm{r}-\bm{r}^{\prime}|)\rho(\bm{r}^{\prime})}{\int\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{{\rm
u},1}(|\bm{r}-\bm{r}^{\prime}|)},$ (S34)

and

$2a(\bm{r})\approx 2a=k_{\rm B}T\int\\!\mathrm{d}\bm{r}\,c^{(2)}_{{\rm
u},1}(r).$ (S35)

As an interesting aside, for simple fluids like a Lennard–Jones fluid, one
might consider using the random phase approximation (RPA), [3]

$c^{(2)}_{{\rm u}}(|\bm{r}-\bm{r}^{\prime}|)\approx c^{(2)}_{{\rm
u},0}(|\bm{r}-\bm{r}^{\prime}|)-\beta u_{1}(|\bm{r}-\bm{r}^{\prime}|),$ (S36)

where $u_{1}$ denotes the long-ranged attractive part of the interparticle
potential (through Weeks–Chandler–Anderson splitting [4]). In this case, our
coarse-graining reduces to

$\overline{\rho}(\bm{r})\approx\frac{\int\\!\mathrm{d}\bm{r}^{\prime}\,u_{1}(|\bm{r}-\bm{r}^{\prime}|)\rho(\bm{r}^{\prime})}{\int\\!\mathrm{d}\bm{r}\,u_{1}(r)},$
(S37)

with

$2a=\int\\!\mathrm{d}\bm{r}\,u_{1}(r).$ (S38)

This is identical to the coarse-graining procedure specified by Weeks in the
context of local molecular field theory [2].

For complex fluids such as water, making a reasonable range separation
described by Eq. S17 is not straightforward as for simple liquids. Therefore,
for practical purposes, we use a Gaussian weight of width $\lambda$ to coarse-
grain the density,

$\overline{\rho}(\bm{r})\approx\frac{1}{(2\pi\lambda^{2})^{3/2}}\int\\!\mathrm{d}\bm{r}^{\prime}\,\exp\left(-\frac{|\bm{r}-\bm{r}^{\prime}|^{2}}{2\lambda^{2}}\right)\rho(\bm{r}^{\prime}),$
(S39)

where $\lambda$ and $a$ are parameters that need to be determined. We will
provide a basis to justify these parameters in Section S4.

##  S2 Square-gradient functional for the slowly-varying density field

###  S2.1 Parameterization

The form of $\omega_{\mathrm{coex}}(\rho_{\mathrm{s}})$ specified by Eq. S23
is approximate. This is evident from its symmetric form; the curvatures at
$\rho_{\rm v}$ and $\rho_{\rm l}$, as determined by $C$, are identical. Or, in
other words, the compressibilities of the vapor and liquid phases are
identical. For our purposes, we are more concerned with the compressibility in
the liquid phase, and we choose $C$ accordingly. To be precise, denoting the
bulk compressibility as $\chi_{\rm u}$, we have

$\chi_{\rm u}=\left(\frac{\delta\rho_{\rm
u}}{\delta\mu}\right)_{T}=\frac{1}{\omega^{\prime\prime}(\rho_{\mathrm{l}})}=\frac{1}{C(\rho_{\rm
l}-\rho_{\rm v})}.$ (S40)

Following the standard compressibility relation [3]
$\chi_{\mathrm{u}}=\beta\rho_{\mathrm{u}}S(0)$, where $S(0)$ is the
$k\rightarrow 0$ limit of the liquid structure factor $S(k)$, we can write

$C=\frac{k_{\mathrm{B}}T}{\rho_{\rm
u}S(0)(\rho_{\mathrm{l}}-\rho_{\mathrm{v}})^{2}}.$ (S41)

In practice, to calculate $C$ using Eq. S41, for each water model at a
particular $T$, we perform a simulation of the bulk liquid at density
$\rho_{\mathrm{u}}=\rho_{\mathrm{l}}$ and extract $S(0)$ as described in
Section S3. This ensures consistency in $\chi_{\rm u}$ between our rapidly-
and slowly-varying functionals. We note that using an asymmetric form of
square gradient theory would allow one to distinguish the compressibilities of
the liquid and vapor phases, as well as differences in density gradients on
either side of the free liquid-vapor interface [5]. Initial investigations
along these lines, however, proved challenging, and we therefore settled on
the simple form described above.

If we set $D=0$, then Eq. S23 reduces to a standard approximate form
encountered in textbooks [3]. If we then choose $m$ to attain the surface
tension of the free liquid-vapor interface, we have found that the resulting
density profiles compare poorly to results from molecular simulations. While
$\omega_{\rm coex}$ describes properties of the bulk fluid, in order to
simultaneously obtain $\gamma$ and reasonable density profiles, we have taken
the pragmatic approach of introducing the second term in Eq. S23. We then
parameterize $D$ and $m$ by the following procedure:

  1. (i)

We start with $D=0$ and find a value for $m$ such that the resulting
functional gives the correct liquid–vapor surface tension $\gamma$.

  2. (ii)

If the resulting density profile for the free liquid–vapor interface is too
sharp compared to molecular simulation, this means that the contribution to
$\varOmega_{\mathrm{s}}$ from the local term is too large relative to that
from the gradient term. We therefore reduce the relative contribution from
$\omega_{\mathrm{coex}}(\rho_{\mathrm{s}})$ by guessing a negative value of
$D$ to decrease the barrier height in the region
$\rho_{\mathrm{v}}<\rho_{\mathrm{s}}<\rho_{\mathrm{l}}$. If the resulting
density profile is too wide, we instead guess a positive value of $D$.

  3. (iii)

With this new value of $D$, we find a new value for $m$ that gives the correct
surface tension $\gamma$.

  4. (iv)

We return to step (ii) if the resulting density profile still deviates
significantly from the simulation data.

We note that, since the interfacial width is more sensitive to $D$ than $m$,
we do not need to adjust these parameters too much. As a remark, were we to
know the equation of state for the water models that we investigate, then we
would be able to determine $\omega_{\rm coex}$ exactly. We would then focus
our efforts entirely on parameterizing $m(\rho_{\rm s})$ on properties of the
free liquid-vapor interface; in the absence of such data, the approach we
describe above provides a practical means to parameterizing $\omega_{\rm
coex}$ for our purposes.

In the main paper, we showed results using three different water models: SPC/E
[6], RPBE-D3 [7, 8] and mW [9] (further details in Section S3). To
parameterize the functional for each water model, we use data from simulations
of the free liquid–vapor interface and the uniform fluid at the liquid
coexistence density at $300\,\mathrm{K}$ for each water model. For the mW
model, we also investigated $T=426\,\mathrm{K}$. For SPC/E water, the
liquid–vapor coexistence properties at $300\,\mathrm{K}$ have been determined
by Vega and Miguel from direct coexistence molecular dynamics (MD) simulations
[10]. For water described with the electronic density functional RPBE-D3, the
liquid–vapor coexistence curve (binodal) has been determined by Schienbein and
Marx [11] through _ab initio_ Gibbs ensemble Monte Carlo simulations. However,
the surface tension was not computed. Wohlfalhrt _et al._ trained a neural
network potential [12] using energies and forces computed at the RPBE-D3 level
and were able to reproduce the coexistence curve from Schienbein and Marx, in
addition to determining the surface tensions at different temperatures using
direct coexistence MD simulations. We therefore use the liquid–vapor
coexistence properties at $300\,\mathrm{K}$ from Wohlfahrt _et al._ [12] For
the mW model, we use liquid–vapor coexistence properties determined by Coe _et
al._ from grand canonical Monte Carlo (GCMC) simulations [13]. Since the
critical temperature of mW ($T_{\mathrm{c}}=917\,\mathrm{K}$) is much higher
than that of real water ($T_{\mathrm{c}}=647\,\mathrm{K}$), we investigate
$T=426\,\mathrm{K}$ for mW as it gives the same scaled temperature as ambient
water $T/T_{\mathrm{c}}\approx 0.46$ (argued in Ref. 13). We have also
considered $T=300\,\mathrm{K}$ where the surface tension for mW is in better
agreement with the other water models. In instances where the interfacial
density profile is not given in the above references, we obtained it with
direct coexistence MD simulations of our own (in such cases, we checked for
consistency, e.g., by computing the surface tension). The bulk simulations at
constant density $\rho_{\rm u}=\rho_{\rm l}$ to extract the bulk
compressibility $\chi_{\rm u}$ are described in Section S3. We summarized the
parameters used in Table S1.

Water model | SPC/E | RPBE-D3 | mW | mW
---|---|---|---|---
$T$ [K] | 300 | 300 | 426 | 300
$\rho_{\rm l}$ [$\mathrm{nm^{-3}}$] | 33.234a | 30.06b | 32.203c | 33.405c
$\rho_{\rm v}$ [$\mathrm{nm^{-3}}$] |  $4.747\times 10^{-4}$a |  $1.04\times 10^{-3}$b |  $1.652\times 10^{-3}$c |  $1.337\times 10^{-5}$c
$\gamma$ [$\mathrm{mJ\,m^{-2}}$] | 63.6a |  $68\pm 7$b |  $55.6\pm 2$c | 65.8395c
$k_{\mathrm{B}}T\rho_{\mathrm{u}}^{-1}\chi_{\mathrm{u}}=S(0)$ | 0.0630 | 0.0598 | 0.0456 | 0.0265
$C$ [$\mathrm{J\,mol^{-2}\,nm^{9}}$] | 1.079 | 1.537 | 2.296 | 2.523
$D$ [$\mathrm{J\,mol^{-2}\,nm^{15}}$] | 30 | 90 | -15 | -15
$m$ [$\mathrm{kJ\,mol^{-2}\,cm^{3}\,\text{\AA}^{2}}$] | 1255 | 1500 | 1100 | 1148
Table S1: Parameters for the slowly varying functional. For the SPC/E, RPBE-D3
and mW water models, square gradient theory is parameterized based on
liquid–vapor coexistence at temperature $T=300\,\mathrm{K}$ (and
$T=426\,\mathrm{K}$ for mW). The first part of the table gives thermodynamic
properties of water. The second part gives parameters for the square–gradient
functional.
aDetermined by Vega and Miguel in Ref. 10.
bDetermined by Wohlfahrt _et al._ in Ref. 12.
cDetermined by Coe in Ref. 13. We have used data for $T=420\,\mathrm{K}$ for
$T=426\,\mathrm{K}$ since it is the closest state point tabulated, with
indicated uncertainties accounting for any discrepancy.

###  S2.2 Free liquid–vapor interface

The resulting local grand potential density and liquid–vapor interface
profiles from square gradient theory parameterized for each water model are
shown in Fig. S1.

Figure S1: Liquid–vapor interface from square gradient theory. The local grand
potential density at coexistence $\omega_{\mathrm{coex}}(\rho_{\rm s})$ given
in Eq. S22, and corresponding liquid–vapor interface profiles $\rho(z)$ given
in Eq. S23 are shown in (a) and (b), respectively. The square gradient theory
is parameterized for the water models and state points indicated in the legend
in (b). The function $\omega_{\mathrm{coex}}(\rho_{\rm s})$ has two minima at
$\rho_{\mathrm{v}}$ and $\rho_{\mathrm{l}}$ corresponding to the vapor and
liquid phases, respectively. The density profiles for SPC/E and RPBE-D3 are
taken from Ref. 12 and Ref. 14, respectively. For mW, the density profiles are
from our own direct coexistence simulations. We note that the density of the
liquid phase for RPBE-D3 deviates the most from the experimental liquid
density ($33.3\,\mathrm{nm^{-3}}$), which is a shortcoming of this electronic
density functional well-documented in the literature [15, 11, 16].

###  S2.3 Bulk compressibility

Figure S2: Liquid bulk compressibility from square gradient theory. From the
simple square gradient form given in Eq. S22, an increase in
$\delta\mu=\mu-\mu_{\mathrm{coex}}$ stabilizes the liquid phase further and
shifts the liquid minimum to a higher density in a linear fashion, as shown in
(a) for SPC/E water at $300\,\mathrm{K}$. The circles indicate the global
minimum of each curve, locating the stable liquid density at
$\rho_{\mathrm{u}}$. Note that for $\beta\delta\mu\approx 10^{-3}$ the
deviation from coexistence is relatively small. The linearity between
$\rho_{\rm u}$ and $\delta\mu$ is depicted in the inset of (b); the gradient
of this line gives the bulk compressibility $\chi_{\mathrm{u}}$, as determined
by the square gradient functional. As seen in (b), as we move away from
coexistence, this estimate for $\chi_{\mathrm{u}}$ deviates more from values
obtained from simulations (detailed in Section S3), which are indicated by the
circles (error bars smaller than the symbols). In the region that is shaded
gray ($\rho_{\mathrm{u}}<\rho_{\mathrm{l}}$), the liquid phase is metastable.
The domain of applicability of the theory ($0<\beta\delta\mu<0.1$) is shaded
green.

The ensemble to naturally work in for cDFT is the grand canonical ensemble.
Water under ambient conditions lies close to its liquid–vapor coexistence [17]
with a chemical potential deviation from coexistence $\beta\delta\mu\approx
10^{-3}$. For a bulk fluid at chemical potential $\mu$, the corresponding
uniform density $\rho_{\mathrm{u}}$ is determined by the global minimum of the
local free energy density $\omega(\rho_{\rm s})$. For the form that we
adopted, as given in Eq. S22, the term linear in $\rho_{\rm s}$ adds a bias to
the double-well form of $\omega_{\mathrm{coex}}(\rho_{\rm s})$, making the
liquid phase more stable at higher $\rho_{\rm s}$. This behavior is
demonstrated in Fig. S2(a) for the square gradient theory parameterized on
SPC/E water; note that we are using large values of $\beta\delta\mu$ to
exaggerate the behavior. In Fig. S2(b) we plot the bulk compressibility of the
bulk liquid determined from simulations, and from square gradient theory.
Owing to the approximate form of $\omega_{\rm coex}(\rho_{\mathrm{s}})$,
agreement with the simulation results is limited to $0<\beta\delta\mu\lesssim
0.1$. In Section S7, we also show that this range of $\delta\mu$ is where the
contact theorem (a statistical mechanical sum rule) for a fluid at a planar
hard wall is also obeyed.

In the following, we determine the uniform liquid density $\rho_{\mathrm{u}}$
from the square gradient form corresponding to a chemical potential deviation
of $\beta\delta\mu\approx 10^{-3}$ and the bulk compressibility
$\chi_{\mathrm{u}}$ for each water model, as summarized in Table S2.

Water model | SPC/E | RPBE-D3 | mW | mW
---|---|---|---|---
$T$ [K] | 300 | 300 | 426 | 300
$\rho_{\rm u}$ [$\mathrm{nm^{-3}}$] | 33.236 | 30.079 | 32.204 | 33.406
Table S2:  Bulk liquid density and compressibility from square gradient
theory. The table gives the stable liquid density $\rho_{\mathrm{u}}$ that
corresponds $\beta\delta\mu\approx 10^{-3}$ determined from the minimum of
$\omega(\rho_{\mathrm{s}})$ for each water model and state point considered.

###  S2.4 Adapting the reference slowly-varying density in the crossover
regime

By construction, the theory captures the solvation behavior in the limits of
small solutes and macroscopic planar surfaces most accurately, with a gradual
crossover between these limits. The exact behavior in the crossover regime,
however, is inevitably more challenging to capture. To explore this aspect, we
consider the solvation of hard spheres of increasing radius $R$ in SPC/E water
at 300 K where simulation data up to $R=1\,\mathrm{nm}$ are available [18], as
shown in Fig. S3.

Figure S3: Adapting the reference slowly-varying density in the crossover
regime. In (a), $\varOmega_{\mathrm{solv}}/4\pi R^{2}$ for solvation of hard
sphere solutes in SPC/E water at 300 K and $\beta\delta\mu=10^{-3}$ predicted
from the theory using a constant $m$ (dotted line) and using
$m(\rho^{*}_{\mathrm{s}})$ (solid line) where $\rho^{*}_{\mathrm{s}}$ is a
measure of how depleted the density profile is. The linear form of
$m(\rho^{*}_{\mathrm{s}})$ employed is shown in the inset. For both
approaches, $\varOmega_{\rm solv}/4\pi R^{2}\sim\gamma$, indicated by the open
circle. In (b), $\rho(r)$ and $\rho_{\mathrm{s}}(r)$ for various solute sizes
are shown in solid and dashed lines respectively. The case of the solvent in
contact with a hard wall is shown in the inset. The values for
$\rho^{*}_{\mathrm{s}}$ are marked with open circles.

Similar to previous studies with LCW theory [1, 19, 20, 14], the functional
using the square gradient theory as parameterized in the previous subsection
with a constant value for $m$ predicts that $\varOmega_{\mathrm{solv}}/4\pi
R^{2}$ first increases up to a value higher than the liquid–vapor surface
tension $\gamma$ at $R\approx 0.7\,\mathrm{nm}$, before slowly decreasing to
approach $\gamma$ from above as $R\to\infty$. (More precisely,
$\varOmega_{\mathrm{solv}}/4\pi R^{2}\sim\gamma+\gamma_{\rm wv}$, where
$\gamma_{\rm wv}$ is the surface tension between the vapor phase and the wall.
Owing the small vapor pressure of water, it is common to neglect this
contribution.) Simulations, however, predict that
$\varOmega_{\mathrm{solv}}/4\pi R^{2}$ does not exceed the surface tension up
to $R=1\,\mathrm{nm}$. This observation suggests that in the crossover regime,
the best slowly-varying reference density slightly deviates from what is given
by $\varOmega_{\mathrm{s}}[\rho_{\mathrm{s}}]$ parameterized with constant
$m$. In principle, the coefficient of the square-gradient term could depend on
density, i.e., $m(\rho_{\rm s})$. [Note: this would introduce a term
proportional to $m^{\prime}(\rho_{\rm s})|\nabla\rho_{s}|^{2}$ upon
minimization.] Such a parameterization is not straightforward. For practical
purposes, we instead treat $m$ as a function of a characteristic feature of
$\rho_{\rm s}^{\ast}$. Specifically, $\rho^{*}_{\mathrm{s}}$ is the minimum
value of the density at positions around the solute where $\rho(\bm{r})$
intercepts $\rho_{\mathrm{s}}(\bm{r})$. Since $\rho^{*}_{\mathrm{s}}$ gives an
indication of how depleted the density profile is, adapting a constant $m$ to
$m(\rho^{*}_{\mathrm{s}})$ allows us to quantitatively describe the simulation
data. Here, we employ a linear dependence for $m(\rho^{*}_{\mathrm{s}})$,

$m(\rho^{*}_{\mathrm{s}})=m+t\rho^{*}_{\mathrm{s}}.$ (S42)

This enforces the constant $m$ value is recovered when $\rho_{\rm
s}^{\ast}\approx\rho_{\rm v}$. For SPC/E water, we have used
$t=-20323\,\mathrm{kJ\,mol^{-2}\,cm^{3}\,\text{\AA}^{5}}$. Note that the
solvation behaviors for spheres with $R<0.4\,\mathrm{nm}$ and at a planar wall
are not sensitive to this change. For hard spheres, $\rho^{*}_{\mathrm{s}}$
amounts to the contact density; the definition we describe facilitates its
application to soft-core solutes.

We emphasize that even with a constant $m$, the results for both the solvation
free energy and density are still semi-quantitative. Therefore, in cases where
we lack simulation data in the crossover regime, we use a constant value of
$m$ parameterized on the free liquid-vapor interface.

##  S3 Simulation details for correlation functions

In this section, we describe the procedure followed to obtain the spherically
symmetric direct correlation function of a bulk fluid from simulation.

###  S3.1 Simulation details

For each interatomic potential, MD simulations of bulk liquid water at ambient
conditions were performed. All simulations were carried out with the LAMMPS
simulation package [21]. Cubic simulation boxes with periodic boundary
conditions were used. For the structure factor calculation, the systems
contained $13824$ water molecules for SPC/E and RPBE-D3, $15000$ water
molecules in the case of the mW model. Dynamics were propagated using the
velocity Verlet algorithm with a time-step of $1\text{\,}\mathrm{fs}$. The
temperature was maintained at $T=$300\text{\,}\mathrm{K}$$ or
$T=$426\text{\,}\mathrm{K}$$ using the CSVR thermostat [22]. The system was
equilibrated for $200\,\mathrm{ps}$ and production runs were performed for at
least $2\,\mathrm{ns}$ in the NVT ensemble.

###  S3.2 Interatomic potentials

Three different interatomic potentials for water with different levels of
complexity are considered:

  1. (i)

SPC/E: This simple point charge water model [6] contains a Lennard-Jones
centre on the oxygen atom and partial charges on the hydrogen and oxygen
atoms. The geometries of water molecules were constrained using the RATTLE
algorithm [23]. All Lennard-Jones interactions were truncated and shifted at 1
nm. Long-ranged electrostatic interactions were evaluated using
particle–particle particle–mesh Ewald summation [24] (using a 1 nm cut off in
real space) such that the RMSE in the forces was a factor of $10^{5}$ smaller
than the force between two unit charges separated by a distance of 0.1 nm
[25].

  2. (ii)

mW: The monatomic water potential [9] is a coarse-grained model of water in
which there are no explicit hydrogens. It takes the form of a Stillinger–Weber
potential [26] and the tetrahedrality is achieved by three-body terms. For
this single-site water model, the fluid density is only described by a singlet
particle density so no angular component is neglected in our theory.

  3. (iii)

RPBE-D3: An ab initio approach to model water can be achieved through quantum
mechanical calculations. For the underlying electronic structure, the Revised
Perdew–Burke–Ernzerhof (RPBE) generalized gradient approximation for the
exchange–correlation functional [7] supplemented by Grimme’s D3 dispersion
correction [8] is used. Since thorough sampling of the liquid state’s phase
space at this level of theory is prohibitively expensive, we perform
simulations using a surrogate neural network potential [27] trained on forces
and energies from ab initio simulation trajectories as given in Ref. 12, using
LAMMPS [21] with an interface to the n2p2 package[28].

###  S3.3 Workflow for extracting the direct correlation function

We denote the 3D spatial Fourier transform of function $f(\bm{r})$ as

$\hat{f}(\bm{k})=\int\\!\\!\mathrm{d}\bm{r}\,f(\bm{r})\exp(-i\bm{k}\cdot\bm{r}),$
(S43)

and the inverse Fourier transform of $\hat{f}(\bm{k})$ as

$f(\bm{r})=\frac{1}{(2\pi)^{3}}\int\\!\\!\mathrm{d}\bm{k}\,\hat{f}(\bm{k})\exp(i\bm{k}\cdot\bm{r}).$
(S44)

For a radially symmetric function $f(r)$, the corresponding expressions become

$\hat{f}(k)=4\pi\int^{\infty}_{0}\\!\\!\mathrm{d}r\,r^{2}\frac{\sin(kr)}{kr}f(r),$
(S45)

and

$f(r)=\frac{1}{(2\pi)^{3}}\int^{\infty}_{0}\\!\\!\mathrm{d}k\,\frac{4\pi
k}{r}\sin(kr)\hat{f}(k).$ (S46)

For a bulk system containing $N$ water molecules at density
$\rho_{\mathrm{u}}$, we extracted the radial distribution function $g(r)$

$g(r)=\frac{1}{N^{2}}\sum^{N}_{i,j}\left\langle\delta(r-|\bm{r}_{i}-\bm{r}_{j}|)\right\rangle,$
(S47)

by generating a histogram of intermolecular distances and appropriate
normalization. Note that we considered the positions of the oxygen as the
positions of the water molecule in the case of SPC/E and RPBE-D3 water. The
resulting radial distribution functions for the different water models are
shown in Fig. S4(a). The high-wavevector part of the structure factor
($k>1\,\text{\AA}^{-1}$) are computed from [3]

$S(k)=1+4\pi\rho_{\mathrm{u}}\int^{\infty}_{0}\\!\\!\mathrm{d}r\,\frac{\sin(kr)}{k}\left[g(r)-1\right].$
(S48)

Figure S4: Bulk structure of different water models from simulation. The
radial distribution functions $g(r)$ are shown in (a), with the structure
factor $S(k)$ given in the inset. The direct correlation function determined
from simulation to be used as input to cDFT is shown in (b). The water models
and temperatures considered are denoted in the legend in (b).

In the low-wavevector limit ($k<1\,\text{\AA}^{-1}$), we also calculate the
structure factor directly from the simulation trajectories using[29]

$\begin{split}S(k)&=\frac{1}{N}\sum^{N}_{i,j}\langle\exp[-i\bm{k}\cdot(\bm{r}_{i}-\bm{r}_{j})]\rangle\\\
&=\frac{1}{N}\left\langle\left[\sum^{N}_{i}\sin(\bm{k}\cdot\bm{r}_{i})\right]^{2}\right\rangle+\frac{1}{N}\left\langle\left[\sum^{N}_{i}\cos(\bm{k}\cdot\bm{r}_{i})\right]^{2}\right\rangle.\end{split}$
(S49)

The smallest wavevector accessible using this method from our simulations is
$k\approx 0.08\,\text{\AA}^{-1}$. Additionally, we also determine the
isothermal compressibility,

$\kappa_{T}=-\frac{1}{V}\left(\frac{\partial V}{\partial P}\right)_{T},$ (S50)

where $V$ is the volume and $P$ is the pressure of the fluid, to calculate the
$k\rightarrow 0$ limit of the structure factor, given by the relation [3]

$S(0)=\rho_{\mathrm{u}}k_{\mathrm{B}}T\kappa_{T}.$ (S51)

We determine the isothermal compressibility by finite difference [30]

$\kappa_{T}=\frac{1}{\rho}\left(\frac{\partial\rho}{\partial
P}\right)_{T}\approx\frac{\ln(\rho_{1}/\rho_{0})}{P_{1}-P_{0}}.$ (S52)

To evaluate this expression, the system is simulated in the NVT ensemble at
densities $\rho_{0,1}=\rho_{\mathrm{u}}\pm 1\,\mathrm{nm^{-3}}$ and the
resulting pressures $P_{0,1}$ are sampled.

The full structure factor is obtained by combining data obtained from Eqs.
S48, S49 and S51. Data at low $k$ were interpolated using Scipy’s
UnivariateSpline routine [31] and appropriate Savitzky–Golay smoothing filter
[32] to be as fine as $\delta k=0.01\,\text{\AA}^{-1}$. The resulting
structure factors for the different water models are shown in the inset of
Fig. S4(a).

We then obtained the direct correlation function using the Ornstein–Zernike
relation [33, 3]

$1-\rho_{\mathrm{u}}\hat{c}^{(2)}(k)=\frac{1}{1+\rho_{\mathrm{u}}\hat{h}(k)}=\frac{1}{S(k)},$
(S53)

where $\hat{c}(k)$ and $\hat{h}(k)$ are the Fourier transforms of the direct
correlation function $c^{\mathrm{(2)}}(r)$ and the total correlation function
$h(r)=g(r)-1$. The resulting direct correlation functions in real space for
the different water models are shown in Fig. S4(b).

###  S3.4 Extracting solvation free energies of hard spheres

When the solutes of interest are small hard spheres with radius $R\lesssim
0.4\,\mathrm{nm}$, the solvation free energies can be determined from
simulations of the bulk liquid as

$\varOmega_{\mathrm{solv}}=-k_{\mathrm{B}}T\ln p(0;v),$ (S54)

where $p(0;v)$ is the probability of forming a spherical cavity of volume
$v=4\pi R^{3}/3$ containing $N_{v}=0$ solvent molecules [34]. For all water
models and state points considered, we computed $\varOmega_{\mathrm{solv}}$
for hard spheres of radii up to $\sim 0.4\,\mathrm{nm}$.

##  S4 Coarse-graining procedure

As derived in Section S1, the coarse-graining procedure emerges from the
theory as a consequence of splitting the density field into a slowly-varying
and rapidly varying component. While the definition of the coarse-graining
given in Eq. S31 involves the slowly-varying component of the direct
correlation function for the auxiliary system, $c_{\mathrm{r},1}^{(2)}$, for
all practical purposes, we use a Gaussian weight as in Eq. S39. This leaves
$a$ and $\lambda$ as parameters, for which we have used $a=200$ kJ cm3 mol-2
and $\lambda=0.08\,\text{nm}$ for SPC/E water. Here, we provide an analysis,
for SPC/E water at $300$ K, that explains how we arrived at these values.
Similar values are used for the other water models, which are given in Section
S6.

###  S4.1 Hard sphere splitting

For complex fluids, it is generally difficult to separate out the part of the
correlation function that is slowly-varying in space. To make progress, we
therefore assume that the short-ranged correlations in the uniform fluid can
be approximated by that of a uniform hard-sphere fluid at the same density,

$c^{(2)}_{\mathrm{u}}(r)\approx
c_{\mathrm{hs}}^{(2)}(r)+c_{\mathrm{u},1}^{(2)}(r),$ (S55)

where $c^{(2)}_{\mathrm{hs}}$ is the two-body direct correlation function of
the hard sphere fluid. The coarse-graining procedure becomes

$\overline{\rho}^{({\mathrm{hs}})}(\bm{r})=\frac{\\!\int\\!\mathrm{d}\bm{r}^{\prime}\,[c^{(2)}_{\mathrm{u}}(|\bm{r}-\bm{r}^{\prime}|)-c_{\mathrm{hs}}^{(2)}(|\bm{r}-\bm{r}^{\prime}|)]\,\rho(\bm{r}^{\prime})}{\\!\int\\!\mathrm{d}\bm{r}^{\prime}\,[c^{(2)}_{\mathrm{u}}(|\bm{r}-\bm{r}^{\prime}|)-c_{\mathrm{hs}}^{(2)}(|\bm{r}-\bm{r}^{\prime}|)]},$
(S56)

and

$2a^{({\mathrm{hs}})}=k_{\mathrm{B}}T\\!\int\\!\mathrm{d}\bm{r}\,[c^{(2)}_{\mathrm{u}}(r)-c_{\mathrm{hs}}^{(2)}(r)].$
(S57)

This is effectively treating water as a simple van der Waals fluid.

Here, we have used $c^{(2)}_{\mathrm{u}}$ obtained from simulation of bulk
SPC/E water and the Percus–Yevick expression [35] for $c_{\mathrm{hs}}^{(2)}$,
which is given by[3]

$c^{(2)}_{\mathrm{hs}}(r)=\begin{cases}-\sigma_{1}-\frac{6\eta\sigma_{2}}{d}r-\frac{\eta\sigma_{1}}{2d^{3}}r^{3}&\text{if
}r\leq d\\\ 0&\text{if }r>d,\\\ \end{cases}$ (S58)

where $d$ is the hard sphere diameter, $\eta=\frac{\pi d^{3}\rho}{6}$ is the
packing fraction with $\rho=\rho_{\mathrm{u}}=33.2\,\mathrm{nm^{-3}}$ and

$\sigma_{1}=\frac{(1+2\eta)^{2}}{(1-\eta)^{4}},\quad\sigma_{2}=-\frac{(2+\eta)^{2}}{4(1-\eta)^{4}}.$
(S59)

It is reasonable to choose $d$ to be at least the molecular diameter, and
acceptable values could range from $0.275$–$0.295\,\text{nm}$, corresponding
to the range of the first peak in $g(r)$. We then find
$a^{({\mathrm{hs}})}\approx 40$–$500$ kJ cm3 mol-2.

The hard sphere splitting of the direct correlation of SPC/E is shown in Fig.
S5(a) using $d=0.286$ nm. While such a splitting allows one to separate out
some of the short-ranged correlations, we see that
$c^{(2)}_{\mathrm{u}}(r)-c_{\mathrm{hs}}^{(2)}(r)$ is not completely slowly-
varying. This reflects the fact that water is not a simple fluid, and that the
hard-sphere fluid is not a suitable reference system to describe its behavior
quantitatively. We can circumvent the problem of finding an appropriate
splitting for the direct correlation function by carrying out the coarse
graining using a Gaussian weight, as given in Eq. S39. This is equivalent to
approximating the component of the direct correlation function responsible for
slowly-varying density $c^{(2)}_{\mathrm{u},1}$ with the form

$c^{(2)}_{\mathrm{u},1}(|\bm{r}-\bm{r}^{\prime}|)\approx\frac{2\beta
a}{(2\pi\lambda^{2})^{3/2}}\exp\left(-\frac{|\bm{r}-\bm{r}^{\prime}|^{2}}{2\lambda^{2}}\right).$
(S60)

With an appropriate choice for $a$ and $\lambda$, as will be rationalized in
the following subsection, this form of $c^{(2)}_{\mathrm{u},1}$ resembles the
slowly-varying behavior of $c^{(2)}_{\mathrm{u}}(r)-c_{\mathrm{hs}}^{(2)}(r)$,
as shown in the inset of Fig. S5(a).

We determined the free energies of solvating hard sphere solutes using this
hard sphere splitting. The results are shown in Fig. S5(b), together with the
HNCA and our theory using a Gaussian coarse graining. While the hard sphere
splitting is able to remedy the HNCA’s failure to capture the hydrophobic
crossover, $\varOmega_{\mathrm{solv}}$ for solutes with $R>0.4\,\mathrm{nm}$
is still overestimated in comparison to using a Gaussian coarse graining (and
results from molecular simulations, as described in the main article).

Figure S5: Hard sphere splitting of the direct correlation function (a) The
bulk direct correlation function of SPC/E water
($\rho_{\mathrm{u}}=33.2\,\mathrm{nm^{-3}}$, $300\,\mathrm{K}$) is split into
a hard-sphere component $c^{(2)}_{\mathrm{u},0}$ (Percus–Yevick
approximation[35, 3] with hard sphere diameter $d=0.286\,\text{nm}$, density
$\rho=\rho_{\mathrm{u}}$) and the remaining component
$c^{(2)}_{\mathrm{u},1}$. While $c^{(2)}_{\mathrm{u},1}$ by definition should
be a function that is slowly-varying over the molecular length scale, the hard
sphere splitting still leaves rapidly-varying oscillations in
$c^{(2)}_{\mathrm{u},1}$. Instead, by using Eq. S60 ($a=200$ kJ cm3 mol-2,
$\lambda=0.08$ nm), these rapidly-varying oscillations are smoothed out, as
shown in the inset of (a) in reciprocal space. In (b), we show
$\varOmega_{\mathrm{solv}}/4\pi R^{2}$ with the HNCA (dotted purple), our
theory with the hard-sphere splitting (dashed pink) and our theory using a
Gaussian weight (solid blue).

###  S4.2 Coarse-graining parameters

To determine the coarse-graining parameters, we first use
$c^{(2)}_{\mathrm{u}}(r)-c_{\mathrm{hs}}^{(2)}(r)$ from the hard sphere
splitting described above as a guide. For plausible values of $d$ based on the
position of water’s oxygen–oxygen radial distribution function ($0.275\lesssim
d/\text{nm}\lesssim 0.295$), we estimate a broad range for $a$
($a^{({\mathrm{hs}})}\approx 40$–$500$ kJ cm3 mol-2). As we need to perform a
simulation of the bulk fluid to parameterize $c^{(2)}_{\rm u}$, we can narrow
down this estimated range at negligible further computational cost by
comparing $\varOmega_{\mathrm{solv}}(a,\lambda)$ obtained from theory to the
result from molecular simulations for a solute on the verge of the crossover
regime. In Fig. S6(a) we show $\varOmega_{\mathrm{solv}}(a,\lambda)$ for
$R=0.32\,\mathrm{nm}$, identifying a narrow range of the coarse-graining
parameters $a\approx 100$–$300$ kJ cm3 mol-2 and $\lambda\approx
0.06$–$0.11\,\text{nm}$. Similar to $m(\rho^{*}_{\rm s})$ discussed in Section
S2, the range of acceptable $a$ and $\lambda$ gives an extra flexibility in
describing the crossover regime, as shown in Fig. S6(b). For SPC/E water, we
settled on $a=200$ kJ cm3 mol-2 and $\lambda=0.08$ nm.

Figure S6: Determining the coarse-graining parameters. The solvation free
energy of a hard sphere with $R=0.32\,\mathrm{nm}$ associated with different
sets of parameters $(a,\lambda)$ for the coarse-graining procedure is shown in
(a). The four sets of parameters selected all give $\varOmega_{\mathrm{solv}}$
in good agreement with the value obtained from simulation, which is indicated
by a marker on the colorbar. The parameters for each set are: ($a=200$ kJ cm3
mol-2, $\lambda=0.8\,\text{\AA}$) for A; ($a=300$ kJ cm3 mol-2,
$\lambda=0.8\,\text{\AA}$) for B; ($a=200$ kJ cm3 mol-2,
$\lambda=1.0\,\text{\AA}$) for C; and ($a=300$ kJ cm3 mol-2,
$\lambda=1.0\,\text{\AA}$) for D. In (b), set A is shown to give the best
agreement for $\varOmega_{\mathrm{solv}}/4\pi R^{2}$ both for $R\leq 0.4$ nm
(shown in the inset, simulation data are our own) and across the crossover
regime at $R\approx 1$ nm (simulation data from Ref. 18). The arrow indices
$\gamma=63.6\,\mathrm{mJ\,m^{-1}}$ for SPC/E water [10]. In all cases, we have
used the same parameterization for $m(\rho_{s}^{\ast})$ as described in
Section S2.

##  S5 Numerical procedure for LCW-style cDFT minimisation

###  S5.1 Self-consistent cycle

Eqs. S32 and S29 form a self-consistent pair of equations for the fluid
density that minimize the functional in Eqs. S18. For clarity, we rewrite
them:

$\displaystyle\omega^{\prime}(\rho_{\rm s}(\bm{r}))=m\nabla^{2}\rho_{\rm
s}(\bm{r})+2a[\overline{\rho}(\bm{r})-\overline{\rho}_{\rm s}(\bm{r})],$
(S61a) $\displaystyle\rho(\bm{r})=\rho_{\rm
s}(\bm{r})\exp\left[-\beta\phi(\bm{r})+\int\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\rm
s}(\bm{r},\bm{r}^{\prime})\delta_{\rm s}\rho(\bm{r}^{\prime})\right].$ (S61b)

The self-consistent cycle consists of the following steps:

  1. (i)

Make an initial guess for the slowly-varying density field, e.g.,
$\rho_{\mathrm{s}}(\bm{r})=\rho_{\mathrm{u}}$ and the fluid density, i.e.,
$\rho(\bm{r})=\rho_{\mathrm{u}}\exp\left[-\beta\phi(\bm{r})\right]$.

  2. (ii)

Fix $\rho_{\mathrm{s}}(\bm{r})$ and solve for a new fluid density
$\rho(\bm{r})$ from Eq. S61b using Picard iteration:

    1. (a)

Calculate the right hand side of Eq. S61b using the current
$\rho^{(i)}(\bm{r})$; call this $\tilde{\rho}^{(i)}(\bm{r})$.

    2. (b)

Calculate the new full density
$\rho^{(i+1)}(\bm{r})=(1-\alpha)\rho^{(i)}(\bm{r})+\alpha\tilde{\rho}^{(i)}(\bm{r})$.

    3. (c)

Set the new density $\rho^{(i+1)}(\bm{r})$ as the current density
$\rho^{(i)}(\bm{r})$. Go to step (a) until the solution for $\rho(\bm{r})$ is
converged.

  3. (iii)

Fix $\rho(\bm{r})$ and solve for a new slowly-varying density
$\rho_{\mathrm{s}}(\bm{r})$ from Eq. S61a using Picard iteration. We will
recast Eq. S61a in an iterative form in the next subsection.

    1. (a)

Calculate the right hand side of the iterative equation using the current
$\rho_{\mathrm{s}}^{(i)}(\bm{r})$; call this
$\tilde{\rho}_{\mathrm{s}}^{(i)}(\bm{r})$.

    2. (b)

Calculate the new slowly-varying density
$\rho_{\mathrm{s}}^{(i+1)}(\bm{r})=(1-\alpha_{\mathrm{s}})\rho^{(i)}(\bm{r})+\alpha_{\mathrm{s}}\tilde{\rho}_{\mathrm{s}}^{(i)}(\bm{r})$.

    3. (c)

Set the new density $\rho_{\mathrm{s}}^{(i+1)}(\bm{r})$ as the current density
$\rho_{\mathrm{s}}^{(i)}(\bm{r})$. Go to step (a) until the solution for
$\rho_{\mathrm{s}}(\bm{r})$ is converged.

  4. (iv)

Repeat step (ii) and (iii) until both solutions for
$\rho_{\mathrm{s}}(\bm{r})$ and $\rho(\bm{r})$ are converged.

The values of the mixing parameters $\alpha$ and $\alpha_{\mathrm{s}}$ used
depend on the system in question and typically $\alpha\in[0.01,0.05]$ and
$\alpha_{\mathrm{s}}\in[0.05,0.2]$. Solutions are considered converged when
difference between the current and new density profiles is less than the
tolerance ($\mathrm{rtol}\,=10^{-5}$, $\mathrm{atol}\,=10^{-8}$) using Numpy’s
allclose routine [36].

###  S5.2 Solving for the slowly-varying density field

To recast Eq. S61a in an iterative form, we can rewrite the term involving
$\nabla^{2}\rho_{\mathrm{s}}$ in terms of the coarse-grained density field
$\overline{\rho}_{\mathrm{s}}$[1]. By a second order Taylor expansion of
$\overline{\rho}(\bm{r})$ about $\overline{\rho}(\bm{r}^{\prime})$ in Eq. S39,
we can write

$\begin{split}\bar{\rho}(\bm{r})&\approx\int\\!\mathrm{d}\bm{r}^{\prime}\,\frac{1}{(2\pi\lambda^{2})^{3/2}}\exp\left(-\frac{|\bm{r}-\bm{r}^{\prime}|}{2\lambda^{2}}\right)\left[\rho(\bm{r})+\frac{1}{2}\nabla^{2}_{\bm{r}^{\prime}}\rho(\bm{r}^{\prime})|_{\bm{r}}(\bm{r}-\bm{r}^{\prime})^{2}+...\right]\\\
&=\rho(\bm{r})+\frac{1}{2}\lambda^{2}\nabla^{2}\rho(\bm{r}).\end{split}$ (S62)

Therefore, we can re-express the Laplacian of the slowly-varying density as

$\nabla^{2}\rho_{\mathrm{s}}=\frac{2}{\lambda^{2}}[\overline{\rho}_{\mathrm{s}}(\bm{r})-\rho_{\mathrm{s}}(\bm{r})].$
(S63)

By substitution of Eq. S63 into Eq. S61a, we reach an iterative form of Eq.
S61a

$\rho_{\mathrm{s}}(\bm{r})=\overline{\rho}_{\mathrm{s}}(\bm{r})-\frac{\lambda^{2}}{2m}\omega^{\prime}(\rho_{\mathrm{s}})+\frac{a\lambda^{2}}{m}\left[\overline{\rho}(\bm{r})-\overline{\rho}_{\mathrm{s}}(\bm{r})\right],$
(S64)

which is solved numerically by Picard iteration. We also note that for the
coarse-graining procedure, the Gaussian convolution is done in Fourier space,
i.e

$\hat{\bar{\rho}}(\bm{k})=\hat{\rho}(\bm{k})\exp\left(\frac{1}{2}|\bm{k}|^{2}\lambda^{2}\right).$
(S65)

###  S5.3 Solving for the rapidly-varying density field

By using the approximation in Eq. S30 for $c^{(2)}_{\mathrm{s}}$, Eq. S61b
becomes

$\rho(\bm{r})=\rho_{\mathrm{s}}(\bm{r})\exp\left[-\beta\phi(\bm{r})+\frac{\rho_{\mathrm{s}}(\bm{r})}{2\rho_{\mathrm{u}}^{2}}\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{u}}(|\bm{r}-\bm{r}^{\prime}|)\rho_{\mathrm{s}}(\bm{r}^{\prime})[\rho(\bm{r}^{\prime})-\rho_{\mathrm{s}}(\bm{r}^{\prime})]\right],$
(S66)

where the integral

$f_{1}(\bm{r})=\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{u}}(|\bm{r}-\bm{r}^{\prime}|)\rho_{\mathrm{s}}(\bm{r}^{\prime})[\rho(\bm{r}^{\prime})-\rho_{\mathrm{s}}(\bm{r}^{\prime})],$
(S67)

can be performed in Fourier space. By defining

$f_{2}(\bm{r})=\rho_{\mathrm{s}}(\bm{r})[\rho(\bm{r})-\rho_{\mathrm{s}}(\bm{r})],$
(S68)

the convolution theorem allows us to compute $\hat{f}_{1}$ as a product in
Fourier space

$\hat{f}_{1}(\bm{k})=\hat{c}^{(2)}_{\mathrm{u}}(|\bm{k}|)\hat{f}_{2}(\bm{k}).$
(S69)

###  S5.4 External potentials

We have considered spherical solutes with three different external potential
forms:

  1. (i)

Hard sphere solutes: For an ideal hydrophobic solute, the external potential
is

$\phi(r)=\begin{cases}\infty&\quad\quad\quad\text{for $r\leq R$}\\\
0&\quad\quad\quad\text{for $r>R$},\end{cases}$ (S70)

where $R$ is the hard sphere radius.

  2. (ii)

Lennard-Jones solutes: For solvation of a single Lennard-Jones particle with a
fixed well-depth of $\epsilon_{\mathrm{s}}=0.5\,\mathrm{kJ\,mol^{-1}}$ of
increasing size $\sigma_{\mathrm{s}}$ in SPC/E water (with Lennard–Jones
parameter as $\epsilon_{\mathrm{w}}=0.65\,\mathrm{kJ\,mol^{-1}}$,
$\sigma_{\mathrm{w}}=0.3166\,\mathrm{nm}$), as considered in Ref. 37, the
solute–solvent interaction is given by

$\phi(r)=4\epsilon_{\mathrm{sw}}\left[\left(\frac{R}{r}\right)^{12}-\left(\frac{R}{r}\right)^{6}\right],$
(S71)

where the effective strength is
$\epsilon_{\mathrm{sw}}=\sqrt{\epsilon_{\mathrm{s}}\epsilon_{\mathrm{w}}}$ and
the effective radius is
$R=\frac{1}{2}(\sigma_{\mathrm{s}}+\sigma_{\mathrm{w}})$.

  3. (iii)

For solvation of solutes with attractive tails following Coe _et al._ [38,
39], the external potential is

$\phi(r)=\begin{cases}\infty&\quad\quad\quad\text{for $r\leq R$}\\\
\epsilon_{\rm
sf}\bigg{[}\frac{2\sigma_{\mathrm{s}}^{9}}{15}\left(\frac{1}{r_{-}^{9}}-\frac{1}{r_{+}^{9}}\right)+\frac{3\sigma_{\mathrm{s}}^{9}}{20(r+r_{\rm
min})}\left(\frac{1}{r_{+}^{8}}-\frac{1}{r_{-}^{8}}\right)\\\
+\sigma_{\mathrm{s}}^{3}\left(\frac{1}{r_{+}^{3}}-\frac{1}{r_{-}^{3}}\right)+\frac{3\sigma_{\mathrm{s}}^{3}}{2(r+r_{\rm
min})}\left(\frac{1}{r_{-}^{2}}-\frac{1}{r_{+}^{2}}\right)\bigg{]}&\quad\quad\quad\text{for
$r>R$},\end{cases}$ (S72)

where $\epsilon_{\rm sf}$ is the effective solute–fluid attraction strength,
$R$ is the solute hard sphere radius, $r_{\rm min}$ is the location of the
attractive tail and $r_{\pm}=r+r_{\rm min}\pm R$. We considered solutes with
$R$ and $\epsilon_{\rm sf}$ as a multiple of the size and attractive strength
of a water “molecule” in the mW model
($\sigma_{\mathrm{mw}}=0.23925\,\mathrm{nm}$,
$\epsilon_{\mathrm{mw}}=25.895\,\mathrm{kJ\,mol^{-1}}$).

##  S6 Comparison to mDFT

###  S6.1 Theory

In mDFT [40, 41, 42, 43, 44], the excess intrinsic Helmholtz free energy
functional is written as an expansion about a bulk fluid of uniform density
$\rho_{\mathrm{u}}$ at chemical potential $\mu$

$\mathcal{F}^{(\rm ex)}_{\rm intr}[\rho]=\mathcal{F}^{(\rm ex)}_{\rm
u,intr}[\rho_{\rm u}]-k_{\rm B}T\int\\!\\!\mathrm{d}\bm{r}\,c_{\rm
u}^{(1)}(\bm{r})\delta_{\rm u}\rho(\bm{r})-\frac{k_{\rm
B}T}{2}\int\\!\\!\mathrm{d}\bm{r}\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,\delta_{\rm
u}\rho(\bm{r})c_{\mathrm{u}}^{(2)}(|\bm{r}-\bm{r}^{\prime}|)\delta_{\rm
u}\rho(\bm{r}^{\prime})+\mathcal{F}_{\rm bridge}[\rho],$ (S73)

where $\mathcal{F}^{(\rm ex)}_{\rm u,intr}$ is the intrinsic excess Helmholtz
free energy, $c_{\rm u}^{(1)}$ and $c_{\rm u}^{(2)}$ are the one-body and two-
body direct correlation functions of the uniform fluid and $\delta_{\rm
u}\rho=\rho-\rho_{\rm u}$. The “bridge” functional, $\mathcal{F}_{\rm
bridge}$, accounts for contributions to $\mathcal{F}^{(\rm ex)}_{\rm intr}$
beyond quadratic order in $\delta_{\rm u}\rho(\bm{r})$. Neglecting
$\mathcal{F}_{\rm bridge}$ amounts to the hypernetted-chain approximation
(HNCA) of integral equation theories [3]. Upon combining with the ideal
contribution in Eq. S3 and relating $c^{(1)}_{\mathrm{u}}$ to the excess part
of the intrinsic chemical potential,

$-k_{\mathrm{B}}Tc^{(1)}_{\mathrm{u}}(\bm{r})=\mu^{\mathrm{ex}}=\mu-
k_{\mathrm{B}}T\ln(\Lambda^{3}\rho_{\mathrm{u}}),$ (S74)

the intrinsic Helmholtz free energy functional becomes

$\begin{split}\mathcal{F}_{\mathrm{intr}}[\rho]=&\,k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\rho(\bm{r})[\ln(\Lambda^{3}\rho(\bm{r}))-1]+\mathcal{F}_{\mathrm{u,intr}}^{\mathrm{(ex)}}[\rho_{\mathrm{u}}]+\mu\\!\int\\!\\!\mathrm{d}\bm{r}\,\delta_{\mathrm{u}}\rho(\bm{r})-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\delta_{\mathrm{u}}\rho(\bm{r})\ln(\Lambda^{3}\rho_{\mathrm{u}})\\\
&-\frac{k_{\mathrm{B}}T}{2}\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{u}}(|\bm{r}-\bm{r}^{\prime}|)\,\delta_{\mathrm{u}}\rho(\bm{r})\,\delta_{\mathrm{u}}\rho(\bm{r}^{\prime})+\mathcal{F}_{\rm
bridge}[\rho]\\\
&+\mathcal{F}_{\mathrm{u,intr}}^{\mathrm{(id)}}[\rho_{\mathrm{u}}]-k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\rho_{\mathrm{u}}[\ln(\Lambda^{3}\rho_{\mathrm{u}})-1].\end{split}$
(S75)

Note that with the final two terms, we have simply added and subtracted the
intrinsic Helmholtz free energy for an ideal gas of uniform density
$\rho_{\mathrm{u}}$. Upon rearrangement, we arrive at

$\begin{split}\mathcal{F}_{\mathrm{intr}}[\rho]=&\mathcal{F}_{\mathrm{u,intr}}[\rho_{\mathrm{u}}]+(\mu-
k_{\mathrm{B}}T)\\!\int\\!\\!\mathrm{d}\bm{r}\,\delta_{\mathrm{u}}\rho(\bm{r})\\\
&+k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\rho(\bm{r})\ln\left(\frac{\rho(\bm{r})}{\rho_{\mathrm{u}}}\right)-\frac{k_{\mathrm{B}}T}{2}\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{u}}(|\bm{r}-\bm{r}^{\prime}|)\,\delta_{\mathrm{u}}\rho(\bm{r})\,\delta_{\mathrm{u}}\rho(\bm{r}^{\prime})+\mathcal{F}_{\rm
bridge}[\rho].\end{split}$ (S76)

The grand potential functional in mDFT is then given by

$\begin{split}\varOmega_{\phi}[\rho]=&\varOmega_{0}[\rho_{\mathrm{u}}]+\\!\int\\!\\!\mathrm{d}\bm{r}\,\phi(\bm{r})\rho(\bm{r})+k_{\mathrm{B}}T\\!\int\\!\\!\mathrm{d}\bm{r}\,\left[\rho(\bm{r})\ln\left(\frac{\rho(\bm{r})}{\rho_{\mathrm{u}}}\right)-\delta_{\mathrm{u}}\rho(\bm{r})\right]\\\
&-\frac{k_{\mathrm{B}}T}{2}\\!\int\\!\\!\mathrm{d}\bm{r}\\!\\!\int\\!\\!\mathrm{d}\bm{r}^{\prime}\,c^{(2)}_{\mathrm{u}}(|\bm{r}-\bm{r}^{\prime}|)\,\delta_{\mathrm{u}}\rho(\bm{r})\,\delta_{\mathrm{u}}\rho(\bm{r}^{\prime})+\mathcal{F}_{\rm
bridge}[\rho],\end{split}$ (S77)

where

$\varOmega_{0}[\rho_{\mathrm{u}}]=\mathcal{F}_{\mathrm{u,intr}}[\rho_{\mathrm{u}}]-\mu\\!\int\\!\\!\mathrm{d}\bm{r}\rho_{\mathrm{u}}.$
(S78)

Recently, Borgis et al. proposed a bridge functional based on the weighted
density approximation (WDA) [45, 46] parameterized with thermodynamic
properties of the bulk solvent [44, 47]. The functional was demonstrated to
accurately predict solvation free energies for a database of solutes with a
variety of shape and size. We shall refer to the theory as the HNC+WDA. In
this ansatz, the bridge functional takes a cubic form [47]

$\mathcal{F}_{\mathrm{bridge}}[\rho]=k_{\mathrm{B}}T\int\\!\mathrm{d}\bm{r}\,\frac{a_{\mathrm{b}}}{\rho_{\mathrm{u}}^{2}}\left(\overline{\rho}(\bm{r})-\rho_{\mathrm{u}}\right)^{3},$
(S79)

where $\bar{\rho}(\bm{r})$ is a weighted density given by a convolution of the
solvent density with a Gaussian function

$\overline{\rho}(\bm{r})=\int\\!\mathrm{d}\bm{r}^{\prime}\,\frac{1}{(2\pi\lambda_{\mathrm{b}}^{2})^{3/2}}\exp\left(-\frac{|\bm{r}-\bm{r}^{\prime}|^{2}}{2\lambda_{\mathrm{b}}^{2}}\right)\rho(\bm{r}^{\prime}),$
(S80)

where $\lambda_{\mathrm{b}}$ is the coarse-graining length.

The two parameters in this functional, $a_{\mathrm{b}}$ and
$\lambda_{\mathrm{b}}$, need to be determined. The parameter $a_{\mathrm{b}}$
is fixed by imposing the liquid–vapor coexistence condition, ensuring that

$\varOmega_{0}[\rho_{\mathrm{l}}]=\varOmega_{0}[\rho_{\mathrm{v}}].$ (S81)

If we assume that $\rho_{\mathrm{l}}\approx\rho_{\mathrm{u}}$ and
$\rho_{\mathrm{v}}\approx 0$, then one finds

$a_{\mathrm{b}}=1-\frac{1}{2}\rho_{\mathrm{u}}\hat{c}^{(2)}_{\mathrm{u}}(0),$
(S82)

where
$\hat{c}^{(2)}_{\mathrm{u}}(0)=\int\mathrm{d}\bm{r}c_{\mathrm{u}}^{(2)}(r)$.
Using the standard compressibility relation [3], $a_{\mathrm{b}}$ can be
directly related to the bulk compressibility $\chi_{\mathrm{u}}$ (or the
$k\rightarrow 0$ limit of the structure factor)

$a_{\mathrm{b}}=\frac{1}{2}\left(1+\frac{1}{k_{\mathrm{B}}T\rho_{\mathrm{u}}^{-1}\chi_{\mathrm{u}}}\right)=\frac{1}{2}\left(1+\frac{1}{S(0)}\right),$
(S83)

leaving $\lambda_{\mathrm{b}}$ as the only parameter to be determined.

One approach to determine $\lambda_{\rm b}$ is to require that
$\varOmega_{\mathrm{solv}}/4\pi R^{2}\sim\gamma$ for hard sphere solutes with
large $R$. This approach, however, does not guarantee the behaviour for small
solutes. Alternatively, one can choose $\lambda_{\mathrm{b}}$ such that
$\varOmega_{\mathrm{solv}}$ for small hard sphere solutes
($R<0.4\,\mathrm{nm}$) agrees well with estimates from molecular simulation
(which can be computed from bulk simulations, see Section S3). With this
second approach, however, it is not guaranteed that
$\varOmega_{\mathrm{solv}}/4\pi R^{2}\sim\gamma$ for large solutes.

###  S6.2 Numerical results

We now proceed to show our results using the HNC+WDA functional for the
different water models we have investigated, at the state points considered in
Section S2. Here, we have parameterized $\lambda_{\mathrm{b}}$ using
simulation data for $\varOmega_{\mathrm{solv}}$ of hard spheres with
$R<0.4\,\mathrm{nm}$ with the corresponding water models and state points. The
parameters used are given in Table S3.

Water model | SPC/E | RPBE-D3 | mW | mW
---|---|---|---|---
$T$ [K] | 300 | 300 | 426 | 300
$\rho_{\rm u}$ [$\mathrm{nm^{-3}}$] | 33.236 | 30.079 | 32.204 | 33.406
$a_{\mathrm{b}}$ | 8.5640 | 8.7297 | 11.3269 | 18.9831
$\lambda_{\mathrm{b}}$ [Å] | 1.0 | 1.0 | 0.95 | 1.0
Table S3: Parameters of the HNC+WDA functional from mDFT. For each water model
and temperature considered, the density $\rho_{\rm u}$ has been determined to
have a small chemical potential deviation from liquid–vapor coexistence from
square gradient theory. The parameter $a_{\mathrm{b}}$ is determined by the
bulk compressibility, as given in Eq. S83 by imposing the liquid–vapour
coexistence condition. The stated values of $\lambda_{\rm b}$ give
$\varOmega_{\rm solv}$ in best agreement with results from molecular
simulations for $R<0.4$ nm.

The results from the HNC+WDA functional for solvation of hard spheres are
shown in Fig. S7. We also show the results obtained from the cDFT approach
that we have presented. By construction, the HNC+WDA functional agrees very
well with simulation data for $R<0.4\,\mathrm{nm}$. As the solute size
increases, the HNC+WDA functional begins to disagree with our LCW-style
theory. In the case of SPC/E and RPBE-D3 water as seen in Figs. S7(a) and (b),
with $\lambda_{\mathrm{b}}\approx 1\,\text{\AA}$, while the HNC+WDA functional
can capture “hydrophobic crossover” similar to our approach, $\varOmega_{\rm
solv}/4\pi R^{2}$ is not guaranteed to converge to $\gamma$ at the
$R\rightarrow\infty$ limit. The results for the mW model of water, shown in
Figs. S7(c) and (d), are dramatic; HNC+WDA fails, even qualitatively, to
describe the expected behavior at large $R$. This highlights the difficulty in
capturing faithfully the behavior at each length scale when only a single
functional is constructed. In contrast, the theory from this work, which
solves for two functionals self-consistently, gives results that are
physically reasonable at both length scales.

Figure S7: Assessing the robustness of mDFT and our LCW-style cDFT to
different interatomic potentials for water. We show
$\varOmega_{\mathrm{solv}}/4\pi R^{2}$ for hard sphere solutes of increasing
radius $R$, predicted from the HNC+WDA approximation in the top panels (using
various values of the coarse-graining parameter $\lambda_{\mathrm{b}}$
indicated in the legend) and our LCW-style cDFT in the lower panels [dotted
lines indicate our theory parameterized with constant $m$, while solid lines
indicate our theory with $m(\rho_{\rm s}^{\ast})$]. Results are shown for (a)
the SPC/E water model at 300 K, (b) the RPBE-D3 electronic functional at 300
K, (c) the mW water model at 426 K, and (d) the mW water model at 300 K. In
each panel, we indicate $\gamma$ with an arrow, and simulation data with
circles.

The parameters we have used for the LCW-style cDFT to give the results, with a
constant $m$ and with $m(\rho^{*})$, are given in Table S4.

Water model | SPC/E | RPBE-D3 | mW | mW
---|---|---|---|---
$T$ [K] | 300 | 300 | 426 | 300
$a$ [kJ cm3 mol-2] | 200 | 300 | 300 | 300
$\lambda$ [nm] | 0.08 | 0.08 | 0.11 | 0.08
$m$ [$\mathrm{kJ\,mol^{-2}\,cm^{3}\,\text{\AA}^{2}}$] | 1255 | 1550 | 1100 | 1148
$t$ [$\mathrm{kJ\,mol^{-2}\,cm^{3}\,\text{\AA}^{5}}$] | -20323 | -21115 | -21115 | -18370
Table S4: Parameters of LCW-style cDFT for different water models. This table
gives the parameters used to give the results presented in Fig. S7.

##  S7 Measure of hydrophobicity and critical drying

###  S7.1 Local compressibility

A measure of the hydrophobicity, or more generally solvophobicity, of a solute
or substrate in the grand canonical ensemble is the local compressibility, as
proposed by Evans and Stewart [48]

$\chi(\bm{r})=\left(\frac{\partial\rho(\bm{r})}{\partial\mu}\right)_{T},$
(S84)

i.e., as the derivative of the local density $\rho(\bm{r})$ with respect to
the chemical potential $\mu$ at a fixed temperature $T$. In this work, we also
consider the contribution to $\chi(\bm{r})$ from the slowly-varying density
$\rho_{\mathrm{s}}(\bm{r})$, which is given as

$\chi_{\mathrm{s}}(\bm{r})=\left(\frac{\partial\rho_{\mathrm{s}}(\bm{r})}{\partial\mu}\right)_{T}.$
(S85)

In practice, we compute $\chi(\bm{r})$ through a finite difference

$\chi(\bm{r})\approx\frac{\rho(\bm{r};\mu+\Delta\mu)-\rho(\bm{r};\mu)}{\Delta\mu},$
(S86)

and similarly for $\chi_{\mathrm{s}}(\bm{r})$

$\chi_{\mathrm{s}}(\bm{r})\approx\frac{\rho_{\mathrm{s}}(\bm{r};\mu+\Delta\mu)-\rho_{\mathrm{s}}(\bm{r};\mu)}{\Delta\mu},$
(S87)

where $\Delta\mu$ is a small finite change in chemical potential. Since there
are available GCMC results from Refs. 38, 39, we will focus on mW water at
$426\,$K ($T/T_{\mathrm{c}}=0.46$).

###  S7.2 Water in contact with a planar hard wall

Figure S8: Critical drying for water at a planar hard wall from LCW-style
cDFT.  Density and local compressibility profiles for different chemical
potentials from cDFT. The fluid modelled mW at $426\,\mathrm{K}$ is in contact
with a planar hard wall. The expected values at contact from the respective
contact theorem for $\rho(0^{+})$ and $\chi(0^{+})$ are marked as open
circles.

To show the critical drying behavior prominently, we determine the density
profile $\rho(z)$ and local compressibility profile $\chi(z)$ of the fluid in
contact with a planar hard wall at various chemical potentials approaching the
coexistence point, as shown in Fig. S8. The equilibrium $\rho(z)$ at each
$\beta\delta\mu$ is determined carefully by starting with different initial
guesses with the vapor layer intruding up to various points away from the
wall, and the solution with the lowest grand potential is selected. We
determine $\chi(z)$ with a finite difference with $\Delta\mu$ at least two
times smaller than $\delta\mu$.

As $\beta\delta\mu\rightarrow 0$, an intruding vapor layer between the wall
and the liquid develops that becomes thicker with $\delta\mu$. This behaviour
is in agreement with what is expected as one moves away from the critical
drying point [38, 48, 49].

For a fluid in contact with a hard wall, the contact theorem [48] for the
density profile is

$\rho(0^{+})=\beta P,$ (S88)

where $P$ is the bulk pressure. In our square gradient functional, the
equivalent relation the contact density should follow is

$\rho(0^{+})=\beta\delta\mu\rho_{\mathrm{u}}-\beta\omega_{\mathrm{coex}}(\rho_{\mathrm{u}}),$
(S89)

consistent with the common tangent construction. For the local
compressibility, the contact theorem [48] is

$\frac{\chi(0^{+})}{\chi_{\mathrm{u}}}=\frac{1}{S(0)},$ (S90)

Our results for both $\rho(z)$ and $\chi(z)$ obey these relationship well for
$\beta\delta\mu\leq 10^{-1}$, establishing the limit for which the theory is
applicable away from coexistence.

We define the equilibrium vapor layer thickness $\ell_{\mathrm{eq}}$ as the
distance from the wall to where the density first approaches
$\rho(z)/\rho_{\mathrm{u}}=0.5$. From binding potential analysis as laid out
in Ref. 50, the divergence behaviors expected approaching the drying point are

$\ell_{\rm eq}\sim-\ln\delta\mu,$ (S91)

and

$\chi(\ell_{\rm eq})\sim\delta\mu^{-1}.$ (S92)

As in Ref. 50, we evaluate the local compressibility at $z=\ell_{\mathrm{eq}}$
using

$\chi(\ell_{\rm
eq})=\left(\frac{\partial\rho(z)}{\partial\mu}\right)_{T,\,z=\ell_{\rm
eq}}=-\rho^{\prime}(\ell_{\rm eq})\left(\frac{\partial\ell_{\rm
eq}}{\partial\delta\mu}\right)_{T},$ (S93)

where the prime denotes differentiation with respect to $z$. For a range of
$10^{-5}\leq\beta\delta\mu\leq 7\times 10^{-1}$, we verified that the
predictions for $\ell_{\mathrm{eq}}$ and $\chi(\ell_{\rm eq})$ from the theory
follow the scaling expected for critical drying behavior excellently, as shown
in the main paper.

###  S7.3 Water around hydrophobic spherical solutes

Having established that the LCW-style cDFT obey the scaling relations for
critical drying, we proceed to compute the local compressibility $\chi(r)$ of
the fluid around spherical hydrophobic solutes. We consider two cases. First
we investigate how $\chi$ changes as we increase the radius $R$ of a hard
sphere solute. Second, we consider how $\chi$ changes as we vary the
attractive interaction strength between a large solute and water (see Eq.
S72). The resulting compressibility and density profiles from our cDFT
approach are shown in Fig. S9. We have also shown the GCMC simulations results
from Coe _et al._ from Ref. 39. Overall, our results are in good agreement the
molecular simulations, indicating that the LCW-style cDFT approach contains
the essential physics of critical drying.

Figure S9: Measuring hydrophobicity across length scales with cDFT. The
compressibility $\chi(r)$ (top panels) and density profiles $\rho(r)$ (bottom
panels) are shown for hard sphere solutes of different radii $R$ [panel (a)
shows results from cDFT, panel (b) shows results from GCMC simulations]. We
also show results obtained with a solute with a hard repulsive core
($R=4.1\,\mathrm{nm}$) but with different solute–solvent attractive strengths,
indicated by $\epsilon_{\mathrm{sf}}/\epsilon_{\mathrm{mw}}$ [panel (c) shows
results from cDFT, panel (d) shows results from GCMC simulations]. The results
from cDFT are in overall good agreement with GCMC data from Ref. 39.

##  S8 Temperature dependence of hydrophobic solvation

###  S8.1 Parameterization

In this section, we assess how well our theory captures the temperature
dependence of hydrophobic solvation, paying particular attention to the
“entropic crossover” present in water but not in simple liquids. We will
consider solvation in RPBE-D3 and mW at state points neighbouring the triple
point up to the critical point. The fluid is in the liquid phase at
coexistence with density $\rho_{\mathrm{u}}=\rho_{\mathrm{l}}$. We show the
available data points from simulations for these water models in Fig. S10.

Figure S10: Temperature dependence of inputs of theory from simulation. For
the binodal curves in (a) and liquid–vapor surface tensions in (b), simulation
data for RPBE-D3 and mW are in taken from Refs. 12, 13, respectively. For the
bulk compressibilities in (c), we have obtained these values with our own
simulations of the bulk fluid at the liquid coexistence density
$\rho_{\mathrm{u}}=\rho_{\mathrm{l}}$.

For each water model at each temperature considered, we parameterized the
slowly-varying functional as described in Section S2 and summarized the
parameters for the theory in Table S5.

|  $T$ [K] | $T/T_{\mathrm{c}}$ |  $\rho_{\rm l}$ [nm-3] |  $\rho_{\rm v}$ [nm-3] |  $\gamma$ [$\mathrm{mJ\,m}^{-2}$] | $S(0)$ | $C$ | $D$ | $m$
---|---|---|---|---|---|---|---|---|---
RPBE-D3 | 300 | 0.47 | 30.0517 | 0.001003 | 68 | 0.0598 | 1.537 | 90 | 1500
350 | 0.55 | 29.2829 | 0.006686 | 57 | 0.0757 | 1.532 | 50 | 1750
400 | 0.63 | 27.9458 | 0.047133 | 46 | 0.1070 | 1.429 | 20 | 2100
450 | 0.71 | 26.3412 | 0.157111 | 36 | 0.1489 | 1.391 | -20 | 2600
500 | 0.79 | 24.3356 | 0.424535 | 28 | 0.2324 | 1.286 | -30 | 2800
550 | 0.87 | 21.9622 | 1.042952 | 15 | 0.3951 | 1.204 | -100 | 2400
mW | 300 | 0.33 | 33.4053 | 0.000013 | 65.8395 | 0.0265 | 2.523 | -15 | 1148
400 | 0.44 | 32.5038 | 0.000895 | 57.5166 | 0.0401 | 2.414 | -25 | 1215
500 | 0.54 | 31.1015 | 0.011361 | 47.7113 | 0.0712 | 1.943 | -30 | 1480
600 | 0.65 | 29.3987 | 0.064951 | 36.4794 | 0.1362 | 1.448 | -30 | 1680
700 | 0.76 | 27.1951 | 0.265254 | 23.4635 | 0.2956 | 0.998 | -35 | 1850
800 | 0.87 | 24.1901 | 0.951576 | 11.6765 | 0.8246 | 0.618 | -50 | 1950
Table S5: Parameters of the functional at different temperatures. The liquid
and vapor densities at coexistence, $\rho_{\mathrm{l}}$ and
$\rho_{\mathrm{v}}$, and surface tensions $\gamma$ have been determined with
simulations for RPBE-D3 in Ref. 12 and for mW in Ref. 13. The critical
temperature is at $T_{\mathrm{c}}=632\,$K for RPBE-D3 and
$T_{\mathrm{c}}=917.6\,$K for mW. We have determined the compressibilities of
the bulk liquid phase with our own simulations, given as the $k\rightarrow 0$
limit of the structure factors $S(k)$. The final three columns of the table
give the parameters for the slowly-varying functional in our theory: $C$ given
in $\mathrm{J\,mol^{-2}\,nm^{9}}$, $D$ given in
$\mathrm{J\,mol^{-2}\,nm^{15}}$ and $m$ given in
$\mathrm{kJ\,mol^{-2}\,cm^{3}\,\AA^{2}}$.

The direct correlation functions of the uniform liquid phases are obtained
from bulk simulations via the Ornstein–Zernike relation as described in
Section S3. We show the temperature dependence of these functions as well as
the square gradient theory of the free liquid–vapor interface for each water
model in Fig. S11.

Figure S11: Comparison of the bulk liquid local structure and liquid–vapor
interface of RPBE-D3 and mW. The direct correlation functions of the bulk
liquid at different temperatures are extracted from our own bulk simulations
and are shown in (a) and (d). The local grand potential densities
$\omega_{\mathrm{coex}}(\rho_{\mathrm{s}})$ parameterized based on bulk and
coexistence properties for each water model are shown in (b) and (e). The free
liquid–vapor interfaces for each water model at different temperatures are
shown in (c) and (d), with the solid lines being predictions from the theory
and open circles data from direct coexistence simulations. For RPBE-D3, the
density profiles are taken from Ref. 12. For mW, the density profiles are from
our own direct coexistence simulations.

###  S8.2 Entropic crossover

With the functions and parameters determined, we then determine the solvation
free energy of hard spheres of increasing radii $R$ at different temperatures
from the theory, as shown in Fig. S12. For both RPBE-D3 and mW water, the
theory predicts an “entropic crossover”, whereby $\varOmega_{\mathrm{solv}}$
increases with $T$ for small solutes with $R\lesssim 0.3\,\mathrm{nm}$ and
decreases with $T$ for larger solutes with $R\gtrsim 0.3\,\mathrm{nm}$, in
good agreement with data from simulations.

Figure S12: Capturing the entropic crossover with cDFT. The solvation free
energy at different temperatures are shown for hard sphere solutes of
different radii $R$ in RPBE-D3 and mW water. Panels (a) and (c) show the
results from cDFT while panels (b) and (c) show the results from simulations.
The results from cDFT are in good agreement with simulation, capturing the
entropic crossover in both water models.

## References

  * Lum, Chandler, and Weeks [1999] K. Lum, D. Chandler, and J. D. Weeks, “Hydrophobicity at small and large length scales,” The Journal of Physical Chemistry B 103, 4570–4577 (1999).
  * Weeks [2002] J. D. Weeks, “Connecting local structure to interface formation: A molecular scale van der Waals theory of nonuniform liquids,” Annual Review of Physical Chemistry 53, 533–562 (2002).
  * Hansen and McDonald [2013] J. Hansen and I. McDonald,  _Theory of Simple Liquids: with Applications to Soft Matter_ (Elsevier Science, 2013).
  * Weeks, Chandler, and Andersen [1971] J. D. Weeks, D. Chandler, and H. C. Andersen, “Role of repulsive forces in determining the equilibrium structure of simple liquids,” The Journal of Chemical Physics 54, 5237–5247 (1971).
  * Parry, Rascón, and Evans [2016] A. O. Parry, C. Rascón, and R. Evans, “The local structure factor near an interface; beyond extended capillary-wave models,” Journal of Physics: Condensed Matter 28, 244013 (2016).
  * Berendsen, Grigera, and Straatsma [1987] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, “The missing term in effective pair potentials,” The Journal of Physical Chemistry 91, 6269–6271 (1987).
  * Hammer, Hansen, and Nørskov [1999] B. Hammer, L. B. Hansen, and J. K. Nørskov, “Improved adsorption energetics within density-functional theory using revised Perdew-Burke-Ernzerhof functionals,” Phys. Rev. B 59, 7413–7421 (1999).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu,” The Journal of Chemical Physics 132, 154104 (2010).
  * Molinero and Moore [2009] V. Molinero and E. B. Moore, “Water modeled as an intermediate element between carbon and silicon,” The Journal of Physical Chemistry B 113, 4008–4016 (2009).
  * Vega and de Miguel [2007] C. Vega and E. de Miguel, “Surface tension of the most popular models of water by using the test-area simulation method,” The Journal of Chemical Physics 126, 154707 (2007).
  * Schienbein and Marx [2018] P. Schienbein and D. Marx, “Liquid–vapor phase diagram of RPBE-D3 water: Electronic properties along the coexistence curve and in the supercritical phase,” The Journal of Physical Chemistry B 122, 3318–3329 (2018).
  * Wohlfahrt, Dellago, and Sega [2020] O. Wohlfahrt, C. Dellago, and M. Sega, “Ab initio structure and thermodynamics of the RPBE-D3 water/vapor interface by neural-network molecular dynamics,” The Journal of Chemical Physics 153, 144710 (2020).
  * Coe, Evans, and Wilding [2022a] M. K. Coe, R. Evans, and N. B. Wilding, “The coexistence curve and surface tension of a monatomic water model,” The Journal of Chemical Physics 156, 154505 (2022a).
  * Huang and Chandler [2002] D. M. Huang and D. Chandler, “The hydrophobic effect and the influence of solute-solvent attractions,” The Journal of Physical Chemistry B 106, 2047–2053 (2002).
  * Morawietz  _et al._ [2016] T. Morawietz, A. Singraber, C. Dellago, and J. Behler, “How van der Waals interactions determine the unique properties of water,” Proceedings of the National Academy of Sciences 113, 8368–8373 (2016).
  * Schmidt  _et al._ [2009] J. Schmidt, J. VandeVondele, I.-F. W. Kuo, D. Sebastiani, J. I. Siepmann, J. Hutter, and C. J. Mundy, “Isobaric-isothermal molecular dynamics simulations utilizing density functional theory: An assessment of the structure and density of water at near-ambient conditions,” The Journal of Physical Chemistry B 113, 11959–11964 (2009).
  * Cerdeiriña  _et al._ [2011] C. A. Cerdeiriña, P. G. Debenedetti, P. J. Rossky, and N. Giovambattista, “Evaporation length scales of confined water and some common organic liquids,” The Journal of Physical Chemistry Letters 2, 1000–1003 (2011).
  * Huang, Geissler, and Chandler [2001] D. M. Huang, P. L. Geissler, and D. Chandler, “Scaling of hydrophobic solvation free energies,” The Journal of Physical Chemistry B 105, 6704–6709 (2001).
  * Huang and Chandler [2000a] D. M. Huang and D. Chandler, “Cavity formation and the drying transition in the Lennard-Jones fluid,” Phys. Rev. E 61, 1501–1506 (2000a).
  * Huang and Chandler [2000b] D. M. Huang and D. Chandler, “Temperature and length scale dependence of hydrophobic effects and their possible implications for protein folding,” Proceedings of the National Academy of Sciences 97, 8324–8327 (2000b).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Computer Physics Communications 271, 108171 (2022).
  * Bussi, Donadio, and Parrinello [2007] G. Bussi, D. Donadio, and M. Parrinello, “Canonical sampling through velocity rescaling,” The Journal of Chemical Physics 126, 014101 (2007).
  * Andersen [1983] H. C. Andersen, “Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations,” Journal of Computational Physics 52, 24–34 (1983).
  * Hockney and Eastwood [1988] R. Hockney and J. Eastwood,  _Computer Simulation Using Particles_ (Adam-Hilger, 1988).
  * Kolafa and Perram [1992] J. Kolafa and J. W. Perram, “Cutoff errors in the Ewald summation formulae for point charge systems,” Molecular Simulation 9, 351–368 (1992).
  * Stillinger and Weber [1985] F. H. Stillinger and T. A. Weber, “Computer simulation of local order in condensed phases of silicon,” Phys. Rev. B 31, 5262–5271 (1985).
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, “Generalized neural-network representation of high-dimensional potential-energy surfaces,” Phys. Rev. Lett. 98, 146401 (2007).
  * Singraber, Behler, and Dellago [2019] A. Singraber, J. Behler, and C. Dellago, “Library-based LAMMPS implementation of high-dimensional neural network potentials,” Journal of Chemical Theory and Computation 15, 1827–1840 (2019).
  * Sedlmeier, Horinek, and Netz [2011] F. Sedlmeier, D. Horinek, and R. R. Netz, “Spatial correlations of density and structural fluctuations in liquid water: A comparative simulation study,” Journal of the American Chemical Society 133, 1391–1398 (2011).
  * Motakabbir and Berkowitz [1990] K. A. Motakabbir and M. Berkowitz, “Isothermal compressibility of SPC/E water,” The Journal of Physical Chemistry 94, 8359–8362 (1990).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods 17, 261–272 (2020).
  * Savitzky and Golay [1964] A. Savitzky and M. J. E. Golay, “Smoothing and differentiation of data by simplified least squares procedures.” Analytical Chemistry 36, 1627–1639 (1964).
  * Ornstein and Zernike [1914] L. S. Ornstein and F. Zernike, “Accidental deviations of density and opalescence at the critical point of a single substance,” Royal Netherlands Academy of Arts and Sciences 17, 793–806 (1914).
  * Hummer  _et al._ [1996] G. Hummer, S. Garde, A. E. García, A. Pohorille, and L. R. Pratt, “An information theory model of hydrophobic interactions.” Proceedings of the National Academy of Sciences 93, 8951–8955 (1996).
  * Percus and Yevick [1958] J. K. Percus and G. J. Yevick, “Analysis of classical statistical mechanics by means of collective coordinates,” Phys. Rev. 110, 1–13 (1958).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” Nature 585, 357–362 (2020).
  * Fujita and Yamamoto [2017] T. Fujita and T. Yamamoto, “Assessing the accuracy of integral equation theories for nano-sized hydrophobic solutes in water,” The Journal of Chemical Physics 147, 014110 (2017).
  * Coe, Evans, and Wilding [2022b] M. K. Coe, R. Evans, and N. B. Wilding, “Density depletion and enhanced fluctuations in water near hydrophobic solutes: Identifying the underlying physics,” Phys. Rev. Lett. 128, 045501 (2022b).
  * Coe, Evans, and Wilding [2023] M. K. Coe, R. Evans, and N. B. Wilding, “Understanding the physics of hydrophobic solvation,” The Journal of Chemical Physics 158, 034508 (2023).
  * Zhao  _et al._ [2011] S. Zhao, R. Ramirez, R. Vuilleumier, and D. Borgis, “Molecular density functional theory of solvation: From polar solvents to water,” The Journal of Chemical Physics 134, 194102 (2011).
  * Jeanmairet  _et al._ [2013] G. Jeanmairet, M. Levesque, R. Vuilleumier, and D. Borgis, “Molecular density functional theory of water,” The Journal of Physical Chemistry Letters 4, 619–624 (2013).
  * Jeanmairet, Levesque, and Borgis [2013] G. Jeanmairet, M. Levesque, and D. Borgis, “Molecular density functional theory of water describing hydrophobicity at short and long length scales,” The Journal of Chemical Physics 139, 154101 (2013).
  * Jeanmairet  _et al._ [2016] G. Jeanmairet, N. Levy, M. Levesque, and D. Borgis, “Molecular density functional theory of water including density–polarization coupling,” Journal of Physics: Condensed Matter 28, 244005 (2016).
  * Borgis  _et al._ [2020] D. Borgis, S. Luukkonen, L. Belloni, and G. Jeanmairet, “Simple parameter-free bridge functionals for molecular density functional theory. Application to hydrophobic solvation,” The Journal of Physical Chemistry B 124, 6885–6893 (2020).
  * Tarazona [1985] P. Tarazona, “Free-energy density functional for hard spheres,” Phys. Rev. A 31, 2672–2679 (1985).
  * Curtin and Ashcroft [1985] W. A. Curtin and N. W. Ashcroft, “Weighted-density-functional theory of inhomogeneous liquids and the freezing transition,” Phys. Rev. A 32, 2909–2919 (1985).
  * Borgis  _et al._ [2021] D. Borgis, S. Luukkonen, L. Belloni, and G. Jeanmairet, “Accurate prediction of hydration free energies and solvation structures using molecular density functional theory with a simple bridge functional,” The Journal of Chemical Physics 155, 024117 (2021).
  * Evans and Stewart [2015] R. Evans and M. C. Stewart, “The local compressibility of liquids near non-adsorbing substrates: a useful measure of solvophobicity and hydrophobicity?” Journal of Physics: Condensed Matter 27, 194111 (2015).
  * Evans and Wilding [2015] R. Evans and N. B. Wilding, “Quantifying density fluctuations in water at a hydrophobic surface: Evidence for critical drying,” Phys. Rev. Lett. 115, 016103 (2015).
  * Evans, Stewart, and Wilding [2017] R. Evans, M. C. Stewart, and N. B. Wilding, “Drying and wetting transitions of a Lennard-Jones fluid: Simulations and density functional theory,” The Journal of Chemical Physics 147, 044701 (2017).