Performance of Minnesota functionals on predicting core-level binding energies of molecules containing main-group elements

Here we explored the performance of M06, M06-L, M11, and M11-L Minnesota functionals on predicting core-level 1s binding energies (BEs) and BE shifts (ΔBEs) for a set of 20 organic molecules containing main-group elements C → F (39 core levels in total). The broadly used Hartree–Fock (HF) and Becke–Lee–Yang–Parr (B3LYP) methods have also been studied for comparison. A statistical analysis comparing with X-ray photoelectron spectroscopy (XPS) experimental values shows that overall BEs estimations only deviate a small percentage from the experimental values, yet the absolute deviations are generally too large, with the different methods over/underestimating the reported values. However, taking the contribution of relativistic effects of BEs into account leads to larger differences. Overall, the performance of the explored Minnesota functionals is not satisfactory, with errors of up to 1 eV, except for the M06-L meta-GGA functional. In this case, the mean absolute deviation is below 0.1 eV and thus within XPS chemical resolution. Hence, M06-L poses itself as a rather accurate and computational expense-wise method for estimating BEs of organic molecules. Nevertheless, the observed deviations almost cancel when considering ΔBEs with respect to some arbitrary reference, with errors within 0.2–0.3 eV, indicating that these are largely systematic, which in turn implies that the corresponding methods have room for improvement.


Introduction
The electron Spectroscopy for Chemical Analysis (ESCA), more commonly known as X-ray Photoelectron Spectroscopy (XPS), is an experimental technique widely used in many materials and surface science laboratories and facilities, either from research, or in applied industry.It is mainly used for the elemental analysis of bulk materials, especially for surfaces, given its surface sensitivity [1,2], but it has also been applied to the detection of gas phase molecules [3].Furthermore, XPS is currently used to observe in situ the evolution of an heterogeneously catalyzed reaction, allowing for the characterization of reactants, intermediates, and products, thus serving as a powerful tool to determine the reaction mechanism [4][5][6].
The XPS usefulness hangs on the elemental analysis by measurements of core level electron Binding Energies (BEs).The BEs are characteristic of a given element or, more specifically, of a given element in a given chemical environment and electronic state.Thus, BEs provide qualitative information the species present in a sample, quantitative information of their quantity, plus qualitative/quantitative information of the different chemical environments and electronic states of a given detected element [1][2][3][4][5][6].Such a detailed information enables BEs from XPS experiments to be used as chemical fingerprints in condensed phase systems, reflecting the chemical properties and bonding between the species in it [7].The small variations for a given element in different chemical environments are often weighted by means of BEs shifts (ΔBEs), and normally allow for distinguishing structural as well as oxidation state of the atom [8].
However the assignment of a given XPS peak to a given element chemical environment is by no means straightforward, and, actually, ab initio estimations come at hand for such purpose.Indeed, core level BEs and their ΔBEs have been accurately predicted from ab initio calculations at the well-known Hartree-Fock (HF) level of theory [9][10][11], but system size limitations have driven the research endeavors in using Density Functional Theory (DFT) based methods for such purpose [12][13][14].
A milestone study in this aspect is the one of Takahata and Delano [15], who studied 35 small organic molecules containing B→F atoms, by estimating a total of 59 1s core electron BEs, and tackled the effect of using increasingly in size basis sets and a sequence of 21 different DFT exchange-correlation (xc) functionals, either within the Local Density Approximation (LDA) or the Generalized Gradient Approximation (GGA), but also considering meta-GGA and hybrid xc functionals.There the Voorhis-Scuseria (VS98) [16] metaGGA and the Becke-Lee-Yang-Parr (BLYP) [17,18] GGA where found to be the best tested choices for computing BEs, with Mean Absolute Errors (MEA) of ~0.2 and ~0.3 eV, respectively.
However, this previous study did not include some popular xc functionals used for main group element molecules, these are the hybrid functional, which incorporate a percentage of HF exchange in their expression, and known to improve the thermochemistry of main-group molecules.A clear representative of them is the threeparameters BLYP hybrid xc (B3LYP) [19], which has been recently found to be acutely suited in estimating BEs for a series of N-containing molecules in gas phase [20,21].
Another family of meta-GGA and hybrid xc functionals that have not been explored in this matter are the so-called Minnesota functionals, developed by the working group of Prof. Donald Truhlar, a series of continuously improved xc functionals aiming to target an overall description of any chemical system [22].
In the present work we want to bridge this gap by evaluating the performance in estimating 1s core level BEs for the B3LYP functional and a series of Minnesota xc functionals.In particular, we chose the M06-L meta-GGA functional [23], given that such xc has proven to deliver an accurate description even in very complex systems [24], like the known CO adsorption on Pt(111) puzzle [25].The M11-L meta-GGA functional [26] represents an update with improved accuracy with respect M06-L.For any of these two variants, a hybrid xc functional can be constructed, by adding 27% and 42.8% of HF exchange to their ansatz to the M06-L and M11-L xc functionals, respectively, forming the so named M06 [27] and M11 [28] hybrids.For comparison with previous works [15], HF method has been also tested.
A possible way to obtain the energy values in Eq. 1 is to make use of separate Self-Consistent-Field (SCF) calculations.The resulting procedure is usually referred as ΔSCF [9,11] and Eq. 1 usually rewritten as where the subindex i in  !indicates the ionized core whereas E N SCF and E i N-1  are the variationally optimized energy for the initial system with N electrons and the final systems with (N-1) electrons and the corresponding i core hole.
The ΔSCF calculations have been performed on the molecular set as proposed by Takahata and Delano [15] [30,31].
For the core hole state the occupied orbitals are selected using an overlap criterion instead of the usual Aufbau approach.All calculations have been carried out in a spin restricted fashion and are non-relativistic.
Note by passing by that relativistic effects are different for different core levels and increase with the atomic number of the ionized atom.In order to discuss the accuracy of the different methods in predicting BEs it is convenient to have a reliable estimate of the contribution of the relativistic effects.To this end, results from relativistic and non-relativistic calculations for the C→F isolated atoms at the HF level of theory provided by Bagus [32] have been used.These relativistic calculations were carried out with the DIRAC program [33] and the non-relativistic calculations were carried out with the CLIPS code [34].The wave functions were based on the average of configurations and do not take into account the multiplet splittings for these open shell atoms [35].We compared fully relativistic four-component Dirac HF wave functions and energies with non-relativistic HF wavefunctions and energies for the C→F isolated atoms.The basis sets used for these calculations were the same as for the other calculations.It is worth pointing out that previous works [20,21] validated the GAMESS results for the core hole states by comparing to results obtained with CLIPS.

