Structural and Dynamic Properties of Solvated Hydroxide and Hydronium Ions in Water from Ab Initio Modeling

Both $\rm OH^{-}$(aq) and $\rm H_{3}O^{+}$(aq) ions are amphiphilic in nature as compared to the water molecules. Fig. 1 shows the hydrophobicity and hydrophilicity of different ends in both ions, as well as the number of HBs for each end. The oxygen site of $\rm OH^{-}$(aq) is hydrophilic and the hydrogen site is hydrophobic; on the contrary, the oxygen site of $\rm H_{3}O^{+}$(aq) is hydrophobic while the hydrogen site is hydrophilic. Fig. 2 illustrates the radial distribution functions $\rm g(r)$ of solvated hydroxide and hydronium ions, as well as those of the ambient liquid water from AIMD simulations using both SCAN and PBE0-TS XC functionals. Note that the superscripts $\rm*$ and $\rm w$ denote an atom from an ion and a water molecule, respectively. The error bars of each curve shown in Fig. 2 represent the standard errors obtained by dividing the MD trajectory into several segments with the length of each segment being 2.419 ps (10,000 snapshots); these standard errors of sampling are relatively small compared to the amplitude of the radial distribution functions. Therefore, the errors do not influence the discussion below. On the one hand, the two ions are hydrophobic, as respectively illustrated in Figs. 2 (a) and (f), $\rm g_{H^{*}O^{w}}(r)$ of hydroxide and $\rm g_{O^{*}H^{w}}(r)$ of hydronium both exhibit the reduced HBs of the first solvation shell as compared to ambient liquid water. Specifically, the hydroxide (hydronium) ion shows a strong repulsion against water molecules on its hydrogen (oxygen) site by largely reducing the donated (accepted) HBs with neighboring atoms. On the other hand, the two ions are hydrophilic, which can be seen in Figs. 2(b) and (e) as the first peak of $\rm g(r)$ shifts to the left when compared to its counterpart in pure liquid water. In detail, for the first high peak of $\rm g_{O^{*}H^{w}}(r)$ shown in Fig. 2(b), the oxygen site of the hydroxide ion attracts multiple neighboring water molecules via forming HBs, while the first high peak of $\rm g_{H^{*}O^{w}}(r)$ shown in Fig. 2(e) implies that the hydrogen sites of the hydronium stably donate three protons to neighboring water molecules and form HBs. As will be explained below, the amphiphilic nature of the two ions leads to their unique structural and dynamic properties.

Unlike some GGA functionals that yield qualitatively different radial distribution functions for hydroxide, Tuckerman, Chandra, and Marx (2006) the SCAN functional yields a $\rm g(r)$ similar to that predicted by PBE0-TS, as shown in Fig. 2. Regarding the computational costs, SCAN takes the advantage of being much more efficient than the hybrid functional. Chen et al. (2017) Therefore, it is valuable to compare properties of hydroxide and hydronium in detail as obtained from the two functionals and we have the following findings. First, the hydrophobic sites of the two ions are investigated. We observe a less hydrophobic end in hydroxide (hydrogen site) as compared to that in hydronium (oxygen site), because the first and the second peaks of $\rm g_{H^{*}O^{w}}(r)$ for hydroxide in Fig. 2(a) are considerably higher than those in the $\rm g_{O^{*}H^{w}}(r)$ for hydronium shown in Fig. 2(f), suggesting more water molecules locate around the hydrogen site of hydroxide than the oxygen site of hydronium. When comparing the results from SCAN to those from PBE0-TS, the hydrophobic site of hydroxide, i.e., the hydrogen of hydroxide, becomes less hydrophobic by donating slightly more HBs to its neighbors. This is evidenced as SCAN yields a higher first peak of $\rm g_{H^{*}O^{w}}(r)$ than PBE0-TS in Fig. 2(a). The result is also consistent with the HB data shown in Fig. 1, where the averaged number of donating HBs in hydroxide is 0.577 and 0.471 from SCAN and PBE0-TS, respectively. On the contrary, the hydrophobic site of hydronium, i.e., the oxygen site, becomes more hydrophobic by adopting the SCAN functional, which is evidenced by a lower first peak of $\rm g_{O^{*}H^{w}}(r)$ located around 3 Å in Fig. 2(f). Moreover, the peak locates further away from the central atom than the one from PBE0-TS, suggesting that it would be more difficult for hydronium to accept HBs in SCAN (0.138) than in PBE0-TS (0.140), as listed in Fig. 1.

Second, we focus on the hydrophilic sites of the two ions, where multiple HBs form between the ion and water molecules. Notably, the symmetric picture for the twin topological defects, i.e., hydroxide and hydronium ions, is invalid from both XC functionals if we inspect the coordination number of the two ions. Specifically, the height of the first peak at the hydrophilic end of hydroxide (Fig. 2(b)) is higher than that of hydronium (Fig. 2(e)), indicating a larger coordination number. For example, according to the HB analysis of the SCAN trajectory in Fig. 1, hydroxide accepts 4.101 HBs on average, which is larger than hydronium that donates 2.990 H-bonds. Meanwhile, the oxygen site of hydroxide described by PBE0-TS donates 3.922 HBs, which exhibits less hydrophilic property as compared to the one by SCAN. Furthermore, we observe more non-H-bonded water molecules in the second shell of hydroxide from the SCAN functional, which can be seen by a higher second peak of $\rm g_{O^{*}H^{w}}(r)$ in Fig. 2(b). The above structural features of hydroxide imply that SCAN gathers more H-bonded water molecules around both ends of hydroxide than PBE0-TS.

