Mechanism of Efficient Proton Conduction in ... - ACS Publications

Dec 3, 2015 - Max Planck Institute for Solid State Research, Heisenbergstrasse 1, D-70569 .... Cláudio M. Lousada , Ricardo M. F. Fernandes , Nadezda...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCB

Mechanism of Efficient Proton Conduction in Diphosphoric Acid Elucidated via First-Principles Simulation and NMR Rachel A. Krueger,†,§ Linas Vilčiauskas,‡,∥ Jan-Patrick Melchior,† Gabriel Bester,†,⊥ and Klaus-Dieter Kreuer*,† †

Max Planck Institute for Solid State Research, Heisenbergstrasse 1, D-70569 Stuttgart, Germany Department of Chemistry, New York University, New York, New York 10003, United States § Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, United States ∥ Department of Chemical Engineering, University of Texas, Austin, Texas 78712, United States ⊥ Department of Chemistry, University of Hamburg, D-20146 Hamburg, Germany ‡

ABSTRACT: Diphosphoric acid (H4P2O7) is the first condensation product of phosphoric acid (H3PO4), the compound with the highest intrinsic proton conductivity in the liquid state. It exists at higher temperature (T > 200 °C) and lower relative humidity (RH ≈ 0.01%) and shows significant ionic conductivity under these conditions. In this work, ab initio molecular dynamics simulations of a pure H4P2O7 model system and NMR spectroscopy on nominal H4P2O7 (which contains significant amounts of ortho- and triphosphoric acid in thermodynamic equilibrium) were performed to reveal the nature and underlying mechanisms of the ionic conductivity. The central oxygen of the molecule is found to be excluded from any hydrogen bonding, which has two interesting consequences: (i) compared to H3PO4, the acidity of H4P2O7 is severely increased, and (ii) the condensation reaction only leads to a minor decrease in hydrogen bond network frustration, which is thought to be one of the features enabling high proton conductivity. A topological analysis of diphosphoric acid’s hydrogen bond network shows remarkable similarities to that of phosphonic acid (H3PO3). The hydrogen bonding facilitates protonic polarization fluctuations (Zundel polarization) extending over several molecules (Grotthuss chains), the other important ingredient for efficient structural diffusion of protons. At T = 160 °C, this is estimated to make a conductivity contribution of about 0.1 S/cm, which accounts for half of the total ionic conductivity (σ ≈ 0.2 S/cm). The other half is suggested to result from diffusion of charged phosphate species (vehicle mechanism) that are present in high concentration, resembling conduction in ionic liquids.

1. INTRODUCTION Proton transfer in hydrogen-bonded systems plays a fundamental role in problems spanning physics, chemistry, and biology.1−6 Experimental7 and theoretical8−17 studies have provided an increasingly detailed picture of proton motion in water. Recent ab initio work has begun to establish a similarly detailed picture of proton transport in phosphoric acid (H3PO4).18 The remarkably high intrinsic conductivity of H3PO4 may be explained by the contributions of a structural diffusion mechanism. Quick, local rearrangements in the hydrogen bond network serve to displace the proton over a short distance within its hydrogen bond, creating a pair of oppositely charged defects. A series of protonic displacements create a polarized Grotthuss chain, with positive and negative defects separated by two or more molecules.18 So that the charge transfer is not reversed along the chain, depolarization must occur through solvent reorganization. Depolarization occurs efficiently because the hydrogen bond network is intrinsically frustrated: the ratio of protons to unprotonated oxygens is 3:1. Frustrated hydrogen bond networks, as observed in H3PO4,18 are characterized by a population of © 2015 American Chemical Society

f rustrated protons: protons covalently bonded to hydrogenbond-accepting oxygens that cannot themselves participate in hydrogen bonds because of a scarcity of acceptors. Frustrated protons readily form new hydrogen bonds that stabilize excess charge to produce irreversible charge separation and thus conduction. One of the remarkable properties of phosphorus oxoacid systems is that they comprise a large family of stable compounds that allow systematic study of the effects of hydrogen network topology, proton density, and interactions on longrange proton transport. One “degree of freedom” in this family is provided by the ability of phosphoric acid to condense and form oligomeric or polymeric species with a general formula of Hn+2PnO3n+1. We are primarily interested in the range of compositions around stoichiometric H3PO4, where H3PO4 and H4P2O7 (Figure 1) are most prevalent. The second “degree of freedom” comes from the ability of the phosphorus atom to Received: October 3, 2015 Revised: December 1, 2015 Published: December 3, 2015 15866

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B

2. COMPUTATIONAL AND EXPERIMENTAL METHODOLOGY The AIMD trajectories for liquid diphosphoric acid were generated using the Car−Parrinello method42 as implemented in the CPMD v3.17 code.43 A cubic simulation box with a side length of 15.68 Å containing 27 H4P2O7 molecules was used. Simulations of an H3PO4 system containing the same number of phosphorus tetrahedra have yielded dynamic and structural data in good agreement with experiment,18,20 indicating that the number of molecules included in the simulation was sufficient to minimize periodicity artifacts. The resulting density of ρ = 2.07 g/cm3 corresponds to a molar volume vH4P2O7 = 2vH3PO4 − vH2O, where vH3PO4 is calculated from an experimental density value of ρ = 1.885 g/cm3 for phosphoric acid just above its melting point. The electronic structure was represented using Kohn− Sham44 single-particle orbitals expanded at the Γ-point in a plane wave basis, with the kinetic energy cutoff set to 100 Ry. The Becke−Lee−Yang−Parr (BLYP)45,46 exchange-correlation functional was used, and core electronic states were treated using norm-conserving Troullier−Martins47 pseudopotentials. Hydrogen atoms were assigned the mass of deuterium to improve adiabatic decoupling, suppress nuclear quantum effects that are beyond the scope of this study, and allow the use of a longer time step (4 au ≈ 0.1 fs). Electrons were assigned a fictitious mass μ of 500 au. A 2 ns classical pre-equilibration was carried out using GROMACS v4.6.348 with modified GROMOS 87 force field parameters.49 This configuration served as the starting point for two AIMD canonical (NVT) equilibration runs. The kinetic temperature was controlled using massive50 Nosé−Hoover chain thermostating.51 The two resulting sets of positions and velocities were then used as starting points for two NVE production runs: a primary 50 ps run with temperature averaging 520 K, and an 18 ps run with an average temperature of 400 K. Both temperatures are well above the melting point of diphosphoric acid (∼70 °C). The higher temperature trajectory allowed observation of a greater number of proton transfer events and faster diffusion, while the lower temperature trajectory allowed direct comparison of the proton dynamics to the previously generated AIMD trajectories for other phosphorus oxoacids. Unless otherwise noted, the results presented are based on analysis of the 50 ps 520 K simulation. A batch of a sample with composition close to H4P2O7 was prepared by reacting P2O5 (Sigma-Aldrich, 99.998%) with an adequate amount of bidistilled water in a closed vessel to avoid material loss during the exothermic reaction. The composition was controlled through recording the intensities of the three well-resolved 31P NMR lines for H3PO4 (39.4 mol % corresponding to 23.3% intensity), H4P2O7 (51.5 mol % corresponding to 60.8% intensity), H5P3O10 (9.0 mol % corresponding to 15.9% intensity), and H6P4O13 (∼0.01 mol %) at T = 400 K. This corresponds to the total composition P2O5· 2.17H2O, which is close to the nominal composition of H4P2O7 (P2O5·2H2O). Samples from the same batch were used for ac impedance and NMR spectroscopy using sample holders without any free volume after filling to avoid compositional changes during heating. For impedance spectroscopy, a conductivity cell with two platinum electrodes of 5 mm diameter facing each other at a separation of about 10 mm was used. The T-shape of the cell allowed filling through the single opening without producing bubbles within the sample space, with a thermocouple positioned in the open leg of the T. The geometrical factor

