The Grotthuss mechanism for bifunctional proton transfer in poly(benzimidazole)

Poly(benzimidazole) (PBI) has received considerable attention as an effective high-temperature polymer electrolyte membrane for fuel cells. In this work, the Grotthuss mechanism for bifunctional proton transfer in PBI membranes was studied using density functional theory and transition state theory. This study focused on the reaction paths and kinetics for bifunctional proton transfer scenarios in neutral ([PBI]2), single (H+[PBI]2) and double-protonated (H2+[PBI]2) dimers. The theoretical results showed that the energy barriers and strength for H-bonds are sensitive to the local dielectric environment. For [PBI]2 with ε = 1, the uphill potential energy curve is attributed to extraordinarily strong ion-pair H-bonds in the transition structure, regarded as a ‘dipolar energy trap’. For ε = 23, the ion-pair charges are partially neutralized, leading to a reduction in the electrostatic attraction in the transition structure. The dipolar energy trap appears to prohibit interconversion between the precursor, transition and proton-transferred structures, which rules out the possibility for [PBI]2 to be involved in the Grotthuss mechanism. For H+[PBI]2 and H2+[PBI]2 with ε = 1, the interconversion involves a low energy barrier, and the increase in the energy barrier for ε = 23 can be attributed to an increase in the strength of the protonated H-bonds in the transition structure: the local dielectric environment enhances the donor–acceptor interaction of the protonated H-bonds. Analysis of the rate constants confirmed that the quantum effect is not negligible for the N–H+ … N H-bond especially at low temperatures. Agreement between the theoretical and experimental data leads to the conclusion that the concerted bifunctional proton transfer in H2+[PBI]2 in a high local dielectric environment is ‘the rate-determining scenario’. Therefore, a low local dielectric environment can be one of the required conditions for effective proton conduction in acid-doped PBI membranes. These theoretical results provide insights into the Grotthuss mechanism, which can be used as guidelines for understanding the fundamentals of proton transfers in other bifunctional H-bond systems.


