Dynamics of structural diffusion in phosphoric acid hydrogen-bond clusters

Abstract

The dynamics and mechanism of structural diffusion in H₃PO₄ hydrogen-bond (H-bond) clusters were studied under excess proton conditions using ab initio calculations at the RIMP2/TZVP level and Born–Oppenheimer molecular dynamics (BOMD) simulations. The outstanding feature of our theoretical approach was the use of combined vibrational and NMR spectroscopic techniques that enable local dynamical motions and their activation energies to be studied in detail. Based on the concept of presolvation, H⁺(H₃PO₄)_n (n = 2–5) clusters were selected as model systems, in which the rate-determining intermediate complexes in the structural diffusion process were studied in detail. The potential energy curves for proton displacement suggested that Zundel-like complexes are preferentially formed in linear H-bond chains in low local-dielectric environments, whereas Eigen-like complexes are stabilized in high local-dielectric environments. The RIMP2/TZVP results revealed that the Eigen–Zundel–Eigen mechanism, which leads to a positive charge displacement, is a short-range process (n = 2–3), and structural diffusion is enhanced by fluctuations in the H-bond chain length and local-dielectric environment. The vibrational and ¹H NMR spectra obtained from BOMD simulations confirmed the proposed mechanism and suggested that for protonated H₃PO₄ clusters, structural diffusion can proceed without the reorientation of H₃PO₄ molecules (frustrated rotation), as in the case of neat liquid H₃PO₄. Because the H⁺ species cannot be directly probed in experiments, this theoretical study suggested for the first time characteristic properties of exchanging protons in H₃PO₄ systems under excess proton conditions, which can be used to assist experimental investigations.

---

Introduction

The dynamics and mechanism of proton transfer have been topics of interest not only in fundamental physical chemistry but also in the development of proton-conducting materials, which can be effectively used in electrochemical devices such as fuel cells. Because the rate of proton transfer from anode to cathode partly determines the efficiency of a fuel cell, theoretical and experimental studies have been conducted in recent decades to develop electrolytes that possess high proton conductivity and excellent thermal, mechanical and chemical stabilities.^1,2 This is because commercially available fuel cells with hydrated proton-exchange membranes such as Nafion® can only be effectively operated at temperatures lower than the boiling point of water.²

Possessing high proton-solvating and self-ionization abilities, as well as the ability to form extensive hydrogen-bond (H-bond) networks in the condensed phase, phosphoric acid (H₃PO₄) exhibits remarkably high intrinsic proton conductivity in the pure state and in aqueous solutions^1,3 and, therefore, has been used effectively as an electrolyte in distributed fuel-cell power generators.⁴ Due to the presence of various charged species and difficulties in directly probing proton exchange, the dynamics and mechanism of proton transfer in liquid and solid H₃PO₄ are complicated and not well understood;^1,2,5 proton transfer in liquid H₃PO₄ has been suggested to occur almost exclusively through a structural diffusion process, also regarded as the protonic-chain conduction mechanism.⁶

Nuclear magnetic resonance (NMR) spectroscopy has been effectively used to obtain fundamental information on the molecular structures and dynamics of atoms in different local chemical environments.⁷ For a prototypical system such as protonated water clusters, ¹H NMR experiments suggested the chemical shift of the shared proton (“naked proton”) in the Zundel complex (H₅O₂⁺) to be 21.3 ppm,⁸ whereas for neutral water clusters, the chemical shifts in polar and non-polar solvents were 4.6–5.5 ppm.⁹ In regards to H₃PO₄ systems, ¹H and ³¹P pulsed gradient Hahn spin-echo (PGSE) NMR spectroscopy has been used to study mobile species in 85 wt% (14.6 M) H₃PO₄ solution.¹⁰ The self-diffusion coefficients (D) obtained based on the analysis of the ¹H NMR peak at 9.8 ppm and the approximation that H⁺ and ¹H diffuse at the same rate suggested that protons are transferred much faster than phosphorus-carrying species. It was concluded in ref. 10 that the activation energy obtained from the slope of the line width-temperature dependence in conventional ¹H NMR experiments and that derived from a simple Arrhenius-type equation in ¹H PGSE NMR experiments are not the same due to the different time and length scales probed by these two techniques, providing energies of 11 and 25 kJ mol⁻¹, respectively.

¹H and ³¹P pulsed field gradient (PFG)-NMR spectroscopy was used to study the dynamics of mobile species in liquid anhydrous H₃PO₄ over the temperature range of 290–420 K.¹¹ The average self-diffusion coefficients (D) of all ¹H atoms were measured and compared with those obtained from high resolution quasi-elastic neutron scattering (QNS)¹² and conductivity measurements.^6,13 It appeared that the PFG-NMR and QNS methods yielded similar proton self-diffusion coefficients, whereas that obtained from conductivity measurements is significantly smaller because the fast diffusion process was not probed in the conductivity measurements.¹¹ Analysis of the experimental data led to the conclusion that in the glassy state, the self-diffusion coefficients obtained from ¹H and ³¹P PFG-NMR measurements cannot be described by the simple Arrhenius equation.¹¹

The mechanism of proton transfer in “neat” liquid H₃PO₄ was studied using Car–Parrinello molecular dynamics (CPMD) simulations on 54 H₃PO₄ molecules.¹⁴ Based on the analysis of proton coupling correlation functions, the proton transfer under neutral conditions was concluded to involve structural diffusion, which is associated with specific H-bond rearrangements in the surrounding environment. It was anticipated that to prevent the proton from returning to the original H₃PO₄ and to ensure that the negatively and positively charged species are fully separated, the linking H-bond chain should be broken. This can be accomplished through the rotation of an H₃PO₄ molecule due to “frustrated” configurations in the H-bond network. It was also concluded that in neat liquid H₃PO₄, strong, polarizable H-bonds produce coupled proton motion and a pronounced dielectric response in the medium, leading to the formation of extended, polarized H-bond chains.

In our previous work,¹⁵ ab initio calculations and Born–Oppenheimer molecular dynamics (BOMD) simulations were used to study the dynamics and mechanism of proton dissociation and structural diffusion in hydrated H₃PO₄ clusters. This theoretical study focused on the effects of the introduction of water into liquid H₃PO₄, which leads to an increase in the proton conductivity of the liquid.¹³ Because the effects of electron correlation play important roles in many-electron systems and BOMD simulations based on quantum chemical methods that include the effects of electron correlation are very CPU time consuming, the applicability and performance of second-order Møller–Plesset perturbation theory with the resolution of the identity approximation (RIMP2) were assessed in this work. Based on the concept of presolvation and the observation that the intermediate complexes are preferentially formed under excess proton conditions, H⁺(H₃PO₄)(H₂O)_n (n = 1–4) clusters were selected as model systems. The static results revealed that the rate-determining intermediate complex for proton dissociation is preferentially formed in low local-dielectric environments (e.g., ε = 1), whereas proton transfer from the first to the second hydration shell is driven by fluctuations in the number of water molecules in high local-dielectric environments (e.g., ε = 78). The BOMD results over the temperature range of 298–430 K validated the proposed mechanism and rate-determining process and iterated that insights into the dynamics and mechanism of structural diffusion in microhydrated H-bond systems can be obtained based on a carefully selected presolvation model.

