Considering increasing global energy and environmental pollution crises caused by fossil fuel depletion, the development of renewable energy technologies plays a pivotal role in sustainable applications [1,2]. Among various electrochemical energy conversion technologies, water splitting systems, fuel cells and metal-air batteries have received tremendous interests owing to their high energy density and environmentally benign [3–5]. Importantly, the hydrogen evolution reaction (HER), oxygen evolution reaction (OER) and oxygen reduction reaction (ORR) are essential reactions in these electrochemical processes. The HER at the cathode and the OER at the anode sides of water splitting electrolysis can produce clean hydrogen and oxygen gas, respectively [6,7]. The reversible OER and ORR reactions play critical roles in determining the performance of fuel cells and metal-air batteries [8–10]. Meanwhile, these electrochemical reactions are kinetically sluggish and need to overcome large overpotential to occur. Currently, Pt- and Ir/Ru-based materials are the benchmark catalysts for the HER/ORR and OER, respectively [9,11]. However, their high cost and low abundance limit their large-scale application [3]. Therefore, searching for alternative catalysts with efficient activity and low cost is highly desirable.

Recently, bifunctional electrocatalysts for water splitting, fuel cells and metal-air batteries have attracted tremendous attention due to their intrinsic advantages [4]. For instance, bifunctional electrocatalysts for both HER and OER could simultaneously generate hydrogen and oxygen that achieve overall water splitting and help further understand the performance trend of overall electrocatalysts [12–15]. For the metal-air batteries that require the combination of OER and ORR, using bifunctional electrocatalysts could exhibit better performance than two separate electrocatalysts owing to the best working conditions of two separate electrocatalysts being different [16,17]. Moreover, the use of bifunctional electrocatalysts could lower the cost and simplify the procedures due to their working conditions are the same. Hence, it is of great significance to develop bifunctional electrocatalysts for the HER, OER and ORR. In addition to bifunctional electrocatalysts, single-atom catalyst (SAC) is another kind of catalyst that possesses high performance and low cost, owing to the exposed homogeneous metal active sites and maximum use of metal atoms [18,19]. Given this, the rational design of single-atom bifunctional electrocatalysts for the HER, OER and ORR is of paramount importance and economic interest.

Two-dimensional (2D) covalent organic frameworks (COF) are a class of layer-stacked materials consisting of uniformly exposed metal sites and covalent organic ligands [20,21]. 2D COF materials possess the large surface area, higher efficiency in using metal atoms than the counterpart bulk structures and tunable chemical functionality, which have attracted widespread attention in electrocatalysis [22]. Recently, a new stable dioxin-linked metallophthalocyanine 2D COF was successfully synthesized with the hosted metal atom coordinated by four in-plane nitrogen atoms (namely MN₄@COF) [23]. In this work, we choose Fe as the metal active site as it has been widely studied in the metal-nitrogen-carbon (M-N-C) catalysts for the ORR and is recognized as the most possible alternative to noble-metal-based catalysts [24,25]. However, these Fe-N-C catalysts generally exhibit inferior HER and OER activity. To achieve the goal of designing single-atom bifunctional catalysts for overall water splitting, fuel cells and metal-air batteries, the catalytic activity of these electrochemical reactions need to be improved. It should be pointed out that the axial coordination of the metal active site has become a promising frontier to tailor the electronic structures and activity of M-N-C catalysts [26–28]. Hence, a series of FeN₄@COF and FeN₄-X@COF (X represents the axially coordinated ligands, X = -SCN, -I, -H, -SH, -NO₂, -Br, -ClO, -Cl, -HCO₃, -NO, -ClO₂, -OH, -CN and -F) were constructed, and their HER, OER and ORR catalytic activity were systematically investigated by considering the effect of axial coordination through density functional theory (DFT) calculations. The calculated results suggest that all the designed catalysts are thermodynamic and electrochemical stable, and the axial ligands could enhance the corresponding activity of FeN₄@COF. FeN₄-SH@COF and FeN₄—OH@COF are screened out as efficient bifunctional catalysts for HER and OER, FeN₄—OH@COF and FeN₄-I@COF are promising bifunctional catalysts for OER and ORR. Moreover, given the tremendous speed of progress in the experimental work in this field, it is highly desired that the theoretical results in our work can inspire experimental groups to synthesize the best catalysts identified from our calculations.