Results and Discussion
Let us first analyze the absolute BEs -BE(ΔSCF) -results, which are encompassed in Table 1.At a first glance, any functional explored, and also HF, are targeting well the experimental BE.Indeed, the statistical analysis, evaluated in terms of

Mean Error (ME), Mean Absolute Error (MAE), and the Mean Absolute Percentage
Error (MAPE) show it so: Any functional here explored, even the HF method, excellently targets the experimental value with a deviation below 0.17%.When one goes into detail, HF tends to underestimate the experimental BEs by ~0.3 eV, probably due to the known localization of HF exchange, which limits the electronic relaxation once the core electrons is suppressed.This could be at the origin of the underestimation of BEs of B3LYP by the same amount.
However, when one deals with Minnesota functionals, the panorama changes.
Indeed, M06 hybrid xc is found to underestimate BEs by ~0.7 eV, although one has to keep in mind that related meta-GGA M06-L functional already underestimates them by ~0.3 eV.Thus, the underestimation seems to be reside in the functional by construction, although these can simply come from relativistic effects, as shall be commented below.
On the other hand, the M11 and M11-L reformatted xc functionals tend to overestimate the BEs by ~0.3 and ~0.7 eV, respectively.However, it is clear that by adding the corresponding % HF to M06-L or M11-L meta-GGA xc functionals, the BEs lower by ~0.45 eV, as expected.By further analyzing the data using MAE one can remove sign cancellation errors; this way, M11 hybrid and M06-L meta-GGA seem to be the best xc functionals to carry out BE estimations, with a mean error for both near ~0.3 eV, closely followed by B3LYP with a mean absolute error of ~0.4 eV.
An elemental analysis for C→F is visually presented in Figure 1 where excellent linear relationships are drawn for calculated BEs versus the experimental ones.Several conclusions can be withdrawn from this analysis.The first one is that C 1s are rather accurately described by any of the inspected methods, with just very slightly deviations from ideality.Thus, the over/underestimations commented in This agreement is also reflected when plotting, itemized for the different studied elements, the estimated ΔBEs with respect the experimental ones, as shown in Figure 2.
The agreement in between methods is remarkable, and barely one can only highlight a sensibly larger slope for HF, most acute for C 1s and O 1s cases, slightly followed by M11, and faint value dispersion for F 1s cases.Actually, one can inspect the data gained from the linear fittings in this matter, as shown in Table 3.It is clear that linear fitting are excellent for ΔBEs, with R regression coefficients larger than 0.992.One can estimate, for each inspected element, and method, the mean offset for the obtained range of ΔBEs.By doing so, it is clear that the larger the slope, i.e. the more deviation with the ideal value of 1.0, the larger the offset, as observed for C 1s and N 1s ΔBEs, with slopes of 1.12 and 1.27, and offsets of 0.67 and 0.87 eV, respectively.This is also observed for M11 but to a lower degree, e.g. for N 1s , with a slope of 1.11 and an offset of 0.49 eV.Aside from these two cases, the overall offsets range 0-0.2 eV, really close to the ESCA chemical accuracy of 0.1 eV.It is mandatory to highlight that the lowest offset deviation is found for M06-L meta-GGA, below 0.1 eV regardless of the studied method.
Such offset analysis can be also applied to BEs.Indeed, Table 3 contains the linear fittings for the BEs plots in Figure 1.Again, the linear adjustments are excellent, with regression coefficients above 0.99.Here, however, the offsets can be larger despite a slope near one, reflecting the over/underestimation of the chosen xc functional or HF.
Indeed, the deviations can be as high as 1 eV, as observed for F 1s cases.Nevertheless, one has to keep in mind that these results refer to non-relativistic calculations.We indeed considered the influence of relativistic effect on the 1s BEs on the studied elements.Relativistic (Rel.) and non-relativistic (Non-Rel.)calculations have been carried out for the C→F atoms, and the relativistic contributions and the comparison with the non-relativistic ones are shown in Table 4.The relativistic change in the final state, including the relaxation of the electrons with a core-hole, Diff., is defined as: It has been reported [21] that the relativistic changes on the initial and final states are very close with each other, which means that the relativistic changes are dominant, already at the initial state.Since we want to establish the relativistic contribution to the absolute core level BEs, the changes were studied for the BE(ΔSCF) values.The Diff. values in Table 4 show that the relativistic 1s core BEs are larger than the non-relativistic values, although the relativistic effects lead to a very small increase in the core level BEs for these light atoms as expected; from 0.13 eV for B to 0.75 eV for F, thus increasing along the C→F series, as expected.
When properly accounting for the relativistic effects on 1s core orbitals, a clear distinction in between HF and the other xc meta-GGA and hybrid functionals pops up.worsens the situation, with an overestimation of more than 1 eV.However, despite the overall malfunction of Minnesota functionals in estimating 1s core level BEs, the M06-L meta-GGA is performing excellent.By including the relativistic contributions, i.e., achieving proper estimates, the deviation with respect the experimental values is, in average terms, below 0.1 eV, this is, within the XPS chemical accuracy.Thus, M06-L, known to be a rather accurate xc meta-GGA functional, even suited to simulate particularly complicated systems [24,25], is here reinforced in the sense that is also well-suited in obtaining an excellent description of 1s core electron energies of main group elements.

