Hydrogen Bonds and Vibrations of Water on - American Chemical

Jul 8, 2009 - David J. Wesolowski,. ⊥. David Cole,. ⊥ and Jorge O. Sofo*,†,#. Department of Physics, Department of Geosciences, and Materials Re...
0 downloads 0 Views 1MB Size
13732

J. Phys. Chem. C 2009, 113, 13732–13740

Hydrogen Bonds and Vibrations of Water on (110) Rutile Nitin Kumar,† Sanghamitra Neogi,† Paul R. C. Kent,‡ Andrei V. Bandura,§ James D. Kubicki,| David J. Wesolowski,⊥ David Cole,⊥ and Jorge O. Sofo*,†,# Department of Physics, Department of Geosciences, and Materials Research Institute, The PennsylVania State UniVersity, UniVersity Park, PennsylVania 16802, Center for Nanophase Materials Sciences and Chemical Sciences DiVision, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37830, and St. Petersburg State UniVersity, St. Petersburg 198504, Russia ReceiVed: February 23, 2009; ReVised Manuscript ReceiVed: May 15, 2009

We study the relation between the hydrogen bonding and the vibrational frequency spectra of water on the (110) surface of rutile (R-TiO2) with three structural layers of adsorbed water. Using ab initio molecular dynamics simulations at 280, 300, and 320 K, we find strong, crystallographically controlled adsorption sites, in general agreement with synchrotron X-ray and classical molecular dynamics simulations. We demonstrate that these sites are produced by strong hydrogen bonds formed between the surface oxygen atoms and the sorbed water molecules. The strength of these bonds is manifested by substantial broadening of the stretching mode vibrational band. The overall vibrational spectrum obtained from our simulations is in good agreement with inelastic neutron scattering experiments. We correlate the vibrational spectrum with different bonds at the surface to transform these vibrational measurements into a spectroscopy of surface interactions. 1. Introduction 1

Langmuir proved more than 90 years ago that diverse insulating substances will have between two and five layers of water on the surface when exposed to ambient air. These layers are mainly held by hydrogen bonds because the interaction between the electric field generated by the insulating surface and the dipole moment of the water molecule decays rapidly away from the surface. Ewing and co-workers estimated that only the first layer of water molecules is strongly affected by the electric field of the substrate.2 As a consequence, hydrogen bonds not only hold water molecules together but also are a good part of the interaction of water with surfaces beyond the first layer. There are technological advances that depend on our scientific understanding of the static and dynamical aspects of hydrogen bonding at surfaces. These applications range from the understanding of biological membranes to fluid transport in geologic media to self-cleaning windows. A relatively low formation energy, ranging from 2 to 10 kcal/mol, makes the formation and breaking of these bonds a possible process at normal temperatures.3 As a consequence, the study of structure and stability of systems held together by this type of bond requires careful consideration of entropic effects associated with their dynamics. For the computational study of these phenomena, this characteristic is a challenge and an advantage at the same time. Dynamical simulations are computationally intensive and require unique facilities. However, at the same time, they offer unique insight in a very relevant subject. To study these effects we will exploit the relation between hydrogen bond strength and hydrogen vibrational frequencies associated with stretching modes of water. In particular, we have chosen the (110) surface * Corresponding author. E-mail: [email protected]. † Department of Physics, The Pennsylvania State University. ‡ Center for Nanophase Materials Sciences, Oak Ridge National Laboratory. § St. Petersburg State University. | Department of Geosciences, The Pennsylvania State University. ⊥ Chemical Sciences Division, Oak Ridge National Laboratory. # Materials Research Institute, The Pennsylvania State University.

of rutile (R-TiO2) for our study, because of the extensive availability of experimental information on this system, the technological importance of the material, and the fascinating richness of the effects of water interaction with this surface. The surface of crystalline rutile (R-TiO2, tetragonal crystal structure) is one of the most investigated metal-oxide surfaces.4,5 However, it still remains the source of many controversies, and many of these are related with its interaction with water. For example, a fascinating photoinduced superhydrophilic effect was discovered in 19976 with important implications to the fabrication of self-cleaning surfaces. However, more than ten years later, the microscopic origin of this effect is a matter of debate.7,8 The strength of the interaction between water and rutile is at the center of this discussion. At a microscopic level, this interaction is the result of subtle forces produced by hydrogen bonds and dynamical effects related to vibrational modes and their associated entropies.9 One of the essential unsolved dilemmas in this area is the dissociation state of water on the surface. The molecular or associated form of adsorption corresponds to the formation of a chemical bond between the oxygen of a water molecule and the undercoordinated Ti(IV) atom at the surface (referred to as the terminal oxygen (TO) site). The so-called BO (a structural surface oxygen bonded to two fully coordinated Ti atoms in the (110) surface) can also form a covalent bond with a hydrogen atom and can induce the dissociation of the absorbed water molecule. Many static studies of the geometry and configuration of this chemisorbed monolayer have been done in the past.10-14 Thus, the two extremes of (110) surface termination, which have been explored by static density functional theory (DFT) and classical molecular dynamics (CMD) simulations with force-field parameters optimized on the basis of DFT results,13-21 are that (a) all bridging and terminal sites are hydroxyl groups (the dissociative surface) and (b) the surface is terminated by chemisorbed terminal water molecules and unhydroxylated BO atoms (the associative surface). In fact, pH-dependent protonation of these surface sites can be directly linked to the surface charge in aqueous media.22

10.1021/jp901665e CCC: $40.75  2009 American Chemical Society Published on Web 07/08/2009

