Nonlocal density functional theory of water taking into account many-body dipole correlations: Binodal and surface tension of 'liquid-vapour' interface
Yu. A. Budkov A. L. Kolesnikov
cond-mat.soft





[]ybudkov@hse.ru



School of Applied Mathematics, Tikhonov Moscow Institute of Electronics and Mathematics, National Research University Higher School of Economics, Tallinskaya st. 34, 123458 Moscow, Russia
G.A. Krestov Institute of Solution Chemistry of the Russian Academy of Sciences, Academicheskaya st., 1, 153045 Ivanovo, Russia

[]kolesnikov@inc.uni-leipzig.de
Institut für Nichtklassische Chemie e.V., Permoserstr. 15, 04318 Leipzig, Germany

In this paper we formulate a nonlocal density functional theory of inhomogeneous water. We model a water molecule as a couple of oppositely charged sites. The negatively charged sites interact with each other through the Lennard-Jones potential (steric and dispersion interactions), square-well potential (short-range specific interactions due to electron charge transfer), and Coulomb potential, whereas the positively charged sites interact with all types of sites by applying the Coulomb potential only. Taking into account the nonlocal packing effects via the fundamental measure theory (FMT), dispersion and specific interactions in the mean-field approximation, and electrostatic interactions at the many-body level through the random phase approximation, we describe the liquid-vapour interface. We demonstrate that our model without explicit account of the association of water molecules due to hydrogen bonding and with explicit account of the many-body electrostatic interactions at the many-body level is able to describe the liquid-vapour coexistence curve and the surface tension at the ambient pressures and temperatures. We obtain very good agreement with available in the literature MD simulation results for density profile of liquid-vapour interface at ambient state parameters. The formulated theory can be used as a theoretical background for describing of the capillary phenomena, occurring in micro- and mesoporous materials.     

Nonlocal density functional theory of water taking into account many-body dipole correlations: Binodal and surface tension of 'liquid-vapour' interface
    A.L. Kolesnikov
    December 30, 2023. .../latex/6he-sdm-2/6he-sdm-2-3.tex
=======================================================================================================================================================



§ INTRODUCTION

Theoretical description of water adsorption on micro- and mesoporous materials is a challenge problem for modern physical chemists and chemical engineers. Its great importance is due to numerous industrial applications and fundamental significance of describing water in a nanoconfinement. The examples of technological applications, where the description of water adsorption is highly relevant, are: characterizing of micro- and mesoporous materials at ambient conditions (i.e. determination of pore size distribution, surface area, accessible volume) <cit.>, modelling of water purification from toxic compounds (e.g., ions of heavy metals <cit.>), description of mechanical stability of construction porous materials (concrete, wood, paper, etc.) during capillary condensation/evaporation cycles<cit.> , etc. All these examples have one common feature - the influence of inhomogeneity on the process.


Obtaining a correct description of inhomogeneous water requires a reliable theoretical model based on the first principles of statistical mechanics. Such a theoretical model must account for the electrostatic interactions (including short-range specific interactions, attributed to electron charge transfer) between water molecules and must be based on the nonlocal functional theory to describe correctly the liquid-vapour interface and temperature behavior of the surface tension. Despite the fact that up to now several theoretical models have been formulated <cit.>, none of them satisfy the requirements formulated above.


In paper <cit.> the authors formulated a density functional theory (DFT) taking account of the molecular structure of water within the TIP4P model <cit.> for describing the water liquid-vapour interface. The authors took into consideration the universal intermolecular interactions by the Lennard-Jones pair potential, whereas the electrostatic interactions between sites – through the Coulomb potentials. To take into account for the electrostatic interactions in the total free energy functional, the authors used the mean-field approximation, expanding the anisotropic potential into the multipole series, truncated by the fifth order. They applied the local density approximation with the Weeks Chandler Andersen (WCA) procedure to the dispersion interactions <cit.>. Though the authors obtained rather satisfactory agreement with the experimental liquid-vapour coexistence curve (binodal), the values of the surface tension at all the temperatures were  highly overestimated. Such a discrepancy can be explained in two ways. Firstly, the authors used the local DFT and, thus, did not take into account the nonlocal packing effects which should be important for the water molecules, situated at the interface. Secondly, the authors did not consider the many-body electrostatic correlations of the water molecules which must be significant in condensed liquid phase. Indeed, while for a vapour phase the electrostatic correlations manifest themselves through the effective Keesom pairwise interactions, for the liquid condensed phase it is necessary to take into account the higher electrostatic correlations <cit.>.  Nevertheless, from the result of the manuscript<cit.>
one can conclude that the main electrostatic contribution comes from the short-range part of the Coulombic potential. The author shows that short-range attractive and repulsive interactions play the crucial role in the structure properties of polar and associating pure fluids. Later in the work <cit.>
the authors pointed out that not accounting for the long-range electrostatic interactions leads to errors in the system with non-uniform geometries, however mostly errors arise in electrostatic properties.




