The interactions between molecules and electrode surfaces play a key role in electrochemical processes and are a subject of extensive research, both experimental and theoretical. In this paper, we address the water dissociation reaction on a Pd(111) electrode surface, modeled as a slab embedded in an external electric field. We aim at unraveling the relationship between surface charge and zero-point energy in aiding or hindering this reaction. We calculate the energy barriers with dispersion-corrected density-functional theory and an efficient parallel implementation of the nudged-elastic-band method. We show that the lowest dissociation barrier and consequently the highest reaction rate take place when the field reaches a strength where two different geometries of the water molecule in the reactant state are equally stable. The zero-point energy contributions to this reaction, on the other hand, remain nearly constant across a wide range of electric field strengths, despite significant changes in the reactant state. Interestingly, we show that the application of electric fields that induce a negative charge on the surface can make nuclear tunneling more significant for these reactions.

Surface-mediated reactions involving water appear in many technologically relevant processes that lead to hydrogen production or material corrosion,1 for example. In particular, in electrocatalysis—a key process to achieve a sustainable economy—the chemistry involved in the water splitting reaction can be modified by the application of a potential bias.2 In such a system, the structure of water at the interface, the collective properties of the surface, and the electrolyte composition affect the reactions substantially and present a challenge for both theory and experiments, even when considering idealized electrodes.3–5 

From a theoretical perspective, there are at least three ingredients that must be accounted for when addressing water splitting at (charged) metallic surfaces: charge transfer and reorganization at the interface, (screened) van der Waals interactions between molecules and metal, and the effect of nuclear fluctuations and dynamics on the process of interest. It should be noted that the latter is now accepted to be adequately described only by considering the nuclei as quantum particles.6–11 While the structure of water at metallic surfaces has been addressed by theoretical models numerous times in the past,12–14 even including potential biases in the simulations,14–17 the surface-mediated dissociation of water has received comparatively less attention from theoretical work despite its undeniable importance. For example, under alkaline conditions, the hydrogen dissociation of a water molecule (the Volmer step) is the first step of the hydrogen evolution reaction and is likely to be rate-determining under these conditions.18 Several groups reported studies of the splitting of water molecules on various transition metals,19–25 but to the best of our knowledge, they never accounted, at the same time, for potential biases, long-range dispersion interactions, and nuclear quantum effects (NQEs).

Specifically, accounting for a potential bias in first-principles simulations is not trivial. A number of methods to include a potential bias in the slab model were proposed, including various degrees of approximation,26–33 without the adoption of a common “default” approach. A fully ab initio approach based on the non-equilibrium Green’s function formalism34 presents a very high computational hurdle for the simulation of rare reactive events. Other recent approaches, such as potentiostats that can be coupled to ab initio simulations31 or grand-canonical self-consistent field (SCF) techniques,30 could provide a good compromise between efficiency and accuracy but still present an elevated computational cost. Effective 2D-periodic models, such as the effective screening medium (ESM),26 have allowed the simulation of ab initio molecular dynamics of water at charged metal interfaces17,35 more than 15 years ago. With the ESM, the water dissociation on a platinum surface under a potential bias could be studied by the ab initio molecular dynamics (AIMD),25,36 which presents a remarkable achievement. However, only relatively short molecular dynamics (MD) trajectories (of a few picoseconds) could be simulated, which considerably limits the predictive power of such simulations, given that the reaction is a rare event on this timescale. Interestingly, it has recently been shown from classical MD simulations that the results obtained by simulations with a constant applied potential between two electrodes in a 2D-periodic model were indistinguishable from those obtained from a 3D-periodic model with a large metallic slab at the center of the simulation cell and an applied electric field37—a possibility that can also simplify the ab initio simulations of electrochemical interfaces.

In this paper, we employ an applied electric-field setup to address the water dissociation reaction directly and understand the role of different surface charges and nuclear quantum effects in modifying the reaction dynamics. We study a model system of a water monomer adsorbed on a Pd(111) surface, modeled by the density-functional theory. Different charges on the surface of the metal are realized by the application of an electric field on a thick Pd(111) slab, which mimics a potential bias. We report the changes in reaction barriers and rates and the effects of zero-point energy at different field strengths. These results serve as the groundwork for more complex system setups and simulation protocols. They also give valuable insights into how far the NQEs and electric fields can impact this reaction.

We simulated a Pd(111) surface using a slab of seven atomic layers in a periodic cell, including a large vacuum layer of 64 Å, using the FHI-aims code.40 A dipole correction41 was applied to compensate for the spurious interactions between periodic images. We employed the PBE exchange-correlation functional42 augmented by dispersion interactions as described by the screened vdWsurf model,43 with the coefficients for Pd taken from Ref. 44. Default light settings of FHI-aims were used for basis sets, and modified light settings with double radial density were used for numerical real-space grids. The lateral dimensions of the Pd slab were chosen as 4 × 4 unit cells, simulated with a 3 × 3 k-point mesh. We placed the metallic slab at the center of the unit cell.

The potential bias is mimicked by an external electric field permeating the system, which is implemented as a saw-tooth potential with the potential step in the vacuum region (see Fig. 1). The polarization induced by the electric field causes the surfaces of the slab to effectively become capacitor plates with opposite charges. We ensured that the slab is thick enough to screen the surfaces from each other, by analyzing the charge density distribution at self-consistency under different field strengths. Of course, the surfaces remain coupled by the charge conservation constraint. We observed a linear dependence of the surface charge on the applied electric field in a range from −10 to +10 V/Å. An increase in the field by 1 V/Å induced a surface charge of 0.0364 electrons per Pd atom, which is equivalent to 8.734 C/cm2. The surfaces at zero electric field are slightly negatively charged, which offsets the linear dependence by −0.012 e per atom. We will call a field “negative” or “positive” according to the charge that it induces on the surface of interest.

FIG. 1.

