Article pubs.acs.org/JPCC
Viscosity of Ultrathin Water Films Confined between Aluminol Surfaces of Kaolinite: Ab Initio Simulations Peter J. Feibelman* Sandia National Laboratories, Albuquerque, New Mexico 87185-1415, United States ABSTRACT: Ab initio molecular dynamics simulations of water confined between kaolinite walls were conducted in an effort to make contact with experiments implying dramatic viscosity enhancement for water in nanometer-scale, hydrophilic channels. An earlier ground-state structural optimization of a single water layer on a flat kaolinite(001) surface had yielded a molecular arrangement, which, by appearing hydrophobic to subsequent layers, suggested the possibility of very low flow resistance. Well above the freezing point, however, and under shear, the surface became hydrophilic, as a percentage of first-layer water molecules flipped to expose dangling hydroxyls. This led to simulated steady-state velocity fields consistent with a “no-slip” boundary condition and viscous flow. The magnitude of the viscosity derived from the simulations, only a few times that of bulk water, does not lend theoretical weight to the notion of dramatic enhancement under nanoconfinement.
I. INTRODUCTION Transport of water in narrow channels is of great importance in geophysics, in biological and artificial membranes, and in sensor technology. Gaining a grasp of how confining walls’ atomic arrangements, morphologies, and chemical identities affect water flow can thus be expected to have broad impact. A review of the experimental literature reveals a wide range of shear viscosities reported for water in narrow channels.1 The impetus for the work reported in this article is the suggestion, based on interfacial force microscope2 (IFM) experiments, that in channels a few nanometers wide, and bounded by hydrophilic surfaces, water behaves as though its viscosity were roughly a million times enhanced relative to that of bulk water.3−8 Is this surprising inference credible? Can it be reconciled with surface balance and surface forces apparatus studies of water between mica walls,9−12 in which no such enhancement has been observed, or with MD simulations of water narrowly confined by walls of amorphous silica functionalized with hydrophilic self-assembled monolayers?13,14 Undertaking a theoretical effort to answer these questions unavoidably raises scientific and technical issues: Is viscosity indeed, is hydrodynamicsa meaningful way to describe flow in a channel only a few water layers thick? Do classical force fields describe flow in narrow channels faithfully? Would ab initio molecular dynamics simulations represent an improvement? Do classical and ab initio approaches to water transport in narrow channels yield similar or very different results? Experimental results help to justify using hydrodynamic constructs to describe transport in ultranarrow channels. For instance, in an IFM “drainage” experiment, where a tip was caused to approach a sample at various speeds, the force required to expel the water initially between tip and sample was proportional to the speed of approach, as hydrodynamics requires, for tip−sample separations between 1 and 2 nm.15 In addition, a solution of the Navier−Stokes equations for the force vs tip−sample separation correctly described the observed © 2013 American Chemical Society
quadratic increase in resistance to approach as the tip−sample distance diminished.4 The task of accounting theoretically for a million-fold viscosity enhancement appears inherently forgivingat first. Even an order-of-magnitude error resulting from one or another approximation should not be especially serious when what is at issue is an observed increase of 6 orders of magnitude. Nonetheless, there are concerns. A straightforward way to obtain viscosity values from simulations is to constrain theoretical upper and lower confining walls to move in opposite directions, to integrate Newton’s equations until the system comes to steady state, and then to compute the average drag forces that the upper and lower wall atoms experience. The viscosity is obtained by dividing the imposed shear rate into the calculated drag. A technical obstacle to this protocol is that if the walls move too slowly, bringing the system to steady state can involve very many time steps. In the IFM and related experiments, probe tip speeds were typically of order ∼200 nm/s. This is 10−100 million times slower than the 2−20 m/s speeds typically used in classical molecular dynamics simulation of confined water. For example, in their classical study of water confined between kaolinite walls (conducted concurrently with this one), Haria et al. simulated a system containing ∼2000 molecules and still needed 80 million, 1 fs time-steps to obtain a converged viscosity for a 20 m/s shear.16 For 2 m/s, even runs of this length did not provide a reliable estimate of the viscosity of nanoconfined water. Thus, as discussed further at the end of this article, one needs to consider the implications of extrapolating theoretically convenient wall speeds down by factors of 10−7 and smaller still to compare to experiment. Received: December 10, 2012 Revised: March 1, 2013 Published: March 6, 2013 6088
dx.doi.org/10.1021/jp312152h | J. Phys. Chem. C 2013, 117, 6088−6095
The Journal of Physical Chemistry C
Article
surfaces in the simulation were perfectly flat. “No-slip” and roughness are not inevitably related, in other words. In what follows, section II details the numerical methods used in this study: section II.A reviews the structural properties of the kaolinite (001) crystal face and section II.B the DFT approximate methods used. Section III reports the results of the DFT simulations of water viscosity in narrow, kaolinitebounded channels.
Also possibly impeding an explanation of the observed large viscosity enhancements is that given the intermolecular force law adopted, the channel walls must be hydrophilic. But, examples of hydrophilic surfaces in density functional theory17,18 (DFT) calculations are not plentiful. Even Pt(111), a precious metal surface known experimentally to support complete wetting at 0 K, is not wetting according to DFT energy calculations19 based on using the typically satisfactory Perdew, Burke, Ernzerhof (PBE) approximate functional.20 The periodic ground-state water arrangement on Pt(111) was deciphered in ref 19 on the basis of scanning tunneling micrographs and shown to correspond to a DFT/ PBE relative minimum energy. Nonetheless, the DFT/PBE binding of water molecules is stronger in a three-dimensional ice crystal. With that in mind, the ab initio MD simulations reported here were conceived to take advantage of Hu and Michaelides’ discovery that water adopts a 2-dimensional wetting arrangement on kaolinite(001) in the DFT universe defined by use of the PBE functional.21−23 Thus, PBE kaolinite(001) can stand as the basis for an ab initio study of resistance to water shear in a hydrophilic nanochannel. In general, the advantage of ab initio as opposed to classical force-field calculations is their freedom from semiempirical parameters. This should be important in studies of water adjacent to metals, where reliable classical force fields are not commonly available. By contrast, for water sheared between clay surfaces like kaolinite’s, force fields are well developed.24 Thus, the present study can be viewed as groundwork for ab initio shear simulations with more challenging wall materials. Of great interest is whether the compute-intensive ab initio approach produces meaningful results despite simulation cells containing an order of magnitude fewer molecules than used in classical MD, despite integrating for no more than a few picoseconds, as opposed to nanoseconds, and despite the systematic errors inherent in an unavoidably approximate functional like PBE.20 On the last of these concerns, one should understand that despite the imposing notion of an ab initio study, the approximate functionals used in current-day DFT offer far from a perfect description of water behavior. For a cogent example, although PBE is widely used and broadly satisfactory in studies of the structure of condensed matter, PBE water is structured (glassy) at 300 K,25 where it should easily flow. One must therefore set the computational thermostat to ∼350 K as a way of estimating what happens at room temperature in nature. [This would be disastrous for studies aimed at biochannels, needless to say.] A main outcome of the simulations reported here is a kaolinite-confined water viscosity enhanced by less than an order of magnitude relative to bulk water’s. This is far short of a million but agrees with the small enhancement also found in the classical simulations of ref 16. A second result of the present study is insight into the “noslip” boundary condition typically imposed in the hydrodynamic interpretation of flow between hydrophilic walls. Although there are no “dangling” OH’s in the PBE first-layer structure of T = 0 K water on kaolinite,21,22 at temperature and under shear a percentage are, and H-bonding between them and adjacent water molecules imposes “no-slip.” That is, the water molecules immediately adjacent to the walls move with them. This illustration of the physical basis of the commonly used boundary condition is noteworthy in that the kaolinite
II. NUMERICAL METHODS A. Kaolinite Channel Walls. IFM observations of enhanced water viscosity under confinement refer to motion of a probe tip relative to a sample, both hydrophilically functionalized. At a minimum, to simulate such observations using DFT, it is necessary to identify surfaces that are hydrophilic in the “DFT universe” corresponding to the choice of a particular approximate functional. As noted in the previous section, this is not as easy as it may appear. For example, water has been known for a decade and a half to form a periodic, 2dimensional layer on the Pt(111) crystal face.26 And recently, STM imaging as well as other experiments have led to a molecular-level picture of this wetting layer.19,27,28 DFT calculations confirmed that the proposed molecular arrangement corresponds to a relative minimum energy19 and that it accounts for corresponding (measured) vibration spectra27 as well as the topography observed through thermally programmed desorption of Kr atoms deposited on it.28 Still, the DFT/PBE binding energy/H2O of the Pt wetting structure optimizes to a value of only 0.57 eV,28 some 0.07 eV less than the PBE water molecule’s binding energy in bulk ice.29 Thus, thermodynamically, at 0 K in the DFT/PBE universe, water should not adsorb as a complete wetting layer on Pt(111) but should form 3-dimensional ice clusters instead. Given this situation, Hu and Michaelides’ discovery that water forms a 0 K, 2-dimensional wetting layer on kaolinite(001) in the DFT/PBE world makes this material useful for an effort to simulate the viscosity of water confined between hydrophilic surfaces. The source of kaolinite(001)’s attractiveness to water is its intrinsic, dangling, and easily rotated OH ligands, arranged hexagonally, with a spacing close to commensurate with the basal plane of ice Ih.17 Bulk kaolinite is a stack of silica−alumina bilayers (see Figure 1). The OH ligands hanging off the Al atoms of the nth bilayer H-bond to the O atoms of the silica in the n + 1st. Cleaving kaolinite(001) exposes the O−H ligands of the last alumina layer. Upon water deposition, they form hydrogen bonds to the
Figure 1. Two bilayers of a bulk kaolinite crystal. Gold, gray, red, and white spheres represent Si, Al, O, and H atoms. Note that the nth bilayer is H-bonded to the n + 1st by the OH hanging off its alumina layer. 6089
dx.doi.org/10.1021/jp312152h | J. Phys. Chem. C 2013, 117, 6088−6095
The Journal of Physical Chemistry C
Article
Si atoms of upper and lower films to displace in the −z and +z directions, continuing to integrate the Newton’s equations of motion and thus preserving, to a degree, information already gained on the approach to steady state. Because PBE water is glassy even at temperatures where it ought to melt,25 I ran the ab initio simulations at 350 K, maintaining that temperature by use of a Nosé−Hoover (N− H) thermostat.38,39 Narrow confinement of the water prevents evaporation. The thermostat amounts to an added fictitious mass in the dynamical equations, which by dissipating energy allows the real particles to equilibrate. In the thermostat’s original version, the N−H mass couples to all three components of the real-particle velocities. In a reformulation introduced for the sake of evaluating shear response,40 I allowed the N−H mass to couple only to the y-components of particle velocities, where the z-coordinate runs across the waterfilled channel, and the shear is imposed in the x-direction. Thus, neither the calculated shear rate, d⟨vx⟩/dz, nor the velocity of the fluid normal to the channel walls, ⟨vz⟩, is directly influenced by the thermostat. Energy is dissipated in the plane normal to the channel walls and transverse to the imposed shear. Because this mathematical artifice could introduce anisotropy into the fluid behavior, it is important to confirm (see below) that interatomic forces restore isotropy, allowing energy equipartition to occur. Besides checking equipartition, the drag forces on the upper and lower channel walls, ⟨Fx(upper)⟩ and ⟨Fx(lower)⟩, must be computed to evaluate a viscosity. With the upper wall moving in the +x direction, ⟨Fx(upper)⟩ should be negative and ⟨Fx(lower)⟩ equal to it in magnitude and opposite in sign. The shear viscosity, γ, is then given by
O atoms of the H2O molecules, while what might have been dangling O−H ligands of the waters form hydrogen bonds to the O atoms of the alumina surface. The result is that with a layer of water on the surface there are no dangling O−H’s and the binding is strong. Whether or not this arrangement renders the water-covered surface hydrophobic at T = 0 K, at finite T and with other water molecules present, as in a water-filled channel, one can expect a fraction of the first-layer H2O’s to flip O−H’s and bond to neighboring water molecules instead of the crystal surface. Thus, the first water layer might remain anchored to the kaolinite wall, while at the same time binding to a fraction of adjacent water molecules and thus enforcing a no-slip boundary condition. Accordingly, kaolinite(001) walls oriented with O− H’s dangling into the channel from top and bottom are a reasonable choice for DFT simulations of the viscosity of narrowly confined water. B. Details of the ab Initio Simulations. DFT MD simulations were conducted using the VASP plane-wave code,30−33 the Perdew−Burke−Ernzerhof (PBE) version of the generalized gradient approximation,20 and the projectoraugmented wave treatment of electron−nuclear interactions.34,35 As of May 2011, VASP makes available a suite of constrained dynamics options particularly convenient for the present study.36 One is the SHAKE algorithm.37 I took advantage of it by requiring all nearest O−D separations to equal 1 Å, leaving the D−O−D bond angle free to bend. (Substituting D for H slows the bending vibration, making Newton’s equations somewhat less stiff.) A second helpful option is to force certain atoms to displace by a specified vector in each time step. That makes it easy to drag channel walls at a fixed velocity. At the beginning of each simulation step, a subset of atoms of the kaolinite channel walls is moved by the desired wall velocity times the integration timestep, Δt. VASP’s displacement constraint is also convenient for “graceful” channel expansion or compression, i.e., changing the channel width with minimal impact on molecular density and velocity distributions. The kaolinite wall on either side of the water-filled channel was represented by eight surface unit cells (4 × 2) of a single silica−alumina bilayer, with the alumina sheet and its dangling hydroxyls facing inward. Lattice repeat distances of lx = 20.784 and ly = 18.042 Å were taken from a DFT/PBE bulk kaolinite optimization. I populated this cell with 160 or 224 D2O molecules, initially in layers of 32, corresponding to the static, T = 0 wetting arrangement optimized in refs 21 and 22. To achieve reasonable numerical accuracy given short O−H bonds, I used a 400 eV plane-wave basis cutoff. The surface Brillouin zone was sampled with a single point at Γ̅ . Although integrating the Newton’s equations of a manyparticle system in contact with a thermal reservoir must eventually lead to equilibrium under static conditions and to steady-state under a constant driving force, efficient use of computer resources means starting as close as possible to the thermal end results. For a channel whose upper and lower boundaries are essentially identical, I therefore started from an arrangement with equal numbers of up and down pointing O− H bonds. To impose shear on the water layer, the Si atoms of the upper wall were forced to move at a constant speed in the +x direction, and of the lower wall in the −x direction, while their relative positions were fixed at values drawn from the bulk kaolinite optimization. To compress the water film, I forced the
γ = [⟨Fx(lower)⟩ − ⟨Fx(upper)⟩]/[A yz d⟨vx⟩/dz]
(1)
where Ayz is the area of the channel unit cell perpendicular to the direction of drag (that is, to the x-direction) and d⟨vx⟩/dz is the shear rate. For this equation to make sense, ⟨vx⟩ should be a linear function of z. Because integration times in the calculations were never more than a few picoseconds, it was necessary to verify that the computed average drag forces on the walls had achieved a good approximation to steady state when the integration stopped. For that purpose, I recomputed forces dropping the contributions from the first 1 to 3 ps of integration to see how that changed the results. Such comparisons yielded error estimates of 45° or 45° to the Kaolinite Wall, and the Percentage Oriented at an Angle 45°
% with an angle
