H-Bond Features of Fully Hydroxylated Surfaces of Crystalline Silica

Sep 18, 2009 - How Does Silica Catalyze the Amide Bond Formation under Dry ... Todd R. Zeitler , Jeffery A. Greathouse , Julian D. Gale , and Randall ...
6 downloads 0 Views 3MB Size
Download PDF
17876

J. Phys. Chem. C 2009, 113, 17876–17884

H-Bond Features of Fully Hydroxylated Surfaces of Crystalline Silica Polymorphs: A Periodic B3LYP Study Federico Musso,†,‡ Mariona Sodupe,*,† Marta Corno,‡ and Piero Ugliengo*,‡ Departament de Quimica, UniVersitat Auto`noma de Barcelona, Bellaterra 08193, Spain, and Dipartimento di Chimica IFM and NIS - Nanostructured Interfaces and Surfaces - Centre of Excellence, UniVersita` degli Studi di Torino, Via P. Giuria 7, 10125 Torino, Italy ReceiVed: June 6, 2009; ReVised Manuscript ReceiVed: August 6, 2009

Different models of hydroxylated surfaces of quartz, cristobalite, and tridymite have been studied with the hybrid B3LYP functional using the Gaussian basis set of polarized double-ζ quality with periodic boundary conditions. Starting from the optimized bulk structures of the polymorphs, 2D slabs exhibiting low (hkl) crystallographic planes have been cut, dangling bonds healed by hydroxyl groups, and the final structures fully optimized. The H-bond pattern at a given surface depends on the (hkl) plane and on the OH group density, exhibiting isolated, weakly interacting pairs, short chains, or strings extending through the whole surface. Cases in which no H-bonds are present envisage either a very low OH density or slab structural rigidity which hinders the OH groups to establish H-bond contacts. The thermodynamics of surface hydroxylation of the considered polymorphs has been shown to correlate with the strength of the H-bonds formed at the surfaces measured by the bathochromic shift of the ν(OH) stretching frequency with respect to the value for a free surface OH group. Simulation of the vibrational spectra in the OH stretching region for all surfaces of each polymorph showed a general good agreement with the experimental spectra recorded on polycrystalline powdered samples validating the present surface models for further studies on molecular adsorption. 1. Introduction Crystalline silica is one of the most abundant species on the earth’s crust, and its structural properties are involved in many different reactive and catalytic processes that are not yet completely clarified.1-3 For instance, in medicine, the toxic behavior of crystalline silica is well-known. Prolonged inhalation of micrometric dust of quartz and its polymorphs of major density, cristobalite and tridymite, provokes serious lung diseases like pulmonary fibrosis and silicosis.4-6 Moreover, it seems that the inhalation of crystalline silica dust could also provoke lung cancer, even if this aspect is still discussed because of the difficulties in interpreting experimental data. The pathogenic effect of silica on lung tissue is still unclear, but one of the hypotheses proposed by the scientists is that the adsorption of biomolecules on the surface could trigger the toxicity.7 Indeed, in many fields like materials science, pharmaceutics, chromatography, heterogeneous catalysis, or drug delivery, the surface properties more than the bulk features determine the behavior of silica. Surface properties are determined by the nature, the density, and the distribution of surface silanols (Si-OH) and siloxanes (Si-O-Si), the chemical functionalities that impart, respectively, the hydrophilic and hydrophobic character to the mineral. Silanols, which represent the most common termination of a silica surface, are generally involved in the adsorption of molecules which could be bound by H-bonds with their hydrogen acceptor sites, while siloxanes act in the adsorption of apolar molecules via weak dispersive forces. Different types * Corresponding author. E-mail: [email protected] (P.U.); [email protected] (M.S.). † Universitat Auto`noma de Barcelona. ‡ Universita` degli Studi di Torino.

of silanols can be found on a silica surface: vicinals or geminals, for example, that could be isolated or interacting through a H-bond (see Figure 1). Their distribution and density vary depending on the silica polymorph and surface, and also on the external factors and conditions like the previous treatment at which the material has been submitted or the temperature. An increase of temperature provokes the condensation of surface hydroxyls interacting via a hydrogen bond with liberation of a water molecule.8 Surface properties of silica have been investigated experimentally with different techniques, mainly with microcalorimetry and IR and Raman spectroscopy.9-16 The density and the distribution of the OH terminations determine different hydrogen bond patterns at the surface that are difficult to describe experimentally because of the presence of broad bands in the region of the spectrum corresponding to the OH stretching mode. Therefore, theoretical calculations within the cluster or the periodic approach are a powerful tool to describe the properties of these surfaces. Unfortunately, ab initio methods that include a systematic description of the electron correlation can only be used with relatively small systems because they are strongly time demanding. On the other hand, considering the classical approach, some computational studies have been done in order to investigate silica surfaces with potential functions,17 but most of the force fields developed for silica structures are parametrized only for the bulk properties and are not able to accurately describe H-bond interactions. In this field, an improvement has been made recently by some of us with the development of a new force field for silica (FFSiOH) that shows an excellent agreement with B3LYP data in the description of geometries, frequencies, and energetics of both bulk and surface properties of silica.18

10.1021/jp905325m CCC: $40.75  2009 American Chemical Society Published on Web 09/18/2009

Hydroxylated Surfaces of Crystalline Silica Polymorphs

J. Phys. Chem. C, Vol. 113, No. 41, 2009 17877 localized basis sets: a 6-31G*, very similar to the one we chose, and a larger one. They observe that, even if the structure is quite consistent between GGA and B3LYP, there are some differences in the description of H-bonds, with the H · · · O distance being slightly shorter or longer in the case of the hybrid functional, and in the computed stability of the surface with respect to the quartz bulk. H-bonds and energies seem to be dependent from the basis set as well. The comparison with their results allows validation of our structures and energetic stabilities and will be carried out later on in the paper. All of the calculations have been performed with the periodic CRYSTAL03/06 code and the B3LYP hybrid functional associated with the Gaussian basis set, which has been successfully tested in solid state simulations.34 2. Computational Details

Figure 1. Different types of surface silanols exhibited by the considered surfaces.

