Long-range proton and hydroxide ion transfer dynamics at the water/CeO2 interface in the nanosecond regime: reactive molecular dynamics simulations and kinetic analysis

The structural properties, dynamical behaviors, and ion transport phenomena at the interface between water and cerium oxide are investigated by reactive molecular dynamics (MD) simulations employing neural network potentials (NNPs). The NNPs are trained to reproduce density functional theory (DFT) results, and DFT-based MD (DFT-MD) simulations with enhanced sampling techniques and refinement schemes are employed to efficiently and systematically acquire training data that include diverse hydrogen-bonding configurations caused by proton hopping events. The water interfaces with two low-index surfaces of (111) and (110) are explored with these NNPs, and the structure and long-range proton and hydroxide ion transfer dynamics are examined with unprecedented system sizes and long simulation times. Various types of proton hopping events at the interface are categorized and analyzed in detail. Furthermore, in order to decipher the proton and hydroxide ion transport phenomena along the surface, a counting analysis based on the semi-Markov process is formulated and applied to the MD trajectories to obtain reaction rates by considering the transport as stochastic jump processes. Through this model, the coupling between hopping events, vibrational motions, and hydrogen bond networks at the interface are quantitatively examined, and the high activity and ion transport phenomena at the water/CeO2 interface are unequivocally revealed in the nanosecond regime.


S.I.
Biased MD and Temperature Accelerated MD

S.I.1 Biased MD
The Biased MD (BMD) is introduced by adding a bias potential, which is a function of collective variables (CVs), to the potential energy of the system.These CVs are chosen so that the corresponding energy barrier can be easily overcome within the timescales of DFT-MD simulations.In this study, BMD was employed to promote proton hopping between water molecules and surface O atoms.For this purpose, the minimum distance between H atoms and surface O atoms was taken as CV, where it was defined as   ln exp Here, si is the distance between the reference atom and atom i in the system.In this study, an approximate free energy curve that was a function of the minimum distance between a particular O atom of the CeO2(111) surface and H atoms of water molecules was first obtained by the umbrella sampling prior to the BMD simulation, and this energy curve with the opposite sign was used to design a proper bias potential.The parameter β was set to 500 Å.
This bias potential flattens the barrier of proton hopping between water molecules and surface O atoms.
Since the approximate free energy curve does not need to be the precise true free energy curve of the system as long as the curve flatten the barrier, we roughly estimated the function by using a smaller simulation cell without dispersion correction.The free energy curve was obtained by the umbrella sampling method employing a CeO2(111) slab model that had a surface area of p(3×3) and three O-Ce-O tri-layers, and 50 water molecules were placed over the surface.The DFT-MD simulation of 10 ps was performed for each umbrella window, and a total of 10 windows were used, where a quadratic potential centered at predefined points was applied with a force constant k ranging from 2 to 60 eV/Å 2 depending on windows (see Table S1).During the umbrella sampling, the bottom O−Ce−O tri-layers were fixed, the mass of hydrogen was replaced with that of deuterium, and the timestep and temperature were set to 1.0 fs and 360 K, respectively.The biased probability distributions from different umbrella windows were reweighted and patched using the weighted histogram analysis method (WHAM).[1] The computed free energy curve is shown in Figure S1, along with the histograms of the minimum distance for each window.
Table S1.Center of restraints s (k) and force constant k.After estimating the free energy curve, a bias function that roughly cancels the energy barrier was designed to smooth out unnecessary detailed profiles of the curve (see Figure S1).This potential was applied to all surface O atoms in the BMD simulations.The same bias potential was used for the BMD simulations at the water/CeO2(110) interface.The PLUMED package [2] interfaced with CP2K was used for the enhanced sampling.