H-Bonds and Vibrations of Water on (110) Rutile To make things more interesting, in real samples there is always a hydration layer of structurally perturbed and dynamically hindered water molecules covering these adsorbed molecules.23 The dissociation state is not static. It is the result of a dynamical equilibrium dominated by hydrogen bonding among the surface atoms, the adsorbed water molecules, and the hydration layer. Protons are shuttled at the surface via the Grotthuss mechanism with an efficiency that brings the surface to its dynamical equilibrium. The dynamical aspect of the equilibrium between dissociative and associative absorption is at the core of the long-standing discrepancy both between theory and experiment and between different theories about the amount of water that is dissociated at the surface.9,12,24-27 This long controversy is evidence to the fact that the energies involved are all small and of the order of the energy of different arrangements of the hydrogen bond network around the adsorption sites.28 We study these interactions at an atomistic level with an emphasis on the dynamics of the system by relating the strength of the interactions to the vibrational density of states (VDOS) of the system. This provides both a tool to analyze our quantum mechanical molecular dynamics simulations and a method to relate our analysis to measurable quantities accessible in current experiments. The interplay among hydrogen bonds, dissociation, and vibrational frequencies is a local effect that dominates the macroscopic behavior of the material. At the local level, the stretching mode frequency of an OH bond depends strongly on the surrounding environment. The red shift of the OH vibrational frequency upon the formation of hydrogen bonds is a wellknown effect. When isolated, the vibrational frequency of this mode appears as a well-defined peak at about 105 THz (3500 cm-1 or 434 meV). In the liquid, the proximity to oxygen atoms from other molecules creates a strong attraction force that weakens this covalent bond. This hydrogen bonding with other water molecules lowers the vibrational frequencies of the stretching band; this band can extend as low as 85 THz (2800 cm-1 or 350 meV).2 The strong attraction produced by hydrogen bonds, evidenced by the lower vibrational frequency, leads the way to the dynamical process of water dissociation and proton transfer in acidic environments. The so-called Grotthuss mechanism can be seen as the diffusion of an excitation on the hydrogen bond network. A surface with exposed oxygen atoms in different bonding situations, such as (110) rutile, can generate hydrogen bond centers of diverse strength, including stronger bonds compared with those offered in liquid water. All of these subtle differences are reflected in the vibrational energies of the system and can be measured by inelastic neutron scattering (INS) due to the large incoherent scattering cross section of hydrogen.29 Through our simulation we are able to deconstruct the contributions of different atoms to the vibrational spectrum and elucidate the local interactions at play in this system. To explore the above-mentioned effects in a computer simulation that combines a realistic description of the interactions with the dynamical aspects of the water coverage, we performed extended-time ab initio molecular dynamics (aiMD) simulations of the system consisting of three water molecules per Ti2O4 surface unit on the (110) surface of rutile at three different temperatures. The trajectories obtainedsfor 288 atoms and more than 15 pssare used to calculate the VDOS of the system. We focus on the water contribution to the VDOS. This contribution is dominant at frequencies between 10 and 120 THz. The deconvolution of the VDOS into contributions of individual atoms provides a unique tool to study the local hydrogen bond dynamics at the surface and analyze the signature

J. Phys. Chem. C, Vol. 113, No. 31, 2009 13733

Figure 1. Left panel: the simulation cell. Oxygen is depicted in red, hydrogen in white, and titanium in gray. The light blue dashed line marks hydrogen bonds. Right panel: a histogram of the position of oxygen atoms (red) and hydrogen atoms (white) along the direction perpendicular to the (110) surface. F (z) corresponds to the distribution function of the oxygen atom along the z direction of the cell. This histogram was averaged over 10 ps of simulation at 300 K. BO denotes the layer of surface BO atoms. L1 corresponds to the water layer on top of the undercoordinated titanium atom. L2 and L3 denote other water layers. Two peaks in the hydrogen distribution are highlighted. H1 is formed by hydrogen atoms bridging between surface oxygen atoms and between hydrogenated BO and L2 (L1-L1, L1-BO, and L2-BH), while H2 are hydrogen atoms bridging L1 and L2 and between BO and L2.

of the dissociation state on measured vibrational properties at the interface. 2. Methods The aiMD simulation cell consists of a periodic slab geometry composed of alternating rutile slabs and vacuum. The rutile slabs are oriented with the surface in the (110) direction and terminated at the bridging oxygen (BO) atoms, resulting in an anhydrous surface that exposes repeating rows of undercoordinated BO atoms protruding above the plane of Ti atoms, fully coordinated surface oxygens (SO) in the Ti plane, and undercoordinated surface Ti(IV) atoms. The dimensions of the cell in the directions parallel to the surface are 11.836 Å along the surface rows and 12.994 Å across these rows (constituting eight Ti2O4 repeating surface units). The slab contains 48 Ti atoms in three planes (two surface planes and a central plane that is fixed at the bulk rutile structure). Perpendicular to the slabs, the repeat length is 24 Å; this allows a gap of approximately 14.7 Å between BO atoms. Numerous ab initio studies of the TiO2 (110) surface showed that a vacuum gap of 6-10 Å is sufficient for reproducing correctly the surface properties in the 3D slab model. Thus, Zhang and Lindan21 used a gap around 7 Å for the study of two-layer water adsorption on TiO2 (110). They concluded that such a gap is sufficient to converge a typical adsorption energy to within 1% with respect to the gap. Taking into account these previous experiences, we consider that the value we use of 14.7 Å is enough for the study of three-layer water adsorption. In this case (as it can be seen in Figure 1), the average distance between the oxygen atoms of third-layer water molecules at the top and bottom of the simulation cell is about 6 Å. This is also sufficient to prevent the strong interaction between molecules of different adsorption layers. Each of the two surfaces of the slab contains eight BO atoms and eight undercoordinated (five-fold coordinated) titanium atoms. On each surface, free of defects, we arrange 24 water molecules, making a total of 48 water molecules in the simulation cell. This amount of water was chosen to mimic as closely as possible the hydration situation determined experimentally for the rutile

13734

J. Phys. Chem. C, Vol. 113, No. 31, 2009

