N-Doped CrS2 Monolayer as a Highly-Efficient Catalyst for Oxygen Reduction Reaction: A Computational Study.

Searching for low-cost and highly-efficient oxygen reduction reaction (ORR) catalysts is crucial to the large-scale application of fuel cells. Herein, by means of density functional theory (DFT) computations, we proposed a new class of ORR catalysts by doping the CrS2 monolayer with non-metal atoms (X@CrS2, X = B, C, N, O, Si, P, Cl, As, Se, and Br). Our results revealed that most of the X@CrS2 candidates exhibit negative formation energy and large binding energy, thus ensuring their high stability and offering great promise for experimental synthesis. Moreover, based on the computed free energy profiles, we predicted that N@CrS2 exhibits the best ORR catalytic activity among all considered candidates due to its lowest overpotential (0.41 V), which is even lower than that of the state-of-the-art Pt catalyst (0.45 V). Remarkably, the excellent catalytic performance of N@CrS2 for ORR can be ascribed to its optimal binding strength with the oxygenated intermediates, according to the computed linear scaling relationships and volcano plot, which can be well verified by the analysis of the p-band center as well as the charge transfer between oxygenated species and catalysts. Therefore, by carefully modulating the incorporated non-metal dopants, the CrS2 monolayer can be utilized as a promising ORR catalyst, which may offer a new strategy to further develop eligible electrocatalysts in fuel cells.

1. Introduction

The rising energy demand and depletion of fossil fuels have attracted increasing interest in alternative energy sources. This is because the traditional fossil fuels are not renewable and will produce harmful combustion products, such as CO, CO2, NO, and SO2, which thus pose a serious challenge to human health and environmental protection [1,2,3,4,5,6]. In this regard, electrochemical energy storage and conversion technologies, such as fuel cells, were regarded as efficient and clean devices for reducing the global energy shortage [7,8]. However, the efficiency of these energy-related devices was greatly determined by the oxygen reduction reaction (ORR) due to its sluggish kinetics and high overpotential [9,10,11,12,13]. At present, Pt-based materials represent the most promising ORR catalysts [14,15,16,17]. However, their practical application is severely limited by their high price and low abundance. Thus, the search for highly efficient, low-cost alternative electrocatalysts is urgent to boost the ORR.

In recent years, single-layer transition metal dichalcogenides (TMDs) have become a hotspot of theoretical and experimental research owing to their large surface area, stability, and peculiar electronic structure [18,19,20,21,22,23]. Thus, TMD-based materials have been widely employed in nanoelectronics, nanophotonics, the absorber layer in solar cells, anode materials, field-effect transistors, and so on [24,25,26,27,28]. More interestingly, the heteroatom doped TMDs have attracted great interest because they have more active sites, low cost, high stability, and high efficiency as potential electrocatalytic catalysts [29,30,31,32,33,34,35]. For example, Xiong et al. reported that co-doped MoS2 possesses high catalytic activity for oxygen evolution reaction [29]. Moreover, Lv et al. proposed that N-doped MoS2 displayed a high faradaic efficiency and low onset overpotential for CO2 electroreduction to CO [30]. Theoretically, Li et al. demonstrated that B-doping causes the VS2 monolayer to exhibit outstanding catalytic activity towards nitrogen reduction reactions [31]. Moreover, Singh et al. reported that the ORR activity of N-doped WS2 monolayers should originate from the introduction of spin density caused by N doping [32]. In addition, Tian et al. showed that Co doped 1T-TiS2 exhibits significantly enhanced performance toward the ORR [33].

As a representative of 2D TMDs, the chromium disulfide (CrS2) monolayer was synthesized via the chemical vapor deposition (CVD) method in 2019 [21]. Interestingly, due to its unique properties, the CrS2 monolayer has wide applications in spintronic devices [36,37]. For example, Chen et al. suggested that CrS2 has the most diverse electronic and magnetic properties: antiferromagnetic (AFM) metallic, non-magnetic (NM) semiconductor, and ferromagnetic (FM) semiconductor with a Curie temperature of ~1000 K [36]. Moreover, Zhang et al. reported that the magnetic properties of CrC2 monolayer can be effectively tuned by doping metal atoms [37].

