ARTICLE pubs.acs.org/JPCC
Periodic Density Functional Theory Study of Water Adsorption on the r-Quartz (101) Surface A. V. Bandura,† J. D. Kubicki,*,‡ and J. O. Sofo§ †
St. Petersburg State University, St. Petersburg, Russia Department of Geosciences and the Earth & Environmental Systems Institute and §Department of Physics, The Pennsylvania State University, University Park, Pennsylvania 16802, United States
‡
bS Supporting Information ABSTRACT: Plane wave density functional theory (DFT) calculations have been performed to study the atomic structure, preferred H2O adsorption sites, adsorption energies, and vibrational frequencies for water adsorption on the R-quartz (101) surface. Surface energies and atomic displacements on the vacuum-reconstructed, hydrolyzed, and solvated surfaces have been calculated and compared with available experimental and theoretical data. By considering different initial positions of H2O molecules, the most stable structures of water adsorption at different coverages have been determined. Calculated H2O adsorption energies are in the range -55 to -65 kJ/mol, consistent with experimental data. The lowest and the highest O-H stretching vibrational bands may be attributed to different states of silanol groups on the watercovered surface. The dissociation energy of the silanol group on the surface covered by the adsorption monolayer is estimated to be þ80 kJ/mol. The metastable states for the protonated surface bridging O atoms (Obr), which may lead to hydrolysis of siloxane bonds, have been investigated. The calculated formation energy of a Q2 center from a Q3 center on the (101) surface with 2/3 dense monolayer coverage is equal to þ70 kJ/mol which is in the range of experimental activation energies for quartz dissolution.
1. INTRODUCTION There have been numerous experimental and theoretical1,2 studies of water interactions with silica surfaces because of the common occurrence of each of these materials in both man-made and environmental systems. Molecular orbital and density functional theory (DFT) calculations at various levels of approximation have been applied to calculate the molecular interaction between different silica species and water. The surface relaxation and reconstruction of silica minerals,3-9 the interaction of their surface functional groups with water and ions,10-14 as well as the mechanism of silica dissolution15-20 have been extensively explored by quantum mechanical methods. Rapid interaction of Si dangling bonds with H2O molecules and formation of silanol surface groups, in which a hydroxylated O atom is bound to only one Si atom, have been predicted previously.5,10 Efforts have been made to develop classical force fields that can be applied to SiO2-water systems. Classical force fields have been used extensively to model quartz-water interfaces using molecular dynamics simulations.21-34 However, the development of reactive force fields that allow bond-making and bondbreaking events such as Hþ transfer has been fewer in number.30-32 Simulations such as these can provide estimates of potential energy surfaces for bond-making and bond-breaking events (e.g., ref 35) that can be used to parametrize reactive force fields. The chemical reactions between reactive water and dangling bonds on a freshly cut silica surface were analyzed32 by using the SERIALREAX force field. These simulations predicted that water molecules penetrate a silica r 2011 American Chemical Society
film through a proton-transfer process called “hydrogen hopping”, which is similar to the Grotthuss mechanism. Computer simulations of water interaction with R-quartz (R-SiO2) surfaces using quantum mechanical methods developed for periodic systems have begun in the past few years. To date, there have been reports investigating a number of reconstructed5,7,9 quartz surfaces. Using first-principles molecular dynamics, Rignanese et al.10 investigated the interactions of water with the (001) R-quartz surface. Both the cleaved (unreconstructed) surface, which presents nonbridging oxygen atoms, and a surface, which is characterized by three-membered and six-membered rings with siloxane bonds, were considered. The cleaved surface was found to be strongly hydrophilic. When the surface is partially hydrated, the formation of H-bond chains was predicted between H2O molecules and adjacent silanol groups. A low-energy electron diffraction (LEED) study of the (001),36 (100), and (101)37 surface structure modifications induced by thermal treatments showed that reconstruction under air occurs at temperatures greater than 500-600 C. The changes in bonding pattern without increasing the unit cell (UC) size on pristine R-quartz low-index surfaces was also found in periodic DFT calculations.5 A DFT molecular dynamics investigation9 of the (001) surface at 300 K showed a densification of the two uppermost Received: November 8, 2010 Revised: February 11, 2011 Published: March 07, 2011 5756
dx.doi.org/10.1021/jp1106636 | J. Phys. Chem. C 2011, 115, 5756–5766
The Journal of Physical Chemistry C
ARTICLE
Table 1. Calculated and Experimental Lattice Parameters of r-Quartz parameter
a
experimenta
direct optimization, 800 eV cutoff
equation of state, 400 eV cutoff
a, Å
4.914
5.03
5.04
c, Å
5.405
5.52
5.52
c/a
1.100
1.097
1.096
Si-O bond length, Å
1.608-1.610
1.625-1.628
1.625-1.628
Si-O-Si angle, deg
143.7
147.4
147.4
K, GPa
37.247
-
31.0
K0
6.047
-
4.4
Ecoh, eV Egap, eV
19.148 9.349
19.0 5.67
19.1 5.67
Geometrical parameters from ref 46.
layers of SiO2 tetrahedral units, with three-membered and sixmembered rings in 1 1 geometry. An experimental X-ray reflectivity and atomic force microscopy study38 of the quartz (100)- and (101)-water interfaces showed that most of the surface silanol groups are single sites (i.e., Q3 = 3 Si-O-Si linkages and 1 Si-OH) with a small percentage of geminal (i.e., Q2) silanediol groups found on the (100) surface. Periodic DFT calculations of the hydroxylated (101) and (001) surfaces4,5 and atomistic calculations of the hydroxylated (001), (100), (101), (101),22 and (011)23 surfaces showed the presence of intersite silanol H-bonds. In this work, periodic DFT calculations of pristine vacuumreconstructed (R-SiO2 only), hydroxylated (SiOH surfaces), and solvated (SiOH þ H2O) (101) surfaces of R-quartz were performed. This surface was selected because it has been previously studied via sum-frequency generation (SFG) spectroscopy39 and X-ray reflectivity,38 so it is possible to make detailed comparisons between DFT results and experimental data. The results are also intended to serve as a benchmark for SiO2-H2O reactive force fields. H-bond networks between the silanol Si(OH) groups have been characterized for the hydroxylated surface. Structure and energetics of adsorbed water layers at different coverage have been analyzed to facilitate the understanding of mechanisms of water adsorption on silica surfaces.
2. COMPUTATIONAL DETAILS Periodic DFT calculations were performed with the VASP simulation package.40,41 The gradient-corrected exchange correlation functional of PBE42 was used. Electronic wave functions were modeled with the plane wave (PW) basis set using projectoraugmented wave (PAW) pseudopotentials.43,44 The energy cutoff for the plane wave basis set was fixed at 400 eV in most surface simulations. The bulk UC of R-quartz has hexagonal symmetry, corresponding to the P3221 space group. Prior to the surface calculations, the bulk UC parameters were optimized. Two different cell optimization methodologies were used: (1) direct optimization of shape, volume, and ionic positions using an 800 eV energy cutoff, and (2) optimization of cell shape and ionic positions at constant volume using a 400 eV energy cutoff. The results of these two tests are reported in Table 1. The Brillouin zone was sampled with the Γ-centered 6 6 6 k-points using a Monkhorst-Pack45 mesh. Both procedures produced similar values of cell parameters (see Table 1). The set of volume-energy data obtained in the second case was used
to fit the Murnaghan50 equation of state. As a result, the bulk modulus K and its first derivative K0 were obtained. Calculated properties are compared with experimental data in Table 1. To estimate the cohesive energy, Ecoh, the O and Si atomic energies were calculated at the Γ-point using a 10 10 10 Å cubic box in the appropriate spin state and stoichiometrically subtracted from the calculated R-quartz lattice energy. Agreement between theoretical and experimental values of Ecoh is excellent (Table 1). As expected from previous work,51 the electronic band gap using PBE is underestimated relative to the experimental value. Properties of bulk R-quartz have been calculated using quantum mechanical methods in previous studies.11,52-55 Results produced by different DFT techniques are extensively compared in refs 53-55. The PBE exchange-correlation functional leads to overestimation of the cell parameters by approximately 1-2% depending on the basis set and pseudopotential used.53-55 This result is consistent with our optimization results. Also, it has been shown53 that the 800 eV energy cutoff and a 6 6 6 k-point mesh are sufficient to obtain convergence for geometrical parameters within 0.1% if a PW basis set is used. All surface slab calculations were carried out with constant cell parameters obtained in the previous R-quartz crystal structure optimization. In most calculations, the minimum (primitive) surface UC has been used for the (101) R-quartz surface (i.e., a = 5.03, b = 7.467, γ = 109.682). This primitive cell corresponds to the conventional base-centered rectangular cell for space group C2. Models with three, four, and five layers of Si3O6 units (including 29, 38, and 47 Si and O atomic planes, accordingly) were considered; the third vector c was set to 20, 25, and 30 Å, correspondingly. The minimum slab is approximately 12 Å thick and has a vacuum gap of 8 Å. A 6 4 1 k-mesh was applied to these slab calculations. A number of calculations were performed using an extended surface cell. For this purpose, we prepared a slab consisting of a 2 2 supercell with dimensions 10.060 14.934 Å. In this case, the Brillouin zone was sampled with the Γcentered 3 2 1 k-points mesh. Whenever possible, symmetry constraints (space group C2) were used to create equivalence on both slab sides. Ionic relaxation was performed with the convergence criterion for the forces on atoms set to 0.01 eV/Å. For selected systems where practical, vibrational frequencies were calculated. In these calculations, the energy cutoff was fixed at 800 eV, and the HIGH precision level was used for other VASP options. The convergence criterion for the electronic cycle was set at 10-8 eV, and the positions of all atoms were optimized until the forces on atoms were smaller than 0.002 eV/Å. To calculate the Hessian matrix, finite central differences were used with 0.015 Å step size. 5757
dx.doi.org/10.1021/jp1106636 |J. Phys. Chem. C 2011, 115, 5756–5766
The Journal of Physical Chemistry C
ARTICLE
Table 2. Atomic Displacements (Å) Normal to the (101) Plane Relative to the Bulk Crystal Structure for the Two Upper Si3O6 Layers on a Symmetrically Reconstructed Surface Calculated Using the Five-Layer Slab Model (400 eV Cutoff)a first Si3O5 layer atom O2c Si1c
number of planeb -
second Si3O6 layer z-shift -0.27 -0.39
2 3
atom
number of planeb
z-shift
O7
10
0.17
O8 Si4
11 12
0.04 0.14
O3
4
-0.11
O9
13
-0.12
O4
5
-0.14
O10
14
-0.11
Si2c
6
-0.30
Si5
15
0.10
O5c
7
-1.02
O11
16
0.46
O6
8
0.45
O12
17
0.06
Si3
9
0.06
Si6
18
0.01
a
Displacements toward the vacuum are positive and toward the bulk are negative. b Sequence number of atomic plane in the cleaved unrelaxed slab starting from the top. O and Si atoms are numbered in the same order. c Atoms exhibiting large lateral displacements (>0.5 Å) during the surface 1 1 reconstruction due to formation of a twomembered Si ring (Figure 1b).
Figure 1. (a) Cleaved R-quartz (101) surface. (b) SiO2-terminated (101) relaxed surface for symmetric cleavage of R-quartz. (c) Adsorption of water molecule on the SiO2-terminated (101) surface of Rquartz. (View along the [010] direction). Red balls, O atoms; gray balls, Si atoms; blue balls, H atoms; black lines, UC boundaries.
To calculate H2O adsorption energies, the energy of an isolated H2O molecule was calculated in a 10 10 10 Å cubic cell. In these calculations, the Brillouin zone sampling was restricted to the Γ-point. For calculation of the isolated H2O vibrational frequencies, an 800 eV cutoff and 2 2 2 k-point mesh were used.
3. RECONSTRUCTION OF PRISTINE (101) SURFACE To estimate the formation energy of the hydroxylated surface, one needs to know the energy of the nonhydroxylated surface. For this purpose, calculations of the pristine (101) R-quartz surface were made. Pristine silica surfaces are unstable and will not survive in ambient conditions. Here, we consider the relaxation of pristine silica surfaces only for comparison to hydroxylated surfaces. Homolytic cleavage of R-quartz can produce two types of (101) terminations, symmetrical (equal composition of both created surface layers) and asymmetrical (different composition of both created surface layers). In the first case, the silyl Q3(þ) and siloxy Q3(O-) sites are produced in equal amounts, whereas in the second case, only one type of center is created on each surface. For both terminations, the pristine surface is unstable and undergoes fast reconstruction. Reconstruction of the symmetrically terminated surface, which is completed within 2 ps as shown by classical computer simulations,24 can lead to the formation of two-membered rings R2, [-O-Si(
