Water at Hydrophobic Substrates: Curvature, Pressure, and

and Roland R. Netz. Langmuir 0 (proofing), ... and Colin Nuckolls. Journal of the American Chemical Society 0 (proofing), ..... A.V. Nguyen. Advan...
0 downloads 0 Views 168KB Size
Download PDF
4756

Langmuir 2004, 20, 4756-4763

Water at Hydrophobic Substrates: Curvature, Pressure, and Temperature Effects Shavkat I. Mamatkulov and Pulat K. Khabibullaev Heat Physics Department of Uzbekistan Academy of Sciences, 28 Katartalstr., 700135 Tashkent, Uzbekistan

Roland R. Netz* Sektion Physik and Center for Nanoscience, Theresienstr. 37, Ludwig-Maximilians-University Munich, 80333 Munich, Germany Received October 29, 2003. In Final Form: January 27, 2004 We studied the water density profile close to spherical and planar hydrophobic objects using molecular dynamics (MD) simulations. For normal pressure and room temperature, the depletion layer thickness of a planar substrate is ∼2.5 Å. Even for quite large spherical solutes with a radius of R ) 18 Å, the depletion layer thickness is reduced by 30%, which shows that substrate curvature and roughness is an experimentally important factor. Rising temperature leads to a substantial increase of the depletion layer thickness. The compressibility of the depletion layer is found to be surprisingly small and only ∼5 times higher than that of bulk water. A high electrostatic surface potential of 0.5 V is found, which presumably plays an important role in the presence of charged solutes, since it can promote adsorption into the interfacial layer.

I. Introduction The reaction of water to the presence of a hydrophobic object depends crucially on the size (i.e., the radius of curvature for the spherical shape) of the solute: For small nonpolar objects (with a radius of no more than a few angstroms), the hydration is characterized by the clathratelike geometry of the water shell surrounding the solute, with water molecules oriented such that hydrogen bonding groups (donors and acceptors) are parallel to the solute surface.1,2 The water density close to the hydrophobic surface is increased relative to the bulk density in this case. For very large nonpolar objects or in the limiting case of a planar hydrophobic substrate in contact with water, on the other hand, the water density is reduced at the hydrophobic surface.3,4 This crossover of the hydration properties as a function of solute size has been extensively studied in the past using analytical theory5,6 and computer simulations.7-9 As a consequence of this subtle crossover, the effective interaction between two solutes embedded in liquid water depends on the solute size, being (at short distances) repulsive between small hydrophobic objects10 and attractive between large hydrophobic objects.11-14 These effects constitute the driving force behind the self(1) Pratt, L. R.; Chandler, D. J. Chem. Phys. 1977, 67, 3683. (2) Rapaport, D. C.; Scheraga, H. A. J. Phys. Chem. 1982, 86, 873. (3) Stillinger, F. H. J. Solution Chem. 1973, 2, 141. (4) Lee, C. Y.; McCammon, J. A.; Rossky, P. J. J. Chem. Phys. 1984, 80, 4448. (5) Lum, K.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570. (6) Huang, D. M.; Geissler, P. L.; Chandler, D. J. Phys. Chem. B 2001, 105, 6704. (7) Wallqvist, A.; Berne, B. J. J. Phys. Chem. 1995, 99, 2885. (8) Hummer, G.; Garde, S. Phys. Rev. Lett. 1998, 80, 4193. (9) Wallqvist, A.; Gallicchio, E.; Levy, R. M. J. Phys. Chem. B 2001, 105, 6745. (10) Pangali, C.; Rao, M.; Berne, B. J. J. Chem. Phys. 1979, 71, 2975. (11) Wallqvist, A.; Berne, B. J. J. Phys. Chem. 1995, 99, 2893. (12) Spohr, E. J. Chem. Phys. 1997, 106, 388. (13) Koga, K. J. Chem. Phys. 2002, 116, 10882. (14) Pertsin, A. J.; Hayashi, T.; Grunze, M. J. Phys. Chem. B 2002, 106, 12274.

