Dynamics of proton transfer in imidazole hydrogen-bond chains

Abstract

The dynamics and mechanism of proton transfer in the imidazole (Im) hydrogen-bond (H-bond) chain were studied using a unit cell of the Im crystal structure (H⁺(Im)_n, n = 2–4) as a model system and B3LYP/TZVP calculations and Born–Oppenheimer molecular dynamics (BOMD) simulations as model calculations. The B3LYP/TZVP results suggested that only linear H-bond structures are involved in proton transfer and H⁺(Im)₂ is the smallest, most active Zundel-like intermediate complex, which is preferentially formed in a low local dielectric environment. The potential energy curves for proton displacement confirmed the Eigen–Zundel–Eigen scenario that consists of breaking and forming H-bonds in the Im H-bond chain, and because the energy barrier for the reorientation of Im molecule is high, the ring-flip process is ruled out from the proton-transfer mechanism. The BOMD results over the temperature range of 298 to 500 K confirmed that proton transfer in the Im H-bond chain is a local (short-range) process by showing that the activation energies for proton displacement in H⁺(Im)₂ and H⁺(Im)₄ are nearly the same and comparable to experimental and theoretical values. The proton transfer profiles and vibrational spectra suggested that at low temperatures, the N–N vibration, transferring proton and librational motion in the protonated H-bond are synchronized (coherent), resulting in effective structural diffusion process. The ¹H NMR results confirmed these findings and further revealed that the dynamics of proton transfer at low and high temperatures are different due to the interferences of vibrational and librational motions and increases in the oscillatory shuttling motion at elevated temperatures. These theoretical results lead to the conclusion that the rate-determining process of proton transfer in the Im system is the oscillatory shuttling motion in the Zundel-like intermediate complex and does not necessarily involve reorientation of Im molecule as a key process.

---

Introduction

The dynamics and mechanism of proton transfer in the condensed phase are important for understanding the mechanisms of charge transfer in biological systems¹ and for the development of electrolyte membranes² in electrochemical devices such as fuel cells.³ In aqueous solutions, theoretical and experimental information has been accumulated, from which mechanisms of proton transfer have been proposed.^4–6 In addition to the diffusion of hydrated protons through the “vehicle mechanism” by which protons bound to water molecules travel through the electrolyte via mutual diffusion,⁷ Eigen and de Maeyer⁸ suggested that “structural diffusion”, also known as the “Grotthuss mechanism”,⁹ is responsible for the fast proton transfer in aqueous solutions; the structural diffusion is characterized by proton shuttling in the hydrogen-bond (H-bond) network of water before the actual transfer.

For water-based fuel cells, due to the limitation of the operating temperature which has to be above the boiling point of water, attempts have been made to replace Nafion®, a perfluorinated polymer with trifluoromethane sulfonic (triflic) acid functional groups, with water-free proton-conducting materials.¹⁰ Possessing high proton conductivity and thermal and mechanical stabilities, aromatic heterocyclic compounds have received considerable interest as effective proton solvents,² among which imidazole (Im) has been frequently selected as a representative model.^11–15 Liquid and solid Im have received special attention because of their ability to transport protons at temperatures above the boiling point of water; the boiling and melting temperatures of Im are 529 and 363 K, respectively.¹⁶ Consisting of two N atoms that can act simultaneously as a proton donor and acceptor, Im possesses high amphotericity and the ability to form extensive H-bond chains in the liquid and solid states.¹⁷ Due to the presence of the aromatic heterocyclic ring, the mechanism of proton transfer in the Im H-bond chain is expected to be more complex than in water. Experimental and theoretical studies have anticipated that proton transfer in the Im H-bond chain consists of a Grotthuss-like mechanism and reorientation of the Im ring;^11–16 proton transfer in liquid and solid Im was anticipated to involve rotation of the Im ring about one of the two possible axes.¹ However, because the reorientation time scale of the Im H-bond chain is beyond the ¹⁵N NMR detectable range,¹⁵ there has been controversy concerning whether the rotation of the Im ring (ring-flip mechanism) is indeed the rate-determining process. Outstanding results pertaining to the dynamic processes in liquid and solid Im are summarized as follows.

To study the dynamics of protons in the solid state, proton conductivity in Im single crystals was measured as a function of temperature.¹⁶ Because the unit cell at room temperature is monoclinic with four Im molecules forming a linear H-bond chain in the c crystallographic axis,¹⁸ the activation energy obtained from conductivity measurements (approximately 164 kJ mol⁻¹ or 1.7 eV) is considered to be too high for the Grotthuss-like mechanism (breaking and forming of the Im H-bond chain) to be the rate-determining step. In contrast, the activation energies associated with the lifetime of the H-bond network formed from Im tethered to oligomeric aliphatic backbones in glassy and crystalline phases were predicted by constant-pressure molecular dynamics (NPT-MD) simulations to be only 9.3 and 13.3 kJ mol⁻¹, respectively, whereas the activation energy for the Im liquid under neutral conditions is 8.7 kJ mol⁻¹.¹⁹ Although the dynamics of Im molecules in the linear H-bond chain are not well understood, molecular reorientation was anticipated in ref. 16 to be the rate-determining step of proton transfer in the Im crystalline solid. The ring-reorientation mechanism was concluded based on solid-state ¹⁵N NMR experiment conducted on imidazolium methyl sulfonate to be the rate-determining step of proton transfer in Im-based materials with an activation energy of 63.7 kJ mol⁻¹.¹¹

The pronounced anisotropy of proton conductivity in the Im crystalline solid was studied using solid-state ¹⁵N NMR spectroscopy.¹² In contrast to the mechanism proposed in ref. 16, two-dimensional ¹⁵N exchange results revealed that the proton-transfer mechanism is not governed by the reorientation of the Im ring and that the charge transfer along the c crystallographic direction takes place through the tunneling of proton in the linear H-bond structure. These findings are in accord with MD simulations performed using a reactive MS-EVB (multistate empirical valence bond) model on an excess proton solvated in the Im liquid, suggesting that the proton is tightly associated with Im and that proton transfer in the Im H-bond chain is mediated by a Zundel-like complex (Im⋯H⁺⋯Im) through the Eigen–Zundel–Eigen scenario with a free energy barrier of approximately 28 kJ mol⁻¹.¹⁴ Based on the observation that proton transfer in the Im H-bond chain involves approximately four Im molecules (a unit cell of the Im crystal structure), the dynamics of the proton-transfer process were concluded to be a local event with very short spatial/temporal correlation described by a few shuttling or recurrent proton-transfer events.¹⁴

