Origin of the catalytic activity of phosphorus doped MoS2 for oxygen reduction reaction (ORR) in alkaline solution: a theoretical study

Phosphorus doped MoS2 nanosheets (P-doped MoS2) have been reported as excellent oxygen reduction reaction (ORR) catalysts with four-electron selectivity in alkaline solution. By performing density functional theory (DFT) calculations, we revealed the detailed reaction mechanism and the key reaction sites on surface of P-doped MoS2 for ORR catalysis. The double P-doped MoS2 (2P-MoS2) is calculated to be more stable than the single P-doped MoS2 (P-MoS2), and the configuration with two P atoms in neighboring sites exhibits the highest stability. The surface of P-doped MoS2 is found highly active for dissociation of O2. Comparative calculations reveal that P-MoS2 is unsuitable as ORR catalyst due to the high dissociation barrier of H2O (1.19 and 2.06 eV for the first and second adsorbed H2O), while the 2P-MoS2 shows good ORR catalytic activity with much lower dissociation barrier of H2O (0.62 eV). Furthermore, we elucidated that the ORR catalytic activity in 2P-MoS2 originates from the activated S2 atom, which provides an extra adsorption site for the first H2O and the following OH group benefited from the enhanced hydrogen bond interaction. Our results illustrate the mechanisms of doped MoS2 based catalysts and provide rational way for designing ORR catalysts with high activity.

experimentally. The experimental data indicate that N-and P-doped single layered MoS 2 possess wonderful ORR catalytic performance in acidic and alkaline environment, respectively 36,37 . Soon after, the mechanism about the ORR catalytic effect on N-and P-doped MoS 2 single layer is studied by simulations 38 . However, the catalytic mechanism on the experimental results of P-doped 2D materials as ORR catalysts is not studied, especially the reported generation of extra ORR activity on P-doped MoS 2 nanosheet with four-electron selectivity in alkaline solution, which have been successfully synthesized via the method of pyrolysis 37 . More importantly, Huang and the co-worker have found through the mapping analysis that, P doping on the MoS 2 monolayer can be dispersed on the surface evenly, and with the increase of P-doping amount, the activity of P-doped MoS 2 on ORR can be improved obviously, which reveals that the active site should be related with the P atoms.
In this work, by density functional theory (DFT) calculations, we explored the origin of catalytic activity in the P-doped MoS 2 monolayer as the ORR electrocatalyst. Previous theoretical and experimental studies 37 have indicated that, the P-doped MoS 2 monolayer prefers to work as ORR catalyst in alkaline solution rather than acid solution. According to our simulation results, the P-doped MoS 2 with doping concentration of 5.5% is optimal for ORR catalytic activity (the optimal value observed experimentally is 4.7% 37 ), since it results in stable adsorption of H 2 O molecule and the extra OH group through the activated S site, and facilitates the following ORR steps to proceed with a lower barrier in alkaline media. Our work reveals the origin of catalytic activity of P-doped MoS 2 in catalyzing ORR in alkaline solution and provides rational proposal for designing novel ORR catalysts.