assembly of organic substances in water, such as protein folding, enzymatic reactions, micellization, etc. Beyond that, it seems fair to say that, without a proper characterization of the behavior of water at hydrophobic surfaces, no true understanding of the properties of such surfaces and their interactions will be possible. Most likely, many of the features interpreted as being inherent to surfaces themselves might in fact reflect properties of the interfacial water layer instead (e.g., protein adsorption resistivity,15 zeta potentials,16 surface potentials,17 and polymeradsorption energies,18 just to mention a few). In this paper, we primarily concentrate on the water density profiles close to hydrophobic curved and planar surfaces using molecular dynamics (MD) simulations. Our work is motivated by recent scattering experiments where the water density depletion at planar nonpolar substrates was determined. As the main result of those experiments, it was shown that the effective depletion thickness (defined as the thickness of a steplike depletion layer consisting of a vacuum with the same integrated depleted amount as the rounded and smeared-out depletion profiles found in experiments) is ∼2-3 Å on hydrophobic polystyrene substrates19 and ∼5 Å on hydrophobic self-assembled monolayers20 using neutron reflectivity methods. In X-ray reflectivity measurements, a value of ∼1 Å on paraffin substrates was found.21 Finally, ellipsometry measure(15) Harder, P.; Grunze, M.; Dahint, R.; Whitesides, G. M.; Laibinis, P. E. J. Phys. Chem. B 1998, 102, 426. (16) Chan, Y.-H. M.; Schweiss, R.; Werner, C.; Grunze, M. Langmuir 2003, 19, 7380. (17) Kreuzer, H. J.; Wang, R. L. C.; Grunze, M. J. Am. Chem. Soc. 2003, 125, 8384. (18) Friedsam, C.; Del Campo Becares, A.; Jonas, U.; Seitz, M.; Gaub, H. E. New. J. Phys. 2004, 6, 9. (19) Steitz, R.; Gutberlet, T.; Hauss, T.; Klo¨sgen, B.; Krastev, R.; Schemmel, S.; Simonson, A. C.; Findenegg, G. H. Langmuir 2003, 19, 2409. (20) Schwendel, D.; Hayashi, T.; Dahint, R.; Pertsin, A. J.; Grunze, M.; Steitz, R.; Schreiber, F. Langmuir 2003, 19, 2284. (21) Jensen, T. R.; Ostergaard Jensen, M.; Reitzel, N.; Balashev, K.; Peters, G. H.; Kjaer, K.; Bjornholm, T. Phys. Rev. Lett. 2003, 90, 086101.

10.1021/la036036x CCC: $27.50 © 2004 American Chemical Society Published on Web 04/27/2004

Water at Hydrophobic Substrates

ments point to a depletion thickness of ∼2-3 Å.22 The reason for the discrepancies among different experiments is not well understood, but with the aforementioned dependence of water hydration on the radius of curvature of the solutes, it is clear that surface roughness is one important factor (among many others, e.g., small traces of attractive interactions between wall and water molecules) and will, if present, reduce the depleted amount. Our simulations are done for two different geometries: First, we consider a single spherical solute with a varying radius inside a cubic box filled with water (modeled as SPC/E water) at constant pressure by using periodic boundary conditions. In the second part, we look at a planar hydrophobic substrate (formed by a layer of parallel alkane molecules oriented normal to the interface) in contact with a water slab. In contrast to previous numerical implementations of the same problem, we fix the lateral area of the interface and apply a pressure normal to the interface, which minimizes spurious effects due to the interfacial tension between the hydrophobic substrate and the water phase.23 We clearly observe the crossover from water adsorption at small nonpolar solutes (with a radius of 4 Å) to depletion for the largest spherical cavity considered (with a radius of 18 Å) and for the planar substrate. We consider two different definitions of the depletion thickness. The first is defined with respect to the density of the hydrophobic wall, which is the experimentally relevant number; here, we find a thickness of ∼2.5 Å in the planar case. The second is defined with respect to the effective potential felt by the water molecules (as will be explained in more detail later on), and this is a definition which was used in previous theories and many simulations; in this case, the depletion thickness is larger and of the order of 3 Å. A main problem in the analysis of the scattering experiments described above is the possibility of adsorption of gas or other organic contaminants at the hydrophobic substrate, which would lead to a scattering signature quite similar to the one caused by a depletion layer consisting of a vacuum. To furnish an independent way of characterizing the depletion layer, we studied the temperature and pressure dependence of the hydration at spherical cavities and planar substrates. Quite surprisingly, the depletion layer of planar substrates has a rather low compressibility and is only slightly deformed for pressures up to 1000 atm. Temperature increase on the other hand leads to an enhanced depletion layer, in agreement with previous simulations,21 which might be an experimentally viable route for distinguishing depletion layers from surface layers of contaminants. Finally, we also calculate the temperature dependence of the electrostatic potential profile across the interface. The top water layer is oriented such that the hydrogen atoms predominantly point away from the aqueous halfspace.21,24,25 This way, the image charge repulsion is minimized, as will be explained further below. The surface potential due to the orientation is of the order of 0.5 V and turns out to be almost independent of temperature. This water surface potential has consequences on the distribution of charged species (such as salt ions) close to the surface and leads to charging effects even on neutral substrates.26 (22) Granick, S. To be submitted for publication. (23) Marrink, S.-J.; Marcelja, S. Langmuir 2001, 17, 7929. (24) van Buuren, A. R.; Marrink, S.-J.; Berendsen, H. J. C. Colloids Surf., A 1995, 102, 143. (25) Brodskaya, E. N.; Zakharov, V. V.; Laaksonen, A. Langmuir 2001, 17, 4443. (26) Jungwirth, P.; Tobias, D. J. Phys. Chem. B 2001, 105, 10468.

