Solvent Induced Proton Hopping at a Water-Oxide Interface

Despite widespread interest, a detailed understanding of the dynamics of proton transfer at interfaces is lacking. Here we use ab initio molecular dynamics to unravel the connection between interfacial water structure and proton transfer for the widely studied and experimentally well-characterized water-ZnO$(10\bar{1}0)$ interface. We find that upon going from a single layer of adsorbed water to a liquid multilayer changes in the structure are accompanied by a dramatic increase in the proton transfer rate at the surface. We show how hydrogen bonding and rather specific hydrogen bond fluctuations at the interface are responsible for the change in the structure and proton transfer dynamics. The implications of this for the chemical reactivity and for the modelling of complex wet oxide interfaces in general are also discussed.

the adsorption energy by at most +15 meV/H 2 O. VASP is a plane wave code that uses the projector augmented wave method to treat the core electrons. We used six valence electrons for oxygen (2s 2 2p 4 ) and twelve for zinc (3d 10 4s 2 ) and a cut-off of 500 eV for the plane wave expansion. We have employed a k-point mesh density of 6 × 4 × 1 per surface unit cell with the Monkhorst-Pack scheme. S3 Because it is only possible to sample k-space at the Γ point in CP2K a 6 × 4 unit cell has been used to reproduce the same k-point density as in the VASP calculations. A dipole correction has been used in VASP while in CP2K the periodicity along the direction of the vacuum has been explicitly removed.
The clean ZnO(1010) is modelled using a 4 bilayer slab with the bottom bilayer kept fixed at its bulk position, using a range of surface unit cells and a vacuum region of at least 15Å. The structure of a water monomer adsorbed on the surface is shown in Fig. S1(a). It is modelled as an H 2 O molecule adsorbed on the surface at 1/6 ML (i.e. using a (3 × 2) surface unit cell). It has been shown in Ref. S4 that at 1/6 ML the interaction between neighboring water molecules is negligible. Also going from a 3 to a 4 bilayer slab gave rise to a change in the adsorption energy of the monomer of just 11 meV. We compared the relative stability of water adsorbed at ML coverage in the molecular (M) and partially dissociated (PD) states, as shown in Figs S1(b) and (c), respectively.

Monomer and ML water adsorption
The adsorption energies of water at low coverage (1/6 ML) and at ML coverage obtained with the two codes using several functionals are shown in Table S1. The table shows that water becomes more stable upon increasing the coverage from 1/6 to 1 ML due to hydrogen bonding and the PD state is more stable than the M state (as noticed in Ref. S4,S5 ). This holds regardless of the code used (whether it is CP2K or VASP), or if dispersion is accounted for (using PBE-D2 or optPBE-vdW and optB86b-vdW S6,S7 ) and also if a fraction of exact exchange is used with the HSE06 functional. S8 The adsorption energy computed with CP2K is at most 10 meV/H 2 O more stable than the one obtained from VASP, indicating that the basis set superposition error in our CP2K calculations is small. Adding a dispersion correction with any of the schemes used or the inclusion of 1/4 exact exchange with the HSE06 functional yields to further stabilization compared to PBE. On the whole, the energy difference between the PD and the M state is ≈ 0.1 eV with all functionals, with PBE being closer to the experimental value in Ref., S5 even if zero-point energy is accounted for. Table S1: Adsorption energy of water on ZnO(1010) predicted according to different functionals, in the case of monomer adsorption (at a coverage θ = 1/6) and at monolayer coverage in the case of molecular (θ M = 1) and partially dissociated (θ P D = 1) states. Within brackets are adsorption energies corrected for zero point energy effects, obtained in the harmonic approximation from frequency calculations. Unless specified all reported values are obtained with VASP. All quantities are in eV/H 2 O.

S3
AIMD simulations Set-up AIMD simulations of a LF and of the water ML on ZnO(1010) have been performed with the CP2K/QUICKSTEP package using the set-up described earlier. The slab of the AIMD simulations is made of a (6 × 3) unit cell with 18 water molecules in the case of the water ML and 144 in the LF simulation. The initial structure of the AIMD of the water ML was for water in the optimized M state, whereas the initial structure of the LF simulation has been obtained by running a 2 ns long classical MD simulation using a shell model for ZnO (1010), S9 a polarizable SPC water model S10 and by using a Lennard-Jones potential to model the water-ZnO interactions. The parameters for the Lennard-Jones potential have been obtained by a fitting to our DFT-PBE data for water monomer adsorption. The AIMD simulations are performed within the Born-Oppenheimer approximation in the canonical ensemble at a target temperature of 360 K. 360 K has been used to partially compensate for the lower diffusivity of water using the PBE functional. S11 The AIMD trajectories are 40 ps long with a 1 fs timestep, and we use deuterium masses for the hydrogens. This enables stable AIMD simulations to be performed with a 1 fs timestep. Specifically, the energy drift for the AIMD of liquid water on ZnO(1010) is < 0.01 meV/(atom·ps). The first 5 ps of the simulation are used for equilibration and the remaining 35 for analysis.
Equilibration in all the AIMD simulations has been performed for the initial 3 ps using the Bussi-Donadio-Parrinello thermostat, S12 while for the remaining time the target temperature has been maintained with the Nosé-Hoover chain thermostat. S13

