Reliable energy level alignment at physisorbed molecule-metal interfaces from density functional theory.

A key quantity for molecule-metal interfaces is the energy level alignment of molecular electronic states with the metallic Fermi level. We develop and apply an efficient theoretical method, based on density functional theory (DFT) that can yield quantitatively accurate energy level alignment information for physisorbed metal-molecule interfaces. The method builds on the "DFT+Σ" approach, grounded in many-body perturbation theory, which introduces an approximate electron self-energy that corrects the level alignment obtained from conventional DFT for missing exchange and correlation effects associated with the gas-phase molecule and substrate polarization. Here, we extend the DFT+Σ approach in two important ways: first, we employ optimally tuned range-separated hybrid functionals to compute the gas-phase term, rather than rely on GW or total energy differences as in prior work; second, we use a nonclassical DFT-determined image-charge plane of the metallic surface to compute the substrate polarization term, rather than the classical DFT-derived image plane used previously. We validate this new approach by a detailed comparison with experimental and theoretical reference data for several prototypical molecule-metal interfaces, where excellent agreement with experiment is achieved: benzene on graphite (0001), and 1,4-benzenediamine, Cu-phthalocyanine, and 3,4,9,10-perylene-tetracarboxylic-dianhydride on Au(111). In particular, we show that the method correctly captures level alignment trends across chemical systems and that it retains its accuracy even for molecules for which conventional DFT suffers from severe self-interaction errors.

Interfaces between molecules and metals are central components in emerging technologies such as organic and molecular electronics and third-generation photovoltaic cells. These interfaces can often strongly influence the electrical and optical characteristics of a specific structure or device.[1,2] A key physical quantity of interest in this context is the energy level alignment between frontier electronic states in the molecule and the Fermi energy of the substrate, EF.[3] It quantifies the energy barrier for electrons or holes crossing between molecule and metal and thus determines the efficiency of charge-carrier injection. As an example, in metal–molecule–metal junctions energy level alignment dictates the current at a certain bias and thus directly influences the charge-transport characteristics of nanoscale molecular devices.[4,5] To understand the inter-relation of structure, chemical composition, and molecule–metal interactions with energy level alignment in full microscopic detail is therefore tremendously important and as such lies at the focus of current scientific efforts in the fields of organic and molecular electronics; see, for example, refs (6−18).

Energy level alignment can be assessed experimentally with high accuracy using photoemission spectroscopy (PES).[19−22] To calculate it reliably from first-principles is still, however, a major challenge for modern electronic-structure methods.[23,24] Many-body perturbation theory, typically using the GW approximation (where G is the Green’s function and W the dynamically screened Coulomb interaction) to Dyson’s equation, in principle allows for a rigorous determination of excitation energies[25,26] and, therefore, energy level alignment. In practice, GW has indeed been successfully applied to a variety of molecules and metals separately. Calculations of the combined molecule–metal interface, which is inherently a complex heterogeneous system, remain, however, scarce within present-day GW. This is because GW calculations are highly demanding computationally, can strongly depend on the underlying starting point,[27,28] and are challenging to converge numerically.[29,30] Density functional theory (DFT)[31] offers an excellent trade-off between accuracy and efficiency and therefore became the workhorse methodology of contemporary molecule–metal interface calculations. However, DFT in its various approximations has been repeatedly shown to not allow for quantitative predictions of energy level alignment.[16,29,32−37] In principle, DFT eigenvalues (other than the energy of the highest-occupied state, an issue discussed below[38]) are not equivalent to quasi-particle excitation energies, even within exact DFT.[39] However, it has been shown theoretically that DFT eigenvalues can be reasonable, often quantitatively useful approximations to quasi-particle excitation.[40−42] Nevertheless, conventional approximations to DFT do not fulfill exact physical conditions for the asymptotic decay of the exchange-correlation (xc) potential of either molecule or metal. As elaborated below, these are essential for achieving reliable energy level alignment in DFT.

