Molecular Understanding of the Catalytic Consequence of Ketene Intermediates under Confinement

Neutral ketene is a crucial intermediate during zeolite carbonylation reactions. In this work, the roles of ketene and its derivates (viz., acylium ion and surface acetyl) associated with direct C–C bond coupling during the carbonylation reaction have been theoretically investigated under realistic reaction conditions and further validated by synchrotron radiation X-ray diffraction (SR-XRD) and Fourier transformed infrared (FT-IR) studies. It has been demonstrated that the zeolite confinement effect has significant influence on the formation, stability, and further transformation of ketene. Thus, the evolution and the role of reactive and inhibitive intermediates depend strongly on the framework structure and pore architecture of the zeolite catalysts. Inside side pockets of mordenite (MOR), rapid protonation of ketene occurs to form a metastable acylium ion exclusively, which is favorable toward methyl acetate (MA) and acetic acid (AcOH) formation. By contrast, in 12MR channels of MOR, a relatively longer lifetime was observed for ketene, which tends to accelerate deactivation of zeolite due to coke formation by the dimerization of ketene and further dissociation to diene and alkyne. Thus, we resolve, for the first time, a long-standing debate regarding the genuine role of ketene in zeolite catalysis. It is a paradigm to demonstrate the confinement effect on the formation, fate, and catalytic consequence of the active intermediates in zeolite catalysis.