Langmuir, Vol. 20, No. 11, 2004 4757 Table 1. Parameters for the Different Water Models Considereda ww model (kJ/mol) SPC/E TIP3P TIP4P

0.6500 0.6364 0.6480

σww (Å)

H-O-H bond angle (deg)

O-H distance (Å)

q (e)

3.166 00 -0.8476 1.0000 3.150 61 -0.8340 0.9572 3.153 65 -1.0400 0.9572 O-M: 0.1500

109.47 109.47 104.52

a Here, q is the partial charge on the oxygen atom, σ ww is the Lennard-Jones radius, and ww is the Lennard-Jones interaction parameter. In the TIP4P model, the oxygen charge is slightly displaced from the oxygen atom toward the hydrogen atoms by a distance of O-M.

Table 2. Buckingham Parameters for the Interaction between the Water Molecules and the Spherical Cavitiesa

aR

R (Å)

A

B (nm-1)

4 10 18

2.45 × 108 1.6 × 1015 1.154 × 1026

48 35 33.3

denotes the cavity radius.

II. Computer Simulation Details For the water-water interactions within our MD simulations, we used the standard Lennard-Jones (LJ) potential

[( ) ( ) ]

qRqβ σww Uww(r) ) ΣRβ + 4ww rRβ r

12

-

σww r

6

(1)

augmented by Coulombic interactions between all charged sites. The LJ interaction is centered on the oxygen atoms, with r denoting their mutual distance and ww and σww being oxygen-oxygen interaction parameters. As a test, we performed bulk simulations of three different standard water models, which will be discussed below, for which we list the respective parameters in Table 1. For all our simulations in the presence of hydrophobic solutes or interfaces, we chose the SPC/E model. For the interactions between the spherical hydrophobic objects and the water molecules, we used the purely repulsive part of the Buckingham potential

Usw(r)/kBT ) A exp(-Br)

(2)

where A and B are the Buckingham potential parameters, r is the solute water separation, and kBT is the thermal energy. The radius, R, of the spherical solutes is defined as the solute-water distance at which the repulsive interaction is of the order of the thermal energy, that is, Usw(R)/kBT ) 1. Buckingham potential parameters for the different cavities are presented in Table 2. The hydrophobic sphere radius was varied from 4 to 10 to 18 Å. The simulation cells were cubic boxes containing 1024, 2013, and 3315 water molecules, respectively, and we employed periodic boundary conditions using the NPT (constant particle number, pressure, temperature) ensemble. For the case of the planar hydrophobic substrate, we used an inhomogeneous ensemble; namely, we fixed the lateral area (in x-y directions), while the box size in the z direction was controlled by a constant pressure. The simulation cell thus was a rectangular prism, and again, we employed periodic boundary conditions in all three directions. At the beginning of our simulation, the cell dimensions were Lx ) Ly ) 39.02 Å and Lz ) 84 Å, containing N ) 2781 SPC/E water molecules and a hydrophobic wall. The hydrophobic wall was built using

4758

Langmuir, Vol. 20, No. 11, 2004

Mamatkulov et al.

64 alkane molecules (CH3-(CH2)18-CH3) with an initial length of l ) 24.345 Å each. In the first phase of equilibration, the alkane molecules were placed along the z direction at the center of the simulation box with the terminal groups on both ends fixed on a square lattice. Then, all methylene groups were equilibrated for 20 ps while keeping the terminal groups facing one interface fixed and allowing the other terminal groups only to move along the z direction. For a second period of 20 ps, the roles of the two interfaces were exchanged. After that, all methylene groups were only allowed to move along the z direction. After thermalization over 100 ps at normal pressure and a temperature of T ) 300 K, the cell contained water with a density of 0.995 g/cm3. The density of the alkane wall was 0.82 g/cm3. The interactions involving methylene groups were also modeled by a LJ interaction potential, which in a slightly different notation is written as

Uij(r) )

C(12) ij r12

-

C(6) ij r6

(3)

