A central aim of multiscale modeling is to use results from the Schrödinger equation to predict phenomenology on length scales that far exceed those of typical molecular correlations. In this work, we present a new approach rooted in classical density functional theory (cDFT) that allows us to accurately describe the solvation of apolar solutes across length scales. Our approach builds on the Lum–Chandler–Weeks (LCW) theory of hydrophobicity [K. Lum *et al.*, J. Phys. Chem. B **103**, 4570 (1999)] by constructing a free energy functional that uses a slowly varying component of the density field as a reference. From a practical viewpoint, the theory we present is numerically simpler and generalizes to solutes with soft-core repulsion more easily than LCW theory. Furthermore, by assessing the local compressibility and its critical scaling behavior, we demonstrate that our LCW-style cDFT approach contains the physics of critical drying, which has been emphasized as an essential aspect of hydrophobicity by recent theories. As our approach is parameterized on the two-body direct correlation function of the uniform fluid and the liquid–vapor surface tension, it straightforwardly captures the temperature dependence of solvation. Moreover, we use our theory to describe solvation at a first-principles level on length scales that vastly exceed what is accessible to molecular simulations.

## I. INTRODUCTION

Many of the most fundamental processes in nature, including protein folding, crystallization, and self-assembly, occur in solution. Far from being an innocent bystander, the solvent often plays a vital role in determining the static and dynamic behaviors of these complex processes,^{1–4} owing to the delicate balance of solute–solute, solute–solvent, and solvent–solvent interactions. This provides a strong motivation to faithfully describe solvation behavior across a broad range of fields, from biological and chemical to physical and materials sciences. Solutes of interest can range in length scale from microscopic species^{5} and nano-particles^{6} to macromolecules^{7} and extended surfaces.^{8} Solvation is a multiscale problem.

The solvent, of course, comprises individual molecules. Molecular simulations, therefore, provide a natural approach to describe solvation and bestow fine details at time and length scales that can be challenging to access with experimental approaches alone. Depending upon the approximations made in describing the intermolecular interactions, molecular simulations provide one of the most accurate means to estimate solvation free energies. Yet the relatively high computational cost associated with molecular simulations makes their routine use for solutes much larger than small organic compounds cumbersome and inefficient.^{9}

Implicit solvation models alleviate the computational burden by describing the solvent degrees of freedom as a structureless continuum.^{10–13} Such approaches are numerically efficient, making it possible to routinely handle large macromolecules in solution. The major drawback of implicit solvation models, however, is their failure to account for essential solvent correlations that may hold a prominent role in the process under investigation (see, e.g., Ref. 2). They also often make assumptions concerning the validity of macroscopic laws applied to the microscopic domain, which can lead to inconsistent results.^{14} For these reasons, approaches that provide a coarse-grained description of solvation while retaining information about essential molecular correlations become very appealing. Such approaches should, for example, capture changes in the average equilibrium density field upon changing the solute–solvent interaction, as depicted in Fig. 1, without resorting to explicitly averaging microscopic degrees of freedom.

Motivated by both the seminal Lum–Chandler–Weeks (LCW) theory^{15} and more recent developments in molecular density functional theory (mDFT),^{16–21} in this article we present a classical density functional theory (cDFT)^{22–24} for the solvation of apolar solutes. Although our approach differs from the original LCW theory, it retains the essential feature of appropriately accounting for both slowly and rapidly varying components of the density field, as shown in Fig. 1(b). From a practical viewpoint, the theory we present lends itself more readily to numerical evaluation than LCW theory, including application to solutes with soft repulsive cores and attractive tails. Moreover, as we will discuss, our approach also offers conceptual advantages.

## II. RELEVANT BACKGROUND THEORY: cDFT FOR SOLVATION

*ρ*(

**) is the average of a microscopic density field of the fluid, whose chemical potential is**

*r**μ*and $Fintr$ is the intrinsic Helmholtz free energy functional independent of the external potential

*ϕ*

