Solvation in atomic liquids: connection between Gaussian field theory and density functional theory

For the problem of molecular solvation, formulated as a liquid submitted to the external potential field created by a molecular solute of arbitrary shape dissolved in that solvent, we draw a connection between the Gaussian field theory derived by David Chandler [Phys. Rev. E, 1993, 48, 2898] and classical density functional theory. We show that Chandler's results concerning the solvation of a hard core of arbitrary shape can be recovered by either minimising a linearised HNC functional using an auxiliary Lagrange multiplier field to impose a vanishing density inside the core, or by minimising this functional directly outside the core --- indeed a simpler procedure. Those equivalent approaches are compared to two other variants of DFT, either in the HNC, or partially linearised HNC approximation, for the solvation of a Lennard-Jones solute of increasing size in a Lennard-Jones solvent. Compared to Monte-Carlo simulations, all those theories give acceptable results for the inhomogeneous solvent structure, but are completely out-of-range for the solvation free-energies. This can be fixed in DFT by adding a hard-sphere bridge correction to the HNC functional.


Introduction
In a world of hard-core numerical simulations on huge computers where most problems in solution chemistry are formulated in terms of molecular dynamics simulations and subsequent data analysis, it is wise to keep simpler methods that make it possible to derive analytical results or to perform the calculations with reasonable computer resources. Such methods rely on the statistical mechanics of atomic and molecular liquids that has been developed in the second half of the last century and are found by now in classical textbooks [1][2][3]. Along this vein is the beautiful and appealing recent theoretical work of Dung Di Caprio and Jean-Pierre Badiali who were able to formulate the description of classical fluids at equilibrium as a formally exact field theory [4][5][6][7][8]; this formalism was applied to model atomic and molecular fluids at solid interfaces [9][10][11][12]. Other more traditional approaches include molecular integral equation theories in the reference interaction site (RISM) [13][14][15][16], molecular [17][18][19][20][21][22][23][24], or mixed [25,26] picture and the density functional theory (DFT) in its atomic [27][28][29] or molecular version [30][31][32][33][34][35][36][37][38].
The basic theoretical principles of classical DFT can be found in the seminal paper by Evans [27] and subsequent excellent reviews by him [27][28][29] and other authors [39]. The advent in the late 1980's of a quasi-exact DFT for inhomogeneous hard sphere mixtures, the fundamental measure theory (FMT) [40][41][42][43][44][45], has recently promoted a great deal of applications to atomic-like fluids in bulk or confined conditions or at interfaces. Classical "atomic" DFT can be nowadays considered as a method of choice for many chemical engineering problems [46,47]. Much less applications exist for molecular fluids for which solvent orientations should be considered. The description has been generally limited to generic dipolar solvents or dipolar solvent/ions mixtures [32][33][34][35]; such an approach may be already considered as "civilized" compared to primitive continuum models [34]. We have proposed an extension of molecular DFT to arbitrary fluid/solvents (the so-called MDFT method) with the goal of describing the solvation of three-dimensional molecular object in those solvents [36][37][38][48][49][50][51][52][53][54][55][56]. Note that a 3D-version of the RISM equations [57][58][59][60][61][62], as well as a RISM-based DFT approach [63,64] have also been recently developed with the same goal.
In this paper, we also elaborate on a field theoretical approach that is different from the one by di Caprio and Badiali -and certainly starts from a less fundamental ground. We refer to the Gaussian field theory (GFT) of fluids developed by Chandler and collaborators [65][66][67][68]. Our main focus will be to draw a connection between the GFT approach of Chandler and our favorite classical DFT in the context of molecular solvation, i.e., a liquid submitted to an external potential field v(r) created by a molecular solute of arbitrary shape dissolved at infinite dilution in it. For simplicity, we restrict the discussion to atomic or pseudo-atomic solvents (such as CCl 4 ) modelled by spherical Lennard-Jones particles for which only the position r matters.

