Proton dissociation and transfer in a phosphoric acid doped imidazole system

Jittima Thisuwan and Kritsana Sagarik*
School of Chemistry, Institute of Science, Suranaree University of Technology, Nakhon Ratchasima 30000, Thailand. E-mail: kritsana@sut.ac.th; Fax: +66 44224635; Tel: +66 44224635

Received 6th August 2014 , Accepted 10th November 2014

First published on 11th November 2014


Abstract

The dynamics and mechanisms of proton exchange in a phosphoric acid (H3PO4) doped imidazole (Im) system were studied using a quantum chemical method at the B3LYP/TZVP level and Born–Oppenheimer molecular dynamics (BOMD) simulations. The theoretical studies began with selecting the appropriate presolvation models for proton dissociation and transfer, which were represented by embedded and terminal hydrogen bond (H-bond) structures of H+(H3PO4)(Im)n(n = 1–4), respectively. B3LYP/TZVP calculations confirmed that excess proton conditions are required to promote proton exchange, and the intermediate complexes are preferentially formed in a low local-dielectric environment. In contrast, a high local-dielectric environment is required to stabilize the positive charge and prevent the proton from returning to the original Im molecule. The static results also revealed that the embedded structure with n = 2 represents the smallest, most active intermediate complex for proton dissociation, whereas the terminal structure with n = 3 favors proton transfer in the Im H-bond chain. BOMD simulations confirmed the static results and further suggested that the fluctuations of the H-bond chain lengths and the local-dielectric environment must be included in the proton dissociation and transfer mechanisms. Analyses of the time evolutions of torsional angles in the Im H-bond chain showed characteristic solvent structure reorganization, regarded as helical-rotational motion, which drives the proton away from the H3PO4 dopant. The current theoretical results showed in detail for the first time the interplay among the “key molecular motions” that fundamentally underlie the dynamics and mechanisms of proton exchange in the H3PO4 doped Im H-bond system.


Introduction

Proton exchange membrane fuel cells (PEMFCs) have received considerable attention due to their potential as effective power generators.1–3 Commercially available PEMFCs employ hydrated Nafion® as a polymer electrolyte membrane and therefore have serious limitations on the operating temperature, as Nafion® based fuel cells cannot be operated above the boiling point of water. Over the past decade, attempts have been made to develop electrolyte membranes with high proton conductivity and thermal and mechanical stabilities at elevated temperatures. Aromatic heterocyclic compounds such as benzimidazole (BI), pyrazole (Pz) and imidazole (Im), as well as polymers containing amide functional groups such as poly(oxodiazole)(POD) and poly(benzimidazole)(PBI), have shown high potential as alternative proton solvents.4 Due to the autoprotolysis property, a high amphotericity and being a good proton solvent, as well as the ability to form an extensive H-bond network, these aromatic heterocyclic materials possess high proton conductivities in the condensed phase.4

PBI is an aromatic heterocyclic polymer that exhibits excellent thermal, mechanical and chemical stabilities.4 Although PBI is an intrinsic proton conductor, experiments have revealed that blank PBI is a proton insulator,5–7 and a sample preparation in which the PBI membrane is dipped into an acid solution represents an essential step for the conductivity measurements.7,8 Doping a PBI membrane with acids such as H2SO4, H3PO4 and HClO4 can considerably increase the proton conductivity, among which PBI doped with H3PO4 yields the highest proton conductivity.7 In the pure state, because all hydrogen atoms are capable of forming extensive H-bonds, H3PO4 shows the highest intrinsic proton conductivity.9 However, it was demonstrated that the proton conductivity of H3PO4 decreases considerably when strong electrolytes such as H2SO4 were added.10 In contrast, the conductivity was found to increase significantly with water content. The latter is attributed to the promotion of self-ionization and an increase in the number of protons.8

Given the high proton conductivity with respect to the c crystallographic direction11,12 and properties as an excellent proton solvent, Im has been frequently employed as a prototypical molecule in the study of proton transfer in aromatic heterocyclic materials.11–16 Experiments have shown that Im can form extensive H-bond chains in the condensed phase, along which protons can be conducted by the Grotthuss mechanism.11,12 It was suggested that to complete the proton transfer process, the Grotthuss mechanism must be followed by a cooperative ring reorientation (ring-flip) as an important rate-determining step.16 However, solid-state 15N NMR experiments demonstrated that for pure Im, the Grotthuss mechanism need not involve the reorientation of the Im ring.11

The influence of Im and 1-methyl Im on the conductivity of 99 wt% H3PO4 solutions and H3PO4 doped PBI membranes was studied over the temperature range of 350 to 473 K under various humidity conditions.17 It appeared that the conductivity of the solutions decreased with an increasing amount of Im and increased with temperature. The former is attributed to a reduction in the number of proton charge carriers as a consequence of the acid–base reaction, and the viscosity of the solutions was suggested to play an important role in the proton conduction, especially at temperatures below 373 K. Because proton conduction in pure H3PO4 is perturbed by the addition of Im,17 to investigate the mechanism underlying this effect, ab initio molecular dynamics (AIMD) simulations on a H3PO4-Im 2[thin space (1/6-em)]:[thin space (1/6-em)]1 system (36 H3PO4 and 18 Im molecules) were performed at 400 K.18 Due to the acid–base reaction, stable ion-pair complexes between H2PO4 and ImH+ and Zundel-like ions, characterized by proton shuttling between two Im molecules (H+(Im)2), were observed in AIMD simulations. Analyses of the local structures and energetic and dynamic results suggested that, due to the proton transfer from H3PO4 to Im, the proton density in the H3PO4 phase is reduced and the formation of the ion-pair complexes leads to significantly shorter H-bond chains compared to those in the pure state. These reduce the cooperative effects of the proton transport mechanism and the effectiveness of proton transfer in the Im H-bond chain.18

Insights into proton dissociation and transfer in an aqueous solution were obtained in our previous work19 in which the dynamics and mechanisms of proton exchange in hydrated H3PO4 under excess proton conditions were studied using H+(H3PO4)(H2O)n(n = 1–4) as the model system with ab initio calculations and Born–Oppenheimer molecular dynamics (BOMD) simulations at the RIMP2/TZVP level as the model calculations. Based on the static results in the gas phase and continuum aqueous solvent, mechanisms of proton dissociation and transfer that involve fluctuations of the number of water molecules in the H-bond chain and the local-dielectric environment were proposed. It was concluded that the smallest, most active intermediate complex for proton dissociation (n = 1) is formed in a low local-dielectric environment (e.g., ε = 1) and that the proton transfer from the first hydration shell is driven by fluctuations in the number of water molecules in a high local-dielectric environment (e.g., ε = 78) through the Zundel complex, which connects the first and second hydration shells. BOMD simulations over a temperature range of 298 to 430 K confirmed that only the H-bond chains with asymmetric hydration at H4PO4+ are involved in the proton dissociation and transfer. This finding validated the proposed mechanisms by showing that an agreement between the theoretical and experimental activation energies can be achieved based on the proposed rate-determining processes.19

In the work presented here, the dynamics and mechanisms of proton dissociation and transfer in a H3PO4 doped Im H-bond system were studied using the quantum chemical method at the B3LYP/TZVP level and BOMD simulations. Based on the concept of presolvation20,21 and the observation that a unit cell of the Im crystal structure is monoclinic with four molecules forming a linear H-bond chain along the c crystallographic axis, H+(H3PO4)(Im)n(n = 1–4) were selected as model systems. The characteristic properties of the exchange protons in the intermediate complexes in low and high local-dielectric environments were investigated and analyzed in detail. The elementary steps of proton dissociation and transfer were proposed based on the static results and confirmed by BOMD simulations in the temperature range between 298 and 423 K. The interplay among the key molecular motions that fundamentally underlie the dynamics and mechanisms of proton exchange in the H3PO4 doped Im H-bond system were analyzed and discussed.

Computational methods

Presolvation models

