Electrostatic potential-derived charge: a universal OER performance descriptor for MOFs

Metal–organic frameworks (MOFs) provide opportunities for the design of high-efficiency catalysts attributed to their high compositional and structural tunability. Meanwhile, the huge number of MOFs poses a great challenge to experimental-intensive development of high-performance functional applications. By taking the computationally feasible and structurally representative trigonal prismatic secondary building units (SBUs) of MOFs as the entry point, we introduce a descriptor-based approach for designing high-performance MOFs for the oxygen evolution reaction (OER). The electrostatic potential-derived charge (ESPC) is identified as a robust and universal OER performance descriptor of MOFs, showing a distinct linear relationship with the onset potentials of OER elemental steps. Importantly, we establish an ESPC-based physical pattern of active site–intermediate binding strength, which interprets the rationality of ESPC as an OER performance descriptor. We further reveal that the SBUs with Ni/Cu as active site atoms while Mn/Fe/Co/Ni as spectator atoms have excellent OER activity through the variation pattern of ESPC along with metal composition. The universal correlation between ESPC and OER activity provides a rational rule for designing high-performance MOF-based OER electrocatalysts and can be easily extended to design functional MOFs for a rich variety of catalytic applications.


Introduction
Metal-organic frameworks (MOFs) are a class of materials assembled from SBUs and organic ligands according to a certain topological structure. 1,2 Owing to their nature of large specic surface area, penetrating pore structures, abundant and atomically dispersed coordination unsaturated active sites (CUSs), exible tunability of the composition and structure, etc., these materials have been extensively studied in molecular storage/sieving, 3-7 catalysis, 8-10 sensors, [11][12][13] batteries, 14,15 and phase transitions. [16][17][18] During the past six years, over 100 000 MOFs have been collected in the Cambridge Structural Database (CSD) 19,20 and the number keeps on increasing dramatically. This, however, makes it challenging for conventional trial-and-error experiments to identify promising MOFs for high-performance functional applications.
The development of high-throughput computational methods has provided a solution to screen high-performance MOFs efficiently and effectively. [21][22][23][24][25][26][27][28][29][30][31] Many studies have focused on the screening of gas storage/sieving, [28][29][30] using computationally cheap characterization for the relevant geometry (pore size distributions, pore diameters, surface area, pore volume, etc.) and molecular-level Grand Canonical Monte Carlo (GCMC)/molecular dynamics (MD) simulations. These efforts have greatly contributed to the development of high-performance MOFs by screening promising materials from tens of thousands or even half a million MOF structures. However, high-throughput computational screening is also called "bruteforce screening," 32 meaning that it can only handle low-cost molecular-level simulations. Even with high-performance supercomputers, it is still difficult to deal with expensive electronic-level calculations. There are few electronic-level screening efforts involving specic catalytic reaction such as photocatalysis 27 or thermocatalysis. [22][23][24][25][26] Due to the costly quantum mechanics-based computational approach, the electronic structure or reaction barrier can only be calculated for a few tens to hundreds of MOFs even using an optimized workow.
Herein, in order to break the barrier of traditional "bruteforce screening" and achieve rational design of MOFs with little computational cost, we develop a promising descriptor-based pathway. Specically, the pivotal reaction driving the transformation of energy structure, like the OER, was chosen as the target reaction, and the secondary building units (SBUs) that play a decisive role in catalytic performance were selected as a rational simplied scheme for MOFs. In particular, we selected the widespread trigonal prismatic SBUs possessing a small number of atoms as a case, whose characteristics such as universal atomic 5-coordination conguration and 3-metal coupling endow them with excellent transferability and compositional tunability, making them an ideal platform for mining descriptive models. With monometallic SBUs, we identify ESPC as a highly accurate, instructive, transferable, and readily available OER performance descriptor. The descriptor was then employed to predict the OER activity of SBUs with bimetallic composition, and the accuracy of the descriptive model was improved using DFT-derived data. Furthermore, by correlating ESPC with OER performance of SBUs with different structures and compositions, a universal linear scaling relationship between ESPC and OER performance was identied, illustrating the universality of ESPC as an OER performance descriptor. For the purpose of using ESPC to guide the actual experimental synthesis, we deeply analyzed the relationship between ESPC and bond order as well as band structure, established an ESPC-based physical pattern of active siteintermediate binding strength. Moreover, the variation regularity of ESPC with metal composition was revealed, which demonstrated that the MOFs with Ni/Cu as active sites and Mn/ Fe/Co/Ni as spectator atoms have excellent OER activity.

