Density Functional Theory simulations were performed with the SIESTA package. We used the van der Waals density functional by Dion et al, with the modified exchange correlation by Klimeš, Bowler and Michaelides. The valence-electron wave functions are expanded using a double-ζ plus polarization (DZP) basis set with the orbital radii extension defined using a 30 meV energy shift, 6 while the core electrons were replaced by normconserving Troullier-Martins pseudopotentials. The fineness of the real space grid used for integration was defined by setting an energy cutoff of 300 Ry. smearing of the electronic occupations was determined using a Fermi-Dirac distribution, and 300 K of electronic temperature. Self-consistency was considered to be achieved when the changes of the density matrix elements were smaller than 10 -6 as well as smaller than 10 -5 eV for the Hamiltonian matrix elements. We used periodic boundary conditions with sufficiently large unit cells (at least 20 Å of vacuum space along all directions was used for determining the simulation cell for all systems, in order to avoid interaction between replicas in neighbouring cells). Geometry optimizations were performed using the conjugate gradient method until all forces were lower than 10 meV/Å. Constant-height STM simulations were performed in the Tersoff-Hamann approximation 11 using the local density of states, taken 4 Å above the molecular plane to account for the tip-sample distance and the locally different orbital decay into the vacuum. The STM post-processing tool distributed with SIESTA was used for an accurate extrapolation of the local density of states at such high distance over the molecular plane, avoiding spurious effects associated with the finite range of the numerical atomic orbitals used as a basis set.