Dehydroxylation and Silanization of the Surfaces of β-Cristobalite

J. Phys. Chem. B 2001, 105, 8007-8013

8007

Dehydroxylation and Silanization of the Surfaces of β-Cristobalite Silica: An ab Initio Simulation S. Iarlori,† D. Ceresoli,‡,§ M. Bernasconi,‡,* D. Donadio,‡ and M. Parrinello| Pirelli CaVi e Sistemi S.p.a., Viale Sarca 222, I-20126 Milano, Italy, Istituto Nazionale per la Fisica della Materia and Dipartimento di Scienza dei Materiali, UniVersita´ di Milano-Bicocca, Via Cozzi 53, I-20126 Milano, Italy, and Max-Planck-Institut fu¨ r Festko¨ rperforschung, Heisenbergstr. 1, D-70569 Stuttgart, Germany ReceiVed: March 1, 2001

Dehydroxylation and silanization processes on the silica surface are studied by ab initio molecular dynamics. The (100) and (111) surfaces of β-cristobalite are used as two possible models of the hydroxylated amorphous surface. The activation energy and latent heat for the dehydroxylation reactions of the (100) surface computed by constrained ab initio molecular dynamics are in reasonable agreement with experimental data on the amorphous surface. Adhesion reactions of silanes are simulated aiming at elucidating the binding mechanism of organosilanes used for instance as silica-polymer coupling agents. The simulation have provided insights on the occurrence of multiple silica-silane bonds and on the role of hydrolization of silane by physisorbed water in the adhesion on the wet surface.

I. Introduction Many technological applications of amorphous silica rely on its specific surface properties.1,2 A freshly formed amorphous silica surface will include a distribution of reactive sites which are known to react rapidly with atmosferic moisture leading to the formation of surface silanol groups (Si-OH). On the other hand, sol-gel silicas, grown in water solution, display a variety of surface silanols depending on sinthesis conditions. The concentration, distribution, and nature of silanols largely determine the technologically relevant properties of hydroxylated silica surface. To quote just one example, an important aspect in the tire manufacturing technology is represented by the adhesion properties of polymers on the surface of hydroxylated silica, added to reinforce the mechanical properties of the polymeric blend. The silica-polymer adhesion properties can be improved with the adsorption of organosilanes coupling agents which simultaneously bind to the silica and the polymer. At present, the TESPT ((EtO3-Si-(CH2)3-S2)2) molecule is the best performing silica-polymer coupling agent in tires technology.3,4 A detailed microscopic knowledge of the adhesion mechanism of the organosilanes is obviously of paramount importance for the design of new, better performing coupling agents and for the optimization of the technological process. Here we report on the first attempt to model the silicaorganosilane adhesion via ab initio molecular dynamics (MD) simulations.5,6 The first step in the elucidation of the reaction mechanism necessarily requires a full characterization of the surface silanols which are the expected reactive sites for the organosilane adhesion. Although a detailed microscopic description of the * To whom correspondence should be addressed. E-mail: [email protected]. † Pirelli Cavi e Sistemi S.p.a. ‡ Universita ´ di Milano-Bicocca. § Present address: International School for Advanced Studies (SISSA), Via Beirut 4, I-34014 Trieste, Italy. | Max-Planck-Institut fu ¨ r Festko¨rperforschung.

hydroxylated silica surface is still lacking, various experimental techniques including NMR,7-9 infrared7,10-15 and Raman16,17 spectroscopy, as well as chemical probes have provided many data on the properties of surface silanols which are often rationalized by modeling the surface as an alternation of patches of the hydroxylated (100) and (111) surfaces of β-cristobalite.1,7 This is the crystalline phase of silica with density and refractive index closest to those of amorphous silica. Moreover, the two main faces of β-cristobalite can sustain the two types of silanol groups identified experimentally on the amorphous silica surface, namely the single silanols (a single hydroxyl attached to a surface Si) typical of the (111) surface and the “geminal” silanols (two hydroxyls attached to the same surface Si), which are typical of the (100) surface. To study the silanization processes, we have adopted this β-cristobalite model as well and first studied the (100) surface which together with our previous calculations18 on the (111) surface complete the description of the clean surfaces needed to address their chemical activity. As expected on simple geometrical arguments, we have found that the vicinal, geminal silanols on the (100) surface are H-bonded as opposed to the single, isolated silanols on the ideal (111) surface18 which are not. The calculated vibrational frequencies of H-bonded and isolated silanols are in reasonable agreement with experimental data10,19 on amorphous surfaces. Then, as a simpler and test case of surface chemical reaction, we have studied the dehydroxylation process on the (100) surface of β-cristobalite. The dehydroxylation reaction leading to the formation of five-membered rings has been simulated by forcing the condensation of vicinal H-bonded silanols via constrained ab initio MD.20,21 The calculated dehydroxylation energy and energy barrier for the reaction are in good agreement with available experimental data on amorphous silicas surface1,22,23 and give us confidence on the reliability of our framework to address the central problem of our modeling efforts: the silanization process. We have studied the adhesion reaction on both (100) and (111) silica surfaces of tetraethoxysilane (Si(EtO)4, TEOS) and