Introduction
Fuel cells have been accepted as environmentally friendly electrochemical devices that can effectively transform chemical energy into applicable electrical energy. However, because the most commonly used polymer electrolyte membrane (PEM) in fuel cells is hydrated Nafion ® , the operating temperature has been a serious limitation; water-based fuel cells cannot be operated at and above the boiling point of water. As an aromatic heterocyclic polymer with extraordinary thermal stability and high proton conductivity at relatively high temperatures (400-600 K), poly(benzimidazole) (PBI) has received considerable attention as an effective PEM in fuel cells [1]. Although pristine PBI is not a proton-conducting polymer (σ = 10 −12 S cm −1 ) and cannot be directly used as a PEM [2], acid-doped PBI membranes are excellent proton-conducting polymers [3]. Because the chemistry of acid-doped PBI membranes has been discussed in detail in several review articles [4], only information relevant to the present study is briefly summarized.
Based on the Scatchard method [5], both imide groups in PBI (figure 1) were reported to be preferentially protonated with a considerably high protonation constant compared with sites with lower affinity. The room temperature acid dissociation constants (K a ) for these two protonated sites are 5.4 × 10 −4 and 3.6 × 10 −2 , respectively. It was found that when both imide groups are protonated, proton transfers take place through hydrogen bonds (H-bonds) between the protonated and nonprotonated imino nitrogen groups (the N-H + … N H-bond) [6]. Experiments revealed that in anhydrous states in the temperature range between 298 and 433 K, proton conductivities for PBI membranes with acid doping levels of two range from 10 −9 to 10 −5 S cm −1 [7].
Acid-doped PBI membranes can be prepared directly from H 3 PO 4 and H 2 SO 4 solutions [8], for which the conductivity depends upon the acid concentration, temperature and duration of immersion [8]. The H 3 PO 4 -doped PBI membrane exhibits excellent conductivity. For example, at 403 K, σ ranges between 5.0 × 10 −3 and 2.0 × 10 −2 S cm −1 and is increased to 3.5 × 10 −2 S cm −1 at 463 K [3]. In the high conductivity states, the H 3 PO 4 -doped PBI membrane was observed to be very swollen, soft and flexible with a decreased mechanical strength. Analysis of the conductivity in the temperature range of 298 to 403 K revealed that the conductivity of H 3 PO 4 -and H 2 SO 4 -doped PBI membranes can be explained by using the Arrhenius equation [9], for which the maximum degree of protonation was concluded to be two and one for the acids, respectively, implying that the imide groups in the benzimidazole ring can be fully protonated by H 3 PO 4 . It was shown that the high proton conductivity in acid-doped PBI membranes can be explained using the Grotthuss mechanism, not the weak electrolyte theory [9], and σ is approximately 10 −7 S cm −1 at 303 K [9].
Fontanella et al. [10], based on high-pressure conductivity measurements on 600 mol% H 3 PO 4 -doped PBI membrane (mediated by the H-bond networks between polymer chains) and 85% H 3 PO 4 solution (mediated by the H-bond networks between H 3 PO 4 ), suggested that at room temperature, the conductivity decreases with increasing pressure, resulting from an increased viscosity. However, as the temperature is increased, the acid-doped PBI membranes behave like a polymer electrolyte, in which ion transport is enhanced by polymer segmental motions, and, at 348 K, the conductivity increases with increasing pressure. These results indicate that the ion transport mechanisms in a 600 mol% H 3 PO 4 -doped PBI membrane and 85% H 3 PO 4 solution are different, especially at high temperature and pressure.
Temperature-and pressure-dependent dielectric spectroscopic experiments also revealed that the PBI molecule is fully protonated after blending with H 3 PO 4 [11]. The results showed that the temperature dependence of the conductivity can be described using the Arrhenius equation, for which the activation energy does not change with H 3 PO 4 concentration. Based on the variation in activation volume with temperature and acid concentration, proton transfer between phosphate and imidazole moieties through H-bonds and the self-diffusion of phosphate are anticipated to be the two most important mechanisms. This rules out the proton transport mechanism through the polymer segmental motion proposed in [10].
In our previous study [12], proton transfer in a unit cell of the imidazole (Im) crystal structure (H + [Im] n , n = 2-4) was studied using the density functional theory (DFT) method with the Becke, 3-Parameter, Lee-Yang-Parr (B3LYP) hybrid functional and TZVP basis set. The B3LYP/TZVP results revealed that H + [Im] 2 is the smallest, most active Zundel-like complex. Potential energy curves for proton transfer in N-H + … N H-bonds indicated that a single-well potential is favourable in a low local dielectric environment (ε = 1), whereas in a high local dielectric environment (ε = 23), a double-well potential dominates, and fluctuation of the local dielectric environment helps promote the Grotthuss mechanism through the Eigen-Zundel-Eigen-like scenario. Born-Oppenheimer molecular dynamics (BOMD) simulations confirmed that the interconversion between single-and double-well potentials (the Eigen-Zundel-Eigen-like scenario) results from the fluctuation of the number of Im molecules in the H-bond chain, and the rate-determining process for proton transfer is a local (short-range) process in which the N-N vibration, oscillatory shuttling motion of the H-bond proton and librational motion of the Im H-bond chain are coherent. Our theoretical studies of the Grotthuss mechanism also revealed that the dynamics of proton transfer can be affected by the motions of aromatic rings, e.g. in phosphonic acid-functionalized polymer ( poly(vinyl-phosphonic acid)) [13] and H 3 PO 4 -doped imidazole (Im) systems [14]; these motions can be compared with the polymer segmental motion in [10].
According to the above literature review, the following remarks can be made: (i) Experiments studied the temperature dependence of the conductivity of acid-doped PBI membranes using the Arrhenius equation and suggested that proton hopping in acid-acid (e.g. -P-O-H + … O=P-) and acid-PBI Hbond (e.g. N-H + … O) networks dominate especially below 373 K, whereas the information for 'bifunctional proton transfer' in the PBI-PBI (N-H + … N) H-bond chain is restricted; the Nheterocycles in PBI can act as a proton solvent and contribute to proton conduction in the PBI membrane [11]. (ii) There has been controversy as to whether or not proton transfer in a PBI royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 211168 membrane can be described by the Grotthuss mechanism, for which the Arrhenius equation can be used. (iii) Because the size of the PBI systems is large in view of quantum chemical methods, theoretical studies of proton transfer in a PBI membrane are restricted and the quantum effect has never been previously taken into account; the transportation of small particles such as H + has been shown in many cases to be governed by the quantum effect ( proton tunnelling), for which the rate constants are increased at low temperatures and the tunnelling of protons in the Im H-bond chain along with the c crystallographic direction is anticipated based on solid-state 15 N NMR spectroscopy [15].
In this study, the Grotthuss mechanism for bifunctional proton transfer in H-bonds between PBI molecules was investigated using the DFT method with the DZP and TZP basis sets. Based on our previous theoretical studies and reported experimental data, neutral ([PBI] 2 ), single (H + [PBI] 2 ) and double-protonated (H 2+ [PBI] 2 ) dimers were chosen as model systems, for which the effect of the size of the basis set and basis set superposition error (BSSE) was initially evaluated. The equilibrium and transition structures and energetics associated with the optimized proton transfer paths in low and high local dielectric environments were computed and analysed in detail; ε = 1 and ε = 23 were included in the model calculations to simulate the lowest and highest effects of the electric field, respectively, for the surrounding PBI molecules on the bifunctional proton transfer process. The kinetics and proton conductivities for the PBI membrane were studied over the temperature range of 200-500 K using transition state theory (TST), and the quantum effect is discussed based on the rate constants obtained with second-order and full Wigner corrections.