_{i}. The grand potential functional $\Omega \varphi i$ is minimized by the corresponding equilibrium density

*ρ*

_{i}, which satisfies

*β*= 1/(

*k*

_{B}

*T*),

*k*

_{B}is the Boltzmann constant,

*T*is the temperature, and $ci(1)$ is the one-body direct correlation function,

*i*for the quantities in Eq. (2).

^{25–27}Moreover, in cases where hard spheres act as a suitable reference fluid, attractive interactions can be reasonably treated in a mean-field fashion.

^{23}While approaches based on fundamental measure theory may capture some of the essential physics of more complex liquids such as water, it is unreasonable to expect quantitative agreement. Instead, in such cases, the mDFT approach

^{16–21}has shown great promise. The essential idea behind mDFT is that one can use the two-body direct correlation function,

*ρ*

_{u}to parameterize the grand potential functional

*δ*

_{u}

*ρ*(

**) =**

*r**ρ*(

**) −**

*r**ρ*

_{u}. The “bridge” functional, $Fbridge$, accounts for contributions to the excess part of $Fintr$ beyond quadratic order in

*δ*

_{u}

*ρ*(

**). The solvation free energy is simply**

*r*Neglecting $Fbridge$ amounts to the hypernetted-chain approximation (HNCA) of integral equation theories.^{28} Furthermore, upon linearizing $\Delta uFintr(id)$, mDFT within the HNCA is equivalent, up to a small correction factor, to Chandler’s Gaussian field theory for solvation.^{29,30} Aside from being numerically tractable, the HNCA provides reasonable accuracy for small solutes. However, as we have emphasized, solvation is a multiscale problem. To see this, consider that the solute simply excludes solvent density within a radius *R* of its center. Whereas within the HNCA, *Ω*_{solv} scales indefinitely with the solute’s volume,^{29–32} for large enough solutes we know from macroscopic theory that it scales with the solute’s surface area, i.e., *Ω*_{solv} ∼ 4*πγR*^{2}, where *γ* is the liquid–vapor surface tension. Many previous studies have shown that for water under ambient conditions, this “hydrophobic crossover” occurs at *R* ≈ 1 nm.^{15,33–36}

The failure of the HNCA for large solutes arises because it cannot describe water’s proximity to its liquid–vapor coexistence at ambient conditions.^{23,37} A natural progression, then, is to attempt to encode the physics of coexistence through $Fbridge$, while maintaining reference to the homogeneous fluid with density *ρ*_{u}. Such an approach has been used with some success within the mDFT framework for simple point charge water models.^{20,21} However, as we show in the supplementary material, existing bridge functionals of this kind are generally not robust to the choice of water model. We also note that attempts to use the “hard-sphere bridge functional” for water^{18,38,39} have proved problematic, failing to describe simultaneously the pressure and surface tension.^{20} In this article, we will instead follow more closely the key idea underlying LCW theory: the fluid’s density can be separated into a slowly varying component that can sustain interfaces and liquid–vapor coexistence and a rapidly varying component that describes the local structure on microscopic length scales.^{15,40}

In the following, we will outline how these ideas from LCW theory can be used to develop a cDFT approach to describe the solvation of apolar solutes across both small and large length scales. We will validate our theory against available simulation data for the solvation of both hard- and soft-core solutes. We will show that the physics of critical drying, which is essential for a faithful description of solvophobicity on large length scales, is well-described. We will also demonstrate that our approach captures a distinguishing feature of solvophobicity in complex liquids such as water—the “entropic crossover”—that is absent in simple liquids. We will use our theory to describe the hydrophobic effect at an *ab initio* level on length scales inaccessible to molecular simulations.

## III. A cDFT BUILT ON SEPARATION OF LENGTH SCALES

### A. Expansion about an inhomogenous density

*ρ*

_{s}(

**) that acts as a suitable reference. In line with Eq. (1),**

*r**ρ*

_{s}minimizes the grand potential functional $\Omega \varphi s$ prescribed by a slowly varying external potential

