Practical estimation of XPS binding energies using widely available quantum chemistry software

Chemical shifts observed in high‐resolution X‐ray photoelectron spectroscopy (XPS) spectra are normally used to determine the chemical state of the elements of interest. Often, these shifts are small, or an element is present in several oxidation states in the same sample, so that interpretation of the spectra is difficult without good reference data on binding energies of the likely constituents. In many cases, reference spectra taken from pure reference samples of the chemical components can aid the peak fitting procedure. However, reference materials are not always available, so that it becomes necessary to estimate the binding energies of likely components through quantum chemical calculations. In principle, such calculations have become much easier than in the past, due to the availability of powerful personal computers and excellent software. In practice, though, care needs to be taken in the approximations, assumptions, and settings used in applying such software to calculate binding energies.


| INTRODUCTION
With the discovery of new materials and new applications for those known, the need for surface and interface analysis has notably increased. Examples of the application of surface analysis are biology and bioengineering (characterization of immobilized antibodies, 1,2 interaction between replacement joints with bones and tissues 3,4 ); catalyst (characterization and behavior understanding 5,6 ); coatings; 7 corrosion; 8,9 defect analysis; 10 fouling; 11 and many others. Consequentially, the importance of X-ray photoelectron spectroscopy (XPS), its spectra interpretation, and the quantification obtainable from it has grown dramatically.
At NEXUS (the UK National EPSRC XPS Users' Service), we often need to make comparisons between XPS spectra from our users' samples and spectra we record from pure standard materials. In many cases, though, reference materials are not available or are not available in pure form or in known compositions, and in these cases spectra from advanced materials can be difficult to interpret. Typically, the interpretation of the results includes peak-fitting, the results of which can depend very sensitively on the quality of one's initial knowledge of the binding energies of the states included in the fit.
As is well known, XPS employs X-rays to eject core electrons of the elements present on the surface of the samples of interest. The kinetic energy of these electrons is measured and, because the energy of the incident x-ray is known, the binding energy of the electrons is obtained. The result of this analysis is a spectrum in which measured photoelectron intensity is plotted as function of the binding energy.
Each electron, from a given orbital of a given atom, will have a specific binding energy which allows the analyst to assign the peaks to specific elements and quantify those from the area of the peaks. The binding energies of the core electrons are affected (even though very mildly if compared with valance electrons) by the surrounding atoms. The consequent chemical shift is detectable in XPS and can be used to determine the chemical states of elements, although the shifts are normally quite small, and often peak fitting routines need to be used because they overlap in binding-energy range.
As a qualitative guide in the interpretation of the spectra, one can bear in mind that the binding energy of the core electrons generally increases when the valence electrons are involved in bonds and are therefore slightly more distant from the nucleus. The small amount of screening provided by the valence electrons is reduced further due to these bonds, and it is the core electron binding energies which are most important in XPS, so that the indirect effect of chemical bonding is seen in small binding energy shifts. This effect can be expected to increase with increasing electronegativity of the surrounding atoms. When the element is present only in a few chemical states, or the atomic environments give rise to very different binding energies, this process of distinguishing different functionalities is quite straightforward. However, in many cases, one has to deal with samples which contain elements in many different oxidation states as well as similar chemical environments which present similar binding energies. In these cases, the interpretation of the spectra is challenging, and there is the need to use reference spectra and data for the correct assignment of peak components. [12][13][14][15] Sometimes, reference spectra for the systems of interest are not available or, often, the literature is contradictory, and, in these cases, calculations which can predict XPS binding energies provide a good support for the spectral interpretation. In the first 2 decades in which XPS was widely used, quantum chemistry approaches were very time consuming and required expensive computers which were available only to a few.
Nowadays, with the rapid development of faster processors, these calculations require less time and can often be performed (depending on the size of the system of interest) on relatively inexpensive computer clusters or even on multiple-core processors in PCs or laptops.
Since the 1970s, many works have been published on the prediction of XPS binding energies using different approaches, typically using Koopmans' theorem or the so-called ΔSCF method, and other DFT methods. Quantum computational methods have successively been applied to calculations of Auger energies, 16 the understanding, and use of shake-up satellites structures 17,18 as well as multiplet splitting. 19 Besides using many different methods, the works in the literature employ many different software packages as well as home-made ones to implement those methods. Parallel with the development of these methods, there has been a trend towards a larger proportion of computational chemistry studies using a smaller number of comprehensive commercial software packages, for example GAUSSIAN. GAUSSIAN with its Gaussview interface is commonly used also by less experienced users and by those (like many busy surface analysts) who are not experts in computational chemistry. Most of the existing published work on binding-energy estimation, in contrast, is addressed to an expert audience rather than to surface analysts, who are the most concerned with these BE matters, and does not explain in simple terms how to perform calculations of the electron BE.
This work has the aim to guide the reader on the methods employed for the prediction of XPS binding energy and the ideas behind these methods. We illustrate this with 2 very effective methods, which can be run on PCs of only slightly above-average performance, using the popular software GAUSSIAN. Binding energy calculations for a list of molecules (containing elements of the second row of the periodic table) were carried out. These were compared with the experimental binding energies showing the accuracy of the methods.

| METHODS FOR THE CALCULATION OF ELECTRON BINDING ENERGIES
Calculating an electron binding energy means calculating the ionization energy from the respective core level. In order to do this, a quantum chemical approach should take into account the ejection of the core electron from its orbital (with energy E), the orbital relaxation upon ionization (R), and the electron correlation which will be changed after the ionization (ΔC). The binding energy can thus be expressed approximately as follows: There are other terms which can be added and which will be mentioned later on.
Depending on the complexity of the approach used, the BE can be approximated to different degrees of accuracy. It is important to emphasize that these calculations are normally applied to solid samples. This means that even the most accurate quantum chemical approach will not give exact absolute values of the BEs (because the BE in gas phase is higher than that in solid phase). However, very reliable results can be obtained for chemical shifts if this effect is constant across a set of cases of interest.

| Orbital energy
A simple way to approximate the binding energy is by taking BE = − E, as in several published works. 20,21 This first level of approximation here consists of applying Koopmans' theorem which states that the ionization energy of a molecular system is equal to the negative of the orbital energy. This approximation, though very crude, is quite useful and fits well enough with experimental data, when there are no significant final state effects (or rather those final state effects are constant). This means that this method can give a useful guide to evaluate the trends in initial state effects which are of inductive type and caused by the electron density over the atomic site where core ionization occurs (and therefore of analytical value in interpreting the surface chemistry). To obtain the energies of all the orbitals of interest of a given molecule, it is essential to first perform a geometry optimization of the molecule of interest. In geometry optimization, the arrangement in space of the atoms corresponding to a minimum of potential surface energy is found. Because a conformational change from a minimum to another is not expected to influence the binding energy of core electrons, it is generally not necessary to make sure that the minimum found is the global minimum.
Once this is done, we then need to consider individual orbitals, and therefore we need a way to identify them. A way to identify the orbitals easily is by using natural bond orbital (NBO) analysis. Employing NBO analysis, the orbitals which are associated almost exclusively with a single atom (for example core orbitals) are localized as natural atomic orbitals (NAOs); involving bonding or antibonding between 2 atoms is localized by using only the basis set atomic orbitals (AOs) of these atoms; and the Rydberg like orbitals are made orthogonal to one another. 22 Using a small basis set, it is also a good way to start and identify the molecular orbital. Once the calculation is run, orbitals can be visualized by opening the Gaussian checkpoint file (.chk) in Gaussview.
Alternatively, the output file (.out), if population analysis was required (pop = full), orbitals can be recognized on the output text file.
Depending on the accuracy of the method (Hartree-Fock or DFT), and even the choice of basis-set and the functional (in the case of DFT), different levels of accuracy may be expected. What needs to be taken into account is that, with the increase of the basis set, there is an increase in computational effort and time, and eventually this will become prohibitive.

| Relaxation energy
A better approximation can be obtained by taking BE = −E + R, thus including the relaxation effect, the response of the electrons to the formation of core hole. This can be taken into account by calculating the binding energy as the difference between the total energy of the ion and that of the neutral molecule through self-consistent field calculation (ΔSCF). 21 In this way, the BEs calculated account for both initial and final state effects. The latter are mainly caused by the relaxation.

It follows that a comparison between the results obtained by
Koopmans' theorem and ΔSCF will help us understand if the BE is affected only by the electron density of an atomic site or is also influenced by the capability of the electrons to screen the core hole. It is very important to know this when the BE is, for example, employed to probe the electron density of complexes. 23 As mentioned earlier, ΔSCF consists in calculating the binding energy of a specific electron of a specific atom in a molecule as the difference between the total energy of the ion (obtained by removing that specific electron from the molecule) and that of the neutral molecule. If, for example, we were aiming to calculate the C1s and O1s binding energy of methanol (CH 3 OH), we would need to perform a geometry optimization of the neutral molecule to find its total energy (at the optimized geometry) and then perform a single point calculation (on the same geometry optimized for the neutral species) for the molecule with a core hole in the C1s and a single point calculation for the methanol with a hole in the O1s. From the difference between each of the ionized molecules and the neutral one, the 2 binding energies are found. This means that, for every molecule of interest, n + 1 calculations are needed, where n is the number of electrons of which we wish to know the binding energy of. The way to calculate the energy of a molecule with a core hole will be explained in the computational method section.
ΔSCF gives reliable results and can be considered a good approximation even if we neglect electron correlation. This is because the relaxation energies upon core ionization are much larger than the change in correlation energies derived from them.

| Electron correlation
If we wish to take into account the electron correlation (BE = −E + R+ ΔC), then DFT methods need to be used. One approach is to use the unrestricted generalized transition state (uGTS) model. [24][25][26][27] The ionization energy is related to the energy of the Kohn-Sham orbitals 28 (which are in DFT the analogues of molecular orbitals in Hartee-Fock calculations) optimized for a transition state with fractional occupation number. Similar results are obtained through ΔKS. [29][30][31] Although, as mentioned, these methods give a treatment of both relaxation and correlation, they seem to be very functional dependent, so the "due diligence" required in ensuring the functional is fit-for-purpose (given the large literature for work on similar systems) does require some effort. Typically, the improved accuracy of these methods (compared with ΔSCF) is limited and not rewarding for organic compounds. However, the analysis of more complex systems containing metallic elements usually needs them.
Finally, to further improve the accuracy of the calculation of the binding energy, 2 factors can be taken into account: the relativistic correction and the zero-point vibrational energy correction.
Although these corrections could be made by means of calculation or approximations, 24,25,30,32,33 they are very small. While they could make a difference in high precision KS type calculations, they will typically not influence the Koopmans's theorem or the ΔSCF results.

| COMPUTATIONAL DETAILS
The binding energies of a set of 14 structures were calculated, using Koopman's theorem and the ΔSCF methods. ΔKS was also applied to some of the structures. All calculations were performed using the The first line of the input is the name of the check-up file which will be created (or from which input is read). The second line specifies the type of calculation and method selected and from where the starting geometry will be read (in this case from the coordinates given at the end of the input).
After the space, the next line is a comment line in which we are free to write information on the calculation. After the space, the 2 digits indicate, respectively, the total charge of the molecule and its spin multiplicity. After the molecule coordination needs to be introduced.
From the checkpoint and output file, we can distinguish the core orbitals of interest and see their energy.

| How to perform ΔSCF
To obtain the binding energy of a specific atom in a molecule from ΔSCF, 2 calculations are needed. The first is the one set out above. It gives us an optimized geometry, the total energy of the neutral species, and the opportunity to recognize and identify the positions of the orbitals of interest. The second calculation (the calculation of the energy of the core ionized molecule) requires 2 steps when performed with GAUSSIAN. The first step entails the swap of the core orbital of interest with the highest occupied molecular orbital (HOMO). The second step consists in the ionization; the calculation of the energy of a positively charged molecule. Because the HOMO has been swapped with the core orbital of interest, the electron will be taken from the latter. The steps are illustrated in Figure 1. An example of an input file is provided below.  The check point file from the first calculation was copied in a new file named xxx_Step2.chk. The geometry of the molecule will be read now from this (geom = check). The keyword "guess = only" requests that calculation terminate once the initial guess is computed, while "save" saves the generated guess in the checkpoint file. The keyword "alter" allows the exchange of orbitals. The orbitals to exchange will be indicated at the end of the input file for both α and β orbitals. In the example, the α orbitals will stay the same, while for the β, the orbital 1 (the lowest energy C1s) will be exchanged for the HOMO which is the orbital 17. The orbitals to select can be identified by the output file of the first calculation. It is also possible to visualize them in Gaussview.
In the second step of the second calculation, the molecule is ionized; thus, the charge is changed from 0 to 1 and the spin multiplicity from 1 to 2. Because there might be problems with the convergence/stability of a calculation with a molecule with a core hole, the keyword scf = DM, which calls for a direct minimization, is employed.
ΔKS calculations were explored, by means of the same approach using the functionals B3LYP and WB97XD. 35 The binding energies of the elements constituting the chemical species, illustrated in Figure 3, calculated by means of Koopmans's theorem and ΔSCF. The values were reduced by 5.47 eV (representing the difference between the free molecule and the solid state). The experimental binding energies for the related polymer form are also provided. 12 An approximation of the relaxation energies is given from the difference between the Koopmans's theorem values and those of ΔSCF. Note that not all the C1s species present in the model molecule are present in the polymer (and thus in the experimental data). Also, the C1s in the chlorine containing molecules was not always stable during ΔSCF calculation; thus, some theoretical values are missing An alternative approach for ΔSCF (not employed for this work) could be obtained by calculating the energy of the ionized species by using the equivalent core approximation. 36 An equivalent core species is used instead of the one we wish to ionize. For example, C1s ionized state of CO can be approximated as NO + .

| RESULTS AND DISCUSSION
Binding energy predictions were performed on the structures illustrated in Figure 2. These were used as models of the polymers from which the experimental energy was measured by Beamson  In Figure 3, the results obtained from both Koopmans's theorem Although the values obtained via Koopmans's theorem are higher than the others, because they exclude relaxation effects, the trend is the same as that observed experimentally and via ΔSCF. Table 1 also shows an approximate value for the relaxation energies given by the difference between the ΔSCF and the Koopmans's theorem's result.
Consistent results are obtained for each element, and the relaxation energy seems also the increase linearly with the increase in binding energy.