Structure and Dynamics of Proton Transfer in Liquid Imidazole. A

Oct 1, 2012 - Ailin Li†, Zhen Cao†, Yao Li†, Tianying Yan†*, and Panwen Shen†. † Institute of New Energy Material Chemistry, Tianjin Key Laboratory of...
0 downloads 0 Views 1MB Size
Subscriber access provided by University of Guelph Library

Article

The Structure and Dynamic of Proton Transfer in Liquid Imidazole. A Molecular Dynamics Simulation Ailin Li, Zhen Cao, Yao Li, Tianying Yan, and Panwen Shen J. Phys. Chem. B, Just Accepted Manuscript • Publication Date (Web): 01 Oct 2012 Downloaded from http://pubs.acs.org on October 2, 2012

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 B 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 28

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

The Structure and Dynamic of Proton Transfer in Liquid Imidazole. A Molecular Dynamics Simulation Ailin Li,† Zhen Cao,†, ‡ Yao Li,† Tianying Yan,†,* and Panwen Shen† †

Institute of New Energy Material Chemistry, Tianjin Key Laboratory of Metal- and Molecule-Based Material Chemistry, Nankai University, Tianjin 300071, China ‡

Current Address: Department of Chemistry, University of Chicago, Chicago, Illinois 60637

*

Corresponding Author: [email protected]

RECEIVED DATE

TITLE RUNNING HEAD: MD Simulation of PT in Liquid Imidazole

ACS Paragon Plus Environment

1

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 2 of 28

Abstract.

Proton transfer (PT) via the Grotthuss mechanism in liquid imidazole (im) at 393 K is studied with molecular dynamics simulation using a reactive multistate empirical valence bond (MS-EVB) model. It is found that the proton is tightly binded to an imidazole to form an imidazolium (imH+), which is solvated in a distorted Eigen-like complex (im-imH+-im), while the Zundel-like complex (im-H+-im) is rare. PT occurs via an Eigen-Zundel-Eigen scenario for switching the identity of imH+ from an Eigen-like complex to another, intermediated by a Zundel-like complex. Structural and dynamical analyses demonstrate that PT in imidazole can be considered as a local event with very short spatial/temporal correlation, characterized by few “rattling” or recurrent PT events. At long time scale, the trend of PT correlation function may be recast with diffusion model of reversible geminate recombination toward the power-law decay. The formation of the hydrogen bonds (HBs) for the imidazole molecules between the first and second solvation shell of imH+ is crucial to pave the PT pathway. The above features may be understood by the flexibility of the HBs in liquid imidazole, as a stable HB network is essential for the Grotthuss mechanism.

Keywords. proton transfer; imidazole; Grotthuss mechanism; Eigen-like complex; Zundel-like complex

ACS Paragon Plus Environment

2

Page 3 of 28

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

I. Introduction. Imidazole (im) and im-based compounds have recently received considerable attention for their roles in the proton exchange membrane (PEM) fuel cells as proton transfer (PT) medium at the intermediate high temperature.1-2 Unlike the volatilization of water when the temperature is over 80 ℃, the boiling point of imidazole is as high as 256 ℃ and thus can be an alternative proton medium in the application of PEM fuel cells.1,3 PT between imidazole moieties occurs principally via the Grotthuss mechanism.4-5 The two nitrogen atoms on the imidazole ring can either donate (from the hydrogen-bearing nitrogen, NHB) or accept (by the non-hydrogen-bearing nitrogen, NNB) a proton. Indeed, the hopping of proton between donor and acceptor is fast once they are aligned along the PT pathway, as demonstrated by Car-Parrinello molecular dynamics (CPMD) simulations on the PT in crystalline imidazole and im-based compounds.6-8 In contrast to the relatively low ionic conductivity in crystalline imidazole or solid im-based polymers, the ionic conductivity of liquid imidazole is found to be several order of magnitude higher (~10-3 S/cm at melting point Tm = 90 ℃).1 The pioneer experimental study by Kreuer et al.3 shows that liquid imidazole displays excellent PT rate at the temperature of 120 ℃, which is comparable to the PT rate in water at room temperature, with the temperatures of both systems ~30 ℃ higher than their individual melting points. A reactive multistate empirical valence bond (MS-EVB) model for PT in liquid imidazole was developed,9 and MD simulation based on the MS-EVB model shows that the self-diffusion coefficient of the protonated center is 0.20 ± 0.01 Å2 ps-1 at 120 ℃, which is about two times higher than that of a vehicular imidazolium (imH+), in good agreement with the experimental measurements.3 This extra high proton conductivity is reasonable because in the Brønsted acid system charge migration can be achieved without large motion of molecules within the hydrogen bond (HB) network, such as proton conduction in phosphoric acid. According to the CPMD simulation, Vilčiauskas et al. found PT in phosphoric acid was driven by the formation of extended HB chains and subsequent solvent reorganization, which is very close to the Grotthuss mechanism.10 On the other hand, it was ACS Paragon Plus Environment

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 4 of 28