*ϕ*

_{s}(

**). Performing the expansion around**

*r**ρ*

_{s}gives

*δ*

_{s}

*ρ*(

**) =**

*r**ρ*(

**) −**

*r**ρ*

_{s}(

**), and the one- and two-body direct correlation functions are defined in an analogous manner to Eqs. (2) and (4). Upon addition of $Fintr(id)$ to Eq. (8) and substitution of $cs(1)$ according to its definition [Eq. (2)], we arrive at a slightly modified version of Eq. (5),**

*r**ϕ*

_{s}were known [and assuming that $cs(2)$ can be reasonably approximated], the reference density

*ρ*

_{s}would result from minimization of $\Omega \varphi s$ and, in turn,

*ρ*from minimization of Eq. (9). For example, with

*ϕ*

_{s}(

**) = 0, the HNCA is recovered. In the general context of solvation, however, the appropriate choice for**

*r**ϕ*

_{s}is unknown.

### B. Specifying an appropriate inhomogeneous reference using coexistence solutions

*ϕ*

_{s}is challenging. As we are interested in liquid water at thermodynamic states close to coexistence, we adopt a strategy that exploits the fact that, at coexistence, inhomogeneous density fields minimize the grand potential when subject to appropriate boundary conditions. For example, in a planar geometry, solutions corresponding to the free liquid-vapor interface are obtained by minimizing the grand potential at coexistence, subject to the conditions

*ρ*

_{v}and

*ρ*

_{l}are the vapor and liquid densities, respectively, at coexistence.

*ρ*

_{s}. However, the above-mentioned example demonstrates the associated challenges. Imagine a liquid in contact with a planar hard wall. Close to coexistence, the average solvent density profile will resemble that of a free liquid–vapor interface with a well-defined separation from the wall. In contrast, in the above-mentioned example, any density profile corresponding to free translation of the interface is equally plausible. We therefore seek a procedure that allows us to “pick” the appropriate coexistence solution. In fact, with the approximations that we will make, the coexistence solutions will ultimately only enter implicitly in our approach. To this end, we introduce a second expansion around an “auxiliary” reference density,

*ρ*

_{r}, which is an equilibrium solution at coexistence. Expanding the slowly varying reference density

*ρ*

_{s}around

*ρ*

_{r}, we write

*δ*

_{r}

*ρ*

_{s}(

**) =**

*r**ρ*

_{s}(

**) −**

*r**ρ*

_{r}(

**) and again, the one- and two-body direct correlation functions are defined in an analogous manner to Eqs. (2) and (4). Substituting Eq. (10) into Eq. (8) and keeping both expansions to second order, we obtain**

*r**ρ*(

**) and**

*r**ρ*

_{s}(

**) with differences between**

*r**ρ*

_{s}(

**) and**

*r**ρ*

_{r}(

**). As before, we will add $Fintr(id)$ to Eq. (11). We will also substitute $cr(1)$ according to Eq. (2), with**

*r**ϕ*

_{r}=

*δμ*=

*μ*−

*μ*

_{coex}, where

*μ*

_{coex}is the chemical potential at coexistence. The grand potential then reads

*ρ*

_{s}; provided that we have some procedure for specifying

*ρ*

_{s}to suit our needs, we do not need to know

*ϕ*

_{s}. Instead, we have transferred the problem to specifying the appropriate boundary conditions for

*ρ*

_{r}.

In the context of solvation, how to choose boundary conditions on *ρ*_{r} is not immediately obvious. Moreover, we appear to have complicated matters, as we now need to deal with three density fields (*ρ*, *ρ*_{s}, and *ρ*_{r}). We will now outline a series of approximations that enable us to consider a theory expressed explicitly in terms of *ρ* and *ρ*_{s} only, as well as a procedure to “pick” an appropriate slowly varying reference density.

*ρ*

_{s}that keeps it sufficiently close to

*ρ*while maintaining a slowly varying nature, we introduce the following “pseudofunctional,”

