Geometry optimizations were performed using the BP86 generalized-gradient-approximation (GGA) exchange-correlation functional in combination with the Ahlrichs def2svp basis set. Computations were density-fitted using the auxiliary basis set by Weigend et al., and corrected for dispersion using Grimme's D3 method with Becke-Johnson damping. Vibrational frequencies of optimized geometries were verified to be real and final electronic wavefunctions were checked for stability. Energies and molecular orbitals of final junction structures were calculated using the range-separated and dispersioncorrected ωB97XD hybrid GGA functional in combination with the Ahlrichs def2tzvpp basis set. We verified that relativistic corrections make minor modifications to φ b and do not change any qualitative trends shown in this work. The energies of the optimized junction (E junction ), the same junction geometry minus the molecule (E electrodes ), and the optimized isolated molecule (E molecule ) were used to calculate the binding energy of the junction,