Understanding the Nature of Transition States in the Confined

Subscriber access provided by YORK UNIV

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Understanding the Nature of Transition States in the Confined Nanospace of Different Acidic Zeolites on Desulfurization Mechanism of Thiophene Xingxiang Guo, Xinfeng Mao, Chao Deng, Yingxin Sun, and Sheng Han J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b09924 • Publication Date (Web): 19 Dec 2018 Downloaded from http://pubs.acs.org on December 22, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Understanding the Nature of Transition States in the Confined Nanospace of Different Acidic Zeolites on Desulfurization Mechanism of Thiophene

Xingxiang Guo†, Xinfeng Mao†, Chao Deng††, Yingxin Sun*†, Sheng Han*†

†School

of Chemical and Environmental Engineering, Shanghai Institute of Technology, Shanghai 201418, China

††Jiangsu

Key Laboratory of Pesticide Science, College of Sciences, Nanjing Agricultural University, Nanjing 210095, China

*Corresponding authors: Yingxin Sun and Sheng Han Tel: +86-21-6087-7214

E-mail address: [email protected] and [email protected]

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT We have presented a two-layer ONIOM study on the unimolecular and bimolecular hydrodesulfurization (HDS) mechanisms of thiophene catalyzed by H-Beta, H-MCM-22, H-ZSM-5, and H-FAU zeolites. The rate-determining step for HDS is the hydrogenation step. The bimolecular desulfurization mechanism is more favorable to occur due to the lower free energy barriers, which is good accord with experimental observations. Hammond’s postulate was validated by analyzing the thermodynamics and geometrical structures of transition states (TSs). The difference charge density (DCD), reduced density gradient (RDG), localized orbital locator (LOL), and Fukui functions were used to understand the nature of the TSs. For the dimerization step of two thiophene molecules, the DCD analysis indicates a dramatic electron transfer from the thiophene fragment to the hydrothiophene fragment in the TS. High LOL values between two α-C atoms of two thiophene rings indicate a strong interaction and a high degree of covalence between these two C atoms. For the step of H2S formation, a DCD analysis suggests a clear electron transfer from the H2S fragment to the organic carbon chain in the TS; the plots of RDG and LOL indicate a weak interaction between the C atom of organic chain and H2S fragment. The confined nanospace in zeolite has important stabilization effects on TSs and intermediates. The calculated results suggest that almost all the reaction steps in the unimolecular and bimolecular mechanisms on four zeolites are entropy-decreased processes. The stabilization effect of ZSM-5 framework mainly originates from the VDW attractive interactions. The classic electrostatic interactions favor the formation of all the TSs on MCM-22. In the bimolecular mechanism, the behaviors of TSs in Beta are contrary to those in ZSM-5; the electrostatic interactions favor the formation of TSs over FAU but no clear variation trend of VDW interactions is found.


Page 2 of 48

Page 3 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. INTRODUCTION The removal techniques of sulfur in gasoline have attracted increased attention because the sulfur compound can cause SOx emissions and affects the activity of automobile catalytic converters.1,2 With the rapid growth of global energy consumption and increasingly strict environmental legislation, the demand of improving low-quality fuels increases.3 There are many traditional techniques about the removal of sulfur compound, such as hydrotreatment of fluidized catalytic cracking (FCC) feed, hydrodesulfurization (HDS), and adsorption desulfurization. However, these strategies have some disadvantages. Hydrotreatment needs high operating temperature and economic costs. The HDS process requires an H2 atmosphere, often at high pressures, which can lead to a significant loss of octane number. Although adsorption desulfurization avoids the drawbacks above, the low adsorption capacities and selectivity of materials limit its applications.4 Many studies have showed that acidic zeolites alone can catalyze thiophene desulfurization.5-10 Although the traditional HDS catalysts include metal sulfides such as W-Ni or Ni-Mo dispersed over γ-Al or acidic zeolite support,11 zeolites are very efficient catalysts because of high internal surface areas, good thermal stabilities, and high catalytic activities for many important reactions such as the HDS process. For example, Chica et al. studied the adsorption dynamics and stoichiometry of thiophene in different carriers on H-ZSM-5 and H-Y zeolites. They found that H-Y had a great amount of adsorption.12 Welters et al. studied H(x)/NaY-supported metal sulfide catalysts and found that acidic zeolite can promote the HDS conversion and product selectivity; the acid sites can remove thiophene sulfur under reaction conditions.10 Furthermore, Shan et al. reported that the cracking of C–S bond and hydrogen transfer are two essential steps for thiophene and alkylthiophene desulfurization over USY zeolite; low temperatures are favorable for the hydrogen transfer reactions while the cracking step of C–S bond is dominant at high temperatures.13 The similar conclusions are obtained by the following several research groups. Yu et al. studied the adsorption, desorption, and reactions of thiophene on H-ZSM-5 and Co/H-ZSM-5. Their experimental results showed that thiophene interacts with acidic OH groups in H-ZSM-5 through hydrogen bonding; the hydrogen-bonded thiophene undergoes ring-opening or oligomerization on acidic sites in H-ZSM-5.9 Jaimes et al. postulated that thiophene conversion in ZSM-5 occurs according to the selective adsorption with the formation of carbenium ions as transition states followed by (a) thiophene ring-opening and bimolecular reactions or (b) alkylation reactions.14



Page 4 of 48

Although a great deal of experimental work has been done about thiophene desulfurization over the past few decades, a detailed mechanism on the thiophene cracking remains unclear from a theoretical standpoint due to the complexity of desulfurization process. Little theoretical work has been reported. Saintigny et al. used a

3T

zeolitic

cluster

model

(H3SiOHAl(OH)2OSiH3)

to

study

thiophene

desulfurization by a density functional theory (DFT) method; they found that the presence of hydrogen has no effect on activation barriers of thiophene cracking.15 Li et al. reported a mechanism of thiophene cracking on ZSM-5 zeolite. There are two main steps: protonation of thiophene and C–S bond dissociation.16 Rozanska et al. performed a periodic DFT study to investigate the C–S bond cracking step of thiophene and its derivatives over MOR zeolite. They found that zeolite framework could not stabilize the transition states (TSs) and the activation energies of thiophenic derivative cracking are in the range of 270−320 kJ mol-1.17 Dang investigated the loading dependence of the adsorption mechanism of thiophene in FAU zeolite by grand canonical ensemble Monte Carlo (GCMC) simulations. They revealed that in the whole loading range, the thiophene adsorption mechanism could be divided into two parts: ideal adsorption and insertion adsorption.18 In our recently published studies, two-layer ONIOM calculations were performed on C–S bond cracking of thiophenic derivative over H-Beta19 and HDS of thiophene over H-FAU zeolite.20 The calculated results suggested that a bimolecular HDS mechanism is more favorable due to the lower activation barriers. So far, the relationship between the characteristic of HDS and the different topological structures of zeolites still remains unresolved. It is well-known that the bare cluster model can not explain the influence of zeolite framework, and the periodic DFT approach is too expensive to treat large zeolite framework. A recently developed quantum mechanical/molecular mechanical (QM/MM) method as well as the more general ONIOM (Our-own-N-layered Integrated molecular Orbital + molecular Mechanics) scheme have been extensively used to explore the catalytic mechanisms in zeolites.19-22 The scheme combines the advantages of the high accuracy of quantum mechanics calculations and the high efficiency of molecular mechanics force field. We have investigated successfully several reaction mechanisms on desulfurization process of thiophene catalyzed by zeolites using QM/MM method in previous studies as mentioned above.19,20 In this work, the HDS reaction mechanism of thiophene in a series of acidic zeolites like H-Beta, H-FAU, H-MCM-22, H-ZSM-5 was investigated by a two-layer ONIOM study. These four zeolites are typical catalysts and have wide applications in


Page 5 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the petrochemical industry. An important problem is whether the conclusions from the HDS process are the same over the four different zeolites. Therefore, the following two interesting theoretical questions will be thoroughly considered. (1) What is the rate-determining step for the HDS process over the four different zeolites? (2) What is the relationship between the characteristic of thiophene desulfurization and the topological structure of zeolites? The elucidation of the reaction mechanisms provides useful insights into the elemental steps of reactions and can help to optimize the reaction conditions and design new catalysts for industrial production. 2. COMPUTATIONAL MODELING AND METHODS Four nanoclusters, 188T H-Beta,23 204T H-MCM-22,24 208T H-ZSM-5,25 and 156T H-FAU26 were taken from their lattice structures. As shown in Figure 1, the H-Beta model covers the perpendicular cavity between the channels A and B that are represented by a 12-membered ring (12MR). One silicon atom is replaced with an aluminum atom, and a proton is attached to the bridging oxygen atom (O1) bonded directly to the Al atom (Figures S1, Supporting Information).23 For H-MCM-22, the 204T cluster model consists of three intersecting channels A, B, and C; the aluminum atom prefers the T4 site and the charge-balancing proton is attached to the O1 atom (Figures S2).24 The H-ZSM-5 contains the straight channel A and the zigzag channel B; a silicon atom was substituted with an aluminum atom at T12 and a charge-balancing proton favors bonding with O1 site (Figures S3).25 The H-FAU cluster model covers an area of two supercages connected to each other through the 12MR window; the charge-balancing proton is attached the bridging oxygen atom O1 bonded directly to the Al atom (Figures S4).26 The terminal H atoms of the cluster model were used to saturate the dangling bonds of the Si atoms and the Si–H bonds were fixed at 1.47 Å. A two-layer ONIOM approach was employed in this work, the total energy of the whole system is calculated from the following formula:27 EONIOM = EQM(inner) + EMM(real) – EMM(inner)

(1)

Herein, the real region contains an inner region and a large outer region. Both QM and MM calculations need to be performed for the inner region. The real region is calculated only at the MM level. The inner region consists of the active Brønsted acidic site and reactant; the outer region is the extended zeolite framework. In this study, we employed the nonbonded terms of UFF force field28 to treat the van der Waals (VDW) and classic electrostatic interaction energies. The VDW energies are



Page 6 of 48

calculated from a Lennard-Jones 12-6 type expression:

Evdw

6   r 0 12  rij0   ij    ij     2    r     rij  i j  ij   

(2)

