Identifying the active nickel sites and hydrogen species in the phosphosulfide overlayer on the Ni₂P hydrotreating catalyst: a DFT-D3 study

Abstract

Density functional theory (DFT) calculations and an atomistic thermodynamic approach were employed to identify the active nickel sites and determine the chemical nature of the hydrogen species in the nickel phosphosulfide overlayer formed on the (0001) and (101̄0) surfaces of the Ni₂P hydrotreating catalyst. Our results showed that Ni(1) sites are less resistant to sulfur than Ni(2) sites, in agreement with experimental observations. Under HDS conditions, the fully sulfided (0001) surface of Ni₂P was found to be structurally similar to the (111) surface of the low-activity sulfide Ni₃S₂, where closely coordinated Ni atoms are prone to sulfidation. Non-hydrogenated surfaces were always the most stable, suggesting that NiH and SH groups are unlikely to form on these surfaces. In contrast, the most stable (101̄0) surface was identified as a hydrogenated surface with coordinatively unsaturated Ni(2) atoms surrounded by SH groups. Depending on the electronegativity of their ligands (S or Ni), the hydrogen species exhibit either protonic (H^δ+) or hydridic (H^δ−) character. Our results suggest that protons within SH groups are the most likely reactive hydrogen species on the nickel phosphosulfide overlayer under reaction conditions, providing atomic-level insight into the origins of HYD activity in Ni₂P catalysts.

---

1. Introduction

Nickel phosphide (Ni₂P) has emerged as a promising candidate for a new generation of industrial hydrodesulfurization (HDS) catalysts.^1–7 According to the literature, the active phase of the Ni₂P catalyst is a nickel phosphosulfide overlayer that forms during the HDS reaction.^1,8–12 This overlayer can also be generated by pretreatment with sulfidation agents, such as H₂S or CS₂ in the presence of H₂, or by reducing sulfur-containing precursors, such as nickel hexathiodiphosphate (Ni₂P₂S₆).^13–16

The active sites in HDS catalysts are primarily linked to coordinatively unsaturated metal sites (CUSs), sulfhydryl (SH) groups, and metal hydrides (MH).^17,18 CUSs are S²⁻ anion vacancies located at the edges or corners of catalyst nanoparticles. SH groups, on the other hand, result from the dissociative adsorption of H₂ onto surface sulfur atoms. These groups behave as Brønsted acid sites.^19,20 MH species are formed when H₂ dissociatively adsorbs onto CUSs. Consequently, it is essential to understand the structural changes of Ni₂P under reaction conditions and the specific catalytic functions of the active sites.

In an early study, Oyama et al. investigated the catalytic activity of Ni₂P supported on low- and high-surface-area siliceous supports (SiO₂-L, SiO₂-H, and MCM-41) for the HDS of 4,6-dimethyldibenzothiophene (4,6-DMDBT).²¹ Extended X-ray absorption fine structure (EXAFS) analysis confirmed the presence of two types of nickel sites: Ni(1), with tetrahedral coordination to P, and Ni(2), with square pyramidal coordination to P. The analysis further revealed both Ni(1) and Ni(2) sites are present in large crystallites, whereas Ni(2) sites are more numerous on the small crystallites. Catalytic selectivity measurements showed that small crystallites exhibit higher selectivity toward HYD products, while large crystallites have higher selectivity toward DDS products. Based on product selectivity, it was concluded that the Ni(1) sites are primarily responsible for the DDS pathway, while the Ni(2) sites are highly active for the HYD pathway.

In a later study, Oyama et al. investigated the sulfur-resistant properties of the Ni₂P catalyst in deep hydrodesulfurization.³ EXAFS analysis of spent catalysts revealed that Ni(2) sites are bound to sulfur with a lower Ni–S coordination as the particle size decreased. Based on EXAFS results, the authors concluded that the active Ni(2) sites on highly dispersed Ni₂P catalysts are much more tolerant of sulfur than the tetrahedral Ni(1) sites. Furthermore, the higher selectivity for HYD products correlates with an increased number of surface Ni(2) sites, which predominate in smaller Ni₂P particles.

To elucidate the catalytic mechanism of the HDS reaction on the Ni₂P catalyst, it is essential to understand the state of adsorbed hydrogen on the nickel phosphosulfide overlayer. To the best of our knowledge, there are no experimental studies on the chemical nature of hydrogen species adsorbed on Ni₂P as an HDS catalyst. Only Nelson et al., using density functional theory (DFT) calculations, have studied the presence of SH groups on the (0001) surface of the Ni₂P catalyst.²² The authors calculated the free energy changes associated with H₂S adsorption and its dissociation into –SH and sulfur on the (0001) Ni₃P₂-terminated surface as a function of temperature (100–1000 K) under H₂ and H₂S partial pressures of 5.0 and 0.1 MPa, respectively. Their results indicated that the relative stability of the adsorbed sulfur species follows the order H₂S < SH < S. These results suggest that hydrogen bound as an SH group might not be thermodynamically stable under HDS reaction conditions.

On the other hand, electrocatalysis experiments have shown that incorporating sulfur improves the performance of Ni₂P in the hydrogen evolution reaction (HER).^23–25 Although the HER and HDS occur under different conditions, both processes involve hydrogen atoms adsorbed on the surface.²⁶ It was found that sulfur incorporation accelerates the adsorption rate of reaction intermediates on the catalyst surface, thereby improving HER kinetics. Therefore, sulfur incorporated into the catalyst may act as an activity promoter or selectivity modifier.²⁷

In this work, we combined DFT calculations with an atomistic thermodynamic analysis to investigate the stability of hydrogen species on the phosphosulfide overlayer formed on the (0001) and (101̄0) surfaces of the Ni₂P catalyst under HDS conditions. Our goal is to elucidate the nature of the active nickel sites and the chemical characteristics of the hydrogen species, which are key to understanding the hydrogenation functionality of the phosphosulfide overlayer.

2. Computational details

2.1. Computational methods

DFT calculations were performed using the Vienna ab initio simulation package (VASP).^28,29 The Perdew–Burke–Ernzerhof exchange–correlation functional with a generalized gradient approximation was used.³⁰ The projector augmented wave (PAW) method was used to describe the electron–ion interactions.³¹ Standard PAW potentials with a valence electron distribution of Ni 3d⁸4s², P 3s²3p³, S 3s²3p⁴, and H 1s¹ were used. A plane-wave kinetic energy cut-off of 400 eV was employed. The energy and atomic force convergence criteria were set to 10⁻⁵ eV and 0.01 eV Å⁻¹, respectively. The Gaussian smearing method with a width of 0.05 eV was used to describe the partial occupancies of the electronic states. The DFT-D3 correction developed by Grimme et al. was used to capture dispersion interactions, which are important to the energetics of catalytic transformations.^32–35 Gamma-centered k-point meshes of 5 × 5 × 9 and 7 × 7 × 7 were used to sample the Brillouin zone for the optimization of bulk Ni₂P and Ni₃S₂, respectively. The calculated lattice parameters for bulk Ni₂P are a = b = 5.836 Å and c = 3.307 Å, while for bulk Ni₃S₂ (rhombohedral representation), a = b = c = 4.034 Å. These values are in good agreement with the reported experimental values for Ni₂P (a = b = 5.859 Å, c = 3.382 Å)³⁶ and Ni₃S₂ (a = b = c = 4.082 Å).³⁷ All surface calculations were performed using symmetric slabs. The (0001) Ni₃P₂ surface termination was modeled using a nine-atomic-layer (2 × 2) slab with a 3 × 3 × 1 gamma-centered k-point mesh (Fig. S1a). The (101̄0) AB Ni₂P surface termination was modeled with ten-atomic-layer (2 × 2) slabs with a 3 × 5 × 1 Monkhorst–Pack k-point mesh (Fig. S1b). The (111) surface of Ni₃S₂ was modeled with eight-atomic-layer (2 × 2) slab with a 3 × 3 × 1 gamma-centered k-point mesh (Fig. S1c). A vacuum space of 20 Å was set to avoid interactions between periodic slabs, and all atoms were fully relaxed during the calculations. Bader charge³⁸ and electron density difference (EDD) analyses were performed to assess the charge distribution and the chemical nature of the hydrogen species adsorbed on the surface. The EDD was calculated from the following equation:

