Density Functional Theory Investigation into Structure and Reactivity of

J. Phys. Chem. C 2008, 112, 12653–12662

12653

Density Functional Theory Investigation into Structure and Reactivity of Prenucleation Silica Species Charles L. Schaffer and Kendall T. Thomson* School of Chemical Engineering, Purdue UniVersity, West Lafayette, Indiana 47907 ReceiVed: October 4, 2006; ReVised Manuscript ReceiVed: April 16, 2007

We investigated the condensation and reverse hydrolysis reactions of several silica clusters ranging in size up to the octamer cage. Using density functional theory (DFT) at the B3LYP/6-31G(d,p) level of theory, we found that under neutral conditions the reactions proceeded through a single-step, SN2-type mechanism with formation of a pentacoordinated transition state. Under acidic conditions, the reactions of the smaller species followed a two-step mechanism with formation of a stable pentacoordinated intermediate, while the larger species proceeded through a single-step mechanism similar to the neutral route. In vacuo energy calculations showed a decreased activation barrier for formation of the larger ringed oligomers compared to the barriers of the smaller species. Calculation also showed an increasing kinetic and thermodynamic barrier for the silica clusters to break back into smaller less-condensed constituent species as they oligomerized to form larger species. To study the influence of solvation on these reactions, we used a hybrid implicit/explicit hydration model that explicitly accounted for water in the calculations. Results on the silica dimer cluster revealed a marked change in both the mechanism and energetics of the reactions. The results suggested the presence of excess water (1) elongated the hydroxyl bonds of the cluster, (2) elongated the siloxane bonds, (3) stabilized both the reactive intermediates and transition states, and (4) acted as a proton acceptor to facilitate the reactions. The calculated rate-limiting activation barrier for the condensation reaction to form the dimer species was in excellent agreement with previous experimental and theoretical results. Application of the hybrid solvation model to calculate energy of solvation for several silica oligomers showed a decrease in solvation energy as the silica clusters grew in size. I. Introduction Crystalline tectosilicates, such as zeolites, are among some of the most commercially important materials. Along with their traditional uses for catalysis, ion exchange, and separations, these materials continue to enjoy widespread popularity due to the potential to synthesize ‘designer’ materials for specific applications.1 To this end, much attention has been paid to elucidating the synthesis mechanism from prenucleation to crystallization. Despite this effort, there have only been a few mechanisms proposed for the formation of such complex solids.2–8 Even those lack the detail necessary to move this field from heuristic methods to a rational and well-defined synthesis route necessary for precise control of material properties. Of all of the stages in the synthesis mechanism, the one most elusive to researchers could be the events occurring within the first few minutes (or seconds) of particle formation. The early stages of synthesis are characterized by rapid hydrolysis of a soluble silica source and subsequent condensation of the resulting silica species. The exact nature of the condensation reactions is still not well understood, but it is believed that this process could have a significant influence on the structural properties of the resulting crystalline solid.9 Therefore, understanding these reactions at a fundamental level is crucial to controlling the synthesis of these materials. Among the tasks undertaken by researchers in this area have been to identify prenucleation silica species, provide reaction path analysis, propose mechanism, and calculate reaction energies (i.e., * To whom correspondence should be addressed. E-mail: thomsonk@ ecn.purdue.edu. Phone: (765) 496-6706. Fax: (765) 494-0805.

activation barriers) for a few of the lower weight oligomers. These studies have been performed using a combination of experimental and computational techniques. The starting point for zeolite synthesis is the sol-gel solution. The identification of silica species present in solution is not a trivial task. Due in part to the number of different silica oligomers present in solution, it is difficult for experimental techniques to identify individual cluster properties and almost impossible to track individual reactions. Instead, most experiments are only able to identify concentrations and bulk structural features (i.e., the number of Si-O bonds) of species in solution. One technique that is particularly well-suited for such investigations is 29Si nuclear magnetic resonance (NMR) spectroscopy. Knight et al. used this technique to study condensed silica species present in aqueous solutions. In a series of papers,10–14 they identified species containing Qn sites with n ) 1, 2, and 3 (In Qn notation, where n represents the number of bridging oxygens bonded to a silicon). A major result of their work was the identification of a Si8O20H8 species, the cubic octamer. Subsequent investigations have identified structures containing tetrahedrally coordinated silicon sites believed to be derivatives of the cubic octamer.15,16 This finding is significant because the Q4 site is the principle building unit of the final zeolite product. Recently, Pelster et al.9 were able to track the evolution of different silica species undergoing condensation reactions in aqueous solutions. Using electrospray mass spectrometry (ESI MS) in conjunction with three different reactor systems, they were able to follow the growth of silica oligomers in solution. They showed that using different ammonium salts as templating agents lead to different products, a result supported

10.1021/jp066534p CCC: $40.75  2008 American Chemical Society Published on Web 07/31/2008

12654 J. Phys. Chem. C, Vol. 112, No. 33, 2008

Schaffer and Thomson

SCHEME 1: SN2-Attack Mechanism of Silica Condensation and Reverse Hydrolysis Reaction under Neutral Conditions

SCHEME 2: “Flank-Side” Attack Mechanism of Silica Condensation and Reverse Hydrolysis Reaction under Neutral Conditions

SCHEME 3: Mechanism of Silica Condensation and Reverse Hydrolysis Reaction under Acidic Conditions

by other investigations of similar systems.4–6 They identified the cubic octamer, prismatic hexamer, and cyclic-trimer as the predominate products when using TMA+, TEA+, and TPA+, respectively. They also identified Q4 species previously char-

Figure 1. Schematic representation of reactions building up to larger oligomers from the monomer species.

