The Activity Improvement of the TM3(hexaiminotriphenylene)2 Monolayer for Oxygen Reduction Electrocatalysis: A Density Functional Theory Study.

Polymer electrolyte membrane fuel cells (PEMFCs) are one of the most prominent clean energy technologies designed to achieve hydrogen utilization and solve problems such as low efficiency and high pollution associated with fossil fuel combustion. In order to bring about PEMFC commercialization, especially for automobile applications, developing high-activity and -selectivity catalysts for the oxygen reduction reaction (ORR) is of critical importance. Based on the density functional theory, the catalytic activity of the conductive, two-dimensional metal-organic frameworks TM3(HITP)2 monolayer (where HITP = hexaiminotriphenylene; TM = Ni, Co, Fe, Pd, Rh, Ru, Pt, Ir, and Os) for ORR has been investigated systematically. Furthermore, the classical volcano curves of the ORR activity, as a function of the OH binding, are found where the Ni, Pd, and Pt located at the weak binding side suffer from the sluggish *OOH formation and prefer the inefficient 2e - mechanism, while for other elements belonging to the strong binding side, the reactions are hindered by the poison due to ORR intermediates. Based on the free energy profiles, the corresponding overpotentials μORR exhibit the inverted volcano curve as a function of the atomic number of the 3d/4d/5d TM active center in the same period. Based on the μORR data, ORR activity decreases in the order of Ir > Co ≈ Rh > Ni ≈ Pd > Pt ≈ Fe > Ru > Os. Herein, the Co, Rh, and Ir central atoms exhibit enhanced catalytic activity in combination with the desirable selectivity of the O2 reduction to H2O. This systematic work may open new avenues for the development of high-performance non-PGM catalysts for practical applications of ORR.

Introduction

Hydrogen is a potential candidate for future energy provision, as stated in the concept of the hydrogen economy, so as to solve the issues of the rising global energy demands, depletion of fossil fuel reserves, and associated environmental pollution issues. Due to its high efficiency, ease of operation, and low emission, the polymer electrolyte membrane fuel cell (PEMFCs) is the most prominent technology to derive benefit from the proposed hydrogen cycle, which leads to the production of electricity from the electrochemical oxidation of hydrogen, with water as its endproduct (Colić and Bandarenka, 2016; Xia Z. et al., 2016; Chen et al., 2017, 2018). However, a critical obstacle to its commercialization is the dominant voltage loss associated with the sluggish oxygen reduction reaction (ORR), even when catalyzed by the noble Pt nanoparticle (Nørskov et al., 2004). In this regard, a significant amount of Pt is required to achieve the desirable power density, making the overall cost prohibitively high (Debe, 2012). In order to overcome the economic bottleneck, design, and application of earth-abundant alternatives for ORR electrocatalysis are at the heart of PEMFC research (Xia W. et al., 2016).

Inspired by the pioneering work on cobalt phthalocyanines acting as cathode catalysts (Jasinski, 1964), a tremendous amount of research is being performed on the TM–Nx carbon materials for ORR (TM denotes transition metals), especially TM–N4 active motif. The intrinsic active characteristic of TM in Nx carbon materials is experimentally supported, where the presence of the TM atom would boost the ORR performance compared with the metal-free counterpart (Peng et al., 2013; Yin et al., 2016; Yang et al., 2017); the corresponding activity would be suppressed by adding the SCN− ions and CO molecule (Zhang et al., 2016; Yang et al., 2017). Linear correlations between the content of TM–Nx and ORR activity have been observed (Yang et al., 2017). Furthermore, ORR activity shows its dependence on the TM center atom, and this is supported by the density functional theory (DFT), a theoretical work of Rossmeisl et al., where the ORR activity of the TM–N4 embedded graphene has been systematically investigated, and where the Fe, Ir, Mn, Ru, and Rh doping is identified as boosting the ORR (Calle-Vallejo et al., 2011). Additionally, the ORR mechanism is sensitive to the TM–N4 active motif. Liu et al. have synthesized the Fe–Nx and Co–Nx doped carbon nanofiber and realized that the Fe–N4 promotes 4e− ORR in comparison with the 2e− pathway of the Co–N4 one (Liu et al., 2016). The same conclusion has been achieved by Kattel et al., where the O–O bond scission is presented and the efficient 4e− reduction of O2 to 2H2O is preferred on the Fe–N4 sites, but the reduction of O2 to H2O2 is found enhanced on the Co–N4 motif (Kattel et al., 2013, 2014).

