CO adsorption on the GaPd( % 1 % 1 % 1) surface: a comparative DFT study using diﬀerent functionals

CO adsorption on the polar ( % 1 % 1 % 1) surface of the intermetallic compound GaPd is examined within ab initio methods using an all-electron full-potential electronic structure approach. Comparison between the PW-LDA, GGA-PBE, GGA-RPBE, GGA-revPBE, and hybrid HSE06 functionals is considered through bulk, clean surface and CO adsorption calculations. The choice of the functional is found to have a strong influence in the description of single CO adsorption on the surface model proposed in literature. As expected from the so called ‘‘CO adsorption puzzle’’, diﬀerences in the obtained results demonstrate that classic LDA and PBE functionals can only partially describe the complex CO adsorption bonding scenario on a surface containing transition metal elements (in this case Pd atoms), where the energies of the substrate–adsorbate electronic states are shifted, yielding important differences in the absolute values of the adsorption energies, vibrational frequencies and surface–adsorbate interaction. So far the hybrid functional HSE06 correctly retrieves all the tendencies observed experimentally as confirmed by comparing our first-principles results to experimental findings.


I. Introduction
The intermetallic compound GaPd is experiencing great interest due to its catalytic properties in the semi-hydrogenation of acetylene, an essential cleaning step for the polyethylene production. [1][2][3] Recently, different studies about the structure determination and chemical reactivity of the polar surface GaPd(% 1% 1% 1) have been published. [4][5][6][7][8][9] The work of Rosenthal et al. is an experimental study that joins a variety of surface science techniques. 4 In particular, results of thermal desorption spectroscopy (TDS) of adsorbed CO molecules indicated that different surface terminations are obtained depending on the preparation temperature (670 K or 870 K). Scanning tunneling microscopy (STM) images revealed flat terraces separated by a unique step height. Lowenergy electron diffraction (LEED) together with X-ray photoelectron diffraction (XPD) showed no signs of surface reconstruction or segregation. It was proposed that the surface prepared at 870 K can be explained with an atomistic model derived from a bulk truncation containing in the topmost atomic layer a single protruding Pd atom per surface unit cell. For the surface prepared at 670 K, a denser surface model that could contain either Pd-rich or Ga-rich terminations was proposed. Another study, a joint application of experimental LEED, STM, and plane-wave quantum chemical calculations using the PBE-GGA functional, 10 was performed by Prinz et al. 5 Calculated relative surface energies and qualitative comparison between experimental and simulated STM images in addition to LEED results, lead the authors to the conclusion that the GaPd(% 1% 1% 1) surface can be explained with a model containing only one protruding Pd atom per surface unit cell, located in the center of 3 Ga atoms, in agreement with the previous findings. 4 However, a different surface model was proposed in the work of Krajčí et al. 6 In this study, the determination of the structure of the GaPd(% 1% 1% 1) surface was performed exclusively by plane-wave quantum chemical calculations using the PW91-GGA functional. As the previous publication by Prinz et al. it includes relative surface energy calculations and qualitative comparison of the electronic structure with STM images and photoemission spectroscopy. However, their conclusion points towards a Ga-terminated surface model and not to a Pd-terminated model as in the previous works. The publication of Krajčí et al. includes a study of the CO adsorption in different surface models, 6 but as in previous ref. [7][8][9] calculations were performed at the GGA level of theory, thus not addressing the ''CO adsorption puzzle''. [11][12][13][14][15] For understanding processes occurring at solid surfaces, the atomic structure determination must be complemented by the analysis of the electronic structure, with both depending on the interplay between the chemical species present at the surface at their corresponding coverages. Thus, to comprehend and predict the behavior of solid surfaces, an atomistic approach is vital. The density functional theory (DFT), offering good accuracy at low computational cost, has become the main tool in the last decades for understanding properties of molecules, aggregates, and materials on the atomic scale. 16,17 DFT is in principle exact. 18 However, since the exact form of the exchange-correlation density functional is not known, approximations must be introduced. 19 The most widely used approximations are the local density approximation (LDA) and the generalized gradient approximation (GGA). In complex and more realistic systems, the choice of the right functional describing the exchange and correlation interactions has a strong impact on the accuracy of the obtained results, sometimes leading to considerable errors in the predicted energies, misleading final interpretations. [11][12][13][14][15][20][21][22][23][24][25][26] Particularly, DFT within the Perdew-Wang parametrization for the LDA (PW-LDA), 27 the Perdew-Burke-Ernzerhof (PBE) functional for the GGA, 10 and the modified PBE functional, namely the RPBE, 28 fail to predict the correct adsorption site for the CO molecule on most of the surfaces of transition-and noble-metal elements. [11][12][13][14][15] The underlying problem, known as the ''CO-adsorption puzzle'', is due to the inherent lack of accuracy of both LDA and GGA to describe at the same time the electronic structure of small molecules, such as CO, the electronic surface structure of metals containing d (and/or f) states, as well as the interactions of the surface and the molecule. This leads to a wrong positioning of the orbitals of the CO molecule upon surface adsorption, causing an overestimation of the adsorption energies. 11,14,15 Additionally, the error also depends on the atomic structure of the adsorption site, which leads to wrong sitepreference, i.e. wrong energetic ordering of the different adsorption sites. 12,14,15,29,30 Thus, the quest for an appropriate DFT method that includes a treatable exchange-correlation term to solve the problems linked to the ''CO adsorption puzzle'' for intermetallic compounds containing transition-metals at their surfaces must be faced before any comparison with available experimental data. In the present work, a comparative DFT study of the CO adsorption on the Pd 1 -terminated model 4,5 of the (% 1% 1% 1) surface of the intermetallic compound GaPd is conducted using different exchange-correlation density-functional approximations. This case study constitutes an excellent opportunity to unveil the effects coming from both surface and sub-surface transition-metal elements on the adsorption of the small, yet very important CO molecule.