10.1021/jp010800b CCC: $20.00 © 2001 American Chemical Society Published on Web 07/26/2001

8008 J. Phys. Chem. B, Vol. 105, No. 33, 2001 dimethyldiethoxysilane (Si(CH3)2(EtO)2). These molecules have the same functional groups for adhesion to silica (the ethoxy groups, EtO) of TESPT. The silanization mechanism of the surface of precipitated silica with TESPT3,4,24,25 and other organosilanes1,26 has been studied experimentally mainly by NMR. These measurements have shown that (i) TESPT reacts with the surface hydroxyls releasing ethanol, (ii) multiple bonding of TESPT on the surface is possible, and (iii) by increasing the quantity of water physiadsorbed on the surface, and by working in an acid or basic environment, the efficiency of the TESPT reaction improves.3,4,24 The experimental studies, however, have not provided a complete microscopic description of the process, which is a necessary step before alternatives to TESPT can be designed. Our aims are, in particular, a better understanding of two aspects of the adhesion process. First, from the NMR study, it has not been possible to tell whether the silane reacts with the surface by forming more bonds or it reacts primarily with a single surface hydroxyl and policondensates. Second, it is unclear whether the ethoxy groups of the silane can directly react with the surface hydroxyl or the silane must get first hydrolyzed by the physiadsorbed water in order to make the silane-silica adhesion possible.24 Anticipating our results, we find that the geometry of silanols on the (100) surface provides a simple mechanism for the formation of two silicasilane bonds. Moreover, the calculation of activation energies for the adhesion reactions suggests that organosilanes must be first hydrolyzed by the physiadsorbed water before reacting with surface hydroxyls. This intermediate reaction seems necessary in order to reproduce the adhesion rate measured experimentally.3,4,24 The paper is organized as follows: after a brief description of our theoretical framework and its benchmark on the properties of bulk crystalline silica (section II), our results on the structural and vibrational properties of hydroxylated (100) surface of β-cristobalite are reported in section III. The simulation of the dehydroxylation and silanization reactions are reported in sections IV and V, respectively.

II. Computational Method We performed ab initio MD simulations using the CarParrinello scheme.5,6 Our approach is based on density functional theory (DFT) within the local density approximation of the exchange-correlation functional augmented by the gradient corrections proposed by Becke27 and Lee, Yang, and Parr28 (BLYP). Gradient corrections to the local density approximation (LDA) are necessary to reproduce the correct energy hierarchy of the different crystalline polymorphs of silica29 as shown by calculations using others exchange and correlation energy functionals.29 We have chosen the BLYP functional because it provides a good description of the hydrogen bond.30,31 Only valence electrons are treated explicitly, and electron-ion interactions are described by norm-conserving pseudopotentials generated according to the Troullier and Martins procedure.32 Kohn-Sham orbitals are expanded in plane waves up to a kinetic energy cutoff of 70 Ry. This energy cutoff is large enough to correctly reproduce structural properties at equilibrium. In the calculation of the phonon frequencies and IR spectra, a larger energy cutoff of 80 Ry turned out to be necessary.33 Phonon frequencies ωn and their corresponding normalized eigenmodes ξn are derived by diagonalization of the dynamical matrix calculated by taking finite differences of the atomic forces. The IR spectrum is given by the imaginary part

Iarlori et al. of the transverse dielectric response function in the long wavelength limit as

Im ⊥(ω) )

FnR )

2π2

|Fn|2

∑ 3Ω n

∑I ∑β

δ(ω - ωn)

ω I ZRβ

I ξnβ

xmI

where I runs over the atoms, R and β run over Cartesian axis, n runs over the normal modes, and mI are atomic masses. The transverse effective charge tensors ZRβ are computed from finite difference of the total polarization P with respect to atomic displacement uIβ as I ZRβ )Ω

∂PR ∂uIβ