Because the proton exchange processes in H-bonds are complicated, care must be exercised in selecting the model systems. In the present case, the model systems for proton dissociation and transfer were chosen based on the concept of presolvation;20,21 the local energy fluctuations and dynamics in solution lead to the weakening or breaking of some H-bonds in the first and second solvation shells, resulting in a reduction in the solvent coordination number such that the proton-accepting species possesses a solvation structure corresponding to the species into which it will be transformed to complete the structural diffusion process. It was shown in our previous work that the dynamics of the transferring protons can be studied reasonably well because the moderate sizes of the presolvation models help reduce the complexity in the H-bond systems and allow theoretical calculations to be made using quantum chemical methods at higher level of accuracy, provided that an appropriate presolvation model is selected and that all of the important rate-determining processes are included in the model calculations.22–26

In the present work, because a unit cell of the Im crystal structure is monoclinic with four Im molecules forming a linear H-bond chain along the c crystallographic axis, H+(H3PO4)(Im)n(n = 1–4) were selected as model systems. The cluster size (n = 4) is consistent with the average number of Im molecules obtained from MD simulations with excess proton solvated in liquid Im, using a reactive multistate empirical valence bond (MS-EVB) model.27 Based on the concept of presolvation and some test calculations, the presolvation models for proton dissociation and transfer were proposed and are shown in Fig. 1, regarded hereafter as the “embedded” and “terminal” linear H-bond structures, respectively. Experimental and theoretical studies which support the linear H-bond arrangements in the presolvation models are summarized as follows.


image file: c4ra08198f-f1.tif
Fig. 1 (a and b) Embedded and terminal presolvation models employed in the studies of proton dissociation and transfer in H+(H3PO4)(Im)n(n = 1–4), respectively. 1st = first solvation shell; 2nd = part of second solvation shell; (1)–(7) = labels of H-bonds; (I)–(VII) = labels of Im molecules.

In condensed phase, due to the amphoteric character, computer simulations28 and experiments12 showed that Im molecules form polymeric chains, along which excess proton can be relayed between two Im molecules through the Grotthuss shuttling mechanism. Possessing only one proton donor (N–H) and one proton acceptor (N) in the same molecule, the linear H-bond chains in liquid Im are not as extensive as the H-bond networks in water (two proton donors and two proton acceptors in the same molecule). Because the proton donor and acceptor are used in the H-bond chain formation, the interactions between parallel chains could be the C–H⋯π H-bond and π⋯π interactions, which are considerably weaker than the interactions in the protonated Im H-bond chains; for the benzene dimers, the interaction energies for the C–H⋯π H-bond in the T-shaped structure, and the π⋯π interaction in the sandwich structure, computed at the MP2 level, are only −8.8 and −3.6 kJ mol−1,29 respectively, compared to the interaction energy of the protonated Im dimer of −116 kJ mol−1.30

Our theoretical study on proton transfer in protonated Im H-bond clusters30 revealed that the intermediate complexes (the H-bond structure with a shared-proton) exist only in the linear H-bond chains (see Table S2), which is in accordance with the protonic-chain conduction mechanism;10,31 B3LYP/TZVP geometry optimizations suggested that the linear H-bond arrangements are the global and local minimum energy geometries in the protonated pure,30 and doped Im systems. It should be noted that, although the globular H-bond structures could possess larger contact area and smaller free energy, based on the characteristic properties of the transferring proton (e.g., asymmetric stretching coordinates and vibrational frequencies), the H-bonds in the globular structures are not susceptible to proton transfer.30 Because we are interested in the dynamics of the transferring protons in the intermediate complexes, linear Im H-bond chains were selected as the presolvation models; the lifetime of the intermediate complex has been shown in many H-bond systems to represent the rate-determining processes in proton transfer mechanisms.

To investigate solvent effects on proton dissociation and transfer, a continuum solvent model was included in the model calculations. To account for the electrostatic effects introduced by the solvent, a conductor-like screening model (COSMO), with the dielectric constant of liquid Im (ε = 23), was employed. The application of the continuum solvent model is justified by the following studies.

In liquid phase, analyses of a MS-EVB model revealed that the first solvation shell of the ImH+ ion is highly ordered due to the formation of the strong protonated H-bond (ImH+⋯Im), whereas the second solvation shell is highly disordered.32 These suggest that the interactions connecting the Im molecules in the first and second solvation shells are weak. As the first and second solvation shells are loosely connected, mechanical and chemical interactions between the two solvation shells could be approximately neglected in the model systems, and it is reasonable to incorporate the electrostatic effects introduced by the weak-polar organic solvent by a continuum solvent model.

MD simulations also suggested that the proton transfer in liquid Im can be considered as a “local event”,16,27 with very short spatial/temporal correlation27 (involving only Im molecules close to ImH+) and driven by local rather than long-range cooperative motions (reorientations) of the Im molecules.16 Based on Car-Parrinello molecular dynamics (CPMD) simulations,16 the time scale for the long-range cooperative reorientation process is about 30 ps, which is considerably longer than the proton transfer step of about 0.3 ps. Therefore, the effects of the fluctuations of bulk solvent structures on the conformations of the solute (ImH+⋯Im) could be neglected, and the approximation of the electrostatic effects of the local-dielectric environment by a continuum solvent model is reasonable and justifiable.

Quantum chemical calculations

Three basic steps used in our previous work24–26,33,34 were employed: (1) searching for equilibrium structures and intermediate complexes in the proton dissociation and transfer pathways using the Test-particle model (T-model) potentials;35–39 (2) refining the T-model optimized structures using the quantum chemical method; and (3) performing BOMD simulations on the candidate intermediate complexes using the refined structures as the starting configurations.

It is worth remarking on the quantum chemical method employed in this study. As discussed in the previous work,24 because of the ability to approximate the effects of electron correlation with reasonable accuracy and computational resources, density functional theory (DFT) methods have been used extensively in studies of interacting systems. However, the accuracy of the DFT methods depends on the chemical systems under investigation; the DFT methods tend to underestimate the interaction energy (ΔE) for the systems with strong effects of the electron correlation. The T-shaped and parallel-displaced (PD) structures of a phenol-benzene dimer are a typical example in which the interplay between the electrostatic and dispersion interactions was examined in detail;40 the stability of the former is dominated by the electrostatic interaction and the latter by the dispersion interaction. However, B3LYP/6-311++G(d,p) calculations predicted the PD structure to be unstable, with positive interaction energy (ΔE), whereas ΔE of the T-shaped structure with the O–H⋯π interaction is considerably higher than the experimental value, −8.4 kJ mol−1 and −16.7 kJ mol−1,41 respectively. Because the effects of the electron correlation could play an important role in the present model systems, the applicability and performance of the B3LYP/TZVP calculations were examined by benchmarking with the second-order Møller-Plesset perturbation theory (MP2) method, with the resolution of the identity (RI) approximation and the TZVP basis set, with abbreviated RIMP2/TZVP calculations. The correlations between the key properties, e.g., the H-bond structures, interaction energies and vibrational frequencies, are illustrated in Fig. S1. The methods were in agreement and obtained correlation coefficient (R2) for ΔdDA of 0.97. Therefore, to ensure reasonable computational efforts, the B3LYP/TZVP calculations were employed in both the static and dynamic calculations. All of the quantum chemical calculations and BOMD simulations were performed using the TURBOMOLE 6.5 software package.42

Static calculations

To obtain the fundamental information for the discussion of the dynamics and mechanisms of proton dissociation and transfer in the H3PO4 doped Im system under an excess proton condition, the structural, energetic and spectroscopic properties of H+(H3PO4)(Im)n(n = 1–4) were computed using B3LYP/TZVP calculations, both in the gas phase and the continuum solvent (COSMO, ε = 23). The interaction energies (ΔE) of the H-bond complexes were computed from ΔE = E(H+(H3PO4)(Im)n) − [E(H4PO4+) + nE(Im)]; E(H+(H3PO4)(Im)n) is the total energy of the equilibrium structure of H+(H3PO4)(Im)n; E(H4PO4+) and E(Im) are the total energies of the isolated species at their equilibrium structures. The solvation energies (ΔEsol) were approximated from ΔEsol = E(H+(H3PO4)(Im)n)COSMO − E(H+(H3PO4)(Im)n); E(H+(H3PO4)(Im)n)COSMO and E(H+(H3PO4)(Im)n) are the total energies of the equilibrium structures obtained with and without COSMO, respectively.

