Bias-dependent local structure of water molecules at a metallic interface

Abstract

Understanding the local structure of water at the interfaces of metallic electrodes is a key issue in aqueous-based electrochemistry. Nevertheless a realistic simulation of such a setup is challenging, particularly when the electrodes are maintained at different potentials. To correctly compute the effect of an external bias potential applied to truly semi-infinite surfaces, we combine Density Functional Theory (DFT) and Non-Equilibrium Green’s Function (NEGF) methods. This framework allows for the out-of-equilibrium calculation of forces and dynamics, and directly correlates to the chemical potential of the electrodes, which is introduced experimentally. In this work, we apply this methodology to study the electronic properties and atomic forces of a water molecule at the interface of a gold surface. We find that the water molecule tends to align its dipole moment with the electric field, and it is either repelled or attracted to the metal depending on the sign and magnitude of the applied bias, in an asymmetric fashion.

---

1 Introduction

Following the need for new and renewable sources of energy worldwide, fuel cells using electrocatalysts can be thought of as viable options.^1,2 The interface between metal (electrode) and water in these systems is the electrochemical central point, since it is the region where charge transfer can take place. A better understanding of the metal–water interface is also an essential requisite for predicting the correspondence between the macroscopic voltage and the microscopic interfacial charge distribution in electrochemical fuel cells. This reactivity is governed by the explicit atomic and electronic structures built at the interface as a response to external conditions, such as an applied potential.^3–5

The advance in experimental techniques for studying surfaces in the last decades started to provide important results concerning the local structure of water at interfaces, revealing a bias-dependent behavior.^6–8 Notably, in a more realistic system applied to catalysis, the metal will be at a given potential. From a theoretical perspective this makes the task of simulating this setup difficult.⁹ In fact, accurate calculation of the electrostatic potential at electrically-biased metal–electrolyte interfaces is a challenge for ab initio simulations with periodic boundary conditions.^10,11

One possible way to simulate this electrochemical cell under an explicit bias is to account for the polarization of the metal by charging each atom on the electrode, enforcing a constant potential and using the image charges method.^12–14 Although this can provide interesting insights into the problem, the characterization is limited to the use of empirical models. Neurock’s studies of water/metal interfaces in the presence of an applied potential^15,16 were among the first ones to use Density Functional Theory (DFT)^17,18 to address this problem. The resulting electrode potential – which is related to the Fermi energy of the system – is compared to an internal reference potential by artificially inserting a vacuum layer into the center of the solution region. To be able to fully compute all of these energies when the system is charged, the water molecules in the center of the liquid layer are usually held fixed during the optimization of the charged systems. More recently, N. Bonnet et al. proposed a methodology where an external potentiostat was added to the system¹⁹ following a previous work where the charges at the surface were controlled by including a medium with a given permittivity.²⁰ The idea in the later work was to use the potential energy of a fictitious system, akin to the fictitious mass in Car–Parrinello first-principles molecular dynamics. Therefore, current methodologies²¹ attempt to simulate the effect of finite bias at the metal by altering their charge (adding/subtracting electrons), whereas in experiments the potential is the quantity that is controlled.

The electrochemical cell can be thought of as two metallic electrodes which act as charge reservoirs, with the two metal plates separated by a solution (mostly composed of water). This arrangement is analogous to the one encountered in simulations of electronic transport: a central scattering region coupled to electrodes.^22–24 Thus, we propose in this work, as an alternative, to use open boundary conditions by employing the non-equilibrium Green’s function (NEGF) method combined with DFT to properly compute the effect of an external bias potential applied to electrodes. While standard DFT implementations are not suited to treating extended systems under an external bias, NEGF has been designed to treat out-of-equilibrium situations. Most notably it allows for the inclusion of truly semi-infinite metallic electrodes, which set the correct chemical potential for the metal, and a clear reference potential. Their combination has been developed over the past decade to describe the current–voltage characteristics of nanoscopic systems.^23–25 It treats an open system under the influence of an external bias, and although dynamics – or forces – are typically ignored in such systems, they can be incorporated into the methodology. In this work, we apply this framework to a system consisting of a single water molecule between two Au(111) surfaces, at different configurations and as a function of an external voltage.

2 Methodology

2.1 General methodology

