{'fraction_non_alphanumeric': 0.04181889487481165, 'fraction_numerical': 0.008999525254401718, 'mean_word_length': 4.668421668421669, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Application of a two-length scale field theory to the solvation of charged molecules: I. Hydrophobic effect revisited

On a basis of a two-length scale description of hydrophobic interactions we develop a continuous self-consistent theory of solute-water interactions which allows to determine a hydrophobic layer of a solute molecules of any geometry with explicit account of solvent structure described by its correlation function. We compute the mean solvent density profile n(r) surrounding the spherical solute molecule as well as its solvation free energy. We compare the two-length scale theory to the numerical data of Monte-Carlo simulations found in the literature and discuss the possibility of a self-consistent adjustment of the free parameters of the theory. In the frameworks of the discussed approach we compute also the solvation free energies of alkane molecules and the free energy of interaction of two spheres separated by some distance. We describe the general setting of a self-consistent account of electrostatic interactions in the frameworks of the model where the water is considered not as a continuous media, but as a gas of dipoles. We analyze the limiting cases where the proposed theory coincides with the electrostatics of a continuous media.

On a basis of a two-length scale description of hydrophobic interactions we develop a continuous self-consistent theory of solute-water interactions which allows to determine a hydrophobic layer of a solute molecules of any geometry with explicit account of solvent structure described by its correlation function. We compute the mean solvent density profile n(r) surrounding the spherical solute molecule as well as its solvation free energy ∆G(r). We compare the two-length scale theory to the numerical data of Monte-Carlo simulations found in the literature and discuss the possibility of a self-consistent adjustment of the free parameters of the theory. In the frameworks of the discussed approach we compute also the solvation free energies of alkane molecules and the free energy of interaction ∆G(D) of two spheres of radius R separated by the distance D. We describe the general setting of the self-consistent account of electrostatic interactions in the frameworks of the model where the water is considered not as a continuous media, but as a gas of dipoles. We analyze the limiting cases where the proposed theory coincides with the electrostatics of a continuous media.

I. INTRODUCTION
We consider the current work as a step towards the development of a continuous thermodynamic theory of solute-solvent interactions. The construction of such a kind of a theory is of extreme importance for biochemical applications connected with the reliable and fast determination of the protein-ligand binding constants in solvent (water) with the precision comparable to the accuracy of explicit molecular-dynamical simulations.
There are few kinds of approaches which take into account the effect of fluctuating media on interactions between solvated molecules. Tentatively we can split them into the following three groups: i) "Explicit approaches". Generally such methods are based on numerical simulations taking into account motion of many water molecules. Such approaches are very time consuming because they demand huge computational resources.
ii) "Implicit approaches". These: a) either exploit the heuristic concepts of the "surface area accessible by the solvent" and of the "width of the hydrophobic layer", or b) make some effective renormalization of "bare" (i.e. vacuum) parameters of interacting molecules: in the frameworks of these approaches the effect of water on thermodynamic properties of dissolved substances is computed on the basis of some additional data on solubility, evaporation etcsee, for example, [1,2,3].
iii) "Intermediate approaches". The works belonging to that group (including our one), offer attempts which, on one hand, get rid of explicit numerical simulation of water molecules, and, on the other hand, preserve an effective account for water. The development of these methods on the basis of a clear understanding of a physical origin of solute-water interactions could, at least, give a better control for the parameters used in implicit approaches.
Most commonly, in the "implicit approaches" it is assumed [4,5] that the free energy of solvation is proportional to the solvent accessible surface area of a solute molecule. It is obviously true for solutes much larger than a water molecule. However, for the solutes of sizes comparable to the size of water molecules, as well as for the solutes of complex geometry, the question of the determination of the corresponding hydrophobic layer, being of great importance, is often resolved by means of purely empiric conjectures. For example, it is supposed usually that the hydrophobic layer is equal to the surface covered by a center of a water molecule "rolling" around the solute surface. Some modifications of the method deal with the so-called "molecular surface area" -see, for example, [6]. Another method of implicit account for the solvent in the so-called "Gaussian approximation" has been proposed by [7].
Among the theoretical attempts to develop a constructive theory of solvent-solute interactions, the following groups of works (to our point of view) deserve special attention: 1. The detailed computation of the solvent density-density correlation function in the presence of dissolved substances (with or without electrostatics): a) either introducing weighted densities and explicitly designing the form of the density functional which is entirely dictated by physical constraints [8], or b) by the construction of an appropriate "bridge functional" (i.e. closure relation) for the correlation function satisfying the hierarchy of integral equations [9]. On the basis of the correlation function all the thermodynamic quantities can be easily determined.
2. The computation of the averaged profile of the solvent density described by the fluctuating field in presence of solute molecules (with or without electrostatic interactions). The solvent structure on small scales is taken into account by its bulk correlation function, while the effect of solute solvation is considered by forcing the total solvent density to be zero inside the solute molecule (see, for example, [10] and references therein). Minimization of the corresponding free energy allows one to compute the density and the solvation free energy.
3. The "Scaled Particle Theory"(SPT) [11,12] and the "Information Theory" (IT), developed in the series of papers [13]. Both these approaches, SPT and IT, exploit the relation between the probability to find a cavity of a given shape containing some number of particles, with the solvation free energy of a solute molecule of that shape. The probability of a cavity formation in an ensemble of fluctuating particles (water molecules) is computed on the basis of semi-empirical combinatorial arguments.
Our solvation model borrows its ideas in the second group from the list above, namely from the series of works [14,15,16] where a two-length scale description of hydrophobic interactions was developed. In this approach the solvent density is decomposed in two components, the slowly varying field describing the mean solvent (water) density, and the short-ranged density fluctuations. We believe that such an approach is optimal from different points of view: on one hand it is physically clear, being "semi-microscopic", and on other hand it can be used as a constructive computational tool of account for water, much faster than the corresponding explicit approaches but without the loss of the precision. So, our aim is to construct a theory of solute-water interactions which would allow to determine the hydrophobic layer of solute molecules of any geometry with account of water structure described by its correlation function. The parameters of such a theory will be determined by testing them for solutes with simple geometries and comparing these cases with the experimental data. Briefly, the main ingredients of our work are as follows: i) We consider water as a continuous inhomogeneous media and describe it by a fluctuating continuous density field; the discreteness of the water structure on small scales is taken into account through the fluctuations controlled by the water correlation function χ in the bulk; ii) The solvation free energy of the substance situated in water is obtained in the selfconsistent mean-field approximation based on slightly modified Ginzburg-Landau-Chandler functional [14,15].
Our work is devoted to the calculation of the solvation free energy and of the free energy of interactions of solute molecules of any arbitrary shapes.
Specifically, we start with a two-scale Ginzburg-Landau Hamiltonian and perform an analytical integration over the short-ranged fluctuations. This provides a single-variable Hamiltonian yielding an integro-differential equation. The solution of this equation ensures the equilibrium mean density profile of the water, and the corresponding value of the Hamiltonian gives the solvation free energy. We pay a special attention to the physical clarity and consistency of all steps of derivation of our model. We discuss the hidden dangers on the way and compare our approach to the Monte-Carlo numerical simulations carried out in [10]. We determine in a self-consistent way the width and the geometry of a hydrophobic layer surrounding the solute molecules. So, we believe that a very important auxiliary goal of our work concerns the prediction of the relevant number of water shells (the "solvent accessible surface area") around a solute molecule.
Summarizing the said above, the advantages of the fluctuational approach with respect to other phenomenological theories of "implicit water" are as follows: a) There is no difference in computations of the solvation free energy and of the free energy of interactions of solute molecules; b) The detailed structure of the solvent is taken into account on small scales; c) We do not need any special determination of the "hydrophobic layer" by "rolling" the water molecule around the solute surface. The width and the geometry of a hydrophobic layer is computed on the basis of: a) the true geometry of a solute molecule, and b) the supposition about the solute-selvent and solvent-solvent interactions.
d) The theory contains small number of adjustable parameters, which in turn can be determined from the independent experiments (either numerical, or real).
Our paper is organized as follows. In the Section II we describe the model and derive the basic equations for the equilibrium and the full densities, and for the solvation free energy, as well as we solve numerically the obtained equations, compare them to the numerical Monte-Carlo simulations of the work [10] and propose the way of tuning the adjustable parameters of the Hamiltonian; in the Secion III we compute the solvation free energy of alkanes, as well as the free energy of interaction of two spherical solute molecules; the results of the work are summarized in the Section IV, where we also derive the equations for electrostatics in fluctuating dipolar environment.

