A theoretical investigation on the hydrodesulphurisation mechanism of hydrogenated thiophene over Cu–Mo-modified FAU zeolite

ABSTRACT We have theoretically investigated the hydrodesulphurisation (HDS) mechanism of hydrogenated thiophene over Cu–Mo-modified FAU zeolite using a two-layer ONIOM (our Own N-layered Integrated molecular Orbital and molecular Mechanics) study. The thiophene is hydrogenated to 2,3-dihydrothiophene (2,3-DHT), 2,5-dihydrothiophene (2,5-DHT) and tetrahydrothiophene (THT) due to moderate free energy barriers. Hydrogenolysis desulphurisation (HYD) and concerted direct desulphurisation (DDS) are discussed. Ring-opening, hydrogen transfer and C−S bond cleavage steps of thiophene derivatives are involved in the HYD process. The rate-determining steps are the hydrogen transfer step for 2,5-DHT and C−S bond cracking step for 2,3-DHT and THT. The concerted DDS pathway is probably more favourable than the HYD pathway in the desulphurisation of 2,5-DHT. The difference charge density (DCD) analysis reveals that for the ring-opening process, the electrons are migrated from the organic chain to the Cu–Mo catalytic centre. The reduced density gradient (RDG) plots indicate that both steric hindrance and a weak van der Waals (VDW) interaction exist between organic fragment and catalytic centre for all transition states (TSs). The localised orbital locator (LOL) maps suggest that there are strong covalent interactions and weak VDW interactions between the atoms in the forming chemical bonds.


Introduction
Organic sulphides in fossil fuels are mainly thiophene, benzothiophene (BT), dibenzothiophene (DBT) and 4,6-dimethyldibenzothiophene (4,6-DMDBT), which are harmful to the environment [1]. The combustion of these sulphides promotes the generation of SO x gas and acid rain. Some organic sulphides lead to the deactivation of catalysts and corrode refining equipment and internal combustion engine compressors [2]. Therefore, the control of sulphur content in fuel becomes stricter. Some desulphurisation technologies such as hydrodesulphurisation (HDS), oxidative desulphurisation (ODS) and adsorption desulphurisation have been used [2][3][4]. As the demand for ultra-low sulphur fuels increases, the HDS techniques require larger reactors and higher operating pressures. The ODS easily oxidises the sulphide in the fuel into advanced sulphur composition, but the cost is too high and the technology still needs to be improved. Although adsorption desulphurisation avoids the above-mentioned drawbacks, the low adsorption capacity and selectivity of the materials limit its application [5].
In recent years, transition metal-modified zeolite catalysts are widely used in the industrial desulphurisation of thiophene. The metal-modified zeolites can be employed in both adsorptive desulphurisation and HDS of thiophene. For the adsorptive desulphurisation, the metal zeolites can enhance the adsorption ability for the thiophene through the π-complexation between metal atoms and thiophene molecules. The acidic zeolites have low adsorption abilities because of the weak interaction between acidic protons and thiophene. For the HDS process, metal zeolites are more suitable because the acidic zeolites often lead to various by-products through hydrocarbon conversion reactions, and thus decrease the octane number of gasoline. Many transition metal-modified zeolites such as Cu-Na-Y [6], Cu(I)-Y [7], AgY [8], nickel (II)-Y [9], Ni-Mo-MCM-41 [10], CuCeY [11] are widely used for the adsorption desulphurisation and HDS process. The metal atoms such as Cu, Ag, Ni, Mo and Ce are the effective desulphurisation sites for thiophene, BT, DBT and 4,6-DMDBT.
Although much experimental work has been done on HDS process on transition metal-supported zeolites, a detailed mechanism on the HDS over metal-supported zeolites remains unclear from a theoretical standpoint due to the complexity of HDS reaction system. Much related theoretical work has also been reported. Jayaraman pointed out that the presence of organic nitrogen compounds and polycyclic aromatic hydrocarbons on Cu(I)Y zeolite could not remove dibenzothiophene well through π-complexation [12]. In the work of Wang et al., the Cu(I)-Y had the best adsorption capacity on cationexchanged Y zeolites among the thiophene and BT [13]. Bian et al. investigated the oxidative desulphurisation mechanism of aromatic sulphides on Ti-MWW zeolite and found that thiophene had the highest oxidation activation energy than BT, DBT and 4,6-DMDBT [3]. If the Ce atoms are added into the Cu-Y zeolite, the Cu I Ce IV Y zeolite has a higher sulphur adsorption capacity and high selectivity for sulphur compounds such as BT, similar to Ce IV Y, as reported by Song et al. [14]. Jiao et al. investigated the HDS mechanism of BT on Mo 3 S 9 cluster catalyst. Liu et al. used the density functional theory (DFT) method to study the mechanism of thiophene HDS on a sulphide Co-Mo catalyst supported zeolite L [15]. In the work of Jiao and Liu et al., the desulphurisation mechanisms of BT and thiophene are similar but the energy barriers are different in the Mo 3 S 9 cluster and Co-Mo-L zeolite.
FAU zeolite is a typical industrial catalyst and widely used in HDS of thiophene. In this work, the thiophene HDS mechanism catalysed by Cu-Mo-supported FAU zeolite will be studied through a two-layer ONIOM method to understand the nano-confinement effect of the zeolite channel and the nature of the transition states. The Cu and Mo atoms were chosen because of their good performance in experimental HDS process. In this study, we investigated two interesting theoretical questions: (1) the rate-determining step for the HDS process over the Cu-Mo-supported FAU zeolite and (2) the relationship between the characteristic of thiophene desulphurisation and the topological structure of zeolite. It is expected that the theoretical results will provide useful insights into the experimental design of more efficient desulphurisation catalysts.

