Methanol-to-hydrocarbon initiation reactions over a zeolite catalyst: reactive molecular dynamics simulations

Abstract

Environmental concerns and energy security drive the need for renewable alternatives to fossil fuels. Methanol is an ideal candidate as it enables the synthesis of olefins and fuels from biomass using zeolite-catalyzed processes, which inspired the concept of “Methanol Economy” since it ideally implies an oil-independent scenario. Although MTH industrial units already exist, many fundamental aspects remain unknown. Therefore, in this contribution, we developed a ReaxFF reactive force field to study the dynamic features of the equilibration stage of the MTH process in an H-ZSM-5 zeolite. Simulations were run from 600 to 1200 K during 1000 ps, using a constant number of molecules, pressure, and temperature (NPT ensemble). Methanol conversion increases from 800 to 1000 K to form water and the crucial intermediate surface methoxy species (SMS), whose production diminishes at 1200 K because of the prevalence of undesired methane. Generally, temperatures above 1200 K can lead to questionable reactions due to entropy effects. Humidity at 800 K modifies the nature of the zeolite acidity from static to dynamic, embodied in hydronium ions, which enhances methanol conversion via hydrogen transfer reactions and framework activation, namely, water protonation leaves a negatively charged framework that eventually facilitates the dissociation of protonated methanol in water and an SMS. Cation diffusion was pervasive, and it is hypothesized that this relieves the entropic penalties of several relevant reactions. This phenomenon showcases a critical manifestation of dynamic effects that complements experimental and theoretical research mainly conducted by static density functional theory methods. Overall, dynamic effects involving diffusing cations and modification of the zeolite interior by water are complex and call for more extensive investigations.

---

1. Introduction

The world faces the challenge of transitioning from constrained fossil fuels to renewable energy resources.¹ Among the alternatives, solar and wind power generation have significantly expanded worldwide and are expected to experience substantial growth from 2025 to 2035.² However, these technologies necessitate long-duration storage devices to manage extended unfavorable seasonal changes.³ Hydrogen is a promising energy carrier since it can be produced from the electrolysis of water using any renewable energy. The downside of hydrogen-based storage systems is that they must overcome risks related to transportation and handling, which requires optimization of over twenty properties before they can serve as an adequate substitute for liquid fossil fuels.⁴ In any case, hydrogen will undoubtedly be an essential part of the energy economy. Alternatively, fuel property predictions have shown methanol as an ideal candidate to improve engine combustion systems.⁵

Methanol is a liquid under environmental conditions, which enables its distribution through pipeline networks and allows for large-scale storage in tanks without the inherent hazards associated with hydrogen.⁶ Over 90% of methanol comes from natural gas, even though significant efforts have focused on its production from renewable biomass and waste via gasification,^7,8 as well as from the reaction of CO₂ with renewable hydrogen.⁹ Although methanol and its blends have been directly used as fuels in internal combustion engines,¹⁰ over 70% is used to produce other chemicals. In particular, the synthesis of hydrocarbons like olefins and gasoline through processes catalyzed by zeolites is the subject of comprehensive research because this approach offers an energy scenario independent of oil, the so-called “Methanol Economy.”¹¹

In 2010, the Dalian Institute of Chemical Physics of China developed, constructed, and commissioned the first industrial unit of methanol conversion towards olefins, that is, the methanol-to-olefin (MTO) process, after almost thirty years of research and development.¹² Despite the initial industrial success, many fundamental aspects of zeolite-catalyzed processes remain unknown and are the subject of intense investigations due to the complexity of their reaction mechanisms. Such complexity leads to the production of not only olefins but also alkanes and aromatics through the more generally denoted methanol-to-hydrocarbons (MTH) process.

Although the MTH process displays intricated reaction networks, its chemistry has a well-known general classification grouped into stages. The first stage is the equilibration, where the partial dehydration of methanol produces a mixture of dimethyl ether (DME), water, and surface methoxy species (SMS), along with unconverted methanol and zeolite Brønsted acid sites (BASs).¹³ The next stage is the induction period, during which the first carbon–carbon bond formation occurs in a transient fashion. This stage also generates the first olefins, namely ethene and propene, as well as hydrocarbon pool (HPC) species.¹⁴ The mechanism underlying the first carbon–carbon bond has been a challenging and controversial research topic. Over the past decades, more than twenty direct reaction pathways have been proposed; however, most involve prohibitively high energy barriers, rendering them unlikely.¹⁵ Among these studies, the reaction between an SMS and carbon monoxide,¹⁶ a direct product of methanol decomposition, has gained significant acceptance, although alternative pathways based on the methane-formaldehyde mechanism also constitute the backbone of recent postulates.^17,18 Following the production of the first olefins, these compounds can undergo oligomerization and cyclization reactions to yield larger olefins and poly methylbenzenes (MB), which are part of the HCP species. The mechanistic details of this process are still elusive, nonetheless.^19,20 The subsequent stage produces the primary products, ethene and propene, in a steady fashion through a dual cycle in which HPC species serve as co-catalysts. The aromatic-based cycle predominantly yields ethene via poly MB compounds, whereas the alkene-based cycle primarily generates propene from larger olefins. The composition of the HPC and the relative contribution of each cycle are determined by the internal structure of the zeolite. The core reactions in these mechanisms are methylation and either cracking or dealkylation, depending on whether the active species are alkenes or poly MB. The mechanism of this stage is nevertheless more sophisticated and involves different types of elementary steps, including those connecting both cycles.²¹ In the final stage, HCP species further react and condense to form coke, which can accumulate and block the active sites of the zeolites, causing deactivation. Interestingly, reports suggest that water addition inhibits coke formation and prolongs the catalyst lifetime.²²

Studies of zeolite-catalyzed processes are abundant and have addressed specific aspects related to each stage of the MTH chemistry. For instance, researchers have analyzed the reactions responsible for the formation of deactivating species precursors during the induction period,^14–16 identified the roles of propene and methanol molecules in the steady formation of olefins via the HCP mechanism, and established rules to prevent catalyst deactivation caused by polyaromatic hydrocarbons that block or poison the zeolite cages.^23,24 From a theoretical perspective, research on the first stages of the MTH process has focused on aspects of framework methylation,^25,26 methanol dehydration to DME,^27,28 and primary olefin formation.^29,30 A principal feature of these studies is their reliance on static density functional theory (DFT) and metadynamics simulations. While these approaches have provided valuable mechanistic insights and accurate electronic structure data,^31,32 they exhibit limited predictive capability because they rely on the intuition of the chemical modeler. Other DFT limitations include the restricted ability to account for dynamic effects such as ion diffusion, framework flexibility, or entropy constraints.^33,34

A complementary tool to DFT and metadynamics methodologies that circumvents these limitations is reactive molecular dynamics (MD) simulations, using the ReaxFF reactive force field. Several studies have used this approach to determine the effects of metal dopants on the hydrothermal stability of zeolites,^35,36 the acidity depletion induced by water,³⁷ and zeolite dealumination.³⁴ Mechanistic aspects of the MTH process have also been addressed. Still, they are limited to high temperatures, reducing the sampling to entropically favored steps, or do not include in the training set relevant findings brought to light in the last decade.^33,38 This contribution presents the development of a new ReaxFF force field that describes the initiation reactions of MTH processes. It utilizes a prior parametrization for confined water in zeolite–clay composites³⁹ and, in addition, incorporates new data regarding C/O bonds, angles, and torsional energetics, as well as DFT profiles of the first steps of methanol conversion. The new force field examined various conditions of temperature and initial humidity, revealing important dynamic features that clarify the roles of cations and water in the initial reactions of methanol towards species such as DME and SMS.

2. Computational methodology

2.1. Force field

We used a ReaxFF force field to study the equilibration reactions of the MTH mechanism. This force field models reactive events by incorporating the bond order (BO) concept of Tersoff⁴⁰ and a polarizable charge description.⁴¹ The ReaxFF force field also has covalent and non-bonded interaction terms that facilitate the study of diverse materials and systems involving different phases.⁴² The reader can refer to the original publication of this force field for further details.^43,44

The C/Si/O/Al/Ca ReaxFF force field developed by Pitman and van Duin³⁹ was the starting point of the parametrization performed in this work. These authors used a combination of previously reported H/O/Si,⁴⁵ H/O/Al,⁴⁶ Ca/O, and Ca/H ReaxFF descriptions⁴⁷ with updated parameters trained against extensive DFT data containing H₂O clusters binding energies; H₂O, O₂, and H₂ bond dissociation energies; hydrogen transfer reactions involving H₂O, OH⁻/H₂O, and H₃O⁺/H₂O clusters; and ice crystals equations of state, among other H/O related information.

It is relevant to highlight that the H/O/Si⁴⁵ and H/O/Al⁴⁶ parameters have been validated and accurately describe dry and hydrated silica and aluminum clusters. More specifically, the reference force field³⁹ can reproduce the experimental details of clay interlayer hydration, silanol formation on exposed surfaces, the geometries of large-scale clay–zeolite composites, and water diffusion through the different sections of such composites. Besides, several studies have evaluated zeolite molecular models using this force field. For instance, Tian et al.⁴⁸ obtained a notable agreement between simulated and experimental structures of fused silica, Zheng et al.³⁸ confirmed the crystallinity of a H-ZSM-5 zeolite up to 2000 K using reactive MD simulations, and Chen et al.⁴⁹ evaluated ZSM-5, hydrated ZSM-5, and hydrated H-ZSM-5 molecular models that could reproduce relevant values of bond lengths and valence angles reported in the literature. While several of these investigations have focused on the conversion of various hydrocarbons, the reference force field³⁹ does not explicitly include any C/O parameterization. Thus, this crucial aspect requires additional attention.

Since interactions of hydrocarbons with the zeolite and relevant intermediates, such as water, are crucial in the MTH process, the ReaxFF C/O parameters were the object of additional training. To this end, the corresponding fitting, using the successive one-parameter search technique,⁵⁰ considered single, double, and triple bond dissociation energies, as well as valence angles and torsions of oxygenated hydrocarbon species. The new version of the force field also contains parameters for the C–O–Si and C–O–Al valence angles. This part of the training included profiles equivalent to approximately 700 energy points. The chemistry of the intermediate carbon monoxide was particularly challenging to describe, and its refinement required DFT data regarding dimerization reactions and its bonding with the zeolite framework. The C/O parametrization also used DFT results associated with the MTH mechanism. It covered methanol and DME adsorption on a BAS of the zeolite, the subsequent protonation reactions, the first carbon–carbon bond formation reaction between the SMS and carbon monoxide,¹⁶ surface acetyl species (SAS) adsorption, and SAS reactions with water, methanol, DME, and molecular hydrogen. The training set for these MTH reactions incorporated over 270 additional energy points and consisted of profiles that included reactants, products, saddle points, and reaction path structures.

All the DFT calculations of this contribution model T2 and T5 aluminosilicate clusters at the ωB97X-D/6-311+G(2df,2p)//ωB97X-D/6-311G(d,p) level of theory, which has shown satisfactory performance in the study of phenomena occurring in confined environments like those of zeolites.⁵¹ The parametric set for O/H, H/O/Si, H/O/Al, and Ca/O/H was unchanged and used as defined by Pitman and van Duin.³⁹

In addition to the validation reported in the literature, further testing shows that the force field enhanced in this work successfully reproduces lattice parameters of different zeolite topologies. For the FER, IFR, MFI, and TON frameworks, the maximum deviations in cell lengths are 2.36%, 2.36%, 1.91%, and 0.52%, respectively. Meanwhile, the discrepancies in cell angles remain consistent at 0.00% across all cases. This information is detailed in Table S1 of the SI. These errors fall within the accepted ranges for lattice parameter estimations of materials such as zeolites and perovskites using DFT and force field methods.^52,53 The role of water in zeolite-catalyzed processes is significant since humidity can modify the acidity of a porous environment.²² Therefore, we also evaluated the energetics of water-assisted proton transport in the first coordination sphere of the aluminum atom based on the system studied by Ryder et al.⁵⁴ and Joshi et al.³⁷ Fig. S1 of the SI shows the structures of the corresponding reactant, transition state, and product, and Fig. S2 displays the potential energy diagram. Overall, the force field estimates an energy barrier of 7.63 kcal mol⁻¹, which aligns well with the DFT value of 7.28 kcal mol⁻¹ computed in this work.

