Density-Functional Tight-Binding Molecular Dynamics Simulations of

Dec 11, 2017 - (46) Potential energy surface exploration by means of the quantum mechanics/molecular mechanics (QM/MM) method has shown that the diffe...
1 downloads 17 Views 3MB Size
Subscriber access provided by RMIT University Library

Article

Density-Functional Tight-Binding Molecular Dynamics Simulations of Excess Proton Diffusion in Ice I, Ice I, Ice III, and Melted Ice VI Phases h

c

Aditya Wibawa Sakti, Yoshifumi Nishimura, Chien-Pin Chou, and Hiromi Nakai J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b10664 • Publication Date (Web): 11 Dec 2017 Downloaded from http://pubs.acs.org on December 13, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Density-Functional Tight-Binding Molecular Dynamics Simulations of Excess Proton Diffusion in Ice Ih, Ice Ic, Ice III, and Melted Ice VI Phases

Aditya Wibawa Sakti,a Yoshifumi Nishimura,b Chien-Pin Chou,b and Hiromi Nakaia,b,c,d,*

a

Department of Chemistry and Biochemistry, School of Advanced Science and Engineering,

Waseda University, Tokyo 169-8555, Japan b

c

Research Institute for Science and Engineering, Waseda University, Tokyo 169-8555, Japan

CREST, Japan Science and Technology Agency, Tokyo 102-0075, Japan

d

ESICB, Kyoto University, Kyoto 615-8520, Japan

1 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT The structural, dynamical, and energetic properties of the excess proton in ice were studied using density-functional tight-binding molecular dynamics simulations. The ice systems investigated herein consisted of low-density hexagonal and cubic crystalline variants (ice Ih and Ic) and high-density structures (ice III and melted ice VI). Analysis of the temperature dependence of radial distribution function and bond order parameters served to characterize the distribution and configuration of hundreds of water molecules in a unit cell. We confirmed that ice Ih and Ic possess higher hexagonal symmetries than ice III and melted ice VI. The estimated Grotthuss shuttling diffusion coefficients in ice were larger than that of liquid water, indicating slower proton diffusion process in high-density structures than in low-density systems. The energy barriers calculated based on the Arrhenius plot of diffusion coefficients were in reasonable agreement with experimental measurement for ice Ih. Furthermore, the energy barriers for high-density structures were several times larger than those of low-density systems. The simulation results were likely related to the suppression of proton transfer in disordered water configurations, in particular, ice with low hexagonal symmetry.

2 Environment ACS Paragon Plus

Page 2 of 43

Page 3 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1. INTRODUCTION An understanding of the physical properties of the excess proton in ice is valuable for further investigation of the fundamental proton transfer process in various novel materials.1–3 Introducing an excess proton to the ice structure is thought to cause an ionic defect due to violation of the Bernal–Fowler rules.4,5 The ionic defect induces continuous proton abstraction by adjacent water molecules via the Grotthuss mechanism.6,7 In low-density ordinary hexagonal ice (ice Ih), fast proton diffusivity has been experimentally observed using fluorescence and pump-probe spectroscopy techniques.4,8–10 The time-correlated single-photon counting measurements have also revealed that the proton diffusivity in bulk ice Ih is ten times larger than in bulk water at the same temperature.4,11,12 There have been several experimental studies focusing on proton transfer in ice systems other than ice Ih. For example, cubic ice (ice Ic), a low-density metastable structure13 containing merely chair-form hexamer arrangements of water, has shown similar proton diffusion behavior to ice Ih, despite lacking the boat hexamer configuration.14 Spectroscopic monitoring of the hydrogen-deuterium exchange reaction in amorphous ice and ice film has suggested a proton transfer mechanism based on the coupling of proton migration relay and water molecular reorientation.15–17 The simulations of ice structures has been performed for different objectives, including the assessment of theoretical and computational approaches for describing hydrogen bonding interactions,18–23 prediction and interpretation of structural and vibrational properties,21,24–31 elucidation of nucleation and crystal grow processes,32–45 and modeling of proton transfer.46–48 The recent atomistic simulations using the multistate empirical valence bond method have shown that the diffusivity of the excess proton in ice Ih at 190 K is 1.27 and ~2800 times faster than that

3 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

in liquid water at 298 K and on the surface, respectively.46 Potential energy surface exploration by means of the quantum mechanics/molecular mechanics (QM/MM) method has shown that the difference between the proton transfer mechanisms in ice Ih and in liquid water6,49,50 arises from destabilization of the Zundel (H5O2+) form due to a constraint on the distance between the two oxygens.47 The QM/MM study has further suggested that the existence of repulsive forces makes the Eigen (H3O+)10,50–53 form in ice less stable than the case of liquid water.47 Although these research efforts have provided valuable insights into the dynamical properties of proton transfer, comparison of proton diffusivity between low- and high-density ice structures at the molecular level remains incomplete. In the present study, we examine structural, dynamical, and energetic properties of the excess proton in low-density ice Ih and ice Ic, as well as ice III, with dense tetragonal crystalline form, and melted ice VI, as a representative of high-density liquid structure. The systems containing approximately one thousand atoms have been treated based on the density-functional tight-binding (DFTB)54–56 method. The DFTB-based simulations have been successfully applied to bulk aqueous solutions57–67 and water clusters,57,60,68–70 yet calculation of the solid states of water has been limited to benchmarking20 and comparison71 purposes. The present applications help compare the density dependence of proton diffusivity in ice in a more direct manner and further reveal the performance of DFTB in proton diffusion through an ordered/disordered hydrogen bond network. The remainder of this paper is organized as follows. In Section 2, computational methodology and simulation details are described. Section 3 presents the simulation results and associated discussions. Finally, we provide concluding remarks in Section 4.

