Molecular simulation on mechanism of thiophene hydrodesulfurization on surface of Ni 2 P

Abstract

Hydrodesulfurization reaction, as the last step of hydrothermal cracking reaction, is of great significance for the reduction of viscosity and desulfurization of heavy oil. Based on Density Functional Theory and using Dmol3 module of Materials Studio, this research simulated the adsorption and hydrodesulfurization of thiophene on Ni₂P (001) surface, and discussed the hydrodesulfurization reaction mechanism of thiophene on Ni₂P (001) surface. It was found that the direct hydrodesulfurization of thiophene had more advantages than the indirect hydrodesulfurization of thiophene. Finally, the optimal reaction path was determined: C₄H₄S+H₂→C₄H₆.

Keywords

density functional theory thiophene adsorption

Introduction

With the increasing demand for oil resources, thin oil field exploitation of China has entered the late stage of development, and the production of thin oil is becoming increasingly scarce. Therefore, the exploitation of heavy oil has become the main means to make up for the shortage of production. However, the viscosity of heavy oil is very high, and there are a large number of heteroatoms in its structure, which makes it difficult to exploit and transport heavy oil. The hydrothermal cracking reaction that occurs in the process of steam injection recovery of heavy oil can crack a part of the polymer components to produce light oil and improve the quality and reduce the viscosity, which has attracted the attention of researchers. The last and most critical step of the hydrothermal cracking reaction is hydrodesulfurization. The hydrodesulfurization reaction is to break the C–S bond and then form hydrogen sulfide with the hydrogen in the system, which reduces the viscosity and removes the sulfur element. With the development of technology, the industrial demand for heavy oil desulfurization is becoming more and more strict, so it is particularly important to further reduce the sulfur content by deep hydrodesulfurization (HDS) (Bataille et al., 2000) of heavy oil components. In order to achieve the effect of deep desulfurization of heavy oil (Angelici, 2006), it is necessary to reduce the steric hindrance effect of some sulfur-containing compounds by hydrogenation reaction, but the traditional catalyst has poor catalytic effect and weak hydrogenation ability, so it is difficult to achieve the effect of deep desulfurization (Sun and Prins, 2008). Therefore, researchers from various countries are committed to developing an HDS catalyst with excellent desulfurization effect (Oyama, 2003) to achieve the purpose of deep desulfurization of heavy oil components.

At present, transition metal elements and their sulfides have been widely used, such as Mo, Ni, Co, etc. (Chu et al., 2013; Li et al., 2015; Song et al., 2014), which are common desulfurization catalysts nowadays. Recently, carbides, nitrides, and phosphides have appeared in the literature because of their good physical properties and excellent hydrodesulfurization characteristics. Such catalysts have not been widely used in industry, but their application has broad prospect. Okamoto et al. (Okamoto et al., 2003) prepared precious metal catalysts such as Pd and Pt, and their experimental studies showed that such catalysts had better hydrodesulfurization effects than traditional transition metal catalysts. Vasudevan et al. (Jayamurthy and Vasudevan, 1994) found that the Pt catalyst had better hydrodesulfurization activity than the Co-Mo/Al₂O₃ catalyst system. Diana et al. (Stephanie et al., 2003) synthesized a supported catalyst and found that the hydrodesulfurization activity of a nickel phosphide catalyst supported on the surface of silica was better than that of a conventional sulfide catalyst. Maccra et al. (McCrea et al., 1997) synthesized a β-Mo₂C catalyst supported on the surface of alumina (A1₂O₃), and found that its effect on thiophene hydrodesulfurization was better than that of traditional transition metal catalysts. Puello-Polo et al. studied the catalytic effect of Ni-Mo carbide alloy catalyst in the hydrodesulfurization of thiophene and found that its catalytic performance was better. In the deep hydrodesulfurization reaction, the stability and activity of transition metal phosphides are exceptionally high (Clark et al., 2002; Oyama et al., 2004; Phillips et al., 2002). Moreover, compared with commercial catalysts, transition metal phosphides have lower oxygen consumption and better hydrogenation selectivity. Therefore, transition metal phosphide is one of the research hotspots of new catalytic materials.

There are a variety of organic sulfur compounds in heavy oil. Because of the typical structural characteristics of thiophene and the difficulty of desulfurization, thiophene is generally used as a model molecule for research (Zhu et al., 2006). The development of industrial hydrodesulfurization (HDS) catalysts has roughly experienced transition metal sulfides, noble metals, transition metal carbides and nitrides, and transition metal phosphides. Among them, the transition metal phosphides, especially nickel phosphide (Ni₂P), have excellent hydrodesulfurization activity, sulfur poisoning resistance and stability. Besides, transition metal phosphides have better hydrogenation selectivity and lower oxygen consumption than commercial catalysts, which is an important direction of hydrodesulfurization research at home and abroad (Song et al., 2012). Because the nickel phosphide crystal does not have a layered structure like ordinary sulfides, the activation edge and activation angle on the surface of the nickel phosphide crystal have a larger pore space, which can provide enough space for the reaction to promote the reaction (Hua et al., 2013). In addition, nickel phosphide has the characteristics of an approximate ellipsoidal structure, which makes it easier for nickel phosphide to form isotropic crystals and expose more active centers in the reaction process to promote the forward reaction (Nag et al., 1979). However, few studies are focused on exploring the viscosity reduction reaction mechanism by transition metal at the atomic. Therefore, the phosphide Ni₂P is chosen as the catalyst in this study.