2.2. Molecular model

The structure used in the simulations has an MFI topology, as specified by the International Zeolite Association (IZA).⁵⁵ The initial catalytic system consisted of a periodic 2 × 2 × 2 supercell, whose lattice constants were a = 40.2 Å, b = 39.8 Å, and c = 26.8 Å, with aluminum atoms in T12 position replacing the original tetrahedral silicon atoms. The density of the aluminum atoms was four per unit cell or, equivalently, of a Si/Al ratio of 23, which is below the limit of 27 determined experimentally for the associated H-ZSM-5 zeolite.⁵⁶ Several theoretical studies have considered different active tetrahedral positions, but T12 is representative because it lies in the intersection of the sinusoidal and straight 10-ring channels, a kinetically favored position that also permits the interaction of molecules with considerable size.¹³ Balancing protons were initially placed at the O26 position. Still, these protons constantly hop around the first coordination sphere of the tetrahedral aluminum atom before reacting under the conditions of the simulations used in this work. The modifications applied to the initial MFI unit cell, necessary to build the initial catalytic system, were done using the VESTA software.⁵⁷

The total number of methanol molecules added to the system was 32, the same quantity of BASs in the zeolite supercell. All the simulations included this acidity and the amount of initial precursor. To analyze the effect of moisture in the initial stages of the MTH mechanism, 8, 16, and 32 water molecules were added to the system, representing substoichiometric and stoichiometric quantities relative to the BASs in the first two cases and the last case, respectively. Experimental studies and theoretical calculations indicate that water concentration in zeolite can be significantly higher.^58,59 However, these conditions could damage the zeolite structure through dealumination and negatively impact pertinent reactivity in the studied systems.⁶⁰ Moreover, packing both gaseous and solid phases using Packmol software⁶¹ presents convergence problems at high water loads.

2.3. Molecular dynamics simulation details

This work ran duplicates of three sets of simulations. The first set evaluated the effect of the temperatures 600, 800, 1000, and 1200 K on the initial reactions of the MTH process, starting from a system consisting of an acidic zeolite with stoichiometric quantities of methanol. The second set of simulations examined the same reactions but considered the effect of water loads, as outlined in Section 2.2, rather than temperature, which remained constant at 800 K for all initial configurations. The third set of simulations was similar to the second set, but, in this case, the temperature was 1000 K. The 600–800 K temperature range is common in methanol conversion experiments on zeolites of industrial significance, such as H-SAPO-34 and H-ZSM-5.^62,63 Higher temperatures are not of practical interest, but they are instrumental to obtaining theoretical insight into initiation reactions of methanol conversion that could be transferable to low temperatures.^64,65 Note that 1000 and 1200 K are significantly lower than conditions implemented in similar reactive MD studies.³⁸

The platform used in this work was LAMMPS.⁶⁶ All the simulations started with an energy minimization by adjusting the positions of atoms to eliminate excessive forces. This minimization or relaxation avoids unreasonable behavior in the initial simulation steps and was completed when either the energy or force tolerance criteria, that is, 1 × 10⁻⁶ and 1 × 10⁻⁶ kcal mol⁻¹ Å⁻¹, respectively, were satisfied. In the next stage, the temperature was set to 100 K via a Gaussian distribution of the atoms' velocities and zeroing the angular momentum of the simulation box, a condition under which a two-step equilibration was carried out. In the first equilibration stage, the system ran for 100 ps in the NVT ensemble, and in the second one, the trajectories were integrated for another 100 ps in the NPT ensemble at 1.0 atm. A heating rate of 10 K s⁻¹, also implementing the NPT ensemble, increased the temperature from 100 K to 600, 800, 1000, or 1200 K, depending on the conditions to evaluate. After setting the final temperature, the isobaric-isothermal production of the MD simulations started and lasted 1000 ps in all cases. This time was sufficient to distinguish relevant dynamic features of the initial reactions of the MTH process.

The time step of every simulation stage was 0.1 fs, which is sufficiently small to ensure energy conservation and correctly describe the vibrational modes associated with different reaction coordinates. The damping parameters for the thermostat and barostat were 50 ps and 100 ps, respectively, and both ensembles, NVT and NPT, utilized the Nosé–Hoover^67,68 time integration to calculate the atoms’ velocities and positions at each time step. The volume of the simulation box in the NPT ensemble varies, which is ideal to account for the framework flexibility and its influence on the mechanism of methanol conversion.^34,69 In fact, analysis of the trajectories with 8 initial water molecules at 800 K revealed volume changes between +1.9 and −4.7%, closely aligning with the +1.6 and −3.8% changes reported by Wispelaere et al.⁶⁹ These authors implemented a different methodology by modeling an H-SAPO-34 zeolite with defined quantities of water and methanol molecules at 623 K and using ab initio molecular dynamics at the revPBD-D3/DZVP-GTH level of theory with pseudopotentials. Nonetheless, the marked agreement in the volume variation supports the validity of the ensemble used here to account for the framework flexibility dynamically. Finally, a bond order cutoff of 0.4 was the criterion of atom connectivity, as established in a previous work considering a similar zeolite system.³⁴ This bond order cutoff only defines the molecule recognition – this does not affect the ReaxFF energies or forces.

3. Results and discussion

3.1. Parametrization

Examples of species participating in the initiation reactions of the MTH mechanism include the framework BASs and oxygenated organic species adsorbed or diffusing through the zeolite pores.⁷⁰ More specifically, studies have experimentally identified or inferred the existence of DME, methane, carbon monoxide, formaldehyde, and SMS, among others.²¹ Given their relevance at the initial stages of the MTH process, the parametrization of the ReaxFF force field performed in this work considered these molecules, as shown in Fig. 1 for several cases. The panels (a) and (b) of this figure present the C–O single and triple bond length energies, respectively, calculated by DFT (blue lines) and molecular mechanics (MM) using the trained ReaxFF force field (red lines), while panels (c) and (d) display, in each case, the H–O–C and C–C–O valence angle profiles. Methanol, carbon monoxide, and ethanol were model structures for these calculations. ReaxFF well reproduces the DFT profiles corresponding to the C–O single bond length and C–C–O valence angle energies. Even as the evaluated lengths and angles become small and large, the force field and DFT estimations remain in close agreement. Carbon monoxide can be present at the initial stages of the MTH process as a product of the initial decomposition of methanol. Its ReaxFF C–O triple bond length energy curve well represents the DFT data, with some deviation in regions away from the minimum, a commonly observed trend in force field training.^64,71 Despite some slight deviations, the H–O–C valence angle energies also exhibit a good agreement between DFT and ReaxFF values.

Fig. 1 Comparison between distortion energies computed via density functional theory (DFT) at the ωB97X-D/6-311+G(2df,2p)//ωB97X-D/6-311G(d,p) level of theory and the ReaxFF force field developed in this work for (a) C–O single bond, (b) C≡O triple bond, (c) H–C–O angle, and (d) C–C–O angle. White, dark gray, and red spheres represent hydrogen, carbon, and oxygen atoms, respectively.

The panels (a) and (b) in Fig. 2 show the potential energy diagrams (PEDs) of the initiation reactions of the MTH process and an alternative path starting from DME, respectively. Capital letters denote stationary points, as explicitly defined in the caption. The first stage of the MTH mechanism consists of a well-defined equilibration, where methanol (B) dehydrates to form DME (K) and water. The first barrier of Fig. 2(a) portrays the first step of this process, which proceeds through the protonation of methanol by a zeolite's BAS. Eventually, the protonated molecule releases water, and the remanent methyl group can form an SMS (F), an unequivocally identified intermediate that plays a significant role in the first carbon–carbon coupling reactions.^72,73 The reverse reaction across the first barrier of Fig. 2(b) illustrates how an SMS and additional methanol associate to produce DME and recover the original zeolite's BAS (A). All previously described steps are reversible, indicating the coexistence of methanol, DME, SMS, water, and BASs during the equilibration. Letters C, L, D, and M signify the methanol and DME protonation saddle points in the first two cases, which are submerged regarding the reactants, and the subsequent saddle structures leading to the SMSs from the corresponding protonated molecules in the last two cases.

Fig. 2 Comparison between the potential energy diagrams (PEDs) computed via density functional theory (DFT) at the ωB97X-D/6-311+G(2df,2p)//ωB97X-D/6-311G(d,p) level of theory and the ReaxFF force field developed in this work for the initiation reaction of the (a) methanol-to-hydrocarbon (MTH) process and (b) an alternative path starting from dimethyl ether (DME). The PEDs display reactions concerning protonation, zeolite methylation or surface methoxy species (SMS) formation, first carbon–carbon coupling, and surface acetyl species (SAS) reaction involving methanol and DME. White, dark gray, red, brown, and dark pink spheres represent hydrogen, carbon, oxygen, silicon, and aluminum atoms, respectively. Each capital letter detonates a stationary point: A is the zeolite Brønsted acid site (BAS), B is methanol, C is the intermediate complex of methanol protonation, D is the intermediate complex of zeolite surface methylation from a methanol molecule, E is carbon monoxide, F is SMS, G is the intermediate complex of the first carbon–carbon coupling, H is the SAS, I is the intermediate complex of the reaction between SAS and methanol, J is ethyl acetate, K is DME, L is the activated complex of DME protonation, M is the activated complex of zeolite methylations from a protonated DME species (or DME formation from methanol and a SMS), and N is the intermediate complex of the reaction between a SAS and DME.

C and L ReaxFF energies agree with the DFT values, and those of D and M qualitatively reproduce the theoretical trends with some underestimations. The discrepancies in the SMS (F) formation are an artifact introduced to the force field, which, in practical terms, emulates the application of a biased potential used to accelerate the dynamic trajectories and obtain reactivity within feasible simulation times.⁷⁴ The systems showed negligible conversion without the lowering of this energy barrier. Besides, because SMS is the very first reactive intermediate and there are no competing products at this stage of the MTH process,^21,26 it is unlikely that such a modification will cause mechanistic misinterpretations or distort the reaction dynamics. This fact is especially valid as we concentrate our study on analyzing the dynamic aspects of the equilibration stage alone. The option of raising the temperature to speed up the simulations is not suitable due to the potential significant alteration of the free energy profiles, which could affect the likelihood of meaningless reaction pathways.⁷⁵

Experimental and theoretical studies of the MTH process have demonstrated that SMS participate in the steps concerning the first carbon–carbon bond formation.⁷⁶ Nonetheless, details of the reactions involved have been elusive and a subject of debate for years, and, as a result, relevant contributions have proposed mechanisms like carbene insertion,⁷⁷ methane-formaldehyde,⁷⁸ methoxymethyl cation,¹⁸ and many others. In an endeavor aimed at clarifying the feasibility of the mechanisms describing the direct carbon–carbon coupling, Lesthaeghe et al.¹⁵ found, using DFT, that most of them were unrealistic by virtue of their abnormally high energy barriers. More recently, Liu et al.¹⁶ proved that carbon monoxide is a decomposition product of methanol over inert zeolite surfaces that could react with SMS to form the first carbon–carbon bond. The authors used DFT to estimate an energy barrier of 19.1 kcal mol⁻¹, a value 60% lower than that of previously proposed mechanisms. This reaction was included in the force field training, as shown in Fig. 2.

It is necessary to highlight that the methane-formaldehyde mechanism and its variations have been popular to explain the first carbon–carbon bond formation and the presence of methane in the first stages of MTH processes;⁵¹ however, the free energy barriers of these processes are high, rendering them unfeasible in comparison with the mentioned carbonylation of SMS proposed by Liu et al.¹⁶ Regardless of the diverse studies covering the topic, the picture of the first carbon–carbon bond formation is still under debate, and several authors have recently proposed variations of mechanisms based on the methane-formaldehyde scheme that involve more favorable free energy landscapes.¹⁸ Under this scenario, the selection of the reaction or steps leading to the first carbon–carbon bond to train the force field is challenging. Still, the mechanism proposed by Liu et al.¹⁶ is selected here because it has solid experimental evidence and the lowest energy barrier.