Methods
Spin-polarized density functional theory (DFT) calculations were performed by using the linear combination of atomic orbitals (LCAO) method embedded in DMol 3 . A polarized basis set of double numerical polarization (DNP) 33 was employed to ensure the computational accuracy as well as the descriptions of the hydrogen bond. The electronic exchange-correlation effect was described by the Perdew-Burke-Ernzerhof (PBE) 34 functional of general gradient approximation (GGA) level. The large orbital cutoff of 6Å was chosen for good convergence, and TS 35 method was adopted for a long-range van der Waals (vdW) correction. For description of the electron-ion interaction, the DFT-based relativistic semi-core pseudopotentials (DSPPs) 36 were adopted. For reliable adsorption energy, the basis set superposition error (BSSE) correction 37 was adopted, by specifying the catalysts as Counterpoise 1 and adsorbed O-contained species as Counterpoise 2. During the geometric optimizations, all atoms were allowed to relax until the models reach the ground state without any imaginary frequency. The convergence criteria for self-consistent eld (SCF) iteration and the max force on each atom were 10 −8 Ha and 0.001 HaÅ −1 , respectively.
The method for bond order (BO) calculation used in this work was given by Mayer in 1986: 38 where BO AB is the bond order between atom A and atom B, P a and P b are the density matrices for spin a and b, and S is the overlap matrix. The BOs obtained by this method are close to the corresponding classical values, and unlike Mulliken BOs, Mayer quantities are less dependent on the basis set choice and they are transferable, so they can be used to describe similar systems.
The widely studied four-electron mechanism is used to simulate the entire OER process: 39 where the * represents the bare catalyst while the HO*, O*, and HOO* denote different intermediates adsorbed on the catalyst, respectively. For each elemental step, DG i (i = 1, 2, 3, 4) is the change of Gibbs energy and can be expressed as follows: 41 where DE is the reaction energy calculated by DFT, DZPE represents the change of zero-point energy under the harmonic approximation, and TDS represents the entropy correction. The computational hydrogen electrode (CHE) model (G(H + + e − ) = 1/ 2G(H 2 )) proposed by  was adopted to calculate the Gibbs energy of proton-electron pairs (H + + e − ). We used the G(O 2 ) = 4.92 + 2G(H 2 O) − 2G(H 2 ) to derive the free energy of triplet oxygen molecule since the energy calculated by DFT is inaccurate. 45 Based on the above Gibbs energy results, the Gibbs energy difference of the potential-limiting step is dened as follows:

Construction of an ESPC-based OER model
Due to the limitations imposed by the large number of atoms in MOFs for periodic DFT calculations, SBUs, which play a decisive role in the catalytic process, are usually selected as an alternative to MOFs for calculation. 46,47 Among them, trigonal prismatic SBUs are a common class of SBUs ( Fig. 1), whose few atoms make themselves feasible to be implemented in highthroughput calculations. In addition, the strong interaction of three metal atoms connected by a m 3 oxygen atom facilitates the ne-tuning of electronic structure by component orthogonalization. Moreover, the 6-coordinated metal atoms can readily desorb solvent molecules upon activation to form 5coordinated unsaturated active sites that are widely present in MOFs, making the calculation-derived conclusions of such SBUs highly transferable. Therefore, we chose the trigonal prismatic SBUs as an entry point to study the general principles for the design of MOFs with high OER performance. In this study, the SBUs with monometallic components are labeled as M a (M a = V, Cr, Mn, Fe, Co, and Ni), and the SBUs with bimetallic components are labeled as M a 2M b 1 (M a = V, Cr, Mn, Fe, Co and Ni; M b = V, Cr, Mn, Fe, Co, Ni, and Cu). It is important to emphasize that we do not conduct specic screening, but rather pay attention to the descriptors of the OER onset potential by a proper amount of calculations. An adequate linear tting requires data points with large variability, so we use M b as the active site here since in M a 2M b 1 systems, the ESP-derived charge (ESPC) of M b is expected to be more variable than that of M a relative to the monometallic counterpart (due to the relatively low contents of M b , more like the dopant).
To obtain an accurate OER descriptive model, we comprehensively considered three most promising descriptors. They are electrostatic potential-derived charge (ESPC), density of states (DOS), and bond order (BO), which represents the spatial distribution of electrons, the energy distribution of electrons, and the chemical bond strength, respectively. Among them, we calculated the d-band center from the DOS for the convenience  of quantication. For simplicity, we consider only one descriptor to correlate OER onset potentials. Additionally, we focus on the accuracy of the descriptors, while taking into account transferability, instructiveness, and computational simplicity. By correlating the Gibbs energy difference of the four primitive steps (target quantities) with the above candidates, we nd that both ESPC and d-band center exhibit a linear scaling relationship with the target quantities with different accuracy ( Fig. 2a and b). However, the BO shows an exponential scaling relationship with the target quantities, causing its predictive capacity drop sharply at larger values ( Fig. 2c and S1 †). Furthermore, the BO describing the binding strength of the adsorbed molecule to the active site cannot indicate the intrinsic properties of the material, leading to its poor transferability and instructiveness. Therefore, the BO was rst excluded. For the d-band center, it requires expensive DOS calculation and may not be generalized well for MOFs possessing both metallic and semiconducting nature due to the common underestimation of the bandgap, 48 and thus, it is not satisfactory in terms of transferability and computational simplicity. In contrast, ESPC can be obtained by tting the charge distribution near the atom, which is computationally feasible. In addition, it is a physical quantity describing the properties of the atom with desirable transferability and instructiveness. Therefore, ESPC was nally chosen as the OER performance descriptor, and the ESPC-based OER descriptive model derived from six monometallic trigonal prismatic SBUs (Fig. 2a) is as follows: 1, 2, 3, 4) is the Gibbs energy difference of the elemental step in unit eV, and ESPC is the value of electrostatic potential-derived charge in unit e.

Model testing and optimization
The OER performance of the SBUs with bimetallic composition was predicted by the proposed model and the results were validated by DFT calculations. The minor difference between the predicted values and the validated ones indicates the great accuracy (R 2 z 0.87) of the descriptive model (Fig. 3a).
Furthermore, the improved model shows similar predictive ability (R 2 z 0.90) aer being fed with the DFT data of SBUs with bimetallic composition (Fig. 3b), proving the effectiveness of tting scaling relations using monometallic SBUs. In addition, the ESPC provides an excellent description (average R 2 > 0.90) of the OER performance in the case of changing the active atoms while keeping the spectator atoms consistent (Fig. 4), which is of great signicance for studying the doping and active site implanting systems. Moreover, the ESPC could also rationalise the experimentally tested OER activity well ( Fig. S3 and S4 †). The correlations between them are similar to the predictive model raised here, indicating the potential for practical applications.
Based on the proposed model, the strong linear scaling relationships between DG i (i = 1, 2, 3, 4) and the ESPC are shown in Fig. 5 and Table S1, † where DG 1 and DG 2 are negatively correlated with the ESPC, while DG 3 and DG 4 show the opposite trends. This means that there is an onset potential minimum near the intersection of the four prediction lines. When the ESPC is less than 1.1, the height of the nearly coincident DG 1 and DG 2 prediction lines is higher than that of DG 3 and DG 4 , indicating that the potential limiting step (PLS) is step 1 or 2 (* + H 2 O / HO* + H + + e − or HO* + H 2 O / O* + H + + e − ), and when the ESPC is larger than 1.1, the height of the prediction line of DG 3 is the highest, implying that the potential limiting step is step 3 (O* + H 2 O / HOO* + H + + e − ). ,and represent DG 1 , DG 2 , DG 3, and DG 4 , respectively. The shaded red region represents a 95% prediction interval for the linear model.
Transferability is another important criterion for evaluating a descriptive model, especially for MOFs with diverse compositions and structures, as designing material oen involves adjusting composition and structure. Therefore, the ESPC was further employed to predict the OER performance of the SBUs with typical 3-, 4-, and 5-coordinated congurations and representative metals (Fig. 6a-c). Excitingly, the ESPC presents a near-perfect prediction (R 2 = 0.97) of the OER performance of the catalysts (Fig. 6d and Table S2 †), demonstrating the general applicability of the ESPC-based descriptive model. It is pointed out that the prediction accuracy of ESPC for OER activity of typical SBUs appears to be higher than that of trigonal prismatic SBUs. This is because we introduced component orthogonality in trigonal prismatic SBUs to accurately capture the effect of bimetallic coupling on the OER activity. And as can be seen in Fig. 4, in the case of xed M a metals, the prediction accuracy of the model for OER activity of the trigonal prismatic SBUs is similar to that of the typical SBUs.
In addition, our goal is to nd an optimal solution over a vast material space. To validate the proposed model, a 3D volcano plot of OER activity is constructed (Fig. 7, for more details see Fig. S2 †), which shows the evolution trend of the OER performance under the ESPC guidance (blue arrows). It can be clearly seen that the predicted line just passes the top of the mountain, indicating that tuning the ESPC is a desirable option for designing MOFs with high OER performance.

