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
First published on 11th November 2014
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.
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:
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.
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.
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.
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
The behaviors of proton dissociation and transfer were primarily evaluated from the asymmetric stretching coordinate (ΔdDA), computed from ΔdDA = |dA–H − dB⋯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.
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 10000 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
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.
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.
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.
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).
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).
![]() | ||
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.
![]() | ||
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.
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.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4ra08198f |
This journal is © The Royal Society of Chemistry 2014 |