In this Letter, we develop and apply a first-principles method within the DFT framework that allows for a quantitatively reliable calculation of energy level alignment in physisorbed molecule–metal systems. The method extends the “DFT+Σ” approach, grounded in many-body perturbation theory, which introduces approximate gas-phase and image charge self-energy terms to correct the level alignment obtained from conventional density functional theory.[34,43−46] Specifically, our current approach rests within this framework but is based on fulfilling exact properties of the xc-potential, directly affecting the ionization potential (IP) of the gas-phase molecule and the image-plane position of the metal. For the molecule, we use an optimally tuned range-separated hybrid (OT-RSH) functional to achieve highly accurate ionization energies and, more generally, outer-valence excitation energies and affinities in the gas phase. We then determine a nonclassical image plane position from the DFT calculation of the metal, from which the substrate-induced polarization energy is deduced using an electrostatic model. By comparison with experimental and theoretical literature data, we show that the combination of these two physically motivated calculations result in a computationally inexpensive DFT+Σ approach, which (for reasons explained below) we call DFT+Σaxc, that yields quantitatively useful energy level alignment information.

To illustrate the limitations of conventional DFT and test our new approach, we rely on a set of several prototypical molecule–metal interfaces for which comprehensive experimental and theoretical data are already available in the literature. These are the interfaces of 1,4-benzenediamine (BDA),[29,34,36,47] Cu–phthalocyanine (CuPc),[48,49] and 3,4,9,10-perylene-tetracarboxylic-dianhydride (PTCDA)[50−55] with the Au(111) surface, as well as the interface between benzene (Bz) and the semimetal graphite(0001).[33,36] These interfaces, the structures of which are shown in Figure 1,[56] are well-suited test cases for our approach as they are physisorbed systems, that is, they are weakly coupled to the substrate. This means that for these systems, complex system-dependent chemical interactions, including orbital hybridization and/or charge transfer between substrate and molecule, are largely absent. This allows for a well-defined distinction between metallic and molecular states, thereby facilitating our assessment of the theoretical methodology. Moreover, for the three molecule–Au(111) interfaces, we use geometries that were previously determined using dispersion corrected DFT and found to be in excellent agreement with experiment.[47,49,54] This way, we can exclude discrepancies arising from inadequate structural information.

Figure 1

Schematic structural representation of the studied systems: (a) benzene (Bz) on graphite(0001), and (b) BDA, (c) CuPc, and (d) PTCDA on Au(111). Thin black lines indicate the surface unit-cell.

Figure 2 shows the density of states projected onto the molecule (PDOS) for all four systems, calculated within the Perdew–Burke–Ernzerhof (PBE) form[57] of the generalized gradient approximation (GGA), using the plane-wave code VASP[58] (see Supporting Information for full technical details of the band-structure calculations). Comparing the PBE–PDOS with experimental and GW reference data, it can be seen that the agreement of the calculated energy level alignment with the reference data is erratic and strongly system-dependent: Misalignment of the PBE–PDOS with the mean value of the reference data ranges from being virtually zero in the case of CuPc–Au(111) to a devastating 1.3 eV for BDA–Au(111) and 1.2 eV for Bz–graphite(0001), with the difference for PTCDA–Au(111) lying within this rather large error bar. These results are in excellent agreement with results reported in the literature,[29,33,36,49,53] where the quantitative failures of conventional local and semilocal functionals for energy level alignment are well documented. Here, we emphasize that even the chemical trend in energy level alignment extracted from PBE-calculated PDOS is entirely different than that found in experiment: PBE predicts the highest-occupied energy level of BDA to be much closer to the Fermi-energy than the highest-occupied states of CuPc and PTCDA, which is in strong contrast to the trend observed in PES experiments (see Figure 2). Combined with the large quantitative deviations, this strongly impedes any reliable PBE-based prediction of energy level alignment.

Figure 2