demonstrated in the MS-EVB simulation that protonated center highly localized on a specific imidazole,9 in distinct contrast to the solvated proton in liquid water.11-14 This is reasonable because the Brønsted acid system is generally characterized by the weak HB network.15 In the simulation of PT in liquid 2,4,5-trifluoroimidazole by Deng et al., it was found that PT was significantly retarded, due to the substitutions of small hydrogen atom by the bulker fluorine atom, which effectively prevents the proton donor and acceptor from approaching to each other.16 Therefore, PT rate is limited by the translational motion in liquid imidazole and its derivatives, and PT occurs spontaneously without barrier9,17-18 when the proton donor and acceptor approach to each other. The theoretical study by Mangiatordi et al. showed that the frustrated rotation of the protonated imidazole and the fluctuation of the polymer backbone played an important role in the PT process of poly(4vinyl-imidazole) (P4VI).19 The above studies highlight the importance of the surrounding environment and dynamics of the imH+ for PT in liquid imidazole. Compared to the less understood microscopic solvation structure and dynamics of PT in the nonaqueous solutions, the studies of PT in aqueous solution have a quite long history.20-22 Though there are still debates on the mechanism of PT in water, a clear scenario of it will improve our understanding on PT in the nonaqueous medium such as liquid imidazole. For this purpose, we give a brief summary for the recent progress on the mechanism of PT in water herein. The fast mobility of proton in water is attributed to the Grotthuss mechanism, which enhances the diffusion coefficient of proton by a factor of about 4 comparing to that of water.4,13 PT between water is barrierless,13,23,24 and PT in water follows Eigen-Zundel-Eigen (EZE) scenario between two protonated complexes, i.e., Eigen (H9O4+) and Zundel (H5O2+), with a distorted Eigen as the main hydration conformer.25-27 Due to the delocalization feature of the proton hydration in water, water molecules beyond the second solvation shell are involved in PT, as suggested by the experimental Infrared28-29 and dielectric30 spectroscopies, as well as the computer simulations.23,27,31-33 Consequently, PT in water is characterized by significant “rattling”, i.e., the fast “ballistic” motion of the proton between two water molecules which makes proton return to the original water molecule after two successive ACS Paragon Plus Environment

4

Page 5 of 28

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

PT,13,25,34-37 or the recurrence of PT back to the original carrier over longer time scale of several picoseconds,13,25 due to the strong spatial/temporal correlation of the HB network surrounding the charge defect created by the excess proton, with a few successful attempts to break the cage.13,25,34,36 A recent study38 demonstrated that hydrophobic solvent depresses PT in water by stabilizing the solvation structure of H3O+. The purpose of this study is aimed on the mechanism of PT in liquid imidazole, a Brønsted acid system, via MD simulation with an MS-EVB model.9 The comparisons with PT in water are made when they are appropriate. We find in this study that proton is mainly solvated in a distorted Eigenlike complex (im-imH+-im),6 as depicted in Figure 1(a), and PT occurs via the Grotthuss mechanism with an EZE scenario for switching the identity of imH+ from an Eigen-like complex to another, intermediated by a Zundel-like complex (im-H+-im),6-7 as depicted in Figure 1(b). The above finding is similar to PT in water. On the contrary, PT in imidazole can be considered as a local event with very short spatial/temporal correlation, characterized by few “rattling” or recurrent PT events which do not extend beyond the second solvation shell of imH+, comparing to the fact that more water molecules beyond the second solvation shell of H3O+ are involved for PT in water.31,39 More importantly, the formation of the HB for the imidazole molecules between the first and the second solvation shell of imH+ is crucial to pave the PT pathway. The above finding is in contrast to PT in bulk water, in which the breakage of HB of one water molecule in the first solvation shell is required for PT.4,40-41 On the other hand, a proton wire does exist in the confined environment of nano-pore4244

as well as in the water mediated acid-base PT reactions,45 in which PT along the proton wire is

ultrafast. We found in this study that such proton wire is also critical in understanding PT in imidazole. The above features may be understood by the flexibility of the HBs in liquid imidazole, as a stable HB network is essential for the Grotthuss mechanism. II. Methods. Details on the development of PT in liquid imidazole with the MS-EVB method as well as the force field parameters were described in the previous study.9 Briefly, MD simulations with MS-EVB ACS Paragon Plus Environment

5

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 6 of 28

model were performed for an excess proton solvated in 216 imidazole molecules with periodic boundary conditions (PBC) corresponding to the liquid density of 1.03 g/cm3. The production runs were performed with NVE ensemble at 393 K, and a total length of 22.7 ns trajectory was simulated, consisting of five independent runs. Another trajectory of 60 ns was simulated in order to analyze the time correlation functions, C(t), of PT events. We reuse the same trajectory of 60 ns with multiple time-origins to calculate C(t) for better statistics. In the MS-EVB method the system is represented by the linear combinations of the EVB state |i>, i.e., |Ψ>=∑i ci|i>, in which ci is the probability amplitude of the ith EVB state, and the adiabatic state is obtained by solving Schrödinger equation H|Ψ>=E|Ψ>. A two-state partition for the Zundel-like complex is depicted in Figure 1(b). For the above treatment, we define the center of excess charge rCEC as the weighted average of the center of charges (COCs) of individual EVB states,

N

rCEC = ∑ ci2 ri ,COC

(1)

i =1

in which ci2 is the probability and ri,COC is the COC of the imH+, respectively, of the ith EVB state. Classically, the imH+ may be viewed as the one in the EVB state of the highest probability c12, because the system likely collapses into the main EVB state of c12 ~ 0.97,9 as the imH+ depicted in Figure 1(a). III. Results and Discussion. III.1 The solvation of the proton in imidazole. In order to illustrate the solvation structure of proton in liquid imidazole, we analyze the two-dimensional distribution function P(δ,rNN), with δ the reaction coordinate defined by δ =|rCEC – r1,COC|-|rCEC – r2,COC|, in which rCEC is the CEC, r1,COC and r2,COC are the COC positions of the imH+ identified by the two EVB states9 of the highest probability (imH+) and the second highest probability, i.e., the tighter ligand imidazole (1x), and rNN is the distance between NHB on imH+ and NNB on imidazole (1x), respectively. The result is depicted in ACS Paragon Plus Environment

6