4 Environment ACS Paragon Plus

Page 4 of 43

Page 5 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

2. COMPUTATIONAL METHODS AND DETAILS 2.1. DFTB Method The DFTB method is an approximated density-functional theory (DFT) method,54–56 which is roughly three orders of magnitude faster than DFT methods using double-zeta basis sets72 owing to the pre-computed parameters either calculated directly from or fitted to DFT reference values. The second- (DFTB2)54 and third-order (DFTB3)55 extensions of the DFTB method, which were based on the corresponding order of Taylor expansions of the DFT total energy expression with respect to the charge density fluctuations, have improved the description of dipole moments and hydrogen bonding. The total energy expression of the DFTB3 model is written as follows:

E DFTB3 =

1 1 1 rep 0 U AB + ∑ Dµν H µν + ∑ γ AB ∆QA ∆QB + ∑ ΓAB ∆QA2 ∆QB . ∑ 2 AB 2 AB 3 AB µν

(1)

The first term of eq. 1 denotes a distance-dependent pair repulsive potential between atoms A and 0 B. The second term consists of Dµν and H µν , which are the density matrix and Hamiltonian

matrix depending on only reference density, respectively. The γ AB and ΓAB terms are distancedependent functions54,55 that describe the interaction between two induced Mulliken charges ∆Q A and ∆QB . The γ AB term between hydrogen and other atomic species was refined with an

additional damping function for resolving issues of hydrogen bond interactions.55,73,74 By neglecting the fourth term of eq. 1, it becomes the so-called second order DFTB (DFTB2). See refs. 56 and 72 for more theoretical details of the DFTB models.

5 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.2. DFTB-MD Simulations of Excess Proton Diffusion in Ice In this study, four different ice structures (ice Ih, ice Ic, ice III, and melted ice VI) were first simulated classically. The number of water molecules in the three-dimensional periodic unit cell is 360, 216, 324, and 270 for ice Ih, ice Ic, ice III, and melted ice VI, respectively. The preliminary equilibration, which involves an initial energy minimization under the NVE ensemble, thermalization under the NVT ensemble, and equilibration under the NPT ensemble, was performed using the TIP4P/ice75 and TIP4P/200576 force fields for low-density variants (ice Ih and ice Ic) and high-density systems (ice III and melted ice VI), respectively. In each step, a single 1 ns trajectory was calculated using a time step of 1.0 fs. The temperature and pressure during NVT and NPT simulations were maintained using the Langevin thermostat and barostat as implemented in the NAMD program.77 The target temperatures were set to 50, 70, 120, 230, 250, and 270 K for ice Ih and ice Ic, while the latter three temperatures were taken into account for ice III and melted ice VI. Standard pressure (0.101325 MPa) was chosen for adjusting the density of ice Ih and ice Ic. On the other hand, higher pressures were used in the equilibration of highdensity systems (280, 290, and 300 MPa for ice III and 1100 MPa for melted ice VI). The details of the simulation box (unit cell size and density) are shown in Tables S1–S4 in Supporting Information. To investigate proton diffusion in bulk ice, one excess proton was added to each classically equilibrated structure. Figure 1 shows a comparison of the investigated ice structures. The system was equilibrated quantum mechanically for 10 ps under the NVT ensemble followed by a 100 ps production run under the NVE ensemble. The DFTB potential energy was calculated at the DFTB3 level using the 3OBw parameter set,63 which was optimized for reproducing experimental oxygen-hydrogen radial distribution function (RDF) of bulk water. The

6 Environment ACS Paragon Plus

Page 6 of 43

Page 7 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

DFTB3/3OBw method was shown to improve the description of neutral and protonated bulk water properties with a minor drawback in describing the relative energy of protonated water clusters due to the energy differences between isomers might be smaller than the DFTB3 model systematic error of ca. 3 kcal/mol.63 All DFTB-MD simulations were carried out using a time step of 0.5 fs with the in-house DC-DFTB-K program.78

2.3. Evaluation of Structural, Dynamical, and Energetic Properties of Ice and the Excess Proton The calculated DFTB-MD trajectories were analyzed to extract the physical properties of ice and the excess proton in it. The RDF and bond order parameters were evaluated as structural properties. This subsection begins with a brief explanation of how to calculate the latter quantity. The bond order parameter distinguishes different crystal structures in molecular simulations. The quantity was originally developed by Steinhard and coworkers,79 then the modified approach was proposed by Lechner and Dellago.80 Throughout this study, the configuration of oxygen atoms is the focus. The orientational order parameter is given by

qlm (i) =

1 N (i ) ∑Ylm (rij ), N (i) j =1

(2)

where N(i) is the number of nearest neighbors of the i-th oxygen atom, l is an input integer parameter, m is an integer that varies from –l to +l, and Ylm (rij ) denotes spherical harmonics of bond vector rij . When qlm(i) is averaged over the neighbors of i, the averaged form of orientational order parameter can be defined as 1 qlm (i ) = N (i )

N (i )

∑q

lm

( k ),

k =0

7 Environment ACS Paragon Plus

(3)

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 43

where N (i ) is the average number of nearest neighbors over the number of oxygen atoms. Eq. 3 is used to calculate the averaged local bond order parameter as follows:80  4π l ql ( i ) =  ∑ qlm (i)  2l + 1 m = − l

1/ 2

2

  . 

(4)