Quantum chemical calculations
Our theoretical studies of the dynamics and mechanisms for proton transfer showed that the DFT method with the B3LYP hybrid functionals was applicable to heterocyclic aromatic compounds [12]. This is due to the ability to approximate H-bond interaction energies with reasonable accuracy and computational resources. It should be noted that the PBI system is a π system, in which the π−π interaction could be important. However, because the dipole moments of the neutral and protonated monomers are relatively high, the dipole-dipole interaction is mainly responsible for the bifunctional H-bond interactions in the coplanar dimers; for example, for the neutral monomer, the dipole moments (μ) in ε = 1 and 23 are 4.1 and 5.9 D, respectively, and those for the single-protonated monomer are 16.3 and 19.8 D, respectively. In this work, the monomer of the commercially available PBI membrane, [2,2 0 -(m-phenylene)-5,5 0bibenzimidazole], shown in figure 1 [16], was used as a model molecule. Because the PBI systems studied in this work are large, the B3LYP method was applied primarily with the DZP basis set (abbreviated B3LYP/DZP), for which the effect of the size of the basis set was systematically investigated using the B3LYP/TZP method and the counterpoise correction of BSSE.

Equilibrium structures and reaction path optimizations
Our previous studies revealed that for protonated heterocyclic aromatic systems such as Im [12] and H 3 PO 4doped Im [14], the elementary processes of proton transfer can be studied reasonably well using presolvation models, in which only important precursor, transition state and proton transferred molecules are taken into account [17]. For the protonated H-bond in the Im system [12], our results also showed that the ratedetermining process for proton transfer is characterized by shuttling of the H-bond proton with a vibrational frequency strongly redshifted compared with the neutral H-bond. The theoretical results suggested that for the Im and H 3 PO 4 -doped Im systems, the H-bonds susceptible to proton transfer are strong with extraordinarily short distances, and the smallest, most-active intermediates are protonated dimers.
In this work, because experiments showed that proton transfer in PBI membranes occurs effectively under acidic conditions [2] and theories suggest that the fluctuation of the local dielectric environment can affect the proton transfer potential energy curves, [PBI] 2 , H + [PBI] 2 and H 2+ [PBI] 2 in ε = 1 and ε = 23 were chosen as model systems; the thermal energy agitations generally lead to fluctuations of the dipole orientation and local dielectric environment in the vicinity of molecules. The dielectric royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 211168 constants in the gas phase and bulk were used to simulate (approximate) low and high local dielectric environments, because the most populated value at/in the vicinity of the bifunctional H-bonds and the range over which it fluctuates at a particular temperature are not known.
The equilibrium geometries for the PBI dimers were optimized using the B3LYP/DZP method, and the conductor-like screening model (COSMO) [18] was used to simulate the effect of the local dielectric environment. The strength of the PBI-local dielectric environment interaction was approximated using the solvation energy (ΔE Solv ), computed from the difference between the total energies (E Total ) obtained with and without COSMO. In this study, emphasis was placed on H 2+ [PBI] 2 because experiments suggested that both imine groups in the Im rings (figure 1) play important roles in proton transfer in acid-doped PBI membranes, regarded as bifunctional proton transfer [19].
The H-bond interaction energies (ΔE H-bond ) for ε = 1 and ε = 23 were computed from ΔE H-bond = E N-H…N -(E N-H + E N ), where E N-H…N is the total energy of the N-H … N H-bond system and E N-H and E N are the total energies of the proton donor and acceptor moieties inside the supermolecule, respectively. To study the effect of BSSE, counterpoise correction [20] was applied, for which ΔE H-bond/CP = E N-H…N -(E N-H(A) + E N(D) ); (A) and (D) denote the total energies computed with the 'ghost' basis sets (without electrons and nuclei) for the proton acceptor and donor moieties, respectively. The likelihood of proton transfer in H-bonds was anticipated using the asymmetric stretching coordinate (Δd DA ), which was computed from Δd DA = |R N-H -R N…H |; R N-H and R N … H are the N-H and N … H distances in the N-H … N H-bond, respectively. In this study, the H-bond susceptible to proton transfer is characterized by Δd DA close to zero; the H-bond proton located at/or close to the centre of the H-bond is generally governed by a barrierless potential energy [21]. All of the DFT calculations were performed using the TURBOMOLE 7.50 software package [22].
Based on the equilibrium structures and energetics of the PBI dimers, six proton transfer scenarios were hypothesized for ε = 1 and ε = 23. The hypothesized paths were optimized using the nudged elastic band (NEB) method with the Limited-Memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimizer included in the ChemShell software package [23]. In the reaction path optimizations, 15 images connecting the equilibrium structures of the precursor and product structures were optimized. The relative energies with respect to the precursor (ΔE Rel ), ΔE H-bond and ΔE H-bond/CP along with the optimized proton transfer paths were plotted against the mass-weighted Cartesians.