The amphiphilic nature of the two ions can be further analyzed by the HB statistics of the hydroxide and hydronium ions. As found by previous studiesTuckerman, Chandra, and Marx (2006); Chen et al. (2018), the hydronium cation stably donates three HBs but the hydroxide ion has two H-bonded solvation structures, i.e., a nearly tetrahedral structure with three acceptor HBs and a hyper-coordinated structure with four or more acceptor HBs. While the PBE functional predicts more tetrahedral-like structures of hydroxide, the inclusion of van der Waals interactions and self-interaction correction helps to soften and stabilize the hypercoordination of hydroxide. Chen et al. (2018) Importantly, both SCAN and PBE0-TS functionals include the van der Waals interactions and predict a majority of hydroxide ions that accept four or more HBs to form the so-called hyper-coordinated structure, which are shown in Figs. 3(a) and (e). For instance, 68.3% and 74.7% of hydroxide ions with four accepted HBs are found in the SCAN and PBE0-TS trajectories, respectively. Interestingly, the SCAN functional predicts substantially more hydroxide ions that accept 5 HBs (21.2%) than the PBE0-TS functional (9.0%). Additionally, SCAN and PBE0-TS give rise to 10.2% and 16.1% of hydroxide ions that accept 3 HBs. The above analysis suggests that although both hydroxide ions predicted by SCAN and PBE0-TS tend to form the stable hyper-coordinated solvation structure, the distributions of formed HBs still exhibit some deviations.

Early AIMD simulations predicted the solvation structure of the hydroxide to be planar-like.Tuckerman, Chandra, and Marx (2006) However, a pot-like solvation structure of hydroxide was suggested by the neutron diffraction experiment coupled with empirical potential structure refinement. Botti et al. (2003) Fig. 4 illustrates the probability density plot of the spatial distributions of the first solvation shell oxygen atoms that are H-bonded to the hydroxide ion. The results obtained from the SCAN and PBE0-TS functionals are shown in Figs. 4(a) and (b), respectively. From the figure, we can see that the usage of both functionals leads to a pot-like shape of the solvation structure, which is close to the experimental data. Botti et al. (2003)

In comparison with the pot-like shape of the solvation structure from PBE0-TS, SCAN also provides a pot-like shape of the solvation structure but with a thinner bottom, which is closely related to the number of accepted H-bonds. As previously discussed, SCAN yields a larger portion of hydroxide ions which accept 5 HBs as compared to PBE0-TS. In this regard, we examine the solvation structure features of hydroxide with different coordination numbers as described by the two methods, in order to clarify how different XC functionals affect the solvation structure in detail. Figs. 4(a) and (b) plot the solvation structures of hydroxide in terms of accepting 4 (cyan) or 5 (grey) HBs in liquid water. For both functionals, we find that the pot-like solvation structure with 5 HBs occupies a larger volume as compared to the 4 HBs counterparts, especially on the top and bottom areas of the pot-like solvation structure. This indicates that the accepted fifth oxygen atom locates further away from the hydroxide ion. Interestingly, when accepting its fifth H-bond, the hydroxide described by SCAN tends to locate the HB at the bottom of the ‘pot’, while PBE0-TS tends to place the HB on the top. It is worth mentioning that while the excessive HB does not influence the distribution of donated HB in the SCAN trajectory, it substantially affects the distribution of donated HB in the PBE0-TS trajectory, as Fig. 4(b) shows that the lid of the pot becomes thinner when the hydroxide accepts 5 HBs as compared to 4 HBs.

We further compare the planarity order parameter $\rm p(r)$ to understand the local solvation structure of hydroxide that accepts 4 HBs. The parameter $\rm p(r)$ is defined as the distance from the fourth accepted oxygen atom of hydroxide to the plane formed by the other three oxygen atoms. Chen et al. (2018) Fig. 5 illustrates the $\rm p(r)$ of hydroxide from the SCAN and PBE0-TS trajectories. The two trajectories lead to substantial differences in $\rm p(r)$. For example, $\rm p(r)$ of SCAN exhibits a notable distribution from zero to around 2.0 Å, while $\rm p(r)$ of PBE0-TS owns a peak located at a larger planarity. Essentially, the results suggest that the hydroxide ion described by SCAN gives rise to a more planar distribution of accepted O atoms, while the PBE0-TS softens the directional HB strength by mitigating self-interaction error and forms a pot-like solvation structure. The analysis is in accordance with the above comparison made between the cyan isosurfaces in Figs. 4(a) and (b). Notably, the pot-like solvation structure by SCAN relies more on the relatively abundant 5-coordinated HB structure, which yields a solvation structure with a thinner bottom as compared to PBE0-TS. In regard to the relationship between solvation structure and diffusion, a more planar geometry of the solvation structure of hydroxide deviates more from a tetrahedral HB network of water than PBE0-TS, resulting in a more stable asleep mode of hydroxide.

