Langmuir 2008, 24, 135-141
135
A Computer Simulation Study of Stick-Slip Transitions in Water Films Confined between Model Hydrophilic Surfaces. 1. Monolayer Films Alexander Pertsin*,†,‡ and Michael Grunze† Angewandte Physikalische Chemie, UniVersita¨t Heidelberg, Im Neuenheimer Feld 253, 69120 Heidelberg, Germany, and Institute of Organo-Element Compounds, Russian Academy of Sciences, 28 VaViloV Str., Moscow 119991, Russia ReceiVed July 21, 2007. In Final Form: September 24, 2007 The shear behavior of monolayer water films confined in a slit-like pore between hydrophilic walls is simulated in the quasistatic regime using the grand canonical Monte Carlo technique. Each wall is represented as a hexagonal lattice of force sites that interact with water through an orientation-dependent hydrogen-bonding potential. When the walls are in registry, the water oxygen atoms form either a crystal- or fluid-like structure, depending on the period of the wall’s lattice. In both cases, however, the monolayer structure is orientationally disordered. Both the crystaland fluid-like monolayers prove to be capable of experiencing well-defined stick-slip transitions, with the largest yield stress occurring in the crystal-like case. Beyond the yield point, the crystal-like monolayers “melt”, but their structure and molecular motion differ in many respects from those characteristic of normal fluids.
1. Introduction There are many physical situations, ranging from nanotribology to biology, where water plays the role of a lubricant confined between closely spaced surfaces. When the confined water film becomes molecularly thin, the lubricating properties of water exhibit dramatic changes associated with the effects of confinement and interfacial interactions on the structure and rheology of water. In general, the frictional resistance of fluid films can be divided into three regimes of film thickness.1 Films thicker than 5-10 molecular dimensions flow like a bulk fluid. Below this thickness range, the effective viscosity of confined fluid rises, and the shear response of the system exhibits an elastic component in addition to a viscous one. As the film thickness is reduced below 2-4 molecular dimensions, the effective viscosity diverges and the system behaves like a solid. In this regime, sliding over a macroscopic distance cannot be accomplished until a critical shear stress (yield stress) is exceeded. This is the so-called “stick” phase of the shear response. Beyond the yield point, the shear stress abruptly declines, and the system enters the “slip” phase. The stick-slip transition has been observed in nearly all fluids studied, including oils and organic solvents. As to water and aqueous solutions, their shear behavior in the ultrathin-film regime still remains in many respects unclear. Thus, Granick and his colleagues1-3 studied the shear response of aqueous electrolyte films using a surface force apparatus (SFA) modified to measure oscillatory shear forces. The films were confined between mica sheets to a thickness D ) 6 ( 2 Å. The shear response was found to be dependent on the shear rate. The stick-slip transition was only observed when the films were deformed faster than their intrinsic relaxation time (∼0.01 s). At low deformation rates (in the linear response regime), the viscoelastic properties of the nanoconfined films were characteristic of fluids, although the effective viscosity was orders of magnitude higher than that of † ‡
Universita¨t Heidelberg. Russian Academy of Sciences.
(1) Granick, S. Phys. Today 1999, 52, 26. (2) Dhinojwala, A.; Granick, J. Amer. Chem. Soc. 1997, 119, 241. (3) Zhu, Y.; Granick, S. Phys. ReV. Lett. 2001, 87, 096104.
bulk water. Similar results were reported by Klein et al.,4-6 except that confined water and aqueous salt solutions retained most of their fluidity. (At D < 10 Å, the increase in viscosity was within a factor of 3). The high fluidity of water nanoconfined between mica surfaces was also observed in molecular dynamics (MD) simulations for three-layer and thicker films, whereas bilayer films (D ≈ 6 Å) showed a significant increase in shear viscosity.7 On the other hand, there have been a few computer simulations, albeit not concerned directly with the shear behavior, that clearly showed the ability of water to crystallize when confined to monoand bilayer thickness.8-13 The crystallization occurred both in the case of structureless nonorienting hydrophobic walls8-10 and with atomically structured walls bearing a lattice of LennardJones (non-hydrogen-bonding) force sites.11-13 Stable crystallike water monolayers were also observed experimentally in a neutron diffraction study of the structure of water confined in a slit-shaped organic nanopore.14 In all the cases,8-14 the confined water films would thus be expected to exhibit the slip-stick transition if their shear behavior would be examined. The apparent contradiction between the ability of nanoconfined water films to crystallize8-14 and their ability to retain fluid-like shear response when confined between mica sheets1-7 originates, in our view, from the high specificity of the water-mica system. It is known that mica surfaces, when in contact with water, acquire a negative charge, because the surface K+ ions desorb and dissolve in water.15 When the negatively charged mica sheets are brought (4) Raviv, U.; Laurat, P.; Klein, J. Nature 2001, 51, 413. (5) Raviv, U.; Klein, J. Science 2002, 297, 1540. (6) Klein, J.; Raviv, U.; Perkin, S.; Kampf, N.; Chai, L.; Giasson, S. J. Phys.: Condens. Matter 2004, 16, S5437. (7) Leng, Y.; Cummings, P. T. J. Chem. Phys. 2006, 94, 026101. Note that shear rates accessible in MD simulations are several orders of magnitude higher than those typical of SFA measurements. (8) Koga, K.; Tanaka, H.; Zeng, X. C. Nature 2000, 408, 564. (9) Koga, K. J. Chem. Phys. 2002, 116, 10882. (10) Slovak, J.; Tanaka, H.; Koga, K.; Zeng, X. C. Physica A 2003, 319, 163. (11) Zangi, R. J. Phys.: Condens. Matter 2004, 16, S5371. (12) Zangi, R.; Mark, A. E. J. Chem. Phys. 2003, 119, 1694. (13) Zangi, R.; Mark, A. E. Phys. ReV. Lett. 2003, 91, 025502. (14) Janiak, C.; Scharmann, T. G.; Mason, S. A. J. Am. Chem. Soc. 2002, 124, 14010. (15) McGuiggan, P. M.; Israelachvili, J. N. J. Mater. Res. 1990, 5, 2232.
10.1021/la702209g CCC: $40.75 © 2008 American Chemical Society Published on Web 11/30/2007
136 Langmuir, Vol. 24, No. 1, 2008
to the spacing of one to two water diameters, the counterions rush into the intersheet space to neutralize the surface charge. The resulting local concentration of counterions between the sheets may reach as high as 15 M, as estimated by Zhu and Granick.3 It can therefore be conjectured that it is the surface charge and huge counterion concentration that prevent confined water from experiencingsat least, at low shear ratessthe slipstick transition between mica sheets. The objective of our work is to support the above conjecture by demonstrating that pure ionless water confined between electrically neutral walls can undergo the stick-slip transition. To this end, we simulate the shear behavior of pure water confined between model walls bearing no electric charges.16 To avoid conceptual problems involved in the treatment of stick-slip transitions using nonequilibrium MD (see ref 17 for details), we resort to the “quasistatic” approach suggested and developed by Schoen, Diestler, and their colleagues.17-23 In this approach, which corresponds to an infinitely small shear rate, the thermodynamic state of the water film passes through a succession of equilibrium states, each distinguished by a different lateral alignment of the confining walls. In our simulation study, the properties of the film in each of these equilibrium states are calculated using the grand canonical Monte Carlo (GCMC) technique at fixed temperature, volume, and chemical potential, while the number of water molecules is allowed to fluctuate. The choice of the GCMC technique is dictated by the conditions of SFA experiments, where the confined water is open to the bulk water reservoir. To avoid capillary evaporation, we consider hydrophilic walls, such that the water-wall binding energy is the same as the interaction energy of water molecules between themselves. The lateral coupling of the wall to confined water, which is necessary for inducing shear stress in the water film, is ensured by a pattern of discrete force sites arranged on the walls. The shear behavior of confined water is analyzed in terms of the shear stress-strain dependence, in parallel with an analysis of the structure of confined water. In this paper, we restrict ourselves to water films of monolayer thickness, while leaving bilayer and thicker films for a separate paper. 2. Methods 2.1. Interaction Potentials. In calculations of the potential energy of the system, the interactions between water molecules are described with the TIP4P model,24 which represents the water molecule as a Lennard-Jones force site at the oxygen atom, two positive partial charges on the hydrogens, and a compensating negative charge moved off the oxygen to a point 0.15 Å along the bisector of the H-O-H angle. The Lennard-Jones part of the water-water potential is cut off at 8.2 Å, while the electrostatic Coulombic interactions are smoothly damped in the range between 8 and 8.5 Å using a switching function suggested by Lee and Rossky.25 No correction for the long(16) Similar models, where surface charges and counterions were deliberately omitted, were used in earlier studies of water confined between silicate surfaces (see, for example, Mulla, D. J. J. Colloid Interface Sci. 1984, 100, 576). (17) Bordarier, P.; Schoen, M.; Fuchs, A. H. Phys. ReV. E 1998, 57, 1621. (18) Schoen, M.; Cushman, J. H.; Diestler, D. J. Mol. Phys. 1994, 81, 475. (19) Schoen, M.; Hess, S.; Diestler, D. J. Phys. ReV. E 1995, 52, 2587. (20) Schoen, M.; Diestler, D. J.; Cushman, J. H. J. Chem. Phys. 1994, 100, 7707. (21) (a) Cushman, J. H. Nature 1990, 347, 227. (b) Curry, J. E.; Zhang, F.; Cushman, J. H.; Schoen, M.; Diestler, D. J. J. Chem. Phys. 1994, 101, 10824. (c) Curry, J. E.; Cushman, J. H. Mol. Phys. 1995, 85, 173. (22) Bock, H.; Schoen, M. J. Phys.: Condens. Matter 2000, 12, 1545. (23) (a) Schoen, M.; Bock, H. J. Phys.: Condens. Matter 2000, 12, A333. (b) Bock, H.; Schoen, M. J. Phys.: Condens. Matter 2000, 12, 1569. (c) Bock, H.; Diestler, D. J.; Schoen, M. Phys. ReV. E 2001, 64, 046124. (d) Bock, H.; Diestler, D. J.; Schoen, M. J. Phys.: Condens. Matter 2001, 13, 4697. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (25) Lee, S. H.; Rossky, P. J. J. Chem. Phys. 1994, 100, 3334.
Pertsin and Grunze Table 1. Parameters of the Wall-Water Potential in Eq 1 interaction
d, Å
, kcal/mol
m
n
O···A H···A
3.85 2.14
0.26 1.58
6 8
12 12
range electrostatic interactions is made, because the original TIP4P model was calibrated without such a correction. (Our earlier GCMC simulations of bulk TIP4P water26 showed that the inclusion of the long-range electrostatic contribution worsened the agreement with the experimental water density under ambient conditions, as compared to the original TIP4P model.) The wall-water potential used in our simulations mimics the interaction of water with a hydrophilic proton-acceptor surface. The potential is made orientation dependent to reflect the trend of hydrogen bonds to the linear configuration A···H-O, where A denotes a protonacceptor force site (or an “atom”) on the confining wall. The force sites are arranged on a two-dimensional hexagonal lattice with a period l. The interactions of a force site with the oxygen and hydrogen atoms of a water molecule are described by inverse power (m - n) potentials of the form φ(r) ) [n(r0/r)n - m(r0/r)m]/(m - n)
(m < n)
(1)
where and r0 are the depth and position of the potential well. The interaction energy of a lattice site a and a water molecule w is calculated as the sum φaw ) ∑k (w)φak(rak), where k runs over the oxygen and two hydrogens of molecule w. To adapt φaw to modeling hydrogen bonds A···H-O, we followed a simple expedient similar to that used in an early work by McGuire et al.27 in describing hydrogen bonding. In the model suggested by the authors,27 there are no explicit angular terms, and the preference of the linear A·· ·O-H configuration is achieved by making the A···H interaction much stronger and somewhat shorter ranged than the A···O one. With a judicious choice of potential parameters, the region of the A···H potential minimum is overlapped with the repulsive wing of the A···O potential, and it is this overlap that makes the linear A·· ·H-O configuration most favorable: Any deviation from linearity increases the overlap and makes the A···H-O interaction more repulsive. The specific potential parameters used in our simulations are presented in Table 1. The parameters are normalized so that the binding energy of the A···H-O bond is exactly -1 kcal mol-1. The equilibrium A···O separation is calibrated to be 3.4 Å. The angular dependence of the associated potential φaw is similar to that of the water-graphite carbon potential suggested in our earlier work (see Figure 2 in ref 28). In simulations of wall-water interactions, the A···H and A···O potential well depths were multiplied by a scaling factor η such that the binding energy of a single water molecule to the wall with a particular lattice period l was equal to that of the lowest energy TIP4P water dimer (-6.24 kcal mol-1). In other words, a water molecule saw the wall as being as favorable for binding as another water molecule. In our choice of the lattice period l, we did not seek to model a specific solid surface. This gave us the freedom to vary the lattice period l and to explore the effect of l on the shear behavior of confined water films. Three lattice periods were tried: l ) 2.76 (the shortest O···O distance in ice Ih), 3.0, and 3.2 Å. The corresponding values of η were found to be 1.65, 1.88, and 2.03, respectively. 2.2. Simulation Procedure. In GCMC simulations of open confined systems, the chemical equilibrium between the system and the (fictitious) bulk reservoir is maintained by allowing density fluctuations through particle insertions and deletions. The basic weakness of GCMC is a too low probability of performing successful insertions and deletions at densities typical of fluids. To improve the acceptance of insertion and deletion attempts, we employed the (26) Pertsin, A.; Platonov, D.; Grunze, M. J. Chem. Phys. 2005, 122, 244708. (27) McGuire, R. F.; Momany, F. A.; Scheraga, H. A. J. Phys. Chem. 1972, 76, 375. (28) Pertsin, A.; Grunze, M. J. Phys. Chem. B 2004, 108, 1357.
Water Films Confined between Hydrophilic Surfaces
Langmuir, Vol. 24, No. 1, 2008 137
excluded volume mapping29 and Swendsen-Wang filtering30 techniques, which allowed us to reject improbable positions for insertion or deletion of a water molecule based on a computationally inexpensive predictor (for details, see our previous paper).31 Further improvement in the sampling efficiency was achieved due to implementation of a rotational bias procedure,32 which made it possible to reject improbable molecular orientations. For computational convenience, the unit cell of the hexagonal lattice of the wall’s force sites was taken to be orthogonal, with periods a ) x3l, b ) l. The unit cell contained two force sites, one at the origin and the other at the center of the unit cell. The simulation box was a rectangular prism of height h (the slit width) with lateral dimensions Lx ) nxa and Ly ) nyb. The integers nx and ny were taken equal to 7 and 11, respectively, so that Lx and Ly were roughly 30 Å. The increase of nx to 11 and ny to 18 (Lx ≈ Ly ≈ 50 Å) kept the simulation results practically unchanged, which suggested that system-size effects were insignificant. A typical length of the GCMC run was (2-4) × 106 GCMC cycles, each comprising N0 moves, where N0 was the initial number of molecules in a given cycle. The total number of distinct configurations accepted in a single GCMC run came to 108. (For a rough comparison, the number of configurations in a 100 ns MD trajectory with a typical time step of 4 fs is 2.5 × 107.) The simulations were carried out at room temperature. To ensure the equilibrium of the confined system with water bulk, the density-corrected chemical potential33 of the confined system was equated to the value corresponding to bulk water at ambient conditions, µ′′ ) -6.204 kcal mol-1, as determined in separate simulations of bulk TIP4P water. The key quantity calculated in our simulations was the stress tensor, TRβ (R, β ) x, y, z). Its diagonal component Tzz determined the normal pressure, PN, whereas the nondiagonal components Tzx and Tzy represented the x- and y-components of the shear stress. The explicit formulas for TRβ can be found in ref 28. By plotting TRβ as a function of their conjugate strains, we could determine the respective yield stress. In interpreting the observed shear behavior, standard distribution functions and order parameters characterizing the water film structure were calculated. The mobility of water molecules was assessed in terms of the mean-square displacement (MSD), d, and its components d ) dx + dy + dz dχ ) dχ(n) ) dχξ ) dχ + dξ
1
N
∑(χ
N i)1
i 0
- χin)2
(2)
(χ, ξ ) x, y, z)
where χ0i and χni are, respectively, the χth coordinates of the oxygen atom of molecule i in the beginning and at the end of the Monte Carlo run of length n. Note that dxy ≡ d| and dz ≡ d⊥ represent the lateral and normal components of d. The calculations of MSD were performed within the canonical (NVT ) const) ensemble, with the number of water molecules N fixed at the respective grand canonical ensemble average, 〈N〉. The positions of the force sites in the upper and lower walls were related by the equations x′ ) x, y′ ) y + βb, and z′ ) z + h, where the prime refers to the upper wall. That is, we considered only shear strains in the y direction. Note that for β ) (k, the walls are in registry, and they are out of registry for β ) ((k + 1)/2 (k being an integer). That is, β can well be used as a quantitative measure of registry.
3. Results and Discussion We began our simulations with the scanning of a fairly wide range of wall-to-wall separations, 6.0 e h e 12.4 Å, to explore (29) Stapleton, M. R.; Panagiotopoulos, A. J. Chem. Phys. 1990, 92, 1285. (30) Swendsen, R. H.; Wang, J.-S. Phys. ReV. Lett. 1987, 58, 86. (31) Pertsin, A.; Grunze, M. J. Phys. Chem. B 2004, 108, 16533. (32) Shelley, J. C.; Patey, G. N. J. Chem. Phys. 1995, 102, 7656. (33) Pertsin, A.; Grunze, M. J. Chem. Phys. 2006, 125, 114707.
Figure 1. The dependence of the normal pressure PN on the slit width h for different wall lattice periods l at β ) 0.
Figure 2. The relative change in the normal pressure PN as a function of registry β for different l at h ) 6.4 Å.
the general thermodynamic behavior of the system when the confining walls are in registry. At all separations h tried, the confined water film proved to be thermodynamically stable: Even when the starting configuration of the system included only one water molecule, the GCMC run led to capillary condensation and the formation of a continuous water film. As expected, the film had a layered structure and the normal pressure PN showed pronounced oscillations (Figure 1). At compressive pressures (PN > 0), water monolayers were stable up to h ) 6.5 Å, water bilayers were formed at 7.8 e h e 9.3 Å, and the range of water trilayers (not shown in Figure 1) extended from 10.4 to 12.4 Å. The behavior of PN(h) at different l’s was essentially similar. The main difference was that PN(h) slightly shifted to shorter h as l was increased. The shear was simulated by displacing the upper wall relative to the lower one along the y direction by a fixed strain βb and then allowing the system to equilibrate. Thereafter, β was incremented, the simulation was repeated, and so on. Considering that PN(β) is symmetric and Tzy(β) is antisymmetric with respect to β ) 0 and 0.5, the registry β was varied only in the symmetrically independent range, 0 e β e 0.5. At h ) 6.4 Å, which corresponds to a monolayer film at all l tried (Figure 1), PN monotonously increased with increasing registry β. To make allowance for the differences in PN(0) at different l, we illustrate in Figure 2 the behavior of PN(β) in terms of its relative change, ∆pN ≡ [PN(β) - PN(0)]/PN(0). The largest increase in ∆pN was found at l ) 3 Å. In absolute units, it resulted from the change in PN from 2.5 to 8 kbar. The same trend was observed for the shear stress Tzy, which showed the strongest dependence on registry at l ) 3 Å (Figure 3). For all three lattice periods tried, the yield point, βyd, occurred at β ) 0.2 (to an accuracy determined
138 Langmuir, Vol. 24, No. 1, 2008
Pertsin and Grunze
Figure 3. The shear stress along y as a function of registry β for different l at h ) 6.4 Å.
Figure 4. The radial distribution function for different l and β at h ) 6.4 Å.
by the increment size, ∆β ) 0.1), and the associated magnitudes of the yield stress, S ) -Tzy(βyd), were 1.5, 2.2, and 1.3 at l ) 2.76, 3.0, and 3.2 Å, respectively. It is of interest to represent the different stress-strain curves (Figure 3) in terms of reduced quantities, β˜ ) β/βyd, T˜ zy ) Tzy/S. Such a renormalization was suggested by Bock and Schoen22 based on a theory of corresponding states. It was applied to ultrathin atomic fluid films that were confined between chemically corrugated substrates composed of an alternating sequence of weakly and strongly adsorbing strips. It was shown that at β < βyd the function T˜ zy(β˜ ) was fairly insensitive to the relative width of strongly adsorbing strips; i.e., T˜ zy(β˜ ) was unique for different chemical corrugations. In our case, we have only one value of registry, β ) 0.1 (β˜ ) 0.5), to check the sensitivity of T˜ zy(β˜ ) to the model parameter, l. The calculated values of T˜ zy(0.5) proved to be scattered in a narrow range from 0.65 to 0.75, which is comparable to the scatter of T˜ zy(0.5) in the model system studied by Bock and Schoen.22 Thus, it appears that the approach suggested by the authors22 is applicable to our model system as well. The simulated dependence of the yield stress on the applied load, S(PN), was well-described by a linear relation
Table 2. Coefficients of Eq 3 as a Function of the Wall’s Lattice Period
S ) S0 + CPN
(3)
which is generally used in interpreting experimental shear force results.34,35 The second term in the right-hand side of eq 3 represents the contribution to the friction force due to the external load. At high loads, such that CPN . S0, the slope parameter C takes the meaning of the effective friction coefficient. The first term, which is responsible for the resistance to shear at zero load, is usually regarded as an internal contribution to S arising from the adhesion force between the surfaces. Note that in our model system the direct interaction between the confining walls is not considered, and S0 comes solely from the wall-water interaction. The simulation results for the coefficients of eq 3, which are presented in Table 2, show that the confined water possesses a substantial resistance to shear. Thus, the calculated values of S0 are 2 orders of magnitude larger than the experimental values typical of classical boundary lubricants, such as fatty acid monolayers,36 and about 1 order of magnitude larger than the data for monolayers of octamethylcyclotetrasiloxane (OMCTS) (34) Kumacheva, E.; Klein, J. J. Chem. Phys. 1998, 108, 7010. (35) Gee, M. L.; McGuinnan, P. M.; Israelachvili, J. N.; Homola, A. M. J. Chem. Phys. 1990, 93, 1895. (36) Tabor, D. Proc. Inst. Mech. Eng. 1991, 205, 365.
l, Å
S0, kbar
C
2.76 3.0 3.2
1.1 1.7 0.97
0.08 0.11 0.07
and simple hydrocarbons like cyclohexane.34,35 For the simulation model used, the large adhesion contribution S0 is hardly surprising in view of the presence of strong wall-water interactions modeling hydrogen bonding. The results for the slope parameter C are close to 0.1 (Table 2). This is an order of magnitude smaller than the experimental values of C for simple hydrocarbons and OMCTS34,35 but larger than for fatty acid monolayers (0.030.04).36 To understand the origin of the above-discussed shear behavior, we now turn to the structure and molecular mobility of water in confined films. We begin with the films that are formed at l ) 3 Å and show the most pronounced stick-slip transition. A remarkable property of these films was that their constituent water oxygens formed a perfect hexagonal lattice identical to the hexagonal lattice of force sites on the confining walls. The ensemble average number of water molecules in the simulation cell, 〈N〉, was found to be 154, which was exactly equal to the number of force sites on each of the walls (nx × ny × 2 ) 7 × 11 × 2 ) 154). The ideal commensurability (matching) of the wall and water lattices was so advantageous that 〈N〉 remained constant in a fairly wide range of film thickness, from 6.2 to 7.0 Å. The monolayer’s in-plane radial distribution function, g(r), displayed long-range correlations characteristic of a hexagonal lattice with a period l ) 3 Å (Figure 4). During the simulation run, MSD showed only slight oscillations around ∼0.33 Å2 (Figure 5). This is noticeably larger than the values observed in ice by X-ray and neutron diffractometries (∼0.04 Å2).37 Nevertheless, based on the form of g(r) and the independence of MSD on the length of the simulation run, the water monolayer formed at l ) 3 Å can definitely be regarded as a two-dimensional crystal. Note that, even in this crystal-like monolayer, the normal component of MSD was much smaller than the lateral component: d⊥ ) 0.06 Å,2 d| ) 0.27 Å.2 In all other monolayers to be discussed in the following, d⊥ remained practically unchanged, whereas the lateral mobility of water molecules substantially increased. As a consequence, the contribution of d⊥ to d became negligible. (37) Goto, A.; Hondoh, T.; Mae, S. J. Chem. Phys. 1990, 93, 1412.
Water Films Confined between Hydrophilic Surfaces
Figure 5. The mean-square displacement of water oxygens as a function of the length of the simulation run, n, for different l at h ) 6.4 Å and β ) 0.
Figure 6. A typical configuration of the water monolayer film at h ) 6.4 Å, l ) 3 Å, and β ) 0.
Figure 7. The cosine distribution of the orientations of OH bonds for different separations from the midplane at l ) 3 Å, h ) 6.4 Å, and β ) 0.
Despite the high long-range translational order in the arrangement of water oxygens, the monolayer structure at l ) 3 Å was orientationally disordered. This can well be appreciated from the snapshot of the system in Figure 6. There were, however, some preferences in molecular orientations. In the middle of the water monolayer, most O-H bonds assumed orientations close to 90° with respect to the z-axis (Figure 7). These orientations were obviously associated with hydrogen bonds between water molecules and the walls. Somewhat closer to the walls, the most preferred orientation of the O-H bonds changed to 0 and 180° for the upper and lower walls, respectively, which pointed to a growth in the proportion of hydrogen bonds formed by water with the wall atoms. An analysis of the hydrogen-bonding network showed that all water molecules in the monolayer formed fairly strong hydrogen bonds with, at least, one of the walls. The change in the wall lattice period l by only 7-8% (from 3.0 to 2.76 and 3.2 Å) broke down the favorable wall-water matching. At h ) 6.4 Å, 〈N〉 changed from the initial 154 molecules at l ) 3 Å to 126 and 171 molecules at l ) 2.76 and
Langmuir, Vol. 24, No. 1, 2008 139
Figure 8. The mean-square displacement of water oxygens as a function of the length of the simulation run, n, for different β at h ) 6.4 Å and l ) 3 Å.
3.2 Å, respectively. The corresponding g(r)’s showed no signs of long-range order and were typical of fluids (Figure 4). Further evidence in favor of the fluid-like nature of the water film was provided by the behavior of MSD (Figure 5), which could be described by a linear function of the number of MC cycles, n, similar to the time dependence of MSD for the Brownian selfdiffusion.38 The fluid-like nature of the monolayer water films at l ) 2.76 and 3.2 Å may appear to contradict the observed ability of these films to sustain a considerable shear stress, only moderately lower than that of the crystal-like film at l ) 3.0 Å (see Figure 3 and Table 2). A similar shear behavior was reported earlier by Schoen et al.19 for simple atomic monolayers confined between atomically structured walls. The authors19 explained the apparent contradiction between the observed fluidity and the resistance to shear in terms of a so-called “pinning” mechanism. According to this mechanism, the film atoms were trapped in effective holes formed by their nearest-neighbor wall atoms, thereby resisting transversal displacement of the opposing walls relative to each other. The resistance of atomic fluid films to shear was also observed with chemically corrugated substrates composed of weakly and strongly adsorbing strips.22,23 Under favorable thermodynamic conditions, the confined fluid formed a “bridge” phase, in which high-density portions of the fluid stabilized between strongly adsorbing strips alternated with low-density portions between weakly adsorbing strips. It was the bridge phase that was responsible for the resistance of confined fluid to shear stress.22,23 In our model system, the confining walls are chemically homogeneous and the wall-water interaction cannot be reduced to steric “ball-in-hole” effects. However, the observed resistance to shear can still be explained in terms of the pinning mechanism except that in our case the function of “pins” is served by bridge hydrogen bonds A‚‚‚H-O-H‚‚‚A′, where the force sites A and A′ belong to different walls. More complicated bridges, which involve more than one water molecule, can also be conceived. Now we return to the l ) 3 Å monolayer and consider the effect of shear on the structure and molecular mobility of confined water. As the registry β was increased, the equilibrium number of water molecules between the walls remained unchanged, except at β ) 0.4 and 0.5, when the monolayer lost one of the initial (38) In a simulation study of a simple monoatomic liquid confined in a prototypal slit-like pore to a monolayer thickness, Schoen et al.18 observed deviations from linearity in the d(n) dependence. Unfortunately, in the more complicated case of water, the convergence of d(n) was substantially worse to be able to discern such deviations.
140 Langmuir, Vol. 24, No. 1, 2008
Pertsin and Grunze
crystalline configuration in Figure 10a. It can be seen that the sheared system maintains a long-range regularity in the arrangement of the close-packed rows of water oxygens along the x-axis. By contrast, the regularity in the y direction is lost. When coupled with the simulation results for the x and y components of MSD, this suggests that the shear-induced melting occurs along the shear axis (y), whereas the system retains a long-range order in the x-direction. As to the large oscillations and hops observed in dy at β ) 0.4 and 0.5, they may well originate from collective migrations of the individual segments of the closepacked rows along the y-axis. Such migrations can be promoted by the defects present in the rows (Figure 10b).
4. Conclusions Figure 9. The x and y components of the mean-square displacement of water oxygens as a function of the length of the simulation run, n, for different β at h ) 6.4 Å and l ) 3 Å. (The y-component of MSD at β ) 0.5 is not shown because it practically coincides with the relevant total MSD in Figure 8).
Figure 10. A schematic representation of the in-plane monolayer structure when the confining walls are in registry (a) and out of registry (b). The oxygen atoms in different close-packed rows are shown as black and white only to highlight individual rows.
154 water molecules. As can be judged from the behavior of g(r) at β ) 0.5 (Figure 4), the long-range order weakened, but to a much lesser extent than was observed on going from l ) 3.0 to 2.76 and 3.2 Å. The most interesting was the effect of β on MSD (Figure 8). Up to the yield point at β ) 0.2, the water film remained crystallike, in the sense that MSD remained independent of the length of the simulation chain, n. The only difference from the initial system (β ) 0) was the appearance of a slight anisotropy in the lateral component of MSD, d|: dx ) 0.13 Å,2 dy ) 0.17 Å.2 At β ) 0.3, MSD showed a linear dependence on n, as in the fluidlike films formed at l ) 2.76 and 3.02 Å (Figure 5). In addition, the anisotropy of lateral displacements substantially increased: The slope of dy(n) was larger than that of dx(n) by a factor of 2.5 (Figure 9). Further increase of β to 0.4 led to a drastic change in the behavior of MSD (Figure 8). Aside from a substantial increase in the magnitude of MSD, the dependence d(n) became unstable and very noisy. As seen from Figure 9, the bizarre behavior of d came from dy, whereas dx remained a linear function of n, as in a normal fluid. Note that the increase of β from 0.3 to 0.4 had practically no effect on dx. By contrast, further increase of β to 0.5 returned dx to its “solid-like” behavior, as observed at β e 0.2. This unexpected result was confirmed in a series of independent simulations using different starting configurations and different sequences of random numbers. To understand the origin of the erratic behavior of dy at β ) 0.4 and 0.5, we analyzed snapshots of the relevant monolayer configurations. A typical configuration at β ) 0.5 is shown schematically in Figure 10b in comparison with the ideal
The conclusions following from the above-described simulation results can be formulated as follows: (1) The experimentally observed fluidity of water confined between mica sheets down to a monolayer thickness1-6 can hardly be regarded as a general feature of water. In common with mica, the model walls studied in this work are hydrophilic, but in contrast, they are not charged and do not create a huge concentration of cations in the confined space. The water monolayers confined between the model walls exhibit a welldefined stick-slip transition characteristic of other liquids studied by now. (2) The epitaxial effect of the confining walls does enhance the stick-slip transition in confined water monolayers, but it is not a necessary condition for the transition to occur. The transition also takes place when the lattice of confining walls is incommensurate with the confined water film. (3) The appearance of resistance to shear in confined water monolayers is not necessarily associated with the formation of a crystal-like structure. The confined water films sustain a substantial shear stress while remaining fluid from the viewpoint of translational order and the magnitude of molecular displacements. (4) Fluid water monolayers formed beyond the yield point exhibit unusual features that distinguish them from normal fluids. One is a substantial anisotropy in the in-plane component of MSD, which implies an anisotropy in the coefficient of in-plane self-diffusion. Moreover, when the confining walls are out of registry (β ) 0.5), the model system turns out to be fluid-like in the direction of shear, while retaining a long-range translational order in the perpendicular direction. Last, quite unusual is the erratic behavior of MSD at β g 0.4, which seems to reflect collective motions of water molecules in the direction of shear. Acknowledgment. This research was supported by the Office for Naval Research and the Verband der Chemischen Industrie. We are grateful to Prof. S. Granick for helpful remarks on the manuscript.
Appendix Aside from the above-described simulation work, a few comparative simulations were also performed to get an idea of the dependence of shear stress and normal pressure on temperature, wall rigidity, and the choice of the water model. The comparison was made for a single set of system parameters: l ) 3 Å, h ) 6.4 Å, and β ) 0.2. The effect of temperature was evaluated in a simulation at a temperature of 343 K. To appreciate the effect of wall rigidity assumed in the original model, we considered a flexible wall in which each atom was attached to its ideal lattice site with a spring. The spring constant was adjusted so that the MSD of the wall atom was about 0.2 Å.2 As an
Water Films Confined between Hydrophilic Surfaces
Langmuir, Vol. 24, No. 1, 2008 141
Table 3. Normal Pressure and Shear Stress (kbar) at l ) 3.0 Å, h ) 6.4 Å, and β ) 0.2, Calculated with: (a, b) Original Model at T ) 298 and 343 K, Respectively; (c) Flexible Wall; (d) SPC Water original model 298 K 343 K PN -Tzy
4.5 2.2
5.2 2.3
flexible wall
SPC water
6.8 2.1
6.5 1.2
alternative to the TIP4P water-water potential, the simpler, threesite SPC model39 was used, while the water-wall potential was kept unchanged. The simulation results are presented in Table 3. It can be seen that the increase in temperature and the inclusion of thermal motion for wall atoms lead, as expected, to a noticeable increase in PN. By contrast, Tzy is affected only slightly. The small influence of the thermal vibrations of wall atoms on Tzy seems to result from the competition of two opposite trends. On one hand, the normal component of the thermal vibrations enhances the effective roughness of the wall, which should lead to an increase in the shear stress. On the other hand, the lateral components of the thermal vibrations introduce disorder in the lateral arrangement of the wall atoms, which should weaken the epitaxial effect of the wall and hence the lateral coupling of the wall to confined fluid. (39) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht 1981; p 331.
The replacement of the TIP4P water model by the SPC one also leads to a perceptible increase in PN. This can be associated with the known trend of SPC to underestimate water density at a fixed pressure,24 which should result in an overestimated pressure at a fixed density. Considering that the replacement of TIP4P by SPC kept the density of confined water unchanged (〈N〉) 154 in both cases), the observed increase in PN does not seem unexpected. Unlike the effects of temperature and the thermal motion of wall atoms, the change of the water model gave rise to a substantial change in Tzy. A comparison of the structure of the confined SPC and TIP4P monolayers revealed perceptible differences in the orientational distributions of water molecules. Inasmuch as both simulations used exactly the same water-wall potential, the observed differences originate solely from water-water interactions, more precisely, from the orientation-dependent, electrostatic part of the water-water potential.40 In this respect, the TIP4P model is more adequate. Due to the shift of the negative charge off the oxygen, it better reproduces the molecular quadrupole moment, as well as the water structure and properties.24,41 LA702209G
(40) The Lennard-Jones O‚‚‚O potentials of the TIP4P and SPC models are practically identical. (41) Jorgensen, W. L.; Tirado-Rives, J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6665.