II. THE MODEL
In this Section we present the self-consistent description of solvation which takes into account simultaneously: i) the smoothly varying average density field n(r), and ii) the water structure on small scales described by the density-density correlation function χ(r, r ′ ).
Following the general scheme of the works [16], we describe the solvent by two fields: n(r) -the average density varying on long distances, and ω(r) -the small-scale fluctuating field controlled by the correlation function of the water χ(r, r ′ ) = ω(r)ω(r ′ ) in the bulk. The solvent cannot penetrate into the volume occupied by the solute molecule hence the total solvent density ρ(r) = n(r) + ω(r) is nullified inside the solute: ρ(r) = 0 for all r ∈ v in , where v in is the volume occupied by the solute molecule. The fact that we force the total solvent density ρ(r) = n(r) + ω(r) to be zero inside the solute, results in an effective interactions between n(r) and ω(r). In general, one can permit also the direct coupling between n(r) and ω(r) everywhere in the solute.
The thermodynamic properties of the solvent (water) in the bulk are obtained using the self-consistent Ginzburg-Landau (GL) free energy functional which describes the liquidvapor phase transition. The average solvent density n(r) plays the role of an order parameter. The interactions between long-and small-scale fields n(r) and ω(r) can enforce the "drying", shifting the liquid-vapor transition. Minimizing the corresponding free energy functional, we compute the desired density profile and the solvation free energy.
It should be noted that the fully microscopic description of the solute-solvent interactions certainly does not require introducing of any scale separation-everything is described by the single microscopic field ρ(r). However in such a description the interactions between water molecules should also be described on the microscopic level. Such kind of theory has been developed in the works [8]. It seems to us that due to the heavy machinery, the application of such approach to the solvation and interactions of objects of complex architecture in water is yet rather restrictive for practical purposes. So, the nature of two different scales in the approach chosen in our work, is the consequence of the mean-field (GL) description of the water-water interactions.
Before going into details let us point out what is different in our consideration with respect to the original two-length scale theory of Chandler et al. [10,14,15,16]. i) We do not divide the space in "cells" for the large-scale field n(r) (as in [10]). We describe the system by two fully continuous densities. The equilibrium density is determined selfconsistently without additional averaging over small-scale fluctuations. We believe that this would result in better description of solvation of objects of rough shape and of the free energy of interactions between solutes located close to each other.
ii) We can tune the potential W (n(r)). Currently W (n(r)) describes the behavior of the system exactly at the vapor-liquid phase transition, but our consideration can be straightforwardly generalized to other situations.
iii) We can consider the charged objects in the frameworks of the "discrete electrostatics", where the density n(r) becomes the density of water dipoles. The set of basic equations for the density field n(r) and the electric field φ(r) is obtained in Section III, while in the forthcoming paper [17] we numerically solve the derived system of equations. iv) We can incorporate in the theory the direct Van-der-Waals solute-solvent potentials and consider the smooth solute-solvent boundary (see the Section III for more details).