In our previous work,²⁰ proton dissociation and transfer in a phosphoric acid (H₃PO₄)-doped Im system were studied using a quantum chemical method at the B3LYP/TZVP level and Born–Oppenheimer molecular dynamics (BOMD) simulations. This work focused on the dynamics of protons in the intermediate complexes and the effects of fluctuations in the local dielectric environment and H-bond chain length. The theoretical results showed that the intermediate complexes are preferentially formed in short H-bond chains in a low local dielectric environment and a high local dielectric environment is required to stabilize the positive charge and prevent protons from undergoing recurrent transfer. The BOMD results obtained over the temperature range of 298 to 423 K confirmed that the fluctuations in the local dielectric environment and H-bond chain length must be included in the proton-dissociation and proton-transfer mechanisms. Most importantly, a helical-rotational motion was suggested to drive protons away from the H₃PO₄ dopant.

Although fundamental information regarding the H₃PO₄-doped Im system has been obtained,²⁰ due to smaller electrostatic effects, the dynamics and mechanism of proton transfer in the undoped Im H-bond chain are expected to be different. To study proton transfer in heterocyclic aromatic compounds and to prove whether the ring-flip process is the rate-determining step, proton transfer in the Im H-bond chain was investigated in this work. The theoretical methods applied successfully in our previous work²⁰ were used with an emphasis on the dynamics in the protonated H-bond and the reorientational motions (ring-flip and helical-rotational motions) in the Im H-bond chain. Because some spectroscopic results on protonated aromatic heterocyclic systems have been documented, characteristic vibrational and ¹H NMR spectra associated with the motions in the Im H-bond chain are herein discussed in detail.

Theoretical methods

Quantum chemical calculations

Although the DFT method at the B3LYP/TZVP level has been tested and extensively used in our previous study,²⁰ to ensure reliable results, they were examined again in this work using static and spectroscopic properties obtained from RIMP2/TZVP calculations²¹ as benchmarks. Three basic steps applied successfully in our previous work^20,22,23 were employed in this work. All of the quantum chemical calculations and BOMD simulations were carried out using the TURBOMOLE 7.02 software package.²⁴

Presolvation model

Our experience^20,22,25 has shown that proton transfer in protonated H-bond systems can be studied reasonably well based on the concept of presolvation,²⁶ in which the weakening or breaking of some H-bonds in the first and second solvation shells is required to form intermediate complexes. Because a unit cell of the Im crystal structure consists of four Im molecules forming a linear H-bond chain in the c crystallographic direction¹⁸ and excess proton conditions are required to promote proton transfer,²⁰ H⁺(Im)_n (n = 2–4) were selected as model systems. The applicability of these model systems is justified by the observation that proton transfer in the Im H-bond chain is a local process, as discussed in detail in ref. 20. In addition, because proton transfer in protonated H-bond systems is sensitive to the local dielectric environment,^20,22,25 B3LYP/TZVP calculations and BOMD simulations were performed both in the gas phase (ε = 1) and in continuum solvent. In this work, a conductor-like screening model (COSMO) with the dielectric constant of the Im liquid (ε = 23)²⁷ was included in the static and dynamic calculations.

Static calculations

Structures and energetics. To provide fundamental information for the discussion of the dynamic results, the structural, energetic and spectroscopic properties of H⁺(Im)_n (n = 2–4) were computed using B3LYP/TZVP calculations. The likelihood of proton transfer in the Im H-bond chain was anticipated based on the asymmetric stretching coordinate (Δd_DA), computed from Δd_DA = |d_N–H − d_N⋯H| (Fig. 1a);²⁰ the intermediate complex of proton transfer is generally characterized by a Δd_DA value that is close to zero.²² The relative orientation of the Im rings in the Im H-bond chain was characterized by the dihedral angles (ω); for example, in Fig. 1b, the dihedral angle of the two Im rings in H⁺(Im)₂ is ω₁ = ∠C₁–N₅–N₁₄–C₁₁. In addition, to visualize the relative orientation of the Im rings, reference vectors (V_n) coinciding with the C₂ axes of the H⁺(Im) molecules were defined, and projections of these vectors similar to the Newman projection were used to characterize the ring orientations, e.g., projection of V₁ and V₂ in Fig. 2a. The line of sight for the projections in Fig. 2 was defined lengthwise down the Im H-bond chain from the right-hand side to the left-hand side.

Fig. 1 (a) Velocity vectors (v^NH,MD₁, v^NH,MD₂ and v^NN,MD₁₂) used in the calculations of vibrational spectra of N–H⋯N H-bond; Q₁ = symmetric N–H stretch; Q₂ = asymmetric N–H stretch and; Q₃ = N–N vibration. (b) Intermolecular velocity vectors (v₁ and v₂) used in the study of torsional motion of N–H⋯N H-bond (ω₁ = ∠C₁–N₅–N₁₄–C₁₁); Q₄ = torsional rotation.

Fig. 2 (a and b) Reference vectors (V₁–V₃) and Newman projections used to describe relative orientations of Im molecules in H⁺(Im)₂ and H⁺(Im)₃, respectively. θ = rotational angle of V₂ reference vector used in the discussion of the ring-flip mechanism.

In this work, the interaction (ΔE) and solvation (ΔE^sol) energies of H⁺(Im)_n (n = 2–4) were computed using B3LYP/TZVP calculations,²⁰ and the potential energy curves for proton displacement and for the reorientation (ω) of the Im molecule in the Im H-bond chain were constructed using a freeze-scan method. They were represented by plots of the relative energy (E^rel) as a function of Δd_DA and ω (Fig. 1b), respectively.

Vibrational and NMR spectra. The vibrational frequencies of H-bonds in H⁺(Im)_n (n = 2–4) were computed based on the harmonic approximation with a scaling factor of 0.9614; this scaling factor has been proved to be appropriate for the B3LYP/TZVP method.²⁸ The asymmetric N–H stretching mode (Q₂) defined in Fig. 1a was used in the discussion of the dynamic behavior of the transferring proton.²⁰ As an effective probe for H-bond formation,²⁹ the isotropic shielding constants (σ_H⁺) were computed at the B3LYP/TZVP level using the gauge including atomic orbital (GIAO) method.^30,31 The ¹H NMR chemical shifts of the H-bond protons (δ_H⁺) were computed from σ_H⁺ with respect to tetramethylsilane (TMS), which was determined at the B3LYP/6-311+G** level to be 31.97 ppm.³² Based on the conclusion that σ_H⁺ of strong, protonated H-bonds depend only on the H-bond distance and not the neighboring solvent molecules,³³ and the observation that the effects of continuum solvent are small,^33–35 only σ_H⁺ of the intermediate complexes in the gas phase were computed and discussed.²²