In this paper, the elementary reaction of heavy oil model molecules on catalyst surface HDS was studied theoretically at the atomic by combining the quantum chemistry method in computer simulation technology. Firstly, the MS (Materials Studio) software was used to simulate the heavy oil model molecules and the Ni₂P catalytic system, then the thiophene adsorption and HDS processes on the Ni₂P (001) surface were obtained. Secondly, the most likely reaction path for thiophene hydrodesulfurization on the Ni₂P (001) surface was determined by calculating the structural parameters, adsorption energy and activation energy of the adsorption system. This study will more accurately get the theoretical reaction mechanism of thiophene hydrodesulfurization on Ni₂P (001) surface and contribute to the further development of deep HDS.

Calculation method

In this paper, combining with the quantum chemistry method in the computer simulation technology, the Dmol3 module in the Materials Studio software package was used for the calculation. Based on the density functional theory (DFT) numerical calculation method, the thiophene adsorption and hydrodesulfurization reactions were simulated. The Generalized Gradient Approximation (GGA) and PW91 density functional method (GGA-PW91) were used to calculate the exchange correlation energy. In the calculation, the double numerical basis set plus orbital polarization function (DNP) was used to process the valence electron wave function, and the electronic structure was analyzed by Density of States (DOS). And LST/QST tools were used to search for transition states. The Density Functional Semi-Core Pseudopotential (DSPP) was used for the electrons in the core of a metal atom. While for other atoms, such as hydrogen, carbon, nitrogen, and sulfur atoms, All Electrons were used. All calculations used spin polarization. Considering the calculation efficiency and calculation accuracy comprehensively, the system accuracy was set to Medium accuracy, the track thermal occupation generally adopted the default 0.005 Hartree, and the Density Mixing was set as 0.17. The thermal occupancy value was compatible with the selection of the Density Mixing damping coefficient and they were suitable for the value calculated this time. The iterative convergence criteria for energy, gradient, displacement, and Self-Consistent Field (SCF) under Medium accuracy were 22 × 10⁻⁵ Hartree, 4 × 10⁻³ Hartree/Å, 5 × 10⁻³Å and 1 × 10⁻⁵ Hartree respectively.

Surface energy was a concept to characterize the degree of chemical bond destruction between molecules when creating material surfaces. In the theory of solid physics, atoms on the surface of a substance had more energy than atoms inside it. According to the principle of lowest energy, atoms would spontaneously move from the surface of a substance to the interior of a substance. Therefore, another definition of surface energy was the extra energy on the surface of a material relative to the inside of the material. Its definition formula was equation (1)

E_{surf} = (E_{slab} - E_{bulk}) / 2 S

(1)

where

E_{slab}

was the energy of the surface model,

E_{bulk}

was the energy of a volume model with the same number of atoms,

S

was the surface area,

E_{surf}

was the surface energy.

The smaller the surface energy calculated by the formula, the more stable the corresponding surface was.

Different adsorption methods corresponded to different adsorption energies. The definition of adsorption energy was shown in formula equation (2)

E_{ads} = E_{adsorbate} + E_{surf} - E_{adsorbate / surf}

(2)

where

E_{adsorbate}

the energy of the molecules adsorbed in the gas phase,

E_{surf}

was the energy of the base, which was the energy

E_{slab}

of the surface model mentioned in equation (1),

E_{adsorbate / surf}

was the total energy of the surface adsorption system, Ha；The conversion factor was 1Ha = 627.5087 kcal/mol.

According to the above definition, a positive adsorption energy indicated stable adsorption, and the higher the adsorption energy value, the more stable the adsorption.

Simulation content

Construction and optimization of model molecules

Ni₂P was a metal-rich phosphide, and its crystal had a hexagonal structure with unit cell parameters were $P \overline{6} 2 m$ , a = b = 0.5859 nm, c = 0.3382 nm (Ankudinov et al., 1998; Rundqvist, 1961). Figure 1(a) showed a molecular configuration constructed based on the thiophene structural composition, and Figure 1(b) showed a periodic structure model constructed based on the Ni₂P unit cell parameters. The structure of thiophene molecule and Ni₂P single cell was optimized by the above method. The calculated parameters of the Ni₂P lattice were a = b = 5.9507 Å and c = 3.3972 Å(1 nm = 10Å). The optimized bond length was shown in Table 1. The energy of the optimized thiophene was E_thiophene = –15048.33 eV, and the energy of nickel phosphide was E_bulk = –274184.59 eV.

Figure 1.

Schematic of thiophene and Ni₂P.

Table 1.

The band lengths of Ni₂P and thiophene model molecules before and after optimization.

Value	Bond type	Exptl/Å	GGA + PW91/Å	Rel.error (%)
Ni₂P	Ni(1)-P(1)	2.226	2.299	0.13%
	Ni(1)-P(2)	2.209	2.240	1.40%
	Ni(2)-P(1)	2.369	2.401	1.35%
	Ni(2)-P(2)	2.457	2.488	1.26%
