Jacob’s Ladder as Sketched by Escher: Assessing the Performance of Broadly-Used Density Functionals on Transition Metal Surface Properties

The present work surveys the performance of various widely-used density functional theory (DFT) exchange-correlation (xc) functionals in describing surface observable properties of a total of 27 transition metals with face-centered cubic (fcc), body-centered cubic (bcc), or hexagonal close-packed (hcp) crystallographic structures. A total of 81 low Miller index surfaces were considered employing slab models. Exemplary xc functionals within the three first rungs of Jacob’s ladder were considered including the Vosko-Wilk-Nussair (VWN) xc within the local density approximation (LDA), the Perdew-Burke-Ernzerhof (PBE) within the generalized gradient approximation (GGA), and the Tao-Perdew-Staroverov-Scuseria (TPSS) as a metaGGA. Hybrids were excluded in the survey as known to fail in properly describing metallic systems. In addition, two variants of PBE were considered, the PBE adapted for solids (PBEsol), and the revised PBE (RPBE) aimed at improving adsorption energies. Interlayer atomic distances, surface energies, and surface work functions are chosen as the scrutinized properties. A comparison to available experimental data including single-crystal and polycrystalline values shows that no xc is best at describing all the surface properties, although in statistical mean terms the PBEsol xc functional is advised, yet PBE is recommended when considering both bulk and surface properties. Based on the present results a discussion on adapting GGA functionals to the treatment of metallic surfaces in an alternative way to metaGGA or hybrids is provided. *Corresponding author: francesc.vines@ub.edu

leads to unphysically localized TM bands disrupting the proper description of transition metal properties. 21 Furthermore, a very recent study shows that it is possible to describe the thermochemistry of 3d TMs without having to rely on hybrids approaches. 22 The above mentioned previous studies pointed to xc functionals within GGA approximation as the best choice for describing TM bulk properties. In particular, on average the Perdew-Burke-Ernzerhof (PBE) 23 xc was found to be the most accurate functional while the Tao-Perdew-Staroverov-Scuseria (TPSS) 24 would be the next best suited and broadly used from another Jacob's ladder rung, often considered necessary when dealing with main group element molecular systems, as is the case for adsorbates on surfaces. 25 Note in passing by that there have been recent improvements on semilocal metaGGAs, 26 including intermediate-range description of dispersive forces, which, in principle, can lead to a significantly description upgrading on a diversity of chemical systems. 27 Finally, the Vosko-Wilk-Nussair parameterization of LDA would be the best adapted in that rung. 28  The valence electron density was expanded in a plane-wave basis set with a 415 eV cutoff for the kinetic energy, while the effect of the atomic cores into the valence electron density was described through the projector augmented wave (PAW) method. 37 A tetrahedron smearing method of 0.2 eV width was used to speed up the electronic convergence, 38 yet final energy values were corrected to 0 K. An electronic convergence criterion of 10 -6 eV was used, and ionic relaxation was considered converged when forces acting on atoms were smaller than 0.01 eV/Å. The electronic structure calculations were by default non-spin polarized, except for magnetic Fe, Co, and Ni bulks and surfaces. An optimal Monkhorst-Pack 39 grid of 7×7×7 k-points dimensions was found to be sufficient for accurate bulk total energy calculations in most stringent metals, and therefore used in all cases. 19,20 In the case of slab calculations, a 7×7×1 k-points grid was used.
The surface energies, γ, are calculated following Eq. 1, where E slab is the energy of the optimized slab, E bulk the energy of a metal atom in its bulk environment, and N the number of atoms in the employed slab model. The A term is the exposed surface area exposed in each of the two exposed facets of the slab model.
The work function, φ, is defined as the minimum energy needed to remove an electron from a solid, located in its Fermi level (E F ), and place it in the vacuum energy level, V, then; where, in order to acquire V the electron electrostatic potential energy is averaged for each surface along the normal to surface direction, until a constant value is found in the vacuum region, see Figure 2. E F level was obtained from the total density of states (DOS) sampled by circa 10,000 points, ensuring a numerical precision of 0.001 eV.
Finally, the surface relaxation is evaluated as the layer contraction/expansion percentage, Δ ij , given with respect bulk environment. Note that for the studied TMs the interlayer spacing in between two vicinal layers i and j is constant for a given crystallographic direction in bulk, and so Δ ij is given as where !" !"#$ is the bulk interlayer distance along the direction perpendicular to the surface whereas !" is the equivalent interlayer distance in the slab model. In the bulk, the ij index is a constant whereas in the slab refers to a pair of consecutive atomic layers. Hence, i = 1 corresponds to surface layer, j = 2 to the first subsurface layer. Thus !" would refer to the interlayer distance between the surface and the subsurface layers in the slab model. Within this definition, negative values of Δ ij imply an interlayer distance shortening, whereas positive values refer to interlayer distance lengthening.

Results and Discussion
For convenience, results concerning surface energies are presented first, followed by a discussion on the results for surface work functions. The analysis of the structural relaxation is then reported at the end of this section, followed by an overall discusison.