where r stands for the separation between methylene groups themselves or between methylene groups and water molecules. The coefficients appearing in eq 3 are related to the LJ interaction parameters in eq 1 by C(n) ij ) 4ijσnij. To model the interaction between methylene groups, we chose the same values as those employed by Lee et al.,4 namely, ss ) 1.197 kJ/mol and σss ) 3.76 Å. For the interaction between methylene and water, the LJ parameters were determined by the combination rule C(n) ij (n) 1/2 1/2 ) (C(n) ) 0.882 kJ/mol. i Cj ) , which gives sw ) (sw) All production and parametrization simulations as well as simulation analyses were performed with the Gromacs package.27 The LINCS method and a leapfrog algorithm with a time step of 2 fs were employed.28 Temperature was kept constant by the Berendsen thermostat, and pressure coupling was done using Berendsen and Parrinello-Rahman methods.29,30 We used the particle mesh Ewald (PME) method for Coulombic interactions, which is known to reduce finite-size effects substantially for dipolar systems such as water. For all other (nonCoulombic) interactions, we introduced a cutoff. Since this cutoff also applies to the interaction between the cavity and the water molecules, it has to be chosen larger than the cavity radius; otherwise, the water molecules would not feel the repulsive potential sufficiently far away from the cavity center. To treat water-water interactions correctly, it should also not fall below a certain minimal value. For the simulations of a spherical solute immersed in water, the cutoff distance was always larger than the hydrophobic object radius and chosen as rc ) 9, 12 , and 20 Å for a cavity radius of R ) 4, 10, and 18 Å, respectively. For the planar case, the LJ interaction cutoff radius was chosen as rc ) 12 Å. Simulation results where we checked that our simulation results do not depend on the choice of the cutoff radius will be presented further below. To obtain a feeling for the differences between the various standard water models used in the literature, we (27) van der Spoel, D.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tielemann, D. P.; Sijbers, A. L. T. M.; Hess, B.; Feenstra, K. A.; van Drunen, R.; Lindahl, E.; Berendsen, H. J. C. Gromacs User Manual, version 3.1.1, University of Groningen: Groningen, The Netherlands, 2002; www.gromacs.org. (28) Allen, T. D.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (29) Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (30) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182-7190.

Figure 1. Temperature dependence of the bulk density of a few different water models compared with the experimental values. The simulation box contained 1024 water molecules in a periodic cubic box. The system was thermalized for 100 ps and averaged for 1 ns. Table 3. Bulk Water Density as a Function of Temperature and Pressure for the SPC/E Water Model At a Fixed Pressure of p ) 1 bar T (K)

F (g/cm3)

250 265 273 300 320 340 360

1.022 54 1.017 40 1.012 81 0.999 79 0.987 34 0.971 69 0.954 08

At a Fixed Temperature of T ) 300 K p (bar)

F (g/cm3)

1 500 1000 2000

0.999 79 1.001 86 1.023 42 1.059 99

first performed simulations of bulk water, using 1024 water molecules in a periodic cubic box. In Figure 1, the obtained water densities for three different water models are compared with the experimental values. As can be seen, the SPC/E model performs best for a wide temperature range; we therefore have chosen this model to produce all further results in this paper. Table 3 shows the resulting densities for the SPC/E model as a function of temperature (for fixed pressure) and as a function of pressure (for fixed temperature). An issue of much concern in such simulations is the possibly quite large relaxation time as one looks at water structures around large cavities or solutes. The same is true for water at planar hydrophobic substrates. To estimate the relaxation times in our simulation study, we followed two strategies: First, we performed simulations where we watched the relaxation of a global quantity, namely, the system volume, as one parameter is changed abruptly. The solid line in Figure 2a shows the behavior of the total volume of 3315 water molecules as the radius of the spherical solute is increased from R ) 15 Å to R ) 18 Å at time t ) 0; the simulations were performed at a fixed pressure of p ) 1 bar and at a fixed temperature of T ) 300 K. As one can see, the system equilibrates over a few picoseconds. The broken line in Figure 2a shows the volume as the temperature is abruptly increased from T ) 300 K to T ) 340 K at time t ) 0 with a fixed cavity radius of 18 Å. Here, one sees that the relaxation is such a fast process that it is hardly resolved. In Figure 2b, the

Water at Hydrophobic Substrates

Langmuir, Vol. 20, No. 11, 2004 4759

Figure 2. (a) Behavior of the total system volume when the radius, R, of a spherical solute is increased from 15 to 18 Å at time t ) 0 for a fixed temperature of 300 K (solid curve) and when the temperature increases from 300 to 340 K for a fixed cavity radius of 18 Å (broken curve). The simulation box contained 3315 SPC/E water molecules, and all simulations were performed at a constant pressure of p ) 1 bar. The relaxation time is a few picoseconds. (b) Autocorrelation function of the dipole moment in a spherical shell around a cavity with a radius of 18 Å at a temperature of T ) 300 K and a pressure of p ) 1 bar. The total simulation time was 100 ns. The relaxation time is a few picoseconds. (c) Water density profile around a cavity with a radius of 10 Å at a fixed temperature of T ) 300 K and pressure of p ) 1 bar for three different values of the LJ cutoff radius. The results agree within the numerical inaccuracy.