Periodic DFT Caculations
The topologies of mordenite (MOR) were extracted from the official website of the International Zeolite Association (IZA). [S1] The cell parameters of SSZ-13 and MOR zeolites are (a = b = 13.67 Å, c = 14.77 Å) and (a = 18.26 Å, b = 20.53 Å, c = 15.08 Å), respectively. For the model of SSZ-13, one of silicon or phosphorus in T1-O2-T1 sites was replaced by aluminum or silicon to introduce the surface methoxy species (SMS) on O2, as shown in Figure S1. On the other hand, one of silicon in T4-O10-T4 (12MR) and T3-O8-T3 (8MR) were replaced by aluminum to introduce the SMS for MOR, and denoted as MOR-12MR and MOR-8MR, respectively. To account for zeolite framework flexibility and guest molecule dynamics in full at designated reaction temperatures, ab initio molecular dynamic (AIMD) simulations were performed by periodic density functional theory (DFT). This is carried out by first performing equilibrium AIMD simulations of reactant state species for 5 ps at the isothermal-isobaric (NPT) ensemble, in which the amount of substance (N), pressure (P) and temperature (T) are conserved, to relax the cell parameters and atom positions.
Next, the metadynamics (MTD) method was exploited to explore the occurrence of carbonylation in MOR-12MR as well as MOR-8MR during at least 35 ps AIMD simulations. Finally, the thermal stability of various intermediates (viz., CH2CO, Zeo-OCOCH3, and CH3CO + ) were sampled during 50 ps canonical (NVT, i.e., moles, volume, and temperature) AIMD simulations. Notably, the free energy profile in Figure 1c and the critical distance evolutions of direct C-C bond coupling in Figure 2e were reproduced based on our pervious work. [S2] All AIMD simulations were employed by CP2K software, [S3] and the linked PLUMED code [S4] was used to carry out the MTD simulations. The Perdew-Burke-Ernzerhof (PBE) functional [S5] with consideration of Grimme D3 dispersion corrections, [S6] i.e., the PBE-D3 functional, was chosen for the DFT calculations. Limited by large computational cost of triple-zeta basis set in AIMD simulations, the double-zeta (ζ) valence polarized (DZVP) basis set [S7] together with the Goedecker-Teter-Hutter (GTH) pseudopotential [S8] were used for the system. We have tested the convergence of DZVP basis set with SZV, TZVP, and TZVPP basis sets by the intrinsic activation energy (ΔE ‡ , 0 K) of C-C bond coupling of SMS and CO in MOR-12MR, MOR-8MR and SSZ-13 as listed in Table S1. Compared to the qualitative error on ΔE ‡ by single-zeta basis set (SZV), double-zeta plus polarized basis set (DZVP) can quantificationally and quantitatively describe the ΔE ‡ with the error smaller than 11 kJ/mol relative to triple-zeta plus polarized basis set (TZVP) and double polarized basis set (TZVPP). Notably, PBE-D3 functional used in this work, usually underestimates the free energy barriers in zeolite chemistry according previous benchmark studies, [S9, S10] our static DFT calculations based on cluster model also confirmed this point by the free energy barriers of PBE-D3/dgdzvp and B3LYP-D3/dgdzvp in Table S2. However, the static periodic calculations and MTD-AIMD simulations based on PBE-D3/DZVP carried out the free energy barriers highly closed to that by the hybrid functional of HSE06, B3LYP and wB97XD. [S11, S12] Therefore, the PBE-D3/DZVP method used here is reasonable to qualitatively and quantitatively compare the barriers of C-C bond formation between SMS and CO in different zeolites. During the self-consistent field (SCF) procedure, a 360 Ry density CUTOFF criterion with finest grid level was employed along with multi-grids number 4 (NGRID 4 and REL CUTOFF 70). The temperature was controlled by a chain of five Nosé-Hoover thermostats, [S13] and the integration time step was set to 0.5 fs during the AIMD simulations. Figure S1. Periodic models and cluster models of CHA and MOR zeolites used under periodic boundary conditions (PBC) during the AIMD simulations. Cluster models A were used for energy decomposition analysis (EDA) and prediction of orbital energies, while models B were adopted for the intrinsic bond orbital (IBO) analyses with intrinsic reaction coordinate (IRC) pathways.
Metadynamics (MTD) simulation [S14] was employed to explore the CC bond coupling between SMS and carbon monoxide (CO). This is done by first optimizing the geometrical structures and cell parameters of the configuration by periodic DFT followed by a 5 ps NPT molecular dynamics simulation to relax the structure and cell parameters under experimental conditions such as temperature and pressure. Accordingly, MTD simulations were performed at 473 K in the NVT ensemble; an advanced sampling technique employed for enhancing the probability of sampling chemical reactions or rare events provided that a limited number of collective variables (CV) describing the reaction coordinates were well-defined. Nonetheless, MTD simulation is normally biased by regularly spawning Gaussian hills along the chosen CVs, which was defined by coordination numbers (CN): where rij represents the interatomic distance between atoms i and j and r0 is the reference distance.
In this work, parameters d0, n, and m were set to 0, 6, and 12, respectively. Accordingly, CVs may be defined in terms of CNs, viz., CV1 = CN(CSMSOZeo); CV2 = CN(CSMSCCO), which were used to describe the proton transfer process, as shown in Scheme S2. By setting a reference distance of 2.1 Å and 1.5 Å for CN(CSMSOZeo) and CN(CSMSCCO), respectively, and a Gaussian hill width and height of 0.035 and 2 kJ/mol for both CV1 and CV2, and a spawning every 50 fs, MTD simulations were allowed to continue till the height of additional hills no longer affect the resultant free energy profile. Based on the sum of spawned Gaussian hills, the 2D Gibbs free energy profile of the reaction may be reconstructed. Accordingly, a lowest free energy path (LFEP) was derived by means of the MEPSA software. [S15] Scheme S2. Assignments of collective variables (CV1 and CV2) associated with the CC bond coupling between SMS and CO during MTD simulations.
Umbrella sampling method, as the other free energy sampling techiques, was employed to explore the mobility of ketene and acylium ion in MOR. To describe the diffusion path of ketene and acylium ion between 8MR and 12MR channels, we define the distance (d) between two center of mass (COM) as the collective variable (CV), one COM is the ketene or acylium, the other is the 8MR window of side pocket oppsite to BAS, eight windows are respectively the independent simulation that represents a point along the CV from 7 to 0. The additional potential introduced by a simple harmonic term: Where x is the real-time value of CV, x0 is the value of specificed value of CV, k is set to 30 kJ/mol for ketene from 12MR to 8MR channel. All biased umbrella sampling simulations were consisted 10 4 steps at NVT ensemble, and the setting to k lead to the reasonable sampling indicated by the inter-overlaping of different windows. The free energy profile was generated by the sampling obtained in each window based on the weighted histogram anaylsis method. [S16-S17]