Despite these encouraging research works, such TM–N4 carbon materials generally suffer from low activity caused by the relatively few catalytic sites as well as the experimental challenge of the well-controlled active motif (Liang et al., 2013; Peng et al., 2013; Palaniselvam et al., 2016), reducing the competition with the state-of-the-art Pt/C catalysts. Paying attention to the high TM–N4 density in combination with the electronic conductivity, Miner et al. have developed the attractive Ni3(HITP)2 as ORR electrocatalyst (Miner et al., 2016). However, the performance is far from satisfying expectations, with its incomplete oxidation of O2 and the predominant production of H2O2 under the working potential (Miner et al., 2016). In this regard, utilizing the aforementioned information, the development of the efficient TM3(HITP)2 catalysts can be achieved via the variation of metal nodes (Choi et al., 2015; Zhang et al., 2015). It has naturally raised the interest to search for the optimal combinations of the TM3(HITP)2, possessing superior ORR activity as well as selectivity.

In the manuscript, DFT calculations are used within an electrochemical framework to analyze the ORR electrocatalysis over the TM3(HITP)2monolayer. In particular, the stability of the ORR intermediates is calculated, which allows to evaluate the thermodynamic ORR free energy and its overpotentials. The data provide the fundamental understanding of the mechanism of Ni3(HITP)2 and further identify optimal candidates as catalysts. According the d–partial density of states, an atomistic insight of the activity origin has been provided by a thorough comparison among the considered systems. Herein, for simplification, our attention is mainly focused on the monolayer structure shown in Figure 1 as a representative model due to the weak interaction between the interlayers of the 2D-layered materials (Sheberla et al., 2014; Chen et al., 2015; Miner et al., 2016).

Figure 1

The schematic structure of the TM3(HITP)2 monolayer. Green, blue, gray, and white denote the TM, N, C, and H atoms.

Methods

All calculations have been performed within the DFT framework, as implemented in the DMol3 code (Delley, 1990, 2000). The generalized gradient approximation with the Perdew–Burke–Ernzerhof (PBE) functional is employed to describe exchange and correlation effects (Perdew et al., 1996). The DFT semicore pseudopotentials (DSPP) core treat method is implemented for relativistic effects, which replaces core electrons by a single effective potential and introduces some degree of relativistic correction into the core (Delley, 2002). The double numerical atomic orbital augmented by a polarization function (DNP) is chosen as the basis set (Delley, 1990). A smearing of 0.005 Ha (1 Ha = 27.21 eV) to the orbital occupation is applied to achieve accurate electronic convergence. In order to ensure high-quality results, the real-space global orbital cutoff radius is set as high as 5.2 Å. In order to accurately describe the long-range electrostatic interactions of the ORR intermediates with catalysts, the PBE-D method with the TS van der Waals (vdW) correction is employed. In the geometry optimization of structures, the convergence tolerances of energy, maximum force, and displacement are 1.0 × 10−5 Ha, 0.002 Ha/Å, and 0.005 Å, respectively. The spin-unrestricted method is used for all calculations. A conductor-like screening model (COSMO) was used to simulate a H2O solvent environment for all calculations (Todorova and Delley, 2008), which is a continuum model, where the solute molecule forms a cavity within the dielectric continuum. The DMol3/COSMO method has been generalized to periodic boundary cases. The dielectric constant is set as 78.54 for H2O. Some previous results have shown that this implicit solvation model is an effective method to describe solvation (Sha et al., 2011; Zhang et al., 2015). The 15 Å-thick vacuum is added to avoid the artificial interactions between the TM3(HITP)2 monolayer and its images. The proposed structure of the TM3(HITP)2 monolayer is shown in Figure 1, where the atoms in the red square are fixed in all of the structure-optimization calculations.

The adsorption energy (Eads) of the reactant O2 is calculated by the following,

The adsorption energy (Eads) of the ORR intermediates are calculated relative to H2O and H2 (Calle-Vallejo et al., 2011), through,

where Esys, Esubstrate, EH2O, and EH2 are the total energy of the adsorption systems, the TM3(HITP)2 monolayer, H2O molecule, and H2 molecule, respectively. The Eads < 0 corresponds to an exothermic adsorption process.