Schematic of the simulation setup. A periodic cell is permeated by an external electrostatic potential φext with a discontinuity in vacuum far away from the slab.

FIG. 1.

Schematic of the simulation setup. A periodic cell is permeated by an external electrostatic potential φext with a discontinuity in vacuum far away from the slab.

Close modal

The adsorption geometries of a water molecule on a Pd(111) surface were obtained by relaxation from the high-symmetry sites of the surface: top, bridge, and fcc/hcp hollow sites. If more than one geometry was found after the relaxation, the lowest-energy structure was chosen.

For the phonon calculations, a modified version45 of the Phonopy package46 was used and either only the water molecule or the molecule and the first layer of the surface were included to build the Hessian. The deeper layers of the surface were included as a rigid environment. By doing so, we lose the coupling Hessian elements between the atoms of the molecule and the surface beyond the first layer, which can affect the low-frequency molecule-surface phonons. Since the molecule is small and light, we assume that long-wave collective vibrations of a surface do not play a significant role in water reactions because the chemically relevant frequencies of the H–O–H bending and O–H stretching are not substantially affected.

The reaction paths were obtained using the climbing-image nudged elastic band (CI-NEB) method.47 We used a new implementation of the CI-NEB in the i-PI48 code, in which several instances of the FHI-aims code could connect simultaneously to the i-PI and efficiently calculate the forces on each replica in parallel. This implementation is described in detail in Ref. 49 and is available through the main repository of the code. We used a NEB path with nine intermediate replicas with the spring strength of 20–40 eV/Å. We used the FIRE50 algorithm for optimization. The tolerance for forces was set to 0.05 eV/Å for the optimization of the NEB path and to 0.01 eV/Å for the climbing-image optimization.

We have found the preferred adsorption geometries of a single water molecule at the Pd(111) surface for different electric field strengths and report them in Fig. 2.

FIG. 2.

Adsorption geometries of a single water molecule on the Pd(111) surface under −0.44 V/Å (left), no field (center), and +0.44 V/Å (right). The large green, white, and blue spheres denote the 1st, 2nd, and 3rd layers of the Pd(111) surface, respectively.

FIG. 2.

Adsorption geometries of a single water molecule on the Pd(111) surface under −0.44 V/Å (left), no field (center), and +0.44 V/Å (right). The large green, white, and blue spheres denote the 1st, 2nd, and 3rd layers of the Pd(111) surface, respectively.

Close modal

We found that the water molecule prefers a flat orientation if no field is applied, which agrees with previous studies.19–21 In a negative electric field, the flat orientation remains up to −0.37 V/Å, after which the molecule “stands up,” pointing both hydrogen atoms toward the surface. In a positive field, the molecule deviates gradually from the flat geometry and hydrogen atoms turn more and more toward the vacuum. These observations are also in qualitative agreement with the earlier work, for example, regarding a water monomer on Au(111) treated with non-equilibrium Green’s functions.34 We did not observe any abrupt changes in the range between −0.37 and +1 V/Å, as evidenced by the angles between the plane defined by the water molecule and the surface plane, reported in Fig. 3(a).

FIG. 3.

a) Empty black circles, left-hand axis: adsorption energy of a single water molecule on the Pd(111) surface as a function of the applied electric field. The values are calculated with PBE + vdWsurf. Green diamonds, right–hand axis: angle θ between the water molecule and the Pd surface (see the inset; the hydrogen atoms are equivalent). b) Charge induced on the adsorbed water molecule depending on the applied electric field, measured by the Mulliken (filled points) and Hirshfeld (empty points) analysis. The black triangles show the total charge on the molecule, and the blue squares (red circles) show H (O) contributions to it.

FIG. 3.

a) Empty black circles, left-hand axis: adsorption energy of a single water molecule on the Pd(111) surface as a function of the applied electric field. The values are calculated with PBE + vdWsurf. Green diamonds, right–hand axis: angle θ between the water molecule and the Pd surface (see the inset; the hydrogen atoms are equivalent). b) Charge induced on the adsorbed water molecule depending on the applied electric field, measured by the Mulliken (filled points) and Hirshfeld (empty points) analysis. The black triangles show the total charge on the molecule, and the blue squares (red circles) show H (O) contributions to it.

Close modal

We calculated the adsorption energies of water Eads in the presence of an electric field. The references to calculate Eads at each field strength were taken as the isolated subsystems optimized in the presence of the respective electric field embedding. In addition, we report the values of the adsorption energy related to the zero-field reference in the SI. The adsorption energies are shown in Fig. 3(a) and summarized in Table I. The highest adsorption energy is observed when no electric field or a slight positive field of 0.07 V/Å is applied. The value obtained without an electric field is 0.48 eV, which is larger than the values of 0.22–0.31 eV reported in previous studies.20–22 The reason for this discrepancy is that no van der Waals (vdW) interactions were included in these previous studies, as we show in Table I. The inclusion of the vdW interactions in Ref. 12 results in adsorption energies very similar to the ones we obtain. When increasing the field toward positive values, Eads gradually decreases (i.e., the adsorbed system becomes less stable). When a negative field is applied, Eads also gradually decreases down to a field of −0.44 V/Å, after which it remains approximately constant. This change in the trend corresponds to the flip of the preferred adsorption orientation discussed above.

TABLE I.

Adsorption energy Eads, potential energy barrier Ea, ZPE-corrected barrier of dissociation, and the corresponding quasi-harmonic quantum TST dissociation rates kd of a water monomer on a Pd(111) surface, depending on the applied electric field. The electric field of 1 V/Å corresponds to the surface charge of 0.0364 electrons per Pd atom (8.734 C/cm2) in this work. The values are calculated with PBE + vdWsurf, unless specified otherwise.