acterized using 29Si NMR spectroscopy. In total, they were able to identify the possibility of 38 species present at some time in solution. It is generally argued that the hydrolysis reaction proceeds by direct nucleophilic attack of silicon by the oxygen of water.17–19 The reverse water producing condensation reaction proceeds via a reaction in which an oxygen of a hydroxyl group associates with either another silicon of the same cluster (ringforming reaction) or a silicon of a different cluster (chainlengthening reaction).20 A variety of mechanisms for the hydrolysis reaction have been discussed with two being generally accepted as most probable. The first involves a SN2-type backside attack in which the nucleophile, the oxygen of water in the case of hydrolysis, attacks from the side opposite the leaving group (see Scheme 1). Pereira et al.21 have studied the condensation reaction of two silicic acid monomers using a combined density functional/conductor-like screening model (COSMO22) to simulate hydration effects. The investigation supports a SN2-type mechanism on the basis of energetic and statistical analysis of the reaction intermediates. The second mechanism involves a “flank” or lateral attack by the nucleophile (see Scheme 2). Elanany et al.23 argued in favor of a lateral attack mechanism based on analysis of the hydrolysis reaction of tetramethoxysilane, Si(OCH3)4. The findings were based on

Structure and Reactivity of Prenucleation Silica Species

J. Phys. Chem. C, Vol. 112, No. 33, 2008 12655 Because these reactions occur in solvated environments, it is also of interest to explore the effects solvation has on these reactions. After a discussion of our in vacuo results, we present our solvation model and results obtained for the dimer molecule using this model. From here on, we can use the terms hydrolysis and condensation somewhat synonymously as we consider hydrolysis simply the reverse reaction of the water producing condensation reaction (see Schemes 1–3). Therefore, activation energies for each reaction can be obtained simply by judicious examination of the reaction energy diagram of one reaction or the other.

Figure 2. Hybrid implicit/explicit solvation model showing solute surrounded by a first shell of explicit solvent molecules and an outer layer of implicit solvent.

molecular dynamics (MD) simulations using a tight-binding quantum chemical routine “Colors”.24 The polymerization reaction of silica species has been widely studied using MD simulations.21,23,25–28 Recently, Gelb et al.26 performed simulations on systems containing 729 silicic acid monomers in aqueous solutions with varying densities. Using potentials developed by Fueston and Garofalini,28 the polymerization details, including structural features and activation energies to form the dimer molecule, were investigated and found to be in good agreement with experimental results20 (Eact ) ∼12-15 kcal/mol). However, MD simulations, such as the ones performed in this study, can only be used to analyze details occurring on relatively short time-scales and, as the author warns, provide little insight into the individual reactions of the silica species. Ab initio techniques have also been used to study the stability and reactions of silica species.29–35 Catlow et al.34 used density functional theory (DFT) with varying levels of theory to study the structural and energetic properties of various silica clusters containing up to eight silicons. They reported in vacuo condensation energies to form each oligomer directly from the monomer and analyzed the structural details of the most stable conformers. Pereira et al.21 conducted an ab initio investigation of the reaction mechanism and energetics of the formation of a dimer species from two monomers in a simulated hydrated environment. Their calculated activation energy to form the silica dimer was in good agreement with previous experimental results. Trinh et al.35 performed a similar study on the oligomerization reaction of several silica species under neutral and anionic conditions using the COSMO22 model to account for solvent effects. We investigate, via electronic density functional theory, the reactions of several silica clusters involved in prenucleation and elucidate the mechanism and energetics of their condensation and reverse hydrolysis. We analyzed seven silica clusters ranging in size up to the octamer cage. Figure 1 shows a schematic of how each oligomer could be built up from the monomer. The clusters analyzed in this work are a good, albeit not inclusive, representation of the types of clusters found in sol-gel solutions. We examined two linear species, the dimer and trimer, three single ringed species, the cyclic-trimer, branched cyclic-trimer, and cyclic-tetramer, and two multiringed species, the prismatic hexamer (prism) and cubic octamer. Note in Figure 1 we show the linear tetramer for coherence although the condensation reaction to form the linear tetramer was not studied in this work. Beginning with gas phase (in vacuum) calculations, we present a detailed analysis of the mechanism and energetics of the reactions under neutral and simulated acidic conditions.

II. Computational Details Structural optimizations and reaction energies for each of the clusters were obtained using density functional theory within the Gaussian 0336 software package. In vacuo calculations for the condensation and reverse hydrolysis reactions were performed using the B3LYP hybrid density functional37 with a double-ζ basis set, 6-31G(d,p)38 with added polarization functions. Zhao et al.39 tested the performance of 44 DFT methods with three basis sets (6-31G(d,p), aug-cc-pVTZ, and MG3S) for nonbonded interactions by developing four benchmark databases meant to test: hydrogen bonding, charge transfer, dipole interaction, and weak interactions. After averaging the results of each DFT method for all four databases, the best performing method was found to have an 11% error when compared to the benchmark values, while the B3LYP functional yielded an 18% error. Previous researchers have used the B3LYP functional to study similar silica systems.21,31,34,35 Nicholas et al. tested the sensitivity of basis sets by optimizing the geometry of various analogs of disiloxane and found best agreement with experiment using large basis sets with added polarization functions.40,41 We have tested the performance of geometry optimization on disiloxane (H3Si-O-SiH3) using both the 6-311G and 6-31G(d,p) basis sets with the B3LYP functional. The tests confirmed the inadequacies of the unpolarized basis sets to capture the structural features of the molecule. The calculated Si-O-Si bond angle, Si-O bond lengths, and Si-H bond lengths were found to be 179.9°, 1.67 Å, and 1.49 Å, respectively, for the 6-311G basis set and 152.7°, 1.64 Å, and 1.48 Å, respectively, for the 6-31G(d,p) basis set. Experimental values for the Si-O-Si bond angle, Si-O bond lengths, and Si-H bond lengths are 144.1-151.2°, 1.634 Å, and 1.486 Å, respectively.42 The mechanisms and energies of each reaction were studied using a two-part strategy. First, approximate reaction paths were obtained using a constrained optimization procedure. During the constrained optimizations, we successively decreased the distance between two atoms performing energy minimizations at each step, thus mapping out an energy profile along the bonding pathway. The geometries obtained from the constrained optimizations were then used as initial guesses to locate transition states and stable intermediates allowing all degrees of freedom to fully relax. Transition states and intermediates were obtained using the synchronous transit-guided quasi-Newton (STQN) method43,44 and Berny optimization,45 respectively. Frequency calculations were performed for each intermediate and transition state. For each optimized intermediate, we verified that all frequencies were positive and for transition states we verified the existence of only one imaginary frequency. All reported electronic energy activation barriers, denoted ∆Eact in the text, were corrected for zero-point vibrational (ZPV) energies. The frequency calculations also provide thermochemical analysis of the reactions at 1 atm pressure and a temperature