The behaviors of proton dissociation and transfer were primarily evaluated from the asymmetric stretching coordinate (ΔdDA), computed from ΔdDA = |dA–HdB⋯H|, where dA–H and dB⋯H are the A–H and B⋯H distances in the A–H⋯B H-bond, respectively. The H-bond that is susceptible to proton transfer is represented by a ΔdDA value that is close to zero as for a proton located near or at the center of the H-bond. As an effective probe for the dynamics of proton transfer in the H-bond, this study employed vibrational spectroscopy. The vibrational frequencies of protons in H+(H3PO4)(Im)n(n = 1–4) were computed using the B3LYP/TZVP method with the harmonic approximation. They were used as guidelines to assign the vibrational modes predicted by the BOMD simulations. As the anharmonicity is neglected, the vibrational frequencies obtained from quantum chemical methods are generally overestimated. A scaling factor of 0.9614 proved to be appropriate for the B3LYP/TZVP calculations43 and was therefore employed in this work. The dynamics of the proton dissociation and transfer were discussed based on the asymmetric O–H and N–H stretching frequencies, νOH and νNH, respectively24,26,44,45 Additionally, because the motions of protons moving in single- and double-well potentials are different, H-bonds that are susceptible and not susceptible to transfer were distinguished by the threshold asymmetric N–H stretching frequencies (νNH*) that were calculated from the plots of νNH and ΔdDA.33 The relationship between νNH and ΔdDA is described by a reflected normal distribution function, for which the second derivative is equal to zero at νNH = νNH*. The H-bond proton is considered to be susceptible to transfer when νNH < νNH*.33 To provide energetic information for the discussion of the solvent effects, two-dimensional potential energy surfaces (2D-PES) of the transferring protons in the smallest, most active intermediate complexes were constructed in both the gas phase and the continuum solvent by calculating ΔE at various RN–O and RN–H distances. As in our previous studies,23 equally spaced grid points were generated from ΔE and the 2D-PES were plotted as a function of RN–O and ΔdDA using the SURFER program.46 The interaction energy paths connecting the single- and double-well potentials, as well as the path with the highest energy barriers, were identified and discussed in detail.

Dynamic calculations

The dynamics and mechanisms of proton transfer in H-bonds were successfully studied in our previous work by performing NVT BOMD simulations on the smallest, most active intermediate complexes in the presolvation models.19,23–26 In the present work, we are interested in the dynamics of the transferring protons (e.g., the oscillatory shuttling and structural diffusion motions) in the intermediate complexes, which occur locally in a very short time scale; according to the MS-EVB MD simulations,27 these involve linear H-bond chains with few Im molecules. As the intermediate complexes with the shared-proton structures are repeatedly formed and relaxed, and BOMD simulations involve only limited number of Im molecules, the periodic boundary conditions are not necessary in this case; the periodic boundary conditions are generally employed in computer simulations of condensed phase to solve the problem of surface and finite-size effects, and therefore important in the study of the long-range effects in bulk systems.

To study proton dissociation and transfer in H+(H3PO4)(Im)n(n = 1–4), NVT BOMD simulations were performed at the B3LYP/TZVP level. Two thousand steps of 1 fs each were used in the equilibration and 10[thin space (1/6-em)]000 steps were used in the property calculations, corresponding to 2 and 10 ps, respectively. The NVT ensemble was concluded to be appropriate for the dynamic processes that take place on relatively short time scales and in the continuum solvent,47 and the time step proved to be small enough to determine reliable vibrational spectra, activation energies and dynamic properties of the transferring protons.22–26,48,49

Because realistic dynamics at constant temperatures are desirable in this work, the BOMD simulation conditions were carefully chosen, especially the temperature control algorithm. The Nosé–Hoover thermostat was selected based on the conclusion that the transport properties and thermodynamic distributions that are statistically indistinguishable from the NVE ensemble can be achieved when globally applied to the system using any coupling strength;47 the Nosé–Hoover thermostat adjusts velocities by coupling the equations of motion to extended dynamic variables and gives an oscillatory relaxation to the target temperature. To ensure the minimum perturbation on the dynamics of the transferring protons, the energy conservation and stability of the numerical integration in the NVT BOMD simulations were assessed at 323 K, using the smallest, most active intermediate complex as a representative. To test the stability of the numerical integration, a drift in property A (δA) over a time interval [0, tmax] was computed from δA = bAtmax/σA;50 bA is obtained from a linear least-square fit to the data, A(t) = aA + bAt + ε(t), and σA is the standard deviation of ε(t). Therefore, for an ideal, isolated-energy conserved system, bA is 0, and hence, | δA| = 0. The drift in property A is significant if | δA| ∼ 1 and the result obtained by averaging A(t) over the simulation time tmax is considered to be unreliable. Based on this definition, the drifts in energies (Etot, Epot and Ekin), the H-bond structures (RN–O, RN–H and RO–H) and the velocity of the transferring proton (vH+) were calculated over the integration period of tmax = 2 ps and included in Table S1. The energy conservation and velocity plots over the entire simulation length (10 ps) are illustrated in Fig. S2. Additionally, to confirm that the dynamics are not perturbed by the Nosé–Hoover thermostat, the mean-square displacement (MSD) of the transferring proton was computed and included in Fig. S2. Because the trends of the energy conservation and velocity plots showed no systematic drifts and a general linear relationship was observed for the MSD plot, the dynamic properties computed under the present NVT BOMD simulation conditions, including the Nosé–Hoover thermostat, were determined to be well represented. Additionally, as the energy and velocity drifts are only 0.22, the numerical integration in the NVT BOMD simulations is confirmed to be stable.

The dynamic behaviors of proton dissociation and transfer in the H3PO4 doped Im system were primarily discussed using proton transfer profiles, with variations of RA–B, RA–H and RB–H with the simulation time.49 Proton transfer reactions are generally not concerted and can take place in both the forward and reverse directions, with the oscillatory shuttling motion being the rate-determining process.49 This is because the transfer reactions involve quasi-dynamic equilibrium between the close-contact A–H⋯B and the shared-proton A⋯H⋯B H-bonds.49 It is therefore not straightforward to determine the direction of the proton displacement in the Im H-bond chain. Because a stable covalent bond must be formed in the direction of the proton displacement to complete proton hopping in a single H-bond, it is reasonable to correlate the probability of finding a proton moving in a specific direction with the probability of covalent bond formation. Hence, to estimate the probability of finding a proton transfer from atom B to A (PBA), the distribution of the A–H covalent bond distance (RA–H) is first computed from the BOMD simulations. Then, based on the criterion for the N–H covalent bond formation (RN–H ≤ 1.3 Å), PBA is approximated by the number of counts of the N–H covalent bond divided by the total number of samples. In this work, a proton is considered to have favorably transferred from B to A when PBA ≥ 0.5. Likewise, the probability of finding the intermediate complex with the shared-proton A⋯H⋯B H-bond (PA⋯H⋯B) can be approximated by the number of counts of RA–H in the range between 1.3 and 1.5 Å, divided by the total number of samples.

The dynamics and mechanisms of proton dissociation and transfer were finally confirmed using the vibrational spectra obtained from Fourier transformations of the velocity autocorrelation functions (VACF) of the atoms in H-bonds.24,25 Only the characteristic N–O and N–N vibrational frequencies (νNO,MD and νNN,MD, respectively), as well as the asymmetric N–H stretching frequencies (νNH,MD), were reported and discussed in detail. The activation energies of the rate-determining processes in the proton dissociation and transfer pathways were determined by performing BOMD simulations over the temperature range of 298 to 423 K. The first-order rate constants (k) for proton exchanges were approximated from the exponential relaxation behavior of the envelope of VACF of the N–O and N–N vibrations, which are based on the assumption that the decay process follows the first-order kinetics: k = ln(2)/τ1/2, where τ1/2 is the half-life of the exponential decay function. The Arrhenius plots were constructed and the activation energies (ΔE‡,Arr) of proton dissociation and transfer were approximated from the slopes of the linear relationship between ln(k) and 1000/T.23

Results and discussion

Static results

