Elucidating the Role of Tetraethylammonium in the Silicate Condensation Reaction from Ab Initio Molecular Dynamics Simulations.

The understanding of the formation of silicate oligomers in the initial stage of zeolite synthesis is important. The use of organic structure-directing agents (OSDAs) is known to be a key factor in the formation of different silicate species and the final zeolite structure. For example, tetraethylammonium ion (TEA+) is a commonly used organic template for zeolite synthesis. In this study, ab initio molecular dynamics (AIMD) simulation is used to provide an understanding of the role of TEA+ in the formation of various silicate oligomers, ranging from dimer to 4-ring. Calculated free-energy profiles of the reaction pathways show that the formation of a 4-ring structure has the highest energy barrier (97 kJ/mol). The formation of smaller oligomers such as dimer, trimer, and 3-ring has lower activation barriers. The TEA+ ion plays an important role in regulating the predominant species in solution via its coordination with silicate structures during the condensation process. The kinetics and thermodynamics of the oligomerization reaction indicate a more favorable formation of the 3-ring over the 4-ring structure. The results from AIMD simulations are in line with the experimental observation that TEA+ favors the 3-ring and double 3-ring in solution. The results of this study imply that the role of OSDAs is not only important for the host-guest interaction but also crucial for controlling the reactivity of different silicate oligomers during the initial stage of zeolite formation.

Introduction

Zeolites are nanoporous aluminosilicate materials widely used in various industrial applications making use of their catalytic and separation properties.[1] Zeolites are typically synthesized from aqueous gel solutions containing various heteroatomic compounds, with inorganic and/or organic cations acting as directing agents of the structure and mobilizing agents (hydroxyl or fluoride anions). Numerous experimental[2−10] studies have focused on the nature and structure of the silicate oligomers in solution, as understanding the formation of silicate oligomers in the initial stage is key to zeolite synthesis.[11,12] The elementary steps for Si(OH)4 oligomerization were extensively studied in computational studies[11,13−27] using a continuum or explicit model of water. A common pathway of the oligomerization reaction is a two-step mechanism with an initial formation of a penta-coordinated intermediate, followed by a water removal stage.[20,21,28−31] Earlier studies (e.g., refs[20, 21, 32]) have shown that it is crucial to include the effect of thermal motion and the presence of explicit water molecules when modeling aqueous chemical reactions that involve solvent molecules that strongly bind to the reagents, or actively participate in the reaction mechanism. The overall picture of free-energy profiles and mechanism could change significantly with dynamic and explicit treatment of solvent.

Cationic organic structure-directing agents (OSDAs) are known to be important as inferred from various experimental studies.[33−36] The use of different organic templates such as tetramethylammonium (TMA+), tetraethylammonium (TEA+), and tetrapropylammonium (TPA+) leads to distinct dominant structures.[8,37,38] For example, TEA+ is crucial in the crystallization of zeolite β-like structures.[39−42] The interaction between TEA+ and zeolite framework was identified as the main driving force for the stability.[43] In the very first stage of silicate oligomer formation in solution, double 3-ring (D3R) and double 4-ring (D4R) structures were observed.[8,38,44] With excess of organic cation, the structures of D4R.8TMA+ and D3R.6TEA+ are the most stable species in solution.[34,35] Interestingly, the experimental study conducted by Inagaki et al.[45] stated that the D3R structure in the presence of TEA+ could be transformed into other structures such as secondary building units (SBUs) of type “4–2”. Thus, the formation of D3R is critical for the initial stage of zeolite formation, but it might not be essential for the final stable zeolite beta structure to be crystallized. Caratzoulas et al. proposed, on the basis of classical molecular dynamics simulations, that the D3R and D4R structures have different stabilities when interacting with 6TMA+ or 8TMA+.[46][46] The effect of organic templates has been investigated to determine the interaction between the silicate oligomer and organic template.[25,33,46,47]Ab initio molecular dynamics studies have addressed the silicate dimerization mechanism in the presence of ODSAs such as TPA+ and TMA+. These studies showed that the activation barrier of dimerization increases in the presence of those cations. However, a comprehensive picture of the role of TEA+ in the formation of silicate oligomers is still lacking, in particular on a molecular level.