Density functional theory and HNC approximation
We begin by recalling the basis of the density functional theory of liquids submitted to an external potential field v(r). The grand potential density functional for a fluid having an inhomogeneous density ρ(r) in the presence of an external field v(r) can be defined as [27,28] is the Helmholtz free energy functional and µ s is the chemical potential. The grand potential can be evaluated relatively to a reference homogeneous fluid having the same chemical potential µ s and particle density ρ 0 Following the general theoretical scheme introduced by Evans [27,28], the density functional F [ρ] can be split into three contributions: an ideal term, an external potential term and an excess free-energy term accounting for the intrinsic interactions within the fluid, with the following expressions of the first two terms There are several ways of arriving at an exact expression of the excess free-energy, i.e., using an adiabatic perturbation of the pair potential (the so-called adiabatic connection route in electronic DFT), of the external potential, or of the density itself. A conventional approximation is to express the excess term as an expansion around the homogeneous density ρ 0 The first term is the (two-body) direct correlation function (DCF) of the homogeneous solvent, that depends on r 12 = |r 2 − r 1 |, and can be thus denoted as c(r 12 ; ρ 0 ). We define the so-called bridge 33005-2 functional in terms of the higher-order direct correlation functions which thus starts with a cubic term in ∆ρ. Setting F B [ρ] = 0 corresponds to the so-called homogeneous reference fluid (HRF) approximation. It can be shown to be equivalent to the hypernetted chain (HNC) approximation in integral equation theories [29]. The input of the theory is thus a direct correlation function of the pure solvent, which can be extracted from simulation or experimental data by measuring the total correlation function h(r) = g(r) − 1 and solving subsequently the Ornstein-Zernike equation, i.e., in Fourier space: (2.8) χ(r) is the structure factor, or the density susceptibility, measuring density-density correlations at a given distance in the fluid. The excess free energy can thus be also expressed in terms of the inverse susceptibility Minimization of equation (2.3) with respect to ρ gives the equilibrium density

Chandler's Gaussian field theory
Along the same lines as above, Chandler considered the case of a liquid of density ρ 0 , characterised by its intrinsic density susceptibility χ(r), containing a solute creating an external potential v(r) outside a hard core that defines an inside volume V in where the density ρ(r) is zero and where by convention v(r) = 0. Chandler writes a gaussian field Hamiltonian for the pure fluid and the partition function of the fluid + solute system as a field integral where the product of delta-functions imposes the constraint of zero-density inside the core. Performing the Gaussian integral exactly, Chandler arrives at the expression of the solvation free energy and, by functional differentiation with respect to the external potential, at the one-particle equilibrium density
One of the main results in Chandler's paper is that the susceptibility of the medium, defined as χ(r 1 , r 2 ) = δ ρ(r 1 ) /δv(r 2 ), is altered by the presence of the hard core and changed from χ(r 1 , r 2 ) = χ(|r 1 − r 2 |) for the infinite medium to an effective susceptibility that is not translationally invariant anymore.

