Role of the Binding Motifs in the Energy Level Alignment and Conductance of Amine-Gold Linked Molecular Junctions within DFT and DFT + Σ

: We investigate, using density functional theory (DFT), the electronic and conducting properties of benzenediamine connected to gold electrodes via different tip structures. We examine a series of binding motifs to the electrodes and calculate the junction spectral properties. We consider corrections to the position of molecular resonances at the junction and discuss different approaches to the calculation of these shifts. We relate the magnitude of these corrections to resonance energies to the atomistic structure of the tip. Benzenediamine DFT-based transmission spectra can be well approximated by a Lorentzian model involving only the highest occupied molecular orbital (HOMO). We show how benzenediamine calculated conductance values in quantitative agreement with previous experiments can be achieved from the combination of DFT-based spectra and corrections to the DFT-based HOMO energy and an accessible Lorentzian model.


Introduction
Understanding and controlling charge transport in single molecule junctions have been important since it was first suggested that a single molecule might function as an active electronic component [1], and therefore, that metal-organic interfaces are building blocks for the next generation of electronic devices [2,3]. Charge transport properties are mainly measured using scanning-probe methods, in particular, scanning tunneling microscopy break junctions (STM-BJ) [4,5], and mechanically-controlled microscopy break junctions (MC-BJ) [6,7]. Such measurements are typically carried out in solution and at room temperature, and data from thousands of measurements are compiled to generate conductance histograms [8][9][10]. However, in those experiments, the geometry of the interface on the atomic scale is not known, as it cannot be measured in situ, and it is also changing during the experiment, or from sample to sample. Therefore, theoretical methods to understand and guide experiments are extremely valuable [3,11].
Density functional theory (DFT) is the fundamental approach to calculate the electronic properties and optimize the geometry of the junction [12,13]. Similarly, electron transport calculations are normally performed within the non-equilibrium Green's function (NEGF) formalism. DFT-NEGF, even at zero bias (DFT-Landauer), has been able to correctly reproduce trends in conductance with the correct physical picture [14][15][16][17][18][19][20]. However, despite the significant progress that DFT simulations have offered, their limitations are well known. Among them, perhaps the most relevant for electron transport is the underestimation of the fundamental energy gap at the metal/molecule interface [21][22][23][24]. To go beyond the semi-quantitative picture provided by DFT, it would be desirable to make quantitative comparisons with experimentally measured conductance. For this, corrections to the DFT electronic structure at the interface need to be made. The most widely used approach is the DFT + Σ method [25][26][27][28][29][30][31], where corrections to conducting orbital(s) are calculated externally and added onto the converged system Hamiltonian after the DFT cycle. The correction is a self-energy composed of two terms: one that corrects the molecular gas-phase and one that accounts for the nonlocal polarization due the metal electrodes [25,32,33].
A typical molecular junction is composed of a single molecule bonded to electrodes on either side using chemical linker groups. Electrodes are generally metallic, and Au has been the most commonly used electrode material in scanning-probe studies because of its inertness, which enables consistent and reproducible measurements over a wide range of conditions [5,34]. Linker groups connect the molecule to the electrodes mechanically and electronically [35]. Depending on their chemical nature, linker groups bind to the electrodes either forming donor-acceptor bonds with surface asperities (for example, amine groups -NH 2 ) [25,36], geometry-dependent binding using pyridine groups (−NC 5 H 5 ) [37,38], covalent bonding (such as thiol groups −SH) [39][40][41] or through carboxylic groups (−COOH) [42,43]. The atomistic details of the metal-molecule interface have been shown to tune the electronic and conducting properties of the junction [44][45][46][47]. Experimentally, variations in conductance strongly depend on the nature of the linker group. While for thiolate-Au linkers, they are very large [9,[47][48][49][50], for amine-terminated molecules, these are significantly smaller [25,26,33,51]. The benzenediamine (BDA)-Au interface is a prototypical metal-molecule structure that has been extensively investigated experimentally [25,33,51] and theoretically [26,27,52,53], and is, therefore, an excellent benchmark. STM-BJ studies of BDA between Au electrodes yield a conductance peak at 6.4 × 10 −3 G 0 [25,51,52].
In this work, we focus on the simulation of Au-BDA-Au junctions and investigate how sensitive the electronic and conducting properties are with respect to the atomistic termination of the Au tip structures between molecules and electrodes. We carry out this analysis using both DFT and DFT + Σ formalisms, where, for the latter, we also discuss in detail the magnitude of the necessary corrections to DFT orbital energies at the junction. This study illustrates, for a range of BDA interface geometries, how simple post-processing corrections to DFT-based transmission properties can achieve very good agreement with measured conductance values.

Methods
We employed the well-known SIESTA and TranSIESTA codes [13,[54][55][56]. We constructed the junction geometries from the knowledge that the amine groups bind selectively to undercoordinated Au sites on the metal (111) surface [5,25]. We modeled the N-Au contact using three possible motifs consisting of one, three or four Au atoms (corresponding to adatom, trimer or pyramidal tips, respectively) and considered all possible combinations of these motifs. The resulting six structures are shown in Figure 1: adatom-adatom, adatom-trimer, trimer-trimer, adatom-pyramid, trimer-pyramid and pyramid-pyramid. We relaxed the positions of the atoms in the molecule and Au tip atoms until the residual Hellman-Feynman forces fell below 0.02 eV/Å. We used an exchange correlation functional which accounts for van der Waals (vdW) interactions [57]. For structure optimization, the real-space grid was defined with an equivalent energy cut-off of 250 Ry, while the reciprocal space was sampled using a 2 × 2 × 1 Monkhorst-pack mesh. For calculations of density of states, we used a 5 × 5 × 1 Monkhorst-pack grid. Subsequent calculations of charge transport were performed at zero bias for optimized geometries using the Landauer formalism [58] as implemented by TranSIESTA [13,56]. Au atoms were described using a single-ζ polarized basis set, while a double-ζ polarized basis was used for molecular atoms. In these calculations, 5 × 5 × 1 and 15 × 15 × 1 Monkhorst-pack grids were used in reciprocal space for calculating the Green's function and transmission spectra, respectively.
We also employed the DFT + Σ method [25][26][27]29,30], which we implemented in SI-ESTA. In this approach, corrections are added explicitly into the system Hamiltonian (H→H + Σ) by an orbital dependent operator of the form, where Σn is the self-energy correction for the nth molecular orbital, and mol denotes the wavefunction states of the molecule. These states are calculated from the Hamiltonian of the molecular subspace Hmol, contained into the Hamiltonian of the total system (H). The correction operator (Equation (1)) acts only on the molecular subspace Hmol ⊆ H, by construction [26,33,59,60]. The self-energy operator Σ can be constructed from a separate calculation of the relaxed isolated molecule [27,60] or from the Hamiltonian of the molecular subspace Hmol directly cropped from the converged ground-state Hamiltonian of the junction. Either way, the correction operator is introduced into the total Hamiltonian, which is diagonalized again in order to obtain the corrected electronic properties of the system. The self-energy correction term consists of two parts that address: (1) the underestimated gap of the isolated molecule in conventional DFT (Σ 1 n, gas phase correction) and (2) the lack of renormalization due to the metallic electrodes (Σ 2 n, polarization due to metallic surface). The total correction is the sum of these two contributions which, as discussed below, have opposite signs, i.e., Σn = Σ 1 n + Σ 2 n. In principle, this correction could be calculated for every molecular orbital n. However, since we are interested in conductance, we focused on the region around the Fermi level, and calculated it only for the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). We applied ΣHOMO to all occupied states and ΣLUMO to all unoccupied states [26,53].
The first contribution to the self-energy, the gas phase correction Σ 1 n, is defined as the difference between the DFT energy level and the quasiparticle energy level. We calculated the ionization potential (IP), i.e., the energy required to remove an electron from the ground state, and the electron affinity (EA), i.e., the energy required to add an electron to Subsequent calculations of charge transport were performed at zero bias for optimized geometries using the Landauer formalism [58] as implemented by TranSIESTA [13,56]. Au atoms were described using a single-ζ polarized basis set, while a double-ζ polarized basis was used for molecular atoms. In these calculations, 5 × 5 × 1 and 15 × 15 × 1 Monkhorst-pack grids were used in reciprocal space for calculating the Green's function and transmission spectra, respectively.
We also employed the DFT + Σ method [25][26][27]29,30], which we implemented in SIESTA. In this approach, corrections are added explicitly into the system Hamiltonian (H→H + Σ) by an orbital dependent operator of the form, where Σ n is the self-energy correction for the nth molecular orbital, and ψ mol n denotes the wavefunction states of the molecule. These states are calculated from the Hamiltonian of the molecular subspace H mol , contained into the Hamiltonian of the total system (H). The correction operator (Equation (1)) acts only on the molecular subspace H mol ⊆ H, by construction [26,33,59,60]. The self-energy operator Σ can be constructed from a separate calculation of the relaxed isolated molecule [27,60] or from the Hamiltonian of the molecular subspace H mol directly cropped from the converged ground-state Hamiltonian of the junction. Either way, the correction operator is introduced into the total Hamiltonian, which is diagonalized again in order to obtain the corrected electronic properties of the system.
The self-energy correction term consists of two parts that address: (1) the underestimated gap of the isolated molecule in conventional DFT (Σ 1 n , gas phase correction) and (2) the lack of renormalization due to the metallic electrodes (Σ 2 n , polarization due to metallic surface). The total correction is the sum of these two contributions which, as discussed below, have opposite signs, i.e., Σn = Σ 1 n + Σ 2 n . In principle, this correction could be calculated for every molecular orbital n. However, since we are interested in conductance, we focused on the region around the Fermi level, and calculated it only for the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). We applied Σ HOMO to all occupied states and Σ LUMO to all unoccupied states [26,53].
The first contribution to the self-energy, the gas phase correction Σ 1 n , is defined as the difference between the DFT energy level and the quasiparticle energy level. We calculated the ionization potential (IP), i.e., the energy required to remove an electron from the ground state, and the electron affinity (EA), i.e., the energy required to add an electron to the ground state. These quantities are defined in terms of (DFT) total energies as IP = E(N) − E(N − 1), and EA = E(N + 1) − E(N), where E(N i ) is the total energy of the system with N i electrons [28,29,61]. Furthermore, the gas-phase correction to the HOMO (LUMO) is calculated as a difference between the DFT energy position of the isolated molecule E DFT HOMO (E DFT LUMO ) and IP (EA), The second term, accounting for the polarization due to the metallic surface Σ 2 n , is approximated by a classical image charge model [32,61]. It is modeled as the potential energy of a point charge distribution between two image planes. We used a point charge q n = 1 for the n = HOMO or LUMO, corresponding to each atom in the molecule. Each point charge q n,i is located at a vertical distance z i from the image plane z top (z bottom ), corresponding to the top (bottom) electrode. Therefore, the self-energy term is calculated as