Kinetics and proton conductivity calculations
Based on harmonic TST [24,25], the precursor, transition state and product structures obtained from the reaction path optimizations were used in the study of the kinetics of the bifunctional proton transfer process. Because the six proton transfer scenarios hypothesized in this work involve covalent bond dissociation and formation only inside the PBI dimers, unimolecular rate constants were calculated. The classical (k Class ) and quantized-vibrational (k Q-vib ) rate constants were computed over the temperature range of 200-500 K; these temperatures are lower than the glass transition temperature of PBI (713 K) [26]. k Class and k Q-vib were computed to study the possibility to use these semiclassical rate constants in the discussion of proton transfer in H-bond systems. For a one-dimensional energy profile, k Class is given by [27] where Q ‡ and Q R are the partition functions of the transition state and precursor equilibrium structures obtained from the NEB method and geometry optimizations, respectively, and ΔE ‡ is the energy barrier. k B and h are the Boltzmann and Planck constants, respectively. For k Q-vib , the energy barrier obtained with the zero-point vibrational energy (DE z ZPE ) is used Q z ZPE and Q R ZPE are the partition functions of the transition state and precursor equilibrium structures obtained with the zero-point vibrational energies, respectively.
To approximate the effect of quantum mechanical tunnelling, the crossover temperature (T c ) (the temperature below which the transition state is dominated by quantum mechanical tunnelling) was computed [28,29] royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 211168 V z is the imaginary frequency of the transition structure. To assess the effect of quantum mechanical tunnelling, the rate constants with quantized vibrations and second-order Wigner correction (k S-Wig ) [28,29] were computed based on the assumption that proton tunnelling occurs at the top of the barrier. The Wigner correction to the rate constant was carried out using the Wigner transmission coefficient (κ S-Wig ) ð2:4Þ and the Wigner-corrected rate constant is given by is equal to 1 at the classical limit (h = 0). In addition, the rate constants with full Wigner correction (k F-Wig ) [28,29] were also computed to assess all of the rate constants calculated in this work; although the quantum effect (the proton tunnelling), which generally leads to reduction of the energy barrier, could be studied more accurately using the quantum instanton method, our preliminary reaction path optimizations using this method required extensive computational resources and therefore are not applicable for the present PBI system. All of the kinetics and thermodynamics calculations were performed using the DL-FIND program [30] included in the ChemShell package [23]. To correlate the theoretical results obtained using the present model systems with the experimental data [6], proton conductivities (σ) were computed using the Arrhenius equation where A, T and ΔE ‡ are the pre-exponential term, temperature and energy barrier, respectively. σ was computed only for the rate-determining processes obeying the Arrhenius equation.

