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;Chen et al., 2017Chen et al., , 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-N x carbon materials for ORR (TM denotes transition metals), especially TM-N 4 active motif. The intrinsic active characteristic of TM in N x 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 Yang et al., 2017). Linear correlations between the content of TM-N x 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-N 4 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-N 4 active motif. Liu et al. have synthesized the Fe-N x and Co-N x doped carbon nanofiber and realized that the Fe-N 4 promotes 4e − ORR in comparison with the 2e − pathway of the Co-N 4 one . The same conclusion has been achieved by Kattel et al., where the O-O bond scission is presented and the efficient 4e − reduction of O 2 to 2H 2 O is preferred on the Fe-N 4 sites, but the reduction of O 2 to H 2 O 2 is found enhanced on the Co-N 4 motif (Kattel et al., 2013(Kattel et al., , 2014. Despite these encouraging research works, such TM-N 4 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 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-N 4 density in combination with the electronic conductivity, Miner et al. have developed the attractive Ni 3 (HITP) 2 as ORR electrocatalyst (Miner et al., 2016). However, the performance is far from satisfying expectations, with its incomplete oxidation of O 2 and the predominant production of H 2 O 2 under the working potential (Miner et al., 2016). In this regard, utilizing the aforementioned information, the development of the efficient TM 3 (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 TM 3 (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 TM 3 (HITP) 2 monolayer. 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 Ni 3 (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).

METHODS
All calculations have been performed within the DFT framework, as implemented in the DMol 3 code (Delley, 1990(Delley, , 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 FIGURE 1 | The schematic structure of the TM 3 (HITP) 2 monolayer. Green, blue, gray, and white denote the TM, N, C, and H atoms.
1.0×10 −5 Ha, 0.002 Ha/Å, and 0.005 Å, respectively. The spinunrestricted method is used for all calculations. A conductor-like screening model (COSMO) was used to simulate a H 2 O 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 DMol 3 /COSMO method has been generalized to periodic boundary cases. The dielectric constant is set as 78.54 for H 2 O. 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 TM 3 (HITP) 2 monolayer and its images. The proposed structure of the TM 3 (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 (E ads ) of the reactant O 2 is calculated by the following, where E sys , E substrate , E H2O , and E H2 are the total energy of the adsorption systems, the TM 3 (HITP) 2 monolayer, H 2 O molecule, and H 2 molecule, respectively. The E ads < 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 H 2 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 G pH and G U are the free energy contributions due to variation in pH value and electrode potential U, respectively. G U = -eU, in which U is the potential related to the standard hydrogen electrode. G pH = -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 TM 3 (HITP) 2 monolayer is fully constrained. The G < 0 corresponds to an exothermic adsorption process. The free energy G of O 2 is derived as  (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 Dmol 3 code has been widely employed for evaluating the potential of the TM-based carbon electrocatalysts (Wang et al., 2015(Wang et al., , 2016bHou 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 Ni 3 (HITP) 2 , the energetics of the O 2 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 E ads 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.
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 E ads data are given. From the figures, the E ads of the reactant O 2 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 O 2 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 O 2 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 O 2 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 dpartial 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 E ads (O 2 )/E ads (OOH)/E ads (O) and E ads (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, E ads (OOH) = 0.90E ads (OH) + 3.02, E ads (O) = 1.88E ads (OH) + 1.02.
To further evaluate the ORR activity of TM 3 (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 i 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 (R 1 ) and * O formation (R 2a ) are observed even at U = 0 V. The corresponding values of the free energy change G i 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 i of 0.68 and 0.48 eV, respectively. Based on the information, R 1 is determined as the rate-determining step (RDS), pointing to the fact that the Ni 3 (HITP) 2 monolayer suffers from insufficient O 2 activation (Greeley et al., 2009), in accordance with the endothermic capture of the O 2 molecule as shown in Figure 2A.
To clarify the selectivity of the 2/4e − mechanism, the comparison between R 2a and R 2b is considered. In the alternative 2e − pathway, the G 2b of the H 2 O 2 formation (R 2b ) 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 G 2a , indicating the slightly energy favorable condition of R 2b . Herein, the predominant production of H 2 O 2 would be expected on the Ni 3 (HITP) 2 monolayer, which is in consistence with the experiment results of the 3D Ni 3 (HITP) 2 . In summary, the Ni 3 (HITP) 2 monolayer slightly prefers the 2e − mechanism, with the RDS located at the * OOH formation.
Owing to the different adsorption abilities of TM 3 (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 i of the elementary steps R i at U = 0.4 V as a function of E ads (OH) is plotted in Figure 4B. As the weakening of the E ads (OH) occurs with the values ranging from −0.05 to 1.68 eV, G 1 and G 2a increase from the negative to the positive values, while the opposite tendency is found for G 3 and G 4 , which decrease from the positive to the negative values. That is, the steps R 1 and R 2a change from the exothermic to endothermic process, while the steps R 3 and R 4 become more thermodynamically favorable. Obviously, the RDS steps are identified as R 4 for group 8 and group 9 with the exception of Os, where the RDS step is R 3 , while the RDS step is located at R 1 for group 10.
Frontiers in Chemistry | www.frontiersin.org On the basis of the assumption that the activation barrier for the RDS is equal to G max , the activity variation of TM 3 (HITP) 2 referred to as Ni 3 (HITP) 2 is estimated by the rate constant k via the simple Arrhenius equation in the following, where A is the prefactor, k B is the Boltzmann constant, and T is the temperature (298.15 K). The classic volcano curve is observed in Figure 5A, where the In(k TM /k Ni ) as a function of the E ads (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 O 2 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., 2017Gao et al., , 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 Ni 3 (HITP) 2 monolayer (0.68 eV), group 9 possesses superior catalytic performance with the smaller G max 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 G max 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 TM 3 (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) , the prominent improvement of the ORR activity on the Ir 3 (HITP) 2 monolayer is confirmed. Besides the ORR activity, the ORR selectivity is characterized by the ratio of the corresponding reaction rate k 2a (k 2b ) for step R 2a (R 2b ), via the following equation, The data are plotted in Figure 5B. As shown, In(k 2a /k 2b ) shows a linear relationship with E ads (OH). For the weak side, the values of G 2a and G 2b 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(k 2a /k 2b ) are −2.72, −7.78, and −0.39, indicating the energy favorability of R 2b compared with that of R 2a for group 10. Furthermore, considering the In(k 2a /k 2b ) data, the 2e − reduction of O 2 to H 2 O 2 is prevalent at the Ni and Pd centers, while the mixing of 2/4e − ORR is reasonable for the Pt center. By contrast, R 2a is preferred for other elements indicated by the positive In(k 2a /k 2b ). That is, the 4e − mechanism is dominant for group 8 and group 9. Obviously, as the E ads (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 O 2 activation favors the H 2 O 2 formation; the latter poisoned by the O-containing intermediates generally prefers the H 2 O formation. The formation of H 2 O 2 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 longterm interaction, in consistence with our calculation data. For the Ir 3 (HITP) 2 monolayer, the E ads are slightly increased to −0.39, 3.49, 2.11, and 0.61 eV for O 2 , 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 G max from 0.22 to 0.30 eV. Despite the activity decrease, the 4e − reduction pathway remains, supported by the values of G 2a with −0.23 eV and G 2b with 1.26 eV. Herein, the same phenomenon is found for Co 3 (HITP) 2 and Rh 3 (HITP) 2 with the G max 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.
Despite the presence of the Ni 3 (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 Cr 3 (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 . In this regard, although the theoretical candidates have been rationally predicted, the ORR performance of the optimum TM 3 (HITP) 2 with the Co/Rh/Ir active sites crucially needs further experimentally verification.

CONCLUSION
Based on the DFT, the ORR mechanisms on TM 3 (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 O 2 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 TM 3 (HITP) 2 . The materials of the Co 3 (HITP) 2 , Rh 3 (HITP) 2 , and Ir 3 (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.