II. Computational details
Calculations are conducted using the all-electron full-potential code FHI-aims. 31 Effects of different approximations in exchangecorrelation functional are explored by comparing results obtained with (i) the local density approximation (LDA), (ii) the generalized gradient approximation (GGA) and (iii) the hybrid Hartree-Fock/ GGA approximation. Atomic relaxations, electronic structure, and energetics are calculated using the Perdew-Wang parametrization for the LDA (PW-LDA), 27 the GGA Perdew-Burke-Ernzerhof (PBE) functional, 10 the modified PBE functional RPBE, 28 the revised PBE functional revPBE, 32 and the hybrid Hartree-Fock/GGA-PBE Heyd-Scuseria-Ernzerhof functional (HSE06). 33 In the latter, the fraction of the Hartree-Fock exchange a and the screening parameter o are set to 0.25 and 0.11 bohr À1 , respectively, as recommended. 33 FHI-aims employs numeric atomic orbital (NAO) basis sets. Final results are obtained using ''tight'' numerical settings, basis set level up to tier 1 for Ga and Pd, and up to tier 2 for C and O species. 31,34,35 Relativistic effects are considered for both core and valence electrons employing the scaled zeroorder regular approximation (ZORA). 31,36,37 Spin polarization is considered in both bulk and surface calculations, where a fixed quantization (collinear) axis is used to define the spin density. Atomic relaxations are conducted by calculating first energyderivatives (forces) with respect to the nuclear coordinates, and their minimization using the trust radius method of the Broyden-Fletcher-Goldfarb-Shanno BFGS algorithm. 38,39 For the bulk material, also unit cell size/shape relaxation is considered. In all our calculations, forces are converged to values better than 10 À4 eV Å À1 , and total energies and charge densities are converged to values better than 10 À6 eV and 10 À5 e Å À3 , respectively. A well-converged G-centred mesh of 20 Â 20 Â 20 k-points is used to sample the Brillouin zone in bulk calculations. For the slab calculations an 8 Â 8 Â 1 k-point mesh was used. Tests in bulk calculations -for all here tested functionals -using an 8 Â 8 Â 8 mesh and a 20 Â 20 Â 20 mesh result in energy differences smaller than 10 meV, while no significant differences are found in the calculated bulk and surface atomic and electronic structures.
PBE and HSE06 calculations including both the Tkatchenko-Scheffler (TS) correction 40 for van der Waals interactions and the modified Tkatchenko-Scheffler (TS surf ) correction that includes short-range screening effects, particularly developed for the study of surface-adsorbate systems, [41][42][43] have been considered. The results obtained with the HSE06 (PBE) functional show that the dispersion corrections of both TS and TS surf contribute with E0.1 eV (E0.2 eV) to the final adsorption energies of the studied adsorbate-substrate system. Furthermore, their contribution to the CO adsorption energy differences between the most favoured adsorption sites T 1 and H 1 is o10 meV (o50 meV), for the results obtained with the HSE06 (PBE) functional. Differences between TS-and TS surf -corrected adsorption energies for the studied adsorbate-substrate system are o10 meV (o20 meV) for the results obtained with the HSE06 (PBE) functional. Since for all the here tested functionals, the CO adsorption energy difference between the sites T 1 and H 1 is E1 eV (see Table 4), van der Waals dispersion corrections are found to play only a minor role in the CO adsorption on the GaPd(% 1% 1% 1) surface and are not further considered in the present study. These results are in line with the fact that CO is a very small and stable closed-shell molecule and that Pd is the element that presents the smallest correction for the many-body collective response (screening) of the C 6 coefficients of different metallic surfaces. [41][42][43] The C-O stretching vibrational frequency of the free molecule (gas phase) together with the vibrational frequencies for CO adsorbed at different surface sites are calculated in the harmonic approximation by diagonalizing the Hessian matrix obtained by finite differences of analytic forces, using ''tight'' numerical settings and basis set level up to tier 1 for all atomic species. The motion of the C and O atoms, as well as the two topmost atomic layers Pd 1 and Ga 3 is considered, with displacement steps of 0.002 Å for each atom along each Cartesian coordinate. Lower-lying atomic layers of the slab are kept frozen at their corresponding relaxed positions.
Chemical bonding analysis in real space is carried out using the concept of electron localizability indicator (ELI). 44 The ELI in the ELI-D representation 45,46 is computed by the help of an interface to the FHI-aims code. For both the electron density (ED) and the ELI-D, the attractor locations, basins and the number of electrons contained inside the basins are determined by the DGrid package. 47