where Ω is the cell volume. The change in polarization is computed within the Berry phase approach of refs 34 and 35. We have checked the reliability of the BLYP functional by computing the structural parameters of R-quartz and β-cristobalite in the bulk. The unit cell of R-quartz has nine atoms in the unit cell with space group P312.36 We optimized the internal parameters at the experimental density (a ) 4.91 Å and c ) 5.401 Å36). We used only one special point37 in the Brillouin zone (BZ) integration which is enough to achieve convergence in the structural properties (within 1.5%) as checked in ref 38. The computed internal atomic coordinates in the notation of ref 36 are u ) 0.4647 for Si and x ) 0.4103, y ) 0.2768, and z ) 0.1103 for O atoms. Experimental values are39 u ) 0.4697, x ) 0.4135, y ) 0.2669, and z ) 0.1191. The crystalline β-cristobalite phase can be obtained by heating R-cristobalite above 270 °C at normal pressure.40 The space group Fd3m assigned experimentally to β-cristobalite has been interpreted in the past as an average structure composed of small domains with six possible orientations of the structure with I4h2d symmetry.40 More recent experimental data have shown that the coexistence of different domains of lower symmetry is in fact dynamical and not static as due to domains with fixed orientation.41 Although thermodynamically stable only above 270 °C, the low symmetry β-cristobalite structure is locally stable at low temperatures (0-300 K) in the small simulation cells with periodic boundary conditions used here. We first modeled β-cristobalite in a cubic supercell containing 8 formula units with the I4h2d symmetry. We optimized the internal coordinates at the experimental density of 2.18 g/cm3. Only the Γ point has been used in the integration of the BZ. The internal structural parameters of β-cristobalite are assigned by the SiO ˆ Si angle, two OSˆ iO angles, and the Si-O bond length. The calculated values are 142.8° (SiO ˆ Si), 107.4° and 113.6° (OSˆ iO), and 1.628 Å (Si-O) in good agreement with the experimental values 146.7°, 107.8° and 112.8°, and 1.611 Å40 and with previous LDA calculations.42 The surfaces are modeled by a slab geometry with 3D periodic boundary conditions. The initial geometry was obtained by cutting the bulk and saturating the surface oxygen atoms with hydrogen. In the surface plane, the supercells contain two surface unit cells. The slabs were 11-15 atomic layers thick and separated by 6 Å of vacuum. Geometry optimization were performed using both a simulated annealing procedure and a Broyden-Fletcher-Goldfarb-Shanno algorithm.43 In the Car-

Dehydroxylation of β-Cristobalite Silica

J. Phys. Chem. B, Vol. 105, No. 33, 2001 8009

Parrinello MD simulations, we used a time step of 0.12 fs and a fictitious electronic mass of 800 au. In the MD and geometry optimizations, the BZ integration has been restricted to the supercell Γ point. Finally, the simulation of the dehydroxylation reaction via condensation of vicinal hydroxyls on the (100) surface has been performed by the method of constraints, originally developed for classical MD.20 Constrained ab initio MD has already been used to simulate other chemical reactions in solution21 and at surfaces.44 We performed several simulations by constraining the reaction coordinate at different values. The reaction coordinate is the distance between the Si and O atoms belonging to different vicinal silanols. The simulation allowed us to compute the energy profile as a function of the reaction coordinate which provides the energy barrier for the reaction and the dehydroxylation energy. Although in principle the constrained dynamics provides an estimate not only of energy but also of the free energy differences between initial, final, and transition states,20,21,44 the latter have not been calculated because of insufficient statistical accuracy of our short simulation runs. The same technique has been used to compute the activation energy for the silanol-organosilanes reaction reported in section V. III. Hydroxylated (100) Surface The optimized geometry of the (100) slab is reported in Figure 1. The surface unit cell contains two geminal silanols on each surface of the slab. The orientation of the surface hydroxyls has been obtained as follows: we first performed a simulation of a thinner slab containing only three silicon layers. The initial orientation of the hydroxyl groups added to saturate undercoordinated atoms has been chosen not too far from the equilibrium geometry found in previous Hartree-Fock (HF) calculations.45 We broke the symmetry between the two surfaces of the slab by choosing slightly different H arrangements. After simulated annealing and geometry optimization, a difference in the orientations of the hydroxyls on the two surfaces appear. The two surface configurations are hereafter referred to as surface A and surface B (see Figure 1). Then we increase the slab thickness to four silicon layers. Geometry optimization of the thicker slab did not produce any appreciable change in the surface structure. On both surfaces A and B, the geminal silanols form H-bonded chains along the [100] direction, with the hydrogen bonds being 1.68-1.71 Å long. However, although on surface B the two geminal silanols in the unit cell are equivalent, each of them accepting and donating one hydrogen bond, the two geminal silanols on surface B are inequivalent. In fact, in the latter surface, one geminal silanol donates two H-bonds, whereas the other accepts two H-bonds. To obtain the surface energies of the two (A and B, cf. Figure 1) configurations, we have computed the total energy of a relaxed slab with both surfaces in configuration B. It turns out that surface B is slightly favored with respect to surface A by 7.2 kJ/mol per geminal silanol. This energy difference is comparable to the thermal energy at the R-β transition temperature in cristobalite (270 °C). Therefore, the two silanols configurations are expected to be both present at the surface of crystalline β-cristobalite. The orientation of the hydroxyl can be assigned by the dihedral angle formed by the OH group with the plane defined by the SiO bond and the normal to the surface. This dihedral angle (here called β as in ref 45) has values of +144° and -44° for one geminal silanol and +143° and -50° for the other silanol on surface B. On surface A, the dihedral angles are 143° and 132° for one geminal and -59° and -61° for the

