An effective scheme to determine surface energy and its relation with adsorption energy

Surface energy is fundamental in controlling surface properties and surface-driven processes like heterogeneous catalysis, as adsorption energy is. It is thus crucial to establish an effective scheme to determine surface energy and its relation with adsorption energy. Herein, we propose a model to quantify the effects of the intrinsic characteristics of materials on the material-dependent property and anisotropy of surface energy, based on the period number and group number of bulk atoms, and the valence-electron number, electronegativity and coordination of surface atoms. Our scheme holds for elemental crystals in both solid and liquid phases, body-centered-tetragonal intermetallics, fluorite-structure intermetallics, face-centered-cubic intermetallics, Mg-based surface alloys and semiconductor compounds, which further identifies a quantitative relation between surface energy and adsorption energy and rationalizes the material-dependent error of first-principle methods in calculating the two quantities. This model is predictive with easily accessible parameters and thus allows the rapid screening of materials for targeted properties.


Introduction
Surface energy is fundamental and dominant in controlling surface structure, reconstruction, roughening, nanoparticles' size, and crystal's shape for solids [1][2][3][4] . As surfaces are the region where materials interact with media, surface energy is also of great importance in determining surface-driven processes such as heterogeneous catalysis, gas sensing and biomedical applications [5][6][7][8] , as adsorption energy is. Thus, it is crucial to establish an effective scheme to determine surface energy and its relation with adsorption energy, particularly by means of the easily accessible intrinsic characteristics of materials. However, the electronic and geometric factors that dictate the change of surface energy from one material to the next (materialdependent property) and control the anisotropy of surface energy, still remain elusive so far. On the other band, although it is generally assumed a positive correlation between surface energy and adsorption energy 9,10 , the quantitative relationship as well as the underlying physical picture still remains debated.
Many (semi-)empirical models have been proposed for determining surface energy of solids. Broken-bond models [11][12][13] correlate surface energy with the energy of broken chemical bonds but are only applicable into the lowindex surfaces of elemental crystals. The Miedma model 14 and Stefan model 15 predict surface energy with the experimental energy of vaporization and are inconvenient for the application in many cases like surface alloys. The Friedel model 16 describes surface energy with the d-band width and thus is only applicable into transition metals (TMs) and TM alloys. Although these models have been applied with some success, they are not related to the easily accessible intrinsic properties of materials, perform with limited universality and effectiveness, and can hardly reveal the connection between surface energy and adsorption energy.
First-principle methods have been a workhorse in characterizing surface properties. Compared to experiments, the widely used semi-local functionals underestimate surface energy and overestimate adsorption energy on Pt and Rh, but underestimate both two quantities on Ag and Au, exhibiting an unexpected material-dependent property 9,10,[17][18][19][20] . This drawback poses an obstacle to the reliable predictions of surface reconstructions. Random phase approximation (RPA) 21 calculations can resolve the dilemma on Pt and Rh and produce improved results, which was attributed to the better description of electronic structure by RPA than by the semi-local functionals 10 . However, RPA calculations are too expensive to be widely used for material screening and still generate materialdependent error in describing surface energy and adsorption energy 10,19 . These methods clearly indicate a strong correlation between surface energy and adsorption energy, however, their numerical characteristics prohibit the understanding of the underlying physical framework. Hence, it is urgently needed to set out from the intrinsic property of materials to understand the correlation between surface energy and adsorption energy.
Here we propose a universal picture to determine surface energy and its correlation with adsorption energy, by using the period number and group number of bulk atoms, and the valence-electron number, electronegativity and coordination number of surface atoms. This model is predictive for a variety of materials covering elemental crystals in both solid and liquid phases, alloys, and semiconductor compounds. Furthermore, our scheme builds a quantitative relation between surface energy and adsorption energy, which allows the estimation of adsorption energy with surface energy and rationalizes the material-dependent error of first-principle methods in calculating surface energy and adsorption energy.