Computational methods
FAU zeolite is a faujasite-type zeolite that has two typical supercages interconnected by a 12-membered ring (12MR) window in the crystal unit cell. We constructed a 156 T FAU nanocluster from the original crystal structure [16]. A silicon (Si) atom is replaced by an aluminium (Al) atom. A molybdenum (Mo) atom connects with the Al atom through two bridging oxygen (O) atoms. A copper (Cu) atom connects with the Mo atom through two bridging sulphur (S) atoms. The Cu and Mo atoms with two bridging S atoms are placed around the [(≡SiO) 3 Al(OSi≡)] cluster in the Cu-Mo-supported FAU zeolite model. The corresponding zeolite model is displayed in Figure 1.
We carried out a two-layer ONIOM study. In the ONIOM scheme, the total energy of the entire system is calculated according to the following equation [17]: Here, the whole region contains an inner region and a larger outer region. High-level quantum mechanics (QM) calculations are needed for the inner region and reaction molecules. Molecular mechanics (MM) calculations are needed for the rest of zeolite framework. The MM energy term is used to treat the van der Waals and classical electrostatic interactions. The MM electrostatic interaction energy is obtained from the Coulomb's law, i.e. the classic charge-charge interaction. We used the DFT method with the M06 hybrid functional and the universal force field (UFF) to describe the inner and outer regions, respectively [18,19]. It is well-known that the B3LYP functional often underestimates the activation barrier height of the transition state and is not accurate enough to treat the medium-range VDW interactions [20]. The M06 functional developed by Zhao et al. exhibits a remarkable accuracy for predicting the main group thermochemistry, kinetics, non-covalent interaction, electronic excitation energy and aromatic-aromatic stacking interaction [21][22][23]. In particular, M06 functional can be used to deal with the problem of multi-reference rearrangement involving reactions in which both organic and transition metal bonds are formed or broken.
A zeolite model called ONIOM (M06 (qm1 (6-31G (d,p)): qm2 (3-21G)): UFF) is employed in order to save the computational time. The inner QM region can be divided into two parts: qm1 and qm2 subregions. The active centres of zeolite and reactants are included in the qm1 subregion and calculated by the 6-31G (d,p) basis set; the qm2 subregion includes the rest of inner QM region and treated by the 3-21G basis set. We employed the 6-31G (d,p) basis set for Si, Al, O, C, S and H atoms and the LANL2DZ effective core potential for Mo and Cu atoms. The use of the mixed basis sets can be found in many published works [24][25][26][27]. For FAU, the qm1 and qm2 subregions contain 12 and 6 T atoms, respectively. From the previous theoretical works, the 18 T QM size is reliable enough in the ONIOM model for studying the desulphurisation reaction network in the bulk catalytic system [28].
The atoms in the catalytic centre Cu-S-Mo cluster and reactant molecules were relaxed, and the zeolite framework was rigid in the geometry optimisation. The intrinsic reaction coordinates (IRCs) were calculated to verify the transition states [29]. The typical experimental temperature of 673.15 K is used for energy thermal corrections. In order to obtain more accurate energies, we carried out single-point energy calculations for QM region at M06 (qm1(6-311+G(2df,2p)):qm2 (6-31G(d,p))) level.
The classical electrostatic interaction energy in the force field calculations is calculated by the Coulomb's law, i.e. the classic charge-charge interaction that is denoted as follows: where the q i is the force field atomic charge of atom i. The van der Waals energy in the UFF force field is obtained according to the Lennard-Jones 12-6-type expression as [19]: where 1 ij is the well depth and r 0 ij is the VDW bond length. 1 ij and r 0 ij are calculated from a geometric mean combination rule as follows: 1 i and r 0 i parameters for atom i can be taken from the UFF force field.
The difference in charge densities (DCDs) [30,31] was employed to describe the migration of electrons in a transition state structure as the following equation: where ρ zeolite-organicfragment , ρ organicfragment and ρ zeolite are the electron densities of the TS, organic fragment and isolated zeolite in the TS, respectively. The isosurface plots of reduced density gradient (RDG) and Ω function are employed to distinguish and visualise different types of interactions of important intermediates and transition states in the zeolite system [32]. These two functions are given as follows: V(r) = Sign(l 2 (r))r(r) where ρ(r) is the total electron density, and λ 2 is the second largest eigenvalue of the Hessian matrix of electron density.
The electron density is further analysed by the local orbital locator (LOL). The LOL is defined as the following equation [33]: where t LSDA s , t exact s and ν σ are the LOL function, non-interacting Kohn-Sham kinetic energy density and local spin density approximation (LSDA), respectively.
In the catalytic reactions, the organic fragments in transition states in the zeolite channel are cationic and the zeolitic framework is anionic in nature. It is important to describe accurately the electrostatic interaction between the cationic transition state and the anionic framework of zeolite. In present study, the QM atomic charges were obtained from the DFT calculations and the MM charges were taken from our previous work [34]. Gaussian 09 and Multiwfn software packages were employed for the ONIOM calculations and the analysis of DCDs, RDGs and LOLs functions, respectively [35,36].

