Proton Transport in Triflic Acid Hydrates Studied via Path Integral Car

Dec 7, 2009 - Xiaolu Cheng , Jonathan D. Herr , and Ryan P. Steele ... Ali A. Hassanali , Jérôme Cuny , Michele Ceriotti , Chris J. Pickard , and Mi...
0 downloads 0 Views 1MB Size
16574

J. Phys. Chem. B 2009, 113, 16574–16589

Proton Transport in Triflic Acid Hydrates Studied via Path Integral Car-Parrinello Molecular Dynamics Robin L. Hayes,† Stephen J. Paddison,*,‡ and Mark E. Tuckerman*,¶ Department of Chemistry, Courant Institute of Mathematical Sciences, New York UniVersity, New York, New York, and Department of Chemical and Biomolecular Engineering, UniVersity of Tennessee, KnoxVille, Tennessee ReceiVed: August 13, 2009; ReVised Manuscript ReceiVed: October 29, 2009

The mono-, di-, and tetrahydrates of trifluoromethanesulfonic acid, which contain characteristic H3O+, H5O2+, and H9O4+ structures, provide model systems for understanding proton transport in materials with high perfluorosulfonic acid density such as perfluorosulfonic acid membranes commonly employed in hydrogen fuel cells. Ab initio molecular dynamics simulations indicate that protons in these solids are predisposed to transfer to the water most strongly bound to sulfonate groups via a Grotthuss-type mechanism, but quickly return to the most solvated defect structure either due to the lack of a nearby species to stabilize the new defect or a preference for the proton to be maximally hydrated. Path integral molecular dynamics of the mono- and dihydrate reveal significant quantum effects that facilitate proton transfer to the “presolvated” water or SO3- in the first solvation shell and increase the Zundel character of all the defects. These trends are quantified in free energy profiles for each bonding environment. Hydrogen bonding criteria for HOH-OH2 and HOH-O3S are extracted from the two-dimensional potential of mean force. The quantum radial distribution function, radius of gyration, and root-mean-square displacement position correlation function show that the protonic charge is distributed over two or more water molecules. Metastable structural defects with one excess proton shared between two sulfonate groups and another Zundel or Eigen type cation defect are found for the mono- and dihydrate but not for the tetrahydrate crystal. Results for the tetrahydrate native crystal exhibit minor differences at 210 and 250 K. IR spectra are calculated for all native and stable defect structures. Graph theory techniques are used to characterize the chain lengths and ring sizes in the hydrogen bond network. Low conductivities when limited water is present may be attributable to trapping of protons between SO3groups and the increased probability that protons transfer to waters bridging two different sulfonate groups. Introduction The proton exchange membrane (PEM) is an essential component in fuel cells deemed suitable for transportation and stationary applications. These energy converters have been touted as a key ingredient in a sustainable and carbon neutral energy portfolio.1 The role of the PEM is to facilitate the transport of protons from anode to cathode, while excluding all other species (i.e., the reactant gases). Nafion, a perfluorosulfonic acid (PFSA) ionomer, has been the PEM of choice for the past couple of decades due to its high proton conductivity when fully hydrated accompanied by substantial chemical and mechanical stability at elevated temperatures (∼90 °C) and under the typically corrosive environment of an operating fuel cell. The Achilles heel that prevents widespread adoption of this technology in vehicles today is the poor proton conductivity when the water content drops below 10-15 λ ) H2O/SO3H.2 In practice, either lower than optimal temperatures or complicated schemes are required to maintain sufficient humidity during operation. The new quest is to find a robust, inexpensive PEM material that exhibits high proton conductivity (>0.1 S cm-1) when either little or no water is present so that higher operating temperatures and less expensive catalysts may be used. * To whom correspondence should be addressed. E-mail: (S.J.P.) [email protected]; (M.E.T.) [email protected]. † Department of Chemistry, New York University. ‡ University of Tennessee. ¶ Department of Chemistry, Courant Institute of Mathematical Sciences, New York University.