While the force field reproduced the energy profile of this reaction, our simulations indicated that the free energy of this type of process remains too high^31,32 to be observed within reasonable simulation times. More specifically, systems consisting of methylated zeolite and different concentrations of carbon monoxide were completely inert even after 5 ns of simulation. Significant reactivity was observed only after implementing accelerated dynamics using the tracking method⁷⁹ (data not shown). These results imply that rigorous reactive MD studies of reactions involving carbon–carbon bond formation, or the production of methane and formaldehyde, require efforts that deserve a separate study, and the force field presented in this work constitutes a good starting point for those endeavors.

The product of the SMS + CO reaction is an acetyl group that can bond the zeolite to form a SAS and later react with water, remaining methanol, DME, or any other molecules present at this stage. The rightmost part of the diagrams in Fig. 2 illustrates the energy barriers for the reaction of SAS with methanol and DME in panels (a) and (b), respectively. As observed, the ReaxFF force field can also qualitatively reproduce DFT profiles for these cases. The induction period of the MTH process covers all the reactions from the equilibration stage to the formation of the first olefins, ethane and propene, and the HCP.¹⁴ The associated induction time is of the order of minutes, and the theoretical study of this stage represents a particular challenge for atomistic simulations because of the timescale and complexities of the transformations involved.⁸⁰

3.2. Temperature effect on species profiles

Fig. 3 and 4 show the species concentration profiles of the most relevant neutral and charged species, respectively, identified in the initiation mechanism of the MTH process as a function of the temperature. Starting with the profiles presented in Fig. 3(a) for the feedstock methanol, it is evident that the temperature accelerates its consumption since simulations performed at 600, 800, 1000, and 1200 K report a total of 30, 25, 19, and 14 molecules, respectively, after 1000 ps of trajectory. The low reactivity observed at 600 K is understandable for the reactive MD trajectories carried out at the short timescale of this study, given that the first stages of the MTH process occur in the order of minutes;⁸¹ however, it is still notable that we even observe limited reactivity, without the necessity to enforce excessively high temperatures.³⁸ The chemistry of methanol decomposition displayed by the simulations goes beyond the profiles of the force field training presented in Section 3.1, showing SMS, protonated methanol, and methyl cations as the most prominent products. Of these intermediates, SMS and the methyl cations are reactive and most likely represent incipient precursors of the HCP and short olefins.⁸²

Fig. 3 Number profiles of (a) methanol, (b) water, (c) surface methoxy species (SMS), and (d) methane as a function of time of simulations conducted at temperatures in the 600–1200 K range.

Fig. 4 Number profiles of (a) hydronium ions and (b) methyl cation species as a function of time of simulations conducted at temperatures in the 600–1200 K range.

Methanol protonation is an essential step in the MTH initiation mechanism. Its product reacts fast and is present in the trajectories only in low quantities. The simulations showed a steady maximum number of protonated methanol molecules at different conditions, three at 600 K and two at the remaining temperatures (data not shown). Conversely, their approximated lifetimes displayed considerable differences. At 600 K, the life span lay in the 40–120 ps range and decreased to the 5–20 ps range at 800 K. As observed, the temperature favors the reactivity, and at 1000 and 1200 K, protonated methanol remained within the zeolite pores for less than one picosecond on average.

The water production observed in Fig. 3(b) is proportional to methanol consumption, as the temperature also promotes it. After 1000 ps, the number of produced water molecules at 600, 800, 1000, and 1200 K is 1, 2, 12, and 16, respectively. However, at intermediate steps, the simulations at 600 and 800 K report a corresponding maximum of 2 and 4. The simulations also show that the water, a product of protonated methanol decomposition, is not a simple spectator inside the zeolite since this species interacts with the BASs to also get protonated and diffuse as a hydronium ion, whose profile is illustrated in Fig. 4(a). The diffusion of this ion is not a new phenomenon and has been documented in studies related to proton mobility.⁸³

Fig. 3(c) demonstrates that increasing temperature in the 600–1000 K range enhances SMS formation, as the number of such SMS increases from 1 to 9 by the end of the simulations. Nonetheless, the SMS number decreases to 5 at 1200 K due to a competing pathway that produces methane from the methyl cation precursor. This trend is further supported by Fig. 3(d), which shows a substantial increase in methane production at 1200 K relative to lower temperatures. Specifically, the number of methane molecules reaches 9 at 1200 K, decreases to 3 at 1000 K, and is negligible at lower temperatures.

Although the trends of methanol consumption and water production correlate well, their relationship with methyl cation profiles is more complex and can be roughly linked to the temperature. As shown in Fig. 4(b), elevated temperatures favor the diffusion of this cation through the zeolite pores, but its total number remains relatively low during the trajectories; for example, they are a maximum of 2, 4, 4, and 5 species at 600, 800, 1000, and 1200 K, correspondingly. These results indicate that the diffusion of cations like H₃O⁺, CH₃⁺, and protonated methanol through the zeolite pores is recurrent in the initiation mechanism of the MTH process. Although notable and pertinent static DFT studies do not explicitly discuss this phenomenon,^31,32 it is not an entirely new topic. In fact, ion exchange was the first property systematically studied in zeolites, and its origins date to the 1850s.⁸⁴

The following sections discuss aspects observed during the trajectories relevant to the initiation reactions of the MTO process: kinetics of methanol consumption and mechanistic aspects of SMS, methane, and DME formation.

3.2.1. Methanol consumption kinetics. To semiquantitatively assess the reliability of the methanol decomposition behavior, the corresponding apparent activation energy was computed by considering a first-order reaction and plotting the natural logarithm of the rate constant against the inverse of the temperature, as presented in Fig. 5. The obtained value of 5.30 kcal mol⁻¹ is in reasonable agreement with that estimated by Pérez-Uriarte et al.⁸⁵ of 8.08 kcal mol⁻¹. The authors studied the effect of replacing methanol with DME as the initial reactant to produce olefins over an H-ZSM-5 zeolite embedded in a reactor free of heat and mass transfer limitations and considering a 598–673 K temperature range. The kinetic mechanism of the study followed a lumped methodology where methanol mainly produced C₂–C₄ olefins, and whose parameters originated from an empirical fitting with experimental profiles. The apparent activation energy shown in Fig. 5 also agrees well with the value of 6.28 kcal mol⁻¹ reported by Kaarsholm et al.,⁸⁶ who proposed a kinetic model based on the HCP mechanism and fitted its parameters to results obtained in a fluidized bed reactor, using a phosphorous-modified H-ZSM-5 catalysis. In this case, methanol theoretically produced high-order hydrocarbons resembling the HCP. The experiments covered temperatures between 673 and 823 K and a co-feed with hexene simulated recirculation conditions. Curiously, Chen et al.⁸⁷ calculated apparent activation energies of around 9.08 and 6.45 kcal mol⁻¹ to produce methane and C₂–C₃ olefins, respectively, from the methanol–DME equilibrium by also fitting a lumped mechanism with experimental results obtained in an oscillating microbalance reactor, but with a SAPO-34 catalyst instead. Overall, the apparent activation energy of Fig. 5 is consistent with several experimental-based results of studies that used zeolites as a catalyst and methanol as a source of short olefins or intermediate large hydrocarbons. However, other models may provide different results depending on the research focus, assumptions, and methodologies.⁸⁸

Fig. 5 Arrhenius plot for methanol conversion. The apparent activation energy (E_a) was calculated using first-order kinetics and a linear fit between the data of the natural logarithm of the rate constant and the inverse of the temperature. The equation is y = −2666.92·x + 22.62, and the associated correlation coefficient is 0.98.

3.2.2. Mechanistic aspects of the initiation steps: water, SMS, methane, and DME. Water concentration is proportional to methanol consumption because water originates from two reaction steps consisting of (i) methanol protonation and (ii) dissociation of this protonated methanol molecule to produce water and a methyl cation. Later, this methyl cation can (i) form SMS, (ii) form methane, or (iii) diffuse through the zeolite pores. Fig. 6 shows the sequence of methanol protonation from panels (a) to (c), how a diffusing protonated methanol molecule forms water and a methyl cation in panels (d) and (e), and the formation of the relevant intermediate SMS from a diffusing methyl cation in panels (f) and (e).

Fig. 6 Snapshots of the first steps of the methanol-to-hydrocarbon (MTH) mechanism, as shown by the reactive molecular dynamics (MD) trajectories. (a) Methanol approaches a Brønsted acid site (BAS) of the zeolite to (b) form an intermediate complex that leads to a (c) protonated methanol cation and a negatively charged framework. The charged methanol species (d) diffuses through the zeolite channels until (e) it dissociates in water and a methyl cation near a negatively charged oxygen belonging to the zeolite lattice. This proximity facilitates the (f) encounter between the methyl cation and bridging Al–(O⁻)–Si group in the framework to (g) form a surface methoxy species (SMS). White, dark gray, red, brown, and dark pink spheres represent hydrogen, carbon, oxygen, silicon, and aluminum atoms, respectively.

Mechanistic studies based on DFT calculations have established a stepwise route for SMS formation. In the initial step, methanol approaches the BAS of the zeolite and gets adsorbed, eventually producing protonated methanol. Later, the protonated methanol rotates around the initial BAS region until the methyl group faces the oxygens of the zeolite framework to form the bond of the SMS and release water.⁸⁹ The critical difference between DFT calculations of the references and the reactive MD simulations of this work is that the protonated methanol does not rotate around the initial acid center of the zeolite but diffuses through the pores for periods that depend on the temperature, as explained previously, before forming the SMS and water, as presented in Fig. 6. Although the rotational energy barrier is low, approximately 4.6 kcal mol⁻¹,²⁸ entropy may favor diffusion to a certain degree.⁹⁰ The diffusion of cations, including the hydronium ion, is an exhibition of dynamic effects in zeolite-based processes.

Methane formation is dominant in the initial stages of methanol conversion over acidic zeolites and sustained throughout the catalysis lifetime.^91,92 Liu et al.⁹³ presented concrete numbers. In experiments conducted using an H-ZSM-5 zeolite at 748 K, the authors measured selectivities of 50% for methane, 25% for formaldehyde, and 25% for CO and CO₂ species at a methanol conversion of 0.24%. The methane formation mechanism is more debated, nonetheless. Performing deuterium tracing experiments in a SAPO-34 at 673 K, Wei et al.¹⁷ concluded that methane forms at the beginning of the MTH process, mainly following the methane-formaldehyde mechanism.^51,78 In the same work, the authors discarded the possibility of methane formation from the hydrogens of the framework since experiments with deuterated zeolite, Z–OD, produced non-deuterated methanol, CH₄. In this work, methane formation followed such a forbidden pathway, that is, methane formed through the reaction Z–OH + CH₃⁺ → Z–O⁺ + CH₄, where the methyl cation is a diffusing species and the hydride source is the hydrogen attached to the zeolite.

Hydride transfer reactions in zeolite catalysts occur during the conversion of alkanes to fine chemicals, particularly through the carbonium ion cracking cycle.⁹⁴ Studies indicate that low Si/Al ratios and confinement effects significantly increase the rate of hydride-dependent chemistry.^95,96 Similarly, strong acidity, confinement, and elevated temperatures promote the transport of charged species within zeolite pores. This effect is demonstrated by proton mobilization across both narrow and wide regions of the zeolite framework, as well as by hydrogen–deuterium (H/D) isotope exchange between BASs of the zeolite and methane.⁹⁷ Because the system under study involves considerable acidity, molecular-sized pores and channels, and high temperatures, it is possible to hypothesize that these factors likely create favorable conditions for the methane formation pathway observed in the trajectories, namely, Z–OH + CH₃⁺ → Z–O⁺ + CH₄. However, this observation is inherently speculative and cautions against conducting simulations of zeolite-catalyzed processes at elevated temperatures to accelerate the dynamics. Further research at temperatures above 1000 K is necessary to confirm this methane formation pathway. Consequently, the following sections on the effect of water present only results obtained at 800 K and 1000 K.

