The local electron attachment energy and the electrostatic potential as descriptors of surfaceadsorbate interactions

Two local reactivity descriptors computed by Kohn–Sham density functional theory (DFT) are used to predict and rationalize interactions of nucleophilic molecules (exemplified by CO and H2O) with transition metal (TM) and oxide surfaces. The descriptors are the electrostatic potential, VS(r), and the local electron attachment energy, ES(r), evaluated on surfaces defined by the 0.001 e Bohr 3 isodensity contour. These descriptors have previously shown excellent abilities to predict regioselectivity and rank molecular as well as nanoparticle reactivities and interaction affinities. In this study, we generalize the descriptors to fit into the framework of periodic DFT computations. We also demonstrate their capabilities to predict local surface propensity for interaction with Lewis bases. It is shown that ES(r) and VS(r) can rationalize the interaction behavior of TM oxides and of fcc TM surfaces, including low-index, stepped and kinked surfaces spanning a wide range of interaction sites with varied coordination environments. Broad future applicability in surface science is envisaged for the descriptors, including heterogeneous catalysis and electrochemistry.


Introduction
The ability to predict and rationalize the interactions and reactions of extended surfaces is of great interest in numerous research fields spanning from heterogeneous catalysis, electrochemistry and corrosion, to nanotoxicity, biocompatibility and biodegradability. In the present paper we discuss an alternative direction to assess surface reactivity as compared to those traditionally used in the materials science and surface chemistry communities. We will show that two local properties, the local electron attachment energy [E(r)] and the electrostatic potential [V(r)], can be used for characterization of site-resolved Lewis acidity of transition metal and oxide surfaces when evaluated on contours of constant electron density, i.e. isodensity surfaces. The predictions made by these properties are compared to experimental and computed adsorption energies of two nucleophilic probe molecules, i.e. the Lewis bases CO and H 2 O. Both E(r) and V(r) are obtained at the 0.001 a.u. isodensity surfaces, an approach which has been found useful for predictions and understanding of the interactive behavior of isolated molecules as well as transition metal and oxide nanoparticles. [1][2][3][4][5][6][7][8] Through the analysis of these properties, we have previously been able to rationalize e.g. the enhanced catalytic behavior of Au nanoparticles; undercoordinated sites (e.g. corner sites) are associated with a depletion of electron density (known as a s-hole) and a corresponding local maximum in V(r). 4 The s-holes are more pronounced at the corner sites than at terrace sites; this explains the stronger interactions at corners, which are necessary for an increased catalytic activity compared to the inert, perfect surfaces of Au. We have found that the interactions at undercoordinated Au sites are similar in origin to those of e.g. hydrogen and halogen 9 bonding. A similar behavior is also found for Ag and Cu nanoparticles. 7 In light of these observations, we have introduced a new type of bonds, regium bonds, 7 that take place between the s-hole of a neutral Au, Ag or Cu site and a Lewis base, e.g. CO or H 2 O. Striving to explore similarities and differences within chemical interaction behavior over the periodic table, we will herein investigate if the above concepts are more generally applicable and can be used to also explain interaction tendencies of extended surface as well as of other metals and materials.
The use of DFT-based descriptors to characterize the reactive properties of surfaces is a common practice in the surface science community. An important example is the d-band model of Hammer and Nørskov, 10 which have been successfully employed to rationalize trends in adsorption energies on transition metal substrates and thereby, in the extension, to predict e.g. catalytic activity of metals and alloys. 11,12 Compared to the traditional descriptors, including the d-band model, both E(r) and V(r) have the advantage of being local in their nature. This means that they can not only be used to analyze the difference in reactivity/activity between different materials, but also the variation in reactivity over the same surface, i.e. the regioselectivity. In addition, we have found that the E(r) and V(r) are able to predict the local adsorption behavior for the interaction of Lewis bases with oxide surfaces; a process for which there are no generally applicable descriptor, including the d-band model and its extensions. [13][14][15] In the following we will demonstrate that local minima in E(r), E S,min , and local maxima in V(r), V S,max , on isodensity contours indicate sites susceptible towards interaction with electron donating molecules. We will also show that the relative magnitude of E S,min and V S,max closely reflects trends in adsorption energies.