It is also necessary to mention the phenomenological density functional theories of water which do not take into account the electrostatic interactions explicitly <cit.>. In paper <cit.> the authors formulated a molecular DFT of water, which allowed authors to predict with good accuracy the  temperature of freezing at atmospheric pressure. To perform the numerical calculations, the authors used the experimental pair correlation functions oxygen-oxygen, oxygen-hydrogen, and hydrogen-hydrogen. Thus, despite the success of formulated molecular DFT, its application to different molecular systems requires the external experimental data regarding the site-site correlation functions (from X-ray scattering or full-atomistic computer simulations). In paper <cit.> the authors used a similar molecular phenomenological DFT approach to describe the liquid-vapour interface of water at a temperature of 298 K. Despite the fact that the authors obtained a good fitting for experimental values of densities of the coexisting phases and surface tension, it remained unclear how this approach could describe these quantities at the other state parameters. Moreover, the formulated theory deals with the polynomial approximation for the excess free energy, whose phenomenological coefficients are not related to any statistical theory. In paper <cit.> the authors developed a phenomenological DFT which operates with the position-orientation number density of structured fluids. Despite the fact that this theory is based on the bulk equation of state, describing the anomalous behavior of water below 4^0 C and taking explicit account of the hydrogen bonding between the water molecules, whereas the free energy functional accounts for the nonlocality through the gradient term, the surface tension obtained is twice its experimental value. As the authors mentioned <cit.>, this discrepancy is most probably due to the simplicity of the gradient correction used. In paper <cit.> a self-consistent DFT, taking explicit account of hydrogen bonding through the Statistical Associating Fluid Theory (SAFT) <cit.>, is applied to investigating the phase behavior and surface tensions of water and aliphatic alcohols. The authors showed that for the bulk phases, their theory is reduced to an equation of state that provides an accurate description of saturation pressures as well as vapor-liquid phase diagrams. Near the critical region, the long-range fluctuations were taken into account using a renormalization group theory. It is worth mentioning the similar SAFT-based density functional approach <cit.>, where the authors not only described the saturation pressure and surface tension, but also analyzed the effective interactions between the hydrophobic hard rods, immersed in liquid water. In the manuscript<cit.>, the authors used the SAFT-VR density functional theory in order to describe the vapor-liquid interface of associating and non-associating molecules, including water. The functional treats short-range repulsion, chain and association contributions in the local density approximation. A good description of both binodal and surface tension was achieved by including interfacial data in the optimization scheme <cit.>. In the work <cit.>, the authors developed classical density-functional theory of rigid-molecular fluid and applied it to calculations of thermodynamic and structural properties of water. The used functional contains the hard-spherical contribution (White Bear mark II), contribution of the attractive Van der Waals interactions (on the weighted-density approximation level), and electrostatic contribution, described within the mean-field approximation. The using hard-sphere diameter as a fitting parameter allowed authors to describe surface tension with good accuracy. 


Despite the fact that in the mentioned papers the authors obtained a very good fitting for the saturation pressure and liquid-vapour surface tension, these are only a few theories accounting explicitly the electrostatic interactions between water molecules. The existing SAFT-based density functional theories can be considered as good tools for chemical engineering applications, but cannot answer the physically reasonable question: How big the contribution of the electrostatic interactions between water molecules to the surface tension of the liquid-vapour interface? Moreover, in principle, the SAFT-based theories in their present form cannot elucidate the role of the electrostatic correlations in the capillary phenomena taking place with confined water. 

One of the simplest water models is the SPC/E  - a 3-site model. Each site carries a point-like charge and additionally, the Lennard-Jones interaction potential between oxygen atoms is applied. Despite the fact that SPC/E model does not take into account hydrogen bonding explicitly, it can reproduce some thermodynamic and structural properties with sufficient accuracy <cit.>. Inspiring the fact that such a simple model allowed the authors to describe equilibrium water in MD simulations, we will formulate a nonlocal simple DFT approach for the description of inhomogeneous liquid water. Describing water molecules as spherically symmetric dipolar particles, having two oppositely charged sites and taking into account the many-body electrostatic correlations, universal inter molecular interactions, and short-range specific interactions, related to electron charge transfer, we will describe with good accuracy the binodal and the surface tension at ambient conditions. In the framework of the formulated theory, we will elucidate the role of different inter molecular interactions in the chemical potential of water at the liquid-vapour interface.




§ THEORY OF BULK WATER

Let us formulate a simple statistical theory of liquid (vapour) water in the bulk phase. We describe each water molecule as two sites with charges ± q, separated by distance l, so that the dipole moment is p=ql=1.85 D.  We assume that sites with charges -q<0 interact with each other through the Coulomb potential, the Lennard-Jones (LJ) potential, and the attractive square-well short-range potential (see the description below). The sites with a charge +q>0 interact with all types of sites through the Coulomb potential only. We would like to note that the LJ potential describes the excluded volume and dispersion interactions, while the square-well potential takes into account at the primitive level the short-range specific interactions between water molecules. In other words, in this study we do not take explicit account of the association between water molecules, resulting in the formation of a hydrogen bond network. Instead, we replace the asymmetric chemical interactions by the effective spherically symmetric short-range square-well potential. Such a simplification is quite acceptable, because in this study we do not consider the fine effects, related to the network of the hydrogen bonds (see, for instance, <cit.>). As we will show below, such a primitive description of the effect of chemical interactions will allow us to describe the liquid-vapour coexistence curve and surface tension of water at the normal conditions that are the goal of the present study. Note that explicit account of the association effect was made in the general context in papers <cit.>. We also neglect the effect of static electronic polarizability of water molecules, as it is much less than the orientation polarizability, related to the permanent dipole moment at the normal temperature. Due to the fact that our model is the first statistical theory of liquid water taking into account dipole correlations of molecules at the many-body level within the random phase approximation, we have neglected for simplicity the orientation nonlinear effects, such as formation of chain-like clusters in accordance with "head-to-tail" mechanism taking place usually in magnetic and ferroelectric fluids <cit.>. As we will show below, such an assumption allows us to describe successfully the phase coexistence curve and surface tension of liquid-vapour interface of water at the ambient conditions.


Bearing in mind all the above written model assumptions, we can write the density of Helmholtz free energy of water as follows:

    f(ρ,T) = f_id(ρ,T)+f_ex(ρ,T),

where 

    f_id(ρ,T)=ρ k_BT(ln(ρλ^3)-1)

is the free energy of the ideal gas.
The excess free energy can be written in the following form

    f_ex(ρ,T)=f_sr(ρ,T) + f_el(ρ,T),