In this work, the dynamics and mechanism of structural diffusion in H₃PO₄ H-bonds were studied under excess proton conditions using the RIMP2 method and BOMD simulations. Because the H⁺ species cannot be directly probed in experiments¹⁰ and previous theoretical studies focused on proton transfer in liquid H₃PO₄ under neutral conditions,^14,16 to provide fundamental information that can be used as a guideline in experimental investigations, the present theoretical study focused on the dynamics of exchanging protons (H⁺) in the rate-determining intermediate complexes and the reorientation of the H₃PO₄ H-bond chain (frustrated rotation), which has been suggested to be the determining factor in the rate of proton transfer in neat liquid H₃PO₄.¹⁴ Based on the concept of presolvation, H⁺(H₃PO₄)_n (n = 2–5) clusters were selected as model systems, for which characteristic vibrational and ¹H NMR spectra of the H-bonds in the intermediate complexes were computed. The potential energy curves for proton displacement in ε = 1 and 61 (ref. 17) were constructed and used to suggest the structural diffusion mechanism. The dynamic results obtained from BOMD simulations were used to confirm the proposed mechanism and are discussed in comparison with the available theoretical and experimental data.

Computational methods

Quantum chemical calculations

To consider the effects of electron correlation, the RIMP2 method was employed with the TZVP basis set, abbreviated RIMP2/TZVP. This ab initio method has been tested and used successfully in our previous studies on protonated H-bond systems, e.g., protonated water clusters,¹⁸ H₃PO₄-doped imidazole (Im) systems,¹⁹ and microhydrated H₃PO₄ clusters.¹⁵ Because several of our previous studies involved density functional theory (DFT),^19–22 to investigate the performance of the DFT method, the RIMP2/TZVP results were compared with those obtained from B3LYP/TZVP calculations; although the DFT method has been applied effectively in various H-bond systems, there have been cases in which the DFT method fails to account for the equilibrium structures and interaction energies, especially for systems dominated by dispersion interactions.²³ Three basic steps that were applied successfully in our previous studies^24–26 were used in this work; (1) searching for equilibrium structures of the intermediate complexes in the structural diffusion process using the Test-particle model (T-model) potentials;^15,24,25 (2) refining the equilibrium structures using RIMP2/TZVP and B3LYP/TZVP geometry optimizations and; (3) performing BOMD simulations using the refined intermediate structures as the starting configurations. All of the quantum chemical calculations and BOMD simulations were performed using the TURBOMOLE 7.02 software package.²⁷

Presolvation model

The concept of presolvation,⁶ which was used successfully in our previous studies on liquids and aqueous solutions,^{15,18,19,21,28} was adopted in this work; “the local-energy fluctuations and dynamics in aqueous solutions lead to the weakening or breaking of some H-bonds in the first and second hydration shells, resulting in a reduction in the water coordination number such that the proton-accepting species possesses a hydration structure corresponding to the species into which it will be transformed to complete the structural diffusion process”.¹⁵ Because experiments have shown that H₄PO₄⁺ plays an important role in structural diffusion in liquid H₃PO₄ (ref. 10) and there is only a 10% probability for quasi-coherent proton transfer in a four H-bond chain,¹⁴ to apply the concept of presolvation, H⁺(H₃PO₄)_n (n = 2–5) clusters were selected as model systems. This choice is justified by the conclusion that the dynamics and energetics of the structural diffusion process in protonated H-bonds can be studied reasonably well using the smallest intermediate complexes as representatives.^21,22,28,29

In addition, because the fluctuations of local-dielectric environment have been shown to play an important role in the formation and stabilization of the intermediate complexes,^{15,18,19,21,28} the polarization effects of the surrounding H₃PO₄ molecules were taken into account;^14,30 our previous studies showed that for the intermediate complex in a low local-dielectric environment, the shape of the potential energy curve for proton displacement is characterized by a single-well potential, whereas in a high local-dielectric environment, a double-well potential dominates,^15,31 which is in accordance with the localized-charge solvation (LCS) concept.³² In this work, to account for the dielectric response of the medium, a conductor-like screening model (COSMO)³³ with the dielectric constant of liquid H₃PO₄ (ε = 61)¹⁷ was included in the static and dynamic calculations.

Static calculations

Intermediate complexes. The equilibrium structures of H⁺(H₃PO₄)_n (n = 2–5) clusters were obtained from RIMP2/TZVP and B3LYP/TZVP geometry optimizations, in which the likelihood of proton displacement was anticipated from the asymmetric stretching coordinate (Δd_DA in Fig. 1) and the rate-determining intermediate complexes in the structural diffusion mechanism were characterized by the Δd_DA values being close to zero. The total interaction (ΔE^Tot) and solvation (ΔE^sol) energies were computed using both RIMP2/TZVP and B3LYP/TZVP methods.²⁸ The potential energy curves for proton displacement and for rotation of the H-bond (e.g., ω(1) in Fig. 1) were constructed using RIMP2/TZVP calculations with the freeze- and relax-scan methods, respectively;¹⁸ in the freeze-scan method, all of the equilibrium geometrical parameters, including the R_O–O distance, were kept constant and the R_O–H distance was varied, for example, in the range of 0.7–1.7 Å, whereas in the relax-scan method, all of the equilibrium geometrical parameters were allowed to change when the dihedral angle ω(1) was varied in the range of 0–360°.

Fig. 1 Velocity vectors (v₁^OH,MD, v₂^OH,MD, v₁^OO,MD and v₂^OO,MD) used in the calculations of vibrational spectra of O–H⋯O H-bond; Q₁ = symmetric O–H stretch; Q₂ = asymmetric O–H stretch and; Q₃ = O–O vibration; Δd_DA = asymmetric stretching coordinate; torsional rotation of H-bond (1) (ω(1) = ∠P₁–O₂–O₁₁–P₁₀).

Vibrational and NMR spectra. Characteristics of exchanging protons in the intermediate complexes were discussed based on vibrational frequencies and ¹H NMR chemical shifts obtained using the RIMP2/TZVP and B3LYP/TZVP methods. In this work, the asymmetric O–H stretching frequencies (ν^OH) were computed using the harmonic approximation with a scaling factor of 0.9434.²⁷ The likelihood of the formation of the rate-determining intermediate complex was also anticipated from the threshold asymmetric O–H stretching frequency (ν^OH*) obtained from the plot of ν^OH versus Δd_DA; the relationship between ν^OH and Δd_DA is represented by a reflected normal distribution function, of which the second derivative is zero at ν^OH = ν^OH* (inflection point). ν^OH* is the vibrational frequency at which a Zundel-like complex with shared-proton structure is transformed into an Eigen-like complex with close-contact structure.²⁰ In this work, Zundel-like and Eigen-like complexes were characterized by the shape of the potential energy curve by which the dynamics of exchanging protons are governed, which are single-well and double-well potentials, respectively.

The isotropic shielding constants of exchanging protons were calculated using the gauge-including atomic orbital (GIAO) method;^34,35 the ¹H NMR shielding constants at the electron correlated level were approximated from

where and are the isotropic shielding constants computed at the Hartree–Fock (HF) and RIMP2 levels, respectively.³⁶ The ¹H NMR chemical shifts of the exchanging protons were derived from with respect to tetramethylsilane (TMS), which was 31.97 ppm at the MP2/TZVP level.³⁷ Because the isotropic shielding constants of protons in strong, protonated H-bonds depend only on the H-bond distance and not on the neighboring solvent molecules,³⁸ and because the effects of the continuum solvent (reaction field) on the H-bond structure are generally small,^37–39 it is reasonable to discuss only the ¹H NMR results in the gas phase.

Dynamics calculations

BOMD simulations. In BOMD simulations,²¹ classical equations of motions of nuclei were integrated, whereas forces on nuclei were computed from quantum gradients with the molecular orbitals (MO) updated by solving Schrödinger equations based on the Born–Oppenheimer approximation; within the framework of BOMD simulations, nuclei undergo classical Newtonian dynamics on a quantum potential hypersurface. BOMD simulations are therefore more accurate, as well as more CPU time consuming, compared to classical MD simulations, in which forces on nuclei are determined from predefined empirical or quantum pair potentials.