Mechanistic studies have established well-defined stages for the overall conversion in the MTH process. The first stage is the equilibration among methanol, DME, water, and SMS molecules, where DME is the gaseous dehydration product of methanol.¹³ The simulations of this work showed little DME, with a maximum number of two for this species between 800 and 1200 K (data not shown). The simulations at 600 K did not display any DME molecule.

Theoretical studies since the 1990s have defined two possible pathways for DME formation in the equilibration stage: one dissociative (or stepwise) and the other associative (or concerted).⁸⁹ The dissociative pathway consists of the formation of DME through an SMS, which is a product of methanol protonation and dehydration. Then, a second methanol molecule approaches the SMS to react and form a protonated DME, which neutralizes by giving the proton to the framework. The associative pathway involves the simultaneous adsorption of two methanol molecules in a BAS of the zeolite with the rearrangement necessary to facilitate interactions between intermolecular methyl and hydroxyl groups aimed at forming a DME species.⁹⁸ Our simulation showed a variation of the dissociative pathway for the DME formation, consisting of a reaction between methanol and a diffusing and short-lived methyl cation initially belonging to a protonated methanol molecule or an SMS, as sequentially shown from panel (a) to (e) of Fig. 7. The transition state of the associative pathway has a significant stabilization enthalpy. Still, its entropy penalty is also considerable, which must favor the dissociative route of DME formation or its variation at the high temperatures implemented in the simulations of this work, that is, between 800 and 1200 K.^27,99

Fig. 7 Snapshots of dimethyl ether (DME) formation through a variation of the dissociative pathway involving a methyl intermediate,^26,27 as shown by the reactive molecular dynamics (MD) trajectories. (a) A methanol molecule in the gas phase encounters a diffusing methyl cation to (b) form a protonated DME species. This cation (c) finds a negatively charged oxygen in the framework to (d) get adsorbed and (e) transfer a proton to recover a Brønsted acid site (BAS) of the zeolite by releasing a neutral DME molecule. White, dark gray, red, brown, and dark pink spheres represent hydrogen, carbon, oxygen, silicon, and aluminum atoms, respectively.

A relevant aspect of the reactive MD results discussed in the previous paragraphs involves the diffusion of short-lived cations. This behavior contrasts with static DFT calculations, where prevailing assumptions suggest that these cations react quickly at the nearest available reactive site. This discrepancy is likely due to the inherent limitations of the static methodology.³³ Cation diffusion within the zeolite pores is a well-documented phenomenon that determines critical aspects of the zeolite chemistry, such as acidity via water-mediated proton transport,⁸³ cation exchange applications like detergency and nuclear water removal, and reduction of NO_X pollutants mediated by NH₃.¹⁰⁰

3.3. DFT analysis on the effect of temperature on initiation reactions of the MTH process

Section 3.1 explained the difficulties in observing the reaction that forms the first carbon–carbon bond within feasible simulation times (less than 5 ns). A standard way to accelerate reactivity is to increase the temperature, but results in Section 3.2 show that this approach can promote pathways not reported in the literature, such as methane formation via the reaction Z–OH + CH₃⁺ → Z–O⁺ + CH₄. To shed light on the complications that arise from setting high temperatures in reactive MD simulations of zeolite systems, static DFT calculations were performed for reactions (R1)–(R3).

Z–OCH₃ + CO ↔ Z–OC(=O)CH₃ (R1)

Z–OCH₃ + CH₃OH ↔ Z–OH + CH₄ + CH₂O (R2)

Z–OH + CH₃OH ↔ Z–OCH₃ + H₂O (R3)

Reaction (R1) involves the first carbon–carbon bond formation between SMS and carbon monoxide to form the SAS product. Reaction (R2) represents a competing pathway that becomes significant when the methanol concentration is high, such as during the equilibration stage of the MTH process. This alternative route leads to the production of undesired formaldehyde and methane. The final reaction (R3) involves the reversible generation of SMS from methanol and the zeolite BAS. These calculations employed a 44T finite cluster model, placing the aluminum atom at the T12 position and bonding the balancing proton to the O26 oxygen atom. The reactive MD trajectories simulate a structure with this configuration, as detailed in Section 2.2. Hydrogen atoms are bonded to the outermost silicon atoms at fixed positions to maintain the tetrahedral orientation and prevent unrealistic deformations. Geometry optimization and vibrational frequency calculations employed the ωB97X-D/6-311+G(2df,2p)//ωB97X-D/6-311G(d,p) level of theory,⁵¹ consistent with the methodology used to generate the force field training set. Intrinsic reaction coordinate (IRC) calculations verified the connectivity between each saddle point, reactants, and products. All vibrational frequencies were treated within the harmonic approximation, and the reference state for non-zeolite species was the gas phase. All electronic structure calculations were performed using the Gaussian16 software package.¹⁰¹

Fig. 8(a) presents the free energy barriers (ΔG^≠) for the competing reactions (R1) and (R2), as well as the enthalpy barrier (ΔH^≠) for reaction (R1), and Fig. 8(b) shows the free energy and enthalpy barriers for both the forward and reverse steps of reaction (R3). The results obtained in this study align closely with those reported by Plessow and Studt,³² who employed a periodic model of the MFI framework using the PBE-D3 method and the projector-augmented-wave (PAW) approach. For example, these authors reported free energy barriers at 673 K for reactions (R1) and (R3) of 42.0 and 40.0 kcal mol⁻¹, respectively. Following the methodology presented here, the corresponding values are 47.0 and 40.0 kcal mol⁻¹.

Fig. 8 (a) Gibbs free energy barrier (ΔG^≠) as a function of the temperature of the competing reactions (R1) and (R2), including the enthalpy barrier (ΔH^≠) of (R1). These two properties help visualize the entropy contribution (−T·ΔS^≠) in this reaction. (b) Gibbs free energy and enthalpy barriers of the forward and reverse steps of reaction (R3) as a function of the temperature. An increase in temperature narrows the free energy barrier of the forward and reverse steps, reaching the turnover point at approximately 1300 K. Above this temperature, the reverse step of (R3) is more prominent. The reactions (R1), (R2), and (R3) are Z–OCH₃ + CO ↔ Z–OC(=O)CH₃, Z–OCH₃ + CH₃OH ↔ ZOH + CH₄ + CH₂O, and Z–OH + CH₃OH ↔ Z–OCH₃ + H₂O, respectively.

Fig. 8(a) shows that the free energy barriers of reaction (R1) are lower than those of reaction (R2) in the 250–2000 K range, with the difference increasing at higher temperatures. This trend is positive because the formation of the first carbon–carbon bond prevails over the production of methane, apparently pushing the MTH process forward. Although higher temperatures favor reaction (R1), the associated increase in the free energy barrier results in reduced overall reactivity. This trend is a consequence of the dominance of the entropy term on the free energy barrier, as indicated by the gap between the solid and dashed blue lines in Fig. 8(a) that quantifies the term −T·ΔS^≠. Fig. 8(b) illustrates another example of the adverse effect of elevated temperatures on reactivity. The reverse conversion of the BAS and methanol to water and SMS, the most relevant methylating agent in the MTH process, dominates above 1300 K. The depletion of SMS towards the initial reactants would occur despite the relation of the enthalpies of reaction (dashed lines) remaining almost unchanged. The turnover of the free energy barriers at 1300 K is small. However, it fundamentally indicates that reactions over zeolites could exhibit this trend and, in some instances, have a significant impact on the overall system, leading to speculative pathways at elevated temperatures. The studies by Gounder and Iglesia^75,102 provide in-depth insights into the roles of entropy and temperature in zeolite catalysts.

The free energy profiles shown in Fig. 8(a) clarify the observed lack of reactivity between SMS and carbon monoxide (R1) or methanol (R2) in the reactive MD simulations, as described in Section 3.1. At 800 K, the free energy barriers for reactions (R1) and (R2) are 50.7 and 56.5 kcal mol⁻¹, respectively, which are too high to be overcome within the practical simulation times (less than 5 ns). This result is unsurprising, as the induction period, including reactions such as (R1) and (R2), occurs over minutes, a timescale inaccessible to the trajectories and methodology implemented here. From this perspective, the simulations behave well. To address this limitation, it is necessary to develop techniques that accelerate rare events, thereby enabling exploration of a broader range of the MTH reaction space. For instance, such methods would facilitate the formation of primary products such as ethene and propene from methanol and the zeolite BASs. Achieving this objective is particularly challenging because the entropy constraints are significant in zeolite systems.

3.4. Effect of water loads at 800 K

Water can affect zeolite-based processes through different mechanisms. Literature has reported that water can compete for adsorption on the acid centers of the zeolites, either BAS or Lewis acid sites (LASs), inhibiting reactions like ethene oligomerization. This effect is relevant since such a reaction is associated with coke formation and the subsequent catalyst deactivation.²⁵ Water also participates in the solvation and stabilization of reactants or transition states, hindering or promoting a given process, respectively, depending, among other factors, on the nature of the involved species.¹⁰³ As presented in Section 3.2, water can modify the nature of the BASs of a zeolite from static to dynamic because of the formation of hydronium ions and protonated water clusters.¹⁰⁴ Due to water's rich and diverse role in zeolite-catalyzed reactions, it is pertinent to study the effect of humidity on the initial stage of the MTH process.

Fig. 9 and 10 present the species profiles of the most relevant stable and charged species, respectively, identified in the initial reactions of the MTH process at different initial water loads. The temperature of the simulations was 800 K, and the initial number of methanol molecules was 32, while those of water were 0 (red), 8 (green), 16 (blue), and 32 (brown). Fig. 9(a) shows that starting with 8 water molecules slightly promotes the methanol consumption during most of the simulated time. The most noticeable discrepancy with the reference system (no initial water) is between approximately 800 and 900 ps, where water facilitates the consumption of up to 8 additional methanol molecules. However, by the end of the trajectories, both simulations exhibit a comparable alcohol concentration. This variation is likely a result of the time and size scale limitations inherent in this type of reactive MD simulations. At higher water loads, that is, 16 and 32, promoting effects also manifest and affect the methanol profiles through the simulations. The systems with an initial water load of 16 and 32 increase methanol consumption up to 10% and 17%, mainly after 400 and 300 ps, respectively. In any case, these systems also end up with a similar number of methanol molecules as the reference configuration.

Fig. 9 Number profiles of (a) methanol, (b) water, (c) surface methoxy species (SMS), and (d) methane as a function of time of simulations conducted at 800 K and the initial water loads of 0, 8, 16, and 32.

Fig. 10 Number profiles of the (a) hydronium ion and (b) methyl cation species as a function of time of simulations conducted at 800 K and the initial water loads of 0, 8, 16, and 32. Panel (b) displays smoothed profiles due to the low values of the data (between 0 and 6) and their rapid variation through the trajectories. The smoothing consisted of averages in windows of 4 ps.

Fig. 9(b), corresponding to the number of water molecules, shows numbers that differ from the loads specified at the beginning of the production period. This behavior results from the formation of hydronium ions during heating after the equilibration. The water loads specified in the legend result from the summation of species of Fig. 9(b) and 10(a) at zero time. Despite being initially depleted due to the formation of hydronium ions, the overall trends of Fig. 9(b) show an increase in water concentration for all the scenarios studied. For the reference simulation and the initial water loads of 8, 16, and 32, the dynamics present an increase of 3, 4, 3, and 7 in the number of water molecules, respectively, after 1000 ps of trajectory. This behavior is a consequence of methanol dehydration. In any case, water never reaches its initial concentrations, and its final numbers fall by 1, 5, and 8 for the initial water loads of 8, 16, and 32, respectively.

The profiles illustrated in Fig. 9(c) indicate that the initial presence of water encourages SMS production since the lines corresponding to initial 8, 16, and 32 water loads are almost always above that of the reference after 200 ps. In particular, an initial water load of 8 significantly favors the formation of SMS after 650 ps of simulation, ending with around 4 (80%) more SMS than the reference at 1000 ps. For the initial water load of 16, significant promoting effects arise between 600 and 950 ps. High water concentrations accelerate the formation of SMS because the simulation with an initial water load of 32 significantly increases the SMS formation after 300 ps; however, the overall trend is intricated and mixed with profiles starting with a lower number of water molecules.