Electric field (V/Å)Eads (eV)Ea (eV)ZPE-corrected Ea (eV)kd300 K (s−1)Note
+0.74 0.40 1.32 1.12 4.3 × 10−5  
+0.44 0.44 1.18 0.98 8.6 × 10−3  
+0.15 0.47 1.08 0.86 5.3 × 10−1  
+0.07 0.48 1.06 ⋯ ⋯  
No field 0.48 1.03 0.83 4.7 This work 
 0.31 ⋯ ⋯ ⋯ This work, no vdW 
 0.33 ⋯ ⋯ ⋯ Ref. 38, PW91, no vdW 
 ⋯ 1.12 ⋯ ⋯ Ref. 19, PBE, no vdW 
 0.22 1.09 0.87 ⋯ Ref. 20, PW91, no vdW 
 0.30 1.05 ⋯ ⋯ Ref. 21, PW91, no vdW 
 0.31 1.15 0.96 ⋯ Ref. 22, PW91, no vdW 
 0.24 ⋯ ⋯ ⋯ Ref. 12, PBE, no vdW 
 0.42 ⋯ ⋯ ⋯ Ref. 12, revPBE + vdW-DFa 
 0.46 ⋯ ⋯ ⋯ Ref. 12, optPBE + vdW-DFa 
 0.53 ⋯ ⋯ ⋯ Ref. 12, optB88 + vdW-DFa 
−0.07 0.45 1.02 ⋯ ⋯  
−0.15 0.41 1.00 0.80 1.4 × 101  
−0.29 0.34 0.98 0.78 1.9 × 101  
−0.44 0.31 1.02 0.85 2.4 × 10−1  
−0.74 0.30 1.01 0.90 1.1 × 10−1  
−1.00 0.30 ⋯ ⋯ ⋯  
Electric field (V/Å)Eads (eV)Ea (eV)ZPE-corrected Ea (eV)kd300 K (s−1)Note
+0.74 0.40 1.32 1.12 4.3 × 10−5  
+0.44 0.44 1.18 0.98 8.6 × 10−3  
+0.15 0.47 1.08 0.86 5.3 × 10−1  
+0.07 0.48 1.06 ⋯ ⋯  
No field 0.48 1.03 0.83 4.7 This work 
 0.31 ⋯ ⋯ ⋯ This work, no vdW 
 0.33 ⋯ ⋯ ⋯ Ref. 38, PW91, no vdW 
 ⋯ 1.12 ⋯ ⋯ Ref. 19, PBE, no vdW 
 0.22 1.09 0.87 ⋯ Ref. 20, PW91, no vdW 
 0.30 1.05 ⋯ ⋯ Ref. 21, PW91, no vdW 
 0.31 1.15 0.96 ⋯ Ref. 22, PW91, no vdW 
 0.24 ⋯ ⋯ ⋯ Ref. 12, PBE, no vdW 
 0.42 ⋯ ⋯ ⋯ Ref. 12, revPBE + vdW-DFa 
 0.46 ⋯ ⋯ ⋯ Ref. 12, optPBE + vdW-DFa 
 0.53 ⋯ ⋯ ⋯ Ref. 12, optB88 + vdW-DFa 
−0.07 0.45 1.02 ⋯ ⋯  
−0.15 0.41 1.00 0.80 1.4 × 101  
−0.29 0.34 0.98 0.78 1.9 × 101  
−0.44 0.31 1.02 0.85 2.4 × 10−1  
−0.74 0.30 1.01 0.90 1.1 × 10−1  
−1.00 0.30 ⋯ ⋯ ⋯  
a

See Ref. 39 for details about the vdW-DF with different functionals.

We investigated how the character of the bond between the molecule and the surface changes with an electric field by projecting the total electronic DOS onto the atoms. We show such a projected DOS in Fig. 4. By comparing the positions of the molecular orbitals in the gas-phase to those on the surface, we observe in Fig. 4 that what used to be the HOMO level of the water molecule in isolation becomes substantially broader for field values ranging from +0.74 to −0.15 V/Å, which points toward a hybridization with the Pd d-bands. For more negative fields, the hybridization is drastically reduced, which is consistent with the “standing” orientation of the water molecule and the character of the HOMO level of water (centered on oxygen). This is also reflected on the Hirshfeld and Mulliken analyses of the charge accumulated on the molecule, shown in Fig. 3(b). On a positively charged surface, the Hirshfeld and Mulliken analyses agree for the total charge, even if they predict different values for individual atoms. The molecule donates roughly 0.2 e to the surface, indicating orbital hybridization and electron transfer from the molecule to the surface. At an electric field below −0.44 V/Å, when the molecule “stands,” the Mulliken analysis shows, in contrast, a negligibly small charge, while the Hirshfeld method shows about −0.15 e. We interpret the discrepancy between the Hirshfeld and Mulliken charges as a specific feature of the Hirshfeld definition of charge, which assigns local changes of the electron density to atoms proportionally to their free-atom electron density at the considered point. Therefore, a change of the “tail” of surface density can be assigned to the molecule even if there is no hybridization of the orbitals of the surface and the molecule.

FIG. 4.

Electronic density of states projected on the atomic species at different values of an external electric field. The black lines denote Pd, and the red lines represent the sum of O and H contributions. The green lines on the zero-field plot show the DOS of an isolated water molecule, aligned so that the 1s orbital of its oxygen matches the 1s orbital of the adsorbed water molecule, after which all values are shifted such that zero corresponds to the Fermi level of the calculations that include Pd(111).

FIG. 4.

Electronic density of states projected on the atomic species at different values of an external electric field. The black lines denote Pd, and the red lines represent the sum of O and H contributions. The green lines on the zero-field plot show the DOS of an isolated water molecule, aligned so that the 1s orbital of its oxygen matches the 1s orbital of the adsorbed water molecule, after which all values are shifted such that zero corresponds to the Fermi level of the calculations that include Pd(111).

Close modal