where

    f_sr(ρ,T)=ρ k_B T (-ln(1 - η) + 3η/1 - η + 3 η^2/2(1 - η)^2)+1/2B ρ^2,

is the contribution of the short-range interactions, including the excluded volume interactions, dispersion interactions, and short-range specific interactions; B is the parameter of attractive interactions accumulating the contributions from WCA attraction tail

    V_WCA(|𝐫|)=
    {-ϵ    if    r < 2^1/6σ
    4ϵ[(σ/r)^12 - (σ/r)^6]    if    2^1/6σ < r < r_c
    0    if    r > r_c
    }

and 

    V_spc(|𝐫|)=-ϵ_swΘ(σ_sw/2-|𝐫|),

where r_c is the cutoff of LJ pairwise potential of interactions; Θ(r) is the Heviside step-function; σ and ϵ are, respectively, the size and energy parameters of LJ potential; ϵ_sw and σ_sw are, respectively, the energy and size parameters of the square-well potential. The latter describes the contribution of short-range specific interactions <cit.> which in our case is related to the electron charge transfer of the water molecules. We assume for our calculations that r_c=5σ. Such an assumption allows us to make our DFT calculations (see the next section) less time-consuming <cit.>. The first term in (<ref>) describes the contribution of the excluded volume interactions of the hard spheres within the Percus-Yewick approximation; η= π d_BH^3ρ/6 is the packing fraction of hard spheres with the effective Barker-Henderson diameter, determined by the following Pade approximation <cit.>

    d_BH = σ1.068 ϵ/k_B T + 0.3837/ϵ/k_B T + 0.4293.

 The second term in (<ref>) describes the total contribution of the attractive interactions, thus B can be defined as:
 
    B =∫ d𝐫(V_WCA(|𝐫|) + V_spc(|𝐫|) ) = -32√(2)/9πϵσ^3 + 16/3πϵσ^3[ (σ/r_c)^3 - 1/3(σ/r_c)^9] -πϵ_swσ_sw^3/6.


The contribution of the electrostatic interactions, which in our case are reduced to the short-range dipole-dipole interactions, can be described by the free energy of the dipolar hard spheres (see Appendix 1):

    f_el = -k_B T/l^3(1-3/4α)σ(y),

where l is the dipole length and y = ρ p^2/3 ϵ_0 k_B T, p is the dipole moment. The auxiliary function

    σ(y) = √(6)/4π[2(1+y)^3/2-2-3y]

is also introduced; the compressibility factor α of the hard spheres (see Appendix 1) in Percus-Yevick approximation has the following form

    α=η(4-η)(2+η^2)/(1+2η)^2.

 Before that, we talk only about the thermodynamic properties, however, it is interesting to study also the structural ones. Here, we compare calculated structure factor and pair correlation function of the homogeneous water at 298 K with presented literature experimental data and DFT calculations. In general, there are two methods to calculate pair correlation function within DFT - test particle method and method base on the Ornstein-Zernike (OZ) equation. Recently <cit.>, it was shown that the first one gives more accurate results and especially enforce the correct behavior at small distances. However, we unable to use it due to the lack of explicit known form of the effective electrostatic pair potential, thus we will use OZ-based method (see  Appendix IV). The Fourier transform of structure factor in the homogeneous limit can be calculated by the standard expression <cit.>:

    S(k) = 1/1 - ρ c^(2)(k),

where k=|k|; the pair correlation function can be expressed through structure factor as:

    g(r) = 1 + 1/2 π^2 r ρ∫_0^∞ dk ksin(kr)(S(k) - 1).







Figures <ref> and <ref> represent the comparison between structure factor calculated within this model, molecular density functional theory <cit.>, and X-Ray scattering experiments <cit.>. We obtain only partial qualitative agreement between our model and X-Ray obtained structure factor. Our results does not demonstrate the split of first pick <cit.> and the next picks, on bigger k, are shifted to the left. The first pick approximately corresponds to the average position of the first experimental ones and the amplitudes of picks are in agreement with experimental data. By means of Eq.(<ref>) we calculated the pair correlation function. The Fig. <ref> represents X-Ray obtained <cit.> and calculated within cDFT <cit.> g_OO,  and g(r)Θ(r - d_HS) calculated in the present study. We plotted only part of the whole pair correlation function (cutted approximately on the distance of effective hard-sphere diameter), due to the numerical artifacts on the low distances. The pair correlation function does not go to zero in the core region. Such behavior was already previously reported in the literature <cit.> and can be attributed to the accuracy of the used method. Nevertheless, the obtained pair correlation function is in a good agreement with the results of molecular density functional theory except the position of the first maximum. The latter is shifted to the higher values by approximately 0.3 . However, we would like to note that both theories fail to reproduce the real structure of water and demonstrates a significant shift of the second and third picks.




§ NONLOCAL DENSITY FUNCTIONAL THEORY OF INHOMOGENEOUS WATER

Based on the bulk theory which was formulated in the previous section, let us formulate a density functional theory of inhomogeneous water. We start from the grand thermodynamic potential of inhomogeneous water in external potential field with the potential energy V_ext(r), which can be written in the following form

    Ω[ρ(𝐫)] = F_id[ρ(𝐫)] + F_ex[ρ(𝐫)]+∫_V d𝐫ρ(𝐫) V_ext(𝐫) - μ∫_V d𝐫ρ(𝐫),

where V is system volume, F_id[ρ(𝐫)] is the free energy of the ideal gas, F_ex[ρ(𝐫)] is the excess free energy of water, μ is the chemical potential and ρ(r) is the single-particle density. Note that within this consideration ρ(r) is the average density of the centers of mass of the water molecules. Thus, we construct the effective nonlocal DFT in terms of the simple fluid theory <cit.> for water that is a molecular liquid in its nature. Such a simplification can be considered as an example of coarse-graining. The ideal gas free energy is

    F_id[ρ(𝐫)] = k_B T ∫_V d𝐫ρ(𝐫)[ln(λ^3ρ(𝐫)) - 1]