In this work, ab initio molecular dynamics (AIMD) simulations were performed to study the formation of silicate oligomers in the presence of TEA+ in aqueous solution, incorporating the water molecules explicitly. The free-energy profiles of the formation pathways of different silicate oligomers were obtained. The results show that the pathway for the 4-ring formation is favorable over that of the 3-ring formation in the presence of TEA+. More interestingly, during the reaction, TEA+ molecules prefer to form a complex with most of the silicate structures (dimer, 3-ring, 4-ring, linear and branched tetramer), except the trimer species. In contrast to the case of TMA+, the presence of TEA+ shows that the formation of the 3-ring structure is more favorable than that of the 4-ring. This trend is in excellent agreement with experimental observation.[8,38] From the present study, we can infer that the formation of D3R or D4R can be controlled by the single ring formation step.

Simulation Method

Our computational setup was similar to that of earlier computational studies of silicate oligomerization reactions in aqueous solution.[20,21,28] The ab initio molecular dynamics simulation is based on a density functional theory (DFT) description of the electronic structure. We employed the CP2K package[48] to carry out the molecular dynamics simulation using a Born–Oppenheimer approach as implemented in the Quickstep module.[49] Here, Goedecker–Teter–Hutter (GTH) pseudopotentials[50,51] are used to account for the interactions between the electrons and the atomic nuclei. The BLYP exchange–correlation functional[52,53] was used, together with Grimme’s type D2[54] to account for the long-range van der Waals interactions. The double-zeta valence DZVP-MOLOPT basis set complemented with polarization functions[55] was employed for all atom types. An energy cutoff of 400 Ry was chosen for the auxiliary plane-wave basis set. The molecular dynamics trajectories were generated with a time step of 0.5 fs. We applied a velocity rescaling thermostat[56] with a time constant of 300 fs to impose the temperature of the system set at 350 K.

The simulation cell was a periodic orthorhombic box (12 × 12 × 25 Å3) with a density of about 1 g/cm3, similar to that of the experimental value at ambient conditions. The initial geometry of silicate oligomer and organic cation TEA+ was first optimized in the gas phase. This structure was then solvated with around 140 water molecules evenly distributed in the simulation box. Subsequently, to generate a representative configuration, a 20 ps equilibration run was performed in the NVT ensemble. The total number of atoms in the system was approximately 450 atoms. Due to the high simulation cost of ab initio MD, we did not consider with systems with a higher concentration of silicate and TEA+ in the present study.

Reaction pathways were obtained by tracing a proper reaction coordinate using the method of constraints.[57,58] For each value of the reaction coordinate, the initial configuration was taken from the last configuration of the simulation at the previous value of the reaction coordinate. After 1 ps of equilibration, a 10 ps trajectory was generated to gather data. The total trajectory of the simulations of a reaction pathway was around 200 ps, distributed typically over 20 values of the reaction coordinate.

The free-energy (ΔG) profiles of the oligomerization reactions were obtained by thermodynamic integration using eq where F is the calculated constraint force and r is a fixed value of the reaction coordinate. The intergral is evaluated numerically on the basis of the calculated values of the constraint force at each of the (∼20) reaction coordinate values. The errors of the constraint force are typically below 10–5 Hartree/Bohr in a 10 ps production run. This approach has generic applicability, and it is extensively used in earlier studies to calculate free-energy barrier reactions in solution.[20,59−62]

A common two-step mechanism of silicate condensation reaction in the basic conditions[28] is described in Scheme . The first step is to form a fivefold coordinated intermediate with OSi–O bonding. For this stage of the reaction pathway, the distance between O3 atom and Si2 was taken as a reaction coordinate. Here, O3 atom is defined as the reactive oxygen. The second stage consists of a water removal step, where the distance between Si2 and O4 was taken as the reaction coordinate. For ring closure reactions, a similar mechanism has been considered.[28,63] Note that in the ring closure reaction the silicon and oxygen atom in the first reaction step are of the same oligomer molecule. We investigated six oligomerization reactions, from the formation of a dimer up to a 4-ring structure. A schematic process of the reactions is provided in Figure . For each reaction, ab initio MD with an explicit water model was used to calculate the free-energy profile and elucidate the reaction mechanism.