HF is overestimating
In light of the above-presented discussion, the obtained results for the six explored methods, with and without relativistic corrections, are overall good when computing ΔBEs of main group organic molecules, and predict absolute core level BEs with a small percentage deviation.However, when relativistic effects are considered, the meta-GGA M06-L xc functional stands out in between its fellows, showing chemical accuracy in both BEs and ΔBEs, with the concomitant advantage of being computationally more affordable than other tested hybrid functionals, and so poses itself as excellent for describing at the same time the thermochemistry of organic molecules and their core state energies.

Conclusions
Here we have explored the performance of HF, M06-L and M11-L meta-GGA Minnesota xc functionals, and B3LYP hybrid functinal plus M06 and M11 Minnesota xc hybrid functionals in predicting 1s core level BEs of a set of 20 molecules (39 core levels explored in total) extracted from the previous study of Takahata and Delano [15] containing samples of molecules containing C→F main group elements.The obtained results using ΔSCF methodology have been compared to the reported experimental references.This has been carried out in a non-relativistic fashion, yet the relativistic effects have been explicitly considered on isolated atoms, and, since are known to be independent of the particular chemical environment, added a posteriori on the obtained estimates.
The analysis yields that computed absolute core level BEs are, overall, and regardless of the method, close to the reported experimental values, with small percentage deviations, but sensible when considered in an absolute deviation, due to the inherent different construction of the explored functionals and methods.The BEs values have also been analyzed in terms of shifts, ΔBEs, where such building differences are cancelled, showing in overall terms a good performance; actually B3LYP and M06 hybrids, or M06-L and M11-L meta-GGAs could be used to get ΔBEs since their deviations, of the order of 0.2-0.3eV, are close to the XPS chemical resolution of 0.1 eV.
However, when relativistic effects are included, i.e., carrying out a properly realistic description of the core level relaxation, all methods show strong over/underestimations of the BEs, which can be as large as 1 eV.There is however one exception to this: the M06-L meta-GGA Minnesota xc functional, which, after applying the pertinent relativistic corrections, displays an overall mean deviation of the estimated BEs with respect the experimental measurements below 0.1 eV, this is, well within the XPS chemical resolution, and so, poses itself as an excellent choice when modeling core level BEs of main group organic molecules, plus being computationally less demanding than hybrid functionals, which in turn display a lower accuracy in BEs.