The methane profile in Fig. 9(d) manifests a negligible effect of initial water loads, showing only one molecule after 500 ps. Methane formation from methyl cations, whose quantities are equally low due to the enhanced production of the competing SMS, rationalizes the low number of methane molecules detected in the simulations. The presence of extra water molecules and the formation of hydronium ions also explain the low methane concentration, as this process depletes hydrogen from the framework, which, together with methyl cations, constitutes the source of methane production, as explained in Section 3.2.2.

Fig. 10(a) presents the concentration profiles of hydronium ions. The results evidence a marked increase in this species compared to the reference system without initial water molecules and corresponding simulations evaluating the temperature effect shown in Fig. 4(a), where the maximum number of hydronium ions detected was 4 at 1000 K. In Fig. 10(a), the reference system and those with initial water loads of 8 and 16 increase the number of hydronium ions through the trajectories by 5, 6, and 3, respectively, to end up with 10, 12, and 14 molecules of this cation in each case. On the other hand, the number of this species oscillates and remains, on average, the same, around a value of 15 for the initial water load of 32. The most relevant implication of these results is the change in the zeolite acidity. Although these phenomena have been reported previously in the context of water and zeolite BASs interactions from theoretical and experimental perspectives,^37,105 water protonation in the context of an actual process like the MTH studied in this work is scarce. It is relevant to highlight that altering the acid site nature during methanol conversion processes may be restricted to specific initial water concentrations.

The presence of water has minimal impact on the profiles of methyl cations, as depicted in Fig. 10(b), likely because of its reactivity and for being an intermediate species between the dehydration of methanol and the formation of SMS, two processes favored under the studied conditions. An exception is the simulation of reference with no initial water molecules (red line), which shows a visible number of 5 methyl cations by the end of the simulations. This anomaly is likely a result of the absence of the mentioned promoting effect of water in this case. Fig. 10(b) presents a maximum of 4 methyl cations for the reference simulations, rather than 5, due to the smoothing process described in the legend.

Fig. S3 and S4 in the SI illustrate the effects of water on the profiles of the most relevant species involved in the MTH initiation reactions during simulations conducted at 1000 K. The overall trends at this temperature are similar to those observed in the simulation at 800 K for several profiles. For instance, water continues to enhance methanol consumption in a significant fraction of the trajectories after 300 ps. The difference relies on the insensitivity of the methanol conversion to water loads since, in these cases, the methanol profiles look almost indistinguishable. Water also seems to promote the formation of SMS, but only at the initial water load of 32 between 300 and 800 ps.

On the other hand, the rise in temperature also induces marked differences. The simulations run at 1000 K increase methanol consumption up to 8 molecules at the initial water load of 16, compared with the simulations performed at 800 K, and for the reference simulations and those with initial water loads of 8 and 32, an additional 4, 3, and 6 water molecules are consumed at the end of the trajectories. These results indicate that the combined action of temperature and initial water promotes methanol conversion. Regarding the SMS, their final number of molecules is more uniform at the end of the trajectories at 1000 K than at 800 K, being 7 for all the cases with initial water loads and 10 for the reference system. The most dramatic effect of temperature increase is on the water, hydronium ions, and methane profiles. The number of water molecules steadily increases for all cases, while the number of hydronium ions correspondingly decreases. The number of these ions is around 5 by the end of simulations run with initial water loads of 16 and 32, significantly lower than the 14 molecules observed for the equivalent trajectories at 800 K. An increased methanol conversion and reduced hydronium ion formation indicate that carbon and protons are available to react, leading to a significant production of methane. In other words, temperature favors methane formation, a behavior that agrees with experimental evidence, as presented in the work of Zhao et al.¹⁰⁶

Methane is an undesirable byproduct because it is inert and challenging to convert into valuable olefins and hydrocarbons. These characteristics and our simulations suggest that a temperature of 800 K is optimal for the initial stage of methanol conversion when the process considers initial water loads. This finding is consistent with experimental research that examines how temperature impacts catalytic performance. Wang et al.¹⁰⁷ studied the ethene-to-aromatic process in an H-ZSM-5 zeolite. The authors found that raising the temperature from 573 to 623 K shifts the production of deactivating coke toward alkylbenzenes and their cations, namely, toward beneficial HCP species. Similarly, Pérez-Uriarte et al.⁶³ establish an optimal temperature of 648 K for olefin selectivity starting from DME in an H-ZSM-5 zeolite. Zhao et al.¹⁰⁸ achieved maximum catalyst lifetime at 698 K when using a mixture of methanol and molecular hydrogen to synthesize olefins with H-SAPO-34 zeolite.

3.5. Discussion of the water effect on the initial stages of the MTH mechanism

The following subsections discuss several implications of the effect of water on the initiation reactions of the MTH process: modification of the zeolite acidity, methanol conversion under modified acidity conditions, hydrogen transfer reactions and framework activation induced by water, and chemistry of hydronium ions and SMS on the MTH process.

3.5.1. Modification of the zeolite acidity. The most significant implication of Fig. 10(a) is the widespread presence of hydronium ions, which alters the nature of zeolite acidity. Generally, this shift changes acidity from static to dynamic, significantly impacting reactions such as those occurring in the MTH process.¹⁰⁹ Several studies have highlighted that this dynamic acidity influences reactivity because of its considerable proton affinity compared with the BASs of the lattice and its influence on the electrostatic environment within the zeolite pores,²² which affects adsorption and diffusion properties of molecules, particularly those of hydrophobic characteristics.¹¹⁰

The most likely structures of protonated water are hydronium ions (H₃O⁺), Zundel cations (H₅O₂⁺), and protonated clusters with three or more water molecules, which tend to form at low temperatures and high water concentrations.^104,105 The simulations in this study did not show significant protonated water clusters due to consistently lower water concentrations than those of the BASs, resulting in an insufficient condition to form structures with high proton affinity. Besides, thermal effects at 800 K strongly disfavor the stability of protonated water clusters compared to room conditions.³⁷ Fig. S5 of the SI shows the profiles of the Zundel cation at 800 K and different water loads. The figure indicates that simulations with a water load of 32 present a maximum of two molecules of this cation through the trajectories, and only one molecule is detected in the other conditions.

It is relevant to mention that the work of Wang et al.¹⁰⁵ provides convincing evidence of a pronounced dynamic formation of hydronium ions, even at low water concentrations like those representing one water molecule per acid site. Franke and Simon¹¹¹ also supported the existence of hydronium ions at temperatures above 500 K, as long as water is at sufficient concentration within the zeolite interior. The results of these authors are significant because the hydronium ion profiles presented in Fig. 10(a) are from simulations run at representative water concentrations and high temperatures. In other words, the theoretical and experimental observation of acidity modification and hydronium ions formation agree with the reactive MD simulation performed in this work, which corresponds to an MTH process with comparable methanol, BASs, and water concentrations.

3.5.2. Methanol conversion under modified acidity conditions. Water can affect methanol conversion via competition for the acid center of a zeolite, solvation of methanol and transition states, shifting of chemical equilibrium, and modification of the electrostatic environment inside the zeolite.¹¹² Regarding the competition for adsorption on a zeolite BAS, Van der Mynsbrugge et al.²⁶ showed that water cannot compete with methanol at 670 K based on DFT calculations at the revPBE-D3/DZVP-GTH level of theory involving an MFI framework. These results contrast with the theoretical study of De Wispelaere et al.,¹¹³ who found a marked adsorption competition between methanol and water at 603 K, leading to lower methanol protonation. However, the zeolite used in this case had a CHA topology. The profile of methanol molecules presented in Fig. 9(a) qualitatively agrees with the conclusions of Van der Mynsbrugge et al.²⁶ from the theoretical viewpoint because water addition did not harm methanol protonation or conversion. Besides, adsorption of water and methanol was a rare event, likely due to the higher temperatures employed, that is, 800 K, which favored the formation of hydronium ions and protonated methanol molecules, as explained in Sections 3.2 and 3.5.1.

Water can also inhibit methanol reactivity via solvation. Nonetheless, recent research on ethanol dehydration on acidic zeolites (beta, TON, MFI, etc. frameworks) indicates that water concentration needs to be higher than that of the alcohol (up to seven times) to prevent reactivity.^103,114 Similar conditions should also apply to methanol. In the simulations conducted in this study, water is consistently present in substoichiometric amounts of methanol, making the formation of stable clusters impossible. Furthermore, the temperature of 800 K used in our simulations complicates the creation of such assembled structures, as in the case of water.³⁷

3.5.3. Hydrogen transfer reactions and framework activation induced by water. Competition for adsorption on a BAS is one of the most discussed mechanisms in the literature regarding how water can significantly influence methanol conversion.¹¹⁵ Since 800 K is a temperature that disfavors particular physical interactions in the systems under investigation, the adsorption competition can shift towards a protonation one. Fig. 10(a) shows that water significantly consumes zeolite protons, yet it does not hinder methanol protonation but encourages it. Fig. 11 presents the profiles of protonated methanol molecules as a function of initial water loads. The simulations of reference and that of an initial load of 32 host a maximum of 2 and 4 protonated methanol molecules, respectively, while the trajectories with intermediate initial water loads report up to 3 units of this cation. It should be noted that Fig. 11 presents mainly the trends, as the smoothing process described in the legend affects the absolute values.

Fig. 11 Number profiles of protonated methanol molecules as a function of time of simulations conducted at 800 K and the initial water loads of 0, 8, 16, and 32. This figure displays smoothed profiles due to the low values of the data (between 0 and 4) and their rapid variation through the trajectories. The smoothing consisted of averages in windows of 1 ps.

The reason behind the promoting effect of water to protonate methanol is linked with hydrogen transfer reactions, which also connects with the increased SMS production via polarization in the zeolite environment. The phenomenology of these processes can be described as follows: when water reacts with a BAS, it forms a hydronium ion that diffuses through the zeolite pores. Additionally, this reaction generates a negatively charged oxygen ion in the zeolite, effectively activating the framework. Panels (a) to (d) of Fig. 12 illustrate this stage. Experimental and theoretical evidence support the existence of a stabilizing electrostatic environment within the zeolite and relatively stable hydronium ions at high temperatures.^25,105,111 The diffusing hydronium ion can easily protonate methanol via a hydrogen transfer reaction, as presented in panels (e) to (g) of Fig. 12. Protonated methanol will also diffuse through the zeolite pores, and the initially activated framework will facilitate its dissociation to form an SMS and one water molecule, as presented for a similar process in panels (d) to (f) of Fig. 6. The overall reaction of zeolite methylation, i.e., CH₃OH + Z–OH ↔ H₂O + Z–OCH₃, remains unchanged. The difference relies on the presence of an extra water molecule, which promotes methanol protonation and its dissociation via a hydrogen transfer reaction and framework activation, respectively.

Fig. 12 Snapshots of water-mediated methanol protonation as shown by the reactive molecular dynamics (MD) trajectories. (a) A water molecule moves towards a Brønsted acid site (BAS) of the zeolite to (b) interact and (c) form a hydronium ion and leave the oxygen of the lattice negatively charged, namely, activate the framework. (d) The hydronium ion diffuses through the zeolite channels (e) until it encounters a methanol molecule. (f) Both species react (g) to recover the initial water and produce a protonated methanol species. The fate of this cation is to form water and surface methoxy species (SMS), where the initially activated framework plays a role. Panels (d) to (g) of Fig. 6 illustrate these subsequent steps.

Equilibrium arguments dictate that water addition would shift the equilibrium of the methanol dehydration reaction, namely, 2CH₃OH + Z–OH ↔ CH₃OCH₃ + H₂O + Z–OH, towards the reactants. De Wispelaere et al.¹¹³ suggested that such a shift could affect the equilibration stage of the MTO process and reduce DME formation. Similarly, Pérez-Uriarte et al.⁶³ conducted experimental studies on the conversion of DME-to-olefins (DTO) using an H-ZSM-5 zeolite. Their findings indicated that water addition reduced both DME conversion and the production of SMS intermediates, which are crucial for hydrocarbon production. Although this reasoning is correct, the simulation results of this work show that water at substoichiometric quantities does not play a thermodynamic role in the equilibration stage but acts as a “catalyst” that facilitates the protonation of methanol and SMS formation via hydronium ions and framework activation.