TABLE 1: In Vacuo Energies (kcal/mol) for the Condensation and Reverse Hydrolysis Reaction4

1 Corrected for zero-point vibrational energies. 2 We only show the schematic representation of the neutral hydrolysis and condensation reaction. See Scheme 3 for the mechanis of the acidic reaction. Lines represent Si-O-Si linkages and circles represent the silicon monomer, Si(OH)4. 3 The inclusion of ZPV energy corrections lowered the calculated energies resulting in a negative activation barrier. 4 Energies are calculated at the B3LYP/6-31G(d,p) level of theory.

of 298.15 K. We report the free energy of activation, ∆Gact, which measures the change in Gibbs free energy from each reactive intermediate to the subsequent transition state. The details of the procedures for calculating thermochemical data including the use of partition functions can be found in the book by McQuarrie and Simon.46 To incorporate solvation effects in to the calculations, we used a hybrid implicit/explicit solvation model to investigate the condensation and reverse hydrolysis reaction of the silica dimer. We incorporated hydration effects into our calculations by encasing the solute within a coordinating shell of water molecules, surrounding each silica fragment. The entire system was then fully optimized at the B3LYP/6-31G(d,p) level of theory. After optimization, additional water molecules were placed around the fragment and another complete optimization was performed. The process was continued until there were no structural changes observed within the solute. The use of explicit solvent molecules makes it difficult to ensure we are at a global

energy minimum. To overcome this, we have tested many starting geometries of the solvent and attempted perturbations on the optimized system. For the ringed species, we attempted several starting geometries with explicit solvent molecules located inside the ring. However, the lowest energy configuration of the explicit solvent was found to be outside the ring where hydrogen bonding was strongest. The use of an implicit solvent layer further corrects for areas of the solute that are not solvated with explicit solvent. Previous researches have used similar solvation techniques which combine discrete incorporation of solvent molecules with continuum descriptions of the solvent.47–51 Once the hydration sphere was fully optimized, the intermediate and transition state geometries of the reaction were obtained using the same constrained optimization technique outlined above but allowing the water molecules to re-equilibrate at each step in the optimization. Each of the optimized solute/explicit solvent intermediates and transition states were then subjected to a single-point energy calculation using a polarizable con-


J. Phys. Chem. C, Vol. 112, No. 33, 2008 12657

Figure 4. In vacuo reaction pathway for the branched cyclic-trimer hydrolysis and reverse condensation reaction at the Q1, Q2, and Q3 silicon sites. Transition states of each reaction are shown in the figure. See Figure 3 for pictures of the reactant and product for each reaction.

Figure 3. Reactants and products of the hydrolysis and reverse condensation of a branched cyclic-trimer at different silicon sites. Corresponding silicons and oxygen from attacking water are numbered in each reaction. Hydrogens of the cluster are omitted for clarity.

tinuum model (PCM).52–54 Thus every intermediate and transition state was solvated with a first layer of explicit solvent molecules and a surrounding continuum layer of implicit solvent (see Figure 2). The advantage of such a solvation model is that local interactions of the solvent with the solute (i.e., hydrogen bonding) are accounted for by the use of explicit solvent, whereas the long-range interactions are accounted for by the implicit solvent. We refer to this technique in the text as the “hydration sphere model” (HSM). To model the implicit solvent layer, we used the PCM method available in Gaussian 03 with the default parameters for water. To build the molecular cavity within the polarizable continuum, we used radii from the UFF force field and included spheres on individual hydrogen atoms. We tested our hybrid solvation model using the UFF radii and found good agreement with experimentally reported52 solvation free energies for several neutral compounds. Using the HSM technique, we calculated the free energy of solvation according to

∆Gsolv ) Gsolution - Gsolvent - Gsolute where Gsolution is the energy of the entire system, Gsolvent is the energy of the solvent only, and Gsolute is the in vacuo energy of the solute. To obtain the energy of the pure solvent, we removed the solute from the optimized hydration sphere and allowed the explicit solvent molecules to fully relax. We use ∆G solv to denote the free energy of solvation because, similar to the in

vacuo study, the zero-point energy and free energy corrections were added to the electronic energy at the end of the calculations once a frequency analysis was performed. We should be clear here that the solvation free energies we calculate are not intended to be precise measures, as our calculations of free energies did not include configurational entropies of the extended solvent. However, since differences in free energies rather than absolute free energies are reported, we believe that our prescription, although neglecting differences in solvent configuration free energies between pure solvent and solvated silica species, nevertheless provides a good estimate of solvation free energy by including the enthalpic contributions arising from solvent-solute interactions. We operate under the assumption that the dominate contribution to the free energy differences are enthalpic and entropic contributions can be partially included via vibrational modes. This, however, may not include important entropically derived phenomena (including the hydrophobic effect). III. Results and Discussion A. In Vacuo Analysis. With the exception of the cyclictrimer and cyclic-tetramer, the optimized structures for each cluster are in good agreement with previous studies32,33 conducted at similar levels of theory. The lowest energy conformers for each of the single ring structures were found to be nearly planar. However, using a local density functional and modest basis-set, Periera et al. report a chairlike and crown conformation for the geometries of the cyclic-trimer and cyclic-tetramer, respectively. Detailed geometries for each structure, optimized at the B3LYP/6-31G(d,p) level of theory, can be found in the Supporting Information. The sol-gel processes can occur under varying ionic conditions (acidities). For our in vacuo study, we first considered the gas phase hydrolysis reaction in the presence of a hydronium