Figure 1. Optimized geometry of the hydroxylated (100) surface. The orientation of the hydroxyl groups are different on the two opposite surfaces of the slab. Top panel: side view. Central panel: top view of the top surface (surface A). Bottom panel: top view of the bottom surface (surface B). Only surface atoms are shown in the top views. The surface supercell (50.8 Å2) periodically repeated in the plane is shown in the top views. The supercell contains two geminal silanols on each surface and a total number of 36 atoms with in-plane periodic boundary conditions. Dotted lines indicate H bonds. Oxygen, silicon, and hydrogen atoms are represented by black, gray, and white spheres, respectively.

other. The B configuration is close to the geometry found in previous HF calculation45 (β equal to +125° and -50°) where only one geminal silanol per surface unit cell has been considered. The rotational dynamics of the hydroxyls on surface A has been analyzed in a MD simulation at room temperature (in the larger slab).46 The dihedral angle β fluctuates by as much as 45°, but the H-bond length never exceeds 2 Å. Therefore, no H-bond breaking or interconversion between configurations A and B has been seen in our simulation run 0.5 ps long. This result implies the presence of an energy barrier for the breaking of a H-bond higher than the thermal energy. Previous HF calculation45 gives an energy barrier of 33 kJ/mol for the simultaneous breaking of two H-bonds induced by the rotation of the hydroxyls in the B configuration. We did not attempt to estimate this energy barrier better.

8010 J. Phys. Chem. B, Vol. 105, No. 33, 2001

Iarlori et al.

We calculated the IR active modes in the OH stretching regions as described in section II. In our slab, we have eight hydroxyls (four on each surface of the slab) and, therefore, eight OH stretching modes. The four modes corresponding to the stretching of the hydroxyls which are H-bond donors are in the range of 3190-3373 cm-1, whereas the four stretching modes of the hydroxls which do not donate a H bond are in the range of 3608-3614 cm-1. These frequencies have been computed at the lower 70 Ry cutoff which overestimates the OH stretching modes of the H2O molecule by ∼40 cm-1 with respect to the frequencies at full convergence with the cutoff energy (80 Ry) and by ∼70 cm-1 with respect to experiments.33 By assuming a constant error in the stretching frequencies of silanols, our results are in satisfactory agreement with the experimental IR spectra1,10,19 of amorphous silica surface which show a very broad peak around 3530 cm-1, attributed to H-bonded silanols, and several sharp peaks in the range of 3740-3750 cm-1, attributed to non-H-bonded silanols. Calculation of the IR spectrum of the (111) surface18 gives the OH stretching frequency of isolated silanols at 3715 cm-1 (at 80 Ry cutoff) still in good agreement with the aforementioned experimental data. IV. Dehydroxylation Reactions As emerged from IR7,10-15,19, Raman,16,17 and NMR47,8 measurements, the amorphous silica surface undergoes dehydroxylation reactions via condensation of vicinal silanols upon heating under vacuum. H-bonded silanols condense first at temperatures above 190 °C; the stronger the hydrogen bonding the easier the hydroxyls removal. Only isolated (non-H-bonded) silanols remain on a silica surface which has been evacuated above 500 °C (see refs 1 and 7 for a review of experimental data). Further dehydroxylation above 500 °C eliminates a large fraction of isolated hydroxyls and produces highly reactive sites able to chemisorb ammonia, methanol, and other chemical probes.48-51,11-15 In the current literature,48-51,14,15 these latter sites are customarily, even if rather tentatively, assigned to strongly metastable two-membered silicon rings. Our previous ab initio simulations of the dehydroxylation reaction of the (111) surface of β-cristobalite provided18 compelling evidence of the existence of two-membered silicon rings on the dehydroxylated surface. Here, we study the dehydroxylation of the (100) surface leading to the larger five-membered silicon ring which are expected to be formed in the early stages of dehydroxylation at lower temperatures. We computed the energy barrier for the condensation of two vicinal silanols on the (100) surface by performing constrained MD.20,21 The reaction coordinate is assumed to be the distance between the silicon atom of a geminal silanol and the nearest oxygen atom of the vicinal silanol (see Figure 2a). We performed constrained MD run at 300 K and then geometry optimization at several values of the reaction coordinate. The minimum energy as a function of the reaction coordinate (Si-O distance) is shown in Figure 3. The geometry of the transition state is shown in Figure 2b; it corresponds to the configuration of maximum energy in Figure 3. The dehydroxylation reaction leading to the final configuration shown in Figure 2c is endothermic, i.e., energy must be supplied to remove a water molecule via condensation of the two silanols. The calculated dehydroxylation energy is 42 kJ/mol, much lower than the experimental average energy of dehydroxylation (∼81 kJ/mol for temperatures in the range of 350-650 °C16). Because the condensation of vicinal silanols on the (100) surface generates a low-strained silicon ring, it is conceivable that the theoretical dehydroxylation energy is lower than the experi-

