Interface dipoles of organic molecules on Ag(111) in hybrid density-functional theory

We investigate the molecular acceptors 3,4,9,10-perylene-tetracarboxylic acid dianhydride (PTCDA), 2,3,5,6-tetra uoro-7,7,8,8-tetracyanoquinodimethane (F4TCNQ), and 4,5,9,10-pyrenetetraone (PYTON) on Ag(111) using densityfunctional theory. For two groups of the HSE(\alpha, \omega) family of exchange-correlation functionals (\omega = 0 and \omega = 0.2\AA) we study the isolated components as well as the combined systems as a function of the amount of exact-exchange (\alpha). We find that hybrid functionals favour electron transfer to the adsorbate. Comparing to experimental work-function data, we report for (\alpha) ca. 0.25 a notable but small improvement over (semi)local functionals for the interface dipole. Although Kohn-Sham eigenvalues are only approximate representations of ionization energies, incidentally, at this value also the density of states agrees well with the photoelectron spectra. However, increasing (\alpha) to values for which the energy of the lowest unoccupied molecular orbital matches the experimental electron affinity in the gas phase worsens both the interface dipole and the density of states. Our results imply that semi-local DFT calculations may often be adequate for conjugated organic molecules on metal surfaces and that the much more computationally demanding hybrid functionals yield only small improvements.


Introduction
A reliable theoretical description of the interaction between inorganic substrates and conjugated organic molecules is of fundamental importance for a variety of fields, including organic photovoltaics and (opto)electronics. For instance, the adsorption of strong electron acceptors creates an interface dipole modifying the work function of the inorganic substrate. Thus, a layer of acceptors can be inserted at a metal/organic interface to tune the charge-injection barriers into the organic material [1][2][3]. Here we focus on these so-called charge injection layers and investigate by means of density-functional theory (DFT) how molecular adsorbates affect the interface dipole formation.
In the framework of effective single-particle theories such as Kohn-Sham (KS) DFT, electron transfer from the substrate to the organic acceptor occurs through (partially) filling the lowest unoccupied molecular orbital (LUMO), which becomes pinned at the substrate's Fermi energy. The dipole that builds up as a result of the charge transfer increases the work function. However, most common exchange-correlation functionals, in particular (semi) local functionals, suffer from noticeable electron self-interaction errors and the absence of the derivative discontinuity in the exchange-correlation potential. Subsequently, they place unoccupied orbitals too low and occupied orbitals too high in energy, which could result in a spurious charge transfer. This deficiency can be largely reduced by employing hybrid functionals that include a fraction α of exact exchange. However, many hybrid functionals worsen the bulk properties of metals [4], which could adversely affect the description of the adsorption process. We address this conundrum and investigate how the transition from (semi) local to hybrid functionals affects metal/organic interfaces for two groups of the Heyd-Scuseria-Ernzerhof (HSE) [5] family. HSE builds on the Perdew-Burke-Ernzerhof (PBE) [6] generalized gradient functional and the corresponding hybrid functional PBEh [6][7][8].
In PBEh a fraction α of PBE exchange is replaced by exact (Hartree-Fock) exchange. A common choice is α = 0.25, more commonly known as PBE0 [6][7][8]. In HSE an additional range-separation parameter ω is introduced that limits the exact exchange to the short range. ω = 0 then corresponds to PBEh. Another common choice is the HSE06 parameterization that corresponds to α = 0. 25 and ω = 0.2 Å −1 [5,9]. In HSE06 exact-exchange is thus active in a range of ≈5 Å, which is of the order of typical adsorption distances of molecules on Ag (111).
We have recently shown for pyridine on the non-polar ZnO(1010) surface that the final work function depends strongly on the mixing parameter α, while the interface dipole and therefore the adsorption-induced work function change do not [10]. However, there the interface dipole originates almost exclusively from the formation of a covalent bond. We here extend our previous study to charge transfer systems. As test cases, we study perylene-tetracarboxylic acid dianhydride (PTCDA), 2,3,5,6-tetrafluoro-7,7,8,8-tetracyanoquinodimethane (F4TCNQ) and 4,5,9,10-pyrenetetraone (PYTON) (structures and full names are given in figure 3) adsorbed on Ag (111), since for all three cases, experimental unit cells and work function (changes) as well as (semi) local DFT results are available [11][12][13][14]. In addition, the silver lattice constant differs by less than 0.2% between PBE, HSE06 and PBE0 [15] and therefore we can keep the lattice constant fixed for all studies reported below.