Density of states projected onto the molecule (PDOS), as calculated using the PBE (green thin solid line), HSE (red thin dashed line) and DFT+Σaxc (blue thick solid line) approximations, for all studied systems. Shaded areas denote reported experimental PES data (from ultraviolet PES[34,49,52] in blue and from resonant X-ray PES[34] in yellow) and theoretical data (from GW calculations[33] in green) for the highest-occupied state and indicate their peak widths or uncertainties, respectively.

Figure 2 also shows the PDOS calculated with the screened range-separated hybrid (RSH) functional suggested by Heyd, Scuseria, and Ernzerhof (HSE).[59,60] The fraction of exact exchange included in hybrid functionals mitigates differences in orbital self-interaction errors.[61,62] Therefore, HSE is known to improve the description of electronic properties of molecules[63−67] while retaining similar accuracy as PBE for metals.[36,68] Indeed, the agreement with the reference data in Figure 2 is somewhat improved over PBE-DFT calculations. For Au(111)–BDA, this is in agreement with ref (36). However, as with PBE also in HSE the lowest electronic energy barrier among the gold-molecule systems is predicted for BDA, which is in strong disagreement with experiment. Clearly, even when applying computationally demanding hybrid DFT functionals, chemical trends in energy level alignment can be incorrect.

The discrepancies with experiments are due to two separate issues in DFT.[33,36] First, for the unknown exact DFT xc-potential it can be shown that the energy of the highest occupied electronic state equals the negative of the IP,[69,70] which is known as the IP-theorem. In approximate DFT functionals such as PBE and HSE this criterion is not obeyed. We demonstrate this well-known issue in Table 1, which reports the PBE and HSE calculated energy of the highest occupied molecular orbital (HOMO) of Bz, BDA, CuPc, and PTCDA in the gas phase, again compared to high-level PES and GW literature data. Clearly, the HOMO energy of each system calculated with either PBE and HSE strongly underestimates the IP, confirming a well-documented deficiency of many approximate DFT functionals.[40] Furthermore, the result is not even qualitatively correct, in the sense that the predicted ordering of the IPs is inconsistent with experiment.

Table 1

Energy of the Highest-Occupied Molecular Orbital (HOMO) of Bz, BDA, CuPc, and PTCDA in the Gas Phase, Calculated with PBE, HSE, OT-RSH, and GW (Taken from References (30, 88, 97 and 107)), Compared with Experimental Values (Taken from References (61 and 112−114)) for the Ionization Potential (IP)a

	PBE	HSE	OT-RSH	GW	experimental IP
Bz	–6.3	–6.9	–9.3	–9.4	9.3
BDA	–4.2	–4.8	–7.1	–7.3	7.3
CuPc	–4.9	–5.0	–6.2	–6.2	6.4
PTCDA	–6.1	–6.5	–8.1	–8.0	8.2

All quantities are given in units of eV.

Second, upon molecule–metal contact the molecular ionization energy and, ergo, the energy level alignment will change because of the polarization response in the metallic substrate, an effect known as substrate-induced renormalization of the molecular energy levels.[33,71−73] As a function of molecule–metal distance, the energy change of the molecular levels relative to EF follows a classical imagelike potential over a wide spatial range,[33] which is why renormalization is sometimes referred to as an “image-charge effect”. The importance of incorporating such effects into theoretical calculations of conductance has been repeatedly demonstrated,[16,43−46,74−76] and the impact of renormalization on conductance in molecular-electronic devices has recently also been studied experimentally.[77] In DFT functionals such as PBE and HSE, which are based on local correlation only, renormalization must be absent because from both a DFT and a GW perspective it is a nonlocal correlation effect.[33,36,71−73]

Taken together, these deficiencies impede quantitatively useful energy level alignment calculations within DFT.[33,36] However, it is often observed that despite their shortcomings with respect to orbital energies, the same classes of functionals do predict bare-substrate work functions and adsorbate-induced dipole effects accurately.[35,49,53,78,79] This suggests that the calculations surveyed so far provide for an electron density distribution that is sufficiently accurate to allow for nonself-consistent correction schemes built upon it. The central aspect of the new DFT+Σ approach suggested here is to correct the energy level alignment errors in a two-step procedure without any input from higher-level theory or experiment.

