Proton-transfer dynamics in ionized water chains using real-time Time Dependent Density Functional Theory

In density functional-theoretic studies of photoionized water-based systems, the role of charge localization in proton-transfer dynamics is not well understood. This is due to the inherent complexity in extracting the contributions of coupled electron-nuclear non-adiabatic dynamics in the presence of exchange and correlation functional errors. In this work, we address this problem by simulating a model system of ionized linear H-bonded water clusters using real-time Time Dependent Density Functional Theory (rt-TDDFT)-based Ehrenfest dynamics. Our aim is to understand how self-interaction error in semilocal exchange and correlation functionals affects the probability of proton transfer. In particular, we show that for H-bonded (H$_2$O)$_n^+$ chains (with $n>3$), the proton-transfer probability attains a maximum, becoming comparable to that predicted by hybrid functionals. This is because the formation of hemibonded-type geometries is largely suppressed in extended H-bonded structures. We also show how the degree of localization of the initial photo-hole is connected to the probability of a proton-transfer reaction, as well as to the separation between electronic and nuclear charge. These results are compared to those obtained with adiabatic dynamics where the initial wavefunction is allowed to relax to the ground state of the ion cluster, explaining why different functionals and dynamical approaches lead to quantitatively different results.


I. INTRODUCTION
Understanding charge transfer dynamics in aqueous solvated semiconductor surfaces for enhanced photocatalytic activity is a complex undertaking. Elucidating the structure of these hydrated surfaces has been a subject of extensive research not only because of the ubiquitous nature of water but also to better control and improve the (photo)catalytic property of the adherent semiconductor materials. Experimentally, even the characterization of the structural motifs and physisorbed chemical species on the semiconductor surfaces is both difficult and uncertain; for example, the identification of surface species involved in photocatalysis is often based on indirect evidence 1 . In particular, the distinction among oxygen atoms belonging to a (surface-bound) hydroxyl group, an oxidated surface and a peroxide compound is, at best, qualitative. Such species are intermediates in the half water-splitting oxidation reaction 2 , and their identification is linked to a mechanistic understanding of how photo-excited holes are transferred to the surfaceadsorbed species and how protons move away from the reaction sites. In general for water-oxide interfaces, the identification of reactive species -electrons and holesand their distribution in the bulk versus that at the surface, pose an uphill challenge for uncovering the underlying reaction mechanisms 3,4 .
From a theoretical viewpoint, ab initio molecular dynamics studies have led to a significant understanding of * vidushi.sharma@stonybrook.edu † maria.fernandez-serra@stonybrook.edu dynamical processes at the (photo)catalytic interface [5][6][7][8][9] . These include a description of the state of adsorbed water -molecular or dissociated -favored on a TiO 2 anatase surface 10,11 , as well as the facet-dependent behavior of excess electrons that makes one surface-termination better suited for a reduction reaction and the other better suited for oxidation 12 . There remain even more unresolved puzzles. For instance, whether or not the transfer of a (photogenerated) hole to surface-water species leads to a fast separation of the proton remains an open question. Furthermore, the time scale of the proton-transfer reaction and whether it ought to be studied at the adiabatic or non-adiabatic level of theory [13][14][15][16][17] is also still an open question. These questions are difficult to attack largely due to the inherent complexity of the systems involved, in particular due to large system sizes, the necessity for long simulation times, and potential contributions from non-adiabatic terms.
A principal aim of our work is to develop physical intuition for a more tractable system comprising of water molecules and examine proton-transfer probability (P(PT)) as computed by rt-TDDFT-based 18, 19 Ehrenfest dynamics 20-23 upon "photoionization" of such systems. The main idea behind this study is to identify a practical and reliable method to study more complex systems. In the condensed phase, water molecules are organized in large cluster structures, strongly connected by hydrogenbond networks 24,25 . Frequently (aqueous) solvated semiconductor surfaces sustain and enhance the formation of such networks 11,[26][27][28][29] , resulting in the ordered formation of surface chains of water molecules. Thus the very welldeveloped network of H-bonds in water plays a fundamental role in not only defining its structural and dynam-ical properties, but also the properties of the material it interacts with 30,31 .
We also want to understand the effect of hydrogen bonds and their cooperativity -if any -on hole localization and proton transfer rate, and simultaneously evaluate the reliability of a commonly used semilocal generalized-gradient approximate (GGA) exchange and correlation (XC) density functional, namely PBE 32 , on describing these processes. We focus on openended chains of water molecules such that upon photoionization, the hole is always localized at the first H-bond donor oxygen of the chain.
Since removing an electron from a system comprising of n H 2 O molecules makes it extremely prone to self-interaction error (SIE) by conventional standards 33 , we first address this issue before conducting a (TD)DFT-based dynamics of ionized water clusters. We compute the groundstate static binding energies of H 2 O in ionized water clusters, (H 2 O) + n , using various approximate (semilocal and hybrid) density-functionals and higher accuracy wavefunction-based methods. A comparison of these independent binding energies yields an estimate of the underlying SIE inherent to the XC functional used. SIE manifests in the deviation of the binding energy computed using a given XC functional from its value using a more accurate method. We show that SIE in functionals like PBE (GGA) become less severe for condensed phases. This is a crucial observation because most assessments of SIE in PBE among other XC functionals have targeted smaller molecular systems or systems with non-interacting components [34][35][36] . The first part of our work, outlined in Sec. II, addresses this gap by examining the error for an interacting system, (H 2 O) + n , which is closer in nature to the bulk phase.
In Sec. IV, we present simulations of the dynamics of (H 2 O) + n systems using (TD)DFT techniques. Our aim is to understand the proton-transfer mechanism and identify the relevant time scales. We present an rt-TDDFT/PBE-based Ehrenfest trajectory statistics of the proton-transfer in H-bonded ionized-water chains of size n = (2 − 5). We also compare the non-adiabatic Ehrenfest dynamics with adiabatic Born-Oppenheimer molecular dynamics, and further refine the analysis by distinguishing between PBE and the hybrid functional PBE0 37,38 .
Finally in Sec. IV B, we explore the dynamical evolution of the photo-hole in the ionized water cluster. The rt-TDDFT/Ehrenfest approach involves a simultaneous real-time evolution of the electronic subsystem in accordance with the time-dependent Kohn-Sham equations and a classical evolution of nuclei. We investigate explicitly the evolution of hole densities and establish a connection between hole-localization and initiation of proton-transfer dynamics in ionized water chains. We also compare the behavior of three density functionals (with increasing fractions of exact exchange) for calculations of the dynamic hole densities on a particular H 2 O unit of the ionized water chain. Our results indicate that for the larger clusters, PBE simulations are not substantially different from those obtained with hybrid functionals. We explain why PBE produces significantly different results in smaller clusters, with our findings supporting this observation as well as the use of other semilocal XC functionals in rt-TDDFT/Ehrenfest simulations of larger, condensed phase systems.