S.I.2 Temperature Accelerated MD
In the temperature accelerated MD (TAMD), which is also called driven adiabatic free energy dynamics (d-AFED), [3][4] the following Lagrangian is introduced, Here, L0 is the unbiased Lagrangian of the system, and a set of auxiliary variables {sα} is introduced that is coupled to the CVs {Sα} by a harmonic potential with force constant kα.The auxiliary subsystem {sα} is thermostatted to high-temperature T h , while the real system is thermostatted to a much lower system temperature T 0 , which was set at 360 K in this study.The mass of auxiliary variables {μ α } is set to much larger than the atomic mass of the system so that the adiabatic decoupling of the auxiliary variables from the degrees of freedom in the system is guaranteed.
In this work, we introduced the coordination number as CVs.The ordinary coordination number between two atom types is defined as where NX is a number of atoms with atom type X, and cij is a so-called switching function and represented as where d ij is the distance between atom i and j belonging to the atom types X and Y, respectively, and p and q are parameters to define the shape of the switching function.
0 XY R is a cutoff bond distance between atom types X and Y. Based on this coordination number, an associated coordination number is defined as where the indices i, j, and k run over all atoms of atom types X, Y, and Z, respectively.We used the associated coordination number defined as CN[Ce-O-H] to promote proton hopping of water molecules or hydroxide ions coordinated to surface Ce atoms.Here we introduced independent CVs for each of the surface Ce atoms, resulting S5 in a total of 16 CVs for TAMD simulations for each surface.The parameters of p = q = 6 were employed, and cutoff bond distances between Ce−O and O−H were set to 2.8 Å and 1.2 Å, respectively.The TAMD simulations were performed with μ α = 50.0amu, k α = 2.0 a.u., and the temperature for the auxiliary variables was set to T h = 2,000 K.

S.II. Construction and Refinement of NNP
The size of the embedding network (DNN for symmetric function) was set to (10,20,40,80).The fitting network had hidden layers of size (240, 240, 240).The activation function was tanh.The number of training batches was set to 10 6 to ensure sufficient convergence of the loss function, and the batch size was set to 1.The learning rate started from 5×10 −3 and finally decayed to 1.1×10 −8 .
In Figure S2, the improvements in NNP accuracy and oxygen density distribution deviation during the refinement procedure are shown.Here, the maximum sample standard deviation (SSD) of the forces is defined as

S.III. Adsorption Structure of a Water Molecule on the CeO2 Surface
The adsorption structures of a single water molecule were calculated by geometry optimization for two lowindex surfaces.In calculating the energies of gas-phase molecules, the cubic simulation cell of a = b = c = 15.0Å was used.For the calculation of the energies of slab models, the simulation cell was modified to have c = 30 Å from that described in the main text to ensure a vacuum space of at least 15 Å.For geometry optimization, the forces on all atoms were minimized to less than 0.02 eV/Å (4.5×10 −4 hartree/bohr).The adsorption energy of a water molecule E A to the surface was calculated according to where Ewater+surf is the total electronic energy of the water-surface system, while Esurf and Ewater are the energies of the pristine CeO2 slab model and an isolated water molecule, respectively.In this definition, the more negative value of adsorption energy indicates a stronger binding to the surface.The adsorption energies are given in units of eV in Table S2, and the adsorption structures are shown in Figure S3.

S.VI. Radial Distribution Functions
The radial distribution functions (RDFs) between the surface Ce atoms and water O atoms calculated from the NNP-MD simulations for p(8×8) surfaces are shown in Figure S6.The adsorption threshold for each system is also shown.Note that this threshold could affect the values of the proton/hydroxide transfer numbers between surface species.

S.VII. Decomposed Density Profiles
Figure S7 shows the decomposed profiles of the surface species that are hydrogen bonded to (a-d) surface O atoms or (e) adlayer O atoms.In Figure S7(e), the threshold to determine the bulk region for each system is also given.Figure S8 presents the profiles of the ion species that were not coordinated to the surface Ce atoms, which were generated only transiently during the proton/hydroxide transfer processes.Note that H3O + ions coordinated to the surface Ce atoms were generated quite rarely, and their density profiles are less than 10 −6 Å −3 .

S.VIII. Supplementary Figures on Proton Hopping Mechanisms
In Figure S9

