A computational study of the formation of surface methoxy species in H-SSZ-13 and H-SAPO-34 frameworks

Abstract

The methanol-to-hydrocarbons (MTH) reaction on zeolites is vital for the production of higher-order hydrocarbons from sustainable C₁ feedstocks. The formation of the first C–C bond is a key initiation step in the MTH reaction, and surface methoxy species (SMS) are important in many of the proposed pathways to form C₂ species, but the reaction steps that form SMS in zeotype frameworks remain uncertain. Therefore, we have investigated the reaction energies and activation barriers for SMS formation pathways in zeotype frameworks using accurate ab initio simulations, considering isostructural aluminosilicate and aluminophospate frameworks to allow scrutiny of how catalyst composition affects reaction steps. The SMS precursors dimethyl ether (DME) and trimethyl oxonium (TMO) are found to form directly from methanol with relatively low barriers (62 kJ mol⁻¹ and 94 kJ mol⁻¹, respectively); and the protonated forms of DME and CH₃OH, as well as TMO, form SMS with low kinetic barriers (36–48 kJ mol⁻¹). The activation barriers for processes occurring on H-SAPO-34 are consistently 10–23 kJ mol⁻¹ higher than reactions on the isostructural H-SSZ-13 framework, indicating that only kinetic differences exist between aluminosilicates and aluminophosphates for SMS formation. Significant differences are identified in the activation energies for reactions that proceed through the front-side attack S_N2 when compared to the back-side attack S_N2 mechanisms, with reduced electron donation to the carbocation intermediate leading to instability of the front-side attack S_N2 intermediate. Overall, the direct framework methylation step via protonated methanol has the lowest kinetic barrier, which agrees with experimental observations of direct SMS formation, and this result provides a foundation for further mechanistic investigations.

---

1 Introduction

The methanol-to-hydrocarbons (MTH) process is an economically viable route to convert methanol feedstocks to a variety of value added products such as gasoline (methanol-to-gasoline, MTG) and precursors for polymer materials (methanol-to-olefins, MTO).¹ In combination with biomass processing techniques such as biogas-to-liquid (BtL),² the MTH process may be integrated into carbon neutral synthetic routes to fuels and commercial chemicals.³ However, despite extensive research efforts, the underlying mechanism of the MTH process remains uncertain.⁴ The initial formation of the first C–C bond is especially challenging to describe due to the low concentrations of initiating species and the short time-scales of the constituent elementary steps. Further complexity is added by the variation in catalytic activity induced by the underlying structure and chemical composition of the catalyst, which can affect turn-over frequencies, product distributions, and deactivation rates.^5,6 Over the last two decades, the formation of the first C–C bond has been proposed via: indirect mechanisms, where aromatic impurities within the feedstock undergo methylation followed by olefin elimination; and direct mechanisms, where primary olefins are formed from methanol through a range of precursors, including surface methoxy species (SMS),^7,8 trimethyloxonium (TMO),^9,10 carbene,^11,12 dimethyl ether (DME),^7,13 and CO.^14,15 Consensus oscillates between the two categories of mechanistic routes to the initial C–C bond, with recent spectroscopic and theoretical evidence favouring direct mechanisms.⁸

DME has been studied as a key intermediate for the formation of C–C bonds in direct MTH mechanisms, but its importance is debated by time-resolved FTIR (Fourier-transformed infrared) spectroscopy studies.^8,12 These studies demonstrate that CH₃OH methylates the Brønsted acid site before characteristic DME signals appear. This reaction is followed by the formation of ethene via the coupling of surface methoxy species (SMS) through a carbene-like intermediate. However, the coupling of SMS requires that two framework CH₃ groups are brought into close contact. Minova et al.⁸ suggest that the SMS groups may shuttle across the zeolite by hopping between adjacent framework oxygen sites, but density functional theory (DFT) studies show migration of SMS from the Brønsted acid site to be highly endothermic.¹¹ Alternatively, an intermediate, short-lived methylation agent may allow an activated methyl group to shuttle across the catalyst framework to nearby SMS. Growing evidence supports trimethyl oxonium (TMO) as a candidate methylation agent, which has been detected by solid state ¹³C NMR spectroscopy.¹⁶ Computational studies also proposed that TMO may form the first C–C bond with comparable free energy barriers to methylation via SMS.¹⁰

The chemical composition of the zeolitic framework may also affect the formation of C₂ species in the induction period. Experimental studies show that aluminosilicate (zeolite) and aluminophosphate (ALPO) frameworks display different product selectivity and catalyst lifetimes.^5,6 For example, the higher acidic strength of the zeolite H-SSZ-13 compared to the ALPO H-SAPO-34 leads to a higher initial turn-over frequency in the MTH process, but the catalyst lifetime is reduced by increasing the rate that inactive polycyclic species are formed.⁶ Furthermore, DFT studies have demonstrated that the activation barriers for MTH initiation reactions are on average 19 kJ mol⁻¹ higher on H-SAPO-34 compared to the isostructural H-SSZ-13.¹⁷

To complement previous research,^11,13,17,18 additional DFT studies would be valuable to understand the effect of framework composition on the formation of SMS. In the present work, we will analyse SMS formation reactions proceeding from CH₃OH, DME, and TMO in a common framework to help achieve a consensus on the relative reactivity of each precursor. Here, we systematically assess the different routes to form the SMS in MTH reactions using ab initio simulations (Fig. 1). Five potential mechanisms for forming SMS from methanol are considered and grouped into routes proceeding from: methanol (reactions (I) and (II)); DME (reactions (III)–(V)); and TMO (reactions (VI) and (VII)). Activation barriers and reaction energies are calculated using periodic DFT simulations and QM/MM embedded-cluster models. We compare reactivity on the chabazite (CHA)-based, isostructural H-SSZ-13 aluminosilicate and H-SAPO-34 aluminophosphate frameworks, building on the work of Plessow et al.¹⁷ In addition, we study the mobility of SMS species on the H-SSZ-13 zeolite through migration reactions between adjacent framework oxygens. Comparison with experimental results provides further evidence for the most likely initial C₂ precursor, informing how the induction period may be guided towards specific product formation.

Fig. 1 Pathways to SMS formation on zeolitic and aluminophosphate frameworks through direct or indirect routes (I)–(VII), and the migration of SMS from an oxygen site adjacent to the Brønsted acid site to a neighbouring framework oxygen (VIII).

2 Methodology

2.1 Systems

The modeled reactions were performed on the H-SSZ-13 and H-SAPO-34 microporous frameworks. Both frameworks were generated from the CHA framework structure, which has 108 atoms in its unit cell (a = 13.675 Å, b = 13.675 Å, c = 14.767 Å, α = 90°, β = 90°, γ = 120°).¹⁹ The CHA framework has one distinct tetrahedral (T-)site, with 4 non-equivalent oxygen atoms that haven been labelled O1 to O4. The silicate framework is composed of repeating [SiO₄]⁴⁻ units, and the aluminophosphate is composed of alternating and repeating [AlO₄]⁵⁻/[PO₄]³⁻ units. One Brønsted acid site was introduced per unit cell of the H-SSZ-13 and H-SAPO-34 frameworks by replacing the T-site atom (Si or P, respectively) with a species bearing one less atomic charge (Al or Si, respectively); charge neutrality was achieved by placing a charge compensating proton on an adjacent oxygen. The O1 site is selected as the main reaction site, which is generally the most stable site for the proton in both H-SSZ-13 and H-SAPO-34 frameworks.²⁰