*ψ*is not an external potential; it depends parametrically on both

*ρ*(

**) and**

*r**ρ*

_{s}(

**), and penalizes differences between**

*r**ρ*and

*ρ*

_{s}, while ensuring that

*ρ*

_{s}remains slowly varying. (We will specify a form for

*ψ*below.) In this sense, it is not a functional in the usual cDFT sense, and we have introduced it simply as a computational tool to obtain an appropriate slowly varying reference.

*ρ*

_{s}is slowly varying, $Fintr[\rho s]$ can be reasonably approximated by the square gradient form

*ω*(

*ρ*

_{s}) is the local grand potential density. Assuming that

*m*is independent of

*ρ*

_{s}, minimizing $\Omega \u0303\psi $ with respect to

*ρ*

_{s}gives

*ω*′ indicates partial differentiation with respect to

*ρ*

_{s}.

*ρ*and

*ρ*

_{s}that do not feature

*ρ*

_{r}explicitly. To calculate the grand potential given by Eq. (12), however, still requires

*ρ*

_{r}, owing to the ideal terms and the final double integral. Considering the latter, as

*ρ*=

*ρ*

_{s}far from the solute, nonzero contributions to the integral are localized to regions close to the solute. Moreover, close to coexistence, we anticipate that

*δ*

_{r}

*ρ*

_{s}is small. In practice, then, we simply ignore the final term in Eq. (12) when computing the free energy. By similar reasoning, we also approximate $\Delta rFintr(id)[\rho ]\u2212\Delta rFintr(id)[\rho s]\u2248\Delta sFintr(id)[\rho ]$. When calculating the grand potential, we, therefore, adopt the simpler approximate form

*ψ*(

**; [**

*r**ρ*,

*ρ*

_{s}]) needs to be specified. In this work, we will not explore ways to systematically optimize

*ψ*. While, in principle, the fact that

*ψ*does not enter the grand potential explicitly offers some flexibility in specifying its form, it is not completely arbitrary; as mentioned previously, it ought to penalize differences between

*ρ*(

**) and**

*r**ρ*

_{s}(

**) while ensuring that**

*r**ρ*

_{s}(

**) remains slowly varying. (In addition, see below where we discuss our approach in the context of LCW theory.) With these considerations in mind, our choice for**

*r**ψ*is guided by the final term in Eq. (12). Specifically, we assume that $cr(2)$ can be separated as

**−**

*r***′| >**

*r**ℓ*, where

*ℓ*is a molecular length scale, accounting for rapid variations. Conversely, $cr,1(1)$ is slowly varying over length scales comparable to

*ℓ*. In keeping with the requirements for

*ψ*, we, therefore, prescribe

### C. Connecting to LCW theory for a practical cDFT scheme

*ψ*given by Eq. (21), together with Eqs. (15) and (18), the approximate cDFT that we have derived bears a striking resemblance to LCW theory. Indeed,

*ψ*plays the role of the “unbalancing potential” in Ref. 15. The formal similarity can be made more apparent if we introduce the coarse-graining procedure,

^{40}in which $cr,1(2)$ is treated in a random phase approximation. [Although also note that Eqs. (15), (18), and (19) do not constitute a random phase approximation.] In this case, there is a notion of “unbalanced attractive interactions,” and in the original LCW theory,

^{15}it was argued that such unbalanced forces in water could be described by a similar coarse graining over an appropriate length scale that describes the range of attractive interactions.

*ψ*in terms of the coarse-grained density fields as in Eq. (23), therefore, also serves a practical purpose, as we can estimate its value with an approximate coarse-graining, e.g.,

*λ*, we now treat

*a*as a parameter that needs to be determined. As detailed in the supplementary material, we estimate their values by using a hard sphere reference fluid for which $cr,0(2)$ is known. For liquid water at 300 K and

*βδμ*≈ 10