II. SELF-INTERACTION ERROR IN IONIZED WATER CLUSTERS
Local and semilocal generalized-gradient XC functionals (LDA; GGA) suffer from what is known as selfinteraction error (SIE) 39 . The approximation to the exact exchange in these functionals prevents the exact cancellation between the self-Coulomb and self-exchange for all one-electron densities. In order to minimize the total energy, single electron (Kohn-Sham) orbitals tend to over-delocalize their associated electron density. Therefore, for a charged system, the SIE spreads the electron (or hole) artificially over the fragments, yielding too low energies for the delocalized states. The terms "delocalization error" and "self-interaction error" are often used interchangeably 40 , the former typically signifying the physical aspect of the error 41 .
In this section, we explore the SIE in different charged water clusters, (H 2 O) + n , (n = 2 − 5). Previous investigations of SIE in DFT have focused on how the chargedelocalization (and hence, SIE) is affected by the size of the system 42,43 . However, most studies have considered non-interacting molecular systems where the size is tuned by simply repeating the non-interacting units separated by some finite distance. In such cases, the delocalization error worsens with increase in system-size, as ionization would result in the removal of a fraction of electron charge from all the molecules in the system.
In our study of ionized (H 2 O) + n clusters, we also consider the effect of electrostatic interactions among the H 2 O units on the size-dependent charge delocalization or SIE. Here n controls the system size, and Fig 1 illustrates the SIE for two prominent spatial configurations of (H 2 O) + n -hemibonded (HB) and proton-transferred (PT) for each (n = 2 − 5). The hole (i.e. highest occupied molecular orbital (HOMO) of the un-ionized system) has a very different electronic configuration in these two structures. In the hemibonded structure the hole arises from the antibonding combination of the 1b 1 orbitals (i.e. oxygen p z orbitals) of each of the two water molecules. By contrast, in the proton-transferred structure the hole is largely localized on the 1b 1 orbital of the H-donating molecule. This gives rise to varying amounts of charge delocalization and SIE for the same system. We compare different methods used for predicting the energies of the hemibonded and protontransferred (H 2 O) + n clusters. These include the semilocal PBE functional, hybrid functionals which combine different amounts of exact Hartree-Fock (HF) and DFT For the water dimer cation, (H 2 O) + 2 , Sodupe et al. 44 compared different DFT and post-Hartree-Fock methods to determine its ground state structure. Their observation is that GGA functionals overestimate the energies of the delocalized hole, thus incorrectly favoring the hemibonded configuration as the preferred ground state of the dimer ion. However a hybrid functional combining equal fractions of GGA and exact (Hartree-Fock) exchange -BHLYP -improves the hole-localization, and correctly predicts the proton-transferred geometry as the groundstate, in accordance with more accurate wavefunctionbased approaches. In the first row of Table I, we compare the energies obtained at different levels of theory for the hemibonded and proton-transferred dimer ions; our results are in agreement with previous studies 15, [44][45][46] , but these studies have solely focused on the water dimer cation.
The ground-state DFT calculations for all the systems are performed with the 6-311++G** basis set. The binding energy per water monomer (E b ) is evaluated as for all n = (2 − 5) ionized water clusters. In Eq. (1), E (H2O) + n is the total energy of an ionized (+1) cluster with n H 2 O molecules, E H2O is the total energy of a neutral H 2 O monomer and E H2O + is the total energy of a charged (+1) H 2 O monomer.
A crucial feature of PBE when applied to ionized water clusters is that it favors the H 2 O· · · OH 2 bonding interaction, yielding rather low energies for structures containing H 2 O monomers bonded via the O's. We refer to these as hemibonded-type geometries. The binding energies are shown for the simplest case of a water dimer ion (H 2 O) + 2 in panel (a) of Fig. 1. For the hemibonded structure, the XC functionals tested in this work provide very different binding energies ranging from -1.3 eV for PBE to -0.92 eV for BHLYP. PBE distributes the hole density evenly over the two water units over-stabilizing the hemibond configuration. As shown, the energies become larger (less negative) upon increasing the fraction of exact exchange in the XC functional. The exact exchange in the hybrid functionals compensates for the artificial hole delocalization introduced by the approximated exchange, bringing the energies closer to the more accurate MP2 and CCSD values. While the binding energies computed by CCSD and PBE differ by 0.48 eV (see Fig.  1(a)), BHLYP performs much better with a binding energy difference of 0.09 eV compared to CCSD.
In contrast, for the proton-transferred [H 3 O + − · OH] structure, all the density functionals provide very similar binding energies, which are fairly consistent with the energies from higher accuracy theories. The energy computed by PBE is 0.09 eV lower than that given by CCSD. Nonetheless, PBE wrongly selects the hemibonded (HB) structure as the ground-state representation of the ionized water dimer due to the underlying SIE. In Table I, we compare the difference in binding energies of the HB and PT structures to this effect for different methods. A negative energy difference for any given method implies that the HB structure is the preferred ground-state. It should be noted that PBE yields negative energy differences for all n (as expected), thus favoring a hemibond over a proton-transfer. However, (E b,HB − E b,PT ) PBE approaches zero from the left as n increases, implying a reduction in the SIE. In other words, for larger ionized water clusters PBE provides better estimates of binding energies of HB and PT structures, but it still misidentifies the HB geometry as its ground-state configuration. In the following sections, we will show that this is the main reason why there is a large disagreement between PBE and hybrid XC functionals in the description of the dynamical evolution of the ionized dimer. This error decreases with increasing chain length (n), as the occurrence of the hemibonded structure is largely suppressed as soon as extended H-bond chains form. Both CCSD and MP2 predict the proton-transferred geometries as the ground-state of the respective (H 2 O) + n . We remark that the XC functional that closely mimics the trend of wavefunction-based methods is BHLYP.
The binding energy differences between PBE and CCSD methods (∆E CCSD→PBE ) are highlighted for (a) Dimer, (b) Trimer, (c) Tetramer, and (d) Pentamer ionized water clusters in Fig. 1 for the two structural configurations -hemibonded and proton-transferred. The respective molecular structures for all the ionized water clusters 47 are also shown in the insets in Fig. 1. It is observed that ∆E CCSD→PBE is greater for the hemibonded geometries as compared to the proton-transferred ones for a given system size (n). ∆E CCSD→PBE has a smaller spread for the proton-transfer structures. More importantly, ∆E CCSD→PBE , which is indicative of the size-dependent SIE, decreases significantly with an increasing (H 2 O) + n system size. Specifically, it varies from 0.48 eV for n = 2 to 0.21 eV for n = 5 in the hemibonded structures, and from 0.09 eV for n = 2 to 0.04 eV for n = 5 in the proton-transferred structures. This apparent size-dependence of ∆E CCSD→PBE (or SIE), though similar in behavior, has different origins for the two fundamentally different molecular arrangements of the ionized water structures. In hemibonded structures, an increase in n results in a systematic smearingout of the hole over more H 2 O units (with PBE), thereby lowering the "localized" hole density on individual fragments and reducing their contributions to the overall selfinteraction present in the system. On the other hand, increasing n in proton-transferred structures counteracts the "delocalizing bias" 40 of PBE by effectively localizing the added hole in the ionized system. The mostly linear and open-network arrangement of H 2 O molecules in PT clusters enables the hole to be selectively localized over specific H 2 O units. At the microscopic level, this is driven by the cooperative behavior of H-bonds in openchain water geometries. The unidirectional H-bonds in these finite n-chains strengthen each other such that the HOMO of the neutral system -(H 2 O) n -is localized on the oxygen of the H 2 O which exclusively donates an Hbond (and does not accept one). With increasing n the H-bond cooperativity 48,49 becomes relevant. This manifests in a reduction of the SIE. This crucial feature guides much of our work on the excited-state dynamics of ionized water clusters that is described below.
In the following sections, we discuss both adiabatic dynamics (using Born-Oppenheimer approximation) and non-adiabatic dynamics (using mean-field Ehrenfest approach) of the ionized water chains, that is, (H 2 O) + n (for n = 2 − 5). Our aim is to capture the ultrafast processes that drive proton-transfer phenomena in ionized water clusters comprising of mostly linear H-bonds.