autocorrelation function of a local quantity, namely, the dipole moment inside a thin spherical shell around a cavity with a radius of R ) 18 Å, is shown as a function of time. The autocorrelation function falls off to zero over a few picoseconds, which means that there are no slow modes present which would hamper simulation of these systems. It follows that the relaxation is fast enough so that simulations over a few nanoseconds should suffice in order to obtain equilibrium properties. In Figure 2c, we address how the Lennard-Jones cutoff might influence the results of our simulations. As discussed before, we use periodic boundary conditions (implemented by particle mesh Ewald techniques) for the long-ranged Coulombic interactions but cut off all nonCoulombic interactions at a distance of rc. In Figure 2c, we show water density profiles around a cavity with a radius of 10 Å at a fixed temperature of T ) 300 K and pressure of p ) 1 bar for three different values of the LJ cutoff radius, as indicated in the figure. The results agree within the numerical inaccuracy, which shows that the cutoffs used are large enough that the results do not depend on the presence of this cutoff. We performed similar careful checks for the different cavity sizes used and for the planar geometry and got the same result. III. Hydration of Spherical Solutes In the first part, we look at the density distribution of water molecules around a single spherical solute. The interaction between the solute and the water molecules is given by eq 2 and thus is purely repulsive and exponentially decaying. It should therefore reflect the properties of an extremely hydrophobic substance. Denoting the number of water molecules in a spherical shell of volume ∆V ) 4πr2∆r by ∆N, where r is the distance from the center of the spherical solute and ∆r is the thickness of the shell, the density distribution (which is proportional to the radial distribution function) is given by F(r) ) ∆N/ ∆V. In Figure 3, we plot the radial density distribution normalized by the bulk density for different solute radii of R ) 4, 10, and 18 Å and for various temperatures and pressures. In the same plot, we also show the solutewater potential in units of kBT, as given by eq 2. As explained before, the radius of the solute sphere is defined by the distance where the rescaled interaction is unity. The density distribution as a function of the distance from the center of the excluded volume region directly illustrates the change of hydration as the radius increases. For the smallest sphere considered, that is, R ) 4 Å (see

Figure 3a and b), the density exibits oscillations manifesting the microscopic granularity of liquid water and is increased close to the sphere as compared to the bulk value. This is due to some highly ordered, clathratelike structure of water which engulfs the cavity. For R ) 18 Å, on the other hand (see Figure 3e and f), the density depletion is evident. In a spherical shell with a thickness of ∼2 Å around the solute, the water density is almost zero. As we will show in the next section, this behavior for the largest spherical solute considered by us is quite similar to the one at a planar hydrophobic substrate. The density profiles for the solute radius R ) 10 Å (see Figure 3c and d) are intermediate between density enhancement and depression and are close to the radius where the hydration behavior crosses over from adsorption to depletion. Our results for the temperature dependent hydration show that a temperature increase in general decreases the density close to the solute particle. This means that depletion effects are stronger at higher temperatures. Our results for the pressure dependence show that the depletion layer for the largest particle with R ) 18 Å is quite stable against applied pressure. In other words, the depletion layer is not very compressible, as will be discussed in more detail later on. IV. Hydration of the Hydrophobic Wall Although numerous studies of water-solid24,25 and water-vapor23,31,32 interfaces have shown that, for nonpolar purely repulsive interactions between the substrate and water, water is depleted from the interface, the precise length over which water is depleted is still a matter of debate. This is linked to the fact that the ensemble within which the interface is studied is crucial: Somehow the constant pressure property has to be put in, and thus, the system volume or the particle number has to fluctuate. The lateral system size is coupled to the interfacial tension (which is quite large), and thus, for a finite system, the effective lateral pressure is different from the externally imposed pressure. This has been explicitly demonstrated in a recent simulation.23 Grand-canonical simulations, where the system size stays constant and the constant pressure is produced by varying the particle number, suffer from low relaxation rates, since the probability for particle (31) Townsend, R. M.; Rice, S. A. J. Chem. Phys. 1991, 94, 2207. (32) Taylor, R. S.; Dang, L. X.; Garrett, B. C. J. Phys. Chem. 1996, 100, 11720.

4760

Langmuir, Vol. 20, No. 11, 2004

Mamatkulov et al.

Figure 3. Temperature dependence of the water radial distribution function for a spherical nonpolar solute with a radius of (a) 4, (c) 10, and (e) 18 Å at a fixed pressure of p ) 1 bar. Pressure dependence of the water radial distribution function for a solute with a radius of (b) 4, (d) 10, and (f) 18 Å at a fixed temperature of T ) 300 K. The simulation box contained 1024, 2013, and 3315 SPC/E water molecules for a radius of 4, 10, and 18 Å, respectively. The broken line denotes the interaction potential between the solute and water molecules, given by eq 2, rescaled by the thermal energy. Each system was thermalized for 100 ps and averaged for 2 ns.

Figure 4. Snapshot of the MD simulation of a planar hydrophobic slab (made up of 64 alkane molecules) in contact with a water slab consisting of 2781 SPC/E water molecules.

