ARTICLE pubs.acs.org/JPCC
Adsorption of Zn2þ on the (110) Surface of TiO2 (Rutile): A Density Functional Molecular Dynamics Study A. V. Bandura,† J. O. Sofo,‡ and J. D. Kubicki§,* †
Department of Quantum Chemistry, St. Petersburg State University, St. Petersburg, Russia Department of Physics and §Department of Geosciences and the Earth & Environmental Systems Institute, The Pennsylvania State University, University Park, Pennsylvania 16802, United States
‡
ABSTRACT: Adsorption of Zn2þ at the rutile TiO2 (110)-aqueous interface was studied with BornOppenheimer molecular dynamics at 300 K. Simulations were carried out using the periodically repeated slab model with vacuum gap of 15 Å filled with 72 H2O molecules. Two possible adsorption sites, monodentate above bridging oxygen (TiOTi or Obr) and bidentate above terminal oxygens (TiO), were investigated. Sites with different coordination environment for adsorbed Zn2þ differ from each other by the position of Zn2þ above surface level and by characteristic ZnO distances. Obtained results gave evidence that 4-fold coordination of adsorbed Zn2þ is more probable than the 6-fold coordination found for aqueous species. The hydrolysis of H2O molecules was observed in the first coordination shell of adsorbed ion, resulting in formation of OH groups attached to Zn2þ. Calculated energies favor the tetrahedral bidentate structure of hydrated Zn2þ on the rutile surface. The model structures are compared to observed positions of Zn2þ above the rutile (110) surface using X-ray scattering techniques.
1. INTRODUCTION Adsorption of cations to oxide surfaces is an important process in catalysis and environmental transport of contaminants. For decades, conceptualization of this phenomenon has been guided by various constructions of electrical double-layer theory.14 Models of the oxidewater interface have been constructed based on observations of macroscopic experiments such as potentiometry and electrophoresis.5,6 An atomic-level picture of the oxidewater interface and adsorbed species has been elusive due to the difficulty in obtaining structural information on the interface without interference from the bulk oxide and aqueous phases. Spectroscopic methods, such as ATR-FTIR spectroscopy, have been applied to understand adsorption of molecular species,7,8 but such techniques are not extremely useful when applied to adsorbed metal cations, such as Zn2þ because such species generally exhibit broad, low-frequency bands that are difficult to assign to specific bonding environments. Recently, the use of X-ray standing wave techniques (XSW)913 has shown promise in determining the atomic-level structure of adsorbed ions. However, even when combined with extended X-ray adsorption fine structure (EXAFS) spectroscopy,14 the analytical data can be difficult to interpret and result in different interpretations.13,15 Consequently, quantum mechanical calculations have been applied to the same adsorption reactions to help interpret subtle discrepancies between spectroscopic interpretations.13 Although energy minimizations performed using 3D periodic density functional theory (DFT) can produce configurations of adsorbates consistent with experimental observations, they are limited in the configurational space that can be efficiently r 2011 American Chemical Society
sampled and by the fact that they represent 0 K rather than the temperature of the experiment. These limitations have led to the present study where the DFT method is applied to molecular dynamics (MD) simulations to model the system at a finite temperature over a period of time. Because DFT calculations on the number of atoms included in this study are computationally demanding, the level of theory employed here is not as rigorous as for the energy minimization calculations. Even with this short cut, simulation times are on the order of 10 ps when full equilibration could require an order of magnitude longer duration. Our strategy is the following: we generate the initial structures for the DFT-based MD simulations using more efficient, classical molecular dynamics. This step brings the configuration into thermodynamical equilibrium. The resulting structures are used as initial structures to the DFT-MD that, despite its relaxed accuracy, brings the system into local chemical equilibrium because it introduces chemical reactivity into the simulation. DFT-MD is then used to generate input atomic structures to the more accurate atomic relaxations performed with a more accurate level of the theory. Ultimately, the comparison of these results against observed XSW and EXAFS data are used to determine whether the predictions are accurate.
2. COMPUTATIONAL DETAILS All calculations were performed using the plane-wave pseudopotential implementation of density-functional theory (DFT) Received: January 15, 2011 Revised: February 28, 2011 Published: April 22, 2011 9608
dx.doi.org/10.1021/jp200432p | J. Phys. Chem. C 2011, 115, 9608–9614
The Journal of Physical Chemistry C as written into the VASP code.16,17 Exchange and correlation were treated within the PerdewBurkeErnzerhof (PBE) functional.18 The accuracy of this functional was tested for its ability to reproduce observed binding energies and geometries for the water dimer and small water clusters. Ionic cores were described by the projector-augmented wave (PAW) method,19,20 which improves transferability and reduces the number of plane waves required in the expansion of the KohnSham orbitals in order to improve computational efficiency. A Ti-atom pseudopotential, in which the semicore states (3s and 3p) were treated as valence state, was chosen for this study. Similarly, the valence shell of Zn includes the 3d states. Because of the relatively large system size, an energy cutoff (Ecut), 212 eV, has been applied to make MD simulations of thousands of time steps practical. The applicability of this cutoff was tested by performing several static calculations on previously studied systems. Also, a DFT-MD test run with Ecut = 400 eV was made for one of the considered systems, which confirmed the validity of using 212 eV cutoff (below). All MD simulations were performed with single Γ-point sampling of the Brillouin zone. Convergence criterion in the electronic cycle (i.e., computation of the electron energy and density at fixed geometry) was set to 106 eV to provide sufficient accuracy for the energy conservation. The Nose thermostat21 was used to simulate the canonical ensemble at constant temperature and volume. Most of the calculations were conducted at 300 K. The Nose-mass parameter of thermostat corresponds to the period of 20 fs of induced temperature fluctuations with a time step of 0.5 fs. Thermal equilibration was achieved in 2 ps in most cases. Total simulation times were approximately 10 ps. A 0.25 ps period was applied for obtaining the intermediate (partial) averages.
3. PREPARATION OF INITIAL STRUCTURES Because of the mentioned computational restrictions, a TiO2 slab with the minimal reasonable thickness of three Ti-layers parallel to (110) of the R-TiO2 (rutile) crystal was accepted for simulations. Although such slabs could result in an inaccurate ratio for associated and dissociated adsorbed molecules, the structures and relative energies between configurations are not expected to be a strong function of slab thickness if the number of dissociated and associated forms of water molecules attached to the TiO2 surface is not changed. One can also suppose that interaction of water with Zn2þ ion is not influenced considerably by the slab thickness. The 4 2 surface supercell was constructed in the (110) plane. The dimensions for the cell were a = 11.836 Å and b = 12.994 Å, in x and y directions correspondingly, and c = 24 Å in the z direction. A gap between two periodic slabs of approximately 15 Å was filled with 72 H2O molecules and 2 Zn2þ ions, one on each side of the slab. Additionally, a pair of terminal hydroxyl groups was added on each surface of the slab to ensure the charge-neutrality. The resulting total number of atoms in the simulation cell was 370. Positions of middle-layer atoms in the TiO2 slabs were fixed at bulk rutile geometry in all systems. As has been shown in ref 22, the internal constraints, consisting of freezing some internal layers to their atomic bulk positions, help to stabilize the surface energy for relatively thin systems that allow working on larger surfaces. A view of the simulation cell is shown in Figure 1. The number of H2O molecules was estimated from the condition that the average liquid density should be approximately 1.0 g cm3. Obtained mean water density in the
ARTICLE
Figure 1. Simulation cell general view. Initial geometry for the B6I system is shown.
central layer with thickness 4 Å is in the range or 1.051.07 g cm3 for all simulated systems. Initial positions of the H2O molecules were prepared via classical MD simulations using a previously derived force field23 as described below. One goal of this investigation was to estimate the most probable position and structure of the hydrated Zn2þ on the R-TiO2 (110) surface. Experimental data12,13 and our previous static calculations13 showed that there are two most probable positions: monodentate above the bridging oxygen and bidentate above the two neighboring terminal OH groups. Moreover, DFT static calculations gave evidence that the total coordination number of Zn2þ is lower than 6, and one or more coordinated H2O molecules dissociates to produce OH groups attached to Zn. These previous results were taken into account when the initial structures were created for DFMD simulations. Accordingly, three different types of systems have been considered: (B6) the monodentate 6-fold coordinated Zn2þ above the bridging oxygen; (B4) the monodentate 4-fold coordinated Zn2þ above the bridging oxygen; and (T4) the bidentate Zn2þ above the two neighboring terminal OH groups. Because of the large size of the model systems, the proper choice of initial coordinates is critical for providing the reasonably short equilibration times. This means that the water subsystem must be near the equilibrium state before the simulation starts. To satisfy this requirement, two different strategies of preparing the initial particle positions were applied. The first utilized the previously derived force field.23 A classical MD simulation of approximately 100 ps was performed to preequilibrate the B6 and B4 configurations. In the latter case, the Zn2þ complexes were fixed at the previously determined static DFT structure13 because force field parameters for ZnOH 9609
dx.doi.org/10.1021/jp200432p |J. Phys. Chem. C 2011, 115, 9608–9614
The Journal of Physical Chemistry C
ARTICLE
Figure 3. Structures obtained in B6 systems at the end of DFMD runs. (a) B6I, top Zn; (b) B6I, bottom Zn; (c) B6II, top Zn; (d) B6II, bottom Zn.
Figure 2. Running distances (Å) between Zn2þ and the reference (110) Ti-atomic plane in initially monodentate Zn2þ surface complex: (a) B6I starting geometry, average values for last 3.25 ps: 3.91 Å (top Zn), 3.26 Å (bottom Zn); (b) B6II starting geometry, average values for last 4.00 ps: 3.60 Å (top Zn), 3.43 Å (bottom Zn); (c) B4 starting geometry, average values for 5.5 ps: 3.24 Å (top Zn) and 3.30 Å (bottom Zn).
interactions were not available. The classical force field was developed for aqueous Zn2þ in 6-fold coordination; and, hence, it cannot reproduce the tetrahedral environment found in static DFT calculations. Several low-energy states were chosen from the classical trajectory and optimized by applying the force field potentials. Finally, the state with the lowest energy was optimized by DFT (to maximum force on atom of 0.05 eV/Å) prior to the start of the DFT-MD simulation. The second approach is based on the idea that increasing the simulation temperature can considerably speed up the translational diffusion in model systems, producing the H2O distribution that is sufficiently different from the initial state and is not far from equilibrium, if the starting structure was previously equilibrated at low temperature (i.e., simulated annealing). Again, a low-energy
state may be chosen from the end of the DFT-MD trajectory as the starting structure for further simulation at 300 K. Nevertheless, optimization of this state may be necessary to generate a final configuration representative of the lower temperature state. If the optimized energy of the chosen state noticeably differs from energies of other structures chosen, the results may not be comparable due to relatively short simulation times (i.e., ∼10 ps). Simulations of several picoseconds cannot reproduce the translational diffusion of H2O because the corresponding characteristic times of the rotational and translational components of diffusion motion are larger at least by an order of magnitude.2426 Hence, simulated DFT-MD trajectories sample a relatively small configuration space and mainly rotational diffusion only is reflected. The second scheme was applied for constructing of T4 structure from the previously equilibrated B6 structure. Details of these procedures will be discussed further below in context of the structural results.
4. DESCRIPTION OF MODELS: INITIAL AND CALCULATED STRUCTURES 4.1. Octahedral Monodentate Site above Bridging Oxygen. As stated above, the initial positions of H2O molecules and
Zn2þ ions in the first octahedral monodentate structure B6I were prepared using classical MD simulation. The starting structure is shown in Figure 1. The geometry of adsorption complexes on both sides of the TiO2 slab was not equal because symmetry restrictions were not introduced for the MD simulations. The B6I system was equilibrated within 2 ps. After t = 1 ps, one of the adsorbed Zn2þ approached the surface and became 4-fold coordinated. The distance between this Zn2þ and Obr lowered to approximately 3.25 Å. The second Zn2þ moved farther away from the surface positioning approximately 4 Å from the 9610
dx.doi.org/10.1021/jp200432p |J. Phys. Chem. C 2011, 115, 9608–9614
The Journal of Physical Chemistry C
Figure 4. 4-fold coordinated complexes of Zn2þ on (110) TiO2 surface obtained at the end of DFT-MD runs. B4 system: (a) top Zn, (b) bottom Zn; T4 system: (c) top Zn, (d) bottom Zn.
corresponding Obr and remained in octahedral coordination. This general configuration did not change during the rest of the simulation. In part a of Figure 2, the time-dependence of the distances between Zn2þ ions and the nearest Ti atomic (110) plane in unrelaxed (bulk geometry) slab is shown. There are clearly two separate Zn2þ-surface distances. Final configurations of these inner- and outer-sphere adsorbed Zn2þ ions are illustrated in parts a and b of Figure 3. System B6II was also constructed using the same procedure as the previous one, but B6I and B6II initial structures are separated by 200 ps of classical MD trajectory. In this case, after 1 ps pre-equilibration period, the positions of both Zn2þ ions were slightly translated manually so they were in equivalent positions relative to the cell boundaries at a distance of approximately 2.1 Å from the corresponding bridging oxygens (this distance is close to the optimal value). During the next 2 ps of DFT-MD simulation, the Zn2þ positions were kept fixed. The constraint was then released and both ions were allowed to move freely at 300 K. The time-dependence of the distances between Zn2þ ions and the reference (110) plane is given in part b of Figure 2. Although both Zn2þ ions remain bonded to the bridging O atom during all simulation time, one of them (at the bottom position) approaches the surface and reduces its average coordination number to 5. The final structures of the Zn2þ complexes are shown in parts c and d of Figure 3. Thus, for both regarded B6 systems considered, there is a clear tendency for reducing the coordination number of the adsorbed Zn2þ during the progress of DFT-MD trajectory. The last structure of the B6II run and corresponding particle velocities were taken as the initial configuration for next B6III test run with Ecut = 400 eV. The average Zn2þ coordination number was approximately 5.4 during all 1.4 ps of this simulation, and the average distances between Zn2þ ions and the reference
ARTICLE
(110) plane remain practically the same as in the previous run with B6II starting geometry. Thus, increasing of energy cutoff up to 400 eV does not noticeably influence the results obtained with B6II starting geometry using Ecut = 212 eV. At the same time, the simulation with the fixed Zn positions was continued at 800 K for another check of the stability of octahedral complex. This calculation led to a 4-fold coordinated structure (below). 4.2. Tetrahedral Monodentate Site above Bridging Oxygen. The tetrahedral monodentate structure (B4) was prepared via classical force field MD simulation as described above. No change of the Zn2þ total coordination number was found in the B4 system during the 7.5 ps DFT-MD simulation. The final geometry of the 4-fold coordinated complex in B4 system is presented in parts a and b of Figure 4. One of two Zn2þ in the B6I structure that converted to B4 type spontaneously during the equilibration period also represents a similar state. Comparison of the time-dependence of vertical Zn2þ displacements shown in parts a and c of Figure 2 confirm this conclusion. The formation of tetrahedral complexes from octahedral also was observed when the B6II system was heated for short time (0.5 ps) to 800 K and then gradually over 0.25 ps cooled to 300 K keeping the positions of Zn ions fixed. Recent experimental studies15,27,28 have shown that Zn2þ often changes its coordination environment from a 6-fold coordination in aqueous solutions to a 4-fold coordination on adsorption. The conversion of 6-fold structures to 4-fold structures provides evidence that this might also happen at the rutile (110) aqueous solution interface. The reason for the observed phenomenon may be water hydrolysis. The number of hydroxyls is thought to be a controlling factor driving the 6- to 4-fold coordination change because hydrolysis drives a similar Zn2þ coordination change in aqueous solution. The deprotonation of a water molecule from Zn(H2O)62þ resulted in a shorter ZnOH bond and release of at least one H2O. All of the hydrolysis species had a coordination number of less than 6.29 However, the number of hydroxyl groups directly attached to Zn2þ is different in the various models. The complex originated from the B6 system has H2O molecules as the only ligands, whereas the number of OH groups oscillates between 0 and 1 for one Zn2þ in B4 system, and it is exactly 2 for the other. Another explanation is that there is not enough room for 6 ligands near the Zn2þ when it is attached to a bridging oxygen atom on the surface because the neighboring OH groups or water molecules at the terminal positions (above Ti(V)) are too close to the Zn2þ-bonded ligands. As a result, stabilization occurs by reducing the coordination number and shortening the ZnObr distance. Accordingly, in the tetrahedral structure (as can be seen from part c of Figure 2) the Zn2þ is located on average more closely to the reference Ti-atomic surface plane than in B6 systems with octahedral coordination. Thus, the averaged distance from Zn2þ to the reference plane is about 3.50 Å for octahedral complexes whereas it is approximately 3.25 Å for the tetrahedral complexes. 4.3. Tetrahedral Bidentate Site above Terminal Oxygens. A bidentate structure above the two hydroxide groups was prepared from the state obtained in high-temperature DFTMD simulation. A low-energy state was selected from the 800 K trajectory, then simulated for 0.5 ps at the same temperature after manual correction of the position of one of the Zn2þ ions to provide approximate equivalence of the two sites on the opposite slab surfaces. The positions of both Zn2þ ions were 9611
dx.doi.org/10.1021/jp200432p |J. Phys. Chem. C 2011, 115, 9608–9614
The Journal of Physical Chemistry C
ARTICLE
Table 1. Relative Energies (kJ/mol) of Zn2þH2O Complexes on the (110) TiO2 Surface Per One Zn2þ Obtained by Static Optimization of Selected Structures Zn position
system
5 H2O, 2 1 cell
72 H2O, 4 2 cella initial
final
above the bridging oxygen (monodentate site)
B6I
15
79
6-fold coordinated Zn2þ above the bridging oxygen (monodentate site)
B6II
0
0
115
4-fold coordinated Zn2þ above the bridging oxygen (monodentate site)
B4
59
3
144
4-fold coordinated Zn2þ above the two terminal hydroxyls (bidentate site)
T4
115
72
158
2þ
6-fold coordinated Zn
a
Difference between the energy of given structure and energy of the initial (optimized) structure B6II. Final energy was obtained via the optimization of selected low-energy state at the end of simulated trajectory.
fixed during the period of 0.5 ps. Then, the low-energy structure from the final stage of simulation was optimized. As can be seen from the Table 1, this initial structure has the lowest energy among the considered starting geometries. The simulation was continued up to 6.25 ps after releasing of constraints. During the entire simulation time, the coordination number of Zn remains equal to 4 with two H 2O molecules bonded to the cation. No dissociation of H 2O molecules was observed during the entire simulation time. In this structure, the Zn2þ was located approximately at the same distance (on average 3.25 Å) from the reference Ti-atomic surface plane as in B4 system. The structure of surface complexes at the end of simulation period is presented in parts c and d of Figure 4.
capacity at constant volume by the following equation:
5. RESULTS AND DISCUSSION
where Eint internal (potential plus kinetic) energy. The approximate calculation of Cv from available experimental data for water, rutile,30 and Zn2þ(aq) and OH(aq) ions31 gives the value of 7.4 kJ K1 (per mole of simulation cells). We have obtained the value of Cv in a range of 9.7 to 9.9 kJ K1mole1 for the model systems. This value is not far from the experimental estimation taking into account the uncertainty of the simulated quantities (because of the drift of the calculated energies mentioned before), neglect of quantum effects for high-frequency water vibration, and disregarding the surface effects in experimental values. Hence, we can conclude that our sampling of canonical ensemble was sufficiently correct. 5.2. Distribution Functions and Interatomic Distances for Zn2þO Interactions. Our simulations show that the height of the ion at the B4 and T4 sites (3.25 Å), above the reference Ti plane, is significantly lower than that of the ion at the B6 site (3.50 Å). The first value agrees well with the experimental result13 for the dominant A1 monodentate site. However, XSW and EXAFS data analysis derived a Zn2þ ion height of 2.50 Å at the second (minor) A2 site that is significantly lower than we found by our simulations. Therefore, we can conclude that experimentally derived A2/A4 site does not exactly correspond neither to B4 nor to T4 simulated system. Useful structural information can be obtained from the radial distribution functions (RDF) g(r). In spite of the overall low quality of statistics, some pair correlations can be calculated with reasonable precision. In Figure 5, we show the RDF for ZnO interactions calculated for the regarded models. The difference in Zn environment is fairly apparent in those quantities. The most probable Zn water O atom distance in octahedral monodentate complex is rather large, about 2.1 Å, whereas the same distance is smaller (2.0 Å) in tetrahedral monodentate and bidentate complexes. As expected, the covalent ZnOH bond is also shorter 1.9 Å. The values obtained are in excellent
5.1. Energies. The optimized energies of the selected structures are given in Table 1. Energies for all initial structures prepared by classical force-field MD are close to each other. The T4 initial structure based on high-temperature simulations has a remarkably lower energy. Although this may be due to the difference between the classical and DFT potentials for H2OH2O interactions, we consider that the nature of the sites considered is the determining factor. The energy difference between the systems considered (Table 1) was estimated by static optimization of the selected low-energy states at the end of the simulated trajectories. The difference is not as pronounced as was found in static calculations using a 2 1 cell13 because of the damping influence of the large array of H2O molecules. In this case, the difference between the energies of B4 and T4 structures (∼15 kJ/mol) is much smaller than that obtained for 2 1 cell (∼75 kJ/mol). Accordingly, the results for the 4 2 cell agree better with the experimental observations12,13 by XSW and EXAFS that show the partitioning of Zn2þ between those two sites. Nevertheless, static energies favor the tetrahedral bidentate structure (T4) of hydrated Zn2þ on the rutile surface. The average energies obtained in these DFT-MD simulations are not well-defined. In spite of the good control of the kinetic energy and reasonable fluctuations (below), both the total (for system plus thermostat) and the potential energies have negative drift (∼0.1 eV/ps) in all simulations. This drift is possibly due to insufficient energy conservation provided by the attainable accuracy of electronic energy calculations or (this is less probable) because of relatively short equilibration time. To check the quality of sampling, we have estimated the contribution ckin of kinetic energy Ekin to the dimensionless heat
ckin ¼ ðÆE2kin æ ÆEkin æ2 Þ=ðNmov k2 T 2 Þ
ð1Þ
where Nmov is the number of movable atoms (338 for our systems), and k the Boltzmann constant. Values of ckin obtained for simulated systems lay in an interval from 1.42 up to 1.57 (including a test run with 400 eV energy cutoff), which is acceptable compared to the theoretical value of 1.5 for canonical ensembles. Standard deviation of temperature from 300 K (estimated using block averages over 0.25 ps) does not exceed 0.35 K in all simulations. The total heat capacity at constant volume Cv may be calculated via the fluctuation formula: Cv ¼ ðÆE2int æ ÆEint æ2 Þ=kT 2
9612
ð2Þ
dx.doi.org/10.1021/jp200432p |J. Phys. Chem. C 2011, 115, 9608–9614
The Journal of Physical Chemistry C
ARTICLE
complex is sharper and its maximum at 1.95 Å is again close to corresponding maximum in Zn2þO(H2) RDF. Lastly, the Zn2þ Oterm distribution for the bidentate complex is close to Zn2þOH distribution in the tetrahedral monodentate complex.
Figure 5. Radial distribution functions gZnO(r) and cumulative coordination number NZnO(r) for the Zn2þ ligand oxygen atom interactions calculated for different surface complex structure (surface oxygens are excluded from the distributions). Dashed lines represent fitting by the asymmetric double sigmoidal function.
Figure 6. Radial distribution functions gZnO(r) for the Zn2þ nearest surface oxygen atom interactions calculated for different surface complex structure. Dashed lines represent fitted curves.
agreement with the known ZnO bond lengths of 1.96 and 2.10 Å for the four- and 6-fold coordinated aqueous species.15,3234 Also in Figure 5, we compare results obtained with the B6II starting geometry (4.0 ps run using 212 eV energy cutoff) and with the B6III starting geometry (1.4 ps run using 400 eV energy cutoff). In spite of the noticeably low statistics in the second case, one can see a good agreement between the calculated RDF and cumulative coordination number. This gives additional evidence that using the rather small 212 eV energy cutoff does not lead to incorrect structures for our systems. The RDF’s for the interactions of Zn2þ with the adsorption site O are given in Figure 6. These quantities are more characteristic than the Zn2þ ligand RDF for the given adsorption site. The monodentate octahedral site has a broad distribution with maximum at R(ZnObr) = 2.1 Å, which is close to the position of maximum on the Zn2þ O(H2) RDF shown in Figure 5. The analogous RDF for the tetrahedral monodentate
6. CONCLUSIONS Adsorption of Zn2þ at the rutile TiO2 (110)-aqueous solution interface was investigated by DFT molecular dynamics. Two distinct adsorption sites were considered based on previous experimental and modeling work: the first is a monodentate position above the bridging oxygen, and the second is a bidentate position above the two terminal hydroxyls. In the first case, the possibility of octahedral and tetrahedral coordination environments was investigated. In the second case, the tetrahedral coordination was assumed on the basis of previous static calculations. Although adsorption sites for Zn2þ with different coordination environments differ by the position of the ion above the surface level and by the ZnO distances, these properties are quite similar for the two sites with tetrahedral Zn2þ coordination. The values we obtained are in excellent agreement with the experimental ZnO bond lengths for the 4- and 6-fold coordinated aqueous species. Our simulations show that tetrahedral zinc coordination is more stable than octahedral zinc coordination above the bridging oxygen atoms, and this is probably so for the other adsorption sites. We observe the hydrolysis of H2O molecules in the first coordination shell of the adsorbed ion, resulting in formation of OH groups bonded to Zn2þ. The analysis of MD trajectories shows that, on average, one OH group is attached to Zn2þ in the tetrahedral monodentate site, but zero and two are also possible. Calculated energies for selected states at the end of simulated trajectories favor the tetrahedral bidentate structure of hydrated Zn2þ on the rutile surface. Neither adsorption site agrees with the experimentally derived A2/A4 site13 between bridging and terminal groups, and trial runs failed to determine the low-lying energy states in this area. However, high-temperature simulations show that Zn2þ above the bridging oxygens tends to migrate over a region between two Obr and between one Obr and one Oterm. Nevertheless, further theoretical studies are needed to identify the experimentally observed adsorption site A2/A4.13 ’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT This research was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, under the project “Nanoscale Complexity at the Oxide-Water Interface”. Computation was supported by the Center for Environmental Kinetics Analysis (CEKA), an NSF Environmental Molecular Science Institute and by the Materials Simulation Center, a Penn State MRSEC and MRI facility. ’ REFERENCES (1) Hiemstra, T.; Van Riemsdijk, W. H. On the Relationship between Surface Structure and Ion Complexation of Oxide-Solution 9613
dx.doi.org/10.1021/jp200432p |J. Phys. Chem. C 2011, 115, 9608–9614
The Journal of Physical Chemistry C
ARTICLE
Interfaces. In Encyclopedia of Colloid and Interface Science; Hubbard, A. T., Ed. Marcel Dekker, Inc.: New York, 2002, pp 37733799. (2) Bourikas, K.; Hiemstra, T.; Van Riemsdijk, W. H. Langmuir 2001, 17, 749. (3) Sverjensky, D. A. Geochim. Cosmochim. Acta 2001, 65, 3643. (4) Machesky, M. L.; Wesolowski, D. J.; Palmer, D. A.; Ridley, M. K. J. Colloid Interface Sci. 2001, 239, 314. (5) Bourikas, K.; Hiemstra, T.; Van Riemsdijk, W. H. J. Phys. Chem. B 2001, 105, 2393. (6) Fedkin, M. V.; Zhou, X. Y. Y.; Kubicki, J. D.; Bandura, A. V.; Lvov, S. N.; Machesky, M. L.; Wesolowski, D. J. Langmuir 2003, 19, 3797. (7) Potapova, E.; Grahn, M.; Holmgren, A; Hedlund, J. J. Colloid Interface Sci. 2010, 345, 96. (8) Swedlund, P. J.; Miskelly, G. M.; McQuillan, A. J. Langmuir 2010, 26, 3394. (9) Bedzyk, M. J.; Cheng, L. W. Reviews in Mineralogy and Geochemistry. In X-ray Standing Wave Studies of Minerals and Mineral Surfaces: Principles and Applications.Fenter, P., Rivers, M. L., Sturchio, N. C., Sutton, S. R., Eds.; 2002. Vol. 49 Applications of Synchrotron Radiation in Low-Temperature Geochemistry and Environmental Sciences. Mineralogical Soc America: WA, pp 221266. (10) Zhang, Z.; Fenter, P.; Cheng, L.; Sturchio, N. C.; Bedzyk, M. J.; Machesky, M. L.; Wesolowski, D. J. Surf. Sci. 2004, 554, L95. (11) Zhang, Z.; Fenter, P.; Cheng, L.; Sturchio, N. C.; Bedzyk, M. J.; Predota, M.; Bandura, A.; Kubicki, J. D.; Lvov, S. N.; Cummings, P. T.; Chialvo, A. A.; Ridley, M. K.; Benezeth, P.; Anovitz, L.; Palmer, D. A.; Machesky, M. L.; Wesolowski, D. J. Langmuir 2004, 20, 4954. (12) Zhang, Z.; Fenter, P.; Cheng, L.; Sturchio, N. C.; Bedzyk, M .J.; Machesky, M. L.; Anovitz, L. M.; Wesolowski, D. J. J. Colloid Interface Sci. 2006, 295, 50. (13) Zhang, Z.; Fenter, P.; Kelly, S. D.; Catalano, J. G.; Bandura, A. V.; Kubicki, J. D.; Sofo, J. O.; Wesolowski, D. J.; Machesky, M. L.; Sturchio, N. C.; Bedzyk, M. J. Geochim. Cosmochim. Acta 2006, 70, 4039. (14) Stohr, J.; Jaeger, R.; Brennan, S. Surf. Sci. 1982, 117, 503. (15) Roberts, D. R.; Ford, R. G.; Sparks, D. L. J. Colloid Interface Sci. 2003, 263, 364. (16) Kresse, G.; Furthm€uller, J. Phys. Rev. B 1996, 54, 11169. (17) Kresse, G.; Furthm€uller, J. Vasp the Guide; Institut f€ur Materialphysik, Universit€at Wien: Vienna, Austria, 2003. (18) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (19) Kresse, G.; Joubert, J. Phys. Rev. B 1999, 59, 1758. (20) Bl€ochl, P. E. Phys. Rev. B 1994, 50, 17953. (21) Nose, S. Prog. Theor. Phys. Suppl. 1991, 103, 1. (22) Perron, H.; Vandenborre, J.; Domain, C; Drot, R.; Roques, J.; Simoni, E.; Ehrhardt, J.-J.; Catalette, H. Surf. Sci. 2007, 601, 518. (23) Bandura, A. V; Kubicki, J. D. J. Phys. Chem. 2003, 107, 11072. (24) Mamontov, E.; Wesolowski, D. J.; Vlcek, L.; Cummings, P. T.; Rosenqvist, J.; Wang, W.; Cole, D. R. J. Phys. Chem. C 2008, 112, 12334. (25) Mamontov, E.; Vlcek, L.; Wesolowski, D. J.; Cummings, P. T.; Rosenqvist, J.; Wang, W.; Cole, D. R.; Anovitz, L. M.; Gasparovic, G. Phys. Rev. E 2009, 79, 051504. (26) Takamuku, T.; Hirano, T.; Yamaguchi, T.; Wakita, H. J. Phys. Chem. 1992, 96, 9487. (27) Kosmulski, M.; Maczka, E. Langmuir 2004, 20, 2320. (28) Trainor, Th. P.; Brown, G. E., Jr.; Parks, G. A. J. Colloid Interface Sci. 2000, 231, 359. (29) Zhu, M.; Pan, G. J. Phys. Chem. A 2005, 109, 7648. (30) Chase, M. W., Jr. NIST-JANAF Themochemical Tables, Fourth Edition; J. Phys. Chem. Ref. Data, Monograph 9, 1998, 1-1951. (31) Criss, C. M.; Millero, F. J. J. Sol. Chem. 1999, 28, 849. (32) Christensen, A. N. Acta Chem. Scand. 1969, 23, 2016. (33) Bellotto, M.; Rebours, B.; Clause, O.; Lynch, J.; Bazin, D.; Elkaim, E. J. Phys. Chem. 1996, 100, 8527. (34) McDonald, W. S.; Cruickshank, D. W. Z. Kristallogr. 1967, 124, 180. 9614
dx.doi.org/10.1021/jp200432p |J. Phys. Chem. C 2011, 115, 9608–9614