A. The Hamiltonian
We describe the solvent in absence of any solute molecules by the following two-length scale Hamiltonian where n(r) is the smoothly changing (average) solvent density (also considered as the "order parameter" of the theory); ω(r) is the field corresponding to the short-ranged density fluctuations; χ(|r − r ′ |) is the solvent correlation function in the bulk; a is the phenomenological parameter which requires the microscopic determination-see the discussions below; and the self-consistent potential W (n(r)) is chosen in the common form of the standard Ginzburg-Landau (GL) expansion for the order parameter n(r) as the fourth-order polynomial: where n 1 and n 2 are the values of the order parameter n in the vapor and liquid phases correspondingly (below we set, if not specified, n 1 = 0) and b is the coupling constant which in combination with the parameter a defines the surface tension (see (10) below).
The Hamiltonian (1) consists of two (still decoupled) parts. The first term describes the non-local small-scale Gaussian fluctuations of the solvent. These fluctuations extend to the distances controlled by the solvent density correlation function in the bulk χ(|r − r ′ |), The second term describes the fluctuations of the self-consistent potential (2) in the system which could exhibit a phase transition. In our case this is a liquid-vapor (i.e. "drying") phase transition.
The choice of a functional form of the potential (2) is rather arbitrary being to some extent a question of a taste. However it should obey some basic requirements which are imposed mainly by physical reasons. As it has been mentioned already, the Hamiltonian (1) with the potential (2) should describe the vapor-liquid (i.e. "drying") phase transitions. This should be valid only in the vicinity of the solvated object, thus manifesting the nature of a hydrophobic effect. The presence of a drying transition near the solute surface shifts effectively the parameters of the bulk media to the region of the liquid-vapor coexistence [15]. The general form of the Hamiltonian describing the liquid-water coexistence is well known -it is the two-well potential usually parameterized by the GL expansion-see Fig.1. The parameterization being rather arbitrary, might be taken not only in the polynomial form as in the GL theory, but also in any form that correctly describes the region of phase coexistence under conditions close to normal. We believe that the GL expansion is well enough to reflect the major features of the physical picture. The choice of a more sophisticated form of the potential (2) entering in the large-scale part of the Hamiltonian (1) is the possible way of the refinement of the model in future within the frameworks of the mean-field approach.
In Fig.1 for generality we show broader class of potentials W (n) than described by Eq.(2), allowing the nonsymmetric shapes (see Fig.1b). The profiles in Fig.1a,b are drawn for n 1 = 0. The potential W (n) in the GL expansion depends upon four parameters (a, b, n 1 , n 2 ), two of them (the densities n 1 and n 2 ) are fixed in each phase. The parameters a and b are considered as the free adjustable coefficients of the theory. Supposing slightly more general form of the potential 2,W (n) = W (n) + µ(n − n 1 ), shown in Fig.1b, we may consider the relative difference (RD) between the energy densities in each phase, i.e. the difference between the values at the minima ofW (n) (see Fig.1b) as an extra adjustable parameter. It is easy to see that this difference is closely connected to the asymptotic behavior of the solvation energy at large solute sizes. Actually, the RD contributes to the volume-dependent part of the solvation energy. It is obvious that the term which scales linearly with the volume, tends to 1 3 P V solute as the solute volume V solute increases [11] (here P is the liquid pressure). However the simple quantitative estimation shows that on the nanometer scale the volume contribution is negligible in comparison with the surface contribution to the solvation energy [10]. Thus we may neglect the difference in minima of the potential shown in Fig.1b and consider everywhere further the potential symmetric upon the phase change-see Fig.1a. Just such functional form of the potential is proposed above in (2).
We can find the equilibrium densities minimizing the Hamiltonian (1) with respect to the (still decoupled) fields n(r) and ρ(r): In absence of any interactions between the fields ω and n, the solution of (4a) postulates the homogeneity of the liquid ω(r) = const, while Eq.(4b) describes: a) the plain vapor-liquid interface supplied with the boundary conditions n(−∞) ≡ n 1 = 0 (n(∞) = n 2 ) in the onedimensional case, and b) the density profile near the surface of extended macroscopic object in the two-and three-dimensional case with boundary conditions n = 0 at the solute surface and ∇n = 0 for |r| → ∞. Everywhere further in this Section for simplicity we consider the spherical objects only (the solvation of objects of other shapes is discussed at lengths of the Section III).
Qualitatively these solutions can be analyzed by computing the first integral of motion (as in the classical mechanics [18]). Writing the first term in Eq.(4b) in 3D as ∆n(r) = n ′′ (r) + 2 r n ′ (r), multiplying (4b) by n ′ (r) ≡ dn(r) dr and integrating it from r in (the radius of the sphere) to ∞, we obtain: Hence, we get from (5) the following relation: Having the equation (6) and the exact functional form of the potential (2) in hand, we know the value of n ′ (r) at the point r = r in which determines the qualitative behavior of the density profile. In our case the second integrand vanishes, hence the density monotonically increases from n 1 at the point r in to n 2 for |r| = ∞. In practice, the liquid density reaches quickly the bulk value at some distance outside the solute surface. The analysis of an asymmetric potentials with RD can be performed within the same formalism.
The sense of the parameters a and b becomes transparent in the limit of large spherical solute molecules r in ≫ δ, where δ is the characteristic size of the liquid-vapor interface (see below). Substituting (2) into (4b) and replacing the interaction of the field n(r) with the solute molecule by the boundary condition, we get where n avr = n 2 + n 1 2 . Supposing in (7) that r ≫ δ, we can expand the 3D Laplacian and obtain the effective 1D equation whose solution reads It is clear that the symmetric potential of GL form leads to the one-parameter theory. This parameter δ fixes the width of the phase separation interface. The solvation free energyG per the spherical solute area 4πr 2 in in the limit of large molecules can now be expressed as followsG where the coefficient has a sense of a surface tension. Note that for 0 < δ r in ≪ 1 in the leading order of the expansion (10) the true boundary condition n(r in − δ) = 0 cannot be distinguished from the approximate one, n(r in ) = n avr .