The hydrogen bond is largely comprised of electrostatic attractions, which is determined by the electronic structures from DFT calculations. In AIMD, the local electronic structure associated with chemical bonds can be conveniently described by the maximally localized Wannier functions (MLWFs) Marzari and Vanderbilt (1997); Marzari et al. (2012), which are obtained by unitary transformation of the DFT eigenfunctions. To be specific, in a molecular dynamics trajectory, the information for the electronic environment for water molecules can be deduced from the statistical analysis of the MLWF centers. Figs. 6(a), (b), and (c) show the distributions of centers of MLWFs representing electronic information in $\rm OH^{-}$(aq), $\rm H_{2}O$, and $\rm H_{3}O^{+}$(aq) molecules in ambient liquid water, respectively. The four Wannier function centers surrounding a water molecule or water ion are divided into two peaks according to their distances from the central oxygen atom. The peak located closer to the O atom represents the lone pairs of electrons while the peak further from the O atom represents the bonding pairs of electrons constituting covalent bonds. The centers of MLWFs can be utilized to compute the molecular dipole moments with the following formula Sharma, Resta, and Car (2005); DiStasio et al. (2014); Zheng et al. (2018):

where $\mathbf{R_{H_{1}}}$, $\mathbf{R_{H_{2}}}$ and $\mathbf{R_{O}}$ are the coordinates of the three atoms in a water molecule, and $\mathbf{R_{W_{i}}}$ with $i=1,2,3,4$ are coordinates of the four MLWF centers. In this work, we calculate the distributions of molecular dipole moment magnitudes of water molecules involved as donor (acceptor) on the hydrophilic end of $\mathrm{OH^{-}}$ ($\mathrm{H_{3}O^{+}}$). Fig. 7 illustrates the results from both SCAN and PBE0-TS functionals, as well as the dipole moments for pure liquid water obtained from other works. DiStasio et al. (2014); Zheng et al. (2018) Note that for the hydrophobic sites of the two ions, we observe that the neighboring water molecules are less affected by the presence of water ions and their dipole moments are close to those in bulk water.

The positions of the lone pairs and bonding pairs of electrons determine the amphiphilic nature of molecules. Taking the SCAN trajectories as an example, Figs. 6(a) and (b) illustrate that the length of covalent bond in hydroxide and ambient liquid water, i.e., the distance between the oxygen and the proton, is 0.97 and 0.98 $\rm\AA$, respectively. Meanwhile, the distance between the oxygen atom and the covalent bonding pairs in hydroxide and in liquid water is 0.52 and 0.49 $\rm\AA$, respectively. Therefore, the distance between the proton and the electron bonding pairs is smaller in hydroxide (0.45 $\rm\AA$) than in liquid water (0.49 $\rm\AA$), leading to a less positively charged environment of hydrogen in hydroxide. The above explains the origin of the hydrophobic site of the hydroxide anion, where the formation of HBs is hindered due to the weak electrostatic attraction force. On the other hand, the lone pair of hydroxide locates further away from the oxygen atom (0.35 $\rm\AA$) than the lone pairs of liquid water molecules (0.33 $\rm\AA$), which brings the hydroxide ion a more negatively charged local environment around its oxygen site. Consequently, this leads to the hydrophilic property of hydroxide anion, where acceptance of HBs becomes easier than water molecules. As illustrated in Fig. 7(a), in correspondence to the more negatively charged oxygen site, the water molecules accepted by hydroxide become more polarized with an average dipole moment of 3.17 D, which is substantially larger than the average dipole moment of 2.94 D in bulk water. Meanwhile, the average dipole moment of the H-bonded water on the hydrophobic end is only 2.93 D, which is close to the average value of 2.94 D from pure liquid water. In stark contrast, the electron bonding pairs of hydronium cation stay further away from the protons, creating a more positively charged (hydrophilic) environment to stably form three HBs, as illustrated in Fig. 6(c); the oxygen site of hydronium stays closer to the lone pair, leading to a less negatively charged environment that exhibits hydrophobic properties. Induced by the strong positively charged hydrogen site, the dipole moment of water molecules accepting HBs from hydronium also increases significantly. As shown in Fig. 7(b), the average dipole moment of water molecules accepting hydrogen bonds from hydronium increses from 2.94 D in bulk water to the value of 3.17 D. Meanwhile, the average dipole moment of the water molecule that accepts a hydrogen bond from the oxygen site of hydronium reaches 2.96 D, similar to the value in bulk water.