Δρ = ρ(surface + adsorbate) − ρ(surface) − ρ(adsorbate) (1)

where ρ(surface + adsorbate) is the electron density of the whole system, ρ(surface) and ρ(adsorbate) are the electron densities of the surface, and the adsorbate calculated at the adsorption positions, respectively. The VESTA program was used to visualize the electron density difference isosurfaces.³⁹

2.2. Surface models

To construct slab models of the possible structures of the nickel phosphosulfide overlayer on the Ni₂P catalyst, we selected the (0001) and (101̄0) surfaces. Papawassiliou et al.,⁴⁰ using ³¹P solid-state nuclear magnetic resonance (ssNMR), high-resolution transmission electron microscopy (HRTEM), X-ray diffraction (XRD), and DFT calculations, determined that Ni₂P nanoparticles synthesized under excess phosphorus predominantly expose these surfaces. For the (0001) surface, the slab model was constructed from a (2 × 2) supercell with a Ni₃P₂ termination. This surface termination has been reported as the most thermodynamically stable.^41,42 The surface is described by a hexagonal supercell (a = b = 11.672 Å) and contains eight P atoms and twelve Ni atoms in the outermost atomic layer, Ni₁₂P₈ (Fig. S1a). These Ni atoms are referred to as Ni(1). On the surface, the Ni(1) atoms form four three-fold adsorption sites (Ni₃ sites), which are labeled as T sites, as illustrated in Fig. S1. For the (101̄0) surface, the slab model surface was constructed from a (2 × 2) supercell with the AB-type Ni₂P termination, which has also been reported as the most thermodynamically stable surface termination.^41,42 The surface is described by a rectangular supercell (a = 11.672 Å and b = 6.614 Å) and contains four P atoms and eight Ni atoms in the outermost atomic layer, Ni₈P₄ (Fig. S1b). This surface contains two distinct types of Ni atoms: four Ni(1) and Ni(2). In Fig. S1b the Ni(2) nickel atoms are represented by light blue spheres to distinguish them from the Ni(1) nickel atoms. On the surface, the Ni atoms adopt a zigzag arrangement along the y-axis. (Fig. S1). This zigzag configuration exhibits two types of Ni₃ sites: four T1 sites and four T2 sites. Each T1 site consists of two Ni(1) atoms and one Ni(2) atom, while the T2 site consists of two Ni(2) atoms and one Ni(1) atom. It has been shown that the S atoms are preferentially located at the T1 sites.⁴² At 50%S, the S atoms are preferentially adsorbed at the T1 sites in an alternating configuration. The adsorption energy difference between T1 and T2 sites was around 1 eV. For the (111) surface of the Ni₃S₂ catalyst (rhombohedral representation), the slab model was constructed from a (2 × 2) supercell with the Ni₃S₂ termination. The surface is described by a hexagonal supercell (a = b = 11.360 Å), with eight S atoms and twelve Ni atoms in the outermost atomic layer, Ni₁₂S₈ (Fig. S1c). The surface Ni atoms also form four Ni₃ sites, which are also labeled as T sites.

The slab models of the nickel phosphosulfide overlayer on the (101̄0) and (0001) surfaces were constructed following the surface modification procedure previously reported.⁴² These slab models take into account the processes of sulfur adsorption and the replacement of surface phosphorus by sulfur at the outermost atomic layer. For the (0001)-Ni₁₂P₈ surface, the replacement of P by S was evaluated at 0%, 50%, and 100% of S, where 0%S corresponds to no replacement, 50%S corresponds to replacement with four S atoms, and 100%S corresponds to replacement with eight S atoms. Additionally, sulfur adsorption was evaluated at coverages of 75%S and 100%S, which correspond to the adsorption of three and four S atoms, respectively. Surface terminations are labeled according to the elemental composition of the outermost layer and the number of S atoms adsorbed. For the (101̄0)-Ni₈P₄ surface, the replacement of P by S was evaluated at 0%, 50%, and 100% of S, where 0%S corresponds to no replacement, 50%S corresponds to replacement with two S atoms, and 100%S corresponds to replacement with four S atoms. This surface has two types of adsorption sites for sulfur: four T1 sites and four T2 sites (Fig. S1b). For the slab models of the hydrogenated surfaces, we evaluated the dissociative adsorption of one or two H₂ molecules, corresponding to hydrogen coverages of 50% and 100%, respectively.

2.3. Thermodynamic approach

Since HDS operating conditions are sulforeductive conditions in which H₂ is more abundant than H₂S in the gas phase, PH₃ could be produced when P atoms are replaced by S atoms in the presence of H₂. Therefore, the possibility of complete sulfidation of bulk Ni₂P to bulk Ni₃S₂ under relevant conditions was evaluated by the following reaction:

Ni₂P(bulk) + 4/3H₂S + 3/2H₂ → 2/3Ni₃S₂(bulk) + PH₃ + 4/3H₂ (2)

and the corresponding Gibbs free energy, , is then given bywhere and are the DFT total energies of bulk Ni₃S₂ and Ni₂P phases, respectively. The chemical potential of phosphorus, μ_P(T, p), is determined by the chemical equilibrium with the gas-phase mixture of PH₃ and H₂,

μ_P(T, p) = μ_PH3(T, p) − (3/2)μ_H2(T, p) (4)

μ_PH3(T, p) and μ_H2(T, p) are the chemical potentials of PH₃ and H₂, respectively, evaluated assuming ideal gas behavior. The calculations of chemical potentials are more detailed in the SI.

On the other hand, the chemical potential of sulfur, μ_S(T, p), is determined by the chemical equilibrium with the gas-phase mixture of H₂ and H₂S,

μ_S(T, p) = μ_H2S(T, p) − μ_H2(T, p) (5)

where μ_H2S(T, p) and μ_H2(T, p) are the chemical potentials of H₂S and H₂, respectively, also evaluated assuming ideal gas behavior. Under typical HDS conditions, the ratio ranges from 0.01 to 0.05, and the temperature varies between 573 and 700 K. In this study, we selected 573 K as the working temperature because several experimental studies using Ni₂P catalyst have performed HDS measurements at this temperature.^1,2,10,12,13

The possibility of replacing the surface P atoms with S atoms, as well as the adsorption of S atoms under relevant conditions, can be described as

Ni_xP_y + (z + a)H₂S + (3/2z)H₂ → Ni_xP_(y−z)S_z_aS + zPH₃ + (z + a)H₂ (6)

