A classical density functional theory for solvation across length scales

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⁵ and nano-particles⁶ to macromolecules⁷ and extended surfaces.⁸ 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.⁹ 

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.¹⁴ 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¹⁵ 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.

Neglecting $Fbridge$ amounts to the hypernetted-chain approximation (HNCA) of integral equation theories.²⁸ Furthermore, upon linearizing $Δ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², 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.²⁰ 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.

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.

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.

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⁻³). 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⁴³), 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,³³ with Ω_solv increasing with volume for small solutes before crossing over to Ω _solv/4πR² ∼ γ 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 $ρs*$⁠. Specifically, $ρ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(ρ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.³⁰ 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⁻¹. Using the same parameterization for $m(ρs∗)$ as for hard spheres, our results for Ω_solv are in good agreement with available simulation data.⁴² 

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.

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⁵⁴ at T = 426 K.⁵⁵ We see good agreement between the results from our theory compared to grand canonical Monte Carlo (GCMC) simulations by Coe et al.⁵¹ 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 ρ(r) with ρ_s(r) in Eq. (29). As one might expect, as 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.⁵³ Specifically, for a fluid in contact with a planar hard wall, it can be shown that, close to coexistence, χ(ℓ_eq) ∼ δμ⁻¹, 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.

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⁵⁶ 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.⁵⁷ Here, we will demonstrate that our cDFT approach captures this entropic crossover. Moreover, we will do it from first principles.

Factors of k_BT 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⁵⁸ with dispersion corrections.⁵⁹ 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).

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,⁵¹ 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,⁵⁶ 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;¹⁹ 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.⁷¹ 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.

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.⁷² To maintain the temperature, we used the CSVR thermostat.⁷³ For simulations of SPC/E water,⁴³ long-ranged electrostatic interactions were evaluated using the particle–particle particle–mesh method,⁷⁴ such that the RMS error in the forces was a factor of 10⁵ smaller than the force between two unit charges separated by 0.1 nm.⁷⁵ The rigid geometry of the SPC/E water molecules was imposed with the RATTLE algorithm.⁷⁶ For simulations of the neural network surrogate model⁶⁰ of RPBE-D3,^58,59 we used the n2p2 package interface⁷⁷ 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.

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.

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.

The authors have no conflicts to disclose.

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).

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.