A. Adiabatic Born-Oppenheimer molecular dynamics
We choose small clusters of (H 2 O) n (n = 2 − 5) as our model systems. For the water dimer (H 2 O) 2 , we adopt the optimized (un-ionized) H-bonded geometry as shown in Fig. 2a, which corresponds to the energy min-  48 . These structures do not correspond to a minima in the potential energy surface of (H 2 O) n (n > 2) water clusters, which have a closed H-bond network, optimizing the formation of four H-bonds per water molecule 5 . This emphasizes the importance of H-bond cooperativity effects in studying the non-adiabatic dynamics of ionized (H 2 O) + n clusters. More importantly, this choice ensures that the photo-generated hole is always localized on the first molecule of the chain, which always forms a single donor H-bond. We first perform a Born-Oppenheimer molecular dynamics (BOMD) simulation for the neutral (H 2 O) n (n = 2 − 5) system, wherein the nuclei evolve on a single potential energy surface and the electronic structure is solved self-consistently using the DFT module of NWChem package 51 . We use a 6-31++G** basis set and the BHLYP hybrid functional. The adiabatic ground-state dynamics for each of the neutral (H 2 O) n systems is performed at a temperature of 200K using a Langevin thermostat for 100 ps, with a time step of 0.5 fs. Uncorrelated snapshots are chosen from this single adiabatic trajectory (simulated for each of the n = 2 − 5 systems) at random time intervals, to be further used as initial geometries for the excited-state non-adiabatic simulations upon HOMO-ionization at t = 0. All the structures derived from the adiabatic simulations maintain mostly linear one-dimensional H-bonds connecting each water to its neighbor(s).

