Spin-polarized first-principle calculations were implemented by Vienna ab initio simulation package (VASP) with an energy cutoff of 500 eV, the Perdew-Burke-Ernzerhof (PBE) functional within the generalized gradient approximation (GGA) was adopted to address the exchange-correlation interactions . The P and TM doped 5×5×1 graphenes were modeled with a vacuum thickness of 15Å and optimized with a Γ centered 2×2×1 k-point sampling, specially, the 5×5×1 k-point was employed to calculate the density of states (DOS). During the optimized process, all atoms were relaxed until the convergence criterions of energy (1×10-6 eV) and force (0.02 eV Å -1 ) are reached. DFT-D3 developed by Grimme was utilized for the van der Waals correction . During the DFT calculations, the implicit solvent model was employed to consider the solvation effect . In addition, AIMD simulations were carried out to explore the thermodynamic stability of TM-P4-C catalysts at 300K. The Phonopy was adopted to calculate the harmonic phonon dispersion.