Surface Energies
In order to determine the best suited xc functional for this property we first compare the computed surface energies with the zero Kelvin extrapolated experimental data. Notice that such extrapolation is typically done from polycrystalline samples near melting temperatures, and likely correspond to an admixture of different exposed surfaces. 30 However, it is unclear whether the zero Kelvin extrapolated data would belong to the most stable surface, or still to an admixture of competing surfaces. For that purpose we gained Wulff constructions such as those described in Figure 3. For each transition metal, the Wulff constructions were obtained from the calculated surface energies obtained for each xc functional. For fcc and bcc structures these were built either using the visualization for electronic and structural analysis (VESTA) package 40 whereas the WinXMorph suite was used for the hcp ones. 41 There, for each TM and each xc, the constructed Wulff shapes provides the percentage of area exposed by each exposed surface. The specific computed surface energies, !"#! for each particular surface and each employed xc are found in Table   S1 of the Supporting Information, alongside with the experimental surface energies, !"# .
The obtained zero Kelvin Wulff contributions are listed in Table S2.
The performance of the different xc functionals in describing surface energies is analyzed in Figure 4 by comparing the computed estimates to experimental values. Notice that two comparisons are here made, either using those !"#! of most stable surfaces, consistently evaluated for each TM surface and xc, or Wulff shape averaged, thus considering contribution fractions of all surfaces of a given TM at a given xc level, !"#$$ !"#! . The latter provides a better comparison since, as mentioned, a large number of experimental values correspond to an average over exposed surfaces. For clarity, only linear fittings are plotted; although dispersion is evaluated by the linear regression factor, R, see Table S3 in Supporting Information, and the error analysis listed in Table 1. Note here that quantitative analysis of accuracy is made based on mean errors (ME), mean absolute errors (MAE), and mean absolute percentage errors (MAPE).
A close inspection on Figure 4 shows that the approximation of comparing calculated to measured surface energies assuming that experiments correspond to the most stable facets -!"#! case-can be justified since experimental trends are duly followed, especially when considering VWN and PBEsol functionals, which display linear regression with coefficients R over 0.90, slopes near unity, and intercepts of only -0.01 (VWN) and -0.12 (PBEsol) J/m 2 .
The performance of PBE and RPBE in this point can be considered as quite good, whereas metaGGA TPSS displays a poor adjustment, with an R of only 0.51, see Table S2.
Considering an admixture of surfaces, here simulated by Wulff averaged surface case-one realizes that the a better fitting is obtained, see Figure 4. This is accompanied by slightly better R coefficients, slopes, and intercepts, see Table S2.
The only exception to this fact is RPBE; The Wulff averaged calculated value leads to a slope increase from 0.85 up to 0.95, and the intercept changes from -0.34 J/m 2 to 0.18 J/m 2 . This apparent better fit is indeed not, as R drops from 0.88 to 0.56. As a matter of fact, the poor description of surface energies at RPBE does not only apply to the most stable surfaces, but also to others, which scales to the admixture of all contemplated surfaces, see Table S1. This unrealistic treatment of surfaces on an equal foot upraises the Wulff overall surface energies in a quite unbalanced fashion, leading to a higher dispersion. Therefore, the apparently better fitting is an artifact stemming from the poor performance of RPBE in describing surface energies, rather than a proper description based on a realistic physical description.
A quantitative analysis based on the statistical errors is shown in

Work functions
Next, we discuss the performance of the different functionals on predicting surface work functions. The complete set of calculated values is reported in Table S5 in the Supporting Information. Following the same procedure as in the previous section, we rely on the Wulff constructions to obtain an averaged value of φ to be compared with the experimental values as obtained from polycrystalline samples. 31 However, for this observable, there is also a significant number of surface specific measurements obtained from TM single crystals cut in the desired studied directions. Therefore, for this subset of data, a direct comparison is possible and the list of experimental values is reported in functions change with temperature by between 10 -4 to 10 -6 eV/K. 43,44 Therefore, the disagreement due to thermal effects at standard conditions cannot exceed 0.04 eV.
Consequently, this small deviation is neglected in the oncoming analysis and discussion.
Two sets of values are used for comparison to polycrystalline values, !"#$$ !"#! , using Wulff-averaged values, or single crystal data well defined surfaces, !"#$%& !"#! . The trends from linear fittings are reported in Figure 5, whereas the corresponding statistical analysis is summarized in Table 2. As one can see, a better correlation is obtained when using the one-  Figure 5, that, at variance to trends for calculated surface energies, the trends on workfunctions appear to show a somewhat overestimation of large φ values, and an underestimation of small φ ones, with transient points located at ~4 and ~3 eV for single-crystal and polycrystalline cases, respectively. For the oncoming analysis, we decided to better rely on the single-crystal data, given that slab models provide a realistic representation of these systems. Thus, we neglect other higher order Miller index surfaces, and other low-coordinated sites on the nanoparticle approximation, which can affect the simulated values. From this analysis, PBE appears as the best performing xc functional.
Considering the linear fittings, some caveats ought to be commented. The slope (a), intercept (b), and R factors of Figure 5 linear fittings are reported on Table S6. These This said, a further close inspection to the crystallographic itemized statistics, shown in Table  S7, indicates a general pattern in the sense that work function of fcc metal surfaces are better described than those of hcp and bcc ones, and for this situation, TPSS performs better.
However, bcc are better described when comparing to polycrystalline data, although there the best description comes from the PBE functional. At variance with the above commented case of surface energies, the work function experimental accuracy, below ±0.3%, barely helps in concealing estimates and measurements.