The Gibbs free energy changes (ΔG) of the ORR elemental steps have been calculated according to the computational hydrogen electrode (CHE) model developed by Nørskov et al., where the chemical potential of the proton/electron (H+ + e−) in solution is equal to the half of the chemical potential of a gas-phase H2 at the reference relative hydrogen electrode (RHE) potential (Nørskov et al., 2004). Herein, G(H++e−) = 1/2G(H2). The variance of the proton–electron pair free energy with potential is simply determined using the linear free energy dependence of the electron energy on potential, shifting the electron energy –eU, where e is the elementary positive charge and U is the electrode potential of interest on the RHE scale (Nørskov et al., 2004; Nie et al., 2014). The total chemical potential of a proton–electron pair at the potential U is written as follows:

Therefore, for a general electrochemical reaction, the free energy change ΔG for every elemental step can be determined as following:

where ΔE is the electronic energy difference based on DFT calculations, ΔZPE is the change in zero point energy, T is the temperature (equal to 298.15 K here), ΔS is the change in the entropy, and ΔGpH and ΔGU are the free energy contributions due to variation in pH value and electrode potential U, respectively. ΔGU = –eU, in which U is the potential related to the standard hydrogen electrode. ΔGpH = –kTln10 x pH, which is the corrected free energy of H+ ions depending on the concentration. According to the previous works, pH is set as 0 in acid medium and 14 in alkaline medium. The zero-point energies and entropies of the ORR intermediates are calculated from the vibrational frequencies according to standard methods. Following the suggestion of Wilcox et al. (Lim and Wilcox, 2012), in order to reduce the calculation, the TM3(HITP)2 monolayer is fully constrained. The ΔG < 0 corresponds to an exothermic adsorption process. The free energy G of O2 is derived as G(O2) = 4.92+2G(H2O)−2G(H2) by utilizing OER equilibrium at the standard conditions; the G of H2O2 is derived similarly considering that the H2O2/O2 standard equilibrium potential is 0.682 V vs. SHE (Sun et al., 2014). The CHE model has been successfully applied for developing the novel electrocatalysts with prominent ORR performances, where the DFT calculations are in line with the experimental results (Nørskov et al., 2004; Greeley et al., 2009; Favaro et al., 2015; Lang et al., 2015; Jia et al., 2016; Tang et al., 2016; Liu et al., 2017; Li et al., 2018; Xu et al., 2018). Furthermore, the PBE/DNP method in Dmol3 code has been widely employed for evaluating the potential of the TM-based carbon electrocatalysts (Wang et al., 2015, 2016b; Hou et al., 2016; Xiao et al., 2017). Therefore, the reliability of the mentioned approach is confirmed.

Results and discussion

Prior to the investigation of the activity, the essential step is to determine the adsorption behavior of the ORR intermediates. For Ni3(HITP)2, the energetics of the O2 adsorption is endothermic with the value of 0.13 eV, caused by the structure deformation as shown in the inset of Figure 2A, indicating no such adsorption has occurred on the Ni center atom. Besides the reactant, the corresponding Eads of ORR intermediates are 4.08, 3.57, and 1.17 eV for OOH, O, and OH, respectively. In comparison with the Pt(111) (3.55 eV for OOH, 1.38 eV for O, and 0.64 eV for OH) (Xiao et al., 2016), the obviously weak adsorption ability of the center Ni atom has been observed.

Figure 2

The adsorption energies of the O2 molecules (A) and the ORR intermediates (B). Inset: the representative O2 adsorption structure on Ni, Co, and Fe active centers.