where  ij is the well depth and rij0 is the VDW bond length. The  ij and rij0 are calculated from a geometric mean combination rule:

 ij   i   j

(3)

rij0  ri 0  rj0

(4)

The  i and ri 0 values for atom i can be taken from the UFF parameters.28 The MM electrostatic interaction energies are obtained from the Coulomb’s law, i.e., the classic charge-charge interactions:

Eel  

qi q j

i j

rij

(5)

where the qi is the force field atomic charge of atom i . The DFT method with the M06-2X hybrid meta generalized gradient approximation (meta-GGA) functional29,30 and universal force field (UFF)28 was employed to describe the inner and outer regions, respectively. It is well known that the B3LYP functional undervalues the activation barriers heights of transition states and is not accurate enough to describe medium-range van der Waals interactions.31 Recently, Zhao et al. developed new density functions, called M06-class functions. M06-2X function has excellent performance for accurately predicting main-group thermochemistry, kinetics, noncovalent interactions, electronic excitation energies, and aromatic-aromatic stacking interactions.31 In order to save computational time, we employed the zeolite model called ONIOM(M06-2X(qm1(6-31G(d,p)):qm2(3-21G)):UFF) to divide the inner QM region into two parts. The qm1 subregion includes the active center of zeolite and reactant and is treated with the 6-31G (d,p) basis set; the qm2 subregion is treated with the 3-21G basis set. For H-Beta, H-MCM-22, H-ZSM-5, and H-FAU, the qm1 subregion contains 8T, 9T, 8T, and 12T atoms; the qm2 subregion contains 8T, 23T, 26T, and 6T atoms, respectively. During the geometry optimizations, only the 5T active region [(≡SiO)3Al(OH)Si≡] and the reacting molecule were allowed to relax, whereas the rest of the structure was fixed at the crystallographic coordinates. To


Page 7 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