Figure 2. Initial (top panel), transition (central panel), and final (bottom panel) states of the dehydroxylation reaction of the (100) surface via condensation of vicinal silanols which leads to the formation of a fivemembered silicon ring. The reaction coordinate is the Si-O distance shown in the top panel. The color code is the same as that in Figure 1.

Figure 3. Total energy of the slab in Figure 2 as a function of the reaction coordinate in the condensation reaction of the two vicinal silanols. The reaction coordinate is the Si-O distance shown in Figure 2. The energy of the final state is obtained by summing the energy of the final configuration in Figure 2 to the energy of the isolated released water molecule. The highest energy corresponds to the transition state. The final configuration has been optimized without any constrained on the reaction coordinate.

mental one which results from an average over reactions forming four-, three-,17 and two-membered rings.48-51 The dehydroxyl-

Dehydroxylation of β-Cristobalite Silica

Figure 4. Optimized structure of the dimethyldiethoxysilane reacted with a hydroxyl of the (100) surface of β-cristobalite (cf. Figure 1). Top panel: side view. Bottom panel: top view. Only the first few surface layers are shown in the bottom panel. The simulation cell contains 134 atoms with in-plane periodic boundary conditions. Oxygen, silicon, carbon, and hydrogen atoms are represented by black, dark gray, light gray, and white spheres, respectively.

ation energy found here is much lower than the values obtained in previous HF45 and semiempirical calculations52 where only limited relaxations on small systems were performed. In contrast, our calculated dehydroxylation energy is very consistent with the value of 50 kJ/mol obtained in a recent DFT calculation that uses a Gaussian basis set and a larger surface supercell.53 From the difference between the energy of the initial and the transition states in Figure 3, we also obtain an activation energy for the reaction of 313 kJ/mol. Experimentally the activation energy of the dehydroxylation reaction increases by decreasing the silanol coverage.1,22,23 At high silanol coverage, the dehydroxylation proceeds via condensation of H-bonded vicinal silanols with an activation energy in the range of 80-200 kJ/ mol.1,22,23 When vicinal silanols are consumed, further condensation requires the migration of hydroxyls which implies a higher activation energy in the range of 200-300 kJ/mol.1,22,23 Our computed activation energy is larger than the experimental values. However, because of the small size of the supercell and the use of periodic boundary conditions, we are forced to simulate not an isolated geminal pair but an infinite row. This is likely to provide an overestimation of the activation energy.54 V. Silanization of the Silica Surface As a first example of the silanization of the silica surface, we have studied the adhesion of the silane Si(CH3)2(EtO)2 on the hydroxylated (100) surface of β-cristobalite. From the NMR studies, it has not been possible to tell whether the silane reacts with the surface by forming more bonds or it reacts primarily with a single surface hydroxyl and policondensates. It turns out that the geometry of the geminal silanols on the (100) surface

J. Phys. Chem. B, Vol. 105, No. 33, 2001 8011

Figure 5. Same as Figure 4 for the reaction of the silane with two hydroxyls of the (100) surface of β-cristobalite. The two reacted hydroxyls belong to two neighboring, H-bonded, geminal silanols. The simulation cell contains 127 atoms. The color code is the same as that in Figure 4.