Considering the low affinity of the ORR intermediates to the Ni center, it is expected that the ORR activity would be boosted by selecting suitable TMs. Indeed, the change of the TM center definitely exhibits different adsorption behaviors, as shown in Figure 2, where the corresponding Eads data are given. From the figures, the Eads of the reactant O2 and the corresponding ORR intermediates decrease monotonically in the order of Ru ≈ Os > Fe > Co > Rh ≈ Ir > Ni > Pd ≈ Pt. It is obvious that the adsorption ability of the elements is weakened from group 8 to group 10. Carefully reviewing the O2 adsorption as shown in the inset of Figure 2A, similar to Ni, no such adsorption behaviors have occurred at the Pd and Pt active centers due to the endothermic adsorption energies. The unfavorable adsorption behavior of group 10 implies the underlying mechanism of the OOH formation, which originates from the long-range electron transfer from the catalyst to O2 molecules at the outer Helmholtz plane (Choi et al., 2014). On the contrary, group 9 binds to the reactant suitably with the corresponding end-on structure and group 8 possesses a too strong O2 adsorption with the side-on adsorption configuration. Herein, the mentioned adsorption capability is in line with the general prediction of the d-band model (Hammer and Nørskov, 2000). The corresponding evidences are provided by the d–partial density of states (d–PDOS) of the central TM atom plotted in Figure 3. As shown, the d band moves away from the Fermi energy as the TM atom changes from group 8 to group 10. The TM atom with the higher d states possesses stronger adsorption ability. Furthermore, due to the relationship between the adsorption behavior and the d band, the linear scaling relations between Eads(O2)/Eads(OOH)/Eads(O) and Eads(OH) are observed, in accordance with the previous reports of other functionalized carbon materials (Calle-Vallejo et al., 2011; Baran et al., 2014). That is,

Figure 3

The d-PDOS of the central TM atom.

To further evaluate the ORR activity of TM3(HITP)2, according to the experimental condition (Miner et al., 2016), the OOH associative mechanisms in the alkaline solution are taken into consideration, with the elemental steps R listed as follows (Wang et al., 2016a), where asterisks denote active TM sites. Due to the small barrier of proton transfer, which can be ignored at high applied voltages, our attention is only focused on the reaction energies (Nørskov et al., 2004; Calle-Vallejo et al., 2011; Zhang et al., 2015).

Analyzing the free energy plots of the complete 4e− ORR pathway in Figure 4A, the endothermic processes of *OOH formation (R1) and *O formation (R2a) are observed even at U = 0 V. The corresponding values of the free energy change ΔG are 0.28 and 0.08 eV, respectively. The positive values indicate that the mentioned steps are thermodynamically unfavorable (Nørskov et al., 2004). When the ideal potential of 0.4 eV is applied, the mentioned reaction steps are even unfavorable with the increased ΔG of 0.68 and 0.48 eV, respectively. Based on the information, R1 is determined as the rate-determining step (RDS), pointing to the fact that the Ni3(HITP)2 monolayer suffers from insufficient O2 activation (Greeley et al., 2009), in accordance with the endothermic capture of the O2 molecule as shown in Figure 2A. To clarify the selectivity of the 2/4e− mechanism, the comparison between R2a and R2b is considered. In the alternative 2e− pathway, the ΔG2b of the H2O2 formation (R2b) is increased from 0.01 to 0.41 eV as the U ranged from 0 to 0.4 V, being 0.07 eV lower than ΔG2a, indicating the slightly energy favorable condition of R2b. Herein, the predominant production of H2O2 would be expected on the Ni3(HITP)2 monolayer, which is in consistence with the experiment results of the 3D Ni3(HITP)2. In summary, the Ni3(HITP)2 monolayer slightly prefers the 2e− mechanism, with the RDS located at the *OOH formation.

Figure 4

The free energy G diagrams of Ni3(HITP)2 in the alkaline medium (A) and the free energy change ΔGi of the elemental step Ri at the ideal potential of 0.4 eV (B).

Owing to the different adsorption abilities of TM3(HITP)2, the suitable TM would boost the ORR activity as well as the selectivity of the 4e− mechanism. In order to characterize the relationship between the ORR activity and the TM center atom, the reactive free energy change ΔG of the elementary steps R at U = 0.4 V as a function of Eads(OH) is plotted in Figure 4B. As the weakening of the Eads(OH) occurs with the values ranging from −0.05 to 1.68 eV, ΔG1 and ΔG2a increase from the negative to the positive values, while the opposite tendency is found for ΔG3 and ΔG4, which decrease from the positive to the negative values. That is, the steps R1 and R2a change from the exothermic to endothermic process, while the steps R3 and R4 become more thermodynamically favorable. Obviously, the RDS steps are identified as R4 for group 8 and group 9 with the exception of Os, where the RDS step is R3, while the RDS step is located at R1 for group 10.

On the basis of the assumption that the activation barrier for the RDS is equal to ΔGmax, the activity variation of TM3(HITP)2 referred to as Ni3(HITP)2 is estimated by the rate constant k via the simple Arrhenius equation in the following,