Undoubtedly the success in such a pursuit will be facilitated through the securing of a molecular-level understanding of proton transport in existing PEMs. The transfer and long-range transport of protons features importantly in other systems including the enzymatic function of proteins and membrane biophysics.3 Common among these biological systems with hydrated PEMs is confinement of water and the consequent constraint of proton transport occurring in an aqueous medium that is quite different from that of bulk water.4-8 Extensive computational investigations of the mobility of an excessive proton in bulk water has provided a molecular description of the Grotthuss or structural diffusion mechanism where the hydrated proton exists in a continuum of structures including the limiting cases of the Eigen (H9O4+) and Zundel (H5O2+) cations.9-11 Although there is still some controversy concerning a complete molecular description, critical to longrange proton transport in bulk water is dynamical hydrogen bonding among the water molecules of the second and possibly third hydration shells of the excess proton,12 and “presolvation”,9 which posits that the water molecule to which the proton is transferred must first possess a solvation pattern characteristic of hydronium (i.e., accepting one, not two, hydrogen bonds). Clearly the circumstances of hydrated protons in PEMs are more complex than bulk water resulting in different structural and dynamic properties of the absorbed water.6 Undoubtedly, this impacts the mechanism by which protons are transported in these materials, particularly at low levels of hydration.13

10.1021/jp907853p  2009 American Chemical Society Published on Web 12/07/2009

Proton Transport in Triflic Acid Hydrates Substantial experimental and theoretical studies have been undertaken to understand structure-function relationships in PFSA membranes, specifically the connections between the hydrated morphology and transport properties. We refer the reader to several recent reviews2,14-16 and provide only a brief overview here. PFSA ionomers consist of a poly(tetrafluoroethylene) (PTFE) backbone with perfluoroether side chains terminating with sulfonic acid groups. The hydrophobic PTFE backbone when combined with the very hydrophilic and highly acidic perfluorosulfonic acid groups in the macromolecules of PFSA membranes exhibits separation when exposed to water. This separation occurs on a very fine scale (typically a few nanometers) between crystalline regions consisting of aggregated PTFE backbones and domains made up of the hydrated protons and tethered sulfonate groups.17 The shape, size, and connectivity of the water and ion regions remain uncertain despite extensive experimental studies. It is recognized, however, that increasing the density of the sulfonic acid groups through increasing the ion exchange capacity of the membrane does increase the proton conductivity provided that the mechanical stability of the membrane is not compromised (i.e., that it does not dissolve when fully hydrated nor become too brittle when either dry or minimally hydrated).18 Hence, understanding proton transport in PFSA systems with only a little water should help facilitate the development of PEMs suitable for application under low humidity conditions. Theoretical calculations and simulations can offer important insights not readily available from experiments, but it is clear that these must examine and connect chemical and physical phenomena across many time and length scales.19 The large system sizes on the order of tens of nanometers required to capture the experimental morphology and long times necessary to follow long-range proton transport preclude a straightforward application of quantum mechanical methods capable of accurately describing protons. Hence, a variety of simulation techniques that focus on specific spatial and temporal regimes have been utilized. Several different coarse-grained modeling schemes have been employed to investigate the features of PEMs on length scales up to tens of nanometers.20-23 Recently, dissipative particle dynamics (DPD) simulations of PFSA membranes have revealed the differences in hydrated morphologies as a function of side chain chemistry, equivalent weight, molecular weight, and hydration.24 There has been a relative abundance of classical molecular dynamics (MD) simulations of PFSA ionomers.25-30 As these simulations assume a priori the hydrated structure of the proton as a hydronium ion they are unable to capture structure diffusion or hopping of the proton. Their limited system sizes of only a few thousand atoms also make the hydrated morphology of the simulated system questionable. Typically “agreement” has been observed in computed hydronium ion diffusion coefficients and experimental proton self-diffusion coefficients for only minimally hydrated PFSA membranes and despite suggestions to the contrary, this does not imply that there is no hopping or shuttling of the proton at low water contents (i.e., λ e 6 H2O/SO3H). Although several empirical valence bond (EVB) schemes have been employed in MD simulations of PEMs31-33 in an effort to capture the transport of protons via structural diffusion, they have all failed to obtain agreement with experimentally measured proton diffusion coefficients. These simulations suggest that protons can be trapped at the water-sulfonate interface via a contact ion pair (CIP)/solvent separated ion pair (SSIP) defect configuration where the defect shuttles between an Zundel cation with one water in the first