Theory
We shall begin by readopting the formulism of the descriptors to fit the solid-state framework. For a polyatomic non-periodic, isolated system (e.g. a gas phase molecule), V(r) at position r is rigorously defined by where R A and Z A are the position and charge of the A:th atomic nuclei, and r(r 0 ) is the electron density function. The first of these terms comes from the ionic contribution to the electrostatic potential and the second term is the local Hartree potential. In the non-periodic framework, the vacuum potential of a free electron in its ground state is zero on an absolute scale (V abs vac = 0). Due to the periodic representation of standard planewave codes, infinity is implied in three dimensions. Therefore the electrostatic potential [V per (r)] obtained from a periodic calculation is ill-defined and will be arbitrarily shifted by a constant DV relative to the absolute scale: 16 V per (r) = V abs (r) + DV. (2) V per (r) will, furthermore, be shifted to different relative scales depending on the considered system; different atomic compositions, structures, orientations and vacuum distances results in different shifts. 17,18 For the comparison of V(r) [and E(r)] from one surface to another, a prerequisite is the usage of a common scale. Herein we will (approximately) shift V(r) back to absolute numbers by comparison to the potential V per vac at the mid-point of the vacuum region between two slabs (illustrated in Fig. 1a). [17][18][19][20] In summary: Hence V(r) will here be determined as: Turning next to the E(r) property, we define it as (based on its definition for isolated systems 2 ) to fit the solid-state framework. Here r i (r) and e i are the density function and eigenvalue of the i:th electronic state, including contributions from all virtual (unoccupied) one electron states going from the Fermi level, E F , up to a certain energy off-set, E 0 . On an absolute scale, and within exact generalized Kohn-Sham DFT (GKS-DFT), the logical choice of E 0 is zero. 2 This follows directly from Janak's theorem, 21 which stipulates that the addition of an (partial) electron to any one-electron state above e i = 0 corresponds to a net increase in energy. This would thus yield an unbound state that does not contribute to an attractive chemical interaction. Using the same reasoning as for the electrostatic potential, E 0 can be identified as the energy of an electron in the vacuum region of a periodic calculation, i.e. E 0 = E per vac . Since E = qV, and since the charge of an electron is À1 in atomic units, it follows that E per vac = ÀV per vac . In a periodic framework the evaluation of E(r) should also run over the entire k-space. Thus eqn (5) adapts to the form where n kpt is the number of k-points and w k is a symmetryweighting factor accounting for the number of irreducible representations of the studied k-point (see also Fig. 1b). For a spin-polarized calculation, the above summation runs over both the a and b spin-states. ‡ V(r) and E(r) are complementary properties to a certain degree; whereas the V(r) provides information on the electrostatic tendencies upon interaction, the E(r) property also contains components related to charge-transfer and polarization. E(r) can be partitioned into (shown here for the single k-point case) where t i (r) is the local orbital kinetic energy density of the i:th electronic state and defined as t i (r) = À1/2c i *(r)r 2 c i (r), 22 and V XC (r) is the exchange-correlation potential. While V XC (r) is approximately constant on an isodensity surface, t i (r) contains additional information as compared to V(r) and is associated with the energy change of the system upon addition of electron density (i.e. charge-transfer). For the nucleophilic analog to E(r), the average local ionization energy % I(r) of Sjoberg et al., 23,24 we propose the following expression for % I(r) in the framework of periodic Kohn-Sham DFT (and Hartree-Fock):

Surface models
To evaluate the descriptors, surface models where generated for different fcc metal substrates and for the transition metal oxides Cu 2 O and TiO 2 in the anatase and rutile phases. The structures were based on crystal parameters from Xu et al. 25 for the TiO 2 structures, from Yu et al. 26 for Cu 2 O and from Wellendorff et al. 27 for the fcc metals. For TiO 2 four-layer symmetric (1 Â 1) unit cell slabs were used, separated by a 28 Å vacuum distance, with the outermost layer being allowed to relax. A range of surface facets were studied including the anatase (100), (001) and (101), as well as rutile (100), (110)  For comparison of site-specific interactions with V(r) and E(r), and in the cases when experimental data was not available, CO adsorption energies where determined for the different structures using (2 Â 2) unit cells for the low index surfaces. The adsorption energies were calculated by: where E *CO , E * , and E CO are the electronic energies of the adsorbate-surface, free surface, and gaseous CO, respectively. CO was studied in a 15 Â 15 Â 15 Å 3 unit cell.