B. The correlation function
The important component of our model is the bulk correlation function, χ(|r − r ′ |), of a solvent. It appears as an input into the theory. Recall some basic facts concerning the function χ(|r − r ′ |). We have by definition: where n is the mean solvent density in the bulk, h(r − r ′ ) = g(r − r ′ ) − 1 and the function g(|r − r ′ |) has very clear physical meaning: The function g(r) may be determined from experimental data, or numerical simulations, or, say, from some selfconsistent integral equation with corresponding bridge functional [19].
The following relation is the direct consequence of the definition of χ(r − r ′ ): The inverse correlation function Let us pay attention that the integration in (14) is spread only to the volume v in occupied by the solute and hence the function χ −1 in (r, r ′ ) is not translationally invariant. We have used in our work two different correlation functions: i) the correlation function of water obtained from the numerical molecular-dynamical simulations [20], and ii) the correlation function of hard spheres constructed on the basis of solution of self-consistent integral equation completed by the Percus-Yevic bridge functional [22,23]. The functions h(|r−r ′ |) = g(|r−r ′ |)−1 of water and of hard sphere liquid are shown in Fig.2. Both functions correspond to the dimensionless bulk density n 2 = 0.7. Some approximate equations of the theory of liquids in the statistical thermodynamics of classical liquid systems one can find, for instance, in the review article [19].
Working simultaneously with two different correlation functions (of water and of hard spheres) we able able to check the sensitivity of the physical observable quantities (as the The correlation function defines the natural scale, ξ, in the theory. In what follows we fix ξ = 2.78Å as the distance till the first peak in the correlation function of water.

C. The partition function
Now we turn to the solute-solvent interactions. To take into account the effect of the solute solvation, we demand the total density ρ to be zero inside the volume v in of the solute: where r ∈ v in . The direct interaction between the short-and the long-range density fluctuations in the solvent can be written in the simplest possible way where c is the coupling constant. In what follows all the analytic derivations are performed for generic c, however the numerical solutions of the corresponding equations are presented for c = 0 only. The case c = 0 we leave for the future consideration.
The solute-solvent partition function is obtained as follows where Let us rewrite the constraints imposed by the Dirac δ-function in (17) using the functional Fourier transform where N is the number of points in the product where Ξ is some physically irrelevant constant which defines the normalization of the partition function Z.
The partition function Z with the action S[ω, n, ψ] defined in (20) is the key object for our computations. Such quantities as the total, ρ(r), and the averaged, n(r), densities, the solvation free energy, ∆G, the free energy of interactions of different solute molecules in the solvent, etc. can be calculated on the basis of the generating function Z.

D. The mean density
To define the mean (large scale) density n(r), let us require In (21) all other fields except n(r) are supposed to be already integrated over. This leads to the functional equation which can be written explicitly as a set of two coupled integrodifferential equations relating the inner (r ∈ v in ) and outer (r ∈ v out ) regions of the solute molecule: and −a∆n(r) + δW (n(r)) In the current paper we pay attention to the case c = 0 only, i.e. when the direct coupling between the short-and long-ranged fluctuation is absent. In that case the profiles of the mean density n(r) for spherical solute molecules of different sizes are given by the solutions of the simplified integro-differential equation We solve (23a)-(23b) numerically by using the simple technical trick which allows to increase essentially the speed of the computations. Namely, we have observed that the most CPU time is spent for the calculation of χ −1 in (r, r ′ ) on the basis of (14). To get rid of this part of computations we can multiply the (23a) by χ(|r − r ′ |) and integrate over the whole space v tot (v tot = v in ∪ v out ), with the subsequent application of (14). Finally we arrive at the following set of equations: The generic form of the mean density profile n(r) for the spherical solute molecule is shown in Fig.3 for different (still arbitrary) values of the constants a and b of Hamiltonian. The densities of the vapor (n 1 ) and bulk (n 2 ) phases are as follows: n 1 = 0; n 2 = 0.7. The numerical solutions of (24) are obtained with the correlation functions of hard spheres (dashed line) and of water (solid line). The typical plots are shown in Fig.3.
The presence of the regions with the negative values of the mean densities in the vapor phase should not confuse-remember that only the full density ρ(r) is a conserved value and should be zero inside the solute molecule, while the average density n(r) is determined selfconsistently. The existence of the regions with negative values of the long-ranged averaged field n(r) signifies that in these regions the fluctuations of the short-ranged field ρ(r) are not symmetric, such that locally ω = 0. Such regions exist only inside the solute volume v in and in the "physical" volume v out (i.e. outside the solute) the mean density is always positive. The full density ρ(r) = n(r) + ω(r) which minimizes the functional (20) satisfies the following set of equation The equations (25) with c = 0 again can be rewritten in the form which excludes the function χ −1 in (r, r ′ ) from the corresponding expressions.

E. The solvation free energy
The free energy of solvation ∆G is defined as the energy necessary to transport the solute molecule from its environment in the solvent to the vacuum. The partition functions of solvent samples of sufficiently large volume V solvent containing the solute molecule inside can be straightforwardly written on the basis of (17): where H is given by (18a). Conventionally the solvation free energy ∆G is written as follows where Z solvent only is the partition function of the pure solvent (without the solute molecules). At first glance it seems naturally to write Z solvent only simply as However the expression (28) being used in (27) leads to the divergence of the solvation free energy ∆G in practical computations based on the approach developed in (17)- (24). The formal reason for such a divergence deals with the occurrence of uncompensated infinite product of Gaussian integrals (20) in the ratio Z solvent−solute /Z solvent only in (27). The physical origin of such a divergence is due to the mixture of different statistical ensembles associated with the partition functions (26) and (28). Namely, when writing Z solvent−solute as in (26) and imposing the δ-function constraint, we fix some particular value of the field ρ within the solute volume v in what means that (26) is the partition function of the canonical ensemble (with respect to the density ρ). At the same time, for Z solvent only written in the form (28), we allow any density of the field ρ inside the solute volume. Hence, (28) is the partition function of the grand canonical ensemble.
The regularization of eq.(28) is based on the probabilistic consideration of the solvation free energy [15]. Rewrite ∆G in (27) as follows where Z v (N) is a partition function of N molecules inside the volume v in . The continuous analog of the last expression appears when passing from N to the average number of moleculesñ = N/v in . Simultaneously we require the net density ρ to be equalñ inside v in and replace the summation over N by the integration overñ. Now we can rewrite the normalization partition function Z solvent only in the following "regularized" form Instead of this form we use an approximation obtained from (30) by the stationary phase integral calculation technic. Thus, the definition (27) with the partition functions given by (26) and (31) is very natural, justified physically and does not contain any divergences. It is important to notice that the computation of the equilibrium density based on Z solvent only given by (30) is the bulk density, for which the value of corresponding effective Hamiltonian is strictly zero.
The solvation free energy ∆G for any solute molecule is computed on the basis of the equilibrium density profile n eq (r) ≡ n(r) calculated from Eqs.(22a)-(22b). Substituting n eq (r) back into the Hamiltonian (18a) and taking into account (14), we can rewrite the solvation free energy G ≡ H{n eq } as a sum of three terms: where ∆G corr = 1 2 v in n eq (r)χ −1 in (r, r ′ )n eq (r ′ ) drdr ′ ∆G surf = a 2 (∇n eq (r)) 2 dr ∆G int = W (n eq (r)) dr (33) As in the computations of the mean density n(r), here we consider the case c = 0 only. In that case (32) can be rewritten in the following form The contributions to the normalized solvation free energy of a spherical solute molecule of radius R are shown in Fig.4 for each term {∆G corr , ∆G surf , ∆G int } = {∆G corr , ∆G surf , ∆G int }/(4πR 2 in ) separately, as well as for the sum ∆G solv = ∆G corr + ∆G surf + ∆G int . The diameter of the solvent molecule in Fig.4 is set to 1. The solvation free energy of a spherical solute molecule in absence of electrostatic interactions [14,15,16] is proportional (as expected) to the volume of solute molecule for sufficiently small sizes R in (of order of the correlation length in the solvent) and tends to be proportional to the surface area for large R in . We clearly see the non-monotonic behavior of the solvation energy upon the size of a solute molecule [39] reported in some works (see, for example, [10]).

F. Adjustment of the parameters of the Hamiltonian
The free parameters a and b of the Hamiltonian (1)-(2) are chosen from two basic requirements: i) The theoretically computed normalized solvation free energy ∆G solv = ∆G solv /(4πR 2 ) of a spherical solute molecule of a radius R (see (35)) reproduces the corresponding dependence ∆G solv (R) obtained in the Monte-Carlo simulations of [10]; ii) The theoretically computed full density ρ(r|R) for few sizes R of spherical solute molecules (25) reproduces the behavior g(r) of the correlation function in vicinity of the spherical solute of radius R extracted from the Monte-Carlo simulations [10].
The corresponding results are shown in Fig.5 for the following choice of the parameters: a = 9; b = 30. The dashed and solid lines represent the results obtained with the correlation functions of hard spheres and of the water correspondingly. The solvation energy ∆G solv is represented in Fig.5 in dimensionless units ∆G solv /γ, where γ is the surface tension of the flat vapor-water interface. The numerical value γ ≈ 72 mJ/m 2 is taken from [10].
The solvation energies ∆G solv and ∆G solv computed directly using (34) and (35) These transfer coefficients remain unchanged in all further computations. In Fig.4 and Fig.5 the radius of a solute molecule is measured in Angstroms.