2.2 DFT parameters

Energies, forces, and stresses were calculated with the FHI-aims density functional theory (DFT) all-electron full-potential software package (Version: 231118).²¹ Periodic DFT calculations were performed with the PBE exchange–correlation density functional.²² Non-bonding dispersion interactions were incorporated with a non-local many-body-dispersion (MBD-NL) correction parameterised for the PBE density functional (PBE-MBD).^23,24 A converged Brillouin-zone²⁵ sampling was achieved with a (4 × 4 × 4) k-point grid. An FHI-aims default “light” basis²¹ was used for geometry relaxation tasks; total energy calculations were subsequently performed with a “tight” basis set for the converged structures to reduce errors driven by basis incompleteness. Systems were evaluated with a closed-shell (spin-paired) electronic structure, and relativistic corrections are applied to all atoms with the scalar atomic ZORA (zero-order regular approximation) scheme.^21,26 Additional numerical convergence parameters used in the SCF procedures are included in the ESI† (Section S1).

2.3 Geometry relaxation and nudged elastic band

Geometry optimisation was performed with a BFGS algorithm implemented in the Atomic Simulation Environment (ASE) Python package (Version: 3.23).²⁷ The total energy is minimised with respect to the lattice vectors and atomic positions simultaneously using the deformation gradient approach of Tadmor et al.²⁸ Convergence for the geometry optimisation was achieved when the maximum atomic force on all atoms was below 0.01 eV Å⁻¹.

Nudged elastic band (NEB) simulations were performed with the machine learning accelerated NEB (ML-NEB)²⁹ algorithm implemented in the Catlearn Python package.³⁰ The reaction pathway between reactants and products was represented with 6 interpolating images, i.e., 8 images in total. The initial path is generated using the image dependent pair potential (IDPP) interpolation³¹ between reactants and products. The NEB convergence criteria were set to 0.05 eV Å⁻¹ and 0.03 eV for the maximum atomic forces and the uncertainty of the energy, respectively. The highest energy image is then used as the starting point for a dimer method calculation, which further optimises the transition state structure towards a saddle point with a tighter convergence threshold (0.02 eV Å⁻¹). Vibrational frequency calculations were performed to ensure that the reactant and product structures were optimised to a local minimum (zero imaginary frequencies) and the transition state structure converged to a saddle point (one imaginary frequency).

The large number of degrees of freedom in the zeolite model significantly hindered the convergence for transition state calculations; therefore, atoms beyond the fourth nearest neighbour of the adsorption site were constrained for transition state calculations. To ensure the interpolated internal coordinates of the NEB reaction pathway are consistently optimized, the same constraints were applied to geometry optimisations of the reactant and product structures. The constraints introduce a maximum calculated error of ±1 kJ mol⁻¹ to the activation energies and ±2 kJ mol⁻¹ to the reaction energies (Section S2, ESI†).

2.4 QM/MM embedding simulations

Hybrid QM/MM (quantum mechanics/molecular mechanics) simulations were performed on converged ground state and transition state structures with higher-level hybrid-DFT using the Py-ChemShell software (Version: 23.0.0).³² Calculations were performed using the embedded-cluster model, where a small region of chemical interest is described at the high-accuracy QM level of theory, while the extended environment is evaluated with a computationally inexpensive classical force-field. In addition, the total energy of the embedded-cluster is evaluated under open-boundary conditions. QM/MM embedded-cluster simulations therefore measure energetic properties at the zero concentration limit, unlike periodic DFT where the simulation concentration depends on the dimensions of the simulation cell.

In the QM/MM embedding workflow, the QM calculations were performed with FHI-aims²¹ and MM calculations with the GULP software package (Version: 5.2.0).³³ The use of FHI-aims is advantageous because an identical numerical framework is applied for periodic and cluster calculations; energetics derived from periodic and aperiodic DFT should therefore be commensurate. Single-point energy evaluations were performed with the same numerical settings specified in Section 2.2 and the hybrid M06-2X^34,35 exchange–correlation density functional, which accurately predicts both adsorption and activation energies for reactions in zeolites.³⁶ The MM region was represented with a classical force-field derived from Hill and Sauer³⁷ with fixed charges,³⁸ which is available for zeolites only.

For each reaction on the zeolite H-SSZ-13, QM/MM simulations were performed using the optimised structures of the reactant, product, and transition state from periodic simulations. The atomic positions of the molecular species in the periodic unit cell (up to the fifth nearest neighbour) were transferred into an expanded, pristine CHA cluster. The QM region was defined using the delta-cluster method of Migues et al.,³⁹ where atoms described at the QM level are defined by overlapping, atom-centered spheres of radius δ from the constituent atoms of the adsorbate and Brønsted acid T-site. Converged energetics with an error bar of <4 kJ mol⁻¹ was achieved with δ = 11a₀ (see Section S3, ESI†).

3 Results and discussion

3.1 Formation of SMS on H-SSZ-13

The two mechanisms for the direct formation of SMS from methanol (reactions (I), (IIa), and (IIb)) are considered in Fig. 2, with the comparison and discussion herein focused on H-SSZ-13 to simplify the comparisons between reactions. The methylation of the zeolite framework from adsorbed methanol (reaction (I)) has a high barrier (255 kJ mol⁻¹) and is endothermic (63 kJ mol⁻¹). Alternatively, SMS may be formed by the protonation and rotation of methanol at the Brønsted acid site (BAS), followed by the methylation of the deprotonated BAS (reaction (IIb)). Although the formation of the reactant complex for reaction (IIb) through reaction (IIa) is endothermic (108 kJ mol⁻¹), the barrier to framework methylation through reaction (IIb) (37 kJ mol⁻¹) is significantly lower than in reaction (I). FTIR experiments have shown that methanol loading correlates with higher concentrations of protonated methanol on aluminosilicate zeolites,^40,41 and the low barriers for framework methylation from [CH₃OH₂]⁺ indicates that reactions (IIa) and (IIb) are a more feasible route towards forming SMS at these higher loadings. The large difference in activation barriers between reactions (I) and (II) are explored further in Section 3.3.

Fig. 2 Energy profile for SMS formation proceeding via CH₃OH on H-SSZ-13. The energies are shown for the reactants (R), products (P), and transition states (TS) of reactions (I) and (II). The route with the highest barrier is shown in orange, and the lowest barrier in blue. Energies are referenced to the isolated Brønsted acid site of H-SSZ-13. All calculations were performed with PBE-MBD. The annotated transition states structures are shown for reactions (I) and (IIb).

