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.


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, CO 2 , NO, and SO 2 , 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 Nanomaterials 2022, 12, 3012 2 of 11 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 MoS 2 possesses high catalytic activity for oxygen evolution reaction [29]. Moreover, Lv et al. proposed that N-doped MoS 2 displayed a high faradaic efficiency and low onset overpotential for CO 2 electroreduction to CO [30]. Theoretically, Li et al. demonstrated that B-doping causes the VS 2 monolayer to exhibit outstanding catalytic activity towards nitrogen reduction reactions [31]. Moreover, Singh et al. reported that the ORR activity of N-doped WS 2 monolayers should originate from the introduction of spin density caused by N doping [32]. In addition, Tian et al. showed that Co doped 1T-TiS 2 exhibits significantly enhanced performance toward the ORR [33].
As a representative of 2D TMDs, the chromium disulfide (CrS 2 ) monolayer was synthesized via the chemical vapor deposition (CVD) method in 2019 [21]. Interestingly, due to its unique properties, the CrS 2 monolayer has wide applications in spintronic devices [36,37]. For example, Chen et al. suggested that CrS 2 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 CrC 2 monolayer can be effectively tuned by doping metal atoms [37].
Inspired by the interesting properties of the CrS 2 monolayer, in this work, we explored the potential of several non-metal-doped CrS 2 monolayers (X@CrS 2 , 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@CrS 2 systems, we suggested that these doped CrS 2 candidates are very likely to be synthesized in experiments and possess extremely high stability. Moreover, the N@CrS 2 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.

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@CrS 2 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 + ∆G U , 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. ∆G U = −eU, where U is the applied potential. Since DFT generally fails to obtain the energy of the O 2 molecule, the free energy of O 2 (G O2 ) will be computed via the energies of H 2 O and H 2 , namely G O2 = G H2O − 2G H2 + 4.92 eV. Furthermore, the ORR's catalytic activity was evaluated by computing the corresponding overpotential (η) according to the following equations: η = max{∆G 1 , ∆G 2 , ∆G 3 , ∆G 4 }/e + U 0 , 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 U 0 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.

Structures, Stabilities, and Properties of X@CrS 2 Catalysts
First, we investigated the structures and electronic properties of the pristine CrS 2 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 CrS 2 monolayer is 3.04 Å, while the lengths of the formed Cr-S bonds are 2.28 Å. Moreover, the CrS 2 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 CrS 2 monolayer are consistent with previous theoretical reports [50], indicating the accuracy of our computational methods to describe the behavior of the CrS 2 monolayer. 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 + U 0 , 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 U 0 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.

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. 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 Based on the optimized CrS 2 monolayer, various non-metal atoms were introduced to substitute one of the S atoms within the CrS 2 monolayer to construct the doped CrS 2 system. After fully geometrical relaxation, the structures of the X@CrS 2 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 CrS 2 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 CrS 2 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. 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.  To assess the experimental feasibility of these doped CrS2 monolayers, we computed their formation energies (Ef) based on the following equation: [51,52], where the EX@CrS2 and ECrS2 are the total energies of doped and pristine CrS2 monolayers, respectively. μX and μS are the chemical potentials of a doping atom and S atom. Notably, μX can be derived from the energy per atom in their bulk X materials. In contrast, μS 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 μS should meet the relationship from the thermodynamic equilibrium condition: μCrS2 = μCr + 2μS, 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, μS can be obtained by: μS = (μCrS2 − μCr)/2. On the other hand, under S-rich conditions, the μS 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,  To assess the experimental feasibility of these doped CrS 2 monolayers, we computed their formation energies (E f ) based on the following equation: [51,52], where the E X@CrS2 and E CrS2 are the total energies of doped and pristine CrS 2 monolayers, respectively. µ X and µ S are the chemical potentials of a doping atom and S atom. Notably, µ X can be derived from the energy per atom in their bulk X materials. In contrast, µ S 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 µ S should meet the relationship from the thermodynamic equilibrium condition: µ CrS2 = µ Cr + 2µ S , where µ CrS2 represents the total energy of the pristine CrS 2 monolayer per formula unit. Under Cr-rich conditions, the µ Cr was taken from the bulk Cr. Thus, µ S can be obtained by: µ S = (µ CrS2 − µ Cr )/2. On the other hand, under S-rich conditions, the µ S can be obtained from its bulk S 8 state, and then µ Cr can be determined by: µ Cr = µ CrS2 − 2µ S . The computed formation energies of these X@CrS 2 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 CrS 2 monolayer is difficult. On the contrary, when O, Si, Cl, As, Se, and Br dopants were introduced into the CrS 2 monolayer, their E f values were always negative, implying notable feasibility for their experimental synthesis. Notably, the CrS 2 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 CrS 2 monolayers can be easily fabricated in experiments, except for B-doping.
To further explore the structural stability of X@CrS 2 , we computed the binding energies (E bind ) of the introduced dopants on the CrS 2 monolayer according to the following equation: , and E X represent the total electronic energies of doped CrS 2 monolayers, the defective CrS 2 substrate, and isolated X atoms in their most stable phase. It can be observed that the E bind values of these dopants on the CrS 2 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@CrS 2 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.  [55]. Therefore, we strongly believe that the as-designed doped CrS 2 monolayer holds great promise for synthesis.
Next, we explored the magnetic and electronic properties of these X@CrS 2 systems. Our results showed that the pristine CrS 2 monolayer is a nonmagnetic material. After introducing these dopants into the CrS 2 monolayer, however, different magnetic behaviors can be observed: (1) C-, O-, and Se-doped CrS 2 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 µ B , while the Si-doped CrS 2 monolayer has a magnetic moment of 2.00 µ B . Moreover, we noted that these dopants will loot different amounts of electrons (0.04~0.95 |e − |) from the CrS 2 monolayer, except for the Si dopant, which can donate about 0.32 electrons to the CrS 2 substrate. As a result, the band gaps of the doped CrS 2 monolayer are reduced to different degrees due to the introduction of impurity levels. For example, the N-doped CrS 2 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 MoS 2 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 CrS 2 monolayers, which consist of 25 Cr and 50 S atoms. We will mainly focus on the trend of the band gaps of the CrS 2 monolayer after non-metal doping.

ORR Catalytic Activity
After knowing that these non-metal doped CrS 2 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 O 2 molecule can be reduced to H 2 O along a four-electron (4e − ) pathway ( Figure 3a): (1) * + O 2 (g) + H + + e − 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 (Ebind) of the introduced dopants on the CrS2 monolayer according to the following equation: Ebind = 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 Ebind 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 μB, while the Si-doped CrS2 monolayer has a magnetic moment of 2.00 μB. 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.

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 OOH * ; (2) OOH * + H + + e − As, Se, and Br dopants were introduced into the CrS2 monolayer, their Ef values always negative, implying notable feasibility for their experimental synthesis. No the CrS2 monolayers doped by C, N, and P atoms are more likely to form under C conditions due to their negative formation energies. Obviously, by carefully tunin reaction conditions, these doped CrS2 monolayers can be easily fabricated in experim except for B-doping.
To further explore the structural stability of X@CrS2, we computed the binding gies (Ebind) of the introduced dopants on the CrS2 monolayer according to the follo equation: Ebind = EX@CrS2 − ECrS2 − EX, where EX@CrS2, ECrS2, and EX represent the total elec energies of doped CrS2 monolayers, the defective CrS2 substrate, and isolated X ato their most stable phase. It can be observed that the Ebind values of these dopants o CrS2 substrate range from −2.52 eV of Br to −7.16 eV of C, indicative of the strong bin strength, thus ensuring their high stability. The stability of X@CrS2 was further eval by using AIMD simulations, where N doping was taken as a representative. As s from Figure S1, there is still no significant distortion of the geometric structure at 5 indicative of its high thermodynamic stability. It should be noted that incorporati various non-metal atoms into 2D TMDs has been already realized experimentally. F ample, Li et al. reported a simple, facile, and effective strategy to fabricate an N-d MoS2 nanosheet [54], while Xin et al. demonstrated the fabrication of a P-doped nanosheet by the one-step hydrothermal method [55]. Therefore, we strongly believ the as-designed doped CrS2 monolayer holds great promise for synthesis.
Next, we explored the magnetic and electronic properties of these X@CrS2 sys Our results showed that the pristine CrS2 monolayer is a nonmagnetic material. Aft troducing these dopants into the CrS2 monolayer, however, different magnetic beha can be observed: (1) C-, O-, and Se-doped CrS2 systems still retain their nonmagnetic due to the absence of unpaired electrons; (2) as for B, N, P, As, or Br-doped system total magnetic moment is close to 1.00 μB, while the Si-doped CrS2 monolayer has a netic moment of 2.00 μB.
Moreover, we noted that these dopants will loot diff amounts of electrons (0.04~0.95 |e − |) from the CrS2 monolayer, except for the Si do which can donate about 0.32 electrons to the CrS2 substrate. As a result, the band ga the doped CrS2 monolayer are reduced to different degrees due to the introduction o purity levels. For example, the N-doped CrS2 system exhibits a smaller band gap o eV than that of the pristine one (0.99 eV), suggesting the enhanced electrical conduct which may facilitate the rapid charge transfer in electrocatalysis [56]. It is noted th GGA method usually underestimates band gaps, whereas the hybrid functional me such as Heyd-Scuseria-Ernzerhof (HSE06), can rectify the band gap [57]. For exampl computed band gaps of the well-established MoS2 monolayer using HSE06 and PBE m ods are 2.20 and 1.79 eV [58], respectively. Despite the fact that a more accurate ban value can be predicted by means of the HSE06 method, the corresponding computa cost is extremely high for the 10 doped CrS2 monolayers, which consist of 25 Cr and atoms. We will mainly focus on the trend of the band gaps of the CrS2 monolayer non-metal doping.

ORR Catalytic Activity
After knowing that these non-metal doped CrS2 candidates exhibit great potent experimental synthesis, high stability, diverse magnetic moment, and enhanced co tivity, we further explored their catalytic activities towards ORR.
Typically, the O2 molecule can be reduced to H2O along a four-electron (4e − ) pat ( and Br dopants were introduced into the CrS2 monolayer, their Ef values were negative, implying notable feasibility for their experimental synthesis. Notably, 2 monolayers doped by C, N, and P atoms are more likely to form under Cr-rich ns due to their negative formation energies. Obviously, by carefully tuning the conditions, these doped CrS2 monolayers can be easily fabricated in experiments, or B-doping. further explore the structural stability of X@CrS2, we computed the binding enerind) of the introduced dopants on the CrS2 monolayer according to the following n: Ebind = EX@CrS2 − ECrS2 − EX, where EX@CrS2, ECrS2, and EX represent the total electronic s of doped CrS2 monolayers, the defective CrS2 substrate, and isolated X atoms in ost stable phase. It can be observed that the Ebind values of these dopants on the bstrate range from −2.52 eV of Br to −7.16 eV of C, indicative of the strong binding , thus ensuring their high stability. The stability of X@CrS2 was further evaluated g AIMD simulations, where N doping was taken as a representative. As shown gure S1, there is still no significant distortion of the geometric structure at 500 K, ve of its high thermodynamic stability. It should be noted that incorporation of non-metal atoms into 2D TMDs has been already realized experimentally. For ex-Li et al. reported a simple, facile, and effective strategy to fabricate an N-doped anosheet [54], while Xin et al. demonstrated the fabrication of a P-doped MoS2 eet by the one-step hydrothermal method [55]. Therefore, we strongly believe that esigned doped CrS2 monolayer holds great promise for synthesis. xt, we explored the magnetic and electronic properties of these X@CrS2 systems. ults showed that the pristine CrS2 monolayer is a nonmagnetic material. After inng these dopants into the CrS2 monolayer, however, different magnetic behaviors bserved: (1) C-, O-, and Se-doped CrS2 systems still retain their nonmagnetic states he absence of unpaired electrons; (2) as for B, N, P, As, or Br-doped systems, the gnetic moment is close to 1.00 μB, while the Si-doped CrS2 monolayer has a magoment of 2.00 μB.
Moreover, we noted that these dopants will loot different s of electrons (0.04~0.95 |e − |) from the CrS2 monolayer, except for the Si dopant, an donate about 0.32 electrons to the CrS2 substrate. As a result, the band gaps of ed CrS2 monolayer are reduced to different degrees due to the introduction of imevels. For example, the N-doped CrS2 system exhibits a smaller band gap of 0.25 that of the pristine one (0.99 eV), suggesting the enhanced electrical conductivity, ay facilitate the rapid charge transfer in electrocatalysis [56]. It is noted that the ethod usually underestimates band gaps, whereas the hybrid functional method, Heyd-Scuseria-Ernzerhof (HSE06), can rectify the band gap [57]. For example, the ed band gaps of the well-established MoS2 monolayer using HSE06 and PBE meth-2.20 and 1.79 eV [58], respectively. Despite the fact that a more accurate band gap n be predicted by means of the HSE06 method, the corresponding computational xtremely high for the 10 doped CrS2 monolayers, which consist of 25 Cr and 50 S We will mainly focus on the trend of the band gaps of the CrS2 monolayer after tal doping. 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.

R Catalytic Activity
To further explore the structural stability of X@CrS2, we computed the binding energies (Ebind) of the introduced dopants on the CrS2 monolayer according to the following equation: Ebind = 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 Ebind 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 μB, while the Si-doped CrS2 monolayer has a magnetic moment of 2.00 μB.
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.

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) [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 [59,60]. We again took the N@CrS 2 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@CrS 2 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 O 2 molecule (1.23 Å). Subsequently, the approach of a second hydrogen induces the dissociation of OOH * into (O * + H 2 O) or the formation of H 2 O 2 . Remarkably, the OOH * formation is exothermic by 2.43 eV, which is much larger than that of H 2 O 2 formation (0.16 eV), suggesting that N@CrS 2 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 H 2 O 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. Nanomaterials 2022, 12, 3012 6 of 11 (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. 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]. As for other X@CrS 2 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@CrS 2 (0.97 V) < C@CrS 2 (1.20 V) < O@CrS 2 (1.53 V) < Se@CrS 2 (1.61 V) < Br@CrS 2 (1.74 V) ≈ Cl@CrS 2 (1.78 V) < P@CrS 2 (2.18 V) < B@CrS 2 (2.42 V) < Si@CrS 2 (2.82 V), as shown in Figure S2, which are all higher than that of N@CrS 2 (0.41 V). In particular, the overpotential of N@VS 2 (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@CrS 2 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@CrS 2 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 H 2 O molecules on the adsorbed oxygenated species ( Figure S4). The results showed that the computed η value for ORR on the N@CrS 2 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].

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 * (∆G OOH* ), O * (∆G O* ) and OH * (∆G OH* ) on different X@CrS 2 systems.
As shown in Figure 4a, obvious linear scaling relationships can be obtained between ∆G OOH* and ∆G OH* by ∆G OOH* = 0.72 ∆G OH* + 3.49 (R 2 = 0.97) as well as ∆G O* and ∆G OH* by ∆G O* = 1.10 ∆G OH* + 0.72 (R 2 = 0.89). Thus, ∆G OH* can be utilized as an eligible descriptor to describe the catalytic trend of these considered X@CrS 2 catalysts. Furthermore, a volcano plot of ORR activity (η ORR ) for X@CrS 2 with the variation in ∆G OH* is obtained (Figure 4b), in which either a strong Si@CrS 2 or a weak binding strength O@CrS 2 with OH * species will induce poor ORR activity. On the contrary, the N@CrS 2 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.
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 (R 2 = 0.97) as well as ΔGO* and ΔGOH* by ΔGO* = 1.10 ΔGOH* + 0.72 (R 2 = 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. 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. To gain deep insight into the remarkable difference of OH * adsorption on X@CrS 2 , we turn to exploring the p-band center (ε p ) model of the non-metal active sites [68], in which Si@CrS 2 , N@CrS 2 , and O@CrS 2 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@CrS 2 suggests its optimal interaction with the oxygenated species, which is responsible for its superior catalytic performance.
In addition, the charge density differences in Si@CrS 2 , N@CrS 2 , and O@CrS 2 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@CrS 2 , N@CrS 2 , and O@CrS 2 , respectively, which is consistent with the binding strengths between them. Thus, the moderate adsorption strength of N@CrS 2 endows it with its high ORR catalytic activity. 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.

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  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.

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

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 CrS 2 monolayers. Our results demonstrated that most of the doped CrS 2 materials generally possess high stability and hold great promise for experimental synthesis. As expected, depending on the kinds of the incorporated dopants, these CrS 2 monolayers exhibit different magnetic and electronic properties. Based on the computed free energy changes in all elementary steps during ORR, the N@CrS 2 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@CrS 2 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 CrS 2 -based materials.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/nano12173012/s1, Figure S1: Total potential temperature and energy fluctuations after equilibration for N@CrS 2 in the AIMD simulation; Figure S2: The computed free energy diagrams for ORR on different X@CrS 2 systems; Figure S3: The computed free energy diagrams for ORR on N@CrS 2 with solvent effect. Figure S4: The computed free energy diagrams and the corresponding intermediates configurations for ORR on N@CrS2 with the explicit solvent models; Tables S1: The computed the overpotentials (η ORR , V) at various energy cutoff (eV) or k-points; Tables S2: The computed the overpotential (η ORR ) and the free adsorption energies (∆G, eV) for various X@CrS2 materials.