The equilibrium structures and the interaction (ΔE) and solvation (ΔEsol) energies obtained from the B3LYP/TZVP calculations in the gas phase and the continuum solvent are listed in Table S2, together with the H-bond distances (RN–N, RN–O and ΔdDA) and the asymmetric N–H stretching frequencies (νNH). The H3PO4 doped Im H-bond structures in Table S2 are labelled with a four character code, e.g., Gn-[m](I–VII) and Cn-[m](I–VII), where G = equilibrium structure in the gas phase; C = equilibrium structure in continuum solvent; and n = number of Im molecules. Different H-bond structures with the same number of Im molecules are distinguished by [m]. Positively charged Im molecules are designated by a Roman numeral in parenthesis, from (I) to (VII). As an example, G3-[2](I) and G3-[2](II) are the same H-bond structures (m = 2 and n = 3) with different positively charged Im molecules, at Im (I) and Im (II), respectively.

Excess proton conditions and tendencies of proton dissociation and transfer

Investigation of the embedded structures and their H-bond properties in Table S2 reveals outstanding acid–base interaction in the H3PO4 doped Im system, specifically two protons are readily dissociated from H4PO4+ and protonate the two neighboring Im molecules (Im (I) and Im (II)), forming the ImH+-H2PO4-ImH+ complex, e.g., in the gas phase, structures G2-[1](I,II) and G3-[1](I,II), and in the continuum solvent, structures C2-[1](I,II) and C3-[1](I,II). Because the B3LYP/TZVP calculations (not shown here) suggest no proton dissociation from the corresponding unprotonated (neutral) embedded structures, one can conclude that proton dissociation and transfer in the Im system takes place preferentially under excess proton conditions,8,9 with the Im-H4PO4+-Im complexes as the precursors. This is in agreement with experiments in which dipping the PBI film into a H3PO4 solution is a necessary step for the conductivity measurement.8 Therefore, to obtain realistic proton dissociation and transfer mechanisms, excess proton conditions must be included in the presolvation models.19,24,25,33

The tendencies of proton dissociation and transfer are primarily discussed using the relationships between νNH and ΔdDA in Fig. 2. The plot of νNH and ΔdDA in Fig. 2a suggests that the threshold asymmetric N–H stretching frequency (νNH*) for the proton dissociation in the gas phase is at νNH* = 2158 cm−1 and Δd*DA = 0.37 Å, approximately 150 cm−1 lower than in aqueous solution (νOH* = 2306 cm−1),19 reflecting a stronger acid–base interaction in the H3PO4 doped Im system. In the continuum solvent, the νNH values are all above 2500 cm−1 without an inflection point (νNH*), suggesting that the close-contact N–H⋯O H-bond with localized positive charge is more favorable than the shared-proton structure (N⋯H⋯O) with delocalized positive charge in COSMO, which is in accordance with the localized charge solvent (LCS) concept: “a polar solvent solvates species with localized charges more efficiently than species with more delocalized charges”.51 Similar results were obtained for proton transfer in the Im H-bond chains, for which the plot of νNH and ΔdDA in Fig. 2b yielded the threshold frequency only in a low local-dielectric environment, in the gas phase at νNH* = 1925 cm−1 and Δd*DA = 0.44 Å. The threshold frequency is approximately 70 cm−1 lower than that of the undoped Im H-bond chains (1992 cm−1),30 reflecting the strong polarization effect of the dopant and a higher probability of finding the shared-proton structure and delocalized protons in the H3PO4 doped Im system; the shared-proton structure is an important intermediate in the proton transfer pathway.


image file: c4ra08198f-f2.tif
Fig. 2 Plots of asymmetric N–H stretching frequencies (νNH) as a function of asymmetric stretching coordinates (ΔdDA), obtained from B3LYP/TZVP calculations. The dashed lines separate the protons that are susceptible and not susceptible to exchange. (a) N–H⋯O H-bonds in the first solvation shell. (b) N–H⋯N H-bonds in Im H-bond chains.

Two-dimensional potential energy surfaces and solvent effects

The two-dimensional potential energy surfaces (2D-PES) for proton dissociation in the smallest embedded structure in the gas phase and the continuum solvent are shown in Fig. 3. The absolute minimum on the 2D-PES of structure G2-[1](I,II) is represented by an asymmetric single-well potential at RN–O = 2.54 Å and ΔdDA = 0.33 Å. At longer RN–O distance (RN–O ≈ 2.67 Å), a double-well potential emerges. On the 2D-PES, a low-interaction energy, barrier-less path connecting the single- and double-well potentials spans from A to C, regarded hereafter as the “dissociation path” (or path AC), whereas the “association path” spanning from A to B (or path AB) is uphill and therefore unfavorable in this H-bond complex. The shapes of the cross-sectional plots in Fig. 3a suggest at least three N–H stretching modes associated with proton dissociation, namely a vibrational mode associated with a proton moving in the single-well potential and two vibrational modes with a proton moving on the dissociation and association paths. Because conventional quantum chemical calculations on vibrational spectra are based on the harmonic approximation, which requires well-optimized geometry, only the νNH value of a proton moving in the single-well potential is obtainable from the static calculations; for proton dissociation in structure G2-[1](I,II), νNH = 2033 cm−1. Therefore, BOMD simulations must be performed to study the other two characteristic modes that involve H-bond structures distorted from the absolute minimum energy geometry.
image file: c4ra08198f-f3.tif
Fig. 3 Two-dimensional potential energy surfaces (2D-PES) of proton dissociation in the gas phase and continuum solvent. (a) Proton dissociation in H-bond (1) of the embedded structure G2-[1](I,II). (b) Proton dissociation in H-bond (1) of the embedded structure C2-[1](I,II). (a), (b) and (c) are the cross-sectional plots of the single- and double-well potentials, respectively.

A comparison of the cross-sectional plots in Fig. 3a and b reveals the effects of the continuum solvent on the shapes of the 2D-PES, namely the single- and double-well potentials are shifted closer to the Im molecules, suggesting stronger N–H covalent bonds and weaker N–H⋯O H-bonds in structure C2-[1](I,II). In other words, due to the partial neutralization of positive charges by the continuum solvent, protons are localized nearer to the nitrogen atoms, resulting in a lower probability of finding the shared-proton structure and the oscillatory shuttling motion. These findings confirm that the intermediate complexes of proton dissociation and transfer in the H3PO4 doped Im system are favorably formed in a low local-dielectric environment and that a high local-dielectric environment is required to stabilize the transferred proton before the successive transfer.

Mechanisms of proton transfer in the Im H-bond chains

Apart from the excess proton condition and the effects of the local-dielectric environment, to propose the elementary steps of proton transfer in the H3PO4 doped Im system, the effects of fluctuations of the H-bond chain lengths are discussed in detail. In this case, ΔdDA of the embedded structures G2-[1](I,II), G3-[1](I,II) and G4-[2](I,II) are analyzed and the effects of the H-bond chain extension are shown schematically in Fig. 4a. Starting from the smallest, most active embedded structure G2-[1](I,II), an extension of the Im H-bond chain at Im (II) leads to structure G3-[1](I,II), with a decrease in ΔdDA of H-bond (1) from 0.33 to 0.26 Å and an increase in ΔdDA of the H-bond (2) from 0.36 to 0.46 Å. These findings suggest that an increase in the Im H-bond chain length at one side of the embedded structure results in a stronger N–H⋯O H-bond on the other side. Further extension of the Im H-bond chain at Im (IV) of structure G3-[1](I,II) yields structure G4-[2](I,II), with a decrease in the ΔdDA of H-bond (1) from 0.26 to 0.23 Å and an increase in the ΔdDA of H-bond (2) from 0.46 to 0.50 Å. Because the two N–H⋯O H-bonds in the embedded structures are simultaneously stabilized and destabilized by the H-bond chain extension, it is reasonable to anticipate that terminal structures, which are precursors for proton transfer in the Im H-bond chain, can be formed in a low local-dielectric environment. In a high local-dielectric environment, a comparison of the ΔdDA of H-bonds (1) and (2) in structures C2-[1](I,II), C3-[1](I,II) and C4-[1](I,II) reveals different results. Because the protonated H-bonds are weakened in a continuum solvent, the ΔdDA values are not sensitive to the H-bond chain extension as in the gas phase. Therefore, the successive addition of Im molecules does not lead to significant changes in the ΔdDA. However, structures C2-[2](III), C3-[2](V) and C4-[3](VII) in Table S2 suggest that the outstanding effect of the continuum solvent is to stabilize the positive charge at the end of the Im H-bond chains.
image file: c4ra08198f-f4.tif
Fig. 4 (a) The effects of Im H-bond chain extension in low local-dielectric environment. (b) Examples of elementary steps of proton transfer in Im H-bond chains. step 1 = Terminal structure formation. step 2 = Proton transfer in Im H-bond chains.