DFT seems to be the most appropriate choice to describe the surface properties of silica, taking into account the balance between the computational cost and the accuracy of the results. Previous DFT analysis of hydroxylated quartz and amorphous silica surfaces with periodic models has been performed with the local density approximation19 and, mainly, with BLYP20-23 and PW9124-27 GGA functionals, as they are common to many periodic codes. For β-cristobalite silica, the properties of hydroxylated surfaces have been investigated, within a periodic approach, with empirical potential functions28 as well as ab initio molecular dynamics with pure GGA functionals.29,30 In our opinion, the hybrid B3LYP functional shows a very good accuracy in the description of geometries, energies, and vibrational properties as well, while pure GGA functionals use to be less accurate in the description of the frequency spectrum in the OH stretching region when hydroxyls are noninteracting or involved in hydrogen bonds.31 The aim of this work is to describe structural and vibrational features of different fully hydroxylated surfaces of the three dense crystalline silica polymorphs, namely, quartz, cristobalite, and tridymite. The description of the distribution and density of surface hydroxyl groups as well as the description of H-bond interactions and their spectral signature is the first step to understand the adsorption properties of silica surfaces. The structure and stability of the (001) surface of R-quartz has been recently investigated by Goumans et al.32 with the B3LYP functional and an atomistic basis set, considering three periodic models of different thickness. In this work, the authors compare the results of the hybrid functional with those obtained with PW91 and PBE33 functionals and a plane wave basis set. They also compare the B3LYP results obtained with two

All calculations for bulk structures and surfaces have been performed with the ab initio CRYSTAL06 code.35 This code implements the Hartree-Fock and Kohn-Sham self-consistent field method for the study of periodic systems,36 and in the present version, it allows one to perform full geometry optimization, both internal coordinates and cell size,37 as well as to compute the phonon spectrum at the Γ point of molecules, slabs, and crystals.38 To define the surfaces, different 2D infinite slab models of variable thickness, cell parameters, and number of atoms inside the cell have been cut out from the bulk. It is worth noting that CRYSTAL03 models true 2D systems at variance with the approach followed by plane wave based codes in which the slab is artificially replicated through infinity also in the direction perpendicular to the slab by including a large amount of empty space. The unsatisfied Si or SiO valences of the surface derived from the cut were filled with OH groups or H atoms, respectively, using the molecular graphic program MOLDRAW,39 also used for the visualization and the geometrical manipulation of all the structures. 2.1. Basis Set. The multielectron wave function is described by linear combination of crystalline orbitals (COs) which, in turn, are expanded in terms of Gaussian type basis sets. The GTO basis set consisted of a 66-21G(d) (Rsp ) 0.13, Rd ) 0.5), 6-31G(d) (Rsp ) 0.2742, Rd ) 0.538), and 6-31G(p) (Rs ) 0.1613, Rp ) 1.1) for Si, O, and H, respectively (outer shell exponent in bohr-2). 2.2. Hamiltonian and Computational Parameters. For all calculations, Becke’s three-parameter (B3) hybrid exchange functional,40 in combination with the gradient-corrected correlation functional of Lee, Yang, and Parr,21 has been adopted. The DFT implementation is such that the exchange-correlation contribution is the result of a numerical integration of the electron density and its gradient, performed over a grid of points. In CRYSTAL, the Gauss-Legendre quadrature and Lebedev schemes are used to generate angular and radial points of the grid, respectively. A good compromise between accuracy and cost of the calculation for geometry optimization and vibrational frequencies has been shown by using a pruned grid consisting of 75 radial points and 1 subinterval with 974 angular points (75, 974).35,41 Default values of the tolerances that control the Coulomb (6 6 6 6) and exchange (12) series have been adopted so that when the overlap between two atomic orbitals is smaller than 10-6 (Coulomb) or 10-12 (exchange) the integral is disregarded.35 The Hamiltonian matrix has been diagonalized42 in 10 reciprocal lattice points (k-points), corresponding to a shrinking factor of 4.35

17878

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

2.3. Central Zone Phonon Frequencies. Complete phonon frequency calculation of quartz, cristobalite, and tridymite fully hydroxilated surfaces has been performed, with no restrictions to the OH fragments. Frequencies have been calculated as the eigenvalues obtained by diagonalizing the weighted Hessian matrix W associated to the Γ point (point k ) 0 of the first Brillouin zone, called the central zone, IR and RAMAN spectra refer to this point):

Wij(k ) 0) )

∑ G

Hij0G

√MiMj

Musso et al. TABLE 1: Structural Data and Relative Stabilities of r-Quartz, r-Cristobalite, and r-Tridymitea Geometry R-Quartz Expb B3LYP space group a (Å) b (Å) c (Å) V (Å3) 〈Si-O〉 (Å) 〈Si-O-Si〉 (deg)

P3121 4.902 4.902 5.399 37.4 1.612 142.6

P3121 4.971 4.971 5.467 39.0 1.632 142.7

R-Cristobalite

R-Tridymite

Expc

B3LYP

Expd

B3LYP

P41212 4.971 4.971 6.928 42.8 1.603 146.7

P41212 5.024 5.024 6.989 44.1 1.629 144.5

C2221 8.756 5.024 8.213 45.2 1.559 179.7

C2221 8.869 4.968 8.267 45.5 1.626 176.8

Relative Stability

In this expression, H0G ij is the second derivative of the electronic and nuclear repulsion energy E at equilibrium u ) 0 with respect to the displacement ui ) xi - xIi of atom A in cell 0 and displacement uj ) xj - xIj of atom B in cell G from their equilibrium positions xIi , xIj . In CRYSTAL, energy first derivatives Φj of E with respect to atomic displacements, Φj ) ∑G νGj ) ∑G (∂E/∂uGj ), are calculated analytically,43,44 while second derivatives are calculated numerically using a three-point formula:

[ ] ∂Φj ∂uj0

0

≈

Φj(0, ... , ui0, ... , 0) - Φj(0, ... , -ui0, ... , 0) 2ui i ) 1, ... , 3N; j ) 1, ... , 3N