Scheme 1

Representation of a Two-Step Mechanism of Silicate Condensation Reaction in the Presence of TEA+. R, TS1, I, TS2, and P Refer to Reactant, Transition State 1, Intermediate, Transition State 2, and Product, Respectively

Figure 1

Left: scheme of the oligomerization reactions considered in this work forming from dimer to 4-ring formation. Right: snapshot of ab initio MD simulations for dimerization reaction. This system consists of a dimer silicate, TEA+, and 132 water molecules.

Results and Discussion

Radial Distribution Functions

We performed an unconstrained 20 ps AIMD simulation of the TEA+–silicate system and analyzed the radial distribution function (RDFs) to collect basic structural information for the water, silicate, and organic template in aqueous solution. The RDFs of oxygen–oxygen (Ow–Ow) and oxygen–hydrogen (Ow–Hw) among water molecules and those of silicon–oxygen (Si–Ow) and nitrogen–oxygen (N–Ow) are shown in Figure . For the Ow–Ow RDF, the first peak is at 2.8 Å, which is in good agreement with the experimental data and earlier simulations.[63,64] This implies that the presence of TEA+ and silicate has no effect on the water structure. The first peak of Si–Ow RDF is located at 3.75 Å for TEA+–water, which is very similar to that of TMA+–water. However, the location of N–Ow RDF peak for TEA+–water is at 5.8 Å, which is further than for TMA+–water (4.7 Å).[63] This is not unexpected, as TEA+ is larger than TMA+, yielding a larger distance water–nitrogen distance in the case of TEA+.

Figure 2

Radial distribution function (RDF) for Ow–Ow, Ow–Hw, Si–Ow, and N–Ow for a TEA+ and silicate solvated in water, as obtained by unconstrained ab initio MD simulations.

Formation of Linear Silicate Oligomers

Figure presents representative snapshots of the silicate dimerization reaction in close contact with TEA+. The reaction mechanism of silicate condensation is very similar to that reported in earlier studies.[28,63] As described in Scheme , the first reaction step is to form SiO–Si bond resulting in a fivefold silicate intermediate. The second step is the removal of the water to form the dimer product. The calculated free energies of the transition states, intermediate, and product are listed in Table . Taking the reactant as the reference, the free-energy value of the first reaction barrier is 63 kJ/mol. The second activation barrier, calculated as the free-energy difference between the transition state TS2 and the intermediate, is 24 kJ/mol. The resulting overall activation barrier of the dimer formation is 69 kJ/mol, which is in a similar range to that reported in previous computational studies.[18,21,28] It is known that the presence of counterions would increase the overall reaction barrier.[21] In this case, TEA+ raises the total activation barrier of dimerization by only 8 kJ/mol compared with the case without cation. This effect is minor compared to the effect of the presence of TMA+, which yields an increase of 15 kJ/mol.

Figure 3

Representative snapshot from ab initio MD simulations of the dimerization reaction with the mechanism in Scheme . The organic template TEA+ stays close to the silicate structure along the whole reaction pathway.

Table 1

Calculated Free-Energy (kJ/mol) Profiles along the Silicate Formation in the Presence of TEA+ Obtained by Ab Initio MD

free energy	reactant	TS1	intermediate	TS2	product
dimer	0	63	45	69	1
trimer	0	69	55	83	19
linear tetramer	0	63	49	85	18
3-ring	0	67	55	82	20
4-ring	0	75	68	97	53
branched tetramer	0	67	51	80	23

When looking at the free-energy barriers of the second step for the formation dimer, trimer, and linear tetramer, we see that these are in the range of 24–36 kJ/mol. These values, as well as the reaction mechanism, are very are similar to those observed in simulations of systems with other cations such as Li+, NH4+,[21] and Na.+.[28] The leaving hydroxyl group forms well-defined hydrogen bonds with water molecules. It is protonated either directly by another silicate hydroxyl group or via a proton transfer chain mediated by one or more water molecules.[63,65] The presence of ODSAs such as TEA+, TMA+, or TPA+ seems to have a minor effect on the water removal step.[30,63] This is in consistent with previous studies showing that this step is more determined by the hydrogen-bonding network between silicate and water.[65]