Because the intermediate complexes (strong H-bonds with shared-proton structures) are favorably formed in a low local-dielectric environment and a high local-dielectric environment is required to stabilize the positive charge and drive protons in the Im H-bond chain, the fluctuations of the local-dielectric environment must be included as a part of the mechanisms of proton dissociation and transfer. Therefore, based on the analyses of the static results, examples of the elementary steps, which involve the fluctuations of the H-bond chain lengths (n = 2 to 4) and the local-dielectric environment (ε = 1 and 23) are proposed in Fig. 4b. The elementary steps show the terminal structure formation in step 1 and proton displacement to the end of the Im H-bond chains in step 2.

Similar effects of fluctuations of the H-bond chain lengths and dopant were reported by Sadeghi et al.52 in which BOMD simulations revealed that proton transfer in the protonated linear H-bond chains of water molecules (n = 6) occurs through a series of semi-collective motions, resulting from the fluctuations of the H-bond chain lengths and reorganization of the water molecules. Additionally, based on the BOMD results on the protonated linear H-bond chains of water molecules (n = 5) with an ammonia molecule (NH3) attached to one end of the chain, it was demonstrated that the doping effect is described by proton transfer throughout the H-bond chain until it reaches the NH3 dopant.52 Because the dopant in the present system is an acid, proton transfer takes place in the opposite direction, from the proton donor (H4PO4+) to the end of the Im H-bond chains, as shown in Fig. 4b.

It should be noted that due to the thermal energy fluctuations and dynamics, as well as the local geometrical requirements, H-bonds in liquids can be considered as “forming” and “breaking”.53 For liquid Im in excess proton conditions, MD simulations based on an MS-EVB model suggested continuous forming and breaking of protonated H-bond chains in which the average H-bond cluster size is approximately four (H+(Im)4). The proton transfer in the Im H-bond chain is therefore considered to be a local event, characterized by the fluctuations of the H-bond chain lengths between the Eigen-like (n = 3) and Zundel-like (n = 2) complexes.27 Because the fluctuations of the Im H-bond chain lengths and the structure reorganization lead to fluctuations of the local polarization that determine the local-dielectric environment, and because the dielectric spectra of imidazolium compounds suggested that the dipolar relaxation times can be quite short54 and could synchronize with the lifetimes of the intermediate complexes in the proton transfer pathways, it is justifiable to include the effects of the local-dielectric environment in the proton dissociation and transfer mechanisms.

Effects of Im H-bond chain structure reorganization

To discuss the effects of Im H-bond chain structure reorganization, the relative energy diagrams and potential energy curves for proton displacement in H+(H3PO4)(Im)4 are constructed and shown in Fig. 5a. The relative energy diagram shows that structures G4-[3](I), G4-[3](III) and G4-[3](V) are close in energy and are accessible through the Im H-bond chain structure reorganization, e.g., through the rotations of the dihedral angles ω3 and ω5. It appears that in a low local-dielectric environment (ε = 1), the orientations of the Im molecules, especially Im (III) and Im (V), determine the relative stabilities and the shapes of the potential energy curves that define the energy barriers for proton displacement; e.g., to convert structure G4-[3](I) to G4-[3](III), ω5 decreases from 107 to 93 degrees. The rotation of ω5 results in the modification of the shape of the potential energy curve of H-bond (3) in such a way that the lowest minimum of the double-well potential shifts from Im (I) to Im (III). As a consequence, the positive charge (the proton in H-bond (3)) displaces from Im (I) to Im (III). Similarly, to further drive the proton away from the H3PO4 dopant (from Im (III) to Im (V)), structure G4-[3](III) must be changed to structure G4-[3](V) by the rotation of ω3 from 87 to 101 degrees. It should be recognized that in a low local-dielectric environment, although the rotations of ω5 and ω3 modify the shapes of the potential energy curves, the energy barriers (ΔE) in the double-well potentials are not high enough to prohibit protons from returning to the original Im molecules, e.g., in H-bond (5), ΔE is only approximately 3.3 kJ mol−1, compared with the thermal energy fluctuation at room temperature of approximately 2.5 kJ mol−1. To visualize the roles played by the continuum solvent, the relative energy diagram and potential energy curves for proton displacement in structures C4-[3](I), C4-[3](III), C4-[3](V) and C4-[3](VII), which are shown in Fig. 5b, are discussed. It appears that in general, in a high local-dielectric environment (when COSMO is switched on), the energy barriers are substantially higher than in the gas phase. It is therefore reasonable to conclude that a high local-dielectric environment stabilizes the positive charge and prohibits protons from shuttling back to the original Im molecule.
image file: c4ra08198f-f5.tif
Fig. 5 Effects of Im H-bond chain structure reorganization on the energetics of proton transfer. Potential energy curves of proton displacement in H-bonds (3) and (5) are included. (a) Proton displacement in H+(H3PO4)(Im)4 in low local-dielectric environment (ε = 1). (b) Proton displacement in H+(H3PO4)(Im)4 in high local-dielectric environment (COSMO, ε = 23).

Dynamic results

BOMD simulations on the presolvation models reveal that in the continuum solvent (COSMO, ε = 23), a proton spends the majority of the time at the bottom of the double-well potential (a close-contact A–H⋯B H-bond structure) and there is no intermediate complex with the oscillatory shuttling motion (the shared-proton A⋯H⋯B H-bond structure) observed at any temperature. This is in accordance with the static results and confirms that the intermediate complexes are effectively formed in a low local-dielectric environment. Because the quasi-dynamic equilibrium between the oscillatory shuttling and structural diffusion motions has been shown in our previous studies to be a rate-determining process,24,26,44,45,49 the dynamics of proton dissociation and transfer are discussed based only on the BOMD results on the intermediate complexes in a low local-dielectric environment (ε = 1).

Dynamics of proton dissociation

To verify the elementary steps of proton dissociation and transfer suggested by the static calculations, the BOMD results on structures G2-[1](I,II), G3-[1](I,II) and G3-[2](III) at 353 K are analyzed in detail. Structures G2-[1](I,II) and G3-[2](III) are, according to νNH (νNH < νNH*), intermediate complexes for the proton transfer in the embedded and terminal presolvation models, respectively. Because H-bonds (1) and (2) in structure G2-[1](I,II) are equivalent, only the dynamics in H-bond (1) are discussed. The proton transfer profile (variations of the RN–O, RN–H and RO–H distances) in Fig. 6a is dominated by a close-contact N–H⋯O H-bond, with an occasional shared-proton N⋯H⋯O H-bond, examples of which are shown in panels P1 and P2. Due to the strong acid–base interaction, a stable Im+-H2PO4-Im+ structure is formed, evident from the probability of proton transfer from O to N, PON = 0.87, computed from the distribution of the N–H covalent bond distance (RN–H) in Fig. S3a.
image file: c4ra08198f-f6.tif
Fig. 6 (a and b) Variations of the RN–O and RN–H distances in H-bond (1) of structure G2-[1](I,II), obtained from BOMD simulations at 353 K. (c) Characteristic vibrational spectra of proton dissociation in H-bond (1) of structure G2-[1](I,II), obtained from BOMD simulations at 353 K. (d) Arrhenius plot of proton dissociation in H-bond (1) of structure G2-[1](I,II), obtained from BOMD simulations over the temperature range of 298 to 423 K. S and L = small- and large-amplitude N–O vibrations, respectively. ΔE‡,Arr = Arrhenius activation energy of proton dissociation.