Moreover, a value of uj ) 0.001 Å has been chosen as the displacement for each atomic coordinate and the tolerance of the convergence of the SCF cycles has been reduced from 10-6 (default value) to 10-11 Hartree. For all considered cases, it has been checked that the computed harmonic frequencies were all positive, ensuring that the slab models are all true minima on the B3LYP potential energy surface. 2.4. Geometry. Lattice constants and internal coordinates have been simultaneously optimized within the same run, using analytical gradients and upgrading the numerical Hessian with the Broyden-Fletcher-Goldfarb-Shanno algorithm.45-48 Convergence criteria are particularly stiff, since the localization of the minimum on the energy surface is crucial for calculating frequencies. For this reason, the tolerance for the maximum allowed gradient and displacement have been set to 0.0003 hartree Bohr-1 and 0.0012 Bohr, respectively. 3. Results and Discussion 3.1. Bulk Structures. Bulk structures of R-quartz49 (ICSD 93974), R-cristobalite50 (ICSD 77452), and R-tridymite51 with primitive hexagonal P3121, primitive tetragonal P41212, and centered orthorhombic C2221 space groups, respectively, have been obtained from the ICSD database52 and have been fully optimized, relaxing both internal coordinates and cell parameters. Vibrational frequencies of the bulk structures have been computed to obtain the values of enthalpy and compare the relative stabilities. Geometrical data and relative stabilities, both compared with some experimental values, are shown in Table 1. As expected from previous experience, the B3LYP hybrid functional slightly overestimates cell parameters and 〈Si-O〉 bonds while giving the same trend for the three polymorphs resulting from the experiment. Overall, the agreement between experimental and computational data for the three polymorphs is excellent. Even the B3LYP 〈Si-O-Si〉 angle inside the cell, whose value is

∆H (kcal/mol)

Expe

B3LYP

Expe

B3LYP

Expe

B3LYP

0.0

0.0

2.6

2.5

2.9

3.0

Cell parameters as a, b, and c, cell volume as V, 〈Si-O〉 and 〈Si-O-Si〉 average values. b Reference 49, ICSD# 93974. c Reference 50, ICSD# 77452. d Reference 51. e Reference 53. a

very sensitive to small perturbations due to its high flexibility, is in good agreement with experiments (largest deviation being 2.9°). As for energetic, calculated, and calorimetric53 ∆H values associated to the relative stability of the silica polymorphs, they are in excellent agreement with each other, validating the B3LYP functional as a reliable method to study silica materials. 3.2. Surface Models. All of the quartz, cristobalite, and tridymite surfaces studied in this work have been described with a two-dimensional slab model, by cutting a slab of a given thickness along a chosen crystallographic plane from the corresponding optimized bulk structure. This way, the slab model is defined by two surfaces that represent realistic silica surfaces. For quartz, the crystallographic planes have been chosen according to crystallographic information on the natural habit of the crystal structure. In the case of cristobalite and tridymite, however, due to the lack of precise experimental data on crystal habit, the studied planes correspond to those with lower Miller indexes, since the corresponding surfaces should, in principle, be the most developed ones due to their slower rate growth. To distinguish the various models, a capital letter specifying the polymorph (Q for quartz, C for cristobalite, and T for tridymite) together with the Miller indexes for the considered face have been adopted. As a whole, five surface models have been considered for quartz, Q(001), Q(100), Q(010), Q(011), and Q(101), four models for cristobalite, C(001), C(100), C(101), and C(110), and three for tridymite, T(001), T(100), and T(110), respectively. Cutting a silica structure means to break, heterolitically or homolitically, Si-O bonds, leaving highly energetic dangling bonds that must be saturated. Since the process of grinding silica materials occurs in a water rich atmosphere, the most convenient surface model is a fully hydroxylated slab. Depending on the Si-O bond cut, different types of surface silanols are generated. The most common types of SiOH are those described in Figure 1, i.e., geminals, vicinals, isolated and interacting via H-bond silanols. 3.3. Geometrical Features. All of the periodic slabs have been constructed maintaining the same hydroxyl pattern at the top and the bottom surface avoiding any symmetry constraints (P1 layer group) to ensure the maximum degrees of freedom during geometry optimization, except for the Q(010) slab, for which the P21 layer group symmetry has been enforced. Full

Hydroxylated Surfaces of Crystalline Silica Polymorphs

J. Phys. Chem. C, Vol. 113, No. 41, 2009 17879

TABLE 2: Deformation Data on Thickness (τ, Å), Area (A, Å2), and Volume (V, Å3) of the Slab Cellsa τb τa Ab Aa ∆A |∆A| (%) Vb Va ∆V |∆V| (%) a

Q(001)

Q(010)

Q(100)

Q(011)

Q(101)

C(001)

C(100)

C(101)

C(110)

T(001)

T(100)

T(110)

11.0 11.0 21.4 21.1 –0.28 1.3 235.4 232.3 –3.1 1.3

8.6 8.6 27.2 26.9 –0.24 0.9 233.7 231.6 –2.1 0.9

8.9 8.9 27.2 26.7 –0.52 1.9 241.9 237.2 –4.6 1.9

8.1 8.0 34.6 34.3 –0.30 0.9 280.2 274.3 –5.9 2.1

7.6 7.5 34.6 34.5 –0.10 0.3 263.0 258.7 –4.2 1.6

20.0 20.0 25.2 24.9 –0.38 1.5 504.8 497.1 –7.7 1.5

8.1 8.5 35.1 30.8 –4.4 12.4 284.4 261.3 –23.1 8.1

9.5 9.3 43.2 41.7 –1.5 3.5 410.8 388.2 –22.6 5.5

7.9 7.6 49.7 46.2 –3.5 7.0 392.3 350.8 –41.5 10.6

9.3 9.3 22.0 22.0 –0.04 0.2 204.9 204.5 –0.4 0.2

10.9 10.7 42.0 40.3 –1.77 4.2 458.0 430.7 –27.3 6.0

10.9 11.4 72.9 63.9 –9.0 12.4 794.2 727.9 –66.3 8.3

Subindexes b and a indicate geometrical values before and after the geometry optimization.