Results and discussion
The B3LYP/DZP and B3LYP/TZP methods suggested that the equilibrium structures of the [PBI], H + [PBI] and H 2+ [PBI] monomers take on approximately the same planar structure (C-form in figure 1), and the structures of the dimers formed from these monomers are also an approximately planar structure (electronic supplementary material, table S1); the high proton conductivity of the PBI membrane was attributed to the linear PBI chains with a high tendency for coplanar aromatic rings [19]. Therefore, the discussion, henceforth, is focused on these C-form dimers. In this study, the Grotthuss mechanism was hypothesized to occur through the bifunctional proton transfer process in the coplanar H-bond networks, e.g. in H

The DFT methods and equilibrium structures
The total (E Total ) and H-bond interaction energies (ΔE H-bond and ΔE H-bond/CP ) obtained from B3LYP/DZP and B3LYP/TZP calculations are listed in the electronic supplementary material, table S1, together with the H-bond distances (R N-N ) and asymmetric stretching coordinates (Δd DA ). Both B3LYP/DZP and B3LYP/TZP geometry optimizations yield the same dimer structures. For ε = 1, the R N-N distances obtained from the B3LYP/DZP method are slightly shorter with stronger H-bonds (more negative ΔE H-bond ). Although the stability of CG-[1] 2+ is higher (more negative E Total ) than CG-[1] + and CG- 2) kJ mol −1 . It appears that with the counterpoise correction, the discrepancy between the Hbond interaction energies obtained from the two DFT methods is significantly decreased (as already observed for ε = 1), e.g. for CC-[1] 2+ , ΔΔE H-bond/CP = 4.7 kJ mol −1 . Therefore, one can conclude that ΔE H-bond/CP obtained based on the B3LYP/DZP and B3LYP/TZP methods with the counterpoise correction are approximately the same, and it is reasonable to use the B3LYP/DZP results with the counterpoise correction in further discussion.
Because the local dielectric environment affects charges on atoms involved in H-bonds, to study this effect on the strength of H-bonds, the net stabilization energies (ΔE NSE ), defined as the difference between ΔE H-bond/CP for ε = 23 and 1, were calculated, and are given in the electronic supplementary material, To complete the evaluation of the DFT methods, the optimized proton transfer paths for bifunctional proton transfer in CG- [1] 0 are discussed as an example. Because the structure of the potential energy curve governs the kinetics of proton transfer, the relative total energies (ΔE Rel ) obtained based on the B3LYP/DZP and B3LYP/ TZP methods are compared in the electronic supplementary material, figure S1. Additionally, the relative H-bond interaction energies with respect to the precursor CG-[1] 0 (ΔE Rel,H-bond and ΔE Rel,H-bond/CP ) are included in the electronic supplementary material, figure S1; ΔE H-bond and ΔE H-bond/CP for CG-[1] 0 obtained from the B3LYP/DZP and B3LYP/TZP methods are set to 0 kJ mol −1 . It appears that the ΔE Rel values obtained from both methods are approximately identical and ΔE Rel,H-bond and ΔE Rel,H-bond/CP are virtually the same. These results confirm the applicability of the B3LYP/DZP method with the counterpoise correction, and the discussions, henceforth, are based only on these methods.  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 211168 For ε = 23, a similar potential energy curve is shown in figure 2b, in which the same non-concerted, uphill bifunctional proton transfer with lower energy barriers and weaker H-bonds are observed; ΔE ‡ , ΔE H-bond/CP and ΔE Solv for the transition structure CC- [2] 0, ‡ are 63.3, −232.1 and −140.8 kJ mol −1 , respectively. The uphill bifunctional proton transfer for ε = 1 and 23 can be attributed to the strength of the ion-pair associations in the transition structures (CG- [2] 0, ‡ and CC- [2] 0, ‡ ), which can be regarded as 'dipolar energy traps' [32]; because the ion-pair association in CC- [2] 0, ‡ is strongly destabilized by the electric field of the local dielectric environment (ε = 23), ΔE NSE = 324.6 kJ mol −1 (electronic supplementary material, table S2), the energy barrier is lower than that for ε = 1. Because the Grotthuss mechanism is characterized by the interconversion between precursor, transition state and proton transferred structures (e.g. the Eigen-Zundel-Eigen-like scenario) and bifunctional proton transfer in [PBI] 2 is governed by an uphill potential energy curve, the ion-pair transition structure of [PBI] 2 can be excluded from the Grotthuss mechanism for the PBI membrane.