provides a simple mechanism for the formation of two silicasilane bonds. In fact, the silane can react either with one or simultaneously with two surface hydroxyls, releasing one or two ethanol molecules, respectively. The two hydroxyls involved in the adhesion belong to two neighboring H-bonded, geminal silanols. The energy released by the reaction is roughly 19 kJ/ mol larger when the silane forms two bonds with the surface, with respect to the configuration with a single bond. Therefore, there are surface sites where multiple bonds are energetically favored. The configurations with one and two bonds are shown in Figures 4 and 5, respectively. Because the fraction of geminal silanols observed by NMR is 1/21 for Fisher or Baker silica gel,47 the bonding configuration we have simulated is likely to occur but cannot probably account for all of the multiple siloxane bonds seen experimentally by NMR in silica-silane adhesion.3,4,25,26 On the other hand, to address the role of hydrolysis in the silanization process, we have simulated the adhesion of TEOS with a single silanol on the (111) surface of β-cristobalite, given the fact that the single hydroxyl is the most abundant type of silanols present on the amorphous silica surface. In fact, at the hydroxylated (111) surface the concentration of single silanols (4.5 OH/nm2) is very close to the experimental value on the hydroxylated amorphous surface (4.9 OH/nm2, including geminal silanols47). The (111) surface is modeled by a silicon bilayer as shown in Figure 6. In the simulation, TEOS has been forced to react directly with the surface hydroxyl, releasing an ethanol molecule. We performed constrained ab initio MD by using as reaction coordinate the distance between the silicon atom of the silane and the oxygen of the surface hydroxyl. In this way, we imposed that the bond broken in the adhesion is the Si-O bond of the silane in analogy with the description of the hydrolyzation reaction of TEOS in water.55 The energy as a function of the chosen reaction coordinate in the direct adhesion

8012 J. Phys. Chem. B, Vol. 105, No. 33, 2001

Iarlori et al. assessment of the role of hydrolyzation in the adhesion process could come from a time-resolved simultaneous measurements of the ethanol release and of the intensity of the Q3 peak in 29Si NMR. The time-resolved NMR data would measure the consumption of surface hydroxyls and would therefore provide the rate of the real adhesion reaction. Our calculations provide further motivations for such an experiment which has not been performed so far. VI. Conclusions

Figure 6. Transition state configuration for the simulated reaction of TEOS with a hydroxyl on the (111) surface of β-cristobalite. A very thin slab is used to model the surface. The silanols on the bottom surface have been held fixed during the simulation. The transition state corresponds to the configuration of maximum energy in Figure 7.

Figure 7. Total energy for the reaction of TEOS with a single, isolated silanol on the (111) surface of β-cristobalite as a function of the reaction coordinate. The reaction coordinate is the distance between the silicon atom of the silane and the oxygen of the surface silanol (cf. Figure 6).

of TEOS with the surface hydroxyl is shown in Figure 7. The transition state is shown in Figure 6. The activation energy is as large as 154 kJ/mol. Experimentally, the reaction rate of TESPT has been measured24 by monitoring the ethanol release. However, this analysis cannot discriminate between the direct adhesion of the organosilane, as in the reaction sketched in Figure 6, from the hydrolyzation of silane and the successive condensation of the surface and silane hydroxyls, because both reactions lead to the same ethanol release. The activation energy for the reaction estimated experimentally (from the ethanol release) is ∼48 kJ/mol,24 too low with respect to our value for the direct TEOS-hydroxyl reaction (154 kJ/mol). Solvation effects, here neglected, might decrease the activation energy for the direct TEOS-hydroxyl reaction. However, the need of physisorbed water to decrease the activation energy naturally suggests that water would eventually trigger the other reaction channel, i.e., the hydrolyzation of the silane itself. Therefore, the latter reaction is probably what is measured experimentally from the ethanol release.24 The occurrence of a previous hydrolyzation is also in agreement with the observed strong dependence of the ethanol release on the pH and on the moisture content.24 However, the reaction rate measured experimentally for TESPT (k ) 5 × 10-5 s-1, at 50 °C, pH ) 6.3, and moisture content 11.1 wt %)24 is lower than the rate of hydrolyzation of TEOS in water at the same pH conditions (k ) 1.4 × 10-2 s-1, at 20 °C).56 This low rate might be due to an insufficient hydration of the surface at the experimental conditions at which TESPT-silica reaction rates have been measured.24 A definitive