J. Phys. Chem. B, Vol. 113, No. 52, 2009 16575 solvation shell of the sulfonate (CIP) and a more stable Eigen cation with the central H3O+ separated by one H2O from the sulfonate (SSIP). Ab initio methods, which allow bonds to break and form during a simulation, provide increased insight into proton transfer, but are limited to even smaller system sizes than classical MD simulations. Electronic structure calculations of trifluoromethanesulfonic (triflic) acid (CF3SO3H)34 and gas phase oligomeric fragments35 of the short-side chain PFSA membrane with one to seven water molecules reveal that the proton dissociates as a function of the number of water molecules hydrating a sulfonic acid and the separation of the -SO3H groups. Similar to the results from EVB simulations at higher water contents, the H3O+ prefers to be separated from the SO3- by one layer of water. Limited water precluded the formation of a second layer. Without a linear-scaling algorithm, first principles condensed phase calculations are either limited to a small unit cell where finite size effects bias the results or a large cell where the thermodynamic properties cannot be converged. The notable exceptions are the hydrates of triflic acid where the atomic positions are known from experiment.36-40 Ab initio molecular dynamics (AIMD) simulations of the monohydrate of triflic acid (H3O+-CF3SO3-)4 identified a defect where one proton in the unit cell is shared between two sulfonate groups and a second proton is a member of a Zundel ion.41 The next level of theoretical accuracy, where both the electrons and protons are treated quantum mechanically, has not been previously reported. Trifluoromethanesulfonic acid mono-, di-, and tetrahydrate crystals (Figure 1) provide a well-characterized model system to study the effects of SO3- interactions on proton transport as a function of water (λ ) 1, 2, and 4 water/SO3-, respectively) and density of protogenic groups. Figure 2 shows a schematic of the hydrogen bonding network, including labels for each unique atom, for each crystal hydrate. All the systems contain a prototypical proton defect surrounded by sulfonates and/or water; the monohydrate has a hydronium (H3O+), the dihydrate has a Zundel cation (H5O2+), and the tetrahydrate has an Eigen cation (H9O4+). Like Nafion, there are separate hydrophobic (fluorocarbon) and hydrophillic (SO3+, HxO) regions and extended hydrogen bonding networks. In the mono- and dihydrates, segregation takes the form of layers with isolated hydrogen bonded networks interrupted by sulfonate oxygens (Os). In the tetrahydrate, there is a continuous one-dimensional (1D) water wire between oxygen in HxO (Ow) surrounded by interdigitated CF3 groups. The hydrogen bond network must reorient, however, before long-range proton transport can occur. Defects in the native crystal are needed before long-range proton transport can occur. Moreover, constraints due to the crystal structure and periodic boundary conditions completely inhibit vehicular diffusion of H3O+. Hence, this work aims to investigate how the local environment influences defect structure stability and a small number of structural diffusion events. More specifically, the goals of this study are 3-fold: (1) it aims to explore proton transfer events in the different local environments provided by the three hydrates of triflic acid; (2) it seeks to quantify any quantum effects accompanying the proton transfer events; and (3) it proposes to find additional stable protonic defects that may provide insight into how protons move or become trapped under circumstances of low hydration and high PFSA density. Ultimately, we hope to sketch a rough map of how Nafion fails at low hydration in order to point the way to the design of better PEMs. The remainder of this paper is organized as follows. The simulation cells and methodology are first described for each

16576

J. Phys. Chem. B, Vol. 113, No. 52, 2009

Hayes et al.

Figure 1. Average configuration of the (a) monohydrate [H3O+-CF3SO3-], (b) dihydrate [H5O2+-CF3SO3-], and tetrahydrate [H9O4+-CF3SO3-] native crystals. Features are identified by the following scheme: H (silver spheres), hydronium ions (blue spheres), Ow in water (red spheres), Os in SO3- (red rod), S (yellow rods), C (aqua rods), F (green rods), and H-bonds (thin red lines). Two unit cells are shown in each direction. The top row displays the characteristic proton configuration for each crystal, hydronium, Zundel, and Eigen cations, respectively. The middle row (top view) places the proton defect, circumscribed by an ellipse, within the crystal. The bottom row (side view) reveals that the mono- and dihydrates exhibit water layers separated by CF3 groups, while the tetrahydrate has connected water layers with interdigitated CF3 groups.

