Introduction

A detailed understanding of the atomic and electronic structures of electrified surfaces and electrochemical interfaces has critical implications for the development of next-generation semiconductor devices as well as energy conversion and storage devices1,2. Toward this goal, significant advances have been recently made based on operando experimental characterization techniques and corresponding nonequilibrium first-principles simulations3,4,5. For the latter, however, including the electrode potential effects within first-principles simulations remains a significant challenge2. Due to the difficulty of incorporating the nonequilibrium constant electrochemical potential condition, which is the actual control parameter in experiments, the majority of works have adopted approaches that apply an electric field6,7,8 [Fig. 1a, ➀] or assign charges on metal (or effective metal) electrodes4,9,10,11,12,13 [Fig. 1a, ➁]. One promising way to faithfully simulate constant electrochemical potential conditions [Fig. 1a, ➂] is adopting the nonequilibrium Green’s function (NEGF) formalism combined with density functional theory (DFT), which has been established as the standard approach for first-principles finite-bias quantum transport calculations14,15. However, this method suffers from an ill-defined total energy, and accordingly, fully first-principles schemes for the constant-potential simulation of electrified interfaces need to be further developed.

Fig. 1: Formulation and validation of the electric enthalpy in electrified interface within the MS-DFT formalism.
figure 1

a To describe the electrode potential Φ-dependent interaction between a working electrode L and a molecule C with the dipole moment \({{\boldsymbol{p}}}_{C}^{{V}_{b}}\) (green solid arrow), rather than introducing an electric field ε0 across C (➀) or assigning charges on electrodes (➁), MS-DFT explicitly introduces the bias voltage \({V}_{b}\) between L and the reference electrode R, which are characterized by the Fermi functions f with the electrochemical potentials \({\mu }_{L}\) and \({\mu }_{R}\), respectively (➂). Electrode potential Φ is defined as −\({V}_{b}/2\). b The atomic structure of the Au(111)-single water-Au(111) interface model overlaid with the two-dimensional contour plot of the bias-induced electrostatic potential change \(\Delta {v}_{H}\) (top) and the corresponding plane-averaged charge density difference or Landauer resistivity dipole \(\Delta \bar{\rho }\) (bottom) obtained at Φ = −2.0 V. c The \({{\mathcal{F}}}_{C}^{{V}_{b}}\) (red filled diamonds) is shifted downward from \({E}_{L+C+R}^{{V}_{b}}-{E}_{L+R}^{{V}_{b}}\) (orange open diamonds) by the \({{\boldsymbol{\varepsilon }}}_{0}\cdot {{\boldsymbol{p}}}_{C}^{{V}_{b}}\) contribution. The corresponding \({-}\nabla {{\mathcal{F}}}_{C}^{{V}_{b}}\) (dotted lines) curves match well with the z-component of the atomic forces \({F}_{z}\) (gray dashed line).

Full size image

Herein, based on the multi-space constrained-search DFT (MS-DFT) approach, we develop a first-principles approach for the computation of the electric enthalpy of a molecular system in contact with an electrified electrode and, accordingly, the nonequilibrium electrode-molecule adsorption energy. The MS-DFT formalism is a microcanonical ensemble approach for nonequilibrium open quantum systems we recently established by adopting a viewpoint that maps quantum transport processes to multi energy- and real-space excitation counterparts16,17,18,19. Within MS-DFT, the nonequilibrium state of channel C established by left electrode L and right electrode R maintained at different electrochemical potentials \({\mu }_{L}\) and \({\mu }_{R}\), respectively, is calculated by variationally minimizing the total energy functional within the constraint of an applied finite bias voltage \({V}_{b}=({\mu }_{L}-{\mu }_{R})/e\). The possibility of defining the variationally convergent nonequilibrium total energy for a molecular system C sandwiched between the working electrode L and the reference electrode R could possibly open previously unavailable avenues for the investigation of electrified interfaces17. However, rather than the total energy of the entire L + C + R system, the appropriate thermodynamic potential that characterizes the electrified interface is the electric enthalpy of channel C. Accordingly, in this work, we extend MS-DFT and establish a method to calculate the electric enthalpy of a molecular system or the nonequilibrium electrode-molecule adsorption energy. Considering an electrified interface model in which a single water molecule is placed on the Au(111) and graphene electrodes in various conformations, we find that the two cases show very different trends in that the flat water conformation is stabilized on the Au electrode but not on the graphene counterpart. The origins of this drastic difference will be explained in terms of the configuration-dependent equilibrium electronic structures and their changes under finite electrode potential conditions.