Dynamic calculations

BOMD simulations. In this work, NVT–BOMD simulations with the Nosé–Hoover thermostat^36–38 were performed at the B3LYP/TZVP level using the equilibrium structures of H⁺(Im)_n (n = 2–4) as starting configurations. Two thousand steps of 1 fs each were employed in the equilibration and 10 000 steps were used in the property calculations, corresponding to 2 and 10 ps, respectively.^20,39 The NVT-BOMD simulation conditions with the Nosé–Hoover thermostat were tested and used successfully in our previous study on the H₃PO₄-doped Im system.²⁰

Vibrational and NMR spectra. Because vibrational frequencies obtained from static calculations with the harmonic approximation do not account for coupled motions and interferences, the interplay between vibrational and rotational motions and their effects on proton displacement were studied using the vibrational and rotational spectra obtained from Fourier transformations of velocity autocorrelation functions (VACFs). The velocity vectors associated with the vibrational and rotational motions of interest are defined in Fig. 1a and b and were used in the calculations of the VACFs.⁴⁰ In this work, the asymmetric N–H stretching (ν^asym-NH,MD) and N–N vibrational (ν^NN,MD) frequencies, as well as the torsional rotational frequencies (ν^ω,MD), were computed.²⁰ The activation energy (ΔE^‡,Arr) and the first-order rate constant (k) for proton exchange in the Im H-bond chain were determined by performing NVT-BOMD simulations over the temperature range of 298 to 500 K;²⁰ ΔE^‡,Arr was approximated from the Arrhenius plot, namely the slope of the linear relationship between ln(k) and 1000/T.

The dynamics and energetics of proton transfer in the Im H-bond chain were also discussed based on the ¹H NMR chemical shift spectra . In the calculations of , the BOMD results for the smallest, most active intermediate complex over the temperature range of 298 to 500 K were employed. Statistical samplings of the intermediate structures were performed every five BOMD steps and approximately 2500 intermediate structures were used in the GIAO calculations of the instantaneous ¹H NMR shielding constants .⁴¹ The ¹H NMR chemical shift spectra were represented by Lorentzian peak functions,⁴² from which the activation energy (ΔE^‡,NMR) for the exchange between the shared-proton and close-contact structures was approximated using a line-width analysis. Based on the assumption that the change in the ¹H NMR line width as a function of temperature is correlated with the proton-exchange rate, ΔE^‡,NMR was computed through the Arrhenius equation;⁴³ was approximated as the full width at half maximum (FWHM) of the ¹H NMR chemical shift spectra.⁴⁴ In this study, the effective transverse relaxation times were calculated from using and ΔE^‡,NMR from the plot of and 1000/T.⁴³

Results and discussion

Quantum chemical methods

The equilibrium structures of the intermediate complexes and the corresponding static properties obtained from B3LYP/TZVP and RIMP2/TZVP calculations in the gas phase (ε = 1) and in continuum solvent (ε = 23) are listed in Tables S1 and S2,† respectively. The correlations between the properties obtained from these two methods are illustrated in Fig. S1.† The figure shows that both B3LYP/TZVP and RIMP2/TZVP methods yield nearly the same equilibrium structures, and excellent correlations are achieved for all of the static properties both in the gas phase and in continuum solvent, with correlation coefficients (R²) of 0.98–0.99; RIMP2/TZVP calculations yield slightly stronger and shorter H-bonds due to better representation of the effects of the electron correlations. These findings confirm that for the Im H-bond system, B3LYP/TZVP and RIMP2/TZVP calculations yield similar results and the B3LYP/TZVP method can be applied with confidence.

Static results

The values of Δd_DA obtained from B3LYP/TZVP and RIMP2/TZVP calculations confirm that only the linear H-bond structures shown in Tables S1 and S2† are involved in proton transfer in the Im H-bond system. To facilitate the discussion, the Im H-bond structures are labeled with a four-character code, e.g., Gn-[m]^(ω_i=xx) and Cn-[m]^(ω_i=xx); G = equilibrium structure in the gas phase; C = equilibrium structure in continuum solvent; n = number of Im molecule and; ω_i = dihedral angle of H-bond (i). Different H-bond structures with the same number of Im molecules are distinguished by [m]. For example, G2-[1]^{(ω₁=−94)} and G2-[2]^(ω₁=180) are H⁺(Im)₂ (n = 2) with different relative ring orientations, ω₁ = −94 and 180 degrees, respectively. In addition, Im molecules in H⁺(Im)_n are labeled with Roman numerals in parentheses from (I) to (IV).

Intermediate complexes. B3LYP/TZVP calculations predict that in the gas phase, a perpendicular structure (structure G2-[1]^{(ω₁=−94)}) is the absolute minimum energy geometry of H⁺(Im)₂ with ΔE = −116.7 kJ mol⁻¹, R_N–N = 2.67 and Δd_DA = 0.45 Å, whereas a planar structure (structure G2-[2]^(ω₁=180)) is approximately 5.6 kJ mol⁻¹ less stable with R_N–N = 2.69 and Δd_DA = 0.49 Å. The trends of ΔE, R_N–N and Δd_DA of the perpendicular and planar structures (Zundel-like intermediate complexes) are in accord with the other static properties, e.g., σ_H⁺ and ν^NH are smaller for structure G2-[1]^{(ω₁=−94)}.