Inspired by the interesting properties of the CrS2 monolayer, in this work, we explored the potential of several non-metal-doped CrS2 monolayers (X@CrS2, X = B, C, N, O, Si, P, Cl, As, Se, and Br) as ORR catalysts by performing comprehensive DFT computations. According to the computed formation energies and binding energies of these X@CrS2 systems, we suggested that these doped CrS2 candidates are very likely to be synthesized in experiments and possess extremely high stability. Moreover, the N@CrS2 catalyst was screened out as an eligible ORR catalyst with a rather low overpotential of 0.41 V, which originates from its optimal binding strengths with oxygenated intermediates based on the linear scaling relationships and volcano plot.

2. Materials and Methods

All computations were performed based on spin-polarized density functional theory, as implemented in the Vienna ab initio simulation package (VASP 541) [38,39]. The projector augmented wave (PAW) was adopted to describe the interactions between ions and electrons [40,41]. The generalized gradient approximation (GGA) in the form of Perdew–Burke–Ernzerhof (PBE) was adopted for the exchange-correlation functional [42]. Notably, the PBE functional is, nowadays, the most commonly used functional for solid-state calculations, which remains the best for the solids containing 3d-transition elements [43,44,45]. In particular, our purpose in the present work is to search for the ideal ORR catalysts among various candidates, and we thus mainly focus on the catalytic tendency of these candidates for ORR. Thus, although some more modern and better functionals have been proposed in recent years [46,47], the identified scaling relations of oxygenated species on various electrocatalysts and the derived conclusions (such as the corresponding catalytic activity) will not change. The DFT+D3 method in the Grimme scheme was employed to describe the possible weak interactions of the oxygenated species and the catalysts [48]. An energy cutoff of 500 eV was used for the plane wave (Table S1), and the convergence criteria of force and energy were set to 0.02 eV Å−1 and 10−5 eV, respectively. A k-point of 3 × 3 × 1 was sampled in the Brillouin zones (Table S1), and the vacuum space was set to 15 Å.

To estimate the ORR catalytic performance of these X@CrS2 materials, the change in the Gibbs free-energy change (ΔG) of each elementary reaction step during ORR was computed using the computational hydrogen electrode (CHE) model [15,49]: ΔG = ΔE + ΔZPE − TΔS + ΔGU, in which ΔE is the reaction energy directly obtained from DFT computations, and ΔZPE and ΔS represent the difference of zero-point energy and entropy, respectively, which can be derived from the computations on the vibrational frequencies and the standard thermodynamic data. ΔGU = −eU, where U is the applied potential. Since DFT generally fails to obtain the energy of the O2 molecule, the free energy of O2 (GO2) will be computed via the energies of H2O and H2, namely GO2 = GH2O − 2GH2 + 4.92 eV. Furthermore, the ORR’s catalytic activity was evaluated by computing the corresponding overpotential (η) according to the following equations: η = max{ΔG1, ΔG2, ΔG3, ΔG4}/e + U0, where ΔG represents the free energy change in each elementary reaction, including the formation of OOH*, O*, and OH* and the desorption of OH*, and U0 is the equilibrium potential of ORR, which is equal to 1.23 V. According to the above definition, a smaller η value corresponds to a higher catalytic activity for ORR.

3. Results

3.1. Structures, Stabilities, and Properties of X@CrS2 Catalysts

First, we investigated the structures and electronic properties of the pristine CrS2 monolayer with the 2H phase. As shown in Figure 1a, each Cr center is prismatically coordinated by six surrounding S atoms, with the S atoms in the upper layer lying directly above those of the lower layer. Furthermore, the optimized lattice constant for the CrS2 monolayer is 3.04 Å, while the lengths of the formed Cr-S bonds are 2.28 Å. Moreover, the CrS2 monolayer is a direct semiconductor with the band gap of 0.93 eV, in which the conduction band minimum and the valence band maximum located at the K point mainly originate from the contribution of 3d orbitals of the Cr atom (Figure 1b). Notably, these above results on the structure and properties of a pristine CrS2 monolayer are consistent with previous theoretical reports [50], indicating the accuracy of our computational methods to describe the behavior of the CrS2 monolayer.

