|Home | About | Journals | Submit | Contact Us | Français|
We study the effects of nanoscopic roughness on the interfacial free energy of water confined between solid surfaces. SPC/E water is simulated in confinement between two infinite planar surfaces that differ in their physical topology: one is smooth and the other one is physically rough on a nanometer length scale. The two thermodynamic ensembles considered, with constant pressure either normal or parallel to the walls, correspond to different experimental conditions. We find that molecular-scale surface roughness significantly increases the solid-liquid interfacial free energy compared to the smooth surface. For our surfaces with a water-wall interaction energy minimum of −1.2 kcal/mol, we observe a transition from a hydrophilic surface to a hydrophobic surface at a roughness amplitude of about 3 Å and a wave length of 11.6 Å, with the interfacial free energy changing sign from negative to positive. In agreement with previous studies of water near hydrophobic surfaces, we find an increase in the isothermal compressibility of water with increasing surface roughness. Interestingly, average measures of the water density and hydrogen-bond number do not contain distinct signatures of increased hydrophobicity. In contrast, a local analysis indicates transient dewetting of water in the valleys of the rough surface, together with a significant loss of hydrogen bonds, and a change in the dipole orientation toward the surface. These microscopic changes in the density, hydrogen bonding, and water orientation contribute to the large increase in the interfacial free energy, and the change from a hydrophilic to a hydrophobic character of the surface.
Water plays a central role in many of the biomolecular self-assembly processes, including the folding of proteins and the formation of lipid membranes.1–10 In key steps of these biomolecular processes, water is often highly confined,11 reduced for instance to a few layers of water molecules between the extended surfaces of large macromolecules. Considering the diversity in the chemical and physical topology of these surfaces, the interfacial behavior of water in a cellular environment will depend sensitively on the details of the confining surfaces. The wetting behavior of surfaces and its dependence on the physical and chemical nanostructure are also technologically relevant, for instance in the development of low-friction fluid flow channels.11–13
Quantitative experimental studies of water at interfaces are highly challenging,14–17 and roughness can be a relevant parameter.18 Recent advances in both theory and experiments9 have re-invigorated the interest in the molecular underpinning of the hyydrophobic effect and the behavior of water near surfaces. Water confined between a “Janus interface” of adjoining hydrophobic and hydrophilic surfaces was found to fluctuate significantly during shear deformations.19 This observation raises several interesting questions about water near heterogeneous surfaces and how local patchiness (wetting versus nonwetting) may affect hydrophobicity and microscopic density fluctuations.20 To find answers to at least some of these questions, molecular simulations can complement the laboratory experiments because in the simulations surface details can be controlled precisely and molecular information can be obtained directly.
Previous theoretical and simulation studies have typically focused on the structure, thermodynamics, and dynamics of water confined between idealized smooth surfaces,21–32 between atomistic surfaces,33–37 within carbon nanotubes,12 and between realistic protein-like surfaces.38,39 These studies have vastly improved our molecular-level understanding of the expected changes in the water behavior due to confinement and the presence of hydrophilic versus hydrophobic surfaces.
In a laboratory experiment, one way to characterize surfaces as hydrophilic or hydrophobic is to measure the macroscopic contact angle of small water droplets on these surfaces. Recent molecular simulation studies have used a microscopic analog of the macroscopic contact angle to define the hydrophobicity of surfaces.40–42 Giovambattista et al.41 used contact angle data to show that that water behavior can be a non-trivial function of the surface polarity. Godawat et al.42 studied the behavior of water near self-assembled monolayer structures with a wide range of chemistries and showed that the hydrophobicities of such surfaces as measured by contact angles can be characterized well by measures related to the water density fluctuations. We similarly saw enhanced water-density fluctuations near extended non-polar surfaces,43 associated with transitions between locally wet and dry regions and resulting in a broadened liquid-vapor like interfacial density profile.
It is well known in the surface wetting literature that the physical roughness of the surface can cause the contact angle to change significantly.44 This change is typically described by either the Cassie45 or the Wenzel46 model depending on the wetting properties of the smooth surface. These models are widely used to describe the effects of the mesoscopic (micron-scale) roughness of superhydrophobic surfaces, such as Lotus leafs. Much less is known on how surface roughness of the order of few molecular diameters,47–49 which is relevant for instance to characterize protein surfaces, will affect the interfacial behavior of water.
This poor understanding of surface roughness effects on the interfacial thermodynamics also severely limits the applicability of surface-area based solvation models. Many widely used models for ligand binding and self-assembly rely on the interfacial free energy of surfaces to predict their solvation free energy.50 However, at the molecular scale corrections should be made to the interfacial free energy to account for curvature and physical roughness.51,52 Examples of such curvature corrections rely on the macroscopic concept of the Tolman length,53 or on the microscopic observation of curvature dependent “cavity expulsion” of water.54
To explore how molecular-scale physical roughness affects the interfacial free energy, we perform molecular dynamics simulations of water confined between walls. One of the major difficulties in performing simulations of confined water is to define the thermodynamic state of the system and to make sure that water is present in a thermodynamically stable equilibrium. Here, we consider two thermodynamic ensembles, defined by constant normal and parallel pressures (Fig. 1) at a constant temperature and with a constant number of particles. We also monitor individual components of the pressure tensor to make sure that the system is in a stable state during the entire simulation time.
In our model system, the wall roughness can be changed in a controlled manner by modulating the local wall position with periodic functions of given amplitudes and wavelengths (see Fig. 2). We find that increasing wall roughness, as characterized by increasing roughness amplitude or decreasing roughness wavelength, can significantly increase the solid-liquid interfacial free energy. Roughness can even change the sign of the interfacial free energy from negative (hydrophilic) to positive (hydrophobic). We examine the microscopic origins of this change in wetting behavior in the context of hydration structure near rough walls and associated density and hydrogen bond distributions.
We use a system of Nw = 1500 SPC/E water55 molecules in a rectangular box with periodic boundary conditions applied in two directions x and y, and fixed boundary condition applied in the z-direction. The walls are represented by a coarse-grained 9-3 Lennard-Jones potential,
which is frequently used to represent interactions between a fluid particle and a wall made up of 12-6 Lennard-Jones particles. Here, ε and σ are the water-wall interaction energy and distance parameters, respectively. z is the distance between the water oxygen atom and the wall surface. We use ε = 1.15 kcal/mol and σ = 0.346 nm in this work. The resulting well depth of the wall-water interaction is ≈ −1.21 kcal/mol. To represent a rough wall as shown in Figure 2, we use a periodic function of the following form,
where z(x, y) represents the wall position as a function of x- and y-coordinates. a is the roughness amplitude, and λ is the roughness wavelength, which is an integer fraction of the simulation box length in the x and y direction. The resulting interaction V [zi − z(xi, yi)] of an oxygen molecule and the wall depends on all three oxygen coordinates (xi, yi, zi).
The equilibrium behavior of our confined water system is obtained from molecular dynamics simulations using LAMMPS56 in a canonical ensemble. The temperature is held constant at 300 K by using a Nosé-Hoover thermostat with a time constant of 0.5 ps. The simulation time step is 2 fs. For each state we simulate for a total time of 4 ns. The initial 2 ns are discarded as equilibration time, and the remaining 2 ns of the simulation runs are used for the final analysis. Electrostatic interactions are calculated with particle-particle particle-mesh solver. To account for non-periodicity in the z-direction we extend the simulation box to three times the height H and apply the correction proposed by Yeh and Berkowitz.57
We calculate the solid-liquid interfacial free energy between the confined water and confining walls from the difference between normal and parallel pressures,58
where is the average of the z-dependent parallel pressure pP (z). γl and γu are the interfacial free energies per unit area (of the surface projected onto the xy-plane) between confined water and the lower and upper wall, respectively, assuming separability.
To calculate how the wall roughness may affect the interfacial free energy for the confined water system shown in Fig. 2, we first calculate the pressure equations of state as a function of wall distance H. The equations of state differ for the two ensembles under constant parallel and normal pressures. As our starting smooth wall system (a = 0), we use wall-oxygen interaction parameter ε = 1.15 kcal/mol and confinement length H = 16.75 Å to obtain a system close to a stable thermodynamic equilibrium and with about 4 water layers between the walls. To explore the effects of roughness, we then modulate the wall location of the lower surface according to the periodic function function in eqs. 1 and 2. Results for the interfacial free energy will be compared in states defined by conditions of 1 atm pressure in either the normal direction (PN) or the parallel direction (PP) (see Fig. 1).
Figure 3 (left, insets) shows the normal pressure PN and the parallel pressure PP as a function of H for roughness amplitudes a between 0 and 10 Å (from left to right) of the lower wall. The simulation pressures are fitted accurately by a cubic polynomial in H. The fitted equations of state allow us to find the H values for which PN or PP is equal to 1 atm. We then use the respective H value to evaluate the pressure component that is not kept constant, i.e., PP (PN = 1 atm) and PN (PP = 1 atm). From the cubic fits, we also calculate the isothermal compressibility κT = −H−1(dH/dP) in both the normal and the parallel direction. The results are shown in Fig. 3 (right, bottom) where dP/dH is evaluated at PN = 1 atm and PP = 1 atm, respectively. We find that κT increases with increasing a for both constant PN and PP. The confined water is thus more compressible near rough surfaces. We want to point out that this is the overall behavior of the confined system and not the local behavior near a surface as studied previously by Godawat et al.42 and by us.43 A similar increase in κT for confined water system was found by Giovambattista et al. 41 with decreasing surface polarity and therefore increasing surface hydrophobicity. Also, compressibility is higher for a system at constant PN as compared with at constant PP for a given a.
Thermodynamic stability is a possible concern in simulations of confined systems. For roughness amplitudes a ≥ 7.5 Å, we find that the cubic fits to the P − H data produce multiple roots for P (H) = 1 atm. For a ≥ 7.5 Å, we thus expect phase co-existence (in a large system) between a liquid and vapor phase, or fluctuations (in a small system) between the two “phases.” For these large roughness amplitudes a ≥ 7.5 Å, the systems simulated are thus likely metastable.
The simulation systems can also be subject to mechanical instabilities. With one of the pressure components, PN or PP, fixed at 1 atm, the other one can be negative. As shown in Fig. 3 (left panels), PP (PN = 1 atm) is negative for a ≥ 3 Å, and PN (PP = 1 atm) is negative for a ≤ 1.5 Å. To check whether the system underwent phase separation (through the formation of nanoscale bubbles), we track the full pressure tensor, making sure that the in-plane diagonal pressure components agree, Pxx = Pyy, and that the normal pressures on the lower and upper walls are equal, Pzz,l = Pzz,u. We do not find any signatures of such deviations for all the state points considered in this work, in particular not for conditions close to 1 atm pressure in either normal or parallel direction. Also, the slope of the pressure equations of state (Fig. 3, left) is always negative, P/H < 0, which means that the thermodynamic stability condition is satisfied (κT > 0; Fig. 3 right bottom). These observations together provide a sufficient reason to believe that our systems of interest are in a thermodynamic equilibrium (or at least in a metastable equilibrium) during the entire simulation times.
Figure 4 shows the interfacial free energy γ as a function of wall roughness amplitude a (left) and wavelength λ (right) under the constraints of PN or PP being equal to 1 atm. We find that γ changes sign with increasing wall roughness by either increasing a (with λ = 11.6 Å fixed) or decreasing λ (with a = 3 Å fixed). This change in sign means that the hydrophilic smooth surface (γ < 0) becomes hydrophobic (γ > 0) with increasing molecular scale roughness. We also find that γ has a rather strong dependence on a and this dependence is stronger for systems under the constraint of constant parallel pressure PP. The dependence of γ on the wavelength λ for a fixed a = 3 Å is weaker than the dependence on the amplitude a (Fig. 4, right). With increasing λ, the calculated γ values drop toward the smooth-wall limit expected for λ → ∞, albeit slowly. Importantly, physical roughness of the surfaces on the length scale of a single water molecule can change the wetting behavior of confined water qualitatively from hydrophilic (γ < 0) to hydrophobic (γ > 0).
If we assume that the interfacial free energy of the smooth wall, γu, is independent of the roughness of the lower wall, we can evaluate the interfacial free energy for the rough wall by subtracting the value for the smooth wall from the total free energy, γl(a) ≈ γ (a) − γ (a = 0)/2. We find that the interfacial free energy for the rough wall calculated in this way varies from −8.1 mN/m to +89.1 mN/m for PP = 1 atm and from −3.6 mN/m to +41 mN/m for PN = 1 atm.
By assuming a solid-vapor interfacial free energy of γSV ≈ 0 mN/m, we can use Young’s equation based on mechanical equilibrium to find the contact angle θ of SPC/E water on our rough surfaces,
where γSV, γSL, and γLV are the solid-vapor, solid-liquid, and liquid-vapor interfacial tensions, respectively. With a liquid-vapor interfacial tension of γLV = 61.3 mN/m for SPC/E water at 300 K,59 and with γl(a) assumed to be equal to γSL and γSV ≈ 0, we can use eq. 4 to estimate the contact angle. At constant PN = 1 atm, the contact angle changes from 87° for a flat interface to 132° for a = 10 Å and λ = 11.6 Å. For constant parallel pressure, a contact angle of 180° is reached just above a = 6 Å, starting from θ ≈ 82° for a flat surface. For larger amplitudes a ≥ 7.5 Å, the calculated free energy of the liquid-solid interface exceeds that of the liquid-vapor and vapor-solid interface combined. It is thus thermodynamically favorable to separate the fluid phase from the rough surface by an intervening vapor layer, and eq. 4 does not have a solution for the contact angle θ. This metastability is consistent with the equation of state at PP = 1 atm having multiple roots for a ≥ 7.5 Å (Fig. 3). We note, however, that in our model these drastic changes in the wetting behavior of microscopically rough surfaces result from the combined effects of wall roughness and confinement.
To understand the origin of the increased interfacial free energy of confined water near rough surfaces, we calculate various structural measures that characterize the spatial and orientational water distribution. We only perform these analyses for systems under the constraint of constant normal pressure, i.e., PN = 1 atm.
Figure 5(d) shows the water density profiles near the smooth upper and rough lower walls. The water density is plotted as a function of the local (x- and y-dependent) distance from the wall, zu − z and zl − z, respectively, averaged over the entire system. Remarkably, neither at the smooth nor at the rough wall do we find any dramatic changes in the average density of the first water shell. Even for the restrictive confinement of our simulation setup, the presence of a rough surface does not appear to affect the density distribution near the opposite surface. However, as we will discuss below, the local changes in the density induced by roughness are pronounced.
To describe the hydrogen bonding behavior of water, we use energetic criteria.60 A pair of water molecules is defined to be hydrogen bonded if the intermolecular pair energy, which is the sum of Coulombic and van der Waals energies between all the atoms of the two molecules, is below some minimum value. We compute these pair energies for water molecules whose oxygen atoms are within 6 Å. Figure 5(a) shows the resulting pair energy distributions for several roughness amplitudes a, with fixed λ = 11.6 Å. The peak at the lower energy values is due to molecules that are tightly held together, primarily by hydrogen bonds. All the curves are statistically indistinguishable from each other meaning that, on average, the wall roughness does not change the pair energy distribution. From this distribution, we also construct a histogram of the number NHB of hydrogen bonds per water molecule. Below an interaction energy cut-off of −2 kcal/mol, a molecule pair is considered to be hydrogen-bonded. As for the energy distribution, we find no major changes in the distributions P (NHB) of the number of hydrogen bonds, averaged over the entire system, as the roughness amplitude is changed.
The dipolar orientation of water also depends sensitively on the roughness. Figure 5(c) shows the average of cos θ as a function of distance zu − z normal to the walls (z-direction), where θ is the angle between the water dipole and the negative z-axis, and zu is the location of the upper smooth wall. Near the upper smooth wall (zu − z ≈ 0), the dipole-moment vectors are on average lying in the xy-plane, with cos θ ≈ 0, consistent with the early simulations of Lee, McCammon, and Rossky.61 In contrast, the dipoles of water molecules inside the rough grooves (zu − z > 1.5 nm) are on average oriented partially toward the lower rough wall, with cos θ > 0.
Upon closer inspection of individual structures, we find that the water molecules at the tips of the finger-like protrusions penetrating into the valleys tend to form peculiar hydrogen bonding structures. These water molecules are typically connected to the bulk-like water above by two hydrogen bonds, and have one or both of their OH bonds dangling, pointing toward the rough surface. An example of a water molecule accepting two hydrogen bonds and donating none is shown in Fig. 6 (indicated by a solid circle and a light blue arrow).
The locally-averaged number of hydrogen bonds per water molecule is consistent with this picture. With our definition, we have on average four hydrogen bonds in the bulk phase at the center of the slab (Fig. 7 right). Near the flat upper wall, the average number NHB of hydrogen bonds drops to three. In the valleys of the rough interface, NHB drops even further to about 2.5, 2, 1.5, and 1 for a = 1.5, 3, 4.5, and 6 Å.
This dramatic loss of hydrogen bonding interactions in the valleys of the rough surface is accompanied by a local depletion in water density. Whereas the average water density as a function of the distance from the rough surface is only weakly dependent on roughness (Fig. 5(d)), the local water-oxygen density at the center of the valleys is low. As shown in Fig. 7 (left), the water density directly above the center of the valleys drops only gradually toward the wall, in sharp contrast to the steep density profiles at the flat wall and at the center of the protrusions of the rough interface (not shown). These broadened density profiles in the valleys are reminiscent of the liquid-vapor interface and the interface near extended hydrophobic solutes.43
Consistent with the capillary-wave fluctuations in the density near extended hydrophobic solutes,43 simulation snapshots show that some of the valleys are transiently filled by locally dense water “fingers” (solid yellow circle in Fig. 6), whereas other valleys are transiently empty (dashed yellow circle in Fig. 6). This analysis suggests that the valleys in the rough interface are fluctuating between locally “wet” and “dry” states. This local “dewetting” appears to be at least partially caused by the loss in hydrogen bonding in the fingerlike water structures able to penetrate into the valleys. Put together, our results suggest that the combination of loss in hydrogen bond energy, and local transient dewetting of the valleys are important factors contributing to the increase in interfacial free energy γ with molecular-scale surface roughness.
We performed molecular dynamics simulations of SPC/E water confined between two surfaces, one perfectly flat and the other with sinusoidal roughness. Water has strong attractive interactions with both surfaces (well depth of −1.2 kcal/mol in a 9-3 potential). From the difference between normal and parallel pressures, we determine the interfacial free energies. In the limit of infinite surfaces at infinite separations, these interfacial free energies would correspond to the surface tensions. Here, however, the two surfaces are close on a molecular scale, and we expect our interfacial free energies to be affected also by the confinement.
Despite having equal interaction strengths of water with both surfaces, we find that the interfacial free energy is negative only for smooth surfaces, but becomes positive when the roughness amplitude exceeds 3 Å at a wave length of 11.6 Å. This sign reversal implies that the character of the surface changes from hydrophilic to hydrophobic as a result of relatively modest levels of roughness at molecular scales. In terms of contact angles, we find an increase from 82° to just under 180° upon going from a flat interface to a roughness amplitude of a = 6 Å at constant parallel pressure PP = 1 atm.
At first sight, our results may seem counterintuitive: we start out with a surface with a favorable (negative) interfacial free energy, γ (a = 0) < 0. By making the surface rough but maintaining the interaction strength with water as a function of the height z, one in effect increases the surface area accessible to water at fixed projected area. With the apparently favorable interactions between water and the surface, γ (a = 0) < 0, and the increased accessible area, one might naively expect the surface free energy to become even more negative with increasing roughness. However, this simplistic argument fails because it does not take into account that the interfacial water has to distort its hydrogen bond network to follow the bumpy surface. In our model system, the free energy penalty of the distorted hydrogen bond network appears to outweigh the energetic gain from the larger accessible area. As a consequence, we observe partial and transient dewetting of the valleys, and the rough surface changes its character from hydrophilic to hydrophobic.
This strong dependence of the interfacial free energy on the roughness may help in the design of superhydrophobic surfaces,44 and the coating of high-flux nanofluidic and microfluidic channels.11–13 Our results also have implications for the solvation thermodynamics and hydration-driven interactions of macromolecules, including proteins. The surfaces of these molecules are inherently rough on the molecular length scale. Improved surface area models of solvation could possibly be obtained by accounting for the strong effects of molecular roughness on the apparent interfacial free energy.
For protein molecules, the chemical heterogeneity of the surface adds another level of complexity to the solvation calculations.39 MD simulations of water between partially hydroxylated silica surfaces (chemical patterning) showed that the boundary between a hydrophobic and a hydrophilic surface may become blurred because of chemical patchiness.34 Willard and Chandler employed a coarse-grained lattice-gas model to investigate the interface between water-like solvent and chemically heterogeneous surfaces.62 By varying the overall fraction and size of hydrophilic sites on a surface, they found that the interfacial fluctuations are large and spatially heterogeneous. Further, adding small uniform solvent-surface attractive interactions were found to bring the interface closer to the surface while maintaining large density fluctuations close to the dewetting conditions.62
An important open question is the degree to which confinement affects the interfacial free energies. To explore this question, simulations of larger systems can be performed in which the separation of the two walls is increased. It will also be interesting to see if the contact angle inferred from γ are consistent with simulations of water droplets on the rough surfaces.
J.M. would like to thank Dr. Vincent Shen (NIST) for helpful discussions. J.M. was supported by a start-up grant from the Lehigh University. G.H. was supported by the Intramural Research Program of the NIH, NIDDK. This study utilized the high-performance computational capabilities of the Biowulf PC/Linux cluster at the National Institutes of Health, Bethesda, MD (http://biowulf.nih.gov).