The first step in our DFT+Σ approach is to correct the orbital energies for the gas-phase molecules. Here, we accomplish this by employing an OT-RSH functional. In the OT-RSH approach, a range-separation of the Coulomb operator is performed[80,81] with correct Fock-exchange in the long-range and compatible xc in the short-range. The optimal range-separation parameter, γopt,[82−94] is then found by tuning it per system until the IP theorem is satisfied, that is, εHγ = −IPγ,[87] where εH is the HOMO energy and IP the ionization potential calculated from DFT total-energy differences (see Supporting Information for details of the optimal tuning scheme). This procedure is entirely nonempirical and fully embedded in DFT within the generalized Kohn–Sham framework.[87,95] The HOMO energies of the four molecules in gas phase calculated with this optimally tuned RSH (OT-RSH) approach in the QChem code[96] are also shown in Table 1 (see Supporting Information for full technical details of the gas-phase calculations). In line with previous work on molecules,[82−94] the OT-RSH HOMO energies are in very good agreement with the reference data with a remaining error bar of 0.2 eV, which is on par with GW. Here, we employed a short-range Fock-exchange fraction of 20%,[85] because this can additionally mitigate possibly severe orbital-ordering errors in PBE, associated with treating localized and delocalized orbitals on the same footing and arising from self-interaction errors. This has been shown for various gas-phase molecules in general[40,88,97,98] and specifically by Liu et al., who used a similar correction term in the context of charge transport in metal–molecule junctions.[46]

The second step is to calculate the impact of renormalization on energy level alignment. Our model for that is based on previous DFT+Σ work that used an image-charge expression with a classical image plane obtained from the DFT-determined self-consistent response of the metal slab to an external electric field,[99] as a simplified self-energy that accounts for renormalization.[33,43−46] Here, we use a different and more general route to determine the image plane, inspired by fundamental work on the xc-potential of metallic surfaces,[70,100] the xc-hole at a metal surface,[101] and GW calculations on jellium and metal surfaces.[102−104] Specifically, we know that the exact xc-potential far away from a metallic surface must decay in an image-potential manner,[70] that is, in atomic units Vim(z) = −1/[4(z – z0)], where z is the coordinate perpendicular to the surface and z0 is the image-plane position. We also know that the asymptotic form of the PBE potential decays exponentially instead. An important hint for overcoming this deficiency has been given by Eguiluz et al.[102,103] They have used the Sham-Schlüter equation to show that the xc-potential compatible with the self-energy obtained from a GW calculation is (1) virtually the same as that obtained from a (semi)local functional inside the metal and (2) connects seamlessly to the correct asymptotic image potential well outside the metal. This can be used in a very simple manner to deduce the image-plane position, by tuning it such that the plane-averaged curves of the PBE xc-potential and image potential curves have a common tangent point, as suggested in ref (105). Importantly, this way we effectively realize a range-separated potential for the metal surface that has the correct long-range form but with remaining compatible exchange and correlation in the short-range. Furthermore, as discussed in detail on the basis of GW calculations in refs (102−104), the more general image-plane position thus obtained is that of a “quantum-mechanical electron” and not the one of a “classical point charge” obtained from the self-consistent response procedure employed in previous DFT+Σ work.[106] This means that we intrinsically include the effect of exchange-correlation interactions on the effective image-plane position. These can be significant but are absent in the image-plane position derived previously.[103,104] Matching the xc-potential with the image potential is also computationally simpler: it is based on a single DFT calculation, whereas the self-consistent response technique implies an additional series of DFT calculations with explicit external electric fields or test charges.[99]