Figure 1. Structures of phosphonic acid (H3PO3; left), orthophosphoric acid (H3PO4; center), and diphosphoric acid (H4P2O7; right).

form oxoacids in its different oxidation states, yielding H3PO4, H3PO3, and H3PO2.19 In this “matrix” of systems, H3PO4 shows the highest degree of hydrogen bond network frustration as well as intrinsic proton conductivity. Lower proton:oxygen ratios in H4P2O7, H3PO3, and H3PO2 produce lower levels of frustration in the oxoacids’ respective hydrogen bond networks, decreasing the effectiveness of the structural conduction mechanism and lowering the overall conductivity.19,20 Ab initio molecular dynamics simulations of the phosphorus oxoacid family can elucidate the structural and dynamical properties of these systems and clarify the molecular proton conduction mechanisms. Diphosphoric acid (H4P2O7) is a particularly promising candidate for further study. Formed by the self-condensation of H3PO4, H4P2O7 exists at very low relative humidity (RH ≈ 0.01%).21 A small percentage (∼2% at T = 42 °C)22 that severely increases with temperature21 is present even in nominal phosphoric acid (H3PO4) at thermodynamic equilibrium (RH ≈ 0.01%). Understanding the structure and dynamics of the hydrogen bond network formed in the melt of diphosphoric acid (H4P2O7) is thought to be helpful in predicting how H3PO4-containing systems, such as polybenzimidazole−phosphoric acid adducts used as membranes in hightemperature proton exchange membrane (PEM) fuel cells,23−25 respond to decreased humidity. Especially for temperatures above 130 °C, increased electrolyte resistances have been attributed to the formation of H4P2O7.16,26 Since nominal diphosphoric acid exists at low relative humidity, which can easily be maintained at high temperature (200−250 °C), knowing its ionic conductivity and understanding the underlying conduction mechanism is also relevant in the context of the quest for proton-conducting materials operating at high temperature and low humidity.27,28 In fact, proton-conducting solid metal diphosphates, often doped with lower valence metals,29−35 and La-based metaphosphate glasses36,37 do exist at lower humidification than acidic orthophosphates such as CsH2PO4 already employed in “solid acid” fuel cells.38 Diphosphate proton transfer also plays a key role in biological reaction mechanisms39,40 in which P2O74− is generated during ATP hydrolysis. To unravel the basic mechanism of proton transport in lowwater, high-H4P2O7 systems, we performed ab initio molecular dynamics (AIMD) simulations of a model pure H4P2O7 system. AIMD uses density functional theory to calculate forces on nuclei on the fly, capturing polarization, bond breaking, and bond formation to yield results about hydrogen-bonded systems that agree well with experiment.6,10,18,41 Analyzing the simulation trajectories, we examine the structural details of H4P2O7 as well as the dynamics and energetics of proton transfer. To gain macrolevel insight into proton movement in H4P2O7 and to assess the usefulness of our model, we measured conductivities and H and P diffusion coefficients for nominal H4P2O7. The results throw light on proton transport processes in both established and emerging phosphorus oxoacid conductors. 15867

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B

similarity of their molecular structures, which include phosphorus tetrahedra composed of two hydroxide groups and an unprotonated oxygen. Degrees of protonation provide a second point of comparison for the oxoacids’ hydrogen bond structure. The degree of protonation for an oxygen can be estimated by the sum of the empirical bond order52 contributions from each proton in the system. Most oxygens in the H4P2O7 system (Figure 3) fall into two general populations also present in the other phosphorus oxoacids: oxygens donating in

of the cell was determined by measuring the impedance of 0.1 and 0.01 M aqueous solutions of KCl (Merck Certipur) of known specific conductivity. The same setup (HP 4192A LF impedance analyzer) was then used for measuring the conductivity of liquid P2O5·2.17H2O in the temperature range T = 310−415 K. 1 H and 31P tracer diffusion coefficients were determined by the pulsed field gradient stimulated echo (PFG-STE) technique using a Bruker gradient NMR system (diff60 probe, Avance III spectrometer, Ultrashield 400WB Plus magnet). Echo attenuations were recorded as a function of the magnetic field gradient up to Gmax = 29 T/m using an effective pulse duration of δ = 1.2 ms (sinusoidal shape), an additional spoiler gradient, and a constant diffusion time of Δ = 20 ms. Since the 31P signals of orthophosphoric acid (P1), diphosphoric acid (P2), and triphosphoric acid (P3) were well resolved in the Fourier transform of the attenuated echo, the diffusion coefficients of P1, P2, and P3 could be determined separately in the temperature range T = 342−460 K.