Previous DFT studies have shown that SMS is formed favourably from [CH₃OH₂]⁺ compared to CH₃OH. For the latter reaction, Nastase et al.¹¹ predict that the barrier for SMS formation from adsorbed CH₃OH is large (225 kJ mol⁻¹) in H-ZSM-5, which agrees with our calculated barrier for reaction (I) in H-SSZ-13. Furthermore, Di Iorio et al.⁴² predict that the barrier for SMS formation from [CH₃OH₂]⁺ is 129 kJ mol⁻¹,⁴² which compares well to the combined barrier calculated for CH₃OH protonation, [CH₃OH₂]⁺ rotation, and framework methylation presented in Fig. 2 (145 kJ mol⁻¹).

The activation barrier for the protonation of CH₃OH at the Brønsted acid site (reaction (IIa)) is omitted from Fig. 2. This is because the proton from the [CH₃OH₂]⁺ product transfers to the Brønsted acid site upon geometry optimisation, which suggests a near-barrierless reaction. Using a rigorous path integral molecular dynamics approach, which takes nuclear quantum effects into account, Hunt et al.⁴³ calculated the free energy barrier for proton transfer between [H₃O]⁺ and CH₃OH at the MP2 level of theory as 33 kJ mol⁻¹ at T = 50 K. This indicates that the barrier for proton transfer is smaller than the other steps measured in Fig. 2, and therefore it is considered unlikely that the protonation of CH₃OH is a significant barrier to SMS formation.

SMS may also be formed indirectly from dimethyl ether (DME) via reactions (III)–(V) (Fig. 1). Two routes to DME formation from methanol are commonly considered in the literature (Fig. 3):^42,44 the associative mechanism, where two methanol molecules undergo a dehydration reaction at the Brønsted acid site (reaction (III)); and the dissociative mechanism, where CH₃OH is methylated by SMS to form DME through dissociative pathway A (reversed reaction (V)), or protonated DME through dissociative pathway B (reversed reactions (IVb) and (IVa)). For indirect SMS formation routes to compete with direct routes, DME must be formed directly from CH₃OH rather than SMS produced from CH₃OH or [CH₃OH₂⁺]. Therefore, the associative mechanism of DME formation (reaction (III)) must be competitive with both reactions (I) and (II).

Fig. 3 Energy profiles for the formation of DME in H-SSZ-13 from CH₃OH through the associative pathway (blue), or from SMS through the dissociative pathway A (purple), and dissociative pathway B (orange). Reversed reactions (IVa), (IVb), and (V) are denoted as (rIVa), (rIVb), and (rV) respectively. The energies are shown for the reactants (R), products (P), and transition states (TS). Energies are referenced to the isolated Brønsted acid site of H-SSZ-13. All calculations were performed with PBE-MBD. The annotated transition states structures are shown for reactions (III), (IV), and (V).

The associative pathway (reaction (III)) for the formation of DME has previously been modeled in a variety of zeolites with various density functionals (Table 1). The large variation in the calculated reaction barriers (51 kJ mol⁻¹ to 151 kJ mol⁻¹) reflects the different methanol orientations used to model the pre-reaction complex. For example, the B97-D3 results of Nastase et al.¹³ were calculated using an unfavourable frontside-S_N2 mechanism, where the shuttling methyl group is less stable. In contrast, our simulations consider the more favourable backside S_N2 mechanism, which has an activation barrier that is 89 kJ mol⁻¹ lower. Further analysis of the transition state geometry identifies that the distances between the shuttling methyl carbon and the donor/acceptor oxygens are significantly reduced on H-SSZ-13 compared to the same reaction on H-ZSM-5 (donor: 1.93 Å vs. 2.13 Å and acceptor: 2.08 Å vs. 2.35 Å). Crossley-Lewis et al.⁴⁴ calculate a large activation barrier for the backside-S_N2 mechanism (214 kJ mol⁻¹) on H-ZSM-5, which contrasts with the lower barriers calculated on H-SSZ-13 in the present work (62 kJ mol⁻¹) and by Di Iorio et al. (51 kJ mol⁻¹). As shown by the large variation in activation barriers in different zeolites, the relative stability of the transition state for reaction (III) is sensitive to framework topology and/or the relative positioning of the reactants.

Table 1 Activation energies (ΔE^‡) for the formation of DME through the associative mechanism (Fig. 1, reaction (III)). The activation energies are calculated relative to the total energy of the adsorbed CH₃OH molecule. The activation barriers calculated in this work are shown in bold

Zeolite	Model	Density functional	ΔE^‡/kJ mol^−1
H-SSZ-13	Periodic	PBE-MBD	62
H-SSZ-13^42	Periodic	PBE-D3	51
H-ZSM-5^44	Periodic	PBE-D3	214
H-ZSM-22^45	Periodic	RPBE	112
H-SSZ-13 (217 atom QM cluster)	QM/MM cluster	PBE-MBD	64
H-SSZ-13 (217 atom QM cluster)	QM/MM cluster	M06-2X	73
H-ZSM-5 (74 atom QM cluster)^13	QM/MM cluster	B97-D3	151

---

The dissociative mechanisms to DME (dissociative pathways A and B) are considered through reversal of reactions (IVa), (IVb) and (V) (Fig. 3). The barrier to form protonated DME from methanol and SMS (reaction (IVb), reversed) is significantly lower than the barrier to form non-protonated DME (reaction (V), reversed) (65 kJ mol⁻¹vs. 176 kJ mol⁻¹, respectively). Although the reverse of reaction (V) forms the more thermodynamically stable product compared to reversed reaction (IVb) (−72 kJ mol⁻¹vs. +20 kJ mol⁻¹), the reversed reaction (IVa) reaction leads to a thermodynamically favoured product from protonated DME (−88 kJ mol⁻¹). Overall, we show that the associative mechanism is moderately favoured, as reaction (III) has an activation barrier that is 3 kJ mol⁻¹ smaller than the barrier to form protonated DME from CH₃OH and surface SMS (reaction (IVb), reversed). This conclusion is supported by experimental work, which shows that protonated DME is formed after the first signals of DME.⁸

After the formation of DME, SMS may be formed through the degradation of protonated DME (reaction (IVa) and (IVb)) or DME (reaction (V)), as shown in Fig. 4. Reaction (V) is kinetically disfavoured in H-SSZ-13 due to a large activation barrier (248 kJ mol⁻¹); however, framework methylation via protonated DME (reaction (IVb)) proceeds with a far lower activation barrier (45 kJ mol⁻¹), which agrees with the relatively low barriers calculated for direct SMS formation from protonated methanol (reaction (IIb)) compared to methanol (reaction (I)). Furthermore, the concerted protonation and rotation of adsorbed DME through reaction (VIa) is less endothermic than it is for adsorbed CH₃OH through reaction (IIa) (+88 kJ mol⁻¹ and +108 kJ mol⁻¹, respectively), meaning that the total barrier for SMS formation from protonated DME is 14 kJ mol⁻¹ lower than SMS formation from protonated methanol.