We have studied dehydroxylation and silanization reactions on the hydroxylated (100) and (111) surfaces of β-cristobalite as two possible models of the amorphous silica surface. In agreement with previous theoretical works,45,52,53 we find that the geminal silanols on the (100) surface are H-bonded to vicinal silanols. Condensation reactions of the H-bonded silanols lead to the formation of five-membered silicon rings. The activation energy (313 kJ/mol) and latent heat (42 kJ/mol) of this process, as calculated from constrained ab initio MD, are in reasonable agreement with the experimental data on the amorphous surface. As a first attempt to model the adhesion reaction of organosilane used for instance as silica-polymer coupling agents, we have simulated the reaction of TEOS and dimethyldiethoxysilane with the surface silanols on the (111) and (100) faces, respectively. The H-bonded geminal silanol on the (100) face is a very likely site for multiple silane-silica bonds whose existence has been inferred from NMR measurements. The silane can react with the H-bonded hydroxyls releasing one or two ethanol molecules. The calculated energy released in the formation of two silicasilane bonds is 19 kJ/mol larger than the formation energy of a single silane-silica bond. However, because the most abundant silanol type on the amorphous surface is the single silanol, we have studied the reaction of TEOS with the single hydroxyls on the (111) face of β-cristobalite as well. The calculated activation energy for this latter process is 154 kJ/mol. Experimentally it is deduced an activation energy of 48 kJ/mol24 from the measurement of the rate of ethanol release in the adsorption of silanes on wet silica surfaces. This very low value suggests that in the production of ethanol the preferred channel is not the direct silane-surface adhesion, but instead, the silane is first hydrolyzed by the adsorbed water and only later reacts with the surface hydroxyls. Acknowledgment. We thank L. Garro and F. Negroni for providing us with information on the silica-silane reactions. Discussions with P. Sozzani, R. Simonutti, and E. Tosatti are gratefully acknowledged. This work is partially supported by the INFM Parallel Computing Initiative. D.D. is grateful for a graduate scholarship granted by Pirelli Cavi e Sistemi S.p.a. References and Notes (1) The surface properties of silica; Legrand, A. P., Ed.; Wiley: New York, 1998. (2) Iler, R. K. The Chemistry of Silica; Wiley: New York, 1979. (3) Go¨rl, U.; Hunsche, A. H. Proceeding of the meeting of the Rubber DiVision of the American Chemical Society; American Chemical Society: Washington, D.C., October 1997. (4) Hunsche, A.; et al. Kautsch. Gummi Kunstst. 1997, 12, 881. (5) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (6) CPMD, written by Hutter J. et al., MPI fu¨r Festko¨rperforschung and IBM Research Laboratory 1990-2000. (7) Chuang, I.-S.; Maciel, G. E. J. Phys. Chem. B 1997, 101, 3052. (8) Bronnimann, C. E.; Ziegler, R. C.; Maciel, G. E. J. Am. Chem. Soc. 1988, 110, 2023. (9) Kinney, D. R.; Chuang, I.-S.; Maciel, G. E. J. Am. Chem. Soc. 1993, 115, 6786. (10) Hoffman, P.; Kno¨zinger, E. Surf. Sci. 1987, 188, 181.

