Dynamical aspects of intermolecular proton transfer in liquid water and low-density amorphous ices

PHYSICAL REVIEW E 89, 052130 (2014) Dynamical aspects of intermolecular proton transfer in liquid water and low-density amorphous ices Amani Tahat and Jordi Mart´ ı* Department of Physics and Nuclear Engineering, Technical University of Catalonia-Barcelona Tech, Building B5, Northern Campus UPC. Jordi Girona, 1-3, 08034 Barcelona, Catalonia, Spain (Received 3 February 2014; revised manuscript received 18 March 2014; published 19 May 2014) The microscopic dynamics of an excess proton in water and in low-density amorphous ices has been studied by means of a series of molecular dynamics simulations. Interaction of water with the proton species was modelled using a multistate empirical valence bond Hamiltonian model. The analysis of the effects of low temperatures on proton diffusion and transfer rates has been considered for a temperature range between 100 and 298 K at the constant density of 1 g cm−3. We observed a marked slowdown of proton transfer rates at low temperatures, but some episodes are still seen at 100 K. In a similar fashion, mobility of the lone proton gets significantly reduced when temperature decreases below 273 K. The proton transfer in low-density amorphous ice is an activated process with energy barriers between 1–10 kJ/mol depending of the temperature range considered and eventually showing Arrhenius-like behavior. Spectroscopic data indicated the survival of both Zundel and Eigen structures along the whole temperature range, revealed by significant spectral frequency shifts. DOI: 10.1103/PhysRevE.89.052130 PACS number(s): 61.20.Qg,68.08.De,61.30.Hn,31.15.xv I. INTRODUCTION The phenomenon of proton transfer (PT) is a fundamental process playing a key role in a wide variety of technological, medical, and chemical processes, such as in the energy production in fuel cell membranes [1,2], in fundamental molecular mechanisms occurring in viruses, such as the human immunodeficiency-1 protease (HIV-1PR) [3], in the activation of nicotine as it enters an aqueous environment [4], or as the main component in molecular reactions in aerosols [5] which are the base of atmospheric chemistry, to mention only a few. Hydronium species can be found in most aqueous environments, even in interstellar media [6]. From the fundamental point of view, PT is the main reason for the neutral pH of water due to the process of water autoionization [7] and it plays a key role in energy conversion processes such as in photosynthesis and in cellular respiration [8]. From a general point of view, it would be redundant to mention here the importance of the role of the proton as a regulator of any reactive process that takes place in aqueous media. Despite being an ubiquitous element in the chemical physics of solutions, it was not until recently that the structure of aqueous protons was properly understood [9]. This was largely the result of a series of computer simulation experiments that were able to reveal the structural characteristics of PT. Thus, there is a general consensus to describe the aqueous proton in terms of a “default” within the three-dimensional network of hydrogen bonds of water. Hence, dynamics of the proton is regulated by the control of the hydrogen-bond network. The temperature and, more specifically, the thermodynamical state of the system plays a central role both on structure and location of the proton species and also on dynamics of the proton transfer. Whereas plenty of information about PT in liquid water and in biomolecular systems at ambient conditions is available [9–22], its characteristics at low temperatures, such as in supercooled states, for the wide variety of ice *[email protected] classes and for highand low-density amorphous ices (LDA) is still largely unknown. Given the complexity of the phase diagram of water [23,24], studies of structure and dynamics of aqueous lone protons are usually reported for restricted regions of the diagram. For instance, proton arrangements in ice I were studied by means of dielectric relaxation and infrared spectroscopy and analyzed by Von Hippel [25]. Spectroscopic methods such as Fourier-transform infrared (FTIR) allowed Devlin and coworkers [26–28] to monitor the hydrogendeuterium exchange in water an hydrogen chloride adsorbed in ice surfaces and to obtain activation energies of Bjerrum defect formation. Ohmine and coworkers [29,30] have reported results from ab initio calculations and observed that PT in cubic ice Ic is still fast, but significantly less than in liquid water. Experiments on amorphous ice films by reactive ion scattering and low-energy sputtering have revealed that PT occurs up to temperatures of the order of 100 K [31,32]. A recent study by means of reactive ion scattering and infrared spectroscopy [33] indicates the existence of efficient proton-relay channels for hydronium on amorphous ice surfaces. In a recent paper [34], we introduced some of the methods (from the point of view of pattern recognition software) that will be employed in the present study. It is well known that microscopic structure of water typical of liquid bulk is significantly enhanced in ice [35]. At the environment of an excess proton, we should expect further significant differences with local densities of pure liquids and solids due to the presence of the lone proton. Nevertheless, since the calculation of the phase diagram of the water model employed in this work is out of the scope of this paper, we have considered a series of thermodynamical states from liquid to solid at constant density, focusing on the characteristics of PT. In summary, the main aim of the present paper is to report a thorough analysis of PT in bulk (unconstrained) liquid water and LDA ices from 298 K down to 100 K by means of molecular dynamics (MD) simulations combined with empirical valence bond (EVB) calculations, paying especial attention to the dynamical aspects of PT. The technical details will be outlined in Sec. II, the results of solvation structures 1539-3755/2014/89(5)/052130(10) 052130-1 ©2014 American Physical Society AMANI TAHAT AND JORDI MART´ I PHYSICAL REVIEW E 89, 052130 (2014) of the lone proton and dynamics of PT will be described and discussed in Sec. III and we will give some conclusions in Sec. IV. II. METHOD The computer simulation experiments reported in the present study were performed using a multistate empirical valence bond approach. The implementation of this kind of methodology to study chemical reactivity has been described extensively in the literature [11,12,15,36–48]. For such reason, here we will restrict ourselves to introduce only to the main features of the method, following the ideas outlined in Ref. [49]. It is important to note that quantum effects may be described in a variety of ways, whose choice will depend strongly of the characteristics of the system under study. So, for fully quantum systems such as liquid helium or hydrogen, the choice of methods such as path integral Monte Carlo or ab initio (Car-Parrinello) molecular dynamics is in order. When the system under study consists of a quantum particle in a sea of classical molecules, other semiclassical methods such as EVB are also very appropriate. Here we should distinguish between the effects of the quantum nature of the proton on its dynamics, which will be explicitly considered in the EVB framework and those corresponding to the treatment of electronic structure, which are not explicitly addressed in the EVB methodology, being taken into account by means of pseudopotentials and effective interactions. The systems considered in the present work consisted in two parts: (1) a quantum particle (excess proton) and (2) a classical bath (liquid water or amorphous ice) formed by 125 water molecules. The EVB method assumes that the Born-Oppenheimer potential energy surface 0({R}) driving the dynamics of the nuclei with coordinates {R}can be obtained from the lowest instantaneous eigenvalue of the EVB Hamiltonian: ˆ HEVB({R})=|φihij ({R})φj|,(1) where we have considered (as it will be in most of the forthcoming formulas) the criterion of summation over repeated indexes. The EVB Hamiltonian is represented in terms of the basis set {|φi} of diabatic (localized) VB states. In the case of an excess proton in water, these diabatic states are associated to configurations with the H+located in a particular water oxygen. The ground state |ψ0of ˆ HEVB satisfies ˆ HEVB|ψ0=0({R})|ψ0,(2) and it can be expanded as a linear combination of diabatic states as |ψ0= i ci|φi;(3) leading to the final expression for the potential energy surface 0{R}=cicjhij ({R}).(4) Dynamics of the nuclei of mass Mkis governed by the following Newton’s equation of motion: Mk d2Rk dt2=−cicj∇Rkhij ({R}).(5) In the framework of EVB methods matrix, off-diagonal elements hij can be cast out in terms of nuclear coordinates, achieving an excellent agreement with results from full quantum calculations. The parametrization employed in the present work follows those proposed by Schmitt, Voth et al. [43,45,47], which are applied to different environments and provide good results for aqueous protons in all cases. Diagonal elements hii include contributions from stretching and bending intramolecular interactions within the tagged H3O+and also inside the rest of water molecules, which are modeled using a flexible TIP3P force field [50]. In addition, diagonal elements also include intermolecular interactions such as those between hydronium-solvent and solvent-solvent. Conversely, off-diagonal elements hij introduce the coupling between diabatic states iand jand have been modeled including interatomic contributions within a particular (H5O2)+Zundel water dimer spanned by states |φiand |φjplus Coulomb interactions between the dimer and the rest of solvent. A complete list of parameters is provided in Ref. [45]. Within this framework, Schmitt and Voth were able to reproduce geometries and energies of relevant protonated water clusters (H5O2)+,(H 7O3)+and (H9O4)+, obtained from high-level ab initio calculations. The construction of the EVB Hamiltonian was performed following a series of steps: First, we needed to identify the water closest to the excess proton; this water constituted the initial pivot H3O+and the first diabatic state. From this pivot, the rest of the diabatic states were chosen in a tree-like construction via a hydrogen-bond connectivity pattern. The criterion to establish a hydrogen bond was as follows: The maximum oxygen acceptor-proton donor distance was fixed up to 2.8 ˚ A; in addition, we imposed a minimum threshold value of the H-O-O angle of 30o. Finally, all molecules lying in up to the third solvation shell and showing a connecting path with the original pivot were included in the construction of the L×L EVB Hamiltonian matrix, which was properly diagonalized. We found typically Lof the order of ∼10–20 units for the connectivity pattern. In all cases, fluctuations of the total energy never surpassed 1%. At each step, proton transfer was made possible by re-assigning the pivot oxygen label to the instantaneous state exhibiting the largest c2 icoefficient; from this state, the list of participating VB states was reconstructed using the connectivity branching procedure mentioned above. Once the EVB matrix was formed, ground-state eigenvectors and Hellmann-Feynman forces were computed by means of Fk=−ψ0|∂ˆ HEVB ∂xk |ψ0=−cicj ∂ˆ Hij EVB(x) ∂xk .(6) All simulation experiments corresponded to microcanonical runs at temperatures fluctuating around T= 100,125,150,175,200,225,250,273,and 298 ±20 K. In order to differentiate temperature from density effects, we kept densities constant at the value of ambient conditions (ρw= 1gcm −3). According to the phase diagram of the rigid TIP3P model [24] such conditions would correspond to liquid water states (between 298 and 225 K) and to LDA ices (from 200 to 100 K). However, these assignments are only approximate, since (1) our model includes flexibility of the molecular bonds and (2) the pressure in our system will fluctuate, given that we 052130-2 DYNAMICAL ASPECTS OF INTERMOLECULAR PROTON . . . PHYSICAL REVIEW E 89, 052130 (2014) performed our simulations at the microcanonical ensemble (fixed number of particles, volume, and energy), allowing temperature and pressure to slightly fluctuate. Our time step was set to t =0.5 fs for all simulations. We considered equilibration periods of approximately 20 ps, followed by trajectories of hundreds of ps, used to obtain meaningful statistical properties. Long-range interactions (Coulomb terms) were handled by Ewald sum techniques [51], assuming a uniform neutralizing background. III. RESULTS A. Solvation structure of the proton The presence of an excess proton in water promotes a disruption in its local hydrogen-bond structure. The widest accepted microscopic picture of the proton aqueous solvation involves a series of structures intermediate between two limiting cations: the three-coordinated hydronium (H9O4)+ known as the Eigen cation [52] and the symmetric Zundel dimer (H5O2)+[53]. During simulations continuous interconversions between the two structures generating a hybrid (H9O4)+/(H5O2)+complex are usually observed [54,55]. The typical timescale for those interconversions is of several picoseconds and involves displacements in oxygen-oxygen (O-O) distances and changes of the hydrogen connectivity pattern between the complex at its first end (eventually) second solvation shells. Snapshots of the local water structure at several different temperatures are shown in Fig. 1. There we have only included, for the sake of clarity, those molecules having the largest weighting coefficients cifrom Eq. (4)(ofthe order of 20 molecules), which usually are the water molecules of the first and second solvation shells of the instantaneous hydronium species. From this picture we can have a first direct indication of the fact that the local environment of the proton (first water shell) is essentially the same in all cases, namely a three coordinated Eigen cation, although the general structure of the local cluster may be quite different for the three temperatures considered. We observe, for instance, that as temperature decreases the local environment of the proton becomes more localized, with a closer first solvation shell. We will analyze microscopic details of solvation by considering FIG. 1. (Color online) Snapshots of local configurations around the pivot water (oxygen in blue) at different thermodynamic states (left to right): T=150,225,and 298 K. Only water molecules having largest coefficients Ci(typically 20–30 molecules) are explicitly shown. 0 0.5 1 1.5 2 2.5 ρΟ Ο(r) (Å-3) 298 K 273 K 225 K 150 K 100 K 0 1 2 3 4 5 ρΟ∗Ο(r) (Å-3) 298 K 273 K 225 K 150 K 100 K 234567 r ( Å ) 0 0.5 1 1.5 2 ρΟ∗Η(r) (Å-3) 298 K 273 K 225 K 150 K 100 K 234567 r ( Å ) 0 0.5 1 ρΟ Η(r) (Å-3) 298 K 273 K 225 K 150 K 100 K FIG. 2. (Color online) Oxygen pivot (O∗)-oxygen solvent (O) site-site pair correlation functions at different thermodynamic states (top left). Oxygen pivot (O∗)-hydrogen solvent (H) (bottom left). Oxygen-oxygen pair correlation functions (top right), oxygenhydrogen pair correlation functions (bottom right). the local pivot-water density fields given by ρo∗α(r)=1 4πr2 i δro∗−rα i−r,(7) where rO∗is the coordinate of the instantaneous pivot oxygen and rα idenotes the coordinate of site α(O,H) in the ith solvent molecule. Spatial correlations of oxygen pivot-oxygen water sites are shown in top panel at left side of Fig. 2. Interestingly, structures of the pivot-oxygen profiles are dominated by a first solvation shell located at r=2.5˚ A, including three or four acceptor water molecules, being essentially the same at all temperatures. That similar local structures reveal that the proton is able to promote a considerable extent of solvent clustering in its close vicinity at low temperatures in ice ambients in a close fashion to what is seen at ambient conditions. Moreover, the presence of a second shell located at ∼r=4.5˚ A is also clearly perceptible in all thermodynamic states investigated. These molecules correspond to the second solvation shell of the hydronium, and some of them may act as hydrogen bond acceptors or donors with respect to the inner ones. As temperature decreases, the position of such second shell waters tends to move at lower distances, from ∼4.6˚ A at 298 K to ∼4.2˚ A at 100 K. The local cluster tends to become smaller as temperature goes down. This promotes a larger extent of water localization in a similar fashion as it happens in cubic ice [29,30], since packing and ordering is common for solid-like states, instead of the typical tendency to disorder of liquid-like states. The analysis of the oxygen pivot-hydrogen water profiles is shown at the bottom panel at left side of Fig. 2and it provides complementary information. Here, in all cases, we found main 052130-3 AMANI TAHAT AND JORDI MART´ I PHYSICAL REVIEW E 89, 052130 (2014) peaks located at r=3.15 ˚ A. At low temperatures, these peaks include exclusively the six hydrogen atoms corresponding to the water molecules belonging to the first solvation shell; as temperature increases the number of hydrogen atoms included in the first peak raises to ∼9. An analysis of the connectivity of these new hydrogen atoms reveals that, typically, two of them belong to molecules from the second shell acting as hydrogen bond donors to the inner ones. At 100 K, we observe how the structure of water hydrogens around the pivot oxygen enhances and reveals three maxima at ∼4.3,5.2,and 6˚ A, suggesting a tendency of the system to evolve towards a more ordered solid-like configuration. We finally remark that, regardless the temperature investigated, we found no evidences of pivot acceptor hydrogen bonding of the type O–H···O. As we will discuss in the following sections, these observations concerning the hydrogen bond connectivity may have relevance in determining the mechanisms that drive the transfer of the proton. At the right column of Fig. 2we displayed the oxygen-oxygen ρO−O(r) and oxygen-hydrogen ρO−H(r) radial density fields. The locations of the main peaks are in overall good agreement with the findings of Botti et al. [16] (Fig. 2) obtained from neutron diffraction experiments of HCl dissolved in water, which indicates that the local water ordering due to the presence of the lone proton is not qualitatively different of that produced by protons of strong acids like HCl, where Zundel and Eigen structures solvating a lone proton in water survive to the influence of large anions such as chlorine [18]. Interestingly, up to three-four water layers can be observed in the O-O and O-H profiles for temperatures corresponding to LDA ices, which indicates that in such systems solid-like structures tend to appear. B. Dynamics of proton transfer We will start analyzing the nature of the proton transfer dynamics in liquid water and ices by direct inspection of the time evolution of the pivot oxygen label during 50 ps time intervals, as shown in Fig. 3. Seven representative temperatures (100, 125, 150, 225, 250, 273, and 298 K) are shown, but the effect of temperature on the frequency of proton transfer episodes is directly seen, by simply counting the number of transitions in the four plots: At T=298 K (bottom panel), approximately ∼20–25 water molecules retain the pivot label during time intervals of the order of 0.5 ps or longer, roughly delivering a transfer time of 0.4–0.5 ps−1. That number is at least fivefold smaller as we move to ice-like temperatures, keeping the density constant. A few PT can be still seen at 150 K and even one of them has been captured at 100 K (given the short time interval considered here). We should point out that the predicted rate of transfer at ambient conditions is a factor ∼8 larger than the one inferred from results of NMR experiments [56–58], being that this a well-known deficiency of the semi-classical picture adopted here; moreover, the explicit incorporation of quantum fluctuations in the transferring proton yields a better agreement with the experiments, leading to rates at least twice as large as the semi-classical ones [43]. The overall jump patterns look quite similar in all thermodynamic states and can be regarded as a sequence of episodes 0 50 100 0 50 100 0 50 100 0 50 100 Label of pivot oxygen 0 50 100 0 50 100 0510 15 20 25 30 35 40 45 50 t (ps) 0 50 100 100 K 125 K 150 K 225 K 250 K 273 K 298 K FIG. 3. Time evolution of pivot-oxygen labeling in different aqueous environments, from liquid to LDA ice systems (top to bottom). in which the proton resides in one water during a few ps, interrupted by intervals in which the proton resonates rapidly between two valence bond states, establishing what is usually called a “special” bond [42]. Some few isolated spikes reveal single attempts of aborted transitions in all cases. The crude picture provided above may be substantially improved using time correlation functions. To do this, we will need that in all cases proton transfers were sufficiently frequent to collect statistics without employing special techniques designed to analyze rare events dynamics, such as transition path sampling or other importance sampling schemes [51]. Equilibrium time correlation functions for the population relaxation of different reactant species have been defined previously [47,59]. The general form of these functions is as follows: C(t)=δhi(t).δhi(0) (δhi)2(8) where the difference function δhi(t)=hi(t)−hidescribes the instantaneous fluctuation of the population of ith reactant away from its equilibrium value. The characteristic function hi(t) is unity 1 if the tagged reactant species is present in the system at time tand 0 otherwise. In the present case, population relaxations of the pivot label are the most natural functions to investigate. From previous observations, we can expect that time correlation functions of hi(t) will show at least three different temporal domains: (1) a resonant time τrsn in the subpicosecond scale, associated to the rapid exchange of the pivot label, i.e., the excess proton, along a “special” bond, represented by spikes in the history of 052130-4 DYNAMICAL ASPECTS OF INTERMOLECULAR PROTON . . . PHYSICAL REVIEW E 89, 052130 (2014) 0 0.4 0.8 1.2 1.6 2 t ( ps ) -1 -0.8 -0.6 -0.4 -0.2 0 ln C(t) 298 K 273 K 250 K 225 K 150 K 100 K FIG. 4. (Color online) Logarithm of the population relaxations for the pivot-oxygen label at different thermodynamical states. Linear fits between 1–2 ps are represented by dashed lines. the pivot labels depicted in Fig. 3; (2) a second time scale τprs characterizing the averaged lifetime of the resonance episodes; and (3) the residence time τrsd of the proton when attached to one particular oxygen pivot. Results for the population relaxation of the pivot label are shown in Fig. 4. The presence of more than one relaxation time is clear from the absence of a single clear linear regime in the time interval analyzed. By means of Onsager’s regression hypothesis[60], we can obtain the proton transfer rates kpfrom the long time slopes of C(t) (see Fig. 4): kp=lim t→∞ −dln C(t) dt .(9) The average mean residence time of the proton in a pivot water is τrsd =k−1 p. The results for the thermodynamical states considered in this work are reported in Table I. TABLE I. Dynamical parameters for the aqueous protons at different thermodynamic states: proton transfer rates τ−1 rsd , residence time τrsd, and diffusion coefficient of the lone proton Dp. T(K) τ−1 rsd (ps−1)τrsd(ps) Dp(˚ A2/ps) 100 0.022 45.5 0.17 125 0.043 23.3 0.26 150 0.090 11.1 0.33 175 0.095 10.5 0.46 200 0.098 10.2 0.57 225 0.101 9.9 0.68 250 0.197 5.1 0.81 273 0.225 4.4 0.87 298 0.540 1.9 0.94 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 T-1(K-1) -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 ln kp FIG. 5. (Color online) Logarithm of proton transfer rates kpas a function of inverse temperatures. Straight lines represent the best linear fits for different temperature ranges: overall (red line); between 298 and 225 K (green line); between 225 and 175 K (orange line); and between 175 and 100 K (blue line). The general trend is a systematic slowdown of proton dynamics when the system is cooled down to LDA ice states: proton transfer rates decrease and, equivalently, the estimated residence times τrsd increase. The comparison to other simulation works, reveals a good overall agreement with findings from Day et al. [47], who obtained a value for the proton transfer rate of 0.3 ps−1at room temperature (300 K), for an EVB model different of the one used in the present work. In the case of cubic ice, it was observed [29,30] that the ratio between PT rates of liquid and ice phases is about a factor 2. This fact was attributed to the larger extent of localization and alignment of the O −O∗ pairs, which would be at the basis of the PT mechanism in ices. In our case, when cooling down the system PT rates tend to decrease, following a monotonic behavior and showing a tendency to satisfy an Arrhenius-like dependence (see Fig. 5) with some clear deviations. Hence, assuming the following dependence of the proton transfer rate with temperature: kp∼Ae−Ek kBT,(10) where Ais a proportional factor and kBis Boltzmann’s constant, we can obtain an estimation of the PT activation energy Ek. Given the slope of the straight line shown in the Arrhenius plot of Fig. 5(red line), we get a value of Ek=3.2 kJ/mol for the overall linear fit to the whole set of values. This is in qualitative good agreement (order of magnitude) with the 10 kJ/mol obtained by Moon et al. [61,62] by means of a reactive ion scattering technique for the activation energy of PT at the surface of polycrystalline ice film, prepared at 135 K. Therefore, PT is an activated process at low temperatures and it requires to surmount an energy barrier of quite important magnitude. When the process of PT is mediated by hydroxide ions, the energy barrier has been measured at 9.6 kJ/mol [32]. It should be noted that at 052130-5 AMANI TAHAT AND JORDI MART´ I PHYSICAL REVIEW E 89, 052130 (2014) the different subintervals of temperatures (between 298 and 250 K, from 225 up to 175 K and from 175 to 100 K) different slopes can be obtained, revealing different behavior when the system is at liquid phase, which undergoes the phase transition from liquid to LDA ice and at the LDA ice phase (in all cases at temperature ranges different from experimental data, because of the present water model). A fine tuning of the slopes is also shown at each range. At the liquid phase (green line), we get an estimated activation energy of 12.8 kJ/mol, at the interval between 225 and 175 K the value is just of 0.6 kJ/mol and at the LDA ice phase, we get about 3.4 kJ/mol. We should note that Luz and Meiboom [63] obtained 10 kJ/mol for the activation energy of PT in pure water between 288 and 348 K from proton magnetic relaxation measurements, in excellent agreement with the 12.8 kJ/mol reported here in the range of 250–298 K. In conclusion, the largest energy barrier to surmount for the proton corresponds to the liquid phase and the lowest to the range 225 to 175 K. This suggests that the slowdown of PT at low temperatures is mainly due to the lack of thermal energy, since the energy barriers are lower than at high temperatures, probably due to the higher degree of localization of the proton, as we can note from Fig. 2, where LDA ice has a local structure richer than that of water at ambient conditions, with 2–3 water layers clearly distinguishable. Some years ago [20], the existence of a direct relationship between the likeliness of PT and the distance between oxygens of solvating water molecules was suggested. At small distances, if the lone proton is equally shared by two waters (Zundel dimer) it corresponds to a minimum in the external potential of the proton along the O-O axis. However, if the proton is closer to one of the waters, the potential shows a maximum. Nevertheless, it was indicated [20] that the correct picture should be given by a twodimensional potential depending on at least two variables: RO−Oand the proton displacement coordinate. Indeed the consideration of multidimensional reaction coordinates for water autoionization was already suggested by Geissler et al. [7]. In our case, the connection between O-O distances and activation energies for PT would be in agreement with the simple description pointed out above, since states with lowest activation energies (between 225 and 175 K) are those showing O-O distances smaller than the corresponding ones at liquid states (see Fig. 2, left side). C. Proton diffusion The diffusion coefficient of aqueous protons at ambient conditions is known to be approximately fourfold that observed for neat water. So, the experimental value is of 0.93 ˚ A2/ps [64] for a proton diffusing in water at 298.15 K and at the density of 1 gcm−3, whereas the value of the diffusion coefficient of bulk liquid water is of 0.23 ˚ A2/ps [65]. The main reason of such an enhancement of the diffusion is the Grotthus translocation mechanism [66], in addition to the usual hydrodynamic Stokes mass diffusion. The calculation of diffusion coefficients of aqueous protons Dpin ice ambients and constant densities reveals interesting changes in this scenario, as it happens at other thermodynamic conditions, such as at interfaces or at high temperature and in supercritical states [49]. We obtained Dpfrom long time slopes of mean square displacements of the proton coordinate rp, in the usual way: Dp=1 6lim t→∞ d dt |rp(t)−rp(0)|2,(11) where the proton coordinate was defined as a weighted sum of the coordinates of the Lpivot molecules, ri pvt: rp= L  i c2 iri pvt.(12) Results for the diffusion coefficients are shown in Fig. 6and numerically reported in Table I. Two important features should be discussed: (1) The general trend of the proton mobility is a neat reduction from 0.94 ˚ A2/ps at room temperature to lower values of about one order of magnitude smaller at 100 K. Dpat 298 K is in overall good agreement with previous works [49] and in excellent agreement with the experimental value reported above (this may be fortuitous). (2) The simulation results also predict the qualitative changes observed experimentally in hexagonal ice networks [25], where the mobility of protons is due to collective effects, affecting the diffusive regime of the proton and turning it from a highly mobile solute at ambient conditions into a much slower particle at lower temperatures, with transport properties very much akin to those of a prototypical cation of small size, such as Li+, whose diffusion coefficient is of the order of 0.1 ˚ A2/ps at 298 K [49]. Again proton diffusion can be regarded as an activated process (see linear fits in Fig. 6), but in this case the overall activation energy (red line, ED∼1.5 kJ/mol) is very far from experimental findings reporting activation energies of water self-diffusion, between 14 and 70 kJ/mol for water at ice surfaces and in bulk, respectively [62]. The fits at the two phases reveal a slightly higher activation energy at the liquid phase (blue line, about 2.9 kJ/mol vs. 1.5 kJ/mol for LDA ices, green line). 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 T-1(K-1) -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 ln Dp FIG. 6. (Color online) Logarithm of diffusion coefficient of proton Dpas a function of inverse temperatures. Straight lines represents the best linear fits, for liquid states (blue line), LDA ices (green line), and overall (red line). 052130-6 DYNAMICAL ASPECTS OF INTERMOLECULAR PROTON . . . PHYSICAL REVIEW E 89, 052130 (2014) The qualitative change in the results for the diffusion constants suggests that the role of the Grotthus mechanism as a key factor to determine the abnormal, high conductivity of the proton should become less important at low temperatures. A crude estimate of this contribution is normally obtained using simple random walk arguments: results from the proton transfer rates show that at ambient conditions the proton jumps a characteristic O-O distance, say dOO ∼3.5˚ A, every τ∼2 ps. During this time interval, the center of mass of a water molecule travels typically l∼(6Dwτ)1/2∼2˚ A; so lis comparable to dOO and the proton transfer increases its mobility in a sizable fashion. A similar calculation performed at T=250 K gives an average time for PT of approximately τ∼5 ps and water diffusion is significantly slower (about 0.02 ˚ A2/ps [67]) yielding l∼0.8˚ A and showing that spatial displacements during the proton translocation are negligible compared to that of the center of mass of the water molecules operated by ordinary diffusion. Even though our model EVB Hamiltonian is likely to predict a lower rate of proton transfer and a subestimation of the actual diffusion, we do believe that the differences are sufficiently large to guarantee that the qualitative picture captured by our simulation experiments remains physically meaningful. D. Proton spectroscopy From the experimental point of view, the nature of microscopic vibrations in water is usually investigated by means of Raman or infrared spectroscopy. Infrared spectra are measured through the absorption coefficient, α(ω) or, equivalently, by means of the imaginary part of the dielectric constant, ε(ω) [68]. Such quantum properties can be computed in the EVB framework with the aid of an absorption lineshape function I(ω), i.e., the Fourier transform of the time derivative of the dipole moment ˙μ(t)[42,69]. In the present study, we have used another observable, the velocity autocorrelation function of the lone proton Cp(t)=vp(0) ·vp(t),(13) where the proton velocity vp(t) can be obtained directly from the time derivative of its position rp: vp(t)=drp(t) dt .(14) From Eq. (13), by means of the usual Fourier transform, we can obtain a vibrational density of states Sp(ω)[45]: Sp(ω)=∞ 0 dt Cp(t)eiωt.(15) With this assignment, we have computed Sp(ω) for several thermodynamic states considered along the present work. The length of Cp(t) has been of 0.5 ps, long enough to capture all relevant proton vibrations, but also much shorter than the proton residence time (always larger than 2 ps, see Table I). As a matter of fact, we will be able to obtain relevant modes of vibration of the hydronium H3O+complex. The results are shown in Fig. 7, together with the corresponding Sp(ω) obtained from supplementary simulations of an isolated Zundel dimer (H5O2)+and of an Eigen complex (H9O4)+, i.e., in the gas phase at 298 K, in order to help explaining 0500 1000 1500 2000 2500 3000 3500 4000 4500 ω (cm-1) 0 0.5 1 Sp(ω) 298 K 250 K 200 K 150 K 100 K 0 0.5 1 Sp(ω) Zundel Eigen (A) (B) (C) (A) (B) (C) (D) (E) (E) (F) FIG. 7. (Color online) Vibrational densities of states of the lone proton Sp(ω) (in arbitrary units): (top) proton in isolated Zundel and Eigen complexes in gas phase at 298 K; (bottom) proton in water from liquid (298 K) to LDA ice states (100 K). the physical origin of the bands observed in Sp(ω). We have chosen to show the full frequency range, although the spectral domain where proton vibrations are located is between 1500 and 3600 wave numbers [69]. A common feature in the spectrum of the excess proton is found at all temperatures (bottom of Fig. 7): a series of maxima between 500 and 4500 cm−1is observed, structured into two groups of vibrations, the first one between 500 and 1900 wave numbers and a second one between 1900 and 4500 wave numbers. Since these are frequencies of proton vibrations, they will describe both regular molecular water motions and vibrations due to the particular characteristics of the lone proton/hydronium complex. It has to be pointed out that in the present case the OH stretch of the proton will be coupled to other vibrational degrees of freedom, i.e., the microscopic motion associated with such a wide vibrational band should be regarded as a combination of collective vibrations involving the proton and water molecules nearby. Leaving apart the band with maxima around 700 cm−1, typical of librational modes in water [70,71], we observe maxima centered around 1500, 2300, 3000, and 3900 cm−1. The locations of the maxima associated to proton vibrations are in good qualitative agreement with experimental data available. So, FTIR measures of HCl and NaCl aqueous solutions at different concentrations at room temperature [69], where relevant maxima associated with hydrated protons were found around 1200, 1800, and 2900 cm−1. Headrick et al. [72] reported proton vibrations at 3160 cm−1for a Zundel dimer from photoevaporation of argon in photofragmentation mass spectroscopy [72], which is also in good qualitative agreement with the features reported in the present work. Finally, Kobayashi et al. [30] found a value for the stretching of the proton in cubic ice around 2600 cm−1.This value, about 13% larger than our findings, may be attributed 052130-7 AMANI TAHAT AND JORDI MART´ I PHYSICAL REVIEW E 89, 052130 (2014) to the fact that the proton structure in cubic ice shows an important extent of directionality towards its surrounding oxygens (O −H+−O∗) that would favor fast vibrational motions, instead of the less ordered LDA structures considered here. In order to enlighten further the meaning of the spectral densities of states reported, we can establish a relationship with the data displayed at the top of Fig. 7. There we can observe that in gas phase ambients the vibrational bands associated with the Zundel dimer and the Eigen complex are centered around different frequencies. In particular, if we focus on the relevant range for proton vibrations (between 1500 and 3500 wave numbers in this case) proton vibrational modes at the Zundel dimer occur at different frequencies (about 1880 and 3500 cm−1) than those inside the Eigen complex (around 1400, 2750, and 3650 cm−1). First of all, these values should be compared with experimental data. So, on the one hand, Schwartz [73] reported the finding of a frequency maximum about 2660 cm−1for a H9O+ 4cluster (Eigen complex) from infrared absorption spectra of several water clusters in the gas phase. Such frequency has been attributed [74] to an H-bonded H3O+stretch. On the other hand, a maximum at about 1740 cm−1was reported from experimental measurements of the gas-phase infrared spectrum of the protonated water dimer [75](H 5O+ 2). In summary, there is an overall reasonable agreement of our findings with experimental data. The reported results from computed vibrational density of states by Schmitt and Voth [45] for a different potential model were of 1550 and 2860 cm−1for the two complexes, what indicates again a good agreement with our results. The relationship of the proton vibrations for the excess proton in bulk at different temperatures (bottom of Fig. 7) with the corresponding findings for the proton located at Zundel and Eigen complexes in vacuo (top of Fig. 7) can be explained as follows: (1) The signatures of the Eigen-like bands are found, for our model, about 1400 (A), 2750 (B), and 3650 cm−1(C) (top of Fig. 7). Bands located at corresponding maxima in bulk (bottom of Fig. 7) are observed in all thermodynamical states considered, although the maxima are centered at the values of 1400 (A), 2370 (B), and 4050 cm−1(C) at 298 K. The first of them remains essentially unchanged with temperature, whereas the band at 2370 wave numbers is red-shifted by 70 cm−1at 100 K. Finally, the highest frequency feature is also red-shifted by 150 cm−1at 100 K. The values around 4000 cm−1(more than the corresponding 3650 cm−1obtained in vacuo) should be attributed to shortcomings of the potential model, since hydrogen vibrations in water are always restricted to smaller values [76]. (2) Concerning Zundel-like bands, the lowest frequency one located at 1880 cm−1(D) is absent in the liquid and LDA ice spectra, whereas the band at higher frequency around 3500 cm−1(E) appears to be of low intensity at ambient conditions and located at 3400 cm−1. At lowest temperatures, this mode would red-shift by some 100 cm−1. (3) A band around 2960 cm−1(F) at 298 K and absent for LDA ice states should be, at the light of our results, unassigned with respect of the Zundel or Eigen structures. Together with this description, we can compare our results with those of Vuilleumier and Borgis [42] for a flexible SPC/E model, where the stretching modes of the hydronium complex were found at 2000 and 2650 cm−1and with those obtained by Voth and coworkers [69], who assigned the modes around 1680–1880 and 3250–3400 cm−1to pure Zundel-like vibrations, the modes around 1580–1640 and 2700–2950 to pure Eigen-like vibrations and the bands centered between 3400–3600 and 3650–3720 cm−1to a linear combination of Zundel and Eigen modes. From their interpretations, we can assign the mode (F) to pure Eigen-like vibrations which disappear as the system is cooled down. Finally, the mode (C) should be related with combinations of Zundel-like and Eigen-like vibrations. As a general trend, we observe a tendency of LDA ices to the slowdown of vibrational modes associated to the proton, given the red-shifts obtained for bands (B), (C), and (E) as well as the absence of mode (F) at low temperatures. These fact are probably due to the enhancement of the water structure for low temperature liquid and LDA ices, together with the important reduction of proton delocalization, as observed by the reduction of PT rates and proton diffusion. IV. CONCLUDING REMARKS In the present work a thorough analysis of the structure and dynamics of an excess proton in liquid water and LDA ices has been reported. We employed MD simulations together with a multidimensional empirical valence bond procedure, in order to construct a suitable Hamiltonian for the semi-classical system, formed by a quantum particle (the lone proton) embedded in a sea of classical TIP3P waters. Our findings have revealed the enhancement of the local structure of the proton in LDA ices. At the lowest temperature considered in the present study (100 K), the environment of the proton typical of ambient conditions, consisting of a mixture of Zundel and Eigen-like structures has evolved to a network of water clusters mainly formed by Zundel and Eigen-like complexes whit enhanced directionality, as it can be inferred from our spectroscopical data. However, such extent of alignment between the proton and its surrounding oxygens would be less would be less important that in the case of cubic ice Ic [29,30]. The proton in ice-like environments remains trapped to an hydronium complex for long time intervals, given an averaged transfer time of about 50 ps, whereas at 298 K, the mean time for a proton transfer is of the order of 2 ps. However, as it has been indicated in some recent experimental findings [31,32] PTs still occur. The activation energy for PT has been estimated to be 3.2 kJ/mol, in reasonable overall agreement with experiments [61,62], which reported a value of 10 kJ/mol for proton transfer in surface ices. Diffusion of the proton tends to decrease when the system is cooled down, changing from 0.94 ˚ A2/ps at 298 K up to a factor 6 smaller at 100 K. Being diffusion also an activated process, we found that activation energy barriers for diffusion are lower than those of PT roughly by a factor 2. The analysis reported in the present paper is a first step into the study of PT in two-dimensional systems, such as water drops between hydrophobic graphene plates, currently under progress in our laboratory, which may be closer to the experimental setups employed to obtain the data available that we used for the characterization and validation of our simulations. 052130-8 DYNAMICAL ASPECTS OF INTERMOLECULAR PROTON . . . PHYSICAL REVIEW E 89, 052130 (2014) ACKNOWLEDGMENTS We thank financial support from the Direcci´ o General de Recerca de la Generalitat de Catalunya (Grant No. 2009SGR1003) and the Spanish Ministry of Economy and Competitiveness (Grant No. FIS2012-39443C02-01). [1] C. Wang, M. Waje, X. Wang, J. M. Tang, R. C. Haddon, and Y. Yan, NanoLett. 4,345 (2004). [2] D. B. Spry, A. Goun, K. Glusac, D. E. Moilanen, and M. D. Fayer, J. Am. Chem. Soc. 129,8122 (2007). [3] J. Trylska, P. Grochowski, and J. A. McCammon, Protein Sci. 13,513 (2004). [4] M. P. Gaigeot, A. Cimas, M. Seydou, J. Y. Kim, S. Lee, and J. P. Schermann, J. Am. Chem. Soc. 132,18067 (2010). [5] R. Bianco and J. T. Hynes, Theor. Chem. Acc. 111,182 (2004). [6] J. M. Hollis, E. B. Churchwell, E. Herbst, and F. C. De Lucia, Nature 322,524 (1986). [7] P. L. Geissler, C. Dellago, and D. Chandler, Science 291,2121 (2001). [8] A. M¨ uller, H. Ratajczak, W. Junge, and E. Diemann, Electron and Proton Transfer in Chemistry and Biology (Elsevier Science, Amsterdam/New York, 1992). [9] D. Marx, M. E. Tuckerman, J. Hutter, and M. Parrinello, Nature 397,601 (1999). [10] M. E. Tuckerman, D. Marx, M. L. Klein, and M. Parrinello, Science 275,817 (1997). [11] T. J. F. Day, U. W. Schmitt, and G. A. Voth, J. Am. Chem. Soc. 122,12027 (2000). [12] S. Walbran and A. A. Kornyshev, J. Chem. Phys. 114,10039 (2001). [13] M. E. Tuckerman, D. Marx, and M. Parrinello, Nature 417,925 (2002). [14] H. J. Baker and H.-K. Nienhuys, Science 297,587 (2002). [15] A. A. Kornyshev, A. M. Kuznetsov, E. Spohr, and J. Ultrup, J. Phys. Chem. B 107,3351 (2003). [16] A. Botti, F. Bruni, S. Imberti, M. A. Ricci, and A. K. Soper, J. Chem. Phys. 121,7840 (2004). [17] D. Asthagiri, L. R. Pratt, and J. D. Kress, Proc. Natl. Acad. Sci. USA 102,6704 (2005). [18] A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, J. Chem. Phys. 125,014508 (2006). [19] G. A. Voth, Acc. Chem. Res. 39,143 (2006). [20] D. Marx, Chem. Phys. Chem. 7,1848 (2006). [21] J. M. J. Swanson, C. M. Maupin, H. Chen, M. K. Petersen, J. Xu, Y. Wu, and G. A. Voth, J. Phys. Chem. B 111,4300 (2007). [22] S. H. Lee and J. C. Rasaiah, J. Chem. Phys. 135,124505 (2011). [23] A. K. Soper, Mol. Phys. 106,2053 (2008). [24] C. Vega, J. L. F. Abascal, M. M. Conde, and J. L. Aragones, Faraday Discuss. 141,251 (2009). [25] A. von Hippel, IEEE Transact. Elec. Insul. 23,825 (1988). [26] W. B. Collier, G. Ritzhaupt, and J. P. Devlin, J. Phys. Chem. 88, 363 (1984). [27] M. Fisher and J. P. Devlin, J. Phys. Chem. 99,11584 (1995). [28] J. P. Devlin, N. Uras, J. Sadlej, and V. Buch, Nature 417,269 (2002). [29] C. Kobayashi, S. Saito, and I. Ohmine, J. Chem. Phys. 113,9090 (2000). [30] C. Kobayashi, S. Saito, and I. Ohmine, J. Chem. Phys. 115,4742 (2001). [31] C-W. Lee, P-R. Lee, Y-K. Kim, and H. Kang, J. Chem. Phys. 127,084701 (2007). [32] J. H. Kim, Y-K. Kim and H. Kang, J. Chem. Phys. 131,044705 (2009). [33] E. S. Moon, Y. Kim, S. Shin, and H. Kang, Phys. Rev. Lett. 108, 226103 (2012). [34] A. Tahat, J. Mart´ ı, A. Khwaldeh, and K. Tahat, Chin. Phys. B 23,046101 (2014). [35] D. Eisenberg and W. Kauzmann, The Structure and Properties of Water (Clarendon Press, Oxford, 1969). [36] C. A. Coulson and U. Danielsson, Ark. Fys. 8, 239 (1954). [37] R. S. Mulliken, J. Chem. Phys. 61, 20 (1964). [38] J. Aqvist and A. Warshel, Chem. Rev. 93,2523 (1993). [39] A. Warshel, Computer Modeling of Chemical Reactions in Enzymes and Solutions (Wiley, New York, 1980). [40] R. Vuilleumier and D. Borgis, J. Phys. B. 102, 4261 (1998). [41] R. Vuilleumier and D. Borgis, Chem. Phys. Lett. 284,71 (1998). [42] R. Vuilleumier and D. Borgis, J. Chem. Phys. 111,4251 (1999). [43] J. Lobaugh and G. A. Voth, J. Chem. Phys. 104,2056 (1996). [44] U. W. Schmitt and G. A. Voth, J. Phys. Chem. B 102,5547 (1998). [45] U. W. Schmitt and G. A. Voth, J. Chem. Phys. 111,9361 (1999). [46] M. Cuma, U. W. Schmitt, and G. A. Voth, Chem. Phys. Lett. 258, 187 (2000). [47] T. J. F. Day, A. V. Soudackov, M. Cuma, U. W. Schmitt, and G. A. Voth, J. Chem. Phys. 117,5839 (2002). [48] D. E. Sagnella and M. E. Tuckerman, J. Chem. Phys. 108,2073 (1998). [49] D. Laria, E. Gu` ardia, and J. Mart´ ı, J. Am. Chem. Soc. 126,2125 (2004). [50] L. X. Dang and B. M. Pettitt, J. Phys. Chem. 91,3349 (1987). [51] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic Press, San Diego, 2002). [52] M. Eigen and L. de Maeyer, Proc. R. Soc. (London) A 247,505 (1958). [53] G. Zundel and H. Metzger, Z. Phys. Chem. 244, 456 (1968). [54] M. Tuckerman, K. Laasonen, M. Sprik, and M. Parrinello, J. Chem. Phys. 103,150 (1995). [55] M. Tuckerman, K. Laasonen, M. Sprik, and M. Parrinello, J. Phys. Chem. B 99, 5794 (1995). [56] S. Meiboom, J. Chem. Phys. 34,375 (1961). [57] H. G. Herzt and R. Manner, Z. Phys. Chem. Neue Folge 135,89 (1983). [58] R. Pfeifer and H. G. Hertz, Ber. Bunsen-Ges. Phys. Chem. 94, 1349 (1990). [59] R. Semino and D. Laria, J. Chem. Phys. 136,194503 (2012). [60] D. Chandler, Introduction to Modern Statistical Mechanics (Oxford University Press, New York, 1987). [61] E-S. Moon, J. Yoon, and H. Kang, J. Chem. Phys. 133,044709 (2010). 052130-9