The schematic models and optimized configurations of FeN₄@COF and FeN₄-X@COF with the axial ligand (-SCN, -I, -H, -SH, -NO₂, -Br, -ClO, -Cl, -HCO₃, -NO, -ClO₂, -OH, -CN and -F) coordinated to the Fe active center are shown in Figs. 1a and b and Fig. S1 (Supporting information). It should be pointed out that one of the critical issues in designing efficient electrocatalysts is their structural stability during long-term use. Therefore, we calculated the corresponding formation energy (E_f) and dissolution potential (U_diss) of all the considered catalysts to theoretically evaluate their thermodynamic and electrochemical stabilities [29,30]. The E_f and U_diss are calculated by the following equations: E_f = E_total - E_substrate - μ_Fe - μ_ligand and U_diss = U_diss°(bulk) - E_f/ne, respectively, where E_total and E_substrate refer to the total energies of catalyst and substrate, respectively. μ_Fe is the energy of a Fe metal atom in its most stable bulk structure, and μ_ligand is the coordinating axial ligand, which uses the corresponding hydrides and balances with H⁺ + e^– which equals the energy of 1/2 H₂ under the standard conditions. U_diss°(bulk) is the standard dissolution potential of the Fe bulk metal (U_diss°(Fe, bulk) = −0.45 V), n is the number of electrons involved during the dissolution process (n = 2) [29]. Since μ_Fe is referenced with respect to its bulk metal, the designed catalysts with negative values of E_f are evaluated to be thermodynamically stable against the clustering of Fe atoms. The designed catalysts with positive values of U_diss vs. standard hydrogen electrode (SHE) are considered to be electrochemically stable. The calculated corresponding results of E_f and U_diss are displayed in Figs. 1c and d, which indicates that all the designed catalysts possess good thermodynamic and electrochemical stabilities that could meet the stability criteria for the electrocatalysts.

The distinct electronic properties of catalysts were explored to gain deep insight into their electrocatalytic performance. As displayed in Fig. S2 (Supporting information), all the designed catalysts could ensure efficient electron transfer during the electrochemical reaction process. Moreover, the calculated partial density of states (PDOS) results (Fig. S3 in Supporting information) suggest that the different axial ligands of the designed catalysts exhibit different contribution to the electronic states around the Fermi level and the electronic states around the Fermi level are mainly contributed by the d orbitals of the Fe atoms. Additionally, the hybridization between the d orbitals of the Fe atoms and the p orbitals of the axially coordinated non-metal atoms of the ligands further demonstrate their strong interaction. The Bader charge analysis (Fig. S4 in Supporting information) indicates that a large amount of charge (1.03–1.37 e) transfer from the Fe atoms to other moieties of the designed catalysts, making the Fe atoms positively charged. These positively charged Fe atoms of the catalysts are considered to be active sites in the following electrochemical reaction processes. Obviously, the Fe active sites possess comparable or more positive charge (1.15–1.37 e) in FeN₄-X@COF (X = -NO₂, -Br, -ClO, -Cl, -HCO₃, -NO, -ClO₂, -OH, -CN and -F) than the Fe in FeN₄ (1.14 e), suggesting that these axially coordinated ligands serve as the electron-acceptor which could attract electrons from the Fe atom. On the other hand, the Fe sites coordinated by -SCN (1.10 e), -I (1.03 e), -H (1.05 e) and -SH (1.07 e) carry less positive charge compared to the FeN₄ (1.14 e), indicating that these ligands function as the electron-donor and could donate electrons to the Fe atom. The induced changes in the charge of the Fe active sites indicate that the electronic properties are significantly regulated by the axial coordination ligands, and the catalytic activity of HER, OER and ORR would be affected thereafter.