nanoparticles that have been measured in quasielastic neutron scattering experiments and treated by DFT optimized CMD simulations by Mamontov et al.23 The nanoparticles used in their experiments predominantly exposed the (110) surface, and the hydration level obtained by exposing the nanopowders to laboratory air corresponded to about 3.5 H2O molecules per Ti2O4 formula unit at the surface.20 CMD simulations12 and X-ray crystal truncation rod studies48 confirm that the first sorbed water layer consists of terminal water molecules that bond to the highly hydrophilic five-fold coordinated titanium atoms, at sites equivalent to oxygen lattice sites in the bulk rutile structure; we will call the oxygen atom of this chemisorbed molecule the TO, as is usually done in the literature on the subject. The network of hydrogen bonds between water molecules in the second and third layers is rigid enough to make the relaxation times longer than what is achievable by aiMD. Thus, we need to prepare an initial configuration by a more efficient method and feed these already relaxed atomic positions into the aiMD engine. Hence, we use CMD to generate our initial structures. First, we prepare a bulk water sample by running CMD simulations with the extended simple point charge (SPC/E) water potential.30 The bulk water sample was prepared in a periodic cell with x- and y-dimensions being equal to the corresponding dimension of the rutile (110) slab. The z-direction of the periodic cell was determined to accommodate the total number of water molecules included in the system. This cell was run for 100 ps at a temperature of 2000 K to break any possible hydrogen bond network arbitrarily introduced by our initial positions. A temperature ramp going from 2000 to 300 K was run for another 100 ps, and an additional 100 ps was used to thermalize the system at 300 K. According to relaxation times of SPC/E water,31 this procedure should provide a reasonably equilibrated hydrogen bond structure for liquid water. A slab is formed from this bulk water system by removing the periodic boundary conditions in the z-direction. This slab of the previously described bulk water system with its normal parallel to the z-direction was located on top of the rutile (110) surface already “decorated” with the layer of terminal water molecules. With this system we run CMD at 300 K, with the force field of Bandura and Kubicki,15 for times longer than 1 ns until the total energy drift is negligible. For this run we used the General Utility Lattice Program (GULP).32 We use the resulting relaxed and thermalized atomic positions as the input to aiMD. Forces for the aiMD simulations have been calculated within DFT33,34 with a plane wave basis set as implemented in the Vienna ab initio simulation package (VASP).35-38 Core electrons have been treated with a frozen projector augmented wave scheme.39,40 The titanium 1s, 2s, and 2p orbitals are treated as core, hence leaving 12 electrons per atom to be considered as valence electrons. Although our setup uses a low precision setting with the energy cutoff of 212 eV, we did shorter aiMD runs with an energy cutoff of 400 eV that yielded a similar vibrational spectrum. The effect of a higher cutoff on the vibrational spectrum is a reduction of all frequencies by 5%. As we will see below, the conclusions of this work do not depend on this small reduction. Moreover, as this is a rather big periodic cell, we used only one k-point at Γ for sampling the Brillouin zone. The exchange and correlation potential was treated in the generalized gradient approximation of Perdew et al.41,42 Extensive parallel CMD runs with the ab initio forces were done in a normal volume and temperature ensemble where the temperature was kept constant using a Nose-Hoover thermostat as implemented in VASP. The time step was set to 0.5 fs to

Kumar et al. properly sample the high-frequency vibrations of the hydrogen atoms. The equations of motion were integrated using the Verlet algorithm. The aiMD simulations required approximately one week of central processing unit time for each 15 ps trajectory on 256 Cray X1E vector processors or up to 2048 Opteron processors of the Cray XT4 at the National Center for Computational Sciences, Oak Ridge National Laboratory. The trajectories obtained from the aiMD simulations at 280, 300, and 320 K were used to calculate the VDOS of the system. This is the Fourier transform or power spectrum of the velocity autocorrelation function, which is calculated using memory functions determined from a stochastic model of the velocities, which is similar to the maximum entropy method.43 The full method is implemented in the nMoldyn analysis program,44 which permits the decomposition of the VDOS into the contributions from each individual atom in the simulation. We use the standard definition for hydrogen bonds, wherein the oxygen accepting the hydrogen bonds is considered an acceptor oxygen and the oxygen with which the hydrogen is covalently bonded is called a donor oxygen. A hydrogen bond is considered to exist if the distance between the acceptor oxygen and a hydrogen covalently bonded to a donor oxygen is less than 2.3 Å and the angle among the hydrogen, the donor oxygen, and the acceptor oxygen is less than 30°.45,46 However, we have tested other alternative definitions, and our results are not affected by the particular definition selected.47-49 3. Results 3.1. Sorbed Water Structure. A snapshot of the simulation cell is shown on the left panel of Figure 1 along with a profile of the oxygen and hydrogen positions averaged over 10 ps of production time of the aiMD simulation at 300 K on the left panel. The profile shows that water forms a layered structure on the surface of rutile. The first layer is composed of TO atoms narrowly distributed in the z-direction just next to the very narrow peak of the BO atoms. We designate this layer as L1, as consistent with Mamontov et al.23 The TO bond to the surface is very strong because of the undercoordinated state of the metal atom. Its bond length is only about 0.3 Å longer than latticeequivalent bonds in the bulk crystal and shortens by approximately 0.25 Å when the terminal water dissociates to form a terminal OH group. This change contributes to the broadening of the L1 peak that can be used as an indicator of the degree of dissociation of water on the surface. There is a much smaller change in the BO-Ti bond length when this oxygen forms a covalently bonded surface hydroxyl. The water molecules in L1 form hydrogen bonds with adjacent water molecules in L1, with BO atoms, and with water molecules in the next layer, denoted by L2. This bonding situation is reflected in the hydrogen atom profile. Two peaks are labeled in the profile. The peak labeled H1 collects the hydrogen atoms forming H-bonds between adjacent TOs, between TOs and BOs, and between hydrogenated bridging oxygens (BH) and L2 water molecules. Hydrogen atoms participating in the H2 peak form H-bonds between L1 and L2 and between BO and L2. Molecules in L2 are hydrogen-bonded to TO and BO and to water molecules in the outermost layer (L3). Molecules in L3 are hydrogen-bonded to other water molecules in L2 or L3, but there is also a large density of non-hydrogen-bonded hydrogen atoms in L3. The same layered structure at the rutile-water interface was found in CMD simulations of both bulk SPC/E water12 and SPC/E water at approximately three-layer coverage.20 We recorded the position of every oxygen atom in L2 and folded that position to an irreducible quarter of the surface unit

H-Bonds and Vibrations of Water on (110) Rutile

Figure 2. Density plot of the oxygen atom positions in the L2 layer projected in three orthogonal directions. The most prominent peak in the density of water molecules in L2 is located around the BO atoms. The x-z and y-z panels show that these water molecules, which are hydrogen-bonded to the BO atom, form an umbrella-shaped distribution because of the high strength of this bond. The other prominent site corresponds to water molecules bridging between two terminal water molecules at L1. This density plot has been averaged over 10 ps of simulation time. In the density plot, red corresponds to high density, while yellow corresponds to lower density.