relaxation of cell parameters and atomic coordinates for the majority of the considered slabs brings the surface OH groups in contact via H-bond interactions, as will described in detail below. A critical point in modeling silica surfaces within the slab model is the influence of the slab thickness on the surface properties. Indeed, H-bond interactions between surface hydroxyls may cause relevant deformations of the internal slab structure which, in turn, are not realistic of the behavior of the real material in which surfaces are linked to a macroscopic bulk. In order to evaluate this effect, the deformation of each surface model has been analyzed by comparing the values of thickness, surface area, and volume of the crystalline 2D cells before and after geometry optimization. Results are summarized in Table 2. Considering that the position of the OH groups inside the cell changes because of the formation of hydrogen bonds in the optimized structure, the thickness has been calculated considering the distance within the highest and lowest silicon layer in the slab. The cell volume has been defined as the product of the slab cell surface area and the thickness. The difference between the volume of the cell before and after the optimization is a figure of merit of the structure deformation because it takes into account the variations of both cell surface area and slab thickness. Table 2 shows that for three structures the deformation exceeds the warning value of 5% of the cell surface area and 8% for the volume, while for the remaining slabs the deformation remains well under 3%. C(100) is the most critical case and needs further study due to the remarkable deformations in the internal structure of the model. In this case, the formation of H-bonds between surface OH groups induces a large variation in the Si-O-Si angles with a sizable rearrangement of the structure across the whole layer thickness. In order to evaluate if the flexibility showed by the C(100) model is a real feature of this surface or an artifact due to the layer thickness, slabs of increasing thickness have been optimized while keeping the very same surface OH termination. Figure 2 and Table 3 summarize the results. The considered slab models are labeled as C(100)T8, C(100)T10, and C(100)T13, respectively, where the indexes T8, T10, and T13 indicate the thickness of the slab in Å units. The increase of the thickness changes the features of the optimized surface because the rigidity of the structure is increased moving along the T8 f T10 f T13 series by hindering large deformations of the Si-O-Si angles responsible for the link between SiO4 tetrahedra. As revealed by Figure 2, the increased structural rigidity hinders the OH groups to come close enough and interact via H-bonds. Therefore, as shown in Figure 2, the passage from T8 to T10 results in a weakening of the surface H-bonds (from 1.9 to 2.2 Å, approximately), while the T10 to T13 passage entirely prevents the formation of H-bonds at the C(100) surface.

Figure 2. Hydrogen bond features of the C(100) surface as a function of slab thickness. Slab thickness and H-bond distances expressed in Å.

TABLE 3: Deformation Analysis of C(100)T8, C(100)T10, and C(100)T13 Slab Modelsa τb τa Ab Aa ∆A |∆A| (%) Vb Va ∆V |∆V| (%)

C(100)T8

C(100)T10

C(100)T13

8.1 8.5 35.1 30.8 –4.4 12.4 284.4 261.3 –23.1 8.1

10.6 10.9 35.1 32.9 –2.2 6.2 372.2 359.1 –13.1 3.5

13.1 13.2 35.1 34.5 –0.6 1.7 460.0 455.4 –4.6 1.0

a The thickness (τ, Å), area (A, Å2), and volume (V, Å3) of the slab cells are subindexed by labels b and a to refer to structures before and after geometry optimization, respectively.

The same analysis has been carried out on the other slab models exhibiting a cell volume deformation higher than 8%, i.e., the C(110) and T(110) slabs. The reason for these large deformations is due to larger crystalline unit cells which allow as many as eight H-bonds to be formed within a single surface

17880

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

Musso et al.

Figure 5. B3LYP optimized structures of the top hydroxylated R-tridymite surfaces. Hydrogen bond distances in Å as dotted line.

Figure 3. B3LYP optimized structures of the most representative top faces of hydroxylated R-quartz surfaces. Hydrogen bond distances in Å as dotted line.

Figure 4. B3LYP optimized structures of the top hydroxylated R-cristobalite surfaces. Hydrogen bond distances in Å as dotted line.

unit cell. Since C(110) and T(110) slabs already contain 42 and 72 atoms inside the unit cell, respectively, an increase of the thickness would imply, in both cases, an exceedingly high computational cost. A different strategy is to reoptimize only the atomic coordinates of both slabs, keeping the crystalline cell parameters frozen at their initial values. This approach simulates the rigidness of the bulk structure present below the surface in a real crystal and avoids the adoption of a thicker slab by means of imposing the unit cell resulting from the bulk. For both cases, however, at variance with the C(100) surface, the two models do not show the same sensitivity of the H-bond features on the slab thickness and rigidity. Indeed, only the atoms of the most external layers undergo sizable displacements due to the formation of H-bonds, the more internal atomic layers being almost unaffected. For the case of the C(110) slab, the optimization with the frozen cell parameters gives a final structure very close to the fully optimized one, with the same H-bond pattern envisaging a four-member ring on each face of the slab, although with H-bond distances somewhat larger. The formation energy per cell unit in the constrained system is -132.5 kJ/mol and the surface energy -0.22 J/m2, to be contrasted with the values of -153.5 kJ/mol and -0.27 J/m2 for the fully optimized structure. In conclusion, for C(110), relevant changes of the surface features with the slab thickness are not expected. A very similar behavior has been observed for T(110). The distribution and pattern of OH groups as well as the H-bond features of quartz, cristobalite, and tridymite models are shown in Figures 3, 4, and 5. Only the top surface of the slabs is reported, since the two surfaces are, in each case, almost identical. The geometrical features of the crystalline cells are summarized in Table 4 where the values of H-bond lengths have been averaged using top and bottom H-bond length values (difference between top and bottom is about 0.001 Å). Depending on the distribution of the hydroxyl groups on the surface, the optimized structures show different H-bond patterns: C(100) and T(001) models, with the highest slab rigidity and the latter with the smallest OH density as well, only show isolated surface silanols; Q(010), Q(011), C(001), and T(100)

Hydroxylated Surfaces of Crystalline Silica Polymorphs

J. Phys. Chem. C, Vol. 113, No. 41, 2009 17881