B. Non-adiabatic Molecular Dynamics
The geometries extracted from the BOMD trajectory are ionized by removing an electron from the HOMO of the system, i.e. setting the net charge of the sys-tem to +1. We then perform rt-TDDFT-based Ehrenfest dynamics to simulate the excited-state dynamics, in which the electron dynamics is treated quantum mechanically, and the nuclei are classically propagated on a single mean-field surface given by an average over several electronic states. The TDDFT-based Ehrenfest dynamics of electrons and nuclei is formally given by: where R J denotes the nuclear coordinates and φ i are the time-dependent Kohn-Sham orbitals such that the electronic density is given by: ρ(r, t) = i |φ i (r, t)| 2 . V nn and V en are the nuclear and electron-nuclei interaction terms respectively, where v ext (R, r, t) is the external potential due to moving nuclei, v H [ρ](r, t) = d 3 r ρ(r ,t) |r−r | gives the Hartree potential and v xc [ρ](r, t) is the XC potential (derived from the approximate E XC [ρ]).
The non-adiabatic simulations performed for multiple HOMO-ionized (+1) geometries give rise to an ensemble of Ehrenfest trajectories (ETs). We sample these ETs to obtain the relevant time-scales and mechanisms of the proton-transfer reaction occurring in (H 2 O) + n . The rt-TDDFT/Ehrenfest dynamics simulations are performed using the real-space grid code, Octopus-8.2 52 . The optimal values of grid spacing and size of the simulation box are obtained from energy convergence tests on the systems. We use a default spacing of 0.23Å and a spherical simulation box with a radius of 8Å for the dimer, 10 A for the trimer, and 15Å for all other water chains. The coupled electron-ion dynamics uses an enforced timereversal symmetry (ETRS) propagator 53 with a time step of 1.3 attoseconds (as). This was determined to be the maximum time step that conserves the energy within a suitable range (0.7 eV over 20 fs), and the nuclear dynamics is similar to that obtained with smaller time steps (see supplemental Fig. S2a, S2b). All the ETs are propagated for up to 200 fs in order to simulate explicit dissociation or molecular rearrangement in the excited system.

