Langmuir 1997, 13, 5459-5464
5459
The Origins of Hydration Forces: Monte Carlo Simulations and Density Functional Theory Jan Forsman,*,† Clifford E. Woodward,‡ and Bo Jo¨nsson† Physical Chemistry 2, University of Lund, Chemical Center, P.O. Box 124, S-221 00 Lund, Sweden, and School of Chemistry, University College, University of New South Wales, Australian Defence Force Academy, Canberra ACT 2600, Australia Received April 14, 1997. In Final Form: July 21, 1997X The force between two infinite planar surfaces, interacting with the intervening solvent via a shortrange exponentially decaying attraction and a soft repulsion, has been calculated using Monte Carlo simulations and density functional theory. In the former case, the intervening liquid is chosen to mimic water, whereas the density functional calculations has been applied to a nonpolar liquid. It was found that the most important feature of the surface-solvent potential, in terms of how it governs the solvation force, was its range. When the decay length was larger than about half a molecular diameter, a strong repulsion was found at short separations (1 (20) ze1
while n j (z) is a coarse grained density, given by a weighted density average, taken over the excluded volume, σ3, centered on z
n j (z) )
∫
3 4σ3
R(z,h)
L(z)
dz′n(z′)[σ2 - (z - z′)2]
(21)
where L(z) ) max(z - σ, 0) and R(z,h) ) min(z + σ, h). The bulk pressure and the chemical potential are related to the bulk density, nB, via the simple gvdW equation of state
PB )
([
µ ) -kBT ln
kBTnB 3
1 - nBσ
]
-
16πn2Bσ3 9
(22)
)
32πnBσ3 1 - nBσ3 1 (23) + nB 9 1 - nBσ3
The equilibrium density profile, neq(z), minimizes Gs and thus (26) Nordholm, S.; Haymet, A. D. J. Aust. J. Chem. 1980, 33, 2013. (27) Nordholm, S.; Johnson, M.; Freasier, B. C. Aust. J. Chem. 1980, 33, 2139. (28) Freasier, B. C.; Nordholm, S. J. Chem. Phys. 1983, 79, 4431.
5462 Langmuir, Vol. 13, No. 20, 1997
Forsman et al.
satisfies the following integral equation,
n(z)
(
∫
3 ) exp - 3 1 - n(z)σ 4σ 3
∫
h
0
R(z,h)
L(z)
n(z′)[σ2 - (z - z′)2] dz′ -11-n j (z′)σ3
)
dz′n(z′)φ(|z - z′|) + Φ(z,h) - µ kBT
(24)
which was solved iteratively, until the maximum absolute value of the functional derivative, δGs/δn(z), was less than 5 × 10-6. We used an integration grid of 0.05 in the z-coordinate. Test calculations were performed using smaller grids, and essentially identical results were then obtained. The free energy minimizations runs very fast, and a calculation typically takes a few seconds on an ordinary work station.
Results Simulations. We begin by describing the values we used for the surface potential parameters in eq 17. The soft repulsion parameter A was fixed to prevent the oxygens from approaching closer than about 1 Å to the surfaces. We chose the following parameters for the orienting surface: C ) 36 kJ/mol and τC ) 1.6 Å. This potential has a range and magnitude to significantly order only those solvent molecules adjacent to the surface. Any orientational ordering in other layers occurs mainly via solvent-solvent interactions. For the nonorienting surfaces, C ) 0. We considered several values for the parameters, B and τB, associated with the spatially dependent part of the surface potential. These were chosen about some reference values, Bref and τB,ref. Reasonable reference values were determined as follows. We used 1.6 Å as the reference value for the decay length τB, which was the same as τC. For orienting surfaces, we required that the choice of Bref and value of P| at the starting separation should give reasonable properties for the system at (essentially) infinite separation. The parameters B and C determine the hydrophilicity of the surface which is manifest in the surface tension of the wall-solvent interface. The latter is half of the value of Gs at infinite separation. In our work, the inherent surface energy of the walls is not included in Gs; thus, the limiting value is given by -2γw cos θ, where γw is the air-solvent surface tension (∼70 mJ/m2 for water at 1 atm), and θ is the surface-solvent contact angle, which is less than 90° for a hydrophilic surface. For contact angles less than 0°, the solvent completely wets the surface. The limiting value of Gs is also determined by the bulk pressure, i.e., the value of P| at infinite separation. The value of P| at the starting separation should be such that the liquid solvent does not cavitate as the separation increases. That is, the solvent in the bulk is liquid. Further, we also required the bulk to have properties not too far removed from those of water under STP conditions. We chose to set P| ) 500 bar at a separation of 12 Å and Bref ) -14 kJ/mol. These values ensured that the liquid did not cavitate and gave a reasonable limiting value of Gs of ∼-50 mJ/m2. It should be noted that in previous work on hydrophobic surfaces21 we estimated the surface tension of SPC water at an inert surface to be ∼80 mJ/m2 at ∼1500 bar. Assuming a similar value for the air-solvent surface tension of our bulk fluid, gives a contact angle for the reference surface potential of ∼72°. Other values for B and τB were obtained by transforming parameters at a separation of 12 Å while constant chemical potential was maintained. The latter condition was met by suitably adjusting P| and ensured that the bulk solvent properties remained the same, thus providing us with meaningful comparisons. Hysteresis was checked for by
Figure 1. Free energy of interaction per unit area as a function of the separation, under various water-surface potential conditions, but at constant bulk water conditions: solid, Bref, Cref, τB,ref, τC,ref; circles, Bref, Cref, 1.5τB,ref, τC,ref; diamonds, 1.5Bref, 0, τB,ref, 0; squares, 1.5Bref, Cref, τB,ref, τC,ref.
performing “closed loop” transformations of the potential parameters and was found to be negligible. Given a set of parameters, FEDM simulations were performed in order to generate curves of P| versus h using the methods outlined in the previous section. In order to obtain the free energy of interaction between the surfaces we needed to estimate the bulk pressure; see eq 15. This was done by considering the component of the pressure tensor normal to the surfaces
Ps )
d(P|h) dh
(25)
This quantity appears to approach the bulk pressure more rapidly than does P|. For all the potential parameters, Ps is already quite close to the bulk pressure at separations beyond 12 Å (see Figure 3). We estimated the bulk pressure as the average value of Ps in the interval 13-15 Å and obtained 122 bar, with a standard deviation of ∼13 bars. Given the bulk pressure, the free energy Gs as a function of separation could then be obtained from eq 15. Figure 1 shows the free energies for a selection of surface potential parameters. Varying the parameters of the surface potential has a profound influence the limiting value of the free energy. Clearly, the surface becomes more hydrophilic given increases in the range and the strength of the attractive and orientational components of the surface potential. Using the simple estimate outlined above, the most attractive surface has a contact angle very close to zero. All the free energy curves are oscillatory at short separation, reflecting structuring of the solvent between the surfaces. For the shorter ranged potentials (τB ) 1.6 Å), the interaction between orienting surfaces does not appear to be significantly more repulsive than nonorienting surfaces. Indeed, at this value for τB, none of the surface interactions appear to have any significant repulsion underlying the oscillatory component and thus are inconsistent with experimental measurements. However, in the case of the longest ranged surface potential, τB ) 2.4 Å, there is a clear underlying repulsion. This is due to the direct removal of solvent particles from the energetically favorable regions of the surface potential as the surfaces are squeezed together. Clearly, the range of the surface potential will then be mirrored in the range of the surface interaction. We also considered values for the surface potential parameters that maintained the chemical potential while
Origins of Hydration Forces
Figure 2. Free energy of interaction per unit area as a function of the separation, under various water-surface potential conditions, but at constant bulk water conditions and similar bulk interfacial tensions: solid, Bref, Cref, τB,ref, τC,ref; circles, 0.55Bref, Cref, 1.5τB,ref, τC,ref; diamonds, 1.3Bref, 0, τB,ref, 0; stars, 0.76Bref, 0, 1.5τB,ref, 0.
keeping the hydrophilicity of the surfaces approximately constant. This was accomplished by changing two parameters simultaneously to approximately maintain the limiting value of the surface free energy. For example, we coupled an increase of the decay length, τB, with a decrease in B. The free energy curves corresponding to these parameters are given in Figure 2. The limiting values of the free energy appear to be within 25 mJ/m2 of each other. As observed in Figure 1, the net interaction free energy becomes more repulsive when τB is increased to 2.4 Å. A significant repulsion is obtained for both orienting and nonorienting surfaces. This repulsion appears larger when the surface induces orientational order in the solvent. Even the free energy curves for the shorter ranged surface potentials, τB ) 1.6 Å, show that orientational ordering also leads to a more repulsive net interaction. This is despite the fact that the nonorienting surface appears to be slightly more hydrophilic. Recall that the comparison between the orienting and nonorienting surfaces in Figure 1 did not seem to indicate any significant increase in repulsion for orienting surfaces. In that case, however, the strength of the nonorienting surface potential was constant and relatively large. In the present comparisons, the parameters have been adjusted to approximately maintain surface hydrophilicity. This suggests that hydrophilicty induced by an orienting surface potential leads to a more repulsive surface interaction than if one uses a nonorienting attractive potential. This must be due to the effect that orientational order has on solvent-solvent interactions. The interaction between orientationally ordered dipoles has a longer range than that between (correlated) rotating dipoles. This notwithstanding, it would appear that the most significant influence on the surface interaction is obtained by increasing the range of the surface potential. Certainly, there is no noticeable underlying repulsion between the surfaces at separations greater than 5 Å unless the range of the surface potential is sufficiently large. We anticipate that similar results would have been obtained had we chosen to increase the range of the orienting potential, τC, instead. In Figure 3 we present the corresponding net solvation pressures, Ps(h). The data obtained using twice as many molecules, with the reference system parameters, are also provided. We see that any size-dependent effects were small.
Langmuir, Vol. 13, No. 20, 1997 5463
Figure 3. Net solvation pressure, Ps, as a function of the separation. The notation is identical to the one in Figure 2, with the addition of the results, indicated by plus signs, obtained for the reference system using twice as many molecules (in the separation interval 5 Å < h < 12 Å).
Figure 4. Free energy per unit area as a function of separation. Simulation data are denoted in the same way as in Figure 2, whereas the dashed curves show the results from gvdW calculations with τB ) 1.6, 2.4, and 3.2 Å, giving increasingly stronger repulsion.
The repulsion observed to accompany the increase in range of the nonorienting surface potential should not only occur with solvent molecules containing dipoles. That is, the mechanism responsible for the repulsion is not due to any special orientational structuring in water, for example, “ice cages” or “hydrogen-bonded networks”, as has sometimes been proposed. Thus, similar behavior should be observed for a simple nonpolar solvent between “solvophilic” walls. In order to test this assertion, we performed density functional calculations on a confined liquid composed of nonpolar particles, interacting with each other via a Lennard-Jones potential. The Lennard-Jones parameter σ was set equal to 3.1 Å, making σ3 approximately equal to the mean molecular volume in water, at STP conditions. The values of the amplitudes in the surface-solvent and solvent-solvent interactions were chosen to give the same bulk pressure, limiting interfacial tension and density profiles similar to those in the simulations reported above. Qualitatively, the results are not very sensitive to these parameters. The free energy curves, for different decay lengths, τB, are given in Figure 4. A strong solvation repulsion is also found in the nonpolar system, provided τB is sufficiently large. Comparing these with the simulation results for corresponding values of τB and a
5464 Langmuir, Vol. 13, No. 20, 1997
Figure 5. Comparison of results using fine-grained gvdw calculations, coarse-grained gvdW calculations, and calculations by assuming zero compressibility, respectively, in order of increasing limiting interfacial tension: solid, τB ) 1.6 Å; diamonds, τB ) 3.2 Å.
nonorienting potential, there is almost quantitative agreement, though the curves are shifted slightly with a trivial bulk free energy constant, due to the difficulty of obtaining exactly the same limiting interfacial tension in the simulations. In order to more clearly see the qualitative changes obtained as the range of the surface potential is varied, one can use more approximate approaches, which will reduce or remove the oscillations in the free energy curves. For instance, in the coarse-grained version of the gvdW theory, entropy contributions are entirely local; i.e., n(z) is replaced with n j (z) in eq 19: thus packing effects are not accounted for. An even simpler approach is to set the compressibility of the solvent equal to zero in the gvdW calculations, which means that the solvent density between the surfaces will be constant and equal to its bulk value. The results using these approaches are presented in Figure 5. By using the coarse-grained gvdW version, we obtain free energy curves that should be qualitatively more similar to those expected for a solvent confined between bilayers, where the “softness” of the surfaces will reduce density oscillations. Note that the equilibrium separation grows rapidly as the decay length increases. At sufficiently large separations, the attractive long ranged van der Waals contribution due to the truncation of solvent-solvent interaction caused by the surfaces will always dominate. Conclusions The solvent contribution to the interaction between hydrophilic surfaces will always be repulsive. That is, there is always a free energy increase when such surfaces are brought from infinite separation to contact and thus are completely dehydrated. Another way to view this is that the free energy of interaction, Gs, calculated in this study, will be zero at contact (h ) 0) and have a negative value at infinite separation for hydrophilic surfaces. The question of interest is, what is the range of this interaction and what is it due to? Our results indicate that the range of the interaction between inert hydrophilic surfaces is
Forsman et al.
affected mainly by the range of the surface potential. We use the term “inert” to mean that there are no mobile surface groups. The increase in range of the hydration force is seen in both orienting and nonorienting surfaces. This effect is observable in nonpolar solvents, as well, as it does not rely on any orientational ordering in the solvent. Indeed, such repulsive interactions in nonaqueous solvents have been found experimentally.5,9,29 Bergenståhl and Stenius9 suggested that only systems where the solvent molecules are able to form hydrogen bonds with hydrophilic sites at the surfaces will show repulsive solvation forces. Such hydrogen bonds may provide the necessary strength of interaction with the surface in order to produce repulsive hydration forces. However, our study does not support the conclusion that hydrogen bonds between the solvent molecules explains the range of the interaction. Our observations suggest that the range of repulsive hydration forces measured when water is a solvent is due to the surface potential. Other effects that we have not considered in this study may prove viable candidates for the mechanism behind repulsive forces, for example, the steric interaction between mobile surface groups.10 Thus far, we have not discussed the direct surfacesurface interaction. The aim of this article was to consider the solvent-mediated contribution to surface interactions only, in particular, to consider its dependence on aspects of the surface potential. However, some comments may be in order. It is possible that if the surface-solvent interaction is increased, the direct attraction between the surfaces will grow concomitantly and may even cancel the enhanced solvation repulsion30 asymptotically. However, at close separations, linear response theory (which underpins the asymptotic analysis) may not be valid for the interactions between hydrophilic groups on the different surfaces. Higher order effects may well act to saturate the direct interaction between surfaces at short separation, leading to a weaker attraction than predicted from an asymptotic analysis. Furthermore, surfaces that contain fixed hydrophobic groups, which cannot correlate, may indeed repel one another asymptotically, depending on their initial alignment. Thus, the sign and magnitude of the total interaction between surfaces is intrinsically tied to the way the surface groups themselves are modeled. The general problem requires a more extensive model that explicitly accounts for hydrophobic surface groups. It is only in this way that one can treat surface-solvent and surface-surface interactions consistently. We reiterate that such a model is beyond the scope of the present work. Of course the relevance of our conclusions to reality depends on the reliability of the SPC model. This model posesses only three partial charges and therefore cannot form three-dimensional tetrahedral structures. In this context, it should be noted that, despite this shortcoming, we did find an increased surface repulsion for orienting surfaces. It may be worthwhile, therefore, to carry out similar investigations on other models for water including, perhaps, an explicit treatment of surface groups. LA9703813 (29) Persson, P. K. T.; Bergenståhl, B. A. Biophys. J. 1985, 47, 743. (30) Attard, P.; Kjellander, R.; Mitchell, D. J.; Jo¨nsson, B. J. Chem. Phys. 1988, 89, 1664.
