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.


Technical details of the molecule-metal and bare metal calculations
Our molecule-metal and bare metal calculations rely on density-functional theory (DFT) as implemented in the electronic-structure program VASP. 1 Valence-core interactions were treated within the projector-augmented-wave 2 (PAW) formalism using regular ("normal") PAW potentials for all elements. The valence-electron wavefunctions were expanded in a plane-wave basis (kinetic-energy cutoff: 400 eV). Along the surface normal of the moleculemetal calculations, we inserted a vacuum gap of ~20 Å and applied a dipole correction. We employed a convergence criterion of 10 -4 eV for the total energy per unit cell in the electronic self-consistent cycle.
As mentioned in the main text, the coordinates for the molecule-metal systems are taken from previous computational studies that were based on dispersion-corrected DFT. For the calculations with the HSE06 functional, we performed a three-step procedure to facilitate electronic convergence: (i) We performed a PBE single-point calculation without dipole correction; (ii) using the wavefunctions of step (i) as a starting guess, we carried out an HSE06 calculation without dipole correction, for which the evaluation of the Fock-exchange part in the exchange-correlation potential on the reciprocal grid was reduced in the k x -and k ydirection (each by a factor of two); and (iii) using the wavefunctions of step (ii) we performed HSE06 calculations with dipole corrections using the standard k-point grids (vide supra).
For determining the optimal image-plane position (see Fig. 3 and associated discussion in the main text), we note that it is well-known that in the vacuum the exchange-correlation (xc) part converges much more slowly than the electrostatic part of the total potential. Therefore, when computing the local potential of bare graphite(0001) and Au(111), we employed more strict settings in the band-structure calculations. Specifically, we increased the number of FFT-grid points and enhanced the representation of the charge density, potential and augmentation 3 charges. Furthermore, we have employed an additional support grid for the augmentation charges, which further increases the accuracy of calculated quantities.
To compute the local xc-potential of these two inorganic slabs, we then calculated the total local potential and the electrostatic part of it separately. The final plane-averaged xc-potential we used for the fitting procedure was then obtained by subtracting the plane-averaged electrostatic part of the potential from the plane-averaged total potential of the slabs.
As can be seen in Fig. 3 of the main text, the extracted plane-averaged xc-potential of the surfaces still exhibits some numerical noise in vacuum. Note that this occurs rather far away (ca. 4 Å) from the surface and is therefore not expected to be detrimental to our fitting procedure. To guarantee that the fit of the image-plane potential to the xc-potential is not affected by this numerical noise, we smoothed the extracted plane-averaged xc-potentials of graphite and gold with a cubic-spline interpolation. This allows us to find the "optimal" image-plane position, z 0 , in a well-defined way by sweeping z 0 in the image-charge potential until the value and derivative of the xc-and image-potential coincide at one point. Due to the finite number of points used in sampling the real space grid, the numerical accuracy of this approach is estimated to be ±0.1 Å.

Technical details of the tuning scheme and gas-phase calculations
For the formalism and details of the optimally-tuned range-separated functional (OT-RSH) procedure, we refer the reader to Refs. 5,6. All gas-phase calculations presented in this article were obtained within the QChem code. 7 We used the cc-PVTZ 8 basis functions for BZ, CuPc and PTCDA. For the BDA molecule, we employed the aug-cc-PVQZ basis, for which we found that the cc-PVTZ basis set was not sufficient for convergence, in accordance with ref.
9. For the gas-phase calculations, we used geometries that were previously published: BDA