Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Luana S.
Pedroza
*^{ab},
Pedro
Brandimarte
^{cd},
Alexandre Reily
Rocha
^{e} and
M.-V.
Fernández-Serra
^{fg}
^{a}ICTP South American Institute for Fundamental Research, Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo SP 01140-070, Brazil. E-mail: l.pedroza@ufabc.edu.br
^{b}Centro de Ciências Naturais e Humanas, Universidade Federal do ABC, Santo André, São Paulo, Brazil 09210-170
^{c}Centro de Física de Materiales, 20018 Donostia – San Sebastián, Gipuzkoa, Spain
^{d}Donostia International Physics Center, 20018 Donostia – San Sebastián, Gipuzkoa, Spain
^{e}Instituto de Física Teórica, Universidade Estadual Paulista, São Paulo SP 01140-070, Brazil
^{f}Department of Physics and Astronomy, Stony Brook University, Stony Brook, New York 11794-3800, USA
^{g}Institute for Advanced Computational Sciences, Stony Brook University, Stony Brook, New York 11794-3800, USA

Received
16th May 2017
, Accepted 3rd October 2017

First published on 11th October 2017

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.

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.^{9} 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^{19} following a previous work where the charges at the surface were controlled by including a medium with a given permittivity.^{20} 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^{21} 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.

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:

(1) |

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:

(2) |

ρ^{L/R} = G(E)Γ_{L/R}(E)G^{†}(E) | (3) |

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:^{29} 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:

(4) |

(5) |

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

(6) |

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.

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.

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_{1} 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,^{53} 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.

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_{2}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,^{53} as seen from a Bader analysis^{54} 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) |

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.

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.