In this work, the dynamics of proton exchange in ε = 1 and 61 were studied using BOMD simulations over a temperature range of 298–430 K. The time step used to solve the equations of motion was set to 1 fs. In the BOMD simulations, 2000 steps were used in the equilibration calculations and 20 000 steps were used in the property calculations. We were aware of possible statistical errors due to the restricted number of time steps used in the analysis of time-dependent properties and applied the so-called “overlapped data collection” method⁴⁰ in all our dynamic calculations, e.g., the velocity autocorrelation function (VACF) and mean-square displacement (MSD). To provide extra samples without extending the simulation time, the BOMD trajectory (e.g., in Fig. S1a†) was divided into equal overlapped intervals (nstep) with equally spaced time origins (nskip). The overlapped data collection method improved the quality of the averaged dynamic properties by generating additional sets of samples for statistical analyses.⁴⁰

Because information on the effects of temperature on the dynamics of proton exchange is necessary to calculate the activation energy,^{15,18,19,21,28} BOMD simulations were performed over the canonical ensemble (NVT), in which heat exchange between the model system and solvent environment was allowed using a thermostat bath. In this work, a Nosé–Hoover chain thermostat was applied to each degree of freedom in the model systems.^41–44 Setting the thermostat relaxation to twenty times the time step was proven in our previous work^15,19,31 to generate reasonable local-energy fluctuations in the model systems. In this work, to ensure that the energetics and dynamic properties were well represented, the energy conservation and the stability of the numerical integration, as well as drifts (δ_A) in the energies and velocities of the exchanging proton, were evaluated in the NVT-BOMD simulations using the method discussed in ref. 45 and the smallest intermediate complex as a representative.

Vibrational spectra. To study the dynamics of proton exchange, vibrational spectra of H-bonds were computed from the Fourier transformation of the VACF. Because the dynamics of protons in H-bonds are coupled with various degrees of freedom, velocity vectors associated with the vibrational modes of interest were defined and used in calculations of the VACF.⁴⁶ In this work, the symmetric (ν^sym-OH,MD) and asymmetric O–H stretching (ν^asym-OH,MD) frequencies, as well as the O–O vibrational (ν^OO,MD) frequencies, were computed using the velocity vectors defined in Fig. 1.¹⁹ The activation energy (ΔE^‡,Arr) of the structural diffusion process was approximated from the Arrhenius equation based on NVT-BOMD simulations over a temperature range of 298–430 K, from which the first-order rate constants (k) were determined from the exponential relaxation behavior of the envelope of the VACF of the O–O vibration; based on the assumption that the relaxation process follows first-order kinetics, k = ln(2)/τ_1/2, where τ_1/2 is the half-life of the exponential decay function. ΔE^‡,Arr values were computed from the slope of the linear relationship between ln(k) and 1000/T.¹⁹

Proton self-diffusion coefficients. To discuss the dynamics of structural diffusion of the H⁺ species, the proton self-diffusion coefficients (D^§(T)) were computed from NVT-BOMD simulations using the Einstein relation, from which D^§(T) were determined from the slope of the MSD of the exchanging proton. Based on our experience,²² the calculations of the self-diffusion coefficient of a single particle (H⁺) confined in a short H-bond distance are not straightforward due to the limitation of the allowed displacement;^40,47 to obtain the correct linear relationship in the MSD plot, the time interval over which the MSD are calculated cannot be too large, because the value of the self-diffusion coefficient approaches zero when the MSD exceed the square of the maximum possible displacement. For the exchanging proton in a protonated H-bond,²² we have shown that the MSD cannot exceed the square of the maximum possible shuttling distance. Our preliminary BOMD simulations revealed that for the protonated H₃PO₄ H-bond with R_O–O = 2.38 Å and the covalent bond distance R_O–H = 1.02 Å, the maximum possible shuttling distance (R_O–H^Max) is 0.34 Å (Fig. S1b†), corresponding to the MSD value of 0.12 Ų (the square of 0.34 Å). The value of the MSD is associated with the longest time interval of approximately 500 fs (Fig. S1c†), which is in accordance with the result in ref. 22.

In addition, the simple Arrhenius equation,

was used to approximate the activation energy (ΔE^‡,D§) from the plot of ln(D^§(T)) against 1000/T, while our NVT-BOMD simulations yielded self-diffusion coefficients of the H⁺ species (D^§(T)), ¹H PFG-NMR experiments yielded average self-diffusion coefficients of all of the ¹H species (D*(T)), which were represented by a fitting of eqn (3).k_B is the Boltzmann constant, D₀ is the pre-exponential factor, H is the effective energy, T is the absolute temperature, and T₀ is the temperature at which the free volume vanishes.

NMR spectra. The dynamics and energetics of proton exchange were also discussed based on the ¹H NMR chemical shift spectra obtained from NVT-BOMD simulations over a temperature range of 298–430 K. Statistical samplings of the intermediate structures were performed every five BOMD steps, and approximately 2500 intermediate structures were used in the GIAO calculations of the instantaneous ¹H NMR shielding constants . To plot the ¹H NMR spectra using the instantaneous chemical shifts,¹⁸ a histogram with a bin size of 0.05 ppm was constructed. Based on the structure of the histogram and the characteristic chemical shifts obtained from the static calculations, the centers of the peaks were assigned and used as initial estimates in a multi-peak fitting process using Lorentzian peak functions. The fitted Lorentzian peak functions were combined to generate the total ¹H NMR spectra. The activation energies (ΔE^‡,NMR) for the exchange between the shared-proton and close-contact structures were determined based on a ¹H NMR line-width analysis;⁴⁸ was approximated as the full width at half maximum (FWHM), and the effective transverse relaxation time was calculated from using . ΔE^‡,NMR was estimated from the slope of the linear relationship between and 1000/T.⁴⁸

Remarks should be made on the analysis of NMR spectra. In ¹H NMR experiments, a motional correlation time (τ_c) with respect to the ¹H NMR time scale can be used to study the dynamics of chemical exchange processes, in which either slow or fast exchange comprises the limiting cases.⁴⁹ It was cited in ref. 49 “In the slow exchange limit, the NMR spectrum is a superposition of the resonances arising from each species, while in the fast exchange limit, a single peak representing an average chemical shift is observed” and “Hence, changes in the line-widths as a function of temperature can be correlated with an exchange rate Ω, which in the fast exchange limit follows the Arrhenius equation”. In H₃PO₄ systems, because the transverse relaxation time (T₂) was reported to be much shorter than the longitudinal relaxation time (T₁),¹¹ and because the dynamics of short-life processes (e.g., oscillatory shuttling and structural diffusion motions) were of interest, we decided to discuss the energetics and dynamics of structural diffusion using the effective transverse relaxation time .^48,49

Results and discussion

Static results

All of the static results on H⁺(H₃PO₄)_n (n = 2–5) clusters obtained from RIMP2/TZVP and B3LYP/TZVP calculations are included in Tables S1 and S2,† respectively. Correlations between the static properties obtained from these two methods are illustrated in Fig. S2.† Examples of the static results on neutral H₃PO₄ clusters obtained from RIMP2/TZVP calculations are listed in Table S3.† Tables S1 and S2† show that the RIMP2/TZVP and B3LYP/TZVP methods give approximately the same equilibrium structures, and excellent correlations are obtained for all of the static properties; the highest correlation coefficients (R²) are found for Δd_DA and R_O–O with R² = 0.98, whereas the lowest R² values are for the interaction energy per H-bond (ΔE^H-bond) and , R² = 0.92 and 0.93, respectively. These R² values lead to the conclusion that RIMP2/TZVP and B3LYP/TZVP calculations yield approximately the same results, and when computational resources are restricted, B3LYP/TZVP calculations can be applied to H₃PO₄ H-bond systems with confidence.