of the hydrates. This is followed by a lengthy section describing the results including hydrogen bonding criteria, radial distribution functions, free energy profiles, radius of gyration, rootmean-square displacement, defect species, IR spectra, and characterization of the hydrogen bonding as rings or chains by means of graph theory. Finally, we summarize the main results in the conclusion section. Methodology Simulation Cell. Monoclinic (space group P21/c) lattice constants and positions were taken from experiment, and extrapolated to 15 K below the melting temperature if the thermal expansion coefficients were available, resulting in the following parameters: monohydrate36 at 293 K (H3O+CF3SO3-, a ) 5.9606 Å, b ) 9.9722 Å, c ) 9.7043 Å, β ) 98.685°), dihydrate37 at 252 K (H5O2+CF3SO3-, a ) 10.1532 Å, b ) 8.8704 Å, c ) 7.812 Å, β ) 97.646°), and tetrahydrate39 at 210 K (H9O4+CF3SO3-, a ) 7.6200 Å, b ) 8.6804 Å, c ) 12.8297 Å, β ) 91.940°). To estimate the influence of temperature, the same tetrahydrate crystal was used at 250 K. Each periodic unit cell contains four HxOy+-SO3CF3- formula units (i.e., four proton defects per unit cell).

Computational Methods. First principles density functional based Car-Parrinello molecular dynamics (CPMD)42 and PathIntegral CPMD (PI-CPMD),43 as implemented within the CPMD code,44 were used to study proton defect structures and proton delocalization in both the native and defect crystals. The generalized gradient approximation BLYP exchange correlation functional,45 Troullier-Martins norm conserving pseudopotentials,46 a kinetic energy plane wave cutoff of 180 Ryd, the Γ-point approximation for k-point sampling, a fictitious cp mass of 300 au, and a time step of 0.05 fs ensured that energy differences are reliable to ∼1 kcal/mol. Van der Waal forces, which are necessary to stabilize the native crystals in the direction perpendicular to the CF3 layers (see Supporting Information), were included with dispersion-corrected atomcentered potentials (DCACP)47 on H, O, C, and F. In bulk water, BLYP combined with the complete basis limit48 and dispersion corrections49 has been shown to yield reasonable agreement with experimental properties. Although deuterium replaced hydrogen in all the CPMD simulations (but not PI-CPMD simulations), D will generally be referred to as H. Table 1 summarizes all the simulations. The CPMD simulations started from the geometry-optimized positions, annealed

Proton Transport in Triflic Acid Hydrates

J. Phys. Chem. B, Vol. 113, No. 52, 2009 16577

Figure 2. Schematic of the hydrogen bond network in (a) H3O+-CF3SO3-, (b) H5O2+-CF3SO3-, and H9O4+-CF3SO3- native crystals. Atoms are labeled as follows: blue (hydronium), red (O in water, Ow), orange (O in SO3-, Os), and arrows (H). Arrows point from hydrogen donors to hydrogen acceptors. Numbers identify symmetry unique species. Each unit cell contains four symmetry equivalent HxOy+ - CF3SO3- formula units. The connectivity for both infinite crystals and the periodic unit cells (green box) used in this work are shown. In the monohydrate (a) all the H3O+ are isolated by Os. In the dihydrate (b) there is a 1D hydrogen bond network, but interrupted by Os that will prevent long-range proton transport. In the tetrahydrate (c) there is a 2D hydrogen bond network with a 1D water wire, but the hydrogen bonds have the wrong orientation to facilitate long-range proton transport.

the system to the target temperature, equilibrated the system using the NVT ensemble with Nose´-Hoover chains50 on each ionic degree of freedom for 10 ps, and collected thermodynamic information for at least 10 ps in the NVE ensemble. Native crystals started from the experimental positions while the defect crystals either moved one of the H from the end of the + + corresponding NVE native crystal simulation (d1H3O , d1H5O2 , + + + d2H5O2 , d1,NVTH9O4 ), removed an H (d2H9O4 ) and placed it back + in a different location (d3H9O4 ), or started from a published defect 41 H3O+ ). The purpose of the defect crystal simulations structure (d1 was to find metastable defect structures and promote proton hopping events. Given the limited statistics, the defect crystal

results are suggestive, but not conclusive, of characteristic behavior in these systems. PI-CPMD simulations, which treat both the electrons and nuclei quantum mechanically, can only access quantum statistical quantities, not dynamics. To reduce the estimated relative error in the width of the distance distribution functions51 to less than 3%, the Feynman Path integral was discretized using P ) 32 Trotter slices. Nose´ Hoover chains on both electronic and ionic degrees of freedom and the staging transformation52 were used throughout the NVT simulation. Thermodynamic quantities are calculated after the first 2 ps of the simulation to allow the system to equilibrate.