^{−3}(corresponding to

*ρ*

_{u}≈ 33.234 nm

^{−3}for SPC/E water), we obtain an acceptable range of the coarse-graining parameters

*a*≈ 200–300 kJ cm

^{3}mol

^{−2}and

*λ*≈ 0.08–0.11 nm. Values in this range for

*λ*are much smaller than that used in the original LCW theory (

*λ*≈ 0.38 nm), although they are comparable to the coarse-graining length derived in Ref. 36 for a lattice-based version of LCW theory that emphasizes the importance of capillary wave fluctuations.

*ω*

_{coex}takes a quartic form

*h*(

*ρ*

_{s}) is the Heaviside step function. The curvature at the minima is determined by

*C*, whose value we set to be consistent with the compressibility of the bulk fluid

*χ*

_{u}. Together,

*D*and

*m*determine

*γ*and the shape of the free liquid–vapor interface. Details of our parameterization procedure are provided in the supplementary material. In most of what follows, we will treat

*m*as a constant independent of density; for spherical solutes, this yields the correct limiting behavior for both

*R*→ 0 and

*R*→ ∞. For

*R*in the crossover regime, we will show that the extra flexibility afforded by allowing

*m*to vary in a systematic, yet practical, way with

*ρ*

_{s}yields quantitative agreement for intermediate solute sizes.

*ρ*

_{s}as we have here.

In this section, we have derived a cDFT for solvation that aims to describe both the short-wavelength perturbations induced by effects of excluded volume and the physics of liquid–vapor coexistence relevant for larger solutes. We have also described approximations for the coarse-grained density and inhomogeneous two-body direct correlation function that facilitate numerical evaluation. While our theory is appropriate for arbitrarily complex forms of apolar solute–solvent interaction, for demonstrative purposes, we will focus on spherical solutes and planar walls. These simple geometries, however, are sufficient to make connections with existing theories and highlight the advantages of our approach.

## IV. SOLVATION OF SPHERICAL SOLUTES, BOTH HARD AND SOFT

We first consider the paradigmatic test case for any theory of hydrophobicity: the solvation of a hard-sphere solute in water under ambient conditions (*T* = 300 K, *βδμ* ≈ 10^{−3}). To this end, in Figs. 2(a) and 2(b), we present *ρ*(*r*) and *Ω*_{solv} obtained from self-consistent solutions of Eqs. (15) and (23), parameterized for a simple point charge model of water (SPC/E^{43}), as we increase the solute radius *R*. First treating *m* as a constant (dotted line in Fig. 2), we see that our cDFT approach broadly captures the behavior observed in molecular simulations,^{33} with *Ω*_{solv} increasing with volume for small solutes before crossing over to *Ω* _{solv}/4*πR*^{2} ∼ *γ* for *R* ≳ 1 nm. Like previous treatments rooted in LCW theory, however, we also see that these results overestimate *Ω*_{solv} in the crossover regime.

Empirically, agreement with simulation results in the crossover regime could be attained by reducing the value of *m*. In the context of our theory, this implies that the optimal reference density for finite *R* has a sharper interface than the free liquid-vapor interface. Simply reducing *m* in such a manner, however, would result in an incorrect limiting behavior as *R* → ∞. However, this observation motivates us to vary *m* in a linear fashion with a characteristic feature of the slowly varying density $\rho s*$. Specifically, $\rho s*$ is the value where *ρ*_{s} = *ρ* for the first time as *r* increases. Further details are provided in the supplementary material. As we see in Fig. 2(b) (solid line), with $m(\rho s*)$, our cDFT describes the simulation data quantitatively.