Now we compare the detailed electronic structure differences of the two ions from SCAN and PBE0-TS. First, as shown in Figs. 6(a), (b), and (c), a larger distance between oxygen and proton (the length of O-H covalent bond) suggests that oxygen binds more loosely with the proton, and SCAN is known to strengthen the covalent bonds as compared to PBE Chen et al. (2017); here we find PBE0-TS yields protons that locate closer to the oxygen atom, implying a stronger covalent bond as compared to SCAN. Besides, for the pure water described by the SCAN functional, both lone pairs and bonding pairs of ambient liquid water molecules shown in Fig. 6(b) shift closer to the oxygen atom. This indicates that SCAN provides a less negatively charged environment on the oxygen site and a more positively charged environment on the proton side for pure liquid water. In other words, the oxygen and proton sites of water molecules respectively become more hydrophobic and more hydrophilic in the SCAN trajectory than in the PBE0-TS counterpart. The two competing effects lead to a more pronounced first peak of radial distribution functions as observed in Figs. 2(c) and (d), which implies that SCAN provides a more strengthened HB at a short-range scale when compared to PBE0-TS. In addition, when comparing the electronic structure of pure liquid water, the distribution of MLWF centers shown in Fig. 6(b) becomes more delocalized for the PBE0-TS trajectory than the SCAN trajectory, which is consistent with the effect of mitigating the self-interaction errors as provided by PBE0-TS. However, the differences of MLWF distribution for the two ions, as inferred from Figs. 6(a) and (c) by the two functionals, are not substantial.

Second, for the hydroxide anion illustrated in Fig. 6(a), the location of the electron bonding pair is slightly further away from the central oxygen, as predicted by SCAN and compared to that from PBE0-TS. Nevertheless, the proton locates further away from the bonding pair as described by SCAN, suggesting a stronger hydrophilicity of the originally hydrophobic proton site in hydroxide. This leads to an increased number of donated HBs from PBE0-TS (0.471) to SCAN (0.577), as abovementioned. For the other site of the hydroxide anion (oxygen), when comparing the results from SCAN to those from PBE0-TS, the former provides a shortened distance between the lone pair and the oxygen atom, creating a more hydrophobic oxygen site of hydroxide. However, a competing effect arises due to the usage of SCAN, i.e., the protons of water molecules become more hydrophilic. As a result, the oxygen site of hydroxide in the SCAN trajectory accepts more HBs (4.101) than that of PBE0-TS (3.922), as listed in Fig. 1. The above differences that exist in the electronic structures of hydroxide anion lead to the deviation of the solvation shell structure obtained from the SCAN and PBE0-TS trajectories, as illustrated in Fig. 4. As for the average dipole moment of water molecules that donate H-bonds to the hydroxide ion, Fig. 7(a) shows that the values are 3.17 and 3.19 D from the SCAN and PBE0-TS exchange-correlation functionals, respectively. We conclude that the two functionals yield similar dipole moments for the water molecules at the hydrophilic end of the hydroxide ion.

Third, by analyzing the distribution of MLWF centers of hydronium cation from the SCAN and PBE0-TS trajectories shown in Fig. 6(c), we observe that SCAN provides a more hydrophilic environment for the proton sites of hydronium because the distance between the bonding pair and the protons is larger. In addition, the protons of hydronium interact with the more hydrophobic oxygen of water molecules. As a result, hydronium from both trajectories donates three stable HBs, suggesting that SCAN and PBE0-TS yield similar short-ranged structure information for the proton sites of hydronium. In addition, Fig. 7(b) shows that the two functionals predict similar dipole moments for the water molecules that accept H-bonds from the hydronium ion; Specifically, the value is 3.17 and 3.11 D from SCAN and PBE0-TS, respectively. On the other hand, for the oxygen site of hydronium, the lone pair from SCAN locates slightly closer to the oxygen atom as compared to the one from PBE0-TS, which creates a slightly more hydrophobic environment for hydronium. In this regard, it is consistent that hydronium cation accepts 0.138 HBs from SCAN, which is marginally smaller than the 0.140 HBs from PBE0-TS, as listed in Fig. 1.

In conclusion, for water ions including the hydronium and hydroxide ions, as well as the water molecules, we analyze their electronic structures via the distribution of MLWF centers. In general, compared to the results from PBE0-TS, we find that the SCAN functional yields a closer distance between the lone pairs and the oxygen atom, and a larger distance between the bonding pairs and protons. In addition, the hydrophilic and hydrophobic nature of the two ions are affected by considering the delicate interactions between the ions and neighboring water molecules. Consequently, when comparing the SCAN results with the PBE0-TS results, the hydroxide anion from SCAN becomes more hydrophilic, while the hydronium cation from SCAN becomes slightly more hydrophobic.