The HER performance of the designed FeN₄@COF and FeN₄-X@COF catalysts were subsequently investigated. The Gibbs free energy of adsorbed hydrogen on catalyst (ΔG_H^*) is a key descriptor for evaluating the HER activity of the catalyst [31]. For an ideal HER catalyst, the value of ΔG_H^* should be close to zero. In detail, the positive ΔG_H^* value indicates that the adsorption of hydrogen on the catalyst is kinetically unfavorable, while, the negative ΔG_H^* value suggests that the adsorbed hydrogen on the catalyst is difficult to release, which inhibits the HER activity. The calculated ΔG_H^* value of all the designed catalysts are displayed in (Fig. 2a). The exchange current (i₀) that can reveal the intrinsic rate of electrons in the HER process under equilibrium potential were also studied [31]. The theoretical exchange current was calculated by the equation: , here k₀, k_b and T refers to the reaction rate constant at zero overpotential, the Boltzmann constant and the temperature, respectively. For illustrative purpose, the k₀ value was set to 1. Therefore, as shown in Fig. 2b, a volcano curve can be plotted by using i₀ as a function of the calculated ΔG_H^* to compare the HER performance on the designed stable catalysts. Importantly, the catalytic activity of the HER can be quantitatively evaluated from the position of the calculated ΔG_H^* and i₀ with respect to the plotted volcano peak. The closer the location of the calculated ΔG_H^* is to the volcano peak, the better the HER catalytic activity of one catalyst [31]. Obviously, from Fig. 2a we can see that the calculated ΔG_H^* values of FeN₄-SH@COF, FeN₄—OH@COF and FeN₄—ClO@COF are close to the optimal value zero, thus, their location are around the peak of the plotted volcano curve with large exchange current (Fig. 2b), which suggests that they could serve as promising HER catalysts.