Analysis of the trajectories of the formation of dimer, trimer, and linear tetramer suggests that the active oxygen in the first step (O3 in Scheme ) has no direct contact with counterion TEA+. This can be observed in the snapshots of the dimer formation in Figure . In earlier studies, it has been shown that the direct interaction with this active oxygen in the first step substantially increased the first activation barrier and hence the overall barrier.[21,28] In this case, it is apparent that hydrophobic cation such as TEA+ favors greater distance coordination. Therefore, TEA+ has less effect on the activation barrier than a relatively small inorganic cation does (figure ).

Figure 4

Calculated free-energy profile (kJ/mol) of formation of linear silicate oligomer as a function of reaction coordinate.

When comparing free-energy values for the formation of linear species such as dimer, trimer, and tetramer, as listed in Table , we see that the activation energy is largest for the linear tetramer and lowest for the dimer. The relative stability of the dimer appears to be larger than that of the trimer and the linear tetramer. This indicates that the rate-limiting step for linear growth of silicate, in the presence of TEA+, is the formation of the linear tetramer. This trend is opposite to the case of TMA+, where the formation of the linear tetramer has the lowest overall activation barrier among the linear species.[63]

Table 2

Total Free-Energy Barriers (kJ/mol) Obtained by Ab Initio MD of Silicate Oligomerization Reaction in the Presence of TEA+a

free-energy barrier	TEA+	TMA+	without cation
	(this work)	ref[63]	ref[32]
dimer	69	78	61
trimer	83	94	53
linear tetramer	85	74
3-ring	82	89	72
4-ring	97	80	95
branched tetramer	80	72	101

The energies without cation[32] and with TMA+[63] are added for comparison.

Comparison of the present results with those obtained for systems with inorganic cations[21,28] gives the following picture. The relative barrier heights for dimer and trimer formation in the presence of TEA+ are the opposite of that found in an NH4+ cation system, where the dimer formation barrier is considerably higher than that of the trimer formation. In the presence of Na+, the dimer and trimer barrier are comparable, while the Li+ cation tends to have the same effect as TEA+ on the relative heights of barrier for the dimer and trimer species. Thus, the results indicate a different dominant linear species in silicate growth in the presence of organic and inorganic counterions.

Formation of Ring and Branched Silicate Oligomers

A critical step in zeolite synthesis is the formation of branched or ring structures at the initial stage of silicate production. The creation of initial 3-ring or 4-ring structures in particular is a determining factor for the zeolite’s final ring structure. The mechanism of the formation of 3-ring and 4-ring configurations is comparable to that seen in the simulations of linear structure formation. Details of this ring closure reaction were elucidated in earlier computational studies.[28,63]Figure shows snapshots of the calculated 3-ring reaction pathway in the presence of cation TEA+. In the first step, the active oxygen atom from one end of the linear tetramer binds to the Si atom at the other end, yielding an intermediate ring. In the second step, the removal of water results in the final product. It is interesting that the leaving hydroxyl group is protonated by means of an external transfer mechanism, receiving the proton from another water in the surrounding. This is due to the specific arrangement of the hydrogen-bond network around the silicate when the hydroxyl group moves away from the fivefold Si intermediate.

Figure 5

Representative snapshot of ab initio MD simulations of the 3-ring closure reaction.

The free-energy profiles of the formation of the 3-ring, 4-ring, and branched tetramer are shown in Figure , with the numerical values listed in Table . The branched tetramer has the lowest free-energy barrier (80 kJ/mol), similar to that for the formation of the 3-ring (82 kJ/mol). The 4-ring formation has a significantly higher barrier than the 3-ring formation. This trend is opposite to that observed in simulations of the TMA+ cation system,[63] where the 3-ring formation has a higher barrier than the 4-ring formation. The results clearly indicate that the 4-ring formation in the presence of TEA+ is less favorable, showing a higher free energy of the intermediate, transition state, and product than the 3-ring and branched tetramer pathways.