The potential energy curves shown in Fig. S2a† reveal that the energy barriers for proton displacement in structures G2-[1]^{(ω₁=−94)} and G2-[2]^(ω₁=180) are nearly the same within the thermal energy fluctuations at room temperature (approximately 9.5 and 11.0 kJ mol⁻¹, respectively), with the rotational energy barrier of H-bond (1) in Fig. S2b† approximately equal to 5 kJ mol⁻¹. The former are naturally higher than that obtained from a relaxed-scan method in ref. 45 (approximately 5 kJ mol⁻¹ at the B3LYP/6-311G** level), and the latter suggests almost barrier-less rotation of the N–H⁺⋯N H-bond. Structure G2-[1]^{(ω₁=−94)} is characterized by ν^NH = 1952 cm⁻¹ and δ_H⁺ = 21.7 ppm, whereas for structure G2-[2]^(ω₁=180), ν^NH = 2020 cm⁻¹ and δ_H⁺ = 21.2 ppm. The vibrational frequencies are slightly lower than the ν^NH value of protons in the asymmetric low-barrier H-bond (LBHB) of H⁺(Dih)₂ (Dih = 4,5-dihydro-1H-imidazole or C₃H₆N₂), calculated based on the B3LYP/D95+(d,p) method with the harmonic approximation equal to 2172 cm⁻¹;⁴⁶ moreover, δ_H⁺ are slightly higher than the resonance of the mobile proton in the asymmetric N–H⁺⋯N H-bond in H⁺(Pyd)₂ (Pyd = pyridine), determined experimentally at 250 K to be 20.2 ppm.⁴⁷ The spectroscopic results indicate a stronger protonated H-bond in H⁺(Im)₂ than in H⁺(Dih)₂ and H⁺(Pyd)₂.

The equilibrium structure of H⁺(Im)₃ is represented by an Eigen-like structure (structure G3-[1]^{(ω₁=−93)} in Table S1†), which is characterized by perpendicular arrangements of all of the Im molecules. It appears that as the number of Im molecule increases from two to three, the N–H⁺⋯N H-bond is destabilized and the proton in H-bond (1) is transferred from Im (I) to Im (II). The destabilization effect is manifested by increases in the static properties of the protonated H-bonds in structure G3-[1]^{(ω₁=−93)}, e.g., R_N–N = 2.74, Δd_DA = 0.61 Å, and ν^NH = 2421 cm⁻¹. In addition, to study the energetics of the “ring-flip” process,¹⁵ the potential energy curve for the rotation (θ) of the V₂ reference vector in structure G3-[1]^{(ω₁=−93)} (Fig. 2b) was constructed using a freeze-scan method and is shown in Fig. S2c.† The B3LYP/TZVP results suggest that the highest energy barrier at θ = 90 degrees is approximately 147 kJ mol⁻¹, which is higher than that of the ring-flip process in neutral Im H-bond chains, determined from B3LYP/6-311G(d,p) calculations to be 123 kJ mol⁻¹.¹⁹ These relatively high energy barriers rule out the likelihood that the ring-flip process is directly involved in the proton-transfer mechanism; the ring-flip process requires high energy to simultaneously break two protonated H-bonds in the Im H-bond chain.

The B3LYP/TZVP results in Table S1† and the potential energy curves in Fig. S2a and d† show nearly the same static properties of the protonated H-bonds in H⁺(Im)₂ and H⁺(Im)₄, which suggest similar behavior of proton displacement in H-bond (1) of structure G2-[1]^{(ω₁=−94)} and H-bond (2) of structure G4-[1]^{(ω₁=−92)}. Structure G4-[1]^{(ω₁=−92)} is represented by perpendicular arrangements of all of the Im molecules with the energy barrier for proton displacement in H-bond (2) slightly higher than that of H-bond (1) in structure G2-[1]^{(ω₁=−94)}, approximately 12 kJ mol⁻¹ for the perpendicular arrangement of the Im (II) and Im (III) molecules. These results suggest that structure G2-[1]^{(ω₁=−94)} is the smallest, most active intermediate complex that can be used as a representative model in the dynamic calculations.

Effects of local dielectric environment. Table S1† suggests that although the Im H-bond structures in low and high local dielectric environments are not substantially different, the protonated H-bonds are destabilized in a high local dielectric environment with protons staying closer to the proton donors; for example, for structure C2-[1]^{(ω₁=−94)}, ΔE = −32.6 kJ mol⁻¹, R_N–N = 2.74 and Δd_DA = 0.60 Å, whereas for structure G2-[1]^{(ω₁=−94)}, ΔE = −116.7 kJ mol⁻¹, R_N–N = 2.67 and Δd_DA = 0.45 Å. The destabilization of the protonated H-bonds in COSMO with ε = 23 is also evident from the blue shifts of ν^NH, which confirm that the shared-proton structures (Zundel-like intermediate complexes) are not preferentially formed in a high local dielectric environment.^20,41

Mechanism of proton transfer. Based on the static results, possible proton-transfer paths are proposed in Fig. 3, in which all of the intermediate complexes and the corresponding potential energy curves for proton displacement, as well as the effects of the Im H-bond chain extension, are included; paths (a) and (b) are under neutral and excess proton conditions in a low local dielectric environment, respectively, whereas (c) shows the effects of a high local dielectric environment on the potential energy curves. Because the shapes and the well depths of the potential energy curves on path (a) are characteristics of ordinary H-bonds (asymmetric single-well potential with the energy barriers of approximately 131 kJ mol⁻¹) and the extension of the Im H-bond chain from n = 2 to 4 does not lead to significant changes in the shapes and the well depths, one can rule out the possibility of the formation of the Zundel-like intermediate complex under neutral conditions.

Fig. 3 Proposed mechanisms of proton transfer which involve fluctuations of Im H-bond chain length and local dielectric-environment. The intermediate complexes are illustrated with the potential energy curves for single-proton displacement. (a and b) Under neutral and excess proton conditions (ε = 1), respectively. (c) Under excess proton conditions (ε = 23). E^rel = relative energy; Δd_DA = asymmetric stretching coordinate.

In contrast, under excess proton conditions on path (b), starting from the Zundel-like intermediate complex in a low local dielectric environment (structure G2-[1]^{(ω₁=−94)}), an extension of the Im H-bond chain at Im (II) leads to an Eigen-like complex (structure G3-[1]^{(ω₁=−93)}), in which the proton and the minimum of the potential energy curve of H-bond (1) move from Im (I) to Im (II). Further extension of the Im H-bond chain at Im (III) leads to structure G4-[1]^{(ω₁=−92)} with an increase in the energy barrier of H-bond (1) (R_N–N = 2.77 and Δd_DA = 0.65 Å) and a decrease in the energy barrier of H-bond (2) (R_N–N = 2.68 and Δd_DA = 0.48 Å). Because the potential energy curve of H-bond (2) is nearly identical to that of H-bond (1) in structure G2-[1]^{(ω₁=−94)}, the protonated H-bond formed between Im (II) and Im (III) (H-bond (2)) becomes a new precursor for the next transfer. This finding confirms the Eigen–Zundel–Eigen scenario proposed in ref. 14, in which H-bond breaking and forming were suggested to represent the key processes in proton transfer in the Im H-bond chain.