3. RESULTS AND DISCUSSION We begin by analyzing the structure of the H4P2O7 system, which is described by radial pair distribution functions. The O−O radial distribution function (RDF) determined from the AIMD simulation for H4P2O7 (Figure 2) shows a peak at r ≈ 2.54 Å for oxygen positions of neighboring neutral molecules, indicating very strong intermolecular hydrogen bonding. This short separation is almost equivalent to the peak separation of oxygens in neutral H3PO3 molecules (r ≈ 2.55 Å) and close to the hydrogen bond length for neutral H3PO4 molecules (r ≈ 2.60 Å). The pattern of contraction around charged defects present in H3PO4 is also preserved in H4P2O7, with positively charged species showing stronger hydrogen bonding but undercoordinated first solvation shells and negatively charged species showing weaker contraction but overcoordinated solvation shells with respect to neutral species. However, the effects of single ionization on coordination are weaker in H4P2O7 than in H3PO4 and H3PO3. This may be attributed to the fact that single ionization in H4P2O7 results from proton transfer to or away from an oxygen in one of the two phosphorus tetrahedra. The change in coordination upon ionization is likely to be more pronounced on the side of the molecule where the transfer occurred and weaker on the other side, producing on average a smaller change in the coordination of the molecule. The near equivalence of H3PO3 and H4P2O7 hydrogen bond lengths already suggests that the high degree of similarity in the liquid structure of bulk H4P2O7 and H3PO3 mirrors the

Figure 3. Probabilities of degrees of protonation [P(dp)] for oxygens in H4P2O7 molecules. Probabilities for H3PO4, H3PO3, and H3PO220 are shown for reference. The degree of protonation for a given oxygen is the sum of the empirical bond orders,52 with all protons falling inside the real space cutoff. Note that the narrow peak at vanishing protonation (shaded blue) corresponds to the central oxygen of diphosphoric acid. The population with a degree of protonation >1 (overprotonated) is shaded orange.

a single hydrogen bond and oxygens accepting two hydrogen bonds. It is important to note that overprotonated oxygens that both donate and accept hydrogen bonds fall into the former category, but have degrees of protonation greater than 1. In H4P2O7, though, a distinct population of oxygens are seen to have degrees of protonation near zero (Figure 3). Careful scrutiny of the simulation trajectory confirmed that oxygens at the center of H4P2O7 molecules neither accept nor donate

Figure 2. O−O radial distribution functions [gOO(r)] and integrated coordination numbers [nOO(r)] for charged species, for (from left to right) H3PO3, H3PO4,20 and H4P2O7. Intramolecular contributions are neglected. 15868

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B hydrogen bonds. A similar, smaller population of oxygens that do not participate in hydrogen bonds also exist in H3PO2, where the ratio of unprotonated oxygens to protons capable of forming hydrogen bonds is 1:1. In the energetic ground-state configuration, then, each oxygen likely participates in exactly one hydrogen bond as either a donor or an acceptor. At higher temperature, though, the increase in configurational entropy is associated with the appearance of a population of overprotonated oxygens and a corresponding population of oxygens that do not form hydrogen bonds. Nonparticipation of central oxygens in the H4P2O7 hydrogen bond network yields a proton to participating oxygen ratio identical to that of H3PO3, indicating similar degrees of intrinsic hydrogen bond network frustration. Thermal fluctuations and proton transfer events also contribute to the observed levels of frustration. To gain a more sophisticated understanding of hydrogen bond geometry and participation, a hydrogen bond definition based on potential of mean force (PMF) values for distance−angle pairs was applied (Figure 4). Developed for analysis of hydrogen bonding in water,53 this definition uses the contour crossing the saddle point of the PMF plot to mark out a population of favorable O−H distance and H−O−H angle combinations as hydrogen-bonded configurations. PMF plots reveal further geometric similarities between hydrogen bond configurations in H4P2O7 and H3PO3, where angle dependence is virtually identical. The presence of saddle points in the PMFs around 50° indicates that configurations in the approach to hydrogen bond formation are nonlinear because of structural constraints. The average participation (bonds donated and accepted) for each H4P2O7 molecule was 7.36 bonds (one oxygen nonparticipating, two oxygens just accepting protons, and four oxygens mainly donating protons). This is just under twice the level of participation for H3PO3 molecules (3.88 bonds). The effects of the hydrogen bond donor−acceptor mismatch that marks a frustrated hydrogen bond network may be quantified by calculating the fraction of hydrogen bonds in which the donating oxygen accepts another hydrogen bond (a donor−acceptor oxygen). As expected, the percentages of hydrogen bonds that include donor−acceptor oxygens in H4P2O7 and H3PO3 are nearly identical, 15.7% and 13.5%, respectively, producing the population of overprotonated oxygens shaded in Figure 3. With the addition of an acidic proton in H3PO4, the prevalence of bonds with donor− acceptor oxygens nearly doubles to 30.7%. The importance of hydrogen bond network frustration for solvent reorganization and conduction is also apparent from the kinetics of hydrogen bond breaking and formation. Population correlation functions of the kind used to describe proton transfer kinetics54 may also be used to characterize hydrogen bond kinetics.55 The decay rates of the autocorrelation functions [C(t)] shown in Figure 5 describe how quickly protons not involved in hydrogen bonds become involved in hydrogen bonds (given that a proton is not involved in a hydrogen bond at time t0, what is the probability that it has not formed a hydrogen bond by time t?), and how quickly existing hydrogen bonds break (given that a hydrogen bond exists between two molecules at time t0, what is the probability that the two molecules remain continuously connected by a hydrogen bond at time t?). Protons not involved in hydrogen bonds have been separated into two populations: frustrated protons, which are covalently bonded to hydrogen bond acceptor oxygens, and nonfrustrated protons, which are covalently bonded to oxygens that do not accept hydrogen

Figure 4. Potentials of mean force as functions of O−H separation (rOH) and the O−H−O hydrogen bond angle (α) (−kT ln(g(rOH,α))) in units of kT. α gives the angular deviation from a linear configuration. The basin formed around the global minimum by the contour line shown crossing the saddle point represents hydrogenbonded molecule pairs.

bonds. It is apparent that frustrated protons form hydrogen bonds more quickly, and that hydrogen bond formation is a much faster process than hydrogen bond breaking. Integrating the correlation functions gives an approximate time scale of 32 fs for frustrated protons to form hydrogen bonds and 39 fs for nonfrustrated protons to form hydrogen bonds. In contrast, hydrogen bond breaking occurs on a time scale of 800 fs. Stabilizing a configuration that results in separation of charge always requires the formation of a new hydrogen bond. When a population of frustrated protons exist, they may form new hydrogen bonds without first breaking old ones, allowing depolarization of Grotthuss chains to occur efficiently and 15869

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B