In the present study, ql (i ) values of l = 4 and 6 were estimated for all the investigated ice structures. The oxygen vector positions of last five MD steps of the NVT trajectory were sampled, which included at least 1080 configurations for the smallest system size. Note that an arbitrary choice of sampling points is allowed as long as the system is in equilibrium. Next, we describe the details on the evaluation of proton transfer rates and Grotthuss diffusion coefficients. Similar to the case of excess proton diffusion in bulk water,58,64 the following accumulation function was calculated for each NVE trajectory h (t ) = h (t − δt ) + δh (t ),

(5)

where t and δt are the simulation time and time step of the simulation, respectively. The value of h(t) at t = 0 is set to zero. The second term on the right-hand side of eq. 5 is defined as follows:

0 if no proton jumps  δh(t ) = 1 if proton jumps to a new oxygen atom (forward shuttling) . − 1 if proton jumps to the last oxygen atom (backward shuttling) 

(6)

The proton transfer rate was determined from the slope of eq. 5, which is converted into the 2 Grotthuss diffusion coefficient, DG, by multiplying the factor of l 6 (l is the proton transfer

length estimated as 2.5 Å).6 Finally, calculation of the energy barrier for the Grotthuss diffusion process is accounted for. The target quantity, Ea, was estimated using the Arrhenius relationship between the diffusion coefficients and temperature, namely,

8 Environment ACS Paragon Plus

Page 9 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

E ln DG = ln DGo − a . RT

(7)

In eq. 7, DGo denotes a pre-exponential factor, T is the temperature, and R is the universal gas constant.

9 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4. RESULTS AND DISCUSSION 4.1. Structural Properties of Ice Figure 2 shows the temperature dependence of the O–O pair RDF for the low-density ice systems, ice Ih and ice Ic. The O–O pair RDFs at the same temperature are close to each other for the first three peaks, while the shapes of low-density variants are in sharp discord with one another at distances greater than 5.9 Å. At temperatures less than 120 K, there are nine solvation shells within the distance of 10 Å. The most remarkable difference between the graphs for ice Ic and ice Ih concerns the peaks located at 6.4 and 8.6 Å. At these distances, the peak intensity for ice Ic is more than twice less than that for ice Ih. An increase in temperature to 230 K or higher leads to broader peak widths than those at low temperatures due to easier disruption of the hydrogen bond network. The diminished peaks may be identified from the weak shoulders. This observation suggests that the configuration of water molecules is disordered by slow migration rather than fluctuation around the original position. Figure 3 shows the temperature dependence of the O–O pair RDF for high-density ice structures, namely, ice III pre-equilibrated at 280 MPa and melted ice VI. The corresponding results for ice III pre-equilibrated at 290 and 300 MPa are presented in Figure S1 of the Supporting Information. We have performed calculations using temperatures higher than 230 K due to the lack of stability. For ice III, five solvation shells are found within the distance of 10 Å. The comparison with the results of ice Ih and ice Ic shows two distinctive points in the RDF of ice III: one is the absence of shoulder at 5.3 Å for the second solvation peak and the other is the similar peak height between the second and third solvation peaks. The simulated results suggest that the hydronium ion in bulk ice III is highly stabilized. For melted ice VI, the first solvation peak decreases with an increase in temperature, as observed in other ice structures. In the RDF of

10 Environment ACS Paragon Plus

Page 10 of 43

Page 11 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

melted ice VI at 230 K, there is a shoulder for the first solvation peak and a slightly pronounced minimum at 4.3 Å. In contrast, the shoulder for the first solvation peak disappears at 250 and 270 K and the peak patterns corresponding to second and subsequent solvation shells show close resemblance to each other. The dissociation of hydrogen bonds adjacent to the first solvation shell seems to be induced in the case of temperature higher than 250 K. The average bond order parameters ( ql (i ) in eq. 4) for l = 4 and 6 quantify the regularity of the configuration of water molecules in terms of the presence of cubic and hexagonal symmetries in the crystal structure.80–82 Figure 4 shows a comparison between q4 (i) and q6 (i) at 230, 250, and 270 K. The ice structures investigated are neatly classified and ice Ic appears to possess the highest cubical and hexagonal symmetries. Several spreading outliers may arise from the disordered structures involving an excess proton. Another low-density system, ice Ih, is less cubic-symmetrized than ice Ic, which is consistent with previous force field simulation results.45 For ice Ih and ice Ic, the standard deviation from the mean value increases in both bond order parameters as temperature becomes larger (see Tables S5–S8 in Supporting Information for statistical data). Hence, the entropic effect may contribute to loose cubic and hexagonal symmetries. The structure of ice III is characterized by a similar degree of cubic symmetry and a roughly halved hexagonality with respect to ice Ih. It is to be noted that no significant difference is observed for the pre-equilibrated structures at higher pressures (290 and 300 MPa). The melted ice VI yields a smaller q4 (i) value than ice III, resulting in the most disordered water configuration in this study. In Figure 4(c), the results calculated for bulk water with density of 0.997 g/cm3 are also plotted for comparison. The liquid phase retains intermediate cubic symmetry between ice III and melted ice VI.

11 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In summary, the signatures of the ice structures are obtained with the RDF and bond order parameters. We highlight that low- and high-density ice structures are clearly separated by the degree of hexagonal symmetry in the unit cell.