cell. The projections of the resulting density plots are shown in Figure 2 for an average of 10 ps of simulation time at 300 K. Other temperatures display similar adsorption patterns. The most prominent location corresponds to the water molecules in L2 hydrogen-bonded to the BO at the surface (Figure 2). The presence of this strong hydrogen bond is evidenced in the umbrella-like density distribution shown in the x-z view. This picture shows that most L2 water molecules are strongly hydrogen-bonded to the BO of the surface, moving around it while keeping the bond distance mainly constant. Another, less visited, adsorption site involves L2 water molecules located above and between two terminal L1 sites. L3 molecules show no identifiable structure relatable to the (110) surface structure, though L3 does not have the hydrogenbonding structure of bulk water. This suggests that the effect of surface bonding and the distortions produced by it on the water structure are basically limited to the first and second adsorbed layers. Of course, in the direction perpendicular to the surface there is “layering” going beyond the second layer. This layering is a simple termination effect seen even in weakly interacting liquids. However, the L3 molecules are clearly bound to the surface at the temperatures investigated (280-320 K), since any molecules that escaped from one surface quickly migrated to the other surface in the CMD simulations used to thermalize the aiMD initial configurations. The ordering produced in the direction parallel to the surface, limited to only the second water layer, is the result of weaker and more subtle interactions between the surface and the water molecules. L1 is formed by rows of water molecules in the terminal position. These rows are parallel to the BO rows. It is clear that many of the surface adsorption properties are influenced by the degree of dissociation of the water molecules adsorbed in L1. The structure of water molecules in L2 and L3 does not escape this influence. This was clearly demonstrated in the CMD simulations of Prˇedota et al.17 Because of the nonreactivity of their force field, those simulations were done at a fixed amount

J. Phys. Chem. C, Vol. 113, No. 31, 2009 13735 of dissociation in L1 that is set at the beginning of the simulation. This allowed them to determine the percentage of dissociation at will and study the adsorption properties for a given fixed dissociation. The adsorption sites extracted from these CMD simulations were found to be different in the case of a fully associated configuration compared with a fully dissociated case. In our aiMD simulations, the amount of dissociation is oscillating in thermodynamic equilibrium. Because the amount of dissociation for hydrated (110) rutile is low (∼25% in our simulations), the patterns observed are closer to the patterns corresponding to the fully associated case of the CMD results. An associated configuration provides two donor hydrogen bonds from the water in L1 and two strong acceptor sites provided by the unprotonated BOs. The dissociated surface, when BOs and TOs are OH groups, only provides two donor links because of hindrance blocking the possible accepting sites. Our density plots are in good agreement with the L2 adsorption sites of the CMD simulations at 300 K, for the case of no dissociation.12 However, the L2 absorption site between TO atoms that we see as a secondary site in our simulations is not present in the CMD oxygen density plots.17 In our simulations we see that this site is mainly produced by a water molecule in L2 that accepts hydrogen bonds from the molecules in L1. Water molecules in L1 are strong hydrogen bond donors because the bond with Ti weakens the OH covalent bond. The force field used in Prˇedota et al.’s classical simulations (and in Mamontov et al.23) treats the water molecules in L1 as SPC/E molecules, in the same way as water molecules in the liquid.17 Our aiMD simulations show that, because of bonding with the undercoordinated surface Ti atom, these water molecules have different bonding characteristics than water molecules in the liquid. The bonding between the oxygen and the titanium atoms weakens the OH bonds in the water molecule and donates more charge to the oxygen atom. In spite of this detail, the density contour around the bridging atom shows the same shape in both the CMD and aiMD simulations. X-ray crystal truncation rod experiments done by Zhang et al.50 also indicated two (and only two) well-defined structural layers of sorbed water beyond the BO layer, at the interface between the bulk liquid water and the polished and annealed (110) surfaces of rutile single crystals. These experiments indicate essentially the same L1 configuration as observed in this study, including a shortening of the Ti-O bond when the surface becomes more dissociated because of an increase in the solution pH. The X-ray study also showed that the hydration layer L2 is highly structured with well-defined absorption sites. We see that the position labeled AW1 by Zhang and co-workers, on top of the BO atoms, is at nearly the same vertical position above the Ti surface plane as our main L2 sorption site, though the average lateral position indicated by the experiment is mainly atop the BO, whereas our simulations would suggest that the average position of this sorption site is shifted toward the TO site (Figure 2). In the analysis of the experimental information, two other adsorption sites appear. We do not find in our simulations water molecules located at these sites. The crystal truncation rod experiments are done on the surface of a crystal immersed in a thick layer of water. Our simulations, instead, are designed to mimic the nanoparticles measured in neutron scattering with a much lower level of hydration. This different hydration level may influence the adsorption sites through changes in the dissociation level at the surface or through changes in the hydrogen bond network of the liquid surrounding the surface. Other possible sources of this discrepancy include (a) surface defects and step-edge sorption sites, surface con-

13736

J. Phys. Chem. C, Vol. 113, No. 31, 2009

Figure 3. VDOS as a function of frequency for the three temperatures simulated in this work. In the scale of this plot, mainly hydrogen atom contributions to the VDOS are seen. The peak centered at 25 THz corresponds to librations of the water molecules and that centered at 50 THz to bending modes, and the band of modes centered at 100 THz corresponds to stretching modes of the OH bonds. The sharp peak at 120 THz is produced by vibrations of free hydroxyls. The main temperature effect observed in this range is a slight shift to higher frequencies, as the temperature is increased, of the stretching band and an opposite change in the libration band. Both are the results of a weakening hydrogen bond network.

tamination, and other artifacts associated with the X-ray experiments involving real crystals and aqueous solutions and (b) artifacts in the simulation results associated with the intrinsic approximations (i.e., use of DFT to parametrize the simulations and the treatment of core electrons) as well as boundary conditions imposed to make the simulations practical, including the limited lateral dimension, the thickness of the rutile slab, or the amount of vacuum between surfaces. The cell size we have studied is at the frontier of what is possible today using world-class computing facilities. This interesting and intriguing discrepancy is currently under investigation by our group in collaboration with the experimental group. 3.2. Sorbed Water Dynamics. The hydrogen bond layered structure at the rutile-water interface has a strong influence on the vibrational modes of the water system at the interface. It will be clearly seen from our results that the hydrogen bonds associated with BO and TO sites are in general stronger than hydrogen bonds found in bulk liquid water at equivalent conditions. The higher strength is manifest in the vibrational frequencies of the stretching modes. Mamontov et al.20,49 were able to make a qualitative link between the translational and the rotational dynamics of water sorbed on rutile surfaces on the tens of picoseconds to tens of nanoseconds time scale, obtained from CMD simulations and quasielastic neutron scattering measurements.23 This time scale is mainly inaccessible in aiMD simulations at the system size required in this study. However, as will be shown below, we can make a similar link with INS, which can probe the vibrational dynamics of the system. Calculating the Fourier transform of the velocity autocorrelation function, we obtained the total VDOS shown in Figure 3. The contribution of bonds involving hydrogen atoms is dominant in the energy range shown. Contributions from other bonds (Ti-O) are relevant at frequencies below about 10 THz. The overall shape of the vibrational spectrum is typical of condensed water systems. We observe four characteristic peaks. The peak centered on 25 THz corresponds to a librational movement where the molecules move as rigid bodies under the mutual interaction produced by hydrogen bonding. This mode