III. GaPd crystal structure and bulk calculations
The intermetallic compound GaPd crystallizes in the cubic FeSi type of structure, space group P2 1 3 (No. 198, Pearson symbol cP8), containing 8 atoms per unit cell. 48 At 295 K and atmospheric pressure, the lattice parameter of GaPd is a = 4.89695(6) Å, with a volume V = 117.4 Å 3 . 48 The first coordination shell of the Pd (Ga) atoms consists entirely of Ga (Pd) atoms. The shortest interatomic distance corresponds to a Pd-Ga distance equal to 2.5399(2) Å, located along the threefold axis, and the closest Pd-Pd distance is 3.0084(1) Å. 48 Because of the absence of inversion symmetry, two possible enantiomorphic forms of the GaPd compound exist, namely form A and form B. 4,48 The calculations of the atomic and electronic structure of the GaPd compound are initially based on lattice parameters and atomic coordinates from its enantiomorphic form B, as described in ref. 4 and 48 (but accidentally labelled form A in ref. 48). Fig. 1 shows the structure model of form B. Each enantiomorph possesses a characteristic but inverse stacking sequence of four non-equivalent atomic layers along the polar [111] direction. Following the notation developed by Rosenthal et al.,4 in the case of form B, the sequence of atomic planes is Pd 1 /Ga 3 /Pd 3 /Ga 1 as shown in Fig. 1.
The structural stability of crystals typically correlates with the lattice parameters. 49 The structural stability of solids is characterized by the enthalpy of formation DH f , which is linked to the cohesive energy and the elasticity of a given compound. 50 For the intermetallic compound GaPd, a rough approximation of DH f at 0 K, no pressure, and without considering vibrational contributions -here called ''formation energy'' DE f -can be obtained by calculating the difference between the total energy of bulk GaPd unit cell, and four times the respective bulk energies of elemental Ga and Pd. Table 1 presents the lattice parameter a and the formation energy DE f of GaPd, obtained with different functionals. The relative errors (in %) of the calculated values with respect to the experimentally obtained ones (at room temperature) are given in parenthesis, where a positive (negative) sign indicates an over(under)estimation. To compute DE f , the stable a-modification of Ga and the Cu type of structure of Pd, were considered. [51][52][53] As expected, 49 the LDA calculation underestimates the lattice parameter of GaPd with respect to the experimental value, while GGA calculations overestimate it. Bulk calculations using the hybrid functional HSE06 were also conducted. The corresponding optimized lattice parameter is the closest (+0.6%) to the experimental value. The choice of the functional has also an important effect on the obtained values of DE f . Although the latter is negative in all cases, which is in agreement with the high melting point of GaPd (1338 K), 54 both LDA and GGA calculations overestimate DE f compared to the experimental value of DH f = À5.96 eV. 55 The RPBE functional gives the highest deviation from the experimental value (+16.38%). In general, all GGA-based calculations give similar results, with the LDA value being this time closer to the experimental DH f (+5.89%) than the GGA values. Although the nearest value to the experimental DH f is the one obtained with the hybrid HSE06 functional (À5.87%), this is not much better than the LDA result.
Another important property that can be correlated to the structural stability of a solid is its electronic structure. Fig. 2 presents the bulk density of states (DOS) of GaPd calculated with different functionals, with the lattice parameter relaxed for each functional (see Table 1). Table 2 summarizes the main features of the obtained DOS (defined in Fig. 2), including the opening of the DOS between EÀ5 and EÀ7 eV (position of its center O C and length O L ), the peak P d of the Pd 4d band, the Pd 4d bandwidth W d , and the DOS at the Fermi level n(E F ). In line with previous results, in the valence region the most important partial contribution to the total DOS of GaPd is given by the Pd 4d states, that are shifted below the Fermi level in comparison to pure Pd. 1,[56][57][58][59] Although the GGA functionals perform similarly, e.g., |DP d | r 0.11 eV, |DW d | r 0.18 eV, the differences between the DOS calculated with LDA and GGA functionals are significant, e.g., DW LDA-RPBE d = 0.75 eV, DP LDA-RPBE d = À0.50 eV. However, a PBE calculation at the geometry optimized with the PW-LDA functional (denoted below PBE*) removes almost completely the observed Fig. 1 Projected crystal structure of the enantiomorph B of the intermetallic compound GaPd. 48 The red box depicts the unit cell, gray (green) spheres represent Pd (Ga) atoms.  Here, the obtained differences cannot be explained as mainly a geometric effect. The observed differences are due to the self-interaction error present in the LDA and GGA functionals, which is partly removed in hybrid functionals including the HSE06. 14,33,60-65 As a result, the latter provide a better description of the electronic structure of molecular systems, and have been demonstrated to yield better alignment of the metallic band in solids, but can overestimate its width. 14,33,[60][61][62][63][64][65] This can explain the difference between the HSE06 and LDA/GGA positions of the Pd 4d band in the bulk GaPd DOS. Since the chemical reactivity of a surface depends on the electronic structure of both the adsorbed molecular species and the surface, the obtained results suggest that the predicted chemical properties of the surfaces of GaPd will be sensitive to the approximations used in the employed exchange-correlation functional.
Following the bulk stacking sequence of atomic planes Pd 1 /Ga 3 / Pd 3 /Ga 1 along [% 1% 1% 1], the construction of the Pd 1 -terminated model of the GaPd(% 1% 1% 1) surface is achieved by simple bulk truncation, resulting in a (O3 Â O3)Rot301 surface unit cell. 4,5 Fig. 3 shows a top view of the Pd 1 -terminated model. The topmost atomic layer (Pd 1 ) consists of only one Pd atom per surface unit cell. The second (third) layer Ga 3 (Pd 3 ) consists of three Ga (Pd) atoms per surface unit cell -connected with dotted lines in Fig. 3. The simulation of the Pd 1 -terminated model is achieved through asymmetric slabs. In this configuration, the slab is composed of a certain number N F of atomic layers that are kept fixed at the optimized bulk atomic positions (lying at the bottom of the slab), while another number N R of atomic layers on top of them (including the surface of interest) can fully relax. The thickness of the vacuum region separating two consecutive slabs is set to 20 Å. The area of the surface unit cell calculated with the LDA, PBE, and HSE06 functionals is 40.6 Å 2 , 42.6 Å 2 , and 42.3 Å 2 , respectively. Calculations conducted with the PBE functional and an 11-layer-thick slab consisting of different N R and N F layers -keeping the lattice parameters of the supercell fixed -show that the total energy of the slab changes with the ratio N F /N R , decreasing with N R . The most noticeable differences are between the here tested extreme values N F /N R = 1/10 and N F /N R = 8/3 (DE slab = 104.5 meV). To determine the thickness of the slab sufficient to simulate the Pd 1 -terminated model, surface energy differences and layer-resolved surface density of states (sDOS) calculations were compared for an 11-, a 20-, and a 29-layerthick slab using the same conditions as described above, but (i) keeping constant the number of fixed layers lying at the bottom (N F = 4) or (ii) keeping constant the number of atomic layers   that can fully relax (N R = 7), all slabs having identical bottom termination. Thus, the here calculated surface energy differences contain only contributions from the surface of interest (top), since the surface energy of the ''frozen'' bottom termination is the same for all slabs and cancels out when calculating surface energy differences. The calculations show differences smaller than 1 meV Å À2 (4 meV Å À2 ) for the surface energy difference between the 20-(11-) and 29-(20-)layer-thick slabs. Furthermore, only negligible changes of the calculated sDOS with the slab thickness were found, all slabs reaching bulk-like behavior in the central layers. Therefore, the 11-layer-thick slab with N F = 4 at the bottom and N R = 7 at the top (including the surface of interest) is used in all surface calculations. Fig. 4 shows the sDOS of the Pd 1 topmost layer of the clean surface model calculated with different functionals. As for the bulk calculations, the differences between the sDOS calculated with the PBE, revPBE and RPBE functionals, at the geometry optimized with the PBE functional, are negligible (smaller than 0.04 eV for O C , O L , P d and W d ). However, the PW-LDA calculation at the geometry optimized with the PBE functional presents a small shift (0.3 eV) towards lower energies of the opening of the sDOS and the Pd 4d band, in comparison to the PBE results. More significantly, the PBE sDOS shows a shift to higher energies (|DO C | E |DP d | = 0.9 eV) and a reduction of the Pd 4d bandwidth (|DW d | E 1.0 eV), when compared to the HSE06 sDOS, in line with the bulk calculations. The bottom panel in Fig. 4 presents the partial contributions to the total HSE06 sDOS (valence region) of the Pd 1 topmost layer. As expected from a Pd-terminated surface, the analysis shows that its main structure -between À6 and À1.7 eV -is primarily dominated by the Pd 4d states. The total area below the sDOS in this energy window is 4.46 states per cell, where 95% of it corresponds to Pd 4d states. The area under the sDOS between À7 and À6.8 eV is 0.12 states per cell and is also dominated by Pd 4d states, but to a lesser extent (69% Pd d states, 12.5% Pd p states and 18.5% Pd s states). Between À11.6 and À8.4 eV, the number of states in the sDOS is negligible. Fig. 5 shows the HSE06 sDOS convoluted with a Gaussian function [66][67][68] (De = 0.25 eV) and compared to the He II ultraviolet photoelectron spectra (UPS) of the GaPd(% 1% 1% 1) clean surface recorded at room temperature of the sample prepared at 870 K. 4 To estimate the contribution of the sub-surface atomic layers, the convoluted surface DOS is calculated considering (a) only the topmost atomic layer Pd 1 , (b) both surface and subsurface atomic layers Pd 1 and Ga 3 , and (c) the three topmost atomic layers Pd 1 , Ga 3 and Pd 3 (emphasized by colored areas in Fig. 5). Moving towards lower energies, the structure of the convoluted sDOS presents two small peaks at 0.2 eV (P 1 ) and À0.7 eV (P 2 ), the two main peaks of the Pd 4d band at À3.0 eV (P 3 )   4 In the top panel, the blue curve (the envelope of the white area denoted ''a'') is the convoluted DOS due to only the top-most layer of the slab. The envelope of the red patterned area (area b) is due to the top and the subsurface layers, while the envelope of the blue patterned area (area c) is due to the three top-most layers of the slab. The inset shows the same three sDOS in the region above the Fermi level (E F ). and À4.3 eV (P 4 ), and a shoulder-like peak at À5.6 eV (P 5 ). After an opening at À6.3 eV, a peak at À7.0 eV (P 6 ) appears together with two subsequent shoulder-like peaks. As can be seen in Fig. 5, there is an excellent agreement between the convoluted HSE06 sDOS and the experimental data.