Correlation analysis of ESPC and metal combination
Uncovering the relationship between ESPC of the active site and metal composition is essential to fully utilize the descriptive model to guide the synthesis of high-performance materials. In view of this, we mapped the ESPC of active site onto the metal composition matrix (Fig. 8 and Table S3 †). The result shows that the ESPC varies in a wide range of 1.0-1.9. When keeping the M a metal unchanged, the ESPC varies with the M b metals as following sequence V > Cr > Mn > Fe > Co > Ni > Cu. In addition, when M b is xed, the ESPC changes slightly with the difference of M a metals, and the overall trend shis towards larger values  with the increase of the effective nuclear charge of M a . The matrix shows that the ESPC is in the best interval when the spectator atoms are Mn/Fe/Co/Ni, and the active site atom is Cu or Ni. Among them, Ni atoms have been widely reported as OER active sites in various catalytic materials, 49-51 such as carbonbased materials, transition metal suldes and MOFs. This phenomenon is consistent with the nding in this study. While very few studies report Cu as the active site for OER, thus,  results in this study may shed light on the design of future Cubased OER catalysts. More importantly, by analyzing the related experimental studies, 52-58 we obtained the OER performance ranking for MOFs containing different elements (i.e., Ni > Co > Fe, Fe > V, Fe > Mn), which is consistent with the ndings of our results, demonstrating that the model proposed our work can capture the variation patterns of OER performance along with metal compositions.
The variation trend of ESPC with metal composition can be explained by the difference of the element's electronegativity. For atoms in the same row of the periodic table, the distance between the valence electrons and the nucleus is similar. A bigger effective nuclear charge leads to stronger attraction to valence electrons, and the atom is less likely to lose electrons. 59,60 In the single-metal systems, metals with larger electronegativity have higher electron retention and smaller ESPC ( Fig. 9a and c). In bimetallic composition systems, the valence electrons of three metal atoms are redistributed through the m 3 -O, and the electrons will ow to the atoms with larger electronegativity. For example, for atoms with less effective nuclear charge (like V), when they are combined with atoms with more effective nuclear charge (like Fe), their valence electrons tend to migrate to the counterpart, making the ESPC of Fe atoms in V2Fe1 less than that of Fe in the SBU with single Fe composition. In contrast, the ESPC of V atoms in V2Fe1 is larger than that of the V in the SBU with single V composition (Fig. 9b). It should be noted that this phenomenon is not distinct for systems in which Co, Ni, and Cu as M b , probably due to their "harder" electron cloud stemming from relatively larger electronegativity.