Figure 1

(a) The obtained stable configurations and (b) computed band structure for the pristine CrS2 monolayer.

Based on the optimized CrS2 monolayer, various non-metal atoms were introduced to substitute one of the S atoms within the CrS2 monolayer to construct the doped CrS2 system. After fully geometrical relaxation, the structures of the X@CrS2 monolayer are presented in Figure 2, while their corresponding structural parameters are summarized in Table 1. The results showed that, after the introduction of dopants, the structure of the CrS2 monolayer is changed in various ways, which are highly dependent on the difference in the radius between dopants and S atoms. In detail, for B, C, N, and O dopants with a smaller radius than the S atom, the introduced non-S atoms are concave from the CrS2 plane, whereas the doping of Si, P, Cl, As, Se, and Br induces an unchanged or slightly outward structure due to their comparable or slightly larger radius.

Figure 2

The obtained stable configurations for X@CrS2 monolayers (N@CrS2 was chosen as a representative).

Table 1

The computed bond’s length of X-Cr (dX-Cr, Å), formation energies (E eV), binding energies (E eV), magnetic moment (μ), charge transfer (Q, e−), and band gaps (E, eV) for various X@CrS2 monolayers.

	dX-Cr	Ef		Ebind	μB	Q	Egap
		S-rich	Cr-rich
pristine	2.28	/	/	/	0.00	/	0.93
B	1.97	2.06	0.94	−5.62	1.00	0.09	0.20
C	1.92	0.76	−0.36	−7.16	0.01	0.77	0.90
N	1.87	0.16	−0.96	−5.75	1.00	0.94	0.25
O	1.93	−2.70	−3.83	−6.96	0.01	0.95	0.96
Si	2.47	−3.40	−4.53	−3.93	2.00	0.32	0.27
P	2.35	0.15	−0.97	−4.15	1.00	0.27	0.95
Cl	2.41	−2.48	−3.61	−3.14	1.00	0.47	0.18
As	2.49	−0.02	−1.14	−3.70	0.96	0.04	0.00
Se	2.42	−1.17	−2.30	−4.76	0.01	0.28	0.98
Br	2.56	−0.46	−1.58	−2.52	1.00	0.31	0.19

To assess the experimental feasibility of these doped CrS2 monolayers, we computed their formation energies (E) based on the following equation: E = EX@CrS2 − ECrS2 + μ − μ [51,52], where the EX@CrS2 and ECrS2 are the total energies of doped and pristine CrS2 monolayers, respectively. μ and μ are the chemical potentials of a doping atom and S atom. Notably, μ can be derived from the energy per atom in their bulk X materials. In contrast, μ is dependent on the growth conditions, in which two extreme cases were considered, including Cr-rich and S-rich conditions. According to previous studies [32,53], the μCr and μ should meet the relationship from the thermodynamic equilibrium condition: μCrS2 = μCr + 2μ, where μCrS2 represents the total energy of the pristine CrS2 monolayer per formula unit. Under Cr-rich conditions, the μCr was taken from the bulk Cr. Thus, μ can be obtained by: μS = (μCrS2 − μCr)/2. On the other hand, under S-rich conditions, the μ can be obtained from its bulk S8 state, and then μCr can be determined by: μCr = μCrS2 − 2μS. The computed formation energies of these X@CrS2 materials were listed in Table 1. We can see that the formation energies for B doping are all positive values, indicating that the formation of this kind of doped CrS2 monolayer is difficult. On the contrary, when O, Si, Cl, As, Se, and Br dopants were introduced into the CrS2 monolayer, their Ef values were always negative, implying notable feasibility for their experimental synthesis. Notably, the CrS2 monolayers doped by C, N, and P atoms are more likely to form under Cr-rich conditions due to their negative formation energies. Obviously, by carefully tuning the reaction conditions, these doped CrS2 monolayers can be easily fabricated in experiments, except for B-doping.