V. CO adsorption on the
Based on the results and discussion presented above, the CO adsorption on the Pd 1 -terminated model of the GaPd(% 1% 1% 1) surface was thoroughly studied applying different functionals. Several starting geometries of the molecule were considered taking into account the atomic structure of the three topmost layers of the Pd 1 -terminated surface model. The CO molecule was allowed to fully relax starting from different positions, 3 Å above the initial adsorption site. The initial C-O bond length was set to the calculated value for CO in gas phase (see Table 4). Fifteen crystallographically unique possible adsorption sites can be identified. Fig. 6 summarizes the initial adsorption sites considered in this study: top (T), bridge (B), and hollow (H) sites. In the T i sites (i = 1 to 3), the CO molecule is initially adsorbed atop one Ga or one Pd atom of the substrate. In the case of B j sites ( j = 1 to 7), the CO molecule initially forms a bridge between two atoms of the substrate, while in the case of H k sites (k = 1 to 5) the CO molecule is initially located between three substrate atoms. Total adsorption energies E ads (a b ) for the different adsorption sites a b were calculated using where E slab+CO (a b ) is the total energy of the substrate-adsorbate system, i.e., slab plus the CO molecule adsorbed at the site a b , fully relaxed (a = T, B, H and b = i, j, k, respectively), E slab is the total energy of the slab before CO adsorption (clean surface), and E CO is the total energy of the CO molecule in the gas phase.
In calculations based on the HSE06 functional, the most favoured adsorption site is T 1 with CO adsorption energy of 1.07 eV (103.24 kJ mol À1 ). In this site, the CO molecule relaxes perpendicular to the surface with the C atom pointing towards the topmost Pd atom of the substrate (d(C-Pd 1 ) = 1.98 Å, see Table 4). When the CO molecule is initially above the sites B 3 , B 6 , B 7 , and H 5 it also relaxes to a position above T 1 . In the case of the starting position above T 2 , the CO molecule desorbs upon relaxation. Similarly, the CO molecules initially above the hollow sites H 2 and H 3 desorb. On the other hand, the CO molecule initially above T 3 moves towards the hollow site H 1 , relaxing perpendicular to the surface, with the C atom pointing towards the substrate (d(C-Ga 3 ) = 2.88 Å, see Table 4). This is also the case for the CO molecule initially relaxed from above the bridge sites B 1 , B 2 , B 4 , B 5 , and hollow sites H 1 and H 4 . The site H 1 is composed of a triangle of Ga atoms with a central Pd atom at 2.04 Å below them. In spite of the presence of the accessible Pd atom, the surface Ga atoms make this site unfavourable (E ads = 0.05 eV = 4.82 kJ mol À1 ) for the adsorption of the CO molecule.
Summarizing, the HSE06 functional predicts that the CO molecule preferentially adsorbs only at the site T 1 , perpendicular to the substrate, with the C atom pointing towards the topmost Pd atom. Furthermore, no CO-Ga interaction is observed. These results are in perfect agreement with the experimental results obtained with the Fourier transform infrared spectroscopy (FTIRS) and the reflection-absorption infrared spectroscopy (RAIRS). 1,4,7,[56][57][58][59] The CO adsorption on the GaPd(% 1% 1% 1) surface at room temperature (50 mbar dosage) resulted in the appearance of only one band that was assigned to CO adsorbed on Pd in the on-top position. 1,4,7,[56][57][58][59] Moreover, the smaller adsorption energy for CO at the site T 1 of the Pd 1 -terminated model of the GaPd(% 1% 1% 1) surface in comparison to CO adsorption on the Pd(111) surface (1.8 eV), 69 is also in line with the thermal desorption spectroscopy (TDS) results, where the CO desorption from the Pd(111) surface needs a higher energy (desorption at 450 K) than from the GaPd(% 1% 1% 1) surface (desorption at 215 K). 4 In order to correlate the experimentally obtained desorption temperatures 4 with the computed adsorption energies, an estimation of the activation energies of desorption E d is given in Table 3. Considering the results and conditions of the TDS experiments for CO desorption on the GaPd(% 1% 1% 1) and the Pd(111) surfaces given by Rosenthal et al. 4