Thiophene	C²-S	1.730	1.729	–0.10%
	C²-C³	1.370	1.373	0.30%
	C³-C⁴	1.420	1.422	–0.60%

Selection of Ni₂P model adsorption surface

Ni₂P crystal had a hexagonal crystal structure, its (001) crystal surface had the largest interplanar spacing and the largest surface density. The (001) crystal surface was more easily broken and exposed to the catalyst surface (Fuks et al., 2010; Nelson et al., 2006; Oyama et al., 2009), so the (001) surface was selected as the adsorption surface in this paper. Since Ni₂P was a close-packed form of ABAB, there were two (001) crystal surfaces, as shown in Figure 2. The cut planes were shown in a 3 × 3 periodic model.

Figure 2.

Tangential sketch of Ni₂P (001) crystal surface (a), the structure is composed of Ni₃P (b), and the structure is composed of Ni₃P₂ (b).

In this paper, a periodic 3 × 3 four-layer flat plate model was chosen to simulate two kinds of surfaces in the Ni₂P (001) direction. The surface relaxation was considered in the calculation process, and the bottom atoms were fixed so that the upper three layers could move freely.

The surface relaxation optimization of two crystal surfaces Ni₃P and Ni₃P₂ was carried out respectively to obtain the total energy E_slab of the substrate. Figures 3 and 4 were the optimized energy diagram and the convergence diagram, respectively.

Figure 3.

Energy diagram during optimization process of the section Ni₃P₂ and Ni₃P.

Figure 4.

Convergence curve during optimization process of the section Ni₃P₂ and Ni₃P.

Figure 3 showed that the energy value gradually decreased and tended to be stable during the optimization process. The horizontal line in Figure 4 represented the standard requirements of the iterative convergence. When the optimization was completed, the energy diagram and the convergence indicator diagram appeared when the convergence condition was satisfied. The surface energy of the two surfaces was calculated according to the calculation formula of the surface energy and the data obtained by optimization. The results were shown in Table 2.

Table 2.

Surface energy of two kinds of crystal surfaces on Ni₂P(001).

Crystal surface	E_slab/Ha	E_bulk/Ha	S/Å²	E_surf/Ha/Å²
Ni₃P₂	–182206.6954914	–10076.0903057 × 18	17.8524 × 17.8524	1.313223234
Ni₃P	–182206.8181559	–10076.0903057 × 18	17.8524 × 17.8524	1.313415674

As seen in Table 2, the surface energy of the Ni₃P₂ cut plane was about 0.0001924Ha/Å² smaller than the surface energy of the Ni₃P cut plane. The conversion factor 1Ha = 27.2114 eV was used to convert it into an electron volt eV of about 0.005 ev/A², that is, the Ni₃P₂ cut plane was more stable than the Ni₃P cut plane. Therefore, this paper adopted Ni₃P₂ to represent the crystal surface in the direction of Ni₂P (001) as the research model.

Results and discussion

Adsorption of thiophene on Ni₂P (001) surface

The thiophene molecule was an axisymmetric structure, and the C atoms in the thiophene ring were numbered, as shown in Figure 5. The atom labeling method of the thiophene ring in Figure 5 was also applicable to its hydrogenated derivative.

Figure 5.

Schematic diagram of thiopene modecule. The letter X in C^X represents the position of C atom in the thiopene ring.

For a certain surface, there were usually four adsorption sites for the adsorbed molecules, namely the top site, the bridge site, and the triple vacancy (hcp site and fcc site). Since Ni₂P was a close-packed form of ABAB, its corresponding triple vacancy was only hcp site. There were two types of thiophene adsorption: parallel adsorption (flat mode) and vertical adsorption (end-connection mode).

An adsorption configuration was constructed on the Ni₃P₂ four-layer surface model. According to the relative position of the adsorption position, the adsorption of thiophene on the Ni₂P (001) surface could be divided into 12 adsorption configurations as shown in Figure 6, among which (a)∼(f) belonged to flat adsorption. While (g)∼(l) molecules were perpendicular to the surface, and belonged to the end-connection adsorption. In the calculation process, each surface model was only allowed to adsorb one adsorptive thiophene molecule, so it corresponded to a lower 1/9ML surface coverage. Table 3 showed various adsorption structures, adsorption energy and some important binding structure parameters.

Figure 6.

Adsorption structure of thiophene on the surface of Ni₂P (001).

Table 3.

Adsorption sites, adsorption energies E_ads (eV), and bond lengths (Å) for adsorbed thiophene derivatives on MoP(001).

Sites	E_ads	d(S-Ni)	d(C²-S)	d(C⁵-S)	d(C²-C³)	d(C³-C⁴)	d(C⁴-C⁵)
thiophene	–	–	1.729	1.729	1.373	1.422	1.373
top1	1.00	2.476	1.828	1.825	1.439	1.478	1.437
top2	1.04	2.365	1.840	1.838	1.463	1.470	1.466
top3	1.02	2.481	1.849	1.849	1.462	1.449	1.464
hcp1	1.06	2.233	1.785	1.785	1.413	1.448	1.412
hcp2	1.04	2.395	1.869	1.836	1.473	1.453	1.467
bridge	1.01	2.408	1.747	1.749	1.398	1.433	1.399
V-hcp	0.94	2.542	1.729	1.730	1.372	1.427	1.372