In what has become the usual approach in electronic transport calculations, the system is divided into three regions: electrodes (left and right, L/R) and scattering region (SR).²⁶ We framed our metal–water system into an analogous arrangement. In this case, the interfacial region (metal–water–metal) represented the scattering region. The first few layers of the metal at the interface needed to be considered part of the scattering region as we required that the charge density at the edge of the SR resembled that of the bulk metal. Then, the electronic charge distributions in the electrodes, left and right, corresponded to the bulk phases of the same material to a prescribed numerical accuracy. A representation of the typical arrangement used in our calculations is shown in Fig. 1(a).

Fig. 1 (a) A schematic view of the metal–water system used for the non-equilibrium calculations; the left (right) electrodes (LE/RE) and scattering region (SR) are indicated. (b) A sketch of the effect of a positive bias potential on a parallel plate capacitor; the corresponding charge accumulated in each plate, as well as the bias ramp, is shown.

When a finite voltage is applied to the electrodes the problem becomes a non-equilibrium one. The electrodes are then ascribed different chemical potentials and current, in principle, can flow. The non-equilibrium Green’s function formalism is a general formalism for calculating the properties of systems in out-of-equilibrium situations, and can be used to tackle our problem of the electrochemical cell. In principle it can be used to address problems where inelastic effects are present, and most importantly, it goes far beyond electronic transport (including ballistic transport).

Within the NEGF approach, if the Hamiltonian can be cast in a bilinear form, the entire problem can be treated in terms of a single Green’s function; in our case the retarded Green’s function for the scattering region is:

where ε⁺ = E + iη, S_SR is the overlap matrix and H_SR[n] is a Hamiltonian which is a functional of the charge density, n(r⃑). In this work the Hamiltonian was taken as the Kohn–Sham (KS) Hamiltonian from DFT. The effects of the electrodes were introduced in the form of the self-energies Σ_L/R, which were obtained by integrating out the degrees of freedom of the leads. As the electrodes were considered to be good metals, the effect of the bias on the left and right electrodes corresponds to a rigid shift (±V/2) of the zero-bias self-energies, setting the boundary condition (illustrated in Fig. 1(b) for a positive bias). This means that the self-energies could be obtained from a separate DFT calculation for the bulk metal, and did not need to be updated self-consistently throughout the calculation. This approach also ensured that a clear reference potential was defined (the chemical μ₀ of the bulk metal) as we assumed the electrodes were charge reservoirs in thermodynamic equilibrium throughout the calculations.

Once the Green’s function of the SR was calculated, all of the observables of the system could be recomputed. In particular the density matrix is expressed as:

where the μ and ν indices run over the SR electronic states, f(E) = 1/(1 + e^E/kT) is the Fermi distribution and μ_R and μ_L are the electrochemical potentials of the right and left electrodes (μ_L/R = μ₀ ± V/2) that define the bias: V = μ_L − μ_R. Finally,

ρ^L/R = G(E)Γ_L/R(E)G^†(E) (3)

is the electrode spectral density matrix, obtained from the Green’s function of the SR and the left (right) coupling matrices, . From the density matrix, the KS Hamiltonian could be computed. The procedure was then repeated until self-consistency was achieved.

Within the ground state DFT framework, the computation of forces on the nuclei was theoretically well-founded thanks to the Hellmann–Feynman theorem,^27,28 and the forces were obtained via the derivative of the total energy. Using a set of localized basis functions the force was decomposed into two terms:²⁹ one that originated from the derivative of the energy of the occupied eigenstates (band structure contribution) and a second one that contained the remaining contributions to the energy. For an ion I, the former one was given by:

where D_μν is the density matrix and Ω_μν is the energy density matrix . The situation is more complex out of equilibrium, where the Hellmann–Feynman theorem does not apply.³⁰ Recently, it was shown that the forces can actually be obtained by the time derivative of the expectation value of the ionic momentum operators:^29,30

As it turns out, for steady-state problems, the final form for the force is equivalent to the equilibrium case:

which can be expressed by eqn (4), replacing the ground state density matrix and energy density matrix with the out-of-equilibrium ones, obtained in terms of the retarded Green’s function.

One important point that still remains is how well defined the forces are when current flows in the system.^31–33 In our case, however, the gap in the scattering region – the band gap of water – was large enough (∼8 eV) to ensure that no current would flow through the arrangement. In that sense, it is important to stress that current-induced forces were not present in our problem. This remains true for the simulation of ionic electrolytes, where ionic currents are expected to exist – and can be captured by this method – but no electronic currents exist.