Figure 6. Relative free energy profiles along the proton transfer coordinate [A(δ)] for H4P2O7 at 400 and 520 K. A(δ) values for H3PO4, H3PO3, and H3PO220 are shown for reference.

Figure 5. Population correlation function for hydrogen bond formation and breaking in H4P2O7 at 520 K. C(t) for hydrogen bond formation considers the population of protons not involved in hydrogen bonds at t = 0, separating these protons into a frustrated population, which are covalently bound to oxygens that accept hydrogen bonds, and a nonfrustrated population bound to oxygens that do not participate in hydrogen bonds. For the hydrogen-bondbreaking function, all hydrogen bonds are considered together.

calculated transfer barrier, as is the case for the excess proton in water.10 Nuclear quantum effects have also been shown to contribute to delocalization of excess protonic charge over nearby molecules.3,10 However, irreversible transfer of protonic charge must be achieved through hydrogen bond breaking and formation, processes that are less sensitive to localized quantum effects. Proton−proton coupling was found to play a key role in overcoming the transfer barrier in H3PO4, with an electrostatic “push” on transferring protons provided by one or more protons approaching the phosphate tetrahedron.18 Resolving the gHH(r) RDF in δ (the δ value is taken for a reference hydrogen atom) (Figure 7) to examine the distribution of other

ensuring that separation of charge is not reversed. The population of frustrated protons present in H3PO4 explains the emergence of a structural conduction mechanism not found in water, where hydrogen bond network frustration is essentially absent. While the high degree of structural similarity between H4P2O7 and H3PO3 provides no evidence for major differences in conduction mechanisms, an examination of the energetics of transfer demonstrates the effects of the electronic structure changes produced by replacing a nonacidic hydrogen in H3PO3 with an oxygen. A free energy profile for proton transfer may be computed using the probability distribution P(ROO,δ), where ROO is the donor−acceptor separation and δ = rH···O − rH−O, according to



A(δ) = −kBT ln[ P(R OO , δ) dR OO]

(1)

It is apparent that the free energy profile for proton transfer in H4P2O7 (Figure 6) breaks the trend of the transfer barrier increasing with hydrogen bond length (H3PO2 < H3PO3 < H3PO4).20 Despite the virtual equivalence of hydrogen bond lengths in H3PO3 and H4P2O7, the proton transfer barrier of ∼0.06 eV for H4P2O7 at 400 K is far lower than the proton transfer barrier for H3PO3 at 383 K. The possibility that the central oxygen in H4P2O7 acts as an electron density sink to stabilize configurations with protons located near the centers of their hydrogen bonds (δ ≈ 0) is supported by prior calculations of gas-phase acidities, which showed ΔGacid to be lower for H4P2O7 (296.3 kcal/mol) than for H3PO4 (318.0 kcal/mol),56 and by a lower experimental pKa value for diphosphoric acid (0.91)57compared to orthophosphoric acid (2.15).58 The lower proton transfer barrier for H4P2O7 as compared to other acids at 400 K already suggests that protons are located, on average, closer to the centers of their hydrogen bonds. At 520 K it falls to ∼0.04 eV, which is close to kBT. A low transfer barrier is a prerequisite for efficient structural diffusion in H3PO4.18 Inclusion of nuclear quantum effects might lower the

Figure 7. HH radial distribution functions resolved in δ, the proton transfer coordinate. The time-averaged value of δ (indicated with dashed lines) is the lowest for diphosphoric acid, the acid with the strongest hydrogen bond. The proton−proton separation peak (rHH ≈ 2.2−2.4 Å) in the transition state of proton transfer (δ ≈ 0) is missing for diphosphoric acid; i.e., the proton is not pushed through the hydrogen bond by electrostatic forces (see the text). 15870

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B

To determine the effect of proton−proton coupling on the formation of polarized Grotthuss chains, which are thought to be key for the autoionization events and fast charge mobility, we analyze the coupled proton transfer propagation probability: if a proton is transferred to a molecule, how likely is it that another proton will be transferred away from the molecule, and what is the range of this coupling in terms of the number of hydrogen-bonded molecules in the system? Here, a slightly modified version of the previously developed18 correlation function for proton transfer as a function of separation in the hydrogen bond network, Cpc(n), is used:

protons around transferring protons produces a distinct peak when δ → 0 for H3PO3 and H3PO4 in the 2.2−2.5 Å range. This corresponds to the distance between two protons hydrogen bonded to a single oxygen, or the covalent and hydrogen-bonded protons of a single oxygen. In both configurations, the motion of one proton can promote the transfer of the other. However, this peak appears more weakly in the H4P2O7 gHH(r) functions, particularly for the higher temperature trajectory, suggesting that the lower barrier allows more transfer configurations to develop without the push provided by a nearby proton. N

Cpc(n) = ⟨∑

N



i=1 j=1 j ≠∀ X

N

Dij(t )|Tt = 0 { ∑ ... k1= 1 k1≠∀ X

N



Djk1(t ′)|tt00+− ττres Dk1k 2(t ″)|tt11+− ττres ... Dkn−1kn(t ′... ′)|ttnn−−11+− ττres }⟩t0 ··· tn res res res

kn = 1 k n≠∀ X

(2)