Results and discussion
Generally speaking, the sulphur vacancies on the sulphided Cu-Mo catalyst serve as the active centre for the HDS reaction of thiophene and its derivatives. As a promoter, Cu atoms promote the activation of molecular hydrogen and the formation of sulphur vacancies between Mo and Cu atoms [37,38]. The main hydrogenated species of thiophene are 2,3-dihydrothiophene (2,3-DHT), 2,5-dihydrothiophene (2,5-DHT) and tetrahydrothiophene (THT) [39,40]. According to the previous work [15], hydrogenolysis desulphurisation (HYD) and direct desulphurisation (DDS) are to be discussed (Schemes 1-5). In the HYD pathway, hydrogenation of C = C bond occurs before the C−S bond is broken. In this work, the HDS means HYD and H 2 molecules participate in the desulphurisation process. In the DDS pathway, the concerted desulphurisation of hydrogenated thiophene is investigated and the H 2 molecules do not participate in the desulphurisation process. All transition states (TSs) and intermediates are identified in the elementary steps. All important information of TSs are shown in the Figures 1-12 and Figures S1−S11. The free energy barriers of TSs are given in Tables 1-5. Schemes 1-5 display the HYD and DDS pathways. In the following sections, when we use the word 'FAU', it represents the Cu-Mo-FAU zeolite.

The formation of sulphur vacancies
Scheme 1 shows the pathway for the formation of sulphur vacancy on the FAU zeolite. The free energy profile is given in Figure S1 (Table 1), and the subsequent desulphurisation of hydrogenated thiophene is occurred on the sulphur vacancy. Liu et al. reported that the TS-2 needed to overcome a higher activation barrier than that in TS-1 in zeolite L [15]. In the present study, the calculated results show that the H−H bond cracking step (TS-1) is more difficult to occur than the formation of H 2 S (TS-2) due to the higher activation barrier over Cu-Mo-FAU zeolite.
The rate-determining step is very important for understanding the difficulty of occurrence of a reaction pathway. According to the Murdoch's work [42], we displayed the energy profile of the formation of sulphur vacancy in Figure  S1. When the whole reaction pathway is divided into several different sections, we calculated the energy difference ΔE i between the TS with the highest energy and the initial intermediate in each section i. The maximum (ΔE i ) corresponds to the rate-determining step. An effective method on how to divide the reaction pathway is that we first choose an intermediate INT(i) and then select another intermediate INT(j)  with lower energy than the INT(i) and the same number of atoms as the INT(i). This method is often used several times throughout the whole reaction pathway to obtain several different reaction sections; for the simple pathways, we probably divide the reaction pathway only once. In Figure S1, one section is involved: complex-1 to complex-2. The largest energy difference is obtained for TS-2 with the free energy barrier of 178.64 kJ mol −1 , showing that the rate-determining step is the formation of H 2 S. The free energy barrier of 178.64 kJ mol −1 is obtained from energy difference between TS-2 and complex-1. In addition, one can see that from Figure  S1, the formation of sulphur vacancy is a thermodynamically endothermic reaction and thus high reaction temperature is needed to produce enough sulphur vacancies.
Decomposing the ONIOM total energy barriers into QM and MM contributions is an important way to investigate the effect of zeolite framework on the reaction mechanisms. According to Equation (1), the total ONIOM free-energy barrier for a given reaction step from INT to TS can be calculated from the following equation [43]: The DDE MM value is called MM energy contribution to the total ONIOM free energy barrier and listed in Tables 1-5. The sign of DDE MM indicates that the zeolite framework is favourable (DDE MM . 0) or unfavourable (DDE MM , 0) for the formation of TS. In Table 1, one can see that zeolite framework is favourable for the formation of TS1 and TS2 because the DDE MM values are negative. On the other hand, the free energy barriers are controlled by the DE QM values because of the small DDE MM values.
The enthalpy barriers (ΔH ≠ ), entropy losses (−TΔS ≠ ), free energy barriers (ΔG ≠ ) at 673.15 K for TS1 and TS2 are also given in Table 1. It can be seen that these reaction steps are mainly dominated by enthalpy barriers because all entropy losses are very small. The entropy losses in TS1 and TS2 are negative, which suggests that the formation of sulphur vacancy is an entropy-increased pathway.

The desulphurisation of thiophene
According to the work of Liu et al., the HDS of thiophene (not hydrogenated thiophene) can occur through ring-opening and C−S bond cracking steps on Co-Mo-supported L zeolite [15]. We investigated carefully the HDS pathway of thiophene on Cu-Mo-supported FAU zeolite and did not find the corresponding transition states. Many experimental work has suggested that CuO/SBA-15 [44], Cu/ ZrO 2 [45], CuCeY [46,47], CuNaY [48] and CuY [49] can been used for adsorptive desulphurisation of thiophene. We optimised adsorbed structures of thiophene on Cu-S2-Mo and Cu-S-Mo FAU zeolites. On Cu-S2-Mo FAU, the C = C double bond of thiophene interacts with the Cu atom; the bond length of Cu … C = C is about 2.21 Å ( Figure S4a). The Cu … C = C structure is consistent with the reported experimental prediction for the η 5 adsorption mode of thiophene [50]. On Cu-S-Mo FAU, the sulphur atom of thiophene interacts with the Mo atom; the bond distance of S … Mo is 2.40 Å ( Figure S4b). The S … Mo adsorption model is in accord with the experimental prediction for the η 1 adsorption mode of thiophene [50]. At 0 K, the adsorption energies of thiophene on Cu-S2-Mo and Cu-S-Mo FAU zeolites are −100.21 and −146.19 kJ mol −1 , respectively. At 673.15 K, the adsorption energies of thiophene on Cu-S2-Mo and Cu-S-Mo FAU are −4.19 and −9.76 kJ mol −1 respectively, indicating that even if at high temperatures, thiophene still tends to adsorb on the zeolite. Therefore, we predict that the HDS pathway of thiophene is probably adsorptive desulphurisation in the absence of H 2 gas. In the presence of H 2 , thiophene will be firstly hydrogenated into 2,5-DHT, 2,3-DHT and THT and then desulphurised through hydrodesulphurisation (HDS) pathway. This conclusion can be obtained in the following sections.