Interpretability of ESPC as an OER activity descriptor
Since the DG total of the OER is a constant (4.92 eV), the DG i (i = 1, 2, 3, 4) actually has a trade-off relationship. Taking V, Co2Cu1 and Mn2Cu1 as examples (Fig. 10), high bonding strength of V atoms with OH and O species induces excessive bonding energy, even making DG 1 and DG 2 become negative values. This strong adsorption leads to little activity of O, so the thermodynamic barrier between the intermediate states of V-O and V-OOH is very high, leading to an ultra-high overpotential (4.27 eV). In contrast, for Co2Cu1, the lower bonding strength of Cu atoms with OH and O facilitates the bonding of OH on Cu-O, which makes DG 1 and DG 2 larger and is not desirable for the entire OER process. However, for Mn2Cu1, the adsorption strength of Cu with three oxygen species is moderate, and the balanced thermodynamic barrier is benecial for obtaining excellent OER performance.
The relationship between ESPC and the adsorption strength of oxygen-containing intermediates can be revealed by BO. As shown in Fig. 11    decreases rapidly, which indicates that the adsorption strength of the O intermediate is weak, and the PLS prefers to occur at the rst or second elemental step. When ESPC is more than 1.1, the adsorption of the O intermediate is strong, and thus, the PLS occurs at the 3rd or 4th elemental step.
The relationship between ESPC and the adsorption strength of oxygen-containing intermediates can be understood more deeply from the perspective of band structure. As the decrease of ESPC, the d-band center becomes progressively lower relative to the Fermi level (Fig. 12), reducing the energy of the antibonding state formed between the adsorbed species and the catalyst. As a result, electrons ow more easily into the antibonding state, leading to an increase in the total energy and weakening the binding strength of intermediates on the catalyst. This weakening of the binding strength is particularly pronounced for O intermediates since the double bond nature of M]O is accompanied by more electron transfer during the adsorption. In addition, with the decrease of ESPC (the active site is from V to Cu, in the order of increasing atomic number), the overall energy of the orbitals that overlapping between the metal 3d orbital and the oxygen 2p orbital (green area in Fig. 13a-g) decreases, which improves the electron lling degree of the overlapping orbitals (Fig. 13h). The stability of the metals is gradually enhanced, indicating that the "bonding ability" of the metals is weakened and leads to the reduced adsorption strength of the OER intermediate on the metals. Thus, the ESPC can be used as a simple indicator to quantify the information of the band structure, which is directly related to the OER performance.

Conclusions
In summary, we unravel that ESPC, a universal descriptor, can describe the OER activity of MOFs containing diverse SBUs very well. By analyzing the variation of ESPC with combinations of different metals, a general design principle was further put forward to enhance the OER activity of MOFs, that is, the active sites should be Ni/Cu while the spectator atoms are Mn/Fe/Co/ Ni. This principle is physically meaningful since ESPC indeed is Fig. 12 The projected density of states (PDOS) for the d-band of the active site in trigonal prismatic SBUs. The vertical dashed lines correspond to the Fermi level, which is shifted to zero and the solid lines denote the d-band centers. The occupied orbitals are filled with color. Fig. 13 (a-g) The projected density of states (PDOS) for the 3d-orbital of the active site and the 2p-orbital of related coordinated oxygen atoms in trigonal prismatic SBUs, the vertical dashed lines correspond to the Fermi level, which is shifted to zero and the overlapping orbitals are filled with green. (h) Correlation between −DG Max and −unoccupied electrons/occupied electrons, the dashed lines provide guidance for eyes. a reection of the bonding strength between metal and oxygencontaining intermediates. The appropriate metal conguration can optimize the d-electron distribution of the active site, which reduces the height of the anti-bonding state of the metaloxygen species and weakens the bonding strength between them as well. All these effects synergistically contribute to the high-efficiency OER process. Our work demonstrates that the regulation of metal components of SBUs is a promising strategy for developing MOF-based energy conversion devices.

Data availability
The structures and calculated data supporting the ndings of this study are available in the ESI. †

Conflicts of interest
There are no conicts to declare.