IV. RESULTS AND DISCUSSION
A. rt-TDDFT diabatic trajectories An extensive study focusing on the dynamics of ionized water monomer and dimer was carried by Chalabala et al. 54 . They compared the results from two non-adiabatic approaches, namely surface hopping and Ehrenfest dynamics, and highlighted the importance of using hybrid functionals in the latter method to obtain accurate simulation outcomes. In this work, we consistently employ the Ehrenfest dynamics approach, but focus on a non-hybrid functional (PBE). We also confirm the results obtained by Chalabala et al. 54 for the smallest system, that is, a dimer cation The HOMO of the neutral (H 2 O) 2 system is localized on the H-bond donor molecule in the optimized H-bonded configuration as shown in Fig. 2a. Ionizing such a configuration results in an unpaired electron, or equivalently, a hole in this "occupied" orbital. The timeevolution of a photoionized water dimer performed using rt-TDDFT/PBE produces one of the two reaction channels, shown in Fig. 3b. We choose to simplify the outcome channels into a binary -proton transfer and no proton transfer -classification. The proton transfer channel indicates that the proton from the ionized molecule is transferred to its immediate H-bonded neighbor. This channel can be further divided into separate channels with bound and unbound molecules, see black and red arrows in Fig. 3a. However, we do not make this distinction in our trajectory-based statistics. In order to sample the population of these channels for the simulation period, we average over 53 individual Ehrenfest trajectories, see Fig. 3b. The initial short-time Ehrenfest dynamics (t < 15 fs) is characterized by a rapid proton transfer seen in a majority (fraction of 0.85) of the simulated trajectories. This is almost always followed by a proton bounce-back to the original donor oxygen, signaled by a decrease in the O d · · · H + bond length. This result confirms that in non-adiabatic Ehrenfest simulations, it is imperative to propagate the nuclei on the excited-state surface for longer times to conclusively determine the underlying mechanism and outcome of the simulated dynamics 9 . The black curve in Fig. 3b shows the successful protontransfer trajectory statistics. A proton is considered as transferred when the H + · · · O a bond distance d Oa···H < 1.5Å, (a = acceptor). The second, no proton-transfer channel, results from a molecular rearrangement of the two water units with the net charge (+1) being shared by them. In these trajectories, the two water monomers transition from a short-living proton-transferred geometry to a hemibonded-type configuration. We observe a prevalence of this channel in the trajectory statistics using the semilocal PBE functional. At the end of the simulation period (t = 200 fs), a proton-transferred structure is formed in only 13% of the initiated trajectories. On the other hand, using the BHLYP (hybrid) functional in Ehrenfest dynamics, Chalabala et al. 54 found 95% of their trajectories resulting in proton-transfer events. We associate this difference to the fact that PBE spuriously determines the hemibonded structure to be lower in energy ∆E PT→HB = −0.18 eV (see Table I).
2. Trimer (H2O) + 3 Fig. 4b shows the proton-transfer population statistics for chainlike structures of ionized water trimers. These structures contain an extra H-bond (compared to the dimer), and the HOMO spin density in a corresponding neutral trimer chain is localized on the single Hdonating oxygen, with some weight on the H-bond acceptor, see Fig. 4a. The trajectories supporting a protontransfer suggest the formation of a long-lasting intermediate [H 3 O + · · · · OH] bonded-pair before the fragments dissociate to give mobile hydroxyl radical and a reactive hydronium cation. For the 20 simulated trajectories, the first proton hop is observed to be relatively fast (within 10 fs in a fraction of 0.90 of the trajectories), and 40% of the trajectories result in a proton-transfer at the end of simulation period (t = 200 fs). This increase in proton-transfer population (in comparison to the dimer ion case) correlates with the energies in Table I. The PBE energies listed in Table I indicate that [E HB − E PT ] PBE < 0. However, the magnitude of the energy difference, |E HB − E PT | PBE , reduces in going from n = 2 to 3, making the hemibonded geometry in trimer not as highly favorable as it is for dimer. Moreover, the presence of another H-bonded water (right-most) in this chain makes it less likely for the central water to twist out of the H-bonded geometry into a hemibonded-type one.
In some trajectories (20%), a secondary protontransfer from H 3 O + to the nearby H 2 O occurs. The simulations reveal two key factors that influence such a process: the interatomic O· · · O distances, and the time it takes for the initially nonparticipating H 2 O to break away from the two interacting water-monomers compared to the first (primary) proton-transfer time. If the first proton-transfer takes too long to complete due to excessive proton-rattling between the donor and acceptor oxygens, the nonparticipating H 2 O is driven away by the electronic forces acting on the system, and the likelihood of a secondary proton transfer in the trimer chain is severely reduced. Since a secondary proton transfer event eliminates the probability of return of the proton to its original donor ( · OH), the bump seen in the secondary proton transfer curve (Fig. 4b) at shorter times (t ∼ 30 fs) contributes toward an increase in the overall "first" P(PT).