SCHEME 4: Mechanism of Hydrolysis and Reverse Condensation Reaction of the Dimer Species in the Presence of Explicit Water Showing Initial Deprotonation of the Cluster by the Surrounding Water and Subsequent Reprotonation Prior to Hydrolysis

TABLE 2: Comparison of the Calculated Rate-Limiting Activation Barriers for the Condensation and Reverse Hydrolysis Reactions Using Each of the Three Modelsa rate limiting activation barrier (kcal mol-1) hydrolysis condensation

in vacuo

HSM

PCM

46.9 45.2

14.4 11.4

46.7 44.7

a Energies are calculated at the B3LYP/6-31G(d,p) level of theory.

cation (H3O+). For conciseness we will refer to this as acidic conditions, while the reaction with neutral water (H2O) will be referred to as neutral conditions. We began by placing a hydronium cation (H3O+) at varying distances from the cluster and optimized the entire system allowing all degrees of freedom to fully relax. In all cases, we observed rapid protonation of one of the oxygens of the silica cluster leaving a neutral water molecule and the protonated silica species. Hydrolysis then proceeded through attack of the charged cluster by a neutral water molecule (see Scheme 3). Our calculations show the hydrolysis reaction proceeded through a SN2-type mechanism under both neutral and acidic conditions. We attempted to force the hydrolysis reaction to proceed via the flank-side mechanism by constraining the Si-O bond distance between the silica of the cluster and oxygen of the attacking water molecule for several silica species (dimer, trimer, and cyclic-trimer). In every case, when we removed the constraint and allowed full relaxation of the molecule, the attacking water was found a large distance away from the silica cluster. When we performed the same constrained optimization technique via the SN2 mechanism, we found upon removing the constraint that the molecule relaxed to the hydrolysis product. We found that for the smaller species the condensation and reverse hydrolysis reaction proceeded through a single-step mechanism under neutral conditions and through a two-step mechanism under acidic conditions. Trinh et al.35 found similar results for the condensation reaction under neutral and anionic conditions. We found the three largest species studied in this

work, the cyclic-tetramer, prism, and cubic octamer, proceeded through a single-step mechanism under both neutral and acidic conditions. The single-step neutral mechanism occurred through a pentacoordinated transition state. The transition state species undergoes an intermolecular hydrogen transfer which leads to cleavage of the Si-O bond. For the smaller species, the dimer, trimer, and cyclic-trimer, the condensation and reverse hydrolysis reaction proceeded though two activated processes under acidic conditions. The first step in the condensation reaction corresponded to the formation of a stable pentacoordinated intermediate and the second step corresponded to the removal of water. For the reverse hydrolysis reaction, the first step corresponded to the formation of a pentacoordinated intermediate and the second associated with a hydride transfer, siloxane bond breakage, and an inversion process for the linear species. This mechanism is in agreement with the one proposed by Pereira et al.,21 who observed a stable pentacoordinated intermediate for the silica dimerization reaction under acidic conditions. For the larger species, under acidic conditions, we observed a single-step mechanism similar to that found for the neutral mechanism. Table 1 shows the in vacuo activation energies and overall reaction energies for the condensation and reverse hydrolysis reactions under neutral and acidic conditions for each species. In many of the reactions, the inclusion of zero-point vibrational energy corrections lowered the calculated energy of the transition state and increased the energy of the intermediates, which significantly influenced the size of the activation barrier. The most prominent example of this result can be seen in the calculated activation barriers for the acidic condensation reactions of the dimer, trimer, and cyclic-trimer. For these reactions, the addition of ZPV energy corrections shows the water removal step to be an unactivated process. Under neutral conditions the activation barrier for the formation of the dimer and trimer species is 45.2 kcal/mol, the activation barriers for the formation of the single-ring cyclictrimer and cyclic-tetramer are 24.3 and 24.9 kcal/mol, respec-


J. Phys. Chem. C, Vol. 112, No. 33, 2008 12659

Figure 5. Condensation and reverse hydrolysis reaction energy diagram comparing the in vacuo, PCM, and HSM calculations. All energies are calculated at the B3LYP/6-31G(d,p) level of theory. Pictures shown in the figure are from the in vacuo and PCM calculations. See Scheme 4 for an illustration of each of the transition states and intermediates under solvated conditions using the HSM technique.

TABLE 3: Calculated Solvation Free Energies per Silicon Using the HSM at the B3LYP/6-31G(d,p) Level of Theory cluster

∆Gsolv/Si (kcal mol-1)a

monomer dimer cyclic-trimer trimer cyclic-tetramer prismatic hexamer cubic octamer

-29.11 -16.95 -15.10 -14.99 -12.05 -4.23 -4.82

a Values include free energy correction and ZPV energy corrections.

tively, and the activation barriers for the formation of the prism and cubic octamer are 14.2 and 28.1 kcal/mol, respectively. Under acidic conditions, the activation barriers for the formation of the larger oligomers and cyclic-species are also lower than the barriers for the formation of the smaller linear oligomers, with the exception of the cyclic-tetramer whose activation barrier is approximately equal to that of the dimer species. The activation barriers for the acid-catalyzed formation of the dimer and trimer species are 15.0 and 44.5 kcal/mol,respectively. The activation barriers for the formation of the cyclic-trimer and cyclic-tetramer species under acidic conditions are 6.24 and 15.4 kcal/mol, respectively, and the activation barriers for the formation of the prism and cubic octamer under acidic conditions are 6.40 and 12.7 kcal/mol, respectively. The lower activation barriers for the formation of the larger oligomers relative to the smaller oligomers indicate that the global rate limiting step for formation of the more condensed oligomers is the formation of their smaller constituent species. This trend was observed previously by Trinh et al.;35 however,they found the formation of the ringed species to have higher activation barriers than the dimerization and trimerization reactions. It is also of interest to compare the activation barrier for the condensation reaction to the activation barrier of the reverse hydrolysis reaction for each species. From the results in Table