Self-interaction and hybrid density functionals
Organic/inorganic interfaces typically contain several 100 atoms per unit cell. The firstprinciples' method of choice for systems of this size is KS DFT. (Semi) local DFT functionals are particularly popular, due to their computational efficiency. In a few notable cases manybody perturbation theory [16][17][18][19][20] or quantum Monte Carlo techniques have been used [21], but due to their computational cost they still remain the exception. Furthermore, any perturbative treatment relies on the assumption that the electron density at the interface is well described by the zeroth-order calculation, which is usually KS DFT. Thus any erroneous electron transfer that might occur in KS DFT is difficult to rectify in such a perturbative treatment [22].
The main deficiency of (semi) local DFT functionals for interface calculations (apart from the missing image effects in the KS energies [16,17,23]) is the self-interaction error [24,25]. This interaction of an electron with itself leads to a delocalization of electron density and a total energy that is no longer piecewise linear as a function of the electron number [26][27][28][29][30]. Concomitantly, the error is not only known as delocalization error [27][28][29] but is also referred to as many-electron self-interaction error [24,25,31] or non-Koopmans compliant error [32]. We illustrate this behaviour for PTCDA and the PBE functional in figure 1.
A priori the deviation from straight-line behaviour says little about the ionization energies of a molecule in the gas phase. The ionization potential (IP) and the electron affinity (EA) are given by the total energy difference between the neutral system and the anion or the cation. This is known as the -self-consistent field ( SCF) approach and performs well for atoms and small molecules also for local and semi-local functionals [33][34][35]. Thus, total energy differences are hardly affected by the self-interaction error. It should be noted, however, that vibronic effects are typically not included in SCF and will also not be included in our study.
In KS DFT, the orbital energies are given by the derivative of the total energy with respect to the particle number [36]. The PBE total energy of PTCDA in figure 1 exhibits concave parabolic behaviour. As a result, the corresponding KS eigenstate energy changes linearly with occupation. For the neutral system (N = 0 in figure 1), the PBE KS eigenvalue then considerably underestimates the experimental value. This underestimation can have profound consequences when we bring two different subsystems into contact. In the limit of negligible chemical interaction, this can give rise to an overestimation of charge transfer [37] and underestimation of charge-transfer excitation energies [38]. In the worst case, spurious charge transfer results even at infinite separation [31], if the highest occupied KS state of one subsystem lies above the lowest unoccupied KS state of the other [22,31,39].
Hybrid functionals are a popular approach to mitigate the self-interaction error, as we will show in more detail in section 4.1. In this work, we study the impact of exact exchange for hybrid functionals of the PBE hybrid (PBEh) form (1) E xc is the exchange-correlation energy, E exact x the Hartree-Fock (i.e. exact) exchange energy, α the mixing parameter, E PBE x the PBE exchange and E PBE c the PBE correlation energy. Functionals in the Heyd-Scuseria-Ernzerhof (HSE) family [5] additionally depend on a rangeseparation parameter ω that spatially splits the exchange energy into a long-range (LR) and a short-range (SR) part At short range, the exchange energy is given by PBEh and at long range by PBE, while the correlation energy is always given by PBE. The HSE family of functionals facilitates a smooth transition from the PBE [6] generalized gradient approximation (α = 0, ω = arbitrary) to the PBE0 [6][7][8] hybrid functional (α = 0.25, ω = 0) via the popular HSE06 [9] functional (α = 0.25, ω = 0.2 Å −1 ). The process of limiting exact exchange to the short range can also be viewed as static screening of exact exchange. ω can then be interpreted as the inverse of a screening length. A value of ω = 0.2 Å thus corresponds to a screening length of ≈ 5 Å, which is slightly larger than the typical adsorption distance of molecules on Ag (111) and thus sufficient for our purposes. In principle, one would like to screen the exact exchange dynamically, preferably within the framework of the random phase approximation [40][41][42]. Unfortunately, at present such calculations are too expensive for large periodic systems. Limiting the range of exact exchange significantly reduces the computational time. However, for any non-zero value of ω the potential of any HSE functional will inherit the incorrect exponential asymptotic decay of PBE. In this work, we will focus on the impact of α, but to clarify the impact of range separation, we use two different groups of the HSE functional family, ω = 0 and 0.2 Å −1 . While it was shown for several solids, including Ag, that HSE06 yields better structural properties and atomization energies than PBE (although the improvement over PBE0 is small) [15,43], a survey of the work function of six different transition metals showed no systematic improvement [43]. Also for several molecular properties, such as heterolytic dissociation of dimers, Rydberg excitations or vibrational properties, HSE06 yields results of the same quality as PBEh [44]. It is not a priori clear if this is also the case for systems with strong charge-transfer character. Addressing this question will be the topic of this paper.
Applying the concept of straight-line behaviour introduced earlier in this section, we show in figure 2 the variation of the LUMO energy ( ) of PTCDA in the gas phase with respect to its orbital occupation. A many-electron self-interaction free description is reached when the orbital energy does not depend on its occupation, i.e. when it forms a straight, horizontal line. For HSE(α), this criterion is never reached. Figure 2 shows a strong positive slope of for all α, corresponding to a convex curvature of the energy versus occupation curve. For PBEh(α), on the other hand, an almost perfect horizontal line is observed for α = 0.7. Smaller values of α show positive, larger values a negative slope. The fact that all lines cross close to n = 0.5 shows that the Slater-Janak transition state relation is fulfilled [36,45]. The results also show that the Slater-Janak transition state is practically independent of the choice of the functional.
For PTCDA (and later also for F4TCNQ and PYTON) we therefore conclude that PBEh(α) can be made self-interaction free and that this requires a large value of α = 0.7. In HSE(α), on the other hand, the self-interaction error can never be fully removed (figure 2). We attribute this to the fact that in PBEh(α), exact exchange is not range limited [6], which implies that the potential is closer to the exact 1/r behaviour, where r is the distance from the molecule [6][7][8]. In HSE(α), on the other hand, exact exchange is short range [5,9] and the potential asymptotically follows the incorrect exponential decay of the PBE functional [46]. It should be noted that the correct 1/r behaviour could be restored in hybrid functionals by inverting the range-separation. Long-range hybrid functionals [46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63][64] apply Hartree-Fock in the long range (rather than PBE in HSE) and a semi-local functional in the short range. However, with a few notable exceptions (e.g. [65,66]), they have so far not been applied to bulk materials and are beyond the scope of the present work.
We will, however, not dismiss the HSE(α) family at this stage, because we have not considered the metallic substrate, yet. For our surface/adsorbate system we therefore pose two questions. Firstly, which α is best for the combined system? And secondly, how do different hybrid functionals perform for the same α? In the present study, we attempt to answer these questions by systematically investigating the effect of α for three different molecules adsorbed on Ag (111).

