J. Phys. Chem. B 1999, 103, 2165-2171
2165
Quantum Mechanical ab Initio Characterization of a Simple Periodic Model of the Silica Surface B. Civalleri, S. Casassa, E. Garrone, C. Pisani,* and P. Ugliengo Department of Inorganic, Physical, and Materials Chemistry, UniVersity of Torino, Via Giuria 5, I-10125 Torino, Italy ReceiVed: July 8, 1998; In Final Form: October 23, 1998
The electronic, structural, and vibrational properties of a thin silica film are studied by means of periodic ab initio techniques, both Hartree-Fock and density functional. The film is cut from bulk edingtonite, a tetragonal silica structure with five SiO2 groups per unit cell; the dangling bonds are saturated so as to obtain a fully hydroxylated surface. This model system is proved to possess a number of attractive properties which make it suitable for the study of some properties of silica surfaces: (1) the parent structure is rather stable, its energy per unit SiO2 being comparable to that of common zeolites such as faujasite; (2) the formation of the surface involves low energy, no reconstruction, and scarce relaxation; (3) the simplest two-dimensional structure contains only five SiO2 groups in the unit cell and can therefore be treated with limited computational effort; (4) the building block of the two-dimensional periodic structure is easily related to standard cluster models of silica; (5) the density of surface hydroxyls is comparable to that experimentally observed on dehydrated silica; (6) the stretching frequency of the OH group at the surface is very close to that obtained with the best cluster models and to the experimental value for isolated hydroxyls on amorphous silica; (7) a few variants of the fundamental structure are possible, and different low-energy surfaces can be cut from them, which allow the simulation of a number of local situations within the same basic model.
1. Introduction There exists a variety of silica surfaces both amorphous and crystalline which are involved in important reactive and catalytic processes; the theoretical simulation of these systems is therefore a subject of active research.1 If we limit ourselves to considering quantum mechanical (QM) ab initio approaches, cluster models are by far the favorite technique due to their versatility and simplicity; through the use of standard computer codes of quantum chemistry and the clever design of the molecular clusters, a lot of valuable information has been collected.2 On the other hand, the accurate QM treatment of clusters comprising more than a few SiO2 groups is practically impossible; Si24O60H24 to simulate a sodalite cage is among the largest model clusters that have been considered.3 The need for more realistic models which take into account solid state effects on the chemical and physical properties of the atoms at the surface is evident, at least for the purpose of calibrating the cluster models. Vigne´-Maeder and Sautet (VS) have recently proposed an extended periodic model of hydroxylated and dehydroxylated surfaces of silica which partly responds to this requirement.4 Following a suggestion formulated by Sindorf and Maciel for the interpretation of 29Si NMR experiments on amorphous silica surfaces,5 they considered the (111) and (100) surfaces of β-cristobalite. Slabs parallel to these surfaces were cut from the perfect crystal of that substance, and the resulting twodimensionally periodic structures were studied at the HartreeFock (HF) level by means of the CRYSTAL92 program;6,7 this periodic simulation was supported by molecular cluster calculations. The slabs had to be very thin (5 atomic layers) in order to keep the computational burden within reasonable limits, since the number of atoms in the unit two-dimensional cell is proportional to film thickness. As a consequence, the films
which resulted were formed by a stack of parallel and disjoint polymers of infinite length, a highly unstable structure whose resemblance to the parent crystal is at least dubious. The VS study indicates however a feasible computational scheme for the simulation of silica surfaces: (i) choice of a suitable crystalline structure of silica; (ii) extraction of thin films parallel to low-index crystal faces; (iii) saturation of the dangling surface bonds with oxygen and hydrogen atoms so as to generate a realistic surface situation; (iv) use of the two-dimensional periodic structure of the film to calculate its properties; (v) comparison with molecular cluster results. This procedure has been followed with success for investigating the surface reactivity of different systems, such as zirconia;8 the VS paper shows how difficult it is to apply it in the case of silica. The present study aims at the same objective and uses essentially the same techniques, but the starting point has been different. While trying to design cluster models of silica whose size could naturally and consistently be increased, a simple “building block” was devised which can form repetitive structures in the three orthogonal directions;9 it will be designated in the following as Y, and its structure is shown in Figure 1. Its chemical composition is (SiO2)5, but it is convenient to regard the terminal atoms of the central (SiO)4 ring (six oxygens and one silicon) as symmetrically divided between the two ends of the ring; in this form, it possesses a 4-fold rotoreflection axis S4 along z. As shown in the same figure, Y is the basic unit of three naturally occurring fibrous zeolites: edingtonite (EDI), thompsonite (THO), natrolite (NAT).10,11 Cutting these structures parallel to the (100) or (001) crystal plane allows us to generate a variety of thin films characterized by different surface topology and closely related to the parent three-dimensional crystals and to model clusters obtained from
10.1021/jp9829187 CCC: $18.00 © 1999 American Chemical Society Published on Web 03/09/1999
2166 J. Phys. Chem. B, Vol. 103, No. 12, 1999
Figure 1. Schematic representation of the building block Y, and of the derived structures edingtonite (EDI), thompsonite (THO), and natrolite (NAT). In Y, the silicon atoms which form the four rings are in black; white atoms do not belong to the repeat unit.
the combination of Y blocks. This point will be dealt with in more detail in section 2. In the present work, we will consider only the simplest of these surface structures, a slab cut from edingtonite parallel to the (100) face, to be labeled EDIxSn, x indicating the direction perpendicular to the surface, S standing for slab, and n indicating the number of Y units perpendicular to the slab. The surface oxygens are saturated with hydrogens as discussed in the next section. The number of atomic layers in the slab is then 4n + 3. Most of the calculations will concern the case n ) 1. An important reference for the present study is a recent paper of ours (to be referred to in the following as I) concerning the structure and stability of six silica polymorphs: R- and β-quartz, sodalite, chabazite, faujasite, and edingtonite.12 CRYSTAL95,13,14 a more recent release of the code employed by VS, was used for calculating the equilibrium geometry and total energy of those bulk crystals and for exploring the dependence of the results on the basis set and on the Hamiltonian adopted. As concerns the latter point, the HF calculations were supplemented with a posteriori estimates of the correlation energy;15-17 furthermore, a variety of density functional15,18,19 and hybrid20 schemes were considered by exploiting the new capabilities of the code.21 The most important results of I which are relevant for the present study can be summarized as follows. (1) An accurate comparison of different basis sets has permitted us to select one among them (Si[6-21G(d)], O[6-31G(d)]) that, while relatively simple, gives results close to those obtained with much more extended sets. The two d exponents are 0.5 and 0.6 Bohr-2 for silicon and oxygen, respectively. This set is used here in conjunction with a 31G(p) set for hydrogen (see next section). (2) The relative stability of the different polymorphs is similar, almost independently of the basis set and of the Hamiltonian adopted; R-quartz and edingtonite are at the two extremes of the scale, their stability differing by about 11 kJ/mol per SiO2 unit. For comparison, faujasite is less stable than R-quartz by about 7 kJ/mol. (3) The optimization of the geometrical parameters was performed in two steps: the GULP code developed by J. Gale and based on shell-model ion pair potentials22,23 provided a reference geometry. The final opti-
Civalleri et al. mization was performed by means of a modified Polak-Ribie`re algorithm24 using the QM energies obtained with the CRYSTAL code. In fact, the use of suitable ion pair parameters25,26 made the GULP optimized geometry very close to the final result in all cases; the corresponding energy differed from the optimized value by 0.5 kJ/mol on average. This has encouraged us to use the semiclassical method for optimizing the geometry of the two-dimensional structures here considered, except for the terminal hydroxyl groups whose configuration was determined through a QM optimization. We adopt here the same methodology as in I and take as a reference the results obtained there for bulk edingtonite. We also consider results obtained with cluster models of the surface, either the simple silanol (SIL) molecule, H3SiOH, or more complicated ones, based on the Y building block. In section 2, we present the model systems which have been selected for an explicit treatment and a variety of others which could deserve attention for further investigations. Some aspects of the computational techniques are recalled in section 3. Section 4 presents and discusses the results. The structural effects and the energetic cost associated with the formation of the slab from bulk edingtonite are analyzed in order to verify if the two-dimensional structure may represent an acceptable model of the surface of amorphous silica. The electronic structure of the slab is compared with that of the bulk and of cluster models of the surface. Special attention is given to the electrostatic field above the slab as compared to that in the corresponding region at the clusters, because the difference between the two models comes to the foreground when considering this important observable. The topological analysis by Bader27,28 is used to gain information on the nature of the different Si-O bonds in our model systems. The calculated stretching frequency of the terminal O-H bond is finally discussed with reference to experimental data, as a check of the “normal” character of the surface. 2. Model Structures The Y building block of formula (SiO2)5 on which our model systems are based is shown in Figure 1, both by explicitly representing all its atoms and in a schematic symmetrized representation of its silicon-only skeleton. A few derived structures are shown in the same figure. The Y blocks may be connected to each other along z to form chains (the so-called 4-1 chains29), which are a prominent component of fibrous zeolites.30 In turn, chains can be joined among them along x and y (the two equivalent directions) either by a simple translation (“parallel” interchain bond) or by translation followed by rotation by π /4 (“antiparallel” bond). In edingtonite, all chains are connected through parallel, in natrolite through antiparallel bonds. Thompsonite is intermediate between the two, since parallel bonding occurs in one direction (x, say), antiparallel in the other; as a consequence, thompsonite no longer possesses a 4-fold axis but only two orthogonal symmetry planes. A more complete analysis of the different crystalline structures which can be obtained from a combination of 4-1 chains has been performed by different authors.29,31 By cutting the interchain bonds in a (100) or (010) plane of these crystals, surfaces are formed of the type shown in Figure 1. If the crystals are cut perpendicularly to the chains, that is, (001) surfaces are considered, more complex and probably less stable twodimensional structures are obtained: according to which the atomic layer is left at the top of the slab and after saturation of dangling bonds with hydrogens, one can generate a surface rich in geminal or vicinal hydroxyls, or both. As anticipated in the Introduction, we consider here EDIx Sn slabs, cut from bulk edingtonite parallel to the (100) face
Simple Periodic Model of the Silica Surface
J. Phys. Chem. B, Vol. 103, No. 12, 1999 2167 TABLE 1: Geometry Data for Bulk Edingtonite and for the Slabs EDIxS1 and EDIxS2a system bulk
EDIxS1 EDIxS2
technique
b
c
(Si‚‚‚Si)x
(Si‚‚‚Si)y
GHF HF GB3 B3 GHF GB3 GHF GB3
6.955 6.917 7.001 7.000 6.947 6.991 6.949 6.995
6.474 6.453 6.492 6.492 6.468 6.487 6.474 6.490
3.765 3.731 3.780 3.773 3.794 3.816 3.779 3.798
3.765 3.731 3.780 3.773 3.754 3.768 3.757 3.774
a
The optimization is performed with respect to energies calculated with different techniques: GHF (GULP with “Hartree-Fock” model potential parameters); GB3 (GULP with “B3-LYP” model parameters); HF (CRYSTAL ab initio HF); B3 (CRYSTAL ab initio B3-LYP). b and c are the cell parameters in the y and z directions (for the bulk, a ) b); (Si‚‚‚Si)x and (Si‚‚‚Si)y are the distances between silicon atoms at the opposite ends along x and y of the (SiO)4 ring of the Y block. All entries are in angstroms.
TABLE 2: Geometry Data for the Terminal Si-O-H Group Obtained Using the CRYSTAL Code with ab Initio HF (HF) or B3-LYP (B3) Techniques, as Indicated (Angles in Degrees, Distances in Å) system SIL A B Figure 2. Structure of EDIxS1, EDIxS2, and of two related cluster models, A and B. In the two slabs, the numbers attached to silicon atoms are the NMR chemical shifts [-δTMS(Qn)] estimated according to Sauer’s formula43 (see text section 4.A); the experimental values are reported in italics.
and containing n ) 1 or 2 layers of 4-1 chains perpendicular to the slab (see Figure 2). By choosing to terminate the slab on both sides with oxygens, which are then saturated with hydrogens, a symmetry plane is maintained at its center. The resulting structure presents a density of 2.2 isolated hydroxyls per square nanometer, close to that experimentally observed on amorphous silica dehydrated at 800 K (1-2 O-H per nm2).32 Starting from the same building block, a variety of clusters can be generated. The two which have been considered here are represented in Figure 2, are identified as A (Si5O7H8) and B (Si6O10H8), and are cut from the EDIxS1 slab by saturating the silicon unfilled valences with hydrogens. 3. Computational Techniques The Si[6-21G(d)] and O[6-31G(d)] basis sets adopted for silicon and oxygen are the same as in I (see Introduction). The standard 31G(p) set was used for hydrogen; the exponent coefficient of the p polarization function is 1.1 Bohr-2. The QM codes GAUSSIAN-9433 and CRYSTAL9513 were used for the clusters and the periodic systems, respectively. Both for the clusters and for the periodic structures we considered two kinds of Hamiltonian, Hartree-Fock and hybrid (LYP functional16 plus HF exchange, according to the Becke formula B3-LYP20); the two schemes will be referred to in the following as HF and B3. The slab geometries were optimized by means of the GULP code22,23 using the parameters designed by Sauer and co-workers with reference to HF25 and B326 computations; the corresponding results will be indicated as GHF and GB3. In a final step, only the four parameters defining the geometry of the terminal Si-
EDIxS1 EDIxS2
technique
Si-O
∠SiOH
O-H
HF B3 HF B3 HF B3 HF B3 HF B3
1.639 1.659 1.611 1.629 1.610 1.628 1.611 1.632 1.611 1.632
119.7 116.9 119.8 117.2 120.0 117.3 122.4 116.6 123.5 116.6
0.942 0.963 0.942 0.963 0.942 0.963 0.942 0.964 0.944 0.964
O-H group were optimized by performing periodic HF and B3 calculations and adopting the same techniques as in I. On the contrary, the geometry of the clusters was fully optimized, by imposing Cs symmetry for SIL and A and C2V symmetry for B. The topological analysis of the electron density was performed by means of the program TOPO developed by C. Gatti and interfaced to CRYSTAL.34,35 The anharmonic O-H stretching frequency has been computed at the HF and B3 level by the following procedure. (i) The O-H distance has been assumed as a normal mode decoupled with respect to all other modes, in agreement with previous work.36-38 (ii) This distance was allowed to vary around its equilibrium value and the potential energy calculated for each distance. (iii) A polynomial curve of sixth degree was used to best-fit the energy points, ensuring a root-mean-square error below 10-6 hartree. (iv) the corresponding nuclear Schroedinger equation was solved numerically following the method proposed in ref 39, by means of the program ANHARM.40 4. Results and Discussion A. Geometries. Tables 1 and 2 report data concerning the optimized geometries; the former displays the relaxation of the silica framework following the formation of the surface, the latter, the dependence of the hydroxyl geometry on the underlying silica structure.41 In all cases, different optimization schemes are compared. Consider first the bulk results. The semiclassical optimization is in excellent agreement with the B3 QM result, while the agreement is not quite as good for the HF type calculations,
2168 J. Phys. Chem. B, Vol. 103, No. 12, 1999 possibly because of the different procedure followed for defining the model potential parameters. Indeed, the B3 force field has been tuned on “mechanically embedded” clusters,25 whereas the older HF parametrization resulted from “free-space” molecular clusters.26 As already shown in I, the B3 structure appears slightly more expanded than the HF one, the difference being however well within 1%, in agreement with previous findings about the sodalite structure.25 The comparison of bulk with slab geometry clearly shows that relaxation is almost negligible, an indication that the twodimensional structures are quite stable. On closer inspection, we can notice that the thinner slab exhibits a certain in-plane contraction, while the distance between silicon atoms of the (SiO)4 ring in the orthogonal (x) direction is increased. All these effects are considerably reduced in the thicker EDIxS2 slab. It can be of interest to compare interchain Si‚‚‚Si distances, for instance calculated according to the GB3 scheme. In the bulk, they are given by dint ) b - (Si‚‚‚Si)y ) 3.221 Å. They are practically unchanged in the slabs, 3.223 and 3.221 Å in a direction parallel to the surface and 3.220 Å in the perpendicular direction for the thicker slab. When considering the equilibrium geometry of the Si-O-H group in the y ) 0 plane, two minima are found, either the proton is bent in the positive or in the negative direction of the z axis (see Figure 1). The two minima are practically equivalent, although the second one is slightly more stable; the data of Table 2 are referred to this case. The O-H bond length is unchanged in all cases at both HF and B3 level, while the Si-O bond is appreciably longer in SIL than in the other structures. The SiO-H angle computed at the HF level appears to increase regularly along the series, from 119.7 for SIL to 123.5 for the thicker slab; this contrasts the B3 result where the angle remains almost constant. Again, B3 bond lengths are consistently larger than HF ones, but the difference is within 1-2%. From the knowledge of the calculated equilibrium geometry, predictions can be formulated about the 29Si NMR magnetic shielding constants of the different silicon atoms. Fifteen years ago, Engelhardt and Radeglia42 proposed a semiempirical formula to correlate the chemical shift of Q4 atoms (silicons with four Si-O-Si angles) to the observed values of the four Si-O-Si angles. This relationship has been recently generalized by Sauer and co-workers43 to include Qn atoms with n e 4 (silicons with n Si-O-Si and 4 - n Si-O-H angles). The formula reads
δTMS(Qn) ) -61.23
∑′F(SiOSi) - 60.45 ∑′′F(SiOH) + 6.58
where F(SiOX) ) cos(SiOX)/[cos(SiOX) - 1] δTMS is the chemical shift in ppm with respect to tetramethylsilane, the sums ∑′ and ∑′′ run over the n Si-O-Si and the 4 - n Si-O-H bonds, respectively, and the SiOX angles are referred to the results of HF geometry optimizations rather than to experimental observations. Using this formula, the value of -δTMS for R-quartz silicon, calculated with reference to the HF optimized geometry reported in I, is 105.7, to be compared to the experimental value of 107.4. The calculated values of -δTMS for the different Si atoms of EDIxS1 and EDIxS2, obtained using the optimized GHF angles, are reported in Figure 2. The observed values in amorphous silica are 99 and 109 for Q3 and Q4 atoms, respectively.5 Again, the agreement is excellent. B. Energetics. The total energy data per unit cell of the periodic optimized structures are as follows (the two entries, in
Civalleri et al. hartree, refer to HF and B3 computations in this order): bulk, -2194.40966, -2200.17020; EDIxS1, -2270.42749, -2276.54768; EDIxS2, -4464.83697, -4476.71745. The energies of an isolated water molecule with the same basis set are -76.02311 and -76.38097 hartree at HF and B3-LYP levels, respectively. These data allow us to gain information about the stability of the three-dimensional structure with respect to exfoliation in the presence of water. In the following, we consider only B3 energies and make reference to unit twodimensional cells. The reaction
EDIxS2 + H2Ovapor f EDIxS1 + 8.1 kJ/mol can be interpreted as the process of forming two hydroxylated surfaces by breaking the interlayer bonds in EDIxS2. The reaction is endothermic and corresponds to a surface formation energy of 0.0148 J/m2. The formation of EDIxS2 from the bulk is more endothermic.
2EDIbulk + H2Ovapor f EDIxS2 + 9.2 kJ/mol It should provide a better estimate of the surface formation energy of the semicrystal (0.0168 J/m2). The higher value of this reaction energy with respect to the previous one could indicate that interlayer bonds are weakened in the final phases of the exfoliation. In all events, the energy required for the formation of a unit area of the hydrated surface appears very small, and exfoliation should correspond to an easy process, even at moderate temperatures. More precisely, hydroxylated fibers parallel to z should be formed since the bonds in the x and y directions are equivalent. In fact, in this type of zeolites, the most extended naturally occurring crystal faces are those parallel to the 4-1 chains.10 It is interesting to compare the above results with those concerning the reaction with water of a cluster B-B (two B clusters connected through an interchain-type Si-O-Si bond). From the corresponding B3 energies (-4913.35287 and -2494.88810 hartree for B-B and B, respectively) it turns out that the breaking of the bond is in this case an exothermic process (∆E ) -19.1 kJ/mol). This confirms the hypothesis that interchain bonds are made stronger by cooperative effects in bulk edingtonite. C. Electronic Structure. To characterize the electronic structure of the different systems, it is expedient to consider different observable quantities such as the electrostatic potential and the electron density or its derivatives. As concerns the latter quantity, the trajectories of the density gradient are particularly interesting because they permit identification of the atomic basins and the bond critical points; the value of the density and of the density Laplacian there can provide information on the character of the bonding interaction27,28 (see below). None of these representations reveals significant differences between the two QM approaches, B3 or HF. The increase of interatomic distances in the B3 with respect to the HF equilibrium structure does not bring in qualitatively new features. Unless explicitly stated, we shall make reference to B3 results. Figure 3 shows the electrostatic potential V(r) above EDIxS1 and EDIxS2, and at the corresponding regions of the clusters. The presence of a second layer of 4-1 chains does not alter the potential at the surface, to confirm the general finding that very thin films can adequately reproduce the electrostatic potential at low-index surfaces of ionic crystals.44 On the contrary, the difference between the description of this quantity by periodic and cluster models is striking. This is even more evident in Figure 4, which displays the electrostatic field along
Simple Periodic Model of the Silica Surface
J. Phys. Chem. B, Vol. 103, No. 12, 1999 2169
Figure 4. B3 electrostatic field along the direction of the OH bond, for different systems, as indicated. The field is expressed in atomic units; multiply by 1033 to obtain the field in kJ/(mol debye).
Figure 3. B3 electrostatic potential at the surface of different structures, as indicated. Consecutive lines differ by 0.01 au (0.271 V); lines corresponding to absolute values larger than 10 times this step are not plotted. Continuous, dashed, and dot-dashed curves refer to positive, negative, and zero potential, respectively.
the direction of the OH bond in a region of distances of chemical interest from H. The curves for the two slabs are perfectly superimposed (only that for EDIxS2 is reported for clarity), while the clusters generate electrostatic fields which are stronger by 4-7 kJ/(mol debye) in the range 1.8-2.2 Å. One would expect from the present result that finite cluster models tend to overestimate the interaction with polar or easily polarizable molecules. Instead, the observed heat of adsorption of ammonia on amorphous silica is consistently larger than that obtained with the use of cluster models of the type here considered.45 Proper account must be taken, however, of the fact that the direction explored in Figure 4 is not representative of the complex structure of the electrostatic field; as is easily seen from inspection of Figure 3, the periodic model generates stronger fields in a direction parallel to the surface, at distances which are still accessible to small molecules. Recent work of ours46 concerning the adsorption of ammonia on the different systems here considered supports this
Figure 5. Value of the Laplacian of the B3 electron density at critical points along Si-O bonds as a function of bond length. Squares refer to EDIxS1, circles to B densities; “surface” labels silicon-hydroxyl bonds.
interpretation. Although the molecule comes closer to the hydroxyl of the A cluster (1.771 Å) than to that of EDIxS1 (1.836 Å), the adsorption energy on the latter system is larger by almost 7 kJ/mol. This is seen to be due to the additional interaction that takes place in the periodic model between two hydrogens of the ammonia molecule and the oxygen atom of the hydroxyl group at a neighboring cell; correspondingly, the equilibrium geometry is considerably different for the two systems. Let us now make a few comments on the results of the topological analysis. Rather than concentrating on atomic quantities, such as atomic charge, dipole, or energy obtained from the integration of suitable functionals over the atomic basin and requiring very heavy computations to provide reliable results, we consider here only the critical points CSi-O along
2170 J. Phys. Chem. B, Vol. 103, No. 12, 1999
Civalleri et al.
TABLE 3: Calculated Frequency Data for the O-H Stretching (cm-1)a system SIL A B EDIxS1 EDIxS2 amorphous silica a
calculation method
ωe
ω01
ω02
ωexe
HF B3 HF B3 HF B3 HF B3 HF B3 (exptl47)
4228 3889 4236 3904 4234 3908 4244 3895 4259 3902 3922
4076 3730 4085 3748 4085 3750 4095 3740 4106 3743 3749
7999 7300 8021 7340 8022 7341 8039 7324 8058 7326 7325
76 80 75 78 75 79 75 78 77 80 87
Symbols are explained in the text.
the different Si-O bonds in EDIxS1 and B. Their position, the value of the density and of the density Laplacian ∇2F there permits us to characterize the bond as more or less ionic. In particular, positive and high values of the Laplacian at a critical point of type (3,-1) (a minimum along the bond, a maximum in the two orthogonal directions) are taken as indicative of a markedly ionic character of the bond. Figure 5 reports these values as a function of bond length. An obvious anticorrelation is observed between the two quantities. In particular, the interchain bond has by far the largest value of the Laplacian and the shortest length; for the rest, cluster and slab bonds are concentrated in a narrow range of values of the two quantities, although the cluster generally presents slightly longer and less ionic bonds. The “surface” Si-O(-H) bonds are located above the correlation line and appear more ionic. D. Vibrational Frequencies. The procedure for solving the nuclear Schroedinger equation for the OH vibration has been presented in the previous section. The first two vibrational transitions, namely the fundamental ω01 and the overtone ω02, the harmonic ωe, and the anharmonicity constant ωexe of the O-H mode are obtained from the eigenvalues E0, E1, and E2 following the scheme
ω01 ) E1 - E0 ω02 ) E2 - E0 ωe,xe ) (2ω01 - ω02)/2 ωe ) ω01 + 2ωexe These quantities are reported in Table 3. The ωe and ω01 values for SIL are always the lowest. HF frequencies increase regularly along the series, while at the B3 level they are appreciably lower for the slabs than for the A and B clusters. In all cases, the spread is very low, about 30 cm-1; the anharmonicity, ωexe, is remarkably constant and almost insensitive to electron correlation. A spectacular agreement is shown in Table 3 between B3 results and experimental data from infrared measurements on amorphous silica.47 This can be partly fortuitous, since it is known that calculated frequencies have some dependence on the quality of the basis set adopted. Recent work by some of us48 has shown that the ω01 frequency for SIL decreases by 47 cm-1 when passing from B3-LYP/DZP to B3-LYP/aug-ccpVDZ levels. If this correction is applied to the OH ω01 B3 value for EDIxS2, a value as low as 3694 cm-1 is arrived at, definitely underestimated when compared to the experimental data for amorphous silica.47 In fact, the depencence on basis set quality of energy derived quantities is usually less important
for periodic structures than for small molecular clusters.6 Our computer technology does not allow us at present to obtain an ab initio periodic estimate of ω01 for EDIxS2 at a B3-LYP/ aug-cc-pVDZ level, but it would be interesting to check in the future the present prediction. 5. Conclusions The results of this study show that the simple slab model here proposed has a number of attractive features which make it suitable for the study of some properties of silica surfaces. It is characterized by an appreciable stability, exhibits a topography which can mimic some aspects of the surface of amorphous silica (especially as concerns the distribution of isolated hydroxyls and their vibrational frequencies), and is computationally inexpensive. It can prove useful in the calibration of cluster models. In particular, while very accurate estimates of vibrational frequencies can be obtained with both models, the electrostatic field experienced at the surface of the periodic system is more shortranged and at the same time more structured than at the corresponding sites of the molecular models; if this aspect is not properly taken into account, the latter can provide an inaccurate description of the interaction with polar or easily polarizable molecules, as proved in section 4.C. In this area of research, the periodic model lends itself to being used in combination with appropriate embedding schemes for the study of local chemical processes.49,50 The model here considered is just one in a family of periodic structures which can be constructed from the same Y building block containing five SiO2 units, much in the spirit of the work undertaken by Ozin et al.51 By properly assembling these Y blocks and cutting the resulting structures along different planes, a variety of closely related surface topographies can be generated, for the simulation of different aspects of the surface reactivity of silica. Acknowledgment. The present work is part of a project coordinated by A. Zecchina and cofinanced by the Italian MURST (Cofin98, Area 03). The authors thank Dr. J. Gale for providing the most recent version of the GULP code and Prof. R. Dovesi for useful discussions. Financial support from CINECA is acknowledged. References and Notes (1) Catlow, C. R.; Ackermann, L.; Bell, R. G.; Cora`, F.; Gay, D. H.; Nygren, M. A.; Pereira, J. C.; Sastre, G.; Slater, B.; Sinclair, P. E. Faraday Discuss. 1997, 106, 1. (2) Sauer, J.; Ugliengo, P.; Garrone, E.; Saunders: V. R. Chem. ReV. 1994, 94, 2095. (3) Sauer, J.; Hill, J.-R. J. Phys. Chem. 1994, 98 1238. (4) Vigne´-Maeder, F.; Sautet, P. J. Phys. Chem. B 1997, 101, 8197. (5) Sindorf, D. W.; Maciel, D. W. J. Am. Chem. Soc. 1983, 105, 1487. (6) Pisani, C.; Dovesi, R.; Roetti, C. Hartree-Fock Ab Initio Treatment of Crystalline Solids. Lecture Notes in Chemistry; Springer: Berlin, 1988; Vol. 48. (7) Dovesi, R.; Saunders: V. R.; Roetti, C. CRYSTAL92 User Documentation; Universita` di Torino: Torino, 1996. (8) Orlando, R.; Pisani, C.; Ruiz, E.; Sautet, P. Surf. Sci. 1992, 275, 482. (9) Civalleri, B. Doctorate Thesis, University of Torino, Torino, 1998. (10) Gottardi, G.; Galli, E. Natural Zeolites; Springer: Berlin, 1985. (11) Meier, W. M.; Olson, D. H. Atlas of Zeolite Structure Types; Butterworth: London, 1987. (12) Civalleri, B.; Zicovich-Wilson, C. M.; Ugliengo, P.; Saunders, V. R.; Dovesi, R. Chem. Phys. Lett. 1998, 292, 394. (13) Dovesi, R.; Saunders, V. R.; Roetti, C.; Causa`, M.; Harrison, N. M.; Orlando, R.; Apra`, E. CRYSTAL95 User’s Manual; Universita` di Torino: Torino, 1996.
Simple Periodic Model of the Silica Surface (14) The B3-LYP calculations have been made possible by a new addition to the CRYSTAL95 code by V. R. Saunders, which has been made available to the present authors. (15) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (16) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (17) Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244. (18) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (19) Perdew, J. P. In Electronic Structure of Solids; Ziesche, P., Eschrig, H., Eds.; Akademie Verlag: Berlin, 1991, pp 11-20. (20) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (21) Towler, M. D.; Zupan, A.; Causa`, M. Comput. Phys. Commun. 1996, 98, 181. (22) GULP (General Utility Lattice Progam) written and developed by Gale, J. D., Royal Institution/Imperial College, UK, 1992-1994. (23) Gale, J. D. J. Chem. Soc., Faraday Trans. 1997, 93, 629. (24) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes; University Press: Cambridge, 1989. (25) Sierka, M.; Sauer, J. Faraday Discuss. 1997, 106, 41. (26) Schro¨der, K.-P.; Sauer, J. J. Phys. Chem. 1996, 100, 11043. (27) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Internal Series of Monographs on Chemistry 22; Clarendon Press: Oxford, 1990. (28) Aray, Y.; Bader, R. F. W. Surf. Sci. 1996, 351, 233. (29) Smith, J. V. Z. Kristallogr. 1983, 165, 191. (30) Taylor, W. H.; Meek, C. A.; Jackson, W. W. Z. Kristallogr. 1933, 86, 373. (31) Alberti, A.; Gottardi, G. Neues Jahrb. Mineral. Monatsh. 1975, 396. (32) Zhuravlev, L. T. Langmuir 1997, 3, 316. (33) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeeseman, J. R.; Keith, T. A.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; HeadGordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, revision C.3; Gaussian, Inc.: Pittsburgh, PA, 1995.
J. Phys. Chem. B, Vol. 103, No. 12, 1999 2171 (34) Gatti, C. TOPOND96 Users’s Manual; CNR-CSRSRC: Milano, 1996. (35) Gatti, C.; Saunders, V. R.; Roetti, C. J. Chem. Phys. 1994, 101, 10686. (36) Ugliengo, P.; Saunders, V. R.; Garrone, E. J. Phys. Chem. 1990, 94, 2260. (37) Senchenya, I. N.; Garrone, E.; Ugliengo, P. J. Mol. Struct. (THEOCHEM) 1997, 368, 93. (38) Garrone, E.; Kazansky, V. B.; Kustov, L. M.; Sauer, J.; Senchenya, I. N.; Ugliengo, P. J. Phys. Chem. 1992, 96, 1040 (39) Lindberg, B. J. Chem. Phys. 1988, 88, 3805. (40) Ugliengo, P. ANHARM - A program to solve Monodimensional nuclear Schroedinger equation; 1989, unpublished. (41) The CRYSTAL95 inputs corresponding to the geometries here reported can be sent on request (e-mail:
[email protected]). (42) Engelhardt, G.; Radeglia, R. Chem. Phys. Lett. 1984, 108, 271. (43) Moravetski, V.; Hill., J.-R.; Eichler, U.; Cheetham, A. K.; Sauer, J. J. Am. Chem. Soc. 1996, 118, 13015. (44) Pisani, C.; Causa`, M.; Dovesi, R.; Roetti, C. Prog. Surf. Sci. 1987, 25, 119. (45) Fubini, B.; Bolis, V.; Cavenago, A.; Garrone, E.; Ugliengo, P. Langmuir 1993, 9, 2712. (46) Civalleri, B.; Garrone, E.; Pisani, C.; Ugliengo, P. Presented at the 7th International Conference on Theoretical Aspects of Heterogeneous Catalysis, Cambridge, UK, 1998. (47) Kustov, L. M.; Borokov, V. Yu; Kazansky, V. B. Russ. J. Chem. Phys. (Khim. fiz.) 1982, 12, 1632. (48) Civalleri, B.; Garrone, E.; Ugliengo, P. J. Phys. Chem. B 1998, 102, 2373. (49) Pisani, C.; Cora`, F.; Nada, R.; Orlando, R. Comput. Phys. Commun. 1994, 82, 139. Pisani, C.; Birkenheuer, U. Comput. Phys. Commun. 1996, 96, 152. (50) Pisani, C.; Birkenheuer, U.; Cora`, F.; Nada, R.; Casassa, S. EMBED96 User’s manual; Universita` di Torino: Torino, 1996. (51) Oliver, S.; Kupermann, A.; Ozin, G. A. Angew. Chem., Int. Ed. Engl. 1998, 37, 46.