1, under neutral conditions the activation barriers for the condensation and reverse hydrolysis reaction are nearly equal for each species. Under acidic conditions, however, the difference between the condensation activation barrier and the hydrolysis activation barrier for each species increases with larger oligomers. This trend indicates that as the species oligomerize there is a larger kinetic barrier to hydrolyze back into less-condensed species. The notion of increasing stability with increasing cluster size under acidic conditions is further supported by examining the overall reaction energy for the condensation reaction of each species. Under acidic conditions, the reaction energy for the condensation of the dimer and trimer species is -13.0 and -8.50 kcal/mol, respectively. The reaction energy for the acid-catalyzed condensation reaction of the cyclictrimer, cyclic-tetramer, prism, and cubic octamer is -26.2, -18.4, -27.3, and -19.6 kcal/mol,respectively. This trend shows that under acidic conditions the clusters become more stable as they condense to form larger oligomers. This result supports previous studies theorizing that as silicic acid polymerizes to form higher ordered structures, the stability of the structures increase.55 The results are also in agreement with experimental observation that under acidic conditions the condensation reaction will tend toward larger oligomers in solution.9 We next investigate the reaction energies associated with condensation and reverse hydrolysis occurring at different silica sites by examining the reactions of the branched cyclic-trimer, which contains Q1, Q2, and Q3 silicons. In the case of the branched cyclic-trimer, the Q2 and Q3 silicons are contained in the trimer ring, while the Q1 silicon is contained in the branching monomer chain. There are three distinct reaction pathways available to the branched cyclic-trimer. Hydrolysis via the Q1 silicon produces a monomer and cyclic-trimer species (see Figure 3). Attack of the Q2 and Q3 silicons are both ring-opening reactions and lead to Q13Q31 and Q22Q21 (the linear tetramer) species, respectively (here the subscripts denote the number of Qn silicons contained in the clusters).

12660 J. Phys. Chem. C, Vol. 112, No. 33, 2008 The reaction associated with attack of the Q1 silicon had almost identical activation barriers as the dimer and linear trimer species. Thus, it appears that for the hydrolysis reaction the leaving group has little effect on the size of the activation barrier. The energetics of the ring-opening reaction via attack of the Q2 and Q3 silicons were also similar to those of the unbranched cyclic-trimer species. When we compare the results of the branched cyclic-trimer with the results presented in Table 1, we see that in every case, the condensation reaction of two separate species is a more activated process than an intramolecular condensation reaction. Figure 4 presents the hydrolysis and reverse condensation reaction energy diagram for each of the Qn, n ) 1-3, pathways for the branched cyclic-trimer. B. Solvated Silica Dimer Condensation and Reverse Hydrolysis Reaction. Calculations performed on the silica dimer showed the presence of explicit water molecules had a significant affect on the structure and energetics of the condensation and reverse hydrolysis reaction. The smallest hydration sphere past which additional water molecules had no effect on the structure of the dimer was ∼10 Å in diameter, corresponding to 14 water molecules in the explicit solvent layer. We optimized the structure of the dimer species both in vacuo and using the HSM technique at the B3LYP/631G(d,p) level of theory. The bond lengths between the silica of the dimer cluster and the terminal oxygens were slightly elongated when compared to the in vacuo cluster and the O-H bond distances were significantly longer than those found in the gas phase. The average in vacuo Si-O bond length was 1.647 Å compared to 1.651 Å in solution and the average O-H bond length was 0.966 Å in vacuo compared to 0.990 Å in solution. In fact, it was observed that the fully optimized hydrated dimer was deprotonated by one of the surrounding water molecules leaving Si2O(OH)5O- and H3O+. A detailed structural comparison of the in vacuo and hydrated clusters, optimized at the B3LYP/G-31G(d,p) level of theory, can be found in the Supporting Information. Under solvated conditions, the hydrolysis reaction proceeded via a two-step mechanism. Prior to formation of the pentacoordinated intermediate, the negatively charged cluster was reprotonated by the attacking water leaving a hydroxyl ion, OH-. The hydroxyl ion then formed the stable pentacoordinated intermediate, and the hydrolysis reactions proceeded similar to the in vacuo acidic reaction (see Scheme 4). However, rather than two neutral monomers as the hydrolysis product, the product was a neutral monomer and a negatively charged monomer. The calculated rate-limiting activation barrier for the hydrolysis reaction using the HSM technique was 14.4 kcal mol-1, which was significantly less than our calculated in vacuo barrier of 46.9 kcal mol-1 (see Table 2). Using the HSM, we calculated the rate-limiting activation barrier for the condensation reaction to be 11.4 kcal mol-1. This value is in good agreement with previous theoretical and experimental studies of the dimer condensation reaction which estimate the activation barrier to be 12-15 kcal mol-1.20,21,26,29 We attributed the decrease in activation barrier to several factors occurring in the aqueous environments. The presence of explicit water resulted in (1) elongation the hydroxyl bonds of the cluster, (2) elongation of the siloxane bonds, and (3) greater stabilization of both the reactive intermediates and transition states. Furthermore, the surrounding water molecules acted as proton acceptors to facilitate the reactions. To compare the hybrid implicit/explicit HSM to a continuum (implicit-only) model, we calculated the silicon dimer reaction energies using the polarizable continuum model (PCM)52,56,57