Finally, one notices that the above methodology relies on the Kohn–Sham Hamiltonian being a good description of the single particle excitations for the system, as is the case in different implementations of the NEGF formalism within DFT.^23,24 This leads to known pitfalls, which are associated with the positions of the molecular energy levels and charge transfer between the surfaces and molecules to mention a few.^34,35 Most of these issues however, pertain to approximations in the exchange–correlation functional, and corrections in different forms can be readily incorporated into the formalism.^36–40 Nonetheless, it is important to point out that local and semi-local functionals tend to perform better for forces and structures compared to total energies and single particle energy levels.

The described methodology was implemented in the Smeagol code^24,25 which was bundled with Siesta.^41,42 In the same way that relaxation and ab initio molecular dynamics can be performed within DFT, one can now use the code to do the same for out-of-equilibrium systems.

2.2 Details of calculations

In this work, we used two different gradient-dependent exchange–correlation (XC) functionals: PBE⁴³ and vdW-DF^PBE, which includes van der Waals corrections (vdW). The vdW-DF^PBE functional is a modified version of the original vdW-DF functional,⁴⁴ in which the revPBE local term is replaced by PBE.⁴⁵ The core electrons were described by norm-conserving pseudopotentials in the Troullier–Martins form.⁴⁶ A basis set of numerical atomic orbitals with double-ζ polarization was used to describe the valence electrons. For both the metal and water the basis set was variationally optimized and ensured that our results (Au lattice parameter and water–metal geometry) were in agreement with the plane-wave calculations.

For the non-equilibrium calculations, each metal slab within the scattering region had 3 layers of (111) planes with 12 Au atoms on each plane of a size of 10.29 × 9.89 Å in the plane perpendicular to the transport direction. This size was chosen because periodic boundary conditions are still applied in this plane, and it is necessary to minimize interaction between periodic repetitions of the water molecules in the plane. The water molecule was placed close to one metal surface (the one defined as the left). In order to minimize the interaction between the surfaces, the right and left side were 20 Å apart. The electrodes, connected to the scattering region, consisted of 3 Au layers each (left and right). Fig. 1 shows a schematic view of the system and its components.

3 Results

Before applying a bias at the electrodes it was important to characterize the ground state configuration of the metal–water system. This was initially done using a (111) surface Au slab with 4 layers and a 2 × 2 in-plane supercell within the standard periodic DFT formalism. The relaxed water structures were then used as starting configurations for the larger gold surfaces. All of the atoms were allowed to move during the geometrical optimization using the conjugate gradient algorithm and with a 0.005 eV Å⁻¹ tolerance criteria on the forces. The final configurations were very similar to the ones used as a starting point. Our results show that the most stable configuration of the water molecule on top of the Au(111) surface is the so-called “flat” one (i.e. the molecule dipole moment is almost parallel to the surface plane), in agreement with previously reported calculations.^47,48 The molecular plane was slightly tilted with respect to the metal plane, with α = 3° (α being the angle between the molecular plane and the surface plane) and the distance between Au and O d_Au–O = 2.92 Å. This structure is illustrated in Fig. 2(a). We also relaxed the metal–water–metal structure using the NEGF formalism at zero bias, obtaining α = 6° and 2.79 Å for the Au–O distance. The small differences are attributed to the effect of using a finite representation of Au surface in the standard DFT calculation.

Fig. 2 z-Component of the force at the center of mass (right panels) as a function of the vertical displacement for different water configurations (left panels): (a) “flat”, (b) “up”, (c) “flat-up”, (d) “perpendicular” and (e) “down”. The vertical displacement is given with respect to the Au–O distance in the ground state. The insets of the right panels show the regions in each graph for which F^CM_z = 0.

It is worth mentioning that the potential energy surface (PES) for this system was very flat in the region of the water on top of the Au. For instance, the energy difference between the configuration where the molecule was “flat” and the one where it had the hydrogens pointing down (towards the metal) was only ∼0.06 eV. Therefore, we also considered four other different configurations for the water molecule, corresponding to rigid rotations of the ground state structure (left panels of Fig. 2). The geometries were labeled according to their orientation: “up” (hydrogens pointing away from the metal), “down” (hydrogens pointing towards the metal), “perpendicular” (α = 90° and θ = 90°), and “flat-up” (one hydrogen higher than the flat configuration, with α = 24° and θ = 84°). The angle θ corresponded to the angle between the molecule dipole and surface normal.