Recently, the PT mechanism of the two ions has received increasing attention from simulations. For example, the diffusion coefficients of the two ions are largely determined by the PT phenomenon. In fact, only part of the PT events contributes to the final diffusion coefficients, the rest are usually classified as non-diffusional PTs such as rattling. Chen et al. (2018); Chandra, Tuckerman, and Marx (2007); Fischer, Dunlap, and Gunlycke (2018); Tse, Knight, and Voth (2015); Arntsen et al. (2021) However, defining the rattling events are non-trivial. Fischer, Dunlap, and Gunlycke (2018); Arntsen et al. (2021) Until now, several classification methods of PT events exist due to the complex behavior of the PT phenomenon itself. Chandra, Tuckerman, and Marx (2007); Tse, Knight, and Voth (2015); Fischer, Dunlap, and Gunlycke (2018); Hassanali et al. (2013) For instance, Hassanali et al. Hassanali et al. (2013) proposed to divide the PT events into single and concerted jumps, the latter of which dominate the PT phenomenon for both hydronium and hydroxide ions. Later, the concerted jump mechanism was revised to be dominant only in hydronium by using a more accurate XC functional PBE0-TS in AIMD simulations. Chen et al. (2018) Recently, Arntsen et al. suggested that the concerted jumps of hydronium were in fact dominated by non-diffusional jumps. Arntsen et al. (2021) In addition, it is worth mentioning that the ‘special pair dance’ (SPD) concept, which was proposed for hydronium, Markovitch et al. (2008) refers to the period during which no PT occurs and the closest H-bonded neighbor of the protonated water molecule quickly interchanges at a frequency of tens of femtoseconds. In this regard, the SPD concept has no overlap with the classification of PT. As a typical phenomenon of the non-diffusional PTs, water ions tend to be trapped among a few water molecules and keep rattling between them, as explained by previous studies using either ring distributionHassanali et al. (2013) or accepted HB of hydronium. Tse, Knight, and Voth (2015) We also observe similar trapping phenomena in the SCAN trajectories. Figs. 8(a) and (b) respectively illustrate the representative PT events for hydroxide and hydronium from the SCAN trajectories, where the change of ion indices with respect to simulation time is plotted. The trapping events are identified as blue circles in Fig. 8, from which the proton (hole) tends to hop among neighboring water molecules with a high frequency for more than a few ps but contribute little to the diffusion coefficients.

In order to remove the non-diffusional PTs, Arntsen et al. Arntsen et al. (2021) proposed to apply a filter to the molecular dynamics trajectories. The filter eliminates the short-ranged hopping events that are followed directly by a back hop returning to the original water. In other words, two molecules were used to define the return events. The method can also be applied to trim the “slingshot" effect which involves more than two water molecules but cannot fully eliminate the non-diffusional PTs. Fischer et al. Fischer, Dunlap, and Gunlycke (2018); Fischer and Gunlycke (2019) inherently defined the PT events within the first solvation shell of hydronium as rattling, which typically involves three or four water molecules within the Eigen complex $\mathrm{H_{3}O^{+}[H_{2}O]_{3}}$. We find similar rattling effects in the SCAN trajectory. Interestingly, we observe additional long-ranged return events involving more than four water molecules, which do not have to rely on the concept of Eigen complex. Specifically, the long-ranged return phenomenon refers to a water ion hops for several steps and occasionally revisits its former hosts by following exactly the same route. In the following, we refer to these long-ranged return events that involve three or more water molecules as the ‘revisit’ events.

Fig. 8 illustrates multiple revisit routes (orange arrows) as obtained from the SCAN trajectories, which suggest that non-diffusional events are in fact a complex phenomenon that involves several water molecules. As illustrated in Fig. 8(b), several revisit routes of hydronium involve three water molecules, most of which result from frequent switches of proton rattling partners. For example, in the first revisit event of hydronium, the proton first rattles between water molecules with indices being 3 and 42, then between water 3 and 55, followed by rattling again between water 3 and 42. Another long-ranged revisit event is shown in Fig. 9, where the revisit route is shown in Fig. 9(a) while the changes of ion index and distance from the initial position with respect to time are plotted in Fig. 9(b). We can see that after the proton hop from the water molecule with index of 9 to another one with the index of 37, the proton takes a long journey by visiting five different water molecules with indices being 37, 58, 33, 5 and 7. The time of starting and ending points is respectively marked as red and orange points in Fig. 9(b). After arriving at water with index of 7, the excess proton rapidly goes back to water with index of 37 by following the same route ($7\rightarrow 5\rightarrow 33\rightarrow 58\rightarrow 37$). The return point is marked as a green point in Fig. 9(b). By observing the change of distance in Fig. 9(b) after the revisit event ends, we find these PT events contribute little to the diffusivity of hydronium. In summary, the existence of the above two non-diffusional PT phenomena, i.e., trapping and revisit, demonstrates that the classification of PT events needs more investigations in future.

Importantly, studying both non-diffusional and diffusional PTs from simulation is informative based on the following three considerations. First, non-diffusional PT events are relevant to experimental findings. For instance, the Eigen-Zundel configuration transition of hydronium is explained by fast rattling between two water molecules, Zundel and Metzger (1968); Eigen, De Maeyer, and Bernal (1958); Agmon (1995) and recent infrared spectroscopy experiments provide new evidences to identify which species dominates the PT process. Fournier et al. (2018); Carpenter et al. (2018) Second, both non-diffusional PTs and the SPD phenomenon are the intrinsic properties of the two ions that need to be addressed from a scientific point of view. A recent work found different PT patterns exist in solvated hydronium and hydroxide ions, i.e., hydroxide experiences much less PT events than hydornium by considering both diffusional and non-diffusional PT events. Chen et al. (2018) In this regard, it is worth investigating more microscopic origins to elucidate the different PT patterns of the two ions, which could be useful in designing relevant experiments. Third, Tse et al. Tse, Knight, and Voth (2015) mentioned that the frequency of concerted hops involving three water molecules depends on the exchange-correlation functional. Therefore, in order to verify the accuracy of SCAN functional in describing both solvated hydronium and hydroxide ions, it is a necessity to carry out comprehensive tests on both structural and dynamic properties, including both diffusional and non-diffusional PTs.