TABLE 4: Structural Data for Quartz (Q), Cristobalite (C), and Tridymite (T) Slabsa a (Å) b (Å) γ (deg) A (Å2) τ (Å) 〈Si-O〉 (Å) 〈O-H〉 (Å) F (OH/nm2) a

Q(001)

Q(010)

Q(100)

Q(011)

Q(101)

C(001)

C(100)

C(101)

C(110)

T(001)

T(100)

T(110)

4.868 4.971 119.2 21.1 11.0 1.633 0.980 9.5

4.841 5.565 90.0 26.9 8.6 1.635 0.970 7.4

4.972 5.361 90.0 26.7 8.9 1.632 0.980 7.5

5.000 7.300 110.1 34.3 8.0 1.636 0.968 5.8

10.202 4.983 137.3 34.5 7.4 1.636 0.972 5.8

4.969 5.002 90.0 24.9 20.0 1.631 0.973 8.0

5.019 6.873 89.5 34.5 13.2 1.631 0.965 5.8

8.300 5.030 90.0 41.7 9.3 1.633 0.974 4.8

6.856 6.733 90.0 46.2 7.6 1.637 0.979 8.7

5.078 5.076 121.5 22.0 9.3 1.629 0.964 4.5

5.080 7.923 89.9 40.3 10.7 1.630 0.971 5.0

8.030 7.951 90.1 63.9 11.4 1.632 0.979 6.3

Cell parameters (a, b, γ), cell area (A), slab thickness (τ), 〈Si-O〉 bonds and 〈Si-O-Si〉 angles, and the surface OH density (F).

models exhibit a single H-bond interaction inside the crystalline cell; the C(110) model shows four OH groups on the surface, interacting through a four-membered H-bond ring, whereas the Q(001), Q(100), Q(101), C(110), C(101), and T(110) models present H-bond chains extending through infinity. The present Q(001) structure is in excellent agreement with the model proposed by Goumans et al.,32 showing the same hydrogen bond chain with a short and a long H-bond. The most important difference is the length of the longest H-bond that changes by less than 0.1 Å, the other structural differences being negligible. Experiments capable of assessing the surface structural details at an atomic scale for silica polymorphs are rather scarce despite the enormous progress in the field.54 To our knowledge, only two recent papers are worth mentioning in that respect: the experimental determination of the structure, at atomic resolution, of thin hydrophobic silica film grown on a Mo(112) substrate by scanning tunneling microscopy (STM), infrared reflection adsorption spectroscopy, and X-ray photoelectron spectroscopy55 and the growth of an ordered ice layer on top of the same surface.56 For the latter case, low-energy electron diffraction (LEED) has been used, together with other spectroscopic techniques, showing the silica and ice superimposed lattices. In principle, the same techniques may be adopted to grow other crystal faces of the kind of those studied here, so that results of the modeling may provide excellent structures to be used in the STM and LEED refinement. 3.4. Energetics. In order to evaluate the stability of each slab model, an ideal reaction process of reconstruction of the surface, starting from the bulk plus water, is considered:

aBSiO2 + bH2O S cSSiOH

(1)

where BSiO2 refers to the bulk unit cell, SSiOH to the slab unit cell, and a, b, and c to the stoichiometric coefficients of the reaction. From the reconstruction energetics, the formation energy per unit cell and the surface energy (eqs 2 and 3) are derived: cell cell ∆Ecell ) cEslab - (aEbulk + bEH2O)

∆Esurf )

∆Ecell 2A

(2)

(3)

where A is the surface unit cell area. Values of ∆Ecell and ∆Esurf are reported in Table 5, which shows that almost all of the surface energies ∆Esurf are negative, resulting in a thermodynamically favored surface formation. Exceptions are for T(001) and C(100): for the former, both the

TABLE 5: Formation Energy of Quartz (Q), Cristobalite (C), and Tridymite (T) Slabsa Q(001) Q(010) Q(100) Q(011) Q(101) C(001) C(100) C(101) C(110) T(001) T(100) T(110) a

∆Ecell (kJ/mol)

∆Esurf (J/m2)

–89.3 –38.1 –87.1 –5.0 –24.8 –61.2 –0.3 –30.9 –153.5 2.7 –40.9 –183.7

–0.35 –0.12 –0.27 –0.01 –0.06 –0.21 –0.00 –0.06 –0.28 0.01 –0.08 –0.24

See text for the definition of ∆Ecell and ∆Esurf.

formation energy and the surface energy are positive, whereas, for the latter, both are close to zero. Interestingly, for these two surfaces, no H-bond is present among the surface silanols after the structure relaxation (see discussion above for structural details). That the H-bond between the surface silanols is the key factor driving the energetics of reaction 1 is also shown by data for the Q(011) slab resulting in the least stable case among the quartz surfaces, for which formation and surface energies are -5.0 kJ/mol and -0.01 J/m2, respectively. Indeed, Q(011) only presents a single and relatively weak H-bond (2.053 Å). Therefore, the thermodynamical stability of the slab and the formation of H-bonds between the surface silanols appear to be correlated. The surface energy ∆Esurf is the result of three different factors: (i) the rupture of Si-O-Si bonds of the bulk followed by the formation of surface silanols; (ii) the structural relaxation of the interior of the slab; (iii) the formation of H-bonds between the surface silanols. Process i is common to all considered surfaces and can be treated as a fixed contribution per Si-OH bond to the surface energy. Its local covalent nature can hardly justify any serious dependence of its value as a function of the considered slab. Processess ii and iii are in some way intermingled, but it becomes clear, from the above discussion, that H-bond patterns and strengths differ heavily depending on the nature of the considered silica polymorph and on the particular surface model. One way to study the possible link between H-bond strength and surface energy is to use the bathochromic shift, ∆ω(OH), of the OH stretching frequency, defined as the difference between the harmonic frequency of an OH group belonging to a given surface model with respect to the OH frequency value for an isolated silanol (in this case, the value of 3904 cm-1 of OH groups belonging to the T(001) surface, see Figure 5 and next paragraph), as a measure of the H-bond strength, as it is customary in the specific literature on H-bond features.57,58 Figure 6 shows such a comparison, in which only the largest ∆ω(OH) value among those computed for all OH

17882

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