Results

Formulation and validation of the electric enthalpy in electrified interface within MS-DFT

While the electric enthalpy was a well-established concept in the field of electromechanics20,21,22, it was adopted by the first-principles calculation community only in relatively recent years for the development of “modern theory of polarization” or density polarization functional theory23,24,25,26. These works are concerned with the channel-only system C within periodic boundary condition or the bulk state, and thus the electric enthalpy of C under the macroscopic field \({{\boldsymbol{\varepsilon }}}_{0}\) defined as

$${{\mathcal{H}}}_{C}^{{{\boldsymbol{\varepsilon }}}_{0}}={E}_{C}^{{{\boldsymbol{\varepsilon }}}_{0}}-{{\boldsymbol{\varepsilon }}}_{0}\,\cdot\,{{\boldsymbol{p}}}_{C}^{{{\boldsymbol{\varepsilon }}}_{0}},$$
(1)

where \({E}_{C}^{{{\boldsymbol{\varepsilon }}}_{0}}\) is the (internal) energy of C under \({{\boldsymbol{\varepsilon }}}_{0}\) and \({{\boldsymbol{p}}}_{C}^{{{\boldsymbol{\varepsilon }}}_{0}}\) is the electric dipole moment of C under \({{\boldsymbol{\varepsilon }}}_{0}\).

On the other hand, for the evaluation of the electric enthalpy for electrified interfaces, we need to explicitly include in addition to channel C the electrodes L + R that maintain the bias voltage \({V}_{b}=({\mu }_{L}-{\mu }_{R})/e\) and generate \({\boldsymbol{\varepsilon }}\)0. In this case, defining \({E}_{C}^{{V}_{b}}\) or \({{\mathcal{H}}}_{C}^{{V}_{b}}\) is problematic because the electronic structure of the channel C is determined via the quasi-Fermi levels or only through the explicit reference to L and R electrodes with the electrochemical potentials \({\mu }_{L}\) and \({\mu }_{R}\), respectively17,18. Instead, by performing an additional MS-DFT calculation for the biased L + R electrode-only system, we can define the electric enthalpy of C according to