The role played by the local dielectric environment can also be visualized in Fig. 3. It appears in (c) that although the shapes of the double-well potentials in low and high local dielectric environments are similar, the energy barriers are higher in high local dielectric environment. These findings confirm that high local dielectric environment partially neutralizes the protonic charge in the H-bond, which leads in this case to increases in the strength of the N–H⁺ covalent bond and the energy barrier that prevents protons from returning to the original Im molecules.²⁰

Dynamic results

The BOMD results on H⁺(Im)_n (n = 2–4) with ε = 1 and 23 confirm that only the linear H-bond structures are involved in proton transfer and the Zundel-like intermediate complexes are preferentially formed in low local dielectric environment; the shared-proton structures are formed only in the Im H-bond chains with n = 2 and 4, in accord with the static results and the H₃PO₄-doped Im system.²⁰ Therefore, only the dynamics of protons in H⁺(Im)₂ and H⁺(Im)₄ in low local dielectric environment (ε = 1) are discussed in detail.

Characteristic vibrations in protonated Im H-bond. To study the dynamics of the Eigen–Zundel–Eigen mechanism,¹⁴ the BOMD results on structure G2-[1]^{(ω₁=−94)} near the melting temperature are discussed first. The proton transfer profiles at 350 K in Fig. 4 show two characteristic N–N vibrations, namely, large- and small-amplitude vibrations, labeled L and S in Fig. 4b and c. The average R_N–N for the large-amplitude vibration, for example, in panel P1, is 2.67 ± 0.10 Å, whereas that of the small-amplitude vibration in panel P2 is 2.65 ± 0.05 Å. The variations of R_N–H with respect to the simulation time in Fig. 4 also suggest two outstanding N–H stretching modes associated with the large- and small-amplitude vibrations, namely, the structural diffusion and oscillatory shuttling motions,⁴⁸ respectively; the structural diffusion motion is driven by the large-amplitude vibration, for which R_N–N varies over a wide range, and at the shortest distance, the proton oscillates and transfers, whereas for the small-amplitude vibration, R_N–N varies within a narrow range and the proton only oscillates in the middle of the H-bond with no actual transfer.

Fig. 4 (a–c) Variations of R_N–N and R_N–H distances and time evolutions of dihedral angle (ω₁) of H-bond (1) in structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations at 350 K. L and S = large- and small-amplitude N–N vibrations, respectively.

The vibrational spectra of H-bond (1) in structure G2-[1]^{(ω₁=−94)} obtained from Fourier transformation of the VACF at 350 K are shown in Fig. 5a. Analysis of the proton transfer profiles in Fig. 4 suggests that the N–N vibrational spectra can be divided into two categories, namely, the primary frequency associated with characteristic N–N vibrations ranging between approximately 0 and 200 cm⁻¹ and secondary frequencies resulting from interferences of other motions in the H-bond (>500 cm⁻¹). At 350 K, the primary frequency is located at the lowest wavenumber, ν^NN,MD_P = 186 cm⁻¹, with three outstanding secondary peaks at 955, 1095 and 1537 cm⁻¹, respectively; the secondary peaks found for the protonated H-bond in poly(vinyl imidazole) (PVIm) are located at 943, 1144 and 1545 cm⁻¹.

Fig. 5 (a) Vibrational spectra of H-bond (1) in structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations at 350 K. (b) Arrhenius plot for proton transfer in H-bond (1) of structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations over the temperature range of 298–500 K. (c) Vibrational spectra for torsional oscillation (ω₁) of H-bond (1) in structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations at 350 K. A and B = oscillatory shuttling and structural diffusion peaks, respectively; C = characteristic peak of N–H stretching mode in the imidazolium cation (H⁺(Im)); P = primary vibrational peak; ΔE^‡,Arr = activation energy obtained from the slope of the linear relationship between ln(k) and 1000/T.

For the asymmetric N–H stretching spectra obtained at 350 K, three characteristic peaks are located at ν^asym-NH,MD_A = 824, ν^asym-NH,MD_B = 1902 and ν^asym-NH,MD_C = 2457 cm⁻¹. Peaks A and B are associated with protons undergoing the oscillatory shuttling and structural diffusion motions, respectively, whereas peak C is characteristic of the N–H stretching in the imidazolium cation (H⁺(Im)); for H⁺(Dih)₂, peaks A and B were reported based on explicit quantum calculations on protons and BOMD simulations at 300 K to occur at 938 and 1828 cm⁻¹, respectively,⁴⁶ and peak C was reported based on an IR experiment on acid-doped poly(benzimidazole) to be in the range between 2550 and 2950 cm⁻¹.⁴⁹ Because actual proton transfer occurs nearly instantaneously through structural diffusion motion, the oscillatory shuttling motion, which possesses a longer lifetime, is considered to be the rate-determining process.

Fig. 5b shows that for proton transfer in structure G2-[1]^{(ω₁=−94)}, the expected Arrhenius behavior is established only in low-temperature range (298–350 K), whereas in high-temperature range (380–500 K), the plot of ln(k) versus 1000/T does not show a good linear relationship. These findings indicate complex dynamics in the protonated H-bond in the high-temperature range. The slope of the Arrhenius plot yields the activation energy (ΔE^‡,Arr) for proton exchange in structure G2-[1]^{(ω₁=−94)} of 6.2 kJ mol⁻¹, which is slightly lower than the activation energies associated with the H-bond lifetime in liquid Im under neutral conditions (360–550 K)¹⁹ and the MS-EVB result under excess proton conditions (373–413 K) of 8.7 kJ mol⁻¹.⁴⁵

Interference of librational motions. To discuss the dynamics in the low- and high-temperature ranges, the proton transfer profiles and the time evolutions of the dihedral angle (ω₁) in Fig. 4 and the relaxation behaviors of the envelope of the VACF of the N–N vibration in Fig. 6 are considered. The VACF plots over the temperature range of 298 to 500 K suggest that the deviation of the Arrhenius plot at elevated temperatures results from interferences of motions with the same or similar vibrational frequencies as the primary N–N vibration; at 350 K, interferences appear in the VACF plot (Fig. 6b) at approximately 0.17 ps and become more pronounced at higher temperatures (e.g., Fig. 6c). To study the dynamics of librational motions in H-bond (1) of structure G2-[1]^{(ω₁=−94)} and their interference effects, vibrational spectra of the torsional oscillation (ω₁) in Fig. 5c are discussed, in which the primary librational peak is broad, ranging from 48 to 166 cm⁻¹, with three outstanding secondary peaks at 1115, 912 and 1421 cm⁻¹, respectively. Because the primary frequencies of the N–N vibration and torsional oscillation are similar, they are hypothesized to partially correlate or mutually interfere at this temperature; for example, in panel P1 in Fig. 4b, the primary librational motion (ω₁) is partially correlated with the primary N–N vibration (with the longest wavelength).