Computational details
All calculations were carried out with the Vienna ab initio simulation package (VASP). 28 The PBE xc-functional was used throughout. For the oxides, on-site +U Hubbard corrections were included as proposed by Dudarev et al. 29 The (U-j) parameters were taken from Xu et al. 25 for the TiO 2 structures and from Yu et al. 26 for Cu 2 O. Grimme's D3 dispersion with Becke-Johnson damping was employed for the CO adsorption studies. 30,31 During optimization the core states were represented by standard PBE-PAW potentials. Valence electrons were represented by plane-wave basis sets with an energy cut-off of 400 eV. 4 Â 4 Â 1 k-point G-centered meshes were used for the sampling of the k-space during the geometry optimizations employing the tetrahedron method with Blöchl corrections 32 for the oxides and the first order Methfessel-Paxton smearing method 33 for the metals. For the final energy evaluations and for the descriptor evaluations, 9 Â 9 Â 1 meshes were used for the oxides and 11 Â 11 Â 1 meshes for the metals. Hard PAW potentials and extended valences where used for the final energy and descriptor evaluations using a plane-wave cut-off of 800 eV. Molecular CO was treated by a single G-point. V(r) where generated in VASP (LOCPOT including only the ionic and Hartree potential contributions). E(r) was computed from the optimized wave function according to eqn (6) in the main article. Spin-polarization was considered for Ni(111). The V S (r) and E S (r) were visualized at the 0.001 a.u. (e Bohr À3 ) isodensity ‡ We note that in future applications, for non-zero Kelvin interactions, a Fermi-Dirac smearing of the electron occupation might be necessary to accurately capture the interaction behavior, thus leading to fractional occupation and corresponding weighting of the valence states. surface in VESTA. 34 V S,max and E S,min where evaluated at the same isodensity surface.

Results & discussion
Now to the central question: how do the E(r) and V(r) properties outlined above correlate with the Lewis acidic behavior of extended surfaces? The first feature that we shall explore is the local character of the descriptors. To investigate this, the properties are first mapped on the 0.001 a.u. isodensity contours of Pt(111), Cu 2 O(111) and rutile TiO 2 (110) shown in Fig. 2. These surface structures are suitable tests for the descriptors since the behavior upon adsorption of molecular species onto theses surfaces is experimentally well characterized.
On Pt(111) both the descriptors suggest that the sites on top of Pt atoms are the most susceptible to interactions with Lewis bases. This since E S,min and V S,max (also known as s-holes 9 ) appear and coincide on these sites, suggesting that the on top sites are the preferred sites based on both charge-transfer (i.e. low electron affinity as indicated by E S,min ) and electrostatic (i.e. strongly positive potential as indicated by V S,max ) arguments. Nucleophilic areas of low electrostatic potential are found at bridge and three-fold hollow sites in between atoms. Note again that E S (r) can only be used to characterize local Lewis acidity whereas V S (r) identifies both Lewis acidic (V S,max ) and Lewis basic (V S,min ) sites. The identified adsorption sites on Pt(111) agree well with the experimental findings that the Lewis basic CO adsorbs on top of Pt atoms. [35][36][37] Standard DFT LDA or GGA calculations do, however, notoriously rank the CO affinities of on top and hollow site incorrectly -known as the CO-puzzle. 38 Our results show that the on top site indeed is the preferred site with respect to electrostatics and the accepting of electron density from the sigma-orbital (HOMO) of CO. This in accordance with earlier studies, which attributed the preference for hollow sites to an incorrect leveling of the CO HOMO/LUMO energies compared to the band structure of Pt resulting in an overestimation of the back donation to the CO LUMO orbital at the hollow sites. 39,40 Fig. 2 also includes the examples of the stoichiometric and non-polar Cu 2 O(111) surface and the rutile TiO 2 (110) surfaces. For Cu 2 O(111), the ideal surface contains two possible Cu adsorption sites, the coordinatively unsaturated Cu (Cu CUS ) and coordinatively saturated Cu (Cu CS ) and two O sites (O CUS and O CS ). Among these, the Cu CUS site has been suggested as the favored site of interaction for Lewis bases, 41 which has been confirmed by DFT calculations of e.g. CO, H 2 S, H 2 O and methanol adsorption. [42][43][44] From the E S (r) map of the surface, we find that the Cu CUS sites indeed correspond to local minima, while the Cu CS sites are attributed lower Lewis acidity, i.e. lower electron affinity. The E S (r) indicates no tendency of O sites to be susceptible towards attack of Lewis bases. Hence, the predictions of the E S (r) are in line with the established knowledge. V S (r) outlines the same Lewis acidic surface characteristics as E S (r), but also ranks the O CUS as more Lewis basic (more negative potential) than O CS , again in agreement with the established picture. 44,45 Fig. 3 shows the same Cu 2 O surface but with 1/2 ML O CUS vacancies (V O CUS ), a commonly observed surface reconstruction. 46 This leaves behind a new type of undercoordinated Cu adsorption sites, Cu@V O CUS . Regardless, the most prominent E S,min and V S,max are still located at the same Cu CUS site as on the unreconstructed surface suggesting these sites remain the most favorable adsorption sites. Again, this is in line with previous understanding of e.g. H 2 O and methanol interactions onto this surface. 44,45 From the above, we note that our methods can be used for comparison to e.g. experimental STM (scanning tunneling microscopy) images, as well as in the initial assessment of reactive sites of complex nanostructured surfaces.
Concerning the rutile TiO 2 (110) surface, there are two possible Lewis acidic Ti sites, one six-fold (Ti 6f ) and on fivefold (Ti 5f ) coordinated. Similarly, there are two Lewis basic O sites, one three-fold (O 3f ) and one two-fold coordinated (O 2f ). E S (r) and V S (r) rank the Ti 5f site as more Lewis acidic (i.e. stronger electron affinity and higher positive potential) than the Ti 6f site, whereas the O 2f is more Lewis basic than O 3f , in line with experimental and previous DFT results for e.g. H 2 O, CO and methanol adsorption. [47][48][49][50]