In order to analyze the effect of different bias voltages (magnitude and sign) on the molecule as a function of the distance of the water molecules to the metal, we first performed non-equilibrium calculations for all of the water structures. For each configuration, we started at the ground state Au–O distance (z = 2.79 Å, which corresponds to zero in the plots) at zero bias and increased/decreased this distance by −0.5 to +2.0 Å, and performed a single point calculation. At each point the forces on the atoms were evaluated for a particular applied bias. Since the water molecule was placed close to one metal surface, we observed that the potential in the water molecule closely followed that of the surface. Therefore, the bias indicated in the plots corresponds to V/2, as it corresponds to the potential that effectively acted on the molecule.

The results for the z-component of the force on the center of mass of the molecule, F^CM_z, are shown in Fig. 2 for all orientations. In general, we observed that low bias (−0.5 and +0.5 V) had a small effect on the forces, independent of the water configuration. However, as we increased the applied bias we observed that the forces close to the minimum were modified. In particular, this effect was more evident for the flat molecule (Fig. 2(a)) and for the configurations where the oxygen was facing the metal, due to strong interaction between the oxygen-b₁ orbital of the water molecule and the metal orbitals.^49,50

Fig. 2 indicates that there was a tendency to modify the position of the minimum configuration when the bias was applied. Moreover, this modification was dependent on the magnitude of the bias and it was asymmetric with respect to positive and negative values. This behavior is similar to what was verified when an electric field was applied:^51,52 the molecule tended to get closer to the metal when the bias was negative and moved away from the metal with a positive bias, evidenced by the position at which F^CM_z = 0. A similar trend can be observed for the “up” and “flat-up” configurations shown in Fig. 2(b and c), respectively. For the “perpendicular” and “down” configurations the molecule was essentially unbound (Fig. 2(d and e), respectively).

A similar behavior was also observed when vdW corrections were included for the flat configuration, as shown in Fig. 3. In agreement with previous work,⁵³ we note that the vdW functional did not significantly change the water–Au interaction. We observe, however, that the barrier in all of the cases increased slightly for higher distances; this is true for zero as well as for finite (positive or negative) bias. This means that, although the equilibrium position of the molecule did not depend on the choice of XC functional (specifically for Au–water systems), the restoration force was larger when we included vdW interactions.

Fig. 3 The magnitude of the force at the center of mass as a function of the vertical displacement for the “flat” configuration for the PBE and vdW-DF^PBE exchange–correlation functionals. The vertical displacement is given with respect to the Au–O distance at ground states corresponding to each functional.

Although by rigidly shifting the position of the molecule along z, one can always find a position for which F^CM_z = 0, that is not the case for all directions concomitantly (see ESI†). This is an indication that, as the absolute value of the bias increases, the orientation of the molecule tends to change as well. Thus, in a second step, starting from the “flat” configuration, we allowed the atoms of the water molecule to move using the conjugate gradient algorithm with a bias applied to the system. The minimum configurations obtained for −1.5 and +1.5 V are shown in Fig. 4, where the geometry obtained for the zero bias case is also shown. The asymmetric behaviour with respect to the bias sign was clearly observed. The geometry for +1.5 V had the hydrogen atoms pointing down and the oxygen atom was 2.84 Å away from the metal. In fact, this is the “down” configuration (α = 90° and θ = 180°), presented in Fig. 2(e), which means that the positive bias leads to an unbound molecule. On the other hand, when the negative bias was applied the hydrogen atoms moved slightly upwards (α = 27° and θ = 63°) and the oxygen got closer to the metal (d_O–Au = 2.69 Å) compared to the neutral case. In essence, one notices that even a relatively small bias can lead to significant structural changes on the metal–water interface.

Fig. 4 Relaxed configurations of the water molecule on the Au surface for different bias voltages (+1.5 V → +, 0 V, and −1.5 V → −), the corresponding angles, and the estimated dipole moment (details in ESI†). The angle α is the angle between the molecular plane and the surface plane, θ is the angle between the isolated molecule dipole and the surface normal, and γ is the angle between the isolated molecule dipole and that calculated from the charge density of the combined Au + H₂O system.