Linearised and partially-linearised HNC approximations and connection to Gaussian field theory
The linearised HNC approximation consists in expanding the ideal term in equation (2.4) at dominant order in ∆ρ so that the functional to be minimised becomes In the presence of a solute with a hard repulsive core [very positive values of the potential v(r)], such approximation will obviously fail to give an exponentially vanishing density inside the core. As considered by Chandler above, this approximation should be complemented by constraints imposing ρ(r) = 0 within the inside volume V in . There are two ways to impose those constraints. The first one, not necessarily the easiest one, is to introduce an auxiliary Lagrange multiplier field λ(r) and minimise the following constrained functional with respect to ρ(r) and λ(r) Thus, the minimisation equations are as follows: These equations can be readily solved by linear algebra to give an equilibrium density that is equivalent to the one in equation (3.4). Replacement in equation (4.3) does give the equilibrium free energy of equation (3.3), except the last log-of-determinant term that includes a measure of the fluctuations that is absent in the functional approach. Numerical estimations shows that it can be safely neglected with respect to the other terms. We conclude that the Chandler's Gaussian field approach is, up to a small log-term correction in the energy, equivalent to a DFT approach with a linearised HNC approximation. From a DFT perspective, however, a natural way to account for the constraint is to minimise the functional outside the core only, i.e., for r ∈ V out . The functional can thus be limited to the outside region and written as This functional can be easily numerically minimised on a three-dimensional grid using for example a quasi-Newton minimiser such as L-BFGS [69] to yield the equilibrium density ρ eq and the associated free energy. Since the above functional is bilinear in ρ(r), the formal solution can be also obtained by matrix inversion, i.e., outside the core This solution looks quite different from that in equation (3.4); in the appendix below it is shown that the two formulas are in fact equivalent. Thus, we arrive at the main conclusion of this paper: the rather involved formal solutions of the Gaussian field approach (equivalent to a functional minimization with Lagrange multipliers, as seen above), which involves the necessity to numerically invert the matrix χ in inside the core and then to perform a double multiplication of this matrix with χ, can be replaced by a simple numerical minimisation of the LHNC functional (4.7) outside the hard core. The basic input is the homogeneous bulk inverse susceptibility χ −1 (r 12 ) [or equivalently, the homogeneous bulk DCF c(r 12 ; ρ 0 )], with no interference whatsoever with the introduction of hard-core conditions. The bulk inverse susceptibility applies everywhere, inside and outside the hard core. The fact that, as noted by Chandler, the introduction of such hard-core boundaries modifies the apparent susceptibility of the medium outside the core is a consequence that applies to the LHNC-DFT approach as it does for the GFT one. It should be also valid at a HNC level; this effect can be measured numerically as χ(r 1 , r 2 ) = δ ρ(r 1 ) /δv(r 2 ) -indeed not an easy task on a 3D spatial grid.
In the following we test the HNC, LHNC (equivalent to Gaussian field theory), and PLHC for the solvation of a Lennard-Jones sphere of an increasing diameter in a Lennard-Jones liquid, in comparison with the reference Monte-Carlo generated by Lazaridis [70]. The LJ solvent is characterised by a particle diameter σ 0 and reduced thermodynamic conditions ρ * = 0.85, T * = 0.88. In figure 1, we display the solvent structure for 3 solute diameters, σ/σ 0 = 0.2, 1, and 2, respectively. The DFT results were obtained by direct functional minimisation using a home-made spherical 1D code. The hard-core volume for LHNC was identified to the void region obtained after HNC minimisation [ρ(r) < ρ min , a fixed, very small value]. The first observation is that none of the approximations is either perfect or clearly off. Apart from the smaller solute, it is seen that the HNC approximation tends to underestimate the first-peak position and overestimate its height. The second observation is that, surprisingly, LHNC and PLHNC give undistinguishable results; for both, the first peak appears now too low for the smaller solutes and has a correct height but with a shift in position for the biggest, as in HNC. PLHNC can be qualified as a better theory since the hard core is defined and handled automatically by the functional. The situation gets really worse when going to the solvation free energies. In figure 2, we compare the results of the 3 approximations when increasing progressively σ/σ 0 to the simulation results of Lazaridis. All of them are off by a large factor and in nearly the same way. The problem has been clearly identified [56,71,72]: all those HNC variants give an apparent pressure which is way too high with respect to the exact pressure, P exact , of the LJ fluid, and thus a spurious ∆P∆V contribution where ∆P = P HNC − P exact , and ∆V is the solute partial molar volume -close to, but not identical to the inside volume V in of the solute. This can be corrected by adding an empirical pressure correction, −∆P∆V, to the DFT-HNC (or LHNC, or PLHNC) free energy [56,71]. Herein below we switch to a more fundamental correction for Lennard-Jones that involves a hard-sphere bridge functional.

Hard-sphere bridge correction
Building the thermodynamics of the Lennard-Jones fluid by taking a suitable hard-sphere fluid as a reference is indeed a classic in liquid-state theory and is at the basis of the Van der Waals theory of fluids. A variant of this idea is to approximate the bridge functional in equation (2.6) by a hard sphere bridge (HSB) functional introduced by Rosenfeld as a universal bridge function [73][74][75].
Here, F HS exc [ρ(r)] represents the one-component hard-sphere excess functional which, up to a very good approximation, can be taken as the fundamental measure theory (FMT) functional of Rosenfeld [40] and Kierlik and Rosinberg [41,42]. The fourth term involves the direct correlation function of the HS fluid at the same density, i.e., Note that defined as in equation (5.1), F HB B [ρ(r)] carries an expansion in ∆ρ of the order 3 and higher which corrects the second order expansion of the excess free energy in equation (2.6).
We show in figure 3 that this HNC+HSB theory works much better than the HNC variants for the prediction of solvation properties of dissolved molecular objects. There we again compare the solvation free energy of the growing LJ sphere to the Monte-Carlo results of Lazaridis [70] using different HS diameters, d. It can be seen that the results are extremely sensitive to the choice of d, and that the best agreement is obtained for d = 1.014σ (indeed close to 1, that would be the initial guess value). For that value, we have plotted in figure 4 the solvent density, g(r) = ρ(r)/ρ 0 , obtained for solute of different sizes by direct MD simulations that we have generated, or by DFT in the HNC or HNC+HSB approximation. It can be seen that the addition of the hard-sphere bridge greatly improves the results compared to the HNC approximation and furthermore yields a very good structure.