Figure 6

Calculated free-energy profile of formation of ring and branched silicate oligomer as a function of the reaction coordinate. The profiles calculated for the two stages are connected in the graph.

The observation that reaction free energies are positive in this work complies with previous theoretical studies.[21,30,65] The reason is that the overall reaction produces one extra molecule of water, which generates an entropically unfavorable rearrangement of the structure of water. The free-energy results presented in Table imply that, in the presence of TEA+, the dimerization reaction has the highest rate. In contrast, in the presence of TMA+, the formation of the linear tetramer and branched tetramer is kinetically the most favorable.[63] The formation of 4-ring appears to be unfavorable with the highest in free-energy barrier (97 kJ/mol) compared to the formation of other silicate structures (Table ). Our findings suggest that the formation of 4-rings is suppressed in the presence of TEA+ in the process of zeolite synthesis. This observation is in good agreement with experimental results, where the formation of 3-ring and D3R is the most dominant species with the organic TEA+ template.[8,38]

Interaction between TEA+ and Silicate Oligomers

In earlier computational studies,[30,63] the variation in the height of the barrier was correlated with the relative position of the cation in the reacting species. In this study, we also surveyed the relative position of TEA+ by calculating the distribution of the distance between the TEA+ nitrogen atom and the nearest Si atom of silicate. The averages of this N–Si distance distribution are plotted against the coordinates of the reaction and illustrated in Figure . It is interesting to note that TEA+ stays close to the silicate during dimer and linear tetramer formations, whereas in the trimerization reaction, it departs from the silicate at the beginning of the reaction. To verify this behavior, we analyzed the trajectory of the simulation with the TEA+ and silicate, which were initially located close to each other. After about 8 ps of simulation, the pair TEA+ and reactant state of trimer begins to depart, as seen in Figure . This indicates that during the trimer formation, there is enhanced stability of separately solvated species, which is related to the hydrogen-bond structure. This has also been observed for silicate formation in the presence of TPA+.[30] On the other hand, it seems that binding of TEA+ to the dimer and linear tetramer gives rise to an enhanced stabilization of the intermediate structures. Table shows that the free energy of the intermediate state for the formation of the dimer and linear tetramer is slightly lower than that for the trimer. The TEA+ does not participate actively in the reaction to water removal. Therefore, the connection between the presence of TEA+ in the silicate species’ first coordination shell and the barrier to the water removal reaction may not be important because this part of the reaction involves a mechanism of proton transfer mediated by the silicate hydrogen-bond network.

Figure 7

Shortest distance between nitrogen atom of TEA+ and Si of silicate as a function of reaction coordinate for linear silicate formation. The vertical bars indicate the variation (standard deviation) of the measured distances. Only in the case of trimer, TEA+ separates from the silicate during the reaction.

Figure 8

Shortest distance between TEA+ and trimer in an independent ab initio MD simulation. The two components were initially close together. A separation between TEA+ and silicate was observed around 8 ps of simulation.

Figure shows the shortest distance between the Si atom of the silicate and the N atom of TEA+ as a function of the reaction coordinate in the formation of ring and branched structures. The position of TEA+ is relatively stable during the whole reaction process. TEA+ is positioned near the silicate oligomer. In the case of 3-ring closure reaction, we did not observe a separation between TEA+ and 3-ring. This trend is opposite to the case of TMA+, while TMA+ tends to be dissociated from the silicate during almost the full reaction pathway.[63]

Figure 9

Shortest distance between the nitrogen atom of TEA+ and Si of silicate as a function of reaction coordinate for ring and branched oligomer formation. The vertical bars indicate the variation (standard deviation) of the measured distances. During the reaction, TEA+ remains mostly close to the silicate.

To better assess the stability of binding of TEA+ to the 3-ring and 4-ring, we performed additional simulations starting with a configuration with the TEA+ that is in close contact with the 3-ring and 4-ring structures (not shown here). During a 15 ps NVT simulation, no separation process was observed in both cases. An interesting observation is that TEA+ remains at the face of both 3-ring and 4-ring structures. While in the presence of TMA+, cation TMA+ moves from the facial position to the edge of 3-ring.