The Gibbs free energy associated with these processes (adsorption and replacement) can be calculated as follows:where is the DFT total energy of a slab model (either the (0001). or (101̄0) surface), in which z surface phosphorus atoms have been replaced by sulfur atoms, in addition to the adsorption of a sulfur atoms on the Ni₃ sites. is the DFT total energy of the pristine surface, which is taken as a reference.

On the other hand, the Gibbs free energy change associated only with S atom adsorption (ΔG^S_ads), is given by

where and are the DFT total energy of sulfided and clean surfaces, respectively. a is the number of S atoms adsorbed.

Finally, the dissociative adsorption of H₂ as NiH or SH groups on the nickel phosphosulfide overlayer was evaluated by the following reaction:

Ni_xP_(y−z)S_z_aS + nH₂ → Ni_xP_(y−z)S_z_aS_2nH (9)

and the corresponding Gibbs free energy, , is given bywhere and are the DFT total energy of hydrogenated and non-hydrogenated surfaces, respectively. n is the number of dissociatively adsorbed H₂ molecules, and μ_H(T, p) is the hydrogen chemical potential, where μ_H(T, p) = (½)μ_H2(T, p).

3. Results and discussion

3.1. Bulk Ni₂P stability under the H₂S/H₂ gas mixture

Group 8 transition metal catalysts (nickel, palladium, and platinum) are highly susceptible to sulfur poisoning, with nickel being the most susceptible to sulfide formation.^43–45 However, nickel phosphide-based catalysts have been shown to be highly resistant to sulfur poisoning.^3,7,12,46 Given that Ni₃S₂ is the most stable stoichiometry of nickel sulfide under HDS reaction conditions,^47–49 we investigated the relative stability of bulk Ni₂P over bulk Ni₃S₂ in the presence of the H₂S/H₂ gas mixture. Bulk Ni₂P belongs to the hexagonal crystal system with the P6̄2m space group (Fig. 1a). It is characterized by the presence of two distinct types of Ni atoms, denoted as Ni(1) and Ni(2).^21,50 Ni(1) atoms have a quasi-tetrahedral coordination geometry with the four nearest-neighbor P atoms, while Ni(2) atoms exhibit a square pyramidal coordination geometry with the five nearest-neighbor P atoms. In Fig. 1a the Ni(2) nickel atoms are represented by light blue spheres to distinguish them from the Ni(1) nickel atoms. Bulk Ni₃S₂ belongs to the trigonal crystal system with the R32 space group. It can be described in both rhombohedral (Fig. 1b) and hexagonal (Fig. 1c) lattice representations. This bulk only has one type of Ni atom, which has a quasi-tetrahedral coordination geometry with the four nearest-neighbor S atoms.

Fig. 1 Optimized crystal structure of (a) bulk Ni₂P, (b) bulk Ni₃S₂ in the rhombohedral representation, and (c) bulk Ni₃S₂ in the hexagonal representation. Blue, light blue, pink, and yellow spheres are the Ni(1) nickel, Ni(2) nickel, phosphorus, and sulfur atoms, respectively. (d) Phase diagram for the complete sulfidation of bulk Ni₂P to bulk Ni₃S₂ as a function of μ_P(T, p) and μ_S(T, p). The dependence of the chemical potentials with the gas phase is translated into pressure scales for T = 573 K on the upper x-axis for the H₂S/H₂ mixture and on the right y-axis for PH₃ at p_H2 = 1 and 10 bar. The red dotted line corresponds to typical HDS partial pressures ratio, .

Fig. 1d shows the phase diagram of the free energy change for complete sulfidation of bulk Ni₂P to bulk Ni₃S₂ as a function of μ_P(T, p) and μ_S(T, p). The regions of stability of the bulk phases are represented by the distinct colored areas. The red dotted line in Fig. 1d represents a typical HDS ratio. It is observed that the complete sulfidation of Ni₂P to Ni₃S₂ may occur when μ_P(T, p) is lower than −6.77 eV. At 573 K, this condition is met when p_PH3 is lower than 3 × 10⁻¹⁵ bar with p_H2 = 1 bar or when p_PH3 is lower than 10⁻¹³ bar with p_H2 = 10 bar. PH₃ is a very strong phosphiding agent, and consequently, competition could exist between the sulfidation and phosphidation processes. Our findings suggest that bulk Ni₂P will not undergo complete sulfidation to bulk Ni₃S₂ as long as there is a trace amount of PH₃ in the gas phase. This result correlates with the sulfur resistance observed in Ni₂P catalysts under HDS conditions.^3,7,12,46

3.2. Phosphosulfide overlayer on the (0001) surface

We revisited the adsorption and replacement processes of sulfur atoms on the (0001) surface by calculating the free energy change as a function of the phosphorus and sulfur chemical potentials, μ_P(T, p) and μ_S(T, p). The analysis was carried out considering the thermodynamic limit in which bulk Ni₂P transforms into bulk Ni₃S₂. Fig. 2a presents the corresponding free energy phase diagram. The stability regions of the different surface compositions are indicated by the different shades of green. Fig. 2b–d show the top view of the corresponding surface structures. Surface terminations are labeled according to the elemental composition of the outermost layer and the number of S atoms adsorbed. In Fig. 2a, the red dotted line represents a typical HDS ratio, while the blue dotted line indicates the phase boundary between the bulk Ni₂P and Ni₃S₂. The (0001)-Ni₁₂P₈_4S, (0001)-Ni₁₂P₄S₄_4S, and (0001)-Ni₁₂S₈_4S surfaces are the most stable at a ratio of 0.01. All of these surfaces have all Ni₃ sites covered by S atoms. The findings of this study are consistent with those previously reported.⁴²

Fig. 2 (a) Phase diagram of the most stable (0001)-Ni_xP_y−zS_z_S_z’ surface composition as a function of μ_P(T, p) and μ_S(T, p). The dependence of the chemical potentials with the gas phase is translated into pressure scales for T = 573 K on the upper x-axis for the H₂S/H₂ mixture and on the right y-axis for PH₃ at p_H2 = 1 and 10 bar. The red dotted line corresponds to the typical HDS partial pressure ratio, . The blue dotted line indicates the phase boundary between the bulk Ni₂P and Ni₃S₂. Top view of the optimized structures of the (2 × 2) supercell: (b) (0001)-Ni₁₂P₈_4S, (c) (0001)-Ni₁₂P₄S₄_4S, and (d) (0001)-Ni₁₂S₈_4S surfaces. Blue, pink, and yellow spheres are the nickel, phosphorus, and sulfur atoms, respectively.

