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., refs20, 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 + . 4646 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/cm 3 , 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 The free-energy (ΔG) profiles of the oligomerization reactions were obtained by thermodynamic integration using eq 1 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 1. The first step is to form a fivefold coordinated intermediate with OSi−O bonding. For this stage of the reaction pathway, the distance between O 3 atom and Si 2 was taken as a reaction coordinate. Here, O 3 atom is defined as the reactive oxygen. The second stage consists of a water removal step, where the distance between Si 2 and O 4 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 1. For each reaction, ab initio MD with an explicit water model was used to calculate the free-energy profile and elucidate the reaction mechanism.
■ 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 nitro-gen−oxygen (N−Ow) are shown in Figure 2. 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 + .
Formation of Linear Silicate Oligomers. Figure 3 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 1, 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  1. 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  The Journal of Physical Chemistry B pubs.acs.org/JPCB Article cation. This effect is minor compared to the effect of the presence of TMA + , which yields an increase of 15 kJ/mol. 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 + , NH 4 + , 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 (O 3 in Scheme 1) has no direct contact with counterion TEA + . This can be observed in the snapshots of the dimer formation in Figure 3. 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 4).
When comparing free-energy values for the formation of linear species such as dimer, trimer, and tetramer, as listed in Table 2, 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 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 Journal of Physical Chemistry B pubs.acs.org/JPCB Article the presence of TEA + are the opposite of that found in an NH 4 + 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 5 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.
The free-energy profiles of the formation of the 3-ring, 4ring, and branched tetramer are shown in Figure 6, with the numerical values listed in Table 2. 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.
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 2 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 2). 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 7. It is interesting to note that TEA + stays close to the silicate during dimer and linear tetramer The energies without cation 32 and with TMA +63 are added for comparison. The Journal of Physical Chemistry B pubs.acs.org/JPCB Article 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 8. 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 1 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 9 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 To better assess the stability of binding of TEA+ to the 3ring 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.   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 freeenergy profile provided in Table 1 indicates that the freeenergy 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 1) 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 4ring 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 forcefield 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 freeenergy 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.