Nonlocal density functional theory of water taking into account many-body dipole correlations: Binodal and surface tension of 'liquid-vapour' interface

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.

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.

I. 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) 1-6 , modelling of water purification from toxic compounds (e.g., ions of heavy metals 7 ), description of mechanical stability of construction porous materials (concrete, wood, paper, etc.) during capillary condensation/evaporation cycles [8][9][10][11] , 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 [12][13][14][15][16][17][18][19][20][21] , none of them satisfy the requirements formulated above.
In paper 12 the authors formulated a density functional theory (DFT) taking account of the molecular structure of water within the TIP4P model 22 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 sitesthrough 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 23 . 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 24,25 . Nevertheless, from the result of the manuscript 26 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 27 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 [13][14][15][16][17][18] . In paper 13 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 14 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 15 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 15 , this discrepancy is most probably due to the simplicity of the gradient correction used. In paper 16 a self-consistent DFT, taking explicit account of hydrogen bonding through the Statistical Associating Fluid Theory (SAFT) 28 , 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 17,18 , 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 29 , 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 29 . In the work 30 , 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 31,32 . 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.

II. 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, 33,34 ). 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 20,21,[35][36][37][38] . 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 39,40 . 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: where is the free energy of the ideal gas. The excess free energy can be written in the following where 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 where r c is the cutoff of LJ pairwise potential of interactions; Θ(r) is the Heviside stepfunction; σ 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 41,42 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 43 . The first term in (4) describes the contribution of the excluded volume interactions of the hard spheres within the Percus-Yewick approximation; η = πd 3 BH ρ/6 is the packing fraction of hard spheres with the effective Barker-Henderson diameter, determined by the following Pade approximation 44 The second term in (4) describes the total contribution of the attractive interactions, thus B can be defined as: 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): where l is the dipole length and y = ρp 2 /3 0 k B T , p is the dipole moment. The auxiliary function σ(y) = is also introduced; the compressibility factor α of the hard spheres (see Appendix 1) in Percus-Yevick approximation has the following form 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 45 , 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 46 : where k = |k|; the pair correlation function can be expressed through structure factor as:

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 where V is system volume, F id [ρ(r)] is the free energy of the ideal gas, F ex [ρ(r)] 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 46 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 and the excess free energy, in turn, consists of several parts: where with the Helmholtz free energy F f mt [ρ(r)] of hard spheres with the effective BH diameter d BH within the fundamental measure theory (FMT) 50 , namely where the free energy density is which depends on six weighted densities: where 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 46,51,52 .
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 53 . Thus, following the idea of the weighted density functional theory 51,52 , we treat the electrostatic free energy functional as follows: where φ el (ρ(r)) = f el (ρ)/ρ is the free energy per fluid particle, depending on smoothed density, which, in turn, is determined in the following way: where ω el (|r|) = 3/(4πR 3 w )Θ(R w − |r|) 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 (13), i.e. from the Euler-Lagrange equation where is the one-particle direct correlation function; β = (k B T ) −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 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.

IV. NUMERICAL RESULTS AND DISCUSSIONS
Now we will consider the application of the model to the description of water liquidvapour 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 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 kP a) 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:   be clearly identified by the comparison of the decaying length in the following equation 46 :  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. 8 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 57 .

V. 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 microand 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.

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 where β = (k B T ) −1 is the inverse thermal energy, F 0 [{ρ α }] is the free energy functional of the reference system without Coulomb interactions between the particles; N is the normalization constant which will be specified below; 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 where δρ α (r) = ρ α (r) −ρ α (r) is the fluctuation of the local charge densities; is the inverse structure operator of the reference system for which we adopt the following approximation where c αγ (r − r ) is the matrix of direct correlation functions of the sites that are not bonded to each other and 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 where is the inverse structure operator in the random phase approximation. Further, choosing the normalization constant N 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 where is the electrostatic free energy and is the screening function. Note that we have subtracted from the final expression the electrostatic self-energy of the molecules G αγ (k) are the Fourier-images of the structure factors of the reference system which can be calculated from the following matrix relation 58 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 where 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 where is the characteristic function corresponding to the probability distribution function g(r).
The screening function 24,25,40,58,59 is 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) (50) and 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 where is the correlation function of the hard spheres 46 .
Further, using the following model characteristic function the approximation and taking the integral (41), we arrive at the following relation for the density of the electrostatic free energy where l is the dipole length and the auxiliary functions and are introduced. Here, α = 1 − Z 0 = 1 − ρk B T χ 0 with the isothermal compressibility χ 0 of the hard spheres and y = q 2 l 2 ρ/3ε 0 k B T . 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 In the Percus-Yewick approximation taking into account the equation of state for the reference system and the relation for the compressibility χ 0 = ρ −1 ∂ρ/∂P 0 , one can easily obtain a relation for the "compressibility factor" where η = πd 3 ρ/6 is the packing fraction of the hard spheres; d is the hard sphere diameter; λ is the thermal de Broglie wavelength.

BULK PHASE
The total pressure P and chemical potential µ are the sum of four contributions: and where P id = ρk B T and µ id = k B T ln(ρλ 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, and The electrostatic contribution to the chemical potential takes the following form with the auxiliary functions The electrostatic contribution to the pressure can be calculated by the relation

VIII. 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 where g W CA (z, z ) is Above we used the following notations r m = 2 1/6 σ and ψ(x) is The three independent weighted functions are: where e z 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: where the derivatives are determined as follows The direct correlation function of electrostatic interactions is