16578

J. Phys. Chem. B, Vol. 113, No. 52, 2009

Hayes et al.

TABLE 1: Summary of Trajectories name + n1H3O + d1H3O H 3O + d2 + n1H5O2 H5O2+ N1 + d1H5O2 + d2H5O2 H9O4+ n1 + n2H9O4 H9O2+ N1 + N2H9O4 + d1,NVTH9O4 H9O4+ d1,NVE + d2H9O4 HO4+ d3

crystal +

-

H3O CF3SO3 H3O+ CF3SO3H3O+ CF3SO3H5O2+CF3SO3H5O2+CF3SO3H5O2+CF3SO3H5O2+CF3SO3H9O4+CF3SO3H9O4+CF3SO3H9O4+CF3SO3H9O4+CF3SO3H9O4+CF3SO3H9O4+CF3SO3H9O4+CF3SO3H9O4+CF3SO3-

PI-CPMD

initial structure

temp (K)

ensemble

time (ps)

no no no no yes no no no no yes yes no no no no

native move H to OsH-Os OsH-Os, published41 native native move H to OsH move H to OsH-Os native native native native move H to OsH move H to OsH remove H add H to OsH-Os

293 293 293 252 252 252 252 210 250 210 250 210 210 210 210

NVE NVE NVE NVE NVT NVE NVE NVE NVE NVT NVT NVT NVE NVT NVT

24.5 20.8 10.0 24.5 13.9 16.3 10.0 12.7 16.6 13.3 13.1 10.0 9.0 10.0 10.0

Temperature effects were explored by repeating the simulation at 210 and 250 K. Results and Discussion Hydrogen-Bonding Criteria. Hydrogen-bonding criteria in these highly acidic environments are expected to differ from pure water. Ideally, hydrogen-bonding criteria would be extracted from a large liquid simulation. Given the prohibitive cost of such a calculation, these criteria are estimated from two

+

limiting cases, d1,NVTH9O4 which has the maximum disorder (i.e., more phase space is explored) of all the reported simulations, + and N1H9O4 , which is well equilibrated and exhibits several different HxO-Os environments. Following Kumar et al.,53 the 2D potential of mean force (PMF) is calculated to systematically extract the hydrogen bonding criteria. One dimension contains the distribution function of the angle between OdH-Oa (θ), while the other contains the distribution function of the H-Oa distance (d). (Throughout the paper, Od ) hydrogen donor, Oa

Figure 3. Hydrogen bond criteria extracted from the 2D potential of mean force. The x-axis is the hydrogen bond distance (d, Å) and the y-axis is the angle (θ, deg) from all OwH-Ox triplets in the simulation. Results are smoothed with a first order Savitzky-Golay filter. Contour lines are + in units of kT. Maximum values of d and θ are extracted at kT ) -1 (white line). Two limiting cases are considered, (a,b) use d1,NVTH9O4 to estimate H9O4+ to estimate hydrogen bond criteria (HB2) with quantum effects. hydrogen bond criteria (HB1) from the least constrained crystal; (c,d) use N1 Panels (a) and (c) consider water-water interactions (x ) w); panels (b) and (d) consider water-sulfonate interactions (x ) s). Given the poor statistics, OwH-Os results are also used for OsH-Os and OwH-F hydrogen bonds.

Proton Transport in Triflic Acid Hydrates

J. Phys. Chem. B, Vol. 113, No. 52, 2009 16579

TABLE 2: Hydrogen Bond Criteria Taken from Figure 3 (a) HB1 and (b) HB2 name dOwH (Å) ∠OwH-Ow (deg) dOsH (Å) ∠OsH-Ox (deg) HB1 HB2

2.2 2.05