In the original approach, only a local δ cutoff was used to identify the protons undergoing transfer. However, due to the complicated topology of hydrogen-bonded networks in which the same molecule accepts and donates many hydrogen bonds, which in turn can be involved in several sequences of correlated proton transfer events, this approach likely overestimates the probability of formation of truly one-dimensional Grotthuss chains. One way to improve the accuracy of the estimated chain formation probability is to remember that the hydrogen bond is in fact directional. Hence, every proton transfer can also be assigned a particular direction which is valid even in the δ → 0 limit. In this way, one can only count those transfer events which have compatible directions: if an initial transfer seeding the chain is found, only the “outgoing” protons should be counted for that particular molecule. This allows us to identify the transfer events by the instantaneous reversal of the hydrogen bond direction rather than a δ cutoff. The mathematical form for this particular correlation function, which we term the proton coupling correlation function (PCCF) is not changed from the original expression,18 with the only notable difference being the change of notation where indicator functions are used instead of δ functions. In eq 2, N corresponds to the number of phosphate tetrahedra in the system (2 times the number of molecules in the case of H4P2O7). The chain determination algorithm requires the construction of an adjacency transfer matrix (D) for each frame. Protons are always assumed to form a hydrogen bond between their two nearest neighbor tetrahedra. If tetrahedra i and j are hydrogen bonded at time t with the proton covalently bonded to an oxygen in i, and the proton changes nearest neighbor oxygens to covalently bond to an oxygen in j, Dij is set to 1 at the moment before the change, and Dji is set to −1. Dij is set to 0 for every time at which no transfer occurs. On the chain formation time scale, displacements for phosphorus tetrahedra were found to be small (∼0.5 Å), so relative positions of the tetrahedra were assumed to be fixed when the hydrogen bond direction was determined. The complete set of matrices D may then be searched to identify transfer events occurring in the same direction along a chain starting with bonds adjacent to the initial transfer, where n = 0 is the initial transfer event and n = 1 counts the probability of a transfer event in the adjacent hydrogen bond occurring in the correct direction. The D matrices are searched, and the coupled events are summed up to n = 6, the length of the longest nonrepeating chain in the simulation box. Once the chain of events associated

with every transfer has been recorded, the function is normalized so that Cpc(0) = 1. As was discovered in the case of H3PO4, coupled transfer events may not be instantaneously correlated (quasi-coherent); a certain response time, τres, which is far below the typical hydrogen bond lifetime but long enough to allow for some fast relaxation of the surrounding hydrogen bonds was introduced. A typical value of τres = 50 fs was found to be suitable to capture the relevant events in these studies but may be tuned to limit the portion of the trajectory sampled after each event. The notation T indicates the entire MD trajectory length, and ⟨···⟩t0...tn indicate an ensemble average that depends on the identification of the times t0, ..., tn in T. Figure 8 shows the resulting PCCFs for proton transfers using the directional graph approach for the three phosphorus

Figure 8. Correlation functions for cooperative proton transfer as a function of separation in the hydrogen bond network (n) for H4P2O7, H3PO3, and H3PO4. For H4P2O7, molecules and tetrahedra have alternately been treated as vertices in the graph. Estimated probabilities for random transfer in each system in the absence of coupling are shown with dotted lines. Transfer events are recorded at the moment before a proton changes its nearest neighbor oxygen, and a response time of τres = 50 fs is allowed after transfer in an adjacent bond.

oxoacids. The estimated probabilities of finding a correctdirection transfer in the connected hydrogen bonds at a distance of n away from the original transfer within τres are 15871

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B given. The analytical estimate of chain extension probability based on the expected number of transfer events within τres and the number of hydrogen bonds in which they can occur is shown for comparison (Figure 8). Two functions have been computed for H4P2O7one in which each molecule is treated as a vertex in the graph, and one in which each phosphorus tetrahedron is treated as a vertex. These models represent the extreme cases for intramolecular coupling. The molecule treatment would be accurate in the case in which coupling between the tetrahedra on each side of the molecule is as strong as coupling within the same tetrahedron. The tetrahedron treatment assumes that transfers on one side of the molecule are independent of transfers on the other. It is probable that the level of coupling across the molecule lies somewhere between these extremes. In each of the oxoacids, and in both H4P2O7 models, phosphorus tetrahedra are added to polarized chains significantly more often than they would be in the random case. Interestingly, the probabilities of forming a chain of transfer events within the hydrogen bond network are similar in all of the oxoacids, despite variation in conduction mechanisms and experimental conductivities. It is also possible that some level of intramolecular coupling increases overall coupling in H4P2O7. Hence, the proton−proton coupling at the instant of transfer captured in the gHH(r,δ) RDF plots indicates a more efficient seeding of Grotthuss chains, but not the probability of increasing the length of an existing chain, which is probably governed by a number of other competing effects such as degree of hydrogen bond frustration and dielectric response of the surrounding medium. Although gHH(r,δ) for H4P2O7 does not show a sharp peak to indicate coupled proton motion, this does not indicate that the coupling does not occur. It is likely that the proton configurations associated with transfer in H3PO3 (such as two protons hydrogen bonded to a common oxygen) are also present in H4P2O7. However, the low free energy barrier for transfer in H4P2O7 increases the mobility (Zundel polarizability) of protons within hydrogen bonds, so the peak corresponding to transfer configurations broadens. The structural diffusion mechanism observed in H3PO3 and H3PO4 thus persists intact in H4P2O7. Self-diffusion coefficients calculated from the simulation trajectory and measured experimentally using PFG-NMR (Figure 9) point to a structural contribution to conductivity in H4P2O7. Diffusion coefficients may be calculated from AIMD mean-squared displacements using the Einstein−Smoluchowski relation, and although the results should be cautiously viewed as accurate only in the long-time limit, quantitative agreement between AIMD and PFG-NMR diffusion coefficients for H3PO4 has been reported.20 Diffusion constants for protons and phosphorus atoms in H4P2O7 estimated from the present MD simulation (DH ≈ 3.2 × 10−7 cm2/s and DP ≈ 5.5 × 10−8 cm2/s at 383 K) are an order of magnitude lower than diffusion constants for H3PO4 at the same temperature. In H3PO4, the DH/DP ratio of 5.4 indicates that the proton diffusion is efficiently decoupled from the phosphate diffusion, with structural diffusion of protons between molecules dominating conductivity rather than vehicular diffusion of protons traveling together with the molecules (ions).20 In H4P2O7, decoupling between proton and phosphorus motion occurs to a similar extent despite the lower degree of frustration, with DH/DP ≈ 5.9 when diffusion constants are adjusted to the same temperature as the H3PO4 trajectory, assuming similar temperature dependence.

Figure 9. Proton and phosphorus tracer diffusion coefficients for orthophosphoric acid (H3PO4),59 phosphonic acid (H3PO3),60 and diphosphoric acid (P2O5·2.17H2O ≈ H4P2O7) as obtained by 1H and 31 P PFG-NMR. For diphosphoric acid, DP1, DP2, and DP3 correspond to the diffusion coefficients of ortho-, di-, and triphosphoric acid species which coexist in thermodynamic equilibrium (see the text). Note that the ratio DH/DP (decoupling of H diffusion from P diffusion) is the highest for diphosphoric acid.