and the equation of Redhead for thermal desorption: 70
where R is the gas constant, n 1 is the frequency factor for desorption, T P is the temperature at which the desorption rate is maximum and b H is the heating rate, a reasonable estimation of E d is obtained. We have used n 1 = 1 Â 10 13 s À1 , which is a typical value for the estimation of E d (with n 1 = 1 Â 10 12 s À1 and n 1 = 1 Â 10 14 s À1 , the results are essentially the same). In agreement Fig. 6 Top view of the Pd 1 -terminated model with the initial adsorption sites considered in this study. Top, bridge, and hollow sites are indicated by black (1 to 3), blue (1 to 7) and red (1 to 5) numbers, respectively. Color code is as in Fig. 1 and 3. with the computed CO adsorption energy difference (0.7 eV), the estimated difference between the activation energy of desorption of CO on GaPd(% 1% 1% 1) and on Pd(111) is 0.6 eV, signifying a stronger interaction between CO and the surface for Pd(111) than for GaPd(% 1% 1% 1). The here obtained HSE06 CO adsorption energies for the sites T 1 and H 1 are somewhat smaller than previous values obtained with semi-local exchange-correlation functionals. Adsorption energies of 1.19 eV (PW91-GGA), 6 1.37 eV (PBE-GGA), 7 and 1.59 eV (RPBE-GGA) 8 have been calculated for the T 1 site, and of 0.69 eV (PW91-GGA), 6 0.56 eV (PBE-GGA), 7 and 0.49 eV (RPBE-GGA) 8 for the H 1 site. All the previous calculations were performed at the GGA level of theory using the pseudopotential plane-wave method, 6,8 and a mixed Gaussian and plane-wave approach. 7