Tetramer (H2O) + 4
As in the previous cases, in a neutral H-bonded water tetramer linear chain, the HOMO is primarily localized on the first molecule of the chain, which forms only one donor H-bond but no acceptor H-bonds (see Fig. 5a). Unlike the n = 2, 3 cases, all the tetramer (n = 4) trajectories exhibit a first proton-transfer step followed by fewer bounce-backs (only in 0.30 of the 34 simulated ETs) to the original O-donor, as shown in Fig. 5b. The fluctuations in proton-transfer population last up to 75-100 fs, after which a stable proton-transfer reaction statistics is obtained. A successful proton-transfer is observed for 70% of the trajectories. This increase correlates with the results in Table I, which shows how the energy difference between hemibonded and proton-transferred structures decreases with increasing chain lengths (n) in PBE. That is, the self interaction error, while still present, becomes smaller.
After the primary proton transfer reaction, the H 3 O + cation may lose its proton to the neighboring H 2 O (secondary proton transfer). This phenomenon can continue to propagate down the chain, depending on two factors: the length of H-bonded chain (n), and the time it takes for the linear order of H-bonds to disappear to form a gas phase, which we refer to as the H-bond chain lifetime, 'τ '. For a particular system (and trajectory), τ serves as a cut-off after which the proton ceases to be pushed any further down the "chain" and maintains its position on the latest water unit. In the case of a water dimer ion, frequent proton-hops back and forth (rattling events) between the donor and acceptor units result in a prolonged lifetime τ . At the same time, longer water clusters allow multiple proton transfers from one unit to the next in the chain, thereby also extending the lifetimes. In general, a faster first proton transfer shortens the lifetime. Additionally, the first and secondary proton-transfer events occur earlier with increasing system size n, see Table II for a comparison of all relevant time scales. Therefore, the cooperative H-bonds in the ionized water chains induce a faster proton dynamics.

Pentamer (H2O) + 5
In the pentamer chain, a fast first proton-transfer step is recorded within an average simulation time of 6.5 fs (for 20 simulated trajectories). This is immediately followed by a secondary proton-transfer seen in most trajectories.
94% of the initiated (20) ETs predict proton-transfer reaction at the end of the simulation time (t = 200 fs) as depicted in Fig. 6b. Thus, rt-TDDFT/PBE gives a very high proton-transfer probability for an ionized pentamer chain. The proton-transfer mechanism is governed by the formation of a bonded H 3 O + ion − · OH radical contact pair at fixed O d −O a distance. Subsequently, the two fragments separate after the occurrence of another proton hop from H 3 O + to the neighboring H 2 O in the chain. An important characteristic of this coupled electron-nuclear dynamics is that the separation of H 3 O + − · OH is fundamentally driven by the downhill electrostatic potential for the proton to move along the H-bonded water chain.
With 94% of the (H 2 O) + 5 ETs exhibiting a protontransfer reaction, we expect the ionized pentamer chain (n = 5) to be the PBE saturation limit as far as enhanced proton-transfer dynamics due to cooperative H-bonding interaction in water is concerned. This implies that for these sizes of chains, PBE and hybrid functionals should yield very similar results. For n > 5, the (H 2 O) + n chains are likely to get divided into smaller sub-chains as the nuclei evolve in time, or they would prefer to exist as bifurcated (branched) structures -both of which can be treated as a composite of the four cases (n = 2 − 5) discussed in this work.
We have also considered the proton-transfer dynamics in a branched-pentamer chain. In this case, the alignment of H-bonds is no longer unidirectional, as a bifurcation allows one of the water units to form three H-bonds with its neighbors (instead of utmost two H-bonds per water exclusively studied so far), see Fig. 7a. 76% of the 20 trajectories studied indicate proton-transfer over the simulation period, as shown in Fig. 7b . This falls midway between the tetramer-and pentamer-Ehrenfest statistics, which is not surprising because the geometric arrangement can be viewed as either a tetramer chain with an extra water, or a pentamer chain with a (nonlinear) displaced water. A distinctive feature of this type of chain is that all the trajectories supporting a first proton transfer also show a transient secondary protontransfer event, see Fig. 7b at t ∼ 25 fs. At the end of the simulation window ∼ 50% of the simulated trajectories result in a secondary proton-transfer , similar to the case of linear (H 2 O) + 5 . This is because H 3 O + formed after a secondary proton-transfer is stabilized by donating two H-bonds. This final structure is also found to occur in  the linear pentamer chain.

B. Role of hole localization and dynamics
Our results indicate that for PBE, the P(PT) in HOMO-ionized H-bonded water chains increases with the length of the chain. We also know 54 that hybrid functionals predict a P(PT) for the dimer similar to what PBE predicts for the pentamer chain. Our study has pointed to the SIE associated with an overestimation for hemibonded-type structures in order to explain the differences between the results for hybrid and non-hybrid functionals. However this is not enough to understand all the results. In particular, the role that the localization of the photoionized hole plays on the P(PT) also needs to be considered. This is also important to understand to what point non-adiabatic simulations are necessary to accurately describe the rate of proton transfer upon single-ionization. In their study of the ionized water dimer, Chalabala et al. 54 showed that the rate of proton transfer is different in BOMD and rt-TDDFT simulations, both for hybrid and GGA functionals. Here we evaluate this in an (H 2 O) + 4 chain, comparing PBE and PBE0. As a proof of concept, we choose n = 4 for BOMD among all lengths considered in this study.