For (a)∼(f) adsorption configuration, the sulfur atoms in top1, top2, and top3 were located at the top position, the sulfur atoms in hcp1 and hcp2 were located at the triple vacancy, and the sulfur atoms in the bridge adsorption configuration were located at the bridge position. In (g)∼(l), the sulfur atoms were located at the triple vacancy in the V-hcp adsorption configuration. Compared with the end-connected mode adsorption configuration, the flat mode was more stable. Among the 12 adsorption configurations, the hcp1 configuration (E_ads=1.06 eV) with the molecular ring located above the hcp position in the flat mode was the most stable, followed by the top2 and hcp2 configurations (E_ads=1.04 eV). The adsorption energy of the other adsorption configuration was relatively low, and the adsorption energy of the V-hcp configuration was the lowest.

In the hydrogenation process, the most stable adsorption mode, that is, the hcp1 adsorption configuration was selected as the initial thiophene adsorption state. Hydrogen atoms were selected as the surface hydrogen source. An adsorption system was composed of thiophene adsorbed at the hcp site and H atoms adsorbed at adjacent vacancies, and used as the initial response state. The product structure formed by adding hydrogen atoms to different positions of thiophene was used as the final state of the reaction.

HDS reaction path of thiophene on Ni₂P surface

Six possible thiophene hydrodesulfurization reaction pathways were summarized as shown in Figure 7 (Ni et al., 2013). Tables 4 and 5 detailed the specific reaction steps of the six thiophene hydrodesulfurization reaction paths (* indicated a certain amount of charge, and the product was named according to the position α or β where H was added).

Figure 7.

Different reaction pathways for the hydrogenation of thiophene.

Table 4.

Reaction mechanisms for the hydrogenation of thiophene on Ni₂P (001) surface (A, B, C).

Step	Mechanism A	Mechanism B	Mechanism C
1	C₄H₄S+H→α-C₄H₅S+	C₄H₄S+H→α-C₄H₅S+	C₄H₄S+H→α-C₄H₅S+
2	α-C₄H₅S+H→α,α-C₄H₆S+	α-C₄H₅S+H’→α,β-C₄H₆S+	α-C₄H₅S+H→α,β-C₄H₆S+
3	α,α-C₄H₆S+H→α,α,β-C₄H₇S+	α,β-C₄H₆S+H→α,β,β-C₄H₇S+	α,β-C₄H₆S+H→α,β,α-C₄H₇S+
4	α,α,β-C₄H₇S+H→C₄H₈S+	α,β,β-C₄H₇S+H→C₄H₈S+	α,β,α-C₄H₇S+H→C₄H₈S+
5	C₄H₈S+H→C₄H₉S+	C₄H₈S+H→C₄H₉S+	C₄H₈S+H→C₄H₉S+
6	C₄H₉S+H→C₄H₁₀+S	C₄H₉S+H→C₄H₁₀+S	C₄H₉S+H→C₄H₁₀+S

Table 5.

Reaction mechanisms for the hydrogenation of thiophene on Ni₂P (001) surface (D, E, F).

Step	Mechanism D	Mechanism E	Mechanism F
1	C₄H₄S+H→β-C₄H₅S+	C₄H₄S+H→β-C₄H₅S+	C₄H₄S+H→C₄H₅S+
2	β-C₄H₅S+H→β,α-C₄H₆S+	β-C₄H₅S+H→β,α-C₄H₆S+	C₄H₅S+H→C₄H₆+S
3	β,α-C₄H₆S+H→β,α,β-C₄H₇S+	β,α-C₄H₆S+H→β,α,α-C₄H₇S+	/
4	β,α,β-C₄H₇S+H→C₄H₈S+	β,α,α-C₄H₇S+H→C₄H₈S+	/
5	C₄H₈S+H→C₄H₉S+	C₄H₈S+H→C₄H₉S+	/
6	C₄H₉S+H→C₄H₁₀+S	C₄H₉S+H→C₄H₁₀+S	/

The reaction products were given different names to distinguish different reaction paths, such as α-C hydrogenation and β-C hydrogenation. According to whether the C-S bond in thiophene breaks directly, the hydrodesulfurization reaction of thiophene was divided into two types: direct hydrodesulfurization reaction and indirect hydrodesulfurization reaction. The indirect hydrodesulfurization reaction was A-E, and the direct hydrodesulfurization reaction was F.

The adsorption systems consisting of catalysts, reactants and reaction products of the six reaction mechanisms in the reaction network diagram were optimized respectively. The optimized stable adsorption configurations were used as the initial state (IS) and final state (FS). Table 6 listed the activation energy E_a of each elementary reaction and the energy change ΔE before and after the reaction (Figure 8).

Table 6.

Activation energy (E_a, eV) and reaction energy (ΔE, eV) of each reaction on Ni₂P (001) surface.