Fig. 6 (a–c) VACF plots of N–N vibration in structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations at 298, 350 and 380 K, respectively.

Because the VACF plot in Fig. 6a shows almost no interference, to prove this hypothesis, the proton transfer profiles, the time evolutions of ω₁ and the vibrational spectra at 298 K are discussed in detail. Fig. 7 suggests that due to moderate thermal energy fluctuations, the small-amplitude vibration (S) does not lead to oscillatory shuttling motion (e.g., panel P1 in Fig. 7b). The proton transfer profile in Fig. 7c suggests that the large-amplitude vibration (L) and the primary librational motion (ω₁) dominate in panel P2 and proton exchange in H-bond (1) is synchronous with these two characteristic motions. In other words, the primary N–N vibration, the transferring proton and the primary librational motion are synchronized (coherent), resulting in smoother and more effective proton exchange at 298 K (Fig. 6a) than at 350 K (Fig. 6b). The vibrational spectra in Fig. 8 support these interpretations by showing that at 298 K, the primary N–N vibrational peak is more structured than that at 350 K; ν^NN,MD_P = 156 cm⁻¹ with small, broad secondary peaks at 1058, 940 and 1568 cm⁻¹, respectively. It should be noted that the secondary librational peak at 912 cm⁻¹ in the librational spectra at 350 K (Fig. 5c) becomes shoulders at 298 K (Fig. 8b) due to weaker thermal energy fluctuations. Therefore, one can conclude that the secondary librational modes ranging between 500 and 1000 cm⁻¹ interfere with the dynamics of proton exchange in the protonated Im H-bond especially at elevated temperatures.

Fig. 7 (a–c) Variations of R_N–N and R_N–H distances and time evolutions of dihedral angle (ω₁) of H-bond (1) in structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations at 298 K. L and S = large- and small-amplitude N–N vibrations, respectively.

Fig. 8 (a) Vibrational spectra of H-bond (1) in structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations at 298 K. (b) Vibrational spectra for torsional oscillation (ω₁) of H-bond (1) in structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations at 298 K. A and B = oscillatory shuttling and structural diffusion peaks, respectively; C = characteristic peak of N–H stretching mode in the imidazolium cation (H⁺(Im)); P = primary vibrational peak.

Effects of H-bond chain extension. The BOMD results in Fig. S3 and S4† show that all of the characteristic properties of the protonated H-bond in H⁺(Im)₂ are not significantly changed when n is increased from 2 to 4. Because of an extensive H-bond chain, the large-amplitude vibrational motion dominates at all temperatures, e.g., at 350 K in P1 in Fig. S3a and b.† The primary N–N vibrational peak at 350 K (Fig. S4a)† is at ν^NN,MD_P = 163 cm⁻¹ with only one well-defined secondary interference peak at 943 cm⁻¹. Peaks B and C in the asymmetric N–H stretching spectra are red shifted compared with those in H-bond (1) of structure G2-[1]^{(ω₁=−94)}; ν^asym-NH,MD_A = 846, ν^asym-NH,MD_B = 1500, ν^asym-NH,MD_C = 2021 cm⁻¹. The red shifts of peaks B and C imply that the structural diffusion motion in H⁺(Im)₄ is slightly more feasible than in H⁺(Im)₂ at the same temperature. Because the primary librational peak (Fig. S4b)† is observed in the range between 38 and 156 cm⁻¹ with only one comparable secondary peak nearby at 974 cm⁻¹, the librational interference effects in H⁺(Im)₄ are expected to be weaker than those in H⁺(Im)₂. This is also evident from smoother exponential decay functions in the VACF plots over the entire temperature range (Fig. S5†). The plot of ln(k) versus 1000/T in Fig. S4c† shows a good linear relationship over the entire temperature range, confirming that the librational interference effects in H⁺(Im)₄ are not as strong as those in H⁺(Im)₂; the slope of the Arrhenius plot yields ΔE^‡,Arr = 5.8 kJ mol⁻¹, which is slightly lower than that in H⁺(Im)₂.

Because the time evolutions of the torsional angles in H⁺(Im)₄ show correlations only among ω₂, the N–N vibration and the transferring proton, e.g., at 350 K in panels P1 and P3 in Fig. S3 and S6 and at 298 K in panels P1 and P2 in Fig. S7 and S8,† the H-bond chain structure reorganization through the variations in ω₁ and ω₃ (helical-rotational motion) can be ruled out from the proton-exchange mechanism. The situation is different from that in the H₃PO₄-doped Im system, in which the potential energy curves for proton displacement and BOMD simulations revealed that proton transfer from H₄PO₄⁺ along the Im H-bond chain is driven by a helical-rotational motion of the Im H-bond chain.²⁰

The spectroscopic and energetic results obtained for H⁺(Im)₄ justify the conclusion that H⁺(Im)₂ is the smallest, most active intermediate complex and the use of structure G2-[1]^{(ω₁=−94)} as a representative. Based on this conclusion, the mean-square displacements (MSD) of the exchanging proton in structure G2-[1]^{(ω₁=−94)} were computed. The slope of the MSD plot yields the excess proton self-diffusion coefficient D = 1.8 × 10⁻⁵ cm² s⁻¹ at 350 K, which is in excellent agreement with the values obtained from MS-EVB simulations on 216 Im molecules and one excess proton⁴⁵ and conductivity measurements at 393 K of 2.0 × 10⁻⁵ cm² s⁻¹.⁵⁰

Characteristic NMR spectra. The ¹H NMR chemical shift spectra of the transferring proton in structure G2-[1]^{(ω₁=−94)} obtained over the temperature range of 298 to 500 K are included in Fig. 9. The spectra are characterized by three well-separated peaks, which reflect the three asymmetric N–H stretching modes (also labeled A, B and C in Fig. 9). The ¹H NMR results obtained at 298 K, shown in Table S3,† suggest that peak A, which is associated with the shared-proton structure and cannot be obtained from the static calculations, is centered at , whereas peaks B and C are observed at lower ppm, approximately at 21.0 and 17.6 ppm, respectively. Because peak A is located at higher ppm than that of the mobile proton in the asymmetric N–H⁺⋯N H-bond of H⁺(Pyd)₂ (21.7 ppm),⁴⁷ the efficiency of proton transfer is expected to be higher in Im systems.