Page 7 of 28

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 2(a). The strong peaks of P(δ,rNN) appear at (±4.70 Å, 2.73 Å), which indicates the formation of the Eigen-like complex, in which the proton is tightly binded to the central imidazole to form imH+, as depicted in Figure 1(a). The negligible probability around δ ~ 0, where the CEC in the midway between two COCs for a typical Zundel-like complex as depicted in Figure 1(b), indicates that the Zundel-like complex is rare for proton in liquid imidazole. Figure 2(b) shows the partial radial distribution function (RDF), g(rNN), for NHB on imH+ and NNB on the neighboring imidazole molecules in the first solvation shell.25,46 Since the Eigen-like complex consists of two imidazole molecules around the imH+, we adopt the decomposition manner, inspired by the analyses utilized for PT in liquid water,25 to distinguish the contributions from the NNB, N1x and N1y. The subscripts 1x and 1y denote the tighter ligand imidazole (1x) and the looser ligand imidazole (1y), respectively. The overall peak (black line) appears at 2.73 Å, with the peak of N1x at 2.69 Å and N1y of 2.83 Å. These peak distances are shorter than the rNN peak of 2.9 Å found in liquid imidazole,47 indicating the stronger HB between imH+ and imidazole than that between the imidazole molecules. The stronger HB associated with imH+ can also be found by comparing Figure 3(a) and 3(b), in which with rNN and ∠NNH angle peaked at (2.73 Å, 8.5°), (2.99 Å, 18°), respectively, for imH+-im and im-im. The rNN of the gas-phase optimized Eigen-like and Zundel-like complexes are 2.55 Å and 2.69 Å, respectively, with B3LYP/6-311G** level of theory.9 Therefore, imH+ is solvated in a distorted Eigen-like complex, as depicted in Figure 1(a), with one rNN shorter than the other, as similar to the distorted Eigen found in water.25 Figure 2(c) shows the probability distribution, P(δ) = ∫drNNP(δ,rNN), in which P(δ,rNN) is depicted in Figure 2(a). It can be seen that P(δ) is ~ 2×10-5 around the region of |δ| < 0.5, where the proton is almost equally shared by two neighboring imidazole molecules, which can be attributed to a Zundellike complex as in Figure 1(b). Therefore, the occurrence of the Zundel-like complex is indeed rare event, and this is in contrast to proton solvated in bulk water13,23,25,46 as well as in nanoconfined water chain44,48-50 with apparent Zundel populations. Figure 2(d) shows the free energy profile,

ACS Paragon Plus Environment

7

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

∆F (δ ) = −kbT ln P(δ )

Page 8 of 28

(2)

in which kb is Boltzmann constant and T = 393 K is the temperature, and the zero free energy corresponds to the Eigen-like complex of the maximum P(δ). It may be deduced from Figure 2(d) that proton hops from an Eigen-like complex of δ ~ ± 4.7 Å to another Eigen-like complex with the identity of imH+ switched, via a Zundel-like complex of δ ~ 0. Therefore, PT in imidazole can be attributed to the EZE scenario, similar to PT in water.25-26 The free energy barrier of PT in imidazole is ~6.8 kcal/mol, which is much higher than that of PT in water of ~1 kcal/mol or lower when nuclear quantum effect is taken into account.13,24,51 On the other hand, it is also notable that rNN distribution peaked at 2.69 Å in Figure 2(b), and at such distance PT is essentially barrierless.9 Therefore, the high free energy barrier in Figure 2(d) may be correlated to the dynamics of the solvation structure of neighboring imidazoles. III.2 The dynamics of proton transfer in imidazole. Based on the previous studies, PT process is correlated with the HB dynamics.31 Also PT via the Grotthuss mechanism requires a relatively stable HB network. We next analyze HB correlation function in imidazole,52-53

CijH ( t ) =

pij ( 0 ) h ( 0 ) h ( t ) pij h

(3)

in which h(t) = 1 if a geminate imidazole pair is hydrogen bonded at time t, and h(t) = 0 otherwise. Since the system consists of imidazoles and one proton, a selection function, pij(0), is utilized to distinguish different HB pairs at t = 0, i.e. p01(0) = 1 for HB pair between imH+ and the imidazole in the first solvation shell, p12(0) = 1 for HB pairs within the first and second solvation shell of imH+, and p>2(0) = 1 for HB pairs beyond the second solvation shell of imH+, respectively, with pij(0) = 0 otherwise. A tight HB criterion of (rNN, ∠NNH) within (3.2 Å, 30°) is adopted for p01(0) = 1, while

ACS Paragon Plus Environment

8

Page 9 of 28

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

a loose HB criterion of (rNN, ∠NNH) within (3.2 Å, 45°) for p12(0) = 1 and p>2(0) = 1, respectively, according to the two-dimensional (rNN, ∠NNH) distribution shown in Figure 3. CijH(t) is calculated based on the time segment when the imH+ identity does not switch for ≥30 ps, the average time between PT events.3,9 The resulting C01H(t), C12H(t) and C>2H(t) are depicted in Figure 4. It is notable that C01H(t) decays considerably slower than C12H(t) and C>2H(t), while the latter two correlation functions relax at similar rate. Also, the bumps characterizing the HB librational motion appear at ~0.06 ps in C01H(t) and ~0.2 ps in C12H(t) and C>2H(t), indicating the stronger HB between imH+ and imidazole with shorter time scale for the HB librational motion. The above comparisons show the stronger HB between imH+ and the imidazole in the first solvation shell. The similar relaxation for C12H(t) and C>2H(t) highlights the spatial localization feature of the proton, which does not extend beyond the first solvation shell of imH+. Such localization feature is different from proton solvated in water, in which the charge defect is delocalized over larger water cluster extended beyond the second solvation shell of H3O+.13,25,31-33 In liquid imidazole, the Eigen-like complex is relatively stable due to the stronger HB between imH+ and imidazole with weaker HB beyond the first solvation shell. Apart from that, it can be seen that the HB correlation function exhibits a power-law decay at long time, which can be attributed to the reversible geminate recombination due to the diffusion of the molecules.21,54 Comparing to the average time between PT events, the relaxation of HB correlation functions in Figure 4 is much faster, with CijH(t) decays an order of magnitude faster than the average ~30 ps. This may be the reason for the slow dynamics of PT. Due to the fast changing HB around the cation, it is less probable for the components involved to reach the proper configuration for PT process. Also the fast relaxation of HB may result in the short memory of proton on its original carrier after PT occurs i.e., a characteristics of short spatial/temporal correlation for PT in imidazole. Once a new imidazolium cation is formed, the structure of the solvation shell soon changes and makes it less probable for the back reaction. And this can be verified by the identity change of the imH+ cation, which is shown in Figure 5(a). As is shown, the abruption of the identity line marks a PT event. 9 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