Fig. 4 Energy profile for SMS formation reactions (III), (IV) and (V) proceeding via DME on H-SSZ-13. The energies are shown for the reactants (R), products (P), and transition states (TS). The route with the highest activation barrier is shown in orange, and the lowest barrier in blue. Energies are referenced to the isolated Brønsted acid site of H-SSZ-13. All calculations were performed with PBE-MBD.

Similar to DME, TMO may form either through the methylation of DME via adsorbed methanol (reaction (VI)), or through the methylation of DME by framework SMS (reversed reaction (VII)), as included in Fig. 5. Compared to the methylation of DME by SMS (reverse reaction (VII)), the coupling of DME and CH₃OH at the Brønsted acid site (reaction (VI)) has a larger barrier (71 kJ mol⁻¹vs. 94 kJ mol⁻¹, respectively) and less favourable thermodynamics (+23 kJ mol⁻¹vs. +67 kJ mol⁻¹, respectively). Therefore, we conclude that TMO is more likely to form from DME and SMS already present in the zeolite. In addition, the low reaction barrier for the reverse of reaction (VI) (27 kJ mol⁻¹) means that adsorbed TMO may degrade to DME and methanol with relative ease. This degradation process may explain the low concentrations of TMO detected in zeolites under operating conditions.¹⁶

Fig. 5 Energy profile for SMS formation proceeding via TMO on H-SSZ-13. The energies are shown for the reactants (R), products (P), and transition states (TS) of reactions (III), (VI), and (VII). The barrier for the reverse reactions (VI) and (VII) (denoted (rVI) and (rVII)) are shown with purple arrows. Energies are referenced to the isolated Brønsted acid site of H-SSZ-13. All calculations were performed with PBE-MBD. The annotated transition state structures are shown for reactions (VI) and (VII).

Framework methylation via TMO (reaction (VII)) occurs with low activation barriers (48 kJ mol⁻¹), indicating a feasible three-step route to framework methylation proceeding via DME formation through an associative mechanism, followed by TMO formation with CH₃OH, after which TMO may degrade to DME and SMS. In addition, the formation energy of the pre-reaction complex for reaction (VII) (50 kJ mol⁻¹) is comparable to the combined protonation and rotation energy for DME through reaction (IVa) (88 kJ mol⁻¹). However, the barrier for the decomposition of TMO to CH₃OH and DME (reverse reaction (VI)) is 23 kJ mol⁻¹ lower than the barrier to form the pre-reaction complex of the TMO methylation reaction (reaction (VII)). This means that TMO is more likely to reform CH₃OH and DME than proceed to the final SMS product.

Overall, SMS are found to form on H-SSZ-13 zeolites from the following species with barriers ordered from lowest to highest: [CH₃OH₂]⁺ (37 kJ mol⁻¹) < protonated DME (45 kJ mol⁻¹) < TMO (48 kJ mol⁻¹) < DME (248 kJ mol⁻¹) < CH₃OH (255 kJ mol⁻¹). Because DME and TMO may be formed with relatively low activation barriers (62 kJ mol⁻¹ and 94 kJ mol⁻¹, respectively), they may be considered feasible precursors to SMS via indirect formation routes. However, the enthalpic cost of protonation and rotation are comparable to the barriers of framework methylation for CH₃OH and DME. Therefore, the formation of protonated CH₃OH and DME are most likely the rate determining steps to SMS formation via the reactions presented. As [CH₃OH₂]⁺ has been observed experimentally at room temperature in tandem with the formation of SMS,⁴¹ and reaction (IIb) is the most energetically favourable SMS route in our work, we consider the [CH₃OH₂]⁺ the most likely precursor to framework methylation.

3.2 Migration of SMS on H-SSZ-13

Minova et al. propose that SMS migration to different O-sites in the zeolite framework is necessary for C₂ formation.⁸ To test the feasibility of SMS migration in H-SSZ-13, reaction (VIII) is modeled by considering the transfer of SMS from its original O-site to neighbouring O-sites in the 8-membered ring (8-MR), as shown in Fig. 6. These reactions are modeled as SMS migrations on the oxygen atoms at same T-site (0th-nearest neighour (NN0)), and to T-sites that are 1 and 2 nearest neighbours (NN1 and NN2) away from the original T-site. SMS migration occurring in a H-SSZ-13 unit cell with multiple Al-substituted T-sites in the 8-MR is also considered, motivated by the relatively low Si/Al ratio (~5) for H-SSZ-13, equivalent to 6 Al atoms per CHA unit cell, when synthesized with the method of Fickel et al.⁴⁶ In comparison, the H-SSZ-13 unit cell with one Al substituted T-site has a high Si/Al ratio (35), largely precluding potential SMS migrations between adjacent Brønsted acid sites. SMS migration between paired Al T-sites (O–Al–O–Si–O–Al–O) is therefore also considered; although the concentration of paired Al sites is typically very low under synthesis methods,⁴⁷ synthetic methods that combine the N,N,N-trimethyl-1-adamantyl ammonium (TMAda⁺) structure directing agent with Na⁺ during the crystallization step induce higher concentrations of paired Al sites under a random distribution.⁴⁸ Comparing between the migration of SMS to a paired Al T-site or a close Si T-sites may then reveal mechanistic differences between these frameworks prepared using different synthetic methods.

Fig. 6 Migration of SMS on H-SSZ-13 from the –O–Al–O– site in the 8-MR to: the same T-site (NN0); a T-site one nearest neighbour (NN1) site away; and a T-site two nearest neighbour (NN2) sites away, as labelled. The H-SSZ-13 framework with isolated (1Al) and paired (2Al) Al substituted T-sites are represented in the top and bottom two rows, respectively. The reaction denoted as 2Al NN2 CH₂ proceeds through a carbene-like intermediate. The activation energies for the corresponding reactions are shown in Table 2.

Table 2 shows the reaction energies and activation barriers for SMS migration between different O-sites. The barrier for SMS migration between O-sites on the same T-site (1Al–NN0) is 200 kJ mol⁻¹, but the reaction is only mildly endothermic (10 kJ mol⁻¹). SMS migrations from an –O–Al–O– site to the nearest and second nearest –O–Si–O– site (1Al–NN1 and 1Al–NN2) also have large activation barriers (238 kJ mol⁻¹ and 211 kJ mol⁻¹, respectively). However, SMS formed at –O–Si–O– site are substantially less stable than at the –O–Al–O– site, leading to either a highly endothermic reaction for NN1 (+90 kJ mol⁻¹), or an unstable geometry for NN2 where CH₃ freely dissociates.

Table 2 Reaction energies (ΔE) and activation energies (ΔE^‡) for SMS migration in H-SSZ-13, as shown in Fig. 6. SMS migration is considered between O-sites on the same T-site (NN0), and T-sites one (NN1) and two (NN2) nearest neighbours from the original T-site. The H-SSZ-13 framework is modeled with one or two aluminium substituted T-sites (1Al and 2Al, respectively). The 2Al–NN2 (CH₂) reaction proceeds through a carbene-like transition state and all other reactions through a CH₃ transition state. X denotes an unstable product geometry for which a reaction energy could not be calculated