148 149

2.1 2.15

142 141 +

) hydrogen acceptor). For the path integral simulation, N1H9O4 , the angles (θc) and bonds (dc) are calculated from the centroid. The PMF is defined as W(d, θ) ) -kT ln g(d, θ), where, like the radial distribution function, g(d, θ) is the ratio of the average number of atoms in the given shell to the number if the atoms were noninteracting. Figure 3 shows the results for (a,c) OwH-Ow and (b,d) OwH-Os. Given the inadequate statistics, the OwH-Os results are also used for OsH-Os and OwH-F. The hydrogen bonding criteria, shown in Table 2 are taken as the maximum θ and d along the kT ) -1 contour. Using the kT ) 0 contour identified unphysical hydrogen bonds, such as two hydrogen bonds from a single OwH! Given a less constrained system, it might be possible to extract hydrogen bond criteria closer to the saddle point as was done by Kumar et al. Location of Od or H*. Defects are identified in the following manner. First each H is assigned to the nearest O, then hydrogen bonds, using the criteria in Table 2, are found. All H3Ow+ or HOs- are treated as charge defects. Usually, there are four defects, however roughly 0.01% of the time in the PI-CPMD simulations, there are five due to atypical configurations involving two H2O split into OH- and H3O+ such as H3Ow+-OwH--H3Ow+. This does not imply that H3Ow+-OwHare nearest neighbors but rather that the quantum defect structure cannot always be reduced to H3Ow+. Recent infrared spectroscopy measurements54 have also identified the same extended H7O3+ defect species, not Zundel or Eigen cations, for small water clusters in organic solvents. Regardless of the physical existence of H7O3+, given the infrequency of these events, they do not alter the reported statistics. For each defect in each Trotter slice, a most active bond, OdH*, is identified by choosing the OH pair with the smallest value of a particular reaction coordinate δ, where δ ) ROaH - ROdH and Rxy ) distance between atoms x and y. Long range proton transport is not observed in the native crystal simulations, but “rattling” Grotthuss-type proton hops do occur. The lack of permanent proton transfer events in the native crystal may be due to the inability of nearby species to stabilize the new proton defect structures. Figure 5 shows the atomic positions averaged over all time steps, Trotter slices, and equivalent atoms for the di- and tetra- hydrate crystals with and without quantum nuclei. Each O and H is colored according to the probability that they are part of an OdH*. The actual values appear in Table 3. A quantum treatment of the proton significantly increases the delocalization of the defect structure. OsH occurs >0.1% of the time in the PI-CPMD dihydrate simulation, indicating that quantum effects slightly increase the likelihood that defects will be trapped at sulfonate groups when water is scarce. Quantum effects in the tetrahydrate are also significant. There is a minor enhancement of defects temporarily migrating along the water wire. More striking is the 6% increase in probability that the defect will be located at Ow,1, the water in the first solvation shell of the central H3O+ that is most strongly bound to two SO3-. The propensity for defects to reside on Ow,1 provides a driving force to keep proton defects within the first solvation shell of the SO3- interface, at least when the sulfonate density is high. Whether this phenomenon remains in the triflic acid pentahydrate crystal40 or other systems with more than 4 H2O/SO3H is an open question.

Figure 4. Difference plot of normalized Ow*-Oa distance vs δ when Oa accepts either one or two hydrogen bonds and Ow* accepts zero [Ow*H - - - Oa T Ow* - - - HOa] in each defect and Trotter slice. The proton transfers to Oa above δ ) 0 (dashed line). Each distribution is normalized prior to subtraction. Color indicates the probability that more Oa species accept either one (red) or two hydrogen bonds (blue). Results + from tetrahydrate simulations using (a) CPMD [n1H9O4 ] with Oa ) Ow,2, H9O4+ and PI-CPMD [N1 ] with (b) Oa ) Ow,1, (c) Oa ) Ow,2, and (d) Oa ) Ow,3. Poor statistics precluded the inclusion of the corresponding CPMD Oa ) Ow,1 and Oa ) Ow,3 results. The expected presolvation (accepts one hydrogen bond) correlates with short O-O distances in the proton transfer region (δ ) 0).