To simplify the discussion, the H-bond structures in Tables S1–S3† are labeled with a three-character code, Gn-[m] or Cn-[m], where G = H-bond structure in the gas phase, C = H-bond structure in continuum solvent, and n is the number of H₃PO₄ molecules. Different H-bond structures with the same number of H₃PO₄ molecules are distinguished by [m]. For example, according to the three-character code, G3-[1] and G3-[2] are different H-bond structures ([1] and [2], respectively) with the same number of H₃PO₄ molecules (n = 3) in the gas phase (G). In contrast, G2-[2] and C2-[2] are the same H-bond structure (n = 2, m = 2) in the gas phase (G) and continuum solvent (C), respectively. To discuss the direction of proton displacement, H₃PO₄ in the H-bond chains are labeled with Roman numerals in parentheses from (I) to (V).

Likelihood of proton transfer. The equilibrium structures in Tables S1 and S3† reveal that although H₃PO₄ can form extensive three-dimensional H-bond clusters, based on the values of Δd_DA, the rate-determining intermediate complexes with shared-proton structures are preferentially formed in protonated linear H-bond chains in ε = 1. The latter is evident from the plots of ν^OH and Δd_DA in Fig. S3,† in which the inflection point is found only in the gas phase at ν^OH* = 2036 cm⁻¹ and = 0.29 Å. For the neutral H₃PO₄ clusters (Table S3†), cyclic H-bond structures dominate with ν^OH and Δd_DA above the threshold values. These findings rule out the likelihood that the neutral cyclic H-bond structures are the intermediate complex in the structural diffusion process. The predominance of the cyclic H-bond structures in the neutral H₃PO₄ clusters is in accordance with the B3LYP/6-311G** results on (H₃PO₄)_n (n = 2–6),¹⁶ for which the values of Δd_DA are all above the threshold value (Δd_DA ≥ 0.4). The H-bond structures that are potentially involved in the structural diffusion mechanism (n = 2–4) are included in Fig. 2. Because the static results for short and long H-bond chains are not substantially different, the smallest intermediate complexes (n = 2) in ε = 1 are discussed first.

Fig. 2 Linear H-bond structures potentially involved in proton exchange in H⁺(H₃PO₄)_n (n = 2–5) clusters obtained from RIMP2/TZVP calculations in low and high local-dielectric environments, ε = 1 and 61, respectively.

Characteristics of intermediate complexes. The smallest intermediate complexes in the gas phase are represented by shared-proton structures with trans and cis conformations; G2-[1] and G2-[2] in Fig. 2 are characterized by Δd_DA = 0 and 0.18 Å and by ω(1) = 161 and 20°, respectively. ΔE^H-bond of G2-[1] is −153.9 kJ mol⁻¹, which is approximately 3 kJ mol⁻¹ more stable than G2-[2], with ν^OH = 928 and 1516 cm⁻¹, respectively. Table S1† suggests that the characteristic ¹H NMR chemical shifts of protons in G2-[1] and G2-[2] are slightly different at = 19.9 and 18.8 ppm, respectively. Because the potential energy curves for proton displacement in Fig. 3a show the characteristics of single-well potentials, G2-[1] and G2-[2] can be considered Zundel-like complexes. The values of 19–20 ppm are regarded as ¹H NMR signatures of exchanging protons in the protonated H₃PO₄ clusters, and are slightly smaller than that in the prototypical system (H₅O₂⁺), determined in ¹H NMR experiments to be 21.3 ppm,⁸ and are also consistent with the values for the hydrated H₄PO₄⁺ clusters, obtained based on the same method to be 18–20 ppm;¹⁵ the strong ¹H NMR downfield shifts (14–22 ppm) are associated with the deshielded protons in short, strong H-bonds, with ΔE^H-bond approximately 63–160 kJ mol⁻¹.⁵⁰

Fig. 3 (a) Potential energy curves of proton displacement in H⁺(H₃PO₄)₂ in low and high local-dielectric environments obtained from RIMP2/TZVP calculations with a freeze-scan method. (b) Potential energy curve for the interconversion between the trans and cis conformations in low local-dielectric environment obtained from RIMP2/TZVP calculations with a relax-scan method.

Because the potential energy curve in Fig. 3b suggests the energy barrier (ΔE^Rel) for the transformation from the trans to cis conformation to be approximately 4 kJ mol⁻¹, which is slightly higher than thermal energy fluctuations at room temperature (2.6 kJ mol⁻¹), and because R_O–O and Δd_DA do not change in the relax-scan path (approximately 2.38 and 0 Å, respectively), one can conclude that in ε = 1, rotation of the protonated H-bond (ω(1)) is feasible but does not seem to affect the energetics of structural diffusion under excess proton conditions.

Effects of H-bond chain extension. The static results in Table S1† and the potential energy curves in Fig. 3 and 4 reveal the outstanding effects of H-bond chain extension. In ε = 1, when the number of H₃PO₄ molecules is increased from n = 2 to 3, two equilibrium structures are obtained from the RIMP2/TZVP geometry optimizations, namely, G3-[1] and G3-[2], as shown in Fig. 2. They are characterized by close-contact structures with ω(1) = 2 and 56° and ω(2) = 170 and 67°, respectively. Because the potential energy curves for proton displacement in H-bonds (1) and (2) in Fig. 4a exhibit the characteristics of double-well potentials with minima close to the oxygen atoms in H₃PO₄ (II), G3-[1] and G3-[2] are considered Eigen-like complexes with ν^OH values, given in Table S1,† ranging from 2200 to 3050 cm⁻¹. Table S1† also shows that the values for protons in Eigen-like complexes (close-contact structures) are 10–16 ppm, whereas those for protons in neutral cyclic H-bonds (Table S3†) are 6–12 ppm. The latter confirm that the ¹H NMR peak at 9.8 ppm used in experiments^10,51 is associated with protons in neutral H-bonds, which is in accordance with the observation that the undissociated H₃PO₄ molecules dominate in orthophosphoric acid;¹⁷ H-bonds with ¹H NMR chemical shifts below 14 ppm are categorized as “moderate” or “week” H-bonds with ΔE^H-bond less than 63 kJ mol⁻¹.⁵⁰

Fig. 4 Potential energy curves for proton displacement in low and high local-dielectric environments obtained from RIMP2/TZVP calculations with a freeze-scan method. (a) H-Bonds (1) and (2) in H⁺(H₃PO₄)₃. (b) H-Bonds (1), (2) and (3) in H⁺(H₃PO₄)₄.

Further extension of the H-bond chains in G3-[1] and G3-[2] leads to G4-[1] and G4-[2] in Fig. 4b, respectively, in which the characteristic properties and the shape of the potential energy curves of H-bond (2) are similar to those of H-bond (1) in G2-[1] and G2-[2], implying that H-bond (2), which is formed from H₃PO₄ (II) and H₃PO₄ (III), is the precursor (Zundel-like complex) for the next transfer. Further extension of the H-bond chain at H₃PO₄ (IV) results in another Eigen-like complex (G5-[1]) with characteristic double-well potentials at H-bonds (2) and (3). The formation of Zundel-like and Eigen-like complexes leads to a positive charge displacement from H-bond (1) to (2), as shown in Fig. 5a, which is in accordance with the Eigen–Zundel–Eigen mechanism in protonated imidazole (Im) H-bond chains,⁵² in which breaking and forming H-bonds in the c crystallographic axis of a unit cell represent the most important elementary steps in the structural diffusion process.³¹