- J. K. Norskov, T. Bligaard, J. Rossmeisl and C. H. Christensen, Nat. Chem., 2009, 1, 37 CrossRef CAS PubMed.
- J. Greeley and N. M. Markovic, Energy Environ. Sci., 2012, 5, 9246 CAS.
- A. Hodgson and S. Haq, Surf. Sci. Rep., 2009, 64, 381 CrossRef CAS.
- P. A. Thiel and T. E. Madley, Surf. Sci. Rep., 1987, 7, 211 CrossRef CAS.
- 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 CrossRef PubMed.
- 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 CrossRef CAS.
- J.-J. Velasco-Velez, C. H. Wu, T. A. Pascal, L. F. Wan, J. Guo, D. Prendergast and M. Salmeron, Science, 2014, 346, 831 CrossRef CAS PubMed.
- K. Ataka, T. Yotsuyanagi and M. Osawa, J. Phys. Chem., 1996, 100, 10664 CrossRef CAS.
- M. Nielsen, M. E. Björketun, M. H. Hansen and J. Rossmeisl, Surf. Sci., 2015, 631, 2 CrossRef CAS.
- F. Calle-Vallejo and M. T. M. Koper, Electrochim. Acta, 2012, 84, 3 CrossRef CAS.
- J. Cheng and M. Sprik, Phys. Chem. Chem. Phys., 2012, 14, 11245 RSC.
- J. I. Siepmann and M. Sprik, J. Chem. Phys., 1995, 102, 511 CrossRef CAS.
- A. P. Willard, S. K. Reed, P. A. Madden and D. Chandler, Faraday Discuss., 2009, 141, 423 RSC.
- M. K. Petersen, R. Kumar, H. S. White and G. A. Voth, J. Phys. Chem. C, 2012, 116, 4903 CAS.
- C. D. Taylor, S. A. Wasileski, J.-S. Filhol and M. Neurock, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 165402 CrossRef.
- J.-S. Filhol and M. Neurock, Angew. Chem., Int. Ed., 2006, 45, 402 CrossRef CAS PubMed.
- P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, 864 CrossRef.
- W. Kohn and L. J. Sham, Phys. Rev. A, 1965, 140, 1133 CrossRef.
- N. Bonnet, T. Morishita, O. Sugino and M. Otani, Phys. Rev. Lett., 2012, 109, 266101 CrossRef PubMed.
- M. Otani and O. Sugino, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 115407 CrossRef.
- M. Mamattkulov and J.-S. Filhol, Phys. Chem. Chem. Phys., 2011, 13, 7675 RSC.
- S. Datta, Electronic Transport in Mesoscopic Systems, Cambridge University Press, 1997 Search PubMed.
- M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor and K. Stokbro, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 165401 CrossRef.
- 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 CrossRef.
- 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 CrossRef CAS PubMed.
- C. Caroli, R. Combescot, P. Nozieres and D. Saint-James, J. Phys. C: Solid State Phys., 1971, 4, 916 CrossRef.
- H. Hellmann, Einführung in die Quantenchemie, F. Deuticke, Leipzig, 1937, vol. 54 Search PubMed.
- R. P. Feynman, Phys. Rev., 1939, 56, 340 CrossRef CAS.
- R. Zhang, I. Rungger, S. Sanvito and S. Hou, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 085445 CrossRef.
- M. Di Ventra and S. T. Pantelides, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 61, 16207 CrossRef CAS.
- T. N. Todorov, D. Dundas and E. J. McEniry, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 075416 CrossRef.
- M. Bai, C. S. Cucinotta, Z. Jiang, H. Wang, Y. Wang, I. Rungger, S. Sanvito and S. Hou, Phys. Rev. B, 2016, 94, 035411 CrossRef.
- J.-T. Lu, J.-S. Wang, P. Hedegard and M. Brandbyge, Phys. Rev. B, 2016, 93, 205404 CrossRef.
- J. B. Neaton, M. S. Hybertsen and S. G. Louie, Phys. Rev. Lett., 2006, 97, 216405 CrossRef CAS PubMed.
- J. M. Garcia-Lastra and K. S. Thygesen, Phys. Rev. Lett., 2011, 106, 187402 CrossRef CAS PubMed.
- C. Toher, A. Filippetti, S. Sanvito and K. Burke, Phys. Rev. Lett., 2005, 95, 146402 CrossRef CAS PubMed.
- M. Strange, C. Rostgaard, H. Häkkinen and K. S. Thygesen, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 115108 CrossRef.
- C. D. Pemmaraju, I. Rungger, X. Chen, A. R. Rocha and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 125426 CrossRef.
- A. M. Souza, I. Rungger, C. D. Pemmaraju, U. Schwingenschloegl and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 165112 CrossRef.
- 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 RSC.
- P. Ordejón, E. Artacho and J. M. Soler, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 10441 CrossRef.
- 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 CrossRef CAS.
- J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865 CrossRef CAS PubMed.
- M. Dion, H. Rydberg, E. Schroder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401 CrossRef CAS PubMed.
- J. Wand, G. Roman-Perez, J. Soler, E. Artacho and M. V. Fernandez-Serra, J. Chem. Phys., 2011, 134, 024516 CrossRef PubMed.
- N. Troullier and J. L. Martins, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993 CrossRef CAS.
- A. Michaelides, V. A. Ranea, P. L. de Andres and D. A. King, Phys. Rev. Lett., 2003, 90, 216102 CrossRef CAS PubMed.
- G. Cicero, A. Calzolari, S. Corni and A. Catellani, J. Phys. Chem. Lett., 2011, 2, 2582 CrossRef CAS.
- A. Poissier, S. Ganeshan and M. V. Fernandez-Serra, Phys. Chem. Chem. Phys., 2011, 13, 3375 RSC.
- J. Carrasco, A. Michaelides and M. Scheffler, J. Chem. Phys., 2009, 130, 184707 CrossRef PubMed.
- J. Huzayyin, J. H. Chang, K. Lian and F. Dawson, J. Phys. Chem. C, 2014, 118, 3459 Search PubMed.
- 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.
- L. S. Pedroza, A. Poissier and M.-V. Fernández-Serra, J. Chem. Phys., 2015, 142, 034706 CrossRef PubMed.
- R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, Oxford, UK, 1994 Search PubMed.
- R. J. P. Keijsers, J. Voets, O. I. Shklyarevskii and H. van Kempen, Phys. Rev. Lett., 1996, 76, 1138–1141 CrossRef PubMed.
- Weiss, C. Wagner, C. Kleimann, M. Rohlfing, F. S. Tautz and R. Temirov, Phys. Rev. Lett., 2010, 105, 086103 CrossRef CAS PubMed.
- 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 CrossRef.
- 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 |

This journal is © The Royal Society of Chemistry 2018 |