Results and Discussion
We began our analysis by studying the electronic properties of the junctions within DFT, as well as their variation with respect to the binding motifs ( Figure 1). We focused only on the position of the HOMO peak, as it has been shown in the literature that HOMO dominates the zero-bias conductance for BDA [25,26,51,52]. Figure 2 shows the density of states projected over the atoms of the molecule for all the binding motifs considered: adatom-adatom, adatom-trimer, trimer-trimer, adatom-pyramid, trimer-pyramid and pyramid-pyramid. The figure is centered around the energy of the DFT HOMO peak (−1.0 eV). The inset shows the same data on an extended energy range. The data are offset vertically for clarity. From the figure, the position of the HOMO peak at the DFT level had a small variation across all different structures, ranging from −0.83 to −1.09 eV. The LUMO peak was found between 2.74 and 2.92 eV. From Figure 2, tip structures involving adatoms broadly resulted in sharper PDOS peaks and resonances closer to the Fermi level. In order to implement the DFT + Σ method, it is necessary to calculate the magnitude of the self-energy correction of the appropriate resonance (Equations (2) and (3)). We first address the calculation of the gas-phase self-energy correction of the HOMO peak, Σ 1 HOMO , which involves the calculation of the isolated molecule. Several ways of computing this correction are possible. One option is to optimize the geometry of the molecule. Another possibility is to compute the molecule in the geometry it adopts at the junction. In the case of BDA, the magnitude of both corrections was the same in all cases regardless of which approach was used, −2.99 eV, although we believe that this is due to the reduced conformational flexibility of the amine linker. We anticipated a spread of values for other linkers, such as methyl-sulfide groups, which rotate when adsorbed with respect to their gas-phase geometry. In this case, it would be more consistent to use the geometry the molecule adopts at the junction [20,35].  The total self-energy correction for the HOMO peak is the sum of both (opposing) contributions. Results for all tip structures are presented in Table 2. As expected, the dependence with tip structure is opposite that of Table 1: the largest shift was found for the pyramid-pyramid combination, where the polarization self-energy was lowest, and structures with one pyramid followed. These differences are relevant for the calculation of conductance in subsequent sections. Furthermore, the total self-energy correction for the LUMO peak was 1.90 eV for the adatom-adatom, adatom-trimer and trimer-trimer tip structures, 2.09 eV for the adatom-pyramid and trimer-pyramid tip structures and 2.31 eV for the pyramid-pyramid tip structure. The second term in the self-energy correction, accounting for the polarization due to the metallic surface, Σ 2 HOMO , was approximated using a classical image charge model, as described in Equation (4). A single point charge was positioned at the geometrical center of the molecule [28,29,53]. The image plane position was taken to be 1.0 Å above the outer atomic plane of each Au [62]. Results for all tip structures are presented in Table 1. The calculated shift of the HOMO resonance due to screening is, for both approaches, given in Table 1. We see that the largest values were obtained for "short" tips that protrude the least from the surface, such as adatoms or trimers. When pyramidal tips were considered, image charge screening was reduced due to the larger vertical distance spanned by these tips. The total self-energy correction for the HOMO peak is the sum of both (opposing) contributions. Results for all tip structures are presented in Table 2. As expected, the dependence with tip structure is opposite that of Table 1: the largest shift was found for the pyramid-pyramid combination, where the polarization self-energy was lowest, and structures with one pyramid followed. These differences are relevant for the calculation of conductance in subsequent sections. Furthermore, the total self-energy correction for the LUMO peak was 1.90 eV for the adatom-adatom, adatom-trimer and trimer-trimer tip structures, 2.09 eV for the adatom-pyramid and trimer-pyramid tip structures and 2.31 eV for the pyramid-pyramid tip structure. Having described the spread of the self-energy values with respect to tip structure, we applied the correction to the Hamiltonian (H→H + Σ) as described in Equation (1). The shift was applied to the molecular Hamiltonian. DFT values are given by the eigenstates of the molecular box of the junction Hamiltonian, as described previously. The initial and final (corrected) energies of the HOMO resonance are given in Table 3. Figure 3 shows the corrected density of states projected onto the atoms of the molecule for tip structure combinations. The figure is centered on the range around −3.0 eV, near the position of the corrected HOMO peak. The inset reproduces the same data on a linear scale and over an extended energy range with the curves offset for clarity. The LUMO peak was shifted following the same procedure described here. The position of the corrected HOMO peak ranged between −2.9 and −3.2 eV. Additionally, the position of the corrected LUMO peak ranged between 4.5 and 4.8 eV. Table 3. Highest occupied molecular orbital (HOMO) peak position, calculated using the DFT and DFT + Σ method.