insertion is very small. To circumvent these problems, we fix the lateral system size, that is, the interfacial area, and only allow the system size in the z direction to fluctuate. We also built up the hydrophobic substrate by self-assembled alkane chains, which seems to be an acceptable representation of the substrate structure used in recent experiments.19-21 In Figure 4, a snapshot of the MD simulation is shown, which serves to illustrate the geometry of the system. As explained in section II, the alkane molecules form a compact slab in the middle of the simulation box. They are only allowed to fluctuate in the z direction and thus allow for fast pressure equilibration. The water slab has a thickness of ∼4 nm, which should be large enough such that bulk water properties are reproduced. We therefore interpret our results as being caused by the single hydrophobic substrate-water interface and neglect interactions between the two interfaces through the finite water slab. In Figure 5, we show the normalized densities of the alkane slab (dotted line) and the water layer (solid line)

Figure 5. Normalized density profile of the hydrophobic alkane slab (dotted line, in the middle) and the water layers (solid line, to the left and to the right) at a constant pressure of p ) 1 bar and temperature of T ) 300 K. The broken line denotes U h (z), the laterally averaged Lenard-Jones potential felt by the water molecules. The simulation box contained 64 hydrophobic molecules and 2781 SPC/E water molecules. The system was thermalized for 100 ps and averaged for 2 ns.

at atmospheric pressure and at a temperature of T ) 300 K. It is clearly seen that between the alkane slab and the water layer a region of reduced density appears. The density profiles are calculated using pointlike atomic form factors, and they denote the nuclear density; they therefore correspond to what would be seen in a neutron scattering experiment. The broken line denotes the laterally averaged interaction potential due to the alkane slab; in other words, this is the potential energy felt by the water molecules. The expression for this laterally averaged hydrocarbon slab-water interaction potential, U h (z), is obtained by integrating the Lennard-Jones pair potential, Usw(r),

Water at Hydrophobic Substrates

Langmuir, Vol. 20, No. 11, 2004 4761

Figure 6. (a) Temperature dependence of the water density profile close to a planar hydrophobic wall at a fixed pressure of p ) 1 bar. The Lennard-Jones potential of the hydrophobic layer is denoted by a broken line (to the left). (b) Pressure dependence of the water profile for a fixed temperature of T ) 300 K. (c) Electrostatic potential distribution for a fixed pressure of p ) 1 bar and various temperatures. (d) Electric field distribution for a fixed pressure of p ) 1 bar and temperature of T ) 300 K. In all cases, the simulation box contained 64 hydrophobic molecules and 2781 SPC/E water molecules. The system was thermalized for 100 ps and averaged for 2 ns.

between one methylene group and a water molecule over the volume of the hydrocarbon medium, weighted with the alkane (or hydrocarbon) density profile, Fhc(z).

U h (z) )

∫ dz′ Fhc(z′) ∫ dr 2πrUsw(x(z - z′)2 + r2)

In addition to the density profile, we computed the electrostatic potential and the normal electric field across the water-alkane interface, which follow from the expressions

(4)

∫-∞z dz′ ∫-∞z′

Ψ(z) - Ψ(-∞) ) -

Substituting eq 3 in eq 4 we obtain

U h (z) ) π

(

C(12) sw

C(6) sw

)

∫ dz′ Fhc(z′) 5(z - z′)10 - 2(z - z′)4

(5)

In Figure 5, it is seen that the potential U h (z) is almost purely repulsive (the attractive region is not discerned on the scale of the figure and is therefore irrelevant) and quite close in shape to the exponentially decaying potentials employed by us for the case of spherical solutes. This a priori unexpected potential shape is due to the folding of the continuously decaying alkane density profile with the laterally averaged LJ potential in eq 5. It is clearly visible that the laterally averaged potential, U h (z), is shifted toward the inside of the alkane slab, which can be understood, since the methylene-water potential contains attractive components. This has important consequences for the definition of the depletion thickness: Due to this shift, a definition of the depletion thickness based on the interaction potential (which is done in most theories and in many simulations) will overestimate the numerical value when compared with experimental values (which are based on the density profile of the hydrophobic substrate). In Figure 6a and b, we show the temperature and pressure dependence of the water density profile close to the hydrophobic surface. We also show (as a broken line to the left) the laterally averaged potential between the alkane slab and the water molecules, U h (z), which is more or less the same for all different temperatures and pressures. In analogy to the results for the largest spherical solute considered in the previous section, the depletion zone increases with increasing temperature21 (one notes a small irregularity between T ) 320 K and T ) 340 K, since here, the order is inverted; apparently, the depletion layer width shows a nonmonotonic temperature dependence). In agreement with common expectation, the depletion layer shrinks for very large applied pressures, but as we have observed before for the case of the largest spherical solute, the pressures needed to compress the depletion layer are quite large.

E(z) )

∫-∞z

Fc(z′′) dz′′ 0

Fc(z′) dz′ 0