Conclusion
In this paper we have shown a close connection between the Gaussian field theory of solvation introduced by Chandler in [65] and density functional theory in a linearised HNC approximation. Chandler's formulae for the solvation density around hard solutes and the associated solvation free energies can be recovered by minimising the LHNC functional with constraints imposed through an auxiliary Lagrange multiplier field. A simpler but equivalent formulation arises when minimising the functional outside the hard core only. Both theories share with the full HNC approximation, or the intermediate PLHNC approximation, the same caveat of greatly overestimating the solvation free energy of dissolved objects. Chandler was indeed aware of these limitations and provided further improvements based on the coupling of GFT at the microscopic scale to a lattice gas model having a correct macroscopic behaviour at larger scales [68]. In DFT, improvements can be made by considering a bridge functional beyond the second order expansion in density. For the Lennard-Jones solvent, the natural bridge that emerges is that of a reference hard fluid, whose hard-sphere diameter should be optimised. The extension of such an approach to molecular liquids, such as water, has been proposed with some success [51,63]. This remains to be further explored and improved -since a water molecule is definitely not a spherical entity.
The interlink between density functional theories and other versions of liquid-state field theories, such as those developed by Jean-Pierre Badiali and his Parisian and Ukrainian collaborators along the years, is also a very interesting subject that merits to be explored in depth in the future.

A. Connection between "inner" and "outer" DFT formulations, and
Chandler's GFT

A.1. Notation
Herein below we will use a discrete matrix notation for the fields and associated functionals. Let V be the liquid volume and be decomposed into an inside volume V in occupied by the hard-sphere solute and the remaining volume V out . We define the functions on a finite three-dimensional grid. Let m points lie inside the solute and n points outside. In that case, the one-variable functions, like density, can be represented as vectors of size (m + n) × 1 (e.g. ρ). The two-variables functions, e.g., the susceptibility function χ(r 1 , r 2 ) are represented as matrices (m + n) × (m + n) (e.g X). Then, the convolution can be represented as a matrix multiplication, e.g., where ∆v is the elementary volume which corresponds to each discretisation point. For simplicity, we will take below ∆v = 1. Let the density inside the solute be ρ in (m × 1 vector), the density outside the solute ρ out (n × 1 vector). The free energy functional can be defined as follows: Here, ∆ρ = ρ − ρ 0 , X is a susceptibility matrix where X in is m × m, X inter is m × n, X out is n × n. It is important to note that, for example, The above functional should be minimised with the constraint ρ in = 0, ∆ρ in = −ρ 0 . There are two approaches to do this: Lagrange multiplier minimisation or restrained minimisation in the outer volume. We show herein below that the two approaches are equivalent to each other and give the same results as Chandler's Gaussian field theory in [65].

A.2. Lagrange multipliers minimization
To perform the minimisation using Lagrange multipliers we add −λρ in to the functional: From the necessary minimization conditions ∂F ∂λ in = ρ in = 0 (A.6)

33005-9
and ∂F ∂ρ in = (X −1 ) in ∆ρ in + (X −1 ) inter ∆ρ out = λ in , the last two equations can be rewritten as From this, we find ∆ρ and the relations: Using the first equation we find Inserting this into the second equation: This is exactly Chandler's Gaussian field expression, equation (3.4), in discretised form [with the understanding that χ −1 in = ( χ in ) −1 ]. Injecting this formula into equation (A.5) also gives the same expression as Chandler for the equilibrium solvation free-energy, equation (3.3), except the last logarithm-ofdeterminant term.

A.3. Direct minimization in outer volume (reduced number of variables)
Instead of performing the minimisation with the Lagrange multipliers, we can minimise the reduced functional which depends only on ρ out : Taking the derivative (X −1 ) out ∆ρ out − (X −1 ) T inter ρ 0 in + βv out = 0 (A. 15) and 16) To see that this is the same as (A.13) we need to invert the matrix X. To do it, let us define

33005-10
By the definition of the inverse matrix we have where I is an identity matrix of appropriate size. We have the following equations: Multiplying the first by A −1 : .19) and inserting this into the third equation: (Here, we use A = A T , C = C T , which is true since X is symmetric). Now, from the last equation in (A.18) Inserting here the expression of Y: We can further simplify this expression. We first express the identity matrix I as Inserting this into the expression of Z we have Cancelling B T A −1 B and C −1 C we get Returning to the expression (A.16) Opening the brackets ∆ρ out = −B T A −1 ρ 0 in − Cβv out + B T A −1 Bβv out = B T A −1 (−ρ 0 in + Bβv out ) − Cβv out (A. 31) or, returning to the original definitions: ∆ρ out = X T inter (X in ) −1 (−ρ 0 in + X inter βv out ) − X out βv out (A. 32) which is the same as (A.13). This terminates the proof for the equilibrium density. The same equivalence can be proved for the equilibrium solvation free-energy.