where A is the prefactor, kB is the Boltzmann constant, and T is the temperature (298.15 K). The classic volcano curve is observed in Figure 5A, where the In(kTM/kNi) as a function of the Eads(OH) is plotted. Based on our results, the adsorption ability is weakened as the atomic number increases from group 8 to group 10. The Ni/Pd/Pt active centers suffer from the too weak binding strength, leading to the energetically unfavorable process of the *OOH formation. Meanwhile, for the Fe/Ru/Os active centers, the too strong interaction with the ORR intermediates accounts for the poisoning of the O-containing intermediates. The suitable binding strength of Co/Rh/Ir active centers indicates the balance between the O2 activation and the catalyst recovery. As discussed by the previous reports (Calle-Vallejo et al., 2011; Viswanathan et al., 2012a; Baran et al., 2014; Zhang et al., 2015), the ORR activity depends on the adsorption of the intermediates. Our results are in accordance with the previous reports in that the bond strength should be compromised in the case of the effective ORR catalysts on the basis of Sabatier principle (Greeley et al., 2009; Gao et al., 2017, 2018); this is similar with the previous studies on metal (Stephens et al., 2012), functional graphene (Calle-Vallejo et al., 2011), metal porphyrine (Baran et al., 2014) as well as the 2D MOF (Gao et al., 2017). Herein, in comparison with the Ni3(HITP)2 monolayer (0.68 eV), group 9 possesses superior catalytic performance with the smaller ΔGmax of 0.50, 0.52, and 0.22 eV for the Co, Rh, and Ir center atoms, respectively. By contrast, inferior activities are observed for group 8 and group 10. The corresponding ΔGmax are 0.93, 1.23, and 1.34 eV for Fe, Ru, and Os atoms, and 0.72 and 0.90 eV for Pd and Pt atoms, respectively. Furthermore, for clarification, the data of the relative overpotentials μORR are summarized and collected for the activity comparison, depicted in Figure 6. As the atomic number of TM increased from Ni to Fe, the μORR reduces and then increases, presenting further evidence to the presence of the classical volcano-shaped activity (Nørskov et al., 2004; Calle-Vallejo et al., 2011; Baran et al., 2014; Zheng et al., 2017). Similar situations are found for the 4d/5d TM3(HITP)2 monolayer. Based on the μORR data, ORR activity decreases in the order of Ir > Co ≈ Rh > Ni ≈ Pd > Pt ≈ Fe > Ru > Os. Furthermore, compared with the data of the TM supported on graphene and macrocyclic molecules (the minimum μORR = 0.4 eV) (Xu et al., 2018), the prominent improvement of the ORR activity on the Ir3(HITP)2 monolayer is confirmed.

Figure 5

The activity enhancement factor of TM3(HITP)2 referred to the Ni3(HITP)2 (A) and the selectivity of the 4e− reduction referred to the 2e− pathway (B).

Figure 6

The overpotentials μORR of the TM3(HITP)2 monolayer for the 4e− mechanism.

Besides the ORR activity, the ORR selectivity is characterized by the ratio of the corresponding reaction rate k2a (k2b) for step R2a (R2b), via the following equation,

The data are plotted in Figure 5B. As shown, In(k2a/k2b) shows a linear relationship with Eads(OH). For the weak side, the values of ΔG2a and ΔG2b are 0.48 and 0.41 eV for Ni, 0.57 and 0.37 eV for Pd, and 0.20 and 0.19 eV for Pt, respectively. The corresponding In(k2a/k2b) are −2.72, −7.78, and −0.39, indicating the energy favorability of R2b compared with that of R2a for group 10. Furthermore, considering the In(k2a/k2b) data, the 2e− reduction of O2 to H2O2 is prevalent at the Ni and Pd centers, while the mixing of 2/4e− ORR is reasonable for the Pt center. By contrast, R2a is preferred for other elements indicated by the positive In(k2a/k2b). That is, the 4e− mechanism is dominant for group 8 and group 9. Obviously, as the Eads(OH) is strengthened, the catalytic selectivity is changed from the 2e− to 4e− pathway. The two sides of the catalytic activity profile shown in Figure 5A essentially distinguish the 2e− catalysts (weak binding side) from the 4e− catalysts (strong binding side) (Viswanathan et al., 2012b; Zagal and Koper, 2016). That is, the former suffering from the insufficient O2 activation favors the H2O2 formation; the latter poisoned by the O–containing intermediates generally prefers the H2O formation. The formation of H2O2 not only decreases the ORR efficiency, but also degrades the proton exchange membrane (Tsuneda et al., 2017). In short, the dramatic enhancements in oxygen-reduction rates and product selectivity are achieved by selecting TM centers in group 9, in line with the previous works (Gao et al., 2017; Wannakao et al., 2017).