Reaction	E_a	ΔE
C₄H₄S+H→C₄H₅S+	007	4.50
C₄H₅S+H→C₄H₆+S	0.41	−3.64
C₄H₄S+H₂→α-C₄H₅S+*	0.43	−0.57
C₄H₄S+H₂→β-C₄H₅S+*	0.89	−0.53
α-C₄H₅S+H→α,α-C₄H₆S+	1.92	0.22
α-C₄H₅S+H→α,β-C₄H₆S+	1.47	−0.23
α,β-C₄H₆S+H→α,β,β-C₄H₇S+	1.28	0.88
α,β-C₄H₆S+H→α,β,α-C₄H₇S+	1.14	−0.38
α,β,α-C₄H₇S+H→C₄H₈S+	2.38	0.12
C₄H₈S+H→C₄H₉S+	2.06	−0.39
C₄H₉S+H→C₄H₁₀+S	1.75	−1.25

Figure 8.

Direct hydrodesulfurization pathway of thiophene.

Reaction mechanism F was a direct hydrodesulfurization reaction. For the reaction C₄H₄S*+H₂*→C₄H₆*+S*, two C-S bond breaks occurred directly after hydrogenation, and the products were butadiene and sulfur atoms. When the first C-S bond break occurred at C²-S, the S atom moved toward the adjacent Ni atom and formed a bond. In this transition state (TS1), the C²-S bond was broken. After TS1, the intermediate product C₄H₅S was formed as the final state of the reaction. In this step, as the hydrogen atom approached, the C²-S bond changed from 1.785 Å single bond to 1.813 Å longer double bond, and the C-C double bond became single bond. The final state of the reaction was further elongated relative to the initial C²-S bond. The energy barrier of the reaction in this step was 0.07 eV, and the endotherm was 4.50 eV. The second C-S bond break occurred at the C⁵-S bond, and the reaction product was a co-adsorption system of C₄H₆ intermediate and S atom. In this transition state (TS2), the activation energy barrier of the reaction was 0.41 eV and the reaction exotherm was 3.64 eV. The variation of C-S bond length in the direct hydrodesulfurization reaction was shown in Table 7.

Table 7.

C-S bond lengths (Å) in direct hydrodesulfurization reactions of adsorbed thiophene on Ni₂P (001).

Reactant	product	d(C²-S)	d(C⁵-S)
hcp	hcp	1.785	1.785
hcp	TS1	1.813	1.825
hcp	C₄H₅S	1.829	1.851
hcp	TS2	1.836	1.847
hcp	C₄H₆+S	1.840	1.849

Reactions A∼E belonged to the indirect hydrodesulfurization process, in which the hydrogenation process and the broken process of the C-S bond competed with each other. The hcp configuration in the most stable flat mode was selected as the initial adsorption state of thiophene. The thiophene adsorbed at the hcp position and H adsorbed at adjacent vacancies were selected as the initial state of reaction. The product formed by adding hydrogen atoms to different positions of the thiophene was selected as the final state of the reaction. A series of hydrogenations that thiophene was hydrogenated to 2-MHT and 3-MHT at position 2 (same as position 5, collectively referred to as alpha position) and position 3 (same as position 4, collectively referred to as beta position) respectively. Then 2-MHT and 3-MHT were further hydrogenated to THT and then hydrogenated to n-butane and sulfur was studied.

For the reaction C₄H₄S*+H₂→C₄H₆S*, it could be seen from Figure 7 that the main difference between the indirect hydrogenation reactions A-E was the difference of the initial hydrogenation position. Step 1 of the reaction mechanism A-C and step 1 of mechanism D and E were studied.

For the reaction C₄H₄S*+H₂→α-C₄H₅S*+*, firstly, the stable adsorption configurations of C₄H₄S+H and α-C₄H₅S on the surface of Ni₂P (001) were calculated, and used as the initial and final states of the reaction. For the reaction C₄H₄S*+H₂→β-C₄H₅S*+*, the initial reaction model was the same as the step 1 of reaction mechanism A-C. Corresponding simulation calculation and transition state screening were performed on the initial reaction. In the transition state TS3, the H atom on the α-C in the C₄H₄S ring was tilted upward at a certain angle to shift the thiophene ring, and the free H atom gradually approached α-C and formed an α-C-H bond. The activation energy barrier of this reaction was 0.43 eV, and the reaction exotherm was 0.57 eV. In TS4, the thiophene ring did not shift significantly, and the H atom gradually approached β-C in C₄H₄S to form a β-C-H bond. The activation energy barrier of this reaction was 0.89 eV, and the reaction exotherm was 0.53 eV. Figure 9 showed the potential energy surface diagram (PES) of two reactions and the schematic diagram of the related intermediate structure.

Figure 9.

Hydrogenation pathways of thiophene on Ni₂P (001) (Mechanisms A∼C(1),D,E(1)).