Results and Discussion
Characterization of P-doped MoS 2 . Before the researches on ORR catalytic performance of P-doped MoS 2 nanosheets, we first characterized their structural stability, structural and electronic properties, and pristine MoS 2 3 × 3 supercell is taken as comparison. The stability of pristine MoS 2 and P-doped MoS 2 are estimated by computing their formation energies (E f ) as follows, where E total is the total energy of the nanosheet, n Mo , n S , and n P are the numbers of Mo, S, and P atoms in the cell. N is the total number of the atoms. μ Mo , μ S , and μ p are the chemical potentials of Mo, S, and P, respectively. μ P is obtained from black phosphorus, while μ Mo is obtained from bulk molybdenum and μ S is obtained from bulk MoS 2 by subtracting the chemical potential of Mo, respectively. According to the definition, the formation energies of MoS 2 , P-MoS 2 and three types of 2P-MoS 2 are calculated as −0.052, −0.070, −0.103, −0.105, and −0.105 eV/atom, respectively (as listed in Table 1). Compared with the pristine MoS 2 and P-MoS 2 , 2P-MoS 2 monolayers exhibit very close formation energy (around 0.1 eV/atom), indicating that the P-doping concentration can be increased from 3.7 atom% to 5.5 atom% easily on pristine MoS 2 monolayer. The bond energy is another parameter to be compared. The bond energy is calculated as followed: where E total is the total energy of nanosheet, n is the number of P atom, the E P and E D denote the energy of isolated P atom and the energy of MoS 2 sheet with vacancy defect (VS), respectively. As listed in Table 1, all the 2P-doped ones have as twice as lager Mo-P bond energy than single P-doped one, indicating the significantly higher stability of Mo-P bond in 2P-MoS 2 than that in P-MoS 2 . As a comparison, the bond energy of Mo-S bond in pristine MoS 2 is 0.224 eV, which is only slightly larger than those of Mo-P bond in 2P-MoS 2 . Attributing to the neighboring position of P and S in periodic table, the doping of P doesn't affect the initial structure of MoS 2 monolayer much, which can be confirmed by the few relative height of P in doped MoS 2 . By comparing the internal energy of three different 2P-doping MoS 2 monolayers through different softwares (see Table S1), we confirm that the doped P atoms in MoS 2 tend to distribute locally (Fig. 1c) other than separately. Thus the first type of 2P-doping MoS 2 should be the most stable configuration, which is the most possible to be experimentally prepared. This could be further validated by the tendency study. As shown in Fig. S3, we evaluated the tendency of increasing P doping concentration against the energy difference between locally distributed one and separately distributed one. A significant tendency is revealed that with the doping concentration of P increasing, the energy difference between locally distributed one and separately distributed one also increases. Thus it is not accidental that 2P-MoS 2 (1) shows the highest stability among the three configurations, and P in P-doped MoS 2 prefers to distribute locally.
Then the electronic properties of P-doped MoS 2 are evaluated. As displayed in Fig. 1, we plotted the band structures of the P-MoS 2 and 2P-MoS 2 , along with the pristine MoS 2 as comparison. Our DFT results indicate that the doping of P leads to vanishing of energy gap [zero gap for P-MoS 2 , 2P-MoS 2 (2) and 2P-MoS 2 (3), and . This is mostly due to the introduction of impurity band, which exists as a single band across or above the Fermi level that is splitting from valence band edge. Further analysis from partial density of states (PDOS) tells that the splitting of impurity band origins from the orbital hybridization between Mo and P, as shown from PDOS plots in Fig. 1. The small or zero gap induced by doping of P not only greatly enhances the conductivity of MoS 2 , but also facilitates the chemical reactivity. Additionally, the difference of electronegativity between S and P atom also plays important role on the charge population. As shown in Fig. 1, the S atoms (S2-S5) surrounding the doped P atom carry more negative charge (about −0.11 |e|) than that in pristine MoS 2 (S1, −0.099 |e|). Meanwhile P atoms carry positive charges (more than +0.08 |e|), which renders the P atom as suitable adsorption site for the nucleophilic species. The adsorption of O 2 on P-doped MoS 2 . As the very beginning of ORR, adsorption of O 2 is taken into consideration first. By putting O 2 molecule on different sites on P-MoS 2 and comparing their E ads , we determined the most energetically favorable adsorption configuration as shown in Fig. 2a. As we can see, the O 2 molecule chemisorbs to the P-MoS 2 by bonding with P atom. The bond length of P-O bonds is about 1.670 Å. And the calculated adsorption energy of O 2 on P-MoS 2 is up to E ads = −1.11 eV, which is much more favorable than that on pristine MoS 2 monolayer (−0.11 eV). And the charge transferred from P-MoS 2 to O 2 is as much as Q = 0.695 |e| (in ref. 33 , E ads = −0.93 eV and Q = 0.30 |e|). Further results from TS search show that the chemisorbed O 2 can be easily dissociated into two O atoms by overcoming a negligible energy barrier (about 0.05 eV). As shown in Fig. 2b, the two dissociated O atoms adsorb on the top site of P atom and the bridge site between P and Mo atoms, respectively. Similarly, the adsorption of other ORR species (OH and H 2 O) also prefers to adsorb onto the P atom, as displayed in Fig. S1 of Supporting Information.
There are two different occasions for adsorption of O 2 molecule on 2P-MoS 2 . For 2P-MoS 2 (1), the O 2 molecule binds with the two neighboring P sites with bridged configuration. For 2P-MoS 2 (2) and (3), it is found that due to the large distance between the two P atoms, O 2 molecule can only separately bind with one P atom as the situation of P-MoS 2 (Fig. S1). In the following discussion, we only investigated the adsorption performance of O 2 on 2P-MoS 2 (1), not only due to its highest stability among the 2P-MoS 2 , but also for its different binding configuration with O 2 compared with P-MoS 2 . Resulting from the high chemical activity of the two adjacent P atoms, the O 2 molecule will directly dissociated and form two P=O bonds without any energy barrier, as shown in Fig. 2c. The bond length of the two P=O bonds are 1.505 Å, which is shorter than the P=O bond (1.515 Å) formed on P-MoS 2 .
The charge population of P-MoS 2 and 2P-MoS 2 (1) after dissociation of O 2 is also compared. It is found that comparing with the S atoms (S1 and S2) in P-MoS 2 , the S atoms surrounding dopants (S1, S2, and S3 in Fig. 2c) on 2P-MoS 2 carry less negative electron, which may be favorable for the ORR steps as we will discuss in the following parts.
ORR on P-MoS 2 . Generally, the total reaction in alkaline solution can be expressed as a 4-electron evolution (*denotes the adsorbed surface): The elementary reaction steps could be expressed as:  Fig. 3b. The activation barrier for this reaction is up to 2.06 eV, even higher than the dissociation of first H 2 O. Moreover, this reaction is endothermic with reaction energy of 0.73 eV, which is thermodynamically unstable. The high activation barrier for dissociation of H 2 O demonstrates that the catalytic activity of P-MoS 2 is low for ORR, which mainly due to the strong P=O bond and low reaction activity of H 2 O on P-MoS 2 . This is in consistent with the previous theoretical investigations 38 .
ORR on 2P-MoS 2 . As is presumed from the charge population and electrostatic potential (ESP, shown in Fig. S2) distribution, the two adjacent doping P atoms on the 2P-MoS 2 (1) should be the most active sites for ORR.  Table 2. It is found that the decrease of charge assigned on S2 (Fig. 2c) accounts for the improved adsorption performance of H 2 O. As we can see in Fig. 2c, due to the strong electronegativity of O, the presence of dissociated O atoms will decrease the negative charge of the surface S atoms. Interestingly, the S2 in 2O adsorbed 2P-MoS 2 (1) exhibits much less negative charge (−0.068 |e| as marked in Fig. 2c) other than that in 2 O adsorbed P-MoS 2 (−0.098 |e| as marked in Fig. 2b). It is the charge decrease as well as stronger hydrogen bond caused by the two nearby P-O bonds [as depicted in structure (2) of Fig. 4] together enhance the adsorption performance of H 2 O on the S2 site. Combined with the fact that the adsorption of H 2 O is strengthened, we assumed that it is this S2 atom that makes 2P-MoS 2 (1) different, and further simulations confirm our proposal. We found that this S2 site can also serve as an adsorption site for OH group, which makes it active during the reaction. As the structure (2) displayed in Fig. 4, accompanied with the hydrogen bond interaction from two P=O bonds, the first added H 2 O molecule could be easily dissociated into two adsorbed OH groups, one binds with the P and the other one adsorbs onto the S2 site. The corresponding activation barrier for the dissociation of first H 2 O is 0.62 eV, a moderate value for reaction under room temperature. As comparison, both the two dissociated OH groups in P-MoS 2 can only adsorb onto the P atom, leading to significant structural distortion and a high barrier of 1.19 eV.
Following that, the formed OH groups would combine with the electron provided from the electrode and solve into the solution as OH−. And then, the exposed P atom will be the most active site on the surface of 2P-MoS 2 (1). Both of the strong reaction activity of P atom and the negative charge on O atom contribute to the dissociation of the second H 2 O without energy barrier. Also, the possibility of the breaking of P=O bond and the   formation of OH on P atom is strongly related with the assist of the hydrogen-bond. The whole reaction circulation and the corresponding energy profile of reaction pathway are shown in Fig. 4.

Comparison between P-MoS 2 & 2P-MoS 2 .
Comparing 2P-MoS 2 with P-MoS 2 , the most significant difference is the activation of the S nearby the doping P atoms (as shown in Fig. 2c) and the strong hydrogen-bond existed in 2O-2P-MoS 2 , which leads to the enhanced adsorption for H 2 O and the chemisorption of OH on S atom during the dissociation of H 2 O molecule. Resulting from the more stable adsorption of H 2 O as well as the OH group on the activated S site, the reaction energy barrier for the dissociation of the first H 2 O is greatly lowered compared with the P-MoS 2 (0.62 eV for 2P-MoS 2 and 1.19 eV for P-MoS 2 ). Associated with the fact that the formation energy for 2P-MoS 2 is almost the same with the P-MoS 2 (Table 1), it is much easily for 2P-MoS 2 to be prepared experimentally. Meanwhile, among the 2P-MoS 2 nanosheets, type 1 is the most energetically favorable (Table S1). Based on the comparison results, we confirm that the activated S atom in 2P-MoS 2 is the origin of ORR catalytic activity of P-doped MoS 2 in alkaline solution. Notably, it is not accident that the 2P-MoS 2 we modeled shows improved ORR catalytic performance. As we have stated, the doping concentration of 2P-MoS 2 is 5.5%, which is quite close to the optimal doping concentration observed experimentally (4.7%) 37 . It is the special structure of 2P-MoS 2 that activates the S2 site and facilitates the adsorption of H 2 O molecule and the extra OH group, which well explains why the MoS 2 with a P doping concentration of 4.7% exhibits the best ORR performance in experiments.

Conclusions
In summary, we performed comprehensive studies on the origin of the ORR catalytic activity of P-doped MoS 2 in alkaline solution through DFT calculations. The 2P-MoS 2 with doping concentration of 5.5% is calculated to be more stable than the P-MoS 2 with doping concentration of 3.7%, which is consistent with the optimal doping concentration of 4.7% from the experiments 37 . And the configuration with two P atoms in neighboring sites is found to possess the highest stability among the possible configurations of 2P-MoS 2 . The surface of P-doped MoS 2 exhibits highly reactive activation for dissociation of O 2 . Except the doping on MoS 2 plane, we also discuss the influence of phosphorus present on edge of MoS 2 . The detailed structure configuration and energy values have been given in Fig. S6. For the two edge S sites that can be placed by P (Mo edge and S edge as shown in Fig. S6 (a)), the P prefer doping in S edge site rather than Mo edge site considering the relative energy difference up to 0.97 eV. For the O 2 adsorption on the P-doped MoS 2 on the edge site, two adsorption sites are considered (shown in Fig. S6 (b)). The positive adsorption energy on the S site indicated that P-doped on the edge site can't activate the S site around P atom. While on the edge P site, the adsorbed O 2 would dissociated directly, which is same with single P atom doped on MoS 2 plane (Fig. 2 (b)). And the adsorption for H 2 O is −0.13 eV, which is also the same with the P-MoS 2 (−0.129 eV, as shown in Table 2). Thus, based on our calculation results, we concluded that the doping of P in the basal plane and edge sites will exhibit similar adsorption performance for adsorption of H 2 O and dissociation for O 2 .
Further comparative calculations reveal that P-MoS 2 is unsuitable as ORR catalyst due to the strong P=O bond and high dissociation barrier of H 2 O on P-MoS 2 (1.19 eV and 2.06 eV for the first and second adsorbed H 2 O). In comparison, the 2P-MoS 2 shows good ORR catalytic activity. The dissociation barrier of the first H 2 O on 2P-MoS 2 (0.62 eV) is much lower than that of P-MoS 2 . And the dissociation of O 2 and the second H 2 O on 2P-MoS 2 is even spontaneous without any activation barriers. While for the two other types of 2P-doped MoS 2 , not only the stability is less than the neighboring doped 2P-MoS 2 , but also the reaction mechanism is similar with the single P doped MoS 2 . By checking the charge population and adsorption performance of the ORR intermediates, we elucidated that the ORR catalytic activity in 2P-MoS 2 originates from the activated S2 atom with increased charge that affected by the two neighboring P=O bonds, which provides an extra adsorption site for the first H 2 O and the following OH group with the help from the enhanced hydrogen bond interaction. Our results reveal the detailed reaction mechanism of MoS 2 based catalysts and provides instructive suggestions for designing ORR catalysts with high performance.

Methods
All the calculations are carried out by DMol 3 program. The electronic exchange and correlation effects are described by the generalized gradient approximation with Perdew-Burke-Ernzerhof (GGA-PBE) functional 39 . In order to accurately describe the van der Waals forces, dispersion correction (DFT-D) is adopted by Grimme approach 40 . The all-electron double numerical atomic orbital including polarized p-function (DNP) 41 is chosen as the basis set with orbital cutoff of 4.9 Å. The convergence threshold values for energies, gradient and displacement are specified as 1 × 10 −5 Ha, 1 × 10 −3 Ha/Å, and 5 × 10 −3 Å, respectively, while the self-consistent-field (SCF) convergence threshold value is 1 × 10 −6 Ha. To enhance SCF convergence efficiency during optimization, an electron thermal smearing value of 0.002 Ha is employed for all the calculations. The conductor-like screening model (COSMO) 42 is used to simulate the aqueous environment, where the dielectric constant is set as 78.54 (water). All the stably existed species during the ORR process are confirmed no imaginary frequencies by performing frequency analysis. To determine the activation barriers (E b ), transition state (TS) searches are conducted by using the complete linear synchronous transit/quadratic synchronous transit (LST/QST) method. All the obtained TS structures exhibit only one imaginary frequency throughout the potential surface, and are ensured to directly connect corresponding reactants and products by nudgedelastic band (NEB) algorithm 43 .
The unit cell of 2H-MoS 2 is optimized by using 9 × 9 × 1 Monkhorst-Pack k-point grids. The single P doped MoS 2 is then obtained by substituting a surface S atom in MoS 2 3 × 3 supercell into P atom, with a doping concentration of 3.7 atom%. The structures of MoS 2 3 × 3 supercell and single P-doped MoS 2 nanosheet (P-MoS 2 ) are depicted in Fig. 1(a,b). The double P-doped MoS 2 nanosheet (2P-MoS 2 ) is obtained by substituting two surface Scientific REPORTS | (2018) 8:13292 | DOI:10.1038/s41598-018-31354-0 S atoms in MoS 2 2√3 × 2√3 supercell into P atoms, with a doping concentration of 5.5 atom%. Differing from the relative positions between the two doped P atoms, there are three possible configurations that are taken into consideration. That is, neighboring sites (type 1), spaced by a Mo atom (type 2), and spaced by a S atom (type 3), as shown in Fig. 1(c-e). Considering the close atomic volume between P and S atom, the influence of doping on lattice size is neglected. Brillouin zone of MoS 2 3 × 3 supercell and the P-doped MoS 2 structures are sampled by 3 × 3 × 1 Monkhorst-Pack k-point grids during the geometry optimization, which is tested converging for total energy. Population analysis of electron is performed by assigning Hirshfeld charges 44 for the optimized structures.
In this article, the adsorption energy (E ads ) is defined as ads m ol sub m ol sub / Where E mol and E sub represent the total energies of the isolated adsorbate and adsorption substrate, and Emol/ sub represents the total energy of the adsorption system. Here, the negative Eads means exothermic process for adsorption. The more negative Eads is, the more stable the adsorption system is. The triplet is taken as the ground state of O2 during the calculation.
In order to calculate the Gibbs free energy of reactions involving electron/proton transfer, the computational hydrogen electrode (CHE) model introduced by Nørskov and co-workers was used 45 . The reaction free energies (ΔG) can be computed as follows: where the ΔE was the reaction total energy directly obtained from Dmol3 calculation, ΔZPE is the change in the zero-point energy, T is the temperature (298.15 K), and ΔS is the change in entropy. ΔG pH = K b T × ln10 × pH represents the free energy contribution due to the variation in the H concentration, and the pH value in this work was assume to 14 in the alkaline medium. ΔG U = −neU, where n was the number of electrons transferred and U was the applied electrode potential. The entropy values are taken from the physical chemistry table considering H 2 and H 2 O in gaseous form at room temperature and atmospheric pressure whereas the entropy of the adsorbed state is considered to be negligible. ZPE of the free molecules are estimated from our DFT calculations considering vibrational frequencies of the molecules in the harmonic approximation. The details of calculation steps to plot Free energy profile are demonstrated in the supporting information.