Schaffer and Thomson available in Gaussian 03.36 We began with each molecular configuration (i.e., transition states and reactive intermediates) optimized in vacuo and performed single-point energy calculations after embedding the structures within the dielectric continuum. The calculations were performed at the B3LYP/631G(d,p) level of theory with water as the solvent and including individual solvation spheres around each hydrogen. The need for explicit hydration spheres around individual hydrogen atoms arises because of the elongation of the hydroxyl bonds. By default, the PCM calculation only places spheres around heavy atoms which are meant to enclose the hydrogen atoms bonded to the heavier atoms. This method fails when the hydrogen atoms are significantly far away from the heavier atoms or when hydrogen transfer is occurring, both of which occurred in this case. Figure 5 presents the reaction energy diagram for the silica dimer condensation and reverse hydrolysis reaction for the in vacuo, and solvated cases using both the HSM and PCM techniques. Energies obtained from single-point energy calculations using the PCM were slightly different than the in vacuo analysis; however, the water producing condensation rate limiting activation barrier was still considerably higher than the experimentally determined value. Table 2 shows a comparison of the calculated rate-limiting activation barriers for the condensation and reverse hydrolysis reactions using each of the three models. The energies calculated using the implicit-only (PCM) hydration model follow the same trends as the hybrid implicit/explicit (HSM) model, namely, they show a decrease in energy associated with the transition state and an increase in energy associated with the product. However, the calculated activation energies are very different than those obtained using the HSM and from experiment. C. Solvation Energy Calculations. In addition to reaction energies, we used our solvation model to calculate the free energies of solvation of several silica species in a water solvent. The calculated solvation energy per silicon for each cluster is shown in Table 3. The results show a general decrease in solvation energy per silicon with increasing cluster size. Similar trends in solvation energies have been observed using PCM as well as molecular mechanics methods to describe solvated environments.34 This trend is a result of fewer hydroxyl groups available for hydrogen bonding with the surrounding water as the clusters grow from Q1 to Q3 species. IV. Conclusions We have investigated the in vacuo condensation and reverse hydrolysis reaction of seven silica species using electronic density functional theory. Under neutral conditions, the reactions followed a single-step SN2-mechanism with formation of a pentacoordinated transition state species. Under acidic conditions, the three smallest clusterssthe dimer, trimer, and cyclictrimersfollowed a two-step mechanism. The first step in the condensation reaction corresponded to the formation of a stable pentacoordinated intermediate and the second step corresponded to the removal of water. For the reverse hydrolysis reaction, the first step corresponded to the formation of a pentacoordinated intermediate and the second associated with a hydride transfer, siloxane bond breakage, and an inversion process for the linear species. A single-step mechanism similar to the neutral route was found for the larger silica species under acidic conditions. The results showed relatively high activation barriers for formation of the smaller oligomers compared to the barriers for rings and larger oligomer formation. This result supports previous theoretical and experimental observations that larger

Structure and Reactivity of Prenucleation Silica Species more condensed species should form rapidly once their smaller less-condensed constituents are formed. Our results for the acidic reactions showed an increasing kinetic and thermodynamic favorability toward the more condensed species as the clusters grew in size. After examining the reaction of the branched cyclic-trimer, we inferred the intramolecular condensation reaction (i.e., ring formation) to be a less activated process than the polymerization (chain lengthening) reaction. To study the influence of solvation on the reactions, we used a hybrid implicit/explicit solvation model that explicitly accounted for water in the calculations. The model accounts for local interactions of the solute and solvent (i.e., hydrogen bonding) by optimizing a coordinating shell of explicit solvent molecules around the solute. Long range interactions are accounted for by the use of a polarizable continuum model. Results on the silicon dimer revealed a marked change in both the energetics and mechanism of the condensation and reverse hydrolysis reactions. The presence of explicit water resulted in elongation of the hydroxyl and siloxane bonds of the cluster as well as stabilization of both the reactive intermediates and transition states. The surrounding water also acted as a proton acceptor to facilitate the reactions. The calculated activation barrier for the condensation reaction to form the dimer species was in excellent agreement with previous experimental and theoretical results. Free energy of solvation calculations showed a decrease in solvation energy as cluster size increased. Acknowledgment. The authors acknowledge the National Science Foundation for support of this work through Grant No. CTS-0238989-CAREER. This work was partially supported by the National Center for Supercomputing Applications through proposal ESC030001 [cu.ncsa.uiuc.edu] as well as supercomputing resources at Purdue University. Supporting Information Available: Detailed structural information for in vacuo and hydrated clusters presented in this work. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Davis, M. E. Ordered porous materials for emerging applications. Nature 2002, 417 (6891), 813–821. (2) Cundy, C. S.; Cox, P. A. The hydrothermal synthesis of zeolites: Precursors, intermediates and reaction mechanism. Microporous Mesoporous Mater. 2005, 82 (1-2), 1–78. (3) Cundy, C. S.; Cox, P. A. The hydrothermal synthesis of zeolites: History and development from the earliest days to the present time. Chem. ReV. 2003, 103 (3), 663–701. (4) Burkett, S. L.; Davis, M. E. Mechanisms of Structure Direction in the Synthesis of Pure-Silica Zeolites 0.1. Synthesis of Tpa/Si-Zsm-5. Chem. Mater. 1995, 7 (5), 920–928. (5) Burkett, S. L.; Davis, M. E. Mechanism of Structure Direction in the Synthesis of Pure-Silica Zeolites 0.2. Hydrophobic Hydration and Structural Specificity. Chem. Mater. 1995, 7 (8), 1453–1463. (6) Burkett, S. L.; Davis, M. E. Mechanism of Structure Direction in the Synthesis of Si-Zsm-5 - an Investigation by Intermolecular H-1-Si-29 Cp Mas Nmr. J. Phys. Chem. 1994, 98 (17), 4647–4653. (7) Kirschhock, C. E. A.; Ravishankar, R.; Jacobs, P. A.; Martens, J. A. Aggregation mechanism of nanoslabs with zeolite MFI-type structure. J. Phys. Chem. B 1999, 103 (50), 11021–11027. (8) Kirschhock, C. E. A.; Ravishankar, R.; VanLooveren, L.; Jacobs, P. A.; Martens, J. A. Mechanism of transformation of precursors into nanoslabs in the early stages of MFI and MEL zeolite formation from TPAOH-TEOS-H2O and TBAOH-TEOS-H2O mixtures. J. Phys. Chem. B 1999, 103 (24), 4972–4978. (9) Pelster, S.A.S.W.; Schuth, F. Monitoring Temporal Evolution of Silicate Species during Hydrolysis and Condensation of Silicates Using Mass Spectrometry. J. Am. Chem. Soc. 2006, 125 (13), 4310–4317.