Photo-hole localization in adiabatic molecular dynamics
The BOMD simulations were performed using the NWChem package 51 with a time-step of 0.25 fs for the propagation of nuclei. BOMD constrains the evolution of the system on a purely adiabatic potential energy surface (PES), and does not time-evolve the electronic states, which therefore, instantaneously adapt to the moving nuclei, with a parametric dependence on the nuclear coordinates. As shown in Fig. 8  In the case of adiabatic BOMD, PBE does not predict any proton-transfer ( Fig. 8(a)), while in PBE0 more than 50% of the trajectories result in a proton transfer ( Fig. 8(b)). We have shown in the previous section that in non-adiabatic coupled electron-nuclear dynamics (Fig. 5b), PBE predicts a proton-transfer rate of ∼ 70% for the same system. The reason for such a discrepancy between adiabatic PBE and PBE0 on the one hand, and between adiabatic and diabatic PBE results on the other, is attributed to the nature of the wavefunction that the system is initialized in. In non-adiabatic Ehrenfest dynamics, there is a choice between initiating the dynamics (at t = 0) with (i) a relaxed wavefunction for the ionized system, and, (ii) an optimized wavefunction for the corresponding neutral system followed by the removal of an electron from its HOMO. In the latter case, the simulation starts with unrelaxed electronic states. Fig. 9 (top) shows that the localization of the hole in the unrelaxed wavefunction is almost the same for PBE and PBE0, both in the dimer as well as the pentamer chains. The results for all chain lengths and for other functionals are presented in supplemental Fig. S7. Moreover, as the number of water molecules in the H-bond chain increases, the hole-localization on the first molecule of the chain decreases by the same amount for all functionals. This indicates that the SIE in the unrelaxed wavefunction is non-dominant. It has been hypothesized that such a definition of the initial wavefunction might generate the "correct dynamics" even with PBE 54 . Accordingly, rt-TDDFT/PBE benefits from the construction of a "correct" initial state 58 .
On the other hand, when the wavefunction is allowed to relax, the hole-spread increases with increasing chain length (see Fig. 9 (bottom), and supplemental Fig.   S8). The differences between PBE and hybrid functionals are significant. In BOMD, such a relaxed wavefunction is always obtained after the convergence of the initial (t = 0) self-consistency cycle. This explains such a large difference between the trajectories for PBE0 and PBE in Fig. 8. In PBE, the P(PT) is totally suppressed. PBE presents a much larger hole-delocalization than PBE0 at the initial step. These results indicate that P(PT) increases with increasing initial photo-hole localization in H-bonded water chains. However, our results also show how this localization decreases with increasing chain length -in fact, this result is general for all XC functionals. It is instructive to analyze how the actual dynamics of the hole in rt-TDDFT, diabatic, dynamics differs between hybrid and non-hybrid functionals. This might help us understand why the proton transfer probability increases despite the decreasing hole localization.