Table 1
Aside, the B3LYP display a less markedly negative offset, only significant for O 1s and F 1s .Last but not least, M06-L meta-GGA and M11 hybrid display very little deviations with respect ideality.It is remarkable that for O 1s and F 1s a trend of experimental values, whereas HF and M06 display a similarly negative offset.These deviations seem to progressively acute when going to O 1s and F 1s core levels.Let us pay attention to ΔBEs, listed in Table2, instead of absolute BE values.The overall statistical analysis reveals that, as expected, the above-commented over/underestimations cancel when comparing BE shifts, and so, the overall deviations for all the tested methods sensibly drop: The Minnesota functionals and B3LYP have MAE values of ~0.2 eV, whereas HF display a larger value of ~0.4 eV.Aside the marked underestimation observed for HF, B3LYP, and M06 and M06-L vanishes.It is noteworthy how by estimating BE shifts the mean error on DFT xc functionals approach to the XPS chemical accuracy of 0.1 eV.

Table 1 .
BE(ΔSCF) results for the 1s core orbitals analyzed, in underlined bold font, of the molecules described in the overall set.Experimental BE values are also given.A summary of the statistic analysis; ME, MAE, and MAPE, is reported.All values are in eV, except MAPE, in %.

Table 2 .
ΔBE(ΔSCF) results for the analyzed 1s core orbitals of atoms in molecules described in the overall set, and highlighted in underscored bold font.Experimental ΔBE values are taken from an arbitrary reference; CH 4 , NH 3 , H 2 O, and CH 3 F for C, N, O, and F. A summary of the statistic analysis; ME, MAE, and MAPE, is reported.All values are given in eV, except MAPE, in %.

Table 3 .
Summary of the regression adjustments of the 1s core level BEs and ΔBEs for each element and the overall set for the different explored methods, including the slope and the linear equation regression coefficient, R. Aside, offsets for BEs are given, in eV, as well as the offsets corrected for the relativistic contributions, Rel-Offset, as corrected by the values reported in Table3.

Table 4 .
Relativistic contributions (Rel.) and non-relativistic (Non-Rel.)results for BE(ΔSCF) B→F calculations.The difference (Diff.) between Rel. and Non-Rel. is also given.All values are given in eV.