Ow,2 also donates two hydrogen bonds to SO3- groups, but accepts two hydrogen bonds as well. Lack of proper presolvation may be why Ow,2 is least likely to accept H*; Ow,2 typically accepts two hydrogen bonds. Ow,3 has the expected presolvation, but has only a single hydrogen bond to a SO3- group that is much further away than those in Ow,1 (1.94 vs 1.80 Å). Figure 4 shows the likelihood that the Oa in a given tetrahydrate Ow*-Oa pair in each defect and Trotter slice will accept either one or two hydrogen bonds as function of O-O distance and δ. In all cases (Oa ) Ow,1, Ow,2, or Ow,3) there is a strong correlation between small O-O distances and the expected presolvation (Oa accepting one hydrogen bond) in the proton transfer region (δ ) 0). This does not mean that the proton never transfers when Oa accepts two hydrogen bonds, but rather that the transfer is less likely to occur. Figure 6 estimates the Eigen character of the defects by coloring each O and H according to the percent of time OdH* satisfies |δ| g 0.2 Å. The two oxygen, Ow* and Ow,1, in the Zundel cation of the dihydrate (a,b) still have Eigen cation character ∼50% of the time. In contrast, Ow*, the central oxygen

16580

J. Phys. Chem. B, Vol. 113, No. 52, 2009

Hayes et al.

+

+

+

+

Figure 5. Percent of time O or H is in the most active bond, OdH* for (a) n1H5O2 , (b) N1H5O2 , (c) n1H9O4 , and (d) N1H9O4 . For each defect (OwH3+ or OsH-) in each Trotter slice, OdH* is the OH with the smallest δ in all OdH-Oa triplets, where δ ) ROaH - ROdH and Rxy ) distance between atoms x and y. Atomic positions are averaged over the entire trajectory, all beads, and all equivalent atoms. A subset of molecules, including periodic images of the same molecule, are shown to highlight the local environment. Hydrogen bond distances (d) are given in Å, while select angles (θ) are given in degrees. O and H (spheres) are colored according to the percent of time they are OdH*. Gray atoms are active less than 0.1% of the time. Hydrogen bonds (distance given in Å) are shown as red dashed lines. SO and SC bonds are shown in black. Quantum effects (b,d) significantly delocalize the defect (9%). Only the dihydrate (b) not the tetrahydrate (d) crystal exhibits OsH defects. Protons are more likely to move to water molecules most strongly interacting with Os as seen by the number of Os neighbors and the hydrogen bond distance.

TABLE 3: Location and Character of the Active OdH* Bond Averaged over Trajectory and Equivalent Defects +

+

+

+

+

+

atom

n1H5O2

N1H5O2

n1H9O4

N1H9O4

n2H9O4

N2H9O2

Ow* Ow, 1 Ow, 2 Ow, 3 Os, 1 H1 H2 H3 H4 H5

60.4 (56.8) 39.6 (51.2)

56.0 (53.6) 43.8 (48.2)

0.0 (-) 96.4 (52.9) 2.1 (99.9) 0.5 (100) 0.4 (100) 0.6 (100)

0.1 (6.0) 87.3 (46.1) 5.8 (82.3) 2.6 (89.7) 1.6 (88.7) 2.7 (85.7)

99.7 (93.8) 0.3 (1.6) 0.0 (-) 0.0 (-) 0.0 (-) 55.5 (90.0) 13.4 (99.8) 31.1 (96.8) 0.0 (-) 0.0 (-)

90.6 (70.9) 6.2 (20.8) 0.7 (10.8) 2.5 (15.6) 0.0 (-) 49.4 (61.1) 17.6 (76.0) 32.9 (68.0) 0.0 (-) 0.0 (-)

98.8 (92.1) 1.1 (18.1) 0.0 (-) 0.1 (12.7) 0.0 (-) 54.5 (87.0) 15.8 (97.8) 29.8 (95.4) 0.0 (-) 0.0 (-)

89.5 (70.3) 7.3 (23.2) 0.8 (12.4) 2.4 (16.0) 0.0 (-) 51.1 (59.8) 18.1 (74.5) 30.7 (68.2) 0.1 (82.7) 0.1 (85.8)

Values given for the percent of time O or H is active (percent of time the active bond is an Eigen cation [|δ| g 0.2 Å]) shown in Figures 5, 6, and 7. Percentages