The OER and ORR catalytic activity of all the designed catalysts were systematically studied. The Gibbs free energy of intermediates govern the intrinsic potential-determining steps of the OER and ORR which correspondingly determine the distinct catalytic activity of a catalyst as proposed by Nørskov et al. [32]. The calculated Gibbs free energies of corresponding intermediates (ΔG_HO^*, ΔG_O^* and ΔG_HOO^*) are exhibited in Fig. S5 (Supporting information). According to the Sabatier principle [31], the adsorption energy of the adsorbates should be at the optimum value, not too strong nor too weak, for the catalyst to facilitate catalytic performance. Thus, identifying efficient OER/ORR catalysts with moderate interaction strength is vital in the following study. The Gibbs free energy difference between the two adjacent intermediates of an ideal OER/ORR (four-electron transfer steps) catalyst should be equated to 1.23 eV (U = 0), making the OER/ORR occur at the thermodynamic limit with the overpotential η of zero. It should be pointed out that the above free energy difference of a catalyst is usually unequal. The overpotential of OER (η^OER) and ORR (η^ORR) is determined by the maximum and minimum Gibbs free energy difference, respectively. The corresponding Gibbs free energy diagrams of the OER and ORR on all the designed catalysts are calculated and shown in Fig. S5. The potential-determining step for each catalyst is colored in blue (for the OER) and green (for the ORR), and the blue and green values are the potential-determining step values for OER and ORR. Moreover, the calculated OER and ORR overpotentials are summarized in Fig. 3a. Obviously, compared with the FeN₄@COF catalyst, the axially coordinated FeN₄-X@COF catalysts exhibit enhanced OER performance, and most of them show comparable to or even better than that of IrO₂ (110) (η^OER = 0.52 V) except FeN₄—NO@COF and FeN₄—ClO₂@COF [33]. Especially, among all the designed catalysts, the FeN₄—CN@COF exhibits the best OER activity coupled with a low η^OER value of 0.28 V, followed by FeN₄-SH@COF (η^OER = 0.33 V), FeN₄—OH@COF (η^OER = 0.33 V) and FeN₄—NO₂@COF (η^OER = 0.38 V), suggesting their promising OER activity. As the reverse reaction of the OER, the Gibbs free energy diagrams of the ORR are also presented in Fig. S5. Compared with the FeN₄@COF catalyst, the designed FeN₄-X@COF catalysts also exhibit enhanced ORR catalytic activity, most of which are comparable to or even lower ORR overpotentials than that of Pt (111) (η^ORR = 0.48 V) except FeN₄—H@COF (η^ORR = 0.82 V), FeN₄—NO₂@COF (η^ORR = 0.59 V), FeN₄@COF (η^ORR = 0.58 V) and FeN₄-SH@COF (η^ORR = 0.55 V) [34]. Additionally, from Fig. 3a, it can be seen that FeN₄—ClO₂@COF is the best ORR catalyst with the calculated η^ORR value of 0.32 V, followed by FeN₄-I@COF (η^ORR = 0.33 V), FeN₄-Br@COF (η^ORR = 0.35 V), FeN₄—OH@COF (η^ORR = 0.37 V), FeN₄—Cl@COF (η^ORR = 0.39 V), and FeN₄—HCO₃@COF (η^ORR = 0.40 V), and all of them exhibit lower ORR overpotentials than that of Pt (111), indicating their promising ORR activity. Hence, based on the above results, it can be concluded that the most of the considered axially coordinated ligands X could improve the OER/ORR activity of FeN₄@COF catalysts to approach or exceed the OER/ORR performance of IrO₂ (110)/Pt (111). In particular, the axially coordinated FeN₄-X@COF catalysts are the promising OER and ORR bifunctional catalysts because of the correspondingly synchronously decreased OER and ORR overpotentials. As reported in literature [28,35,36], the overvoltage of the OER and ORR (Δη = η^OER + η^ORR) reveals the efficiency loss and is commonly used to evaluate the OER/ORR bifunctional catalyst. A bifunctional OER/ORR catalyst with a low Δη value can serve as a highly active reversible oxygen electrode. Among all the considered FeN₄-X@COF catalysts, the FeN₄—OH@COF catalyst possesses the lowest Δη value of 0.70 V, followed by FeN₄-I@COF (Δη = 0.78 V), FeN₄—CN@COF (Δη = 0.78 V), FeN₄-Br@COF (Δη = 0.81 V), FeN₄—Cl@COF (Δη = 0.84 V), FeN₄—HCO₃@COF (Δη = 0.87 V), FeN₄-SH@COF (Δη = 0.88 V) and FeN₄-F@COF (Δη = 0.90 V), lower than that of the combination of IrO₂ (110) and Pt (111). Interestingly, the best-performing OER/ORR bifunctional catalyst FeN₄—OH@COF exhibits better activity than that of the reported penta-coordinated Fe-N-C catalysts [28]. Moreover, the corresponding optimized structures of the adsorption intermediates on the efficient FeN₄—OH@COF and FeN₄-I@COF catalysts are shown in Fig. S6 (Supporting information).