J. Phys. Chem. C, Vol. 112, No. 33, 2008 12661 (10) Knight, C. T. G.; Kirkpatrick, R. J.; Oldfield, E. Two-Dimensional Si-29 Nuclear Magnetic-Resonance Spectroscopic Study of ChemicalExchange Pathways in Potassium Silicate Solutions. J. Magn. Reson. 1988, 78 (1), 31–40. (11) Harris, R. K.; Knight, C. T. G. Si-29 Nuclear Magnetic-Resonance Studies of Aqueous Silicate Solutions 0.5. 1st-Order Patterns in Potassium Silicate Solutions Enriched with Si-29. J. Chem. Soc.-Faraday Trans. 2 1983, 79, 1525–1538. (12) Harris, R. K.; Knight, C. T. G. Si-29 Nuclear Magnetic-Resonance Studies of Aqueous Silicate Solutions 0.6. 2nd-Order Patterns in Potassium Silicate Solutions Enriched with Si-29. J. Chem. Soc.-Faraday Trans. 2 1983, 79, 1539–1561. (13) Harris, R. K.; Knight, C. T. G.; Hull, W. E. Nature of Species Present in an Aqueous-Solution of Potassium Silicate. J. Am. Chem. Soc. 1981, 103 (6), 1577–1578. (14) Harris, R. K.; Knight, C. T. G.; Hull, W. E. Nmr-Studies of the Chemical-Structure of Silicates in Solution. ACS Symp. Ser. 1982, 194, 79– 93. (15) Kinrade, S. D. JCH; Schach, AS; et al., Two substituted cubic octameric silicate cages in aqueous solution. J. Chem. Soc.-Dalton Trans. 2002, 7, 1250–1252. (16) Cho, H.; Felmy, A. R.; Craciun, R.; Keenum, J. P.; Shah, N.; Dixon, D. A. Solution State Structure Determination of Silicate Oligomers by 29Si NMR Spectroscopy and Molecular Modeling. J. Am. Chem. Soc. 2005, 128, 2324–2335. (17) McNeil, K. J.; Dicaprio, J. A.; Walsh, D. A.; Pratt, R. F. Kinetics and Mechanism of Hydrolysis of a Silicate Triester, Tris(2-Methoxyethoxy)Phenylsilane. J. Am. Chem. Soc. 1980, 102 (6), 1859–1865. (18) Zerda, T. W.; Hoang, G. High-Pressure Raman-Study of the Hydrolysis Reaction in Tetramethylorthosilicate (Tmos). J. Non-Cryst. Solids 1989, 109 (1), 9–17. (19) Brinker, C. J. Hydrolysis and Condensation of Silicates - Effects on Structure. J. Non-Cryst. Solids 1988, 100 (1-3), 31–50. (20) Iler, R. K. The Chemistry of Silica; Wiley: New York, 1979. (21) Pereira, J. C. G.; Catlow, C. R. A.; Price, G. D. Silica condensation reaction: an ab initio study. Chem. Commun. 1998, 13, 1387–1388. (22) Klamt, A.; Schuurmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, 5, 799–805. (23) Elanany, M.; Selvam, P.; Yokosuka, T.; Takami, S.; Kubo, M.; Imamura, A.; Miyamoto, A. A quantum molecular dynamics simulation study of the initial hydrolysis step in sol-gel process. J. Phys. Chem. B 2003, 107 (7), 1518–1524. (24) Yokosuka, T.; Kurokawa, H.; Takami, S.; Kubo, M.; Miyamoto, A.; Imamura, A. Development of New Tight-Binding Molecular Dynamics Program to Simulate Chemical-Mechanical Polishing Processes. Jpn. J. Appl. Phys. 2001, 41, 2410–2413. (25) Martin, G. E.; Garofalini, S. H. Sol-Gel Polymerization - Analysis of Molecular Mechanisms and the Effect of Hydrogen. J. Non-Cryst. Solids 1994, 171 (1), 68–79. (26) Rao, N. Z.; Gelb, L. D. Molecular dynamics simulations of the polymerization of aqueous silicic acid and analysis of the effects of concentration on silica polymorph distributions, growth mechanisms, and reaction kinetics. J. Phys. Chem. B 2004, 108 (33), 12418–12428. (27) Garofalini, S. H.; Martin, G. Molecular Simulations of the Polymerization of Silicic-Acid Molecules and Network Formation. J. Phys. Chem. 1994, 98 (4), 1311–1316. (28) Feuston, B. P.; Garofalini, S. H. Oligomerization in Silica Sols. J. Phys. Chem. 1990, 94 (13), 5351–5356. (29) Gordon, M. S.; Davis, L. P.; Burggraf, L. W. The Structure and Stability of Neutral Penta-Coordinated Silicon-Compounds. Chem. Phys. Lett. 1989, 163 (4-5), 371–374. (30) Lasaga, A. C.; Gibbs, G. V. Abinitio Quantum-Mechanical Calculations of Water-Rock Interactions - Adsorption and Hydrolysis Reactions. Am. J. Sci. 1990, 290 (3), 263–295. (31) Okumoto, S.; Fujita, N.; Yamabe, S. Theoretical study of hydrolysis and condensation of silicon alkoxides. J. Phys. Chem. A 1998, 102 (22), 3991–3998. (32) Pereira, J. C. G.; Catlow, C. R. A.; Price, G. D. Abinitiostudies of silica-based clusters. Part I. Energies and conformations of simple clusters. J. Phys. Chem. A 1999, 103 (17), 3252–3267. (33) Pereira, J. C. G.; Catlow, C. R. A.; Price, G. D. Ab initio studies of silica-based clusters. Part II. Structures and energies of complex clusters. J. Phys. Chem. A 1999, 103 (17), 3268–3284. (34) Catlow, C. R. A.; Coombes, D. S. ; Lewis, D. W. ; Pereira, J. C. G.; Slater, B. Modeling Nucleation and Growth in Zeolites. In Handbook of Zeolite Science and Technology; Auerbach, S. M., Carrado, K. A., Dutta, P. K., Eds.; Dekker: New York, 2003; pp 91-128. (35) Trinh, T. T.; Jansen, A. P. J.; vanSanten, R. A. Mechanism of Oligomerization Reactoins of Silica. J. Phys. Chem. B 2006, 110, 23099– 23106. (36) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.;