The nudged elastic band (NEB) paths are shown in Fig. 5, and the corresponding energy barriers are shown in Fig. 6. We verified the transition states by a normal mode analysis. In all cases, the transition state has exactly one imaginary-frequency mode, which corresponds to a direction of barrier crossing. The geometry of the transition state has little dependence on the surface charge, especially if compared to the reactant state. The activation energy Ea for each bias is shown in Fig. 7 and summarized in Table I.

FIG. 5.

Individual NEB node geometries of the reaction path found by the CI-NEB algorithm (side and top views). The large green spheres denote the 1st layer of the Pd(111) surface, and the white and blue spheres of the 2nd and 3rd Pd layers are visible in the hcp- and fcc-hollow sites of the surface, respectively. The label TS marks the transition state.

FIG. 5.

Individual NEB node geometries of the reaction path found by the CI-NEB algorithm (side and top views). The large green spheres denote the 1st layer of the Pd(111) surface, and the white and blue spheres of the 2nd and 3rd Pd layers are visible in the hcp- and fcc-hollow sites of the surface, respectively. The label TS marks the transition state.

Close modal
FIG. 6.

Potential energy along the MEP found by the CI-NEB algorithm. The points represent the energy of the NEB nodes relative to the reactant, and the x-axis shows the corresponding distance between the oxygen atom and the dissociated hydrogen. The barriers are shown for the electric field of +0.44 V/Å (red diamonds), no field (black circles), −0.29 V/Å (blue triangles), and −0.44 V/Å (green squares).

FIG. 6.

Potential energy along the MEP found by the CI-NEB algorithm. The points represent the energy of the NEB nodes relative to the reactant, and the x-axis shows the corresponding distance between the oxygen atom and the dissociated hydrogen. The barriers are shown for the electric field of +0.44 V/Å (red diamonds), no field (black circles), −0.29 V/Å (blue triangles), and −0.44 V/Å (green squares).

Close modal
FIG. 7.

Energy barrier of water splitting calculated by the CI-NEB algorithm. The black circles show the potential energy barrier, and the red squares show the ZPE-corrected barrier, calculated in the harmonic approximation. The dashed green line marks the border between two orientations of the reactant state, which are shown schematically in the insets.

FIG. 7.

Energy barrier of water splitting calculated by the CI-NEB algorithm. The black circles show the potential energy barrier, and the red squares show the ZPE-corrected barrier, calculated in the harmonic approximation. The dashed green line marks the border between two orientations of the reactant state, which are shown schematically in the insets.

Close modal

For zero bias, our Ea value lies below (by 20–120 meV) previously reported calculations,19–22 which also did not include dispersion interactions. Many previous studies employed the PW91 exchange-correlation (XC) potential, but it was shown on various metals that the water dissociation barrier calculated by PW91 and PBE is almost identical; e.g., see the SI of Ref. 22. Thus, the discrepancy with the literature is consistent with the observation of Litman et al. that the inclusion of vdW interactions decreases the barrier of dissociation of water on metals.24 The lowest barrier is observed at −0.29 V/Å, directly before the water molecule changes its adsorption orientation. The positive electric fields only increase Ea. It is reasonable to expect that the availability of negative charges on the surface stabilizes the adsorbed OH radical and favors the reaction, but there is more to it. Analyzing the MEP, we speculate that the positive field strengths stabilize the O–Pd bond at the reactant, as evidenced by the charge-transfer behavior. This strengthening makes it harder for water to move toward the point where the TS is. The negative field strengths destabilize the bond, matching, in fact, the strength of the O–Pd bond at around the flipping point. Beyond the flipping point, the molecule is physisorbed but has to flip back in order to reach the transition state (see the MEPs in Fig. 5), which again costs energy.

As discussed previously in the literature, the main nuclear quantum effect that affects the water dissociation on metallic surfaces around room temperature is zero-point energy, with tunneling playing a minor role.24 We confirm this observation here by an analysis of the frequency of the TS mode at all electric field strengths, reported in Table II. We estimate tunneling crossover temperatures Tc to be below 185 K for all field strengths, obtaining a very similar value as the one reported in Ref. 20 for this reaction without the application of a field. We remark that the highest Tc is observed for −0.74 V/Å and that the field strengths reported here can change Tc by 80 K. This observation shows that the application of the field and the consequent negative charging of the surface can considerably increase the curvature of the barrier (assuming that its shape around the TS is reasonably approximated by an inverted parabola) and increase the importance of tunneling at higher temperatures in this reaction.

TABLE II.

Imaginary frequency of the unstable mode at the transition state (ωTS) and the tunneling crossover temperature (Tc) for water dissociation on the Pd(111) surface at different electric fields.

Electric field (V/Å)ωTS (i*cm−1)Tc (K)
+0.74 451.7 103.4 
+0.44 564.4 129.2 
+0.15 621.7 142.4 
645.4 147.8 
−0.15 673.1 154.1 
−0.29 640.5 146.7 
−0.44 729.5 167.0 
−0.74 800.5 183.3 
Electric field (V/Å)ωTS (i*cm−1)Tc (K)
+0.74 451.7 103.4 
+0.44 564.4 129.2 
+0.15 621.7 142.4 
645.4 147.8 
−0.15 673.1 154.1 
−0.29 640.5 146.7 
−0.44 729.5 167.0 
−0.74 800.5 183.3 

We calculate the harmonic ZPE contributions to the reactant and transition states and summarize the changes in the dissociation barrier in Fig. 7 and in Table I. We note that a low-frequency vibrational mode of the reactant, corresponding to the rotation of water hydrogens around the oxygen, parallel to the surface plane, becomes progressively lower in frequency and even slightly negative (40 cm−1) toward higher positive field values. We have tightened the optimization thresholds and increased the numerical and basis set accuracies to confirm this observation. First, we can conclude that this model system is, indeed, slightly unstable above +0.44 V/Å  within the harmonic approximation and the molecule could be freely rotating parallel to the surface plane. Second, these modes do not play a significant role in the water dissociation path and, because this is a feature of the model that would disappear if water was embedded in an H-bonded network, for the ZPE analysis, we consistently ignored this mode in all calculations. We should add that because the mode is so low in energy (always below 30 cm−1), its contribution would always be less than 6 meV to the total ZPE.