According to the data in Figure 9, it could be analyzed from the perspective of kinetics that the activation energy of step 1 of reaction A-C was 0.46 eV lower than that of step 1 of D and E, and H introduced in A-E step 1 was more inclined to combine with the α-C in C₄H₄S, which was consistent with the analysis of adsorption behavior. The C-S bond of the C₄H₄S ring in step 1 of reaction A-C was longer than that in step 1 of reaction D and E, which indicated that the C-S bond in step 1 of reaction A-C was more likely to break, and the reaction was more likely to progress to step 1 of reaction A-C. From the perspective of thermodynamics, the energy of the product in step 1 of A-C was 0.57 eV lower than that of the reactant, and the energy of the product in step 1 of D and E was 0.53 eV lower than the energy of the reactant. Both of these reactions were exothermic, and Lowering the temperature could promote the reaction.

α-C₄H₅S would have two different hydrogenation positions in the next step, resulting in two different reaction pathways. one was the hydroisomerization of α-C₄H₅S at α-C, and the other was hydrogenation of α-C₄H₅S at β-C. The two reactions formed dihydrothiophene (MHT).

For α-C₄H₅S*+H*→α,α-C₄H₆S*+* and α-C₄H₅S*+H*→α,β-C₄H₆S*+*, the dominant adsorption configurations of α-C₄H₅S+H and α,α-C₄H₆S, α-C₄H₅S+H and α,β-C₄H₆S on the Ni₂P (001) surface were optimized. The configurations were used as the initial state and the final state to carry out the transition state search. The activation energy barrier of the reaction α-C₄H₅S*+H*→α,α-C₄H₆S* +* was 1.92 eV, and the heat absorption of the reaction was 0.22 eV. The activation energy barrier of the reaction α-C₄H₅S*+H*→α,β-C₄H₆S*+* was 1.47 eV, and the reaction exotherm was 0.23 eV. Figure 10 was a schematic diagram of a reaction potential energy surface map (PES) and related intermediate structures.

Figure 10.

Hydrogenation pathways of thiophene on Ni₂P (001)(Mechanisms A(2),B,C(2)).

As could be seen from the data, the activation energy of step 2 of reaction mechanism B and C was reduced by 0.45 eV compared with step 2 of reaction mechanism A, which indicated that the reaction was easier to proceed along reaction mechanism B and C. The energy difference before and after the step 2 of reaction A was 0.22 eV, which was an endothermic reaction; the energy difference between step 2 of reaction B and C was −0.23 eV, which was an exothermic reaction. Therefore, the second-step hydrogenation reaction of thiophene on the Ni₂P (001) surface was easier to proceed along step 2 of mechanism B and C, and the product was α,β-C₄H₆S.

For the reaction α,β-C₄H₆S*+H₂*→C₄H₇S* further hydrogenation of α,β-C₄H₆S might generate α,β,α-C₄H₇S or α,β,β-C₄H₇S. The IS, TS, FS and activation energy, etc. of reaction mechanism B, C step 3 were listed in Figure 11. For reactions α,β-C₄H₆S*+H*→α,β,β-C₄H₇S*+* and α,β-C₄H₆S*+H*→α,β,α-C₄H₇S*+*, the dominant adsorption configuration was calculated first, and the reaction energy barrier was shown in Figure 11. In TS7, the thiophene ring had a significant displacement, the hydrogen atoms on β-C were tilted away from the adsorption surface, and the free-state H atoms gradually approached β-C to form β-C-H bonds. After the β-C bond was hydrogenated, it further tilted up. The activation energy barrier of this reaction was 1.28 eV, and the heat absorption of the reaction was 0.88 eV. In TS8, the thiophene ring was deformed and the S-C² bond had obvious elongation. The activation energy barrier of this reaction was 1.14 eV and the reaction exotherm was 0.37 eV.

Figure 11.

Hydrogenation pathways of thiophene on Ni₂P (001) (Mechanisms B(3),C(3)).

It could be seen from Figure 11 that the activation energy of step 3 of reaction mechanism C was lower than that of step 3 of reaction mechanism B by 0.14 eV. Meanwhile, the energy difference before and after the step 3 of reaction C was −0.37 eV, and the reaction was exothermic. The energy difference before and after step 3 of reaction B was 0.88 eV, and the reaction is endothermic. Therefore, the hydrodesulfurization of thiophene on the Ni₂P (001) surface was easier to perform according to step 3 of reaction mechanism C.

For the reaction α,β,α-C₄H₇S*+H*→C₄H₈S*+*, the dominant adsorption configurations of α,β,α-C₄H₇S+H and C₄H₈S on the Ni₂P (001) surface were calculated. In TS9, the β-C slowly tilted upward, the steric hindrance effect of thiophene decreased, and at the same time, the free hydrogen atoms gradually approached the β-C to form the β-C-H bond. The product of this reaction was tetrahydrothiophene (THT). The activation energy barrier of the reaction was 2.38 eV, and the reaction endothermic was 0.12 eV (Figure 12).

Figure 12.

Hydrogenation pathway of thiophene on Ni₂P (001) (Mechanism C(4)).