(a) Influence of the exchange-correlation functional on adsorption energies and atomic relaxations
To gain further insight into the influence of the exchangecorrelation functional on the description of CO adsorption on the GaPd(% 1% 1% 1) surface, LDA-and GGA-based adsorption energies were calculated for CO above the sites T 1 and H 1 . In both cases, the CO molecule was allowed to fully relax starting from the final geometries obtained with the HSE06 functional. Table 4 presents the C-O bond length and the interlayer distances of the topmost layers of the Pd 1 surface model before (clean surface, CO in gas phase) and upon CO adsorption on T 1 and H 1 , as well as the respective adsorption energies, obtained with LDA, GGA, and hybrid functionals. The adsorption energies E ads * are calculated at the HSE06 slab geometry, while the adsorption energies E ads are calculated using slabs built with the lattice parameter optimized with each functional (see Table 1). For the here tested LDA and GGA functionals, the differences between E ads and E ads * are smaller than 0.02 eV for CO adsorption on T 1 , while they can account for differences up to 0.08 eV (LDA) for CO adsorption on H 1 . The results obtained with the HSE06 functional show that, the main effect of the atomic relaxation of the clean surface is the inward relaxation of the topmost Pd atom (0.18 Å), while lower lying atomic planes exhibit much smaller relaxations (o0.04 Å). Upon adsorption of CO on top of a Pd atom (site T 1 ), the latter relaxes outward again, and the surface interlayer distances become almost the same as the ones from the initial bulk positions (the differences are smaller than 0.02 Å). This indicates a strong interaction between the CO molecule and the topmost Pd atom. On the other hand, when CO relaxes above the site H 1 , the interlayer distances relax back to the values calculated for the clean surface (before CO adsorption), as if the CO molecule and the surface acted as two independent systems, corroborating the findings discussed above. Similar to the HSE06 results, the LDA and GGA calculations show an inward relaxation of the topmost Pd atom as main effect of the atomic relaxations on the clean surface. However, deviations in the interlayer distances up to 0.08 Å for LDA and 0.04 Å for the GGA functionals are obtained, when compared to the HSE06 results. After CO adsorption above the site T 1 , all tested functionals predict that the CO molecule adsorbs on top of the Pd atom, with the molecule perpendicular to the surface and the C atom pointing towards the surface. The LDA functional predicts a much stronger CO-surface interaction than the GGA and hybrid functionals. Furthermore, while the C-Pd 1 distance d(C-Pd 1 ) calculated with the different GGA functionals is closer to the HSE06 value than the one obtained with the LDA functional, the latter predicts a C-O bond length d C-O in better agreement with the HSE06 result. The obtained differences are more pronounced in the case of CO adsorption above H 1 . Although for all tested functionals the CO molecule relaxed towards the center of the three Ga atoms of the Ga 3 layer (CO perpendicular to the surface, with the C atom pointing towards the surface), the C-Ga 3 distance changes significantly as function of the exchange-correlation functional. Here, the obtained LDA and PBE final geometries indicate that the CO molecule adsorbs on the site H 1 (d(C-Ga 3 ) LDA = 1.55 Å, d(C-Ga 3 ) PBE = 1.68 Å), i.e. CO does interact with the surface at the site H 1 . However, the adsorption energies E ads calculated with these two functionals present a difference of 0.55 eV. On the other hand, both the final geometry and the adsorption energy obtained with the revPBE and RPBE functionals indicate that the CO molecule desorbs (in both cases E ads E 0 eV and d(C-Ga 3 ) = 3.49 Å), i.e. does not interact with the surface at the site H 1 . In line with the ''CO-adsorption puzzle'', 11,14,15 the obtained results show significant differences between the LDA, GGA and hybrid functionals. Even though the difference in the adsorption energies for the sites T 1 and H 1 is similar for all tested functionals (E1 eV), total adsorption energies change significantly with the choice of the functional.