S.IX. Difference Caused by Proton Assignment Procedure
In this section, we introduce the stable state picture (SSP) [9] assignment for comparison with the nearestneighbor (NNB) assignment of the hydrogen-oxygen bond.While SSP criteria for proton hopping chemical species have already been proposed to analyze population correlation functions, [10] they might be too complicated to apply to a relayed transport analysis, because their strict application to a chain hopping process would require tracking all involved molecules until they are considered to be in stable states.In the following, we propose a simple SSP focusing only on proton states.
Initially at time t = 0, all hydrogens are assigned by using the NNB criterion.At time t > 0, if a hydrogen atom has only one O atom within 1.1 Å and no other O atoms within 1.5 Å, the hydrogen atom is regarded to be stably bonded to the O atom.These thresholds were determined based on the O−H RDF of bulk water (see Figure S10).If a hydrogen atom does not satisfy this criterion, the hydrogen atom is considered to be still bonded to the previously bonded O atom.By using this SSP assignment, we can eliminate counting rapid recrossing around transition states of proton hopping.
A total of 401,647 and 85,747 proton hopping events were observed for 4 ns equilibration NNP-MD simulations of the water/CeO2(111) and (110) interfaces, respectively, under the NNB assignment.These numbers were reduced to 72,381 and 8,546 under the SSP assignment.In Table S3 and Table S4, the number of proton/hydroxide ion transfers via the NNB and SSP assignments are shown, omitting event classes that occur only less than twice during 4 ns simulations of p(8×8) surfaces.Note that each transfer mechanism involves multiple proton hopping events.For reference, the number of transfers that are different from the reverse of the previous transfer (New) and the same as the reverse of the previous transfer (Rattling) are also given.The procedure for identifying "New" and "Rattling" transfers is explained in Section S.XI.While the number of transfers labeled "Rattling" significantly reduced under the SSP assignment, the changes in "New" transfers are relatively moderate.This implies that distinguishing "New" transfers from "Rattling" transfers itself causes eliminating recrossing transfers similar to the SSP.number of PHS occupancies of the corresponding acid sites.To examine this environment dependence of the reaction rates, we show the results of two different models.In the first model (csMSM-a), we omit all environment dependencies while maintaining the semi-Markov nature of the model.In the second model (csMSM-b), we consider not only the number of PES around the base sites but also their spatial configurations to the reaction site pairs.Thus, the csMSM-a and csMSM-b are simpler and more complicated variants of the csMSM, respectively.

S.X. Histogram of Lifetime of Transient Species during SPF/SPR/SPT-II Processes
Figure S12 shows the results of the KMC based on the reaction rates calculated by using the csMSM-a and csMSM-b, including the same results in Figure 6.The KMC based on the csMSM-a does not reproduce the hydroxylation ratio of the surface O sites, and there is no significant improvement in the KMC based on the csMSMb compared to that of the csMSM.Therefore, we concluded that the csMSM contained sufficient environment dependence in this regard.