In this work, we study both diffusional and non-diffusional PTs for the two ions. Essentially, a conceptually correct picture of the structural diffusion of $\rm OH^{-}$(aq) and $\rm H_{3}O^{+}$(aq) ions relies on a precise description of the HB network of water, as well as the solvation structures of water ions. Although the issues are well addressed by the PBE0-TS functional, Chen et al. (2018) the computational cost is extremely high. As abovementioned, the number of formed HBs between molecules stands out as a major indicator of hydrophobicity and hydrophilicity. Generally dictated by electronic structures, the HB network of water and ions is mainly determined by the used XC functional. Laasonen, Csajka, and Parrinello (1992); Laasonen et al. (1993); Gillan, Alfè, and Michaelides (2016); Gaiduk, Gygi, and Galli (2015); Schmidt et al. (2009); Wang et al. (2011); Miceli, de Gironcoli, and Pasquarello (2015); Chen et al. (2017); DiStasio et al. (2014); Burke, Perdew, and Wang (1998); Boese et al. (2000); Tuckerman, Chandra, and Marx (2006); Sakti, Nishimura, and Nakai (2020) Therefore, we validate the performance of the SCAN XC functional in providing the frequency of the PT events and compare them to the PBE0-TS results.

In order to quantitatively describe the frequency of PT events in these AIMD trajectories, we categorize the PT events according to the number of transfer events during one burst, as done in Ref. Chen et al., 2018. We define a single jump as a PT event not followed by any PT events in the following 0.5 ps. The 0.5 ps is a characteristic time to compress the water wire, as introduced in Refs. Hassanali et al., 2011. If more than one PT event occurs within 0.5 ps and the proton does not return to any of its former host molecules in these jumps, these PT will be labeled as a concerted jump. The time cutoff of 0.5 ps is picked on an empirical basis to facilitate comparison between XC functionals, while more physically meaningful time cutoff could be deduced from the PT correlation function. Luzar and Chandler (1996); Chandra, Tuckerman, and Marx (2007); Berkelbach, Lee, and Tuckerman (2009) Note that the rest of PT events as labelled as simple rattling events, which are a subset of rattling events from other definitions that involve two Arntsen et al. (2021) or more molecules. Fischer, Dunlap, and Gunlycke (2018) As abovementioned, none of these definitions could fully eliminate the non-diffusional PTs. In addition, most definitions are applied to hydronium only. In this regard, further investigations for defining rattling for both hydronium and hydroxide ions are needed.

Figs. 10(a) and (c) display statistics of the concerted PT events for hydroxide by using the SCAN and PBE0-TS functionals, respectively. Generally, both XC functionals reveal that the total frequency of PT in hydroxide is lower than hydronium. We also find that the two functionals yield the same conclusion that PT events in the hydroxide system are dominated by both single and double jumps, while triple and quadruple jumps are relatively rare. In stark contrast, as illustrated in Figs. 10(b) and (d), where both functionals yield the same conclusion that the concerted PT events for hydronium is dominated by the double jumps, the frequency of which is substantially larger than the counterpart in hydroxide. For instance, the frequency of double jumps from SCAN (PBE0-TS) is 12.0 (22.1) times per 10 ps for hydronium, much larger than the value of 4.1 (3.1) for hydroxide. In addition, the triple- and quadruple-jumps are relatively rare in hydronium.

We define a forward jump of proton as jumping to a neighboring water molecule that is not the host in the previous jump. In this regard, for any PT events that concern two or more jumps, the second jump must be a forward jump before it can be classified into a multiple jump. It is useful to collect the waiting time between the two successive jumps for the PT events of the two ions, and the results are shown in Fig. 11. The data in Fig. 11 are collected from the SCAN trajectories but similar results can be deduced from the PBE0-TS trajectories. As Fig. 11 shows, 76% and 58% of the forward jumps occur within 0.5 ps after the first jumps for hydronium and hydroxide, respectively. Besides the water wire compression time of 0.5 ps as suggested in Ref. Hassanali et al., 2011, the analysis here also implies that 0.5 ps is a reasonable time scale to include most forward jumps for the two solvated ions in liquid water. Additionally, we observe the hydroxide ion tends to wait longer before the second jump than the hydronium ion. For instance, the longest waiting time is 7.1 and 1.7 ps for the hydroxide and hydronium ions, respectively. These results demonstrate the inactive nature of hydroxide compared to hydronium, which agrees with the different PT frequencies of the two ions as illustrated in Fig. 10.

Since only the simplest rattling event is included in this classification, i.e., rattling between two water molecules, the single and multiple PTs defined here could contribute little to diffusion. We find that hydronium undergoes more simple rattling events than hydroxide, which could be attributed to the fast alternate between the Zundel ($\rm H_{5}O_{2}^{+}$) and Eigen ($\rm H_{9}O_{4}^{+}$) configurations. Marx, Chandra, and Tuckerman (2010); Marx (2006) For example, we observe the hydronium ion rattles between two water molecules continuously for at most more than 30 times within a time period of approximately 5 ps; meanwhile, the hydroxide ion in the SCAN trajectory rattle for at most 8 times between two water molecules for no longer than 2 ps. This result demonstrates a major difference in PT pattern between hydroxide and hydronium.