Page 10 of 28

Compared to bulk water, PT in liquid imidazole is much slower and with less recurrent PT events, and it is not as frequent as that in water.25 The recurrent PT event occurs within very short time segment (less than 0.2 ps), as shown in the inset. No frequent exchange of imH+ identity is found, which implies the Zundel-like complex is rare in imidazole. It is notable that the Eigen-like complex is characterized by the stronger HB between imH+ and its neighboring imidazoles than that between the imidazoles, and an extended HB network does not exist in liquid imidazole. In order to detect the relationship between the PT process and the dynamics of the HB network, the time evolution of the cluster size centered on imH+ maintained by the HBs is investigated. In searching the cluster, the loose HB criterion is adopted. The results are shown in Figure 5(b) for the whole cluster size (black line), the size starting from im1 in the inset against the PT direction (blue line), and the size starting from im2 along the PT direction (red line), respectively, for 1 ps before and after PT. The time origin of PT, t = 0 ps, is defined to be the last time when the proton retained on im1. The vertical dashed lines mark 0.2 ps, the characteristic time for the librational motion of the HB between imidazoles shown in Figure 4, before and after PT. It is clearly shown in Figure 5(b) that the elongation of the imidazole chain along the PT direction is accompanied with the decrease of the imidazole chain against the PT direction, with the red line and the blue line being nearly symmetric within statistical uncertainty, while the size of the whole cluster remains to be 3.9 to 4.0 in the PT process. The HBs between the imidazole molecules in the first and the second solvation shell of imH+ easily breaks and forms. Thus the PT events should be correlated to the change of distance between the nitrogen atoms that forms the HBs. PT involves the switching of the identity of imH+ for the transition from one Eigen-like complex to another. As depicted in Figure 6(a), the proton hops from N2 to N3, and the imH+ switches from im1 to im2. Accordingly, im2 or im1 becomes the imidazole in the first solvation shell of imH+, respectively. By applying the time dependent radial distribution functions (TD-RDFs),25,35 we calculated the nitrogen-nitrogen TD-RDFs, g(rNN)’s, during different time segments to show the time evolution of the solvation structure. Here the time segments were chosen to be 0-0.2 ps, 0.2-0.4 ps, ACS Paragon Plus Environment

10

Page 11 of 28

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

and 0.8-1.0 ps, respectively, before and after PT. The time origin of PT, t = 0 ps, is defined to be the last time when the proton retained on im1, which is the same as in Figure 5(b). In Figure 6, the left panels and the right panels represent the g(rNN)’s before and after PT, respectively. Figure 6(b) and 6(c) show g(rNN)’s between proton donor N2 and acceptor N3 (centered on N2 and only N2 and N3 involved), Figure 6(d) and 6(e) are g(rNN)’s against the PT direction (centered on N1), as well as in Figure 6(f) and 6(g) for g(rNN)’s along the PT direction (centered on N4). In Figure 6(d) to 6(g), when centering on the NHBs (N1 and N4), only the TD-RDFs with the NNBs of different imidazoles were calculated. That is, g(rNN) with N5, N3, N6, etc., and g(rNN) with N6, N2, N5 etc, are calculated, when centered on N1 and N4, respectively. It should be noted that in the above assignment, N3 is treated as NNB when centered on N1, and N2 is treated as NNB when centered on N4, no matter PT occurs or not. It can be seen from Figure 6(b) and 6(c) that the N2-N3 peak at shorter distances with increasing intensities before PT, while the above trend reverses after PT. At time segments of 0-0.2 ps before and after PT, the position of the peak is located at ~2.55 Å. Such rNN highlights the formation of Zundel-like complex when PT occurs.9 Also, the position and intensity of the peak for 0.8-1.0 ps before and after PT are in good agreement with the g(rNN) for imH+ with its tighter ligand imidazole (1x) in Figure 2(b), so that the Zundel-like complex quickly switches to the distorted Eigen-like complex. In the above analyses, it becomes apparent that the mechanism of PT in liquid imidazole can be attributed to the EZE scenario, similar to PT in water.13,25,27 On the other hand, as can be seen from Figure 6(a), im1 and im2 remain two-coordinated during the PT process. This is in contrast to PT in water, in which the coordination number of the proton acceptor switches from four to three.4,25,35 Such difference seen in the switching of the coordination number during PT may be attributed to the single lone-pair (LP) in the sp2 hybrid nitrogen atom in imidazole, while two LPs in the sp3 hybrid oxygen atom in water.21 Since PT correlates with the dynamics of the first and the second solvation shell, it is of interest to study the dynamics of the solvation structure around imH+ via TD-RDFs. It can be seen from Figure ACS Paragon Plus Environment

11

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 12 of 28

6(d) and 6(e) that the first peaks of the g(rNN)’s, attributed to the N1-N5 spatial correlation, peak at longer distances with decreasing intensities, as approaching PT from t = -1.0 ps to t = 1.0 ps. It suggests that the N5-carrying imidazole, which was in the first solvation shell of im1, tends to move away from im1, as proton hopping from im1 to im2. On the contrary, the first peaks of the g(rNN)’s in Figure 5(f) and 5(g), attributed to the N4-N6 spatial correlation, peak at shorter distances with increasing intensities, as approaching PT from t = -1.0 ps to t = 1.0 ps. Therefore, the N6-carrying imidazole, which was in the second solvation shell of im1, tends to move toward im2, as proton hopping from im1 to im2. Another feature from Figure 6(d) and 6(g) is that when the first peak of g(rNN) represents imH+ with the imidazole in its first solvation shell, it peaks at ~2.7 Å, which is in agreement with the g(rNN) for imH+ with its tighter ligand imidazole (1x) in Figure 2(b). On the other hand, when the first peak of g(rNN) represents the imidazoles between the first and the second solvation shell of imH+, as shown in Figure 6(e) and 6(f), it peaks at ~3.0 Å. Comparing to the NHBNNB RDF in liquid imidazole,47 in which the first peak appears at ~ 2.9 Å, it can be seen that the effect of the excess proton does not extend beyond the first solvation shell, highlighting the spatial localization feature of PT in imidazole. By taking a closer inspect on Figure 6(b) to 6(g), it can be seen that the intensities and peak positions of the g(rNN)’s at intermediate time (0.2-0.4) ps resemble those at longer time (0.8-1.0 ps) more than those at short time (0-0.2 ps) before and after PT. Thus, the change in the solvation structure is seen to be most drastic within 0.2 ps before and after PT. Previous simulation shows that the PT in aqueous solution needs the HB interaction even extended to the water molecules beyond the first solvation shell around the hydronium cation.31,33 The above analyses show that in liquid imidazole the formation of the HB in the second solvation shell along the PT direction is essential to pave the PT pathway. And due to the flexibility of the HB network in imidazole, PT occurs within very short time. The above observation highlights the temporal localization feature within sub-picosecond of PT in imidazole.