Fig. 9 (a–c) ¹H NMR chemical shift spectra of transferring proton in H-bond (1) of structure G2-[1]^{(ω₁=−94)} obtained from BOMD simulations at 298, 350 and 450 K, respectively.

The structures and intensities in the ¹H NMR chemical shift spectra of the transferring proton in structure G2-[1]^{(ω₁=−94)} in Fig. 9 show two outstanding patterns, which confirm different leading protonated species in the low- and high-temperature ranges. In the low-temperature range (298–350 K), peaks C and B are comparable and dominate peak A, implying that the structural diffusion motion is more feasible; the imidazolium (H⁺(Im)) cation (peak C) and the close-contact (N–H⁺⋯N) structures (peak B) dominate, in accord with the MS-EVB results reported in ref. 14, in which protons are associated with an Im molecule to form H⁺(Im) and the shared-proton structure is rare. However, over the high-temperature range (400–500 K), peaks A (shared-proton structure) and B (close-contact structure) are comparable and dominate peak C, indicating that the oscillatory shuttling motion is increased.

The activation energies (ΔE^‡,NMR) for proton exchange in the low- and high-temperature ranges are computed using the results of the line-width analysis in Table S3† and the plots of versus 1000/T in Fig. 10. The expected Arrhenius behavior, which reflects motional narrowing of the ¹H NMR peak with respect to temperature, is observed for peak A both in the low- and high-temperature ranges (Fig. 10a) with ΔE^‡,NMR = 7.3 and 12.5 kJ mol⁻¹, respectively. The former is in good agreement with the values obtained from the vibrational analysis (ΔE^‡,Arr = 6.2 kJ mol⁻¹) and the MS-EVB result (8.7 kJ mol⁻¹);⁴⁵ the MS-EVB result was obtained from the linear relationship between ln(D) and 1000/T over the temperature range of 373 to 413 K. ΔE^‡,NMR is higher in the high-temperature range because of stronger thermal energy fluctuations, the increase in the oscillatory shuttling motion and the interferences of vibrational and librational motions as discussed above. For the structural diffusion peak (peak B), the Arrhenius plot in Fig. 10b shows the expected behavior only over the low-temperature range, with an ΔE^‡,NMR value of approximately 2 kJ mol⁻¹, confirming that the oscillatory shuttling motion which possesses higher ΔE^‡,NMR is the rate-determining process.

Fig. 10 Plots of the natural log of effective transverse relaxation time as a function of 1000/T. ΔE^‡,NMR = activation energy obtained from ¹H NMR line shape analysis. (a and b) For peaks A and B, respectively.

Conclusion

In this work, the dynamics and mechanism of proton transfer in Im H-bond chains were studied based on the concept of presolvation and the B3LYP/TZVP method and BOMD simulations as model calculations. The advantages of our theoretical treatments are the use of combined vibrational and NMR spectroscopic techniques which allow local dynamical motions and their activation energies to be studied in details. The B3LYP/TZVP results suggest that the linear H-bond structure with perpendicular arrangements of Im molecules represents the most preferential structure, and in a low local dielectric environment, the shared-proton structure with n = 2 is the smallest, most active intermediate complex. Because the potential energy barrier for the reorientation of Im molecule is significantly higher than that for proton displacement in H-bond, the ring-flip process is ruled out from the proton-transfer mechanism.

Based on the static results, the mechanism of proton transfer in the Im H-bond chain was confirmed to be represented by the Eigen–Zundel–Eigen scenario, which involves breaking and forming local H-bonds; the breaking and forming H-bonds are required for the generation of the Zundel-like (n = 2) and Eigen-like complex (n = 3) formations, respectively, and the fourth Im molecule is required to form the Zundel-like complex for the next transfer. In this work, the potential energy curves for proton displacement also indicated that a high local dielectric environment destabilizes the shared-proton structure and stabilizes the N–H⁺ covalent bond, which leads to an increase in the energy barrier that prevents the proton from returning to the original Im molecule.

The BOMD results confirmed that only the linear H-bond structures are involved in proton transfer and that the smallest, most active intermediate complex are preferentially formed in a low local dielectric environment. An analysis of the proton transfer profiles in H⁺(Im)₂ and H⁺(Im)₄ suggested that the N–N vibrational spectra can be divided into two categories, namely, primary frequencies associated with slow N–N vibrations and secondary frequencies resulting from interferences of faster motions in the protonated H-bond. However, three characteristic peaks associated with protons undergoing oscillatory shuttling and structural diffusion motions and N–H stretching in the imidazolium cation (H⁺(Im)) were observed in the asymmetric N–H stretching spectra. The BOMD results further suggested that at low temperatures, the primary N–N vibration, the transferring proton and the primary librational motion are synchronized (coherent), resulting in smoother and more effective proton exchange than at high temperatures.

Because the time evolutions of torsional angles show no evidence for reorganization in the long H-bond chain structure, the helical-rotational motion is ruled out from the proton-exchange mechanism. Moreover, because the activation energies for proton displacement in H⁺(Im)₂ and H⁺(Im)₄ obtained from vibrational analyses are almost the same, H⁺(Im)₂ can be used as a representative in dynamic calculations. The structures and intensities in the ¹H NMR chemical shift spectra showed two outstanding patterns, which confirm different leading protonated species in the low- and high-temperature ranges; the structural diffusion dominates at low temperatures, and the likelihood of the oscillatory shuttling motion is increased at high temperatures. ¹H NMR analyses confirmed these findings and further suggested that the interference between vibrational and librational motions affects the activation energy. These theoretical results lead to the conclusion that the rate-determining process of proton transfer in the Im system is the oscillatory shuttling motion in the Zundel-like intermediate complex and does not necessarily involve reorientation (ring-flip) of Im molecule.

Acknowledgements

Scholarship for Worapong Bua-ngern is from the Office of the Higher Education Commission (OHEC), Thailand under the program Strategic Scholarships for Frontier Research. Prof. Kritsana Sagarik obtained financial supports from the Thailand Research Fund (TRF), the Advanced Research Scholarship (BRG-5580011). Dr Sermsiri Chaiwongwattana received SUT post-doctoral research fellowship. This work was supported in part by the Higher Education Research Promotion and National Research University Project (HERB-NRU), OHEC, Thailand. High-performance computer facilities provided by the following organizations are gratefully acknowledged; School of Mathematics and School of Chemistry, SUT; National e-Science project of the National Electronics and Computer Technology Center (NECTEC), National Science and Technology Development Agency (NSTDA).