We find SCAN yields a substantially fewer number of double, triple, and quadruple PT hops than PBE0-TS, as shown in Figs. 10(b) and (d). For instance, the frequency of double jumps from SCAN (12.0) is about half of that from PBE0-TS (22.1). Consequently, we obtain a relatively lower diffusion coefficient of hydronium from SCAN as compared to those from PBE0-TS. The reason can be analyzed via the three water molecules (form a water wire) involved in the double jump. In fact, the third water molecule most likely comes from the second nearest neighbor of the first water molecule.

In this regard, we plot $\rm g_{O^{*}O^{w}}(r)$ from the two trajectories (SCAN and PBE0-TS) in Fig. 12. First, the two functionals yield similar average dipole moment (3.17 D from SCAN and 3.11 D from PBE0-TS, in Fig. 7) of the first-shell neighboring water molecules of hydronium; the increased dipole moments as compared to water originate from the hydrophilic end of the hydronium ion, forming three strong H-bonds that enable the PT events. Second, we can see that the first peak of $\rm g_{O^{*}O^{w}}(r)$ from PBE0-TS is closer to the central oxygen atom of hydronium as compared to SCAN. This indicates that SCAN gives rise to a slightly more hydrophobic hydronium than PBE0-TS, in agreement with previous analysis. As a result, the hydronium described by PBE0-TS is more likely to donate its protons to its neighboring molecules, resulting in more frequent PT events. Third, the second peak from the SCAN trajectory in Fig. 12 is significantly lower than that of PBE0-TS, suggesting a relatively lower local water density in the second solvation shell of hydronium. The lack of enough water molecules in the second solvation shell of hydronium prevents the protons from jumping two or more times within a burst of PT events.

In the end, we would like to discuss the presolvation structure with the SCAN and PBE0-TS XC functional. Past researches have addressed s strong correlation between the PT phenomenon and the presolvation structure of ions, which refers to the HB structure of an ion or molecule in resemblance to the species into which it is turning after PT. Tuckerman, Chandra, and Marx (2006) For example, hydroxide is predicted to convert from a hyper-coordinated structure to a three-coordinated structure and donate a HB before PT. Tuckerman, Marx, and Parrinello (2002); Tuckerman, Chandra, and Marx (2006); Marx, Chandra, and Tuckerman (2010) We follow the method of Ref. Berkelbach, Lee, and Tuckerman, 2009 and plot the $\mathrm{g_{O^{*}H}(r)}$ and $\mathrm{g_{H^{*}O}(r)}$ in time periods of 530-480 fs, 210-160 fs, 100-50 fs and 50-0 fs before the PT event occurs. On the one hand, the RDF and the integrated coordination numbers of hydroxide obtained from both SCAN and PBE0-TS trajectories are displayed in Fig. 13. We observe the first peak of $\mathrm{g_{O^{*}H}(r)}$ moves leftward and decreases in height when the trajectory time approaches the PT event, resulting in the decrease in the coordination number (Figs. 13(a) and (b)). For example, from the time period of 530-480 fs to 50-0 fs before PT, SCAN (PBE0-TS) predicts the coordination number at 2.5 $\mathrm{\AA}$ decreases monotonically from 5.07 (4.90) to 4.65 (4.44), indicating a significant loss in the accepted HB of hydroxide before PT. The first peak of $\mathrm{g_{H^{*}O}(r)}$, however, does not render significant changes with respect to time by either functional (Figs. 13(c) and (d)). On the other hand, previous works have proposed that hydronium tends to accept a HB before PT. Meanwhile, the nearest neighbor of hydronium which is going to accept the excessive proton tends to break the other accepted HB. Berkelbach, Lee, and Tuckerman (2009); Biswas et al. (2016); Marx, Chandra, and Tuckerman (2010) In this regard, we plot $\mathrm{g_{O^{*}H}}(r)$ and $\mathrm{g_{O^{{}^{\prime}}H}}(r)$ in time periods of 530-480 fs, 210-160 fs, 100-50 fs and 50-0 fs before the PT, which are shown in Fig. 14. Here $\mathrm{O^{{}^{\prime}}}$ is used to denote the oxygen atom of water molecule that accepts the excess proton. We notice that the two functionals slightly diverge on the change of the coordination number of $\mathrm{g_{O^{*}H}}(r)$ from 530-480 fs to 50-0 fs before PT: the PBE0-TS functional predicts an increasing coordination number at 2.5 $\mathrm{\AA}$ from 3.22 to 3.49, while the SCAN functional only predicts an increase from 3.29 to 3.36. Furthermore, both functionals show similar decrease of the accepted coordination number of the water molecule that accepts the excess proton: the coordination number of $\mathrm{g_{O^{{}^{\prime}}H}}(r)$ at 2.5 $\mathrm{\AA}$ from the SCAN (PBE0-TS) trajectory decreases from 3.77 (3.96) to 3.54 (3.79). In summary, we find the presolvation picture is mostly valid for both functionals, except for the donation of HB by hydroxide before PT. By using the SCAN and PBE0-TS functionals, we observe the occurrence of PT between a hyper-coordinated hydroxide and its neighboring water molecule or between a hydronium without accepted HB and a water with two accepted HB. In this regard, it would be worth exploring the more detailed PT mechanisms with different functionals in future works. Importantly, the reorganization of HB network before PT might be informative to the fast dynamics captured by the fast-developing infrared spectroscopy. Roberts et al. (2009, 2011); Thämer et al. (2015); Napoli, Marsalek, and Markland (2018); Carpenter et al. (2018); Yuan et al. (2019) In spite of this, the validation of presolvation structure with SCAN-based AIMD simulations against spectroscopic experiment awaits future studies.