III.3 The kinetics of PT in imidazole. The PT events can be identified by monitoring the distance between rCEC and the center of the pivot imidazole at t = 0 ps. The kinetics of PT events in ACS Paragon Plus Environment

12

Page 13 of 28

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

imidazole is studied with two different time correlation functions, which are continuous PT correlation function Cc(t) and intermittent PT correlation function Ci(t),

Cc ( t ) =

and Ci ( t ) =

h (0) H (t ) h

,

h (0) h (t ) h

(4a)

(4b)

In the above equations, we choose the distance between rCEC and the center of the imidazole ring (CI) of the pivot imidazole, rCI, as the geminate pair following the approach by Agmon et al..32 Accordingly, the Boolean function is defined as h(t) = 1 if rCEC-CI(t) < rcut at time t, and h(t) = 0 otherwise; whereas H(t) = 1 if rCEC-CI(t) retains to be less than rcut continuously up to time t, and H(t) = 0 otherwise. The cutoff radius rcut is set to be 2.9 Å, which is the first minimum of the RDF between rCEC and rCI of the pivot imidazole, as shown in figure S1 in the supporting information. Extending rcut to 3.0 Å makes no difference in Cc(t) and Ci(t), because the appearance of the low probability in this region. When calculating the time correlation functions, all the molecules are allowed to diffuse in the expanded space beyond PBC, since in a finite system of N molecules confined in PBC box, the time correlation functions converge to 1/N rather than to zero, and this artificial result may affect the observation of the asymptotic behavior.32,54 It should be noted that the above definitions of h(t) and H(t) are equivalent to the charge population functions by Tuckerman et al.,34,36 though the selection of the protonated molecule is slightly different, because CEC may not be easily defined in the AIMD simulations. The time correlation functions are shown in Figure 7(a) for Cc(t) and (7b) for Ci(t). For the definition of Cc(t) in Eq. (4a), the reaction may be written as34,36 kd k2  → im L H +  imH + ← → im + H +  ka

ACS Paragon Plus Environment

(5a)

13

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 14 of 28

in which H+ represents CEC, which hops from an imidazole to another across the dividing surface of rcut = 2.9 Å when PT occurs, with kd and ka the intrinsic dissociation and association rate constants, respectively. The original pivot imidazole retains in the first solvation shell of CEC, i.e., im---H+, for some time, and then a complete reaction of rate constant k2 (or k2 + D, in which D denotes diffusion) separates the geminate pair.34,36 The second step is irreversible, as indicated by the one-direction arrow of rate constant k2, in accordance with the definition of H(t). For such reaction, the solution of the coupled linear differential equation shows that Cc(t) exhibits a biexponential decay.34,36 Interestingly, it can be seen that single exponential decay fits reasonably well to Cc(t), as shown in Figure 7(a) of the blue dashed line for fitting Cc(t) with exp(-t/τc), with τc = 29.73 ps. Such time scale is about two order of magnitude longer than that in water.36 Compared to the biexponential decay of PT in water, the single exponential decay in liquid imidazole may imply the dominance of slow process of PT and the ballistic motion of proton in the Zundel-like complex is not prevailing. It can be seen that single exponential decay fits reasonably well to Cc(t), indicating that the slow process dominates PT dynamics with negligible contribution from the fast process. For PT in water, Tuckerman et al. attributed the slow and fast contributions from Eigen and Zundel, respectively.34,36 Here, the lack of fast process in Cc(t) highlights that Zundel-like complex is rare in imidazole, in good agreement with Figure 2. Also, the fitted τc is in good agreement with the average time of ~30 ps for proton hopping estimated experimentally.3 Meanwhile, it is notable that Ci(t) decays much slower than Cc(t) as comparing Figures 7(b) and 7(b). The fitting to Ci(t) with biexponential decay, aexp(-t/τ1)+(1-a)exp(-t/τ2), gives a = 0.15, τ1 = 35.90 ps, τ2 = 80.69 ps, respectively. It can be seen that τ1 is similar to τc, which characterizes the history dependent recombination of the geminate pair, while τ2 may be attributed to the history independent recombination of the geminate pair. The biexponential decay fits Ci(t) reasonably well up to approximately 300 ps, when Ci(t) has decayed to smaller than 2% of its initial value of 1. In the region beyond 300 ps, Ci(t) exhibits a nonexponential decay, which is slower than the biexponential decay for the reaction proposed in Eq. (5a). ACS Paragon Plus Environment

14

Page 15 of 28

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

The difference seen in Cc(t) and Ci(t) is reasonable, for the former does not take into account the history independent recombination as well as the long time recombination of the geminate pair due to the slow diffusion, while the latter does. Based on the above consideration, the reaction of Ci(t) defined in Eq. (4b) may be written as32,54 kd D  → im L H + ← imH + ← → im + H +  ka

(5b)