Fig. 5 (a) Eigen–Zundel–Eigen scenario for proton exchange in H⁺(H₃PO₄)_n (n = 2–5) clusters which involves fluctuations of the H-bond chain length. Zundel-like and Eigen-like complexes are shown with the potential energy curves for single-proton displacement in low local-dielectric environment, obtained from RIMP2/TZVP calculations. (b) Short-range Eigen–Zundel–Eigen mechanism for proton exchange in protonated H₃PO₄ H-bond chain which involves the H-bond structures with n = 2 and 3, proposed based on the NVT-BOMD results.

Effects of local-dielectric environment. The equilibrium structures in Fig. 2 and the static results in Table S1† suggest that in COSMO (ε = 61), only cis conformations are formed, with Δd_DA ≥ 0.34 Å and ν^OH ≥ 1986 cm⁻¹. In general, it appears that when the dielectric response of the solvent is included in the calculations, all of the protonated H-bonds are destabilized; for example, ΔE^H-bond of C2-[2] is −50.5 kJ mol⁻¹, compared with −150.6 kJ mol⁻¹ for G2-[2]. Fig. 3 and 4 show that the destabilization of the protonated H-bond results in significant changes in the shape of the potential energy curve for proton displacement in such a way that double-well potentials dominate, confirming that Eigen-like complexes are promoted in ε = 61. One can therefore conclude that fluctuations in the local-dielectric environment also play an important role in the structural diffusion mechanism; Zundel-like complexes are formed in low local-dielectric environments, whereas the actual proton displacement occurs in high local-dielectric environments through the formation of Eigen-like complexes.

Dynamic results

NVT-BOMD simulations. The energy conservation and velocity plots obtained from NVT-BOMD simulations on H⁺(H₃PO₄)₂ at 298 K are shown in Fig. S4.† The average total (E^Tot), potential (E^Pot) and kinetic (E^Kin) energies, as well as the average velocity (v^H+) of the exchanging proton and temperature (T), are listed in Table S4,† together with the corresponding drifts (δ_A) computed with t_max = 2 ps. Because the highest drift is only 0.12 for v^H+, which is much less than one,⁴⁵ and trends in the energy conservation and velocity plots do not show systematic drifts, the dynamic properties obtained under the selected conditions of the NVT-BOMD simulations, including the Nosé–Hoover thermostat, are considered well represented.

Dynamic behaviors of protons. Due to a stronger H-bond interaction and lower thermal agitation, the characteristic dynamic behaviors of protons are outstanding in ε = 1 at 298 K. To study the dynamic behaviors of protons, the NVT-BOMD results on the smallest intermediate complex (Zundel-like complex with n = 2) are primarily discussed using the proton transfer profiles, the time evolutions of the dihedral angle ω(1) and the vibrational spectra shown in Fig. 6. In ε = 1, both small-amplitude (S) and large-amplitude (L) O–O vibrations and associated oscillatory shuttling and structural diffusion motions²⁰ are observed, for example, in Fig. 6a, with the average R_O–O of 2.41 ± 0.06 and 2.42 ± 0.08 Å, respectively. The time evolution of the dihedral angle ω(1) suggests that at 298 K, the trans conformation (ω(1) = 143 ± 25°) is more preferential than the cis conformation, and there is no correlation between the librational motion (rotational oscillation of H-bond (1)) and proton exchange events over the entire simulation, which confirms that frustrated rotation is not included in the mechanism of structural diffusion in this protonated H-bond system.

Fig. 6 (a) and (b) Examples of proton transfer profiles and time evolutions of the dihedral angle ω(1) of H-bond (1) in H⁺(H₃PO₄)₂ obtained from NVT-BOMD simulations in low local-dielectric environment. (c) and (d) Examples of vibrational spectra of H-bond (1) in H⁺(H₃PO₄)₂ obtained from NVT-BOMD simulations in low local-dielectric environment. R_O–O and R_O–H are in Å and ω(1) in degree. L and S = large-amplitude and small-amplitude O–O vibrations, respectively; A and A′ = oscillatory shuttling peaks; B = structural diffusion peak.

The vibrational spectra in Fig. 6c show two structured peaks associated with the oscillatory shuttling motions (labeled A and A′); the gray curves in all of the vibrational spectra were obtained by fittings of Lorentzian peak functions to the characteristic peaks (red curves) to differentiate the characteristic peaks from vibrational interference peaks. Peaks A and A′ are located at ν_A^asym-OH,MD = 819 cm⁻¹ and ν_A′^asym-OH,MD = 1066 cm⁻¹, whereas a broad peak associated with the motion of structural diffusion (peak B) is seen at ν_B^asym-OH,MD = 1599 cm⁻¹. These characteristic peaks are in accordance with the IR frequencies reviewed in ref. 50, in which the O–H stretching frequencies for complementary acid–base solutions were suggested to be 500–2000 cm⁻¹.

Investigation of the time evolution of the dihedral angle ω(1) in Fig. 6a suggests a “librational dynamic equilibrium” in the trans conformation (ω(1) = 143 ± 25°), which results in a splitting of the oscillatory shuttling peak. Because peaks A and A′ are more structured than peak B, the oscillatory shuttling motions dominate at this temperature. The proton transfer profiles in Fig. 6b reveal that as the temperature is increased to 400 K, the large-amplitude O–O vibration (R_O–O = 2.44 ± 0.08 Å) and the structural diffusion motion dominate with a strong librational motion of H-bond (1) (ω(1) = 99 ± 30°). The predominance of the structural diffusion motion is also manifested in the vibrational spectra in Fig. 6d, in which peak B is prominent at ν_B^asym-OH,MD = 1616 cm⁻¹, whereas peaks A and A′ are considerably smaller and slightly redshifted to ν_A^asym-OH,MD = 791 cm⁻¹ and ν_A′^asym-OH,MD = 993 cm⁻¹.

Dynamics in extended H-bond chains. To discuss the effects of H-bond chain extension, the proton transfer profiles and vibrational spectra of the Zundel-like complexes with n = 2 and 4 obtained from NVT-BOMD simulations at 298 K are compared in Fig. S5.† In ε = 1, although the proton transfer profiles in Fig. S5a and b† show similar characteristics for the exchanging protons, due to the extension of the H-bond chain, the likelihood of oscillatory shuttling in H-bond (2) of H⁺(H₃PO₄)₄ is virtually lower than that in H-bond (1) of H⁺(H₃PO₄)₂; in H⁺(H₃PO₄)₄, the small-amplitude (S) and large-amplitude (L) O–O vibrations are seen with an average R_O–O of 2.42 ± 0.05 and 2.46 ± 0.09 Å, respectively. This is also reflected in the vibrational spectra in Fig. S5d,† in which the structural diffusion peak (peak B) is shifted to a higher wavenumber and dominates the oscillatory shuttling peaks (peaks A and A′); ν_A^asym-OH,MD = 808 cm⁻¹, ν_A′^asym-OH,MD = 1111 cm⁻¹ and ν_B^asym-OH,MD = 1834 cm⁻¹. The predominance of oscillatory shuttling in H-bond (1) in Fig. S5c† suggests that the smallest Zundel-like complex is the rate-determining intermediate complex for structural diffusion in protonated H₃PO₄ H-bond chains.