$${{\mathcal{F}}}_{C}^{{V}_{b}}={E}_{L+C+R}^{{V}_{b}}-{E}_{L+R}^{{V}_{b}}-{{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}\,\cdot\,{{\boldsymbol{p}}}_{C}^{{V}_{b}},$$
(2)

or the non-equilibrium electrode-molecule adsorption energy of C as

$$\Delta {{\mathcal{F}}}_{C}^{{V}_{b}}={{\mathcal{F}}}_{C}^{{V}_{b}}-{E}_{C}^{0}={E}_{L+C+R}^{{V}_{b}}-{E}_{L+R}^{{V}_{b}}{\boldsymbol{-}}{E}_{C}^{0}-{{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}\,\cdot \,{{\boldsymbol{p}}}_{C}^{{V}_{b}},$$
(3)

where \({E}_{L+C+R}^{{V}_{b}}\) (\({E}_{L+R}^{{V}_{b}}\)) is the non-equilibrium total energy of the L + C + R (L + R) junction under \({V}_{b}\). Namely, compared with \({{\mathcal{H}}}_{C}^{{{\boldsymbol{\epsilon}}}_{0}}\) in Eq. (1), we defined \({{\mathcal{F}}}_{C}^{{V}_{b}}\) in Eq. (2) such that it includes the additional energy gain due to the formation of an interface between C and L + R under \({V}_{b}\) (see also Supplementary Note 1).

Note that while the nonequilibrium junction total energy \({E}_{L+C+R}^{{V}_{b}}\) was so far an ill-defined quantity, MS-DFT now provides a formal justification as well as a practical way to calculate it17. Once the MS-DFT calculation for the biased L+R electrode-only system is performed, we can straightforwardly evaluate \({{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}\) according to

$$-\nabla \left[{v}_{H,L+R}^{{V}_{b}}\left({\boldsymbol{r}}\right)-{v}_{H,L+R}^{0}\left({\boldsymbol{r}}\right)\right]={{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}({\boldsymbol{r}}),$$
(4)

where \({v}_{H}^{0}\) (\({v}_{H}^{{V}_{b}}\)) is the L+R equilibrium (nonequilibrium) Hartree electrostatic potential. The dipole moment \({{\boldsymbol{p}}}_{C}^{{V}_{b}}\) can be also obtained by extracting \({\rho }_{C}^{{V}_{b}}\left({\boldsymbol{r}}\right)={\rho }_{L+C+R}^{{V}_{b}}\left({\boldsymbol{r}}\right)-{\rho }_{L+R}^{{V}_{b}}\left({\boldsymbol{r}}\right)\) and calculating

$${{\boldsymbol{p}}}_{C}^{{V}_{b}}=\int d{\boldsymbol{r}}\,\left[{\rho }_{C}^{{V}_{b}}\left({\boldsymbol{r}}\right)-{\rho }_{{atm}}\left({\boldsymbol{r}}\right)\right]\,{{\boldsymbol{r}}}_{C},$$
(5)

where \({\rho }_{{atm}}\) is the summation of atomic valence charge densities.

Utilizing the Au(111)-single water-Au(111) junction model with the electrode-electrode distance of d = 20 Å (see Methods for details), we first demonstrate how the electronic and energetic properties of an interfacial water are calculated within MS-DFT under the applied working electrode potential defined as \(\varPhi =-{V}_{b}/2=-\left({\mu }_{L}-{\mu }_{R}\right)/2e\). For the given electrode−electrode distance d, the channel electric field ε0 is determined according to the electrode L potential \(\varPhi\) and will be specified together with \(\varPhi\). In Fig. 1b top panel, for the interfacial water in the ‘H-down’ configuration located at the O atom distance from the Au surface \(\overline{{OM}}\) of 3.5 Å, we show the contour plot of \(\Delta {v}_{H}\left({\boldsymbol{r}}\right)\equiv {v}_{H,L+C+R}^{{V}_{b}}\left({\boldsymbol{r}}\right)-{v}_{H,L+C+R}^{0}\left({\boldsymbol{r}}\right)\) calculated under the negative potential Φ = −2.0 V. We note that this potential drop corresponds to the electric field ε0 = −0.25 VÅ−1, which is stronger than −0.20 VÅ−1 that is expected from the device geometry of d = 20 Å. This inconsistency results from d being effectively reduced due to the finite atomic radius of surface Au atoms (1.74 Å27). We also show in the bottom panel of Fig. 1b the corresponding plane-averaged bias-induced electron density change or Landauer resistivity dipole \(\Delta \bar{\rho }\left({\boldsymbol{r}}\right)\equiv {\bar{\rho }}^{{V}_{b}}\left({\boldsymbol{r}}\right)-{\bar{\rho }}^{0}\left({\boldsymbol{r}}\right)\). The \(\varPhi\)-induced electron transfer from the negatively charged working electrode to the H atoms of the water molecule (dotted box in the Fig. 1b bottom panel) stabilizes the H-down configuration of water, and the bias-induced electrostatic potential drop and interfacial charge transfer behaviors calculated within MS-DFT are presented in Fig. 1b top and bottom panels, respectively. The resulting energetic stabilization is quantified in terms of the electric enthalpy \({{\mathcal{F}}}_{C}^{{V}_{b}}\) profile. In Fig. 1c, we show \({{\mathcal{F}}}_{C}^{{V}_{b}}\) and \({E}_{L+C+R}^{{V}_{b}}-{E}_{L+R}^{{V}_{b}}\) as a function of \(\overline{{OM}}\). Note that the difference the two represents the \({{\boldsymbol{-}}{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}\,\cdot\,{{\boldsymbol{p}}}_{C}^{{V}_{b}}\) term and the corresponding downshift of \({{\mathcal{F}}}_{C}^{{V}_{b}}\) curve from the \({E}_{L+C+R}^{{V}_{b}}-{E}_{L+R}^{{V}_{b}}\) counterpart indicates the stabilization of the H-down water configuration due to the same directions of \({{\boldsymbol{p}}}_{C}^{{V}_{b}}\) (negative sign in the H-down configuration) and \({{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}\) (negative sign at the negative Φ). Within the MS-DFT formalism, since the non-equilibrium Hamiltonian or total energy is a well-defined quantity, the forces on the nuclei I are straightforwardly calculated according to the Hellmann-Feynman theorem28,29

$${\vec{F}}_{I}=-\frac{\partial \left\langle {\psi }^{L/C/R}\left|\,\hat{H}\right|{\psi }^{L/C/R}\right\rangle }{\partial {\vec{R}}_{I}},$$
(6)

where the \({\psi }^{L/C/R}\) are the spatially-resolved L/C/R Kohn-Sham states that satisfy the bias constraint of \(e{V}_{b}={\mu }_{L}-{\mu }_{R}\). We then confirmed that the non-equilibrium atomic forces on the center of mass of the water molecule projected along the surface-normal z direction \({F}_{z}\) (gray dashed line) almost perfectly matches the derivatives of the \({{\mathcal{F}}}_{C}^{{V}_{b}}\)–\(\overline{{OM}}\) curve (blue filled circles), numerically validating the formulation of the electric enthalpy presented above.

Au-water adsorption energies dependent on the electrode potential and water geometry

We now examine in more detail the non-equilibrium atomic forces \({F}_{z}\) and adsorption energies \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) of a single water molecule on the electrified Au(111) surface as a function of the water geometry and bias voltage. In Fig. 2a‒c middle and right panels, for the three representative flat (Fig. 2a, left panel), H-up (Fig. 2b, left panel), and H-down (Fig. 2c, left panel) water orientations, we provide as functions of the distance \(\overline{{OM}}\) the z-component non-equilibrium atomic forces \({F}_{z}\) and the non-equilibrium adsorption energies \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\), respectively. While the atomic forces \({F}_{z}\), which are also available from DFT-NEGF calculations15 (see also Supplementary note 2 and Supplementary Fig. 1), are limited to the determination of stable non-equilibrium water geometries (Fig. 2, middle panels), the adsorption energies \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) provide quantitative information on the bias-dependent energetic stabilization/destabilization trends of different water configurations (Fig. 2, right panels). For the flat configuration that assumes the optimal \(\overline{{OM}}\) distance z = 3 Å in equilibrium (Fig. 2a), we observe that no significant Φ-induced changes in \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) are obtained. It can be understood that the \({{\boldsymbol{p}}}_{{H}_{2}O}^{{V}_{b}}\) of the flat configuration is almost perpendicular to the \({{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}\), and thus the \({{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}\,\cdot\,{{\boldsymbol{p}}}_{{H}_{2}O}^{{V}_{b}}\) contribution becomes negligible. On the other hand, in the case of H-up (Fig. 2b) and H-down (Fig. 2c) configurations, where the \({{\boldsymbol{p}}}_{{H}_{2}O}^{{V}_{b}}\) is parallel to the \({{\boldsymbol{\varepsilon }}}_{0}^{{V}_{b}}\), we observe that the \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) curves shift with Φ and induce the strongest \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) of −0.206 eV at Φ = +2.0 V for the H-up configuration and −0.239 eV at Φ = −2.0 V for the H-down configuration.

Fig. 2: Non-equilibrium forces and adsorption energies of a single water molecule on metal surfaces.
figure 2

For the (a) ‘flat’, (b) ‘H-up’, and (c) ‘H-down’ water configurations (left panels), the z-component atomic forces (middle panels) and adsorption energies \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) (right panels) are shown as functions of the distance. Red triangles, green circles, and blue squares are the data obtained at Φ = −2.0 V, 0.0 V, and 2.0 V, respectively. d Relaxed configurations of a single water molecule on the Au(111) surface at Φ = −2.0 V (left), 0.0 V (middle), and 2.0 V (right). The angle θ is defined as the angle of the water plane from the Au surface-normal direction. e The \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) variations as functions of θ at different electrode potentials.