Investigation of the proton transfer profile in panel P1 of Fig. 6b reveals two characteristic N–O vibrations associated with proton dissociation, namely small (S) and large (L) amplitude vibrations. In our previous study, these two characteristic modes were shown to underlie the oscillatory shuttling and structural diffusion motions of protons in H-bonds, respectively.49 For structure G2-[1](I,II), the Fourier transformation of VACF of the atoms in panel P1 (between 3500 and 4200 fs) yields the vibrational spectra in Fig. 6c in which the N–O vibrational frequency (νNO,MD) is outstanding at 894 cm−1, with strong interferences of torsional oscillations at frequencies below 500 cm−1 and bending motions in the range between 1000 and 1550 cm−1.30 According to the static results in Table S2, peak B at 2177 cm−1 and peak C at 2322 cm−1 are characteristic asymmetric N–H stretching frequencies (νNH,MD) of protons moving on the dissociation and association paths (see Fig. 3), respectively, whereas peak A at 1633 cm−1 is the vibrational signature of the shared-proton in the N⋯H⋯O H-bond, which is a proton moving in single-well potential. These transient vibrational frequencies can be used to monitor the acid–base proton dissociation in experiments. The Arrhenius plot in Fig. 6d yields the activation energy of proton dissociation (ΔE‡,Arr) in H-bond (1) in structure G2-[1](I,II) to be 8.7 kJ mol−1, which is considerably lower than that in an aqueous solution, ΔE‡,Arr = 13.1 kJ mol−1.19

To confirm the effects of the Im H-bond chain extension suggested by the static results, the proton transfer profiles, distributions of the N–H covalent bond distances (RN–H) and vibrational spectra of protons in H-bonds (1) and (2) of structures G2-[1](I,II) and G3-[1](I,II) are discussed. Extension of the Im H-bond chain at Im (II) of structure G2-[1](I,II) leads to structure G3-[1](I,II). The analyses of the proton transfer profiles in Fig. 6a, 7a and e and the distributions of RN–H in Fig. S3 confirm that H-bond (1) is stabilized and H-bond (2) is destabilized upon the Im H-bond chain extension. The former is evident from a significant increase in the probability of finding the shared-proton structure (RN–H between 1.3 and 1.5 Å) in H-bond (1), from PN⋯H⋯O = 0.10 in structure G2-[1](I,II) to 0.22 in structure G3-[1](I,II) (see Fig. S3a and S3b), whereas in H-bond (2), PN⋯H⋯O decreases from 0.10 to 0.01 (see Fig. S3a and S3c, where H-bond (1) and H-bond (2) in structure G2-[1](I,II) are equivalent). Because the shared-proton structure of H-bond (1) is more stable, the activation energy (ΔE‡,Arr) for proton dissociation in structure G3-[1](I,II) is higher; the Arrhenius plot in Fig. 7d yields ΔE‡,Arr = 10.8 kJ mol−1, compared with 8.7 kJ mol−1 in structure G2-[1](I,II). The stabilization effect of H-bond (1) in structure G3-[1](I,II) is also indicated by the red-shifts of νNH,MD of peaks A, B and C in Fig. 7c compared with those of structure G2-[1](I,II) in Fig. 6c. The destabilization effect of H-bond (2) is manifested by the very rare proton exchange event observed in Fig. 7e and the disappearance of the three characteristic peaks in Fig. 7g. Additionally, in Fig. 7c, as the relative intensity of peak A is higher than peaks B and C, the extension of the Im H-bond chain is confirmed to increase the probability of finding the oscillatory shuttling motion, leading to a longer lifetime of the rate-determining process in H-bond (1) and a higher ΔE‡,Arr in structure G3-[1](I,II) compared with structure G2-[1](I,II).


image file: c4ra08198f-f7.tif
Fig. 7 (a and b) Variations of the RN–O and RN–H distances in H-bond (1) of structure G3-[1](I,II), obtained from BOMD simulations at 353 K. (c) Characteristic vibrational spectra of proton dissociation in H-bond (1) of structure G3-[1](I,II), obtained from BOMD simulations at 353 K. (d) Arrhenius plot of proton dissociation in H-bond (1) of structure G3-[1](I,II), obtained from BOMD simulations over the temperature range of 298 to 423 K. (e and f) Variations of the RN–O and RN–H distances in H-bond (2) of structure G3-[1](I,II), obtained from BOMD simulations at 353 K. (g) Characteristic vibrational spectra of H-bond (2) of structure G3-[1](I,II), obtained from BOMD simulations at 353 K. S and L = small- and large-amplitude N–O vibrations, respectively. ΔE‡,Arr = Arrhenius activation energy of proton dissociation.

The effects of doping on the dynamics of proton transfer

The effects of doping are verified by comparing the proton transfer profiles and vibrational spectra of a proton in H-bond (3) of structure G3-[2](III) with those in H+(Im)4.30 The proton transfer profile in panel P1 of Fig. 8a and b is dominated by a large-amplitude N–N vibration in which RN–N varies over a very wide range, between 2.4 and 2.9 Å. It is shown in panel P1 that the proton exchange in H-bond (3) takes place almost exclusively through the structural diffusion mechanism, without the oscillatory shuttling motion as the rate-determining process.49 In contrast, in the undoped Im H-bond system (H+(Im)4), the proton transfer profile in Fig. S4a is characterized by both large- and small-amplitude N–N vibrations, labeled L and S, respectively, in panel P1 of Fig. S4b. Due to the interference between various modes of vibrations, especially below 2500 cm−1, the oscillatory shuttling and structural diffusion frequencies of a proton in H-bond (3) cannot be precisely determined from the vibrational spectra in Fig. 8c. However, based on the static results in Table S2, the two vibrational frequencies associated with the proton exchange in H-bond (3) can be estimated: peaks B and C at 1901 and 2345 cm−1, respectively. For H+(Im)4,30 the oscillatory shuttling and structural diffusion frequencies can be easily assigned from the vibrational spectra in Fig. S4c, for peaks A, B and C at 1195 cm−1, 2205 and 2508 cm−1, respectively. The red shifts of peaks B and C in structure G3-[2](III), compared with those in H+(Im)4, confirm that H3PO4 helps promote the structural diffusion process by strengthening the Im H-bond chain and enhancing the large-amplitude N–N vibration in the H3PO4 doped Im H-bond system. This results in a lower activation energy (ΔE‡,Arr) and a higher rate of proton transfer in H-bond (3) of the doped system; the Arrhenius plot in Fig. 8d yields ΔE‡,Arr of 1.22 kJ mol−1, compared with ΔE‡,Arr = 3.28 kJ mol−1 in H+(Im)4 (see Fig. S4d).30
image file: c4ra08198f-f8.tif
Fig. 8 (a and b) Variations of the RN–N and RN–H distances in H-bond (3) of structure G3-[2](III), obtained from BOMD simulations at 353 K. (c) Characteristic vibrational spectra of proton transfer in H-bond (3) of structure G3-[2](III), obtained from BOMD simulations at 353 K. (d) Arrhenius plot of proton transfer in H-bond (3) of structure G3-[2](III), obtained from BOMD simulations over the temperature range of 298 to 423 K. ΔE‡,Arr = Arrhenius activation energy of proton transfer.

It should be noted that ΔE‡,Arr is computed during the proton exchange process, e.g., in panel P1 of Fig. 8a and b, which involves the quasi-dynamic equilibrium between the oscillatory shuttling and structural diffusion motions. Therefore, ΔE‡,Arr might not be directly measurable in an experiment in which different experimental techniques generally yield different activation energies and where other effects could play important roles. As the proton transfer profile in Fig. 8a shows both the shared-proton and close-contact structures, the actual proton transfer is expected to take place not very often and a frequency factor must be included in the model calculations to obtain a more accurate first-order rate constant for the proton displacement in the Im H-bond chain. It should be added that although our presolvation models do not take into account a continuous H-bond network, because proton transfer in the Im H-bond system can be considered a local process with a short spatial, temporal correlation,27 it is reasonable to investigate only the rate-determining process (the dynamic behavior of the intermediate complex), which, in this case, is the quasi-dynamic equilibrium between the oscillatory shuttling and structural diffusion motions (the close-contact N–H⋯O and shared-proton N⋯H⋯O H-bonds).

Dynamics of the Im H-bond chain structure reorganization