The water–Au interaction is mostly electrostatic in nature, and there is almost no charge transfer between water and the metal in the neutral case,⁵³ as seen from a Bader analysis⁵⁴ of the cases with and without bias (see Fig. S2 of the ESI†). At the same time, the effect of the bias on configurations can be understood in terms of a combination of increase/decrease in Pauli repulsion and small charging of the surface. Fig. 5(a and b) show the differences in charge density for different biases compared to the zero-bias case for the “flat” configuration. The corresponding insets of Fig. 5 indicate that most of the change in charge on the molecule is located on the oxygen. That transfer is small, however, as seen in both the insets and Fig. 5(d) which shows the change in charge density averaged over planes perpendicular to the surface. At the same time, by calculating the fluctuations in the density differences between the positive and negative bias,

ΣΔρ = Δρ_1.5,0 + Δρ_−1.5,0 (7)

= (ρ_V=1.5 − ρ_V=0) + (ρ_V=−1.5 − ρ_V=0), (8)

presented in Fig. 5(c), we notice that, although small (the value of the isosurface is 5.4 × 10⁻⁵ e per ų), it is asymmetric. Furthermore the final values of the water molecule dipole moments are similar to those of the isolated molecule in the cases with and without bias, and those dipoles tend to align with the field. The estimated magnitude of the dipole moments as a function of the applied bias are presented in the table of Fig. 4 together with the angular deviation from the dipole of an isolated molecule with the same orientation.

Fig. 5 Flat configuration. Relative changes in the charge density in comparison to the zero bias, Δρ_V,0 = ρ_V − ρ₀ for (a) V = +1.5 V and (b) V = −1.5 V (isosurface values = ±8.1 × 10⁻⁴ e per ų). Insets show an xy-sectional plane taken at the water molecule center of mass. (c) Charge density fluctuation at positive and negative applied biases, i.e. ΣΔρ = Δρ_1.5,0 + Δρ_−1.5,0 (isosurface value = ±5.4 × 10⁻⁵ e per ų). In all of the cases, red (blue) indicates an excess (deficiency) of electrons. (d) Laterally averaged difference in the charge density in comparison to the zero bias: black for +1.5 V and red for −1.5 V; the water–metal system indicates the position of the atoms.

Thus, the picture that arises is the following: for the positive bias the left hand side surface becomes negatively charged and tends to repel the negatively charged oxygen. At the same time as a small amount of charge is transferred to the oxygen, the overlap between the oxygen orbitals and gold surface orbitals tends to move the molecule away due to Pauli repulsion. The opposite behavior is expected for the negative bias. This last point is evidenced in Fig. 6 where we show the difference between our Au–water system and the charge density in an isolated capacitor with a ±1.5 V bias applied and an isolated water molecule in an equivalent electric field. In doing this we remove effects from the charge density that would arise solely from the electric field, focusing instead on the effects due to the water–surface interaction, namely the Pauli repulsion. For zero bias (Fig. 6(b)) the signature of Pauli repulsion, namely the “pillow” density of states between the molecule and metallic surface, is already visible.^55,56 For the positive bias, the interaction between molecule and surface increases and the nodes in the density disappear, an indication of decreased Pauli repulsion. The interaction is thus more attractive. On the other hand, for the negative bias there is larger repulsion as indicated by a slightly larger gap between the “pillow” region and the charge density associated with the molecule.

Fig. 6 Differences in the charge density between the Au + water system, a parallel plate capacitor and an isolated H₂O molecule submitted to an equivalent electric field for (a) V = 1.5 V, (b) V = 0, and (c) V = −1.5 V. An isosurface value of ±8.4 × 10⁻³ e per ų was considered in all of the plots, where red (blue) indicates an excess (deficiency) of electrons.

4 Conclusions

In conclusion, the inclusion of electronic effects via DFT in the description of water–metal interactions is important to advance the comprehension of the local structure of water at an electrochemical interface. In this work we showed how DFT combined with NEGF can be used to describe water–metal systems under an external bias potential. This framework allows for a description of a truly semi-infinite metallic electrode that sets a reference chemical potential that can be controlled by applying an external bias without adding/removing additional charge to the system. This allows for a more direct comparison with the experimental setups. This methodology now allows one to properly calculate the forces and therefore perform relaxation or dynamics of water–metal systems out of equilibrium, simulating an electrochemical cell in the sense that ionic currents could be considered now. This can be achieved by performing a molecular dynamics simulation of the system with the bias applied. In principle, our model allows one to consider a more realistic electrochemical cell with a thicker region of water molecules, although it will be computationally more expensive. One possible way to increase the system size in a cheaper way is to consider only the double layer region at the DFT level and the other water molecules at a molecular mechanics level (QM/MM method).^57,58