Full size image

Further extending these analyses, we first performed non-equilibrium geometry optimizations and as shown in Fig. 2d obtained θ = 162.0° (~H-down) at Φ = −2.0 V and θ = 31.6° (~H-up) at Φ = +2.0 V, respectively. Subsequently, we prepared several Au-water interface configurations by varying θ from 10.0° to 170.0° and calculated their \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) at different \(\varPhi\) (\(-2.0{\rm{V}}\le \varPhi \le 2.0{\rm{V}}\)) and present the results in Fig. 2e. In equilibrium, the flat geometry is the energetically most stable conformation with the energy barriers of 0.033 eV and 0.062 eV for reorientations to θ = 162° and θ = 31.6°, respectively (gray star marks). As Φ is negatively increased, H-down conformations become stabilized: At Φ = −1.0 V, the configurations with θ = 114.3°~162.0° became energetically metastable with negligible energy differences, and at Φ = −1.5 V the θ = 162° configuration becomes the energetically most stable configuration. The adsorption energy at this configuration is \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}={-}0.185\) \({\rm{eV}}\), which is further increased to \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}={-}0.232\) eV at Φ = −2.0 V. Similarly, as Φ is positively increased, H-up conformations become stabilized, but, due to the larger reorientation energy barrier along this direction, the \(\theta =31.6\)° configuration becomes the most stable configuration only at Φ = +2.0 V with the smaller adsorption energy of \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}=\)−\(0.188\) eV.