in which H+ represents CEC, as in Eq. (5a), . The difference between Eq. (5a) and (5b) is seen for the second step, in which the contribution from the diffusion of the geminate pair for the long time recombination is taken into account. Such diffusion controlled mechanism of the reversible geminate recombination has been discussed by Eigen et al. due to the high intrinsic rate constant of PT in water.55-56 Agmon et al. has developed a kinetic model based on diffusion equation with appropriate boundary condition at the dividing surface and provided analytical solutions to the free diffusion model.57 It demonstrates a power-law tail that is proportional to t-3/2 for the long time decay of Ci(t).32,58-61 Alternatively, the diffusion model in Eq. (5b) can be solved numerically by the freely available software (SSDP, ver. 2.66)62 for the spherically symmetric diffusion equation, which has been applied successfully on the PT in water,32 as well as the power-law decay of the HB correlation function in water.21,54 For PT in water, Ci(t) shows an initial biexponential decay up to 30 ps, followed by a transient region between 30 ps to 100 ps, and the power-law tail of the t-3/2 becomes evident after 100 ps.32 Here, due to the slow PT rate in liquid imidazole, deviation from the biexponential decay is not clear until ~ 500 ps. As is shown in Figure 7(b), Ci(t) does not provide good statistical results below 10-3, but the trend corresponds well to the solution to the diffusion model calculated by the SSDP software,62 as shown in the red dashed line in Figure (7b). For the diffusion model in Eq. (5b), the parameters for dissociation and association rate, kd and ka/(4πR2) are 0.016 ps-1 and 0.03 Å ps-1, respectively. Herein R represents the dividing surface for PT reaction, and R = 2.9 Å for PT in imidazole as shown in the inset of Figure 7(a), and R = 1.6 Å for PT in water.32 The relaxation time scale kd-1 is 62.5 ps, which is comparable to the τ2 in the biexponential fitting as

ACS Paragon Plus Environment

15

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 16 of 28

the blue dashed line in Figure 7(b). Compared to the rate parameters of PT in liquid water (kd = 1.07 ps-1 and ka/(4πR2) = 1.6 Å ps-1 ) from the same diffusion model in Eq. (5b),32 both the dissociation and association in liquid imidazole are two orders of magnitude slower, it can be deduced that PT in liquid imidazole is much less frequent, but accompanied by much less recurrent PT process, in agreement with previous analyses. For the long-time region the Ci(t), the deviation from the biexponential decay is clear and the power-law tail of t-3/2 becomes evident. Thus, the long-time asymptotic power-law decay of Ci(t) seems to be common for the reversible geminate pairs.

VI. Summary. In this study, we performed MD simulations with MS-EVB model on an excess proton solvated in liquid imidazole. It is found that proton is tightly binded to an imidazole to form an imH+, which is solvated in a distorted Eigen-like complex (im-imH+-im). PT occurs via an EZE scenario for switching the identity of imH+ from an Eigen-like complex to another, intermediated by a Zundellike complex, with a free energy barrier of ~6.8 kcal/mol. Apart from that, PT in liquid imidazole is characterized by fast switching of the imidazoles, especially those around the second solvation shell around imH+, as the HB correlation function relaxes much faster than the PT correlation function by one order of magnitude. The high flexibility of imidazoles around imH+ accounts for much less “rattling” or recurrence of PT, and proton does not retain the memory of its original carrier after 0.2 ps it hops. Therefore, PT in liquid imidazole is a local event with very short spatial/temporal correlation. The above finding highlights the intrinsic high rate constant of PT as imH+ and imidazole approach to each other to a favorite rNN for PT to occur. At long time scale, the trend of PT correlation function may be recast with diffusion model of revisable geminate recombination toward the power-law decay. Based on the above observations, the dynamics of the second solvation shell is the key factor to determine the PT rate. The analyses on the TD-RDF of the second solvation shell around imH+ shows that PT is indeed correlated strongly with the dynamics of the second solvation shell, with the peak along the PT direction enhanced and that against the PT direction depressed immediately before ACS Paragon Plus Environment

16

Page 17 of 28

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

and after PT. In addition, the time evolution of cluster size centered on imH+ shows that the average cluster size is about 3.9 to 4.0, which consists of a relatively stable Eigen-like complex, with continuous formation and breakage of HBs for the imidazole in the second solvation shell. Within sub-picosecond when PT occurs, the cluster size along the PT direction increases and imidazole molecules against the PT direction decreases. Therefore, a PT pathway is paved, and a new Eigenlike complex is formed with the proton located on the next imidazole. The cluster size maintained by HB remains the same during PT process, as can be attributed to the relatively stable Eigen-like complex. The above discussion on the PT mechanism highlights the importance on the formation of HB with the imidazoles in the first solvation shell of imH+. In liquid imidazole, the solvation shell around the imH+ switches very rapidly, and it is difficult for the formation of a proton wire, and this may account for the relatively slow PT rate in imidazole. We find in this study that a stable second solvation shell around imH+ is crucial for PT in liquid imidazole. We hope that the above finding may improve the understanding on the mechanism of PT in Brønsted acid system toward further development on the anhydrous electrolytes for fuel cells.

Acknowledgement. This work is supported by NSFC (21073097), Natural Science Foundation of Tianjin (12JCYBJC13900), and NCET-10-0512. The simulations were performed on TianHe-1(A) at National Supercomputer Center in Tianjin, China. The authors are grateful to an anonymous reviewer for indicating the reversible geminate recombination toward the long time power-law decay of the time correlation functions.

Supporting Information Available. Figure S1 for the radial distribution function between the center of excess charge (CEC) and the center of imidazole ring (CI). This material is available free of charge via the Internet at http://pubs.acs.org.

References.

ACS Paragon Plus Environment

17

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 18 of 28