Because the static results suggest that the shapes of the potential energy curves for proton displacement are sensitive to the Im H-bond chain structure reorganization (see Fig. 5), the effects of a conformation change on the dynamics of proton transfer are discussed using the proton transfer profiles of structure G3-[2](III), which were obtained from the BOMD simulations at 298 K. Variations of RN–N and RN–H in H-bonds (3) and (5), as well as the time evolutions of ω3 and ω5, as shown in Fig. 9, confirm that the Im H-bond chain structure reorganization does affect the dynamics of proton transfer. Investigation of the time evolutions of the torsional angles suggests division of the proton transfer profiles into four panels of 2500 fs each (P1 to P4).
image file: c4ra08198f-f9.tif
Fig. 9 Proton transfer profiles of H-bonds (3) and (5), and time evolutions of ω3 and ω5 in the terminal structure G3-[2](III), obtained from BOMD simulations at 298 K. (a–c) = H-bond (3); (d–f) = H-bond (5).

image file: c4ra08198f-f10.tif
Fig. 10 Conformations of structure G3-[2](III) sampled from panel P1 to P4 and projections of the reference vectors (vref), obtained from the analyses of the time evolutions of ω3 and ω5 in Fig. 9. Angles in degree; (III) = vref of Im (III); (V) = vref of Im (V).

To visualize the effects of the conformation change on the proton displacement in H-bond (3), a reference vector (vref) coinciding with the C2 axis of ImH+ is defined, and projections of vref, similar to the Newman projection, are used to monitor the dynamically changing Im H-bond structures. They are shown in Fig. 10. In this case, the line of sight is defined lengthwise down the Im H-bond chain from the H3PO4 dopant. Based on the assumption that the probability of finding a proton moving in a specific direction is correlated with the probability of the N–H covalent bond formation, the distributions of RN–H in H-bond (3) of Im (III) in panel P1 to P4 are constructed and shown in Fig. 11.


image file: c4ra08198f-f11.tif
Fig. 11 Distributions of the N–H covalent bond distances (RN–H) in H-bond (3) of structure G3-[2](III), obtained from BOMD simulations at 298 K. (a–d) Distributions in panel P1 to P4 in Fig. 9, respectively.

In the initial state in panel P1, ω3 and ω5 only fluctuate in a narrow range, with average values of 〈ω3〉 = 95 and 〈ω5〉 = 84 degrees. In panel P2, the variations in ω3 and ω5 become unsynchronized; while ω3 librates in a narrow range, regarded as a “rotational waiting state”, ω5 changes from 84 to 354 degrees, the latter of which can be regarded as a critical angle. In panel P3, while ω5 is maintained at the critical angle, ω3 changes to 282 degrees, in response to the rotation of ω5. Finally, in panel P4, both are maintained at 〈ω3〉 = 280 and 〈ω5〉 = 270 degrees. As the torsional rotation of ω5 induces the rotation of ω3, with a rotational waiting state of approximately 2500 fs, the rotational motions of ω5 and ω3 can be collectively regarded as “helical-rotational motion”, which underlies the Im H-bond chain structure reorganization. The Im H-bond chain structure reorganization affects the proton transfer in H-bond (3) in such a way that the population of the protons moving away from the dopant is increased from panel P1 (PNN = 0.3) to P2 (PNN = 0.5) and P3 (PNN = 0.8). This finding supports the static results in Fig. 5a in which the rotation of ω5 leads to a modification of the shape of the potential energy curve of H-bond (3) in such a way that the lowest minimum of the double-well potential shifts from Im (I) to Im (III), resulting in a higher probability of finding a proton moving towards Im (III), which confirms the effects of the Im H-bond chain structure reorganization suggested by the static results.

Conclusion

In this work, the dynamics and mechanisms of proton dissociation and transfer in an H3PO4 doped Im H-bond system were theoretically studied using B3LYP/TZVP calculations and BOMD simulations as model calculations and H+(H3PO4)(Im)n(n = 1–4) as model systems. The theoretical studies began with choosing appropriate presolvation models for proton dissociation and transfer, represented by embedded and terminal linear H-bond structures, respectively. The structures, energetic and vibrational frequencies of protons in the H-bonds of the presolvation models were studied in low (ε = 1) and high local-dielectric (ε = 23) environments. The static results confirmed that the excess proton condition is required to promote proton dissociation and transfer, and it must be included in the presolvation model. The 2D-PES and potential energy curves of the proton exchange suggested that the intermediate complexes for proton dissociation and transfer are formed in a low local-dielectric environment, whereas a high local-dielectric environment is required to stabilize the positive charge and prevent a proton from being shared or returning to the original Im molecule, which is in accordance with the localized charge solvent (LCS) concept. The B3LYP/TZVP calculations also revealed that the embedded H-bond structure with n = 2 represents the smallest, most active intermediate complex for proton dissociation, and similarly, the terminal H-bond structure with n = 3 is representative for proton transfer in the Im H-bond chain. The mechanisms of terminal structure formation and proton transfer in the Im H-bond chain, which involve fluctuations of the H-bond chain lengths and the local-dielectric environment, were proposed based on the static results. Additionally, because the shapes of the potential energy curves are sensitive to the torsional rotations of the H-bonds, the Im H-bond chain structure reorganization is included as an important process in the proton transfer mechanism. A comparison of the threshold frequencies for proton transfer in H+(H3PO4)(Im)n(n = 1–4) and in the undoped Im H-bond system (H+(Im)4) suggested that the strong polarization effects of the dopant lead to a higher probability of finding a delocalized proton and improve proton conduction in the H3PO4 doped Im system.

Because the static results showed that the embedded H-bond structure with n = 2 and the terminal H-bond structure with n = 3 represent the smallest, most active intermediate complexes for proton dissociation and transfer in a low local-dielectric constant environment, respectively, BOMD simulations were performed on these two H-bond structures. The vibrational spectra obtained from the BOMD simulations suggested three transient vibrational peaks associated with the acid–base proton dissociation, with the activation energy of the proton dissociation in the smallest, most active intermediate complex of 8.7 kJ mol−1. It appeared that the extension of the embedded H-bond chain increases the probability of finding the oscillatory shuttling motion and the lifetime of the intermediate complex, resulting in a higher activation energy for proton dissociation from the protonated acid. The comparison of the dynamic results on the H3PO4 doped Im system with those on H+(Im)4 indicated that the structural diffusion process in the Im H-bond chain of the doped system is almost exclusively through the large-amplitude N–N vibration, without the oscillatory shuttling motion as the rate-determining process, with an activation energy of only 1.22 kJ mol−1. For the undoped Im H-bond system, the proton transfer process is characterized by both the oscillatory shuttling and structural diffusion motions. Analyses of the proton transfer profiles and the time evolutions of the torsional angles revealed coupled torsional motions (helical-rotational motion), which drive the proton away from the dopant and confirm that the Im H-bond chain structure reorganization plays an important role in the proton transfer in the H3PO4 doped Im H-bond system. The results showed for the first time the interplay among the translational, rotational and vibrational motions that underlie the dynamics and mechanisms of proton transfer in the linear H-bond chain. Because the discussion are based on the results on carefully selected presolvation models, the present theoretical study provides new insights into the dynamics and mechanisms of proton exchange processes in doped, weak-polar organic solvent. These pieces of fundamental information are important and can be used as guidelines for further theoretical and experimental investigations on proton transfer processes in aromatic heterocyclic materials, which could eventually lead to an improvement of the proton transfer processes in non-aqueous, high temperature fuel cells.

Acknowledgements

The authors would like to acknowledge the financial supports provided by Suranaree University of Technology through the SUT-Ph.D. program (SUT-Ph.D./05/2554). This work was supported in part by the Higher Education Research Promotion and National Research University Project (HERB-NRU), the Office of the Higher Education Commission (OHEC), Thailand. Prof. Kritsana Sagarik received financial supports from the Thailand Research Fund (TRF), the Advanced Research Scholarship (BRG-5580011). Super- and high-performance computer facilities provided by the following organizations are gratefully acknowledged; School of Mathematics and School of Chemistry, SUT; National Electronics and Computer Technology Center (NECTEC), National Science and Technology Development Agency (NSTDA); the Thai National Grid Center (THAIGRID), Ministry of Information and Communication Technology (MICT).

