The Activity Improvement of the TM₃(hexaiminotriphenylene)₂ Monolayer for Oxygen Reduction Electrocatalysis: A Density Functional Theory Study

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–N_x carbon materials for ORR (TM denotes transition metals), especially TM–N₄ 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 (Zhang et al., 2016; 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₄ 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₄ active motif. Liu et al. have synthesized the Fe–N_x and Co–N_x doped carbon nanofiber and realized that the Fe–N₄ promotes 4e⁻ ORR in comparison with the 2e⁻ pathway of the Co–N₄ 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 O₂ to 2H₂O is preferred on the Fe–N₄ sites, but the reduction of O₂ to H₂O₂ is found enhanced on the Co–N₄ motif (Kattel et al., 2013, 2014).

Despite these encouraging research works, such TM–N₄ 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–N₄ density in combination with the electronic conductivity, Miner et al. have developed the attractive Ni₃(HITP)₂ as ORR electrocatalyst (Miner et al., 2016). However, the performance is far from satisfying expectations, with its incomplete oxidation of O₂ and the predominant production of H₂O₂ under the working potential (Miner et al., 2016). In this regard, utilizing the aforementioned information, the development of the efficient TM₃(HITP)₂ 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₃(HITP)₂, 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₃(HITP)₂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₃(HITP)₂ 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 TM₃(HITP)₂ 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 DMol³ 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⁻⁵ 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 H₂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³/COSMO method has been generalized to periodic boundary cases. The dielectric constant is set as 78.54 for H₂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₃(HITP)₂ monolayer and its images. The proposed structure of the TM₃(HITP)₂ 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₂ is calculated by the following,

The adsorption energy (E_ads) of the ORR intermediates are calculated relative to H₂O and H₂ (Calle-Vallejo et al., 2011), through,

where E_sys, E_substrate, E_H2O, and E_H2 are the total energy of the adsorption systems, the TM₃(HITP)₂ monolayer, H₂O molecule, and H₂ 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₂ 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₃(HITP)₂ monolayer is fully constrained. The ΔG < 0 corresponds to an exothermic adsorption process. The free energy G of O₂ is derived as G(O₂) = 4.92+2G(H₂O)−2G(H₂) by utilizing OER equilibrium at the standard conditions; the G of H₂O₂ is derived similarly considering that the H₂O₂/O₂ 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 Dmol³ 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 Ni₃(HITP)₂, the energetics of the O₂ 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.

Figure 2. The adsorption energies of the O₂ molecules (A) and the ORR intermediates (B). Inset: the representative O₂ 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 E_ads data are given. From the figures, the E_ads of the reactant O₂ 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₂ 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₂ 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₂ 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 E_ads(O₂)/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,

To further evaluate the ORR activity of TM₃(HITP)₂, 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₁) 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₁ is determined as the rate-determining step (RDS), pointing to the fact that the Ni₃(HITP)₂ monolayer suffers from insufficient O₂ activation (Greeley et al., 2009), in accordance with the endothermic capture of the O₂ 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₂O₂ 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₂O₂ would be expected on the Ni₃(HITP)₂ monolayer, which is in consistence with the experiment results of the 3D Ni₃(HITP)₂. In summary, the Ni₃(HITP)₂ monolayer slightly prefers the 2e⁻ mechanism, with the RDS located at the *OOH formation.

Figure 3. The d-PDOS of the central TM atom.

Figure 4. The free energy G diagrams of Ni₃(HITP)₂ in the alkaline medium (A) and the free energy change ΔG_i of the elemental step R_i at the ideal potential of 0.4 eV (B).

Owing to the different adsorption abilities of TM₃(HITP)₂, 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₁ and ΔG_2a increase from the negative to the positive values, while the opposite tendency is found for ΔG₃ and ΔG₄, which decrease from the positive to the negative values. That is, the steps R₁ and R_2a change from the exothermic to endothermic process, while the steps R₃ and R₄ become more thermodynamically favorable. Obviously, the RDS steps are identified as R₄ for group 8 and group 9 with the exception of Os, where the RDS step is R₃, while the RDS step is located at R₁ for group 10.

On the basis of the assumption that the activation barrier for the RDS is equal to ΔG_max, the activity variation of TM₃(HITP)₂ referred to as Ni₃(HITP)₂ 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₂ 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 Ni₃(HITP)₂ 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₃(HITP)₂ 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 Ir₃(HITP)₂ monolayer is confirmed.

Figure 5. The activity enhancement factor of TM₃(HITP)₂ referred to the Ni₃(HITP)₂ (A) and the selectivity of the 4e⁻ reduction referred to the 2e⁻ pathway (B).

Figure 6. The overpotentials μ_ORR of the TM₃(HITP)₂ monolayer for the 4e⁻ mechanism.

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₂ to H₂O₂ 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₂ activation favors the H₂O₂ formation; the latter poisoned by the O–containing intermediates generally prefers the H₂O formation. The formation of H₂O₂ 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 Ir₃(HITP)₂ monolayer, the E_ads are slightly increased to −0.39, 3.49, 2.11, and 0.61 eV for O₂, 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₃(HITP)₂ and Rh₃(HITP)₂ 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.

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

Despite the presence of the Ni₃(HITP)₂ 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₃(HITP)₂ 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 TM₃(HITP)₂ with the Co/Rh/Ir active sites crucially needs further experimentally verification.

Conclusion

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

Acknowledgments

We acknowledge the support from the National Natural Science Foundation of China (No. 21503097, 51631004) and the Natural Science Foundation of Jiangsu (No.BK20140518).

References