verify that each transition-state structure has only one imaginary frequency and none for the intermediate, the frequency calculations were carried out at the same level of theory. We performed the intrinsic reaction coordinate (IRC) calculations for the TSs in both the forward and reverse directions to determine two relevant minima.32 In order to gain more accurate interaction energies, single-point energy calculations were further carried out at the M06-2X(qm1(6-311+G(2df,2p)):qm2(6-31G(d,p))) level, the model is called ONIOM(M06-2X//M06-2X:UFF). Since the frequency associated with the soft mode contributes the most to the vibrational entropy, it is also most affected by the numerical error. The harmonic approximation inherently breaks down for large-amplitude nuclear motions and numerical noise becomes significant, as expressed by the presence of spurious low frequencies below 100i cm-1 in magnitude. After carefully checking the related nuclear motions to verify that they corresponded largely to translations or rotations of the reactant molecules, these frequencies were replaced by 12 cm-1 because the entropy contribution from a 12 cm-1 harmonic oscillator approximates that from a pseudotranslational/rotational degree of freedom.33,34 It is noteworthy that if the adsorption energy is calculated as the difference between the total energy of the complex (AB) and the sum of the energies of the separated fragments (A and B), the adsorption energy can always be overestimated due to the basis set superposition error (BSSE). In this work, the BSSE was used to correct the interaction energy using the counterpoise (CP) scheme of Boys and Bernardi.35 The CP method allows an estimation of the BSSE defined as the following equation: BSSE = E(A)A – E(A)AB + E(B)B – E(B)AB

(6)

where E(A)AB is the total energy of fragment A at its equilibrium geometry in the complex AB and obtained at a basis set including all basis functions of fragments A and B without electrons and nuclei of fragment B. E(A)A and E(B)B are the total energies of separated A and B molecules, respectively. In the ONION scheme, the BSSE should be calculated at the QM method only for the inner region. The interaction energy between two molecules A and B can be obtained using the following equation: E(interaction) = EONIOM(AB) – EONIOM(A) – EONIOM(B) + BSSE

(7)

where EONIOM(A), EONIOM(B) and EONIOM(AB) are the total ONIOM energies of the A, B, and AB molecules, respectively. We employed the classic transition state theory (TST) to obtain the rate constant



Page 8 of 48

(k) at 673.15 K in the following form:36 k

 G  k BT exp   h  RT

 k BT  S    H    exp exp      h   R   RT 

(8)

where kB is the Boltzmann constant, h is the Planck constant, and ΔG≠, −TΔS≠, and ΔH≠ are the standard molar Gibbs free energy barrier, entropy loss, and enthalpy barrier for the TSs at 673.15 K, respectively. When the energies were discussed, the final ONIOM (M06-2X//M06-2X:UFF) energy values with thermal corrections at 673.15 K are used unless specifically noted. The difference charge densities (DCDs) were used to evaluate the interaction strength between host and guest molecules.37,38 In this study, the DCDs of the transition states are employed to explain the migration of electrons for different elementary steps. In a given transition-state structure, the electrostatic interaction between the zeolite framework and hydrocarbon fragment is defined as the following equation: Δρ = ρzeolite–hydrocarbon – ρzeolite – ρhydrocarbon

(9)

where ρzeolite–hydrocarbon, ρzeolite, and ρhydrocarbon are the electron densities of the TS, isolated zeolite, and hydrocarbon fragment in the TS, respectively. For distinguishing and visualizing different types of interactions, the isosurface plots of the reduced density gradient (RDG) for important intermediates and transition states in zeolite systems were obtained by calculating the RDG and Ω functions.39 These two functions can be expressed as RDG  r  

1

2  3

  r 



2 1/3

 r 

4/3

  r   Sign  2  r     r 

(10) (11)

where ρ(r) is the total electron density and the λ2 is the second largest eigenvalue of Hessian matrix of electron density. The localized orbital locator (LOL) was proposed to analyze the electronic density. The LOL is defined as the following equations:40

t 

 

 LSDA  exact

t 1  t

0 

(12)

 1


(13)

Page 9 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


 LSDA 

2/3 3 6 2  5/3  5

 exact    i

2

(14) (15)

i

where νσ,  LSDA , and  exact are the LOL function, local spin density approximation, and non-interacting Kohn–Sham kinetic energy density, respectively. The Fukui function and dual descriptor were used to study reactive sites. The Fukui function is defined as    r   f r      N 

(16)

where ρ(r), N, and ν is the electron density, the number of electrons and the external potential, respectively.41 The dual descriptor Δf(r) can be obtained from the following equation:42

f  r    Ns 1  r    Ns 1  r 

(17)

where  Ns 1  r  and  Ns 1  r  are the spin densities for the N+1 and N-1 electron systems, respectively. The condensed dual descriptor ΔfA is often used to reveal the possibility that an atom A could act as a reactive site. The ΔfA can be calculated from the following equations: f A  qNA  qNA1

(nucleophilic attack)

(18)

f A  qNA1  qNA

(electrophilic attack)

(19)

qNA1  qNA1 f  2

(radical attack)

(20)

0 A

where qNA , qNA1 , and qNA1 are the charges of the A atom for the N, N+1, and N-1 electron systems, respectively. The ΔfA is defined as follows. f A  f A  f A

(21)

The transition states in the zeolite channel are cationic in nature for the Brønsted acid catalyzed reactions. The electrostatic interactions between the cation-like transition states and anion-like zeolitic frameworks are crucial to ensure the accurate calculations. However, the charge parameters in the UFF force field were not optimized to represent the electrostatic potentials (ESPs) in zeolites.43 In our work, we employed the M06-2X/6-31G(d,p) level to calculate the ESP energies of four



optimized pure silica zeolite clusters and fit the ESP energies to obtain the MM atomic charge values. The QM atomic charge values were obtained directly from the ESP charges by DFT calculations with the M06-2X functional using the ChelpG scheme.44 The mechanical embedding scheme was used to evaluate the electrostatic interactions between the QM and MM regions. All the calculations were performed using the Gaussian 09 software package.45 The values and plots of DCDs, RDGs, LOLs, and Fukui functions are obtained via Multiwfn software.46 3. RESULTS AND DISCUSSION In our previous study,19,20 we have proposed two reaction mechanisms of HDS process of thiophene molecules over acidic zeolites: unimolecular and bimolecular cracking mechanisms. The corresponding pathways are given in Schemes 1 and 2. All intermediates and TSs are identified in the elementary steps. For each intermediate or TS, energy value is obtained relative to the total energy of the reactant (thiophene) and the zeolite at infinitive separation. The corresponding activation barrier of each TS is obtained as the energy difference between the TS and its previous intermediate. The geometrical parameters of TSs are shown in Figures 2, 3, and S7-S17. The energy profiles are given in Figures S5 and S6 and the stacked column charts of energy barriers are summarized in Figure 4. 3.1. Unimolecular Thiophene Cracking Mechanism. As described in Scheme 1, three main steps are involved in the unimolecular mechanism: (1) C−S bond cracking and ring-opening of the thiophene, (2) hydrogenation of intermediate, and (3) removal of sulfur content. The unimolecular mechanism in the acidic zeolites has been extensively studied experimentally and theoretically by several famous research groups such as Saintigny,15 Rozanska,17 Li,16 Shan,13 and Chica.12 A more detailed discussion will be given below. The first step to the HDS process is the capture of one thiophene molecule by the acid site of zeolite. Two adsorption modes, 1(S) and 2(CC), are obtained; the BSSE-corrected adsorption enthalpies of four zeolites at 673.15 K are given in Table 1. For the Beta zeolite, the 1(S) adsorption mode is slightly more stable than the

2(CC) mode, while the 2(CC) mode is more stable for the other zeolites. The calculated results are the same as the reported value over H-Beta19 and different from that over H-FAU zeolite20 by Sun et al. Although all the adsorption enthalpies are negative, all the BSSE-corrected adsorption Gibbs free energies are positive, suggesting that the adsorption of thiophene is entropy-decreased and unfavorable at


Page 10 of 48

Page 11 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the experimental conditions. However, the thiophene can migrate in the zeolitic pores with high kinetic energies at high experimental temperatures and be frequently close to the acidic sites of zeolites, leading to the subsequent C–S bond cracking and HDS process. It should be noted that in the current study, the BSSE-uncorrected free energies are employed due to the very long computational time for the BSSE corrections of all molecules, unless specifically noted. Figure S7 displays the optimized structures of adsorbed thiophene, complex-TH-1. For the 1(S) adsorption mode, the thiophene interacts with the solid-acid proton from zeolite through an S…H bond. In the 2(CC) adsorption mode, the acidic proton interacts with the Cα = Cβ bond of thiophene. The S…H, H…Cα and H…Cβ distances are small different over four zeolites and in the range of 2.16−2.47, 2.10−2.23, and 2.24−2.60 Å, respectively. The O1–H bond distances are slightly lengthened by less than 0.02 Å. The optimized structures are consistent with the previous experimental observations in the acidic H-ZSM-5 and H-FAU zeolites.6,9 After the adsorption of thiophene on the acidic site of zeolite, the acidic H atom protonates the thiophenic sulfur atom which leads to thiophenic C–S bond cracking and ring-opening. The activation free energy barriers for this step are in the range of 290−450 kJ mol-1 (Beta, 442.19 kJ mol-1; MCM-22, 407.18 kJ mol-1; ZSM-5, 292.96 kJ mol-1; FAU, 354.07 kJ mol-1), with rate constants in the range of 6.8×10-22−2.6×10-10 s-1. Rozanska et al. reported that the activation energy barrier is 318 kJ mol-1 for C–S bond cracking of thiophene catalyzed by H-MOR zeolite.17 The enthalpy barrier and entropy loss are two key factors to determine the activation free energy barrier. Two remarks can be obtained from Table 2. First, all these four free energy barriers are controlled by the enthalpy barriers; the positive enthalpy values indicate that the C–S bond cracking is endothermic. Second, the values of entropy losses are positive, suggesting that the ring-opening process of thiophene is entropy-decreased. Both the enthalpy barrier and entropy loss are positive over four zeolites, which thus leads to the high free energy barriers. Clearly, the C–S bond cracking step is experimentally difficult to occur due to the high activation energies. It is worth noting that the ZSM-5 zeolite with medium-sized pores has the lowest activation energy. In current study, the pore diameters of zeolites are in the following order: ZSM-5 (straight channel, 9.16 Å; zigzag channel, 8.85 Å) < Beta (10.67 Å) < MCM-22 (11.40 Å) ≈ FAU (11.41 Å). Therefore, the moderate pore size is favorable for the formation of transition state structures. The optimized structures of transition state TS-TH-1 are shown in Figure S8. In TS-TH-1, the acidic proton H is bonded to



the sulfur atom; the C–S bond of thiophene is breaking and the new C–O bond is forming. The breaking C–S and forming C–O bond distances are in the range of 2.27−2.97 and 1.63−2.15 Å, respectively. The C–S bond cracking leads to the formation of an alkoxy species, INT-TH-1 (Figure S8). After the C–S bond cracking step, one H2 molecule is co-adsorbed with INT-TH-1 to form complex-TH-2 and then reacts with INT-TH-1 through a transition state TS-TH-2 to produce complex-TH-3. The activation free energy barriers for the hydrogenation step are 244.14 (Beta), 271.53 (MCM-22), 214.95 (ZSM-5), and 264.36 kJ mol-1 (FAU), with rate constants of 1.59×10-6, 1.19×10-8, 2.93×10-4, and 4.29×10-8 s-1, respectively. Again, the ZSM-5 zeolite has the lowest activation energy barrier, which can be explained from the calculated data in Table 2. The entropy loss for ZSM-5 is negative (-6.49 kJ mol-1) but positive for other zeolites, leading to the lower energy barrier for ZSM-5. Our calculated energy barrier is consistent with previously theoretical values (Saintigny et al., about 215 kJ mol-1 on a 3T cluster;15 Sun et al., 251.58 kJ mol-1 at 0 K over H-FAU zeolite20). The free energy barrier is related to the structure of transition state. At the TS-TH-2 shown in Figure S9, both the new C–H and H–Ozeo bonds are forming; the H–H bond is being elongated. The forming C–H, H–Ozeo, and breaking H–H bond distances over Beta, MCM-22, and FAU zeolites are small different and in the range of 1.41−1.46, 1.40−1.42, 0.85−0.88 Å, respectively. The H–H bond distance is lengthened by more than 0.11 Å relative to the isolated H2 molecule (0.74 Å). However, the H–H bond distance is 0.75 Å in the TS-TH-2 over ZSM-5 and is slightly lengthened relative to the isolated H2 molecule, which is responsible for the low activation energy barrier over ZSM-5. A hydrogen bond is formed between sulfur atom and the acidic proton in complex-TH-3 (Figure S10) over zeolites, which is favorable for the subsequent desulfurization process. After the complex-TH-3, the acidic proton attacks the S atom through the transition state TS-TH-3. The activation free energies for this desulfurization step are 284.68 (Beta), 276.18 (MCM-22), 255.69 (ZSM-5), and 310.36 kJ mol-1 (FAU), with rate constants of 1.14×10-9, 5.19×10-9, 2.02×10-7, and 1.16×10-11 s-1, respectively. As shown in Table 2, the ZSM-5 zeolite has the lowest activation energy barrier. Therefore, the ZSM-5 is the best catalyst in the unimolecular HDS of thiophene. The order of reactivity of three steps is hydrogenation > desulfurization > C–S bond cracking. The breaking C−S bond and the forming C−O bond distances for TS-TH-3 are in the range of 2.30–2.94 and


Page 12 of 48

Page 13 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.86–2.35 Å, respectively (Figure S11). After the TS-TH-3, the H2S molecule is formed in the INT-TH-2, as shown in Figure S11. Finally, the H2S molecules diffuse through the zeolite channels and sulfur atoms in thiophene molecules are removed. The rate-determining step for the unimolecular HDS pathway can be obtained by dividing the whole reaction sequence into several different sections. Murdoch reported that each section starts with an intermediate and terminates with another intermediate that is more stable in energy than the previous intermediate.47 The energy difference between the transition state with the highest energy and the initial intermediate in each section should be computed; the largest energy difference will contain the rate-determining step. In our study, only one section, complex-TH-1 to INT-TH-2, is involved. It can be seen from Figure S5, the largest energy difference can be obtained from the transition state TS-TH-2 with the highest energy and the corresponding free energy barriers are 509.61 (Beta), 477.83 (MCM-22), 442.43 (ZSM-5), and 525.14 kJ mol-1 (FAU) relative to the complex-TH-1. Clearly, the rate-determining step for the unimolecular cracking mechanism over four zeolites is the hydrogenation step of INT-TH-1 through TS-TH-2. However, as mentioned above, in the individual steps, the order of reactivity is hydrogenation > desulfurization > C–S bond cracking, which indicates that the hydrogenation step is most favorable. Therefore, the rate-determining step should be obtained from the whole reaction sequence, not individual steps. 3.2. Bimolecular Thiophene Cracking Mechanism. The bimolecular cracking mechanism has been proposed by Shan,13 Aksenov,48 and Chica et al.12,49 As described in Scheme 2, four main steps are involved in the bimolecular mechanism: (1) protonation of thiophene on the acidic sites, (2) dimerization of two thiophene molecules, (3) C−S bond cracking and ring-opening of the dimerized thiophene, and (4) removal of sulfur content. As mentioned above, two adsorption modes, 1(S) and 2(CC), are obtained for the adsorbed thiophene, complex-TH-1. The 1(S) mode is more stable over Beta and FAU and 2(CC) mode more stable over MCM-22 and ZSM-5 when the BSSE-uncorrected free energies are used for saving the computational time as discussed in the previous section 3.1. The adsorbed thiophene molecule (complex-TH-1) is firstly protonated to form the INT-2TH-1 through a transition state TS-2TH-1 (Figure S12); the activation free energy barriers for this step are small (Beta, 87.80 kJ mol-1; MCM-22, 66.98 kJ mol-1; ZSM-5, 67.83 kJ mol-1; FAU, 80.63 kJ mol-1). The second thiophene is then co-adsorbed with the first hydrothiophene (complex-2TH-2) and subsequently



undergoes an electrophilic aromatic substitution step to produce the intermediate INT-2TH-2. The corresponding transition state is TS-2TH-2 (Figure 2) and the activation free energy barriers for this step are 77.90 (Beta), 103.13 (MCM-22), 77.19 (ZSM-5), and 95.24 kJ mol-1 (FAU). The entropy losses for the dimerization step are high and account for the range of 15−46% of the total free energies (Table 3), which indicates that molecular dimerization or decomposition will lead to the large entropic fluctuation. Geobaldo et al. experimentally confirmed the formation of oligomeric thiophene species over H-FAU by IR and UV-Vis spectra.50 Our calculated energy barriers for the dimerization step of two thiophene molecules are low and thus consistent with experimental observations. The resulting INT-2TH-2 can be deprotonated through TS-2TH-3 (Figure S14) and complex-2TH-3 is produced. The activation free energies are in the range of 40−65 kJ mol-1 over four zeolites, suggesting that the deprotonation step will occur easily. The entropy losses of this step over four zeolites are not negligible and account for the range of 11−35% of the Gibbs free energy barriers (Table 3). Interestingly, Li et al.16 reported that the complex-2TH-3 could be formed through a concerted dimerization step of two thiophene molecules with an energy barrier of 101 kJ mol-1 in H-ZSM-5 zeolite. Our IRC calculations suggested that the formation of complex-2TH-3 should occur through a stepwise mechanism (two TSs involved, TS-2TH-2 and TS-2TH-3). The C–S bond cracking of complex-2TH-3 occurs through a transition state TS-2TH-4 and the INT-2TH-3 with a −SH group is produced. De Angelis and Appierto experimentally observed the C−S bond cracking of thiophene and the formation of −SH groups over H-FAU zeolite by an infrared approach.51 Our calculated activation free energies are 181.58 (Beta), 153.97 (MCM-22), 133.92 (ZSM-5), and 125.78 kJ mol-1 (FAU), which are obviously lower than those for the C–S bond cracking step in the unimolecular mechanism mentioned above (more than 290 kJ mol-1 over four zeolites). This comparison indicates that thiophene cracking process should occur preferably via bimolecular mechanism. Rozanska et al.17 reported an activation energy barrier of 318 kJ mol-1 for C–S bond cracking of thiophene catalyzed by H-MOR on the basis of unimolecular mechanism; a lower energy barrier could be expected if the bimolecular mechanism is also investigated. It should be noticed that Li et al.16 reported a lower energy barrier of 66 kJ mol-1 over H-ZSM-5 for the bimolecular C–S bond cracking at the moderate B3LYP/6-31G(d) level; the theoretical level might be responsible for the lower activation energy


Page 14 of 48

Page 15 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


barrier. The structure of TS-2TH-4 is displayed in Figure S15; the breaking C–S bond is elongated (Beta, 3.48 Å; MCM-22, 3.64 Å; ZSM-5, 3.09 Å; FAU, 3.41 Å). The intermediate INT-2TH-3 is a mercaptan species; it contains a −SH group and a C–Ozeo bond close to a C=C double bond (Scheme 2). The INT-2TH-3 can easily be converted into different kinds of mercaptan species through the isomerization of C=C double-bond and isomerization steps of carbonium ions, in accord with the experimental and theoretical results.13,16 The reaction will proceed through a hydrogenation step of INT-2TH-3 (TS-2TH-5 in Figure S16). The activation free energies for this step are in the range of 60−110 kJ mol-1 over four zeolites, with rate constants in the range of 1.3×105−8.2×107 s-1, indicating that the hydrogenation process is easy to occur. The calculated energy barriers are consistent with the previously reported values at 0 K over H-FAU zeolite by Sun et al. (101−122 kJ mol-1 for several different hydrogenation steps in bimolecular mechanism).20 Why are the free energy barriers in the bimolecular mechanism (60−110 kJ mol-1) higher than those in the unimolecular mechanism (214−272 kJ mol-1)? If the question is combined with the comparison of structures between TS-TH-2 and TS-2TH-5, an interesting conclusion can be obtained: the sp2-hybridized vinyl carbenium ion in the TS-TH-2 (Scheme 1) has a lower stability than the sp3-hybridized allyl carbenium ion in the TS-2TH-5 (Scheme 2). So a higher energy barrier is needed in the unimolecular hydrogenation process. The H2S molecule is finally produced and the sulfur atom is removed through a desulfurization step (TS-2TH-6) of complex-2TH-5. The activation free energy barriers are calculated to be 125.31 (Beta), 108.83 (MCM-22), 167.64 (ZSM-5), and 135.98 kJ mol-1 (FAU), which are obviously lower than those in the unimolecular desulfurization step mentioned above (more than 250 kJ mol-1 over four zeolites). Therefore, the H2S molecule can be formed more easily through bimolecular mechanism. Sun et al. reported an activation energy barrier of 136.78 kJ mol-1 at 0 K over H-FAU zeolite employing a two-layer ONIOM approach.20 Our calculated values are generally lower than the result of Saintigny et al.15 (167 kJ mol-1) using a bare cluster model, highlighting the importance of different topological structures of zeolites. The structure of TS-2TH-6 is given in Figure 3. The breaking C−S and the forming C−O bond distances for TS-2TH-6 over four zeolites are in the range of 2.19−2.40 and 2.45−2.98 Å, respectively. From the Schemes 1 and 2, one can see a sp3-hybridized allyl carbenium ion in the TS-2TH-6 (Scheme 2) and a less stable



sp2-hybridized vinyl carbenium ion in the TS-TH-3 (Scheme 1), which probably explains the above-mentioned higher activation energy barrier for desulfurization step in the unimolecular mechanism. Figure S6 gives the energy profiles of the bimolecular HDS mechanism over four zeolites. The rate-determining step can be obtained by the method proposed by Murdoch.47 One section in the bimolecular pathway is involved: complex-TH-1 to complex-2TH-6. The largest energy difference can be obtained from the TS-2TH-5 with the highest energy, and the corresponding free energy barriers are 325.07 (Beta), 319.70 (MCM-22), 322.06 (ZSM-5), and 346.84 kJ mol-1 (FAU) relative to the complex-TH-1. Clearly, the rate-determining step for the bimolecular cracking mechanism over four zeolites is the hydrogenation step of complex-2TH-4 through the TS-2TH-5. It is noteworthy that the activation free energy barriers for the hydrogenation step in the bimolecular HDS mechanism (TS-2TH-5, 319−347 kJ mol-1) are still lower than those in unimolecular mechanism (TS-TH-2, 442−526 kJ mol-1) when the global free energy barriers from the whole reaction sequence are used. Figure 4 displays the stacked column charts of the calculated Gibbs activation free energy barriers. One can see that the sum of energy barriers of three transition states, TS-TH-1, TS-TH-2, and TS-TH-3 in the unimolecular mechanism over ZSM-5 is lower than those over other zeolites (Figure 4a). Similarly, the sums of energy barriers of all six transition states (from TS-2TH-1 to TS-2TH-6) in the bimolecular mechanism over ZSM-5 and FAU are lower than those over Beta and MCM-22 zeolites (Figure 4b). Therefore, the ZSM-5 and FAU would be more promising catalysts in the HDS of thiophene if acidic zeolites are used. Now we make a close analysis of energy profiles as shown in Figures S5 and S6. From the point of view of thermodynamics, the whole HDS process is endothermic for both the unimolecular and bimolecular mechanisms. In our previous investigation on the thiophene cracking catalyzed by H-FAU zeolite,20 we reported that from the energy profiles, the whole HDS process is exothermic, which is different from the current conclusion. Actually, the temperature is a key factor for the energy profiles; we selected the temperature of 673.15 K in current study but 0 K in the previous work. The thermal corrections at 673.15 K to the electronic energies at 0 K should be responsible for the difference of energy profiles mentioned above. On the other hand, although the whole HDS process is endothermic, the high experimental temperatures can provide enough energy for the occurrence of reactions. Hammond’s postulate is an important hypothesis in physical organic chemistry


Page 16 of 48

Page 17 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and useful for understanding the relationship between the geometric structures of transition states and thermodynamics of reactions steps.52 According to the postulate, in an exothermic reaction, the transition state should be closer in energy to the reactants than to the products. Therefore, the transition states will be more geometrically similar to the reactants. In contrast, in an endothermic reaction the transition state is closer in energy and thus more geometrically similar to the products than to the reactants. From the energy profiles in Figures S5 and S6, one can see that except for the hydrogenation and deprotonation steps, all the other steps are generally endothermic over four zeolites. So the transition states, TS-TH-2, TS-2TH-3, and TS-2TH-5 should be geometrically similar to the complex-TH-2, INT-2TH-2, and complex-2TH-4, respectively (See Figures S9, S14, and S16). All other transition states will be structurally similar to their corresponding products. For example, in the TS-TH-3 and TS-2TH-6 for the desulfurization step, the H2S molecule is formed like the corresponding products, INT-TH-2 and complex-2TH-6 (Figures S11 and S17). 3.3. Understanding of the Nature of Transition State Species. It has been well accepted that in zeolite catalysis reactions, the transition-state structures are generally cation-like fragments, while their corresponding previous intermediates are neutrally charged species. The charges will flow between different atoms and the electron density will change from the intermediate to transition state, which reflects the nature of the reaction steps. In this section, we focus on the understanding of the formation of transition state species by analyzing the values of RDG, LOL, DCD, and Fukui functions. In the following discussion, one typical zeolite is often selected as an example because the plots over four zeolites are similar. 3.3.1. C–S Bond Cracking. The isosurface plots of the RDG for complex-TH-1 and TS-TH-1 over FAU are shown in Figures 5a and 5b. The several RDG plots over other three zeolites are given in Figure S18. Generally speaking, a red region appears in both complex-TH-1 and TS-TH-1, indicating the presence of steric hindrance in the thiophene ring. The blue region between the C and S atoms of TS-TH-1 indicates strong attractive interactions, showing that the C–S bond is breaking or forming. The RDG plots on other zeolites are similar to those on FAU (Figure S18). The color-filled maps of LOL for complex-TH-1 and TS-TH-1 over FAU are illustrated in Figures 5c and 5d. As we can see in Figure 5, covalent-bond regions possess high LOL values; the blue circles represent the regions with low electron density between valence shell and inner shell. A low LOL region between C and S atoms indicates the weak covalent effect between C and S atoms and suggests that the



Page 18 of 48

C–S bond is cracking or forming, which is in accord with the prediction by the RDG plots. The crescent-like area around the S atom represents the lone pair electrons. Dual descriptor Fukui function is a useful tool to reveal reactive sites.42 The plot for TS-TH-1 over FAU is given in Figure 6 and the corresponding condensed dual descriptor data (f) are listed in Table 4. One can see that the most positive and negative f values are 0.0745 for C atom and -0.1893 for S atom, indicating that the C atom is nucleophilic site and S atom is electrophilic site. Indeed, the plot of dual descriptor Fukui function of the TS-TH-1 exhibit a large light blue region on C atom and a large light green region on S atom, confirming clearly that these two atoms are strong reactive sites and thus there is a strong interaction between them. The similar conclusion can be obtained for other three zeolites (See Figure S20). This result is also supported by the above-mentioned RDG data of TS-TH-1, i.e., a blue RDG region between the C and S in TS-TH-1 indicates strong attractive interactions. 3.3.2. H2S Formation. We calculated the DCDs of the transition states to elucidate the electrons migration for the formation of H2S (TS-TH-3) in the unimolecular mechanisms. The DCD for the formation of H2S is obtained from the following equation: TS    TS   zeolite   HTS2 S  gas   hydrocarbon  no H 2 S 

(22)

TS where  TS is the electron density of the TS.  zeolite and  HTS2 S  gas  are  hydrocarbon  no H 2 S 

the electron densities of the TS without the H2S fragment and isolated gas H2S molecule at their equilibrium structures of the TS, respectively. The behaviors of DCDs for TS-TH-3 in the four zeolites are the same. After TS-TH-3, the INT-TH-2 and H2S are formed. From the INT-TH-2 with H2S to TS-TH-3, a DCD analysis suggests that for TS-TH-3 over H-Beta zeolite (Figures 7a and 7b), the electron density decreases in the H2S fragment but increases in the organic carbon chain, which reflects a clear electron transfer from the H2S fragment to the INT-TH-2. A similar conclusion can be obtained for the H2S formation in TS-2TH-6 (e.g. Figure S30 for H-MCM-22 zeolite). Figure 8 gives the color-filled maps of LOL for TS-TH-3 over FAU zeolite. The plots of organic chain and H2S fragments are shown in Figures 8b and 8c. From the Figure 8d, it is clear that there is a moderate LOL value between the C atom of organic chain and H2S fragment, indicating a weak interaction between them. The crescent-like lone pair electrons of the S atom is migrating towards the positively charged C atom of organic chain. The similar conclusion can be obtained for other


Page 19 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


three zeolites (See Figure S21-S23). The desulfurization step has been further studied by visualizing the isosurfaces of the RDG. As shown in Figure 9, the RDG isosurface of TS-TH-3 exhibits a light color region between C atom of organic chain and S atom of H2S fragment, confirming the weak interaction between them mentioned above from the color-filled maps of LOL for TS-TH-3. 3.3.3. Dimerization of Two Thiophene Molecules. The DCD for the formation of TS-2TH-2 is obtained from the following equation: TS TS    TS   zeolite  thiophene  hydrothiophene no thiophene   gas 

(23)

TS TS where  TS is the electron density of TS.  zeolite and thiophene  hydrocarbon  no thiophene   gas 

are the electron densities of the TS without the thiophene fragment and isolated gas thiophene molecule, respectively. The DCD plots of TS-2TH-2 over H-ZSM-5 zeolite are shown in Figures 10a and 10b; the DCD analysis indicates a dramatic electron transfer from the thiophene fragment to the hydrothiophene fragment. The electron density decreases in the thiophene but increases in the hydrothiophene. A similar conclusion can be obtained over other three zeolites (Figure S24). Figure 10c gives the isosurface plot of the RDG of TS-2TH-2 over H-ZSM-5 zeolite. There is a red circle between the thiophene ring and hydrothiophene ring, indicating a steric hindrance between these two rings. It is interesting that the red circle is small or does not appear in other zeolites (Figure S25), which may be related to the larger pore sizes of other zeolites. The color-filled maps of LOL for TS-2TH-2 over FAU are shown in Figure 11. High LOL values between two α-C atoms of two thiophene rings indicate a strong interaction and a high degree of covalence between the two C atoms in the dimerization step. A similar conclusion can be obtained over other three zeolites (Figure S26-S28). The plots of dual descriptor Fukui function and condensed dual descriptor values for TS-2TH-2 are shown in Figure 12 and Table 4. The f values for C1 (α-C atom of hydrothiophene fragment) and C2 (α-C atom of thiophene fragment) are 0.0766 and -0.0655, indicating that the C1 and C2 atoms are nucleophilic and electrophilic sites, respectively. Clearly, the result from Fukui data is in good accord with the conclusion from the color-filled maps of LOL. A similar conclusion can be obtained over H-MCM-22 and H-ZSM-5 (Figure S29). 3.4. Confinement Effect of the Different Zeolite Cavities. It has been well known that zeolite cavities have important stabilization effects on transition states and intermediates. Now we discuss the effect of the zeolite



environment on the ONIOM energy barriers by analyzing the QM and MM contributions. The MM contributions can be further decomposed in terms of van der Waals and classic electrostatic interaction energies, as listed in Tables 2 and 3. First, all the free energy barriers in the unimolecular and bimolecular mechanisms are mainly attributed to the QM contributions; all the free energy barriers are controlled by the enthalpy barriers. Second, almost all the entropy losses are positive, showing that these reaction steps on all four zeolites are generally entropy-decreased processes. We calculated the ratio of free volumes to the total volume of four zeolites and found that the ZSM-5 has the smallest percentage of free volumes (ZSM-5, 33.85%; MCM-22, 39.64%; Beta, 44.34%; FAU, 54.76%). Are both the VDW and the electrostatic interactions favorable for the formation of transition states in the ZSM-5 zeolite with medium-sized pores? From Table 3, one can see that on ZSM-5 zeolite, the electrostatic interaction energies are generally positive but the VDW interaction energies are generally negative. The different behaviors suggest that the electrostatic repulsive and VDW attraction forces coexist in the medium-sized ZSM-5 pores and compete with each other. The stabilization effect of ZSM-5 zeolite framework mainly originates from the VDW attractive interactions. However, no quantitative variation trend is found for different transition states on ZSM-5 from Table 2 probably due to the small sizes of molecules in the unimolecular mechanism. The case on MCM-22 is slightly more complicated due to its larger free volume and different from that on ZSM-5. In Tables 2 and 3, all the classic electrostatic interaction energies on MCM-22 zeolite are negative, showing that the classic electrostatic interactions favor the formation of transition states. However, the VDW interactions behave differently in the unimolecular and bimolecular mechanisms. The VDW interaction energies are all negative in the unimolecular mechanism but generally positive in the bimolecular mechanism, indicating that molecules with large sizes will lead to larger spatial hindrance and stronger VDW repulsive forces. A careful comparison of MM energy decomposition from Table 3 suggests that the conclusions in Beta zeolite are contrary to those in ZSM-5 zeolite in the bimolecular mechanism. In Beta zeolite, the VDW interaction energies are generally positive but the electrostatic interaction energies are generally negative. The stabilization effect of Beta zeolite framework mainly originates from the classic electrostatic interactions. For FAU with the largest ratio of free volume, the classic electrostatic interactions still generally favor the formation of transition states in the bimolecular mechanism. However, there is no clear variation trend of VDW


Page 20 of 48

Page 21 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interactions found on FAU zeolite probably due to its larger pore sizes and free volume than other three zeolites. Two remarks can be obtained for the whole MM contributions from Tables 2 and 3. First, it can be seen that except for several deprotonation and hydrogenation steps such as TS-2TH-3 (Beta and MCM-22) and TS-2TH-5 (MCM-22, ZSM-5, and FAU), the MM contributions are small relative to the QM energy values and account for less than 10% of the total ONIOM energies. Second, many negative MM energy values in Tables 2 and 3 suggest that the zeolite environment has a stabilization effect on the corresponding transition state structures. Rozanska et al. carried out a periodic DFT study for the thiophenic derivative cracking over H-MOR zeolite and concluded that the zeolite framework appears not to stabilize the transition state complexes.17 So the different topological structures of zeolites have different stabilization effects on the HDS process. Recently, many works on desulfurization of thiophene over metal zeolites have been reported in the literature.53-62 For example, Liu et al. investigated the mechanism for the HDS of thiophene over zeolite L-supported sulfided Co-Mo catalysts.53 They proposed that the favorable reaction pathway included hydrogenation of thiophene and subsequent C-S bond cracking process. Several other research groups such as Wang,54-56 Reddy,57 Song,58 and Ramirez et al.59 pointed out that MCM-41 supported Co/Ni/Mo systems are effective catalysts for HDS of thiophene and DBT in fuels. In current work, we carried out a preliminary investigation on the HDS of thiophene over Cu and Co substituted Beta zeolite. The Al atom in H-Beta is replaced with a Cu or Co atom. The correct spin multiplicity (SM) of the system was determined by single-point energy calculations. The calculated SM numbers are 3 and 5 for Cu-Beta and Co-Beta zeolites, respectively. Table 5 lists the calculated energy barrier data at the ONIOM(UM06/6-31G(d,p):UFF) level for several important transition states. One can see that all the energy barriers are still high and have no significant reduction relative to those in acidic zeolites, which indicates that the substitution of framework Cu and Co atoms would be inefficient for the HDS of thiophene. The framework metal atoms could not be responsible for the high catalytic activities of metal zeolites. The detailed mechanisms of HDS of thiophene over extra-framework transition metal zeolites have been in progress and will be reported in our future work. 4. CONCLUSIONS This work investigates the HDS mechanisms of thiophene catalyzed by H-Beta,



H-MCM-22, H-ZSM-5, and H-FAU zeolites by a two-layer ONIOM method. The nature of transition states and confinement effect of nanospace in different zeolite channels are revealed in detail. Our major conclusions are outlined as follows. (i) The calculated results suggest that the hydrogenation step is the rate-determining step in both unimolecular and bimolecular mechanisms. A bimolecular HDS mechanism is more favorable than unimolecular mechanism due to the lower activation barriers, which is consistent with the experimental observations. From the point of view of thermodynamics, the whole HDS process is endothermic at 673.15 K. However, the high experimental temperatures can provide enough energy for the occurrence of reactions. Hammond’s postulate is verified by analyzing the thermodynamics and geometrical structures of transition states. (ii) For the step of C–S bond cracking, the plots of RDG and LOL predict a clear interaction between C and S and show that the C–S bond is breaking or forming; the dual descriptor Fukui function suggests that the C atom is nucleophilic site and S atom is electrophilic site. For the formation step of H2S, a DCD analysis suggests a clear electron transfer from the H2S fragment to the organic carbon chain; the plots of RDG and LOL indicate a weak interaction between the C atom of organic chain and H2S fragment. For the dimerization step of two thiophene molecules, the DCD analysis indicates a dramatic electron transfer from the thiophene fragment to the hydrothiophene fragment. High LOL values between two α-C atoms of two thiophene rings indicate a strong interaction and a high degree of covalence between the two C atoms in the dimerization step. The result from Fukui data is in good accord with the conclusion from the color-filled maps of LOL. (iii) The confined zeolite cavities have important stabilization effects on transition states and intermediates. Almost all the reaction steps in the unimolecular and bimolecular mechanisms on four zeolites are entropy-decreased processes. The stabilization effect of ZSM-5 zeolite framework mainly originates from the VDW attractive interactions. The classic electrostatic interactions favor the formation of all the transition states on MCM-22; the VDW interactions are favorable for the unimolecular mechanism but generally unfavorable for the bimolecular mechanism on MCM-22. The conclusions in Beta zeolite are contrary to those in ZSM-5 zeolite in the bimolecular mechanism. The classic electrostatic interactions still generally favor the formation of transition states over FAU in the bimolecular mechanism but no clear variation trend of VDW interactions is found probably due to its larger pore sizes and free volume than other three zeolites.


Page 22 of 48

Page 23 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ASSOCIATED CONTENT Supporting Information The force field charges of zeolitic atoms and BSSE-corrected free energies of all complex structures (Tables S1 and S2), nanoclusters of four zeolites (Figures S1-S4), free energy profiles of unimolecular and bimolecular HDS mechanisms (Figures S5 and S6), optimized structures of transition states and intermediates in unimolecular (Figures S7-S11) and bimolecular (Figures S12-S17) mechanisms, and several typical plots of DCDs, RDGs, LOLs, and Fukui functions (Figures S18-S30). The Supporting Information is available free of charge via the Internet at http://pubs.acs.org. ACKNOWLEDGMENTS This work was funded by Development Foundation for young and middle-aged scientific talents of Shanghai Institute of Technology (Project Number ZQ2018-15), National Science Foundation of China (Project Number 21203118, 21606151, 21504057, and 21707092), IIASA Young Scientists Summer Program (Project Number 21411140044), Shanghai Excellent Technology Leaders Program (Project Number 17XD1424900), Shanghai Leading Talent Program (Project Number 017), Science and Technology Commission of Shanghai Municipality Project (Project Number 18090503800), Shanghai Natural Science Foundation of Shanghai (Project Number 17ZR1441700 and 14ZR1440500), Collaborative Innovation Fund of SIT (Project Number XTCX2015-9), Professor of Special Appointment at Shanghai Institutions of Higher Learning (Eastern Scholar) and Shanghai Association for Science and Technology Achievements Transformation Alliance Program (Project Number LM201680).

REFERENCES (1) Shan, J. H.; Chen, L.; Sun, L. B.; Liu, X. Q. Adsorptive Removal of Thiophene by Cu-Modified Mesoporous Silica MCM-48 Derived from Direct Synthesis. Energ. Fuel. 2011, 25, 3093–3099. (2) Duan, L. H.; Gao, X. H.; Meng, X. H.; Zhang, H. T.; Wang, Q.; Qin, Y. C.; Zhang X. T; Song L. J. Adsorption, Co-adsorption, and Reactions of Sulfur Compounds, Aromatics, Olefins over Ce-Exchanged Y Zeolite. J. Phys. Chem. C 2012, 116, 25748–25756. (3) Moses, P. G.; Hinnemann, B.; Topsøe, H.; Nørskov, J. K. The Hydrogenation and



Direct Desulfurization Reaction Pathway in Thiophene Hydrodesulfurization over MoS2, Catalysts at Realistic Conditions: A Density Functional Study. J. Catal. 2007. 248, 188–203. (4) Hernándezmaldonado, A. J.; Yang, R. T. Desulfurization of Liquid Fuels by Adsorption via π Complexation with Cu(I)-Y and Ag-Y Zeolites. Ind. Eng. Chem. Res. 2003, 42, 123–129. (5) Garcia, C. L.; Lercher, J. A. Adsorption and Surface Reactions of Thiophene on ZSM5 Zeolites. J. Phys. Chem. 1992, 96, 2669–2675. (6) Garcia, C. L.; Lercher, J. A. Hydrogen Bonding of Sulfur Containing Compounds Adsorbed on Zeolite H-ZSM5. J. Mol. Struct. 1993, 293, 235–238. (7) Welters, W. J. J.; Vorbeck, G.; Zandbergen, H. W.; Dehaan, J. W.; Debeer, V. H. J.; Santen, R. A. V. HDS Activity and Characterization of Zeolite-Supported Nickel Sulfide Catalysts. J. Catal. 1994, 150, 155–169. (8) Yu, S. Y.; Li, W.; Iglesia, E. Desulfurization of Thiophene via Hydrogen Transfer from Alkanes on Cation-Modified H-ZSM5. J. Catal. 1999, 187, 257–261. (9) Yu, S. Y.; Garcia-martinez, J.; Li, W.; Meitznerb, G. D.; Iglesia, E. Kinetic, Infrared, and X-ray Absorption studies of Adsorption, Desorption, and Reactions of Thiophene on H-ZSM5 and Co/H-ZSM5. Phys. Chem. Chem. Phys. 2002, 4, 1241–1251. (10) Welters, W. J. J.; Debeer, V. H. J.; Santen, R. A. V. Influence of Zeolite Acidity on Thiophene Hydrodesulfurization Activity. Appl. Catal. A-Gen. 1994, 119, 253–269. (11) Michaud, P.; Lemberton, J. L.; Pérot, G. Hydrodesulfurization of Dibenzothiophene and 4,6-dimethyldibenzothiophene: Effect of an Acid Component on the Activity of a Sulfided NiMo on Alumina Catalyst. Appl. Catal. A-Gen. 1998, 169, 343–353. (12) Chica, A.; Strohmaier, K.; Iglesia, E. Adsorption, Desorption, and Conversion of Thiophene on H-ZSM5. Langmuir 2004, 20, 10982–10991. (13) Shan, H. H.; Li, C. Y.; Yang, C. H.; Zhao, H.; Zhao, B. Y.; Zhang, J. F. Mechanistic Studies on Thiophene Species Cracking over USY Zeolite. Catal. Today 2002, 77, 117–126. (14) Jaimes, L.; De Lasa, H. Catalytic Conversion of Thiophene under Mild Conditions over a ZSM-5 Catalyst. A Kinetic Model. Ind. Eng. Chem. Res. 2009, 48, 7505–7516. (15) Saintigny, X.; Santen, R. A. V.; Clémendot, S.; Hutschka, F. A Theoretical Study


Page 24 of 48

Page 25 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of the Solid Acid Catalyzed Desulfurization of Thiophene. J. Catal. 1999, 183, 107–118. (16) Li, B. R.; Guo, W. P.; Yuan, S. P.; Hu, J.; Wang, J. G.; Jiao, H. J. A Theoretical Investigation into the Thiophene-Cracking Mechanism over Pure Brønsted Acidic Zeolites. J. Catal. 2008, 253, 212–220. (17) Rozanska, X.; Santen, R. A. V.; Hutschka, F.; Hafner, J. A Periodic Density Functional Theory Study of Thiophenic Derivative Cracking Catalyzed by Mordenite. J. Catal. 2003, 215, 20–29. (18) Dang, S. Q.; Zhao, L.; Gao, J. S.; Xu, C. M. Loading Dependence of the Adsorption Mechanism of Thiophene in FAU Zeolite. Ind. Eng. Chem. Res. 2016, 55, 11801–11808. (19) Mao, X. F.; Sun, Y. X.; Pei, S. P. A Theoretical Investigation into Thiophenic Derivative Cracking Mechanism over Acidic and Cation-exchanged Beta Zeolites. Comput. Theor. Chem. 2015, 1074, 112–124. (20) Sun, Y. X.; Mao, X. F.; Pei, S. P. A Two-layer ONIOM Study of Thiophene Cracking Catalyzed by Proton- and Cation-exchanged FAU Zeolite. J. Mol. Model. 2016, 22, 1–16. (21) Boronat, M.; Viruela, P. M.; Corma, A. Reaction Intermediates in Acid Catalysis by Zeolites: Prediction of the Relative Tendency to form Alkoxides or Carbocations as a Function of Hydrocarbon Nature and Active Site Structure. J. Am. Chem. Soc. 2004, 126, 3300–3309. (22) Vreven, T.; Morokuma, K. On the Application of the IMOMO (Integrated Molecular Orbital + Molecular Orbital) Method. J. Comput. Chem. 2015, 21, 1419–1432. (23) Newsam, J. M.; Treacy, M. M. J.; Koetsier, W. T.; de Gruyter, C. B. Structural Characterization of Zeolite Beta. Proc. R. Soc. Lond. A 1988, 420, 375–405. (24) Leonowicz, M. E.; Lawton, J. A.; Lawton, S. L.; Rubin, M. K. MCM-22: A Molecular Sieve with Two Independent Multidimensional Channel Systems. Science 1994, 264, 1910–1913. (25) Van Koningsveld, H.; Vanbekkum, H.; Jansen, J. C. On the Location and Disorder of the Tetrapropylammonium (TPA) Ion in Zeolite ZSM-5 with Improved Framework Accuracy. Acta Crystallogr. B 1987, 43, 127–132. (26) Hriljac, J. A.; Eddy, M. M.; Cheetham, A. K.; Donohue, J. A.; Ray, G. J. Powder Neutron Diffraction and 29Si MAS NMR Studies of Siliceous Zeolite-Y. J. Solid State Chem. 1993, 106, 66–72.



(27) Maseras, F.; Morokuma, K. IMOMM: A New Integrated Ab Initio + Molecular Mechanics Geometry Optimization Scheme of Equilibrium Structures and Transition States. J. Comput. Chem. 1995, 16, 1170–1179. (28) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. UFF, A Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024–10035. (29) Zhao, Y.; Truhlar, D. G. Benchmark Data for Interactions in Zeolite Model Complexes and Their Use for Assessment and Validation of Electronic Structure Methods. J. Phys. Chem. C 2008, 112, 6860–6868. (30) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215–241. (31) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157–167. (32) Gonzalez, C.; Schlegel, H. B. Reaction Path Following in Mass-weighted Internal Coordinates. J. Phys. Chem. 1990, 94, 5523–5527. (33) Brogaard, R. Y.; Wang, C. M.; Studt, F. Methanol–Alkene Reactions in Zeotype Acid Catalysts: Insights from a Descriptor-Based Approach and Microkinetic Modeling. ACS Catal. 2014, 4, 4504–4509. (34) Brogaard, R. Y.; Henry, R.; Schuurman, Y.; Medford, A. J.; Moses, P. G.; Beato, P.; Svelle, S.; Nørskov, J. K.; Olsbye, U. Methanol-to-Hydrocarbons Conversion: the Alkene Mmethylation Pathway. J. Catal. 2014, 314, 159–169. (35) Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553–566. (36) Eyring, H. The Activated Complex in Chemical Reactions. J. Chem. Phys. 1935, 3, 107–115. (37) Li, H.; Ji, Y.; Wang, F.; Li, S. F.; Sun, Q.; Jia, Y. Ab Initio Study of Larger Pbn Clusters Stabilized by Pb7 Units Possessing Significant Covalent Bonding. Phys. Rev. B 2011, 83, No. 075429. (38) Chu, Y. Y.; Han, B.; Zheng, A. M.; Yi, X. F.; Deng, F. Pore Selectivity for Olefin Protonation Reactions Confined inside Mordenite Zeolite: A Theoretical Calculation Study. J. Phys. Chem. C 2013, 117, 2194–2202.


Page 26 of 48

Page 27 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(39) Johnson, E. R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, J.; Cohen, A. J.; Yang, W. T. Revealing Noncovalent Interactions. J. Am. Chem. Soc. 2010, 132, 6498–6506. (40) Becke, A. D. Simulation of Delocalized Exchange by Local Density Functionals. J. Chem. Phys. 2000, 112, 4020–4026. (41) Parr, R. G.; Yang, W. Density Functional Approach to the Frontier-Electron Theory of Chemical Reactivity. J. Am. Chem. Soc. 1984, 106, 4049–4050. (42) Morell, C.; Grand, A.; Toro-Labbé, A. New Dual Descriptor for Chemical Reactivity. J. Phys. Chem. A 2005, 109, 205–212. (43) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. UFF, A Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024–10035. (44) Heinz, H.; Suter, U. W. Atomic Charges for Classical Simulations of Polar Systems. J. Phys. Chem. B 2004, 108, 18341–18352. (45) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A., et al. Gaussian 09, Revision A.1; Gaussian, Inc.: Wallingford, CT, 2009. (46) Lu, T.; Chen, F. W. Multiwfn: A Multifunctional Wavefunction Analyzer. J. Comput. Chem. 2012, 33, 580–592. (47) Murdoch, J. R. What Is the Rate-limiting Step of A Multistep Reaction? J. Chem. Educ. 1981, 58, 32–36. (48) Aksenov, D. G.; Klimov, O. V.; Echevskii, G. V.; Paukshtis, E. A.; Budneva, A. A. Thiophene Conversion in the BIMF Process. React. Kinet. Catal. L. 2004, 83, 187–194. (49) Chica, A.; Strohmaier, K. G.; Iglesia, E. Effects of Zeolite Structure and Aluminum Content on Thiophene Adsorption, Desorption, and Surface Reactions. Appl. Catal. B-Environ. 2005, 60, 223–232. (50) Geobaldo, F.; Palomino, G. T.; Bordiga, S.; Zecchina, A.; Areán, C. O. Spectroscopic Study in the UV-Vis, Near and Mid IR of Cationic Species Formed by Interaction of Thiophene, Dithiophene and Terthiophene with the Zeolite H-Y. Phys. Chem. Chem. Phys. 1999, 1, 561–569. (51) De Angelis, B. A.; Appierto, G. Infrared Spectroscopic Study of Thiophene Adsorbed on Zeolites. J. Colloid Interf. Sci. 1975, 53, 14–19. (52) Hammond, G. S. A Correlation of Reaction Rates. J. Am. Chem. Soc. 1955, 77, 334–338.



(53) Liu, B.; Zhao, Z.; Wang, D.; Liu, J.; Chen, Y.; Li, T.; Duan, A. J.; Jiang, J. Y. A Theoretical Study on the Mechanism for Thiophene Hydrodesulfurization over Zeolite L-supported Sulfided Co, Mo Catalysts: insight into the hydrodesulfurization over zeolite-based catalysts. Comput. Theor. Chem. 2015, 1052, 47–57. (54) Wang, A.; Wang, Y.; Kabe, T.; Chen, Y.; Ishihara, A.; Qian, W. Hydrodesulfurization of Dibenzothiophene over Siliceous MCM-41-supported Catalysts: I. Sulfided Co-Mo Catalysts. J. Catal. 2001, 199, 19–29. (55) Wang, A.; Li, X.; Chen, Y.; Han, D.; Wang, Y.; Hu, Y.; Kabe, T. Deep Hydrodesulfurization over W-based Catalysts Supported by Siliceous MCM-41. Chem. Lett. 2001, 30, 474–475. (56) Wang, A.; Wang, Y.; Kabe, T.; Chen, Y.; Ishihara, A.; Qian, W.; Yao, P. Hydrodesulfurization of Dibenzothiophene over Siliceous MCM-41-supported Catalysts: II. Sulfided Ni-Mo Catalysts. J. Catal. 2002, 210, 319–327. (57) Reddy, K. M.; Wei, B.; Song, C. Mesoporous Molecular Sieve MCM-41 Supported Co-Mo Catalyst for Hydrodesulfurization of Petroleum Resids. Catal. Today 1998, 43, 261–272. (58) Song, C.; Reddy, K. M. Mesoporous Molecular Sieve MCM-41 Supported Co-Mo Catalyst for Hydrodesulfurization of Dibenzothiophene in Distillate Fuels. Appl. Catal. A-Gen. 1999, 176, 1–10. (59) Ramirez, J.; Contreras, R.; Castillo, P.; Klimova, T.; Zárate, R.; Luna, R. Characterization and Catalytic Activity of CoMo HDS Catalysts Supported on Alumina-MCM-41. Appl. Catal. A-Gen. 2000, 197, 69–78. (60) Klimova, T.; Calderón, M.; Ramirez, J. Ni and Mo Interaction with Al-containing MCM-41 Support and its Effect on the Catalytic Behavior in DBT Hydrodesulfurization. Appl. Catal. A-Gen. 2003, 240, 29–40. (61) Cui, J.; Yue, Y. H.; Sun, Y.; Dong, W. Y.; Gao, Z. Characterization and Reactivity of Ni, Mo-supported MCM-41 Catalysts for Hydrodesulfurization. Stud. Surf. Sci. Catal. 1997, 105, 687–694. (62) Sampieri, A.; Blanchard, J.; Breysse, M. ; Brunet, S.; Fajerwerg, K.; Louis, C.; Pérot, G. Hydrodesulfurization of Dibenzothiophene on MoS2/MCM-41 and MoS2/SBA-15 Catalysts Prepared by Thermal Spreading of MoO3. Catal. Today 2005, 107–108, 537–544.


Page 28 of 48

Page 29 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Scheme 1. The unimolecular HDS mechanism of thiophene over four zeolites.



Scheme 2. The bimolecular HDS mechanism of thiophene over four zeolites.


Page 30 of 48

Page 31 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Channel B

Channel B

Channel A Channel C

Channel B

Channel A

(a)

(b)

Channel B

Channel A

Supercage

(c)

(d)

Figure 1. Nanocluster models of four zeolites divided into two areas: the inner region, shown by colored balls, and the outer region by sticks. (a) H-Beta, (b) H-MCM-22, (c) H-ZSM-5, and (d) H-FAU.



C7 C8

C4

C6

S

C5 1.83

Page 32 of 48

C3

C3 C4 1.90

C2

S

C6 C7 C5 C8 1.93 H S

C2

C1

H

C1

O2 O1

O1

Al

Al

(a)

(b)

S

Al

O2

O2 H O1 2.08

C8 C7

C5

C6

2.13 C4 S C3 C1 C2

H C3 C2 S C1 C4

C5

H

1.99

C8 C7 C6

Al

(c)

O1

(d)

Figure 2. Optimized structures of TS-2TH-2 over (a) H-Beta, (b) H-MCM-22, (c) H-ZSM-5, and (d) H-FAU.


Page 33 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


H

C

H

S H 2.45

SH 2.61

2.40

C

O2

O1

Al

O2

(a)

H

2.22

O1

(b)

H S 2.98

C

2.19 O2

O1

H S H 2.60 C H 2.23

Al

O2

Al

(c)

(d)

Figure 3. Optimized structures of TS-2TH-6 over (a) H-Beta, (b) H-MCM-22, (c) H-ZSM-5, and (d) H-FAU.



Free energy barriers (kJ mol-1)

1200

TS-TH-1

1000

TS-TH-2

TS-TH-3

800 600 400 200 0 H-Beta

H-M CM -22

H-ZSM -5

H-FAU

(a)

800 -1

Free energy barriers (kJ mol )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 48

TS-2TH-1 TS-2TH-4

700

TS-2TH-2 TS-2TH-5

TS-2TH-3 TS-2TH-6

600 500 400 300 200 100 0 H-Beta

H-MCM-22

H-ZSM-5

H-FAU

(b) Figure 4. Stacked column charts of the calculated Gibbs activation free energy barriers in (a) unimolecular and (b) bimolecular mechanisms.


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


H-Bonding VDW Repulsive

Page 35 of 48

(a)

(b)

(c)

(d)

Figure 5. Isosurface plots of (a,b) RDGs and color-filled maps of (c,d) LOLs. (a,c) complex-TH-1 and (b,d) TS-TH-1 over H-FAU zeolite. VDW represents the van der Waals interactions and the atoms in the black dotted lines represent the C and S atoms.



(a)

Page 36 of 48

(b)

Figure 6. Plot of dual descriptor Fukui function for TS-TH-1 over H-FAU. (a) positive and (b) negative. The atoms in the black dotted lines represent the C and S atoms.


Page 37 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

Figure 7. Plots of DCDs for TS-TH-3 over H-Beta. (a) positive DCD and (b) negative DCD. The positive DCD represents where the charge density increases and the negative DCD represents where the charge density decreases. The molecule in the black dotted lines represents the H2S.



Page 38 of 48

organic chain

H2S

(a)

(b)

(c)

(d)

Figure 8. Color-filled maps of LOL for TS-TH-3 over H-FAU. (a) TS-TH-3 model, (b) plot of organic carbon chain, (c) plot of H2S fragment, and (d) plot of the whole TS-TH-3 fragment.


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

(d)

Figure 9. RDG plots for TS-TH-3 over four zeolites. (a) H-Beta, (b) H-MCM-22, (c) H-ZSM-5, and (d) H-FAU. The molecule in the black dotted lines represents the H2S.



Page 39 of 48


(a)

(b)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 48

(c) Figure 10. Plots of DCDs and RDG for TS-2TH-2 over H-ZSM-5. (a) positive DCD, (b) negative DCD, and (c) RDG. The molecules in the dotted and solid lines represent the thiophene and hydrothiophene fragments, respectively.


Page 41 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


hydrothiophene thiophene

(a)

(b)

(c)

(d)

Figure 11. Color-filled maps of LOL for TS-2TH-2 over H-FAU. (a) TS-2TH-2 model, (b) plot of thiophene fragment, (c) plot of hydrothiophene, and (d) plot of the whole TS-2TH-2 fragment.



Page 42 of 48

C1

C1

C2

C2

(a)

(b)

C1

C1 C2

C2

(c)

(d)

Figure 12. Plots of dual descriptor Fukui function for TS-2TH-2 over (a,b) H-FAU and (c,d) H-Beta zeolites. (a,c) positive values and (b,d) negative values. The C1 and C2 atoms represent α-C atoms in hydrothiophene and thiophene fragments, respectively.


Page 43 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Calculated BSSE-corrected Adsorption Enthalpies (ΔHads) and Gibbs Free Energies (ΔGads) of Thiophene at the ONIOM(M06-2X//M06-2X:UFF) Level at 673.15 Ka modes

energies

H-Beta

H-MCM-22

H-ZSM-5

H-FAU

1(S)

ΔHads ΔGads ΔHads ΔGads

-67.51 33.58 -64.58 29.14

-45.61 52.33 -59.43 47.81

-50.59 54.83 -57.51 58.85

-46.53 43.06 -62.01 41.42

2(CC) aAll

energy values are reported in kilojoules per mole.



Page 44 of 48

Table 2. Calculated Rate Constants (k), Enthalpy Barriers (ΔH≠), Entropy Losses (−TΔS≠), Free Energy Barriers (ΔG≠) at 673.15 K, QM and MM Contributions to the Free Energy Barriers, and van der Waals (VDW) and Electrostatic (elec) Contributions to the MM Energies for TS-TH-1, TS-TH-2, TS-TH-3 in Unimolecular HDS over Four Zeolitesa zeolites Beta

MCM-2 2 ZSM-5

FAU

aAll

species

k

ΔH≠

−TΔS≠

ΔG≠

QM

MM

VDW

elec

TS-TH-1 TS-TH-2 TS-TH-3 TS-TH-1 TS-TH-2 TS-TH-3 TS-TH-1 TS-TH-2 TS-TH-3 TS-TH-1 TS-TH-2 TS-TH-3

6.81×10-22

414.52 213.17 287.95 400.50 237.01 271.67 283.75 221.44 265.19 328.18 238.23 288.02

27.67 30.97 -3.27 6.68 34.52 4.51 9.21 -6.49 -9.50 25.89 26.13 22.34

442.19 244.14 284.68 407.18 271.53 276.18 292.96 214.95 255.69 354.07 264.36 310.36

426.94 235.98 291.51 419.02 284.91 295.32 298.35 210.17 262.41 351.71 264.01 307.29

15.25 8.16 -6.83 -11.84 -13.38 -19.14 -5.39 4.78 -6.72 2.36 0.35 3.07

2.87 2.68 5.80 -6.02 -11.35 -14.44 0.28 4.09 -6.26 5.31 -0.18 -0.68

12.39 5.48

1.59×10-6 1.14×10-9 3.55×10-19 1.19×10-8 5.19×10-9 2.59×10-10 2.93×10-4 2.02×10-7 4.69×10-15 4.29×10-8 1.16×10-11

-12.63

-5.82 -2.03 -4.70 -5.67 0.69 -0.45 -2.95 0.53 3.75

energy values are reported in kilojoules per mole and rate constants in inverse seconds.


Page 45 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3. Calculated Rate Constants (k), Enthalpy Barriers (ΔH≠), Entropy Losses (−TΔS≠), Free Energy Barriers (ΔG≠) at 673.15 K, QM and MM Contributions to the Free Energy Barriers, and van der Waals (VDW) and Electrostatic (elec) Contributions to the MM Energies for TS-2TH-1, TS-2TH-2, TS-2TH-3, TS-2TH-4, TS-2TH-5, TS-2TH-6 in Bimolecular HDS over Four zeolitesa zeolites

species

k

ΔH≠

−TΔS≠

ΔG≠

QM

MM

VDW

elec

Beta

TS-2TH-1

2.16×106

72.38

15.42

87.80

89.07

-1.27

1.49

-2.76

TS-2TH-2

1.26×107

66.02

11.88

77.90

82.22

-4.32

-6.64

2.32

TS-2TH-3

4.84×109

39.70

4.92

44.62

39.20

5.42

4.43

1.00

TS-2TH-4

1.14×101

170.85

10.73

181.58

183.21

-1.63

4.23

-5.86

TS-2TH-5

1.67×105

109.43

-7.31

102.12

104.05

-1.93

14.38

-12.45

TS-2TH-6

2.65×103

115.36

9.95

125.31

119.45

5.86

12.76

-6.90

TS-2TH-1

8.90×107

54.67

12.31

66.98

68.54

-1.56

14.31

-15.86

TS-2TH-2

1.39×105

72.37

30.76

103.13

110.71

-7.58

-3.57

-4.01

TS-2TH-3

9.13×108

39.19

14.76

53.95

68.57

-14.62

-6.74

-7.88

TS-2TH-4

1.58×101

143.08

10.89

153.97

148.11

5.86

15.23

-9.38

TS-2TH-5

1.32×105

76.70

26.75

103.45

91.33

12.12

19.86

-7.74

TS-2TH-6

5.03×104

112.33

-3.50

108.83

109.37

-0.54

10.11

-10.65

TS-2TH-1

7.64×107

71.28

-3.45

67.83

67.20

0.63

-3.94

4.58

TS-2TH-2

1.44×107

41.44

35.75

77.19

77.83

-0.64

3.23

-3.87

TS-2TH-3

1.60×108

50.43

13.27

63.70

63.59

0.11

-0.72

0.83

TS-2TH-4

5.68×102

121.91

12.01

133.92

133.08

0.84

-0.56

1.40

TS-2TH-5

8.15×107

62.33

5.14

67.47

82.45

-14.98

-17.22

2.23

TS-2TH-6

1.37×100

165.78

1.86

167.64

179.68

-12.04

-16.19

4.15

TS-2TH-1

7.76×106

51.45

29.18

80.63

82.29

-1.66

-0.66

-1.01

TS-2TH-2

5.71×105

68.46

26.78

95.24

96.42

-1.18

-0.67

-0.51

TS-3TH-3

5.82×109

28.53

15.05

43.58

41.18

2.40

-0.46

2.85

TS-2TH-4

2.43×103

129.90

-4.12

125.78

119.06

6.72

12.31

-5.59

TS-2TH-5

1.88×106

84.86

3.70

88.56

102.57

-14.01

-10.49

-3.51

TS-2TH-6

4.11×106

122.83

13.15

135.98

137.91

-1.93

0.28

-2.21

MCM-22

ZSM-5

FAU

aAll

energy values are reported in kilojoules per mole and rate constants in inverse seconds.



Page 46 of 48

Table 4. Condensed Dual Descriptor Values over Four Zeolites zeolites

species

atoms

N

N-1

N+1

f-

f+

f

Beta

TS-TH-1

C S C1 C2 C S C1 C2 C S C1 C2 C S C1 C2

0.0950 0.1146 0.0287 -0.0074 0.2056 -0.0074 0.0388 -0.0225 0.1571 0.1356 0.0555 -0.0284 0.0688 -0.1934 0.0478 -0.0284

-0.0241 -0.0491 -0.0204 -0.0190 -0.0201 -0.0896 -0.0215 -0.0291 -0.0284 0.0120 -0.0394 -0.0296 -0.0661 -0.3374 -0.0291 -0.0309

0.1531 0.2868 0.0334 0.0245 0.2369 0.3303 0.0429 0.0202 0.2302 0.3006 0.0540 0.0553 0.1292 0.1399 0.0480 0.0396

-0.1191 -0.1637 -0.0491 -0.0116 -0.2257 -0.0822 -0.0602 -0.0065 -0.1855 -0.1236 -0.0949 -0.0012 -0.1349 -0.1440 -0.0768 -0.0025

-0.0581 -0.1722 -0.0047 -0.0319 -0.0313 -0.3376 -0.0041 -0.0427 -0.0731 -0.1651 0.0015 -0.0837 -0.0604 -0.3333 -0.0003 -0.0680

0.0610 -0.0085 0.0444 -0.0203 0.1944 -0.2554 0.0561 -0.0362 0.1124 -0.0415 0.0964 -0.0825 0.0745 -0.1893 0.0766 -0.0655

TS-2TH-2 MCM-22

TS-TH-1 TS-2TH-2

ZSM-5

TS-TH-1 TS-2TH-2

FAU

TS-TH-1 TS-2TH-2


Page 47 of 48 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 5. Calculated Free Energy Barriers (ΔG≠) at 673.15 K for TS-2TH-1, TS-2TH-2, TS-2TH-4, TS-2TH-5, and TS-2TH-6 in Bimolecular HDS over H-Beta zeolite at the ONIOM(UM06/6-31G(d,p):UFF) Levela energies

TS-2TH-1

TS-2TH-2

TS-2TH-4

TS-2TH-5

TS-2TH-6

Co

74.31

64.24

142.46

141.36

162.28

Cu

158.61

57.52

178.91

155.65

222.36

aAll

energy values are reported in kilojoules per mole.



TOC Graphic kJ mol-1 700

H-Beta H-ZSM-5

H-MCM-22 H-FAU

600

500

400

300

desulfurization

dimerization

Gibbs free energy barriers

200

100

0

TS-2TH-1 TS-2TH-2 TS-2TH-3 TS-2TH-4 TS-2TH-5 TS-2TH-6

We have presented a two-layer ONIOM study on the unimolecular and bimolecular hydrodesulfurization mechanisms of thiophene catalyzed by H-Beta, H-MCM-22, H-ZSM-5, and H-FAU zeolites. The bimolecular desulfurization mechanism is more favorable to occur due to the lower free energy barriers. It includes six reaction steps, protonation of one thiophene, dimerization of two thiophene molecules, deprotonation of dimerized thiophene, C-S bond cracking, hydrogenation of intermediate, and desulfurization. The Gibbs activation free energy barriers over four zeolites are given in the TOC graphic above. The localized orbital locator (LOL) plots of transition states for the dimerization and desulfurization steps on H-FAU as a typical example are also given in the TOC graphic.


Page 48 of 48

Understanding the Nature of Transition States in the Confined

Recommend Documents