In the reaction C₄H₈S*+H₂*→C₄H₁₀*+*, the hydrogenation reaction of tetrahydrothiophene caused the C-S bond to be broken. The initial state, final state, transition state, and reactive energy barrier of steps 5 and 6 of reaction mechanism were listed in Figure 13 respectively. For C₄H₈S*+H*→C₄H₉*+*, the structure optimization calculation of the dominant adsorption configuration was carried out, which was regarded as the initial state and the final state. In TS10, C₄H₈S had a slight displacement, the ring plane tilted upwards, and the free-state H atoms gradually approached α-C and formed the α-C-H bond. The activation energy barrier of this reaction was 2.06 eV and the reaction exotherm was 0.39 eV.

Figure 13.

Hydrogenation pathway of thiophene on Ni₂P (001) (Mechanism C(5),C(6)).

For C₄H₉*+H*→C₄H₁₀*+*, the structure optimization calculation of the dominant adsorption configuration was carried out, which was regarded as the initial state and the final state. In TS11, the α-C-S bond was broken to form C₄H₁₀ and S atoms, and the S atoms are adsorbed at the triple vacancy Ni atom. The activation energy barrier of this reaction was 1.75 eV, and the reaction exotherm was 0.5 eV.

In summary, the potential energy surface diagrams of the direct hydrodesulfurization and indirect hydrodesulfurization of thiophene on the Ni₂P (001) surface were made as shown in Figure 14. As shown in the figure, the energy barrier of direct hydrodesulfurization of thiophene in the first step was significantly smaller than the energy barrier of indirect hydrodesulfurization of thiophene in the first step (E_a=0.07 eV<E_a=0.43 eV), therefore, the direct hydrodesulfurization of thiophene was more dominant on the surface of Ni₂P (001).

Figure 14.

Sketch for potential relative energy of direct hydrodesulfurization and indirect hydrodesulfurization on Ni₂P(001) surface.

Conclusions

When thiophene molecules were adsorbed in the triple vacancy of the Ni₂P (001) surface by a flat mode, the adsorption energy was the largest and the adsorption configuration was the most stable.

The direct hydrodesulfurization of thiophene had more advantages than the indirect hydrodesulfurization of thiophene. The direct hydrodesulfurization reaction was an endothermic reaction, and the increasing temperature could promote the reaction. The indirect hydrodesulfurization was an exothermic reaction, and cooling down could promote the reaction. In the indirect hydrodesulfurization reaction, the step 4 of reaction mechanism C (TS4) required the highest activation energy, which was the rate-limiting step of the entire reaction. Direct hydrodesulfurization reaction step was C₄H₄S+H₂→C₄H₆.

Footnotes

Acknowledgements

The authors would like to acknowledge the financial support provided by the National Science and Technology Major Project of the Ministry of Science and Technology of China (2016ZX05012-002–005), National Natural Science Foundation of China (No. 51874333) and Natural Science Foundation of Shandong Province, China (No. ZR2017MEE030).

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Riyi Lin

References

Angelici

(2006) Hydrodesulfurization & Hydrodenitrogenation, Encyclopedia of Inorganic Chemistry. Hoboken, NJ: John Wiley & Sons, Ltd.

Ankudinov

Ravel

Rehr

, et al. (1998) Real-space multiple-scattering calculation and interpretation of x-ray-absorption near-edge structure. Physical Review B 58(12): 7565–7576.

Bataille

Lemberton

Michaud

, et al. (2000) Alkyldibenzothiophenes hydrodesulfurization promoter effect, reactivity, and reaction mechanism. Journal of Catalysis 191(2): 409–422.

Chu

Feng

, et al. (2013) Synthesis of Ni/Mo/N catalyst and its application in benzene hydrogenation in the presence of thiophene. Chinese Journal of Catalysis 34(1): 159–166.

Clark

Wang

Oyama

(2002) Characterization of silica-supported molybdenum and tungsten phosphide hydroprocessing catalysts by ³¹P nuclear magnetic resonance spectroscopy. Journal of Catalysis 207(2): 256–265.

Fuks

Vingurt

Landau

, et al. (2010) Density functional theory study of sulfur adsorption at the (001) surface of metal-rich nickel phosphides: Effect of the Ni/P ratio. The Journal of Physical Chemistry C 114(31): 13313–13321.

Hua

Jian

Feng

, et al. (2013) Advances in catalytic hydrodesulfurization mechanism of nickel phosphide catalysts. Acta Petrolei Sinica 29(4): 717–725.

Jayamurthy

Vasudevan

(1994) Temperature-programed desorption and surface-reaction of thipphene over Co-Mo/Gamma-Al2O3 hydrodesulfurization catalysts. The Journal of Physical Chemistry 98(27): 6777–6784.

Guan

Wei

, et al. (2015) A potential regularity for enhancing the hydrogenation properties of Ni₂P. The Journal of Physical Chemistry C 119(5): 2557–2565.

10.

McCrea

Logan

Tarbuck

, et al. (1997) Hydrodesulfurization over alumina-supported molybdenum carbide and nitride catalysts: Effect of Mo loading and phase. Journal of Catalysis 171(1): 255–267.

11.

Nag

Sapre

Broderick

, et al. (1979) Hydrodesulfurization of polycyclic aromatics catalyzed by sulfided CoO-MoO₃/γ-Al₂O₃: The relative reactivities. Journal of Catalysis 57(3): 509–512.

12.

Nelson

Sun

Junaid