SMS migration type	ΔE kJ mol^−1	ΔE^‡ kJ mol^−1
1Al–NN0	10	206
1Al–NN1	94	238
1Al–NN2	X	211
2Al–NN1	10	197
2Al–NN2	10	182
2Al–NN2 (CH2)	1	251

---

The reaction energy for the SMS migration to the first nearest neighbour –O–Si–O– site agrees with endothermic values for SMS migration in the H-ZSM-5 zeolite calculated with the B97-D density functional (126 kJ mol⁻¹).¹¹ Nastase et al. also show that the O–CH₃ bond of the methoxy group is weaker on a framework oxygen than at the Brønsted acid site in H-ZSM-5 (−277 kJ mol⁻¹vs. −358 kJ mol⁻¹, respectively).¹¹ Overall, the calculated activation barriers for SMS migration to the –O–Si–O– sites (211–238 kJ mol⁻¹) are similar to the barrier for the direct formation of SMS from CH₃OH (255 kJ mol⁻¹). As the kinetic barrier for SMS migration is large, the methyl species is concluded to be immobile on the H-SSZ-13 zeolite framework for isolated Al sites.

The reactions barriers for SMS migration between paired Al sites are slightly smaller than the corresponding barriers to migration between the –O–Al–O– and –O–Si–O– sites. The activation barrier for paired Al sites is 15 kJ mol⁻¹ lower for the second (NN2) relative to the first (NN1) nearest neighbour T-sites, which is smaller than for the isolated Al site (27 kJ mol⁻¹). Furthermore, the migration between the two –O–Al–O– sites is only slightly endothermic (10 kJ mol⁻¹), reflecting the similarity in binding strength for SMS with each –O–Al–O– site. The reduced kinetic barriers for SMS migration with a paired Al site, compared to a Si T-site, indicates that higher Al concentrations may facilitate the movement of SMS through the framework.

Introducing a second Al to the framework allows a further proposed mechanism for SMS migration through a carbene-like intermediate (Fig. 6, bottom). In this reaction, the shuttling SMS may simultaneously protonate and deprotonate two Brønsted acid sites in close proximity. We therefore modeled this reaction where, in a single concerted step, SMS simultaneously donates and accepts protons to and from the H-SSZ-13 framework. Compared to the migration of the CH₃ group between paired Al sites (2Al–NN2), the carbene-based mechanism is thermodynamically favoured by 9 kJ mol⁻¹, but the kinetic barrier is 69 kJ mol⁻¹ higher (Table 2). A contributing factor to the increased reaction barrier of the carbene-mechanism is strain introduced to the 8-MR by the transition state, which is demonstrated by the reduction of the 8-MR diameter by 0.58 Å for the transition state compared to the reactant. In contrast, for the corresponding SMS migration reaction with a CH₃ intermediate (2Al–NN2), the diameter of the 8-MR for the transition state only reduces by 0.22 Å with respect to the reactant. We consider the carbene-like reaction as unlikely for the current framework due to the relatively large barrier; however, a more full exploration of paired Al sites at different positions and topological frameworks may support the formation of the carbene-like transition state with reduced strain.

3.3 S_N2 frontside and backside attack and frontside attack on H-SSZ-13

Reactions (I)–(VIII) are bimolecular nucleophilic substitution (S_N2) reactions with activation energies (ΔE^‡) that depend on the attack angle of the nucleophile. This attack angle may align with the backside or the frontside of the bond between the alkyl carbon and the leaving group (Fig. 7), which are labelled herein as S_N2-b and S_N2-f reactions, respectively. In agreement with the activation strain model,^49,50 our results show that ΔE^‡ is significantly larger for S_N2-f pathways than S_N2-b pathways (Table 3). For example, the formation of SMS through the frontside attack of CH₃OH (reaction (I)) is 218 kJ mol⁻¹ higher than the backside attack of [CH₃OH₂]⁺ by the framework oxygen (reaction (IIb)).

Fig. 7 Left: Backside attack S_N2 (S_N2-b) and right: frontside attack S_N2 (S_N2-f) reactions demonstrated for DME formation via reactions (III) and (V), respectively.

Table 3 Activation energies (ΔE^‡) for S_N2 backside (S_N2-b) and frontside (S_N2-f) attack shown for comparison. Hirshfeld charges are given for the carbon atom of the methylating agent, as calculated for the transition state of reactions (I)–(VIII)

Reaction type	Reaction	ΔE^‡/kJ mol^−1	C Hirshfeld charge/|e|
SN2-b	II	37	0.03
	III	62	0.03
	IV	45	0.03
	VI	94	0.03
	VII	48	0.05
SN2-f	I	255	0.11
	V	176	0.10
	VIII	238	0.13

---

In S_N2-b mechanisms, orbital overlap is maximized between the methyl group and the nucleophilic and leaving groups (Fig. 7). The increased orbital overlap reduces charge separation in the transition state complex by increasing the electron donation to the alkyl carbon atom. The reduced charge on the alkyl carbon stabilises the S_N2-b transition state relative to the S_N2-f pathways. Table 3 shows the Hirshfeld charge of the alkyl carbon in the transition state for reactions (I)–(VIII). The charge of the methylating carbon is significantly lower for the S_N2-b mechanisms (0.02–0.05 |e|) than S_N2-f mechanisms (0.10–0.13 |e|), with reduced charge in the alkyl carbon correlating with the more stable transition state complexes. Therefore, the distinction between S_N2-f and S_N2-b pathways is a useful heuristic for identifying favourable pathways for reactions in zeolite frameworks.

3.4 Comparison of H-SSZ-13 and H-SAPO-34

H-SAPO-34 is the isoelectronic and isostructural aluminophosphate analogue to the aluminosilicate H-SSZ-13 zeolite. DFT studies show that the deprotonation energy of H-SAPO-34 (1265 kJ mol⁻¹) is greater than that of H-SSZ-13 (1222 kJ mol⁻¹), indicating that aluminophosphate frameworks are less acidic than the corresponding aluminosilicate frameworks.⁶ The weaker acidity of H-SAPO-34 increases the activation barriers of dehydrogenation reactions,⁵¹ which subsequently reduces the rate of paraffin formation and the deactivation rate of the catalyst.⁶

The differences in reaction energies (ΔE) and activation energies (ΔE^‡) for reactions (I)–(VIII) in H-SSZ-13 and H-SAPO-34 are presented in Fig. 8. The corresponding reaction profiles for the SMS formation reactions from MeOH, DME, and TMO are presented in the ESI† (Sections S4.1, S4.2, and S4.3, respectively). The activation energies are 8–22 kJ mol⁻¹ higher in H-SAPO-34 for most reactions that form H₂O (i.e., reactions where the proton transfers from the Brønsted acid site to the reactant). Our results agree with previous theoretical studies, where the activation barriers increased by 10–30 kJ mol⁻¹ for MTH initiation reactions on H-SAPO-34 compared to H-SSZ-13.^17,51 However, ΔE^‡ in our results increases by 52 kJ mol⁻¹ for reaction (III), which reflects the sensitivity of the associative DME formation mechanism to small changes in framework composition and reactant orientation.