DFT
DFT + Σ binations. The figure is centered on the range around −3.0 eV, near the position of the corrected HOMO peak. The inset reproduces the same data on a linear scale and over an extended energy range with the curves offset for clarity. The LUMO peak was shifted following the same procedure described here. The position of the corrected HOMO peak ranged between −2.9 and −3.2 eV. Additionally, the position of the corrected LUMO peak ranged between 4.5 and 4.8 eV. So far, we have focused on the DFT and corrected electronic properties, discussing how to calculate the magnitude of these corrections. In the final section of the paper, we turn to the electron transport properties and how to apply these corrections to DFT-based So far, we have focused on the DFT and corrected electronic properties, discussing how to calculate the magnitude of these corrections. In the final section of the paper, we turn to the electron transport properties and how to apply these corrections to DFTbased conductance calculations. Figure 4 shows the transmission spectra of the different BDA junctions calculated using the DFT-Landauer formalism. These calculations take the DFT-based electronic structure as input, with its well-known errors in resonance position. Figure 4 highlights the energy range below the Fermi level where the HOMO resonance, which defines zero-bias conductance, appears at the DFT level. As before, the inset plots the same data on a linear scale over an extended energy range. Figure 4 shows that lowbias conductance is determined by the tailing of the HOMO resonance into the Fermi level [25,26,51,52]. In Figure 4, the HOMO peak position is found between −0.85 and −1.00 eV over the tip structures considered, obviously following the DFT-based density of states (see Figure 2 and Table 3). This is to be compared to the HOMO transmission peak between −0.94 and −1.27 eV reported in the literature [25,26]. Calculated conductance values, given by the transmission at the Fermi level, are given in Table 4. Furthermore, as is common for DFT-based approaches, calculated conductance significantly exceeds the experimental value (6.4 × 10 −3 G 0 ) [25,51,52], in this case, by about a factor 10.    To improve the agreement between calculated and measured conductance, it is necessary to correct the position of the conducting orbital. The implementation of DFT + Σ in a fully selfconsistent DFT-NEGF cycle, out of equilibrium, is far from simple. However, trends can often be drawn from studies in equilibrium. Since in many molecular junctions, conductance takes place due to non-resonant tunneling [3,11,34], it is illustrative to consider simple models involving only one molecular resonance. Although single level models fail when transport involves several molecular orbitals [63] or in cases of quantum interference [64], they nevertheless provide a good starting point for many representative molecular junctions. Different ways of fitting the relevant parameters and their accuracy have been discussed [9,35,44,59,[63][64][65][66][67]. Here, we consider a single level (HOMO) at zero bias.
First, the DFT-based conducting peak (in this case HOMO) is fitted using a Lorentzian model of the form, where E HOMO and A are the HOMO peak position and amplitude, respectively, and Γ is the full width at half maximum. The position of the resonance is then corrected using the selfenergy previously addressed (Σ HOMO ). This model assumes that width of the Lorentzian (Γ) is unchanged, which is reasonable since DFT captures the electronic coupling well. For the BDA junctions considered, Γ takes three range of values. For adatom-adatom, adatomtrimer and trimer-trimer structures, the widths were 0.51, 0.52, and 0.55 eV, respectively. For adatom-pyramid and trimer-pyramid structures, the widths were 0.43 and 0.45 eV. Finally, for the pyramid-pyramid structure, the width was the lowest, 0.35 eV. This agrees well with previous calculations, where the conductance peak width showed a modest variation between 0.34 to 0.56 eV [25]. The conductance from the Lorentzian fitted curve is presented in Table 4. The agreement with DFT-Landauer values was very good, where the structures with the widest peaks showed a conductance closest to DFT-Landauer. Table 4 also reports the conductance values calculated from the Lorentzian fit using the DFT + Σ resonance positions. Agreement with experiment was significantly improved. For the adatom-pyramid and trimer-pyramid structures, the corrected conductance was 1.0 × 10 −2 G 0 . This value still exceeded the experimental conductance, but only for a factor of about 1.6, a substantial improvement as compared to the previous factor of 10. For the adatom-pyramid and trimer-pyramid structures, the corrected conductance was about 6.3 × 10 −3 G 0 , which perfectly matches with the experimental value [25,51,52]. Furthermore, for the pyramid-pyramid structure, the corrected conductance was 3.5 × 10 −3 G 0 , even falling below the experimental conductance value by a factor of about 1.8. Additionally, the variation between the calculated conductance values among all tip structures was modest, about 18%. This is comparable to the experimental variation of conductance, about 8% [25]. Therefore, the combination of corrections to resonance positions within the DFT + Σ methodology, and a Lorentzian model of conductance, leads to a very good quantitative agreement with measured conductance values of BDA. The approach that we have described here would be of use to calculate the level-corrected conductance for broad classes of molecular junctions from standard DFT-Landauer calculations and rather straightforward corrections to resonance positions, combined with a Lorentzian transport model.

Conclusions
We studied the simulation of electronic and charge transport properties of BDA junctions with different tip structures. As is well known, at the DFT level, the HOMO position is too close to the Fermi level, and we discussed corrections to the DFT-based HOMO energy within the DFT + Σ formalism. We discussed the two contributions to the resonance energy correction, arising from the self-energy of the isolated molecule and from the screening at the metallic interface. This correction shifts the HOMO resonance further form the Fermi energy. We found that, for BDA junctions, the first term was not sensitive to the details of its calculations. However, the contribution due to interface screening, which was approximated here by a classical image charge model, did show a substantial variation of several tenths of an eV with tip structure. The total correction to the HOMO resonance at the interface was close to 2 eV towards more negative values. For the interfaces considered, DFT-based transmission spectra of BDA yielded conductance values that significantly overestimate the measured conductance by as much as an order of magnitude. We found that a Lorentzian model considering the HOMO resonance matched the DFT transmission spectra well. We showed that a DFT + Σ approach where this Lorentzian model was combined with the corrected HOMO energy produced a significant improvement in the calculated conductance values. For the different interface structures considered, this approach resulted in values in quantitative agreement with experiments.

Data Availability Statement:
The data presented in this study are available on reasonable request from the corresponding authors.

Conflicts of Interest:
The authors declare no conflict of interest.