4.2. Dynamical Properties of Excess Proton in Ice In previous studies on bulk water,64,66 the proton diffusion coefficient was calculated by summing the Grotthuss diffusion coefficient, DG, and vehicular contribution was evaluated as a dynamical property. Nevertheless, in this study, we simply approximate the proton diffusion coefficient to DG because the translational motion of water molecules is negligibly small in ice Ih, ice Ic, and ice III. The proton diffusivity of melted ice VI, wherein water molecules can move more easily than in other ice structures, may be underestimated in this approximation, but the results have been left uncorrected to allow for systematic comparison. Figure 5 shows time-course accumulation functions (h(t) in eq. 5) calculated at different temperatures. In ice Ih, the frequency of proton transfer increases as temperature increases. Even at 50 K, the measurable slope seems to be in line with the experimental observations of proton diffusion in the low temperature region.14,83,84 The similarity of Grotthuss shuttling behavior at 250 and 270 K is likely associated with the acquisition of sufficient kinetic energy by water molecules. The tendency of h(t) in ice Ic is analogous to that in ice Ih with the exception of a lower amount of proton transfer at 50 and 250 K. Remarkably, Grotthuss shuttling is almost suppressed to zero at 50 K. Perhaps the higher cubic symmetry makes proton transfer hard to be activated when the amount of kinetic energy is extremely small. On the other hand, less frequent Grotthuss shuttling at 250 K seems to be a manifestation of the sensitivity of the proton transfer process to the disordered hydrogen bond network, considering that other trajectories starting

12 Environment ACS Paragon Plus

Page 12 of 43

Page 13 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

from different initial structures still yielded moderate slopes. The steepness of the slopes for ice III and melted ice VI is generally lower than that for ice Ih and ice Ic at respective temperatures. This indicates the facile proton transfer event in low-density structures, in other words, in ice with high hexagonal symmetry as discussed in Section 4.1. Conversely, the degree of cubicity seems to be less critical to the capability of Grotthuss shuttling due to comparable behavior of h(t) not only between ice Ih and ice Ic, but also between ice III and melted ice VI.

Table 1 shows the calculated DG values for the investigated ice structures. Overall, the diffusion coefficient becomes larger as temperature increases, indicating reasonable agreement with the experiment in ice Ih.84 The estimated values in ice Ih and ice Ic at 270 K are 2.64 and 3.17 times larger than DG in liquid water at 300 K,64 respectively. Although the previous theoretical study using ab initio-based reactive force field reported DG in ice Ih at 190 K as 0.52 Å2/ps,46 the present fast proton diffusivity indicates that the so-called trapping event46 is rather minor in the trajectories obtained with the DFTB method. The trapping event was caused by dissociation of a water molecule from [H3O(H2O)4]+ complex, which led to a relatively stable [H3O(H2O)3]+ complex.46 The later has a symmetric hydration environment, so that it has a long life-time and suppresses the proton transfer rate. The more frequent proton trapping event in the previous theoretical work leads to a smaller proton diffusivity. In contrast, the previous experimental work on the methanol doped ice reported 10 times larger DG values than that in liquid water. In such system, a crystal defect bigger than the present work is expected by introducing methanol moiety into the ice crystal, thus it leads to larger proton diffusivity.

13 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4.3. Energetic Properties of Excess Proton Diffusion Figure 6 shows the Arrhenius plot of diffusion coefficients, DG, in the investigated ice structures and liquid water.64 The plots for ice Ih and melted ice VI are well-fitted by straight lines, while those of ice Ic and ice III show insignificant deviation from a linear relationship. The estimated energy barriers are 1.50, 1.57, 10.30, 24.55, 8.69, and 14.00 kJ/mol for ice Ih, ice Ic (excluding the results at 50 and 250 K from the fit), ice III at 280, 290, and 300 MPa, and melted ice VI, respectively. In comparison to the energy barrier in liquid water (9.68 kJ/mol),64 lowdensity ice structures exhibit considerably smaller values. This is possibly related to the fact that hydrogen bond rearrangement (Zundel-Eigen transition)50,85,86 is not required in ordered crystal structures. The magnitude of the energy barrier in ice Ih agrees with experimental observations at high temperature regions, ~1 kJ/mol.84 The difference of 0.5 kJ/mol from the experimental diffusion barrier may come from the exclusion of the nuclear quantum effect in the present work. Another noticeable point is the larger DG values for high-density systems (except for ice III at 280 MPa) at 250 and 270 K in comparison to liquid water, despite the larger energy barriers of Grotthuss diffusion.

14 Environment ACS Paragon Plus

Page 14 of 43

Page 15 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

5. CONCLUSIONS We investigated proton diffusion in low- and high-density ice systems at temperatures in the range 50–270 K using the DFTB-MD simulation technique. Characterization of the crystal structures showed that the peak patterns of the fourth and subsequent solvation shells in radial distribution function between oxygen atoms can help categorize ice Ih and ice Ic and that these low-density variants possess higher hexagonal symmetry than ice III and melted ice VI, according to the bond order parameter analysis. The calculated diffusion coefficients and energy barriers of Grotthuss-type proton diffusion were compared with those of liquid water and reported values. The enhanced proton diffusivity from liquid water and low energy barrier in ice Ih and ice Ic followed the experimental trends. The suppressed proton transfer event in ice III and melted ice VI compared to ice Ih and ice Ic was rationalized in terms of the reduced hexagonal symmetry of crystal structures.

15 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ASSOCIATED CONTENT Supporting Information Description The Supporting Information is available free of charge on the ACS Publications website at DOI: … Unit cell sizes and densities of ice structures after classic equilibration, temperature dependence of O–O pair RDF for ice III at 290 and 300 MPa, and mean and standard deviation of bond order parameters.

ACKNOWLEDGMENTS This study was supported in part by MEXT as “Priority Issue on Post-K computer” (Development of new fundamental technologies for high-efficiency energy creation, conversion/storage and use) and Grant-in-aid for Scientific Research (A) “KAKENHI Grant Number JP26248009” from the Japan Society for the Promotion of Science (JSPS). Some simulations were partly performed at the Research Center for Computational Science, Okazaki, Japan. One of the authors (A.W.S.) acknowledges financial support from the Yoshida Scholarship Foundation.