Musso et al. width at half-maximum increasing as a function of the H-bond strength is simulated adopting the empirical relation:

Γi ) 0.72∆i + 2.5 cm-1 proposed by Huggins and Pimentel59 which allows one to compute the full width at half-maximum Γi from the B3LYP computed frequency shift, ∆i ) ωi(OH)[free] - ωi(OH)[Hbonded]. For each polymorph, a full spectrum based on the calculation of the OH harmonic stretching frequencies has been computed, by associating a Lorentzian function Li(V) ) Si(Γi/2π)/[(V ωi)2 + (Γi/2)2] with full width at half-maximum Γi, to each frequency value ωi. In this expression, Si is the B3LYP IR intensity associated to the frequency ωi. The IR spectrum for each silica surface R is then defined as Figure 6. Surface energies ∆Esurf vs the largest bathochromic harmonic OH frequency shift ∆ω(OH) calculated for all polymorph surfaces with respect to the value of 3904 cm-1 for the isolated OH groups of the T(001) surface.

groups of a given surface model has been considered, thereby emphasizing the role of the strongest H-bond as the most representative contribution to the surface energy. The correlation is good for quartz surface models except for Q(001) and Q(100); the same holds for cristobalite, with C(001) and C(110) lying out of the main correlation. Considering that geometrical relaxation within a given slab can be a subtle contribution to the final surface energy value, it is hard to assess which point between Q(001) and Q(100) or C(001) and C(110) should be considered the outlier with respect to the main correlation. The stability of Q(001) with respect to the reaction between bulk quartz and water is in excellent agreement with the value computed by Goumans et al.32 (-0.34 J/m2). Furthermore, the trend of the surface stabilities computed here is in agreement with that computed by Murashov25 for hydroxylated quartz surfaces, with periodic plane wave models and the GGA functional. Even in that case, the author found that the order of decreasing stability is Q(001) > Q(100) > Q(101), in line with the present results. 3.5. Frequencies. In order to validate the present quartz, cristobalite, and tridymite surface models, the IR spectra in the OH-stretching frequency region have been reconstructed by using the whole set of OH stretching frequencies computed for each surface. In this way, it is assumed that the experimental spectrum of a given polymorph is the result of the OH vibrating on all of the possible crystal faces exhibited by the polycrystalline powdered sample. A detailed comparison is unfortunately hampered by the rather broad bands observed in the experimental spectra, as resulting from the H-bond interactions between surface silanols. An important characteristic of H-bonds on silica surfaces is that silanols that act as H-bond donors suffer a bathochromic shift of their stretching frequency in comparison with the OH frequency associated to free silanols. The magnitude of the shift depends on the strength of the interaction and is correlated with the H-bond H · · · O length. Another well-known feature occurring on the IR spectra of real H-bonded systems is that OH bands shifted to lower wavenumbers also broaden due to anharmonic coupling with soft intermolecular modes. This effect cannot be computed within the harmonic approximation adopted here so that the full

HR(V) )

∑ Li(V) i

and each HR(V) contributes to the final spectrum G(ν) with equal weight:

G(ν) )

∑ HR(ν) R

In order to account for the deviation of the frequencies due to the harmonic model and limitations of the level of theory used, the B3LYP harmonic frequencies have been multiplied by a 0.959 scaling factor. This brings the calculated frequency value for the free OH vibrating at 3904 cm-1 on the T(001) surface (the one exhibiting isolated OH group) to the value of 3745 cm-1 which is considered the experimental value of a free surface OH vibration frequency. The very same procedure has been successfully used by some of us8 to reconstruct the infrared spectrum in the OH stretching region of models of amorphous silica fully hydroxylated surfaces, providing excellent agreement with the spectrum of a silica aerosol sample outgassed at 400 K. Figure 7 compares the computed spectra with experimental ones for quartz, cristobalite, and tridymite polycrystalline powders, pretreated at 150 °C, showing a general good similarity between the two, and validating the adopted crystalline surface models. The different H-bond patterns exhibited by the different surfaces as a function of the considered polymorph give rise to rather different infrared spectra. For quartz, only two frequencies at 3730/3721 cm-1 (from the terminal OH of H-bonded pairs belonging to top/bottom Q(011) surface, see Figure 3) contribute to the highest band. Two main features appear centered at about 3600 cm-1 (from all surfaces but Q(001)) and at 3350 cm-1 (from Q(001) and Q(100) surfaces) due to H-bonds of different strength. For cristobalite, six frequencies, respectively, at 3733/ 3721 cm-1 (terminal OH of H-bonded pairs belonging to the C(001) surface, see Figure 4), 3736/3712 cm-1, and 3709/3708 cm-1 (OH at the C(100) surface, see Figure 4) are responsible for the marked doublet at the highest frequencies. The broadband at lower frequencies is due to H-bond interactions occurring at the C(101) and C(110) surfaces and is composed of two main features centered at 3450 and 3350 cm-1, respectively. For tridymite, at least three components close to the 3745 cm-1 band assigned to the stretching vibration of a true isolated silanol group (H-bond free) resulted from the calculation: 3736/3715 cm-1 (from terminal OH of H-bonded pairs of the top/bottom

Hydroxylated Surfaces of Crystalline Silica Polymorphs

Figure 7. B3LYP spectra (continuous lines) in the OH stretching region of R-quartz, R-cristobalite, and R-tridymite compared with experimental spectra (dashed lines). B3LYP spectra have been translated to overlap with the experimental ones at 3800 cm-1 and the intensity scaled with respect to the intensity ratio for the highest frequency peak. Wavenumbers in cm-1.

T(100) face) and 3722 cm-1 (from isolated OH groups of the bottom T(001) face). The main broadband centered at 3456 cm-1 is due to the infinite H-bonded chains of the T(110) surface. Indeed, the presence of eight OH groups (top/bottom of the T(110) surface) in rather strong H-bond all almost vibrating around 3456 cm-1 (standard deviation 26 cm-1) resulted in a too intense band which may explain in part the discrepancy with the experiment. However, one should consider this comparison with all the limits of the adopted methodology and of the adopted surface models. Furthermore, an equal distribution of crystal faces is considered here, whereas the faceting of the real crystals is controlled by a variety of factors (thermodynamic, kinetic, impurities, etc.) not addressed by the present models. 4. Summary and Conclusions Twelve hydroxylated surfaces of quartz, cristobalite, and tridymite with low Miller indexes (hkl) have been described