Formation of 3-Ring vs 4-Ring Oligomers

The formation of a single 3-ring and 4-ring as well as the formation of D3R and D4R is an important notion in zeolite synthesis, as indicated before. An early computational study provided evidence that a 3-ring is less favorable than 4-ring in the absence of a counterion.[17] The higher stability of a 4-ring structure was also observed in the presence of TMA+.[63] However, the present paper reveals the interesting finding that the 3-ring is more stable than a 4-ring structure in the presence of an ODSA cation TEA+. In the presence of TEA+, the free-energy profile provided in Table indicates that the free-energy barrier of 3R formation of both TS1 and TS2 is lower than that of 4-ring formation. As a result, the total activation barrier of 3-ring (82 kJ/mol) is lower than that of 4-ring (97 kJ/mol) by 15 kJ/mol. In addition, the calculated thermodynamic reaction of the 3-ring is also 33 kJ/mol more favorable than the 4-ring formation. These values (Table ) are 20 and 53 kJ/mol for 3-ring and 4-ring, respectively. Earlier simulation[46,47] and experimental[44] studies have shown that the cation TMA+ stabilizes the D4R structure over the D3R. However, in the case of TEA+, a higher stability of the 3-ring/D3R over the 4-ring/D4R was observed.[8,38,45] A stable and neutral charged cluster of D3R with excess 6TEA+ can also be formed.[40] The present results clearly show the more favorable formation of 3-ring over 4-ring in the presence of TEA+ and are consistent with these observations.

Conclusions

The formation of silicate oligomers from dimer to 4-ring in the presence of the organic counterion TEA+ has been studied using ab initio molecular dynamics simulations of a model that incorporates explicit water molecules. The results show that the presence of TEA+ increases the free-energy barriers of all reactions compared to the case without a counterion, consistent with earlier computational studies of systems with other cations.[28,63] The formation of the dimer appears to have the lowest free-energy barrier. An important observation is that the formation of 4-ring is the rate-limiting step in silicate growth due to the relatively high free-energy barrier in the presence of TEA+. In contrast to the case of TMA+, the presence of TEA+ favors the formation of 3-ring over 4-ring. The free-energy barrier of the 3-ring is lower than that of 4-ring formation by 15 kJ/mol.

An interesting observation is that the simulations show that TEA+ does not always directly coordinate with the silicates: the electrostatic interaction may also give rise to an association on a longer range, yielding a well-defined hydrogen-bond structure of the water molecules with negatively charged silicate. This is consistent with the results from classical force-field molecular dynamics simulations.[47] More specifically, the organic cation TEA+ coordinates to the silicate oligomers during the formation of the dimer, branched and linear tetramer, 3-ring, and 4-ring. However, in the case of trimer formation, the TEA+ and silicate are dissociated along the full reaction pathway. Further study on this dissociation in aqueous solution is needed to fully understand the interaction between ODSA and silicate oligomers.

The finding that linear tetramer and 4-ring formation has a higher transition state barrier than 3-ring and trimer formation suggests a dominant appearance of 3-ring structures in the presence of TEA+. Moreover, the calculated reaction free energies demonstrate also that 3-ring formation is thermodynamically more stable than 4-ring formation. This may provide the rationale of the experimental observation that D3R is more favorable than the D4R structure.[8,38] A more definite conclusion requires further investigation to evaluate the free-energy profiles of the single- to double-ring structure reaction.

In conclusion, our findings provide important insights into how TEA+ plays a directing role in silicate oligomerization by regulating those reactions’ thermodynamic and kinetic parameters. TEA+ favors the formation of small oligomers by near contact with silicate structures and the modification of hydrogen bonding in the solvation shell. The results provide a solid basis for further study of the assembly double-ring structures in solution from a single ring. It also provides the essential input for a larger-scale simulations study using kinetic Monte Carlo simulation. These can address silicate and TEA+ concentration and will provide a more detailed picture allowing for a better comparison with experimental observations.