Differences between electrified Au-water and graphene-water interfaces

In addition to the availability of the non-equilibrium total energy and electric enthalpy, another important advantage of MS-DFT is the possibility to faithfully model two-dimensional (2D) electrodes such as graphene16,19. Using the same Au reference electrode, we next studied the graphene working electrode-single water system and found that graphene and Au electrodes show very different behaviors in terms of the electrode potential-dependent stabilization of water configurations. We first present in Fig. 3a the extended data set of the water geometry- and electrode potential-dependent water adsorption energies for the Au electrode case. Here, the H-up, flat, and H-down water configurations were adopted from the zero-Φ geometry optimizations shown in Fig. 2d and \(\overline{{OM}}\) distances were varied between 2.4 Å and 5.0 Å with the 0.2 Å interval. As discussed earlier, they show that, while the H-down (blue triangles) and H-up (green circles) water configurations become stabilized at large negative and positive Φ values, respectively, the flat configuration is energetically most preferable in low electrode potential regimes (red squares). The corresponding data obtained for the graphene electrode case are presented in Fig. 3b, and we find that, in contrast to the Au electrode counterpart, as the electrode potential increases from negative to positive values the H-down configuration directly transitions to the H-up configuration without exhibiting the bias regime where the flat configuration is stabilized. We note that because of the usage of the Au(111) reference electrode there exists a built-in electric field of ~+0.06 VÅ−1 induced by the difference between graphene and Au work function values (4.5 eV and 5.4 eV, respectively). However, other than the slight mismatch between zero-Vb and zero-ε0 limits, the presence of a built-in electric field should not affect the important trend of the absence of electrode potential regimes where the flat water geometry is stabilized on the graphene electrode. We expect that this notable difference would be an important ingredient necessary for the full understanding of future experiments on interfacial water3,4,5.

Fig. 3: Comparison between electrode potential- and water configuration-dependent adsorption energies of a single water molecule on Au and graphene electrodes.
figure 3

The bias-dependent adsorption energies of a single water molecule \(\Delta {{\mathcal{F}}}_{{H}_{2}O}^{{V}_{b}}\) on the (a) Au and (b) graphene electrodes calculated for the H-down (blue triangles), H-up (green circles), and flat (red squares) configurations are presented as functions of the \(\overline{{OM}}\) distance. In each panel, the electric field corresponding to the electrode potential is presented, and the curves for energetically less stable geometries are shown translucent. Additionally, the energetically most stable water structures are displayed in the insets of the panels for the lowest electrode potential cases.

Full size image

Electronic structures of the water adsorbed on electrified Au and graphene surfaces

We next examine the non-equilibrium electronic structures of Au-water and graphene-water interface models and explain their drastic difference in terms of the stability of the flat water configuration discussed above. In Fig. 4a, for the zero-Φ optimized flat-configuration water on the Au electrode (top panels) with \(\overline{{OM}}=3.0\) Å and graphene electrode (bottom panels) with \(\overline{{OM}}=3.0\) Å (Fig. 3, red up arrows), we show the water molecule-projected density of states (PDOS) (left and center panels) and local DOS (LDOS) of the 3a1 O lone pair highest occupied molecular orbital (HOMO) and 1b2 OH bonding HOMO-1 states (right panels) obtained at Φ = −2.0 V and +2.0 V MS-DFT calculations. In Fig. 4b, we also show the corresponding data obtained for the H-down-configuration water on the Au electrode (top panels) with \(\overline{{OM}}=3.6\) Å and the graphene electrode (bottom panels) with \(\overline{{OM}}=3.2\) Å (Fig. 3, blue down arrows). We then find that in general water PDOS downshift (upshift) under negative (positive) Φ due to the electron transfer from (to) Au/graphene to (from) water (See also Supplementary Figs. 2 and 3). While this behavior is universally visible for the deeper HOMO−1 states (orange left triangles in the PDOS figures), we also observe that only in the flat water-Au electrode case (Fig. 4a, top panels) the HOMO (cyan left triangles in the PDOS figures) states are held nearly fixed because of their strong hybridizations with Au states. Importantly, such strong hybridization of the HOMO is missing for the flat-configuration water on graphene case (Fig. 4a, bottom panels), explaining why the flat water configuration is not stabilized on the graphene electrode.