Dynamic evolution of the photo-hole in diabatic trajectories
We analyze the time-evolution of the hole generated upon removal of an electron from the HOMO of the system at t = 0. The hole-density at any time instant is estimated by computing the difference between electronic densities of the ground and ionized states for the same geometric configuration. We characterize the spatial holeevolution using: Here, the x -axis is chosen to lie along the H-bond chain such that ∆ρ(x) measures the variation of the hole density at H 2 O positions along the H-bond axis at any given time. The ionized structures and their corresponding TDDFT densities are selected from a single Ehrenfest trajectory. Additionally, the ground state densities are computed for these structures using static DFT, which are then subtracted from the respective TDDFT densities. We first focus on the hole dynamics in (H 2 O) + 5 since this system exhibits a high probability of proton transfer. Fig. 10 shows the results for the time-evolution of ∆ρ(x) for a selected Ehrenfest trajectory using GGA/PBE. The relative positions of the evolving molecular species along x -axis are also indicated. We have confirmed that a majority of the simulated non-adiabatic trajectories describing proton-transfer at the nuclear scale show a similar behavior. The hole-evolution begins at t = 0 when an electron is photo-excited from the HOMO of the system by explicitly changing the occupation number of the state. At the start of the simulation, the hole density is predominantly localized on the H 2 O unit that exclusively donates a H-bond as shown in the isosurface plot of HOMO of the ionized system (inset, Fig. 10 (a)   The time evolution of Q left (t), for different XC functionals is shown in Fig. 11, where the integrated hole density is computed at every hundredth step-interval (0.13 fs). Fig. 11(a) and (b) show two independent Ehrenfest trajectories for (H 2 O) + 3 , starting from different initial geometric configurations. In Fig. 11(a), PBE and PBE0 trajectories complete a first proton-transfer event, while BHLYP provides an additional secondary protontransfer, however in Fig. 11(b), PBE shows no protontransfer while both the hybrids PBE0 and BHLYP suggest up to a secondary proton-transfer. The insets capture the short-time hole dynamics relevant for the first proton transfer which occurs within 4.61-4.87 fs in the successful cases. In Fig. 11(a), the hole charges at the instant of the first proton-transfer are 0.76 (PBE), 0.78 (PBE0) and 0.79 (BHLYP). Additionally, a secondary proton-transfer occurs for BHLYP at 29.74 fs with a hole-charge of 0.86. In contrast, the values for the trajectory shown in (b) are found to be 0.72 (PBE0) and 0.73 (BHLYP) for the first proton-transfer followed by a secondary proton-transfer occurring around 31.58 fs (PBE0 hole =0.77; BHLYP hole =0.82). PBE fails to show any proton-transfer and the associated hole-charge decays to 0.66 in the early-time dynamics. A point of commonality for all the system sizes n, is that initially all the functionals produce nearly identical hole localizations. Furthermore, there is not a large and obvious difference between the Q left (t) of PBE and PBE0.
We observe that even when there is sufficient delocalization, the electronic charge lags behind the nuclear charge, and its evolution is mediated by the nuclear motion over a prolonged period of time (∼50 fs, which is long when compared to the natural time scale of electronic motion). Fig. 12 (a) and (b) illustrate the underlying effect of H-bond cooperativity present in longer chains on hole density evolution. The overlapping curves highlighted in the insets show that the cooperative Hbonds in larger n chains reduce the dependence of both hole and proton transfer on the choice of functional. After an initial proton transfer reaction, the proton often proceeds to move to the next water molecule in the chain. The nature of this reaction is different from the primary one as both the reactants and the products are inherently different from those found in the primary reaction of interest. In the case of the secondary proton hop reaction, this is initiated from an already-formed Zundel complex 59 , where the proton hops from the hydronium to the next water unit. At the nuclear level, a secondary proton-hop reaction is improbable in shorter chains. In these structures (n = 2, 3), the hemibonded configuration is still accessible (and favored by PBE) in the simulation -which suppresses the 'proton-transfer' branch and results in a dissociation of the H 2 O monomers without any proton transfer. On the other hand, we see a significant secondary proton hop dynamics in longer ionized water chains (n = 4, 5), and the probability of such a reaction grows with the number of water molecules in the chain. This explains why, for n = 5, the probability of a first proton transfer is greater than 90%.

V. CONCLUSIONS
We have studied the proton transfer mechanism upon photoionization of H-bonded water molecular chains, as a function of the length of the chain using non-adiabatic rt-TDDFT simulations in the Ehrenfest approximation. The goal of this study was to understand how the self interaction error in semilocal (GGA) density functionals influences the proton transfer probability. We have shown that for PBE, this probability increases from 13% in the dimer trajectories to 94% in the pentamer trajectories. We have also shown that while the probability is largely underestimated in small chains (dimers and trimers) as compared to what hybrid functionals predict, the error is minimal in longer molecular chains. The results indicate that for PBE, the proton transfer is disfavored in cluster sizes of 3 molecules or less, due to the overestimation of the stability of hemibonded geometries over proton transfer geometries, which is in turn due to the self interaction error. In longer chains (n = 4, 5 in (H 2 O) + n ), the increased H-bond cooperativity makes the transition to hemibonded-type geometries less probable. An increasing fraction of the simulated Ehrenfest trajectories exhibit proton-transfer, often along multiple water molecules in the chain.
There is also a clear difference in photohole delocalization in adiabatic, Born-Oppenheimer molecular dynamic simulations, when comparing PBE to hybrid functionals. However, the charge dynamics in rt-TDDFT simulations is quite comparable for all the functionals studied in this work. This is due to the relaxation of the initial wavefunction of the ionized system. The adiabatic scheme (BOMD/PBE) at comparable time scales fails to describe the evolution of the excited system toward a proton-transfer reaction. The inclusion of nonadiabatic effects is therefore essential for capturing the proton transfer dynamics.
The results presented in this paper have important    In the unbound reaction channel, · OH and H 3 O + dissociate (d O···O > 4.5Å) whereas in the bound channel they continue to interact. Hence, the "bound" and "unbound" channels take into account the separation between the · OH radical formed and the evolving H 3 O + ion.    trajectory starts from the same initial water geometry as in Fig. S5 but shows a successful proton-transfer.