Kumar et al. is not present in gaseous systems. A sharp peak around 50 THz is produced by the bending modes of the water molecules. Above 75 THz is the domain of the two stretching modes of water molecules. These modes are strongly affected by the hydrogen bonding between water molecules. The sharp peak at 120 THz is the vibrational frequency in our simulations that corresponds to a free water molecule OH bond with minimal influence from hydrogen-bonded neighboring molecules. These are mainly due to the water molecules in the outer portion of L3, adjacent to the “vapor” space between the two hydrated surfaces. The broad peak between 75 THz and 120 THz corresponds to OH stretching modes softened by the hydrogenbonding interactions between water molecules. We see that as we go from 280 to 320 K, the stretching mode peak shifts to higher frequencies, while the librations peak gets shifted to lower frequencies. These shifts are produced by a decrease in the density of hydrogen bonds in the system as the temperature is increased. In a recent INS study of water sorbed on anatase nanoparticles (another tetragonal polymorph of TiO2), Levchenko et al.29 measured the vibrational spectra of water confined by the surface of anatase nanoparticles. Although anatase is a different crystalline form of titanium, we expect a similar overall behavior of water on its surface, and this is supported by recent INS studies of water sorbed on rutile nanoparticle surfaces with the (110) surface predominant.51 In Figure 2a of their paper,29 the authors show the VDOS extracted from INS, in the range between 150 meV (∼36 THz) and 500 meV (∼121 THz). Note that at the temperature of their experiment (4 K) the spectrum corresponds to that of frozen, amorphous water. However, dominant features of the spectrum are readily relatable to our simulation results. The three peaks we obtain in our simulation for this energy range are clearly identified in the experiment. The bending mode appears in the experiment at 200 meV (∼48 THz) in excellent agreement with the result of our simulation. The experimental figure shows a broad band of the OH stretching modes at 350-470 meV (∼85-114 THz) with a sharp peak at 450 meV (∼109 THz). As noted by the experimental group, the softening of the OH stretching mode is surprisingly large compared with the spectrum of ice or even with liquid water. Our simulations also show this surprising softening of the OH stretching modes as well as the sharp peak at the upper end of the band. By deconvolution of the individual atomic contributions to the spectrum, we can attribute the softening of the band to strong hydrogen bonds between the BO atoms and the water as well as between water molecules in L1 and molecules in L2. We can also determine that the peak at 450 meV in the experimental figure corresponds to stretching vibrations of OH groups not hydrogen-bonded and pointing to the vapor region. As discussed in the Methods section, this peak appears in our simulation at a higher frequency (120 THz). This difference may be due to either the higher temperatures at which the simulations were conducted or to the use of a low energy cutoff in the simulations. Shorter time molecular dynamics simulations with a higher energy cutoff locate this peak at 113 THz, closer to the experimental frequency. Clearly, our molecular dynamics simulations reproduce satisfactorily the observed vibrational spectra. We will now use our simulation results to explore the microscopic explanation of the observations, namely, the large broadening of the stretching band and its relation with hydrogen bonding at the interface. The hydrogen atoms that populate the H1 and H2 peaks in the density profile of Figure 1 were identified above as being associated with hydrogen bonds to BO and TO sites. The group

H-Bonds and Vibrations of Water on (110) Rutile

Figure 4. Contributions to the stretching band of the VDOS from different layers of hydrogen atoms. As indicated in Figure 1, H1 corresponds to atoms that are the closest to the surface. These atoms are forming bonds between L1 oxygen atoms, between L1 and BO, and between BH and L2. H2 is the layer of hydrogen atoms that are forming bonds between L1 and L2 and BO and L2. The group of hydrogen atoms denoted by H3 is the rest, including those at the interface with the vacuum. The picture shows that at all temperatures the lower part of the spectrum is dominated by the atoms in H2 with contributions from H1 and no participation of atoms in H3. This shows that the strongest bonds are formed among surface species. The sharp peak at 120 THz is produced by OH bonds at the interface between water and vacuum included in the H3 group.

H3, not marked in Figure 1, encompasses all other hydrogen atoms in the cell. This group will include those atoms forming typical liquid water hydrogen bonds and non-hydrogen-bonded OH groups at the interface between water and vacuum. The contributions of the H1, H2, and H3 groups to the stretching band of the VDOS are presented in Figure 4 for the three temperatures presented in this paper. The figure shows that H2 dominates the lower frequency end of the band; it includes the contributions to the lowest energy vibrations and correspondingly is where the strongest hydrogen bonds are formed. H1 has also some weight in the lower frequency end of the band; however, its contribution on average is located at slightly higher energies. H3 shows a stretching band contribution similar to that expected from liquid water in addition to some free hydroxyl group oscillations at the interface between water and vacuum. It may seem surprising that H2 shows stronger hydrogen bonds than H1, because we claim that the surface species, such as L1 water molecules and BO atoms, are the strongest hydrogen bond formers in this system and H1 includes those hydrogen atoms forming bonds among them. However, all hydrogen bonds formed by these surface species are weakened because they are forming bonds between species that are also bonded to Ti atoms at the surface and are not free to move. Consequently, the hydrogen bonds are stretched as in the case of bonds between L1 and BO or compressed as in the case of L1-L1 bonds. H2, instead, contains bonds that are formed between a strong hydrogen bond forming center at the surface and water molecules in the liquid that are free to move and adapt to maximize the strength of the hydrogen bond. In this way, H2 dominates the lower end of the spectrum because it contains the strongest hydrogen bonds in the system. The relation between hydrogen bond strength and vibrational frequency of the stretching mode goes beyond the global picture presented so far, and the power of atomistic simulations is to enable elucidation of the microscopic roots of the macroscopically observed behavior. The simulations reveal that there are many different types of hydrogen bonds at this interface. The

J. Phys. Chem. C, Vol. 113, No. 31, 2009 13737

Figure 5. Frequency of the stretching mode of the different hydrogen atoms involved in hydrogen bonds in the simulation as a function of the hydrogen bond distance. The insets on the upper left corner of the figure show the VDOS of three representative cases. In the lower right corner we show the hydrogen bond distance distribution of the same configurations. The figure highlights the correlation between the bond distance and the vibrational frequency of hydrogen bonds. The symbols are defined in Table 1. Empty symbols correspond to surface water hydrogen bonds, while filled symbols denote either surface-surface or water-water hydrogen bonds. The color indicates the temperature with the same code as in Figure 3, black for 280 K, orange for 300 K, and red for 320 K.