(2006) On the structure and composition of the phosphosulfide overlayer on Ni₂P at hydrotreating conditions. Journal of Catalysis 241(1): 180–188.

13.

Shi

Xia

(2013) Theoretical studies on reaction mechanism of hydrodesulfurization of thiophene catalyzed by Au(111) plane. Chemical Journal of Chinese Universities 34(10): 2353–2362.

14.

Okamoto

Ishihara

Kawano

, et al. (2003) Preparation of Co-Mo/Al₂O₃ model sulfide catalysts for hydrodesulfurization and their application to the study of the effects of catalyst preparation. Journal of Catalysis 217(1): 12–22.

15.

Oyama

(2003) Novel catalysts for advanced hydroprocessing: Transition metal phosphides. Journal of Catalysis 216(1–2): 343–352.

16.

Oyama

Gott

Zhao

, et al. (2009) Transition metal phosphide hydroprocessing catalysts: A review. Catalysis Today 143(1–2): 94–107.

17.

Oyama

Wang

Lee

, et al. (2004) Active phase of Ni₂P/SiO₂ in hydroprocessing reactions. Journal of Catalysis 221(2): 263–273.

18.

Phillips

Sawhill

Self

, et al. (2002) Synthesis, characterization, and hydrodesulfurization properties of silica-supported molybdenum phosphide catalysts. Journal of Catalysis 207(2): 266–273.

19.

Rundqvist

(1961) X-ray investigations of Mn₃P, Mn₂P and Ni₂P. Uppsala University, Sweden.

20.

Song

Dai

Guo

, et al. (2012) Preparation of composite TiO₂–Al₂O₃ supported nickel phosphide hydrotreating catalysts and catalytic activity for hydrodesulfurization of dibenzothiophene. Fuel Processing Technology 96: 228–236.

21.

Song

Wang

, et al. (2014) Effect of titanium content on dibenzothiophene HDS performance over Ni₂P/Ti-MCM-41 catalyst. Journal of Catalysis 311: 257–265.

22.

Stephanie

Diana

Mark

(2003) Thiophene hydrodesulfurization over supported nickel phosphide catalysts. Journal of Catalysis 215: 208.

23.

Sun

Prins

(2008) Hydrodesulfurization of 4,6-dimethyldibenzothiophene over noble metals supported on mesoporous zeolites. Angewandte Chemie (International ed. in English) 47(44): 8478–8481.

24.

Zhu

Zhao

, et al. (2006) Research progress of hydrodesulfurization catalyst and reaction mechanism. Journal of Molecular Catalysis 20(4): 372–383.

Step	Mechanism A	Mechanism B	Mechanism C
1	C₄H₄S+H→α-C₄H₅S+	C₄H₄S+H→α-C₄H₅S+	C₄H₄S+H→α-C₄H₅S+
2	α-C₄H₅S+H→α,α-C₄H₆S+	α-C₄H₅S+H’→α,β-C₄H₆S+	α-C₄H₅S+H→α,β-C₄H₆S+
3	α,α-C₄H₆S+H→α,α,β-C₄H₇S+	α,β-C₄H₆S+H→α,β,β-C₄H₇S+	α,β-C₄H₆S+H→α,β,α-C₄H₇S+
4	α,α,β-C₄H₇S+H→C₄H₈S+	α,β,β-C₄H₇S+H→C₄H₈S+	α,β,α-C₄H₇S+H→C₄H₈S+
5	C₄H₈S+H→C₄H₉S+	C₄H₈S+H→C₄H₉S+	C₄H₈S+H→C₄H₉S+
6	C₄H₉S+H→C₄H₁₀+S	C₄H₉S+H→C₄H₁₀+S	C₄H₉S+H→C₄H₁₀+S

Step	Mechanism D	Mechanism E	Mechanism F
1	C₄H₄S+H→β-C₄H₅S+	C₄H₄S+H→β-C₄H₅S+	C₄H₄S+H→C₄H₅S+
2	β-C₄H₅S+H→β,α-C₄H₆S+	β-C₄H₅S+H→β,α-C₄H₆S+	C₄H₅S+H→C₄H₆+S
3	β,α-C₄H₆S+H→β,α,β-C₄H₇S+	β,α-C₄H₆S+H→β,α,α-C₄H₇S+	/
4	β,α,β-C₄H₇S+H→C₄H₈S+	β,α,α-C₄H₇S+H→C₄H₈S+	/
5	C₄H₈S+H→C₄H₉S+	C₄H₈S+H→C₄H₉S+	/
6	C₄H₉S+H→C₄H₁₀+S	C₄H₉S+H→C₄H₁₀+S	/

Molecular simulation on mechanism of thiophene hydrodesulfurization on surface of Ni 2 P

Abstract

Keywords

Introduction

Calculation method

Simulation content

Construction and optimization of model molecules

Selection of Ni2P model adsorption surface

Results and discussion

Adsorption of thiophene on Ni2P (001) surface

HDS reaction path of thiophene on Ni2P surface

Conclusions

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

ORCID iD

References

Selection of Ni₂P model adsorption surface

Adsorption of thiophene on Ni₂P (001) surface

HDS reaction path of thiophene on Ni₂P surface