## 2 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/A. 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\times 2\times 1\) Monkhorst-pack mesh. For calculations of density of states, we used a \(5\times 5\times 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-\(\zeta\) polarized basis set, while a double-\(\zeta\) polarized basis was used for molecular atoms. In these calculations, 5 \(\times\) 5 \(\times\) 1 and 15 \(\times\) 15 \(\times\) 1 Monkhorst-pack grids were used in reciprocal space for calculating the Green's function and transmission spectra, respectively. We also employed the DFT + \(\Sigma\) method [25; 26; 27; 29; 30], which we implemented in SIESTA. In this approach, corrections are added explicitly into the system Hamiltonian (H\(\rightarrow\)H + \(\Sigma\)) by an orbital dependent operator of the form,  
\[\Sigma=\sum\Sigma_{n}\left|\psi_{n}^{\text{mol}}\,\psi_{n}^{\text{mol}}\right|, \tag{1}\]  
where \(\Sigma_{n}\) is the self-energy correction for the \(n\)th molecular orbital, and \(\left|\psi_{n}^{\text{mol}}\right|\) denotes the wavefunction states of the molecule. These states are calculated from the Hamiltonian of the molecular subspace H\({}_{\text{mol}}\), contained into the Hamiltonian of the total system (H). The correction operator (Equation (1)) acts only on the molecular subspace H\({}_{\text{mol}}\subseteq\) H, by construction [26; 33; 59; 60]. The self-energy operator \(\Sigma\) can be constructed from a separate calculation of the relaxed isolated molecule [27; 60] or from the Hamiltonian of the molecular subspace H\({}_{\text{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 (\(\Sigma_{n}^{1}\), gas phase correction) and (2) the lack of renormalization due to the metallic electrodes (\(\Sigma_{n}^{2}\), polarization due to metallic surface). The total correction is the sum of these two contributions which, as discussed below, have opposite signs, i.e., \(\Sigma n=\Sigma_{n}^{1}+\Sigma_{n}^{2}\). 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 \(\Sigma_{\text{HOMO}}\) to all occupied states and \(\Sigma_{\text{LUMO}}\) to all unoccupied states [26; 53]. The first contribution to the self-energy, the gas phase correction \(\Sigma_{n}^{1}\), is defined as the difference between the DFT energy level and the quasiparticle energy level.