Fig. 8 (A) Reaction energies (ΔE) and (B) activation energies (ΔE^‡) for reactions (I)–(VIII) (Fig. 1) in the CHA-based H-SSZ-13 zeolite and ALPO H-SAPO-34 framework. The calculations were performed with PBE-MBD using periodic DFT.

The ΔE^‡ for SMS formation reactions via oxonium-based reactants (reactions (IIb), (IVb), and (VII)) are lower on H-SAPO-34 than on H-SSZ-13 (Section 3.1). The barriers for these reactions are all small (<48 kJ mol⁻¹), and therefore these steps are unlikely to be rate limiting in the induction period for the MTH reaction. Furthermore, the relative magnitudes of the activation barriers are preserved between H-SAPO-34 and H-SSZ-13 (Section 3.1), which implies that there is no change between the frameworks in the preferred pathway for the formation of SMS. The small changes in reaction energetics show that the topological effects of the CHA lattice are consistent for the aluminosilicate and aluminophosphate frameworks.

For reactions (I), (V), and (VI), the overall differences in ΔE between H-SAPO-34 and H-SSZ-13 are relatively small (<15 kJ mol⁻¹), and there is no clear pattern of reaction preference towards either the H-SAPO-34 or H-SSZ-13 catalyst. However, for framework methylation via oxonium-based reactants (reactions (IIb), and (IVb)), the SMS product is more favourably formed on H-SAPO-34 compared to H-SSZ-13, with ΔE being 20 and 21 kJ mol⁻¹ more exothermic, respectively. The greater stability of the methoxy product implies SMS binds more favourably to the aluminophosphate framework than the aluminosilicate framework.

The migration of SMS across the H-SAPO-34 framework is calculated to be unfavourable, with large ΔE (126 kJ mol⁻¹) and ΔE^‡ (260 kJ mol⁻¹) for the 1Al-NN1 migration reaction, similar to H-SSZ-13. We therefore conclude that SMS are unlikely to move between adjacent oxygen sites on the H-SAPO-34 framework.

3.5 Hybrid-DFT calculations for reactions in H-SSZ-13

GGA (generalised gradient approximation) density functionals, such as PBE, incur significant total energy errors from self-interaction and over-delocalisation.⁵² The self-interaction errors are especially large for molecules with non-equilibrium bond lengths,⁵³ which can lead to large errors in the activation barriers. Hybrid-DFT methods and post-HF approaches provide more accurate activation barriers as, unlike the semi-local GGA density functionals, they include non-local contributions to the total energy. However, the unfavourable scaling of hybrid-DFT with respect to system size limits its use to relatively smaller system sizes.

A significant advantage of the FHI-aims software package is the ability to model periodic and open boundary conditions with an identical numerical framework. Therefore, there are no arising energetic discrepancies from the atomic basis formalism between the application of an aperiodic QM/MM embedded-cluster and a periodic approach. Using this advantageous software infrastructure, the differences in reaction energies and activation barriers were calculated for: (i) periodic and cluster models with the PBE-MBD functional (ΔΔE_cluster), (ii) QM/MM embedded-cluster models using the PBE-MBD and M06-2X functionals (ΔΔE_func) and (iii) the combination of (i) and (ii) (ΔΔE_tot):

ΔΔE_cluster = ΔE_PBE-MBD(QM/MM) − ΔE_PBE-MBD(periodic), (1)

ΔΔE_func = ΔE_M06-2X(QM/MM) − ΔE_PBE-MBD(QM/MM), (2)

ΔΔE_tot = ΔE_M06-2X(QM/MM) − ΔE_PBE-MBD(periodic), (3)

where ΔΔE and ΔΔE^‡ are used herein to represent differences in reaction energies and activation barriers, respectively.

The difference in reaction energetics calculated with eqn (1)–(3) are summarised in Fig. 9, with the absolute values for reaction and activation energies compared in Table 4. The differences in calculated activation barriers, ΔΔE^‡, between the QM/MM embedded-cluster model (M06-2X) and periodic model (PBE-MBD) are predominantly due to the change of the density functional. The contribution of ΔΔE^‡_func to ΔΔE^‡_tot ranges from 8 to 44 kJ mol⁻¹; in comparison, ΔΔE^‡_cluster was smaller, ranging from 2–16 kJ mol⁻¹. For most reactions, a cancellation between ΔΔE^‡_cluster and ΔΔE^‡_func results in a decrease in ΔΔE^‡_tot. Nevertheless, ΔΔE^‡_tot exceeds 30 kJ mol⁻¹ for reactions (V), (VI), and (VIII). Overall, the accurate M06-2X hybrid density functional results in larger barriers, which highlights the importance of accurately incorporating non-local effects in DFT calculations. Goncalves et al. also show that GGA functionals with dispersion corrections underestimate the activation barriers for the initiation steps in the MTH process, with a mean absolute error (MAE) for PBE-D3 relative to MP2 for ΔE^‡ of 42 kJ mol⁻¹, whilst the error for M06-2X was smaller at 7 kJ mol⁻¹.³⁶

Fig. 9 The differences between (A) reaction energies (ΔΔE) and (B) activation barriers (ΔΔE^‡) as calculated with periodic and QM/MM embedded-cluster model approaches for reactions (I)–(VIII). The total differences (ΔΔE_tot) are divided into contributions from the change in density functional (ΔΔE_func) and the transfer from a periodic to an embedded-cluster geometry (ΔΔE_cluster). All comparisons are made for the H-SSZ-13 framework only.

Table 4 Reaction energies (ΔE), forward activation barriers (ΔE^‡_f), and backward activation barriers (ΔE^‡_b) for reactions (I)–(VIII), calculated with PBE-MBD using periodic DFT and with M06-2X using an aperiodic QM/MM embedded-cluster model

Reaction	ΔE periodic PBE-MBD	ΔE QM/MM M06-2X	ΔE^‡f periodic PBE-MBD	ΔE^‡f QM/MM M06-2X	ΔE^‡b periodic PBE-MBD	ΔE^‡b QM/MM M06-2X
I	63	43	255	274	192	231
IIb	−39	−63	37	46	76	109
III	−28	−58	62	73	89	131
IVb	−20	−47	45	55	65	102
V	72	61	248	280	176	219
VI	67	60	94	123	27	63
VII	−23	−43	48	57	71	100
VIII	94	116	238	273	144	157

---