Fig. 4: Electronic structures of the single water molecule on the electrified Au and graphene surfaces.
figure 4

The PDOS of water in the (a) flat and (b) H-down configurations on Au (top panels) and graphene (down panels) calculated at Φ = −2.0 V (left panes) and +2.0 V (central panels). For reference, equilibrium PDOS are shown together (shaded curves). Cyan and orange left triangles indicate the HOMO and HOMO−1 states of the water molecule, respectively. In each case, the corresponding equilibrium HOMO (top) and HOMO-1 (bottom) LDOS extracted from the energy window of [−0.1 eV, +0.1 eV] are shown on the right with the isosurface level of 0.01 states·Å−3eV−1.

Full size image

Comparison between MS-DFT and constant-electric field DFT calculations

We finally discuss the widely used constant-electric field DFT method in light of the MS-DFT formalism. In Fig. 5a, we show for the H-down configuration on Au (left) and graphene (right) electrode cases the plane-averaged \(\Delta {\bar{v}}_{{\rm{H}}}\) plot under Φ = −2.0 V. Note that, as shown in Fig. 5a bottom panels, the constant-\({\varepsilon_0}\) approach treats electrified interfaces by employing a half-cell model or adopting one global Fermi level. For the Au electrode case, the \(\Delta {\bar{v}}_{{\rm{H}}}\) drop is 4 eV, but, as mentioned above regarding Fig. 1b, the potential Φ = −2.0 V (Fig. 5a left panels, red lines) generates the electric field ε0 = −0.25 VÅ−1 (Fig. 5a left-bottom panel, blue dashed line), which is larger than ε0 = −0.20 VÅ−1 (Fig. 5a left-bottom panel, green dotted line) naively expected from the device geometry of d = 20 Å. Namely, within the constant-electric field method, one needs to determine ε0 based on an educational guess, highlighting the advantage of the MS-DFT formalism. On the other hand, for the graphene electrode case, the calculated magnitude of ∆vH = 3.6 eV (Fig. 5a right panels, red lines) is smaller than the applied electrochemical potential difference of 4 eV, and this mismatch is caused by the quantum capacitance effect in graphene19. Then, taking into account also the atomic radius of carbon atoms (0.67 Å27), the corresponding ε0 becomes −0.20 VÅ−1 (Fig. 5a right-bottom panel, green dotted line), rather than −0.25 VÅ−1 derived in the Au electrode case (Fig. 5a right-bottom panel, blue dashed line). In Fig. 5b, we compare the plane-averaged \(\varDelta \bar{\rho }\) curves (top panels) and 2D Δρ contour plots (bottom panels) obtained from MS-DFT and constant-\({\varepsilon_0}\) calculations, respectively. We can then confirm that as long as the correct electric field is adopted (ε0 = −0.25 VÅ−1 in the Au electrode case and \({\varepsilon }\)0 = −0.20 VÅ−1 in the graphene electrode case), the constant-electric field method reproduces well the MS-DFT calculation results that include the non-equilibrium adsorption energies (see Supplementary Fig. 4). Accordingly, MS-DFT provides a formal justification for adopting the constant-electric field method, but we emphasize that \({\varepsilon }\)0 needs to be carefully determined within this approach.

Fig. 5: Justification of the constant-ε method based on MS-DFT.
figure 5