The above procedure is demonstrated in Figure 3 (see Supporting Information for full numerical details) for the Au(111) surface. As shown there, the PBE xc-potential decays too rapidly outside the metal. The image-plane position obtained through the common tangent point procedure is found to be 0.9 Å away from the geometrical edge of the surface of Au(111) [cf. Figure 3] and 0.7 Å away from the one of graphite(0001). Note that due to the finite numerical accuracy intrinsic to the involved calculations, the error in the image-plane position obtained from this approach is estimated to be ±0.1 Å.

Figure 3

Plane-averaged exchange-correlation (xc) potential, obtained from PBE calculations for the Au(111) surface (thick black line). The origin of the x-axis is set to the geometric edge of the slab. Curves from a classical image-charge model with image-plane values, z0, of 0.9 (thick dashed line), 0.5 (thin dotted line), and 1.5 Å (thin solid line) are also shown.

In our proposed DFT+Σaxc approach, (where axc denotes “asymptotic exchange-correlation”) we build on the fact that both PBE and HSE describe short-range exchange and correlation reasonably well, but lack the above-discussed asymptotic exchange and nonlocal correlation contributions that are absolutely crucial for accurate energy level alignment calculations. Following the DFT+Σ framework,[33,34,43−46] we augment the PBE–PDOS, obtained for the full metal–molecule interface from a plane-wave based calculation (see Supporting Information), by shifting each molecular resonance individually using a “one-shot” nonself-consistent correction combining the above two steps. This is a model self-energy correction, but extending previous work it is premised entirely on DFT quantities. For occupied orbitals, it is given (in atomic units) byHere, ε is the ith electronic state, d is the molecule–metal distance, and z0 is the image-plane position. The first term on the right-hand side of eq 1 describes the uncorrected molecule-related PBE eigenvalues in the combined system. The second term is the OT-RSH-based correction of the PBE eigenvalues, obtained from a gas-phase calculation (see Supporting Information) of the molecule. The third term is the impact of substrate renormalization on the molecular levels, obtained from the classical image-charge model with z0 tuned as described above.

The results for the gas-phase correction are easily obtained from Table 1 and range from −1.3 eV for CuPc to −3.0 eV for Bz. The results for the renormalization correction are given in Table 2 and span a modest range, from +1.4 eV for BDA on gold to +1.6 eV for PTCDA on gold. This shows that the two ingredients of the DFT+Σaxc correction are opposite in sign and thus have a tendency to partially cancel.[36,49] However, full cancellation is neither expected in principle nor obtained in practice.[36]

Table 2

Molecule–Metal Distance, d, and Calculated Renormalization Correction Energy for the Occupied Levels of the Studied Systems

	d (Å)	renormalization (eV)
Bz–graphite(0001)	3.24	1.4
BDA–Au(111)	3.50	1.4
CuPc–Au(111)	3.21	1.6
PTCDA–Au(111)	3.18	1.6

In Figure 2 we show the DFT+Σaxc results, calculated with eq 1 for all four molecule–metal interfaces, as thick blue lines. The PDOS is found to be in very good agreement with the experimental and theoretical reference data, as the energy level alignment is predicted by DFT+Σaxc either within or very close to the reported peak widths and error bars of the reference data for all molecule–metal systems. Interestingly, for the cases where PBE already performed reasonably well, the DFT+Σaxc corrections are minor. For Bz-graphite(0001) and BDA-Au(111), however, the calculated correction is on the order of 1 eV and therefore absolutely crucial for accurate energy level alignment prediction. This large correction is also in good agreement with previously reported GW calculations, which were feasible for these moderately complex molecule–metal systems.[33,36] These findings show that the DFT+Σaxc approach correctly reflects the highly system-dependent nature of the physical effects that govern energy level alignment. Furthermore, it correctly reproduces the experimentally determined chemical trends for the molecule–gold interfaces, which highlights its capabilities for reliable predictions of energy level alignment. Importantly, because the corrections of eq 1 are applied per orbital, the suggested scheme does more than just provide a rigid shift of the PBE-computed data, generated by a “scissors operator”. In particular, for CuPc, where PBE is known to produce an incorrect character for the HOMO orbital owing to severe self-interaction errors,[107] the OT-RSH correction term restores the level ordering obtained in gas-phase GW.[88] Therefore, while in Figure 2 the difference between PBE and DFT+Σaxc may superficially appear to be small, it in fact reflects not only a quantitative change in level alignment but also a qualitative change in the nature of the aligned level. We further comment that while here we focused on HOMO resonances for the sake of comparison to (direct) PES, the method is general and potentially applicable to unoccupied resonances (in particular that of the lowest unoccupied molecular orbital), which could be studied with, e.g., inverse PES.