3.5.4. Chemistry of hydronium ions and SMS on the MTH process. Water-mediated protonation reactions in acidic zeolites are not a new topic, but no research addresses their impact on methanol protonation at high temperatures to the authors’ knowledge. An analysis of both experimental and theoretical studies suggests that a requirement for water to enhance protonation rates is to keep its concentration at substoichiometric levels,^112,116 a condition replicated by the simulations conducted in this work. A notable example from the literature illustrating the positive effect of water on protonation reactions is the study of Chen et al.¹¹⁷ They observed an increased rate of benzene hydrogen/deuterium exchange at room temperature when using an H-ZSM-5 zeolite with a Si/Al ratio equal to 15 and a water concentration of approximately one molecule per unit cell. This water load is near the value of 8 used in this work, which resulted in the highest production of SMS, as presented in Fig. 9(c). Examples of beneficial effects upon water addition are decreased coke formation and drops of free energy barriers (up to 8 kcal mol⁻¹) of hydrogen transfer reactions associated with ethene production through the aromatic cycle.⁶² However, hydronium ions can also dealuminate the zeolite and facilitate the collapse of the framework.³⁴

The water-induced framework activation observed in this work shares similarities with framework methylation studies conducted by Matam et al. at room temperature.^118,119 The authors found that high methanol concentration favors hydrogen transfer reactions and SMS formation, and it is hypothesized that framework deprotonation promoted at these conditions decreases the activation barrier of the corresponding methylation reactions. Similarly, the framework deprotonation caused by water in the systems studied here may partially lower the activation barrier for zeolite methylation, namely, SMS formation.²⁵ This fact is of paramount importance because SMS play an active role at different stages of methanol conversion. Specifically, SMS is significant during the induction period, the steady formation of olefins, the subsequent development of secondary products, and, ultimately, the deactivation phase.^21,72

In summary, our simulations reveal that for systems resembling the equilibration of the MTH process in a H-ZSM-5 catalyst, the presence of water at substoichiometric quantities primarily promotes the formation of hydronium ions that diffuse through the zeolite pores in virtue of the electrostatic environment within the framework and high operation temperatures (i.e., 800 K). Besides, hydronium ions induce zeolite activation by leaving negatively charged oxygen atoms around the original BASs, favoring the formation of SMS intermediates. Theoretical and experimental evidence indicate that hydronium ions and SMS participate in relevant reactions through the different stages of methanol conversion. Our simulations showed that hydronium ions can protonate methanol, and the literature illustrates how these ions can assist the ethylation, dealkylation, and proton exchange in aromatics, reduce the activation barrier in the aromatic cycle steps leading to ethylene production, and participate in zeolite dealumination reactions. In addition, SMS plays a primary role in the first carbon–carbon coupling, steady olefin production, and deactivating PAHs and coke formation. Therefore, our simulations suggest that water also impacts the MTH process through the chemistry of species such as hydronium ions and SMS, reflecting dynamic effects that have not been fully considered.

4. Conclusions

This contribution presents a newly developed reactive ReaxFF H/C/O/Si/Al force field that describes the initial reactions of methanol conversion in acidic zeolites. The force field is built on the previous ReaxFF development, focused on studying the dynamics of water confined in smectite clay–zeolite composites.³⁹ This work incorporates new training of the ReaxFF C/O parameters concerning bond dissociation energies, valence angles, torsions, and the challenging carbon monoxide–zeolite energetics. The parametrization used DFT data, including reaction barriers and intermediate energies, for the first stages of the methanol-to-hydrocarbon (MTH) reactions. The new ReaxFF force field was utilized to investigate the dynamic aspects of the equilibration reactions involved in the conversion of methanol over an H-ZSM-5 zeolite. The molecular model consisted of a 2 × 2 × 2 supercell with MFI topology and a Si/Al ratio equal to 23. An initial set of simulations examined how temperatures (600, 800, 1000, and 1200 K) affected the methanol conversion mechanism, a second set tested the effect of initial water loads of 8, 16, and 32 at 800 K, and a third set also tested different water loads but at 1000 K. The trajectories were run for 1000 ps in the NPT ensemble to assess the relevant aspects of the framework flexibility.

The temperature increases the conversion of methanol and the concomitant water production, but also promotes the formation of undesired methane instead of the relevant intermediate surface methoxy species (SMS). This behavior qualitatively agrees with experimental works reported in the literature. The kinetics of methanol consumption exhibited an apparent activation energy of 5.30 kcal mol⁻¹, similar to values derived from empirical models that describe methanol conversion in various types of reactors and zeolites, including H-ZSM-5. A notable manifestation of dynamic effect observed in the simulations was the diffusion of charged molecules like hydronium ions, methyl cations, and protonated methanol. Diffusion of cations is relevant because it facilitates the formation of complex intermediates with loose structures that lead to the critical SMS and dimethyl ether (DME). Entropy factors at the tested temperatures could favor reactions via such loose structures, thanks to their associated lower free energy barriers compared to more rigid geometries involving the zeolite surface. For example, DME forms from the interaction of a diffusing methanol molecule and a methyl cation rather than from a diffusing methanol molecule and an SMS. This trend may extend to multiple stages of methanol conversion. However, above 1000 K, the beneficial effects of temperature diminish because the entropy contribution becomes excessively large, which can result in speculative reactions during reactive molecular dynamics (MD) simulations.

Initial water loads at comparable Brønsted acid site (BAS) quantities strongly influence the dynamics of acidic zeolites at 800 K. The principal outcome is the pronounced formation of hydronium ions, which indicates a change of zeolite acidity from static to dynamic. This process also activates the framework by leaving a negatively charged oxygen atom in the zeolite lattice. Theoretical and experimental evidence show that hydronium ions prevail more than protonated water clusters at the initial water loads and temperatures studied. One of the most widespread mechanisms by which water may affect methanol conversion is via competition for adsorption on the BAS of the zeolite or, in the cases presented in this work, competition for protonation. Nonetheless, the simulations performed in this work displayed an enhanced conversion of methanol and SMS formation instead of demonstrating a detrimental effect. Although theoretical investigations have proved that water cannot compete with methanol for the acid centers of an H-ZSM-5 zeolite, the increase in methanol reactivity originates from hydrogen transfer reactions facilitated by hydronium ions, exemplifying another distinctive manifestation of dynamic effects. Furthermore, the initially activated framework (resulting from the reaction between initial water and BASs to produce hydronium ions) encourages the dissociation of protonated methanol to form SMS and additional water. Previous studies show that high methanol/BAS ratios similarly favor SMS formation due to the preceding framework deprotonation. This scenario implies that initial water acts as a “catalyst” to methylate the zeolite in the equilibration stage of the MTH process. From a general viewpoint, water-mediated hydrogenation reactions are widespread steps in the catalytic mechanisms of zeolites, and SMS plays a pivotal role as an active reactant in the equilibration and induction periods, eventual production of olefins, and catalyst deactivation. Our findings suggest that water affects methanol conversion through complex and dynamic phenomena that demand further investigation.

Author contributions

Edwing Grajales-Gonzalez: conceptualization, investigation, formal analysis, methodology, writing – original draft, writing – review and editing. Adri C. T. van Duin: conceptualization, investigation, methodology, supervision, resources, writing – review and editing. S. Mani Sarathy: conceptualization, funding acquisition, resources, supervision, writing – review and editing.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

Data will be made available on request and in the online supplementary information (SI). Supplementary information: Table S1 lists percentage errors of lattice parameters determined using the force field developed in this work. Fig. S1 displays the stationary points for the water-assisted proton transport in the first coordination sphere of an aluminum atom. Fig. S2 illustrates DFT and ReaxFF potential energy diagrams of the water-mediated proton transport as a function of the forming ZO-H bond length. Fig. S3 presents number profiles of methanol, water, surface methoxy species, and methane as a function of time of simulations conducted at 1000 K and initial water loads of 0, 8, 16, and 32. Fig. S4 shows number profiles of hydronium ions and methyl cation species as a function of time of simulations conducted at 1000 K and initial water loads of 0, 8, 16, and 32. Fig. S5 shows the number profiles of the Zundel cation as a function of time of simulations conducted at 800 K and initial water loads of 0, 8, 16, and 32. See DOI: https://doi.org/10.1039/d5cp02704g.

Acknowledgements

E. G.-G. and S. M. S. acknowledge the resources of the Supercomputing Laboratory at KAUST (Shaheen and Ibex supercomputers). This research was funded by KAUST.

References