Potential energy curves for bifunctional proton transfer
For H + [PBI] 2 with ε = 1, the potential energy curve in figure 3a reveals that the interconversion between the precursor, transition state and proton transferred structures (CG- [1] involves a low energy barrier; single-proton transfer in the forward direction is slightly more preferential (almost barrierless) than in the reverse direction, ΔE ‡ = 1.8 and 12.8 kJ mol −1 , respectively. Due to the higher aromaticity in H + [PBI] 2 , ΔE ‡ is lower than that in H + [Im] 2 , ΔE ‡ = 9.5 kJ mol −1 [12]. In conclusion, the role played by the local dielectric environment and the interplay among energies for the bifunctional proton transfer paths, the relative solvation (ΔE Rel,Solv ), total (ΔE Rel,Total ) and H-bond interaction (ΔE Rel,H-bond/CP ) energies with respect to the precursors were computed and are included in figure 5. The trends for ΔE Rel,Solv and ΔE Rel,Total in figure 5a suggest that for the protonated H-bond dimers (H 2+ [PBI] 2 and H + [PBI] 2 ), the increase in ΔE ‡ correlates with the increase in ΔE Rel,Solv for the transition structure; the transition structure is destabilized compared with the precursors and the higher the number of protonated H-bonds, the higher the ΔE ‡ for ε = 23. By contrast, for [PBI] 2 , the increase in ΔE ‡ is associated with the decrease in ΔE Rel,Solv ; the positive and negative charges in the ion-pair transition structure are both increased (stabilized) compared with the precursor, and the lower ΔE Rel,Solv is, the higher ΔE ‡ is for ε = 23. The latter prohibits the ion-pair transition structure from being involved in the Grotthuss mechanism.
The trends for ΔE Rel,Solv and ΔE Rel,H-bond/CP in figure 5b yield the information in the same direction. For H 2+ [PBI] 2 and H + [PBI] 2 , the increase in ΔE Rel,Solv for the transition structure is correlated with the strength of the H-bonds; the stronger the H-bonds with respect to the precursor, the higher the ΔE ‡ . This finding is in accordance with all of our previous studies, in which the rate-determining process for proton transfer is characterized by strong, protonated H-bonds with the oscillatory shuttling motion of the shared proton. For [PBI] 2 with ε = 23, the increase in the positive and negative charges along with the proton transfer path leads to strong electrostatic attraction in the ion-pair transition structure, which prohibits the Grotthuss mechanism; the Grotthuss mechanism takes place in dynamic H-bond networks, in which the formation and cleavage of covalent bonds are prerequisites. Based on the above discussion, one can anticipate that the rate-determining scenario for bifunctional proton transfer in the PBI membrane is most likely to be the concerted proton exchange in H 2+ [PBI] 2 , CC- 2+ . This piece of information was used as a guideline for calculation of the kinetics and proton conductivities.         [12] revealed that the intensity of the 1 H NMR chemical shift associated with the shared proton (24 ppm) increases as the temperature increases from 298 to 500 K and the oscillatory shuttling motion of the shared proton possesses lower vibrational frequency than the asymmetric N-H + … N H-bond.
The values for the rate constants in table 1 indicate that the quantum effect is not negligible, except for [PBI] 2 with ε = 1 and 23 and H 2+ [PBI] 2 with ε = 23, figure 6a and b, respectively, for which all of the rate constants obey the Arrhenius equation. These bifunctional proton transfer scenarios will be discussed in comparison with available experimental data. In table 1, the classical rate constants (k Class ) are the lowest, as expected, whereas the rate constants with quantized vibrations (k Q-vib ) and with the second-order Wigner correction (k S-Wig ) are similar, especially at high temperatures. The temperature-dependent plots in figure 6 reveal that the rate constants with full Wigner correction (k F-Wig ) strongly deviate from linearity at low temperatures, indicating a strong quantum effect on bifunctional proton transfer, especially for H 2+ [PBI] 2 in the forward direction for ε = 1 and for H + [PBI] 2 for ε = 1 and 23, figure 6b and c, respectively.
To complete the assessment of the quantum effect, the correlations of the rate constants obtained based on different methods are included in the electronic supplementary material, figure S3 and shown as examples in figure 7. The correlation plots suggest that for the rate-determining scenarios, the linear relationships between k S-Wig and k Q-vib exist over the studied temperature range with almost the same slopes and intercepts; for example, for H 2+ [PBI] 2 in the forward direction with ε = 23 (CC-     royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 211168 preparation methods, measurement techniques and temperatures. Therefore, the agreement between the theoretical and experimental data is considered to be satisfactory.

Mobility of the phenyl groups
To prove the hypothesis of Fontanella et al. [10] that polymer segmental motion plays an important role in ion conductivity, potential energy curves for the torsional motion of the phenyl groups (ω 1 and ω 2 in figure 8a) were constructed for all of the dimers involved in the rate-determining processes.

Proton conductivity
It should be noted that the reported experimental proton conductivities for pristine PBI are controversial.
Proton conductivity values of σ = 2-8 × 10 −4 S cm −1 were at a relative humidity (RH) of 0-100% [19], whereas considerably lower values (approx. 10 −12 S cm −1 ) [2] were suggested to be more accurate, and a pure PBI membrane was concluded to be inapplicable as a solid electrolyte [16]. In this work, based on the pre-exponential factor determined using conductivity measurements in [ figure 9. It appears that our theoretical results, which are based on the Grotthuss mechanism, agree reasonably well over the temperature range of 330-370 K (10 −4 -10 −5 S cm −1 ), above which the experimental proton conductivities are higher, especially for RH > 0%; proton conduction in acid-doped PBI membranes was suggested to result from the Grotthuss mechanism only at temperatures below the boiling point of water (373 K) [34]. It should be noted that because the conductivities measured in experiments depend upon experimental conditions (pH, RH, doping level, etc.), the agreement between the theoretical and experimental values is considered again to be satisfactory.

Conclusion
In this work, the Grotthuss mechanism for bifunctional proton transfer in PBI membranes was studied using the B3LYP/DZP and B3LYP/TZP methods and transition state theory (TST). This study used [PBI] 2 , H + [PBI] 2 and H 2+ [PBI] 2 as model systems, in which the bifunctional proton transfer paths in low (ε = 1) and high (ε = 23) local dielectric environments were optimized using the NEB method. From analysis of the equilibrium structures, relative H-bond interaction energies and energy barriers showed that the B3LYP/DZP and B3LYP/TZP methods yield approximately the same results. It appeared that the potential energy curves and strength of the H-bonds are sensitive to the local dielectric environment, and the uphill bifunctional proton transfer in [PBI] 2 for ε = 1 can be attributed to extraordinarily strong ion-pair H-bonds in the transition structure, regarded as a 'dipolar energy trap'. A high local dielectric environment (ε = 23) stabilizes [PBI] 2 but destabilizes the ion-pair in the transition structure, leading to a decrease in the energy barrier. Because the bifunctional proton   Figure 9. Plots of the conductivities (σ) of acid-doped PBI membranes as a function of temperature obtained from experiments (300% doping level and RH from 0-30%) [6], compared with those of the double-protonated dimer ( transfer paths for [PBI] 2 involve uphill potential energy curves and the Grotthuss mechanism generally takes place in dynamic H-bond networks, in which the formation and cleavage of covalent bonds are the fundamental steps, a dipolar energy trap can prohibit interconversion between the precursor and proton transferred structures, which rules out the possibility for [PBI] 2 to be involved in the Grotthuss mechanism in the PBI membrane. Different results were observed for H + [PBI] 2 and H 2+ [PBI] 2 , in which the increase in the energy barrier for ε = 23 is due to the stabilization of the protonated H-bonds in the transition structure compared with the precursor; while the protonated H-bonds in the transition structure are stronger than those in the precursor, the solute-local dielectric environment interaction is weaker, especially for H 2+ [PBI] 2 . Because the effect is considerably stronger than that of H + [PBI] 2 , one can anticipate that the rate-determining scenario for bifunctional proton transfer in the PBI membrane is most likely characterized by concerted proton exchanges in H 2+ [PBI] 2 , CC-[1] 2+ O CC- [2] 2+, ‡ O CC- [3] 2+ ; although the segmental (torsional) motion of the PBI chain can promote bifunctional proton transfer by lowering the energy barrier, the coupled motion of proton transfer and torsional rotation has a relatively high rate constant and does not appear to be the rate-determining process. Analysis of the rate constants obtained based on different approximations over the temperature range of 200-500 K confirmed that the quantum effect is not negligible for N-H + … N H-bonds, especially at low temperatures, and the rates for concerted bifunctional proton transfer in H 2+ [PBI] 2 in a high local dielectric environment obey the Arrhenius equation.
An attempt was made to correlate the present results with experimental data, especially for bifunctional proton transfer pathways that obey the Arrhenius equation. Based on the rate constants obtained with the second-order Wigner correction, the equilibrium constant for bifunctional proton transfer in H 2+ [PBI] 2 in a high local dielectric environment is similar to the experimental acid dissociation constant obtained using the Scatchard method, and the pK a obtained based on the steady state approximation is in excellent agreement with the reported experimental value. Based on the energy barriers and experimental preexponential term, the temperature dependence of the proton conductivity was studied and compared with the experimental data for the acid-doped PBI membranes. It appears that the computed conductivities agree well with the experimental values over the temperature range of 330-370 K; the temperature range over which the proton conductivity obeys the Arrhenius equation. Finally, because the rate-determining scenario for the Grotthuss mechanism involves proton exchange in a high local dielectric environment, a low local dielectric environment can be one of the necessary conditions for effective bifunctional proton transfer in acid-doped PBI membranes. These theoretical results provide insights into the Grotthuss mechanism, which can be used as guidelines for understanding the fundamentals of proton transfers in other bifunctional H-bond systems.
Data accessibility. The data are provided in the electronic supplementary material [35]. Authors' contributions. All authors gave final approval for publication and agreed to be held accountable for the work performed therein.
Competing interests. We declare we have no competing interests. Funding. This work was supported by Suranaree University of Technology (SUT) and by Office of the Higher Education Commission under NRU Project of Thailand. The financial support provided by the Thailand Research Fund (TRF) (grant no. MRG6180120) to J.T. is gratefully acknowledged.