Effects of local-dielectric environment. The effects of local-dielectric environment are discussed based on the NVT-BOMD results on the smallest Zundel-like complex. The proton transfer profiles and the time evolution of the dihedral angle ω(1) in Fig. 7a reveal that at 298 K, the large-amplitude O–O vibration, the cis conformation and the structural diffusion motion dominate with an average R_O–O = 2.50 ± 0.13 Å and ω(1) = 88 ± 45°, whereas the likelihood of oscillatory shuttling decreases due to an increase in the O–H⁺ covalent character and a decrease in ΔE^H-bond in ε = 61. These findings are confirmed by the vibrational spectra in Fig. 7c, in which the structural diffusion peak (peak B) is slightly blueshifted, ν_B^asym-OH,MD = 1616 cm⁻¹, whereas the oscillatory shuttling peaks (peaks A and A′) cannot be resolved due to strong vibrational interferences below 1500 cm⁻¹. In addition, because the dynamics of a proton in the high local-dielectric environment are governed by a double-well potential (e.g., Fig. 3a), two broad peaks are observed in the asymmetric O–H stretching spectra, as seen in Fig. 7c with peaks C and C′ at ν_C^asym-OH,MD = 2121 cm⁻¹ and ν_C′^asym-OH,MD = 2575 cm⁻¹.

Fig. 7 (a) and (b) Examples of proton transfer profiles and time evolutions of the dihedral angle ω(1) of H-bond (1) in H⁺(H₃PO₄)₂ obtained from NVT-BOMD simulations in high local-dielectric environment. (c) and (d) Examples of vibrational spectra of H-bond (1) in H⁺(H₃PO₄)₂ obtained from NVT-BOMD simulations in high local-dielectric environment. R_O–O and R_O–H are in Å and ω(1) in degree. L = large-amplitude O–O vibration; B = structural diffusion peak; C and C′ = characteristic peaks of proton in double-well potential.

As the temperature is increased to 400 K, although the proton transfer profiles and the time evolution of the dihedral angle ω(1) in Fig. 7b show the same dynamic behaviors of the exchanging proton as those at 298 K, the number of proton exchange events is significantly decreased due to stronger thermal agitation. The effects of thermal agitation are also reflected in the vibrational spectra in Fig. 7d, in which the structural diffusion (peak B) and close-contact peaks (peaks C and C′) are blueshifted to ν_B^asym-OH,MD = 1683 cm⁻¹, ν_C^asym-OH,MD = 2171 cm⁻¹ and ν_C′^asym-OH,MD = 2642 cm⁻¹ with strong interference peaks over the entire spectra. The vibrational spectra are complex when the H-bond chain length is increased from n = 2 to 4 in ε = 61. Fig. S6d† shows that for C4-[1], vibrational interferences make it difficult to differentiate the characteristic peaks from the interference peaks; at 298 K, peaks B, C and C′ are observed approximately at 2018, 2459 and 2925 cm⁻¹, respectively.

Similar vibrational analyses were applied in the study of proton transfer in solid maleimide (MI),⁵³ in which a triclinic crystal structure including eight MI molecules and one unbounded proton was selected as the model system in CPMD simulations. Although this model system was considerably larger than our presolvation model, the envelopes of the VACF of the acidic proton showed similar exponential decay functions, with two characteristic vibrational frequency ranges, approximately from 800–1200 cm⁻¹ and from 2900–3100 cm⁻¹. Our analyses, which were based on the velocity vectors in Fig. 1 (the Q₂ mode), suggest that the former correspond to the shuttling motion of H⁺ in the shared-proton structures, whereas the latter reflect the stretching vibration of O–H⁺ in the close-contact structures.

Structural diffusion mechanism. Although the static results suggest that the shared-proton structures can be formed in the H-bond chains with n = 2 and 4, NVT-BOMD simulations show that oscillatory shuttling motions are outstanding in the smallest Zundel-like complex in ε = 1; the likelihood of oscillatory shuttling decreases in the extended H-bond chain (n = 4), especially in ε = 61. Therefore, the rate-determining process is hypothesized to involve dynamics in the smallest Zundel-like complex in low local-dielectric environments. This hypothesis implies that structural diffusion in the protonated H₃PO₄ H-bond chain is determined (governed) by a short-range process (n = 2 and 3) in which the rate is associated with the lifetime of the H-bond in the smallest Zundel-like complex, as shown in Fig. 5b. To prove this hypothesis, the energetics and dynamics of proton exchange in the smallest Zundel-like complex are further studied and discussed in comparison with available experimental data.

Energetics of the structural diffusion process. The energetics of proton exchange in the smallest Zundel-like complex are primarily discussed using the Arrhenius plots in Fig. 8. It appears that linear relationships between ln(k) and 1000/T are observed both in ε = 1 and 61 with activation energies (ΔE^‡,Arr), approximated from the exponential relaxation behavior of the VACF envelope of the O–O vibration, of 15.8 and 12.5 kJ mol⁻¹, respectively, in reasonable agreement with the activation energies of 7.8–14 kJ mol⁻¹ obtained from ion conductivity measurements of 85–100 wt% H₃PO₄ solutions at 293 K.⁵⁴ The ΔE^‡,Arr values confirm that an increase in the local-dielectric environment reduces the likelihood of oscillatory shuttling and lowers the activation energy for the interconversion between the shared-proton and close-contact structures. In other words, the Zundel-like complex, which is preferentially formed in low local-dielectric environments, can be easily transformed into an Eigen-like complex in high local-dielectric environments. This leads to an increase in the population of protons with structural diffusion in high local-dielectric environments, which is in accordance with the static results in Fig. 3a, which show that the close-contact structure (Eigen-like complex) is stabilized in ε = 61. Because the smallest Zundel-like complex in ε = 1 possess higher ΔE^‡,Arr, further discussions are focused on this rate-determining intermediate complex.

Fig. 8 (a) and (b) Arrhenius plots obtained from NVT-BOMD simulations over the temperature range of 298–430 K in low- and high local-dielectric environments, respectively. ΔE^‡,Arr = activation energy obtained from the exponential relaxation behavior of the envelope of the VACF of the O–O vibration.

The activation energies for proton exchange in protonated H₃PO₄ clusters in ε = 1 and 61 are higher than those for proton dissociation in hydrated H₄PO₄⁺ clusters, ΔE^‡,Arr = 13.1 and 9.9 kJ mol⁻¹ in ε = 1 and 78, respectively,¹⁵ and higher than for proton exchange in the Zundel complex (H₅O₂⁺), ΔE^‡,Arr = 10.2 kJ mol⁻¹ in ε = 1.¹⁸ The trend of ΔE^‡,Arr supports the results of conductivity measurements,¹³ in which the introduction of water into liquid H₃PO₄ increases proton conductivity, which in this case is by lowering the activation energies for proton exchange upon hydration.

Proton self-diffusion coefficients. The MSD plots obtained at 298 and 430 K are shown as examples in Fig. S7a and b,† respectively. Two linear relationships, which indicate a change in the dynamic behavior of proton, are observed on the MSD plots in Fig. S7,† ranging from 30–100 fs and from 100–500 fs. The former is associated with the self-diffusion coefficient of a proton with the oscillatory shuttling motion (peaks A and A′ in the vibrational spectra in Fig. 6c), whereas the latter reflects the structural diffusion motion (peak B in the vibrational spectra in Fig. 6c). Because the oscillatory shuttling is the rate-determining process, only the proton self-diffusion coefficients (D^§) derived from the MSD plots ranging from 30 to 100 fs over a temperature range of 298–430 K are summarized in Table S5,† together with the proton self-diffusion coefficients (D*) obtained from the ¹H PFG-NMR measurements over a temperature range of 290–420 K,¹¹ and will be discussed in detail. The plot of ln(D^§) and 1000/T in Fig. S8† yields the activation energy for proton exchange in the smallest Zundel-like complex, ΔE^‡,D§ = 13.5 kJ mol⁻¹, which is only 2.3 kJ mol⁻¹ lower than the ΔE^‡,Arr obtained from VACF.