To further explore the structural stability of X@CrS2, we computed the binding energies (E) of the introduced dopants on the CrS2 monolayer according to the following equation: EX@CrS2 − ECrS2 − EX, where EX@CrS2, ECrS2, and EX represent the total electronic energies of doped CrS2 monolayers, the defective CrS2 substrate, and isolated X atoms in their most stable phase. It can be observed that the E values of these dopants on the CrS2 substrate range from −2.52 eV of Br to −7.16 eV of C, indicative of the strong binding strength, thus ensuring their high stability. The stability of X@CrS2 was further evaluated by using AIMD simulations, where N doping was taken as a representative. As shown from Figure S1, there is still no significant distortion of the geometric structure at 500 K, indicative of its high thermodynamic stability. It should be noted that incorporation of various non-metal atoms into 2D TMDs has been already realized experimentally. For example, Li et al. reported a simple, facile, and effective strategy to fabricate an N-doped MoS2 nanosheet [54], while Xin et al. demonstrated the fabrication of a P-doped MoS2 nanosheet by the one-step hydrothermal method [55]. Therefore, we strongly believe that the as-designed doped CrS2 monolayer holds great promise for synthesis.

Next, we explored the magnetic and electronic properties of these X@CrS2 systems. Our results showed that the pristine CrS2 monolayer is a nonmagnetic material. After introducing these dopants into the CrS2 monolayer, however, different magnetic behaviors can be observed: (1) C-, O-, and Se-doped CrS2 systems still retain their nonmagnetic states due to the absence of unpaired electrons; (2) as for B, N, P, As, or Br-doped systems, the total magnetic moment is close to 1.00 μ, while the Si-doped CrS2 monolayer has a magnetic moment of 2.00 μ. Moreover, we noted that these dopants will loot different amounts of electrons (0.04~0.95 |e−|) from the CrS2 monolayer, except for the Si dopant, which can donate about 0.32 electrons to the CrS2 substrate. As a result, the band gaps of the doped CrS2 monolayer are reduced to different degrees due to the introduction of impurity levels. For example, the N-doped CrS2 system exhibits a smaller band gap of 0.25 eV than that of the pristine one (0.99 eV), suggesting the enhanced electrical conductivity, which may facilitate the rapid charge transfer in electrocatalysis [56]. It is noted that the GGA method usually underestimates band gaps, whereas the hybrid functional method, such as Heyd–Scuseria–Ernzerhof (HSE06), can rectify the band gap [57]. For example, the computed band gaps of the well-established MoS2 monolayer using HSE06 and PBE methods are 2.20 and 1.79 eV [58], respectively. Despite the fact that a more accurate band gap value can be predicted by means of the HSE06 method, the corresponding computational cost is extremely high for the 10 doped CrS2 monolayers, which consist of 25 Cr and 50 S atoms. We will mainly focus on the trend of the band gaps of the CrS2 monolayer after non-metal doping.

3.2. ORR Catalytic Activity

After knowing that these non-metal doped CrS2 candidates exhibit great potential for experimental synthesis, high stability, diverse magnetic moment, and enhanced conductivity, we further explored their catalytic activities towards ORR.