Dehydroxylation of β-Cristobalite Silica (11) Morrow, B. A.; Cody, I. A. J. Phys. Chem. 1976, 80, 1995. (12) Morrow, B. A.; Cody, I. A. J. Phys. Chem. 1976, 80, 1998. (13) Morrow, B. A.; Cody, I. A. J. Phys. Chem. 1976, 80, 2761. (14) Bunker, B. C.; Haaland, D. M.; Ward, K. J.; Michalske, T. A.; Smith, W. L.; Binkley, J. S.; Melius, C. F.; Balfe, C. A. Surf. Sci. 1989, 210, 406. (15) Bunker, B. C.; Haaland, D. M.; Michalske, T. A.; Smith, W. L. Surf. Sci. 1989, 222, 406. (16) Brinker, C. J.; Tallant, D. R.; Roth, E. P.; Ashley, C. S. J. NonCryst. Solids 1982, 82, 117. (17) Humbert, B. J. Non-Cryst. Solids 1995, 191, 29. (18) Ceresoli, D.; Bernasconi, M.; Iarlori, S.; Parrinello, M.; Tosatti, E. Phys. ReV. Lett. 2000, 84, 3787. (19) Morrow, B. A.; McFarlan, A. J. J. Phys. Chem. 1992, 96, 1395. (20) Ciccotti, G.; Ferrario, M.; Haynes, J. T.; Kapral, R. Chem. Phys. 1989, 129, 241. (21) Sprik, M.; Ciccotti, G. J. Chem. Phys. 1998, 109, 7737. (22) Zhuravlev, L. T. Colloids Surf., A 1993, 74, 71. (23) Gillis-D’Hamers, I.; Vraken, K. C.; Vansant, E. F.; De Roy, G. J. Chem. Soc. Faraday Trans. 1992, 88, 2047. (24) Hunsche, A.; et al. Kautsch. Gummi Kunstst. 1998, 7-8, 525. (25) Simonutti, R.; Comotti, A.; Negroni, F.; Sozzani, P. Chem. Mater. 1999, 11, 822. (26) Sindorf, D. W.; Maciel, G. E. J. Am. Chem. Soc. 1983, 105, 3763. (27) Becke, A. D. Phys. ReV. A 1998, 38, 3098. (28) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785-789. (29) Hamann, D. R. Phys. ReV. Lett. 1996, 76, 660 and references therein. (30) Sprik, M.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1996, 105, 1142. (31) Rothlisberger, U.; Parrinello, M. J. Chem Phys. 1997, 106, 4658. (32) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993. (33) The calculated vibrational frequencies of the water molecule are 1560, 3612, and 3715 cm-1 (1546, 3581, and 3682 cm-1) at 80 Ry (70 Ry) cutoff, to be compared with the experimental values57 1595, 3657, and 3756 cm-1. (34) King-Smith, R. D.; Vanderbilt, D. Phys. ReV. B 1993, 47, 1651. (35) Resta, R. ReV. Mod. Phys. 1994, 66, 899. (36) Wyckoff, W. G. Crystal Structures, 4th ed.; Interscience: New York, 1974; p 318.

J. Phys. Chem. B, Vol. 105, No. 33, 2001 8013 (37) Monkhorst, H. J.; Pack, J. D. Phys. ReV. B 1976, 13, 5188. (38) Gonze, X.; Allan, D. C.; Teter, M. P. Phys. ReV. Lett. 1992, 68, 3603. (39) Levien, L.; Prewitt, C. T.; Weidner, D. Am. Mineral. 1980, 65, 920. (40) Wright, A. F.; Leadbetter, A. J. Philos. Mag. 1975, 31, 1391. (41) Spearing, D. R.; Farnan, I.; Stebbins, J. F. Phys. Chem. Minerals 1992, 19, 307. (42) Liu, F.; Garofalini, S. H.; King-Smith, R. D.; Vanderbilt, D. Phys. ReV. Lett. 1993, 70, 2750. (43) Press, W. H.; Teukolsky, A. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes; Cambridge Universiry Press: Cambridge, 1992; p 418. (44) Boero, M.; Parrinello, M.; Terakura, K. J. Am. Chem. Soc. 1998, 120, 2746. (45) Vigne´-Maeder, F.; Sautet, P. J. Phys. Chem. B 1997, 101, 8197. (46) In the dynamical simulation, we fixed the silanols of the other surface (including the silicon atoms). (47) Chuang, I-S.; Maciel, G. E. J. Am. Chem. Soc. 1996, 118, 401. (48) Michalske, T. A.; Bunker, B. C. J. Appl. Phys. 1984, 56, 2687. (49) Dubois, L. H.; Zegarski, B. R. J. Phys. Chem. 1993, 97, 1665. (50) Dubois, L. H.; Zegarski, B. R. J. Am. Chem. Soc. 1993, 115, 1190. (51) Grabbe, A.; Michalske, T. A.; Smith, W. L. J. Phys. Chem. 1995, 99, 4648. (52) Branda, M. M.; Castellani, N. J. Surf. Sci. 1997, 393, 171. (53) Mortensen, J. J.; Parrinello, M. J. Phys. Chem. B 2000, 104, 2901. (54) Even in our large slab, the constraints imposed on the bottom layer of the slab prevent a full relaxation of the strain energy of the five-membered ring. In fact, in a separate simulation of the condensation reaction where both surfaces of the slab were free to move, we observed a large rearrengement of the botton surface and a decrease of the activation energy down to 173 kJ/mol. (55) Brinker, C. J.; Scherer, G. W. Sol-Gel Science; Academic Press: Boston, 1996; p 116. (56) Aelion, R.; Loebel, A.; Eirich, F. J. Am. Chem. Soc. 1950, 72, 5705. (57) Kuchitsu, K.; Morino, Y. Bull. Chem. Soc. Jpn. 1965, 38, 805.