Energy Level Alignment at Molecule-Metal Interfaces from an Optimally-Tuned Range-Separated Hybrid Functional
Zhenfei Liu, David A. Egger, Sivan Refaely-Abramson, Leeor Kronik, Jeffrey B. Neaton
EEnergy Level Alignment at Molecule-Metal Interfaces from anOptimally-Tuned Range-Separated Hybrid Functional
Zhen-Fei Liu,
1, 2, a) David A. Egger, a) Sivan Refaely-Abramson, Leeor Kronik, b) and Jeffrey B. Neaton
1, 2, 4, c) Molecular Foundry and Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley,California 94720, USA Department of Physics, University of California, Berkeley, California 94720,USA Department of Materials and Interfaces, Weizmann Institute of Science, Rehovoth 76100,Israel Kavli Energy Nanosciences Institute at Berkeley, Berkeley, California 94720,USA (Dated: 7 March 2017)
The alignment of the frontier orbital energies of an adsorbed molecule with the substrate Fermi level at metal-organic interfaces is a fundamental observable of significant practical importance in nanoscience and beyond.Typical density functional theory calculations, especially those using local and semi-local functionals, oftenunderestimate level alignment leading to inaccurate electronic structure and charge transport properties. Inthis work, we develop a new fully self-consistent predictive scheme to accurately compute level alignment atcertain classes of complex heterogeneous molecule-metal interfaces based on optimally-tuned range-separatedhybrid functionals. Starting from a highly accurate description of the gas-phase electronic structure, ourmethod by construction captures important nonlocal surface polarization effects via tuning of the long-rangescreened exchange in a range-separated hybrid in a non-empirical and system-specific manner. We implementthis functional in a plane-wave code and apply it to several physisorbed and chemisorbed molecule-metalinterface systems. Our results are in quantitative agreement with experiments, both the level alignment andwork function changes. Our approach constitutes a new practical scheme for accurate and efficient calculationsof the electronic structure of molecule-metal interfaces.
I. INTRODUCTION
Interfaces between molecules and metals play a cen-tral role in emerging functional devices in nanoscienceand nanotechnology . When a molecule is adsorbedon a metal surface, several important physical and chem-ical phenomena occur. For example, an interface dipoleforms, altering the work function of the metal surface ;discrete molecular frontier orbitals hybridize with ex-tended metallic states, forming molecular resonances;and substrate screening effects shift orbital energies. Akey physical observable is the energy level alignment be-tween the frontier molecular resonance peak positionsand the Fermi level, E F , of the metal; this alignmentcan be directly linked to the energy barrier and effi-ciency of charge transfer across the interface. For ex-ample, in molecular junctions, highest occupied molec-ular orbital (HOMO) or the lowest unoccupied molec-ular orbital (LUMO) resonance energies relative to E F are central to determining the zero-bias conductance ofthe junction, as shown in Ref. 12. Molecule-metal bond-ing spans a range of binding energies and degrees ofhybridization strengths , from weak physisorption tostrong chemisorption; in many cases where covalent ormolecule-metal interactions associated with significant a) These authors contributed equally. b) Electronic mail: [email protected] c) Electronic mail: [email protected] charge transfer occur, molecular signatures in the inter-face electronic structure are nonetheless observed. Forthese intermediate cases, molecular resonances can beenergetically close to or at the metal Fermi level, whichis generally an indication of strong charge transfer be-tween the two systems and of significant and complexchanges to the interface dipole. Because level alignmentand charge transport are intertwined, its accurate de-scription is of general importance for understanding, con-trolling, and predicting functional properties at a varietyof interfaces, including systems related to energy conver-sion and storage.Experimentally, energy level alignment can be deter-mined by direct photoemission spectroscopy for occupiedmolecular orbitals, and by inverse photoemission spec-troscopy for unoccupied molecular orbitals. Conduc-tance measurements of molecular junctions probecharge transport properties and, thus, include indirectinformation about the level alignment. However, in prin-ciple, distinct binding geometries of molecules in junc-tions can lead to different level alignment and, hence,charge transport properties , and as charge transportmeasurements are usually ensemble averages of multiplegeometries, care must be taken in relating conductanceto level alignment in such measurements .First-principles electronic structure calculations thatcan model individual, well-defined geometries provideadditional information complementary to experiments.From a formal theory viewpoint, molecular levels at in-terfaces are quasiparticle energy levels, i.e., they corre-spond to charged excitations. A rigorous formalism for a r X i v : . [ c ond - m a t . m e s - h a ll ] M a r quasiparticle energies is many-body perturbation theory(MBPT), which in practice is most often used in the GW approximation , where G is the single-particleGreen’s function and W is the screened Coulomb interac-tion. GW has been shown to be an accurate approach fora wide range of molecules and bulk solids . How-ever, several issues have hindered the widespread use of GW calculations for molecule-metal interfaces. Firstly,due to its high computational cost, even with today’scomputing power it is still far challenging to perform GW calculations of several hundred atoms with peri-odic boundary conditions, as is the case for molecule-metal interfaces . Moreover, several benchmarkingstudies, including work on smaller molecule-metal sys-tems, also showed that results from perturbative GW calculations can be challenging and expensive to convergenumerically . Furthermore, it is by now well-knownthat single-shot GW calculations often depend on theunderlying starting point, which adds additionalcomplications for efficient GW calculations of large-scalemolecule-metal interfaces and their functional properties.Balancing accuracy and efficiency is crucial for calcula-tions of molecule-metal interfaces, and density functionaltheory (DFT) is usually the pragmatic choice forfirst-principles calculations at relatively moderate com-putational cost, provided it is accurate enough. For-mally, however, eigenvalues of Kohn-Sham (KS) Hamilto-nians are only zeroth order approximations to quasiparti-cle energies . In fact, there is no theorem guaranteeingthat they are quantitatively accurate, with the impor-tant exception of the HOMO energy, for which the ion-ization potential (IP) theorem holds . For the bandgap of a semiconductor (or, equivalently, the molecularHOMO-LUMO gap), even the exact functional in KSDFT is not necessarily accurate, as it lacks a derivativediscontinuity in the exchange-correlation (XC) poten-tial owing to its strictly local, multiplicative nature .Fig. 1 schematically shows representative energy levelsfor the molecule-metal interface. The right hand side il-lustrates the situation for an isolated gas-phase molecule,for which there is a well-defined quasihole level denotedas HOMO, which is equal to the negative of the IP, and awell-defined quasielectron level denoted as LUMO, whichis equal to the negative of the electron affinity (EA). Forthe molecule adsorbed on the metal surface, in Fig. 1we denote resonances that appear as peaks in the pro-jected density of states (PDOS) - where the DOS hasbeen projected onto the molecular subspace - as HOMOand LUMO, which are now broadened due to hybridiza-tion. Comparing the molecular energy levels of the gasphase and the adsorbed molecule, Fig. 1 illustrates theimportant physical phenomenon of surface-induced gaprenormalization : when the molecule is close to thesurface, electrons in the metal respond to and screensingle-particle excitations in the molecule, i.e., free car-riers in the metal polarize when a hole/electron is addedto the molecule and screen the Coulomb interactions inthe molecule. As a result, when the molecule approaches E E F metal surface gas-phase molecule exact HOMO (= -IP)exact LUMO (= -EA)adsorbed molecule z . . - - . - - . . . vacuum level FIG. 1. Relevant energy levels at a molecule-metal interface.Green lines (shaded areas) represent exact quasiparticle levels(PDOS of the adsorbed molecule). Note that for the case ofthe adsorbed molecule, broadened molecular resonances areshown. Dashed arrows indicate the surface-induced renor-malization of the molecular levels, which depends on themolecule-metal distance, z . The sign of the vacuum levelchange at the interface, a result of the interface dipole, istaken here to be negative. For easy visualization, we alignthe interface-modified vacuum level with that of the gas-phasemolecule. the surface, the HOMO and LUMO energies move closerto E F and the HOMO-LUMO gap is reduced. Forweakly coupled physisorbed molecules, the energy changeof the molecular levels relative to E F as a function of themolecule-metal distance has been shown to follow a clas-sical image-like form , and renormalization is sometimesreferred to as an “image-charge effect”. The effective“image potential” is given by 1 / [4( z − z )], where z is theaverage height of the molecule on the surface and z is theimage plane position. It is an effective one-body poten-tial resulting from a many-body effect, namely non-localcorrelation, rather than from a static perturbation to thesystem. Fig. 1 also shows how the vacuum levels arealigned at the interface and for the gas-phase molecule.The vacuum level can increase or decrease at the inter-face (only one case is shown in Fig. 1), due to formationof the interface dipole.Results from calculations using the local density ap-proximation (LDA) and typical generalized gradi-ent approximations (GGAs) in the Kohn-Sham (KS)DFT scheme, while often accurate for total energy orcharge density related properties of molecule-metal in-terfaces such as work-function changes , typicallystrongly underestimate the band gap of molecules andsemiconductors . We illustrate this with grey dashedlines in Fig. 2, namely that LDA/GGA misplace theHOMO and LUMO levels in the gas phase, such that thegap is underestimated. Moreover, in typical LDA/GGAcalculations of molecule-metal interfaces, one finds thatthe HOMO-LUMO gap for the adsorbed molecule re-mains virtually the same as in the gas phase. In otherwords, LDA and GGA results for the level alignment arenot sensitive to changes in screening environment and fail . . - - . - - . . . E E F metal surface gas-phase molecule exact HOMO (= -IP)exact LUMO (= -EA)adsorbed molecule z OT-RSH HOMO with tuned β OT-RSH LUMO with tuned β LDA/GGA LUMOLDA/GGA HOMOvacuum level
FIG. 2. Energy levels at a molecule-metal interface. Greenlines represent exact quasiparticle levels, grey dashed lines ap-proximate LDA/GGA levels, and red lines OT-RSH levels ascomputed with the method introduced in this paper (see textfor details). Note that for the case of the adsorbed molecule,broadened molecular resonances are shown. Dashed arrowsindicate the surface-induced renormalization of the molecularlevels, which depends on the molecule-metal distance, z . Thesign of the vacuum level change at the interface, a result ofthe interface dipole, is taken here to be negative. For easy vi-sualization, we align the interface-modified vacuum level withthat of the gas-phase molecule. to capture the gap renormalization effect. This is becauseLDA and GGAs consist of local correlation only whereasrenormalization is a nonlocal effect . As we discussbelow, similar findings hold for conventional global andscreened hybrid functionals.The high computational cost and convergence chal-lenges associated with GW , together with the formaland practical issues of popular XC functionals in pre-dicting level alignments, led to the development of severalcorrection schemes following various directions, see Refs.50,51,57,60–69 for examples. Specifically for physisorbedsystems, one can assume weak coupling between themolecule and the metal, and attempt to shift the PDOSpeaks in a non-self-consistent, adjustable-parameter-free a posteriori manner. One example for such a non-self-consistent correction scheme is the DFT+Σ approach ,where Σ denotes a two-step correction to KS eigenval-ues: (1) it uses a gas-phase correction, equal to thedifference between the gas-phase LDA/GGA HOMO orLUMO energy and a value obtained from a more accurateapproach (e.g., GW ) that is closer to the exact quasi-particle energy (cf. Figs. 1 and 2); and (2) a correc-tion accounting for the surface polarization responsiblefor gap renormalization, which is often approximated us-ing a classical image-charge model , 1 / [4( z − z )], where z is the average distance between the adsorbed moleculeand the surface, and z is the image-plane position. Wenote that recently we have augmented this approachby using a non-classical DFT-determined image planeof the metallic surface to compute the surface polariza-tion term, rather than the classical DFT-derived imageplane . The DFT+Σ method has been used to success-fully predict and explain level alignment at physisorbed interfaces and charge transport in specific molecularjunctions where the above-mentioned weak-coupling as-sumption is reasonable .Of major interest are approximate GW -inspired theo-retical approaches that can go beyond the weak-couplingassumption to also treat molecules chemisorbed on met-als, i.e., when a significant amount of charge is trans-ferred between the two subsystems. Examples in-clude systems for which molecular levels are pinnedat E F71,72 , or where significant covalent interactionslead to strong hybridization between molecular andsubstrate orbitals , such as the Au-benzenedithiol-Aujunction . As mentioned above, this situation canlead to the molecular resonances being energetically closeto the Fermi level, which implies a significant chargetransfer. Clearly, the amount of charge transfer will de-termine the distribution of the charge density of the in-terface system, which means that a correction scheme forlevel alignment requires self-consistency. This notion alsoimplies that in case of chemisorption an incorrect descrip-tion of level alignment can yield errors in the density and,at least in principle, result in incorrect interface dipoles, work functions, and total energies of molecule-metal sys-tems. Therefore, we currently seek a theoretical methodthat can accurately and efficiently characterize the in-terface electronic structure in a self-consistent manner,and thus is applicable to physisorbed and, in particular,also to more strongly bound and partially chemisorbedmolecule-metal systems. With such an approach in hand,one is able to improve the prediction of relevant interfaceproperties beyond level alignment, including work func-tion changes and PDOS lineshapes.Here, we present a new method aiming in the veinof reliable and self-consistent determination of the levelalignment at molecule-metal interfaces. Our concep-tual framework is the generalized Kohn-Sham (GKS)scheme , which is still a DFT approach but in con-trast to LDA/GGA builds on an effective potential that isnonlocal. Hybrid functionals, in which (semi)local func-tionals are mixed with a fraction of Fock exchange, area special case of GKS DFT. Popular hybrid functionalswere shown to greatly improve GGA results for both to-tal energy related properties, such as thermochemistry ,as well as quasiparticle energy levels, such as bandgaps . However, conventional hybrid functionals do notin general remedy the energy level alignment problemat interfaces . Notably, conventional hybrids can-not capture the renormalization of the HOMO-LUMOgap. A relatively recent and particularly promising classare optimally-tuned range-separated hybrid (OT-RSH)functionals , in which parameters of the functionalare tuned per system to satisfy important physical con-ditions, without recourse to empirical data or any fitting.OT-RSH was shown to yield accurate frontier orbital en-ergy levels, outer-valence electron spectra, and HOMO-LUMO gaps of gas-phase molecules and molecularcrystals , as well as transport properties . In thiswork, we extend the applicability of this functional toheterogeneous molecule-metal interfaces by proposing aroute for a judicious choice of optimal parameters in theOT-RSH functional. This yields a fully self-consistentOT-RSH scheme which is applicable to both physisorbedand chemisorbed molecule-metal systems and thus ex-tends methods reliant on the weak coupling limit. Usingit for several prototypical test cases, we achieve quantita-tive agreement with experiments for level alignment andwork function changes including those featuring chargetransfer and stronger bonding.The outline of this paper is as follows: In Sec. II, webriefly review the OT-RSH functional as applied to gas-phase molecules and molecular crystals, and then presentour approach for treating molecule-metal interfaces. InSec. III, we apply the method to six molecule-metal in-terfaces that have been well-studied in both theory andexperiment, including two systems where charge transferand Fermi level pinning occur. In Sec. IV, we discusslimitations of the proposed method and outline remain-ing challenges, which is followed by our conclusions inSec. V. II. METHODOLOGYA. OT-RSH for gas-phase molecules and molecularcrystals
In RSH functionals, the Coulomb operator is decom-posed into short-range and long-range components . Inthis approach, proposed by Yanai et al. , the decompo-sition takes the form1 | r − r (cid:48) | = α + β erf ( γ | r − r (cid:48) | ) | r − r (cid:48) | + 1 − [ α + β erf ( γ | r − r (cid:48) | )] | r − r (cid:48) | , (1)where α, β, and γ are parameters and erf ( · ) is the errorfunction. We note that this partition is not unique, butthe choice of the error function is computationally conve-nient. Here, we treat the first term using nonlocal Fockexchange and the second using semi-local exchange. Theexchange-correlation energy can then be written as E XC = αE EXXX , SR + (1 − α ) E GGAX , SR + ( α + β ) E EXXX , LR + (1 − α − β ) E GGAX , LR + E GGAC , (2)where the subscript X (C) denotes exchange (correla-tion), SR (LR) denotes the short-range (long-range) ex-change, and the superscript EXX (GGA) reflects whetherthe corresponding energy component is treated usingFock exchange (GGA exchange or correlation). In thiswork, we follow Ref. 84 and use PBE (Perdew-Burke-Ernzerhof) for the GGA exchange ( ω PBE for the SRpart) and correlation components.For an isolated gas-phase molecule, α is often chosen tobe 0.2, and then β is chosen as 1 − α = 0 . via enforcing full long-range Fock exchange. Here, we will use the notation β = 0 . β value appropriate for isolatedgas-phase molecules. γ is the range-separation parameterand governs the separation between SR and LR. In theOT-RSH approach, the remaining parameter γ istuned separately for each system by minimizing J ( γ ) = | ε γ ( N ) + IP γ ( N ) | + | ε γ ( N + 1) + IP γ ( N + 1) | , (3)where ε are GKS HOMO eigenvalues, the IP is calcu-lated using total energy differences between neutral andcharged systems, and N ( N + 1) indicates the neutralmolecule (anion).Ref. 90 generalized the OT-RSH functional to the caseof molecular crystals. α and γ were chosen according tothe gas-phase values, but β was adapted to account forchanges in the long-range Coulomb screening due to thepresence of the other molecules in the crystal and the dif-ference in dielectric environment compared to an isolatedmolecule. From Eq. (2), we see that α + β governs thefraction of long-range Fock exchange. For a molecularcrystal with an average dielectric constant (cid:15) , Ref. 90 pro-posed to use α + β = 1 /(cid:15) <
1. This choice of β was shownto correctly describe the gap renormalization in molecu-lar crystals compared to a single gas-phase molecule andoptical absorption . B. OT-RSH for molecule-metal interfaces
We first note that for a fixed, z -independent choice of α, β, and γ , the OT-RSH functional cannot capture the1 / [4( z − z )] behavior of the image-charge effect and gaprenormalization of the adsorbate in the weak-couplingphysisorbed limit at a metal surface for all z . To un-derstand why, consider that in Ref. 50, it was shownthat the physical origin of the gap renormalization is thechange in the screened Coulomb interaction, ∆ W , be-tween the isolated molecule and the adsorbate. From aGKS viewpoint, this long-range correlation effect wouldrequire a change in the amount of long-range screenedFock exchange as a function of molecule-metal distance,which Eq. (2) with fixed parameters lacks. In fact, this isalso the case for any standard local, semi-local, or hybridfunctional, as was shown in Ref. 79 using the PBE andHSE (Heyd-Scuseria-Ernzerhof) functionals as exam-ples.As described above, in Ref. 90, β was tuned from β to 1 /(cid:15) − α in order to capture the change in Coulombscreening between an isolated molecule and a molecularcrystal. Inspired by this idea, and given the fact that theimage-charge effect is a long-range effect, we propose toaccount for the change in long-range Fock exchange ata molecule-metal interface by insisting on α + β < β so as to capture the renormal-ization of the orbital resonance energies at the level ofan image-charge model. For α and γ , it is clear that inthe limit of very weak physisorption the optimal value ob-tained in the gas phase is maintained on the surface as themolecular density remains unchanged. For chemisorbedsystems with exchange of charge between molecule andmetal, of course this does not hold anymore, but as weshow below the optimal value for γ is virtually unaffectedby small charge density rearrangements. Thus, for eachmolecule-metal system we choose α = 0 . γ basedon Eq. (3) as obtained for the gas-phase molecule. Theproblem left is how to choose the optimal β for the inter-face such that the renormalization induced by the surfaceis accounted for (see red lines of Fig. 2). Here, we choose β such that the orbital energies renormalize properly ac-cording to 1 / [4( z − z )]. To avoid a complicated func-tional form with explicit treatment of this effect, we hereuse a DFT-based image-charge model: it allows for deter-mining a non-classical image-plane position for each typeof surface in a unique way , from which we determinethe amount of orbital renormalization. If we keep β as asimple scalar parameter, it should implicitly depend on z due to the z -dependence in the image-charge model.Furthermore, we make the assumption that the orbitalenergies of the isolated molecule and the PDOS peaks ofthe adsorbate change by the same amount when β is var-ied from β to β < − α , while keeping α and γ at theirgas-phase values. Therefore we can tune β in the gasphase, which reduces the computational cost, and choose β such that the neutral gas-phase HOMO changes by anamount equivalent to the image-charge energy, such that ε γN ( β ) − ε γN ( β ) = P, (4)where subscript N denotes the neutral gas-phase systemand P = 1 / [4( z − z )] is the image-charge energy deter-mined as in Ref. 70. We note that the substrate pro-vides significant screening, and the β defined in this waycan become negative. Considering that β = − α corre-sponds to the limit of fully screened LR Fock exchange,we choose this value as our lower limit (i.e., β ≥ − . α is chosen to be 0.2). From our experience, ε γN ( β )changes linearly with β for all the molecules we studied.A cartoon for our β -tuning scheme is shown in Fig. 2,where the red line on the right indicates the gas-phaseHOMO calculated using the tuned β .Our tuning procedure, for a given molecule-metal in-terface system with z obtained from a prior DFT calcu-lation, can be summarized as follows:1. Perform a standard gas-phase OT-RSH calculationfor the molecule, i.e., use α = 0 . , β = 1 − α = 0 . γ as in Eq. (3);2. For the metal slab, compute the image-plane posi-tion, z , by matching the long-range XC potentialfrom a local or semi-local functional (PBE in thiswork) to an image potential. To be specific, firstcompute the xy -averaged PBE XC potential for agiven metal slab, V PBEXC ( z ), where z is the variabledescribing the distance away from the metal sur- face, then tune the parameter z such that the twocurves − / [4( z − z )] and V PBEXC ( z ) have a commontangent point (see Refs. 70,101,102 for details);3. Compute the approximate polarization induced bythe surface with a classical image-charge model byevaluating the expression P = 1 / [4( z − z )] (ina.u.), where z is determined from the geometry ofthe interface and is the average distance betweenthe molecule and the top layer of the metal slabalong the surface normal;4. Tune the optimal β value according to Eq. (4),using a series of OT-RSH calculations of the gas-phase molecule, with varying β and fixed α and γ ;5. Compute the electronic structure of the molecule-metal interface in a self-consistent manner usingthis OT-RSH functional with the above-determined α , β , and γ (in general, each molecule-metal systemhas its own set of optimal parameters).In this work, we stop at Step 5. However in principle, onecould re-optimize the geometry of the molecule-metal in-terface, obtain a new z , and iterate the tuning procedureto self-consistency. C. Technical details
All OT-RSH calculations of molecule-metal interfacesare performed using a modified version of QuantumESPRESSO v. 5.2.0. We implemented long-rangescreened Fock exchange based on the existing subroutinesof short-range screened Fock exchange and hybrid func-tionals, using norm-conserving pseudopotentials. This isrealized via the Fourier transform of the first term in Eq.(1), that is α + β erf ( γ | r − r (cid:48) | ) | r − r (cid:48) | F.T.=== ⇒ πq (cid:20) α + β exp (cid:18) − q γ (cid:19)(cid:21) . (5)The Gygi-Baldereschi approach is used to treat theCoulomb potential divergence at small q vectors, as al-ready implemented in Quantum ESPRESSO. A 55 Ryenergy cutoff is used for every system. Fock exchange isevaluated using a reduced k-mesh in our calculations; forsmaller systems, we verified convergence with respect tothe k-mesh. The k-mesh used for each system is specifiedin Sec. III.In the current implementation of hybrid functionals inQuantum ESPRESSO, there is an inner self-consistencyloop which solves the GKS equation with a fixed Fockoperator, and an outer self-consistency loop which up-dates the Fock operator with new GKS orbitals. Un-less noted otherwise, all results shown in this work arefrom “perturbative hybrid” calculations (which are stillself-consistent in a certain sense, see below), in whichthe Fock operator is constructed only once, using pre-converged PBE orbitals, and the resulting GKS equationis then solved self-consistently to obtain eigenvalues andorbitals without updating the Fock operator again. Theresulting eigenvalues, eigenvectors, and charge densityare already different from their PBE counterparts andthis self-consistent solution of GKS equation corrects thelarge errors of PBE in PDOS calculations, in terms ofboth lineshape and peak positions. Our tests show thatthis “perturbative hybrid” approach yields almost iden-tical PDOS around the Fermi level compared to fullyself-consistent hybrid calculations with both inner andouter loops, but is much more efficient computationallybecause we omit the outer self-consistency loop.It is well-known that for transition metals, core-valenceinteractions may have non-negligible numerical effects in GW and exact exchange-based calculations of semi-conductor band structures. In our tests, the semicore sp -states of transition metal atoms are found to affect thelevel alignment by about 0.3 eV. However, explicit inclu-sion of the semicore sp -states significantly increases thetotal number of electrons as well as the required energycutoff. In order to reduce computational costs, we usea pseudopotential (see Ref. 107 for details) with semi-core sp -states (e.g., 4s for Ag) in the valencefor the top layer of the transition metal substrate thatis closest to the adsorbed molecule, and a pseudopoten-tial without semicore sp -states (e.g., 4d for Ag) forthe reminder of the metal slab (see Ref. 108 for a noteregarding the accuracy of mixing pseudopotentials). Inaddition, we find that although total energy convergencerequires a cutoff of about 200 Ry with explicit semicore sp -states, the eigenvalues and energy level alignmentstypically converge at a much lower energy cutoff. The55 Ry cutoff we use is sufficient for the level alignmentof all the systems studied in this work.The geometries of the molecule-metal interfaces areeither adapted from the literature or relaxed usingdispersion-corrected XC functionals (see the results sec-tion for details), as implemented in VASP with theProjector Augmented Wave (PAW) method . Wethen use the relaxed molecular adsorbate geometry in ourgas-phase tuning procedures. Note that this can resultin slightly different ( < . with a cc-pVTZbasis set, and gas-phase tuning based on Eq. (4) is per-formed using NWChem with the same basis set. III. RESULTS
To demonstrate the success of our approach, we reportOT-RSH calculations of six molecule-metal interfaces.We divide them into two categories: systems withoutsignificant charge transfer (i.e., physisorbed); and sys-tems with non-negligible charge transfer (i.e., partiallychemisorbed). The latter systems feature Fermi level pin-ning in the PDOS and a significant overlap between the LUMO resonance peak and E F . As discussed in the intro-duction, the OT-RSH method has repeatedly been foundto be highly accurate for gas-phase molecules. Therefore,here we do not dwell on gas-phase results but instead fo-cus on the molecule-surface interaction.As discussed above, correction schemes that assumeweak coupling and do not update the charge density, suchas the DFT+Σ approach, are known to perform very wellfor interfaces in the first category . But a self-consistentmethod, such as the one presented in this paper, is neededfor accurate calculations of the second category. Wenote that there are specific cases of chemisorption forwhich charge transfer involves not the frontier, but en-ergetically deeper lying molecular resonances. In suchcases, our non-self-consistent scheme may work equallywell . For each system, we compare PBE and OT-RSHresults with literature results from ultraviolet photoemis-sion spectroscopy (UPS) measurements for both energylevel alignment (taken as position of the peak maximum)and work function changes. For the calculated PDOS,a 0.01 Ry broadening is applied to every system, exceptfor the one in Sec. III A 4, where a 0.02 Ry broadening isapplied in order to match the experimental PDOS width. A. Systems without significant charge transfer
In this section, we present results for four systems.To begin, we consider benzene adsorbed on Al(111) and3,4,9,10-perylene-tetracarboxylic-dianhydride (PTCDA)adsorbed on Au(111), two weakly coupled, physisorbedinterfaces for which we expect the standard, perturba-tive DFT+Σ approach to perform well. These systemsserve as a proof of concept for the approach proposed inthis paper. We then consider 1,4-benzenediamine (BDA)adsorbed on Au(111), where charge transfer is negligiblebut hybridization between molecular orbitals and Au d -states is stronger than for the first two systems. In thiscase, a self-consistent approach such as the one presentedhere may change the GGA PDOS lineshape, while apost-processing approach, such as the standard DFT+Σ,does not. Finally, we consider Anthracene-2-selenolate(AntSe) adsorbed on Au(111), an even more complicatedcase as the Se-H bond is replaced by a Se-Au bond uponadsorption, i.e., a new covalent bond is formed.
1. Benzene on Al(111)
Experimental evidence for a rather weak interactionbetween benzene and Al(111) was given in Ref. 114,where UPS in ultrahigh vacuum was used to determinethe energy level alignment. In our calculations of thissystem, we use a 4 × × × × × P DO S ( a r b . un it ) E − E F (eV) UPSPBEOT − RSH 0 2 4 6 8 10 12 14 16 − − − − − − FIG. 3. PDOS of benzene adsorbed on Al(111), whose struc-ture is shown in the inset, as obtained using PBE (grey dashedline) and OT-RSH (red line), compared with results from UPSmeasurements (green) adapted from Ref. 114.
TS approach (i.e., PBE augmented by the Tkatchenko-Scheffler scheme ) including self-consistent screening to capture the impact of the metallic screening on disper-sive interactions for a set of fixed molecular geometriesat various distances. The lowest total energy is achievedfor the molecule lying flat at 3.24 ˚A above the Al(111)surface on a bridge site.Following the tuning protocol for interfaces describedin Sec. II, Eq. (3) yields γ = 0 .
24 bohr − and a tunedHOMO energy of -9.4 eV for benzene. For Al(111), theimage plane is determined to be 1.1 ˚A above the surfaceand, using the optimized distance of our benzene adsor-bate, the image-charge energy is then 1.7 eV. In order toincorporate the polarization energy due to the metal sur-face, we tune β such that the gas-phase HOMO increasesby 1.7 eV, based on Eq. (4). This results in an optimal β value of 0 . E F , in much better agreement with exper-iment (4.0 eV) than PBE (see Fig. 3), which predictsthe resonance at 3.1 eV below E F . However, PBE andOT-RSH yield essentially the same work function changeof -0.3 eV and the same work function of 3.8 eV, sug-gesting our new OT-RSH scheme leads to results on parwith LDA/GGA for predicting accurate work functionsand interface dipoles. We note that the standard,non-self-consistent DFT+Σ yields a HOMO resonancewithin 0.1 eV difference from OT-RSH result, as expectedfor such a weakly coupled interface. We also note thatwhen the vdW-DF2 functional is used to relax the co-ordinates of the molecule and the top layer of Al(111),the benzene molecule is found at about 3.5˚A above theAl(111) surface, which results in a level alignment of 4.2eV based on the OT-RSH scheme. P DO S ( a r b . un it ) E − E F (eV) PBEOT − RSH 0 5 10 15 20 25 30 − − − − FIG. 4. PDOS of PTCDA adsorbed on Au(111), whose struc-ture is shown in the inset, as obtained using PBE (grey dashedline) and OT-RSH (red line), compared with results from UPSmeasurements (shown as green shaded area) from Ref. 121.
2. PTCDA on Au(111)
PTCDA molecules adsorbed on noble metals have beenvery well-studied (see Ref. 118 for a detailed review).It is known that the interaction between PTCDA andAu(111) is rather weak . For consistency withprevious work , our calculations use the same geome-try as in Ref. 122, where dispersion-corrected functionalsyields excellent agreement with experiment for the struc-ture, and 3 layers of Au(111) are used to represent theslab; a 2 × × × × γ is 0.16 bohr − and the OT-RSHHOMO is -8.2 eV. For Au(111), the image plane is 0.9˚A above the top surface, as determined in Ref. 70. Theaverage distance between the molecule and the top layerof the surface is 3.18 ˚A, which yields an image-chargerenormalization energy of 1.6 eV. We find that β = − . E F , which is reasonably close to the experimentalresult that is a shoulder in the UPS spectrum centeredat 1.8 eV . This is a distinct improvement over thePBE result (ca. 1.1 eV), and matches well with previousDFT+Σ calculations , as expected for weakly coupledinterfaces. OT-RSH yields a work function change of -0.7 eV, similar to that of PBE (-0.5 eV) and experiment(-0.5 eV, Ref. 121). For the value of the work function,OT-RSH yields 4.9 eV, in good agreement with PBE andexperiment (4.8 eV).
3. BDA on Au(111)
BDA adsorbed on Au(111) has attracted con-siderable attention in both experimental andtheoretical studies. However, the adsorptiongeometry has only been revealed recently with scanningtunneling microscopy measurements , which suggestthat BDA molecules form self-assembled linear chains onAu(111). Ref. 125 showed that the linear-chain phaseof BDA is energetically favored over isolated monomerphases used in previous calculations . In Ref. 125, stan-dard DFT+Σ calculations were performed to correct thePBE PDOS, without altering the PDOS lineshape. Here,we use the same geometry as in Ref. 125, but calculatethe PDOS of the interface self-consistently with our newOT-RSH approach. We use 3 layers of Au(111) as thesubstrate, and a 4 × × × × γ = 0 .
23 bohr − and HOMO at -7.0eV. As specified above, the image plane of Au(111) is at0.9 ˚A above the surface. The average distance betweenthe molecule and the surface is determined to be 3.66 ˚A.Therefore the polarization due to the substrate is 1.3 eV,based on the image-charge model. Ref. 125 showed thatfor the linear-chain structure, intermolecular polarizationis non-negligible, an effect for the HOMO resonance thatamounts to 0.3 eV. Therefore, when choosing β using Eq.(4), we use P = 1 . β = 0 .
19. We note in passing that for the purpose of tun-ing β , we directly take the intermolecular polarization of0.3 eV from Ref. 125, without attempting to determinethis value from DFT or from classical electrostatic calcu-lations (see Refs. 90,126–129 for such examples).Fig. 5 shows a comparison of our results toexperiment . As can be seen, the OT-RSH HOMO res-onance is at about 1.9 eV below E F , in much betteragreement with experimental results than PBE, whichunderestimates the level alignment by about 1 eV. Wenote that in Ref. 125, DFT+Σ was found to place theHOMO at 1.3 eV below E F , different from the OT-RSHresult in this paper. This is because, firstly, Ref. 125used an image plane position of 1.47 ˚A for Au(111), asdetermined classically from linear response , whichincreases the polarization by 0.3 eV compared to the onedetermined here. In addition, hybridization of molecu-lar orbitals with metallic states can shift the eigenval-ues when taken into account self-consistently, as com-pared to when using a post-processing approach. Impor-tantly, the DFT+Σ calculations in Ref. 125 only rigidlyshifts the PBE PDOS, leaving the lineshape unchanged,in contrast to our OT-RSH approach with its inherentself-consistency when solving the GKS equation. Ac-cordingly, in Fig. 5 one can see that OT-RSH yields aslightly more asymmetric PDOS lineshape, compared toPBE, indicating enhanced hybridization with Au d -statesat lower energy. Lastly, OT-RSH and PBE again yieldabout the same work function change (-1.5 eV) and the P DO S ( a r b . un it ) E − E F (eV) PBEOT − RSH 0 1 2 3 4 5 6 − − − − − FIG. 5. PDOS of linear-chain phase of BDA on Au(111),whose structure is shown in the inset, as obtained using PBE(grey dashed line) and OT-RSH (red line), compared withresults from UPS and x-ray photoemission spectroscopy mea-surements (green and orange shaded areas, respectively, thewidth of which represent the experimental uncertainty) fromRef. 19. same absolute value for the work function (3.8 eV).
4. AntSe on Au(111)
AntSe is known to form an upright-standing self-assembled monolayer (SAM) on Au(111) , and is thusdifferent from the other systems studied here. Adsorptionof AntSe on Au(111) can be assumed to proceed throughcleavage of hydrogen atoms and formation of a covalentAu-Se bond, i.e., it is a chemisorbed system. AntSe onAu(111) has been studied in detail, both experimentallyand theoretically, in Ref. 56; the structure we use herewas identified as the most likely one in that study, whichcompared theory results to experimental data from scan-ning tunneling microscopy. In our calculations, 4 layersof Au(111) serve as the slab, and a 8 × × × × γ = 0 .
17 bohr − and HOMO energy of -7.1 eV. As mentioned before, thelocation of the Au(111) image plane is 0.9 ˚A above thesurface, and the average distance between the moleculeand surface is 7.06 ˚A. This gives an image-charge en-ergy of 0.6 eV. Similar to the BDA/Au(111) discussedabove, due to the small intermolecular distances withinthe SAM, we expect a non-negligible polarization dueto other molecules, in addition to the polarization due tothe metal surface. As mentioned above, we do notattempt to calculate this contribution from DFT alone.To nevertheless determine this additional intermolecularpolarization energy that is responsible for gap differencebetween an isolated molecule and a SAM, we carry out GW calculations (see Ref. 133 for technical details) usingthe BerkeleyGW package for both the H-terminatedgas-phase molecule and the SAM without the metal sub-strate. We find that the gap is reduced by 1.2 eV inthe SAM compared to the gas-phase molecule. However,the convergence of the absolute value of GW quasipar-ticle energies with respect to the vacuum level is veryslow and computationally demanding. Furthermore, thedipole moment of AntSe complicates defining a uniquevacuum level for the SAM. We therefore make the as-sumption that this gap change equals 2 P , which is rea-sonable because the shapes of the HOMO and LUMOare similar. It then follows that the HOMO renormal-izes in the SAM by 0.6 eV, with respect to the gas-phasemolecule. The total renormalization of the HOMO, sum-ming contributions from the substrate and from othermolecules, is therefore 1.2 eV, which gives rise to β = 0 . for this system. OT-RSH predicts a HOMO reso-nance at 1.5 eV below E F , in much better agreement withthe experimental result (1.9 eV) than PBE, which placesthe HOMO resonance at 0.8 eV below E F . Furthermore,in the OT-RSH PDOS the consecutive features at higherbinding energies (lower in energy in Fig. 6) also matchthe UPS data very well, in sharp contrast to PBE. Thegood agreement of OT-RSH shows that the gas-phase ref-erence is physically reasonable, meaning that the HOMOorbital is nearly unaffected by the chemisorption, whichmakes sense given that it is delocalized over the backboneand not localized on the Se linker. Lastly, both OT-RSHand PBE yield excellent agreement with experiment forboth the work function change (about 1.3 eV) and itsabsolute value (3.9 eV in PBE and experiment and 4.0eV in OT-RSH). B. Systems with non-negligible charge transfer
In this section we present two systems involvingthe Ag(111) substrate with the adsorbates PTCDAand 1,4,5,8-naphthalene-tetracarboxylic dianhydride(NTCDA). Because Ag(111) has a lower work functionthan Au(111), the coupling between the moleculeand substrate is stronger and signatures in UPS thatenergetically close to the Fermi level are ascribed to the(formerly unoccupied) LUMO of the molecule. Thissignals non-negligible charge transfer between the metaland the molecule, for which standard DFT+Σ andother non-self-consistent approaches run into difficultiesbecause the weak-coupling assumption of DFT+Σ isviolated. In fact, a naive application of DFT+Σ is P DO S ( a r b . un it ) E − E F (eV) UPSPBEOT − RSH 0 10 20 30 40 50 − − − − − FIG. 6. PDOS of AntSe adsorbed on Au(111), whose struc-ture is shown in the inset, as obtained using PBE (grey dashedline) and OT-RSH (red line), compared with results from UPSmeasurements (green) taken from Ref. 56. challenging for the pinned LUMO, which is of courseunoccupied in the gas phase, but at least partiallyoccupied when the molecule adsorbs on the metalsurface. We will show in the following that our proposedOT-RSH approach also performs well for such cases.
1. PTCDA on Ag(111)
PTCDA interacts more strongly with the Ag(111)than with the Au(111) surface and, as observed inUPS and transport experiments, features a LUMO peakpinned a E F , see, e.g., Refs. 55,57,118,120,121,135–139 for detailed experimental and theoretical discus-sions. This charge transfer is already captured byLDA/GGA functionals . However, we would expectthe LDA/GGA HOMO resonance position (i.e., the sec-ond highest peak in energy in UPS) to be higher in en-ergy than the physical one, just as for the cases discussedabove. The OT-RSH method proposed in this work, aswe will show below, corrects this quantitative error forthe HOMO energy alignment, and at the same time main-tains a correct description of the LUMO pinning effect.The geometry we use is taken from Ref. 122 and has3 layers of Ag(111) as the substrate, and we employ a2 × × × × γ = 0 .
15 bohr − and a HOMO energy of-8.1 eV. Note that the gas-phase HOMO energy and op-timal γ of PTCDA is slightly different from the one re-ported above for the PTCDA/Au(111) case, as we use themolecular coordinates optimized for the interface system,which change due to the different interactions with the0 P DO S ( a r b . un it ) E − E F (eV) UPSPBEOT − RSH 0 5 10 15 20 25 30 35 40 − − − − − FIG. 7. PDOS of PTCDA adsorbed on Ag(111), whose struc-ture is shown in the inset, as obtained using PBE (grey dashedline) and OT-RSH (red line), compared with results from UPSmeasurements (green) taken from Ref. 139.
Au or Ag substrate . For Ag(111), the image plane po-sition is determined to be 1.0 ˚A above the surface. Theaverage distance between PTCDA and the Ag surface is2.87 ˚A, giving rise to an image-charge energy of 1.9 eV.When we tune β using Eq. (4) for this adsorbate, a com-plication arises: with the physically lowest possible β , asspecified above ( β = − α = − . β = − . . The OT-RSH PDOS shows the im-portant feature of LUMO pinning at E F , and on top ofthat corrects the PBE HOMO resonance energy, placingit at 1.6 eV below E F and achieving better agreementwith experiment. OT-RSH and PBE results for the workfunction change are both very small and close to zero,comparable to experimental results. Both OT-RSH andPBE yield a work function of 4.6 eV, in good agreementwith experiment. These findings are highly encouragingand show that our OT-RSH method captures all the rel-evant and highly non-trivial effects at this interface.
2. NTCDA on Ag(111)
Another interesting system that shows strong inter-action is NTCDA adsorbed on Ag(111). Similar toPTCDA on Ag(111), it also features a LUMO peak at E F , as evidenced in UPS measurements . We relaxthe NTCDA-Ag(111) geometry with the PBE+vdW surf method as implemented in VASP, using 3 layers of P DO S ( a r b . un it ) E − E F (eV) UPSPBEOT − RSH 0 5 10 15 20 25 30 − − − − FIG. 8. PDOS of NTCDA adsorbed on Ag(111), whose struc-ture is shown in the inset, as obtained using PBE (grey dashedline) and OT-RSH (red line), compared with results from UPSmeasurements (green) taken from Ref. 141.
Ag(111) and a 4 × × × × × γ for the gas-phase molecule is 0.19 bohr − , according to Eq. (3), andits HOMO energy is -9.7 eV. For Ag(111), the image-plane position is at 1.0 ˚A above the surface, as deter-mined above. The average distance between the opti-mized molecule and the surface is 2.88 ˚A, yielding animage-charge energy of 1.9 eV. This requires β = − . E F . PBE correctly shows aLUMO pinned at E F , describes the HOMO reasonablywell, but misplaces the HOMO-1 resonance (by 0.7 eV)with respect to results from OT-RSH and experiment.This strongly suggests that OT-RSH is not only accu-rate for the frontier resonances, but actually maintainspredictive power within several eVs from E F , ow-ing to the nonlocal exchange operator . OT-RSH andPBE yield similar work function values (4.7 eV and 4.8eV, respectively). Regarding the work function change,OT-RSH and PBE differ slightly more (by 0.4 eV andhave different signs) in this case. However, we are notaware of experimental data for the work function changein this case, and our results suggest this would be a mea-surement of interest. IV. DISCUSSION
In this paper, we have proposed an approach for ac-curate calculation of level alignments at molecule-metal1
TABLE I. Summary of the OT-RSH and PBE results and experimental literature data for the interface systems studied in thiswork. E F − E HOMO is the level alignment for the HOMO resonance, and ∆Φ is the work function change compared to the cleanmetal surface. Φ is the resulting work function of the interface. Where we are not aware of an experimental ∆Φ or Φ value,an “x” sign is used. γ denotes the optimally-tuned range-separation parameter based on Eq. (3) and is in bohr − . ε OT − RSHHOMO isthe gas-phase OT-RSH HOMO energy. z is the DFT-determined image-plane position for the metal surface. P is the surfacepolarization calculated using an image-charge model. For systems with small surface unit cells, such as BDA/Au(111) andAntSe/Au(111), P also includes polarization contribution due to other molecules in the molecular layer. β is tuned according toEq. (4) and is used for the calculation of the interface. The first four systems show negligible charge transfer and are discussedin Sec. III A. The last two systems show significant charge transfer and are discussed in Sec. III B. E F − E HOMO (eV) ∆Φ [Φ] (eV) γ ( a − ) ε OT − RSHHOMO (eV) z (˚A) P (eV) β commentPBE OT-RSH Expt. PBE OT-RSH Expt.Benzene/Al(111) 3.1 4.4 4.0 -0.3 [3.8] -0.3 [3.8] -0.2 [x] -0.5 [4.8] -0.7 [4.9] -0.5 [4.8] -1.5 [3.8] -1.6 [3.8] x [x] 0.23 -7.0 0.9 1.6 0.19 stronger hybridizationAntSe/Au(111) 0.8 1.5 1.9 -1.4 [3.9] -1.5 [4.0] -1.3 [3.9] +0.1 [4.6] 0.0 [4.6] -0.1 [4.8] ,0.1 [4.9] +0.2 [4.8] -0.2 [4.7] x [x] 0.19 -9.7 1.0 1.9 -0.082 Fermi level pinning interfaces based on the OT-RSH functional, and testedit for a set of well-studied and complex systems. Ourresults, summarized in Table I, show that this methodis successful in quantitative predictions of energy levelalignment. To understand its success, we refer again toFig. 2 and compare it to Fig. 1: The essential idea of ourOT-RSH approach is to tune the parameter β and adaptthe long-range Fock exchange to capture the screeningdue to the metal surface. The amount of screeningand the choice of β are determined by P , given in Eq.(4). In Fig. 2, the red lines on the right demonstratethe attempted gas-phase HOMO/LUMO energy with thetuned β value, and the red curve on the left shows theresulting PDOS of the molecule in the interface. In Eq.(4), for the systems studied here the left hand side is al-most linear in β , which means that the tuned β value, aswell as the resulting level alignment at the interface, arenot sensitive to small variations of P .With our method we set out to capture physics inher-ent to GW calculations but absent from currently avail-able DFT approaches for interfaces, namely a sensitivityof XC effects to the dielectric environment specific tomolecule-metal interfaces. The uneven performance ofstandard DFT approaches for interfaces, i.e., the successin describing work functions and interface dipoles butthe failures in predicting level alignments, can be un-derstood from the “nearsightedness principle” of many-electron systems . It states that many important staticphysical quantities, including the density, only dependon changes in the potential at nearby points. This ex-plains the success of local and semi-local XC function-als in predicting the density of molecule-metal systemsand highly relevant quantities that are directly deter-mined from it, notably total energies, work functions,and interface dipoles. However, long-range Coulomb in-teractions are an exception to this principle. They areimportant for molecule-metal interfaces, especially con-sidering the renormalization of molecular resonances dueto the metal surface. Because of the long-range nature ofthe phenomenon, a standard hybrid functional without range separation, or a range-separated hybrid functionalwithout long-range Fock exchange, cannot in general cap-ture the image-charge effect accurately. Therefore, wechoose to adapt the long-range part of the XC poten-tial by tuning β to capture P , thereby implicitly includ-ing the otherwise missing distance-dependent screeningin our approach, which is crucial for its success. At thesame time, we deliberately do not modify the short-rangepart of the XC potential (controlled by α and γ ) whengoing from the gas phase to the surface. Importantly,orbital energies and resonances are expected to respondto β tuning, but the density and quantities that are di-rectly determined from it are expected to be mostly un-changed. This reasoning is strongly supported by ourresults of level alignments and work function changes.Our scheme therefore provides a superior GKS approachthat can predict both energy level alignments and workfunction changes with high accuracy.As we have stressed throughout this paper, one advan-tage of our fully self-consistent approach, as compared toDFT+Σ and other non-self-consistent schemes, is that itcalculates and fully captures the hybridization betweenthe molecule and the metal. Therefore, it can alter thePBE PDOS lineshape and is also applicable to systemswith stronger charge transfer and Fermi level pinning. Inaddition, since it is based on GKS and includes a frac-tion of short-range Fock exchange, it can also yield im-proved relative orbital spacing, as shown in the examplesof AntSe/Au(111) and NTCDA/Ag(111), as well as im-proved orbital ordering, as was demonstrated previouslyfor other systems .In spite of its success, we also would like to point outthat the method proposed here is not a panacea. Inparticular, we focused here on the level alignment andPDOS of the molecule, and did not discuss how the spe-cific choice of α, β, and γ may affect specifics in the elec-tronic structure of the metal substrate. In this context,we would like to mention that full Fock exchange often ar-tificially opens a gap for metals , and that many hybridfunctionals perform worse than local and semi-local func-2tionals for metals . We expect that similar situationsmay occur when one has to choose a large β value follow-ing the tuning procedure. Coincidentally, this is not thecase for the systems studied in this work. The largest β used in this work is about 0.2, and with this amount oflong-range Fock exchange the PDOS of the metal sub-strate stays qualitatively correct. Quantitatively, evenfor a rather small amount of Fock exchange, the theoret-ical description of certain properties of real metals maysuffer .Another limitation is that for PTCDA/Ag(111), andperhaps also other systems, with the physically lowestpossible value for β , i.e., β = − α = − . P by tuning β . For PTCDA/Ag(111), the re-maining difference was small. However, one may imag-ine molecule-metal systems where the deficit is somewhatlarger, and using β = − . . In those cases, the currently proposed OT-RSHapproach may not work as well as for the systems shownin this paper.Lastly, we have used a DFT-based image-charge modelto approximate the otherwise rather complicated changein the Coulomb screening induced by the metal sur-face. As far as this particular approximation is con-cerned, the approach proposed in this work shares thesame advantages and disadvantages as the “standard”DFT+Σ method: The image-charge interaction mimicsstatic (frequency-independent) polarization, and missesdynamical (frequency-dependent) surface polarization ef-fects and the polarization of the molecule due to themetal, which is a higher order effect. Moreover, for casesof non-negligible intermolecular polarization such as po-larizable SAMs, so far we did not attempt to computethis solely from DFT, but for this proof-of-principle studyrelied on GW results.One way of understanding the above limitations is torealize that molecule-metal interfaces are indeed stronglyheterogeneous systems in the sense that the screen-ing perpendicular to the surface is strongly distance-dependent and very different from the screening paral-lel to it. In the OT-RSH approach proposed here, weuse scalar parameters, α, β, and γ , which may not beenough for complicated cases. Work along the lines of“local hybrids” , i.e., spatially dependent parameters α ( r ) , β ( r ) , or even γ ( r ), as well as density-based mixingparameters , may describe the heterogeneous environ-ment more accurately. But the construction of such func-tionals and their implementations may prove to be muchmore difficult. V. CONCLUSION
In this work, we developed an approach to calculatethe energy level alignment at molecule-metal interfaceswith good accuracy, based on an OT-RSH functional.The essential idea is to start with an accurate electronicstructure description in the gas phase and capture thenonlocal surface polarization due to the metal by screen-ing the long-range Fock exchange. We proposed a non-empirical way to tune the long-range Fock exchange permolecule-metal pair based on an image-charge model,and implemented this approach in a plane-wave code.Results from our fully self-consistent approach for sev-eral prototypical, challenging molecule-metal interfacesare in quantitative agreement with experiments, for boththe level alignments and work function changes. Keep-ing its remaining limitations as discussed in this paper inmind, we believe that our OT-RSH approach paves theway for accurate and reliable predictions of energeticsand level alignments at heterogeneous interfaces, espe-cially those related to energy conversion and molecularelectronics. Finally, we note that the efficiency of thesecalculations strongly depends on new algorithm develop-ments for computing nonlocal Fock exchange, such as therecently proposed adaptively compressed exchange oper-ator method . VI. ACKNOWLEDGEMENT
We thank Achim Sch¨oll for providing us the experi-mental data for NTCDA adsorbed on Ag(111) that ispublished in Ref. 141. Work in Berkeley was supportedby the U.S. Department of Energy, Office of Basic EnergySciences, Division of Materials Sciences and Engineeringunder Contract No. DE-AC02-05CH11231. Work per-formed at the Molecular Foundry was also supported bythe Office of Science, Office of Basic Energy Sciences,of the U.S. Department of Energy under the same con-tract number. Work in Rehovoth was supported by theEuropean Research Council, the Israel Science Founda-tion, the United States-Israel Binational Science Foun-dation, the Wolfson Foundation, the Austrian ScienceFund (FWF):J3608-N20, and the Molecular Foundry. Wethank the National Energy Research Scientific Comput-ing center for computational resources. H. Ishii, K. Sugiyama, E. Ito, and K. Seki, Adv. Mater. , 605(1999). A. Kahn, N. Koch, and W. Gao, J. Polym. Sci., Part B: Polym.Phys. , 2529 (2003). N. Koch, Chem. Phys. Chem. , 1438 (2007). N. Ueno and S. Kera, Prog. Surf. Sci. , 490 (2008). J. Hwang, A. Wan, and A. Kahn, Mater. Sci. Eng., R , 1(2009). S. Braun, W. R. Salaneck, and M. Fahlman, Adv. Mater. ,1450 (2009). N. Koch, N. Ueno, and A. T. S. Wee, eds.,
The Molecule-MetalInterface (Wiley-VCH, 2013). W. Liu, A. Tkatchenko, and M. Scheffler, Acc. Chem. Res. ,3369 (2014). M. Willenbockel, D. Luftner, B. Stadtmuller, G. Koller,C. Kumpf, S. Soubatch, P. Puschnig, M. G. Ramsey, and F. S.Tautz, Phys. Chem. Chem. Phys. , 1530 (2015). R. J. Maurer, V. G. Ruiz, J. Camarillo-Cisneros, W. Liu,N. Ferri, K. Reuter, and A. Tkatchenko, Prog. Surf. Sci. , 72(2016). N. D. Lang, Phys. Rev. B , 4234 (1971). S. Y. Quek, L. Venkataraman, H. J. Choi, S. G. Louie, M. S.Hybertsen, and J. B. Neaton, Nano Lett. , 3477 (2007). Y. Xue, S. Datta, and M. A. Ratner, J. Chem. Phys. , 4292(2001). M. A. Reed, C. Zhou, C. J. Muller, T. P. Burgin, and J. M.Tour, Science , 252 (1997). N. J. Tao, Nature Nanotech. , 173 (2006). L. Venkataraman, J. E. Klare, C. Nuckolls, M. S. Hybertsen,and M. L. Steigerwald, Nature , 904 (2006). A. C. D¨urr, N. Koch, M. Kelsch, A. R¨uhm, J. Ghijsen, R. L.Johnson, J.-J. Pireaux, J. Schwartz, F. Schreiber, H. Dosch,and A. Kahn, Phys. Rev. B , 115428 (2003). S. Y. Quek, M. Kamenetska, M. L. Steigerwald, H. J. Choi, S. G.Louie, M. S. Hybertsen, J. B. Neaton, and L. Venkataraman,Nature Nanotechnology , 230 (2009). M. Dell’Angela, G. Kladnik, A. Cossaro, A. Verdini,M. Kamenetska, I. Tamblyn, S. Y. Quek, J. B. Neaton,D. Cvetko, A. Morgante, and L. Venkataraman, Nano Lett. , 2470 (2010). L. Hedin, Phys. Rev. , A796 (1965). M. S. Hybertsen and S. G. Louie, Phys. Rev. B , 5390 (1986). X. Blase, C. Attaccalite, and V. Olevano, Phys. Rev. B ,115103 (2011). S. Sharifzadeh, I. Tamblyn, P. Doak, P. T. Darancet, and J. B.Neaton, Eur. Phys. J. B , 323 (2012). C. Faber, C. Attaccalite, V. Olevano, E. Runge, and X. Blase,Phys. Rev. B , 115123 (2011). C. Faber, P. Boulanger, C. Attaccalite, I. Duchemin, andX. Blase, Philos. Trans. R. Soc. Lond. A–Math. Phys. Eng. Sci. (2014). M. J. van Setten, F. Caruso, S. Sharifzadeh, X. Ren, M. Schef-fler, F. Liu, J. Lischner, L. Lin, J. R. Deslippe, S. G. Louie,C. Yang, F. Weigend, J. B. Neaton, F. Evers, and P. Rinke, J.Chem. Theory Comp. , 5665 (2015). F. Kaplan, M. E. Harding, C. Seiler, F. Weigend, F. Evers, andM. J. van Setten, J. Chem. Theory Comput. , 2528 (2016). M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett. , 1418(1985). W. G. Aulbur, L. J¨onsson, and J. W. Wilkins, in
Solid StatePhysics , Vol. 54, edited by H. Ehrenreich and S. Spaepen (Aca-demic Press, 2000) pp. 1 – 218. P. Rinke, A. Qteish, J. Neugebauer, C. Freysoldt, and M. Schef-fler, New J. Phys. , 126 (2005). Y. Li, D. Lu, S. A. Swanson, J. C. Scott, and G. Galli, J. Phys.Chem. C , 6413 (2008). N. Kharche, J. T. Muckerman, and M. S. Hybertsen, Phys.Rev. Lett. , 176802 (2014). C. Draxl, D. Nabok, and K. Hannewald, Acc. Chem. Res. ,3225 (2014). M. Strange and K. S. Thygesen, Phys. Rev. B , 195121 (2012). M. van Schilfgaarde, T. Kotani, and S. V. Faleev, Phys. Rev.B , 245125 (2006). I. Tamblyn, P. Darancet, S. Y. Quek, S. Bonev, and J. B.Neaton, Phys. Rev. B , 201402 (2011). P. Rinke, A. Qteish, J. Neugebauer, and M. Scheffler, Phys.Status Solidi B , 929 (2008). F. Bruneval and M. A. L. Marques, J. Chem. Theory Comput. , 324 (2013). P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). D. P. Chong, O. V. Gritsenko, and E. J. Baerends, J. Chem.Phys. , 1760 (2002). C.-O. Almbladh and U. von Barth, Phys. Rev. B , 3231(1985). L. Kronik, T. Stein, S. Refaely-Abramson, and R. Baer, J.Chem. Theory Comput. , 1515 (2012). L. Kronik and S. K¨ummel, Top. Curr. Chem. , 137 (2014). J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys.Rev. Lett. , 1691 (1982). J. P. Perdew and M. Levy, Phys. Rev. Lett. , 1884 (1983). L. J. Sham and M. Schl¨uter, Phys. Rev. Lett. , 1888 (1983). S. K¨ummel and L. Kronik, Rev. Mod. Phys. , 3 (2008). J. C. Inkson, J. Phys. C , 1350 (1973). J. B. Neaton, M. S. Hybertsen, and S. G. Louie, Phys. Rev.Lett. , 216405 (2006). J. M. Garcia-Lastra, C. Rostgaard, A. Rubio, and K. S. Thyge-sen, Phys. Rev. B , 245427 (2009). C. Freysoldt, P. Rinke, and M. Scheffler, Phys. Rev. Lett. ,056803 (2009). V. De Renzi, R. Rousseau, D. Marchetto, R. Biagi, S. Scandolo,and U. del Pennino, Phys. Rev. Lett. , 046804 (2005). I. Magid, L. Burstein, O. Seitz, L. Segev, L. Kronik, andY. Rosenwaks, J. Phys. Chem. C , 7145 (2008). L. Romaner, D. Nabok, P. Puschnig, E. Zojer, andC. Ambrosch-Draxl, New J. Phys. , 053010 (2009). A. M. Track, F. Rissner, G. Heimel, L. Romaner, D. K¨afer,A. Bashir, G. M. Rangger, O. T. Hofmann, T. Bu˘cko, G. Witte,and E. Zojer, J. Phys. Chem. C , 2677 (2010). O. T. Hofmann, V. Atalla, N. Moll, P. Rinke, and M. Scheffler,New J. Phys. , 123028 (2013). Y. Huang, E. Wruss, D. Egger, S. Kera, N. Ueno, W. Saidi,T. Bucko, A. Wee, and E. Zojer, Molecules (2014). T. Abu-Husein, S. Schuster, D. A. Egger, M. Kind, T. San-towski, A. Wiesner, R. Chiechi, E. Zojer, A. Terfort, andM. Zharnikov, Adv. Func. Mater. , 3943 (2015). F. Flores, J. Ortega, and H. Vazquez, Phys. Chem. Chem. Phys. , 8658 (2009). M. Rohlfing, Phys. Rev. B , 205127 (2010). A. Greuling, M. Rohlfing, R. Temirov, F. S. Tautz, and F. B.Anders, Phys. Rev. B , 125413 (2011). A. M. Souza, I. Rungger, C. D. Pemmaraju, U. Schwingen-schloegl, and S. Sanvito, Phys. Rev. B , 165112 (2013). M. Yu, P. Doak, I. Tamblyn, and J. B. Neaton, J. Phys. Chem.Lett. , 1701 (2013). A. Migani, D. J. Mowbray, J. Zhao, H. Petek, and A. Rubio,J. Chem. Theory Comput. , 2103 (2014). T. Esat, T. Deilmann, B. Lechtenberg, C. Wagner, P. Kr¨uger,R. Temirov, F. B. Anders, M. Rohlfing, and F. S. Tautz, Phys.Rev. B , 144415 (2015). J. E. Moore and L. Jensen, J. Phys. Chem. C , 5659 (2016). S. Roychoudhury, C. Motta, and S. Sanvito, Phys. Rev. B ,045130 (2016). J. Ma, Z.-F. Liu, J. B. Neaton, and L.-W. Wang, Appl. Phys.Lett. , 262104 (2016). D. A. Egger, Z.-F. Liu, J. B. Neaton, and L. Kronik, Nano Lett. , 2448 (2015). H. V´azquez, R. Oszwaldowski, P. Pou, J. Ortega, R. P´erez,F. Flores, and A. Kahn, Eur. Phys. Lett. , 802 (2004). G. Heimel, S. Duhm, I. Salzmann, A. Gerlach, A. Strozecka,J. Niederhausen, C. B¨urker, T. Hosokai, I. Fernandez-Torrente,G. Schulze, S. Winkler, A. Wilke, R. Schlesinger, J. Frisch,B. Br¨oker, A. Vollmer, B. Detlefs, J. Pflaum, S. Kera, K. J.Franke, N. Ueno, J. I. Pascual, F. Schreiber, and N. Koch, Nat.Chem. , 187 (2013). G. Heimel, L. Romaner, J.-L. Br´edas, and E. Zojer, Phys. Rev.Lett. , 196806 (2006). A. d. M. Souza, I. Rungger, R. B. Pontes, A. R. Rocha,A. J. Roque da Silva, U. Schwingenschloegl, and S. Sanvito,Nanoscale , 14495 (2014). T. Rangel, G.-M. Rignanese, and V. Olevano, Beilstein J. Nan-otechnol. , 1247 (2015). A. Seidl, A. G¨orling, P. Vogl, J. A. Majewski, and M. Levy,Phys. Rev. B , 3764 (1996). A. D. Becke, J. Chem. Phys. , 5648 (1993). J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin, J.Chem. Phys. , 174101 (2005). A. Biller, I. Tamblyn, J. B. Neaton, and L. Kronik, J. Chem.Phys. , 164706 (2011). R. Baer, E. Livshits, and U. Salzner, Annu. Rev. Phys. Chem. , 85 (2010). T. Stein, H. Eisenberg, L. Kronik, and R. Baer, Phys. Rev.Lett. , 266802 (2010). S. Refaely-Abramson, R. Baer, and L. Kronik, Phys. Rev. B , 075144 (2011). T. K¨orzd¨orfer, J. S. Sears, C. Sutton, and J.-L. Br´edas, J.Chem. Phys. , 204107 (2011). S. Refaely-Abramson, S. Sharifzadeh, N. Govind, J. Autschbach,J. B. Neaton, R. Baer, and L. Kronik, Phys. Rev. Lett. ,226405 (2012). D. A. Egger, S. Weissman, S. Refaely-Abramson, S. Sharifzadeh,M. Dauth, R. Baer, S. K¨ummel, J. B. Neaton, E. Zojer, andL. Kronik, J. Chem. Theory Comp. , 1934 (2014). H. Phillips, Z. Zheng, E. Geva, and B. D. Dunietz, Org. Electr. , 1509 (2014). I. Tamblyn, S. Refaely-Abramson, J. B. Neaton, and L. Kronik,J. Phys. Chem. Lett. , 2734 (2014). T. K¨orzd¨orfer and J.-L. Br´edas, Acc. Chem. Res. , 3284(2014). J. Autschbach and M. Srebro, Acc. Chem. Res. , 2592 (2014). S. Refaely-Abramson, S. Sharifzadeh, M. Jain, R. Baer, J. B.Neaton, and L. Kronik, Phys. Rev. B , 081204(R) (2013). L. Kronik and J. B. Neaton, Annu. Rev. Phys. Chem. , 587(2016). Z.-F. Liu, S. Wei, H. Yoon, O. Adak, I. Ponce, Y. Jiang, W.-D. Jang, L. M. Campos, L. Venkataraman, and J. B. Neaton,Nano Lett. , 5365 (2014). A. Yamada, Q. Feng, A. Hoskins, K. D. Fenk, and B. D. Duni-etz, Nano Lett. , 6092 (2016). J. Toulouse, F. Colonna, and A. Savin, Phys. Rev. A , 062505(2004). T. Yanai, D. P. Tew, and N. C. Handy, Chem. Phys. Lett. ,51 (2004). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996); ibid. , , 1396(E) (1997). T. M. Henderson, B. G. Janesko, and G. E. Scuseria, J. Chem.Phys. , 194105 (2008). D. L¨uftner, S. Refaely-Abramson, M. Pachler, R. Resel, M. G.Ramsey, L. Kronik, and P. Puschnig, Phys. Rev. B , 075204(2014). S. Refaely-Abramson, M. Jain, S. Sharifzadeh, J. B. Neaton,and L. Kronik, Phys. Rev. B , 081204 (2015). J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. , 8207 (2003); ibid. , , 219906(E) (2006). A. G. Eguiluz and W. Hanke, Phys. Rev. B , 10433(R) (1989). Y. Li, D. Lu, and G. Galli, J. Chem. Theory Comput. , 881(2009). P. Giannozzi et al. , J. Phys.: Condens. Matter , 395502(2009). F. Gygi and A. Baldereschi, Phys. Rev. B , 4405(R) (1986). A. Marini, G. Onida, and R. Del Sole, Phys. Rev. Lett. ,016403 (2002). S. Sharma, J. K. Dewhurst, and C. Ambrosch-Draxl, Phys.Rev. Lett. , 136402 (2005). We create norm-conserving Ag and Au pseudopotentials usingthe LD1 utility of Quantum ESPRESSO. For Ag with semi-core states, the r c values for 4s, 4p, and 4d orbitals are 1.0,1.8, and 1.2 a.u., respectively. For Ag without semicore states,the r c for 4d, 5s, and 5p orbitals are all 2.4 a.u. For Au withsemicore states, the r c for 5s, 5p, and 5d orbitals are 1.03, 2.3, and 2.3 a.u., respectively. For Au without semicore states, the r c We have checked the effect of mixing these two types of pseu-dopotentials on the level alignment, using the system of 2,3,5,6-tetrafluoro-7,7,8,8-tetracyano-quinodimethane (F4TCNQ) ad-sorbed on Ag(111) surface. The calculation is performed usingHSE functional and the resulting level alignment is in quanti-tative agreement with both PAW calculations in VASP and thepublished results in Ref. 57 using the same functional.
G. Kresse and J. Hafner, Phys. Rev. B , 558 (1993); ibid. , ,14251 (1994); G. Kresse and J. Furthm¨uller, Comput. Mat. Sci. , 15 (1996); Phys. Rev. B , 11169 (1996). P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994). G. Kresse and D. Joubert, Phys. Rev. B , 1758 (1999). Y. Shao et al. , Mol. Phys. , 184 (2014).
M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P.Straatsma, H. J. J. van Dam, D. Wang, J. Nieplocha, E. Apra,T. L. Windus, and W. A. de Jong, Comput. Phys. Commun. , 1477 (2010).
R. Duschek, F. Mittendorfer, R. I. R. Blyth, F. P. Netzer,J. Hafner, and M. G. Ramsey, Chem. Phys. Lett. , 43(2000).
A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. , 073005(2009).
A. Tkatchenko, R. A. DiStasio, Jr., R. Car, and M. Scheffler,Phys. Rev. Lett. , 236402 (2012).
K. Lee, ´E. D. Murray, L. Kong, B. I. Lundqvist, and D. C.Langreth, Phys. Rev. B , 081101(R) (2010). F. Tautz, Prog. Surf. Sci , 479 (2007). H. V´azquez, R. Oszwaldowski, P. Pou, J. Ortega, R. P´erez,F. Flores, and A. Kahn, Europhys. Lett. , 802 (2004). S. K. M. Henze, O. Bauer, T.-L. Lee, M. Sokolowski, andF. Tautz, Surf. Sci. , 1566 (2007).
S. Duhm, A. Gerlach, I. Salzmann, B. Br¨oker, R. L. Johnson,F. Schreiber, and N. Koch, Org. Electron. , 111 (2008). V. G. Ruiz, W. Liu, E. Zojer, M. Scheffler, and A. Tkatchenko,Phys. Rev. Lett. , 146103 (2012).
T. K. Haxton, H. Zhou, I. Tamblyn, D. Eom, Z. Hu, J. B.Neaton, T. F. Heinz, and S. Whitelam, Phys. Rev. Lett. ,265701 (2013).
G. Li, I. Tamblyn, V. R. Cooper, H.-J. Gao, and J. B. Neaton,Phys. Rev. B , 121409(R) (2012). G. Li, T. Rangel, Z.-F. Liu, V. R. Cooper, and J. B. Neaton,Phys. Rev. B , 125429 (2016). L. Romaner, G. Heimel, C. Ambrosch-Draxl, and E. Zojer, Adv.Func. Mater. , 3999 (2008). A. Natan, N. Kuritz, and L. Kronik, Adv. Func. Mater. ,2077 (2010). S. Sharifzadeh, A. Biller, L. Kronik, and J. B. Neaton, Phys.Rev. B , 125307 (2012). D. Vanzo, B. J. Topham, and Z. G. Soos, Adv. Func. Mater. , 2004 (2015). S. C. Lam and R. J. Needs, J. Phys.: Condens. Matter , 2101(1993). A. Bashir, D. K¨afer, J. M¨uller, C. W¨oll, A. Terfort, andG. Witte, Angew. Chem., Int. Ed. , 5250 (2008). N. Sato, K. Seki, and H. Inokuchi, J. Chem. Soc., FaradayTrans. 2 , 1621 (1981). In the GW calculations, 1500 bands and an energy cutoff of 10Ry are used to calculate the polarizability. The dielectric func-tion is calculated using the generalized plasmon pole model.A box truncation and a single q-point (Γ) are used for themolecule. A slab truncation and a 8 × × J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L.Cohen, and S. G. Louie, Computer Phys. Commun. , 1269(2012).
Y. Zou, L. Kilian, A. Sch¨oll, T. Schmidt, R. Fink, and E. Um-bach, Surf. Sci. , 1240 (2006).
A. Kraft, R. Temirov, S. K. M. Henze, S. Soubatch, M. Rohlfing,and F. S. Tautz, Phys. Rev. B , 041402 (2006). M. Rohlfing, R. Temirov, and F. S. Tautz, Phys. Rev. B ,115421 (2007). J. Ziroff, F. Forster, A. Sch¨oll, P. Puschnig, and F. Reinert,Phys. Rev. Lett. , 233004 (2010).
B. Stadtm¨uller, D. L¨uftner, M. Willenbockel, E. M. Reinisch,T. Sueyoshi, G. Koller, S. Soubatch, M. G. Ramsey, P. Puschnig,F. S. Tautz, and C. Kumpf, Nat. Commun. , 3685 (2014). A. Bendounan, F. Forster, A. Sch¨oll, D. Batchelor, J. Ziroff,E. Umbach, and F. Reinert, Surf. Sci. , 4013 (2007).
A. Sch¨oll, L. Kilian, Y. Zou, J. Ziroff, S. Hame, F. Reinert,E. Umbach, and R. H. Fink, Science , 303 (2010).
J. Ziroff, S. Hame, M. Kochler, A. Bendounan, A. Sch¨oll, andF. Reinert, Phys. Rev. B , 161404(R) (2012). T. K¨orzd¨orfer and S. K¨ummel, Phys. Rev. B , 155206 (2010). W. Kohn, Phys. Rev. Lett. , 3168 (1996). N. Dori, M. Menon, L. Kilian, M. Sokolowski, L. Kronik, andE. Umbach, Phys. Rev. B , 195208 (2006). F. Rissner, D. A. Egger, A. Natan, T. K¨orzd¨orfer, S. K¨ummel,L. Kronik, and E. Zojer, J. Am. Chem. Soc. , 18634 (2011).
P. Puschnig, E.-M. Reinisch, T. Ules, G. Koller, S. Soubatch,M. Ostler, L. Romaner, F. S. Tautz, C. Ambrosch-Draxl, andM. G. Ramsey, Phys. Rev. B , 235427 (2011). G. D. Mahan,
Many-Particle Physics , 3rd ed. (Springer, 2000).
J. Paier, M. Marsman, and G. Kresse, J. Chem. Phys. ,024103 (2007).
T. Ules, D. L¨uftner, E. M. Reinisch, G. Koller, P. Puschnig,and M. G. Ramsey, Phys. Rev. B , 155430 (2014). J. Jaramillo, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. , 1068 (2003).
M. A. L. Marques, J. Vidal, M. J. T. Oliveira, L. Reining, andS. Botti, Phys. Rev. B , 035119 (2011). L. Lin, J. Chem. Theory Comput.12