16 Environment ACS Paragon Plus

Page 16 of 43

Page 17 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected]. Tel: +81-3-5286-3452.

Author Contributions The manuscript was written through contributions of all authors. All authors have given approval of the final version of the manuscript.

Notes The authors declare no competing financial interest.

17 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

REFERENCES (1) Chen, J.; Li, X.-Z.; Zhang, Q.; Michaelides, A.; Wang, E. Nature of Proton Transport in a Water-Filled Carbon Nanotube and in Liquid Water. Phys. Chem. Chem. Phys. 2013, 15, 6344–6349. (2) Habenicht, B. F.; Paddison, S. J.; Tuckerman, M. E. Proton Transfer in Perfluorosulfonic Acid Functionalized Carbon Nanotubes. In Polymer Materials and Membranes for Energy Devices; 2010, 1269, 7–12.

(3) Joseph, S.; Aluru, N. R. Why Are Carbon Nanotubes Fast Transporters of Water? Nano Lett. 2008, 8, 452–458. (4) Timmer, R. L. A.; Cox, M. J.; Bakker, H. J. Direct Observation of Proton Transfer in Ice Ih Using Femtosecond Spectroscopy. J. Phys. Chem. A 2010, 114, 2091–2101. (5) Bernal, J. D.; Fowler, R. H. A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl Ions. J. Chem. Phys. 1933, 1, 515–548. (6) Agmon, N. The Grotthuss Mechanism. Chem. Phys. Lett. 1995, 244, 456–462. (7) Cukierman, S. Et Tu, Grotthuss! And Other Unfinished Stories. Biochim. Biophys. Acta, Bioenerg. 2006, 1757, 876–885.

(8) Uritski, A.; Presiado, I.; Huppert, D. Indication of a Very Large Proton Diffusion in Ice Ih. J. Phys. Chem. C 2008, 112, 11991–12002.

(9) Ayrton, W. E.; Perry, J. Ice as an Electrolyte. Proc. Phys. Soc., London 1875, 2, 171–182. (10) Eigen, M. Proton Transfer, Acid-Base Catalysis, and Enzymatic Hydrolysis. Part I: ELEMENTARY PROCESSES. Angew. Chem. Int. Ed. 1964, 3, 1–19. (11) Pines, E.; Huppert, D. Geminate Recombination Proton-Transfer Reactions. Chem. Phys. Lett. 1986, 126, 88–91.

18 Environment ACS Paragon Plus

Page 18 of 43

Page 19 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(12) Pines, E.; Huppert, D.; Agmon, N. Geminate Recombination in Excited‐state Proton‐ transfer Reactions: Numerical Solution of the Debye–Smoluchowski Equation with Backreaction and Comparison with Experimental Results. J. Chem. Phys. 1988, 88, 5620– 5630. (13) Bartels-Rausch, T.; Bergeron, V.; Cartwright, J. H. E.; Escribano, R.; Finney, J. L.; Grothe, H.; Gutiérrez, P. J.; Haapala, J.; Kuhs, W. F.; Pettersson, et al. Ice Structures, Patterns, and Processes: A View across the Icefields. Rev. Mod. Phys. 2012, 84, 885–944. (14) Bove, L. E.; Klotz, S.; Paciaroni, A.; Sacchetti, F. Anomalous Proton Dynamics in Ice at Low Temperatures. Phys. Rev. Lett. 2009, 103, 165901–165904. (15) Collier, W. B.; Ritzhaupt, G.; Devlin, J. P. Spectroscopically Evaluated Rates and Energies for Proton Transfer and Bjerrum Defect Migration in Cubic Ice. J. Chem. Phys. 1984, 88, 363–368. (16) Devlin, J. P. Vibrational Spectra and Point Defect Activities of Icy Solids and Gas Phase Clusters. Int. Rev. Phys. Chem. 1990, 9, 29–65. (17) Lee, C.-W.; Lee, P.-R.; Kim, Y.-K.; Kang, H. Mechanistic Study of Proton Transfer and H/D Exchange in Ice Films at Low Temperatures (100–140 K). J. Chem. Phys. 2007, 127, 084701. (18) Hamada, I. A van der Waals Density Functional Study of Ice Ih. J. Chem. Phys. 2010, 133, 214503. (19) Santra, B.; Klimeš, J.; Tkatchenko, A.; Alfè, D.; Slater, B.; Michaelides, A.; Car, R.; Scheffler, M. On the Accuracy of van der Waals Inclusive Density-Functional Theory

19 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Exchange-Correlation Functionals for Ice at Ambient and High Pressures. J. Chem. Phys. 2013, 139, 154702. (20) Brandenburg, J. G.; Maas, T.; Grimme, S. Benchmarking DFT and Semiempirical Methods on Structures and Lattice Energies for Ten Ice Polymorphs. J. Chem. Phys. 2015, 142, 124104. (21) Liu, Y.; Ojamäe, K. Raman and IR Spectra of Ice Ih and Ice XI with an Assessment of DFT Methods. J. Phys. Chem. B 2016, 120, 11043–11051. (22) Pham, C. H.; Reddy, S. K.; Chen, K.; Knight, C.; Paesani, F. Many-Body Interactions in Ice. J. Chem. Theory Comput. 2017, 13, 1778–1784.