Interlayer distance
Finally, as far as surface relaxation is concerned, a limited amount of experimental data is available, most focused on the surface and first subsurface interlayer distance, δ 12 . In addition, the experimental reference data for the full set of surfaces and TMs is not complete with a sensible heterogeneity of values with different accuracies and experimental precision, sometimes contradictory for a particular TM surface. Because of this we selected as reference data those which display i) higher precision, ii) a larger amount of interlayer distances, and iii) for i) and ii) being on equal foot, those obtained most recently. The full set of data and concomitant references are reported in Table S8 of the Supplementary Information, as well as those used values highlighted in bold.
The calculated data at all the explored xc levels, including the three different interlayer relaxation values (Δ 12 , Δ 23 , and Δ 34 ) are contained in Table S9 of the Supplementary Information. However, the relaxation percentage as usually provided in the experiments is of little use in evaluating xc performance, as it is already a referenced data.
Because of this, it seems more justified to directly compare the xc computed interlayer distance, δ ij , to the experimental values, derived from the experimentally available interlayer relaxations. To obtain !" !"# , the experimental interlayer percentage relaxations were applied to the 0 K extrapolated interlayer distances. 20 The computed values are reported in Table S10 of the Supplementary Information. The performance of the different methods on predicting δ ij has been analyzed statistically and results are reported in Table 3. Graphically, the calculated values nicely fit the experimental set of values, as seen in Figure 6, with very small deviations of typically around 0.15-0.22 Å (Table 3) with little variations in between the different considered xc functionals. This is also reflected in the data from the linear fitting shown in Table S11 in the Supplementary Information, with slopes near ~0.8, intercepts below 0.4 Å, and R 2 circa 0.8. This is in agreement with the rather good description of geometrical structure by LDA, GGA, and meta-GGA xc functionals in general, and on bulk TMs in particular. 19,20 Despite of this, Figure 6 shows that PBE performs slightly better with comparison to experiment superior to that corresponding to other explored xc functionals. This is aligned with the PBE smallest MAPE value in Table 3. However, this result varies with the crystallographic structure, shown in Table S12, with bcc being the most accurate for this property. In addition, it is important to remark that experiments show a wide range of precision, ±1.5 % on average, which actually closes the gap of any DFT calculation, and actually, the differences among functionals seem to be not so acute as the previously found for γ and φ observables.

Overall Functional Assessment
To summarize, a general view of the performance, in terms of MAPE, of the different explored functionals on predicting surface relaxations, surface energy, and surface work function, is provided. Notice that concerning the assessment of the different functionals, results evidence heterogeneity on the properties under study -even crystallographic groups, and specific TMs, see Figure S1. Even so, an overall better performance for specific cases is observed for VWN concerning γ, and for RPBE on φ, yet for the latter property the mean best performance is again for PBE. On average, VWN is best adapted to surface energies, and PBE to work functions and interlayer distances. However, for a general assessment of surface properties, one can add up the obtained MAPEs for each xc functional under study, as shown in Figure 7. A close inspection reveals that, according to this criterion, the most balanced surface properties are provided by PBEsol, despite not being the functional providing better fitting between calculations and experiment for the properties under inspection. Further than that, combining the present analysis with that for TMs bulk properties -cohesive energies, E coh , bulk moduli, B 0 , and shortest interatomic distances δ-as obtained earlier with the same procedure here applied, 19,20 the most accurate xc functional is PBE, yet closely followed by PBEsol and TPSS.
In this sense one would not advice VWN and RPBE when computing TMs bulk or surfaces properties, although VWN is particularly suited when estimating surface energies.
All that said, RPBE has been claimed to represent an improvement regarding the adsorption of main group molecules on TM surfaces, 36 a point of the scope of the present study, although the metaGGA TPSS functional also represents an improvement in the description of the thermochemistry of main group systems, a point that, coupled to present good performance, can place TPSS in a best-compromise situation when studying the interaction of atoms and/or molecules on TM surfaces. 20,24 Moreover, the present study further reveals that, at least for TM extended systems, the development of new xc functionals needs to consider a broad number of properties much as in the same spirit of the G2 and related datasets used in molecular quantum chemistry. 45 This has been exemplified here by assessing the performance of a series of functional on predicting surface (and bulk) properties, which hopefully serves as a spur and guidance for future DFT xc improvements.
Along    Tags for the different surfaces are located nearby the exposed facets.

Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI:   Table S9 and S10 provide the computed interlayer relaxations and distances. Figure S1 shows, color-coded, the best functional for each bulk and surface TM property.