a (upper panels) The plane-averaged electrostatic potential differences \(\Delta {\bar{v}}_{H}\) calculated for the H-down water molecule on the Au (left) and graphene (right) under Φ = −2.0 V within MS-DFT. (lower panels) The \(\Delta {\bar{v}}_{H}\) curves calculated within MS-DFT at Φ = −2.0 V (red solid lines) are compared with those obtained from the constant-electric field method at ε0=−0.25 VÅ−1 (blue dotted lines) and −0.20 VÅ−1 (green dotted lines). b The plane-averaged charge density differences \(\Delta \bar{\rho }\) calculated within MS-DFT at Φ = −2.0 V (red solid lines) and constant-electric field method at ε0 = −0.25 VÅ−1 (blue dashed lines) and −0.20 VÅ−1 (green dotted lines) for the H-down water on Au (left panel) and graphene (right panel). The 2D contour plots of Δρ obtained from MS-DFT calculations at Φ = −2.0 V and the corresponding constant-\({\varepsilon_0}\) DFT calculations (ε0 = −0.25 VÅ−1 for Au and −0.20 VÅ−1 for graphene).

Full size image

Discussion

In summary, based on the MS-DFT framework, we formulated a first-principles theory to calculate the electric enthalpy of electrified interfaces or the non-equilibrium electrode-molecule adsorption energy. Specifically, we demonstrated that the internal energy of the non-equilibrium channel C can be obtained by performing an MS-DFT calculation for the electrodes-only L + R counterpart in addition to that for the full L + C + R junction and identifying the difference between the two total energies. A subsequent DFT calculation for the channel C under the electric field extracted from the L + R electrodes-only MS-DFT calculation then provides the channel C polarization, from which one can determine the electric enthalpy of the channel C and by referencing to the equilibrium channel C total energy the non-equilibrium electrode L-channel C adsorption energy. Applying the developed theory to the Au(111)-single water and graphene-single water interface models, we calculated the electrode potential- and water configuration-dependent adsorption energies and found that Au and graphene electrodes show very different trends: On the Au electrode, the flat-water configuration was found to be the energetically most preferable conformation in the low electrode potential regimes and the H-down and H-up water configurations are stabilized only at large negative and positive electrode potentials, respectively. On the other hand, on the graphene electrode, the flat water-configuration was not stabilized and the direct transition from the H-down to H-up water configuration occurred as the potential sweeps from the negative to positive direction. Analyzing the flat-water configuration electronic structures in the two cases, we finally revealed that the origins of this difference is due to strong (weak) hybridization of the water O lone pair orbitals with the Au (graphene) states and their energetic pinning (shift) with the applied potentials. The non-equilibrium first-principles calculation theory presented here should prove to be a valuable tool in advancing the atomic-level understanding of electrified interfaces, and in the follow-up publications we will report on the simulations on more realistic electrochemical interfaces that are relevant to the development of various energy and electronic devices.

Methods

Computational details

We implemented the developed theory within the SIESTA package30 in combination with the in-house grid operation code31,32. We applied the MS-DFT formalism to an electrode-water interface model where the Au(111) or graphene working electrode L is in contact with a single water molecule C and is balanced by the reference electrode Au(111) R, as shown in Fig. 1b. To avoid the interactions between neighboring images of the water molecule, a supercell with dimensions of 9.993 Å × 8.655 Å for the Au electrode and 12.918 Å × 12.430 Å for the graphene electrode along the surface-parallel direction was used. The distance between L and R electrodes was set to 20 Å. According to a previous study8, the electric field of ~±0.4 VÅ−1 corresponds to the electrode potential of ~±0.3 V with respect to the potential of zero charge of Au(111), making the voltage drops within our computational model reasonably match the electric fields of 0.1−1.0 VÅ−1 within the Helmholtz layer1. The centroid of the electric dipole moment of the water molecule \({{\boldsymbol{p}}}_{{H}_{2}O}^{V}\) was set as the center of mass of the water molecule. All calculations were performed within the van der Waals density functional of optB8833 to obtain accurate electrode-water interfacial properties. Double ζ-plus-polarization-level numerical atomic orbital basis sets were employed together with Troullier-Martins-type norm-conserving pseudopotentials34. A mesh cutoff of 600 Ry for the real-space integration was used, and 5 × 5 × 1 and 8 × 8 × 1 Monkhorst-Pack k-point grids of the Brillouin zone35 were sampled for Au and graphene, respectively. The atomic geometries were optimized until the total residual forces were below 0.005 eVÅ−1 using the conjugate gradient algorithm.