(23) Zhang, L.; Li, W.; Fang, T.; Li, S. Accurate Relative Energies and Binding Energies of Large Ice–Liquid Water Clusters and Periodic Structures. J. Phys. Chem. A 2017, 121, 4030–4038. (24) Kuo, J.-L.; Klein, M. L.; Kuhs, W. F. The Effect of Proton Disorder on the Structure of IceIh: A Theoretical Study. J. Chem. Phys. 2005, 123, 134505. (25) Li, F.; Skinner, J. L. Infrared and Raman Line Shapes for Ice Ih. I. Dilute HOD in H2O and D2O. J. Chem. Phys. 2010, 132, 204505. (26) Li, F.; Skinner, J. L. Infrared and Raman Line Shapes for Ice Ih. II. H2O and D2O. J. Chem. Phys. 2010, 133, 244504.

(27) Murray, E. D.; Galli, G. Dispersion Interactions and Vibrational Effects in Ice as a Function of Pressure: A First Principles Study. Phys. Rev. Lett. 2012, 108, 105502.

20 Environment ACS Paragon Plus

Page 20 of 43

Page 21 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(28) Zhang, J.; Kuo, J.-L.; Iitaka, T. First Principles Molecular Dynamics Study of Filled Ice Hydrogen Hydrate. J. Chem. Phys. 2012, 137, 084505. (29) He, X.; Sode, O.; Xantheas, S. S.; Hirata, S. Second-Order Many-Body Perturbation Study of Ice Ih. J. Chem. Phys. 2012, 137, 204505. (30) Gilliard, K.; Sode, O.; Hirata, S. Second-Order Many-Body Perturbation and CoupledCluster Singles and Doubles Study of Ice VIII. J. Chem. Phys. 2014, 140, 174507. (31) Geiger, P.; Dellago, C.; Macher, M.; Franchini, C.; Kresse, G.; Bernard, J.; Stern, J. N.; Loerting T. Proton Ordering of Cubic Ice Ic: Spectroscopy and Computer Simulations. J. Phys. Chem. C 2014, 118, 10989–10997.

(32) Svishchev, I. M.; Kusalik, P. G. Crystallization of Liquid Water in a Molecular Dynamics Simulation. Phys. Rev. Lett. 1994, 73, 975–978. (33) Matsumoto, M.; Saito, S.; Ohmine, I. Molecular Dynamics Simulation of the Ice Nucleation and Growth Process Leading to Water Freezing. Nature 2002, 416, 409–413. (34) Radhakrishnan, R.; Trout, B. L. Nucleation of Crystalline Phases of Water in Homogeneous and Inhomogeneous Environments. Phys. Rev. Lett. 2003, 90, 158301. (35) Radhakrishnan, R.; Trout, B. L. Nucleation of Hexagonal Ice (Ih) in Liquid Water. J. Am. Chem. Soc. 2003, 125, 7743–7747.

(36) Quigley, D.; Rodger, P. M. Metadynamics Simulations of Ice Nucleation and Growth. J. Chem. Phys. 2008, 128, 154518.

21 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(37) Moore, E. B.; Molinero, V. Ice Crystallization in Water's "No-Man's Land". J. Chem. Phys. 2010, 132, 244504. (38) Weiss, V. C.; Rullich, M.; Köhler, C.; Frauenheim, T. Kinetic Aspects of the Thermostatted Growth of Ice from Supercooled Water in Simulations. J. Chem. Phys. 2011, 135, 034701. (39) Rozmanov, D.; Kusalik, P. G. Temperature Dependence of Crystal Growth of Hexagonal Ice (Ih). Phys. Chem. Chem. Phys. 2011, 13, 15501–15511. (40) Li, T.; Donadio, D.; Russo, G.; Galli, G. Homogeneous Ice Nucleation from Supercooled Water. Phys. Chem. Chem. Phys. 2011, 13, 19807–19813. (41) Moore, E. B.; Molinero, V. Is It Cubic? Ice Crystallization from Deeply Supercooled Water. Phys. Chem. Chem. Phys. 2011, 13, 20008–20016.

(42) Malkin, T. L.; Murray, B. J.; Brukhno, A. V.; Anwar, J.; Salzmann, C. G. Structure of Ice Crystallized from Supercooled Water. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 1041–1045. (43) Johnson, J. C.; Molinero, V. Crystallization, Melting, and Structure of Water Nanoparticles at Atomospherically Relevant Temperatures. J. Am. Chem. Soc. 2012, 134, 6650–6659. (44) Reinhardt, A.; Doye, J. P. K. Free Energy Landscapes for Homogeneous Nucleation of Ice for a Monatomic Water Model. J. Chem. Phys. 2012, 136, 054501. (45) Reinhardt, A.; Doye, J. P. K.; Noya, E. G.; Vega, C. Local Order Parameters for Use in Driving Homogeneous Ice Nucleation with All-Atom Models of Water. J. Chem. Phys. 2012, 137, 194504.

22 Environment ACS Paragon Plus

Page 22 of 43

Page 23 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(46) Park, K.; Lin, W.; Paesani, F. Fast and Slow Proton Transfer in Ice: The Role of the QuasiLiquid Layer and Hydrogen-Bond Network. J. Phys. Chem. B 2014, 118, 8081–8089. (47) Kobayashi, C.; Saito, S.; Ohmine, I. Mechanism of Fast Proton Transfer in Ice: Potential Energy Surface and Reaction Coordinate Analyses. J. Chem. Phys. 2000, 113, 9090–9100. (48) Kobayashi, C.; Saito, S.; Ohmine, I. Mechanism of Fast Proton Transfer in Ice. II. Hydration, Modes, and Transport. J. Chem. Phys. 2000, 115, 4742–4749. (49) Marx, D. Proton Transfer 200 Years after von Grotthuss: Insights from Ab Initio Simulations. ChemPhysChem 2006, 7, 1848–1870. (50) Kirchner, B. Eigen or Zundel Ion: News from Calculated and Experimental Photoelectron Spectroscopy. ChemPhysChem 2007, 8, 41–43. (51) Eigen, M.; de Maeyer, L. Self-Dissociation and Protonic Charge Transport in Water and Ice. Proc. R. Soc. A 1958, 247, 505–533.