(b) Vibrational frequencies
The calculated C-O stretching frequency is 2247 cm À1 for the free molecule and 2180 cm À1 for the CO molecule adsorbed at The calculated red shift with the RPBE functional is 148 cm À1 , much larger than the experimental value (65 cm À1 ). Furthermore, the calculated C-O stretching vibration frequency for CO above H 1 is 2241 cm À1 , applying the HSE06 functional. Thus, the vibrational frequency difference between the C-O stretching mode of the free molecule and the CO molecule above H 1 is only 6 cm À1 , consistent with the nearly-unbounded behavior of the CO molecule above the site H 1 predicted by the HSE06 functional.
(c) Electronic structure of the adsorbate-substrate system The Blyholder model is generally applied to describe CO adsorption on surfaces containing transition metal (TM) elements. [71][72][73][74][75][76][77][78][79][80][81] According to this model, there is (i) a CO-5s interaction with the empty TM states of the surface (molecule-to-surface transfer of electrons), that is compensated by (ii) a p-bonding interaction of the filled TM states of the surface with the CO-2p* antibonding states (surface-to-molecule back-donation). The analysis and discussion of the electronic structure of the CO molecule adsorbed at the GaPd(% 1% 1% 1) surface calculated with the HSE06 functional, in the light of the Blyholder model, is given below. CO-surface interaction on site T 1 . The adsorbate-substrate bonding is captured by the emerging features of the sDOS, mainly found in the topmost atomic layer Pd 1 , as shown in Fig. 7. Particularly, the sDOS of the clean surface is non-uniformly suppressed between À6 and À1.7 eV, and new structures (peaks), located at À11.5, À8.4 and 3.4 eV, are found upon CO adsorption. These new peaks are also present in the partial sDOS of the adsorbed CO molecule. In agreement with experimental photoemission findings, the peaks at À11.5 and À8.4 eV can be attributed to the interaction of the surface with the 4s and 5s + 1p (overlapped) states of the adsorbed CO molecule, respectively. [72][73][74][75] On the other hand, the peak at 3.4 eV above the Fermi level can be interpreted as an indication of the surface/ CO-2p* interaction. 75,76,78 The integration analysis of the partial contributions to the total sDOS of the Pd 1 topmost layer after CO adsorption shows that the area under the Pd 4d band (between À6 and À1.7 eV) corresponds to 4.25 states per cell, which is 0.21 states per cell less than before adsorption. On the other hand, the areas under the sDOS at energy intervals where the peaks due to the surface/CO-4s and surface/CO-(5s + 1p) interactions occur (that were empty before adsorption) correspond to 0.05 states per cell (composed of 73% Pd d states, 8.8% Pd p states and 18.2% Pd s states) and 0.15 states per cell (composed of 80.5% Pd d states, 15.7% Pd p states and 3.8% Pd s states), respectively. Thus, considering that the development of these two peaks comes at the expense of the observed Pd 4d band depletion, the obtained results suggest that a Pd-C s-bond is created upon CO adsorption, in agreement with the Blyholder mechanism. However, differences with this model are also found. Besides the participation of Pd s and p states to the Pd-C s-bond as shown above, the C-O bond length d C-O is almost not affected after the molecule adsorbs on T 1 (differences smaller than 0.01 Å, see Table 4), in line with the small red shift of the CO stretching vibrational frequency upon adsorption (67 cm À1 ). Thus, the C-O bond is only weakly affected by the adsorption, revealing a rather different (weaker) s-bonding interaction than the expected from the Blyholder model, as corroborated by thermal desorption and vibrational spectroscopy. 4,7 Additional confirmation of this behavior is obtained when analyzing the possible Blyholder back-donation process. Upon adsorption, the peak assigned to the surface/CO-2p* interaction is mostly dominated by Pd d states (composed of 78.8% Pd d states, 13.9% Pd p states, and 7.3% Pd s states). Previously, the sDOS of the clean surface was dominated by Pd s states (56.3% of the total area), with the total area below the sDOS in this region (between 2.7 and 3.9 eV) being almost the same before and after CO adsorption (0.25 states per cell). However, this somewhat broad peak is shifted well above the Fermi level, i.e., it is empty. This means that in this case the weakening of the C-O bond [11][12][13][14][15]71 does not occur through the backdonation process, because the latter does not take place. The small red shift of the CO stretching vibrational frequency upon adsorption must then be explained through another mechanism, such as, e.g., charge attraction/repulsion between the positively charged C atom and the negatively charged Pd and O atoms.
CO-surface interaction on site H 1 . Fig. 8 presents the layerresolved sDOS of the Pd 1 model before and after CO adsorption on H 1 . Contrary to the observations of CO adsorption on T 1 , in this case the bonding interaction between the molecule and the topmost atomic layers is almost non-existent, confirming the free-molecule-like behaviour of CO above H 1 . Here, only one peak around À6.4 eV makes the difference between the sDOS calculated before (clean surface) and after CO adsorption above H 1 .
The partial CO-contribution to the total sDOS of the adsorbatesubstrate system clearly shows that its electronic states exhibit the structure as in the gas phase with the 4s, 1p, 5s and 2p* states at À11.7, À8.9, À6.4 and 3.0 eV, respectively.