Computational details
All calculations were performed with the Fritz Haber Institute ab initio molecular simulations (FHI-aims) code [83]. Surfaces and interfaces were modelled by periodic slabs containing five layers of the metal. For surfaces, a region of 30 Å vacuum was inserted between the metal slab and its periodic replica. The height of the unit cell was kept constant when adding the molecule, which causes a reduction of the vacuum layer by less than 5 Å. Polarization through the vacuum was prevented by means of the dipole correction [84]. Unless otherwise noted, all KS energies are reported with respect to the vacuum level (VL). For non-periodic systems, this pertains to the potential energy of an electron at infinite distance. In periodic systems, the VL is given by the electrostatic potential of an electron far above the surface. For practical purposes, here the energy 2 Å beneath the place of the dipole correction is taken. All computational parameters discussed below were converged for α = 0.8 to a threshold of 10 meV for the work function. A 'Tier 2' basis was used for C, N, O and F, which consist of the minimal basis plus one set of basis functions up to an angular momentum of 3 (d-functions) and one set of basis functions up to an angular momentum of 5 (g-functions). The 'Tier 2' basis for H consists of two sets of basis functions up to p and d functions, respectively. For Ag, only a single set of basis functions beyond the minimal basis, up to f -functions, was included. We verified that the removal of the  (111). Only the top Ag-layer (white) is shown for clarity. Grey balls denote carbon atoms, cyan, hydrogen atoms, blue, fluorine atoms, red, oxygen atoms and green balls represent nitrogen atoms.
g-function, which is part of the 'tight' defaults, affects the total work function and the density of states by less than 10 meV, while significantly decreasing the computational demand and the memory requirements of the calculations. To obtain accurate surface dipoles, the cutoff potential of all basis functions was increased from 4 to 6 Å, which causes an increase of the work function of the pristine surface by approximately 50 meV. For integrations, a tightly converged Lebedev grid was used. A 35×35×1 off-k-point grid has been used for the primitive surface unit cell and scaled appropriately for the larger supercells. To account for van der Waals interactions, we employ the vdW surf scheme [85] for the metal-molecule interaction, which is appropriate for metal surfaces and has been shown to yield accurate adsorption distances and energies. Within this method, the van der Waals parameters for Ag are re-evaluated based on the dielectric function of the metal. Throughout this work, the values reported by Ruiz et al [85] were used. The SCF-cycle was converged to a threshold of 10 −4 eV for the total energy and 10 −3 eV for the sum of eigenvalues. The unit cells were deduced from scanning tunnelling microscopy experiments [11,12,14], and are shown in figure 4. All geometries were optimized using the PBE+vdW surf functional until the remaining forces were smaller than 10 −3 eV Å −1 .