Typically, the O2 molecule can be reduced to H2O along a four-electron (4e−) pathway (Figure 3a): (1) * + O2 (g) + H+ + e− ⭢ OOH*; (2) OOH* + H+ + e− ⭢ O* + H2O; (3) O* + H+ + e− ⭢ OH*; (4) OH* + H+ + e− ⭢ H2O (l) [59,60]. We again took the N@CrS2 monolayer as the representative to compute the ∆G values of the above four elementary steps. As shown in Figure 3b, we found that the formation of OOH* species on the N site of the N@CrS2 catalyst is exothermic, with a ∆G value of −0.84 eV. The length of the formed N–O bond is 1.41 Å, and the O–O bond is elongated to 1.48 Å as compared with that of the free O2 molecule (1.23 Å). Subsequently, the approach of a second hydrogen induces the dissociation of OOH* into (O* + H2O) or the formation of H2O2. Remarkably, the OOH* formation is exothermic by 2.43 eV, which is much larger than that of H2O2 formation (0.16 eV), suggesting that N@CrS2 shows a rather high selectivity towards the 4e− pathway by greatly suppressing the competing 2e− one. Once the O* species is formed, it can be further hydrogenated to form OH* and the second H2O with the ∆G values of −0.83 and −0.82 eV, respectively. According to the CHE model, the final step (i.e., OH* desorption) is identified as the potential-determining step (PDS). Thus, the limiting potential is 0.82 V, which is the smallest applied voltage to make the whole reaction still exergonic, corresponding to the overpotential of 0.41 V.

Figure 3

(a) The involved reaction pathway and (b) the computed free energy diagrams for ORR pathways on the N@CrS2 monolayer at zero and applied voltages.

As for other X@CrS2 candidates, the computed free adsorption energies of oxygenated species and the corresponding ηORR are summarized in Table S2. We found that the computed ηORR increases in the order of As@CrS2 (0.97 V) < C@CrS2 (1.20 V) < O@CrS2 (1.53 V) < Se@CrS2 (1.61 V) < Br@CrS2 (1.74 V) ≈ Cl@CrS2 (1.78 V) < P@CrS2 (2.18 V) < B@CrS2 (2.42 V) < Si@CrS2 (2.82 V), as shown in Figure S2, which are all higher than that of N@CrS2 (0.41 V). In particular, the overpotential of N@VS2 (0.41 V) is even lower than that of the well-established Pt catalyst (0.45 V) [15], implying its excellent ORR catalytic activity. As ORR usually occurs in aqueous solution, we also studied the solvent effect on the ORR activity of the N@CrS2 catalyst using the implicit solvation model in VASPsol with a dielectric constant of 80 [61]. It can be seen from Figure S3 that the PDS locates at the last step, and the computed η value is 0.48 V, which is comparable to that of the value without the solvent effect (0.41 V), implying that the solvent effects left the superior ORR catalytic performance of the N@CrS2 monolayer nearly unchanged. Although implicit solvent models can offer fast and inexpensive starting points for estimating the solvation effects [62,63,64,65], they may not fully describe the interactions of ORR adsorbates with water molecules due to the formation of H-bonds. As an alternative to implicit solvent approaches, explicit solvent models may provide a more comprehensive solution to describe the solvation effects on the ORR catalytic performance, which usually requires sampling thousands of solvent configurations, thus resulting in significant computational expense due to the use of classical molecular dynamic (MD) computations based on force fields. To this end, according to previous studies [66], explicit water layers were employed by placing 10 H2O molecules on the adsorbed oxygenated species (Figure S4). The results showed that the computed η value for ORR on the N@CrS2 catalyst using the explicit models is 0.50 V (Figure S4), which is close to our implicit solvent value (0.48 V). Thus, the implicit solvation model can also provide a reasonable estimation of the solvation energy for ORR intermediates, consistent with previous theoretical studies [67].

4. Discussion

Based on the well-accepted Sabatier principle, either too strong or too weak adsorption of reaction intermediates on catalysts can result in poor catalytic activity. This is because adsorption that is too strong will hamper the desorption process, resulting in poisoned catalysts, whereas too weak adsorption will induce insufficient activation of intermediates. Thus, the best catalysts exhibit the optimal adsorption strength with reaction intermediates, which locates at the peak of the volcano plot. Clearly, the ORR’s catalytic activity is intrinsically dependent on the adsorption strength of reaction intermediates with catalysts. To this end, we scaled the adsorption free energies of OOH* (ΔGOOH*), O* (ΔGO*) and OH* (ΔGOH*) on different X@CrS2 systems.