different donors or acceptors in our system and their symbols are as follows (Figure 5, Table 1). BO atoms are represented as Ti2O or as Ti2OH if protonated. Oxygen atoms attached to the five-fold coordinated Ti (L1) are represented as TiHOH if they are in their associated form or as TiOH if they are in their dissociated form as hydroxyl groups. Finally, free water is represented as HOH. The two atoms participating in the hydrogen bond, the accepting oxygen atom, and the donated hydrogen atom are represented in boldface. We have scanned our simulations for the appearance of bonds between these species. During the times that an H-bond is formed, we record the hydrogen bond distance and the contribution to the VDOS of the hydrogen atom involved. The results of this scan are shown in Table 1. Previous molecular dynamics studies performed in smaller unit cells show hydrogen bond distance between water molecules in the hydration layer and water molecules in L1 in the range of 1.5-1.9 Å, in good agreement with our values reported in Table 1.52 Some of the hydrogen bond configurations appear in two different situations depending on the number of hydrogen bonds that the oxygen atom is accepting. These are denoted in the first column of the table as I, for the single bond configuration, and II, for the case of two bonds. We draw this distinction because the properties are quite different, and this is a very interesting effect that can be seen in the table. Hydrogen bond distances are larger, and frequencies are higher in the case of configurations where the oxygen atom is accepting two bonds because when oxygen is accepting two bonds, the individual bonds are weaker. It is important to note that hydrogen bonds show some nonintuitive properties. First, when the hydrogen bond is weaker, the covalent bond between the donor and the hydrogen becomes stronger, and the vibrational frequency goes up. Correspondingly, the stronger hydrogen bonds correspond to the lowest frequencies. This is shown in Figure 5 where the values of Table 1 are displayed. The plot shows a strong correlation between the acceptor-hydrogen distance and the vibrational frequency associated to that hydrogen bond. The insets in the upper-left corner of the figure show the vibrational

13738

J. Phys. Chem. C, Vol. 113, No. 31, 2009

Kumar et al.

TABLE 1: Distances and Characteristic Frequencies of the Different Hydrogen Bonds Present in the (110) Surface of Rutile Exposed to Watera 280 K

300 K

case

P

d

ν