Isolated gas phase molecules
As a first step, it is instructive to study the effects of exact exchange for the isolated subsystems. In principle, it would be best to compare DFT calculations to advanced computational methods, preferably coupled cluster singles double with perturbative triples (CCSD(T)), which is considered the 'gold standard' of chemistry. Unfortunately, CCSD(T) calculations with a converged basis set are currently not available in the literature for charged, spin-polarized molecules of this size. For experimental reference values, on the other hand, we were only able to find gas phase electron affinities for F4TCNQ. Horke et al [86] measured the photoelectron spectrum of F4TCNQ anions at a photon energy resonant with the D 2 ← D 0 transition. By comparing with the spectrum of the unfluorinated derivative TCNQ, they estimated an electron detachment energy of 3.2 eV. However, this analysis is aggravated by the fact that the D 2 state is below the VL for F4TCNQ, but above the VL for TCNQ. Experimental data for gas phase ionization are only available for PTCDA, for which Dori et al [87] determined an IP of 8.2 eV using photoelectron spectroscopy. In neither experiment has extrapolation to T = 0 K or removal of zero-point vibration energies been attempted. It should be mentioned that the experimental ionization threshold may be affected by electron-vibrational coupling or the vibrational fine-structure. These effects are not accounted for in our calculations. Moreover, it is not always clear whether the ionization process in the experiment occurs vertically, i.e. at fixed geometry, or adiabatically. For these reasons, perfect agreement between DFT calculations and experiment must not be expected. For PTCDA and PYTON, electron affinities in the gas phase have not yet been measured to the best of our knowledge.
We first consider the EA in SCF, which is given by the ground-state energies E tot of the charged (N + 1 electrons) and uncharged systems (N electrons): We only consider vertical electron affinities (i.e. the atomic positions were kept fixed in their equilibrium geometry of the neutral molecule) and compare the SCF-EA with the unoccupied LUMO (n = 0) eigenvalues of the neutral molecule, the half (n = 1 2 ) and the completely filled orbital (n = 1) as a function of α for both hybrid functionals. The completely filled orbital corresponds to the highest occupied molecular orbital (HOMO) of the N +1 electron system. In chemistry textbooks, it is commonly referred to as singly occupied molecular orbital (SOMO) [88]. We will adopt this notation here in order to avoid confusion between the HOMO of the N and the N +1 electron system and to clearly distinguish the spin-polarized, partially filled orbital in the gas phase from the partially filled, spin-unpolarized orbital after adsorption (see below). Figure 5 shows that the EA in SCF lies almost exactly half way between the n = 0 and 1 energy and agrees with the n = 1 2 energy. Thus the Slater-Janak transition state relation [36,45] is fulfilled, as expected. We find that upon increasing α, the energy of the unoccupied LUMO strongly increases with respect to the VL, while the energy of the fully occupied orbital decreases. For the fully occupied or fully empty orbitals, the energies are significantly different when the same α but different ω is used. However, PBEh(α) and HSE(α) give virtually identical results for the half-filled orbital and SCF, in agreement with [89], which reflects that the selfinteraction error hardly affects these two approaches.
For atoms and small molecules, SCF typically performs well for ionization potentials and electron affinities already at the PBE level [33][34][35]. However, for F4TCNQ, we find that SCF gives electron affinities that are more than 1 eV lower than in experiment. Also G 0 W 0 calculations [90], that have become the method of choice for electron addition and electron removal energies in solids [91,92], do not provide agreement with experiment, as shown in the appendix. For the closely related molecule TCNQ, CCSD(T) calculations with a small (augcc-pVDZ) basis set were reported [93]. Even they cannot reproduce the experimental results, although a significant improvement over SCF values was reported [93]. The origin of this discrepancy is not yet understood.
Returning to our definition of the many-electron, self-interaction error in section 2, we observe that in our α dependent plots this condition is fulfilled when the energies for n = 0, 1 and SCF-EA (n ∼ 1 2 ) cross at the same point. For our three molecules, we find that for HSE(α) no such crossing point exists. Conversely, for the PBEh(α) functional the lines cross at α ≈ 0.7. In other words, for PBEh(α) a self-interaction-free description of the orbital can be obtained, but not for HSE(α). This is in agreement with our observations in section 2.  [86]. Note that SCF and n = 1 2 are virtually on top of each other. Figure 6 shows a plot equivalent to figure 5 for the ionization potential (IP) of the molecules. In analogy to EA, the -SCF-IP was calculated as Here we find SCF and experiment to agree within 0.25 eV for any α and both ω. Comparing KS eigenvalues with SCF, we encounter a similar situation as for the EA. For HSE(α), the HOMO and the negative of the IP almost never agree, while for PBEh(α), agreement is found close to α ≈ 0.6 throughout. This value of α is somewhat smaller than for the electron attachment energy, i.e. HOMO and LUMO are thus affected differently by selfinteraction. An identical picture is obtained when comparing to perturbative G 0 W 0 energies, see appendix. For the organic molecules in the gas phase, we can now conclude that (a) the KS eigenvalues of the PBEh(α) group are better than those of the HSE(α) group and (b) generally large values of α are needed.

The pristine metal surface
Although it is, in principle, possible to perform SCF calculations for an embedded metal cluster without periodic boundary conditions, it is often challenging to ascertain whether a discrepancy between theory and experiment originates from deficiencies of the DFT functional or from finite size effects. For periodic systems, SCF calculations are aggravated by the fact that the addition of extra charge in the unit cell requires the introduction of a compensating charge to prevent the divergence of the electrostatic energy. The energy contribution from the compensating background then needs to be removed carefully. For periodic surfaces, any homogeneous compensating background will introduce a dipole and thus an additional divergence with respect to an increase in the vacuum separation. Again, correction schemes exist [94][95][96][97], but the physical energy contributions need to be carefully disentangled from those of the artificial, compensating background. We therefore refrain from SCF calculations for periodic surface supercells and determine α differently.
For semiconducting substrates, it has been suggested to make α dependent on the static dielectric constant [98,99], but that is not well defined for metal surfaces. Some of us have recently established a correlation between defect formation energies and the valence band width as a measure of the cohesive energy [100]. This can be used to determine α for semiconductors or insulators, provided accurate reference data are available for the valence bandwidth [100]. For metals, no such relation has been established, yet. Instead, we decided to focus on the work function and the density of states from periodic boundary calculations.
For Ag (111), the evolution of the Fermi level with respect to the VL (i.e. the work function ) is shown in figure 7(a) for a varying fraction of exact exchange. Both HSE(α) and PBEh(α) yield the highest work function for α = 0, i.e. in the PBE limit. Increasing α results and PBEh(α) (circles), compared to experimental results [12-14, 101, 102]. Right: Ag d-band onset as a function of α for HSE(α) (boxes) and PBEh(α) (circles), compared to experimental results [103].
in a significant decrease of the work function of up to 0.4 eV (HSE(α)) and 0.8 eV (PBEh(α)). Stroppa and Kresse [43] reported a similar decrease in work function for several different metals when going from PBE to HSE06, and ascribed it to a change in the surface dipole that results from an electron redistribution in the d bands. Once again, ideally we would like to compare our results with higher-level theory. Unfortunately, coupled cluster calculations for extended surfaces are even less tractable than for molecules in the gas phase [104]. The comparison with the experimental work function is aggravated by the fact that a wide range of values has been reported in the literature. These range from = 4.46 eV [105], and 4.5 eV (Ag (111) prior to F4TCNQ and PYTON deposition [12,14]) to 4.8/4.9 eV (Ag(111) prior to PTCDA deposition [101,102]) as figure 7(a) shows. Despite this considerable spread, all experimental values are larger than the PBE work function. Such an underestimation of the work function is found for many metals [106].
Although there is no requirement that the energetic position of the KS 4d-bands must agree with the corresponding photoemission peak, their incorrect description is often made responsible for artificial hybridizations. It is thus instructive for this work to also briefly study their dependence on α. The onset of the Ag 4d band is shown in the right panel of figure 7. The scale is one order of magnitude larger than for the work function. Between α = 0 and 1, the d band shifts by 5.1 eV (HSE(α)) or 7.0 eV (PBEh(α)). Unlike , the 4d band is too high in energy when no exact exchange is present. This behaviour is well known for late transition metals [107,108]. The calculations coincide with the experimental binding energy of 3.7 eV [103] for relatively low fractions of exact exchange, α = 0.15-0.20 (PBEh(α)/HSE(α)), again in agreement with other studies [15,109]. We can thus conclude that (a) HSE(α) is more appropriate than PBEh(α) for the metal surface and (b) generally small values of α are needed.