An in-depth understanding of the distinct OER and ORR catalytic activities of the studied FeN₄@COF and FeN₄-X@COF catalysts can guild us in designing promising catalysts. As aforementioned, the activity of OER and ORR is governed by the Gibbs free energies of the intermediates on the catalyst. Therefore, it is necessary to establish relationship between the Gibbs free energy of the intermediates and the OER/ORR activity for the rational design of potential catalysts. Fig. 3b shows the comparison of the calculated adsorption Gibbs free energy values of the HO* and HOO* intermediates on the designed catalysts. Obviously, the value of ΔG_HO^* and ΔG_HOO^* are highly correlated and exist a strong linear scaling relationship which mainly originates from their formation of single bond between the Fe and O atoms. ΔG_HOO^* can be expressed as a function of ΔG_HO^* via the equation ΔG_HOO^* = 0.83ΔG_HO^* + 3.21 eV, such linear relationship is consistent with the reported carbon-based catalysts for the OER and ORR [37]. As displayed in the free energy diagrams of the OER (Fig. S5), the potential-determining step generally occurs at the formation of O* intermediate from HO* or the formation of HOO* from O* step, the OER overpotential could be described by using the difference of ΔG_O^*- ΔG_HO^*. This is true confirmed by the volcano curve illustrated in Fig. 3c, where the catalytic activity of all the designed FeN₄@COF and FeN₄-X@COF catalysts toward the OER fall in a line as a function of ΔG_O^* - ΔG_HO^*. Additionally, all of the FeN₄-X@COF catalysts are located at the right side of the FeN₄@COF catalyst, indicating their relatively weaker adsorption of O* intermediate, therefore improving the OER activity. Apparently, the FeN₄—CN@COF, FeN₄-SH@COF and FeN₄—OH@COF catalysts are located around the top of the plotted volcano curve with low overpotentials and stand out to be efficient OER catalysts. For the ORR, the potential- determining mainly occurs at either the first step (from * to HO*) or at the last step (from HOO* to * + O₂), the overpotential could be described by using the ΔG_HO^*. Fig. 3d displays the volcano curve of the ORR overpotential as a function of ΔG_HO^*. It should be noted that the catalysts located at the left leg of the volcano plot suggest the strong HO* adsorption with their performance mainly limited by the step of HO* reduction, and those located at the right leg indicate weak adsorption of the HO* intermediate which could facilitate the HO* desorption and ORR activity. Among all the designed catalysts, FeN₄—ClO₂@COF, FeN₄-I@COF and FeN₄-Br@COF catalysts locate near the apex of the volcano curve, which further demonstrates the above-mentioned high activity. Therefore, it can be concluded that the axial coordinated ligands could tune the intermediate adsorption ability of FeN₄@COF catalysts to the moderate and accordingly enhance the OER and ORR activity. Importantly, FeN₄—OH@COF and FeN₄-I@COF catalysts are both located around the peak of plotted volcano curves and could serve as efficient bifunctional catalysts for the OER and ORR. Note that the four-electron pathway of the ORR from O₂ to H₂O plays an important role in the metal-air batteries, while the two-electron pathway from O₂ to H₂O₂ is undesirable as it may lead to the degradation of the catalyst [38]. For the aforementioned screened promising ORR catalysts that with the overpotential lower than 0.45 V (FeN₄—ClO₂@COF, FeN₄-I@COF, FeN₄-Br@COF, FeN₄—OH@COF, FeN₄—Cl@COF, FeN₄—HCO₃@COF and FeN₄-F@COF), the selectivity for the four-electron pathway is evaluated from the thermodynamic perspective, since the free energies of O* intermediate (2.02, 2.34, 2.24, 2.38, 2.22, 2.16 and 2.44 eV) are lower than 3.52 eV (ΔG_H - ΔG_H) [39,40]. A kinetics studied was further carried out to understand the catalytic selectivity during the ORR process [39,41] and the computation details are shown in Supporting information. Fig. 3e illustrates that all the calculated values of ln(k_O/k_H) are positive, suggesting that the formation of H₂O from O₂ is prioritized on these designed catalysts. The calculated values of ln(k_sys/k_Pt(111)) of the designed FeN₄—ClO₂@COF, FeN₄-I@COF, FeN₄-Br@COF, FeN₄—OH@COF, FeN₄—Cl@COF, FeN₄—HCO₃@COF and FeN₄-F@COF are 8.51, 8.12, 7.35, 6.58, 5.80, 5.42 and 3.86, respectively, illustrating that the reaction rate of these catalysts is faster than that on the Pt (111). Furthermore, the correlation between the OER/ORR catalytic activities and the electronic properties of catalysts was studied. As shown in Fig. S7 (Supporting information), it obviously exhibits volcano-shaped relationship between the Bader charge of the Fe active metal atoms and the OER/ORR overpotential. After axially coordinated by the considered X ligands, both OER and ORR overpotentials are decreased firstly and then increased within the Bader charge changes from 1.03 e to 1.37 e, suggesting that the FeN₄-based catalysts with high activity requires moderate Bader charge of the active site, which still follows the Sabatier principle. Therefore, the Bader charge of the metal atoms can serve as a feasible descriptor to describe the axial coordination effect on the OER and ORR activities of the FeN₄@COF catalysts. Additionally, the above results further demonstrate that the OER and ORR activity of the FeN₄-based catalysts could be modulated using axial coordinated ligands.