Before summarizing, let us also discuss the limitations of our approach. First, as with all computational methods, the numerical accuracy of our computational methodology (e.g., due to the finite basis set for the molecular and molecule–metal calculations) must be monitored for convergence and may limit the expected quantitative agreement of our approach. Second, our electronic structure calculations will generally also depend on the quality of the geometrical structure of the molecule–metal junction. For example, for the CuPc–Au(111) interface the geometry we have used is by ∼0.2 Å closer to the metal surface than that measured in experiment.[49] Moreover, for the same system it has been shown that energy level alignment significantly depends on the lateral coverage.[48] As a second example, for the BDA–Au(111) interface the computed geometry has been found to depend on the theoretical treatment of van der Waals interactions.[108] We note, however, that fine geometrical details are as challenging to capture experimentally as they are theoretically and are thus well beyond the scope of this paper.

Beyond the above standard limitations, we emphasize that all our calculations are performed for flat-lying, physisorbed molecules for which the character of the molecular electronic states in the combined molecule–metal system is similar to the molecular orbitals in gas-phase and there is a well-defined molecular resonance. Therefore, a “one-shot” correction of the molecular resonance energies in the spirit of eq 1 is indeed applicable. Furthermore, a unique image-plane value can be used. Additionally, one could generalize the approach by considering the spatial distribution of the corrected orbitals using, for example, Mulliken charges on each atom, and averaging over the molecule–substrate distance accordingly.[44−46,105] More generally, the polarizability of the molecule may need to be considered so as to obtain a self-consistent screening response,[33,109,110] or indeed, molecular resonances may appreciably overlap the substrate Fermi level, necessitating treatment of dynamic screening beyond the static polarization.[33] Last but not least, generalization to chemisorbed systems presents additional challenges, because the character of the molecular states may be altered, which renders the identification of resonances in the combined system as states arising from specific molecular orbitals more difficult.[111] Projection schemes have previously been used for DFT+ Σ calculations of chemisorbed molecule–metal systems[111] and we believe that projection-based approaches may also be of assistance here. The theoretical framework suggested here is, therefore, already useful for many cases of interest and holds promise as a basis for generalized methods that could provide useful predictions in more complex scenarios.

In summary, we have developed and applied the DFT+Σaxc method for energy level alignment calculations of physisorbed molecule–metal interfaces. In our new approach, both terms of the model GW-like self-energy correction are based on DFT and anchored in the fulfillment of exact asymptotic criteria for the xc-potential of molecule and metal. A central tenet of our approach is that asymptotic exchange and nonlocal correlation components are absolutely crucial for accurate energy level alignment prediction and need to be included in the theoretical treatment of molecule–metal interfaces and molecular junctions. In practice, accounting for these essential contributions was achieved by combining an optimally tuned range-separated hybrid functional in the gas phase with a DFT-based image-charge model for the metal, which restores the asymptotically correct xc-potential. Our detailed analysis and comparison with experimental and theoretical reference data has shown that our approach reproduces energy level alignment with high accuracy and thus correctly captures chemical trends at molecule–metal interfaces that are out of the reach of conventional DFT approximations. Given that the computational cost of our approach is small compared to the underlying DFT calculation of the complete molecule–metal interface, these findings highlight DFT+Σaxc as a reliable and efficient theoretical tool for calculating energy level alignment.