Considering the overall reaction energies (Fig. 9A), the absolute values of ΔΔE_tot are between 7 and 30 kJ mol⁻¹, with a MAE of 20 kJ mol⁻¹. The MAE of the reaction energies is comparable to the MAE of the activation energies (20 kJ mol⁻¹), but the MAE of ΔΔE_func (9 kJ mol⁻¹) is smaller relative to ΔΔE^‡_func (27 kJ mol⁻¹). These results demonstrate that higher accuracy evaluation of non-local quantities are less important for the equilibrium structures of the reactants and products than for the non-equilibrium geometry of the transition states.⁵³ However, the average contribution of ΔΔE_cluster (17 kJ mol⁻¹) to ΔΔE_tot is higher than that of ΔΔE_func (9 kJ mol⁻¹). Furthermore, as ΔΔE_func and ΔΔE_cluster share the same sign, ΔΔE_tot does not benefit from a cancellation of errors, unlike ΔΔE^‡_tot. The increased contribution of ΔΔE_cluster is indicative of a finite cell size error, which is caused by self-interaction between the periodic replicas of the adsorbates. Additional discussion regarding the elimination of the finite cell size error in the QM/MM embedded cluster simulation, which impacts ΔΔE^‡_cluster and ΔΔE_cluster, may be found in ESI,† Section S5.

Overall, the QM/MM embedding simulations significantly improve upon the semi-local periodic simulations by allowing the use of more accurate hybrid-DFT levels of theory, and by eliminating artefacts such as periodic interactions that arise for small unit cells. The potentially large errors associated with transition state barriers especially necessitate the use of density functionals beyond GGAs to ensure accuracy. A workflow for performing QM/MM embedded-cluster simulations for aluminophosphate frameworks would be an important next step in performing accurate transition state calculations at tractable cost for the H-SAPO-34 framework, and is to be investigated in our future work.

4 Conclusions

The reaction and activation energies have been compared for a variety of reactions relating to SMS formation on H-SSZ-13 and H-SAPO-34. SMS formation was considered as proceeding directly from methanol and indirectly with TMO and DME intermediates. The rate determining steps for the formation of DME (62 kJ mol⁻¹ on H-SSZ-13) and TMO (94 kJ mol⁻¹ on H-SSZ-13) have relatively low kinetic barriers, which means indirect routes to SMS formation are energetically viable. Considering the low barriers of framework methylation for [CH₃OH₂]⁺ (37 kJ mol⁻¹), our work supports the reaction chronology of time-resolved FTIR studies,^8,12 where methanol first forms clusters around the Brønsted acid site⁴¹ from which a protonated methanol forms and then this leads to SMS formation on the zeolite framework.

The migration of SMS across the zeolitic framework was studied by modeling the transfer of SMS in the 8-MR of H-SSZ-13 from isolated and paired Al T-sites. The barriers to SMS migration are large (>180 kJ mol⁻¹), with migration to the first-nearest neighbour being kinetically disfavoured for both isolated (206 kJ mol⁻¹) and paired (197 kJ mol⁻¹) Al T-sites. An additional carbene-like mechanism was modeled for the transfer of SMS between paired T-sites, but the barrier (251 kJ mol⁻¹) was higher than SMS migration reactions with a CH₃ transition state. The transition state of this reaction was destabilised by an increase in ring strain, where accommodating the carbene-like transition state reduced the diameter of the 8-MR by 0.58 Å. We propose that other framework topologies and paired T-sites (i.e., on 6-MR and 4-MR structures) may lead to more favourable reaction energetics.

Our periodic DFT calculations were complemented by embedded-cluster QM/MM simulations using accurate hybrid density functionals.^32,34 The barriers calculated with the hybrid density functionals are consistently higher than the corresponding barriers calculated for a periodic unit cell, by an average of 27 kJ mol⁻¹. The embedded-cluster model effectively eliminates the interaction between periodic images, thus correctly modeling reactivity at the dilute concentration limit. In contrast, periodic simulations are demonstrated to measure reaction energetics at concentrations dependent on the size of the unit cell, and this must be considered in energetic analysis alongside experiment.

Comparison of aluminosilicates and aluminophosphates shows a consistent increase in activation barriers for H-SAPO-34 compared to H-SSZ-13, ranging from 8 kJ mol⁻¹ for direct SMS formation from methanol to 52 kJ mol⁻¹ for the associative formation of DME. The trend of higher activation energies for the aluminophosphate frameworks suggest slower formation of SMS and C₂ products, which may be rationalised by the reduced acidity of the H-SAPO-34 framework.^5,40,54 However, the unchanged relative ordering in activation barriers of each constituent reaction is predicted to lead to few changes in the route to C₂ product formation, in agreement with the work of Plessow et al.¹⁷ To calculate reaction barriers and energies for H-SAPO-34 accurately with hybrid-DFT methods and beyond, an analogous QM/MM embedding cluster model for aluminophosphate frameworks is needed. Towards this goal, parameterised force fields for aluminophophates and partitioning schemes that accurately describe bonding at the QM/MM interface are being pursued in continuing work.

Author contributions

Conceptualisation, G. B. and A. J. L.; metholodogy, G. B. and A. J. L.; software, G. B., O. v. V., and A. J. L.; investigation, G. B.; resources, G. B. and A. J. L.; writing – original draft, G. B.; writing – review & editing, G. B., O. v. V. and A. J. L.; visualisation, G. B.; supervision, A. J. L.; funding acquisition, G. B. and A. J. L.

Data availability

Data supporting the DFT simulations of this study are openly available from the NOMAD repository at DOI: https://doi.org/10.17172/NOMAD/2025.01.09-2.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

A. J. L., O. v. V., and G. B. acknowledge funding by the UKRI Future Leaders Fellowship program (MR/T018372/1, MR/Y034279/1). We acknowledge the support of the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via Welsh Government. Via our membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/X035859), this work used the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk).

Notes and references