Fermi-level pinning for molecules on metallic surfaces
Before discussing the adsorbate systems in detail, it is instructive to briefly revisit the concept of Fermi-level pinning. The basic idea behind Fermi-level pinning is illustrated in figure 8. In the limit of large separation and no interaction between the subsystems, the metallic surface is described by its work function (or equivalently, the Fermi energy E F ), whereas the organic adsorbate is characterized by its IP (HOMO) and EA (LUMO). Both systems share a common VL. Once the two are brought into contact, the molecular orbitals broaden (IP and EA in figure 8). Also, the screening of charges near a metal surface allows the adsorbate to be more easily ionized, i.e. IP and EA are reduced [16,17,23] (to IP and EA in figure 8(b)). This effect is commonly referred to as band-gap renormalization or image effect.
If the LUMO was not already below the Fermi level at infinite distance, it might be now. In that case the LUMO becomes partially occupied as indicated in figure 8(b). The electron transfer induces a shift of the adsorbate orbitals proportional to the dipole density [1], µ/A: Since in thermodynamic equilibrium the whole system must share a common Fermi energy, will align the Fermi energy with the LUMO . In the end, the work function of the combined system, , is given by the difference of the LUMO from the VL, see figure 8(c). In real systems the situation is usually more complex, since hybridization and geometry distortions may lead to orbital broadenings which deviate significantly from the ideal Lorentzian shape and induce fractional occupations which pin the LUMO level some tenths of eV away from the peak maximum, as sketched in figure 8(d). Partially filled states at E F associated with the LUMO have been observed in the photoemission experiments for the three molecules discussed here [12,14,101,102].

Level alignment in the non-interacting limit
The previous section illustrates the importance of having a correct description of both the Fermi energy and the molecular frontier orbitals for Fermi-level pinned systems. At the same time, we have established that within HSE(α) and PBEh(α), the metallic surface and the organic adsorbate require very different values of α. The obvious question is whether, and how, a compromise between those two different regimes can be found for a global value of α. Before we examine the molecules adsorbed on the metal surface, however, it is insightful to examine the level alignment in the non-interacting-subsystem regime, in an analogy to figure 8(a). For comparison we combine the data of the isolated molecules from figure 5 and the clean Ag surface from figure 7 into figure 9. For both functional groups, the empty LUMO crosses the Ag Fermi level for PTCDA and PYTON. For F4TCNQ, only PBEh exhibits a crossing of the LUMO and the Fermi energy. For PTCDA and PYTON, the crossing occurs at small α. Similarily, for HSE(α) the filled SOMO is above the Fermi energy for all three molecules for any α. In contrast, PBEh(α) exhibits a crossing for PTCDA and F4TCNQ at large α. We also find that the SCF EA lies significantly above the Fermi energy for both PTCDA and PYTON for any α. Only F4TCNQ shows EAs which are partly below the Ag E F . These results convey the impression that PTCDA and PYTON on Ag(111) might not be Fermi-level pinned systems at all or that there is at least a transition from a charge-transfer to a non-charge-transfer regime, depending on the choice of α. It should be kept in mind, however, that important factors are missing from this picture. On the one hand, the electron density of the more polarizable component is partly displaced, as soon as two systems come into contact [110]. This effect is known as Pauli pushback or the pillow effect and significantly reduces the surface dipole. For Ag (111), reductions of 0.4-0.6 eV are commonly reported [111,112]. On the other hand, PYTON and F4TCNQ undergo noticeable distortions upon adsorption, leading to the formation of an intrinsic molecular dipole, which can affect the relative position of the molecular and surface levels. Even in the absence of such distortions, the π -electrons and the positive nuclei can build up a considerable local dipole moment that shifts all molecular levels relative to the surface [113]. Finally, the screening of the charge on the molecule by the metal electrons is, of course, completely absent in the calculations of the isolated molecule in the gas phase. Indeed, these effects combined are strong enough so that we find Fermi-level pinning for all three molecules on Ag (111), as discussed in the next section.