S.XIII. Environment Dependent Reaction Rates
In this section, we show a part of the raw results of the csMSM analysis.In our cMSM/csMSM models, the reaction rates of reactions involving metal acid sites depend on the number of assigned OH − ions to the sites (Number of hole occupations; #PHO), and those of reactions involving surface oxygen base sites depend on the number of protonated base sites around the sites (Number of surface oxygen protonation, #SOP).In Table S5 and Table S6, the estimated reaction rates are shown, while those of the transfers observed only less than twice during 4 ns simulations of p(8×8) surfaces are omitted.Here, "src." and "dst."represent source and destination of transfers.The values of Figure 6(c) are the averages of raw reaction rates weighted by their frequency of reaction opportunity.
Table S5.Reaction rates of "non-rattling" reactions based on the csMSM for the water/CeO 2 (111) interface.As seen in Section S.XI, these environment dependencies are essential to describe surface hydroxylation.
Mostly, the observed "lateral interaction" was repulsive.For pair formation reactions, denser environments caused lower reaction rates, and the opposite tendency existed for pair recombination reactions.For transfer reactions, if the source site of the proton hole/excess proton was in a dense environment, the interaction promoted the transfer, and vice versa for the transfer destination site.An exception is the SPF reaction of the "nn" pair for the water/CeO2(110) interface.Because a formed OH − could be stabilized by near OSH + , a denser environment appeared to promote pair formation reactions.

S.XIV. Rattling Reaction Rates
Some of the "rattling" part of the reaction rates based on the csMSM are given in Figure S13 and Figure S14.
Although some of the 50% confidence intervals exceed the plot range, they are plotted as they are to show the relative magnitudes between them.

S.XVI. Combination between SSP Assignment and cMSM/csMSM
To investigate the effects of the proton assignment procedure to the reaction rates, we conducted cMSM/csMSM analysis by using the SSP assignment defined in Section S.IX.The results of the KMC simulations for validation are given in Figure S16, and some of the reaction rates based on the csMSM are shown Figure S17.In the case of (c), each curve is plotted thinner as the corresponding survival function becomes less than 20%.

S33
As seen in Figure S16, the hole survival functions significantly changed because holes with short lifetimes caused by rapid recrossing were eliminated by the SSP assignment.However, surface hydroxylation ratios and long-range diffusion properties are similar to those of the NNB assignment.Also, while the short-time behavior of the time-dependent reaction rates was significantly affected, as in Figure S17(c), the long-time behavior and constant terms are relatively similar, as in Figure S17(a), (b), and (d).Therefore, the global tendency of the reaction rate we obtained can be considered essential regardless of the details of the proton assignment.

Figure S17.
The same plots as in Figure 7 but for values obtained by using the SSP assignment.

Figure S1 .
Figure S1.Roughly estimated free energy curve, F(r), as a function of the minimum distance between one of the surface OS atoms and all H atoms of water molecules.The actual bias function with a negative sign and the histograms obtained by the umbrella sampling are also shown with a bin width of 10 −3 Å.

Figure S2 .
Figure S2.(a) Top 5,000 values of the maximum of the SSD of the forces, calculated from the structures of the MD simulation during each refinement step.Ratios of structures exceeding 0.2 eV/Å are also shown.(b) Convergence of the oxygen density distributions.Each density distribution for the initial/refinement step has four curves, which correspond to the four independent trajectories by using different NNPs.The curves are shifted for better visibility, and the density distribution calculated from the product run is also plotted for each shift level.In order to check the convergence in the third refinement for the CeO2(110) system, the results calculated from the fourth refinement step are also shown, while the NNPs from this step are not used in the product run.

Figure
Figure S4 shows the ratio of various surface species as a function of time, starting from undissociated water molecules by NNP-MD simulations for p(8×8) surfaces.The definition of the surface hydroxy groups, O S H + , surface hydroxide ions, O*H − , and adsorbed water molecules, O*H2, are given in the main text.It should be noted that, as our NNPs have been specifically trained to evaluate thermal equilibrium states, the reliability of the relaxation processes far from equilibrium is not guaranteed.

Figure S4 .
Figure S4.Ratio of surface/adlayer species per primitive surface unit cell as a function of time.

Figure
Figure S5 shows snapshots of the trajectories of NNP-MD product run.A water molecule escaped from the water layer into the vacuum region in the snapshot of the water/CeO 2 (110) interface was removed to avoid confusion.

Figure S6 .
Figure S6.The RDFs of the Ce−O pairs.

Figure S7 .
Figure S7.Decomposed density profiles of hydrogen-bonded surface species.Density profiles of total oxygen and hydrogen atoms are also plotted as dotted curves.

Figure S8 .
Figure S8.Density profiles of ion species that are not coordinated surface Ce atoms.Density profiles of total oxygen and hydrogen atoms are also plotted as dotted curves.
Figure S9(a), even if a proton donor/acceptor adlayer molecule desorbs/adsorbs from a Ce site before/after a hopping event, we treat the molecule as an adlayer molecule.Therefore, the case in Figure S9(a) is classified as an APT-I (direct) process rather than an APT-II (solvent-assisted) process.Due to this treatment, while two water molecules cannot adsorb on a single Ce site on the CeO2(111) surface, an APT-I process occurs on a single site (Figure S9(b)).Note that, in the case of the CeO2(110) surface, two water molecules can be coordinated to the same Ce site (see Figure 4(c-ii)).

Figure S10 .
Figure S10.(a) The RDF of the bulk water O−H pair.The bulk water MD trajectory explained in Section 3.5 was used to calculate the RDF.(b) The logarithmic plot of the RDF.

Figure S11 .
Figure S11.Histograms of ΔtPT= tb − ta during SPF/SPR/SPT-II processes.Those during APT-II processes are given in Figure 5, and t a and t b are defined in Figure 3.

Figure S12 .
Figure S12.(a) Probability density distributions of the hydroxylation ratio of the surface O atoms, (b) proton hole survival functions, and (c) mean square displacements (MSD) of the proton holes.

Figure S13 .
Figure S13.Time-dependent reaction rates of "rattling" reactions for the water/CeO2(111) interface.The values are averaged every 100 fs, and environment dependences are also averaged.

Figure S14 .
Figure S14.The same plots as in Figure S13, but for reactions for the water/CeO2(110) interface.

Figure S15 .
Figure S15.Vibration correlation functions of the adlayer and bulk water (a) OH stretching and (b) H 2 O bending.

Figure S16 .
Figure S16.(a) Probability density distributions of the hydroxylation ratio of the surface O atoms, (b) proton hole survival functions, and (c) mean square displacements of the proton holes, obtained based on the SSP assignment.

Table S2 .
Adsorption energies of a water molecule in molecular and dissociative states on the CeO2(111) and (110) surfaces.

Table S3 .
Transfer event frequencies calculated from structures sampled every 5 fs from 4 ns NNP-MD trajectories of the water/CeO2(111) interface.Event classes occurring with a frequency of less than 0.01 ns −1 per primitive surface unit cell (approximately twice during 4 ns simulations of p(8×8) surfaces) are omitted.

Table S4 .
Transfer event frequencies calculated from structures sampled every 5 fs from 4 ns NNP-MD trajectories of the water/CeO2(110) interface.

Table S6 .
Reaction rates of "non-rattling" reactions based on the csMSM for the water/CeO2(110) interface.