Note that the (0001)-Ni₁₂S₈_4S surface is the most stable at very low PH₃ partial pressures. In the context of HDS experiments, the presence of PH₃ in the gas mixture has not been detected, suggesting that only a trace amount of PH₃ might be produced by the replacement of surface P atoms by S atoms. In contrast to the bulk phase, a very low p_PH3 enables the gradual replacement of all surface P atoms with S atoms, resulting in the formation of a phosphosulfide overlayer on the (0001)-Ni₁₂P₈ surface termination. Therefore, our results suggest that the (0001)-Ni₁₂S₈_4S surface could be the most stable surface termination under HDS reaction conditions. In this regard, Oyama et al. used EXAFS technique to study the structural properties of supported Ni₂P catalysts after the HDS reaction of 4,6-dimethyldibenzothiophene.³ According to EXAFS line shape analysis, catalysts with shorter Ni–Ni bond lengths are more likely to form a phosphosulfide phase. It was found that the tetrahedral Ni(1) sites are less resistant to sulfur than the square-pyramidal Ni(2) sites. On the (0001) surface, the Ni(1) atoms are closely coordinated with another two Ni(1) atoms and, therefore, may be prone to sulfidation. For the Ni₂P catalyst, the active sites are believed to be the coordinatively unsaturated surface nickel atoms belonging to the phosphosulfide phase. Our results suggest that the (0001)-Ni₁₂S₈_4S surface could not be active under HDS conditions because the Ni₃ sites, which consist of Ni(1)-type atoms, would be blocked by sulfur atoms.

It is noteworthy that the surface structure of the (0001)-Ni₁₂S₈_4S surface of Ni₂P closely resembles that of the (111)-Ni₁₂S₈_4S surface of Ni₃S₂ (Fig. 3a and b). Both surface terminations share the same composition and feature Ni₃ sites, where each Ni atom is bonded to two surface S atoms. The Ni–Ni bond distances are 2.646 Å and 2.558 Å, respectively. On the (0001)-Ni₁₂S₈_4S surface, the Ni₃ sites are capped above by an S atom and below a P atom. On the (111)-Ni₁₂S₈_4S surface, the Ni₃ sites are capped above and below by an S atom. This results in an almost regular trigonal bipyramid. It is well known that Ni₃S₂ catalysts exhibit very low activity in various hydrotreating reactions.^47,48,51 In this regard, Aray et al. used DFT calculations to investigate the surface stability of the low-index Miller surfaces and the equilibrium morphology of the Ni₃S₂ catalyst.⁵² It was found that only the (111) and (111̄) type surface faces are exposed under HDS conditions. The S-covered surfaces were the most stable in the whole range of the ratio. These findings suggested that the surfaces of the Ni₃S₂ catalyst may not have coordinatively unsaturated metal sites available for the reaction.

Fig. 3 Top and side views of the optimized structures of the (2 × 2) supercell: (a) (0001)-Ni₁₂S₈_4S and (b) (111)-Ni₁₂S₈_4S surfaces. Blue, pink, and yellow spheres are the nickel, phosphorus, and sulfur atoms, respectively. (c) and (d) ΔG^S_ads as a function of μ_S(T, p), respectively. The dependence of μ_S(T, p) with the gas phase H₂S/H₂ mixture is translated into pressure scales for T = 573 K. The red dotted line corresponds to the typical HDS partial pressure ratio, . The clean surfaces are taken as a reference.

For comparison, we calculated the Gibbs free energy change associated with S atoms adsorption (ΔG^S_ads) on the (0001)-Ni₁₂S₈ and the (111)-Ni₁₂S₈ surfaces as a function of μ_S(T, p) (Fig. 3c and d). The clean surface was taken as a reference. Additionally, ΔG^S_ads was calculated for the (0001)-Ni₁₂P₈ and (0001)-Ni₁₂P₄S₄ surfaces (Fig. S2). It is observed that the fully covered surfaces (with four S atoms adsorbed) are the most stable across the entire range of evaluated μ_S(T, p). For example, at T = 573 K and , the energy required to create a sulfur vacancy on the fully covered surface shows the following trend: (111)-Ni₁₂S₈_4S (1.06 eV) ≈ (0001)-Ni₁₂P₈_4S (1.05 eV) > (0001)-Ni₁₂P₄S₄_4S (0.63 eV) > (0001)-Ni₁₂S₈_4S (0.46 eV). For the (0001) surface, vacancy formation becomes energetically more favorable as surface P atoms are replaced by S atoms; this is due to the strain effect exerted by the S atoms.⁴² Within our thermodynamic framework, surfaces with coordinatively unsaturated sites consisting of spatially isolated Ni₃ ensembles may not be stable under HDS conditions.

3.3. Dissociative adsorption of H₂ onto the phosphosulfide overlayer on the (0001) surface

Dissociative adsorption of H₂ onto the nickel phosphosulfide overlayer can lead to the formation of SH or NiH groups. Fig. 4a–d show the optimized structures of the hydrogenated surfaces resulting from the dissociative adsorption of one and two H₂ molecules on the (0001)-Ni₁₂S₈_4S and (111)-Ni₁₂S₈_4S surfaces. Fig. S3 shows the configurations evaluated for the dissociative adsorption of one H₂ molecule on the (0001)-Ni₁₂S₈_4S surface. When the SH group forms, the sulfur atom binds to two of the three Ni(1) atoms at the Ni₃ site. This leaves the third Ni(1) atom coordinatively unsaturated. Fig. 4e and f display as a function of μ_H(T, p). The non-hydrogenated surface was taken as a reference. We also calculated for the (0001)-Ni₁₂P₈_4S and (0001)-Ni₁₂P₄S₄_4S surfaces (Fig. S4). The configurations evaluated for the dissociative adsorption of one H₂ molecule on the (0001)-Ni₁₂P₈_4S and (0001)-Ni₁₂P₄S₄_4S surfaces are shown in Fig. S5 and S6, respectively. It is observed that the weakest hydrogen adsorption occurs on the (111)-Ni₁₂S₈_4S surface. For example, at T = 573 K and p_H2 = 10 bar, the energy required for the dissociative adsorption of one H₂ molecule shows the following trend: (111)-Ni₁₂S₈_2S_2SH (0.97 eV) > (0001)-Ni₁₂P₈_2S_2SH (0.68 eV) > (0001)-Ni₁₂P₄S₄_2S_2SH (0.62 eV) > (0001)-Ni₁₂S₈_2S_2SH (0.39 eV). For the dissociative adsorption of two H₂ molecules on fully sulfur-covered surfaces, the corresponding energy trend is: (111)-Ni₁₂S₈_4SH (1.94 eV) > (0001)-Ni₁₂P₈_4SH (1.44 eV) > (0001)-Ni₁₂P₄S₄_4SH (1.34 eV) > (0001)-Ni₁₂S₈_4SH (1.05 eV). For the (0001) surface, hydrogen dissociative adsorption surface becomes energetically favorable as surface P atoms are replaced by S atoms. Nevertheless, non-hydrogenated surfaces are always the most stable, which suggests that SH groups are unlikely to form on these surfaces.

Fig. 4 Top view of the optimized structures of the (2 × 2) supercell of hydrogenated surfaces: (a) (0001)-Ni₁₂S₈_2S_2SH, (b) (0001)-Ni₁₂S₈_4SH, (c) (111)-Ni₁₂S₈_2S_2SH, and (d) (111)-Ni₁₂S₈_4SH. Blue, yellow, and white spheres are the nickel, sulfur, and hydrogen atoms, respectively. (e) and (f) as a function of μ_H(T, p), respectively. The dependence of μ_H(T, p) with the gas phase H₂ is translated into pressure scales for T = 573 K. The non-hydrogenated surface was taken as a reference. Coordinatively unsaturated nickel atoms are marked by a red asterisk.

3.4. Phosphosulfide overlayer on the (101̄0) surface

