HIGHLIGHTS
• Surface energy and surface Pourbaix diagram reveal that ZrO2($\overline{1} 11$) is the most thermodynamically stable facet and is preferentially occupied by HO* at the equilibrium potential of oxygen evolution reaction (OER).
• Microkinetic modeling analyzed the OER activity of 40 single-metal doped ZrO2 and identified 16 metals exhibit improved catalytic activity, with Rh and Fe dopants showing the remarkable improvement.
• Thermodynamic free energy diagrams, density of states analysis, and ab initio molecular dynamics simulations further confirm that Fe-ZrO2 and Rh-ZrO2 are highly promising catalysts for OER, showcasing low ΔG for the rate-determining step, high conductivity, and exceptional stability.
1 Introduction
Proton exchange membrane water electrolysis offers an attractive means to store and convert renewable energy into sustainable hydrogen [1,2,3]. Water electrolysis comprises two half-reactions: oxygen evolution reaction (OER) and hydrogen evolution reaction (HER). Extensive studies have displayed that OER on the anode suffers from sluggish kinetics and poor stability resulting in a significant efficiency loss, thereby becoming the main obstacle for the practical implementation of electrochemical water splitting [4,5,6,7,8,9,10,11]. Therefore, it is vital to rationally search and design OER electrocatalysts with high activity and durability. Moreover, maintaining a delicate equilibrium between stability and activity poses a significant challenge for various electrocatalysts under OER operating conditions.
Zirconium dioxide (ZrO2), a highly sustainable material, has excellent thermal stability, making it highly desirable in various applications. For instance, nano-ZrO2 displayed remarkable aromatic selectivity and stability during a 120-h test [12]. ZrO2 also boasts exceptional hydrothermal stability during the CO2 methanation process under high temperatures and pressures [13]. Notwithstanding this, the OER activity of ZrO2 was known to be notably low [14,15,16]. Doping, a frequently utilized and remarkably efficient modulation strategy, can enhance the electrical conductivity and catalytic performance of materials [17,18,19,20,21]. For example, Mishra et al. demonstrated that optimal cobalt doping (5%) in CuO improved OER activity significantly with a low overpotential of 120 mV and a small charge transfer resistance of 2.58 Ω [22]. Additionally, Fu and co-workers demonstrated that Ce dopant induced electron redistribution of CoO to stabilize Co-O bonds and displayed an optimal binding strength with OER intermediates, manifesting favorable activity with an overpotential of 261 mV [23]. Therefore, it is possible to leverage the inherent stability of ZrO2 to conduct single-atom doping research, aiming to ensure its stability while enhancing its activity.
Herein, we analyzed the electronic structure and OER performance of ZrO2 by doping a series of metals with flexible valence states, based on spin-polarized density functional theory calculations with van der Waals corrections (DFT-D3). First, we identified the thermodynamically stable ZrO2 surface, the ($\overline{1} 11$) surface, based on surface energy calculations with benchmarking with experimental literature. Subsequently, we analyzed the surface state of the identified stable ZrO2 surface via surface Pourbaix diagram calculations and found that the ZrO2 surface tends to be occupied by HO* under solution conditions at OER equilibrium potential. Finally, via a stratified catalyst design and screening strategy with various potential dopants for ZrO2 (Fig. 1a), we identified Fe as a promising dopant element by considering the effects of potential radiation, abundance, stability, activity, and cost-effectiveness. This work provides experimentalists with potentially superior catalysts worthy of future validation.
Fig. 1 a Proposed strategy for screening metals capable of enhancing the OER activity of ZrO2 via doping. b Calculated surface energies of the three low-index slabs of monoclinic ZrO2. More details can be found in Fig. S1 |
2 Computational Details
The Vienna ab initio simulation package (VASP) was used to perform all spin-polarized DFT calculations [24,25]. The generalized gradient approximation (GGA) in the form of revised Perdew-Burke-Ernzerhof (RPBE) [26] was employed to describe the exchange-correlation potentials. The Hubbard U parameter, representing the strong correlation repulsion energy between electrons with opposite spins, was added to the Fe (U = 4 eV) [27,28], Co (U = 3 eV) [29,30], Ni (U = 6.6 eV) [31], and Mn (U = 4 eV) [27,32] 3d electrons. The electron-ion interaction was described with the projector-augmented wave method [33]. A plane-wave basis set of 480 eV was adopted. ZrO2 generally exists in three polymorphs [34,35,36]: monoclinic at room temperature, tetragonal at 1480-2650 K, and cubic at > 2650 K. Because the high thermodynamic stability is compulsory for conventional electrocatalysis under room temperature, herein, we further discuss the ZrO2 in a monoclinic form. Initial magnetic moments of 5, 0.6, 5, and 5 are set for the magnetic materials Fe, Co, Ni, and Mn, respectively. The zero-damping DFT-D3 method of Grimme was used to correct the van der Waals interactions [37]. The convergence criterion was set to be lower than 0.02 eV Å−1 for the force and 10-6 eV per atom for energy. A 2 × 2 × 1 Monkhorst-Pack was used for k-point sampling to describe the supercell [38]. A vacuum layer of 18 Å was set in the z-direction to avoid the interaction between periodic images. We employed ab initio molecular dynamics (AIMD) simulations in a canonical NVT ensemble to assess the thermodynamic stability of structures. The whole simulation lasted for 10 ps with a time step of 2 fs at 300 K. Furthermore, crystal orbital Hamilton population (COHP) was also calculated to analyze the local chemical bond properties in periodic systems. COHP is obtained via multiplying the Hamiltonian matrix by the corresponding DOS matrix. Combining it with the powerful LOBSTER tool [39], capable of handling plane-wave basis sets, allows for the visualization of COHP diagrams and analysis of bond strength.
The changes in Gibbs free energies for each elementary step were calculated using the following equation [40]:
$\Delta G=\Delta E+\Delta \mathrm{ZPE}-T \Delta S+\int C_{p} d T+\Delta G_{p H}+\Delta G_{U}$
Where ΔE, ΔZPE,ΔS, and Cp are the changes in the electronic energy directly obtained from DFT, zero-point energy, entropy, and heat capacity, respectively. T is the temperature (298.15 K). The entropies of molecules in the gas phase were acquired from the NIST database. ΔGpH is the free energy contribution due to the variations in H concentration, expressed by the Nernst equation as ΔGpH=2.303kBTpH. In this work, we set pH to zero. ΔGU=−neU, where n was the electron transfer number and U was the applied potential.
Under the acidic condition, the thermodynamic potential for the oxidation of H2O to produce O2 involving 4 electron transfers is described as:
$2 \mathrm{H}_{2} \mathrm{O} \rightarrow \mathrm{O}_{2}+4 \mathrm{H}^{+}+4 e^{-} $
The thermodynamic overpotential (η) is the important indicator to judge the catalytic activities of a catalyst, which can be obtained by using the following equation:
$\eta=\frac{\Delta G_{\max }}{e}-1.23 \mathrm{~V} $
Where ΔGmax is the maximum free energy change among the four elementary steps, and 1.23 V is the equilibrium potential. A lower η implies higher catalytic activity, and the η of an ideal catalyst is zero.
$\gamma_{s}=\frac{\left[E_{{slab-unrelax }}-N E_{{bulk }}\right]}{2 A}+\frac{\left[E_{ {slab-relax }}-E_{ {slab-unrelax }}\right]}{A}$
where Eslab−unrelax and Eslab−relax is the total energy of the unrelaxed and relaxed slab, respectively. Ebulk is the bulk energy per atom, N is the number of atoms in the slab, and A is the surface area. The kinetic volcano model was derived based on the potential-dependent kinetic information as a function of current density and adsorption free energies of HO* and O*, with the details and parameters shown in Ref. [43].
$H_{n} O_{m}^{*}+(2 m-n)\left(e^{-}+H^{+}\right) \rightleftarrows *+m H_{2} O $
where m and n are the number of oxygen and hydrogen atoms of the adsorbate, respectively. The differences in free energy were calculated using:
$\begin{aligned} \Delta G= & G_{ {bare }}+m G_{H_{2} \mathrm{O}}-G_{ {total }}-(2 m-n) \\ & \left(\frac{1}{2} H_{2}-U_{S H E}-2.303 k_{B} T * p H\right) \end{aligned}$
Here, USHE refers to the potential relative to the standard hydrogen electrode (SHE), and kB is the Boltzmann constant (8.617343 × 10-5 eV K−1). The surface Pourbaix diagrams were modeled with the CatMath online platform through cloud computing [47].
3 Results and Discussion
3.1 Surface Energy and Surface Pourbaix Analysis
The most stable phase of ZrO2 at low temperature is monoclinic [36] taking the P21/c space group, and the optimized lattice parameters (a = 5.15 Å, b = 5.26 Å, c = 5.30 Å) closely match the experimental values [48]. To determine the most stable exposed facet, we conducted surface energy calculations on three representative low-index planes, namely ($\overline{1} 11$), (111), and (101), all of which are intricate models featuring O-terminated, Zr-terminated, and Zr and O co-terminated interfaces. Figure S1 provides a visual representation of all possible surface structures maintaining a stoichiometric ratio consistent with the bulk, along with their corresponding surface energy values. A lower surface energy indicates a higher propensity for the formation of chemical bonds, thus suggesting a thermodynamic preference for crystal growth along that specific plane. The minimum surface energy values for ($\overline{1} 11$), (111), and (101) planes are 0.725, 0.930, and 1.142 J m−2, respectively (Fig. 1b). This implies that ZrO2 ($\overline{1} 11$) is the most thermodynamically stable exposed facet, aligning consistently with earlier theoretical prediction and experimental observation [49,50]. This finding identifies an important substrate model, ZrO2 ($\overline{1} 11$), which is a prerequisite for OER activity analysis.
As Fig. 2a illustrated, O-terminated ZrO2 configuration holds many various adsorption sites. H*, O*, or HO* species generated through water activation are possible to pre-covered on the surface of catalysts under OER operating conditions, which renders the oxide system may deviate from its pristine stoichiometric form [43] and exerts substantial influences on reaction overpotential. Surface Pourbaix analysis can provide a valuable approximation of the real state under OER conditions [44,45,46]; thus, we calculated the surface Pourbaix diagram (Fig. 2b). The results illustrate that the surface of ZrO2 will be initially occupied by HO* at USHE > 1.13 V vs. reversible hydrogen electrode (RHE), which is lower than the equilibrium potential of OER (1.23 V vs. RHE); that means, under working conditions, the surface is no longer in its pristine state but instead adopts a restructured configuration enriched with HO*. However, given that HO* also serves as the initial adsorption intermediate of OER, we can directly use the original model to analyze the OER activity. Surface Pourbaix diagram provides a detailed analysis of the surface occupation states to help us obtain a theoretical substrate model that is closer to the experiment. To assess the viability of doping various metals onto the ZrO2 surface, we calculated the doping formation energy (Fig. 2c). The most stable doping configurations for each metal are shown in Figs. S2 and S3. Generally, the doping formation energy is calculated to assess the thermodynamic stability of catalysts, wherein the negative value indicates the high stability and vice versa. As Fig. 2c shows, Cu, Ag, Au, and Hg dopants exhibit positive formation energies, which means these dopants are difficult to replace Zr and fix on the ZrO2 stably. Conversely, the remaining metal dopants display negative formation energies, indicating they are thermodynamical feasibility to substitute Zr and form stable catalysts. All feasible metal elements for doping are listed in Fig. 2d, which are subsequently investigated for the OER analysis. Specifically, metals are arranged in the order of the periodic table of elements except Al. The yellow, pink, blue, green, and purple backgrounds correspond to elements from the third, fourth, fifth, sixth periods, and rare earth elements, respectively. This order can identify elements in empty positions, which are subsequently disregarded.
Fig. 2 a Optimized configuration of ZrO2 with top and side views. The purple and red spheres represent Zr and O, respectively. b Calculated surface Pourbaix diagram of the ZrO2 surface. c Formation energy values of different single-metal atoms doped ZrO2. d The feasible metals for doping on the surface of ZrO2. |
3.2 OER Activity Screening Based on the Volcano Activity Model
The GO*-GHO* is generally regarded as a crucial descriptor for evaluating OER activity since it quantifies the capacity of surface oxygen atoms to make and break bonds with hydrogen and oxygen atoms, which is prone to be the most thermodynamically challenging reaction steps for OER [2,43,51]. To achieve more accurate predictions, we employed the microkinetic modeling as the function of GO*-GHO* and potential at various currents [43,52]. The current curves were obtained by numerically solving the rate of HOO formation considering the O→ HOO transition state. Figure 3a illustrates the kinetic OER activity relating potential to GO*-GHO* at a current density of 5 μ A cm−2, featuring the characteristic volcano-shape. On the left side of the peak, decreased GO*-GHO* values correspond to lower OER kinetic activity. Similarly, catalysts exhibit the decaying OER performance along the right arm of the volcano curve. Catalysts close to volcanic peaks, such as Co, Mn, Rh, Fe, and Pt, demonstrate good catalytic activity. Especially, Rh sitting at the volcano peak, with an optimal GO*-GHO* value, requires relatively low overpotentials to achieve high reaction rates in acidic solutions, demonstrating the highest OER activity. In addition, we plot potential as a function of GO*-GHO* at 1 mA cm−2 (Fig. S4) and 1 A cm−2 (Fig. S5), and the results display that the activity order may change as the current increases, but the qualitative order remains constant; that is, a good catalyst is always good at different currents. Figure 3b displays a linear scaling relation between the adsorption energy of O* (EO*) and HO* (EHO*), with the identified slope and intercept close to those of higher-index (where h2 + k2 + l2 > 1) metal oxide surfaces [52]. Figure 3c further demonstrates that the adsorption energy of HOO* (EHOO*) and EHO* exists a good linear correlation and a constant difference between EHOO* and EHO* on effectively doped catalysts, which implies GHOO* – GHO* = (GO* – GHO*) + (GHOO* – GO*) = constant, i.e., GHOO* – GHO* = (GO* – GHO*) + [constant(GO* – GHO*)]. Therefore, GO* – GHO* serves as a distinctive descriptor elucidating the OER activity, further reinforcing the rationale behind the selection of Fig. 3a. Furthermore, EHO* and $\varepsilon_{d^{-}} \text {-up } / \varepsilon_{d} \text {-down }$ demonstrate a robust linear correlation, as depicted in Fig. 3d. Such correlation offers valuable insights into finding optimal HO binding energies by manipulating the $\varepsilon_{d}$ value. Generally, the optimal OER catalyst has surface hydroxide deprotonation free energy value of 1.5 to 1.7 eV [43,53,54,55]. In our work, taking 1.5 eV as the standard fits well with the dynamic trend (Fig. 3a), we use (GO* − GHO* = 1.5 eV) as the descriptor to further embody catalytic behavior. A smaller absolute value of (GO* − GHO* = 1.5 eV), represented by a darker color, suggests better OER performance (Fig. 3e). Figure 3e selects all metals closing to the volcanic peak exhibits improved OER performance compared to the pristine ZrO2, in which Rh has the best catalytic activity with values closest to zero. Theoretical overpotential (η), a fundamental thermodynamic parameter, has emerged as a crucial metric for evaluating and predicting the OER performance of catalysts. Thus, we plot η as a function of GO* − GHO*, which has a volcano relationship (Fig. S6). Rh standing on the peak of the volcano map exhibits the highest OER activity with minimum η, followed by Fe, Pt, etc. This thermodynamic volcano diagram trend is roughly consistent with the dynamic volcano diagram (Fig. 3a), and the gap may be attributed to the omission of transition states in thermodynamics. To acquire a more profound comprehension of the impact of doping on the electronic properties of the catalyst, we calculate the density of states (DOS) of ZrO2 (Fig. S7a). Near the Fermi level, the valence band is predominantly contributed by O orbitals, while the conduction band is primarily composed of Zr orbitals. A wide forbidden region exists between the top of the valence band and the bottom of the conduction band with a bandgap of 2.6 eV, further observed by the band structure diagram (Fig. S7c). Also, the precise HSE functional displayed that ZrO2 is a semiconductor material with a wide bandgap of 4.45 eV (Fig. S7b). Figures S8 and S9 show the DOS of metal doped ZrO2. The bandgap values exhibit a progressively decreasing trend in the following order: Ir (2.54 eV), Ru (2.49 eV), Sn (2.28 eV), Ti (1.82 eV), and Cr (1.71 eV). In, Sm, Pr, Fe, Nd, Ni, Pd, Ga, Co, Pt, Mn, and Rh metals, on the other hand, induce a direct transformation of the catalyst from a semiconductor to a semi-metallic property characterized by a band gap of 0 eV, effectively promoting the electron transfer and thereby improving the conductivity of the ZrO2. In short, the microkinetic volcano diagrams screen out some promising catalysts, the linear relationship obtains the abnormal adsorption points and the correlation between different adsorbates, and DOS uncover enhanced conductivity of catalysts through doping.
Fig. 3 a Kinetic OER activity volcano model as a function of GO*-GHO* at 5 μ A cm−2 (black line). b Scaling relation between EO* and EHO* on feasible metals doped ZrO2. The linear scaling relations of EHO* vs. c EHOO* and d $\varepsilon_{d}$ on elevated metal doped ZrO2. e A heatmap represents the magnitude of values as a color. Values (unit eV) were obtained by GO*-GHO* = 1.5, and darker colors correspond to lower values implying higher activity. |
3.3 Free Energy Analysis
Gibbs free energy diagrams were established to visualize the OER intermediate formations of four doped metals closest to the volcanic peak and with stable intermediates. Figure 4a exhibits the free energy diagram of OER on ZrO2. The first step, the formation of HO*, is thermodynamically feasible. Then, HO* dissociates into H+ and O*, which is thermodynamically challenging with a higher onset potential of 2.15 V. The next proton-electron transfer steps of HOO* formation and the generation of O2 in the last step are energetically supported. Therefore, HO* deprotonation to generate O* on ZrO2 is the most challenging step, which is the potential determining steps (PDS), requiring a η of 0.92 V to make all steps thermodynamically downhill (watermelon red line). Similarly, HO*→ O* is the PDS of Ti-ZrO2 with a η of 0.63 V, Pt-ZrO2 with a η of 0.54 V, Fe-ZrO2 with a η of 0.36 V and Rh-ZrO2 with a η of 0.26 V (Fig. 4b-e), which aligns closely with the microkinetic volcano diagram trend (Fig. 3a). The specific adsorption energies of O* and HO* are calculated in Fig. 4f. The order of adsorption strength for O* is Rh-ZrO2 > Pt-ZrO2 > ZrO2 > Fe-ZrO2 > Ti-ZrO2, and that of HO* is Pt-ZrO2 > Rh-ZrO2 > ZrO2 > Fe-ZrO2 > Ti-ZrO2, which is further explained by the charge density difference map (Fig. S10 and inset in Fig. 4). In short, the free energy diagram reveals that Fe-ZrO2 and Rh-ZrO2 are highly active catalysts with only slightly higher ΔG.
Fig. 4 Gibbs Free energy diagram of OER on a ZrO2, b Ti-ZrO2, c Pt-ZrO2, d Fe-ZrO2, and e Rh-ZrO2. f Adsorption energy values of O* (pink line) and HO* (blue line) on five varied substrates. The inset is the charge density difference of HO*, and the isosurface value was set to 0.001 e Å−3 |
For the comparison of different adsorbate species, we relied on thermodynamic analysis, while the dynamic factor was lacking. Therefore, we performed AIMD simulations on Rh-ZrO2 with different species based on its optimal catalytic activity (Fig. S11a-c). The results show that they are thermodynamically stable. In addition, we also evaluated the thermal stability of ZrO2, Fe-ZrO2, and Rh-ZrO2 (Fig. S11d-f). The results demonstrate that these systems maintain their structural integrity, in which energy and temperature oscillate around equilibrium within thermal perturbations. Thus, they remain thermodynamically stable and catalytic activity for extended duration.
3.4 Interatomic Bond Strength Analysis
COHP can provide insights into the strength and nature of chemical bonds between atoms. Integrated COHP (ICOHP) is a quantitative method for measuring chemical bond strength, and a more negative ICOHP value indicates a stronger bond strength between atoms. As Fig. 5 exhibited, we calculate the COHP and ICOHP of metal-O (O* and HO*) bond to deeply understand the bonding mechanism. The sequence of bonding strength with O* is Rh (− 3.28 eV) > Zr (− 2.07 eV) > Fe (− 1.91 eV) (Fig. 5a–c), and with O in HO* is Rh (− 1.83 eV) > Zr (− 1.66 eV) > Fe (− 1.58 eV) (Fig. 5d–f). The ICOHP values explain the adsorption energy trend well. The strong bonding interaction between Rh and O* and HO* results in a lower energy requirement for stabilizing the OER intermediates, which is the key reason of its optimal catalytic performance. Bader charge and D band center ($\varepsilon_{d}$) are also further quantified (Fig. S12). Compared with Fe dopant, Rh loses more charge to O* and HO* (Fig. S12a) and has a smaller negative $\varepsilon_{d}$ value (Fig. S12b, c), indicating underlying reasons for the strong adsorption capacity of Rh.
Fig. 5 Crystal orbital Hamilton population (COHP) analysis of the interactions between a Rh-O in O*, b Zr-O in O*, c Fe-O in O*, d Rh-O in HO*, f Zr-O in HO*, and e Fe-O in HO* |
Based on comprehensive analysis, both Rh-ZrO2 and Fe-ZrO2 are potential candidates for enhancing OER efficiency. The Rh doping exhibits a lower overpotential, but its shortage limits industrial application. In contrast, Fe, a cost-effective and readily available metal, has not only been successfully doped into ZrO2 in experiments [56,57,58] but also regraded as an ideal dopant for enhancing OER activity [59,60]. Therefore, Fe holds great promise as a modulating agent for industrial-scale OER application.
4 Conclusion
In summary, we have analyzed the effects of 40 different single-metal doping to regulate the OER activity of ZrO2 using a stratified screening process based on spin-polarized DFT-D3, surface Pourbaix analysis, microkinetic modeling, cost-effectiveness analysis, and bonding analysis. We identified 16 metals exhibit improved catalytic activity, with Rh and Fe dopants showing the remarkable improvement. These doped metals reduce the band gap of ZrO2, thereby significantly increasing conductivity. In addition, the thermodynamic free energy diagram shows that the theoretical overpotential of Ti, Pt, Fe, and Rh metals are 0.63, 0.54, 0.36, and 0.26 V, respectively, which aligns well with the kinetic volcano model. Compared with Fe dopant, the smaller negative $\varepsilon_{d}$ value facilitates Rh dopant to stabilize the OER intermediates well and thus exhibits higher catalytic activity, but its high cost and shortage limit industrial-scale application. Therefore, the cost-effectiveness Fe-ZrO2 catalyst holds more promise and practicality. This finding provides valuable insights into the design and development of stable, low-cost, and high-performance OER catalysts for commercial-scale water splitting.
Declarations
Conflict of interest The authors declare that they have no conflict of interest and no competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/s40820-024-01403-7.