To study the temperature dependence of the proton self-diffusion coefficients obtained from the theoretical and experimental methods, the relations between D^§ and D* and 1000/T are illustrated in Fig. 9. It appears that reasonable agreement between the theoretical and experimental values (D^§ and D*) is observed above the melting temperature, whereas at lower temperatures (298–350 K), the proton self-diffusion coefficients obtained from the ¹H PFG-NMR experiments suggest a slower proton diffusion rate. The discrepancy between the theoretical and experimental values at low temperatures occurs because the ¹H PFG-NMR measurements yield the average proton self-diffusion coefficients (D*) of all ¹H in the system,¹¹ including protons that are not directly involved in structural diffusion, whereas D^§ values are obtained from the MSD plots that are constructed based only on the displacement of the exchanging proton (H⁺).

Fig. 9 Plots of proton self-diffusion coefficients (D^§) in low local-dielectric environment as a function of 1000/T. D* = values obtained from ¹H PFG-NMR experiment in ref. 11; D^§ = values obtained from the MSD plots.

Characteristics of NMR spectra. All of the ¹H NMR results obtained from NVT-BOMD simulations over the temperature range of 298–430 K are included in Table S6.† The ¹H NMR chemical shift spectra of the exchanging proton in the smallest Zundel-like complex are approximated by four Lorentzian peak functions and illustrated as examples in Fig. 10; the four Lorentzian peak functions are associated with two oscillatory shuttling modes, one structural diffusion mode and one close-contact mode. The two characteristic oscillatory shuttling peaks (peaks A and A′) appear above 19 ppm, whereas the structural diffusion and close-contact peaks (peaks B and C, respectively) are located at lower ppm. Because the oscillatory shuttling modes are confirmed to determine the rate of structural diffusion in the protonated H₃PO₄ H-bond chains and the structural diffusion peak does not show motional narrowing with respect to temperature, only peaks A and A′ are analyzed and discussed in detail.

Fig. 10 (a)–(d) ¹H NMR chemical shift spectra of the exchanging proton in H-bond (1) in H⁺(H₃PO₄)₂ obtained from NVT-BOMD simulations. A and A′ = oscillatory shuttling peaks; B = structural diffusion peak; C = close-contact peak.

In Fig. 10a, it appears that at 298 K, peak A is outstanding at 19.8 ppm, and peak A′ appears as a shoulder at 19.2 ppm, whereas peak B is centered at 18.4 ppm. Integrations of peaks A and A′ yield the populations of the two oscillatory shuttling modes in arbitrary units of approximately 16 and 8, respectively (2 : 1). Investigation of the structures of peaks A and A′ over a temperature range of 298–430 K reveals the temperature dependence of oscillatory shuttling, as well as the motional narrowing of these two modes ( in Table S6†); although all of the characteristic peaks are located at approximately the same frequency, the shapes and relative intensities are temperature dependent. For example, at 380 K, peak A′ is more well-defined and separated from peak A, indicating a stronger librational motion in the rate-determining Zundel-like structure at elevated temperatures. Integrations of these two peaks yield a ratio for these two oscillatory shuttling modes of approximately 0.8 : 1; this result is in accordance with the proton transfer profiles and the vibrational analyses, in which the likelihood of oscillatory shuttling is decreased at elevated temperatures. The plots of and 1000/T in Fig. 11 suggest that the activation energies (ΔE^‡,NMR) for proton exchange in the oscillatory shuttling modes (A and A′) are 13.8 and 8.6 kJ mol⁻¹, respectively. The values of ΔE^‡,NMR are consistent with ΔE^‡,D§ and ΔE^‡,Arr, which are obtained based on MSD and VACF, respectively.

Fig. 11 Plots of the natural log of effective transverse relaxation time as a function of 1000/T. ΔE^‡,NMR = activation energy obtained from ¹H NMR line-width analysis. (a) and (b) For the oscillatory shuttling motions (peaks A and A′, respectively). (c) For the structural diffusion motion (peak B).

Remarks should be made on the ¹H NMR results discussed in this work and those obtained from ¹H PGSE NMR measurements.⁵¹ To understand the ion conduction mechanisms in H₃PO₄ fuel cells, ¹H, ²H and ³¹P NMR experiments were conducted to measure diffusion coefficients and activation energies of mobile species in water-free H₃PO₄ solutions. Based on the analysis of the ¹H NMR peaks at 10.6 and 12.4 ppm (protons in H-bonds in 105 and 116 wt% H₃PO₄ solutions, respectively) and the approximation that the diffusion rate of ¹H and H⁺ are the same, the activation energies for diffusion obtained from the analysis of these two peaks are 32.4 and 36.1 kJ mol⁻¹, respectively, whereas the activation energies for the correlation times (τ_c) estimated from the ¹H NMR spin-lattice relaxation (T₁) measurements at high temperatures are 21.2 and 21.5 kJ mol⁻¹, respectively. The discrepancies between the experimental activation energies are due to the different methods used in the NMR analyses¹⁰ and the fact that H₃PO₄ solutions are complicated due to the presence of various mobile species that cannot be easily separated (e.g., H₃PO₄, H₆P₂O₇ and H₈P₃O₁₀). Whereas the discrepancies between the theoretical and experimental values are due to the approximation that the diffusion rates of ¹H and H⁺ are the same and the fact that the H⁺ species cannot be directly probed in ¹H NMR experiments. Because the experimental activation energies are systematically higher, they are anticipated to include the dynamics of slower processes (e.g., mobile ¹H species) other than the dynamics of the H⁺ species in the structural diffusion mechanism.

Conclusions

The dynamics and mechanism of structural diffusion in the H₃PO₄ H-bond clusters were studied under excess proton conditions using quantum chemical methods at the RIMP2/TZVP level and NVT-BOMD simulations over a temperature range of 298–430 K. This theoretical study focused on characteristic properties of exchanging protons and on local dynamics that could determine the rate of structural diffusion, as well as on the likelihood that the reorientation of H₃PO₄ molecules plays an important role in the structural diffusion process.

To accomplish our designated goals, based on the concept of presolvation, H⁺(H₃PO₄)_n (n = 2–5) clusters were selected as model systems, for which structures, energetics and spectroscopic properties of potential intermediate complexes in low and high local-dielectric environments were examined. RIMP2/TZVP geometry optimizations suggested that only the linear H-bond structures are involved in the structural diffusion process and that intermediate complexes with shared-proton structures (Zundel-like complexes) are preferentially formed in low local-dielectric environments, whereas close-contact structures (Eigen-like complexes) with stronger O–H⁺ covalent character are stabilized in high local-dielectric environments. The static results also suggested that the positive charge displacement along the H₃PO₄ H-bond chain involves a short-range process, and the fluctuations of the H-bond chain length and local-dielectric environment enhance structural diffusion by transforming the Zundel-like complex (n = 2) to an Eigen-like complex (n = 3), which can be regarded as an Eigen–Zundel–Eigen mechanism. In addition, the potential energy curve for the rotation of the H-bond in the smallest Zundel-like complex revealed that the interconversion between the trans and cis conformations is feasible but not necessarily involved in the structural diffusion process.

Analyses of the vibrational and ¹H NMR spectra obtained from BOMD simulations validated the proposed short-range Eigen–Zundel–Eigen mechanism by showing that the oscillatory shuttling, which is the rate-determining process, is dominant in the smallest Zundel-like complex in low local-dielectric environments. The dynamic results also led to the conclusion that for the protonated H₃PO₄ clusters, in which the Grotthuss mechanism dominates, structural diffusion can proceed without the frustrated rotation, as in the case of neat liquid H₃PO₄. Based on our analyses, the proton self-diffusion coefficients and activation energies obtained from ¹H NMR experiments are anticipated to include the dynamics of slower processes (e.g., mobile ¹H species) other than the diffusion of the H⁺ species. Therefore, the present theoretical study suggested for the first time the energetics, dynamics and spectroscopic properties of the H⁺ species in H₃PO₄ H-bond systems under excess proton conditions, which can be used as guidelines for further experimental investigation.