PBE-D2
A typical structure of the LF on ZnO(1010) from PBE-D2 is shown in Fig. S2(a). The planar averaged density maximum of region 1 in both profiles is much greater than the average density of liquid water simply because of the pronounced layering close to the surface (see e.g. Ref. S14 ).
When comparing PBE and PBE-D2 we find that because of dispersion forces the density in the contact layer obtained from PBE-D2 has seen a slight increase in the first two peaks in region 1 at 2.0 and 2.5Å. This results in an increase in the coverage from 1.16 ± 0.03 to 1.29±0.04 of a ML. The density of PBE-D2 bulk liquid water (region 2) is 0.85±0.14 g/cm 3 , slightly larger than the PBE value of 0.75 ± 0.12 g/cm 3 . Further, the decay at the liquidvacuum interface is more rapid in the PBE-D2 density. Fig. S3 illustrates the percentage of adsorbed Hs as a function of the simulation time for the PBE-D2 LF simulation. the average distance between a surface O atom and a chemisorbed H plus 3 times the standard deviation. The level of dissociation does not change significantly compared to the PBE simulation ( Fig. 3(b) in the main text). Indeed, there is a statistically insignificant increase in the percentage of adsorbed Hs to 56 ± 5% using PBE-D2, compared to 55 ± 5% with PBE.
Above all, it can be seen that the rapid fluctuations in the number of adsorbed Hs is also seen in the PBE-D2 study.
Note that here and in Fig. 3(a) and (b) in the main text we consider as adsorbed Hs as those which are within 1.2Å of a substrate O atom. In the calculation of the number of hops Hs which are covalently bonded to an oxygen atom within the contact layer are defined as being no further than the mean O-H bond length plus one standard deviation (Fig. 3(c)).
We may define a discrete hop even in the LF since a free energy barrier 100 meV is present also in this case and protons are not shared, as opposed to the case of e.g. (see e.g. S15 ).

Dissociation barrier at ML coverage
In the main manuscript we reported results for the water dissociation barrier as a function of intermolecular distance. The energy profiles for each of these dissociation reactions are shown in Fig. S4. Specifically what we report are energy profiles as a function of the O w -H distance for specific values of the O w -O d distances (as defined in the main text). The particular procedure we applied for computing the dissociation barriers was to perform AIMD Figure S3: Time evolution of the percentage of H atoms adsorbed on the surface for the liquid water film simulation shown after 5 ps of equilibration. The inclusion of dispersion does not alter the proton transfer dynamics significantly nor does it produce a statistically relevant increase in the percentage of adsorbed Hs compared to the PBE simulation (see Fig. 3 Fig. S4). It can be seen that using the HSE06 functional only increases the barrier to dissociation by less than 2 meV/H 2 O. Figure S4: Energy pathways along the O w -H bond-length for the transition from molecular water adsorption to partially dissociated adsorption. ∆E is the total energy difference from the initial (

Comparison with previous work
Although we observed an increase in the coverage, the level of dissociation has not changed significantly with the presence of the LF. Accounting for dispersion forces with the empirical PBE-D2 S1 approach also does not alter our conclusion (see Fig. S3). In addition, our low coverage PBE results agree with previous work S4,S5 and the use of other dispersion corrections (with the optPBE-vdW and optB86b-vdW functionals S6,S7 ) or of hybrid functionals (HSE06 S8 ) does not change the stability between the M and the PD state at ML coverage (Table S1). Further, the barrier to dissociation at ML coverage with HSE06 is not increased significantly compared to PBE (see Fig. S4). However, Raymand et al. S16 reported an increase in the level of dissociation to 80% in their reactive force field (ReaxFF) molecular dynamics simulations. The discrepancy has likely arisen from the different methods used.
It is somewhat difficult to establish how reliable the ReaxFF method is also because details of the H and O parameters used have not been reported. S16,S17 We note that Holtaus et al. S18 have compared the level of dissociation of liquid water on ZnO(1210) predicted by DFT-PBE, tight binding and ReaxFF and found a similar discrepancy between DFT-PBE and ReaxFF. They attributed the discrepancy between the two to a change in the charges in the contact layer upon increasing the coverage that is not accounted for in the ReaxFF.