and the excess free energy, in turn, consists of several parts:

    F_ex[ρ(𝐫)] = F_sr[ρ(𝐫)] + F_el[ρ(𝐫)],

where

    F_sr[ρ(𝐫)] = F_fmt[ρ(𝐫)] + F_att[ρ(𝐫)]

with the Helmholtz free energy F_fmt[ρ(𝐫)] of hard spheres with the effective BH diameter d_BH within the fundamental measure theory (FMT) <cit.>, namely

    F_fmt[ρ(𝐫)] = ∫_V d𝐫Φ({n}),

where the free energy density is

    Φ({n}) =k_BT(-n_0ln(1 - n_3) + n_1 n_2 - 𝐧^(v)_1 𝐧^(v)_2/1 - n_3 + n_2^3 - 3 n_2 𝐧^(v)_2 𝐧^(v)_2/24 π (1 - n_3)^2),

which depends on six weighted densities:  

    n_α(𝐫) = ∫_V d𝐫^'ω_α(𝐫 - 𝐫^') ρ(𝐫^'),

where ω_α(|𝐫 - 𝐫^'|) are the weight functions. Four of them are scalar values: ω_2(𝐫_12) = δ(d_BH/2 - |𝐫_12|), ω_3(𝐫_12) = Θ(d_BH/2 - |𝐫_12|), ω_1(𝐫_12) = ω_2(𝐫_12) / (2π d_BH), ω_0(𝐫_12) = ω_2(𝐫_12) / (π d_BH^2),  and two are vector values: ω^(V)_2(𝐫_12) = 𝐫_12/r_12δ(d_BH/2 - |𝐫_12|) and ω^(v)_1(𝐫_12) = ω^(v)_2(𝐫_12)/(2π d_BH). The functional F_att[ρ(𝐫)] is the contribution from the given above attractive potentials V_WCA(r) and V_spc(r):

    F_att[ρ(𝐫)] = 1/2∫∫_V d𝐫_1d𝐫_2ρ(𝐫_1)ρ(𝐫_2)[V_WCA(|𝐫_1 - 𝐫_2|) + V_spc(|𝐫_1 - 𝐫_2|)].



Note that we do not know the real electrostatic free energy functional of dipolar hard spheres. On the other hand, due to the fact that thermodynamic properties of strongly inhomogeneous confined polar fluids must be very different from those are in the bulk phase, we cannot use the local density approximation for the electrostatic free energy functional. Nevertheless, since we know the approximate expression for the electrostatic free energy of the dipolar hard spheres system for the bulk phase, we can construct the phenomenological weighted density functional based on it, using the Curtin-Ashcroft-Tarazona approach <cit.>. We would like to note, that using of the weighted density functional approach can be justified by the fact that electrostatic interactions between water molecules in liquid phase manifest themselves as the effective short-range dipole-dipole interactions [This is related to the fact that the bulk electrostatic free energy of dipolar hard spheres system can be expanded into the power series on density (virial expansion), in contrast to the ionic systems, for which the electrostatic contribution to the total free energy is not an analytic function of density. As is well known, non-analytical behavior of electrostatic free energy of the ionic systems is determined by the long-range electrostatic interactions and translation entropy of ions <cit.>. However, if the ionic groups form the electrically neutral clusters, as in dipolar <cit.> or quadrupolar <cit.> fluids, the effective interactions between particles become short-range and, thereby, analyticity of the electrostatic free energy recovers.]. Thus, following the idea of the weighted density functional theory <cit.>, we treat the electrostatic free energy functional as follows:

    F_el[ρ(𝐫)] = ∫_V drρ(𝐫) ϕ_el(ρ̅(𝐫)),

where ϕ_el(ρ̅(𝐫)) = f_el(ρ̅)/ρ̅ is the free energy per fluid particle, depending on smoothed density, which, in turn, is determined in the following way:

    ρ̅(𝐫) = ∫_V d𝐫^'ρ(𝐫^') ω_el(|𝐫^' - 𝐫|),

where ω_el(|𝐫|) = 3/(4π R_w^3)Θ(R_w -|𝐫|) is the phenomenological weighted function with the phenomenological scale parameter R_w determined the range of smoothing.




The equilibrium density profile is obtained from the minimization of the grand thermodynamic potential (<ref>), i.e. from the Euler-Lagrange equation

    δΩ[ρ(𝐫)]/δρ(𝐫)=0

or

    ρ(𝐫) = ρexp[β(μ_ex - V_ext(𝐫)) + c^(1)(r)],

where

    c^(1)(r)=-δ(β F_ex[ρ])/δρ(r)=c^(1)_fmt(r) + c^(1)_WCA(r) +c^(1)_spc(r) + c^(1)_el(r)

is the one-particle direct correlation function; β=(k_BT)^-1 is the inverse thermal energy; μ_ex=μ_ex(ρ,T) is the bulk excess chemical potential, given in Appendix II. The contributions to the one-particle direct correlation function can be written as

    c^(1)_el(𝐫) =-δ(β F_el[ρ])/δρ(r)= -βϕ_el(ρ̅(𝐫)) - β∫_V d𝐫^'ρ(𝐫^') [ μ_el(ρ̅(𝐫^'))/ρ̅(𝐫^') - ϕ_el(ρ̅(𝐫^'))/ρ̅(𝐫^')] ω_el(|𝐫 - 𝐫^'|),


    c^(1)_WCA(𝐫) = -δ(β F_WCA[ρ])/δρ(r)=-β∫_V d𝐫^'ρ(𝐫^')V_WCA(|𝐫 - 𝐫^'|),


    c^(1)_spc(𝐫) =-δ(β F_spc[ρ])/δρ(r)= -β∫_V d𝐫^'ρ(𝐫^')V_spc(|𝐫 - 𝐫^'|),


    c^(1)_fmt(𝐫) =-δ(β F_fmt[ρ])/δρ(r)= -β∑_α∫_V d𝐫^'∂Φ({n})/∂ n_αω_α(𝐫^' - 𝐫)

Note that the formulated DFT for the bulk phase, where V_ext(r)=0 and, thus, ρ(r)=ρ=const, transforms into the bulk theory discussed in the previous section.





§ NUMERICAL RESULTS AND DISCUSSIONS

Now we will consider the application of the model to the description of water liquid-vapour interface. We used the theory formulated above to fit the experimental values of coexisting densities and surface tension at the ambient temperature T=298 K. In order to calculate the densities of coexisting phases, we used the following system of equations, consisting of mechanical and chemical equilibrium conditions, i.e.

    P(ρ_v,T) = P(ρ_l,T)   and  μ(ρ_v,T) = μ(ρ_l,T),

where ρ_v and ρ_l are, respectively, the densities of the coexisting vapor and liquid phases. To calculate the equilibrium density profile, we used slit geometry with the length H = 30σ - 60σ and discretization step along z axis ≈ 0.02σ. The surface tension was calculated by the following relation

    γ=Ω[ρ(z)]/A + PH,

where Ω[ρ(z)]/A is the part of the grand thermodynamic potential per unit area A corresponding to the liquid-vapour and P is the pressure in the bulk liquid and vapour phases. The expressions for the pressure and chemical potential are given in Appendix II. The detailed description of DFT main equations in slit geometry is presented in Appendix III. Note that the obtained pressure at T=298 K is deviated from the experimental value (≈ 3141.7 kPa) less than 1 %. It is worth noting that the latent heat of vaporization at 298 K is ≈ 40.89 kJ/mol, with experimental value 44 kJ/mol. Thus, fitting the densities ρ_v and ρ_l and the surface tension γ yields the following set of microscopic parameters:


Table I: The obtained set of parameters for water.





 σ (Å)     ϵ (kcal/mol)     σ_sw (Å)     ϵ_sw(kcal/mol)     l (Å)     R_w (Å)
 
 3.01     0.447     6.02     1.68     2.8     3.01 



Using these model microscopic parameters, we calculated the binodal and surface tension at other temperatures. Fig. <ref> shows the comparison between the binodal, calculated within our DFT, molecular simulations of Alejandre et al. <cit.> and Vega et al<cit.>, and the experimental one <cit.>.




As one can see, the theoretical binodal is in very good agreement with both the MD simulation and the experiment at the ambient state parameters. However, our theory slightly underestimates the critical density and overestimates the critical temperature. However, this disagreement between the theoretical and experimental critical points is not very important for applications discussed in the Introduction. Fig. <ref> demonstrates a comparison between the dependences γ vs T, obtained from our DFT  and SAFT-based DFT <cit.>. The surface tension values calculated within our DFT agree very well with the experimental values at the ambient temperatures and deviate significantly only at sufficiently high temperatures that are close to the critical one. The latter is related to the fact that our theory cannot describe correctly the liquid-vapour equilibrium in the critical point vicinity. As is seen from fig. <ref>, the same picture is observed for the SAFT-based DFT <cit.>. We do not show the surface tension calculated in the works <cit.>. The results obtained in them are almost in a perfect agreement with experimental data.  Fig. <ref> demonstrates the comparison between our DFT, molecular simulations <cit.>, and experimental measurements <cit.>. The surface tension obtained by different authors within one water model (SPC/E) are  different, which could be due to differences in the simulation procedures<cit.>. The results of the present DFT is in quantitative agreement with TIP4P72005 water model. Further, we compared the density profile at the liquid-vapour interface calculated within our DFT the MD simulation with profiles available in the literature <cit.> and ab initio Car-Parinello simulations <cit.>. As Figures <ref> and <ref> show, our theory reproduces the simulation density profiles quite well. Note that DFT density profiles were shifted horizontally in order to match the regions of the steep descent of the simulated ones. 


Despite the good consistent of the calculated surface tensions with experimental and simulation data, the density profiles are significantly steeper than simulation ones. That can be clearly identified by the comparison of the decaying length in the following equation <cit.>:

    ρ(z) = 1/2(ρ_v + ρ_l) - 1/2(ρ_l - ρ_v) tanh(z-z_0/d),

where ρ_v and ρ_l are vapor and liquid densities, z_0 is a position of the Gibbs dividing surface and d is the decay length. The same form was used in the work <cit.>, where the authors reported the value of "10-90" thicknesses (t = 2.1972 d) for density profiles calculated in the SPC/E, TIP4P and TIP4P/2005 water models at 300 K. They are 3.39 , 3.67  and 3.22 , respectively. The value obtained in the present model is approximately 2 times lower ≈ 1.76 .
Also, density profiles demonstrate slightly more pronounced oscillation behavior than those are obtained from molecular simulations. 










It is instructive to estimate the contributions of different intermolecular interactions to the excess chemical potential of liquid water predicted by our statistical model. Fig. <ref> shows these contributions as the functions of temperature, calculated along the liquid-phase branch ρ_l=ρ_l(T) of the binodal. As is seen, the contribution of the universal intermolecular interactions (μ_LJ) and electrostatic contribution (μ_el) almost compensate for each other at all the temperatures, so that the total excess chemical potential in the liquid state region within our model is determined by the contribution of short-range specific interactions (μ_spc). This prediction, in principle, could be verified by ab initio Car-Parinello computer simulations <cit.>.  



§ CONCLUDING REMARKS AND PROSPECTS

We have formulated a nonlocal density functional theory of inhomogeneous liquid water. We considered a water molecule as a couple of oppositely charged sites. The negatively charged sites interact with each other via the Lennard-Jones potential, square-well potential, and Coulomb potential, whereas the positively charged sites interact with all types of sites via the Coulomb potential only. Taking into account the nonlocal packing effects in the framework of the fundamental measure theory, dispersion and specific interactions in the mean-field approximation, and electrostatic interactions at the many-body level through the random phase approximation we have described the liquid-vapour interface. Namely, we have shown that our model without explicit account of the association of water molecules and with explicit account of the many-body electrostatic interactions at the many-body level is able to describe the liquid-vapour coexistence curve and the surface tension at ambient state parameters with good accuracy. 


In conclusion, we would like to discuss the prospects of the formulated density functional theory. At first, this theory could be used as a theoretical background for describing the capillary phenomena, such as wetting/dewetting and capillary condensation/evaporation occurring at the solid surfaces of micro- and mesoporous materials. However, we cannot guarantee that application of this theory to the description of confined water will not change the values of the microscopic parameters of water molecules. However, we believe that the obtained microscopic parameters will be quite close to those obtained in the present study. The density functional theory of confined water will allow us to characterise micro- and mesoporous materials using the experimental adsorption isotherms of water vapour. Secondly, the formulated theory can be used for describing thermodynamic properties of other bulk and confined polar fluids, such as dimethylformamide, aliphatic alcohols, etc. However, these issues are the subject of forthcoming publications.  



§ CONFLICTS OF INTEREST

There are no conflicts to declare.


The model development was supported by the RFBR according to research project No 18-31-20015. The project was partially supported by the RFBR according to research project No 18-29-06008. The authors thank Mikhail Kiselev for fruitful discussions and valuable comments.




§ APPENDIX I: DERIVATION OF HELMHOLTZ FREE ENERGY OF A SYSTEM OF DIPOLAR HARD SPHERES

In this appendix we will briefly consider the fluctuation theory of complex fluids, whose electrically neutral molecules can be modelled as a set of clusters of spatially correlated charged centers with charges q_α. A simplest example of such a fluid is a polar fluid, which we will consider in detail. More specifically, we will derive within the random phase approximation (RPA) an analytical expression for the Helmholtz free energy of the dipolar hard spheres system in the bulk.  We start from fluid partition function, which can be written as the following functional integral over the fluctuations of the local number densities ρ_α(r) of the charged centers as follows

    Z=𝒩∫∏_α^𝒟ρ_αexp[-β F_0[{ρ_α}]-β U_cl[{ρ_α}]],

where β=(k_BT)^-1 is the inverse thermal energy, F_0[{ρ_α}] is the free energy functional of the reference system without Coulomb interactions between the particles; 𝒩 is the normalization constant which will be specified below;

    U_cl[{ρ_α}]=1/8πε_0∫ dr∫ dr^'ρ_c(r)ρ_c(r^')/|r-r^'|

is the energy of Coulomb interactions between the charged centers, ρ_c(r)=∑_αq_αρ_α(r) is the local charge density; ε_0 is the vacuum permittivity. 


Further, we expand the free energy of the reference system into the functional power series near the average densities ρ̅_α=ρ of the charged centers

    F_0[{ρ_α}]≈ F_0[{ρ̅_α}]+k_BT/2∫ dr∫ dr^'∑_α,γG_αγ^-1(r-r^')δρ_α(r)δρ_γ(r^'),

where δρ_α(r)=ρ_α(r)-ρ̅_α(r) is the fluctuation of the local charge densities;

    G^-1_αγ(r-r^')=δ^2/δρ_α(r)δρ_γ(r^')(F_0/k_BT)|_ρ=ρ̅

is the inverse structure operator of the reference system for which we adopt the following approximation

    G_αγ^-1(r-r^')=W_αγ^-1(r-r^')-c_αγ(r-r^'),

where c_αγ(r-r^') is the matrix of direct correlation functions of the sites that are not bonded to each other and 

    W_αγ(r-r^')=ρδ_αγδ(r-r^')+ρ(1-δ_αγ)g_αγ(r-r^')

is the matrix of the structure factors of molecules; g_αγ(r-r^') is the probability distribution function of the distance between the α^th and γ^th sites. Thus, we have the following relation

    Z≈exp[-β F_0[{ρ̅_α}]]𝒩∫∏_α^𝒟ρ_αexp[-1/2∫ dr∫ dr^'∑_α,γS_αγ^-1(r-r^')δρ_α(r)δρ_γ(r^')],

where 

    S_αγ^-1(r-r^')=G_αγ^-1(r-r^')+q_αq_γ/4πε_0 k_BT|r-r^'|

is the inverse structure operator in the random phase approximation. Further, choosing the normalization constant 𝒩 from the condition that at q_α=0 the electrostatic contribution to the free energy is equal to zero and calculating the Gaussian functional integral, we arrive at

    F=F_0+F_el,

where 

    F_el=Vk_BT/2∫dk/(2π)^3(ln(1+ϰ^2(k)/k^2)-ϰ^2(k)/k^2)

is the electrostatic free energy and

    ϰ^2(k)=1/ε_0 k_BT∑_αγq_αq_γG_αγ(k)

is the screening function. Note that we have subtracted from the final expression the electrostatic self-energy of the molecules

    E_self=Vk_BT/2∫dk/(2π)^3ϰ^2(k)/k^2.

G_αγ(k) are the Fourier-images of the structure factors of the reference system which can be calculated from the following matrix relation <cit.>

    G^-1_αγ(k)=W_αγ^-1(k)-c_αγ(k).

Now, we will turn to the theory of dipolar fluids, following from the general theory formulated above. In this case, we consider the dipolar particles as pairs of charges ± q.  The relation for the structure factor of the dipolar particles in the RPA takes the following form

    G^-1(k)=W^-1(k)-C(k)=
    [     1/ρ(1-g^2(k))-c_11(k) -g(k)/ρ(1-g^2(k))-c_12(k); -g(k)/ρ(1-g^2(k))-c_12(k)     1/ρ(1-g^2(k))-c_22(k);                           ]
    ,

where

    C(k)=
    [ c_11(k) c_12(k); c_12(k) c_22(k);         ]

is the matrix of the Fourier-images of the direct correlation functions of the reference system, whereas the structure matrix of dipolar molecules has the following form

    W(k)=
    [      ρ ρ g(k); ρ g(k)      ρ;        ],

where

    g(k)=∫ drg(r)e^-ikr

is the characteristic function corresponding to the probability distribution function g(r). The screening function <cit.> is

    ϰ^2(k)=1/ε_0 k_BT∑_αγq_αq_γG_αγ(k)=2ρ q^2/ε_0 k_BT(1-g(k))Q(k),

where

    Q(k)=1-ρ/2(c_11(k)+
    c_22(k)+2c_12(k))(1+g(k))/1-ρ(c_11(k)+
    c_22(k)+2c_12(k)g(k))+ρ^2(1-g^2(k))Δ(k)

and

    Δ(k)=c_11(k)c_22(k)-c_12^2(k).


Let us consider the reference system with c_11(k)=c(k) and c_22(k)=c_12(k)=0, where c(k) is the direct correlation function of the hard spheres. Thus, one of the sites (site 1) is a center of a hard sphere, while the other site (site 2) is a point-like one which does not correlate with the sites of all the species. Essentially, such a reference system describes a set of hard spheres with a grafted point-like particle which can freely penetrate inside the hard spheres. Thus, in this case we obtain

    Q(k)=1+ρ/2h(k)(1-g(k)),

where

    h(k)=c(k)/1-ρ c(k)

is the correlation function of the hard spheres <cit.>.

Further, using the following model characteristic function

    g(k)=1/1+k^2l^2/6,

the approximation

    Q(k)≈ 1+ρ/2h(0)(1-g(k))=1+Z_0-1/2(1-g(k)),

and taking the integral (<ref>), we arrive at the following relation for the density of the electrostatic free energy

    f_el=-k_BT/l^3Λ(y,α),

where l is the dipole length and the auxiliary functions

    Λ(y,α)=(1-3α/4)σ(y)-3√(6)((α+4)(y^2+2y+y√(1+y))+8(1+√(1+y)))/8π(1+√(1+y))


    +3√(6)exp[α y/4(1+√(1+y))^2](y(4+√(1+y)(α+2))+2y^2+4(1+√(1+y)))/4π (1+√(1+y))

and

    σ(y) = √(6)/4π[2(1+y)^3/2-2-3y]

are introduced. Here, α=1-Z_0=1-ρ k_BTχ_0 with the isothermal compressibility χ_0 of the hard spheres and y=q^2l^2ρ/3ε_0 k_BT. We would like to note that with good accuracy one can use the simplified relation for the electrostatic free energy of the dipolar hard spheres fluid

    F_el=-Vk_BT/l^3(1-3/4α)σ(y).

In the Percus-Yewick approximation

    F_0=N k_B T(ln(λ^3ρ)-1)+N k_B T (-ln(1 - η) + 3η/1 - η + 3 η^2/2(1 - η)^2),

taking into account the equation of state for the reference system

    P_0= ρ k_B T 1 + η + η^2/(1 - η)^3

and the relation for the compressibility χ_0=ρ^-1∂ρ/∂P_0,
one can easily obtain a relation for the "compressibility factor"

    α=η(4-η)(2+η^2)/(1+2η)^2,

where η=π d^3ρ /6 is the packing fraction of the hard spheres; d is the hard sphere diameter; λ is the thermal de Broglie wavelength.
































§ APPENDIX II: PRESSURE AND CHEMICAL POTENTIAL IN THE BULK PHASE

The total pressure P and chemical potential μ are the sum of four contributions:

    P = ρ∂ f/∂ρ-f=P_id+P_ex=P_id+ P_sr+ P_el

and 

    μ =∂ f/∂ρ= μ_id+μ_ex=μ_id+μ_sr+μ_el

where P_id=ρ k_BT and μ_id=k_BTln(ρλ^3) are, respectively, the ideal gas pressure and chemical potential; P_ex and μ_ex are the excess pressure and chemical potential, respectively, which, in turn, can be written through the sums of three contributions.


The ideal gas and short-range interaction contributions to the pressure and chemical potential are, respectively,

    P_id+P_sr= ρ k_B T 1 + η + η^2/(1 - η)^3+ 1/2Bρ^2

and

    μ_id+μ_sr = k_B T ln(ρλ^3) +  k_B T (-ln(1 - η) +  η(14 - 13η + 5η^2)/2(1 - η)^3)+Bρ.


The electrostatic contribution to the chemical potential takes the following form

    μ_el = -k_B T/l^3[(1 - 3/4α)Θ_1 + σ(y)Θ_2]

with the auxiliary functions

    Θ_1 =√(6)p^2/4 πϵ_0 k_B T(√(1 + y) - 1),


    Θ_2 = -π d_BH^3/2(1-η)^3(2+η)/(1+2η)^3.

The electrostatic contribution to the pressure can be calculated by the relation

    P_el = μ_elρ - f_el.





§ APPENDIX III

This Appendix presents detailed information about the contributions to the one-particle direct correlation function for the slit geometry. The direct correlation function of the specific interactions is

    c^(1)_spc(z) =  πβϵ_sw∫_z - σ_sw/2^z + σ_sw/2 dz^'ρ(z^') (σ_sw^2/4 - (z - z^')^2)

and the contribution from the WCA potential c^(1)_WCA(z)

    c^(1)_WCA(z) =  ∫_z - r_c^z + r_c dz^'ρ(z^')g_WCA(z,z^'),

where g_WCA(z,z^') is

    g_WCA(z,z^')=
    {4πϵσ^2 ψ(r_m) -πϵ(r_m^2 - (z^' - z)^2)    if    |z - z^'|  < r_m
    4πϵσ^2 ψ(z^' - z)    if    r_m<|z - z^'|< r_c,
    0    if    |z - z^'|>r_c.}


Above we used the following notations r_m = 2^1/6σ and ψ(x) is

    ψ(x) =  -1/5( (σ/r_c)^10 - (σ/x)^10) + 1/2( (σ/r_c)^4 - (σ/x)^4).









The three independent weighted functions are:

    n_3(z) = π∫_z - R^z + R dz^'ρ(z^') (R^2 - (z - z^')^2),


    n_2(z) = 2π R ∫_z - R^z + R dz^'ρ(z^'),


    𝐧^(v)_2(z) = 2 π𝐞_𝐳∫_z - R^z + R dz ρ(z^') (z - z^'),

where 𝐞_𝐳 is the unit vector along the z-axis and R = d_BH/2.
Due to the fact that only three weighted functions are independent, we can express the hard spheres direct correlation function in the FMT approximation in the following way: 

    c^(1)_fmt(z) =  2π∫_z - R^z+R dz^'[1/2∂Φ/∂ n_3(R^2 - (z^' - z)^2) + R∂Φ/∂ n_2 + ∂Φ/∂ n^(v)_2(z^' - z) ],

where the derivatives are determined as follows

    ∂Φ/∂ n_3 = n_2/4π R^2(1 - n_3) + n_2^2 - (n^(v)_2)^2/4π R (1 - n_3)^2 + n_2^3 - 3n_2 (n^(v)_2)^2/12π(1-n_3)^3,


    ∂Φ/∂ n_2 = -ln(1 - n_3)/(4π R) + n_2/2π R(1 - n_3) + n_2^2 - (n^(v)_2)^2/8π(1-n_3)^2,


    ∂Φ/∂ n^(v)_2 = -n^(v)_2/2π R(1 - n_3) - n_2 n^(v)_2/4π(1-n_3)^2.


The direct correlation function of electrostatic interactions is

    c^(1)_el(z) =  -βϕ_el(ρ̅(z)) - 3β/4 R_w^3∫_-∞^∞ dz^'ρ(z^')  [ μ_el(ρ̅(z^'))/ρ̅(z^') - ϕ_el(ρ̅(z^'))/ρ̅(z^')] (R_w^2 - (z - z^')^2).




§ APPENDIX IV: DIRECT CORRELATION FUNCTIONS

In this Appendix, we consider the calculations of the second direct correlation function c^(2)(r,r^') within the nonlocal density functional theory, considered in the main text. 
The second direct correlation function c^(2)(𝐫 ,𝐫^') is the sum of the following contributions:

    c^(2)_el (𝐫 ,𝐫^')    =   -   2β(μ_el(ρ̅(𝐫^'))-ϕ_el(ρ̅(𝐫^')))/ρ̅(𝐫^')ω_el(|𝐫 - 𝐫^'|) 
       -   ∫_V d𝐫^''ρ(𝐫^'') [1/ρ̅( 𝐫^'')∂βμ_el(ρ̅(𝐫^''))/∂ρ̅(𝐫^'')  - 2β(μ_el(ρ̅(𝐫^'')) - ϕ_el(ρ̅(𝐫^''))/ρ̅^2(𝐫^'')]
       ×   ω_el(|𝐫 - 𝐫^''|) ω_el(|𝐫^' - 𝐫^''|),


    c^(2)_WCA (𝐫 ,𝐫^') = -β V_WCA(|𝐫 - 𝐫^'|),


    c^(2)_spc (𝐫 ,𝐫^') = -β V_spc(|𝐫 - 𝐫^'|),


    c^(2)_fmt (𝐫 ,𝐫^') = -β∑_α∑_γ∫_V  d𝐫^''∂^2 Φ({n})/∂ n_α(𝐫^'')∂ n_γ(𝐫^'')ω_α(𝐫^'' - 𝐫)ω_γ(𝐫^'' - 𝐫^').



     
    

The second direct correlation function for the homogeneous water in the Fourier-representation can be written as a sum of several contributions

    ρ c^(2)(k) = ρ (c_py^(2)(k) +  c_spc(k) + c_wca(k) + c_el(k)),

where

    ρ c_spc(k) = 24 ηϵ_swβ(r_sw/d_HS)^3 sin(kr_sw) - kr_rwcos(kr_sw)/(kr_sw)^3,


    ρ c_WCA(k)    =    24 ηϵβ(r_m/d_HS)^3 sin(kr_m) - kr_mcos(kr_m)/(kr_m)^3
       -    96βϵησ/k d_HS^3∫_r_m^r_c dr sin(kr)[ (σ/r)^11 - (σ/r)^5],

and

    ρ c_el(k) = -βρω_el^2(k)∂μ_el/∂ρ - 2βρω_el(k)(1 - ω_el(k))(μ_el - ϕ_el/ρ).

c_py^(2)(k) is the Fourier image of the direct correlation function of the hard spheres system in the Perkus-Yevick approximation (see, for instance, <cit.>), r_rw = σ_rw/2 and the Fourier-image of smoothing function ω_el is:

    ω_el(k) = -3 kR_w cos(kR_w) - sin(kR_w)/(kR_w)^3.