For hard spheres in water, the main advantage of our approach compared to LCW theory is primarily conceptual: *Ω*_{solv} follows directly from the minimization of *Ω*_{ϕ}, as opposed to relying on, e.g., thermodynamic integration.^{34,44,45} From a practical viewpoint, the cDFT approach is also numerically simpler to implement, as discussed in detail in Ref. 29. However, its main advantages compared to LCW theory become apparent once we depart from the solvation of ideal hydrophobes. In LCW theory, it is assumed that the solute–solvent interaction can be separated into repulsive and attractive contributions, *ϕ* = *ϕ*_{rep} + *ϕ*_{att}. While *ϕ*_{att} can be accounted for straightforwardly in the slowly varying part, *ϕ*_{rep} must be approximated by a hard-core potential, and the solvent density is solved subject to the constraint *ρ*(*r*) = 0 inside the solute.^{30} For cDFT approaches such as ours, such an approximation is not needed. Instead, we simply minimize *Ω*_{ϕ} with the appropriate solute–solvent interaction *ϕ*. To illustrate this, in Figs. 2(c) and 2(d), we show results from our theory for the solvation of Lennard-Jones solutes of effective radius *R* (see supplementary material) with a constant well-depth of *ϵ* = 0.5 kJ mol^{−1}. Using the same parameterization for $m(\rho s\u2217)$ as for hard spheres, our results for *Ω*_{solv} are in good agreement with available simulation data.^{42}

## V. THE PHYSICS OF CRITICAL DRYING IN AN LCW-STYLE THEORY

The motivations for our work primarily stem from the desire to develop “semi-implicit” solvation models that retain information on essential solvent correlations. To that end, the results we have presented so far are promising. However, the physics of hydrophobicity is important in its own right. While LCW theory should be considered a seminal contribution in this area, subsequent work from Evans, Wilding, and co-workers^{46–51} emphasizes the central role of the critical drying transition that occurs in the limit of a planar substrate with vanishing, or very weak, attractive interactions with the solvent. The extent to which critical drying is captured by LCW theory (or its subsequent lattice-based derivatives^{2,36,52}) is, however, unclear. Although not equivalent to LCW theory, we can use our cDFT approach to shed light on the extent to which it contains the physics of critical drying.

*χ*(

**) provides one of the most robust indicators of hydrophobicity or, more generally, “solvophobicity.”**

*r*In the specific case of solvation of hard spheres, *χ*(*r*) should exhibit a pronounced maximum that increases in both magnitude and position with increasing *R*. In Fig. 3(a), we present *χ*(*r*) from our theory, parameterized for the coarse-grained mW water model^{54} at *T* = 426 K.^{55} We see good agreement between the results from our theory compared to grand canonical Monte Carlo (GCMC) simulations by Coe *et al.*^{51} In particular, for the largest solute investigated, *R* ≈ 4.1 nm, we see that the maximum value of *χ*(*r*) is over 40 times larger than its bulk counterpart *χ*_{u}. In the supplementary material, we also present results where the strength of attractive solute–solvent interactions is decreased, which also compare favorably to GCMC simulations. In the context of an LCW-style theory such as ours, we can also isolate the slowly varying component of the local compressibility *χ*_{s}(** r**), by replacing

*ρ*(

**) with**

*r**ρ*

_{s}(

**) in Eq. (29). As one might expect, as**

*r**R*increases, so does the significance of

*χ*

_{s}. It is clear, however, that contributions from the rapidly varying density are still important, even for solute sizes that far exceed

*R*= 1 nm, the colloquial hydrophobic crossover point, as seen in the inset.

A more stringent test of critical drying comes not from comparison to molecular simulations but from known scaling behaviors of *χ* from binding potential analyses.^{53} Specifically, for a fluid in contact with a planar hard wall, it can be shown that, close to coexistence, *χ*(*ℓ*_{eq}) ∼ *δμ*^{−1}, where *ℓ*_{eq} is the position of the maximum in *χ*. Moreover, *ℓ*_{eq} ∼ − ln *βδμ*. As seen in Fig. 3(b), the LCW-style cDFT approach that we have derived obeys both of these scaling relations. This is far from a trivial result. At face value, the physics of critical drying and the emphasis placed on liquid–vapor interface formation by LCW theory seem unrelated. By recasting the essential underpinnings of LCW theory in the context of cDFT, we begin to paint a unifying picture of these two different views on hydrophobicity.