Cluster Model Calculations
The intrinsic bond orbitals (IBO) are a non-empirical form of localized molecular orbitals that furnish quantitative interpretation of bonding by assigning electrons in a doubly occupied IBO to the individual atoms. [S18] Together with the intrinsic reaction coordinate (IRC) pathways, IBO evolutions are capable of providing detailed information of orbitals and electrons during bond formation/rupture to ascertain the reaction mechanism. The IBO were obtained with the IboView program [S19] using the B3LYP-D3 and PBE-D3 methods with the 6-31G(d, p) basis set, whilst the IRC pathways were calculated by using the Gaussian 09 program. [S20] To circumvent the computational cost, relevant IBO and IRC calculations were conducted based on simplified cluster model B (cf. Figure S1).
Energy decomposition analysis (EDA) was used to better understand the confinement effect of intermediates in different channels based on static cluster calculations. Based on cluster models clipped from the optimized periodic structures (i.e., model A; see Figure S1), detailed and precise analyses of host-guest and guest-guest interactions may be accomplished based on the ETS-NOCV EDA approach, [S21] which combined the extended transition state (ETS) method with the natural orbital for chemical valence (NOCV) theory by means of the Amsterdam density functional (ADF) program. [S22, S23] The PBE-D3 functional was chosen for the DFT calculations and EDA to warrant the uniformity between dynamic and static calculations. By invoking the frozen core approximation, [S24] a triple-zeta (TZ) Slater type orbital (STO) basis set containing two polarization functions (namely, TZ2P) was adopted for all elements to describe interactions between intermediates and the zeolite framework. Auxiliary STO functions, centered on all nuclei, were used to fit the electron density and to obtain accurate Coulomb potentials in each SCF cycle. For the ETS-NOCV EDA scheme, the interaction energy (ΔEint) between the fragments may further be divided into four components: Where the four energy terms respectively account for the Pauli repulsive interaction among occupied orbitals on the fragments (ΔEPauli), the classical electrostatic interaction between two fragments (ΔEelec), electron distribution of the constituent molecules (ΔEorb), the dispersion interaction (ΔEdisper) due to the use of dispersion corrected PBE-D3 functional.

Sample preparation
Before adsorption of acetyl chloride, H-MOR (Si/Al=10) was shaped in the form of thin wafer (ca. 8.5 mg/cm 2 ) and place in quartz IR cell connected to a high vacuum line for dehydration. The temperature was gradually increased at a rate of 1K/min and the sample was kept at a final temperature of 673 K at a pressure below 10 -3 Pa overnight. After the samples cooling down to ambient temperature, a known amount of acetyl chloride was introduced into the sample, and then the sample glass tube was sealed.

Synchrotron radiation X-ray diffraction
High-resolution SXRD data were collected at beamline I11, Diamond Light Source, Harwell, UK. The energy of the incident X-ray was set at 15 keV. The wavelength (l=0.826576(10) Å) and the 2θ zero point (ZP=0.000361(2)°) were determined by fitting the diffraction data of high-quality silicon powder (SRM640c). Before data collection, sample was loaded into borosilicate glass capillaries (0.5 mm ID) in a glove box. Glass wool was packed on top of the sample. The SR-XRD data were collected in a Debye-Scherrer geometry using multi-analyzer crystals (MAC) detectors in the 2θ range of 0-150 o with 0.001 o data binning. Each data set was collected for 1 h for good statistics.