1. L. R. Brown, The great transition: Shifting from fossil fuels to solar and wind energy, WW Norton & Company, 2015.
2. International Energy Agency, World Energy Outlook 2024, Paris, 2024, https://www.iea.org/reports/world-energy-outlook-2024.
3. T. Brown and J. Hampp, Ultra-long-duration energy storage anywhere: Methanol with carbon cycling, Joule, 2023, 7, 2414–2420.
4. M. Hirscher, V. A. Yartys, M. Baricco, J. Bellosta von Colbe, D. Blanchard and R. C. Bowman, et al., Materials for hydrogen-based energy storage – past, recent progress and future outlook, J. Alloys Compd., 2020, 827, 153548.
5. S. S. Nagaraja, S. M. Sarathy, B. Mohan and J. Chang, Machine learning-driven screening of fuel additives for increased spark-ignition engine efficiency, Proc. Combust. Inst., 2024, 40, 105658.
6. M. Bertau, H. Offermanns, L. Plass, F. Schmidt and H.-J. Wernicke, Methanol: the basic chemical and energy feedstock of the future, Springer, 2014.
7. F. Dalena, A. Senatore, A. Marino, A. Gordano, M. Basile and A. Basile, Methanol Production and Applications: An Overview, in Methanol, ed. A. Basile and F. Dalena, Elsevier, 2018, ch. 1, pp. 3–28.
8. T. J. Deka, A. I. Osman, D. C. Baruah and D. W. Rooney, Methanol fuel production, utilization, and techno-economy: a review, Environ. Chem. Lett., 2022, 20, 3525–3554.
9. G. Olah, A. Goeppert and S. Prakash, The “Methanol Economy”: General Aspects, Beyond Oil and Gas: The Methanol Economy, John Wiley & Sons, Ltd, 2009, pp. 179–184.
10. S. Verhelst, J. W. G. Turner, L. Sileghem and J. Vancoillie, Methanol as a fuel for internal combustion engines, Prog. Energy Combust. Sci., 2019, 70, 43–88.
11. S. S. Tabibian and M. Sharifzadeh, Statistical and analytical investigation of methanol applications, production technologies, value-chain and economy with a special focus on renewable methanol, Renewable Sustainable Energy Rev., 2023, 179, 113281.
12. P. Tian, Y. Wei, M. Ye and Z. Liu, Methanol to Olefins (MTO): From Fundamentals to Commercialization, ACS Catal., 2015, 5, 1922–1938.
13. K. Hemelsoet, J. Van der Mynsbrugge, K. De Wispelaere, M. Waroquier and V. Van Speybroeck, Unraveling the Reaction Mechanisms Governing Methanol-to-Olefins Catalysis by Theory and Experiment, ChemPhysChem, 2013, 14, 1526–1545.
14. L. Yang, T. Yan, C. Wang, W. Dai, G. Wu and M. Hunger, et al., Role of Acetaldehyde in the Roadmap from Initial Carbon–Carbon Bonds to Hydrocarbons during Methanol Conversion, ACS Catal., 2019, 9, 6491–6501.
15. D. Lesthaeghe, V. Van Speybroeck, G. B. Marin and M. Waroquier, Understanding the Failure of Direct C-C Coupling in the Zeolite-Catalyzed Methanol-to-Olefin Process, Angew. Chem., Int. Ed., 2006, 45, 1714–1719.
16. Y. Liu, S. Müller, D. Berger, J. Jelic, K. Reuter and M. Tonigold, et al., Formation Mechanism of the First Carbon–Carbon Bond and the First Olefin in the Methanol Conversion into Hydrocarbons, Angew. Chem., Int. Ed., 2016, 55, 5723–5726.
17. Z. Wei, Y.-Y. Chen, J. Li, P. Wang, B. Jing and Y. He, et al., Methane formation mechanism in the initial methanol-to-olefins process catalyzed by SAPO-34, Catal. Sci. Technol., 2016, 6, 5526–5533.
18. J. Li, Z. Wei, Y. Chen, B. Jing, Y. He and M. Dong, et al., A route to form initial hydrocarbon pool species in methanol conversion to olefins over zeolites, J. Catal., 2014, 317, 277–283.
19. S. Wang, Z. Qin, M. Dong, J. Wang and W. Fan, Recent progress on MTO reaction mechanisms and regulation of acid site distribution in the zeolite framework, Chem Catal., 2022, 2, 1657–1685.
20. M. Vandichel, D. Lesthaeghe, J. Van der Mynsbrugge, M. Waroquier and V. Van Speybroeck, Assembly of cyclic hydrocarbons from ethene and propene in acid zeolite catalysis to produce active catalytic sites for MTO conversion, J. Catal., 2010, 271, 67–78.
21. I. Yarulina, A. D. Chowdhury, F. Meirer, B. M. Weckhuysen and J. Gascon, Recent trends and fundamental insights in the methanol-to-hydrocarbons process, Nat. Catal., 2018, 1, 398–411.
22. K. Stanciakova and B. M. Weckhuysen, Water–active site interactions in zeolites and their relevance in catalysis, Trends Chem., 2021, 3, 456–468.
23. I. M. Dahl and S. Kolboe, On the Reaction Mechanism for Hydrocarbon Formation from Methanol over SAPO-34: 2. Isotopic Labeling Studies of the Co-reaction of Propene and Methanol, J. Catal., 1996, 161, 304–309.
24. M. Guisnet, L. Costa and F. R. Ribeiro, Prevention of zeolite deactivation by coking, J. Mol. Catal. A: Chem., 2009, 305, 69–83.
25. S. A. F. Nastase, P. Cnudde, L. Vanduyfhuys, K. De Wispelaere, V. Van Speybroeck and C. R. A. Catlow, et al., Mechanistic Insight into the Framework Methylation of H-ZSM-5 for Varying Methanol Loadings and Si/Al Ratios Using First-Principles Molecular Dynamics Simulations, ACS Catal., 2020, 10, 8904–8915.
26. J. Van der Mynsbrugge, S. L. C. Moors, K. De Wispelaere and V. Van Speybroeck, Insight into the Formation and Reactivity of Framework-Bound Methoxide Species in H-ZSM-5 from Static and Dynamic Molecular Simulations, ChemCatChem, 2014, 6, 1906–1918.
27. A. Ghorbanpour, J. D. Rimer and L. C. Grabow, Computational Assessment of the Dominant Factors Governing the Mechanism of Methanol Dehydration over H-ZSM-5 with Heterogeneous Aluminum Distribution, ACS Catal., 2016, 6, 2287–2298.
28. M. Fečík, P. N. Plessow and F. Studt, Simple Scheme to Predict Transition-State Energies of Dehydration Reactions in Zeolites with Relevance to Biomass Conversion, J. Phys. Chem. C, 2018, 122, 23062–23067.
29. P. N. Plessow, A. Smith, S. Tischer and F. Studt, Identification of the Reaction Sequence of the MTO Initiation Mechanism Using Ab Initio-Based Kinetics, J. Am. Chem. Soc., 2019, 141, 5908–5915.
30. S. Wang, Y. Chen, Z. Qin, T.-S. Zhao, S. Fan and M. Dong, et al., Origin and evolution of the initial hydrocarbon pool intermediates in the transition period for the conversion of methanol to olefins over H-ZSM-5 zeolite, J. Catal., 2019, 369, 382–395.
31. P. N. Plessow and F. Studt, Unraveling the Mechanism of the Initiation Reaction of the Methanol to Olefins Process Using ab Initio and DFT Calculations, ACS Catal., 2017, 7, 7987–7994.
32. P. N. Plessow and F. Studt, Theoretical Insights into the Effect of the Framework on the Initiation Mechanism of the MTO Process, Catal. Lett., 2018, 148, 1246–1253.
33. C. Bai, L. Liu and H. Sun, Molecular Dynamics Simulations of Methanol to Olefin Reactions in HZSM-5 Zeolite Using a ReaxFF Force Field, J. Phys. Chem. C, 2012, 116, 7029–7039.
34. E. Grajales-González, M. Monge-Palacios, S. Brunialti and S. M. Sarathy, Mechanistic Insights into the Dealumination of an H-ZSM-5 Zeolite Using Reactive Molecular Dynamics Simulations, ACS Catal., 2024, 14, 14110–14126.
35. K. L. Joshi and A. C. T. van Duin, Molecular Dynamics Study on the Influence of Additives on the High-Temperature Structural and Acidic Properties of ZSM-5 Zeolite, Energy Fuels, 2013, 27, 4481–4488.
36. G. M. Psofogiannakis, J. F. McCleerey, E. Jaramillo and A. C. T. van Duin, ReaxFF Reactive Molecular Dynamics Simulation of the Hydration of Cu-SSZ-13 Zeolite and the Formation of Cu Dimers, J. Phys. Chem. C, 2015, 119, 6678–6686.
37. K. L. Joshi, G. Psofogiannakis, A. C. T. van Duin and S. Raman, Reactive molecular simulations of protonation of water clusters and depletion of acidity in H-ZSM-5 zeolite, Phys. Chem. Chem. Phys., 2014, 16, 18433–18441.
38. M. Zheng, X. Li, L. Guo and W. Ge, Dynamic Intermediate Profiles of Zeolite Catalyzed Methanol to Olefins Revealed by Reactive Molecular Dynamics, Energy Fuels, 2021, 35, 1677–1690.
39. M. C. Pitman and A. C. T. van Duin, Dynamics of Confined Reactive Water in Smectite Clay–Zeolite Composites, J. Am. Chem. Soc., 2012, 134, 3042–3053.
40. Tersoff, Empirical interatomic potential for silicon with improved elastic properties, 1988.
41. K. S. Smirnov and B. van de Graaf, Consistent implementation of the electronegativity equalization method in molecular mechanics and molecular dynamics, J. Chem. Soc., Faraday Trans., 1996, 92, 2469–2474.
42. T. P. Senftle, S. Hong, M. M. Islam, S. B. Kylasa, Y. Zheng and Y. K. Shin, et al., The ReaxFF reactive force-field: development, applications and future directions, npj Comput. Mater., 2016, 2, 15011.
43. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, ReaxFF: A Reactive Force Field for Hydrocarbons, J. Phys. Chem. A, 2001, 105, 9396–9409.
44. K. Chenoweth, A. C. T. van Duin and W. A. Goddard, ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation, J. Phys. Chem. A, 2008, 112, 1040–1053.
45. J. C. Fogarty, H. M. Aktulga, A. Y. Grama, A. C. T. van Duin and S. A. Pandit, A reactive molecular dynamics simulation of the silica-water interface, J. Chem. Phys., 2010, 132, 174704.
46. M. F. Russo, R. Li, M. Mench and A. C. T. van Duin, Molecular dynamic simulation of aluminum–water reactions using the ReaxFF reactive force field, Int. J. Hydrogen Energy, 2011, 36, 5828–5835.
47. H. Manzano, S. Moeini, F. Marinelli, A. C. T. van Duin, F.-J. Ulm and R. J.-M. Pellenq, Confined Water Dissociation in Microporous Defective Silicates: Mechanism, Dipole Distribution, and Impact on Substrate Properties, J. Am. Chem. Soc., 2012, 134, 2208–2215.
48. Y. Tian, J. Du, X. Zu, W. Han, X. Yuan and W. Zheng, UV-induced modification of fused silica: Insights from ReaxFF-based molecular dynamics simulations, AIP Adv., 2016, 6, 95312.
49. Z. Chen, P. Zhao, L. Zhao and W. Sun, Molecular Simulation of the Catalytic Cracking of Hexadecane on ZSM-5 Catalysts Based on Reactive Force Field (ReaxFF), Energy Fuels, 2017, 31, 10515–10524.
50. A. C. T. van Duin, J. M. A. Baas and B. van de Graaf, Delft molecular mechanics: a new approach to hydrocarbon force fields. Inclusion of a geometry-dependent charge calculation, J. Chem. Soc., Faraday Trans., 1994, 90, 2881–2895.
51. Z. Wei, Y.-Y. Chen, J. Li, W. Guo, S. Wang and M. Dong, et al., Stability and Reactivity of Intermediates of Methanol Related Reactions and C–C Bond Formation over H-ZSM-5 Acidic Catalyst: A Computational Analysis, J. Phys. Chem. C, 2016, 120, 6075–6087.
52. M. Trachta, M. Rubeš and O. Bludský, Toward accurate ab initio modeling of siliceous zeolite structures, J. Chem. Phys., 2022, 156, 94708.
53. S. F. Yuk, I. Sargin, N. Meyer, J. T. Krogel, S. P. Beckman and V. R. Cooper, Putting error bars on density functional theory, Sci. Rep., 2024, 14, 20219.
54. J. A. Ryder, A. K. Chakraborty and A. T. Bell, Density Functional Theory Study of Proton Mobility in Zeolites: Proton Migration and Hydrogen Exchange in ZSM-5, J. Phys. Chem. B, 2000, 104, 6998–7011.
55. Structure Commission of the International Zeolite Association. Database of Zeolite Structures, 2022.
56. C. Baerlocher, L. B. McCusker and D. H. Olson, MFI – Pnma, in Atlas of Zeolite Framework Types, ed. C. Baerlocher, L. B. McCusker and D. H. Olson, Elsevier Science B.V., Amsterdam, 6th edn, 2007, pp. 212–213.
57. K. Momma and F. Izumi, VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Crystallogr., 2011, 44, 1272–1276.
58. K. Margeta, N. Z. Logar, M. Šiljeg and A. Farkas, Natural Zeolites in Water Treatment – How Effective is Their Use, in Water Treatment, ed. W. Elshorbagy and R. K. Chowdhury, IntechOpen, Rijeka, 2013 DOI:10.5772/50738.
59. S. Esposito, A. Marocco, G. Dell’Agli, B. De Gennaro and M. Pansini, Relationships between the water content of zeolites and their cation population, Microporous Mesoporous Mater., 2015, 202, 36–43.
60. M. Ravi, V. L. Sushkevich and J. A. van Bokhoven, Towards a better understanding of Lewis acidic aluminium in zeolites, Nat. Mater., 2020, 19, 1047–1056.
61. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, PACKMOL: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem., 2009, 30, 2157–2164.
62. C. Zhang, X. Wu, Y. Zhang, W. Zhang, S. Lin and C. Lou, et al., Water-assisted shape-selective production of ethene in methanol-to-olefins reaction on SAPO-34, Chem Catal., 2024, 4, 101025.
63. P. Pérez-Uriarte, A. Ateka, M. Gamero, A. T. Aguayo and J. Bilbao, Effect of the Operating Conditions in the Transformation of DME to olefins over a HZSM-5 Zeolite Catalyst, Ind. Eng. Chem. Res., 2016, 55, 6569–6578.
64. C. Ashraf and A. C. T. van Duin, Extension of the ReaxFF Combustion Force Field toward Syngas Combustion and Initial Oxidation Kinetics, J. Phys. Chem. A, 2017, 121, 1051–1068.
65. E. Grajales-González, M. Monge-Palacios and S. M. Sarathy, Atomistic simulations of syngas oxy-combustion in supercritical CO₂, J. CO₂ Util., 2021, 49, 101554.
66. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19.
67. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81, 511–519.
68. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A:At., Mol., Opt. Phys., 1985, 31, 1695–1697.
69. K. De Wispelaere, B. Ensing, A. Ghysels, E. J. Meijer and V. Van Speybroeck, Complex Reaction Environments and Competing Reaction Mechanisms in Zeolite Catalysis: Insights from Advanced Molecular Dynamics, Chem. – Eur. J., 2015, 21, 9385–9396.
70. N. Nesterenko, J. Aguilhon, P. Bodart, D. Minoux and J.-P. Dath, Methanol to Olefins: An Insight Into Reaction Pathways and Products Formation, in Zeolites and Zeolite-Like Materials, ed. B. F. Sels and L. M. Kustov, Elsevier, Amsterdam, 2016, ch. 5, pp. 189–263.
71. A. C. T. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu and W. A. Goddard, ReaxFFSiO Reactive Force Field for Silicon and Silicon Oxide Systems, J. Phys. Chem. A, 2003, 107, 3803–3811.
72. Y. Jiang, M. Hunger and W. Wang, On the Reactivity of Surface Methoxy Species in Acidic Zeolites, J. Am. Chem. Soc., 2006, 128, 11679–11692.
73. W. Wang and M. Hunger, Reactivity of Surface Alkoxy Species on Acidic Zeolite Catalysts, Acc. Chem. Res., 2008, 41, 895–904.
74. S. Y. Kim, D. Perez and A. F. Voter, Local hyperdynamics, J. Chem. Phys., 2013, 139, 144110.
75. R. Gounder and E. Iglesia, The Roles of Entropy and Enthalpy in Stabilizing Ion-Pairs at Transition States in Zeolite Acid Catalysis, Acc. Chem. Res., 2012, 45, 229–238.
76. X. Wu, S. Xu, W. Zhang, J. Huang, J. Li and B. Yu, et al., Direct Mechanism of the First Carbon–Carbon Bond Formation in the Methanol-to-Hydrocarbons Process, Angew. Chem., Int. Ed., 2017, 56, 9039–9043.
77. W. Wang, A. Buchholz, M. Seiler and M. Hunger, Evidence for an Initiation of the Methanol-to-Olefin Process by Reactive Surface Methoxy Groups on Acidic Zeolite Catalysts, J. Am. Chem. Soc., 2003, 125, 15260–15267.
78. G. J. Hutchings, F. Gottschalk, M. V. M. Hall and R. Hunter, Hydrocarbon formation from methylating agents over the zeolite catalyst ZSM-5. Comments on the mechanism of carbon–carbon bond and methane formation, J. Chem. Soc., Faraday Trans. 1, 1987, 83, 571–583.
79. A. Vashisth, C. Ashraf, W. Zhang, C. E. Bakis and A. C. T. van Duin, Accelerated ReaxFF Simulations for Describing the Reactive Cross-Linking of Polymers, J. Phys. Chem. A, 2018, 122, 6633–6642.
80. E. Braun, J. Gilmer, H. Mayes, D. Mobley, J. Monroe and S. Prasad, et al., Best Practices for Foundations in Molecular Simulations, Living J. Comput. Mol. Sci., 2019, 1, 1–28.
81. L. Qi, Y. Wei, L. Xu and Z. Liu, Reaction Behaviors and Kinetics during Induction Period of Methanol Conversion on HZSM-5 Zeolite, ACS Catal., 2015, 5, 3973–3982.
82. A. D. Chowdhury, A. L. Paioni, K. Houben, G. T. Whiting, M. Baldus and B. M. Weckhuysen, Bridging the Gap between the Direct and Hydrocarbon Pool Mechanisms of the Methanol-to-Hydrocarbons Process, Angew. Chem., Int. Ed., 2018, 57, 8095–8099.
83. P. Losch, H. R. Joshi, O. Vozniuk, A. Grünert, C. Ochoa-Hernández and H. Jabraoui, et al., Proton Mobility, Intrinsic Acid Strength, and Acid Site Location in Zeolites Revealed by Varying Temperature Infrared Spectroscopy and Density Functional Theory Studies, J. Am. Chem. Soc., 2018, 140, 17790–17799.
84. A. Dyer, Ion-exchange properties of zeolites, in Zeolites and Ordered Mesoporous Materials: Progress and Prospects, ed. J. Ĉejka and H. van Bekkum, Elsevier, 2005, pp. 181–204.
85. P. Pérez-Uriarte, A. Ateka, A. T. Aguayo, A. G. Gayubo and J. Bilbao, Kinetic model for the reaction of DME to olefins over a HZSM-5 zeolite catalyst, Chem. Eng. J., 2016, 302, 801–810.
86. M. Kaarsholm, B. Rafii, F. Joensen, R. Cenni, J. Chaouki and G. S. Patience, Kinetic Modeling of Methanol-to-Olefin Reaction over ZSM-5 in Fluid Bed, Ind. Eng. Chem. Res., 2010, 49, 29–38.
87. D. Chen, A. Grønvold, K. Moljord and A. Holmen, Methanol Conversion to Light Olefins over SAPO-34:  Reaction Network and Deactivation Kinetics, Ind. Eng. Chem. Res., 2007, 46, 4116–4123.
88. S. Standl, F. M. Kirchberger, T. Kühlewind, M. Tonigold, M. Sanchez-Sanchez and J. A. Lercher, et al., Single-event kinetic model for methanol-to-olefins (MTO) over ZSM-5: Fundamental kinetics for the olefin co-feed reactivity, Chem. Eng. J., 2020, 402, 126023.
89. S. R. Blaszkowski and R. A. van Santen, The Mechanism of Dimethyl Ether Formation from Methanol Catalyzed by Zeolitic Protons, J. Am. Chem. Soc., 1996, 118, 5152–5153.
90. H. K. Chau, H. D. Mai, A. Gumidyala, T. N. Pham, D.-P. Bui and A. D. D’Amico, et al., Effect of Water on Cumene Dealkylation over H-ZSM-5 Zeolites, ACS Catal., 2023, 13, 9158–9170.
91. W. O. Haag, R. M. Lago and P. G. Rodewald, Aromatics, light olefins and gasoline from methanol: Mechanistic pathways with ZSM-5 zeolite catalyst, J. Mol. Catal., 1982, 17, 161–169.
92. A. Hwang, M. Kumar, J. D. Rimer and A. Bhan, Implications of methanol disproportionation on catalyst lifetime for methanol-to-olefins conversion by HSSZ-13, J. Catal., 2017, 346, 154–160.
93. Y. Liu, F. M. Kirchberger, S. Müller, M. Eder, M. Tonigold and M. Sanchez-Sanchez, et al., Critical role of formaldehyde during methanol conversion to hydrocarbons, Nat. Commun., 2019, 10, 1462.
94. C. Chizallet, C. Bouchy, K. Larmier and G. Pirngruber, Molecular Views on Mechanisms of Brønsted Acid-Catalyzed Reactions in Zeolites, Chem. Rev., 2023, 123, 6107–6196.
95. B. Xu, C. Sievers, S. B. Hong, R. Prins and J. A. van Bokhoven, Catalytic activity of Brønsted acid sites in zeolites: Intrinsic activity, rate-limiting step, and influence of the local structure of the acid sites, J. Catal., 2006, 244, 163–168.
96. S. Li, J. Cao, X. Feng, Y. Du, D. Chen and C. Yang, et al., Insights into the confinement effect on isobutane alkylation with C4 olefin catalyzed by zeolite catalyst: A combined theoretical and experimental study, Chin. J. Chem. Eng., 2022, 47, 174–184.
97. R. Osuga, T. Yokoi, K. Doitomi, H. Hirao and J. N. Kondo, Infrared Investigation of Dynamic Behavior of Brønsted Acid Sites on Zeolites at High Temperatures, J. Phys. Chem. C, 2017, 121, 25411–25420.
98. P. Huber, F. Studt and P. N. Plessow, Reactivity of Surface Lewis and Brønsted Acid Sites in Zeolite Catalysis: A Computational Case Study of DME Synthesis Using H-SSZ-13, J. Phys. Chem. C, 2022, 126, 5896–5905.
99. A. J. Jones and E. Iglesia, Kinetic, Spectroscopic, and Theoretical Assessment of Associative and Dissociative Methanol Dehydration Routes in Zeolites, Angew. Chem., Int. Ed., 2014, 53, 12177–12181.
100. M. Signorile, E. Borfecchia, S. Bordiga and G. Berlier, Influence of ion mobility on the redox and catalytic properties of Cu ions in zeolites, Chem. Sci., 2022, 13, 10238–10250.
101. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, et al., Gaussian16 (Revision A.03), Wallingford, CT, 2016.
102. R. Gounder and E. Iglesia, The catalytic diversity of zeolites: confinement and solvation effects within voids of molecular dimensions, Chem. Commun., 2013, 49, 3491–3509.
103. J. S. Bates, B. C. Bukowski, J. Greeley and R. Gounder, Structure and solvation of confined water and water–ethanol clusters within microporous Brønsted acids and their effects on ethanol dehydration catalysis, Chem. Sci., 2020, 11, 7102–7122.
104. P. Liu and D. Mei, Identifying Free Energy Landscapes of Proton-Transfer Processes between Brønsted Acid Sites and Water Clusters Inside the Zeolite Pores, J. Phys. Chem. C, 2020, 124, 22568–22576.
105. M. Wang, N. R. Jaegers, M.-S. Lee, C. Wan, J. Z. Hu and H. Shi, et al., Genesis and Stability of Hydronium Ions in Zeolite Channels, J. Am. Chem. Soc., 2019, 141, 3444–3455.
106. W. Zhao, B. Zhang, G. Wang and H. Guo, Methane formation route in the conversion of methanol to hydrocarbons, J. Energy Chem., 2014, 23, 201–206.
107. H. Wang, Y. Hou, W. Sun, Q. Hu, H. Xiong and T. Wang, et al., Insight into the Effects of Water on the Ethene to Aromatics Reaction with HZSM-5, ACS Catal., 2020, 10, 5288–5298.
108. X. Zhao, J. Li, P. Tian, L. Wang, X. Li and S. Lin, et al., Achieving a Superlong Lifetime in the Zeolite-Catalyzed MTO Reaction under High Pressure: Synergistic Effect of Hydrogen and Water, ACS Catal., 2019, 9, 3017–3025.
109. C. J. Heard, L. Grajciar and P. Nachtigall, The effect of water on the validity of Löwenstein's rule, Chem. Sci., 2019, 10, 5705–5711.
110. G. Li, B. Wang and D. E. Resasco, Water-Mediated Heterogeneously Catalyzed Reactions, ACS Catal., 2020, 10, 1294–1309.
111. M. E. Franke and U. Simon, Solvate-Supported Proton Transport in Zeolites, ChemPhysChem, 2004, 5, 465–472.
112. Q. Liu and J. A. van Bokhoven, Water structures on acidic zeolites and their roles in catalysis, Chem. Soc. Rev., 2024, 53, 3065–3095.
113. K. De Wispelaere, C. S. Wondergem, B. Ensing, K. Hemelsoet, E. J. Meijer and B. M. Weckhuysen, et al., Insight into the Effect of Water on the Methanol-to-Olefins Conversion in H-SAPO-34 from Molecular Simulations and in Situ Microspectroscopy, ACS Catal., 2016, 6, 1991–2002.
114. J. S. Bates and R. Gounder, Influence of confining environment polarity on ethanol dehydration catalysis by Lewis acid zeolites, J. Catal., 2018, 365, 213–226.
115. K. Stanciakova, J. N. Louwen, B. M. Weckhuysen, R. E. Bulo and F. Göltl, Understanding Water–Zeolite Interactions: On the Accuracy of Density Functionals, J. Phys. Chem. C, 2021, 125, 20261–20274.
116. D. E. Resasco, S. P. Crossley, B. Wang and J. L. White, Interaction of water with zeolites: a review, Catal. Rev., 2021, 63, 302–362.
117. K. Chen, A. Gumidyala, M. Abdolrhamani, C. Villines, S. Crossley and J. L. White, Trace water amounts can increase benzene H/D exchange rates in an acidic zeolite, J. Catal., 2017, 351, 130–135.
118. S. K. Matam, R. F. Howe, A. Thetford and C. R. A. Catlow, Room temperature methoxylation in zeolite H-ZSM-5: an operando DRIFTS/mass spectrometric study, Chem. Commun., 2018, 54, 12875–12878.
119. S. K. Matam, S. A. F. Nastase, A. J. Logsdail and C. R. A. Catlow, Methanol loading dependent methoxylation in zeolite H-ZSM-5, Chem. Sci., 2020, 11, 6805–6814.

---

Footnote

† Present address: Energy and Particle Technology Laboratory, Department of Mechanical and Aerospace Engineering, Carleton University, 1125 Colonel by Drive, Ottawa, ON, K1S 5B6, Canada.

---