(1) Kreuer, K. D.; Paddison, S. J.; Spohr, E.; Schuster, M. Chem. Rev. 2004, 104, 4637-4678. (2) Asensio, J. A.; Sánchez, E. M.; Gómez-Romero, P. Chem. Soc. Rev. 2010, 39, 3210-3239. (3) Kreuer, K. D.; Fuchs, A.; Ise, M.; Spaeth, M.; Maier, J. Electrochim. Acta 1998, 43, 1281-1288. (4) Agmon, N. Chem. Phys. Lett. 1995, 244, 456-462. (5) Agmon, N.; Goldberg, S. Y.; Huppert, D. J. Mol. Liq. 1995, 64, 161-195. (6) Iannuzzi, M. J. Chem. Phys. 2006, 124, 204710. (7) Iannuzzi, M.; Parrinello, M. Phys. Rev. Lett. 2004, 93, 025901. (8) Münch, W.; Kreuer, K. D.; Silvestri, W.; Maier, J.; Seifert, G. Solid State Ionics 2001, 145, 437443. (9) Chen, H.; Yan, T.; Voth, G. A. J. Phys. Chem. A 2009, 113, 4507-4517. (10) Vilčiauskas, L.; Tuckerman, M. E.; Bester, G.; Paddison, S. J.; Kreuer, K. D. Nat. Chem. 2012, 4, 461-466. (11) Voth, G. A. Acc. Chem. Res. 2006, 39, 143-150. (12) Swanson, J. M. J.; Maupin, C. M.; Chen, H.; Petersen, M. K.; Xu, J.; Wu, Y.; Voth, G. A. J. Phys. Chem. B 2007, 111, 4300-4314. (13) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2008, 112, 467-482. (14) Knight, C.; Voth, G. A. Acc. Chem. Res. 2012, 45, 101-109. (15) Schuster, M.; Kreuer, K. D.; Steininger, H.; Maier, J. Solid State Ionics 2008, 179, 523-528. (16) Deng, W.-Q.; Molinero, V.; Goddard III, W. A. J. Am. Chem. Soc. 2004, 126, 15644-15645. (17) Viswanathan, U.; Basak, D.; Venkataraman, D.; Fermann, J. T.; Auerbach, S. M. J. Phys. Chem. A 2011, 115, 5423-5434. (18) Mangiatordi, G. F.; Hermet, J.; Adamo, C. J. Phys. Chem. A 2011, 115, 2627-2634. (19) Mangiatordi, G. F.; Butera, V.; Russo, N.; Laage, D.; Adamo, C. Phys. Chem. Chem. Phys.

2012, 14, 10910-10918. (20) Marx, D. ChemPhysChem 2006, 7, 1848-1870. Addendum: ibid. 1848(1842):1209-1210. (21) Agmon, N. Acc. Chem. Res. 2012, 45, 63-73. ACS Paragon Plus Environment

18

Page 19 of 28

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

(22) Kornyshev, A. A.; Kuznetsov, A. M.; Spohr, E.; Ulstrup, J. J. Phys. Chem. B 2003, 107, 33513366. (23) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. Nature 1999, 397, 601-604. (24) Marx, D.; Tuckerman, M. E.; Parrinello, M. J. Phys.: Condens. Matter 2000, 12, A153-A159. (25) Markovitch, O.; Chen, H.; Izvekov, S.; Paesani, F.; Voth, G. A.; Agmon, N. J. Phys. Chem. B

2008, 112, 9456-9466. (26) Lee, S. H.; Rasaiah, J. C. J. Chem. Phys. 2011, 135, 124505. (27) Markovitch, O.; Agmon, N. J. Phys. Chem. A. 2007, 111, 2253-2256. (28) Stoyanov, E. S.; Stoyanova, I. V.; Reed, C. A. J. Am. Chem. Soc. 2010, 132, 1484-1485. (29) Stoyanov, E. S.; Stoyanova, I. V.; Reed, C. A. Chem. Sci. 2011, 2, 462-472. (30) Tielrooij, K. J.; Timmer, R. L. A.; Bakker, H. J.; Bonn, M. Phys. Rev. Lett. 2009, 102, 198303. (31) Lapid, H.; Agmon, N.; Petersen, M. K.; Voth, G. A. J. Chem. Phys. 2005, 122, 014506. (32) Chen, H.; Voth, G. A.; Agmon, N. J. Phys. Chem. B 2010, 114, 333-339. (33) Day, T. J. F.; Schmitt, U. W.; Voth, G. A. J. Am. Chem. Soc. 2000, 122, 12027-12028. (34) Chandra, A.; Tuckerman, M. E.; Mark, D. Phys. Rev. Lett. 2007, 99, 145901. (35) Berkelbach, T. C.; Lee, H. S.; Tuckerman, M. E. Phys. Rev. Lett. 2009, 103, 238302. (36) Tuckerman, M. E.; Chandra, A.; Mark, D. J. Chem. Phys. 2010, 133, 124108. (37) Woutersen, S.; Bakker, H. J. Phys. Rev. Lett. 2006, 96, 138305. (38) Bonn, M.; Bakker, H. J.; Rago, G.; Pouzy, F.; Siekierzycka, J. R.; Brouwer, A. M.; Bonn, D. J. Am. Chem. Soc. 2009, 131, 17070-17071. (39) Agmon, N. J. Phys. Chem. A 2005, 109, 13-35. (40) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1995, 103, 150-161. (41) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Phys. Chem 1995, 99, 5749-5752. (42) Brewer, M. L.; Schmitt, U. W.; Voth, G. A. Biophys. J. 2001, 80, 1691-1702. (43) Dellago, C.; Naor, M. M.; Hummer, G. Phys. Rev. Lett. 2003, 90, 105902.

ACS Paragon Plus Environment

19

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 20 of 28