The ZPE contribution decreases the barrier by roughly 0.19–0.20 eV at all points except −0.44 V/Å, where the decrease is 0.17 eV, as shown in Fig. 7. This is a considerable contribution in this case, representing 20% of the barrier height. However, such a constant effect is somewhat surprising, given that the geometry of the reactant state changes considerably. In order to resolve the effects of ZPE, we analyzed the vibrational density of states (vDOS) spectra of the reactant and transition states, which are shown in Fig. 8 (top and bottom, respectively). The most pronounced difference appears in the two hindered rotation modes of the reactant, which shift sharply from the range of 400–500 cm−1 down to 200–320 cm−1 when the field reaches −0.44 V/Å (the flipping point). These modes are shown in Fig. 9. It is important to note that hindered rotation modes are expected to be strongly anharmonic and thermally populated, as quantified in Refs. 51; therefore, the harmonic analysis presented here gives only an estimate of the effect.

FIG. 8.

Vibrational density of states of the reactant state (top) and of the transition state (bottom) of the water splitting reaction on a Pd(111) surface at different electric field values. The lowest-frequency mode of the reactant state becomes unstable at +0.44 and +0.74 V/Å; therefore, it is not shown.

FIG. 8.

Vibrational density of states of the reactant state (top) and of the transition state (bottom) of the water splitting reaction on a Pd(111) surface at different electric field values. The lowest-frequency mode of the reactant state becomes unstable at +0.44 and +0.74 V/Å; therefore, it is not shown.

Close modal
FIG. 9.

Hindered rotation modes of a water molecule on a Pd(111) surface at −0.15 V/Å (top) and −0.44 V/Å (bottom).

FIG. 9.

Hindered rotation modes of a water molecule on a Pd(111) surface at −0.15 V/Å (top) and −0.44 V/Å (bottom).

Close modal

Despite multiple shifts in the low-frequency part of the spectrum, the ZPE is dominated by the high-frequency modes, namely the H–O–H bending mode and one of the two O–H stretching modes, which exist in the reactant state and disappear in the transition state. This is the main source of the difference between the reactant and transition state contributions. We present a detailed account of the contribution of particular modes to the total ZPE in the SI. Even though the high-frequency modes also present a strong variation in the reactant state under different applied electric fields, the ZPE contributions are compensated for between the product and reactant states, as an analysis of the cumulative ZPE contribution shows (see SI), resulting in the essentially constant ZPE correction to the reaction barrier.

In order to estimate the extent to which reaction rates can be changed due to the electric field effect, we report reaction rates as calculated from the quantum quasi-harmonic approximation [Eq. (13) in Ref. 24]. We report these rates at a temperature of 300 K in Table I. This is an artificial temperature because the single molecule would have a considerable desorption probability at this temperature and would surely be very mobile. However, a single adsorbed water molecule is already a model system in this context. Therefore, these numbers serve only to build a more quantitative intuition. As expected, the suppression of the reaction rate at positive field values (up to five orders of magnitude) is much more significant than its increase at negative field values (only up to a factor of 4). We predict that the flipping point of the first water layer will also limit the increase in this reaction rate when water is in the liquid state. Because in the liquid, for the first water layer to flip, hydrogen bonds would need to be broken, the flipping point should happen at considerably larger effective field strengths, even if it is difficult to give a quantitative estimate due to the entropic contributions in the liquid. We know from our calculations that the charge on the surface changes by 0.036 e/(V/Å) per surface Pd atom, which puts the flipping point (0.4 V/Å) at an excess surface charge of around −0.014 e/Pd. This excess surface charge can be compared to other studies in the liquid, for example, Ref. 52, in which larger surface excess charges (0.08 e/Pt) were reported to cause a structuring of the first water layer on a Pt surface.

We presented an ab initio study of water dissociation on the Pd(111) surface using a slab model embedded in an electric field in a 3D periodic setup. The electric field effectively induces opposite charges on each surface of the thick metallic slab and mimics the effect of an applied potential bias at the metallic electrode surface.

Analyzing the reaction paths at different field strengths, we concluded that the main feature that affects the reaction barrier is the geometry of the reactant state, which is strongly influenced by the applied field. Between field strengths of −0.3 and −0.4 V/Å, the reactant state abruptly changes its geometry, with the water molecule pointing the hydrogen atoms to the surface. We demonstrated that the lowest barrier (and consequently the highest reaction rate) for the monomer dissociation lies at the electric field strength at which the molecule has these two adsorption geometries equally preferred, i.e., at the flipping point. We rationalized that at this point, the water–surface bond is optimally weakened, making it easier for the system to reach its transition state.

Furthermore, we showed that the application of negative electric fields, which, in the convention used in this paper, causes negative charges to be accumulated at the surface, increases the tunneling crossover temperature significantly, albeit remaining below 200 K for this reaction and the field strengths studied in this work. The NQEs were explicitly evaluated limited to the ZPE in the harmonic approximation, demonstrating that the reaction barrier is reduced by roughly 0.2 eV at all field strengths. The reason for such a weak dependence on the field strength is the dominant role of bond stretching modes in the ZPE of the reaction and the mutual cancellation of contributions of different vibrational modes between the reactant and transition states. The observation that the ZPE contributions are almost constant throughout a large variety of field strengths could simplify more involved simulation protocols involving MD simulations, where one could estimate the NQE corrections at a single field strength and propagate that to others—even if with care—because the effect of NQEs may not be uniform across field strengths for different liquid water properties.53 