Fig. 5a shows the phase diagram of the free energy change for the processes of adsorption and replacement by S atoms on the (101̄0) surface as a function of μ_P(T, p) and μ_S(T, p). The stability regions of the different surface compositions are indicated by the different shades of green. Fig. 5b–g show the top views of the corresponding surface structures. In all of these surfaces, S atoms are preferentially adsorbed at the T1 sites. In Fig. 6a, the red dotted line represents a typical HDS ratio, while the blue dotted line indicates the phase boundary between the bulk Ni₂P and Ni₃S₂. The (101̄0)-Ni₈P₄_2S, (101̄0)-Ni₈P₂S₂_2S, and (101̄0)-Ni₈S₄_2S surfaces are the most stable at a ratio of 0.01. All of these surfaces have the Ni₃ sites partially covered by S atoms in an alternating configuration. The findings of this study are consistent with those previously reported.⁴² It is observed that the (101̄0)-Ni₈S₄_2S surface is the most stable at very low PH₃ partial pressures. Therefore, our results suggest that the (101̄0)-Ni₈S₄_2S surface could be the most stable surface termination under HDS reaction conditions. On the (101̄0)-Ni₈S₄_2S surface, each T1 site consists of two Ni(1) atoms and one Ni(2) atom, with the Ni(1) atoms shared with the next T1 site. As a result, on the uncovered T1 sites, only the Ni(2) atoms are coordinatively unsaturated (Fig. 5f). Our findings are consistent with the EXAFS results by Oyama et al. that showed the square-pyramidal Ni(2) sites bind to sulfur with a lower Ni-S coordination.³ Therefore, the nickel phosphosulfide phase contains coordinatively unsaturated Ni(2)-type atoms available for the HDS reaction.

Fig. 5 (a) Phase diagram of the most stable (101̄0)-Ni_xP_y−zS_z_S_z’ surface composition as a function of μ_P(T, p) and μ_S(T, p). The dependence of the chemical potentials with the gas phase is translated into pressure scales for T = 573 K on the upper x-axis for the H₂S/H₂ mixture and on the right y-axis for PH₃ at p_H2 = 1 and 10 bar. The red dotted line corresponds to the typical HDS partial pressure ratio, . The blue dotted line indicates the phase boundary between the bulk Ni₂P and Ni₃S₂. Top view of the optimized structures of the (2 × 2) supercell: (b) (101̄0)-Ni₈P₄_2S, (c) (101̄0)-Ni₈P₄_4S, (d) (101̄0)-Ni₈P₂S₂_2S, (e) (101̄0)-Ni₈P₂S₂_4S, (f) (101̄0)-Ni₈S₄_2S, and (g) (101̄0)-Ni₈S₄_4S surfaces. Blue, light blue, pink, and yellow spheres are the Ni(1) nickel, Ni(2) nickel, phosphorus, and sulfur atoms, respectively. Coordinatively unsaturated Ni(2) nickel atoms are marked by a red asterisk.

Fig. 6 (a) Phase diagram of the most stable hydrogenated surface composition for the (101̄0)-Ni₈S₄_2S and (101̄0)-Ni₈S₄_4S surfaces as a function of μ_H(T, p) and μ_S(T, p). The dependence of the chemical potentials with the gas phase is translated into pressure scales for T = 573 K on the upper x-axis for the H₂S/H₂ mixture and on the right y-axis for H₂. The red dotted line corresponds to the typical HDS partial pressure ratio, . Top view of the optimized structures of the (2 × 2) supercell: (b) (101̄0)-Ni₈S₄_2SH, (c) (101̄0)-Ni₈S₄_2SH_2NiH, and (d) (101̄0)-Ni₈S₄_4SH hydrogenated surfaces. Blue, light blue, yellow, and white spheres are the Ni(1) nickel, Ni(2) nickel, sulfur, and hydrogen atoms, respectively. Coordinatively unsaturated Ni(1) and Ni(2) nickel atoms are marked by a red asterisk.

3.5. Dissociative adsorption of H₂ onto the phosphosulfide overlayer on the (101̄0) surface

Fig. 6a shows the phase diagram of the free energy change for the dissociative adsorption of one and two H₂ molecules on the (101̄0)-Ni₈S₄_2S and (101̄0)-Ni₈S₄_4S surfaces as a function of μ_H(T, p) and μ_S(T, p). The (101̄0)-Ni₈S₄_2S surface is taken as a reference. For example, the dissociative adsorption of two H₂ molecules in the (101̄0)-Ni₈S₄_4S surface was calculated as . Fig. 6b–d show the optimized structures for the most stable hydrogenated surfaces. Fig. S7 shows the configurations evaluated for the dissociative adsorption of one H₂ molecule on the (101̄0)-Ni₈S₄_2S. When H adsorbs to form an SH group, the S atom binds in a bridging position between the Ni(1) and Ni(2) atoms. Similarly, when the H atom adsorbs to form a NiH species, it occupies a bridging position between the Ni(1) and Ni(2) atoms. For the (101̄0)-Ni₈S₄_2S surface, the dissociative adsorption of one H₂ molecule produces two SH groups, leading to the (101̄0)-Ni₈S₄_2SH surface (Fig. 6b). When the SH groups form on this surface, two Ni(1) atoms become coordinatively unsaturated. On the other hand, the dissociative adsorption of two H₂ molecules produces two SH groups and two NiH groups, leading to the (101̄0)-Ni₈S₄_2SH_2NiH surface (Fig. 6c). For the (101̄0)-Ni₈S₄_4S surface, the dissociative adsorption of two H₂ molecules produces four SH groups, leading to the (101̄0)-Ni₈S₄_4SH surface (Fig. 6d). We also calculated the phase diagram of the free energy change for the dissociative adsorption of H₂ on the (101̄0)-Ni₈P₄_2S, (101̄0)-Ni₈P₄_4S, (101̄0)-Ni₈P₂S₂_2S, and (101̄0)-Ni₈P₂S₂_4S surfaces (Fig. S8). The configurations evaluated for the dissociative adsorption of one H₂ molecule on the (101̄0)-Ni₈P₄_2S and (101̄0)-Ni₈P₂S₂_2S surfaces are shown in Fig. S9 and S10, respectively.

As shown in Fig. 6a, keeping the ratio constant at 0.01 while increasing p_H2, results in two phase transitions at the most stable surface. First, the (101̄0)-Ni₈S₄_2S surface appears, and then the (101̄0)-Ni₈S₄_4SH surface. The hydrogenated surface (101̄0)-Ni₈S₄_4SH is stable at p_H2 above 10 bar. Hydrogenated surfaces resulting from the dissociative adsorption of H₂ on the (101̄0)-Ni₈S₄_2S surface become stable at lower ratios and higher p_H2. For instance, the (101̄0)-Ni₈S₄_2SH surface is stable at p_H2 above 10 bar, while the (101̄0)-Ni₈S₄_2SH_2NiH surface is stable at p_H2 above 100 bar. These results imply that surfaces with coordinatively unsaturated Ni(1) atoms, or with H atoms adsorbed in a bridging position between the Ni(1) and Ni(2) atoms, may only be stable under strong reducing conditions. This trend was also observed on the surfaces with 0% and 50% of replacement with S (Fig. S8). On the other hand, the energy required for the dissociative adsorption of two H₂ molecules on fully sulfur-covered surfaces at T = 573 K and p_H2 = 10 bar follows the trend: (101̄0)-Ni₈P₄_4SH (−0.07 eV) > (101̄0)-Ni₈P₂S₂_4SH (−0.14 eV) > (101̄0)-Ni₈S₄_4SH (−0.42 eV). These results suggest that replacing surface phosphorus atoms with sulfur promotes the H₂ dissociation process.