FT-IR studies
Before adsorption of acetyl chloride, H-MOR (Si/Al=10) was shaped in the form of thin wafer (6 mg·cm -2 ) and place in quartz FT-IR cell connected to a high vacuum line. The sample was in situ thermally treated at 623K under high vacuum (p = 10 -5 mbar) for 1 h, then cooled down to room temperature (250K). The known doses of acetyl chloride (Sigma Aldrich) was adsorbed on the sample which has been tracking by the recording the spectra. The sorption of acetyl chloride was performed till the total or half consumption of the Si(OH)Al band. The FT-IR experiments were performed in rapid-scan mode (80 kHz) using Vertex 70 spectrometer (Bruker). Each spectrum (5 Figure S2. Proposed evolution of ketene and its role in the formation of surface acetate in SSZ-13. Figure S3. Schematic diagram for: (a) ketene diffusion from 12MR channel to 8MR channel, (b) ketene along with the COM distance between ketene and 8MR window (close to BAS) at 473 K.
To explore the mobility of neutral ketene from 12MR channel to 8MR channel in MOR, the diffusion behaviors of ketene in H-MOR have been carried out by umbrella sampling method during the AIMD simulations as displayed in Figure S5a and S5b, (see Computational Details and sampling region of Figure S4 in Supporting Information). The diffusion process of ketene from 12MR channel to 8MR channel is both thermodynamically and kinetically favorable by overcoming a smaller free energy span of 42.2 kJ/mol with 18.4 kJ/mol decreasing on free energy, because ketene has been protonated to be acylium ion with a small free energy span of 11.7 kJ/mol in the path to gradually approach BAS at d = 1.9 Å. It's confirmable that ketene in MOR-12MR can spontaneously diffuse into the MOR-8MR and protonate to be acylium ion at reaction temperature.   Figure S5. 8MR window size of side pocket in MOR and kinetic diameters of CH3COCl, CH2CO and CH3CO + . The window size obtained from IZA website, and the kinetic diameter of CH3COCl are calculated according the method by Mehio et. al. [S28] but the kinetic diameters of CH2CO and CH3CO + are orginated from these of ethene and methane.   As shown the FT-IR spectra of the CH3COCl adsorbed by H-SSZ-13 in Figure S7b and S8a-S8b, the bands in H-SSZ-13 were significantly weaker than these in H-MOR due to constrained diffusion of CH3COCl in the cage-like pore structure of H-SSZ-13. Further, some bands in the C-O stretching vibrations region (2000~2400 cm -1 ) have marginal intensity, if any. The band at 1645~1650 cm -1 was assigned to surface acetyl as consistent with that in H-MOR, (Figure S8b). The intensity of the new band at 2050~2070 cm -1 was enhanced with the increasing loading of CH3COCl, which was attributed to the appearance of ketene species in H-SSZ-13 the most probably perturbed by the interaction with unreacted CH3COCl and HCl. The ketene bands of similar frequencies like 2080 and 2052 cm -1 were observed also in H-MOR. When temperature increased to 323 K, the bands at 2080 and 2052 cm -1 were consumed in favor of the bands at 2170 cm -1 . Notably, the signal of acylium ion was not observed in H-SSZ-13 due to the lack of the confinement effect in CHA topology. Figure S9. The reaction process of ketenes dimerization and further dissociation to CO2 and propadiene in gas phase at 473 K by PBE-D3/6-31G(d, p) method. [S29] In gas phase, ketene easily dimerize to diketene by overcoming a free energy barrier of 93 kJ/mol as a thermodynamic favorable process. The subsequent dissociation of diketene to propadiene and CO2 is kinetically more limited than dimerization process since a free energy barrier of 168 kJ/mol has to be overcome. Notably, whole reaction path, i.e. from ketene to propadiene and CO2, is thermodynamically feasible at 473 K, and these two free energy barriers would be further reduced under the confinement effect of zeolites. As the tautomer of propadiene, propyne is thermodynamically more stable species, it is rapidly formed by the isomerization of propadiene with the participation of Bronsted acid sites (BAS).
According to the isomerization of propadiene over a silica catalyst by Parmentier et. al., [S30] the energy barrier of the isomerization from propadiene to propyne is only 54 kJ/mol under the catalytic effect of silanol. This energy barrier can be further decreased when the isomerization is catalyzed by either HCl or BAS of the higher acidic strength when the CH3COCl is interacting with H-MOR. For this reason, the short-living diketene and propadiene intermediates cannot be observed in FT-IR spectra of CH3COCl adsorbed in the zeolites even if the one FT-IR scan a performed within 0.20 sec. The most probably the higher scanning velocities have to be applied to detect these intermedies on the zeolite surface.