Although we studied a highly idealized system, we expect these findings to be relevant for more complex scenarios. Chemistry under applied electric fields is a booming field in its own right, and extending some of the protocols presented in this paper for methods that can capture nuclear tunneling could open new exciting avenues for situations where such effects play a role.54 In electrochemistry, the change in the orientation of the first water layer under different potential biases has been observed experimentally,5 but the role of these effects in the chemistry that happens at these interfaces is still not fully understood. We believe that further studies in liquid water, taking into account the insights obtained here, namely (i) that the reaction barrier is the lowest at the flipping point of the reactant, (ii) that nuclear tunneling contributions are more significant on negatively charged surfaces, and (iii) that the ZPE has a large impact on the dissociation rate, should lead to a better understanding of the factors determining the water splitting reaction on metallic substrates. Finally, we note that we constrained the total charge of the simulation cell to be zero, while a real electrode would allow charge fluctuations on a larger scale and, possibly, with a lower energy penalty—a possibility that should be explored in future work.

Additional evaluations of adsorption energies, mode resolved contributions to the zero-point energy, and calculations of crossover temperatures are available in the supplementary material.

We acknowledge the fruitful discussions with Alexandre Reily Rocha and Luana Sucupira Pedroza, who have immensely helped us navigate this subject. We further thank Yair Litman and Luana Sucupira Pedroza for reading a preliminary version of this manuscript and giving useful suggestions. We would like to thank the Max Planck Computing and Data Center (MPCDF) for computing time. K.F. acknowledges the support from the International Max Planck Research School for Ultrafast Imaging and Structural Dynamics (IMPRS-UFAST).

The authors have no conflicts to disclose.