Here, we presented how the magnitude and sign of the bias alters the interaction of a prototype system of one water molecule on top of an Au(111) surface. The external bias changes both the position and alignment of the molecule with the surface. In particular, we have showed that a small positive bias leads to an unbound water molecule on a gold surface due to a combination of electrostatic effects and Pauli repulsion. On the other hand, a negative bias increases the oxygen–metal bond, and leads to a slightly rotated water molecule, thus indicating that the introduction of an external bias has a significant influence on the microscopic structure of molecules at a metallic interface.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The work was funded by DOE Early Career Awards No. DE-SC0003871 and DE-FG02-09ER16052. L. S. P. and A. R. R. also acknowledge financial support from ICTP-SAIFR (FAPESP project No. 2011/11973-4) and the ICTP-Simons Associate Scheme. P. B. acknowledges financial support from the FP7 FET-ICT “Planar Atomic and Molecular Scale devices” (PAMS) project (funded by the European Commission under contract No. 610446), the Spanish Agencia Estatal de Investigación (Grant No. MAT2016-78293-C6-4-R) and the Dep. de Educación of the Basque Government and UPV/EHU (Grant No. IT-756-13). This research used computational resources at the Center for Functional Nanomaterials, Brookhaven National Laboratory, which is supported by the US Department of Energy under Contract No. DE-AC02-98CH10886 and the S. Dumont supercomputer at the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil).

References

