{'fraction_non_alphanumeric': 0.04935634705305086, 'fraction_numerical': 0.018299370778733476, 'mean_word_length': 5.038340885325898, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantitative Current-Voltage Characteristics in Molecular Junctions from First Principles

Using self-energy-corrected density functional theory (DFT) and a coherent scattering-state approach, we explain current-voltage (IV) measurements of four pyridine-Au and amine-Au linked molecular junctions with quantitative accuracy. Parameter-free many-electron self-energy corrections to DFT Kohn-Sham eigenvalues are demonstrated to lead to excellent agreement with experiments at finite bias, improving upon order-of-magnitude errors in currents obtained with standard DFT approaches. We further propose an approximate route for prediction of quantitative IV characteristics for both symmetric and asymmetric molecular junctions based on linear response theory and knowledge of the Stark shifts of junction resonance energies. Our work demonstrates that a quantitative, computationally inexpensive description of coherent transport in molecular junctions is readily achievable, enabling new understanding and control of charge transport properties of molecular-scale interfaces at large bias voltages.

Given the diversity of known synthesizable organic molecules, and a lack of intuition connecting transport phenomena to a specific molecule and interface, a rigorous yet pragmatic quantitative theory capable of predicting charge transport properties of molecular junctions with chemical specificity is needed. Charge transport calculations require solution of the electronic structure for an open system out of equilibrium [24][25][26][27][28][29]. The majority of prior theoretical studies have relied on a Landauer approach, simplified to treat electronic interactions at a mean-field level within density functional theory (DFT) using either Green's functions [30][31][32][33] or scattering-states [34]. However, within common approximations to DFT, such as generalized gradient approximations (GGAs) and hybrid functionals, Kohn-Sham orbital energies are known to yield unsatisfactory level alignment [35][36][37][38][39][40], leading to overestimated conductance [24][25][26][27][28][29] and even incorrect trends [41] relative to experiment. Improved treatment of exchange and correlation effects within the junction [29,37,38,[41][42][43][44][45][46][47][48][49][50] has been shown to lead to better agreement with experiment. In particular, recent calculations based on many-body perturbation theory within the GW approximation [15,[42][43][44]51] have been demonstrated to amend junction level alignment and have resulted in conductance in quantitative agreement with measured values [15,[42][43][44]. However, these approaches rely on computationally expensive methods [37,38,[45][46][47][48][49][50], and their accuracy come at the cost of their tractability for junctions relevant to experiments; moreover, direct comparisons to experiment have been limited to low-bias measurements in the linear response regime. The ability to measure, reliably, conductance and nonlinear IV characteristics, and the promise of exploring new phenomena at higher bias, calls for quantitative computational studies of truly nonequilibrium steady-states.
In this Letter, we demonstrate and apply a quantitative framework capable of explaining measured IV characteristics for four molecules -4,4' diaminostilbene (DAS), bis-(4-aminophenyl)acetylene (APA), 1,6-hexanediamine (HDA), and 4,4'-bipyridine (BP)-in contact with gold electrodes for biases as high as 1V. Whereas the magnitude of the currents and shape of the measured IV characteristic are not reproduced by a standard finite-bias DFT approach [30][31][32][33][34], a method based on the GW approximation [51], successful for low-bias conductance and extended here to finite bias, leads to excellent agreement with experiment. Our approach enables inexpensive and accurate calculations of coherent transport properties, which can be used for the design of functional molecular junctions with nonlinear IV characteristics.
Previous studies [15,[42][43][44] have shown amine-and pyridine-bonded molecules can bind preferentially to undercoordinated gold atoms. Building on those studies, we construct three geometries for each junction with trimer (3 gold atoms), trimer-adatom, and adatom binding motifs (see SI for geometries). We relax all junctions using DFT within the GGA of Perdew, Burke, and Ernzenhof (PBE) [52] and a double-ζ-basis set as implemented in SIESTA [53,54]. Details of our DFT calculations are provided in previous work [42][43][44]. Atomic positions are relaxed until Hellmann-Feynman forces are smaller than 0.04 eV/Å. We model our system by 7-layers of 16 gold atoms on both sides, the last 4 layers being constrained to the bulk geometry. Initial trial geometries are adapted from previous works [15, 42-44, 55, 56]. We consider BP junctions in a "high conductance" configuration (HighG), following the notation of Ref. [15].
In what follows, the conductance and IV characteristics are computed from a Landauer-like formula using a coherent scattering-state approach, and include exchange and correlation contributions that correct for zero-bias junction level alignment, as explained below. In this framework, the IV characteristic is expressed in terms of an energy-and bias-dependent transmission function T (ω; V ) as where f is the lead Fermi-Dirac occupation function and ω has units of energy. The transmission function is computed as T (ω; V ) = tr t (ω; V ) t † (ω; V ) , where t and t † are the transmission coefficients of the junction scattering states computed as a function of bias [34,57]. The zero-bias conductance is determined via linear response, , where E F is the junction Fermi level at zero bias.
We obtain the self-consistent steady-state density matrix from DFT-PBE as described in [34,57] using an 8×8 k // -mesh and in a manner equivalent to prior work [30][31][32][33][34]. For finite bias calculations, referred to as DFT(V), the density matrix includes a real-axis integration of the scattering-states in the bias window on an adaptive energy grid [34] with a resolution of up to 10 −7 eV in the vicinity of the molecular resonances. The transmission function is generated in a subsequent step using a 16 × 16k // -mesh.
To correct for inaccuracies associated with DFT-PBE Kohn-Sham eigenvalues for quasiparticle energy level alignment, we employ a physically-motivated electron self-energy correction to the molecular orbital energies in the junction, DFT+Σ, following Ref. [42][43][44]. Formally derivable in the weak coupling limit of the GW approximation [44,58,59], this adjustable-parameter-free model self-energy acts on the molecular subspace and consists of two terms: i) a gas-phase correction accounting for the difference between DFT highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) energies and their gas-phase ionization potential (IP) and electron affinity (EA) (Σ GP ); and ii) an "image charge" term, accounting for the polarization energy associated with static non-local correlations between the electrons (or holes) on the molecule and in the metal that close the gap of the molecule upon absorption (Σ CORR ). In this work, we calculate Σ GP using a ∆SCF method, and following previous studies [35,[42][43][44]51], we approximate Σ CORR with an electrostatic image charge model (see SI for details). Since these self-energy corrections are large compared to those for the metallic bulk and surface Au states [60], especially for states near E F , we neglect similar corrections to the Au states. Values for Σ GP and Σ CORR for each junction appear in the SI. This approach was shown to lead to quantitative agreement with photoemission experiments regarding the level alignment of   I: Measured and calculated zero-bias differential conductances (in units of 10 −4 G0) for the different molecule-gold junctions under study. The computed spread reflects the lowest and highest conductance of the three geometries considered for each junction (see SI for details). Experimental data for BP comes from Ref. [15]. All other experimental data is taken from Ref. [22].
In I, we report computed zero-bias conductances for DAS, APA, HDA, and BP junctions for both DFT-PBE and DFT+Σ, and compare with experiments [15,22]. As with previous work [12,15,[41][42][43][44], we find DFT-PBE overestimates measured low-bias conductances for DAS, APA, and BP junctions by more than an order of magnitude. An exception is the HDA junction, where the overestimate is just a factor of 4, consistent with Ref. [55,56]. In contrast, DFT+Σ improves agreement with experiment significantly, predicting low-bias conductances to within a factor of two and shifting frontier orbital resonances to higher energies of 1.5 eV (BP-LUMO) or more away from E F . For DAS, APA, and BP junctions, the total correction Σ GP + Σ CORR opens the gap between junction HOMO and LUMO molecular resonances by more than 2.5 eV. For HDA, we find that Σ GP and Σ CORR are of the same magnitude but opposite sign, resulting in a modest net correction of 0.15 eV, explaining its relatively accurate DFT-PBE zero-bias conductance by a cancellation of errors.
Calculated IV characteristics for each junction are shown in 1 and compared with experiments. We first generate a DFT-PBE steady-state charge density, including a real-axis integration of the scattering-states in the bias window, at each bias voltage; the current is then determined from an integration of T (ω; V ), subsequently computed on a dense energy grid, as described above. Then, in one approach, which we refer to as DFT-PBE(V), T (ω; V ) is generated from DFT-PBE junction electronic structure, which is known to overestimate the zero-bias conductance; in the other approach, referred to as DFT(V)+Σ, T (ω; V ) is generated from DFT+Σ junction electronic structure, with Σ determined at zero-bias.
As can be seen in 1, DFT(V)+Σ leads to excellent quantitative agreement with experiment, whereas the DFT-PBE(V) currents are too large by an order of magnitude (or more). Interestingly, integration of the zerobias DFT-PBE transmission function (gray curves), assuming a uniform potential drop across the junction, results in a good estimate of the DFT-PBE(V) IV characteristics for DAS, APA, and HDA junctions, suggesting that T (ω; V ) is weakly dependent on applied bias for these junctions. A similar integration results in an overestimate by a factor of 8 to 10 for the BP junctions, which we explain below.
In 2, we plot the evolution of the DFT-PBE charge density, ∆ρ, and electrostatic Hartree potential, ∆U H , with bias. ∆ρ (V ) and ∆U H (V ) are both computed relative to zero-bias densities and potentials, respectively. We compare this evolution to that generated assuming linear responses of the density and the Hartree potential, i.e. ∆ρ ∼ ∂ρ ∂V V =0 × V and ∆U H ∼ ∂U H ∂V V =0 × V , respectively. We find that all 4 junctions remain in the linear response regime for the experimental range of bias (±0.8 V). This finding is consistent with the fact that the DFT+Σ molecular resonances are far from the bias window in these junctions, even for the highest biases achieved, leading to modest changes in occupation numbers. Indeed, a simple estimate of the change in charge density under bias, obtained by fitting the DFT(V=0)+Σ transmission function to a symmetric Lorentzian model for all four junctions (see SI), leads to a maximum expected change in occupation by only 0.007 e − .
Interestingly for BP junctions (2), the DFT-PBE(V) densities and Hartree potentials remain in the linear response regime to 1 V, even though the DFT-PBE level alignment erroneously places the LUMO resonance in the bias window at this voltage. This indicates that, despite the incorrect DFT-PBE energy level alignment, self-consistent DFT-PBE(V) results are nonetheless in line with the DFT+Σ level alignment for the prediction of the densities under bias and do not show significant charging or any (non-physical) deviation from linear response.
In 3, we report the evolution of the molecular resonances energies with bias for two BP junctions. For a symmetric junction, the resonance energies vary as V 2 , as expected from a quadratic Stark effect, with a maximum shift (relative to the average chemical potential) of less than 0.06 eV for DAS and APA, and up to 0.2 eV for HDA and BP (3 (a)). The small magnitude of these shifts explains why a simple integration of the zero-bias T (ω; V = 0) is so effective for these junctions. However, if the junctions are more asymmetric, such as the BP junction shown in 3 (b), the resonance energies may vary more substantially. In fact, for the asymmetric BP junction shown in 3 (b), the LUMO resonance energy shifts linearly with V , with small second-order corrections to the energy. Appreciable first-order corrections to the wavefunction are also evident in this case, via a 50% reduction in the resonant peak maximum. Integration of the zero-bias T (ω; V = 0) is clearly much less effective in this case, as T (ω; V ) is varying significantly with V .
Since the Stark shifts of the LUMO resonance for BP junctions are well-described by a simple function of V , a modified self-energy correction of the form Σ (V ) = Σ GP +Σ CORR +∆E (V ), with ∆E (V ) from a polynomial fit of resonance energy shifts from just a few finite bias calculations, results in predicted currents in close agreement with those computed self-consistently at each V. Since our "Stark-corrected" self-energy Σ (V ) is calculable at very low-bias or from the zero-bias polarizability, it can be used for efficient and quantitative prediction of nonlinear IVs and rectification ratios, at least for offresonance tunneling through non-degenerate levels [61].
Our results for BP junctions, shown in Figures 2 and 3, are consistent with our findings for DAS, APA, and HDA. For each junction, we can draw the following conclusions. First, extrapolations using the zero-bias response functions ∂ρ and ∂U H ∂V V =0 result in accurate densities and potentials at other experimentally achievable biases, obviating the need in these cases to do laborious selfconsistent calculations at many bias voltages. Second, self-consistency at DFT-PBE level is apparently sufficient for these junctions, despite errors in level alignment, provided that the actual resonances are far from the bias window. Third, a straightforward Stark correction ∆E (V ) to the DFT-PBE resonance peak energy at zero-bias E (V = 0), along with a zero-bias self-energy correction, yields accurate IV characteristics, even for asymmetric junctions.
The success of DFT+Σ at finite-bias is noteworthy and deserves further comment. Because of the small polarizabilities of the molecules considered here, modifications of Σ GP in relevant electric fields are just 0.01 eV for DAS and APA, 0.02 eV for HDA, and −0.002 eV for BP, a mere few % correction to the zero-bias self-term [62]. Moreover, using the DFT+Σ frontier orbital energies, the weak densities of states in the bias window for all junctions give rise to negligible changes in occupation under bias, smaller than 0.007 e − , implying equally negligible changes in Σ CORR for these junctions in the off-resonance limit. Thus corrections to Σ (V = 0) at finite bias will be extremely small. We expect the good agreement between the DFT+Σ IV characteristics and experiments in 1 to deteriorate for more polarizable molecules (due to the changes in the Σ GP term), stronger coupling of resonances to the lead states (via a breakdown of the DFT+Σ approximation), and for molecular junctions beyond the linear response regime of the density (associated with changes in Σ CORR term, and high-order terms in V in ∆E (V )), e.g. for molecular resonances closer to the Fermi energy of the metal.
To conclude, we have computed IV characteristics for four different molecular junctions, and compared them with STM-BJ experiments. While state-of-the-art methods based on DFT-PBE largely overestimate measured currents, our pragmatic first-principles approach based on one-shot self-energy corrections that improve level alignment in the junction leads to excellent quantitative agreement. This method opens up new avenues for computationally inexpensive and quantitative modeling of non-equilibrium steady-state properties of molecular junctions: our results suggest that the finite-bias transmission function for both symmetric and asymmetric molecular junctions can be well-approximated by calculating the changes in zero-bias resonance energies under small perturbative fields. These findings offer the possibility to more rapidly understand and predict functional nonlinear properties of junctions, such as the rectification, accurately and efficiently. . Computational resources were provided by NERSC. Experimental methods, level alignment and self-energy corrections, analytical changes in the density of a Lorentzian model, basis sets, geometries, plot of the self-consistent density matrices, Hartree potentials and transmission functions at finite biases are reported in the supporting information. This material is available free of charge via the Internet at http://pubs.acs.org. * Electronic address: pdarancet@lbl.gov † Electronic address: lv2117@columbia.edu ‡ Electronic address: jbneaton@lbl.gov