PFG-NMR diffusion coefficients determined for P2O5· 2.17H2O show nearly identical levels of decoupling: DH/DP = 4.3 for H4P2O7 and 3 for H3PO4 at T = 150 °C (Figure 9. Note that the sample contained significant concentrations of orthoand triphosphoric acid; see the Computational and Experimental Methodology). The decoupling in H4P2O7 is understandable given the limited movement of bulky H4P2O7 molecules, which interact with neighboring molecules via a high number of hydrogen bonds. Despite significantly lowered diffusion coefficients for both H and P, H4P2O7 still displays relatively high conductivity (∼0.2 S/cm at T = 160 °C) under dry conditions (RH = 0.01%). Figure 10 compares the

Figure 10. Total ionic conductivity σ of H3PO4, H3PO3,59 and P2O5· 2.17H2O (≈H4P2O7) as measured by ac impedance spectroscopy. The contributions from structural diffusion are shown as dashed lines.18,60 The dashed line for diphosphoric acid was obtained by assuming a Haven ratio of unity (see the text). 15872

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B

of charged phosphate species) is surprisingly high (σ ≈ 0.2 S/cm at T = 160 °C) when considering the fact that diphosphoric acid exists at a relative humidity of only ∼0.01%. It is therefore suitable for high-temperature, low-humidity operation (e.g., in fuel cells), which is one of the critical issues in fuel cell technology. In this context, it is worth mentioning that other phosphorus acid systems and their acidic salts require thermodynamic stabilization through relatively high water vapor pressure to ensure high proton conductivity. In the case of phosphonic acid-functionalized polymers, reduced hydrogen bond network frustration and local dynamics as a result of condensation are thought to be the major reasons for their limited conductivity under ambient humidification,62 and the very good solid proton-conducting plastic phase of CsH2PO4 needs to be stabilized by a relatively high water partial pressure to operate in a solid acid fuel cell.38 By contrast, phosphoric acid’s condensation product diphosphoric acid exists at low relative humidity and shows fast proton transport. The present work suggests that this unique combination of properties is essentially the consequence of the fact that condensation does not significantly reduce hydrogen bond network frustration.

conductivity of H4P2O7 to the conductivities of the other oxoacids. The contribution from structural proton diffusion (Dstructure) is calculated by subtracting the proton diffusion contribution associated with the diffusion of the phosphate species as a whole (Dvehicle) from the total as-measured H diffusion coefficient (Dtotal): Dstructure ≈ (Dtotal − Dvehicle)

(3) 5

Dvehicle ≈

3I1PDP1 + 2I2PDP2 + 3 I3PDP3 5

3I1P + 2I2P + 3 I3P

IiP

(4)

31

where the are P NMR intensities for phosphoric, diphosphoric, and triphosphoric acid. The IPi are weighted by the ratio of protons to phosphorus atoms in each acid molecule. The Nernst−Einstein relation σstructure =

e 2c H+Dstructure HRT

(5)



is then used to translate this diffusion coefficient into a conductivity by equating the total proton concentration (note that all protons in the system are exchangeable) to the charge carrier concentration (cH+) in the equation and assuming a Haven ratio (H) of unity, corresponding to uncorrelated proton transfer. Structural conductivities calculated thus represent an upper bound on the range of possible structural diffusion contributions. Although the DH/DP ratio shows a high degree of decoupling between proton and phosphorus diffusion, the structural conductivity contribution (≤50%) in H4P2O7 is far lower than in H3PO4 and H3PO3. This observation can only be explained by high conductivity contributions by charged phosphate species (the vehicle mechanism61). Simulation trajectories showed that, at a given time, only 62% of H4P2O7 molecules on average were neutral, as opposed to 75% and 79% in H3PO4 and H3PO3, respectively. Doubly ionized species were virtually absent in H3PO4 and H3PO3, but 2% of H4P2O7 molecules were doubly ionized, with their diffusion making a conductivity contribution 4 times greater than that of diffusion of singly ionized species.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS R.A.K. acknowledges financial support from the Deutscher Akademischer Austausch Dienst (DAAD) Graduate Research fellowship program (Grant A/12/82840). L.V. acknowledges financial support from the DAAD Postdoctoral Research Fellowship program (Grant D/12/41214). We are grateful to R. Usiskin for proofreading the manuscript. We thank the Max Planck Computing and Data Facility for the use of computational resources.



REFERENCES

(1) Kreuer, K. D. Proton Conductivity: Materials and Applications. Chem. Mater. 1996, 8, 610−641. (2) Kreuer, K. D.; Paddison, S. J.; Spohr, E.; Schuster, M. Transport in Proton Conductors for Fuel-Cell Applications: Simulations, Elementary Reactions, and Phenomenology. Chem. Rev. 2004, 104, 4637−4678. (3) Marx, D. Proton Transfer 200 Years after von Grotthuss: Insights from Ab Initio Simulations. ChemPhysChem 2006, 7, 1848−1870. (4) DeCoursey, T. E. Voltage-Gated Proton Channels and Other Proton Transfer Pathways. Physiol. Rev. 2003, 83, 475−579. (5) Wraight, C. A. Chance and Design - Proton Transfer in Water, Channels and Bioenergetic Proteins. Biochim. Biophys. Acta, Bioenerg. 2006, 1757, 886−912. (6) Marx, D.; Chandra, A.; Tuckerman, M. E. Aqueous Basic Solutions: Hydroxide Solvation, Structural Diffusion, and Comparison to the Hydrated Proton. Chem. Rev. 2010, 110, 2174−2216. (7) Woutersen, S.; Bakker, H. J. Ultrafast Vibrational and Structural Dynamics of the Proton in Liquid Water. Phys. Rev. Lett. 2006, 96, 138305. (8) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. AbInitio Simulations Of Water And Water Ions. J. Phys.: Condens. Matter 1994, 6, A93−A100. (9) Agmon, N. The Grotthuss Mechanism. Chem. Phys. Lett. 1995, 244, 456−462.