J. Phys. Chem. C, Vol. 113, No. 41, 2009 17883 within a periodic approach using the B3LYP hybrid functional adopting a double-ζ polarized Gaussian basis set. The surfaces are characterized by different types of OH groups (geminal, vicinals, and isolated) and by different H-bonding patterns that result from the interaction between hydroxyl groups. The considered surfaces are all free from structural defects (strained two- and three-membered (Si-O) rings) associated to glassy or amorphous silica surfaces when outgassed at very high temperatures and imitate surfaces formed in a water rich atmosphere. The H-bond pattern strongly depends on the disposition and on the density of the OH group at the surface. The majority of the considered models exhibit infinitely connected H-bond chains, whereas small local H-bonded rings or simple pairs through completely free OH groups have also been found, as a function of OH density and rigidity of the silica framework. Since the surface models have been obtained by cutting 2D periodic slabs from the bulk structures of the crystalline silica polymorphs, the dependence of the results on the slab thickness has been validated by carefully checking the geometrical deformation induced in the interior of the slab by the surface reorganization. For that, the values of the thickness, the area, and the volume of the cell before and after the geometry optimization were compared. For just the C(100) surface model, the deformations suffered by the creation of the hydroxylated surfaces slab were so large that the surface H-bond pattern evolves from strongly bound chains for the thinnest slab to completely free OH groups for the thickest one. Other cases, for which the deformation was notable, change only very little the surface H-bond features as a function of the slab thickness. The surface energy for each slab is calculated considering the slab formation energy of the chemical reaction involving the silica bulk and water. Results indicate that the surface energy is mainly determined by the H-bond features occurring at the surface. In general, the stronger the H-bond interactions are, the smaller the reaction energy is, so that surface H-bond interactions dictate the thermodynamics of surface formation. Indeed, a correlation between the surface energy and the strength of the surface H-bond represented by the maximum bathochromic shift suffered by the OH groups compared to a free surface silanol has been shown to hold true. Finally, in order to validate the surface models for each polymorph, the IR spectrum in the OH stretching region was reconstructed by using properly scaled harmonic B3LYP frequencies and IR intensities, the empirical full width at halfmaximum for each band, and summing up these contributions for each considered surface. The comparison with the experimental IR spectra of quartz, cristobalite, and tridymite polycrystalline powders was fair, considering the very broad and featureless character of the experimental spectra caused by the H-bond. This result indicates that the proposed slabs constitute a reliable model of crystalline silica surfaces and a good starting point for the DFT simulation of adsorption processes or even surface reactions involving radical terminations, like those that have been suggested to take part in the development of pulmonary diseases provoked by crystalline silica. Acknowledgment. F.M. thanks the Ministerio de Ciencia e Innovacio´n for the FPU fellowship. Financial support from MICINN through the CTQ2008-06381/BQU project is gratefully acknowledged. Experimental spectra for the silica polymorphs have been kindly provided by B. Fubini, I. Fenoglio, and G.

17884

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