12662 J. Phys. Chem. C, Vol. 112, No. 33, 2008 Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (37) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98 (7), 5648–5652. (38) McLean, A. D.; Chandler, G. S. Contracted Gaussian-Basis Sets for Molecular Calculations 0.1. 2nd Row Atoms, Z)11-18. J. Chem. Phys. 1980, 72 (10), 5639–5648. (39) Zhao, Y.; Truhlar, D. G. Benchmark Databases for Nonbonded Interactions and Their Use To Test Density Functional Theory. J. Chem. Theory Comput. 2005, 1, 415–432. (40) Nicholas, J. B.; Winans, R. E.; Harrison, R. J.; Iton, L. E.; Curtiss, L. A.; Hopfinger, A. J. Abinitio Molecular-Orbital Study of the Effects of Basis Set Size on the Calculated Structure and Acidity of Hydroxyl-Groups in Framework Molecular-Sieves. J. Phys. Chem. 1992, 96 (25), 10247– 10257. (41) Nicholas, J. B.; Winans, R. E.; Harrison, R. J.; Iton, L. E.; Curtiss, L. A.; Hopfinger, A. J. An Abinitio Investigation of Disiloxane Using Extended Basis-Sets and Electron Correlation. J. Phys. Chem. 1992, 96 (20), 7958–7965. (42) Almenningen, A.; Bastiansen, O.; Ewing, v.; Hedberg, K.; Tretteberg, M. Acta Chem. Scand. 1963, 17, 2455–2460. (43) Peng, C. Y.; Ayala, P. Y.; Schlegel, H. B.; Frisch, M. J. Using redundant internal coordinates to optimize equilibrium geometries and transition states. J. Comput. Chem. 1996, 17, 49. (44) Peng, C. Y.; Schlegel, H. B. Israel J. Chem 1993, 33, 449. (45) Schlegel, H. B. Optimization of equilibrium geometries and transition structures. J. Comput. Chem. 1982, 3, 214.

Schaffer and Thomson (46) McQuarrieD. A. SimonJ. D. Molecular Thermodynamics; University Science Books: Sausalito, CA, 1999. (47) Aleman, C. Hydration of cytosine using combined discrete/SCRF models: influence of the number of discrete solvent molecules. Chem. Phys. 1999, 244, 151–162. (48) Claverie, P.; Daudey, J. P.; Langlet, J.; Pullman, B.; Piazzola, D. Studies of Solvent Effects. 1. Discrete, Continuum, and Discrete-Continuum Models and Their Comparison for Some Simple Cases: NH4+, CH3OH, and Substituted NH4+. J. Phys. Chem. 1978, 82 (4), 405. (49) Knapp-Mohammady, M.; Jalkanen, K. J.; Nardi, F.; Wade, R. C.; Suhai, S. L-Alanyl-L-alanine in the zwitterionic state: structures determined in the presence of explicit water molecules and with continuum models using density functional theory. Chem. Phys. 1998, 240, 63–77. (50) Bandyopadhyay, P.; Gordon, M. S. Acombined discrete/continuum solvation model: Application to glycine. J. Chem. Phys. 2000, 113 (3), 1104– 1109. (51) Cui, Q. Combining implicit solvation models with hybrid quantum mechanical/molecular mechanical methods: A critical test with glycine. J. Chem. Phys. 2002, 117 (10), 4720–4728. (52) Cossi, M.; Barone, V.; Cammi, R.; Tomasi, J. Ab initio study of solvated molecules: a new implementation of the polarizable continuum model. Chem. Phys. Lett. 1996, 255, 327–335. (53) Miertus, S.; Scrocco, E.; Tomasi, J. Electrostatic interaction of a solute with a continuum. A direct utilizaion of AB initio molecular potentials for the prevision of solvent effects. Chem. Phys. 1981, 55, 117–129. (54) Miertus, S.; Tomasi, J. Approximate evaluations of the electrostatic free energy and internal energy changes in solution processes. Chem. Phys. 1982, 65, 239–245. (55) Brinker, C. J.; Scherer, G. W., Sol-Gel Science: The Physics and Chemistry of Sol-Gel Processing. Academic Press: London, 1990. (56) Cossi, M.; Scalmani, G.; Rega, N.; Barone, V. New developments in the polarizable continuum model for quantum mechanical and classical calculations on molecules in solution. J. Chem. Phys. 2002, 117 (1), 43– 54. (57) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Energies, Structures, and Electronic Properties of Molecules in Solution with the C-PCM Solvation Model. J. Comput. Chem. 2002, 24 (6), 669–681.

JP066534P

Density Functional Theory Investigation into Structure and Reactivity of

Recommend Documents