(52) Amir, W.; Gallot, G.; Hache, F.; Bratos, S.; Leicknam, J.-C.; Vuilleumier, R. TimeResolved Observation of the Eigen Cation in Liquid Water. J. Chem. Phys. 2007, 126, 034511. (53) Kulig, W.; Agmon, N. Both Zundel and Eigen Isomers Contribute to the IR Spectrum of the Gas-Phase H9O4+ Cluster. J. Phys. Chem. B 2014, 118, 278–286. (54) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260–7268.

23 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(55) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge DensityFunctional Tight-Binding Method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931– 948. (56) Seifert, G.; Joswig, J.-O. Density-Functional Tight Binding—an Approximate DensityFunctional Theory Method. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 456–465. (57) Hu, H.; Lu, Z.; Elstner, M.; Hermans, J.; Yang, W. Simulating Water with the SelfConsistent-Charge Density Functional Tight Binding Method: From Molecular Clusters to the Liquid State. J. Phys. Chem. A 2007, 111, 5685–5691. (58) Maupin, C. M.; Aradi, B.; Voth, G. A. The Self-Consistent Charge Density Functional Tight Binding Method Applied to Liquid Water and the Hydrated Excess Proton: Benchmark Simulations. J. Phys. Chem. B 2010, 114, 6922–6931. (59) Zheng, G.; Niklasson, A. M. N.; Karplus, M. Lagrangian Formulation with Dissipation of Born-Oppenheimer Molecular Dynamics Using the Density-Functional Tight-Binding Method. J. Chem. Phys. 2011, 135, 044122. (60) Goyal, P.; Elstner, M.; Cui, Q. Application of the SCC-DFTB Method to Neutral and Protonated Water Clusters and Bulk Water. J. Phys. Chem. B 2011, 115, 6790–6805. (61) Choi, T. H.; Liang, R.; Maupin, C. M.; Voth, G. A. Application of the SCC-DFTB Method to Hydroxide Water Clusters and Aqueous Hydroxide Solutions. J. Phys. Chem. B 2013, 117, 5165–5179.

24 Environment ACS Paragon Plus

Page 24 of 43

Page 25 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(62) Doemer, M.; Liberatore, E.; Knaup, J. M.; Tavernelli, I.; Rothlisberger, U. In Situ Parameterisation of SCC-DFTB Repulsive Potentials by Iterative Boltzmann Inversion. Mol. Phys. 2013, 111, 3595–3607.

(63) Goyal, P.; Qian, H.-J.; Irle, S.; Lu, X.; Roston, D.; Mori, T.; Elstner, M.; Cui, Q. Molecular Simulation of Water and Hydration Effects in Different Environments: Challenges and Developments for DFTB Based Models. J. Phys. Chem. B 2014, 118, 11007–11027. (64) Nakai, H.; Sakti, A. W.; Nishimura, Y. Divide-and-Conquer-Type Density-Functional Tight-Binding Molecular Dynamics Simulations of Proton Diffusion in a Bulk Water System. J. Phys. Chem. B 2016, 120, 217–221. (65) Nakai, H.; Nishimura, Y.; Kaiho, T.; Kubota, T.; Sato, H. Contrasting Mechanisms for CO2 Absorption and Regeneration Processes in Aqueous Amine Solutions: Insights from Density-Functional Tight-Binding Molecular Dynamics Simulations. Chem. Phys. Lett. 2016, 647, 127–131. (66) Sakti, A. W.; Nishimura, Y.; Nakai, H. Divide-and-Conquer-Type Density-Functional Tight-Binding Simulations of Hydroxide Ion Diffusion in Bulk Water. J. Phys. Chem. B 2017, 121, 1362–1371. (67) Sakti, A. W.; Nishimura, Y.; Sato, H.; Nakai, H. Divide-and-Conquer Density-Functional Tight-Binding Molecular Dynamics Study on the Formation of Carbamate Ions during CO2 Chemical Absorption in Aqueous Amine Solution. Bull. Chem. Soc. Jpn. 2017, 90, 1230– 1235.

25 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(68) Yu, H.; Cui, Q. The Vibrational Spectra of Protonated Water Clusters: A Benchmark for Self-Consistent-Charge Density-Functional Tight Binding. J. Chem. Phys. 2007, 127, 234504. (69) Miró, P.; Cramer, C. J. Water Clusters to Nanodrops: A Tight-Binding Density Functional Study. Phys. Chem. Chem. Phys. 2013, 15, 1837–1843. (70) Domínguez, A.; Niehaus, T. A.; Frauenheim, T. Accurate Hydrogen Bond Energies within the Density Functional Tight Binding Method. J. Phys. Chem. A 2015, 119, 3535–3544. (71) Feng, C.; Zhang, R. Q.; Dong, S. L.; Niehaus, T. A.; Frauenheim, T. Signatures in Vibrational Spectra of Ice Nanotubes Revealed by a Density Functional Tight Binding Method. J. Phys. Chem. C 2007, 111, 14131–14138. (72) Gaus, M.; Cui, Q.; Elstner, M. Density Functional Tight Binding: Application to Organic and Biological Molecules. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 49–61. (73) Yang, Y.; Yu, H.; York, D.; Cui, Q.; Elstner, M. Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method:  Third-Order Expansion of the Density Functional Theory Total Energy and Introduction of a Modified Effective Coulomb Interaction. J. Phys. Chem. A 2007, 111, 10861–10873. (74) Elstner, M. SCC-DFTB: What Is the Proper Degree of Self-Consistency? J. Phys. Chem. A 2007, 111, 5614–5621. (75) Abascal, J. L. F.; Sanz, E.; García Fernández, R.; Vega, C. A Potential Model for the Study of Ices and Amorphous Water: TIP4P/Ice. J. Chem. Phys. 2005, 122, 234511.