View Article Online
In addition to rutile(110), we have also investigated the ability of the descriptors to identify CO adsorption sites on other TiO 2 surface facets including the rutile(101) and (100), as well as the anatase(100), (101), and (001) facets. TiO 2 is e.g. an important photocatalyst material, and these facets constitute the most commonly observed surfaces. [51][52][53] Fig. 4 shows the above-mentioned TiO 2 surfaces with CO adsorbed on the most favorable position and with the corresponding 0.001 a.u. E S (r) and V S (r) isocontour maps. Also for these surfaces it is found that the most extreme E S,min and V S,max correspond to the favored adsorption site.
From the above we conclude that E(r) and V(r) show promise as descriptors for the assessment of regioselectivity in surface Lewis acidity. However, the prevailing question is how well E(r) and V(r) can rank trends in adsorption energies for different substrates? Looking first at the series of Cu(111), (110) and (100) surface facets, Wang and Nakamura 54 have determined their H 2 O affinities using GGA DFT as (110) 4 (100) 4 (111). The Cu on top E S,min (V S,max ) for these surfaces are À8.35 eV (0.44 eV), À7.73 (0.34), and À6.53 (0.33) thus reproducing the trends in the adsorption energies. The corresponding surface maps are shown in Fig. 5C-E.
We have previously used the local variations of V S (r) to rationalize the enhanced catalytic activity of Au nanoparticles compared to its ideal crystalline surfaces. 4 A similar analysis is here included for Cu using both V S (r) and E S (r), covering nanoparticles as well as ideal, stepped and kinked extended surfaces (Fig. 5). Analogous to Au, we find that Cu nanoparticles display V S,max (and E S,min ) on top of the exposed atoms, with the most prominent V S,max and E S,min at corner sites, followed by edge and, lastly, terrace sites (Fig. 5B). Similarly, the extended surfaces have the V S,max and E S,min of largest magnitude (i.e. the most positive potential or strongest electron affinity) at corner sites followed by steps/edges, kinks (see definition in Fig. 5) and lastly terrace sites (Fig. 5C-G). Computed CO interaction energies follow the local descriptor values at the site of interaction. The correlations are strongest when comparing adsorption onto different sites for each substrate surface individually, but clear also when comparing adsorption onto all sites of all Cu substrates together (Fig. 5A). We note a tendency towards non-linear behavior going from the undercoordinated to the terrace sites. This could be seen as two different trends for the different classes of sites. We attribute this to a non-negligible effect of neighboring atoms, which could be captured by adding information from a larger coordination sphere to the interaction analysis.
One can also ask how sensitive the descriptors are to the choice of computational method? A test on Cu(100) using the BEEF-vdw 55 DFT xc-functional instead of PBE corroborate the understanding from isolated systems 2 that the V(r) property is rather insensitive to the functional while the absolute values of E(r) shifts depending on the choice of method. We find that V S,max at the on top site is 0.34 eV for both cases, whereas E S,min is À7.73 eV at the PBE level of theory and À7.45 eV using BEEF-vdw. The behavior of the descriptors using different functionals will be further evaluated in a follow-up study.
We next assess the experimental CO adsorption enthalpies of a series of fcc (111) transition metal substrates. The adsorption enthalpies are taken from a database of Wellendorff et al. 27 and includes the Ir, Rh, Ni, Pt, Pd and Cu surfaces. CO adsorbs on top of all but the Pd and Ni surfaces, where the hollow sites are preferred. 27 Using the E S,min value at the on top site from a 0.001 a.u. isodensity contour, we obtain an R 2 of 0.82 compared to the experimental adsorption enthalpies (Fig. 6 top). For CO, we obtain a slightly improved trend (R 2 of 0.89) when comparing to computed on top site adsorption energies at the PBE level § and extending the series with Au and Ag. Fig. 6 bottom include data for computed H 2 O adsorption energies 57 on fcc  (100) surfaces, the adsorption of CO occurs slightly distorted from the on top site, which is well reproduced by the location of the V S,max and E S,min . § The CO adsorption energy bias of the PBE functional with respect to experimental data has here been accounted for by adding the corrections suggested by Petersen et al. 56   View Article Online metals compared to E S,min showing a correlation with an R 2 = 0.82. We note that the main outlier in all the series above is Pt. The corresponding correlation of V S,max versus experimental data is poor R 2 = 0. 16. This suggests that the differences in CO affinity between the metal surfaces is determined by differences in charge-transfer capacity rather than electrostatics, which is in agreement with the success of the d-band model that downplays electrostatics in favor for orbital-band mixing.
Although the examples presented here represents a limited initial study, they demonstrate that the E(r) and V(r) descriptors can be used to predict and rationalize molecular interactions at both metal and oxide surfaces and have great potential for future applications. Results for other adsorbents including e.g. hydrogen atoms, hydrocarbons, water, ammonia and other nitrogen centered adsorbates will be presented in a follow-up study. For future directions, an interesting path is to further explore the Lewis basicity of metal and oxide surfaces. This could, as touched upon in the present study, be evaluated by the negative regions of V(r) in combination with the Lewis basicity analog to E(r), i.e. the % I(r) of Sjoberg et al. 23,24 This will be a focus of a future study.

Conclusion
In conclusion, we have shown that the local electron attachment energy E(r), as well as the electrostatic potential V(r), mapped on isodensity contours can reproduce the regional variations in Lewis acidity of metal and oxide surfaces. The relative magnitude at the site of interaction can, moreover, be used to correlate adsorption energy trends for CO and H 2 O onto fcc metal surfaces. A general applicability of the two properties for various types of surfaces and applications, e.g. in the fields of heterogeneous catalysis and electrochemistry, is envisioned. For instance, the approach shows great promise for the analysis of special sites/motifs on complex nanostructured material surfaces that may correspond to active sites for catalysis or sites that are susceptible to corrosion as well as seeding sites for electrochemical deposition. It is, moreover, likely that the descriptors can be used to screen for new, tailored materials of e.g. high catalytic activity, or for the comparison to experimental STM images as well as in the prediction of nanotoxicity by the differentiation of reactivity of nanoparticles. We furthermore recognize the applicability of our methods as a guiding tool for surface chemistry modeling, providing swift rationalization of the orientation and positioning of interacting molecules.

Conflicts of interest
There are no conflicts to declare.