1. P. Tian, Y. Wei, M. Ye and Z. Liu, ACS Catal., 2015, 5, 1922–1938.
2. N. Dahmen, E. Henrich, E. Dinjus and F. Weirich, Energy Sustainability Soc., 2012, 2, 3.
3. M. Hennig and M. Haase, Fuel Process. Technol., 2021, 216, 106776.
4. T. Omojola, A. J. Logsdail, A. C. V. Veen and S. A. F. Nastase, Phys. Chem. Chem. Phys., 2021, 23, 21437–21469.
5. P. Ferri, C. Li, R. Millán, J. Martínez-Triguero, M. Moliner, M. Boronat and A. Corma, Angew. Chem., 2020, 132, 19876–19883.
6. Z. Shi, M. Neurock and A. Bhan, ACS Catal., 2021, 11, 1222–1232.
7. T. Omojola, Chem. Eng. Sci., 2022, 251, 117424.
8. I. B. Minova, M. Bühl, S. K. Matam, C. R. A. Catlow, M. D. Frogley, G. Cinque, P. A. Wright and R. F. Howe, Catal. Sci. Technol., 2022, 12, 2289–2305.
9. X. Wu, W. Chen, S. Xu, S. Lin, T. Sun, A. Zheng, Y. Wei and Z. Liu, ACS Cent. Sci., 2021, 7, 681–687.
10. T. Sun, W. Chen, S. Xu, A. Zheng, X. Wu, S. Zeng, N. Wang, X. Meng, Y. Wei and Z. Liu, Chem, 2021, 7, 2415–2428.
11. S. A. F. Nastase, A. J. Logsdail and C. R. A. Catlow, Phys. Chem. Chem. Phys., 2021, 23, 17634–17644.
12. I. B. Minova, S. K. Matam, A. Greenaway, C. R. A. Catlow, M. D. Frogley, G. Cinque, P. A. Wright and R. F. Howe, ACS Catal., 2019, 9, 6564–6570.
13. S. A. F. Nastase, C. R. A. Catlow and A. J. Logsdail, Phys. Chem. Chem. Phys., 2021, 23, 2088–2096.
14. Y. Liu, S. Müller, D. Berger, J. Jelic, K. Reuter, M. Tonigold, M. Sanchez-Sanchez and J. A. Lercher, Angew. Chem., 2016, 128, 5817–5820.
15. A. D. Chowdhury, K. Houben, G. T. Whiting, M. Mokhtar, A. M. Asiri, S. A. Al-Thabaiti, S. N. Basahel, M. Baldus and B. M. Weckhuysen, Angew. Chem., Int. Ed., 2016, 55, 15840–15845.
16. X. Wu, S. Xu, W. Zhang, J. Huang, J. Li, B. Yu, Y. Wei and Z. Liu, Angew. Chem., Int. Ed., 2017, 56, 9039–9043.
17. P. N. Plessow and F. Studt, Catal. Lett., 2018, 148, 1246–1253.
18. S. A. F. Nastase, P. Cnudde, L. Vanduyfhuys, K. De Wispelaere, V. Van Speybroeck, C. R. A. Catlow and A. J. Logsdail, ACS Catal., 2020, 10, 8904–8915.
19. L. S. Dent and J. V. Smith, Nature, 1958, 181, 1794–1796.
20. Y. Jeanvoine, J. G. Ángyán, G. Kresse and J. Hafner, J. Phys. Chem. B, 1998, 102, 5573–5580.
21. V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Comput. Phys. Commun., 2009, 180, 2175–2196.
22. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
23. A. Tkatchenko, R. A. DiStasio, R. Car and M. Scheffler, Phys. Rev. Lett., 2012, 108, 236402.
24. A. Tkatchenko, A. Ambrosetti and R. A. DiStasio, J. Chem. Phys., 2013, 138, 074106.
25. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.
26. E. van Lenthe, E. J. Baerends and J. G. Snijders, J. Chem. Phys., 1994, 101, 9783–9792.
27. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, J. Phys.:Condens. Matter, 2017, 29, 273002.
28. E. B. Tadmor, G. S. Smith, N. Bernstein and E. Kaxiras, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 235–245.
29. J. A. Garrido Torres, P. C. Jennings, M. H. Hansen, J. R. Boes and T. Bligaard, Phys. Rev. Lett., 2019, 122, 156001.
30. M. H. Hansen, J. A. G. Torres, P. C. Jennings, Z. Wang, J. R. Boes, O. G. Mamun and T. Bligaard, An Atomistic Machine Learning Package for Surface Science and Catalysis, arXiv, 2019, preprint, arXiv:1904.00904 [physics] DOI:10.48550/arXiv.1904.00904.
31. S. Smidstrup, A. Pedersen, K. Stokbro and H. Jónsson, J. Chem. Phys., 2014, 140, 214106.
32. Y. Lu, M. R. Farrow, P. Fayon, A. J. Logsdail, A. A. Sokol, C. R. A. Catlow, P. Sherwood and T. W. Keal, J. Chem. Theory Comput., 2019, 15, 1317–1328.
33. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
34. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
35. M. Walker, A. J. A. Harvey, A. Sen and C. E. H. Dessent, J. Phys. Chem. A, 2013, 117, 12590–12600.
36. T. J. Goncalves, P. N. Plessow and F. Studt, ChemCatChem, 2019, 11, 4368–4376.
37. J. R. Hill and J. Sauer, J. Phys. Chem., 1994, 98, 1238–1244.
38. P. Sherwood, A. H. D. Vries, S. J. Collins, S. P. Greatbanks, N. A. Burton, M. A. Vincent and I. H. Hillier, Faraday Discuss., 1997, 106, 79–92.
39. A. N. Migues, A. Muskat, S. M. Auerbach, W. Sherman and S. Vaitheeswaran, ACS Catal., 2015, 5, 2859–2865.
40. S. Bordiga, L. Regli, D. Cocina, C. Lamberti, M. Bjørgen and K. P. Lillerud, J. Phys. Chem. B, 2005, 109, 2779–2784.
41. S. K. Matam, S. A. F. Nastase, A. J. Logsdail and C. R. A. Catlow, Chem. Sci., 2020, 11, 6805–6814.
42. J. R. Di Iorio, A. J. Hoffman, C. T. Nimlos, S. Nystrom, D. Hibbitts and R. Gounder, J. Catal., 2019, 380, 161–177.
43. D. Hunt, L. Turi, D. Laria and K. Golobits, Isomerization-assisted Proton Transfers in MeOH-(H2O)2H+, 2025, https://chemrxiv.org/engage/chemrxiv/article-details/67f7b6c281d2151a02d5950a.
44. J. Crossley-Lewis, J. Dunn, I. F. Hickman, F. Jackson, G. J. Sunley, C. Buda, A. J. Mulholland and N. L. Allan, Phys. Chem. Chem. Phys., 2024, 26, 16693–16707.
45. P. G. Moses and J. K. Nørskov, ACS Catal., 2013, 3, 735–745.
46. D. W. Fickel and R. F. Lobo, J. Phys. Chem. C, 2010, 114, 1633–1640.
47. K. Mlekodaj, J. Dedecek, V. Pashkova, E. Tabor, P. Klein, M. Urbanova, R. Karcz, P. Sazama, S. R. Whittleton, H. M. Thomas, A. V. Fishchuk and S. Sklenak, J. Phys. Chem. C, 2019, 123, 7968–7987.
48. J. R. Di Iorio and R. Gounder, Chem. Mater., 2016, 28, 2236–2247.
49. A. P. Bento and F. M. Bickelhaupt, J. Org. Chem., 2008, 73, 7290–7299.
50. I. Fernández and F. M. Bickelhaupt, Chem. Soc. Rev., 2014, 43, 4953–4967.
51. M. DeLuca, C. Janes and D. Hibbitts, ACS Catal., 2020, 10, 4593–4607.
52. N. Mardirossian and M. Head-Gordon, Mol. Phys., 2017, 115, 2315–2372.
53. C. Shahi, P. Bhattarai, K. Wagle, B. Santra, S. Schwalbe, T. Hahn, J. Kortus, K. A. Jackson, J. E. Peralta, K. Trepte, S. Lehtola, N. K. Nepal, H. Myneni, B. Neupane, S. Adhikari, A. Ruzsinszky, Y. Yamamoto, T. Baruah, R. R. Zope and J. P. Perdew, J. Chem. Phys., 2019, 150, 174102.
54. L. Wu and E. J. M. Hensen, Catal. Today, 2014, 235, 160–168.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp00256g

---