To understand the influence of the long-term electrostatic potential on the ORR performance, the vdW correction is further applied for analyzing the optimum Ir/Co/Rh systems, where the corresponding adsorption energy and the free energy changes are listed in Table 1. Intuitively, the adsorption abilities of the mentioned systems are enhanced due to the presence of the long-term interaction, in consistence with our calculation data. For the Ir3(HITP)2 monolayer, the Eads are slightly increased to −0.39, 3.49, 2.11, and 0.61 eV for O2, OOH, O, and OH, respectively, in comparison with the uncorrected values of −0.31, 3.59, 2.16, and 0.68 eV. The perturbation of the binding strength hinders the OH removal from the Ir active site, leading to the uphill of the ΔGmax from 0.22 to 0.30 eV. Despite the activity decrease, the 4e− reduction pathway remains, supported by the values of ΔG2a with −0.23 eV and ΔG2b with 1.26 eV. Herein, the same phenomenon is found for Co3(HITP)2 and Rh3(HITP)2 with the ΔGmax of 0.68 and 0.59 eV, respectively. The corresponding ORR activity with th 4e− reaction mechanism follows the order of Ir > Rh ≈ Co. Despite the numerical variation, the trend is roughly consistent with the results without considering vdW interaction.

Table 1

The Eads and ΔG of ORR intermediates at the potential U = 0.4V with and without the vdW corrections.

	Without vdW			With vdW
	Co	Rh	Ir	Co	Rh	Ir
Eads (O2)	−0.64	−0.41	−0.31	−0.67	−0.50	−0.39
Eads (OOH)	3.42	3.52	3.59	3.20	3.41	3.49
Eads (O)	2.53	2.59	2.16	2.54	2.55	2.11
Eads (OH)	0.46	0.61	0.68	0.28	0.54	0.61
ΔG1	−0.08	−0.12	−0.06	−0.31	−0.23	−0.16
ΔG2a	0.03	0.13	−0.29	0.26	0.20	−0.23
ΔG2b	1.18	1.21	1.15	1.40	1.32	1.26
ΔG3	−0.44	−0.53	0.12	−0.64	−0.56	0.09
ΔG4	0.50	0.52	0.22	0.68	0.59	0.30
ΔGmax	0.50	0.52	0.22	0.68	0.59	0.30

Despite the presence of the Ni3(HITP)2 experimentally, other 2D catalysts are only theoretical models, which need the confirmation of their synthesis. It should be noted that changing metal atoms would significantly modify the corresponding structures. Thus, by replacing the Ni central atom, the dimensionality of the Cr3(HITP)2 could be transferred from 2D to 3D due to the energetic favorable insertion of the spacer linker (Foster et al., 2016). Different structures inevitably lead to distinct catalytic performances (Sun et al., 2014). In this regard, although the theoretical candidates have been rationally predicted, the ORR performance of the optimum TM3(HITP)2 with the Co/Rh/Ir active sites crucially needs further experimentally verification.

Conclusion

Based on the DFT, the ORR mechanisms on TM3(HITP)2 monolayer have been studied. It is found that the selection of central metals affects the adsorption behaviors, tuning ORR activity and its selectivity. It is realized that the adsorption abilities are monotonously enhanced as the d band upshifts from group 10 to group 8. A classical volcano relationship for the predicted ORR activity as a function of calculated OH adsorption energy is observed. From the calculated ΔG data, the ORR activity decreases in the order of Ir > Co ≈ Rh > Ni ≈ Pd > Pt ≈ Fe > Ru > Os. That is, group 9 possesses superior activity compared with other elements. Furthermore, owing to the insufficient O2 activation, the 2e− mechanism is prevalent in group 10, while the desirable 4e− mechanism is dominant for others. These results would throw insights into the nature of the ORR mechanisms of TM3(HITP)2. The materials of the Co3(HITP)2, Rh3(HITP)2, and Ir3(HITP)2 have been screened out and served as the potential candidates for the ORR electrocatalysis.

Author contributions

BX carried out the simulation and wrote the paper. QJ revised the paper. HZ, HL, and XJ entered the discussion. All authors commented on the manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.