As mentioned above, Oyama et al. attribute the high HYD activity observed in small particles to the presence of S-resistant Ni(2) sites.³ Therefore, we used a (2 × 4) supercell to explore the stability of surfaces with SH groups and coordinatively unsaturated Ni(2) atoms. The (2 × 4) supercell size is an expansion along the y-axis. We evaluated the dissociative adsorption of one, two, and four H₂ molecules, corresponding to hydrogen coverages of 25%, 50%, and 100%, respectively. The surface compositions evaluated were the (101̄0)-Ni₁₆S₈_4S, (101̄0)-Ni₁₆S₈_6S, and (101̄0)-Ni₁₆S₈_8S surfaces, which correspond to sulfur coverages of 50%, 75%, and 100%, respectively. The supercell models are shown in Fig. S11. Fig. 7a shows that the (101̄0)-Ni₁₆S₈_2S_4SH and (101̄0)-Ni₁₆S₈_6SH_2NiH surfaces are the most stable hydrogenated surfaces at a ratio of 0.01. The (101̄0)-Ni₁₆S₈_6SH_2NiH surface is only stable at p_H2 above 10⁴ bar. Additionally, it is found that the (101̄0)-Ni₁₆S₈_4SH surface, which has coordinatively unsaturated Ni(1) atoms, is not stable under any reaction conditions (Fig. S12). Surfaces with H atoms adsorbed in a bridging position between the Ni(1) and Ni(2) atoms, such as the (101̄0)-Ni₁₆S₈_4SH_4NiH surface, are only stable under strong reducing conditions. At high ratios, the (101̄0)-Ni₁₆S₈_8SH surface is the most stable hydrogenated surface. These results suggest that hydrogenated surfaces, such as (101̄0)-Ni₁₆S₈_2S_4SH, which have coordinatively unsaturated Ni(2) atoms, may be stable under HDS reaction conditions. Note that the SH groups surround the Ni(2) atoms (Fig. 7b). These findings may help explain the high HYD selectivity observed at Ni(2) sites, as reported by Oyama et al.³ Organosulfur compounds adsorbed on Ni(2) sites could react with H atoms from adjacent SH groups, which may result in the formation of hydrogenated intermediates and subsequent C–S bond cleavage. A detailed kinetic analysis of the HDS reaction mechanism, however, is beyond the scope of the present work.

Fig. 7 (a) Phase diagram of the most stable hydrogenated surface composition for the (101̄0)-Ni₁₆S₈_4S, (101̄0)-Ni₁₆S₈_6S, and (101̄0)-Ni₁₆S₈_8S surfaces as a function of μ_H(T, p) and μ_S(T, p). The dependence of the chemical potentials with the gas phase is translated into pressure scales for T = 573 K on the upper x-axis for the H₂S/H₂ mixture and on the right y-axis for H₂. The red dotted line corresponds to the typical HDS partial pressure ratio, . Top view of the optimized structures of the (2 × 4) supercell of the hydrogenated surfaces: (b) (101̄0)-Ni₁₆S₈_2S_4SH, (c) (101̄0)-Ni₁₆S₈_4SH_4NiH, (d) (101̄0)-Ni₁₆S₈_6SH_2NiH, and (e) (101̄0)-Ni₁₆S₈_8SH hydrogenated surfaces. Blue, light blue, yellow, and white spheres are the Ni(1) nickel, Ni(2) nickel, sulfur, and hydrogen atoms, respectively. Coordinatively unsaturated Ni(2) nickel atoms are marked by a red asterisk.

To gain insight into the chemical nature of the hydrogen species on the surface, we performed Bader charge and electron density difference analyses. The formal charges of the different adsorbed species on the surface are H⁻, SH⁻, and S²⁻ The stabilization of hydrogen as a proton involves an oxidation process, as it is accompanied by the oxidation of S atoms within the SH group. In contrast, the stabilization of hydrogen as a hydride corresponds to a reduction process.⁵³ The Bader charges (q) of the surface atoms for the (101̄0)-Ni₁₆S₈_6S, (101̄0)-Ni₁₆S₈_2S_4SH and (101̄0)-Ni₁₆S₈_6SH_2NiH surfaces are summarized in Table 1. The hydrogen atoms are characterized by a different Bader charge, which reflects their chemical identities. For the (101̄0)-Ni₁₆S₈_2S_4SH surface, when H adsorbs as a proton (H⁺) to form SH groups, each H atom carries a low positive charge of +0.04e. The charge on the S atoms is found to be more negative as S anions (−0.48e) than as SH groups (−0.30e). For the (101̄0)-Ni₁₆S₈_6SH_2NiH surface, when H adsorbs as a hydride (H⁻), each H atom carries a negative charge of −0.17e. Cai et al. reported the Bader charge analysis of hydrogen atoms adsorbed on the RuS₂ catalyst. They found that hydrogen atoms bonded to sulfur atoms carry a positive charge ranging from +0.02e to +0.15e, whereas those adsorbed on ruthenium atoms exhibit a negative charge of approximately −0.20e.¹⁸ The Bader charges obtained in our study are consistent with their results in terms of both magnitude and sign. We also calculated the EDD to illustrate the redistribution of electron density due to the interaction between these species and the surface. Fig. 8a and b show the EDD plots for the (101̄0)-Ni₁₆S₈_2S_4SH and (101̄0)-Ni₁₆S₈_6SH_2NiH surfaces, respectively. The densities of the S, H, and SH species were used as adsorbate references. In these plots, the blue regions represent electron density accumulation, and the red regions represent electron density depletion. For instance, the electron density accumulates on the adsorbed H atoms on the (101̄0)-Ni₁₆S₈_6SH_2NiH surface, which is consistent with a hydridic character. Additionally, the redistribution of electron density differs between the SH groups due to the local chemical environment that each SH group experiences in its immediate surroundings (Fig. 8b). Overall, these hydrogen species exhibit either proton or hydride character, depending on the electronegativity of their ligands (sulfur or metal).¹⁸

Table 1 Bader charges (q) of the species adsorbed on the (101̄0)-Ni₁₆S₈_2S_4SH, and (101̄0)-Ni₁₆S₈_6SH_2NiH surfaces

Surface	q S (e)	q SH (e)	q S(SH) (e)	q H(SH) (e)	q H (e)
a q S charge of the adsorbed S atoms. b q SH charge of the SH groups. c q S(SH) charge of the S atoms in the SH groups. d q H(SH) charge of the H atoms in the SH groups. e q H charge of the H atoms in the NiH especies.
(101̄0)-Ni16S8_2S_4SH	−0.48	−0.26	−0.30	+0.04	—
(101̄0)-Ni16S8_6SH_2NiH	—	−0.26	−0.28	+0.02	−0.17
		−0.26	−0.31	+0.05

---

Fig. 8 EDD plots for the adsorption of the S, H, and SH species on the (a) (101̄0)-Ni₁₆S₈_2S_4SH and (b) (101̄0)-Ni₁₆S₈_6SH_2NiH surfaces. The blue and red regions represent the electron density gain and loss, respectively. The EDD was plotted using an isodensity level of 0.0035e Borh⁻³. Blue, light blue, yellow, and white spheres are the Ni(1) nickel, Ni(2) nickel, sulfur, and hydrogen atoms, respectively.