The hydrogenation of thiophene to 2,5-DHT, 2,3-DHT and THT
In the HYD pathway, the C = C double bond of thiophene molecule is firstly hydrogenated and then the C−S bond is broken [39,50]. Zhu and Hensen et al. pointed out that the hydrogenated intermediates of thiophene are 2,5-DHT, 2,3-DHT and THT [39,40]. The investigation on hydrogenation process from thiophene to 2,5-DHT, 2,3-DHT and THT is crucial to better understand the HDS mechanism of thiophene because hydrogenation of C = C bond of thiophene is probably competitive with the HDS process. Scheme 2 gives the pathway of the hydrogenation of thiophene on the Cu-2S-Mo FAU zeolite.
The initial step is the adsorption of thiophene on the Cu-2S-Mo catalytic centre with a S … Mo adsorption mode (complex-trans-1). Subsequently, the Cu−H is breaking and a new C−H bond in thiophene is forming through the transition state TS-trans-1 (Scheme 2); the bond lengths of Cu−H and C−H are 1.44 and 1.65 Å, respectively (Figure 4(a)). The activation free energy barrier of the hydrogen transfer step is 152.73 kJ mol −1 ( Table 2). The complex-trans-2 undergoes two different hydrogen transfer steps through the transition states TS-trans-2 and TS-trans-4 to form complex-trans-3 and complex-trans-6, respectively. The activation free energy barriers for TS-trans-2 are calculated to be 49.95 kJ mol −1 , In the TS-trans-2 ( Figure  4 that the hydrogenation of thiophene to 2,5-DHT needed to overcome a moderate energy barrier of 93.50 kJ mol −1 on Mo 3 S 9 cluster catalyst [51]. Another hydrogen transfer step through TS-trans-4 needs an activation free energy barrier of 67.57 kJ mol −1 ( Table 2), suggesting that these two hydrogen transfer steps (TS-trans-2 and TS-trans-4) can occur easily on FAU zeolite. In TS-trans-4 ( Figure 4  One H 2 molecule is then chemically co-adsorbed with 2,5-DHT in complex-trans-3 to form complex-trans-4 or coadsorbed with 2,3-DHT in complex-tran-6 to form complextrans-7. Subsequently, THT could be produced from two different pathways. The second pathway for the formation of THT is two successive hydrogen transfer steps in TS-trans-5 and TS-trans-6 from complex-trans-7 and complex-trans-9, respectively (Scheme 2). In the TS-trans-5 (Figure 4(e)), the Cu−H bond is breaking and the Mo−H bond is forming and the 2,3-DHT (complex-trans-7) is converted into 2,5-DHT (complex-trans-8) with a moderate free energy of 134.47 kJ mol −1 ( Table 2). The 2,5-DHT in complex-trans-8 is then desorbed and the 2,3-DHT is adsorbed again (complex-trans-9). The 2,3-DHT in complex-trans-9 is subsequently converted into THT with concerted breaking of Mo−H and S−H bonds (TS-trans-6) and this step needs to overcome a small free energy of 67.84 kJ mol −1 . Therefore, the lower activation barrier from TS-trans-6 (67.84 kJ mol −1 ) than TS-trans-3 (182.97 kJ mol −1 ) indicates that the formation of THT occur more easily through the hydrogenation of 2,3-DHT. Jiao et al. reported that the formation of THT from the conversion of 2,5-DHT needed an energy barrier of 109.20 kJ mol −1 on Mo 3 S 9 cluster catalyst [51]. The corresponding structures of TS-trans-5 and TS-trans-6 are given in Figure 4 Table 2 lists the enthalpy barriers, entropy losses, QM and MM contributions for the free energy barriers at 673.15 K for all TSs. Several interesting viewpoints can be obtained. First, all MM energy values are positive in FAU, suggesting that FAU zeolite framework is unfavourable for the hydrogenation of thiophene to 2,5-DHT, 2,3-DHT and THT. Second, positive entropy losses values for all TSs indicate that the hydrogenation process of thiophene is entropy-decreased. Third, all the free energy barriers are controlled by the QM contributions except TS-trans-6. In TS-trans-6, the entropy loss is larger than the enthalpy barrier, showing the hydrogen transfer step in TS-trans-6 is controlled by entropic effect.

The HYD pathway of 2,5-DHT
In the HYD pathway, the reactant is the hydrogenated thiophene. As mentioned above, one or two C = C bonds of thiophene can be hydrogenated into three different derivatives: 2,3-DHT, 2,5-DHT and THT. Due to the low activation barriers, the hydrogenation process occurs easily on the zeolite. Experimental observations showed that the hydrogen transfer and the cleavage of C−S bond are the main steps in the HYD pathways of the three hydrogenated species [52][53][54].
Scheme 3 displays the HYD pathway of 2,5-DHT on the sulphided Cu-Mo-FAU zeolite. The corresponding energy profile is given in Figure 2. After the sulphur vacancy is formed (Scheme 1), the 2,5-DHT is firstly adsorbed on the sulphur vacancy. Then, one H 2 molecule is co-adsorbed with the 2,5-DHT through a Mo···H 2 interaction on Cu-S-Mo cluster to form a complex-2,5DHT-2. The H−H and C−S bonds in complex-2,5DHT-2 are broken through the transition state TS-2,5DHT-1 ( Figure 5(a)); the Gibbs activation free energy barrier is calculated to be 52.48 kJ mol −1 (Table 3), which indicates that C−S bond cleavage in 2,5-DHT can occur easily on FAU zeolite. Jiao et al. theoretically reported that the activation barrier for the ring-opening step of 2,5-DHT was 138.60 kJ   bond lengths of C−S, S−H and C−H are in the range of 1.58 −2.78 Å . After TS-2,5DHT-3, the S atom of the −SH group replenishes the S vacancy on the zeolite, and the catalyst is regenerated (complex-2,5-DHT-3). The predicted product 2butene was experimentally observed in situ FTIR and XANES studies for thiophene HDS on Ni 2 P/MCM-41 zeolite by Oyama et al. [50].
A concerted direct desulphurisation (DDS) pathway is proposed for the desulphurisation of 2,5-DHT without H 2 molecules (Scheme 3). The two alkyl C−S bonds of 2,5-DHT in the complex-2,5DHT-1 are simultaneously broken through the TS-2,5DHT-4 to form 1,3-butadiene (complex-2,5DHT-4). The activation free energy barrier is 98.54 kJ mol −1 . Jiao et al. studied the formation of 1,3-butadiene on Mo 3 S 9 cluster and the calculated activation barrier was 89.6 kJ mol −1 , which is slightly lower than our result. The structure of TS-2,5DHT-4 is displayed in Figure 5 We calculated the enthalpy barriers, entropy losses, QM and MM contributions for the free energy barriers at 673.15 K for all TSs in Table 3. Several interesting viewpoints can be obtained. First, positive entropy loss values for all TSs indicate that the HDS of 2,5-DHT is entropy-decreased. Second, most of MM energy values except TS-2,5DHT-2 are positive, suggesting that FAU zeolite framework is generally unfavourable for the formation of TSs in the HDS of 2,5-DHT. Third, all the free energy barriers are controlled by the QM contributions. Fourth, all the free energy barriers are controlled by the enthalpy barriers except TS-2,5DHT-3 and TS-2,5DHT-4; In these two TSs, the entropy losses are larger than the enthalpy barriers, showing the C−S bond cleavage step in TS-2,5DHT-3 and concerted elimination step of 2,5-DHT in TS-2,5DHT-4 are controlled by entropic effect.
INT-2,3DHT-1 is converted to 1-butene through hydrogen transfer (TS-2,3DHT-2) and C−S bond cracking steps (TS-2,3DHT-3). The TS-2,3DHT-2 is formed through the hydrogen transfer step of INT-2,3DHT-1. The corresponding activation energy barrier is 141.99 kJ mol −1 , which is nearly equal to the free energy of 141.14 kJ mol −1 for the hydrogen transfer step (TS-2,5DHT-2) in the HYD of 2,5-DHT. At TS-2,3DHT-2 ( Figure 6 The predicted product 1-butene was found in the experimental work in Hensen [40] and Oyama et al. [50]. From TS-2,3DHT-3 ( Figure S5a), we can see that the S−H and C −S bonds are breaking and a new C−H bond is forming. The bond lengths of S−H, C−S and C−H are 3.50, 3.69 and 1.38 Å, respectively. After TS-2,3DHT-3, the S atom of the S −H group replenishes the S vacancy on the zeolite, and the catalyst is regenerated (complex-2,3-DHT-3).
As mentioned above, a concerted DDS elimination pathway is easy to occur for the desulphurisation of 2,5-DHT. Liu et al. reported a single concerted elimination pathway for the desulphurisation of 2,3-DHT to form 1,3-butadiene in the Co-Mo zeolite L [15]. In present study, IRC calculations showed that this concerted transition state did not exist in the Cu-Mo FAU zeolite. In fact, Liu et al. pointed out that the concerted elimination pathway in the Co-Mo zeolite L was unfavourable due to high activation barrier.
The enthalpy barriers, entropy losses, QM and MM contributions for the free energy barriers at 673.15 K for all TSs are listed in Table 4. First, positive entropy loss values for all TSs indicate that the HYD of 2,3-DHT is an entropy-decreased process. Second, most of MM energy values except TS-2,3DHT-2 are positive, suggesting that FAU zeolite framework is generally unfavourable for the formation of TSs in the HYD of 2,3-DHT. This prediction is the same as that for TS-2,5DHT-2. Third, all the free energy barriers are controlled by the QM contributions or the enthalpy barriers.  Species In the energy profile of HYD of 2,3-DHT ( Figure S2), two sections are involved: complex-2,3DHT-2 to INT-2,3DHT-1 (164.72 kJ mol −1 , TS-2,3-DHT-1) and INT-2,3DHT-1 to complex-2,3DHT-3 (253.54 kJ mol −1 , TS-2,3DHT-3). The largest energy difference is obtained in the second section (INT-2,3DHT-1 to complex-2,3DHT-3, 253.54 kJ mol −1 ) which shows that the rate-determining step is not the ring-opening and hydrogen transfer steps but C−S bond cracking of 2,3-DHT (TS-2,3DHT-3). It should be noted that for the same TS-2,3DHT-3, the energy barrier of 253.54 kJ mol −1 is different from the energy barrier of 138.15 kJ mol −1 . The reason is that the former is calculated from the energy difference between TS-2,3DHT-3 and the first intermediate INT-2,3DHT-1 in the second section in the reaction sequence but the latter is calculated from the energy difference between TS-2,3DHT-3 and its previous intermediate INT-2,3DHT-2. The former is an apparent energy barrier for the whole HYD pathway of 2,3-DHT and the latter is only an energy barrier for an individual step, i.e. the C-S bond cleavage. In addition, one can see that from Figure S2, the HYD of 2,3-DHT is a thermodynamically exothermic reaction.
3.6. The HYD pathway of THT Scheme 5 describes the HYD pathway of THT on the Cu-Mo FAU zeolite. The energy profile is given in Figure S3. As shown in Scheme 5, THT is adsorbed on sulphur vacancy (complex-THT-1). Adding four H atoms to thiophene will lead to the complete destruction of the π-bond system of the thiophene ring; all C−S and C−C bonds in THT are elongated. The H 2 and THT molecules are co-adsorbed on the Cu-S-Mo centre to form complex-THT-2. At the TS-THT-1 (Figure 6(c)), the H−H (0.83 Å) and C−S bonds (2.61 Å) are breaking and new C−H bond is forming (1.75 Å). The activation-free energy barrier is found to be 174.24 kJ mol −1 (Table 5), and this is a higher energy barrier than those found in 2,3-DHT (164.72 kJ mol −1 ) and 2,5-DHT (52.48 kJ mol −1 ) which indicates that the ring-opening step is unfavourable for the fully saturated thiophene. Jiao et al. studied the ring-opening of THT on Mo 3 S 9 cluster and reported that the calculated activation barrier was 170.6 kJ mol −1 [51]. Guo et al. studied the HYD of THT on Pt(111) and found that the energy barrier for the ring-opening of THT was 164.99 kJ mol −1 [39]. Their energy barriers are close to our result (174.24 kJ mol −1 ).
TS-THT-2 and TS-THT-3 undergo hydrogen transfer and C−S bond cracking steps, respectively, to produce a butane molecule. From TS-THT-2 ( Figure 6(d)), it can be seen that the Cu−H bond is broken and a new S−H bond is formed; the bond lengths of S−H and Cu−H are 1.52 and 2.18 Å, respectively. The activation-free energy barrier for TS-THT-2 is 111.33 kJ mol −1 . After the C−S bond cracking step from TS-THT-3, the butane molecule is formed and the catalyst is regenerated. The corresponding activation energy barrier is found to be 267.34 kJ mol −1 , which is higher than energy barriers of 138.15 kJ mol −1 in TS-2,3DHT-3 and 38.03 kJ mol −1 in TS-2,5DHT-3, indicating that the C−S bond cracking step is more difficult to occur in the fully saturated thiophene. Guo and Jiao reported different results for the C−S bond cracking of THT on Pt(111) (111.93 kJ mol −1 ) [39] and Mo 3 S 9 cluster (98.8 kJ mol −1 ) [51]. The final product butane was observed in the experimental work in Hensen [50] and Oyama et al. [40]. From TS-THT-3 ( Figure S5b On Co-Mo zeolite L, Liu et al. reported a concerted DDS elimination pathway for the THT but the concerted transition state is unfavourable due to high activation barrier [15]. In the current study, we obtained a concerted transition state in the Cu-Mo FAU zeolite that is structurally similar to the abovementioned concerted TSs in the work of Liu et al. However, the corresponding IRC calculations verified that the obtained TS in the current work did not exist. Therefore, we predict that the desulphurisation of THT should occur by the stepwise HYD pathway.
Two viewpoints can be obtained from the free energy barrier decomposition in Table 5. First, all the free energy barriers are controlled by the enthalpy barriers or QM contributions. Second, the entropy loss and MM energy value in TS-THT-3 are both negative, which indicates that the C−S bond cleavage step is entropy-increased and the FAU framework is favourable for the formation of TS-THT-3.
In the energy profile ( Figure S3 showing that the rate-determining step is the C−S bond cracking of THT. The hydrogen transfer of INT-THT-1 is much easier than the ring-opening and C−S bond cracking steps due to low energy barrier (111.33 kJ mol −1 ). In addition, one can see that from Figure S3, the HYD of THT is a thermodynamically exothermic reaction. In summary, for the HYD pathways of three hydrogenated thiophenes, the rate-determining step is the hydrogen transfer step for 2,5-DHT and C−S bond cracking step for 2,3-DHT and THT. Table 5. Calculated enthalpy barriers (ΔH ≠ ), entropy losses (−TΔS ≠ ), free energy barriers (ΔG ≠ ) at 673.15 K, QM and MM contributions to the free energy barriers for TS-THT-1, TS-THT-2, and TS-THT-3 for the HDS of THT over FAU zeolite. In this work, the organic fragments in the TSs are cationlike and the zeolitic framework is anionlike. The charges will flow among different atoms and the electron density will change when the reaction steps occur, which reflects the nature of the elementary steps. We calculated the difference charge densities (DCDs), the reduced density gradients (RDGs) and localised orbital locators (LOLs) of the transition states to elucidate the migration of electrons for different elementary steps.

The nature for the formation of sulphur vacancy
The DCD plot of TS-1 for the H 2 cleavage is obtained from the following equation: where ρ TS is the electron density of the TS. ρ TS zeolite hydrocarbon (noH2) and ρ TS H2(gas) are the electron densities of the TS without the H 2 fragment and isolated gas H 2 molecule at their equilibrium structures of the TS. After H−H bond cracking, the INT-1 is formed. A DCD analysis (Figure 7(a,b)) suggests that for TS-1, the electron density decreases in the H 2 fragment and increases in the Cu and S atoms, which reflects a clear electron transfer from the H 2 molecule to the Cu and S atoms when the H 2 molecule is cracked. Figure 7(c,d) show the LOL colour-filled map of TS-1 on FAU. As shown in Figure 7(c,d), the blue circle represents the region with low electron density between the valence shell and the inner shell. A yellow region between two hydrogen atoms indicates that there is a weak interaction between them and the H−H bond is forming or breaking, which is consistent with the nature of the transition state TS-1 (i.e. H−H bond cracking step).
The RDG plot for TS-1 over the FAU zeolite is given in Figure 7(e). A red region appears between two hydrogen atoms in the TS-1, which indicates that there is spatial repulsion between these two hydrogen atoms because they are close to each other.
The DCD plots for TS-2 are given in Figure 8(a,b). A careful DCD analysis suggests that the electron density decreases in the H 2 S fragment, reflecting that the electron cloud is migrating from the organic fragment H 2 S to the catalytic centre during the formation of H 2 S molecules. A green LOL region appears between H 2 S and Cu atom for TS-2 in Figure  8(c,d), indicating that there is a weak attraction interaction between them. This weak interaction is consistent with the nature of the transition state for the formation of H 2 S.
The RDG plot for TS-2 over the FAU zeolite is given in Figure 8(e). A red and blue region appears between the H 2 S fragment and catalytic site; the red region represents a spatial repulsion and the blue region indicates a weak Mo … H 2 S coordination attraction interaction. A careful comparison of RDGs between TS-1 (Figure 7(e)) and TS-2 (Figure 8(e)) suggests that there is a wider red region in TS-1 than in TS-2, which indicates a stronger spatial hindrance and thus higher free energy barrier in TS-1 (178.62 kJ mol −1 ) than TS-2 (113.62 kJ mol −1 ). 3.7.2. The nature for the HYD of 2,5-DHT The DCD plots of TS-2,5DHT-1 are given in Figures S6a and S6b. The electron density on the organic framework decreases but the electron density in the catalytic centre increases, which indicates that in the ring-opening process of 2,5-DHT, the electrons migrate from the organic chain to the catalytic centre. The RDG plot for TS-2,5DHT-1 is given in Figure  S6c. There is a red region between the organic chain and catalytic centre, showing a strong ring tension in the TS-2,5DHT-1. Figures S6d to S6f show the colour-filled LOL map of TS-2,5DHT-1. In the LOL map, there is a red region between C atom of the organic fragment and H atom of H 2 fragment, which shows that the C−H covalent bond is nearly formed. A weak yellow region appears between two H atoms, indicating that there is a weak attraction interaction between them.
The DCD plots of TS-2,5DHT-2 are given in Figures S7a and S7b. The electron density on the H atom and the organic fragment decreases and the electron density on the Cu atom increases, which reflect that the electrons migrate from the organic fragment and the H atom to the Cu atom in the hydrogen transfer process. Figures S7c and S7d show the colour-filled LOL map. In the LOL map, a yellowish-green region is formed between the migrating H atom and the Cu atom, showing a weak interaction between H and Cu atoms. A white region appears between the migrating H atom and the sulphur atom in the organic fragment, which indicates that the H−S covalent bond is nearly formed in the TS-2,5DHT-2. The RDG plot for TS-2,5DHT-2 is given in Figure S7e. There is a blue and red region between organic fragment and Cu-S-Mo catalytic centre, showing that both steric hindrance and a weak van der Waals attractive interaction exist in the reactive centre.
The DCD plots of TS-2,5DHT-3 are given in Figure 9(a,b). The electron density of the Cu-S2-Mo centre decreases but the electron density of the thiophene fragment increases, indicating that for the formation of 2-butene, the electron cloud migrates from the Cu-S2-Mo centre of the zeolite to the thiophene fragment. Figure 9(c-e) shows the colour-filled LOL maps of TS-2,5DHT-3. In the LOL map, there is a red LOL region in the H … C atomic pair, which indicates a strong covalent interaction in the H … C region in the TS-2,5DHT-3. The RDG plot is given in Figure 9(f). There is a red region between thiophene fragment and the Cu-S2-Mo catalytic centre, which indicates a strong steric hindrance between them. Although the short bond distances in the S … H … C region lead to steric hindrance, the short S … H … C bond distance increases the degree of covalent interaction for the ones such as the H … C bond, which is consistent with the prediction of LOL.
The DCD plots of TS-2,5DHT-4 are given in Figure 10(a,b). The electron density of Cu-S2-Mo in the catalytic centre of zeolite decreases, and the electron density on the carbon chain atoms of thiophene ring increases, suggesting that the electron cloud migrates from the Cu-S2-Mo centre to the organic carbon chain. Figure 10(c,e,f) shows the colour-filled LOL map of TS-2,5DHT-4. The LOL map describes the atomic interactions of organic fragment and Cu-S2-Mo catalytic centre. The RDG plot for TS-2,5DHT-4 is given in Figure 10 (d). A yellowish-green region appears between the Cu-S2-Mo catalytic active centre and the organic carbon chain atoms, indicating that there is a weak attraction between them. The yellow and green region is accompanied by a little red region, showing that there is weak steric hindrance between the catalytic centre and the organic carbon chain. The weak steric hindrance is consistent with the low energy barrier of TS-2,5DHT-4 (98.54 kJ mol −1 ). Figure 11(a,b) for TS-2,3DHT-1 and Figure 12(a,b) for TS-THT-1 show the DCD plots. We can see that the electron density on the thiophene ring decreases and the electron density on the catalytic centre increases, indicating that the electrons migrate from the thiophene ring to the catalytic centre in the ring-opening process of 2,3DHT and THT. Figure 12(c−e) show the colour-filled LOL map of TS-THT-1. There is a yellow LOL region between Cu atom, H 2 molecule and C atom of organic carbon chain, showing that there is a weak attractive interaction between them.

The nature for the HYD of 2,3-DHT and THT
The RDG plots are given in Figure 11(e) (TS-2,3DHT-1) and 12f (TS-THT-1). In TS-2,3DHT-1 and TS-THT-1, there is a red region between the thiophene ring and catalytic active centre, showing a strong ring tension in these two TSs. A comparison of their RDGs suggests that the RDG plot of TS-2,3DHT-1 is similar to that of TS-THT-1, which is consistent with similar free energy barriers between TS-2,3DHT-1 (164.72 kJ mol −1 ) and TS-THT-1 (174.24 kJ mol −1 ).
The DCD plots predict similar behaviours on the migration of electron for TS-2,3DHT-2 ( Figures S8a and S8b) and TS-THT-2 (Figures S10a and S10b). The electron density on the H atom and the organic fragment decreases, while the electron density on the Cu and Mo atoms increases, which reflects the migration of electrons from the organic fragment to the Cu-Mo catalytic centre. The RDG plots for TS-2,3DHT-2 and TS-THT-2 are given in Figures S8c and S10c, respectively. There are blue and red regions between the organic carbon chain and the Cu-S-Mo catalytic active centre, which verified the existence of steric hindrance and weak van der Waals attraction interaction. The colour-filled LOL maps of TS-2,3DHT-2 ( Figures S8d−S8f) and TS-THT-2 (Figures S10d −S10f) show the similar nature of TSs. The yellow and green region between H atom and Cu atom indicates that there is a weak interaction in the H … Cu atomic pair. A white area appears between the migrating H and sulphur atoms, which shows that the H−S covalent bond is nearly formed in TS-2,3DHT-2 and TS-THT-2.
The DCD plots are given in Figures S9a and S9b for TS-2,3DHT-3 and S11a and S11b for TS-THT-3. The electron density of the Cu-S2-Mo active centre decreases, but the electron density of the organic fragment increases. Therefore, when 1-butene or butane is formed, the electron cloud migrates from the Cu-S2-Mo active centre to the organic fragment. The RDG plots for TS-2,3DHT-3 and TS-THT-3 are given in Figures S9e and S11c, respectively. There is a red region between the organic chain and the Cu-S2-Mo catalytic centre, which indicates that there is a strong steric hindrance between them. The LOL maps in Figures S9c and S9d and Figures S11d−S11f display a red and yellow LOL region in the S … H … C region, which suggests that there is a strong covalent interaction in the S … H … C region for the TS-2,3DHT-3 and TS-THT-3.

Conclusions
In this work, the catalytic ability of the Cu-Mo-supported FAU zeolite on HDS process of the hydrogenated thiophene is investigated by a two-layer ONIOM method. The nature of TSs and the confinement effect of zeolitic nanospace are revealed in detail. The major conclusions are summarised as follows: (1) Thiophene can be hydrogenated to 2,5-DHT, 2,3-DHT and THT with moderate free energy barriers. The sulphur atoms in the hydrogenated thiophene derivatives are removed as H 2 S molecules and 1-butene, 2-butene, 1,3butadiene and butane are produced. The rate-determining steps are the hydrogen transfer step for 2,5-DHT (TS-2,5DHT-2) and C−S bond cracking step for 2,3-DHT (TS-2,3DHT-3) and THT (TS-THT-3). The calculated free energy barriers increase in the following order: (TS-2,5DHT-2, 141.14 kJ mol −1 ) < (TS-2,3DHT-3, 253.54 kJ mol −1 ) < (TS-THT-3, 267.34 kJ mol −1 ). From the energy profile, the HYD of 2,5-DHT, 2,3-DHT and THT are thermodynamically exothermic reactions and the formation of sulphur vacancy is endothermic reaction. The concerted DDS pathway is probably more favourable than the HYD pathway in the desulphurisation of 2,5-DHT. (2) We can obtain useful information from the calculated enthalpy barriers, entropy losses, QM and MM contributions for the free energy barriers. Overall, all the free energy barriers are controlled by the enthalpy barriers and QM contributions except the C−S bond cracking (TS-2,5DHT-3), the formation step of THT (TS-trans-6) and concerted elimination step (TS-2,5DHT-4) of 2,5-DHT. The hydrogenation of thiophene to 2,5-DHT, 2,3-DHT and THT and their HDS processes are entropy-decreased; the formation of sulphur vacancy is entropy-increased. FAU framework is favourable for the formation of sulphur vacancy and generally unfavourable for the hydrogenation of thiophene and HYD process of 2,5-DHT, 2,3-DHT and THT. (3) We studied the nature of transition states by the analysis from RDG, DCD, and LOL functional. The RDG plots indicate that both steric hindrance and a weak van der Waals attractive interaction exist between organic fragment and catalytic centre for all TSs. The steric hindrance is probably responsible for the calculated high activation free energy barriers in the HYD process. The DCD analysis reveals that for the ring-opening process of 2,5-DHT, 2,3-DHT and THT, the electrons are migrated from the organic chain to the Cu-Mo catalytic centre; for the removal of S atom, the electrons migrate in the opposite direction. The LOL maps for all HYD TSs suggest that there are strong covalent interactions and weak van der Waals interactions between the atoms in the forming chemical bonds.

Supplementary Material
Free energy profiles and relative energies for the formation of sulphur vacancy and HYD pathways of 2,3-DHT and THT, Optimised structures of adsorbed thiophene molecule, TS-2,3DHT-3, and TS-THT-3, Plots of DCDs, RDGs, and LOLs