References

1. J. F. Daycock, G. P. Jones, J. R. N. Evans and J. M. Thomas, Nature, 1968, 218, 672.
2. M. F. H. Schuster and W. Meyer, Annu. Rev. Mater. Res., 2003, 33, 233.
3. W. L. Cavalcanti, R. Marschall, P. Tölle, C. Köhler, M. Wark and T. Frauenheim, Fuel Cells, 2008, 8, 244.
4. U. W. Schmitt and G. A. Voth, J. Chem. Phys., 1999, 111, 9361.
5. G. A. Voth, Acc. Chem. Res., 2006, 39, 143.
6. A. Roudgar, S. P. Narasimachary and M. Eikerling, Chem. Phys. Lett., 2008, 457, 337.
7. K.-D. Kreuer, S. J. Paddison, E. Spohr and M. Schuster, Chem. Rev., 2004, 104, 4637.
8. M. Eigen and L. de Maeyer, Proc. R. Soc. London, Ser. A, 1958, 247, 505.
9. C. J. D. von Grotthuss, Ann. Chim., 1806, 58, 54.
10. S. J. Paddison, Annu. Rev. Mater. Res., 2003, 33, 289.
11. I. Fischbach, H. W. Spiess, K. Saalwächter and G. R. Goward, J. Phys. Chem. B, 2004, 108, 18500.
12. B. S. Hickman, M. Mascal, J. J. Titman and I. G. Wood, J. Am. Chem. Soc., 1999, 121, 11486.
13. M. Iannuzzi, J. Chem. Phys., 2006, 124, 204710.
14. A. Li, Z. Cao, Y. Li, T. Yan and P. Shen, J. Phys. Chem. B, 2012, 116, 12793.
15. W. Münch, K. D. Kreuer, W. Silvestri, J. Maier and G. Seifert, Solid State Ionics, 2001, 145, 437.
16. A. Kawada, A. R. McGhie and M. M. Labes, J. Chem. Phys., 1970, 52, 3121.
17. S. Scheiner and M. Yi, J. Phys. Chem., 1996, 100, 9235.
18. G. Will, Z. Kristallogr., 1969, 129, 211.
19. J. A. Harvey, D. Basak, D. Venkataraman and S. M. Auerbach, Mol. Phys., 2012, 110, 957.
20. J. Thisuwan and K. Sagarik, RSC Adv., 2014, 4, 61992.
21. F. Weigend and M. Haeser, Theor. Chem. Acc., 1997, 97, 331.
22. C. Lao-ngam, M. Phonyiem, S. Chaiwongwattana, Y. Kawazoe and K. Sagarik, Chem. Phys., 2013, 420, 50.
23. S. Chaiwongwattana, M. Phonyiem, V. Vchirawongkwin, S. Prueksaaroon and K. Sagarik, J. Comput. Chem., 2012, 33, 175.
24. TURBOMOLE V7.0 2015, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007, available from http://www.turbomole.com.
25. C. Lao-Ngam, P. Asawakun, S. Wannarat and K. Sagarik, Phys. Chem. Chem. Phys., 2011, 13, 4562.
26. T. C. Berkelbach, H.-S. Lee and M. E. Tuckerman, Phys. Rev. Lett., 2009, 103, 238302.
27. See http://krohne.com/en/services/dielectric-constants/.
28. A. P. Scott and L. Radom, J. Phys. Chem., 1996, 100, 16502.
29. R. Konrat, M. Tollinger, G. Kontaxis and B. Kraeutler, Monatsh. Chem., 1999, 130, 961.
30. R. Ditchfield, Mol. Phys., 1974, 27, 789.
31. R. Ditchfield, J. Chem. Phys., 1976, 65, 3123.
32. I. Alkorta and J. Elguero, New J. Chem., 1998, 22, 381.
33. D. B. Chesnut and B. E. Rusiloski, J. Mol. Struct.: THEOCHEM, 1994, 314, 19.
34. D. B. Chesnut, J. Phys. Chem. A, 2005, 109, 11962.
35. S. Taubert, H. Konschin and D. Sundholm, Phys. Chem. Chem. Phys., 2005, 7, 2561.
36. S. Nosé, J. Chem. Phys., 1984, 81, 511.
37. S. Nosé, Progress of Theoretical Physics Supplements, 1991, 103, 1.
38. S. Nosé and M. L. Klein, Mol. Phys., 1983, 50, 1055.
39. M. Phonyiem, S. Chaiwongwattana, C. Lao-ngam and K. Sagarik, Phys. Chem. Chem. Phys., 2011, 13, 10923.
40. P. Bopp, Chem. Phys., 1986, 106, 205.
41. P. Suwannakham, S. Chaiwongwattana and K. Sagarik, Int. J. Quantum Chem., 2015, 115, 486.
42. T. Murakhtina, J. Heuft, E. J. Meijer and D. Sebastiani, ChemPhysChem, 2006, 7, 2578.
43. Y. J. Lee, B. Bingöl, T. Murakhtina, D. Sebastiani, W. H. Meyer, G. Wegner and H. W. Spiess, J. Phys. Chem. B, 2007, 111, 9711.
44. S. R. Narayanan, S. P. Yen, L. Liu and S. G. Greenbaum, J. Phys. Chem. B, 2006, 110, 3942.
45. H. Chen, T. Yan and G. A. Voth, J. Phys. Chem. A, 2009, 113, 4507.
46. T. Lankau and C.-H. Yu, Phys. Chem. Chem. Phys., 2011, 13, 12758.
47. S. Kong, A. O. Borissova, S. B. Lesnichin, M. Hartl, L. L. Daemen, J. Eckert, M. Y. Antipin and I. G. Shenderovich, J. Phys. Chem. A, 2011, 115, 8041.
48. K. Sagarik, M. Phonyiem, C. Lao-ngam and S. Chaiwongwattana, Phys. Chem. Chem. Phys., 2008, 10, 2098.
49. R. Bouchet and E. Siebert, Solid State Ionics, 1999, 118, 287.
50. K. D. Kreuer, A. Fuchs, M. Ise, M. Spaeth and J. Maier, Electrochim. Acta, 1998, 43, 1281.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra17636d

---