4. Conclusions

In this work, we carried out a comprehensive theoretical analysis of the structural and electronic properties of the active nickel and hydrogen species on the phosphosulfide phase formed on the (0001) and (101̄0) surfaces under HDS reaction conditions. The fully sulfided (0001) surface of Ni₂P, composed of Ni(1)-type atoms, closely resembles the (111) surface of the low-activity sulfide Ni₃S₂, suggesting that surfaces with spatially isolated Ni₃ ensembles that are not covered by sulfur are unlikely to be thermodynamically stable. Additionally, hydrogenated (0001) surfaces are not stable because, upon SH group formation, the sulfur atom binds to two of the three Ni(1) atoms within the Ni₃ site, leaving the third Ni(1) atom coordinatively unsaturated.

On the (101̄0) surface, the SH groups were identified as the stable form of hydrogen adsorption. The (101̄0)-Ni₁₆S₈_2S_4SH surface, which has coordinatively unsaturated Ni(2) atoms surrounded by SH groups, was the most stable under HDS conditions. Hydrogen stabilization as a hydride, with H atoms adsorbed in a bridging position between the Ni(1) and Ni(2) atoms, only occurred under strong reducing conditions. Our results indicated that the formation of the proton-hydride pair may not occur under reaction conditions. Instead, only the proton species within SH groups are likely to play an active role in facilitating hydrogenation reactions on the nickel phosphosulfide overlayer formed on the Ni₂P catalyst. It should be emphasized that the surface structures identified in this work correspond to thermodynamically stable surfaces under equilibrium assumptions, rather than a unique description of the dynamic catalytic surface under reaction conditions.

Overall, these findings offer insight into processes occurring at the catalyst surface that may influence HYD selectivity at Ni(2) sites. Organosulfur species adsorbed at these sites may interact with protons from nearby SH groups, potentially forming hydrogenated intermediates and facilitating subsequent C–S bond hydrogenolysis.

Conflicts of interest

There are no conflicts to declare.

Data availability

The supporting data has been provided as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6cp00205f.

Acknowledgements

The DFT calculations were performed using computational resources from Laboratorio de Química Física y Catálisis Computacional (FONACIT Grant 20240PGP66) and Laboratorio de Fisicoquímica Teórica de Materiales (FONACIT Grant 2023PGP320).

References