Electron transfer to acceptor molecules on metallic surfaces
Despite the aforementioned deficiencies of standard (semi) local functionals, remarkable agreement between the computed and the experimental work function has been reported for the adsorption of electron donors or acceptors on metal surfaces [2,71,[114][115][116][117][118][119][120][121]. Previously, this was attributed to a fortuitous (partial) cancellation of errors [115], in particular due to the fact that the HOMO-LUMO underestimation of (semi) local functionals is of the same order of magnitude as the band-gap renormalization of the molecular states due to image effects that are absent from the PBE eigenvalues. Since hybrid functionals correct the former, but not the latter [18], the question arises whether the admixture of exact exchange improves or deteriorates the agreement between theory and experiment for the combined system. We focus first on the adsorption-induced work function changes, , since this is a physical observable depending only on the electron density. Figure 10 shows as a function of α. 3 For all three molecules, we observe a systematic, pronounced increase of the work function change with increasing α. Using the two PTCDA/Ag (111) experiments as measures, we estimate the experimental reproducibility of work functions and work function changes to be of the order of at least ±0.2 eV. Comparing the experimental and the calculated work function changes for PTCDA and F4TCNQ therefore yields a good agreement for small and medium values of α (up to approximately 0.5). For PYTON, the work function change (i.e. the interface dipole ) is more sensitive to α than for the other two molecules, and agreement with the experiment is restricted to small values (approximately 0.25). For the net work function of For PYTON, in contrast, the variation in is much larger than in the other two cases, causing a strong net increase of . In general, especially for small values of α, we find that there is little difference between PBEh(α) and HSE(α). This can be attributed to the observation that the charge rearrangements upon adsorption are mostly localized to the region between the molecule and the first metal layer [12,115] and are thus mostly within the short-range region of HSE. For the sake of brevity and clarity, we will show only the results for the HSE(α) group hereafter.
Since the geometry is kept fixed in our calculations, the change in must be due to a change in hybridization or due to increased electron transfer. To gain more insight, we turn to the density of states after adsorption, which is shown in figure 11 projected onto the respective molecule. Although a KS density of states is only approximately a physical observable, Körzdörfer et al [122] pointed out that it is closely related to the experimental photoelectron spectrum, provided that electron-phonon coupling and disorder can be neglected 4 . If the aforementioned increase of the interface dipole is indicative of a larger electron transfer between metal and molecules, the molecular density of states should be shifted further below the Fermi level. For all three molecules, we find the LUMO to be below the Fermi energy and thus partly filled after adsorption. For α = 0, the LUMO resides too high compared to the experimental results. Increasing the fraction of exact exchange shifts the whole density of states to lower energies. At α ≈ 0.25, the LUMO agrees well with the photoemission experiments. A further increase of α shifts the LUMO to even lower energies. The LUMO therefore behaves qualitatively like the SOMO, rather than the LUMO of the isolated molecule (cf figure 9), except for the fact that it agrees with experiment at much smaller fractions of exact exchange. The observed downward shift causes an increased filling of the LUMO (see below for more details), and thus an increased electron transfer. We attribute the reason for this behaviour to the fact that for all molecules, the LUMO at α = 0 is occupied with more than one electron (see below for details), and the exact exchange shifts orbitals down in energy if they are more than half filled and up if otherwise (cf figure 2). We, therefore, speculate that the general trend of increased electron transfer only holds for a system for which the LUMO is more than 50% filled at the PBE level, and that the reverse (i.e. reduced electron transfer) should be observed for systems with less than half the filling at the PBE level.
The α-induced shift is not the same for every orbital. In fact, the HOMO shifts more than the LUMO , causing an increase in the HOMO -LUMO splitting. Interestingly, Figure 12. Generalized KS gap between LUMO and HOMO on the surface (blue triangles), the gas phase molecule in its neutral state (black square), the singly charged anion (red circles) and the di-anion with a doubly occupied LUMO (orange diamonds). however, figure 11 reveals that both orbitals are simultaneously in excellent agreement with experiment at α ≈ 0.25. At this point, we briefly digress to emphasize another interesting finding. Figure 12 shows the evolution of the gas-phase HOMO-LUMO gap of the neutral molecule, the energy difference between the SOMO and the next lowest orbital in the same spin channel (former HOMO), the energy difference between the doubly occupied LUMO and HOMO of the molecular di-anion and the HOMO -LUMO gap deduced from the adsorbed molecules on the surface (E gap ). To exclude geometry effects, the orbital energies were calculated using the geometry that the molecules assume on the surface after adsorption. The calculation for the di-anion was performed spin restricted, i.e. for all molecules we report the singlet rather than the triplet state. For the surface calculations, we determine the gap as peak-to-peak rather than onset-to-onset, since we used an artificial broadening (see section 3). Furthermore, the experimental broadening-and thus the onset-to-onset gap-is partly determined by the electron-phonon coupling of the states and by disorder, which is not included in our calculations [124][125][126][127].
On the surface, E gap increases much slower than the gap between the HOMO and the LUMO in the gas phase. E gap also does not follow the behaviour of the gas phase di-anion despite the fact that at the surface the LUMO is almost completely filled with two electrons (see next paragraph). In contrast, E gap closely tracks the α-dependence of the radical anion in the gas phase. This result is rather unexpected, not only because the occupation of the LUMO on the surface differs significantly from 1, but also because the gas phase anion is an open shell system and was treated spin polarized, while the adsorbate calculation was not spin polarized. At present, it is not clear whether this result is coincidental or not, as correlation does not necessarily imply causation. Nonetheless, this ambiguity might have contributed to the controversy regarding the contentious issue of integer (a single electron to a fraction of molecules) versus partial charge transfer (fractional charge to all molecules), that is particularly well documented for F4TCNQ [128]. Figure 13. Absolute values of the total charge transfer (green triangles), LUMO occupation (red circles) and charge donation (black squares) for the three test molecules on Ag (111). For PYTON the LUMO is doubly degenerate in the gas phase and we therefore show both states (red circles and blue triangles).
To better understand the evolution of and in terms of the properties of the isolated subsystems, it is important to review the charge-transfer mechanism in more detail. To this aim, we project the density of states onto the molecular orbitals of the hypothetical, free-standing monolayer [129]. By integrating the molecular density of states of each orbital from −∞ to E F we then obtain formal occupations for each orbital. This allows us to quantify donation and back donation separately. Here, we define donation as the charge accumulated by the LUMO, and back donation as the total charge of the molecule minus the donation. Although this definition differs somewhat from the more intuitive definition in which donation is given by the sum over the occupation numbers of all occupied orbitals and the back donation as the sum over all unoccupied levels, it has the advantage of being more easily identifiable with the experimental photoelectron spectra, for which the LUMO occupation is evident. To be consistent with earlier works [12], we employed a Mulliken-like projection. The results are therefore subject to the corresponding shortcomings, such as occupations potentially exceeding the 0-2 range and an increasing ambiguity with diffuse basis sets. However, it had been found earlier that the general trends obtained by this analysis compare well with trends obtained, e.g. by real-space integration of the electron density [130]. Our results for the net charge transfer, as well as donation and back donation are shown in figure 13 (note that the absolute values of all charges are plotted).
For PTCDA, we find a considerable net charge-transfer between PTCDA and Ag (111), in agreement with earlier works [113,115,124,131,132]. Using our charge-partition scheme, the charge transfer amounts to approximately 0.4e for α = 0, which increases gradually to approximately 0.7 e for α = 0.8. The origin of the increased charge transfer can be traced back to increased occupation of the LUMO , which rises from 1.4-1.7e. Charge donation, on the other hand, remains almost constant at approximately one electron. The increased charge transfer then gives rise to an increased interface dipole.
F4TCNQ interacts with Ag(111) through the hybridization of the low-lying CN-orbitals with the Ag d-bands [12] and the filling of the LUMO . Already for α = 0, the LUMO is found almost completely below the Fermi energy, giving rise to an occupation of 1.8e.
Exact exchange further stabilizes the now occupied orbital, shifting it to lower energies. However, since the LUMO is already almost completely filled at the PBE-level, this shift only causes a minor increase of the occupation to 1.9e. At the same time, charge donation remains unaffected at a level of 1.3e. The total amount of charge transfer therefore remains essentially constant, increasing only very little as α increases. As a result, the net work function of F4TCNQ on Ag(111) follows the work function evolution of the pristine Ag surface.
In the gas phase, PYTON has a degenerate LUMO. The adsorption-induced distortion and the interaction with the three-fold symmetric surface lift this degeneracy, albeit only by some tenth of an eV. The hybridization-induced broadening of these orbitals is larger than this split and their densities of states continue to overlap. As a result, both orbitals become occupied upon adsorption. Although both orbitals are π-orbitals, they are not equally localized and thus differently affected by self-interaction. This leads to a reordering of these orbitals on the surface when α increases. While the position of one orbital is essentially unaffected, resulting in a constant occupation of approximately 1.2e, the other orbital shifts strongly with α to more negative energies, and thereby increases its occupation from 1.0 to 1.5e. This molecule therefore accepts up to 2.7 electrons in total. Naturally, the large electron donation must be accompanied by a large electron back donation. Indeed, we find a back donation of approximately 1.5 electrons, originating from the interaction of the carbonyl-states with Ag (111). As for PTCDA and F4TCNQ, the donation is only weakly dependent on the choice of α. Once again, the net result is a strongly increasing charge transfer, which agrees with the reported strong increase of the work function with α.
The overarching question that remains, is why the three different molecules react so differently to changes in α, i.e. why the interface dipole increases much more for PYTON than for PTCDA, and for PTCDA more than for F4TCNQ. We could establish no simple relation between the α-dependent or and the EA, the LUMO or the SOMO, neither for the gas phase equilibrium structure nor the structure including adsorption-induced geometry changes. We thus have to conclude that the generalized KS eigenvalues of the gas phase molecules alone are, in general, not reliable indicators for Fermi-level pinning and cannot be used to predict the size of work function changes. Nonetheless, figure 13 clearly shows that the influence of α on charge back donation is negligible. This in turn implies that any change in the charge transfer, which manifests itself in the interface dipole, must be related to the changes in the LUMO position. We can therefore identify two system-dependent variables. First, the position of the LUMO for α = 0. In contrast to the situation on Cu or Au, where complete filling only occurs at low coverages [128,133], the LUMO is almost completely filled for F4TCNQ even at full coverage on Ag (111). Further α-induced downward shifts can therefore not yield a significant change in charge transfer. PTCDA and PYTON, on the other hand, are only partially filled for α = 0.0 and can still accept more charge. Moreover, the change in the transferred charge is directly proportional to the molecular density of states at E F and therefore depends on the area of the LUMO peak that is shifted below E F . As figure 11 shows, the density of states for any given α is always largest for PYTON and smallest for F4TCNQ. The exceptionally large DOS at E F is a consequence of the near-degeneracy of the two LUMOs and explains the extraordinarily large work function (change) of PYTON compared to the other molecules. This observation is in line with the induced density of interface states model [134], which emphasizes the importance of the adsorption-induced molecular density of states for systems of this kind. A second factor is clearly the sensitivity of the orbital eigenvalues to α, which depends on the orbital's localization and its filling. Figure 5 shows that the PYTON orbitals are much more sensitive to exact exchange than those of PTCDA. This gives a larger downward shift of the newly occupied orbital. Both aspects together form a plausible, if not comprehensive, explanation for the observed variation of and the charge-transfer dependence of the three systems.

Summary and conclusions
All the above findings may confer the impression that the adsorption of Fermi-level pinned, conjugated organic molecules on metal surfaces is already well described at the PBE+vdW level, and that improvements obtained from HSE(α) or PBEh(α) are hardly worth the significantly higher computational effort. Indeed, as a general rule of thumb we expect this to be true. Nonetheless, it is clear that some systems react more sensitively to hybrid functional than others, for example when they exhibit a particularly high density of states at the Fermi energy as the case of PYTON demonstrates. We thus expect molecular adsorbates that have degenerate frontier orbitals to require a beyond PBE+vdW treatment. It should also be kept in mind that we only considered molecules that bind to the metal surface through hybridization and donation/backdonation. It is not yet clear what will happen in the absence of hybridization, especially if charge transfer to the organic layers causes spin polarization. Last but not the least, all our calculations place the same charge on all molecules in the layer, because of the periodic boundary conditions. If, however, charge transfer becomes localized on one out of many molecules, corresponding lattice distortions (polarons) may become important as reported for both polymers and small organic molecules [135].
In summary, we have systematically studied the level alignment of three different conjugated organic molecules on Ag(111) as a function of the exact exchange admixture α in the hybrid functionals PBEh(α) and HSE(α). The orbitals of the isolated molecules in the gas phase require a large fraction of exact exchange (α ≈ 0.7-0.8) to become self-interaction free within the PBEh(α) group, while for HSE(α) this condition is never fulfilled. In contrast, the work function of the clean Ag surface is always underestimated by all hybrid functionals, whereas an accurate position of the d-bands with respect to the Fermi energy is obtained for α ≈ 0.25. In general, we find that screened exchange (HSE(α)) works better than PBEh(α) for Ag.
For the combined inorganic/organic systems, we find that HSE(α) shows qualitatively the same trends as PBEh(α), which we attribute to the fact that the charge rearrangements upon adsorption are fairly short range. HSE(α) thus already contains all the necessary physics at a much smaller computational cost than PBEh(α). The interface dipole upon adsorption is best described at α ≈ 0, whereas the net workfunction after adsorption requires α ≈ 0.25. The photoelectron spectra, in particular the HOMO and LUMO positions, are best reproduced for α ≈ 0.25. For the three systems studied here, more exact exchange leads to a systematic increase in charge transfer, which we attribute to the fact that all molecules pin at the Fermi level and exhibit a large LUMO occupation already at the PBE level. This can be directly related to an increase of the LUMO occupation, whereas charge donation remains almost unaffected. However, overall the difference between α = 0 and 0.25 is relatively small. Given the fact that the reproducibility of organic monolayers on metal surfaces is still low and that experimental results exhibit a certain spread, we conclude that PBE-based results for interface dipoles and densities of states are generally adequate for interfaces between coinage metals and conjugated organic molecules. As the example of PYTON shows, tuning α to match orbital eigenvalues and ionization potentials or electron affinities can even be counterproductive for global and short-range range-separated hybrid functionals. eigenvalues never cross the G 0 W 0 lines, whereas the PBEh(α) HOMO crosses at α ≈ 0.7 and the LUMO at α ≈ 0.8. These values are consistent with our observations in section 4.1. In principle, we could use the intersection between G 0 W 0 and PBEh(α) to design an internally consistent G 0 W 0 scheme [22]. The highest KS eigenvalue of a finite system is given exactly in exact DFT [143,144] and its self-energy correction vanishes. This would be the intersection of the DFT-KS and the G 0 W 0 line. We here rely on the fact that this relation is also fulfilled for approximate DFT functionals and the GW self-energy. The intersection then gives us an optimized α value. Since it also determines an internally consistent starting point for G 0 W 0 , the scheme can be viewed as a simplified self-consistent GW procedure [22]. This internally consistent scheme is, of course, also applicable to solids and interfaces. However, currently G 0 W 0 is not yet implemented for periodic boundary conditions in FHI-aims.