Results
We study surface energies for 45 elemental crystals, 31 alloys, and 12 compounds including TMs, main-group crystals, Mg-based surface alloys and III-V semiconductors, and cleavage energies for 11 body-centered-tetragonal (AB) intermetallics, 11 fluorite-structure (A2B) intermetallics and 13 face-centered-cubic (A3B) intermetallics. Note that surface energy and cleavage energy denote the energy required to cleave the bulk materials for symmetric and asymmetric terminations respectively. These solids exhibit  The material-dependent property of surface energy. We first attempt to describe surface energy and cleavage energy with surface properties, since it is generally accepted that they are directly correlated with the binding energy of surface atoms [11][12][13] . The d-band model and Muffin-Tin-Orbital theory show that the spatial extent of the metal d-orbitals on surfaces is associated with the number of outer electrons 22 , indicating that the binding energy of surface atoms likely depends on the number of valence electrons. In addition, Pauling electronegativity (χ) is directly related to the interatomic binding energy. We thus try to describe surface energy and cleavage energy using valence number and electronegativity of surface atoms (which are known to be descriptive for adsorbatesurface binding on TMs 23 ), with ψ = S v 2 χ β , where Sv signifies the valence number with the maximum value 12 (including both the sp-and d-electrons for TMs). β is a parameter determined by the contribution of d-orbitals and/or sp-orbitals to valence description and electronegativity. The sp-orbitals contribution is 100% with β=1 for main-group crystals (without valence d-electrons), while d-orbitals and s-orbitals have the equal contribution for TMs (each has β = 1/2). β is 1/2 for Ag and Au and is 1 for the other crystals, because Au and Ag exhibit the full-filled d-band and the low position of d-band center relative to the Fermi levels (-3.56 for Au and -4.40 eV for Ag) 24,25 and bind with other atoms mainly via sp states (although they have valence delectrons). All of the ψ values correspond to elemental crystals are listed in Supplementary Tables 1 and 2. Fig.  1a-c and Supplementary Fig. 1a-c plot the DFT-calculated and experimental surface energy γ versus ψ for elemental crystals 14,[26][27][28] . Clearly, γ is approximately linearly scaled with ψ in a broken-line behavior. However, ψ fails to elucidate the trends of surface energy for alkaline metals and alkaline-earth metals (see the insets in Fig. 1a-c and Supplementary Fig. 1a-c).
To address this issue, we introduce another two new bulk parameters to describe surface energy and cleavage energy, the period number (Np) and group number (Ng) of bulk atoms, since bulk properties are also important for surface stability. The underlying mechanism will be explained in the section of "Understanding and progress of the model". We now propose a new descriptor for describing surface energy and cleavage energy as follows, Np ̅̅̅ and Ng ̅̅̅ are the average period number and average group number for all elements, which are constant 4 and 9. Җ is easy to acquire since the involved parameters are available from the periodic table of elements. All of the Җ values for elemental crystals are listed in Supplementary Tables 1 and  2.
We now plot the DFT-calculated and experimental surface energy γ versus Җ for elemental crystals in Fig. 1di and Supplementary Figs. 1-3 14,26-28 . Remarkably, Җ describes the trends of surface energy very well for alkaline metals and alkaline-earth metals as well as other elemental crystals. Meanwhile, the description accuracy with Җ is overall improved compared to that with ψ. These results imply that the period number and group number of bulk atoms, and the valence-electron number and electronegativity of surface atoms together determine surface energy.
The fitted scaling relations of γ versus Җ are as follows, where k and b are the slope and offset of the scaling relation. We identify that the slope k of the linear relation is approximately -0.035 for crystals with Җ > 17 and 0.20 for crystals with Җ < 17. The turning point of the scaling of γ versus Җ is around Cr and Re that exhibit half-occupied dbands, which can be well understood with the d-band occupation. The higher half of d-bands corresponds to antibonding states and the lower half to bonding states, indicating the most stable bonding at half-occupied dbands 11,26 . Notably, the offset b is about 0 for Җ < 17 and 3.8 for Җ > 17 with small differences depending on the surface orientations, which exactly correspond to the anisotropy of surfaces. The physical origin of the slope k likely stems from the ratio term CN/ CN ̅̅̅̅ (where CN and CN ̅̅̅̅ are the usual coordination number and the generalized coordination number 29,30 of surface atoms) with Supplementary Eq. 1. CN/CN ̅̅̅̅ , which is an indicator of the atomic packing density of surface atoms, provides an effective measure to combine together the surfaces across the different crystals such as fcc, bcc and hcp. For the CPSs of fcc, bcc and hcp (that are  31 and has been demonstrated as a reflection of the novel surface electronic states such as the anomalous surface electron-phonon coupling 32,33 .
Since the experimental surface energy for solids is obtained by extrapolating from liquid-phase measurements [34][35][36][37][38][39][40][41] , we also study surface energies for liquid phases of the considered systems, by considering the coordination difference between solid and liquid phases (see the details where CNl and CNs are the coordination number of bulk liquid and solid phases respectively. The correlation between surface energies of liquids and the electronic descriptor Җ is shown in Fig. 1i and Supplementary Tables 1, 2 and 5. The slopes of the fitted linear relation are 0.150 for Җ < 17 and -0.022 for Җ > 17, which are significantly different from those of solids and are in good agreement with the predictions by our scheme. The anisotropy of surface energy for solids. Here we turn to understand the anisotropy of surface energy, by studying three kinds of crystal structures with 37 different surfaces 28 . We find that the surface energies on each solid exhibit a linear relationship with the generalized coordination number CN ̅̅̅̅ (see where the λ and ξ are constant for a given solid. We further demonstrate that λ scales linearly with the electronic descriptor Җ (Fig. 2d), approximately as, Clearly, the predicted prefactors by Eq. (4) are in good agreement with the fits of DFT-calculated results 28 in Fig.  2d. Our scheme allows one to understand deeply the anisotropy of surface energy for different solids. For solids with large or small Җ such as Cu, Au, Ag, Na, K and Ca, the prefactor λ is small, leading to the weak anisotropy of surface energy for these solids. In contrast, for solids with medium Җ such as Rh, Mo, Co and Tc, the large prefactor term λ makes the anisotropy of these solids become greater. All these predictions by Eq. (4) are in accordance with the DFT findings in Fig. 2 and Supplementary Figs. 4 and 5 28 . By correlating ξ with the electronic descriptor Җ, we obtain a linear scaling between ξ and Җ as Eq. (5), A scheme combining material-dependent property and anisotropy of surface energy for solids. We thus propose an entire expression of surface energy for solids based on the electronic descriptor Җ and the geometric descriptor CN ̅̅̅̅ as: where μ is a constant. This equation quantitatively describes the material-dependent nature and anisotropy of surface energy for solids, and can be used to estimate rapidly the trends of surface energy on different materials and surfaces, as the involved parameters are easily accessible (see more details in Supplementary Note 1).

Generalize the model into alloys and semiconductor compounds.
Alloying can generate much complex bulk structures and/or coordination environments than elemental crystals. Body-centered-tetragonal intermetallics is an important family of binary alloys with formula AB, where elements A and B locate at the vertex and the center of the cube. Fluorite-structure intermetallics is one of the important binary alloys with formula A2B. This kind of alloys has a typical fluorite-like structure where the TM element B is arranged in fcc stacking and the tetrahedral interstice is filled with the main-group element A. Facecentered-cubic intermetallics also plays vital role in alloys with formula A3B, in which elements A and B locate at the face center and vertex of the cubic respectively. Mg-based surface alloys are constructed by substituting one Mg atom at the surface with another element M, which can significantly modulate the properties of Mg surfaces. III-V semiconductor compounds are typical Zinc blende structure with cation (such as Al, Ga and In) forming fcc structure and anion (such as P, As, Sb and Bi) filling in the tetrahedral interstice. Each of these four crystal structures exhibits variable coordination environments from one surface orientation to the next. It is encouraging that our established correlation can be readily generalized into AB intermetallics, A2B intermetallics, A3B intermetallics, Mg-based surface alloys, and semiconductor compounds by considering the environment effect of bulk and surface atoms. For bulk atoms, the period number Np and group number Ng are obtained by averaging the period and group number based on the stoichiometric ratio between A (nA) and B (nB), that is: . For surface atoms, ψ is obtained CN ̅̅̅̅ +5.62 for Җ > 50 of A3B intermetallics, implying that our scheme can also characterize the materials-dependent properties and anisotropy of cleavage energy for these alloys.
For Mg-M surface alloys, the proportion of alloying element (η) is a non-negligible factor that controls how the alloys are formed, thereby affects surface energy, as γ = { (6.0η -0.27) × Җ + b 1 , Җ < 5.9 (6.1η -0.91) × Җ + b 2 , Җ > 5.9 for (0001), (101 ̅ 0) and (112 ̅ 0) surfaces, η is 1/10, 1/9 and 1/8 respectively. Overall, Eq. (7) successfully draws the trends of surface energies for Mg-M alloys on different surfaces (Fig. 3d, Supplementary Fig. 10 and Supplementary Table 12) 43 . Fig. 3e shows that the surface energies of semiconductor compounds also exhibit a linear relationship with Җ 44 with the similar scaling rule as that for elemental crystals γ = -kҖ + b, where the prefactor k for semiconductor compounds is -0.24. These results demonstrate that our scheme indeed captures the intrinsic properties of cleavage energy and surface energy and is thus universal in describing the surface stability of elemental crystals, alloys and semiconductor compounds.
Understanding and progress of the model. To deeply understand the electronic descriptor Җ, we analyze the role of the involved parameters (Np, Ng, Sv and χ) by taking TMs as an example. Fig. 4a shows the d-band width that is usually taken as an index of bond energy, with respect to the period number Np and group number Ng 25 . In each period, the d-band width scales with Ng in a broken-line behavior from group 3 to group 11: it first increases until groups 6~8 and then decreases. In each group, the d-band width increases with increasing the period number for groups 3~6, whereas this proportional relationship is broken for groups 8~11. These results indicate the d-band width is linearly related to Ng in each period and is correlated with Np depending on Ng. According to the tightbinding (TB) approximation, the energy gain of forming surface is thus scaled with (Ng) 1/2 because of the downshift of the occupied states 45 , while the Ng-dependent term likely enters the exponential term of Np. This corresponds to the relation of γ ∝N p √Ng , as we found. In addition, Sv and χ (forming ψ) provides a reasonable description of cohesive energy (see Fig. 4b) 46 . As a result, the electronic descriptor Җ, which combines Np, Ng, Sv and χ, can describe surface energy and cleavage energy effectively. Why the early TMs (groups 3-5) are consistent with the main-group crystals instead of the late TMs in surface energy. Compared with the late TMs, the d-bands of the early TMs exhibit the centers above Fermi energy level and the large width, which are close to the sp-bands of the main-group crystals (see Fig. 4c) rather than the d-bands of the late TMs. Looking from practical effect, the early TMs thus behave similarly with the main-group crystals in surface energy.
We find that CN ̅̅̅̅ is a better descriptor in describing the surface-energy anisotropy of solids compared with CN. While CN ̅̅̅̅ and CN are similar in characterizing the surfaceenergy anisotropy of fcc and hcp crystals, CN ̅̅̅̅ performs better than CN for bcc crystals, with R 2 much larger for CN ̅̅̅̅ (>0.90) than for CN (~0.80) (see Fig. 2 and Supplementary  Figs. 4, 5 and 11). CN is the number of the nearest neighbors for a given surface atom, while CN ̅̅̅̅ considers both the first-and second-nearest neighbors 29,30,47 . Therefore, CN ̅̅̅̅ captures more general geometry of surfaces and is more sensitive to the variation of different orientations, thereby providing a more accurate description for the anisotropy of surfaces.
We now do a brief comparison between our model and the conventional models such as the broken-bond models [11][12][13] , Miedma model 14 , Friedel model 16 and Stefan model 15 . First, the conventional models connect surface energy with another energies or d-band width, which require expensive experiments or calculations when being applied into alloys,  Table 13). For the alkaline metals and alkaline-earth metals that have relatively small surface energies, the MAE of our scheme is only 0.12 J/m 2 . When estimating surface energy on 45 elemental crystals, 66 alloys, and 12 semiconductor compounds (totally 884 different values), the MAE of the predictions by our model relative to the DFT calculations is ~0.23 J/m 2 (see the sheet of "Predicted surface energies" in Source Data file), which is about 0.16 eV/atom, smaller than the approximate error of semi-local functionals. Last but not least, the conventional models can hardly reveal the anisotropy of surface energy, while our model can well characterize this effect. Our model thus provides a more explicit and accurate physical picture for the surface stability of materials (see more details in Supplementary Note 3).
Wulff shape plays a vital role in understanding surface properties of crystals especially those of nanoparticles that have relatively large surface areas 48 Table 14).
The quantitative correlation between surface energy and adsorption energy. We now try to uncover the quantitative correlation between surface energy and adsorption energy. In addition to the expression of surface energy with Eq. (6), we had also demonstrated a general expression for adsorption energy Ead on metallic materials as follows 23 , where X and Xm are the actual bonding number and maximum bonding number of the central atom for a given adsorbate. θ is a constant originating from the coupling between the adsorbate states and the substrate-sp states. Adsorption energies of CO, NO and O are shown in Fig. 4d and Supplementary Fig. 12, where the fitted prefactors are consistent with the predictions by Eq. (8) 49,50 (see Supplementary Note 5). Combining Eqs. (6) and (8), we build a quantitative relation between surface energy and adsorption energy, Surface energy is thus correlated with adsorption energy in both positive and negative relationship, dismissing the conventional concept that Ead is positively correlated with γ. Taking CO adsorption as an example, the underlying mechanism can be understood by the different responses of 5σ-and 2π*-metal hybrid orbitals at the different adsorption sites [51][52][53] (see the details in Supplementary Note 6). It is noteworthy that surface energy is approximately linearly scaled with adsorption energy for the late TMs, as their term is approximately constant.
Our scheme provides a simple picture for understanding the connections and distinctions between surface energy and adsorption energy [see Eq. (9)]. While the valence-electron number, electronegativity and coordination of surface atoms are crucial for both surface energy and adsorption energy, the period number and group number of bulk atoms play an additional decisive role for surface energy. For surface energy, a remarkable characteristic is that the effects of electronic and geometric properties are coupled, leading to a coupling term ҖCN ̅̅̅̅ ( (   56 , and 0.23 eV for strongly constrained and appropriately normed (SCAN) 57 meta-generalized gradient approximation (GGA), most of which are smaller than the ±0.2 eV, the approximate error of DFT (semi-)local functionals. These consistencies demonstrate that our established relation between surface energy and adsorption energy is robust and can be applied into the future materials design.

Discussion
We now study the origin of the material-dependent error of first-principle methods in calculating surface energy and adsorption energy.
For surface energy, the difference between the DFTcalculated and experimental results can be expressed as, For simplicity, we first assume that ∆1 = Җ' -Җ and ∆2 = CN ̅̅̅̅ ' -CN ̅̅̅̅ are both constant for a given method, corresponding to the systematical error of first-principle calculations (Җ' and CN ̅̅̅̅ ' denote the calculated electronic and geometric structures). ∆γ is thus dominated by the descriptors Җ for a given method. Eq. (10) deduces that ∆γ should exhibit a material-dependent nature with an approximately linear scaling with Җ and the corresponding slope, mainly determined by the geometric error ∆2, for Җ < 17 is about five times of that for Җ > 17. These deductions have been explicitly demonstrated by the calculations with LDA, PBE, PBEsol, SCAN, SCAN+rVV10 (Vydrov-Van Voorhis 2010 58 ), and RPA (see Fig. 5a-c and Supplementary Fig. 13) 10,19,28 . Notably, RPA calculations with RPA geometries overestimate but RPA with PBE geometries underestimate surface energy, again identifying the key role of geometric error ∆2 (see Fig.  5b and Supplementary Note 8) 10,19 . As these six methods can describe the electronic structures of metals with reasonable accuracy, the intercept for Җ < 17 is small, indicating that the metals with small or large Җ (Җ < 10 or Җ > 60) likely generate smaller ∆γ than those with medium ψ (10 < Җ < 60). This again has been demonstrated by the first-principle methods that ∆γ is much smaller in Al, Sc, Au and Ag than in Pt, Rh and Re (see Fig. 5a-c and Supplementary Fig. 13) 27,28 . These results strongly support the robustness of our model for surface energy.
In the case of adsorption energy, the error of adsorption energy between first-principle calculations and experiments is as follows, . Indeed, the calculations by LDA, Perdew-Wang-91 functional (PW91) 59 , PBE, PBEsol, Bayesian error estimation functional with van der Waals correlation (BEEF-vdW) 60 , Hammer, Hansen, Nørskov modified PBE functional (RPBE) 61 , and RPA fulfill this prediction, generating linear relations between ∆Ead and ψ for CO, NO, and O2 adsorption on TMs (see Fig. 5d-f and Supplementary Fig. 14) 10,17,18 . Compared to LDA, GGAs and meta-GGAs, RPA exhibits a smaller ∆3 (-0.5 versus -0.1 for CO), suggesting that the material-dependent error is small for RPA in describing adsorption energy. For CO adsorption, the constant θ1,2 term in Eq. (11) is most likely larger for GGAs and meta-GGAs than for RPA, as GGAs and meta-GGAs underestimate significantly the energy of 2π* orbits of CO 9,10 . Due to the synergy of 0.1∆3ψ/(Xm+1) and θ1,2, RPBE, PBE and SCAN overestimate the adsorption energy for metals with small ψ such as Pt and Rh but underestimate the adsorption energy for metals with large ψ such as Ag (see Fig. 5d and Supplementary Fig.  14b) 17,20 , whereas RPA underestimate slightly the adsorption energy on all considered TMs (see Fig. 5e) 10 (see more details in Supplementary Note 7).
Overall, our scheme indicates that the observed material-dependent behavior of the first-principle methods in calculating surface energy and adsorption energy is due to the coupling between the constant systematic errors and the electronic properties of materials denoted by Җ and ψ, in which the electronic structure Җ and ψ strongly depend on materials, as shown in Eqs. (10) and (11).
We have proposed an effective model for the accurate determination of surface energy based on the period number and group number of bulk atoms, and the valence, electronegativity and coordination of surface atoms, which holds for main-and transition-group elemental crystals in both solid and liquid phases, AB intermetallics, A2B intermetallics, A3B intermetallics, Mg-based surface alloys and semiconductor compounds. We find that the electronic properties of bulk and surface atoms control the materialdependent property of surface energy, while the electronic properties of bulk and surface atoms and generalized coordination of surface atoms jointly determine the anisotropy of surface energy. This model correlates surface energy with the intrinsic properties of materials, builds the quantitative correlation between surface energy and adsorption energy, and uncovers the origin for the materialdependent error of first-principle methods in calculating surface energy and adsorption energy. Our findings rationalize some theoretical and experimental results, and could be helpful to engineering surface energy and adsorption energy simultaneously for materials design.

Methods
In the study, the density of states (DOS) calculations were performed by CASTEP code 62 with ultrasoft pseudopotentials 63 and PBE functional augmented with TS method 64 . TM surfaces were modeled with four-layer slabs in a unit cell of p(2×2), where the top two layers were fully relaxed and the rest of the layers were constrained in the optimized lattice. A vacuum of 15 Å was adopted to separate the adjacent slabs. We used plane-wave cutoff energy of 450 eV and the Monkhorst-Pack k-point sampling with 8×8×2 meshes for geometry optimization. The conjugate gradient algorithm was utilized with a convergence threshold of 5.0e -6 eV and 0.01 eV/Å in Hellmann-Feynman force on each atom. It is noteworthy that only the DOS of TMs were calculated by our PBE+TS methods, while the rest data such as surface energies, cleavage energies and adsorption energies were cited from literatures 14,[26][27][28][34][35][36][37][38][39][40][41][42][43][44]49,50 . Surface energy and cleavage energy are calculated with the slab models 26 (12) where Eslab is the energy of the slab, Ebulk is the bulk energy per atom, n is the number of the atoms in the slab, and A is the surface area of the slab. Eq. (12) corresponds to surface energy, when cleavage of bulk yields symmetric slabs 26,28,43,44 , e.g. for elemental crystals, which have two identical surface terminations and equal to surface energy. If cleavage of bulk generates asymmetric slabs, e.g. for intermetallic compounds, which have two distinct terminations and different surface energy, Eq. (12) obtains cleavage energy 42 , which is the average of the surface energy of the two different termination. Note that cleavage energy is identical to surface energy only when the slabs are symmetric. Encouragingly, our scheme works well for both surface energy and cleavage energy.
For Mg-based surface alloys, the surface slabs were constructed by substituting one Mg surface atom with another alloying element M on each termination of the slab to ensure the symmetry. The Mgn-2M2 bulk system contains 94 Mg atoms and 2 M atoms (n is large enough) so that the distance between the 2 M atoms is large enough, hence the interaction between the 2 M atoms can be ignored. Then surface energies of Mg-M surface alloys 43