## VI. MULTISCALE SOLVATION FROM FIRST PRINCIPLES

The physics of critical drying discussed so far is an essential component of hydrophobicity at large length scales. Moreover, it is common to both simple and complex fluids that exhibit solvophobicity. A distinguishing feature of solvophobicity in complex fluids such as water, however, is the “entropic crossover.” For small solutes, *Ω*_{solv} increases with increasing *T*, while for larger solutes it decreases. This behavior has implications for the thermodynamics of protein folding^{56} and has been attributed to a competition of microscopic length scales (i.e., solvent reorganization in the first and second solvation shells) that is absent in simple fluids.^{57} Here, we will demonstrate that our cDFT approach captures this entropic crossover. Moreover, we will do it from first principles.

Factors of *k*_{B}*T* aside, temperature dependence in Eqs. (15) and (23) enters implicitly through $cu(2)$ and *γ*. As the approach that we have developed does not assume a simple pairwise additive form for the interatomic potential, the combination of our theory with recent advances in machine-learned interatomic potentials (MLIPs) means we face the exciting prospect of describing, from first principles, the temperature dependence of solvation across the micro-, meso-, and macroscales.

For illustrative purposes, we model water with RPBE-D3, a generalized gradient approximation (electronic) functional^{58} with dispersion corrections.^{59} To obtain $cu(2)$, we have performed our own simulations of bulk water at the liquid density along the coexistence curve, using the MLIP described in Ref. 60, which also provides *γ*(*T*). In Fig. 4(a), we present results from self-consistent solutions of Eqs. (15) and (23) for temperatures ranging from 300 to 550 K. Note that, in the absence of information in the crossover regime, we have simply taken *m* to be independent of density. For all temperatures, we observe the hydrophobic crossover. Importantly, this occurs at progressively smaller values of *R* as *T* increases; we observe the entropic crossover. These results demonstrate that the LCW-style cDFT that we have developed can be used to faithfully coarse-grain solvent effects while maintaining essential molecular correlations to model phenomenology at a first principles level across a broad range of length scales. This approach is summarized schematically in Fig. 4(b).

## VII. CONCLUSIONS AND OUTLOOK

We have presented a cDFT for the solvation of apolar solutes in water, which is accurate at both small and large length scales. Similar to mDFT, we encode information about the liquid’s small length scale fluctuations by parameterizing our functional on the two-body direct correlation function obtained from simulations of the bulk fluid. In contrast to previous approaches, however, the grand potential that we construct is not based on an expansion around the uniform bulk fluid. Rather, from the outset, our theory acknowledges that the perturbations induced in the solvent density field by large solutes are too severe for the uniform fluid to act as a suitable reference. In this work, we therefore establish a self-consistent cDFT framework that permits the use of an inhomogeneous and slowly varying density field as a reference system.

The theory that we have outlined is similar in spirit to, and indeed motivated by, the seminal work on the hydrophobic effect of Ref. 15. By placing the ideas of LCW theory in the context of cDFT, not only do we gain a numerical advantage, but we also provide conceptual insights. For example, the grand potential that we derive [Eq. (12)] suggests a natural form for the “unbalancing potential” that specifies the coarse-graining; the coarse-graining function is the slowly varying part of a two-body direct correlation function of an inhomogeneous density field. This insight, as we explore in the supplementary material, justifies a coarse-graining length scale much smaller than the molecular diameter of a water molecule that one might naively expect.

Our approach also allows us to connect the ideas of LCW theory directly with more recent theoretical descriptions of hydrophobicity from Evans, Wilding, and co-workers.^{48,49} Specifically, the local compressibility in the presence of large hydrophobes obtained with our approach compares favorably to results obtained by GCMC simulations,^{51} and we show that its variation with chemical potential obeys known critical scaling behaviors. We also demonstrate that our approach, similar to previous LCW treatments,^{56} captures the temperature dependence of the hydrophobic effect. It is a curious observation that lattice-based theories of hydrophobicity^{36,52,62} that aim to improve LCW have emphasized the importance of capillary wave fluctuations; these do not enter explicitly in our theory. Nonetheless, our results for the local compressibility and the good agreement with solvation free energies obtained from molecular simulations suggest that we capture the most salient aspects of the interfacial fluctuations necessary to describe hydrophobicity.