The thermodynamic calculations screen out the above-mentioned promising OER and ORR activity. To consider kinetic factors that may be neglected in the above thermodynamic calculations, we constructed a microkinetic model to gain a deeper understanding of the catalytic activity of these catalysts. To visualize the OER and ORR performance of the above-screened catalysts, their corresponding theoretical polarization curves were simulated based on the reversible hydrogen electrode (RHE) and plotted in Fig. 3f with the change in the current density as a function of potential U. The simulated details are shown in the Supporting Information. As displayed in Fig. 3f, for the OER, the simulated polarization curves of FeN₄—CN@COF, FeN₄-SH@COF, FeN₄—OH@COF, FeN₄—NO₂@COF, FeN₄-I@COF, FeN₄—Cl@COF, FeN₄-Br@COF, FeN₄-F@COF and FeN₄—HCO₃@COF exhibit lower onset potentials (at the current density of 10 mA/cm²) of 1.57, 1.68, 1.77, 1.92, 1.92, 1.94, 1.94 and 1.96 V vs. RHE than that of the IrO₂ (110) (2.27 V) catalyst as reported in our previous work [42], respectively, suggesting their high OER activity. For the ORR, the simulated polarization curves of FeN₄—ClO₂@COF, FeN₄-I@COF, FeN₄-Br@COF, FeN₄—OH@COF, FeN₄—Cl@COF, FeN₄—HCO₃@COF and FeN₄-F@COF show higher onset potentials (at the current density of 1 mA/cm²) of 1.03, 1.00, 0.96, 0.93, 0.88, 0.86 and 0.78 V vs. RHE than that of the Pt (111) (0.67 V) catalyst as reported in our previous work [42], respectively, demonstrating their high ORR activity; while the FeN₄—CN@COF, FeN₄-SH@COF and FeN₄—NO₂@COF show lower onset potentials (at the current density of 1 mA/cm²) of 0.66, 0.57 and 0.48 V vs. RHE than that of the Pt (111), respectively. Based on the above thermodynamic and microkinetic results demonstrate that these screened out catalysts are predicted to be promising catalysts for the OER and ORR; the bifunctional catalytic activity of FeN₄-I@COF, FeN₄—CN@COF and FeN₄-Br@COF toward both the OER and ORR is also excellent. Moreover, the FeN₄-SH@COF and FeN₄—OH@COF can serve as efficient bifunctional catalysts for both HER and OER that achieve the overall sustainable water splitting.

In summary, by means of DFT calculations, we systematically studied the HER, OER and ORR catalytic performance of FeN₄@COF electrocatalysts by considering the effect of axial coordination ligands X (X = -SCN, -I, -H, -SH, -NO₂, -Br, -ClO, -Cl, -HCO₃, -NO, -ClO₂, -OH, -CN and -F). The calculated results demonstrate that all the designed catalysts exhibit good thermodynamic and electrochemical stabilities. The change of axially coordinated ligands could modulate the interaction strength between intermediates and catalysts that governs the corresponding catalytic activities. Among the investigated catalysts, the designed FeN₄-SH@COF and FeN₄—OH@COF are proved to be efficient bifunctional catalysts for both HER and OER with calculated ΔG_H^* and η^OER of 0.002 eV, 0.03 eV, 0.33 V and 0.33 V, respectively. FeN₄—OH@COF, FeN₄-I@COF and FeN₄-Br@COF are expected to be promising bifunctional catalysts toward both OER and ORR with calculated overpotential of 0.33, 0.45, 0.46, 0.37, 0.33 and 0.35 V, respectively. The FeN₄-X@COF catalysts generally exhibit enhanced activity toward the OER and ORR compared with that of FeN₄@COF catalysts. Moreover, the Bader charge of catalysts can serve as a feasible activity descriptor for the OER and ORR. This study offers insight into modulating the catalytic activity by control of the axially coordinated ligands and provides theoretical guidance for optimizing iron-nitrogen-carbon catalysts.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

This work was supported by the National Natural Science Foundation of China (Nos. 22102167 and U21A20317). The calculations were performed on the Supercomputing Center of the University of Science and Technology of China and Ningbo Artificial Intelligence & High Performance Computing Center.

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.cclet.2024.109588.