where Ψ(z) denotes the laterally averaged electrostatic potential, E(z) denotes the laterally averaged normal electric field, Fc(z) denotes the charge density, and 0 represents the vacuum permittivity. The results for the potential and electric field distributions due to the water orientation at the interface are shown in Figure 6c and d. Water molecules in the interfacial layer are oriented such that the hydrogen atoms point predominantly away from the water half-space. The mechanism for this is the reduction of the image charge repulsion, as can be understood from the following simple argument: The image charge repulsion of a charge is proportional to the square of the charge strength. The total image charge repulsion of a single charge (i.e., an oxygen atom) is therefore twice as large as the total image charge repulsion of two charges each with half-charge strength (i.e., two hydrogen atoms), assuming both arrangements have the same distance from the interface. Due to this orientation, a huge dipole moment is built up, which leads to a large potential drop of ∼0.5 V, which is in good agreement with previous results for slightly different definitions of the hydrophobic substrate.24,25 This shows that the large surface potential of the water interface is a robust property; the consequences of this, especially the interaction of the dipolar layer with other charged solutes (e.g., salt ions), have only been addressed very recently.23,26 Figure 6c demonstrates that this potential is almost independent of the temperature. To quantify the density depletion layer, we need a prescription to calculate the depletion thickness in an unambiguous manner. In principle, there are two different ways of doing this. In the first definition, which is based on the interaction between the alkane slab and the water layer, we add the normalized water density profile to the reduced average water-alkane slab potential (which is shown as a solid line in Figure 7), calculate the difference of this function to unity (which is the broken line in Figure 7), and integrate over this function. The result is the

4762

Langmuir, Vol. 20, No. 11, 2004

Mamatkulov et al.

half the length obtained in recent neutron scattering experiments19,20 and twice the length obtained with X-rays at hydrophobic substrates;21 possible reasons for the discrepancies between different experiments will be discussed in the next section. V. Conclusions and Discussion

Figure 7. The depletion distance d1 was determined by adding the Lennard-Jones potential of the hydrophobic slab in units of kBT to the water density (solid line) and integrating the difference to unity (broken line) over the vertical coordinate. We here show explicit data for T ) 360 K and p ) 1 atm. Table 4. Hydrophobic Wall-Water Depletion Distance at Various Temperatures and Pressuresa p ) 1 bar T (K)

d1 (Å)

d2 (Å)

273 300 320 340 360

2.880 3.202 3.920 3.540 4.832

2.128 2.565 2.669 2.615 2.946

T ) 300 K p (bar)

d1 (Å)

d2 (Å)

1 500 1000 2000

3.202 2.837 2.760 1.990

2.565 2.310 1.896 1.845

a d is based on the wall-water potential, and d is based on the 1 2 sum of the alkane and water densities.

depletion length with respect to the water-slab potential, which we denote as d1

d1 )

(

U h (z)

F (z)

)

w ∫ dz 1 - kBT - Fw,bulk

where the lower integration boundary is located at the point where the reduced slab potential, U h (z)/kBT, is unity. This length is not measured in experiments but rather corresponds to the length scale extracted in theories where the hydrophobic slab is purely defined by its interaction potential. Similarly, we can calculate a depletion distance on the basis of the sum of the two densities of the hydrocarbon slab and of the water layers, which is closer to what would be measured in a scattering experiment. We denote this depletion distance by d2.

d2 )

(

F (z)

F (z)

)

hc w ∫ dz 1 - Fhc,bulk Fw,bulk

Results for both depletion distances are presented in Table 4 as a function of varying pressure and temperature. As can be seen, the depletion distance d2 based on the density profiles of the alkane slab and the water layers is significantly smaller than the depletion thickness d1 based on the wall potential. At room temperature and normal pressure, we obtain d2 ) 2.565 Å, which is close to one recent neutron scattering experiment19 and to ellipsometry measurements.22 On the other hand, our result is about