Ti2O · · · TiHOH (2) I Ti2O · · · TiHOH (2) II Ti2O · · · HOH (∆) TiOH · · · TiHOH (b) TiOH · · · HOH (O) TiHOH · · · TiHOH (1) Ti2OH · · · HOH (3) TiHOH · · · HOH ()) I TiHOH · · · HOH ()) II HOH · · · HOH (() I HOH · · · HOH (() II

19.5 80.5

1.68 1.77 1.63 1.60 1.70 1.71 1.66 1.57 1.66 1.68 1.73

91.6 100.3 93.9 88.8 94.4 98.6 92.5 84.8 92.4 96.0 101.2

12.9 87.1 25.8 74.2

P

27.6 72.4

320 K

d

ν

1.77 1.65 1.61 1.70 1.72 1.65

97.6 94.4 88.9 95.7 97.7 91.3

1.63 1.69 1.76

91.3 94.3 101.7

P

28.1 71.9

d

ν

1.76 1.66 1.62 1.64 1.75 1.67

97.4 95.4 88.5 94.4 100.1 92.3

1.68 1.70 1.77

94.4 98.6 100.8

a The frequencies (ν) are expressed in THz and the distances (d) in Å. These distances are the peak position of the respective histograms. In some bonding configurations the oxygen may be accepting one or two hydrogen bonds; these situations are labeled with roman numerals I and II besides the configuration case, and the column P displays the relative abundance of these configurations. The symbol between parentheses corresponds to the code used in Figure 5 where these values show a clear linear dependence between the frequency and the hydrogen bond distance. The empty cells correspond to configurations that are not statistically relevant in our simulations; they correspond to configurations whose probabilities are so low that we do not have enough statistical data.

TABLE 2: Characteristics of the Strongest Hydrogen Bonds in the Water-Titanium Interface H-bond

configurationa

d relaxed (Å)a

d range (Å)a

d peak (Å)b

d range (Å)b

frequency (THz)b

Ti2OH-HOH Ti2O-HOH TiHOH-HOH TiOH-HOH

M2(1) M2(2) M2(2) M2(1)

1.74 1.77 1.76

1.6-2.1 1.4-1.9 1.4-1.8

1.66 1.63 1.66 1.70

1.4-1.9 1.4-2.0 1.4-2.0 1.4-2.0

92.5 93.9 92.4 94.4

a Values from Zhang and Lindan, ref 21. b This work. This table lists the four strongest H-bonds found in the water-titanium interface, the configurations where this bond type is found in the previous work by Zhang and Lindan, the values found by them for the distance after relaxation, and their ranges taken from a molecular dynamics simulation of 1.2 ps at 130 K. Also listed are the peaks of the histogram of distances found in our work, their ranges, and the averages of the vibrational frequency of the contribution of the H-atom forming the bond to the stretching mode band of the VDOS.

frequency spectra of three different bonds of diverse strength. The VDOS for the weaker hydrogen bond, shown in the inset that corresponds to the higher vibrational frequency average, also shows a narrower distribution of frequencies compared with the stronger bonds. In the lower-right corner of the figure, we show the bond length histograms of the same hydrogen bonds. The width of the distribution of distances follows a trend opposite to that of the frequency distribution. As expected, the stronger bonds have a narrower distribution of bond lengths compared with that of the weaker bonds. The data presented here confirm the picture that the bonding between the (110) surface of rutile and water is dominated by strong hydrogen bonds produced by the BO atoms and the tightly bound water molecules in L1. To give further support to this picture, we now consider the nature of the bonding leading to the strong adsorption sites shown to appear in L2. As we have discussed before, the water molecules in L2 have two well-defined absorption sites. One forms an “umbrella” above the BO site, and the other lies above and between L1 water molecules (see Figure 2). The interactions that produce these adsorption sites are the strong hydrogen bonds at the interface. The water molecules adsorbed above the BO atom are bonded to it by one strong hydrogen bond, while the molecules adsorbed between the L1 water molecules are supported by two bonds in a bridging position between them. There are four configurations that link L2 with the surface, and those are the strongest hydrogen bonds in the system. The main dynamical features of these configurations are listed in Table 2 along with similar values obtained by Zhang and Lindan.21

Table 2 shows that the same strong bonds that bind the extra water molecule added in the simulations of Zhang and Lindan are the strong bonds appearing in a simulation when a larger water coverage is present. This means that these strong bonds dominate the water surface interaction. As listed in the last column, the average frequency of the stretching mode band is an indication of a strong hydrogen bond. This result validates the idea that was first started by Lindan and co-workers,20 that strong hydrogen bands are responsible for the large red shift observed in the VDOS of these surfaces. Two of these configurations are shown as an inset of Figure 6. In this configuration we can distinguish three types of bonds. Bonds between the BO atom and a water molecule (Ti2O · · · HOH), a typical water-water interaction (HOH · · · HOH), and a bond between the water and the water molecules adsorbed in the L1 layer (TiHOH · · · HOH). From the vibrational contribution of each bond we see that Ti2O · · · HOH and TiHOH · · · HOH are stronger bonds than HOH · · · HOH, that is, that the bonds with the surface are stronger. It is also very interesting to note that although both stretching mode peaks for Ti2O · · · HOH and TiHOH · · · HOH are centered at about the same frequency, the band corresponding to TiHOH · · · HOH is wider. This width reflects the fact that the water molecule in L2 above TiHOH · · · HOH has significant oscillation leading to the weakening and strengthening of the hydrogen bond, which in turn produces a wider variation in the stretching band. The other two configurations, namely, TiOH · · · HOH and Ti2OH · · · HOH, show a similar frequency distribution indicating strong

H-Bonds and Vibrations of Water on (110) Rutile

Figure 6. VDOS contribution of the hydrogen atoms involved in three different types of hydrogen bonds at the interface. The different types of bonds are shown in the upper panel. As a reference, we include the contribution from a regular water-water hydrogen bond. The strength of the TiHOH · · · HOH and Ti2O · · · HOH bonds is apparent from the lower vibrational frequencies of the stretching bands associated with them.

hydrogen bonding. These are the contributions that dominate the red shift of the stretching mode band of the VDOS. 4. Conclusions In this study we have used extended time aiMD simulations, and we have calculated the hydrogen VDOS to study the interactions at the interface of (110) TiO2 with water. The stretching band of the vibrational spectrum of water molecules and OH groups is particularly sensitive to hydrogen bonding. In this heterogeneous interface, there are a large variety of hydrogen bonds between the surface species and the water and among the water molecules at the surface. Each of these types has different characteristic hydrogen bond strengths and, consequently, a characteristic contribution to the stretching mode band of the associated water molecules or hydroxyls. In particular, INS measurements of water confined at the surface of the anatase polymorph of TiO2 showed a surprising softening of this band, extending to frequencies as low as 85 THz (350 meV),29 which we also observe in the VDOS extracted from our simulations. We provide evidence that the softening is produced by strong hydrogen bonding between the surface species and the hydration layer. We identified four very strong hydrogen bonds between the atoms at the surface and the water molecules in the hydration layer. These bonds dominate the contributions to the VDOS in the lower end of the OH stretching band. One of them is the hydrogen bond between BO atoms at the surface and a water molecule in the L2 hydration layer. The other bond forms between a water molecule or a hydroxyl in L1 and a water molecule in L2. These bonds are characterized by a low average frequency of the stretching band due to the strength of the hydrogen bonds they are forming. Additional evidence of this statement is provided by the formation of an epitaxial structure of the water layer in contact with the surface (L2) shown in Figure 2. The two adsorption spots shown in the figure correspond to the strong hydrogen bonds identified. The bond with the BO atom has a better defined vibrational frequency of the stretching mode band. This is because the water

J. Phys. Chem. C, Vol. 113, No. 31, 2009 13739 molecule in its movement keeps the distance to the BO mostly constant. This is clearly seen in the x-z panel of Figure 2 where the density on top of the BO forms an umbrella-shaped cloud. In the other bond, the water molecule in L2 is accepting two bonds from two water molecules in L1. As a consequence, it oscillates between these two configurations, weakening one bond and at the same time strengthening the other. This oscillatory movement broadens the width of the stretching mode band because it generates vibrational frequencies that are higher and lower than the frequency of a single hydrogen bond. It is also interesting to note that these two bonds that we have identified are very strong for different reasons. Interestingly enough, a hydrogen bond can be strong because the acceptor is a strong acceptor or because the donor forms a weak bond with the hydrogen atom involved in the bond. In these two bonds at the surface of rutile, we have almost prototypical examples of these two cases. The BO atom is undercoordinated with respect to other oxygen atoms in the bulk crystal. It has a higher negative charge and attracts the proton strongly, making this a strong hydrogen bond. Opposite to this case, the hydrogen bond between L2 and L1 water molecules is strong because the oxygen atom in L1 water molecules is over-coordinated. It is an oxygen atom in a water molecule that, in addition to its two hydrogen atoms, forms a relatively strong bond with the fivefold coordinated titanium atom at the surface. The weaker covalent bonds in these L1 water molecules form a stronger hydrogen bond with the L2 water molecules of the hydration layer. Because of the weaker covalent bonds formed between oxygen atoms and hydrogen atoms in the L1 layer, these water molecules are prone to dissociate and form a hydroxylated surface. The VDOS can be measured by different experimental techniques, including INS of hydrated nanoparticles. The correlations demonstrated in our work open the possibility to use the vibrational spectrum to study and characterize the hydrogen-bonding interactions at this type of interfaces. Acknowledgment. This work was supported by a grant from the U.S. Department of Energy, Office of Basic Energy Sciences, Geoscience Research Program to Oak Ridge National Laboratory, which is operated by UT Battelle, LLC under Contract No. DE-AC05-00OR22725 for the U.S. Department of Energy. This research used resources of the National Center for Computational Sciences and the Center for Nanophase Materials Sciences at ORNL, which are sponsored by the respective facilities divisions of the DOE Offices of Advanced Scientific Computing Research and Basic Energy Sciences. This work was supported in part by the Materials Simulation Center, a Penn State Center for Nanoscale Science (MRSEC-NSF) and Penn State Materials Research Institute facility. References and Notes (1) Langmuir, I. J. Am. Chem. Soc. 1918, 40, 1361. (2) Ewing, G. E.; Foster, M.; Cantrell, W.; Sadtchenko, V. Thin Film Water on Insulator Surfaces. In Water in Confining Geometries; Bush, V., Devlin, J. P., Eds.; Springer: Berlin, 2003; p 179. (3) Pauling, L. The Nature of the Chemical Bond; Cornell University Press: Ithaca, NY, 1960. (4) Diebold, U. Surf. Sci. Rep. 2003, 48, 53. (5) Henderson, M. A. Surf. Sci. Rep. 2002, 46, 5. (6) Wang, R.; Hashimoto, K.; Fujishima, A.; Chikuni, M.; Kojima, E.; Kitamura, A.; Shimohigoshi, M.; Watanabe, T. Nature 1997, 388, 431. (7) Mills, A.; Crow, M. J. Phys. Chem. C 2007, 111, 6009. (8) Zubkov, T.; Stahl, D.; Thompson, T. L.; Panayotov, D.; Diwald, O.; Yates, J. T. J. Phys. Chem. B 2005, 109, 15454. (9) Petrik, N. G.; Kimmel, G. A. Phys. ReV. Lett. 2007, 99, 196103. (10) Bates, S. P.; Kresse, G.; Gillan, M. J. Surf. Sci. 1998, 409, 336. (11) Goniakowski, J.; Gillan, M. J. Surf. Sci. 1996, 350, 145.

13740

J. Phys. Chem. C, Vol. 113, No. 31, 2009

(12) Harris, L. A.; Quong, A. A. Phys. ReV. Lett. 2004, 93, 086105. (13) Lindan, P. J. D.; Muscat, J.; Bates, S.; Harrison, N. M.; Gillan, M. Faraday Discuss. 1997, 135. (14) Lindan, P. J. D.; Zhang, C. J. Phys. ReV. B 2005, 72, 075439. (15) Bandura, A. V.; Kubicki, J. D. J. Phys. Chem. B 2003, 107, 11072. (16) Bandura, A. V.; Kubicki, J. D.; Sofo, J. O. J. Phys. Chem. B 2008, 112, 11616. (17) Prˇedota, M.; Bandura, A. V.; Cummings, P. T.; Kubicki, J. D.; Wesolowski, D. J.; Chialvo, A. A.; Machesky, M. L. J. Phys. Chem. B 2004, 108, 12049. (18) Casarin, M.; Maccato, C.; Vittadini, A. J. Phys. Chem. B 1998, 102, 10745. (19) Lindan, P. J. D.; Harrison, N. M.; Gillan, M. J. Phys. ReV. Lett. 1998, 80, 762. (20) Lindan, P. J. D.; Harrison, N. M.; Holender, J. M.; Gillan, M. J. Chem. Phys. Lett. 1996, 261, 246. (21) Zhang, C. J.; Lindan, P. J. D. J. Chem. Phys. 2003, 118, 4620. (22) Machesky, M. L.; Prˇedota, M.; Wesolowski, D. J.; Vlcek, L.; Cummings, P. T.; Rosenqvist, J. R.; Ridley, M. K.; Kubicki, J. D.; Bandura, A. V.; Kumar, N.; Sofo, J. O. Langmuir 2008, 24, 12331. (23) Mamontov, E.; Vlcek, L.; Wesolowski, D. J.; Cummings, P. T.; Wang, W.; Anovitz, L. M.; Rosenqvist, J.; Brown, C. M.; Sakai, V. G. J. Phys. Chem. C 2007, 111, 4328. (24) Wendt, S.; Matthiesen, J.; Schaub, R.; Vestergaard, E. K.; Laegsgaard, E.; Besenbacher, F.; Hammer, B. Phys. ReV. Lett. 2006, 96. (25) Allegretti, F.; O’Brien, S.; Polcik, M.; Sayago, D. I.; Woodruff, D. P. Phys. ReV. Lett. 2005, 95, 226104. (26) Harris, L. A.; Quong, A. A. Phys. ReV. Lett. 2005, 95, 029602. (27) Lindan, P. J. D.; Zhang, C. J. Phys. ReV. Lett. 2005, 95, 029601. (28) Beck, T. J.; Klust, A.; Batzill, M.; Diebold, U.; Di Valentin, C.; Tilocca, A.; Selloni, A. Surf. Sci. 2005, 591, L267. (29) Levchenko, A. A.; Kolesnikov, A. I.; Ross, N. L.; Boerio-Goates, J.; Woodfield, B. F.; Li, G.; Navrotsky, A. J. Phys. Chem. A 2007, 111, 12584.

Kumar et al. (30) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (31) Svishchev, I. M.; Kusalik, P. G. J. Phys. Chem. 1994, 98, 728. (32) Gale, J. D.; Rohl, A. L. Mol. Simul. 2003, 29, 291. (33) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133. (34) Hohenberg, P.; Kohn, W. Phys. ReV. 1964, 136, B864. (35) Kresse, G.; Hafner, J. Phys. ReV. B 1994, 49, 14251. (36) Kresse, G.; Hafner, J. Phys. ReV. B 1993, 47, 558. (37) Kresse, G.; Furthmu¨ller, J. Phys. ReV. B 1996, 54, 11169. (38) Kresse, G.; Furthmu¨ller, J. Comput. Mater. Sci. 1996, 6, 15. (39) Kresse, G.; Joubert, D. Phys. ReV. B 1999, 59, 1758. (40) Blo¨chl, P. E. Phys. ReV. B 1994, 50, 17953. (41) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1997, 78, 1396. (42) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (43) Kneller, G. R.; Hinsen, K. J. Chem. Phys. 2001, 115, 11097. (44) Rog, T.; Murzyn, K.; Hinsen, K.; Kneller, G. R. J. Comput. Chem. 2003, 24, 657. (45) Luzar, A.; Chandler, D. Phys. ReV. Lett. 1996, 76, 928. (46) Luzar, A.; Chandler, D. Nature 1996, 379, 55. (47) Hammerich, A. D.; Buch, V. J. Chem. Phys. 2008, 128, 111101. (48) Matsumoto, M. J. Chem. Phys. 2007, 126, 6. (49) Kumar, R.; Schmidt, J. R.; Skinner, J. L. J. Chem. Phys. 2007, 126, 12. (50) Zhang, Z.; Fenter, P.; Sturchio, N. C.; Bedzyk, M. J.; Machesky, M. L.; Wesolowski, D. J. Surf. Sci. 2007, 601, 1129. (51) Levchenko, A. University of California at Davis, Davis, California. Personal communication, 2008. (52) Zhang, C.; Lindan, P. J. D. J. Chem. Phys. 2003, 119, 9183.

JP901665E