(44) Cao, Z.; Peng, Y.; Yan, T.; Li, S.; Li, A.; Voth, G. A. J. Am. Chem. Soc. 2010, 132, 1139511397. (45) Cox, M. J.; Timmer, R. L. A.; Bakker, H. J.; Park, S.; Agmon, N. J. Phys. Chem. A 2009, 113, 6599-6606. (46) Asthagiri, D.; Pratt, L. R.; Kress, J. D. Proc. Nat. Acad. Sci. USA 2005, 102, 6704-6708. (47) Shaik, M. S.; Liem, S. Y.; Yuan, Y.; Popelier, P. L. A. Phys. Chem. Chem. Phys. 2010, 12, 15040-15055. (48) Freier, E.; Wolf, S.; Gerwert, K. Proc. Nat. Acad. Sci. USA 2011, 108, 11435-11439. (49) Kaila, V. R. I.; Hummer, G. Phys. Chem. Chem. Phys. 2011, 13, 13207-13215. (50) Kofinger, J.; Hummer, G.; Dellago, C. Phys. Chem. Chem. Phys. 2011, 13, 15403-15417. (51) Tuckerman, M. E.; Mark, D.; Klein, M. L.; Parrinello, M. Science 1997, 275, 817-820. (52) Luzar, A.; Chandler, D. Nature 1996, 379, 55-57. (53) Luzar, A.; Chandler, D. Phys. Rev. Lett. 1996, 76, 928-931. (54) Markovitch, O.; Agmon, N. J. Chem. Phys. 2008, 129, 084505. (55) Eigen, M.; de Maeyer, L. Proc. R. Soc. London, Ser. A 1958, 247, 505-533. (56) Eigen, M. Angew. Chem. Int. Ed. 1964, 3, 1-19. (57) Agmon, N. J. Chem. Phys. 1984, 81, 2811-2817. (58) Agmon, N.; Pines, E.; Huppert, D. J. Chem. Phys. 1988, 88, 5631-5638. (59) Agmon, N.; Weiss, G. H. J. Chem. Phys. 1989, 91, 6937-6942. (60) Gopich, I. V.; Solntsev, K. M.; Agmon, N. J. Chem. Phys. 1999, 110, 2164-2174. (61) Agmon, N.; Szabo, A. J. Chem. Phys. 1990, 92, 5270-5284. (62) Krissinel', E. B.; Agmon, N. J. Comput. Chem. 1996, 17, 1085-1098.

ACS Paragon Plus Environment

20

Page 21 of 28

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. Illustrations of a distorted Eigen-like complex (a) and a Zundel-like complex (b). In the MS-EVB model, the adiabatic state of the system is represented by the linear combination of the EVB states,9 defined by the identification of the imH+ in different EVB states. A three state partition, |1>, |2>, and |3>, and a two state partition, |1> and |2>, are depicted in (a) and (b), respectively. For PT in liquid imidazole, the main EVB state makes dominant contribution for most of the time,9 and the imH+ can be viewed classically as the proton tightly binded to the central imidazole in (a), with the tighter ligand in the first solvation shell denoted by 1x and the looser ligand denoted by 1y. The blue balls and cyan balls represent nitrogen atoms and carbon atoms. The red ball and orange ball denote the center of excess charge (CEC) and the center of charge (COC), respectively.

ACS Paragon Plus Environment

21

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 22 of 28

Figure 2. (a) Two-dimentional probability distribution function P(δ,rNN), in which δ is the reaction coordinate (see text); (b) The nitrogen-nitrogen RDF g(rNN) for NHB on imH+ and NNB on the neighboring imidazole in the first solvation shell (black line). The red line and blue line represent the contribution from the tighter ligand imidazole (1x) and looser ligand imidazole (1y), respectively, as depicted in Figure 1(a); (c) The distribution of P(δ); and (d) free energy profile of PT along the reaction coordinate δ.

ACS Paragon Plus Environment

22

Page 23 of 28

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. The two-dimensional (rNN, ∠NNH) distribution between the imidazolium (imH+) and the imidazole in the first solvation shell (a) as well as that between two imidazole molecules (b).

ACS Paragon Plus Environment

23

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 24 of 28

Figure 4. HB correlation functions for HBs between imH+ and the imidazoles in the first solvation shell C01H(t) (black line), HBs between the first and second solvation shell of imH+ C12H(t) (red line), and HBs beyond the second solvation shell of imH+ C>2H(t) (blue line), respectively.

ACS Paragon Plus Environment

24

Page 25 of 28

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. (a) The change of the identity of the imH+ cation. The time segment is randomly chosen from the MD trajectory. In the inset, one of the recurrent PT events is amplified, and the dashed lines mark the time length of 0.2 ps. (b) Time evolution of the cluster size centered on imH+ maintained by the HBs (black line). In the inset, the proton transfers from im1 to im2. The red line and blue line represent the cluster size along and against the PT direction, respectively. The vertical black dashed lines mark 0.2 ps, the characteristic time for the librational motion of the HB between the imidazole molecules, before and after the PT events.

ACS Paragon Plus Environment

25

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 26 of 28

Figure 6. The nitrogen-nitrogen time dependent RDFs g(rNN)’s, with different time segments of 00.2 ps, 0.2-0.4 ps, and 0.8-1.0 ps, respectively, before and after PT. (a) Illustration of the structure of solvation structure around the imH+. The proton hops from N2 to N3, and the imH+ switches from im1 to im2. The g(rNN)’s are shown in (b) and (c) for g(rNN)’s between N2 and N3 (centered on N2), in (d) and (e) for g(rNN)’s against the PT direction (centered on N1), and in (f) and (g) for g(rNN)’s along the PT direction (centered on N4).

ACS Paragon Plus Environment

26

Page 27 of 28

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 7. (a) Continuous correlation function Cc(t) (black line). The blue dashed line shows the fitting of Cc(t), with single exponential decay. The inset shows the cut-off radius rcut of 2.9 Å as the threshold of PT. When the distance between rCEC (green ball) and the rCI (red ball) of the pivot imidazole becomes larger than rcut, the Boolean functions shift to zero. (b) The intermittent correlation function Ci(t) up to 5 ns (black line). The blue dashed line is the fitting of Ci(t) with the biexponential decay, and the red dashed line is the diffusion model calculated by the SSDP software.62 The long-time region of the red dashed line shows the power-law decay of t-3/2 behavior. For slow reactions such as PT in liquid imidazole as comparing to PT in water,32 it is more difficult for Ci(t) to reach the asymptotic regime.

ACS Paragon Plus Environment

27

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 28 of 28

Table of Contents Graphic

ACS Paragon Plus Environment

28