Acknowledgements

The authors would like to acknowledge the financial support provided by the Royal Golden Jubilee (RGJ) Ph.D. Program of Thailand Research Fund (TRF) and Suranaree University of Technology (SUT) (Grant No. PHD/0085/2552) to Prof. Kritsana Sagarik and Parichart Suwannakham and the Advanced Research Scholarship (Grant No. BRG-5580011) to Prof. Kritsana Sagarik. The high-performance computer facilities provided by the following organizations are gratefully acknowledged: Schools of Mathematics and School of Chemistry, SUT; National e-Science project of the National Electronics and Computer Technology Centre (NECTEC), the National Science and Technology Development Agency (NSTDA). This work was supported in part by the Higher Education Research Promotion and National Research University (HERP-NRU) Project of the Office of the Higher Education Commission (OHEC), Thailand.

References

1. M. F. H. Schuster and W. Meyer, Annu. Rev. Mater. Res., 2003, 33, 233–261.
2. K. Sundmacher, Ind. Eng. Chem. Res., 2010, 49, 10159.
3. M. Schuster, K. D. Kreuer, H. Steininger and J. Maier, Solid State Ionics, 2008, 179, 523–528.
4. S. J. Peighambardoust, S. Rowshanzamir and M. Amjadi, Int. J. Hydrogen Energy, 2010, 35, 9349–9384.
5. H. Pu, W. H. Meyer and G. Wegner, J. Polym. Sci., Part B: Polym. Phys., 2002, 40, 663.
6. N. N. Greenwood and A. Thompson, J. Chem. Soc., 1959, 3485.
7. J. A. Pople, W. G. Schneider and H. J. Bernstein, High Resolution Nuclear Magnetic Resonance, McGraw-Hill, New York, 1959.
8. H.-H. Limbach, P. M. Tolstoy, N. Pérez-Hernández, J. Guo, I. G. Shenderovich and G. S. Denisov, Isr. J. Chem., 2009, 49, 199–216.
9. M. N. Nakahara and C. Wakai, Chem. Lett., 1992, 21, 809–812.
10. S. H. Chung, S. Bajue and S. G. Greenbaum, J. Chem. Phys., 2000, 112, 8515.
11. T. Dippel, K. D. Kreuer, J. C. Lassègues and D. Rodriguez, Solid State Ionics, 1993, 61, 41–46.
12. J. C. Lassègues and D. Cavagnat, Phys. B, 1992, 180–181, 645.
13. R. A. Munson and M. E. Lazarus, J. Phys. Chem., 1967, 71, 3245–3248.
14. L. Vilčiauskas, M. E. Tuckerman, G. Bester, S. J. Paddison and K. D. Kreuer, Nat. Chem., 2012, 4, 461.
15. P. Suwannakham, S. Chaiwongwattana and K. Sagarik, Int. J. Quantum Chem., 2015, 115, 486–501.
16. L. Vilciauskas, S. J. Paddison and K.-D. Kreuer, J. Phys. Chem. A, 2009, 113, 9193–9201.
17. R. A. Munson, J. Chem. Phys., 1964, 40, 2044.
18. C. Lao-ngam, M. Phonyiem, S. Chaiwongwattana, Y. Kawazoe and K. Sagarik, Chem. Phys., 2013, 420, 50–61.
19. J. Thisuwan and K. Sagarik, RSC Adv., 2014, 4, 61992–62008.
20. C. Lao-ngam, P. Asawakun, S. Wannarat and K. Sagarik, Phys. Chem. Chem. Phys., 2011, 13, 4562–4575.
21. S. Chaiwongwattana, M. Phonyiem, V. Vchirawongkwin, S. Prueksaaroon and K. Sagarik, J. Comput. Chem., 2012, 33, 175–188.
22. K. Sagarik, S. Chaiwongwattana, V. Vchirawongkwin and S. Prueksaaroon, Phys. Chem. Chem. Phys., 2010, 12, 918.
23. K. Kwac, C. Lee, Y. Jung and J. Han, J. Chem. Phys., 2006, 125, 244508.
24. H. J. Boehm and R. Ahlrichs, J. Chem. Phys., 1982, 77, 2028.
25. H. J. Boehm and R. Ahlrichs, Mol. Phys., 1985, 55, 1159.
26. K. P. Sagarik and R. Ahlrichs, J. Chem. Phys., 1987, 86, 5117.
27. Turbomole Tutorial V7.0, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2015, Turbomole GmbH, since 2007, available from http://www.turbomole.com.
28. M. Phonyiem, S. Chaiwongwattana, C. Lao-ngam and K. Sagarik, Phys. Chem. Chem. Phys., 2011, 13, 10923.
29. K. Sagarik, M. Phonyiem, C. Lao-ngam and S. Chaiwongwattana, Phys. Chem. Chem. Phys., 2008, 10, 2098.
30. R. Vuilleumier and D. Borgis, Nat. Chem., 2012, 4, 432–433.
31. W. Bua-ngern, S. Chaiwongwattana, P. Suwannakham and K. Sagarik, RSC Adv., 2016, 6, 99391–99403.
32. B. Koeppe, J. Guo, P. M. Tolstoy, G. S. Denisov and H.-H. Limbach, J. Am. Chem. Soc., 2013, 135, 7553–7566.
33. A. Klamt and G. Schuuermann, J. Chem. Soc., Perkin Trans. 2, 1993, 799.
34. R. Ditchfield, Mol. Phys., 1974, 27, 789.
35. R. Ditchfield, J. Chem. Phys., 1976, 65, 3123.
36. D. B. Chesnut, Chem. Phys. Lett., 1995, 246, 235.
37. S. Taubert, H. Konschin and D. Sundholm, Phys. Chem. Chem. Phys., 2005, 7, 2561–2569.
38. D. B. Chesnut and B. E. Rusiloski, J. Mol. Struct.: THEOCHEM, 1994, 314, 19–30.
39. D. B. Chesnut, J. Phys. Chem. A, 2005, 109, 11962.
40. D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, London, 1995.
41. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695.
42. S. Nosé, J. Chem. Phys., 1984, 81, 511.
43. S. Nosé, Prog. Theor. Phys. Supp., 1991, 103, 1.
44. S. Nosé and M. L. Klein, Mol. Phys., 1983, 50, 1055.
45. R. L. Davidchack, J. Comput. Phys., 2010, 229, 9323–9346.
46. P. Bopp, Chem. Phys., 1986, 106, 205.
47. J. M. Haile, Molecular Dynamics Simulations, John Wiley & Sons Ltd, New York, 1997.
48. Y. J. Lee, B. Bingöl, T. Murakhtina, D. Sebastiani, W. H. Meyer, G. Wegner and H. W. Spiess, J. Phys. Chem. B, 2007, 111, 9711–9721.
49. G. Brunklaus, S. Schauff, D. Markova, M. Klapper, K. Müllen and H. W. Spiess, J. Phys. Chem. B, 2009, 113, 6583–7040.
50. T. Steiner, Angew. Chem., Int. Ed., 2002, 41, 48–76.
51. Y. Aihara, A. Sonai, M. Hattori and K. Hayamizu, J. Phys. Chem. B, 2006, 110, 24999–25006.
52. A. Li, Z. Cao, Y. Li, T. Yan and P. Shen, J. Phys. Chem. B, 2012, 116, 12793–12800.
53. X. Li, L. Yan and B. Yue, RSC Adv., 2015, 5, 80220–80227.
54. R. He, Q. Li, G. Xiao and N. J. Bjerrum, J. Membr. Sci., 2003, 226, 169–184.

---

Footnote

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

---