Martra, Dipartimento di Chimica I.F.M., Universita` degli Studi di Torino. Useful discussion with Prof. C. Pisani is also acknowledged. References and Notes (1) Scott, R. P. W. Silica Gel and Bonded Phases; John Wiley & Sons: Chichester, U.K., 1993. (2) Iler, R. K. The Chemistry of Silica; Wiley-Interscience: New York, 1979. (3) Nawrocki, J. J. Chromatogr. 1997, 779, 29. (4) Fubini, B.; Bolis, V.; Giamello, E.; Volante, M.; Costa, D. “Interaction between the surface of quartz and various substances mimicking the biological environment of the inhaled dust”; Mechanisms in Occupational Lung Deseases, Sebastien, P., Ed.; Colloque INSERM: Paris, 1991; Vol. 203, pp 317. (5) Fubini, B.; Volante, M.; Giamello, E.; Bolis, V. “The role of surface chemistry in the pathogenicity of mineral dust and fibres: some open questions”; Eighth International Conference on Occupational Lung Deseases, 1993, Genova. (6) Fubini, B.; Hubbard, A. Free Radical Biol. Med. 2003, 34, 1507. (7) Fubini, B. Health effects of silica. In The Surface Properties of Silicas; Legrand, A. P., Ed.; John Wiley & Sons: Chichester, England, 1998; pp 415. (8) Ugliengo, P.; Sodupe, M.; Musso, F.; Bush, I. J.; Orlando, R.; Dovesi, R. AdV. Mater. 2008, 20, 1. (9) Knozinger, H. In The Hydrogen Bond; Schuster, P., Zundel, G., Sandorfy, C., Eds.; North Holland: Amsterdam, The Netherlands, 1976; Vol. 3; p 1263. (10) Anedda, A.; Carbonaro, C. M.; Clemente, F.; Corpino, R.; Ricci, P. C. J. Phys. Chem. B 2003, 107, 13661. (11) Bolis, V.; Fubini, B.; Marchese, L.; Martra, G.; Costa, D. J. Chem. Soc., Faraday Trans. 1991, 87, 497. (12) Fubini, B.; Bolis, V.; Cavenago, A.; Ugliengo, P. J. Chem. Soc., Faraday Trans. 1992, 88, 277. (13) Fubini, B.; Bolis, V.; Cavenago, A.; Garrone, E.; Ugliengo, P. Langmuir 1993, 9, 2712. (14) Bolis, V.; Cavenago, A.; Fubini, B. Langmuir 1997, 13, 895. (15) Fenoglio, I.; Ghiazza, M.; Ceschino, R.; Gillio, F.; Martra, G.; Fubini, B. The role of nature and structure of surface sites in the biological response to silica particles. In NATO Science Series II: Mathematics, Physics and Chemistry. Surface Chemistry in Biomedical and EnViromental Science; Blitz, J. P., Gun’ko, V. M., Eds.; Springer: New York, 2006; Vol. 228, pp 287. (16) Fenoglio, I.; Martra, G.; Coluccia, S.; Fubini, B. Chem. Res. Toxicol. 2000, 13, 971. (17) de Leeuw, N. H.; Higgins, F. M.; Parker, S. C. J. Phys. Chem. B 1999, 103, 1270. (18) Pedone, A.; Malavasi, G.; Menziani, M. C.; Segre, U.; Musso, F.; Corno, M.; Civalleri, B.; Ugliengo, P. Chem. Mater. 2008, 20, 2522. (19) Rignanese, G.-M.; De Vita, A.; Charlier, J.-C.; Gonze, X.; Car, R. Phys. ReV. B 2000, 61, 13250. (20) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (21) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (22) Mishler, C.; Horbach, J.; Kob, W.; Binder, K. J. Phys.: Condens. Matter 2005, 17, 4005. (23) Masini, P.; Bernasconi, M. J. Phys.: Condens. Matter 2002, 14, 4133. (24) Perdew, J. P. Unified Theory of Exchange and Correlation Beyond the Local Density Approximation. In Electronic Structure of Solids ′91; Ziesche, P., Eschig, H., Eds.; Akademie Verlag: Berlin, Germany, 1991; p 11. (25) Murashov, V. V. J. Phys. Chem. B 2005, 109, 4144. (26) Yang, J.; Wang, E. G. Phys. ReV. B 2006, 73, 035406.

Musso et al. (27) Rignanese, G.-M.; Charlier, J.-C.; Gonze, X. Phys. Chem. Chem. Phys. 2004, 6, 1920. (28) Nangia, S.; Washton, N. M.; Mueller, K. T.; Kubiki, J. D.; Garrison, B. J. J. Phys. Chem. C 2007, 111, 5169. (29) Cesaroli, D.; Bernasconi, M.; Iarlori, S.; Parrinello, M.; Tosatti, E. Phys. ReV. Lett. 2000, 84, 3887. (30) Iarlori, S.; Cesaroli, D.; Bernasconi, M.; Donadio, D.; Parrinello, M. J. Phys. Chem. B 2001, 105, 8007. (31) Pascale, F.; Zicovich-Wilson, C. M.; Lopez-Gejo, F.; Civalleri, B.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 888. (32) Goumans, T. P. M.; Wander, A.; Brown, W. A.; Catlow, C. R. A. Phys. Chem. Chem. Phys. 2007, 9, 2146. (33) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (34) Dovesi, R.; Civalleri, B.; Orlando, R.; Roetti, C.; Saunders, V. R. Ab Initio Quantum Simulation in Solid State Chemistry. In ReViews in Computational Chemistry; Lipkovitz, K. B., Thomas, R. L., Cundari, T. R., Eds.; Wiley-VCH, Inc: New York, 2005; Vol. 21, pp 1. (35) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M. CRYSTAL06 0.0 User’s Manual; Universita´ degli Studi di Torino: Torino, Italy, 2007. (36) Pisani, C.; Dovesi, R.; Roetti, C. Hartree-Fock ab-initio treatment of crystalline systems. In Lecture notes in Chemistry; Springer: Berlin, Heidelberg, New York, 1988; p 193. (37) Civalleri, B.; D’arco, P.; Orlando, R.; Saunders, V. R.; Dovesi, R. Chem. Phys. Lett. 2001, 228, 131. (38) Zicovich-Wilson, C. M.; Pascale, F.; Roetti, C.; Saunders, V. R.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 1873. (39) Ugliengo, P. http://www.moldraw.unito.it. MOLDRAW: a molecular graphic program to display and manipulate molecular structures; H1 (32-bit) ed. Torino, 2005. (40) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (41) Pascale, F.; Tosoni, S.; Zicovich-Wilson, C. M.; Ugliengo, P.; Orlando, R.; Dovesi, R. Chem. Phys. Lett. 2004, 396, 308. (42) Monkhorst, H. J.; Pack, J. D. Phys. ReV. B 1976, 8, 5188. (43) Doll, K. Comput. Phys. Commun. 2001, 137, 74. (44) Doll, K.; Saunders, V. R.; Harrison, N. M. Int. J. Quantum Chem. 2001, 82, 1. (45) Broyden, C. G. J. Inst. Math. Appl. 1970, 6, 76. (46) Fletcher, R. J. Comput. 1970, 13, 317. (47) Goldfarb, D. Math. Comput. 1970, 24, 23. (48) Shanno, D. F. Math. Comput. 1970, 24, 647. (49) Tucker, M. G.; Keen, D. A.; Dove, M. T. Mineral. Mag. 2001, 65, 489. (50) Schmahl, W. W.; Swainson, I. P.; Dove, M. T.; Graeme-Barber, A. Z. Kristallogr. 1992, 201, 125. (51) Kihara, K.; Matsumoto, T.; Imamura, M. Z. Kristallogr. 1986, 177, 27. (52) NIST. FindIt ICSD Database 1.3.1; National Institute of Standards and Tecnology, 2003. (53) Saxena, S. K.; Chatterjee, N.; Fei, Y.; Shen, G. Thermodynamic Data on Oxides and Silicates; Springer-Verlag: Berlin, 1993. (54) Somorjai, G. A.; Park, J. Y. Surf. Sci. 2009, 603, 1293. (55) Weissenrieder, J.; Kaya, S.; Lu, J.-L.; Gao, H.-J.; Shaikhutdinov, S.; Freund, H.-J.; Sierka, M.; Todorova, T. K.; Sauer, J. Phys. ReV. Lett. 2005, 95, 076103. (56) Kaya, S.; Weissenrieder, J.; Stacchiola, D.; Shaikhutdinov, S.; Freund, H.-J. J. Phys. Chem. C 2007, 111, 759. (57) Jeffrey, G. A. An Introduction to Hydrogen Bonding (Topics in Physical Chemistry), 1st ed.; Oxford University Press: New York, 1997. (58) Gilli, G.; Gilli, P. J. Mol. Struct. 2000, 552, 1. (59) Huggins, C. M.; Pimentel, G. C. J. Phys. Chem. 1956, 60, 1615.

JP905325M