References

  1. T. Koppel, Powering the Future: The ballard fuel cell and the race to change the world, John Wiley & Sons Ltd., New York, 1999 Search PubMed.
  2. J. Larminie and A. Dicks, Fuel Cell Systems, John Wiley & Sons Ltd., Chichester, 2001 Search PubMed.
  3. S. R. Narayanan, S. P. Yen, L. Liu and S. G. Greenbaum, J. Phys. Chem. B, 2006, 110, 3942–3948 CrossRef CAS PubMed.
  4. M. F. H. Schuster and W. Meyer, Annu. Rev. Mater. Res., 2003, 33, 233–261 CrossRef CAS.
  5. H. A. Pohl and R. P. Chartoff, J. Polym. Sci., Polym. Chem. Ed., 1964, 2, 2787–2806 Search PubMed.
  6. X. Glipa, B. Bonnet, B. Mula, D. J. Jones and J. RozieÁre, J. Mater. Chem., 1999, 9, 3045–3049 RSC.
  7. B. Xing and O. Savadogo, J. New Mater. Electrochem. Syst., 1999, 2, 95–101 CAS.
  8. R. He, Q. Li, G. Xiao and N. J. Bjerrum, J. Membr. Sci., 2003, 226, 169–184 CrossRef CAS PubMed.
  9. D. T. Chin and H. H. Chang, J. Appl. Electrochem., 1989, 19, 95–99 CrossRef CAS.
  10. R. A. Munson and M. E. Lazarus, J. Phys. Chem., 1967, 71, 3245–3248 CrossRef CAS.
  11. B. S. Hickman, M. Mascal, J. J. Titman and I. G. Wood, J. Am. Chem. Soc., 1999, 121, 11486–11490 CrossRef CAS.
  12. A. Kawada, A. R. McGhie and M. M. Labes, J. Chem. Phys., 1970, 52, 3121–3125 CrossRef CAS PubMed.
  13. I. Fischbach, H. W. Spiess, K. Saalwächter and G. R. Goward, J. Phys. Chem. B, 2004, 108, 18500–18508 CrossRef CAS.
  14. M. Iannuzzi and M. Parrinello, Phys. Rev. Lett., 2004, 93, 025901 CrossRef CAS.
  15. K. D. Kreuer, A. Fuchs, M. Ise, M. Spaeth and J. Maier, Electrochim. Acta, 1998, 43, 1281–1288 CrossRef CAS.
  16. W. Münch, K. D. Kreuer, W. Silvestri, J. Maier and G. Seifert, Solid State Ionics, 2001, 145, 437–443 CrossRef.
  17. A. Schechter and R. F. Savinell, Solid State Ionics, 2002, 147, 181–187 CrossRef CAS.
  18. L. Vilciauskas, M. E. Tuckerman, J. P. Melchior, G. Bester and K.-D. Kreuer, Solid State Ionics, 2013, 252, 34–39 CrossRef CAS PubMed.
  19. P. Suwannakum and K. Sagarik, submitted.
  20. T. C. Berkelbach, H.-S. Lee and M. E. Tuckerman, Phys. Rev. Lett., 2009, 103, 238302 CrossRef.
  21. D. Marx, A. Chandra and M. E. Tuckerman, Chem. Rev., 2010, 110, 2174–2216 CrossRef CAS PubMed.
  22. C. Lao-ngam, P. Asawakun, S. Wannarat and K. Sagarik, Phys. Chem. Chem. Phys., 2011, 13, 4562 RSC.
  23. C. Lao-ngam, M. Phonyiem, S. Chaiwongwattana, Y. Kawazoe and K. Sagarik, Chem. Phys., 2013, 420, 50–61 CrossRef CAS PubMed.
  24. M. Phonyiem, S. Chaiwongwattana, C. Lao-ngam and K. Sagarik, Phys. Chem. Chem. Phys., 2011, 13, 10923 RSC.
  25. S. Chaiwongwattana, M. Phonyiem, V. Vchirawongkwin, S. Prueksaaroon and K. Sagarik, J. Comput. Chem., 2012, 33, 175–188 CrossRef CAS PubMed.
  26. K. Sagarik, S. Chaiwongwattana, V. Vchirawongkwin and S. Prueksaaroon, Phys. Chem. Chem. Phys., 2010, 12, 918 RSC.
  27. A. Li, Z. Cao, Y. Li, T. Yan and P. Shen, J. Phys. Chem. B, 2012, 116, 12793–12800 CrossRef CAS PubMed.
  28. M. Iannuzzi, J. Chem. Phys., 2006, 124, 204710 CrossRef PubMed.
  29. P. Hobza, H. L. Selzle and E. W. Schlag, Chem. Rev., 1994, 94, 1767 CrossRef CAS.
  30. V. Buagern and K. Sagarik, in preparation.
  31. N. N. Greenwood and A. Thompson, J. Chem. Soc., 1959, 3485 RSC.
  32. H. Chen, T. Yan and G. A. Voth, J. Phys. Chem. A, 2009, 113, 4507–4517 CrossRef CAS PubMed.
  33. C. Lao-ngam, P. Asawakun, S. Wannarat and K. Sagarik, Phys. Chem. Chem. Phys., 2011, 13, 4562–4575 RSC.
  34. C. Lao-ngam, M. Phonyiem, S. Chaiwongwattana, Y. Kawazoe and K. Sagarik, Chem. Phys., 2013, 420, 50–61 CrossRef CAS PubMed.
  35. K. P. Sagarik and R. Ahlrichs, J. Chem. Phys., 1987, 86, 5117 CrossRef CAS PubMed.
  36. K. P. Sagarik and P. Asawakun, Chem. Phys., 1997, 219, 173 CrossRef CAS.
  37. K. P. Sagarik, S. Chaiwongwattana and P. Sisot, Chem. Phys., 2004, 1 CrossRef CAS PubMed.
  38. K. P. Sagarik and S. Dokmaisrijan, J. Mol. Struct., 2005, 718, 31 CrossRef CAS PubMed.
  39. K. P. Sagarik, V. Pongpituk, S. Chiyapongs and P. Sisot, Chem. Phys., 1991, 156, 439 CrossRef CAS.
  40. K. Kwac, C. Lee, Y. Jung and J. Han, J. Chem. Phys., 2006, 125, 244508 CrossRef PubMed.
  41. J. L. Knee, L. R. Knundkar and A. H. Zewail, J. Chem. Phys., 1985, 82, 4715 CrossRef CAS PubMed.
  42. S. J. Paddison, J. New Mater. Electrochem. Syst., 2001, 4, 197–207 CAS.
  43. A. P. Scott and L. Radom, J. Phys. Chem., 1996, 100, 16502 CrossRef CAS.
  44. S. Chaiwongwattana, M. Phonyiem, V. Vchirawongkwin, S. Prueksaaroon and K. Sagarik, J. Comput. Chem., 2012, 33, 175–188 CrossRef CAS PubMed.
  45. C. Lao-ngam, P. Asawakun, S. Wannarat and K. Sagarik, Phys. Chem. Chem. Phys., 2011, 13, 4562 RSC.
  46. SURFER for Window, Golden Software Inc., USA, 6.04 edn, 1997 Search PubMed.
  47. J. E. Basconi and M. R. Shirts, J. Chem. Theory Comput., 2013, 9, 2887–2899 CrossRef CAS.
  48. S. Chaiwongwattana, C. Lao-ngam, M. Phonyiem, V. Vchirawongkwin and K. Sagarik, submitted.
  49. K. Sagarik, M. Phonyiem, C. Lao-Ngam and S. Chaiwongwattana, Phys. Chem. Chem. Phys., 2008, 10, 2098 RSC.
  50. R. L. Davidchack, J. Comput. Phys., 2010, 229, 9323–9346 CrossRef CAS PubMed.
  51. B. Koeppe, J. Guo, P. M. Tolstoy, G. S. Denisov and H.-H. Limbach, J. Am. Chem. Soc., 2013, 135, 7553–7566 CrossRef CAS PubMed.
  52. R. R. Sadeghi and H.-P. Cheng, J. Chem. Phys., 1999, 111, 2086 CrossRef CAS PubMed.
  53. R. O. Watts and I. J. McGee, Liquid State Chemical Physics, John Wiley & Sons, New York, 1976 Search PubMed.
  54. M. Schmollngruber, C. Schröde and O. Steinhauser, Phys. Chem. Chem. Phys., 2014, 16, 10999–11009 RSC.

Footnote

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

This journal is © The Royal Society of Chemistry 2014
Click here to see how this site uses Cookies. View our privacy policy here.