(d) Chemical bonding analysis of the adsorbate-substrate system
Topological analysis of the electron density according to Bader's quantum theory of atoms in molecules (QTAIM) 82 is used in combination with the electron localizability indicator (ELI-D) [44][45][46] to characterize the CO-GaPd(% 1% 1% 1) system. The ELI-D is an efficient quantum chemical tool for the analysis of chemical bonding in real space, particularly, in the understanding of the properties of intermetallic compounds and the modeling of the adsorption of atoms and molecules where both covalent and ionic interactions may occur. 83 The three known ELI-D basins for the free CO molecule in the valence region 83 are retrieved in the surface-adsorbate system: two representing lone pairs -one located at the oxygen atom (4.0 electrons in the free molecule) and one located at the carbon atom (2.5 electrons) -and one basin accounting for the CRO triple bond with an electron count of 3.3 electrons. Upon adsorption of CO on the site T 1 , 0.2 electrons are transferred from the CRO bonding basin to the oxygen lone pair and the lone-pair at the carbon atom turns into a two-center C-Pd bond with an electron count of 2.8. Here, carbon contributes 2.4 electrons, while 0.4 electrons are coming from the palladium atom. In agreement with the analysis of the electronic structure given above, the formation of a C-Pd s-bond is also confirmed by the ELI analysis, as proposed in the Blyholder model. In the same line, the back donation of electrons from the metal into the empty C-O antibonding 2p* states does not take place (see Fig. 7). Instead, the negatively charged palladium atom leads to a shift of the electron density to the oxygen lone pair, thus weakening the C-O bond. In this context, the clean GaPd(% 1% 1% 1) surface resembles a charged egg carton, where the ''dimples'' correspond to the unshielded negatively charged Pd centers, while the positively charged Ga species -shielded by partially dangling bonds -form the ''spacers''. Due to the very stable, closed-shell nature of carbon monoxide, strong Pauli repulsion occurs when the CO molecule approaches a dangling bond. Since the dangling bonds shield most of the surface, only two potential adsorption sites out of 15 crystallographically possible (see Fig. 6) are left for the CO molecule: on-top of the Pd 1 atom (T 1 site) and the hollow site at the center of the triangle formed by the Ga 3 atoms (H 1 site). In addition, the orientation of the CO molecule (carbon oriented to palladium) fits the expectations from the charged egg carton model, resulting from the attractive Coulomb interactions of the positively charged C atom (+1.2) and the negative charge on the palladium (À0.4) according to the QTAIM analysis.
In agreement with the much smaller calculated adsorption energy for the H 1 site, no significant chemical interactions between the substrate and the CO molecule are found by ELI-D analysis upon CO adsorption.

VI. Conclusions
The CO adsorption on the (% 1% 1% 1) surface of the intermetallic compound GaPd has been examined using a quantum-chemical all-electron full-potential approach. Comparison of the GaPd bulk, clean surface, and CO adsorption obtained with LDA-, GGA-, and hybrid Hartree-Fock/GGA-based functionals was considered. The choice of the functional was found to have a decisive influence on the bulk and surface electronic structure of GaPd. In particular, bulk calculations using the different functionals lead to differences in the optimized lattice parameter resulting in different bulk electronic structure. Despite the fact that all the here tested methods find the Pd on-top site as the most preferred one for CO adsorption on the GaPd(% 1% 1% 1) surface when compared to the H 1 site, the obtained results demonstrate that classic LDA and GGA functionals cannot retrieve a correct description of the energies of the substrate-adsorbate electronic states, leading to the wrong description of the CO adsorption on a surface containing Pd atoms, including important differences in adsorption energies and surface-adsorbate interaction as reflected in the calculated vibrational frequencies. On the other hand, experimental findings and the results obtained with the hybrid functional HSE06 are in excellent agreement. This holds for the predicted adsorption sites, adsorption energies, as well as vibrational frequencies. This study thus reveals the computational level necessary to derive reliable predictions from DFT calculations regarding the chemical properties of intermetallic surfaces.