We used the Einstein equation to compute the diffusion coefficient of ions as

where the mean square displacement (MSD) is defined as $\|\mathbf{r}(t)-\mathbf{r}(0)\|^{2}$ with $\mathbf{r}(t)$ depicts the atom position at time $t$. The results are listed in Table 1. In order to obtain a reasonable value of $D$, we split the trajectories into multiple segments with a length of 12 ps for each segment and the initial snapshots were displaced by an interval of 3 ps. Fig. 15 shows the resulting MSDs of solvated hydroxide and hydronium ions, which were taken from these slices of trajectories and averaged. Here we find that the diffusion of hydroxide is substantially slower than hydronium, in agreement with the experimental data and previous analysis. To be specific, the diffusion coefficients of both water ions from SCAN trajectories are relatively smaller than those of PBE0-TS. As listed in Table 1, SCAN (PBE0-TS) predicts $D$ of hydronium and hydroxide to be $\rm 5.7~{}(8.3)\times 10^{-9}m^{2}/s$ and $\rm 2.9~{}(3.7)\times 10^{-9}m^{2}/s$, respectively. Based on previous analysis of the structural, electronic, and dynamic properties of the two ions, we infer that the relatively lower $D$ values from SCAN originate from the different HB networks described by SCAN and PBE0-TS. In detail, the oxygen site of hydroxide described by SCAN is more hydrophilic and the hydroxide anion is more often solvated by a stable hyper-coordinated structure with four or even five HBs connected with neighboring water molecules. The hyper-coordinated structure turns the ion into an asleep mode and hinders the PT events, preventing itself from fast structural diffusion. As a result, the PT of hydroxide prefers the stepwise motion. On the other hand, the solvated hydronium described by SCAN has a relatively smaller hydrophilicity than PBE0-TS. However, the stable three HBs still form between the protons of hydronium and neighboring water molecules. In this regard, the concerted PT character is preserved for hydronium. Notably, a lack of the second neighbors further prevents concerted PTs from occurring, which substantially lower the frequency of PT bursts in hydronium. Consequently, the diffusivity of hydronium is lowered as compared to the result from PBE0-TS. In conclusion, the diffusivities of both hydroxide and hydronium ions are lowered from SCAN as compared to those from PBE0-TS, and both functionals predict that hydroxide diffuses slower than hydronium, in agreement with the experiments.

Furthermore, we calculated the orientation time correlation function (TCF) of the two ions with the SCAN functional, $\mathrm{C_{2}(t)=<P_{2}(\mathbf{u}(0)\cdot\mathbf{u}(t))>}$, where $\mathrm{\mathbf{u}(t)}$ denotes a unit vector fixed on the ion at time $\mathrm{t}$, $\mathrm{P_{2}(x)}$ denotes the second-order Legendre polynomial and $\mathrm{P_{2}(x)=(3x^{2}-1)/2}$. Ma and Tuckerman (2011) The unit vector is picked as the vector along the O-H covalent bond for hydroxide, while it is picked as the normalized sum of unit vectors along all of the three O-H covalent bonds for hydronium. The first 2 ps components of the generated TCF are shown in Fig. 16, together with fitted curves of the biexponential function, $\mathrm{C(t)=ae^{t/\tau_{1}}+(1-a)e^{t/\tau_{2}}}$, where $\mathrm{\tau_{1}}$ and $\mathrm{\tau_{2}}$ denote the time scale of the two exponential decays respectively and $\mathrm{a}$ denotes the proportion of the first component. The orientation TCF of both ions exhibit similar biexponentiality as compared to that of water molecule.Laage et al. (2011) The characteristic time of the first component ($\mathrm{\tau_{1}}$) of both ions are similar in magnitude, indicating that the libration mode of water molecules is largely reserved in the ions. However, the two ions yield different relaxation time for the second component ($\mathrm{\tau_{2}}$). In detail, the TCF of hydronium with a characteristic time $\mathrm{\tau_{2}}$ = 2.22 ps decays substantially faster than that of hydroxide ($\mathrm{\tau_{2}}$ = 3.12 ps), which implies the inactiveness of hydroxide compared to hydronium. Therefore, we confirm that the second relaxation process could be primarily due to the reorientation induced by PT as originally proposed by Ma et al. Ma and Tuckerman (2011) In this way, the difference in the characteristic time $\mathrm{\tau_{2}}$ between the two ions could be explained by the different PT rates between them.