A general theory of solvation should also describe the polarization field induced by charged species such as ions. This is a challenging problem beyond the scope of the present study. There are, however, reasons to be optimistic. For example, the mDFT framework already demonstrates that orientational correlations of the bulk fluid can be used to construct density functionals that describe polarization;^{19} introducing the results from our work should amount to a straightforward modification of the present mDFT framework. Insights from Weeks’ local molecular field theory may also prove useful in developing new functionals^{45,63–67} (see Refs. 45 and 68 for discussions on the relationship between cDFT and local molecular field theory), as could ideas from integral equation theories (see, e.g., Refs. 69 and 70). In a recent study, Sammüller *et al.* adopted an entirely different approach by using deep neural networks to construct the free energy functional.^{71} Whether such an approach can be used for polar fluids remains to be seen, but it shows great promise. In the context of apolar solvation, our results raise the question of whether machine learning can be made easier by first separating the functional into slowly and rapidly varying contributions. Notwithstanding the obvious areas for future development, the results of our work demonstrate a significant step toward efficient, and fully first-principles, multiscale modeling of solvation.

## VIII. METHODS

Here, we provide a brief overview of the methods used; full details are provided in the supplementary material. All molecular simulations have been performed with the LAMMPS simulation package.^{72} To maintain the temperature, we used the CSVR thermostat.^{73} For simulations of SPC/E water,^{43} long-ranged electrostatic interactions were evaluated using the particle–particle particle–mesh method,^{74} such that the RMS error in the forces was a factor of 10^{5} smaller than the force between two unit charges separated by 0.1 nm.^{75} The rigid geometry of the SPC/E water molecules was imposed with the RATTLE algorithm.^{76} For simulations of the neural network surrogate model^{60} of RPBE-D3,^{58,59} we used the n2p2 package interface^{77} with LAMMPS. cDFT calculations were performed with our own bespoke code, which we have made publicly available. Minimization was performed self-consistently using Picard iteration.

## SUPPLEMENTARY MATERIAL

The supplementary material includes further details on parameterizing the functional, simulation details, further details on the numerical implementation of the functional, a comparison to molecular density functional theory, and details on evaluating the local compressibility.

## ACKNOWLEDGMENTS

We acknowledge Robert Jack for many insightful discussions and Bob Evans and Nigel Wilding for comments on an initial draft of the paper. We are grateful for computational support from the UK national high performance computing service, ARCHER2, for which access was obtained via the UKCP consortium and funded by EPSRC Grant No. EP/X035891/1. A.T.B. acknowledges funding from the Oppenheimer Fund and Peterhouse College, University of Cambridge. S.J.C. is a Royal Society University Research Fellow (Grant No. URF\R1\211144) at the University of Cambridge.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Anna T. Bui**: Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Stephen J. Cox**: Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

Code for minimizing the functionals can be accessed at https://github.com/annatbui/LCW-cDFT. Simulation input files and the direct correlation functions obtained are openly available at the University of Cambridge Data Repository, https://doi.org/10.17863/CAM.104278.

## REFERENCES

*Fundamentals of Inhomogeneous Fluids*

*Theory of Simple Liquids: With Applications to Soft Matter*

*T* = 426 K corresponds to *T*/*Tc* ≈ 0.46, with the critical point of mW at *Tc* = 917.6 K (from Ref. 50).

*ab initio*parametrization of density functional dispersion correction (DFT-D) for the 94 elements H–Pu

*Ab initio*structure and thermodynamics of the RPBE-D3 water/vapor interface by neural-network molecular dynamics

*Computer Simulation Using Particles*