III. RESULTS AND DISCUSSIONS: SOLVATION OF OBJECTS OF VARIOUS SHAPES
The current mean-field-type theory of solvation is applied to computation of solvation of neutral molecules (alkanes), as well as to the computation of the free energy of interactions of separated objects (spheres).

A. Solvation of alkanes
The parameters of the theory and the scaling coefficient are adjusted to satisfy the Monte-Carlo simulations of solvation of hard spheres and the corresponding profile of the solvent density near the solvated object (see the previous Section) and are not tuned anymore. In particular, we use the same parameters to compute analytically the solvation free energy of alkanes. As it is shown below, we find very good coincidence of our computed values with the predictions of the "Scaled Particle Theory" which uses the parameters tuned especially to alkanes [25].
The mentioned coincidence needs some elucidations. There is a viewpoint [27,28,29] supported by numerical computations that one can split the interaction between the solute and the solvent into two parts: i) the free energy, ∆G cav , of a "cavity formation", and ii) the dispersion (attractive) part of the Van-der-Waals interactions, ∆G disp . The numerical methods involving the "thermodynamic integrations" allow to compute both the contributions, ∆G cav and ∆G disp separately [29]. These contributions can be also computed separately in the frameworks of the "Scaled Particle Theory" (SPT) developed in [11,12]. The corresponding description of alkanes has been undertaken in [25] where the authors have reported very good quantitative agreement of the sum ∆G cav + ∆G disp with the experimentally measured solvation free energy.
The model described in the previous sections of our work considers the solute as an object bounded by hard walls and hence takes into account only the "cavitation" part of the solvation free energy. The attractive part of the Van-der-Waals interactions is not yet taken into account. However we do not see any principal obstacles in adding the dispersion part of solute-solvent interactions to our model. Moreover, we can "smear" the δ-functional constraint in (19), releasing more realistic form of the Van-der-Waals potential. The corresponding computations are in progress and will be reported in a forthcoming publication [17].
The comparison of our predictions for the cavitation part of the solvation free energy of alkanes with the corresponding contribution ∆G cav extracted from the papers [25] and [26] is shown in the table below.
Let us remind once more that in our model the parameters of the Hamiltonian a = 9, b = 30 and the scaling factors (36) are tuned on the basis of the Monte-Carlo simulations of solvation of hard spheres without any additional adjustment to alkanes. The volume v in occupied by the alkane molecules we compute using the following prescription. We take the coordinates of centers of all atoms from a standard database, then we draw the spheres with the effective radii 3.5Å for all carbon atoms (the length of the C-O bond in water) and 3.05Å for all hydrogen bonds (the length of the O-H bond in water) and find v in as alkanes ∆G cav , kcal/mol (our model) ∆G cav , kcal/mol (from [25]) ∆G cav , kcal/mol (from [   , [30]. The same parameters can be extracted using the force field developed in [31]. As it has already been mentioned above, the dispersion part of the Van-der-Waals (VdW) interaction can be easily taken into account. In the presence of the attractive contribution to the VdW interactions, we should add an extra term U disp to the mean-field potential (1), where U disp = U(r) n(r)dr and U(r) is the attractive part of the VdW interactions.
In the presence of the potential U(r) the equations (23a)-(23b) for the equilibrium mean density are modified in the following simple way: (23a) reminds without changes, while the equation (23b) acquires an additional contribution from the attractive part of the potential U(r) in the region r ∈ v out : −a∆n(r) + δW (n(r)) δn(r) + U(r) = 0 (37)

B. Solvation of cylinders and the free energy of interactions of two separated spheres
In this section we consider the cavitation part of the solvation free energy of molecules of cylindrical geometry, as well as of the free energy of interaction of two hard spheres in the fluctuating media (i.e. in water). Schematically the studied systems are shown in Fig.6. The purpose of this section is two-fold.
1. First of all, we compute the depth of the "skin layer" of the solvated object. The naive computation of the solvation free energy (as well as of the free energy of interactions) of large molecules like ligands and proteins demands essential computational resources. Actually, according to (24) one has to take into account the interactions of the fluctuating density field n(r) with all points inside the solvated object v in . However it is clear that due to the hard wall constraint the fluctuating field n(r) cannot penetrate deep inside the body of a solute. Considering the solvation of cylindric molecules of different widths and lengths, we study the penetration depth of the density field n(r) inside the solute molecule.
The dependence of the solvation free energy ∆G cyl (L|R) upon its length L for two different radii R is shown in Fig.7a. As one sees, from Fig.8b, the slope of the curves in Fig.7 becomes constant approximately for R > R 0 ≈ 1.4ξ = 3.9Å. It means that for R > R 0 the fluctuating density field penetrates inside the body of the solute to some fixed length only (to the depth of the "skin layer"). Hence, we may not account of all internal points of the solute located into its body beyond the skin layer R 0 ≈ 3.9Å.
2. Secondary, we investigate the free energy of interaction ∆G sph = ∆G sph int (D|R) − ∆G sph solv (D → ∞|R) of two spherical solute molecules, each of radius R, upon the distance D between their centers. The free energy of interaction, or "the mean force potential" [32,33] ∆G sph (D|R), with the subtracted solvation energies of two solitary spheres, is computed for the hard-core potential, U(D), where i.e. we again take into account only the cavitation part of the free energy of interactions, neglecting the dispersion (attractive) part of the direct Van-der-Waals interactions between the impenetrable spheres of radius R.
Even for such simplified situation we have found that the behavior of ∆G sph (D|R) on D is qualitatively different for different values of R compared to ξ [40]. Let us stress once more that ∆G sph is the free energy of interaction between two spheres with the subtracted solvation energies of two solitary spheres ∆G sph solv (D → ∞|R). The observed dependencies are shown in Fig.8. For all radii, R = 0.4ξÅ , 0.7ξÅ and 0.9ξÅ we reproduce the non-monotonic dependence ∆G sph (D) on D. This non-monotonicity is the manifestation of the oscillatory behavior of the water correlation function. Such behavior has been found earlier in direct numerical simulations consistent with the predictions of the "Scaled Particle Theory" (see [32] and the references therein), and on the basis of estimation of the solvent accessible surface area [33]. The account of the Van-der-Waals attractive tail in the potential U(D) can essentially change the behavior ∆G sph (D), however the non-monotonic behavior ∆G sph (D) will still hold.
Comparing Fig.8a and Fig.8c, we see that the dependence ∆G sph (D) is sensitive (even qualitatively) to the size R of interacting spheres. The effective attraction between smaller spheres is weaker than that of larger spheres. This effect can be easily understood. For sufficiently large spheres, R ξ, there is an penalty in the free energy of keeping two solute molecules close to each other, which is the "classical manifestation" of the Casimir-type effect. At the same time, the sufficiently small solute molecules (with R ξ) are located in spontaneously created fluctuational cavities of size ξ in water and the penalty for expelling the water from the volume occupied by two solute molecules close to each other is essentially smaller than that for large solutes. Moreover, one can see from Fig.8a that the mean force potential has the attractive (negative) part corresponding to the second water shell only. So, the two sufficiently small spherical solutes can form a bound state separated by one water shell. The found behavior qualitatively coincides with the results of the work [32], where the authors have studied the "potential of the mean force" acting between two spherical solute molecules.

C. Electrostatics in fluctuating dipolar environment
There are different ways to implement the electrostatic interactions into our selfconsistent mean-field description. In many works it is assumed that the dielectric permittivity of the solvent, ε, linearly depends on the solvent density, n(r), i.e., ε(r) = 1 + α n(r). In our approach we get rid of such supposition and follow the self-consistent scheme where the solvent is considered as a "gas" of dipoles. Our method ideologically is very close to the consideration of the Debye screening in electrically neutral plasma [34]. However, to the contrast with plasma, in our case the positive and negative ions are connected by holonomic constraints in pairs, forming the short-ranged dipoles. The screening in the system of extended objects was a subject of many investigations. The description, most appropriate for our goals, has been developed in the appendix of the paper [35].
Below we derive the basic set of equations describing simultaneously electrostatic and hydrophobic interactions in the fluctuating media. The detailed analysis of the obtained equations together with the numerical computation of the free energies of solvation and interactions of charged molecules, will be the subject of the forthcoming publication [17], while here we analyze briefly the limiting case, in which the proposed theory coincides with the electrostatics in the continuous media with some effective density-dependent permittivity ε(r).
Let us introduce the quantities: ρ q mol (r, τ ) -the charge distribution of the solvent molecule with the orientation τ and with the center of mass situated at the point r. Outside of the molecule ρ q mol (r, τ ) = 0; ρ q solute (r) -the charge distribution of the solute; φ eff (r) -the effective mean electrostatic field.
Sticking to the mean-field description developed in [35], we can place each solvent molecule in a self-consistent field electrostatic field, φ(r), averaged over the volume of this solvent molecule. So, we have: The electrostatic potential φ eff (r) is determined by the Poisson equation where µ(r) is the media charge density and ρ(r) = n(r) + ω(r) is the solute density. Neglecting the short-range fluctuations, ω(r) we may write in the mean-field approximation The full Hamiltonian of the system reads now: H[n(r), φ eff (r)] = dr a 2 (∇n(r)) 2 + W (n) + 1 2 v in dr dr ′ n(r)χ −1 in (r, r ′ )n(r ′ ) + λ dr (∇φ eff (r)) 2 + dr ϕ(r) ∆φ eff (r) + 4π ε solute ρ q solute (r) + µ(r) (42) Let us stress that according to [36] we should minimize the functional action of hydrophobic interactions together with the one of electrostatic field, while the Poisson equation (39)) is considered as a constraint and hence enters in the action with the functional Lagrange multiplier ϕ(r).
Minimizing the functional H[n(r), φ eff (r)] with respect to the fields ϕ(r), φ eff (r), n(r), we come to the system of self-consistent Euler equations, fixing the corresponding equilibrium distributions.
To simplify the current consideration, let us regard the case of point dipoles. The molecular density distribution ρ q mol in this case is given by where d is the dipole moment of the molecule. Then φ(r, d) = − dr ′ φ eff (r ′ ) d∇ r ′ δ(r ′ − r) = d∇φ eff (r) this equation together with the linearized version of (40) leads to the following expression for the media charge density µ(r): Inserting the last relation into the Hamiltonian (42) and minimizing the corresponding expression, we get This system should be provided with some physically justified boundary conditions, for example ∂ϕ(r) ∂r → 0 for |r| → ∞; φ eff (r) → 0 for |r| → ∞; n(r) →n for |r| → ∞ (47) The numerical computation of the solvation and interaction free energies of charged extended objects in the dipolar fluctuating media (in the water) based on the set of derived equations (46) is in progress and will be the subject of our forthcoming publication [17].
It is easy to check that in the limit of the weak electrostatic field equations (46) lead to the standard electrostatics in the continuous media with an effective dielectric densitydependent permittivity. Actually, when the field ϕ eff is small, we may assume what means that ∇ϕ ∇φ eff ∼ φ 2 eff (49) We have omitted in (48)-(49) all numerical prefactors, supposing for simplicity that all quantities are dimensionless.
Let us stress that our derivation implies the absence of the correlations between the orientational degrees of freedom of water molecules in the bulk and near the solute surfaceonly in this case we can pre-average the density µ(r) in (41) over the orientational degrees of freedom. We support our supposition by the conclusion made in [37,38] on the basis of extensive numerical simulations: the rearrangement of the orientational degrees of freedom of water molecules near the solute surfaces plays the secondary role in the thermodynamics of solvation (see also [29] for more discussions).
To conclude with, let us repeat that in the current work the most attention is paid to the consideration of the pure hydrophobic (cavitation) effect. In particular, the main achievements are as follows.
1. The method described in the paper provides the constructive basis for the consistent description of solvation effects of hard solute molecules of any shape. One sees that the continuous two-length scale fluctuational approach with only two free parameters tuned to the Monte-Carlo results on solvation of hard spheres, manages to describe the known cavitation contributions to the solvation free energy of alkanes without any additional adjustment of these parameters.
2. In the framework of the same approach we have considered two auxiliary problems: a) the computation of the solvation free energy of solutes of cylindric shape, and b) the computation of the free energy of interactions of two spheres separated by some distance. The consideration of the problem (a) allows one to conclude that the depth of the "skin layer", i.e. the penetration length of the fluctuating field into the body of the solute is of order of 1.4ξÅ, where ξ is the characteristic length scale of the theory (the location of the first maximum in the bulk correlation function of the solvent) while the investigation of the problem (b) permits us to make conjectures about the potential of the mean force acting between hard spheres in the fluctuating media.
3. The developed approach can be extended to take into account the dispersion part of the solute-solvent interactions, as well as the electrostatic contributions of charged solute molecules. We describe the basic steps towards the implementation of both these contributions into the current approach.