As shown in Figure 4a, obvious linear scaling relationships can be obtained between ΔGOOH* and ΔGOH* by ΔGOOH* = 0.72 ΔGOH* + 3.49 (R2 = 0.97) as well as ΔGO* and ΔGOH* by ΔGO* = 1.10 ΔGOH* + 0.72 (R2 = 0.89). Thus, ΔGOH* can be utilized as an eligible descriptor to describe the catalytic trend of these considered X@CrS2 catalysts. Furthermore, a volcano plot of ORR activity (ηORR) for X@CrS2 with the variation in ΔGOH* is obtained (Figure 4b), in which either a strong Si@CrS2 or a weak binding strength O@CrS2 with OH* species will induce poor ORR activity. On the contrary, the N@CrS2 catalyst displays a moderate binding strength with OH* and exhibits high ORR catalytic activity, making it locate at the peak of the volcano curve, and, thus, it becomes the best ORR catalyst among all the studied systems.

Figure 4

(a) Scaling relations between the free adsorption energies of intermediates (ΔGO* vs. ΔGOH* and ΔGOOH* vs. ΔGOH*) and (b) ηORR vs. ΔGOH* on different X@CrS2 systems.

To gain deep insight into the remarkable difference of OH* adsorption on X@CrS2, we turn to exploring the p-band center (εp) model of the non-metal active sites [68], in which Si@CrS2, N@CrS2, and O@CrS2 systems were chosen as the representatives of strong, moderate, and weak OH* adsorption. Notably, according to this model, the position of εp closer to the Fermi level will generally induce a stronger interaction of reaction species with catalysts. Our results demonstrated that the computed εp values of Si, N, and O dopants are −0.98, −2.37, and −3.30 eV, respectively, as shown in Figure 5. The moderate εp value on N@CrS2 suggests its optimal interaction with the oxygenated species, which is responsible for its superior catalytic performance.

Figure 5

The computed partial density of states for Si@CrS2, N@CrS2, and O@CrS2. The black and pink vertical dotted lines represent the Fermi energy level and p-band center (εp), respectively.

In addition, the charge density differences in Si@CrS2, N@CrS2, and O@CrS2 with adsorbed OH* species were computed (Figure 6). Upon OH* adsorption, we found that the charge depletion around the Si, N, and O endures, while the charge accumulation locates at the X–O bonds, indicating the charge transfer from the catalysts to oxygenated intermediates. According to the Bader charges analysis, the charge transfer is 0.71, 0.08, and 0.16 |e−| for OH* adsorption on Si@CrS2, N@CrS2, and O@CrS2, respectively, which is consistent with the binding strengths between them. Thus, the moderate adsorption strength of N@CrS2 endows it with its high ORR catalytic activity.

Figure 6

The charge difference density of OH* adsorption on (a) Si@CrS2N, (b) N@CrS2, and (c) O@CrS2 with the isosurface of 0.003 e/Å3.Cyane and yellow bubbles represent positive and negative charges, respectively.

5. Conclusions

In summary, by performing comprehensive DFT computations, we have systematically investigated the structures, stabilities, and magnetic and electronic properties as well as the ORR catalytic activity of several non-metal doped CrS2 monolayers. Our results demonstrated that most of the doped CrS2 materials generally possess high stability and hold great promise for experimental synthesis. As expected, depending on the kinds of the incorporated dopants, these CrS2 monolayers exhibit different magnetic and electronic properties. Based on the computed free energy changes in all elementary steps during ORR, the N@CrS2 was revealed as a quite promising ORR electrocatalyst due to its lower overpotential (0.41 V) than that of the Pt benchmark (0.45 V). Moreover, obvious scaling linear relationships between oxygenated species can be obtained, which was employed to construct the volcano curve between ORR catalytic activity and OH* binding strength. Understandably, the moderate p-band center and charge transfer from catalyst to oxygenated species render the N@CrS2 catalyst’s optimal interaction with reaction species, thus rationalizing its outstanding ORR catalytic performance. Our findings not only provide a novel strategy for the design of low-cost, highly-efficient ORR electrocatalysts, but also further widen the potential applications of CrS2-based materials.