26 Environment ACS Paragon Plus

Page 26 of 43

Page 27 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(76) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (77) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781–1802.

(78) Nishizawa, H.; Nishimura, Y.; Kobayashi, M.; Irle, S.; Nakai, H. Three Pillars for Achieving Quantum Mechanical Molecular Dynamics Simulations of Huge Systems: Divide-and-Conquer,

Density-Functional

Tight-Binding,

and

Massively

Parallel

Computation. J. Comput. Chem. 2016, 37, 1983–1992. (79) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Bond-Orientational Order in Liquids and Glasses. Phys. Rev. B 1983, 28, 784–805. (80) Lechner, W.; Dellago, C. Accurate Determination of Crystal Structures Based on Averaged Local Bond Order Parameters. J. Chem. Phys. 2008, 129, 114707. (81) Moroni, D.; ten Wolde, P. R.; Bolhuis, P. G. Interplay between Structure and Size in a Critical Crystal Nucleus. Phys. Rev. Lett. 2005, 94, 235703–235706. (82) Coasne, B.; Jain, S. K.; Naamar, L.; Gubbins, K. E. Freezing of Argon in Ordered and Disordered Porous Carbon. Phys. Rev. B 2007, 76, 085416–085425. (83) Uritski, A.; Presiado, I.; Erez, Y.; Gepshtein, R.; Huppert, D. Temperature Dependence of Proton Diffusion in Ih Ice. J. Phys. Chem. C 2009, 113, 10285–10296.

27 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(84) Presiado, I.; Lal, J.; Mamontov, E.; Kolesnikov, A. I.; Huppert, D. Fast Proton Hopping Detection in Ice Ih by Quasi-Elastic Neutron Scattering. J. Phys. Chem. C 2011, 115, 10245–10251. (85) Kobayashi, C.; Iwahashi, K.; Saito, S.; Ohmine, I. Dynamics of Proton Attachment to Water Cluster: Proton Transfer, Evaporation, and Relaxation. J. Chem. Phys. 1996, 105, 6358–6366. (86) Ohmine, I.; Saito, S. Water Dynamics:  Fluctuation, Relaxation, and Chemical Reactions in Hydrogen Bond Network Rearrangement. Acc. Chem. Res. 1999, 32, 741–749.

28 Environment ACS Paragon Plus

Page 28 of 43

Page 29 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 1. Snapshots of the pre-equilibrated unit cell with an excess proton of (a) ice Ih, (b) ice Ic, (c) ice III, and (d) melted ice VI.

29 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. Temperature dependence of the O–O pair RDF for (a) ice Ih and (b) ice Ic. The black dashed line indicates the distance above which the shapes of the low-density variants are distinguishable from each other.

30 Environment ACS Paragon Plus

Page 30 of 43

Page 31 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 3. Temperature dependence of the O–O pair RDF for (a) ice III pre-equilibrated at 280 MPa and (b) melted ice VI.

31 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. Averaged form of cubic ( q4 (i ) ) and hexagonal ( q6 (i ) ) bond order parameters for ice structures and liquid water calculated at (a) 230, (b) 250, and (c) 270 K.

32 Environment ACS Paragon Plus

Page 32 of 43

Page 33 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 5. Temperature dependence of time-course accumulation functions for (a) ice Ih, (b) ice Ic, (c) ice III pre-equilibrated at 280 MPa, (d) ice III pre-equilibrated at 290 MPa, (e) ice III pre-equilibrated at 300 MPa, and (f) melted ice VI. The dashed lines indicate the slopes for estimating the proton transfer rate.

33 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6. Arrhenius plot of Grotthuss diffusion coefficients for ice structures. The previous simulation result for bulk water64 is also given for comparison.

34 Environment ACS Paragon Plus

Page 34 of 43

Page 35 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 1. Calculated Grotthuss diffusion coefficients (DG) (Å2/ps) for ice structures. Ice III T [K]

Ice Ih

Ice III

Ice III

Melted

Ice Ic (280 MPa) (290 MPa) (300 MPa)

ice VI

50

0.10

0.00

70

0.36

0.32

120

0.94

0.72

230

1.84

2.04

0.25

0.23

0.46

0.36

250

1.98

1.14

0.48

1.09

0.76

1.00

270

1.90

2.28

0.54

1.48

0.90

1.21

35 Environment ACS Paragon Plus

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphic

DFTB-MD Low-density DG High-density

36 Environment ACS Paragon Plus

Page 36 of 43

Page 37 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

FIg. 1 184x217mm (96 x 96 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

FIg. 2 191x277mm (96 x 96 DPI)

ACS Paragon Plus Environment

Page 38 of 43

Page 39 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

FIg. 3 134x194mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

FIg. 4 97x100mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 40 of 43

Page 41 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

FIg. 5 112x125mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

FIg. 6 793x396mm (96 x 96 DPI)

ACS Paragon Plus Environment

Page 42 of 43

Page 43 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table of contents 343x209mm (96 x 96 DPI)

ACS Paragon Plus Environment