1. J. K. Norskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Nat. Chem., 2009, 1, 37.
2. J. Greeley and N. M. Markovic, Energy Environ. Sci., 2012, 5, 9246.
3. A. Hodgson and S. Haq, Surf. Sci. Rep., 2009, 64, 381.
4. P. A. Thiel and T. E. Madley, Surf. Sci. Rep., 1987, 7, 211.
5. O. Björneholm, M. H. Hansen, A. Hodgson, L.-M. Liu, D. T. Limmer, A. Michaelides, P. Pedevilla, J. Rossmeisl, H. Shen, G. Tocci, E. Tyrode, M.-M. Walz, J. Werner and H. Bluhm, Chem. Rev., 2016, 116, 7698.
6. M. F. Toney, J. N. Howard, J. Richer, G. L. Borges, J. G. Gordon, O. R. Meiroy, D. G. Wiesler, D. Yee and L. B. Sorensen, Nature, 1994, 368, 444.
7. J.-J. Velasco-Velez, C. H. Wu, T. A. Pascal, L. F. Wan, J. Guo, D. Prendergast and M. Salmeron, Science, 2014, 346, 831.
8. K. Ataka, T. Yotsuyanagi and M. Osawa, J. Phys. Chem., 1996, 100, 10664.
9. M. Nielsen, M. E. Björketun, M. H. Hansen and J. Rossmeisl, Surf. Sci., 2015, 631, 2.
10. F. Calle-Vallejo and M. T. M. Koper, Electrochim. Acta, 2012, 84, 3.
11. J. Cheng and M. Sprik, Phys. Chem. Chem. Phys., 2012, 14, 11245.
12. J. I. Siepmann and M. Sprik, J. Chem. Phys., 1995, 102, 511.
13. A. P. Willard, S. K. Reed, P. A. Madden and D. Chandler, Faraday Discuss., 2009, 141, 423.
14. M. K. Petersen, R. Kumar, H. S. White and G. A. Voth, J. Phys. Chem. C, 2012, 116, 4903.
15. C. D. Taylor, S. A. Wasileski, J.-S. Filhol and M. Neurock, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 165402.
16. J.-S. Filhol and M. Neurock, Angew. Chem., Int. Ed., 2006, 45, 402.
17. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, 864.
18. W. Kohn and L. J. Sham, Phys. Rev. A, 1965, 140, 1133.
19. N. Bonnet, T. Morishita, O. Sugino and M. Otani, Phys. Rev. Lett., 2012, 109, 266101.
20. M. Otani and O. Sugino, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 115407.
21. M. Mamattkulov and J.-S. Filhol, Phys. Chem. Chem. Phys., 2011, 13, 7675.
22. S. Datta, Electronic Transport in Mesoscopic Systems, Cambridge University Press, 1997.
23. M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor and K. Stokbro, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 165401.
24. A. R. Rocha, V. M. García-Suárez, S. Bailey, C. Lambert, J. Ferrer and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 085414.
25. A. R. Rocha, V. M. García-Suárez, S. W. Bailey, C. J. Lambert, J. Ferrer and S. Sanvito, Nat. Mater., 2005, 4, 335.
26. C. Caroli, R. Combescot, P. Nozieres and D. Saint-James, J. Phys. C: Solid State Phys., 1971, 4, 916.
27. H. Hellmann, Einführung in die Quantenchemie, F. Deuticke, Leipzig, 1937, vol. 54.
28. R. P. Feynman, Phys. Rev., 1939, 56, 340.
29. R. Zhang, I. Rungger, S. Sanvito and S. Hou, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 085445.
30. M. Di Ventra and S. T. Pantelides, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 61, 16207.
31. T. N. Todorov, D. Dundas and E. J. McEniry, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 075416.
32. M. Bai, C. S. Cucinotta, Z. Jiang, H. Wang, Y. Wang, I. Rungger, S. Sanvito and S. Hou, Phys. Rev. B, 2016, 94, 035411.
33. J.-T. Lu, J.-S. Wang, P. Hedegard and M. Brandbyge, Phys. Rev. B, 2016, 93, 205404.
34. J. B. Neaton, M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett., 2006, 97, 216405.
35. J. M. Garcia-Lastra and K. S. Thygesen, Phys. Rev. Lett., 2011, 106, 187402.
36. C. Toher, A. Filippetti, S. Sanvito and K. Burke, Phys. Rev. Lett., 2005, 95, 146402.
37. M. Strange, C. Rostgaard, H. Häkkinen and K. S. Thygesen, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 115108.
38. C. D. Pemmaraju, I. Rungger, X. Chen, A. R. Rocha and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 125426.
39. A. M. Souza, I. Rungger, C. D. Pemmaraju, U. Schwingenschloegl and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 165112.
40. A. d. M. Souza, I. Rungger, R. B. Pontes, A. R. Rocha, A. J. R. da Silva, U. Schwingenschloegl and S. Sanvito, Nanoscale, 2014, 6, 14495–14507.
41. P. Ordejón, E. Artacho and J. M. Soler, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 10441.
42. J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón and D. Sánchez-Portal, J. Phys.: Condens. Matter, 2002, 14, 2745.
43. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
44. M. Dion, H. Rydberg, E. Schroder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
45. J. Wand, G. Roman-Perez, J. Soler, E. Artacho and M. V. Fernandez-Serra, J. Chem. Phys., 2011, 134, 024516.
46. N. Troullier and J. L. Martins, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993.
47. A. Michaelides, V. A. Ranea, P. L. de Andres and D. A. King, Phys. Rev. Lett., 2003, 90, 216102.
48. G. Cicero, A. Calzolari, S. Corni and A. Catellani, J. Phys. Chem. Lett., 2011, 2, 2582.
49. A. Poissier, S. Ganeshan and M. V. Fernandez-Serra, Phys. Chem. Chem. Phys., 2011, 13, 3375.
50. J. Carrasco, A. Michaelides and M. Scheffler, J. Chem. Phys., 2009, 130, 184707.
51. J. Huzayyin, J. H. Chang, K. Lian and F. Dawson, J. Phys. Chem. C, 2014, 118, 3459.
52. Note that the positive bias is defined when we have net negative charge on the surface. Therefore, this situation corresponds to the negative electric field in ref. 51.
53. L. S. Pedroza, A. Poissier and M.-V. Fernández-Serra, J. Chem. Phys., 2015, 142, 034706.
54. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, Oxford, UK, 1994.
55. R. J. P. Keijsers, J. Voets, O. I. Shklyarevskii and H. van Kempen, Phys. Rev. Lett., 1996, 76, 1138–1141.
56. Weiss, C. Wagner, C. Kleimann, M. Rohlfing, F. S. Tautz and R. Temirov, Phys. Rev. Lett., 2010, 105, 086103.
57. G. T. Feliciano, C. Sanz-Navarro, M. D. Coutinho-Neto, P. Ordejón, R. H. Scheicher and A. R. Rocha, Phys. Rev. Appl., 2015, 3, 034003.
58. G. T. Feliciano, C. Sanz-Navarro, M. D. Coutinho-Neto, P. Ordejón, R. H. Scheicher and A. R. Rocha, J. Phys. Chem. B, 2017 DOI:10.1021/acs.jpcb.7b03475.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sc02208e

---