4. CONCLUSIONS The molecular structure, hydrogen bonding, and proton transport in diphosphoric acid (H4P2O7) are found to be significantly affected by the fact that the central oxygen of diphosphoric acid is not involved in any hydrogen bonding. This oxygen serves as an electron sink that increases the acidity of the OH groups of the molecule, increases the strength of intermolecular hydrogen bonding, and stabilizes proton transfer by decreasing corresponding transition-state energies. Furthermore, the exclusion of one out of the seven oxygens in the molecule from hydrogen bonding leads to a 2:1 ratio of potential proton donors and acceptors that is identical to that of phosphonic acid (H3PO3), with a correspondingly high level of hydrogen bond network frustration. Both strong hydrogen bonding and hydrogen bond network frustration lead to an efficient decoupling of proton diffusion from phosphate diffusion through a structure diffusion mechanism. From the structure diffusion coefficient, the corresponding conductivity contribution is estimated to account for about 50% of the total conductivity. The total ionic conductivity (comprising both structural proton diffusion and diffusion of a high concentration 15873

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B (10) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. The Nature of the Hydrated Excess Proton in Water. Nature 1999, 397, 601−604. (11) Lobaugh, J.; Voth, G. A. The Quantum Dynamics of an Excess Proton in Water. J. Chem. Phys. 1996, 104, 2056−2069. (12) Schmitt, U. W.; Voth, G. A. The Computer Simulation of Proton Transport in Water. J. Chem. Phys. 1999, 111, 9361−9381. (13) Izvekov, S.; Voth, G. A. Ab Initio Molecular-Dynamics Simulation of Aqueous Proton Solvation and Transport Revisited. J. Chem. Phys. 2005, 123, 044505. (14) Chandra, A.; Tuckerman, M. E.; Marx, D. Connecting Solvation Shell Structure to Proton Transport Kinetics in Hydrogen-Bonded Networks via Population Correlation Functions. Phys. Rev. Lett. 2007, 99, 145901. (15) Berkelbach, T. C.; Lee, H.-S.; Tuckerman, M. E. Concerted Hydrogen-Bond Dynamics in the Transport Mechanism of the Hydrated Proton: A First-Principles Molecular Dynamics Study. Phys. Rev. Lett. 2009, 103, 238302. (16) Chen, H.; Voth, G. A.; Agmon, N. Kinetics of Proton Migration in Liquid Water. J. Phys. Chem. B 2010, 114, 333−339. (17) Hassanali, A. A.; Cuny, J.; Verdolino, V.; Parrinello, M. Aqueous Solutions: State of the Art in Ab Initio Molecular Dynamics. Philos. Trans. R. Soc., A 2014, 372, 20120482. (18) Vilčiauskas, L.; Tuckerman, M. E.; Bester, G.; Paddison, S. J.; Kreuer, K.-D. The Mechanism of Proton Conduction in Phosphoric Acid. Nat. Chem. 2012, 4, 461−466. (19) Vilčiauskas, L.; de Araujo, C. C.; Kreuer, K.-D. Proton Conductivity and Diffusion in Molten Phosphinic Acid (H3PO2): The Last Member of the Phosphorus Oxoacid Proton Conductor Family. Solid State Ionics 2012, 212, 6−9. (20) Vilčiauskas, L. Proton Transport Mechanisms of Phosphoric Acid and Related Phosphorus Oxoacid Systems: A First Principles Molecular Dynamics Study. Ph.D. Thesis, Universität Stuttgart, 2012 (http://elib.uni-stuttgart.de/opus/volltexte/2012/7158/pdf/thesis_ v2.0.pdf). (21) Melchior, J. P. A Multiscale Study of Transport in Model Systems for Proton Conducting Polybenzimidazole Phosphoric Acid Fuel Cell Membranes. Ph.D. Thesis, Universität Stuttgart, 2015. (22) Munson, R. A. Self-Dissociative Equilibria. J. Phys. Chem. 1964, 68, 3374−3377. (23) Zhang, J.; Xie, Z.; Zhang, J.; Tang, Y.; Song, C.; Navessin, T.; Shi, Z.; Song, D.; Wang, H.; Wilkinson, D. P.; et al. High Temperature PEM Fuel Cells. J. Power Sources 2006, 160, 872−891. (24) Ma, Y.; Wainright, J.; Litt, M.; Savinell, R. Conductivity of PBI Membranes for High-Temperature Polymer Electrolyte Fuel Cells. J. Electrochem. Soc. 2004, 151, A8−A16. (25) Jayakody, J. R. P.; Chung, S. H.; Durantino, L.; Zhang, H.; Xiao, L.; Benicewicz, B. C.; Greenbaum, S. G. NMR Studies of Mass Transport in High-Acid-Content Fuel Cell Membranes Based on Phosphoric Acid and Polybenzimidazole. J. Electrochem. Soc. 2007, 154, B242. (26) Lobato, J.; Cañizares, P.; Rodrigo, M. A.; Linares, J. J. PBI-Based Polymer Electrolyte Membranes Fuel Cells. Electrochim. Acta 2007, 52, 3910−3920. (27) Mallant, R. K. a. M. PEMFC Systems: The Need for High Temperature Polymers as a Consequence of PEMFC Water and Heat Management. J. Power Sources 2003, 118, 424−429. (28) Mathias, M. F.; Makharia, R.; Gasteiger, H. A.; Conley, J. J.; Fuller, T. J.; Gittleman, C. J.; Kocha, S. S.; Miller, D. P.; Mittelsteadt, C. K.; Xie, T.; et al. Two Fuel Cell Cars in Every Garage? Electrochem. Soc. Interface 2005, 14, 24−35. (29) Chatzichristodoulou, C.; Hallinder, J.; Lapina, A.; Holtappels, P.; Mogensen, M. Phase Composition and Long-Term Conductivity of Acceptor Doped Ce(PO3)4 and CeP2O7 with Variable P/Metal Ratio and of CeP2O7-KH2PO4 Composite. J. Electrochem. Soc. 2013, 160, F798−F805. (30) Hibino, T.; Kabayashi, K.; Fujita, S. Design of ProtonConducting Sn0.95Al0.05P2O7 with a Mesoporous Structure. J. Mater. Chem. A 2013, 1, 13082.

(31) Jin, Y.; Shen, Y.; Hibino, T. Proton Conduction in Metal Pyrophosphates (MP2O7) at Intermediate Temperatures. J. Mater. Chem. 2010, 20, 6214. (32) Paschos, O.; Kunze, J.; Stimming, U.; Maglia, F. A Review on Phosphate Based, Solid State, Protonic Conductors for Intermediate Temperature Fuel Cells. J. Phys.: Condens. Matter 2011, 23, 234110. (33) Singh, B.; Jeon, S.-Y.; Im, H.-N.; Park, J.-Y.; Song, S.-J. Electrical Conductivity of M2+-Doped (M = Mg, Ca, Sr, Ba) Cerium Pyrophosphate-Based Composite Electrolytes for Low-Temperature Proton Conducting Electrolyte Fuel Cells. J. Alloys Compd. 2013, 578, 279−285. (34) Solomon, J. M.; Adelstein, N.; De Jonghe, L. C.; Asta, M. First Principles Study of Pyrophosphate Defects and Dopant-Defect Interactions in Stronium-Doped Lanthanum Orthophosphate. J. Mater. Chem. A 2014, 2, 1047. (35) Amezawa, K.; Kitajima, Y.; Tomii, Y.; Yamamoto, N.; Widerøe, M.; Norby, T. Protonic Conduction in Acceptor-Doped LaP3O9. Solid State Ionics 2005, 176, 2867−2870. (36) Harley, G.; Yu, R.; De Jonghe, L. C. Proton Transport Paths in Lanthanum Phosphate Electrolytes. Solid State Ionics 2007, 178, 769− 773. (37) Harley, G.; Kreuer, K. D.; Maier, J.; Jonghe, L. C. D. Structural Investigation of Ternary La/Alkaline Earth Phosphate (La(1-x)MxP3Oy) (M = Ba, Ca, Sr) Glasses. J. Non-Cryst. Solids 2009, 355, 932−937. (38) Boysen, D. A.; Uda, T.; Chisholm, C. R. I.; Haile, S. M. HighPerformance Solid Acid Fuel Cells Through Humidity Stabilization. Science 2004, 303, 68−70. (39) Wang, L.; Yu, X.; Hu, P.; Broyde, S.; Zhang, Y. A WaterMediated and Substrate-Assisted Catalytic Mechanism for Sulfolobus solfataricus DNA Polymerase IV. J. Am. Chem. Soc. 2007, 129, 4731− 4737. (40) Lin, I.-C.; Seitsonen, A. P.; Tavernelli, I.; Rothlisberger, U. Structure and Dynamics of Liquid Water from ab Initio Molecular Dynamics−Comparison of BLYP, PBE, and revPBE Density Functionals with and without van der Waals Corrections. J. Chem. Theory Comput. 2012, 8, 3902−3910. (41) Wood, B. C.; Marzari, N. Proton Dynamics in Superprotonic CsHSO4. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 134301. (42) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471−2474. (43) CPMD, version 3.13, Max-Planck-Institut für Festkörperforschung, Stuttgart, and IBM Zurich Research Laboratory, 1995−2010. (44) Kohn, W.; Sham, L. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133. (45) Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098. (46) Lee, C.; Yang, W.; Parr, R. Development of the Colle-Salventti Correlation-Energy Formula into a Functional of the Electron-Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (47) Troullier, N.; Martins, J. Efficient Pseudopotentials for PlaneWave Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 1993−2006. (48) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (49) Spieser, S. A. H.; Leeflang, B. R.; Kroon-Batenburg, L. M. J.; Kroon, J. A Force Field for Phosphoric Acid: Comparison of Simulated with Experimental Data in the Solid and Liquid State. J. Phys. Chem. A 2000, 104, 7333−7338. (50) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. Efficient Molecular Dynamics and Hybrid Monte Carlo Algorithms for Path Integrals. J. Chem. Phys. 1993, 99, 2796. (51) Martyna, G.; Klein, M.; Tuckerman, M. Nosé-Hoover Chains: The Canonical Ensemble via Continuous Dynamics. J. Chem. Phys. 1992, 97, 2635−2643. 15874

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875

Article

The Journal of Physical Chemistry B (52) Pauling, L. Atomic Radii and Enteratomic Distances in Metals. J. Am. Chem. Soc. 1947, 69, 542−553. (53) Kumar, R.; Schmidt, J. R.; Skinner, J. L. Hydrogen Bonding Definitions and Dynamics in Liquid Water. J. Chem. Phys. 2007, 126, 204107. (54) Tuckerman, M. E.; Chandra, A.; Marx, D. A Statistical Mechanical Theory of Proton Transport Kinetics in HydrogenBonded Networks Based on Population Correlation Functions with Applications to Acids and Bases. J. Chem. Phys. 2010, 133, 124108. (55) Luzar, A. Resolving the Hydrogen Bond Dynamics Conundrum. J. Chem. Phys. 2000, 113, 10663−10675. (56) Koppel, I. A.; Burk, P.; Koppel, I.; Leito, I.; Sonoda, T.; Mishima, M. Gas-Phase Acidities of Some Neutral Brønsted Superacids: A DFT and ab Initio Study. J. Am. Chem. Soc. 2000, 122, 5114−5124. (57) Mitra, R. P.; Malhotra, H. C.; Jain, D. V. S. Dissociation Equilibria in Pyrophosphates and Kinetics of Degradation: Part 1Dissociation Constants of Pyrophosphoric Acid. Trans. Faraday Soc. 1966, 62, 167. (58) Ellis, A.; Anderson, D. W. The Effect of Pressure on the First Acid Dissociation Constants of “Sulphurous” and Phosphoric Acids. J. Chem. Soc. 1961, 1765−1767. (59) Dippel, T.; Kreuer, K. D.; Lassègues, J. C.; Rodriguez, D. Proton Conductivity in Fused Phosphoric Acid: A 1H/31P PFG-NMR and QNS Study. Solid State Ionics 1993, 61, 41−46. (60) Schuster, M.; Kreuer, K. D.; Steininger, H.; Maier, J. Proton Conductivity and Diffusion Study of Molten Phosphonic Acid H3PO3. Solid State Ionics 2008, 179, 523−528. (61) Kreuer, K.; Weppner, W.; Rabenau, A. Investigation of ProtonConducting Solids. Solid State Ionics 1981, 3−4, 353−358. (62) Steininger, H.; Schuster, M.; Kreuer, K. D.; Maier, J. Intermediate Temperature Proton Conductors Based on Phosphonic Acid Functionalized Oligosiloxanes. Solid State Ionics 2006, 177, 2457−2462.

15875

DOI: 10.1021/acs.jpcb.5b09684 J. Phys. Chem. B 2015, 119, 15866−15875