Karen Fidanyan: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Software (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Guoyuan Liu: Software (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Mariana Rossi: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Methodology (equal); Project administration (lead); Supervision (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

Phonon calculations for the initial and transition states of the water splitting reaction and the geometry relaxations of the adsorbed reactants and products for different values of an electric field are uploaded to the NOMAD repository. The dataset DOI is 10.17172/NOMAD/2022.09.15–1.

1.
T.
Würger
,
C.
Feiler
,
G. B.
Vonbun-Feldbauer
,
M. L.
Zheludkevich
, and
R. H.
Meißner
, “
A first-principles analysis of the charge transfer in magnesium corrosion
,”
Sci. Rep.
10
,
15006
(
2020
).
2.
H. N.
Nong
,
L. J.
Falling
,
A.
Bergmann
,
M.
Klingenhof
,
H. P.
Tran
,
C.
Spöri
,
R.
Mom
,
J.
Timoshenko
,
G.
Zichittella
,
A.
Knop-Gericke
,
S.
Piccinin
,
J.
Pérez-Ramírez
,
B. R.
Cuenya
,
R.
Schlögl
,
P.
Strasser
,
D.
Teschner
, and
T. E.
Jones
, “
Key role of chemistry versus bias in electrocatalytic oxygen evolution
,”
Nature
587
,
408
413
(
2020
).
3.
Y.-H.
Wang
,
S.
Zheng
,
W.-M.
Yang
,
R.-Y.
Zhou
,
Q.-F.
He
,
P.
Radjenovic
,
J.-C.
Dong
,
S.
Li
,
J.
Zheng
,
Z.-L.
Yang
,
G.
Attard
,
F.
Pan
,
Z.-Q.
Tian
, and
J.-F.
Li
, “
In situ Raman spectroscopy reveals the structure and dissociation of interfacial water
,”
Nature
600
,
81
85
(
2021
).
4.
Y.
Tong
,
F.
Lapointe
,
M.
Thämer
,
M.
Wolf
, and
R. K.
Campen
, “
Hydrophobic water probed experimentally at the gold electrode/aqueous interface
,”
Angew. Chem., Int. Ed.
56
,
4211
4214
(
2017
).
5.
M. F.
Toney
,
J. N.
Howard
,
J.
Richer
,
G. L.
Borges
,
J. G.
Gordon
,
O. R.
Melroy
,
D. G.
Wiesler
,
D.
Yee
, and
L. B.
Sorensen
, “
Voltage-dependent ordering of water molecules at an electrode–electrolyte interface
,”
Nature
368
,
444
446
(
1994
).
6.
K.
Laasonen
,
M.
Sprik
,
M.
Parrinello
, and
R.
Car
, “
‘Ab initio’ liquid water
,”
J. Chem. Phys.
99
,
9080
9089
(
1993
).
7.
D.
Marx
,
M. E.
Tuckerman
,
J.
Hutter
, and
M.
Parrinello
, “
The nature of the hydrated excess proton in water
,”
Nature
397
,
601
604
(
1999
).
8.
D.
Marx
,
M. E.
Tuckerman
, and
M.
Parrinello
, “
Solvated excess protons in water: Quantum effects on the hydration structure
,”
J. Phys.: Condens. Matter
12
,
A153
A159
(
2000
).
9.
X. Z.
Li
,
M. I.
Probert
,
A.
Alavi
, and
A.
Michaelides
, “
Quantum nature of the proton in water-hydroxyl overlayers on metal surfaces
,”
Phys. Rev. Lett.
104
,
066102
(
2010
).
10.
M.
Rossi
,
H.
Liu
,
F.
Paesani
,
J.
Bowman
, and
M.
Ceriotti
, “
Communication: On the consistency of approximate quantum dynamics simulation methods for vibrational spectra in the condensed phase
,”
J. Chem. Phys.
141
,
181101
(
2014
).
11.
M.
Ceriotti
,
W.
Fang
,
P. G.
Kusalik
,
R. H.
McKenzie
,
A.
Michaelides
,
M. A.
Morales
, and
T. E.
Markland
, “
Nuclear quantum effects in water and aqueous systems: Experiment, theory, and current challenges
,”
Chem. Rev.
116
,
7529
7550
(
2016
).
12.
J.
Carrasco
,
J.
Klimeš
, and
A.
Michaelides
, “
The role of van der Waals forces in water adsorption on metals
,”
J. Chem. Phys.
138
,
024708
(
2013
).
13.
S.
Izvekov
and
G. A.
Voth
, “
Ab initio molecular dynamics simulation of the Ag(111)-water interface
,”
J. Chem. Phys.
115
,
7196
7206
(
2001
).
14.
A.
Groß
and
S.
Sakong
, “
Ab initio simulations of water/metal interfaces
,”
Chem. Rev.
122
,
10746
10776
(
2022
).
15.
L. S.
Pedroza
,
A.
Poissier
, and
M.-V.
Fernández-Serra
, “
Local order of liquid water at metallic electrode surfaces
,”
J. Chem. Phys.
142
,
034706
(
2015
).
16.
Z. K.
Goldsmith
,
M. F. C.
Andrade
, and
A.
Selloni
, “
Effects of applied voltage on water at a gold electrode interface from ab initio molecular dynamics
,”
Chem. Sci.
12
,
5865
5873
(
2021
).
17.
O.
Sugino
,
I.
Hamada
,
M.
Otani
,
Y.
Morikawa
,
T.
Ikeshoji
, and
Y.
Okamoto
, “
First-principles molecular dynamics simulation of biased electrode/solution interface
,”
Surf. Sci.
601
,
5237
5240
(
2007
).
18.
N.
Dubouis
and
A.
Grimaud
, “
The hydrogen evolution reaction: From material to interfacial descriptors
,”
Chem. Sci.
10
,
9165
9181
(
2019
).
19.
G.
Wang
,
S.
Tao
, and
X.
Bu
, “
A systematic theoretical study of water dissociation on clean and oxygen-preadsorbed transition metals
,”
J. Catal.
244
,
10
16
(
2006
).
20.
Y.
Cao
and
Z.-X.
Chen
, “
Theoretical studies on the adsorption and decomposition of H2O on Pd(111) surface
,”
Surf. Sci.
600
,
4572
4583
(
2006
).
21.
A. A.
Phatak
,
W. N.
Delgass
,
F. H.
Ribeiro
, and
W. F.
Schneider
, “
Density functional theory comparison of water dissociation steps on Cu, Au, Ni, Pd, and Pt
,”
J. Phys. Chem. C
113
,
7269
7276
(
2009
).
22.
J. L. C.
Fajín
,
M. N. D. S.
Cordeiro
,
F.
Illas
, and
J. R. B.
Gomes
, “
Descriptors controlling the catalytic activity of metallic surfaces toward water splitting
,”
J. Catal.
276
,
92
100
(
2010
).
23.
D.
Donadio
,
L. M.
Ghiringhelli
, and
L.
Delle Site
, “
Autocatalytic and cooperatively stabilized dissociation of water on a stepped platinum surface
,”
J. Am. Chem. Soc.
134
,
19217
19222
(
2012
).
24.
Y.
Litman
,
D.
Donadio
,
M.
Ceriotti
, and
M.
Rossi
, “
Decisive role of nuclear quantum effects on surface mediated water dissociation at finite temperature
,”
J. Chem. Phys.
148
,
102320
(
2018
).
25.
Y.
Qian
,
I.
Hamada
,
M.
Otani
, and
T.
Ikeshoji
, “
Inhibition of water dissociation on a pitted Pt(111) surface: First principles study
,”
Catal. Today
202
,
163
167
(
2013
), part of Special Issue: Electrocatalysis.
26.
M.
Otani
and
O.
Sugino
, “
First-principles calculations of charged surfaces and interfaces: A plane-wave nonrepeated slab approach
,”
Phys. Rev. B
73
,
115407
(
2006
).
27.
N.
Bonnet
,
T.
Morishita
,
O.
Sugino
, and
M.
Otani
, “
First-principles molecular dynamics at a constant electrode potential
,”
Phys. Rev. Lett.
109
,
266101
(
2012
).
28.
R.
Sundararaman
,
W. A.
Goddard
, and
T. A.
Arias
, “
Grand canonical electronic density-functional theory: Algorithms and applications to electrochemistry
,”
J. Chem. Phys.
146
,
114104
(
2017
).
29.
S.
Surendralal
,
M.
Todorova
,
M. W.
Finnis
, and
J.
Neugebauer
, “
First-principles approach to model electrochemical reactions: Understanding the fundamental mechanisms behind Mg corrosion
,”
Phys. Rev. Lett.
120
,
246801
(
2018
).
30.
S.
Hagiwara
,
C.
Hu
,
S.
Nishihara
, and
M.
Otani
, “
Bias-dependent diffusion of a H2O molecule on metal surfaces by the first-principles method under the grand-canonical ensemble
,”
Phys. Rev. Mater.
5
,
065001
(
2021
).
31.
F.
Deißenbeck
,
C.
Freysoldt
,
M.
Todorova
,
J.
Neugebauer
, and
S.
Wippermann
, “
Dielectric properties of nanoconfined water: A canonical thermopotentiostat approach
,”
Phys. Rev. Lett.
126
,
136803
(
2021
).
32.
R.
Khatib
,
A.
Kumar
,
S.
Sanvito
,
M.
Sulpizi
, and
C. S.
Cucinotta
, “
The nanoscale structure of the Pt-water double layer under bias revealed
,”
Electrochim. Acta
391
,
138875
(
2021
).
33.
L. J. V.
Ahrens-Iwers
,
M.
Janssen
,
S. R.
Tee
, and
R. H.
Meißner
, “
ELECTRODE: An electrochemistry package for atomistic simulations
,”
J. Chem. Phys.
157
,
084801
(
2022
).
34.
L. S.
Pedroza
,
P.
Brandimarte
,
A. R.
Rocha
, and
M.-V.
Fernández-Serra
, “
Bias-dependent local structure of water molecules at a metallic interface
,”
Chem. Sci.
9
,
62
69
(
2018
).
35.
M.
Otani
,
I.
Hamada
,
O.
Sugino
,
Y.
Morikawa
,
Y.
Okamoto
, and
T.
Ikeshoji
, “
Structure of the water/platinum interface—a first principles simulation under bias potential
,”
Phys. Chem. Chem. Phys.
10
,
3609
3612
(
2008
).
36.
T.
Ikeshoji
,
M.
Otani
,
I.
Hamada
, and
Y.
Okamoto
, “
Reversible redox reaction and water configuration on a positively charged platinum surface: First principles molecular dynamics simulation
,”
Phys. Chem. Chem. Phys.
13
,
20223
20227
(
2011
).
37.
T.
Dufils
,
G.
Jeanmairet
,
B.
Rotenberg
,
M.
Sprik
, and
M.
Salanne
, “
Simulating electrochemical systems by combining the finite field method with a constant potential electrode
,”
Phys. Rev. Lett.
123
,
195501
(
2019
).
38.
A.
Michaelides
,
V. A.
Ranea
,
P. L.
de Andres
, and
D. A.
King
, “
General model for water monomer adsorption on close-packed transition and noble metal surfaces
,”
Phys. Rev. Lett.
90
,
216102
(
2003
).
39.
J.
Klimeš
,
D. R.
Bowler
, and
A.
Michaelides
, “
Chemical accuracy for the van der Waals density functional
,”
J. Phys.: Condens. Matter
22
,
022201
(
2009
).
40.
V.
Blum
,
R.
Gehrke
,
F.
Hanke
,
P.
Havu
,
V.
Havu
,
X.
Ren
,
K.
Reuter
, and
M.
Scheffler
, “
Ab initio molecular simulations with numeric atom-centered orbitals
,”
Comput. Phys. Commun.
180
,
2175
2196
(
2009
).
41.
J.
Neugebauer
and
M.
Scheffler
, “
Adsorbate-substrate and adsorbate-adsorbate interactions of Na and K adlayers on Al(111)
,”
Phys. Rev. B
46
,
16067
16080
(
1992
).
42.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
3868
(
1996
).
43.
V. G.
Ruiz
,
W.
Liu
,
E.
Zojer
,
M.
Scheffler
, and
A.
Tkatchenko
, “
Density-functional theory with screened van der Waals interactions for the modeling of hybrid inorganic-organic systems
,”
Phys. Rev. Lett.
108
,
146103
(
2012
).
44.
V. G.
Ruiz
,
W.
Liu
, and
A.
Tkatchenko
, “
Density-functional theory with screened van der Waals interactions applied to atomic and molecular adsorbates on close-packed and non-close-packed surfaces
,”
Phys. Rev. B
93
,
035118
(
2016
).
45.
K.
Fidanyan
, “Modified phonopy-FHI-aims interface,” Github. https://github.com/fidanyan/phonopy
46.
A.
Togo
and
I.
Tanaka
, “
First principles phonon calculations in materials science
,”
Scr. Mater.
108
,
1
5
(
2015
).
47.
G.
Henkelman
,
B. P.
Uberuaga
, and
H.
Jónsson
, “
A climbing image nudged elastic band method for finding saddle points and minimum energy paths
,”
J. Chem. Phys.
113
,
9901
9904
(
2000
).
48.
V.
Kapil
,
M.
Rossi
,
O.
Marsalek
,
R.
Petraglia
,
Y.
Litman
,
T.
Spura
,
B.
Cheng
,
A.
Cuzzocrea
,
R. H.
Meißner
,
D. M.
Wilkins
,
B. A.
Helfrecht
,
P.
Juda
,
S. P.
Bienvenue
,
W.
Fang
,
J.
Kessler
,
I.
Poltavsky
,
S.
Vandenbrande
,
J.
Wieme
,
C.
Corminboeuf
,
T. D.
Kühne
,
D. E.
Manolopoulos
,
T. E.
Markland
,
J. O.
Richardson
,
A.
Tkatchenko
,
G. A.
Tribello
,
V.
Van Speybroeck
, and
M.
Ceriotti
, “
i-PI 2.0: A universal force engine for advanced molecular simulations
,”
Comput. Phys. Commun.
236
,
214
223
(
2019
).
49.
K.
Fidanyan
, “
Theoretical investigations of nuclear quantum effects in weakly bonded metal-molecular interfaces
,” Ph.D. thesis,
Humboldt-Universität zu Berlin
(
2022
).
50.
E.
Bitzek
,
P.
Koskinen
,
F.
Gähler
,
M.
Moseler
, and
P.
Gumbsch
, “
Structural relaxation made simple
,”
Phys. Rev. Lett.
97
,
170201
(
2006
).
51.
M.
Rossi
, “
Progress and challenges in ab initio simulations of quantum nuclei in weakly bonded systems
,”
J. Chem. Phys.
154
,
170902
(
2021
).
52.
M.
Otani
,
I.
Hamada
,
O.
Sugino
,
Y.
Morikawa
,
Y.
Okamoto
, and
T.
Ikeshoji
, “
Electrode dynamics from first principles
,”
J. Phys. Soc. Jpn.
77
,
024802
(
2008
).
53.
G.
Cassone
,
J.
Sponer
, and
F.
Saija
, “
Ab initio molecular dynamics studies of the electric-field-induced catalytic effects on liquids
,”
Top. Catal.
65
,
40
58
(
2022
).
54.
O.
Kirshenboim
,
A.
Frenklah
, and
S.
Kozuch
, “
Switch chemistry at cryogenic conditions: Quantum tunnelling under electric fields
,”
Chem. Sci.
12
,
3179
3187
(
2020
).

Supplementary Material