1. K. K. Bando, T. Wada, T. Miyamoto, K. Miyazaki, S. Takakusagi, Y. Koike, Y. Inada, M. Nomura, A. Yamaguchi, T. Gott, S. Ted Oyama and K. Asakura, Combined in situ QXAFS and FTIR analysis of a Ni phosphide catalyst under hydrodesulfurization conditions, J. Catal., 2012, 286, 165–171.
2. L. Yang, X. Li, A. Wang, R. Prins, Y. Wang, Y. Chen and X. Duan, Hydrodesulfurization of 4,6-dimethyldibenzothiophene and its hydrogenated intermediates over bulk Ni₂P, J. Catal., 2014, 317, 144–152.
3. Y.-K. Lee and S. T. Oyama, Sulfur resistant nature of Ni₂P catalyst in deep hydrodesulfurization, Appl. Catal., A, 2017, 548, 103–113.
4. P. J. Topalian, D. R. Liyanage, S. J. Danforth, A. I. d’Aquino, S. L. Brock and M. E. Bussell, Effect of Particle Size on the Deep HDS Properties of Ni₂P Catalysts, J. Phys. Chem. C, 2019, 123, 25701–25711.
5. C. E. Miles, T. R. Carlson, B. J. Morgan, P. J. Topalian, J. R. Schare and M. E. Bussell, Hydrodesulfurization Properties of Nickel Phosphide on Boron-treated Alumina Supports, ChemCatChem, 2020, 12, 4939–4950.
6. X. Zhou, X. Li, R. Prins, A. Wang, L. Wang, S. Liu and Q. Sheng, Desulfurization of 2-phenylcyclohexanethiol over transition-metal phosphides, J. Catal., 2020, 383, 331–342.
7. R. H. Bowker, G. H. Layan Savithra, B. A. Carrillo, K. G. Hubach, T. McDonald, S. L. Brock and M. E. Bussell, Effect of particle size on the sulfur resistance of nickel phosphide hydrodesulfurization catalysts, J. Catal., 2023, 425, 70–79.
8. S. Oyama, Effect of Phosphorus Content in Nickel Phosphide Catalysts Studied by XAFS and Other Techniques, J. Catal., 2002, 210, 207–217.
9. S. T. Oyama, X. Wang, Y.-K. Lee and W.-J. Chun, Active phase of Ni₂P/SiO₂ in hydroprocessing reactions, J. Catal., 2004, 221, 263–273.
10. Y. Shu, Y. Lee and S. Oyama, Structure-sensitivity of hydrodesulfurization of 4,6-dimethyldibenzothiophene over silica-supported nickel phosphide catalysts, J. Catal., 2005, 236, 112–121.
11. T. Kawai, K. Bando, Y. Lee, S. Oyama, W. Chun and K. Asakura, EXAFS measurements of a working catalyst in the liquid phase: an in situ study of a Ni₂P hydrodesulfurization catalyst, J. Catal., 2006, 241, 20–24.
12. S. T. Oyama, T. Gott, K. Asakura, S. Takakusagi, K. Miyazaki, Y. Koike and K. K. Bando, In situ FTIR and XANES studies of thiophene hydrodesulfurization on Ni2P/MCM-41, J. Catal., 2009, 268, 209–222.
13. X. Duan, Y. Teng, A. Wang, V. Kogan, X. Li and Y. Wang, Role of sulfur in hydrotreating catalysis over nickel phosphide, J. Catal., 2009, 261, 232–240.
14. X. Li, Z. Sun, A. Wang, Y. Wang, X. Duan and Y. Chen, XAS study of Ni₂P/MCM-41 passivated by O₂/He and H₂S/H₂, Catal. Commun., 2014, 43, 21–24.
15. S. Tian, X. Li, A. Wang, R. Prins, Y. Chen and Y. Hu, Facile Preparation of Ni₂P with a Sulfur-Containing Surface Layer by Low-Temperature Reduction of Ni₂P₂S₆, Angew. Chem., Int. Ed., 2016, 55, 4030–4034.
16. M. He, T. Li, X. Li, A. Wang, Q. Sheng, S. Shang and Z. Yu, Preparation of Ni₂P with a Surface Nickel Phosphosulfide Layer by Reduction of Mixtures of Na₄P₂S₆ and NiCl₂, ChemCatChem, 2024, e202400823.
17. S. Rangarajan and M. Mavrikakis, DFT Insights into the Competitive Adsorption of Sulfur- and Nitrogen-Containing Compounds and Hydrocarbons on Co-Promoted Molybdenum Sulfide Catalysts, ACS Catal., 2016, 6, 2904–2917.
18. H. Cai, R. Schimmenti, H. Nie, M. Mavrikakis and Y.-H. C. Chin, Mechanistic Role of the Proton–Hydride Pair in Heteroarene Catalytic Hydrogenation, ACS Catal., 2019, 9, 9418–9437.
19. N. Y. Topsoe and H. Topsoe, FTIR Studies of Mo/Al2O3-Based Catalysts, J. Catal., 1993, 139, 641–651.
20. G. Berhault, M. Lacroix, M. Breysse, F. Maugé, J.-C. Lavalley, H. Nie and L. Qu, Characterization of Acidic Sites of Silica-Supported Transition Metal Sulfides by Pyridine and 2,6 Dimethylpyridine Adsorption: Relation to Activity in CH₃SH Condensation, J. Catal., 1998, 178, 555–565.
21. S. Oyama and Y. Lee, The active site of nickel phosphide catalysts for the hydrodesulfurization of 4,6-DMDBT, J. Catal., 2008, 258, 393–400.
22. A. Nelson, M. Sun and A. Junaid, On the structure and composition of the phosphosulfide overlayer on Ni₂P at hydrotreating conditions, J. Catal., 2006, 241, 180–188.
23. L. L. D. Thi, P. N. K. Tuyen and T. Y. Vu, Flower-like S-doped-Ni₂P mesoporous nanosheet-derived self-standing electrocatalytic electrode for boosting hydrogen evolution, Nanotechnology, 2020, 31, 465401.
24. X. Wang, Q. Hu, G. Li, S. Wei, H. Yang and C. He, Regulation of the adsorption sites of Ni₂P by Ru and S co-doping for ultra-efficient alkaline hydrogen evolution, J. Mater. Chem. A, 2021, 9, 15648–15653.
25. Y. Wu, X. Chen, L. Su, Q. Wang and S. Ren, A sulfur-doped Ni₂P electrocatalyst for the hydrogen evolution reaction, New J. Chem., 2022, 46, 7675–7681.
26. M. Šarić, P. G. Moses and J. Rossmeisl, Relation between Hydrogen Evolution and Hydrodesulfurization Catalysis, ChemCatChem, 2016, 8, 3334–3337.
27. A. J. McCue and J. A. Anderson, Sulfur as a catalyst promoter or selectivity modifier in heterogeneous catalysis, Catal. Sci. Technol., 2014, 4, 272–294.
28. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
29. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
30. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
31. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
32. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132, 154104.
33. W. Hujo and S. Grimme, Comparison of the performance of dispersion-corrected density functional theory for weak hydrogen bonds, Phys. Chem. Chem. Phys., 2011, 13, 13942.
34. J. P. P. Ramalho, J. R. B. Gomes and F. Illas, Accounting for van der Waals interactions between adsorbates and surfaces in density functional theory based calculations: selected examples, RSC Adv., 2013, 3, 13085.
35. S. Banerjee, A. Kakekhani, R. B. Wexler and A. M. Rappe, Relationship between the Surface Reconstruction of Nickel Phosphides and Their Activity toward the Hydrogen Evolution Reaction, ACS Catal., 2023, 13, 4611–4621.
36. S. Rundqvist, M. Yhland, R. Dahlbom, J. Sjövall, O. Theander and H. Flood, X-Ray Investigations of Mn₃P, Mn₂P, and Ni₂P, Acta Chem. Scand., 1962, 16, 992–998.
37. M. E. Fleet, The crystal structure of heazlewoodite, and metallic bonds in sulfide minerals, Am. Mineral., 1977, 62, 341–345.
38. R. F. W. Bader, A quantum theory of molecular structure and its applications, Chem. Rev., 1991, 91, 893–928.
39. K. Momma and F. Izumi, VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Crystallogr., 2011, 44, 1272–1276.
40. W. Papawassiliou, J. P. Carvalho, N. Panopoulos, Y. Al Wahedi, V. K. S. Wadi, X. Lu, K. Polychronopoulou, J. B. Lee, S. Lee, C. Y. Kim, H. J. Kim, M. Katsiotis, V. Tzitzios, M. Karagianni, M. Fardis, G. Papavassiliou and A. J. Pell, Crystal and electronic facet analysis of ultrafine Ni₂P particles by solid-state NMR nanocrystallography, Nat. Commun., 2021, 12, 4334.
41. J. He, Á. Morales-García, O. Bludský and P. Nachtigall, The surface stability and equilibrium crystal morphology of Ni₂P nanoparticles and nanowires from an ab initio atomistic thermodynamic approach, CrystEngComm, 2016, 18, 3808–3818.
42. A. B. Vidal, J. L. Peña-Mena, O. Hurtado-Aular, R. Añez and A. Sierraalta, Unraveling the Structure and Surface Chemistry of the Phosphosulfide Phase Formed on Ni₂P under Hydrodesulfurization Reaction Conditions: A DFT Study, J. Phys. Chem. C, 2022, 126, 14187–14200.
43. A. Aguinaga, M. Montes, P. Malet, M. J. Capitán, I. Carrizosa and J. A. Odriozola, Extended X-ray absorption fine structure study of sulphur poisoned Ni/SiO₂ catalysts, Appl. Catal., A, 1994, 110, 197–205.
44. J. R. Rostrup-Nielsen, J. Sehested and J. K. Nørskov, Advances in Catalysis, Elsevier, 2002, vol. 47, pp. 65–139.
45. X. Chen, J. Jiang, F. Yan, K. Li, S. Tian, Y. Gao and H. Zhou, Dry Reforming of Model Biogas on a Ni/SiO₂ Catalyst: Overall Performance and Mechanisms of Sulfur Poisoning and Regeneration, ACS Sustainable Chem. Eng., 2017, 5, 10248–10257.
46. R. Prins and M. E. Bussell, Metal Phosphides: Preparation, Characterization and Catalytic Reactivity, Catal. Lett., 2012, 142, 1413–1436.
47. T. Pecoraro and R. Chianelli, Hydrodesulfurization catalysis by transition metal sulfides, J. Catal., 1981, 67, 430–445.
48. J. Quartararo, S. Mignard and S. Kasztelan, Hydrodesulfurization and hydrogenation activities of alumina-supported transition metal sulfides, J. Catal., 2000, 192, 307–315.
49. G. Kishan, L. Coulier, V. H. J. de Beer, J. A. R. van Veen and J. W. Niemantsverdriet, Sulfidation and Thiophene Hydrodesulfurization Activity of Nickel Tungsten Sulfide Model Catalysts, Prepared without and with Chelating Agents, J. Catal., 2000, 196, 180–189.
50. Y. Lee and S. Oyama, Bifunctional nature of a SiO₂-supported Ni₂P catalyst for hydrotreating: EXAFS and FTIR studies, J. Catal., 2006, 239, 376–389.
51. R. R. Chianelli, T. A. Pecoraro, T. R. Halbert, W. H. Pan and E. I. Stiefel, Transition metal sulfide catalysis: relation of the synergic systems to the periodic trends in hydrodesulfurization, J. Catal., 1984, 86, 226–230.
52. Y. Aray, D. Vega, J. Rodriguez, A. B. Vidal, M. Elena Grillo and S. Coll, First-Principles Study of Low Miller Index Ni₃S₂ Surfaces in Hydrotreating Conditions, J. Phys. Chem. B, 2009, 113, 3058–3070.
53. R. Khare, R. Weindl, A. Jentys, K. Reuter, H. Shi and J. A. Lercher, Di- and Tetrameric Molybdenum Sulfide Clusters Activate and Stabilize Dihydrogen as Hydrides, JACS Au, 2022, 2, 613–622.

---