Using MD simulations, we have investigated the behavior of water at spherical and planar hydrophobic objects as a function of pressure, temperature, and curvature. We have concentrated on effects which are important for the understanding of the recently observed water density depression at planar hydrophobic substrates.19-21 In agreement with previous simulations and theories, we find small cavities (with a radius of 4 Å) to be surrounded by a water layer of increased density, whereas large cavities (with a radius of 18 Å) lead to a depletion layer with a width of ∼2 Å. A cavity with a radius of 10 Å is roughly neutral in that the averaged water density close to the solute is of the order of the bulk density of water. For the planar substrate, we fixed the lateral interfacial area and employed a constant pressure in the normal direction only. This ensemble reduces spurious effects due to the interfacial energy, which can be quite large for small systems. Our main results are: (i) The depletion thickness of a planar hydrophobic substrate (made up of an assembly of mobile alkane chains oriented normal to the interface) is d1 ≈ 3 Å when measured with respect to the reduced substrate potential and d2 ≈ 2.5 Å when measured with respect to the substrate density profile at normal pressure and room temperature. The latter number is relevant when comparing with scattering results, and the former number is typically used in analytical theories. (ii) The depletion layer thickness depends sensitively on the radius of curvature of the solute. For a spherical solute with a radius of R ) 18 Å, we find d1 ≈ 2 Å, and for R ) 10 Å, we find d1 ≈ 0.5 Å. For a real substrate as used in experiments, this means that even small amounts of surface roughness will tend to reduce the depletion layer thickness. (iii) At high temperature, the depletion layer thickness increases substantially. This effect might be used in the experiments to distinguish a vacuum layer of depleted water from a layer of surface-adsorbed contaminants. (iv) Contrary to our naive expectations, the depletion layer shows considerable resistance toward compression. Only at pressures of ∼2000 bar does the depletion layer shrink significantly. (v) The electrostatic surface potential is +0.5 V at a hydrophobic planar substrate and is quite insensitive to temperature variations. It is caused by the water orientation in the topmost surface layer. Our results for the depletion layer thickness on the planar substrate, d2 ≈ 2.5 Å, lie somewhat between the neutron scattering results for self-assembled monolayers, which point to a depletion layer thickness of ∼5 Å,20 and the X-ray scattering results, which give a value of 1 Å.21 However, our finding for the depletion thickness is close to the results of recent elliposometry measurements22 and to the neutron scattering results on polystyrene substrates.19 Bearing in mind that the depletion layer shrinks considerably even for as large of a radius as R ) 18 Å, it is easy to explain smaller values for the depletion layer thickness in experiments: Maybe, here, the substrate shows more intrinsic roughness than in our simulations (one should note that our planar substrate exhibits also

Water at Hydrophobic Substrates

a certain degree of roughness, since it is composed of mobile alkane chains). It is more difficult to interpret the larger depletion layer thicknesses found in one of the two recently performed neutron scattering experiments. Maybe imperfections of the self-assembled monolayer could be a possible reason for the discrepancy among different experiments. A second important point is the dependence of the depletion layer thickness on the degree of hydrophobicity. It has been convincingly demonstrated21 that even small amounts of attractive interactions between the hydrophobic slab and the water molecules reduce the depletion layer width and therefore have to be carefully taken into account in any comparison between experiment and theory. The at first sight astonishingly small compressibility of the depletion layer can be rationalized by a very simple argument. We assume that the relevant free energy scale which is causing the depletion layer is 1 kBT ) 4 × 10-21 J and that this free energy is spent on a single water molecule with a volume of V ) 0.003 nm3. This reflects the fact that water depletion is to some degree driven by entropic effects (a similar argument for the interfacial tension leads to γ ) kBT/V2/3 ≈ 200 mN/m, which correctly predicts the order of magnitude of an experimental value of γ ) 70 mN/m). This is not to say that water in the first interfacial layer is fully oriented; rather, it might be oriented if it was located closer to the substrate at a distance over which water in fact is depleted. The entropic pressure following from that is p ) kBT/V ≈ 1.3 × 109 J/m3 ≈ 1300 bar. This explains the high pressures needed to compress the depletion layer in our simulations, p ) 10002000 bar, which are thus of the order of the entropic pressure. The compressibility of the depletion layer follows from the data in Table 4 to be of the order of κ ≈ 2 × 10-4 bar-1 and is thus only ∼5 times higher than that of bulk water (note that the bulk compressibility which follows from our data in Table 3 is of the order of κ ≈ 3 × 10-5

Langmuir, Vol. 20, No. 11, 2004 4763

bar-1 and thus is 30% smaller than the experimental value). The consequences of the electrostatic potential jump at the water surface can only be vaguely appreciated. A potential of 0.5 V means that a negative ion adsorbing on the surface gains a free energy of ∼20kBT. Adsorption of negative ions should thus be observed at free air-water interfaces. Indeed, it is known that most aqueous interfaces, especially at hydrophobic surfaces, appear to be negatively charged.16 MD simulations, however, show that the interaction between ions and the surface layer of water is quite complicated and dominated by ion-specific effects.26 Many experiments in the past were concerned with the interaction between two macroscopic hydrophobic bodies and found attractions with decay lengths in the whole range from nanometers to micrometers. What is the relation between the depletion layer thickness as defined in our paper and the experimentally measured decay length of the interaction between two hydrophobic surfaces? Since it is now believed that the very long decay lengths observed for the interaction are due to bubble formation and bridging,19 our decay lengths have very little to do with the measurements of interactions between macroscopic bodies, since the water depletion profile changes as two such bodies are brought into contact, leading to cavity and bubble formation. However, the density profile as calculated in this paper might be relevant for adsorption of small hydrophobic solutes (i.e., simple molecules or polymers) on hydrophobic substrates. We are planning to perform simulations on such systems in the future. Acknowledgment. This work was funded by the Alexander-von-Humboldt Stiftung and the DFG through SFB 486 and SFB 563. LA036036X