Proton Solvation and Transport in Realistic Proton Exchange

Jan 26, 2016 - *Phone 773-702-9092; Fax 773-702-9106; e-mail ... Upon equilibration with molecular dynamics, we find that the bundle morphology ...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCC

Proton Solvation and Transport in Realistic Proton Exchange Membrane Morphologies John Savage and Gregory A. Voth* Department of Chemistry, The James Franck Institute, and Computation Institute, The University of Chicago, 5735 S. Ellis Ave., Chicago, Illinois 60637, United States S Supporting Information *

ABSTRACT: Understanding the role of morphology of perfluorosulfonic acid (PFSA) membranes in defining their proton transport behavior is crucial for the development of the next generation of proton exchange membrane fuel cells. In this study, we build large-scale simulation models of three of the most realistic PFSAs morphologies proposed from the results of SAXS experiments and then examine the cation solvation and transport properties of these models. Upon equilibration with molecular dynamics, we find that the bundle morphology immediately flattens into ribbons, in agreement with the locally flat model. The lamellar model, the extreme version of the locally flat model, shows much too fast dynamical properties and too low density. The random model shows acceptable agreement with experiment; however, the bundle model shows the best overall agreement. The addition of sodium as a co-ion to the system makes the water less structured, and while this does not affect water transport it does slows all ionic transport.

■

by Gierke et al.3 Since then many other models have been suggested, including cylindrical water channels,4 polymer bundles, 5,6 and lamellae. 7,8 Furthermore, the issue is complicated by the fact the morphology may change as hydration level (λ), defined as the number of moles of water per mole of sulfonate groups, is changed or with different processing histories.2 Certain trends have emerged, however. First, it has been shown that the persistence length of Nafion is 3−5 nm,9,10 meaning that on the size scales seen in atomic simulations, one should see relatively straight polymer backbones without the hairpins seen in random morphological models. Second, dilute solutions of Nafion (i.e., volume fraction of Nafion >0.6) have been shown to consist of rod-like polymer bundles,11−13 lending further credence to the polymer bundle model, especially for solution-cast samples of Nafion. Finally, recent analysis from Kreuer et al. indicates that the hydrophobic/hydrophilic interface should be locally flat,14 since the negatively charged sulfonate groups naturally repel each other, arguing against the cylindrical water channel model. Most previous molecular dynamics (MD) simulation work on PFSAs has focused on random morphologies (for reviews, see refs 15 and 16). These simulations all show the phase separation of the hydrated PFSA membranes; however, the extremely long relaxation time and high viscosity of the polymer17 preclude the system from relaxing to its true morphology on the time scale of an MD simulation, even when

INTRODUCTION Hydrogen has been proposed as an alternative energy technology, and in that vein proton exchange membrane (PEM) fuel cells are a proven technology, with potential applications from cellphones to automobiles.1 The key component of these fuel cells is the polymer membrane, which keeps the gaseous fuels separate yet allows protons to pass from anode to cathode and complete the electrical circuit. The polymer membrane material used for many years has been Nafion, which has a polytetrafluoroethylene (PTFE) backbone for mechanical strength and a sulfonated side chain, which provides hydrophilic regions in the polymer that allow proton transport when the membrane is hydrated. This phase separation between the hydrophobic and hydrophilic regions is key to the operation of the membrane, yet details about their arrangement within the membrane are poorly understood because of the amorphous nature of the morphology. Most experimental measurements of perfluorosulfonic acid (PFSA) membranes are on micron thick membranes and as such are only able to obtain an average over the very amorphous structure inherent to PFSAs. Scattering experiments can provide general size scales if a morphological model is assumed, and information on proton transport dynamics throughout the whole membrane can be obtained from conductivity experiments.2 In the present simulation study, we will build morphologies that match scattering information to then determine which morphological model matches the dynamical properties best. The morphology of Nafion has been vigorously debated since the first model, the cluster-network model, was proposed © 2016 American Chemical Society

Received: November 14, 2015 Revised: January 25, 2016 Published: January 26, 2016 3176

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C temperature accelerated dynamics is used.18 In order to investigate the effect of morphology in MD simulations, some plausible morphologies must therefore be assumed. Previous work investigating polymer morphologies has, however, shown interesting behavior. Knox et al. showed the rapid interconnection of water domains using large-scale MD systems19 which utilized a nonreactive MD model in which Grotthuss proton hopping was not allowed (i.e., by using simple classical hydronium cations and hence with only vehicular ion transport). On the other hand, reactive MD simulations of Nafion monomers surrounding hydrophilic regions of cylinders, lamellae, and clusters revealed a complex interplay between the attraction of the hydrated excess protons (aka “hydronium cations”; see discussion next paragraph and ref 20) to the hydrophobic/hydrophilic interface and their diffusion in the center of the pore.21−27 The reactive MD simulations allow for Grotthuss hopping transport of the excess protons along with standard vehicular diffusion.28−31 However, to date, there have been no reactive MD simulations that have combined largescale, realistic polymer morphologies with reactive molecular dynamics, and such simulations are the topic of this paper. Given this background, we chose to simulate three different morphological systems of Nafion of EW 1100 g/mol, lamellae (as the extreme version of the locally flat model postulated by Kreuer14), rod-like polymer bundles,6 and a random morphology as a comparison to previous simulation work on PFSAs.25,26 These morphologies were built to length scales outlined in their respective papers, as obtained from SAXS scattering data. Since determining the mechanisms affecting proton transport in PFSA membranes is the goal of our research, using reactive MD to examine this problem is also crucial.21−31 Hydrated excess protons,20 i.e., the H+ ions in acidic solutions associated with water molecules in the limiting and continually interchanging forms H3O+, H9O4+, and H5O2+, are able to transport in water much faster than other similar ions through the Grotthuss (proton hopping) mechanism, the process by which the excess proton charge transports through the hydrogen bond network of water by exchanging covalent and hydrogen bonds.20,32,33 Although ab initio MD can capture this process, it is infeasible to perform simulations of the time and length scale needed to examine proton transport in these large PFSA systems.27 The multistate empirical valence bond (MS-EVB) reactive MD method was developed to capture the bond breaking and forming in the Grotthuss mechanism in an accurate and highly efficient computational manner.28−31 It has previously been used to reveal interesting and important properties of proton transport in various PFSA materials.21−27 In order to probe the underlying dynamics in these PFSA systems further, we also examined the effect of replacing a fraction of the excess protons with sodium cations. Previous work in our group on anion exchange membranes revealed an interesting co-ion effect in simulations and experiments that contained mixtures of chloride and fluoride anions.34 The selfdiffusion constant of fluoride increased by about 140% when the fluoride content was decreased from 100% to 10% (with 90% chloride). The different solvation structure of the two ions resulted in different attraction to the positively charged counterion and different amounts of water “liberated” upon association to the counterion. The extra water shed by the chloride ion was seen to increase the fluoride diffusion constant. Work in other more simple ionic systems with cations such as the cesium ion and an organic ammonium ion showed the same enhancement effect on the ion mobility.35 In

addition to the water liberation hypothesis, the more disordered water structure and faster water dynamics upon addition of chloride suggested to be also responsible for the enhanced ion mobilities. Experimental studies have examined the effect of mixing sodium and protons (as well as other ions such as Li+, Cs+, and Ba2+).36−39 These studies appear to show a drop in conductivity and water transport upon addition of Na+ to the acidic systems. However, this result is complicated by the fact that the addition of Na+ in these studies removed water from the membrane, which would decrease both conductivity and water transport on its own.

■

SIMULATION DETAILS As we will describe later in more detail, the size of the bundle and lamellar simulations are 3 times larger than any previous reactive MS-EVB simulations performed to date. This system size pushes the boundaries of what is presently possible computationally in terms of both memory and CPU time, and we were therefore only able to perform these simulations using the IBM Blue Gene at the Argonne National Laboratory Leadership Computing Facility. Moreover, this also meant that we were limited in the number of simulations we were able to perform; for example, we could not afford to investigate the effect of temperature on the proton transport. As the point of this study is to evaluate the impact of PFSA morphology on the ion transport, we also wanted to use simulation setups that will actually show proton transport. At low water content the morphology of PFSAs changes at around λ = 5, with the hydrophilic domains becoming disconnected and greatly limiting the proton transport. This effect has been described in detail by us and many others16,26,40 and is independent of the morphological model used.41 Therefore, in order to utilize the computational resources available to us most effectively, we concentrated our efforts on higher hydration levels, i.e., for λ = 10 and 15. The random morphology simulation boxes were built using the Monte Carlo chain-growing algorithm developed by Knox and Voth19 for four chains of polymers containing 10 monomers each (the chains are shown in Figure 1). The two straight chain morphologies (lamellar and bundle) were built with 12 straight chains containing 10 monomers each and placed in the z-direction of a box, as outlined below. The chain length of 10 monomers was chosen to compare to the standard MD simulation chain length used in most previous MD simulations, reactive and nonreactive alike.16 In addition, given that 12 chains are needed to build structures that will match scattering experiment sizes, we were limited in the chain length we could use, as these simulations already pushed the boundaries of what is currently possible computationally. Nevertheless, even though these chain lengths are shorter than those in actual PFSA membranes, their end-to-end distances in the equilibrated structures are 14 and 16 nm in the bundle and lamellar morphologies, respectively. These are 3−5 times the persistence length of Nafion9,10 and so should be a reasonable approximation to the full polymer chain, especially when periodic boundary conditions are taken into account as described below. In the lamellar case, six chains were placed on the yz-plane with their side chains pointing the minus x-direction. A total of six more chains were placed adjacent to these chains with their side chains pointing in the positive x-direction. The spacing of the chains was such that they occupied a volume determined 3177

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C

performed for 200 ps in the constant NVT ensemble for equilibration. Each SCI-MS-EVB system was then run for 1 ns using constant NVE dynamics for the production simulations of the reactive trajectories for analysis. The sodium mixture simulations were built by taking the same starting equilibrations used for the constant NVE production runs and replacing 50%, 75%, or 100% of excess protons with sodium ions. The excess protons to be replaced were randomly chosen to ensure a good dispersion of sodium ions throughout the system. Equilibration of these systems was performed for 200 ps in the constant NVT ensemble. Each of these systems was then run for 2 ns using constant NVE dynamics for production simulations of the reactive trajectories for analysis. The longer time for these production runs were to obtain better statistics from the fewer excess protons compared to the 100% H+ systems. The DREIDING force field45 was used for the polymer model parameters with some changes outlined in ref 26. All charges for Nafion were obtained from a Mulliken population analysis46 on the optimized structures at the DFT B3LYP 631** level in Gaussian 0947 and are shown in ref 25. All bond and angle force constants were taken from ref 48, while the equilibrium bond lengths and angles were taken from our DFT optimized structures. For all polymers, the dihedral parameters for the CF backbone were taken from ref 48. The dihedral parameters for the side chains of Nafion were the default DREIDING parameters. All polymer Lennard-Jones (LJ) parameters were taken from ref 48, and the standard Lorentz−Berthelot mixing rules were used for all interactions except for one important interaction between the sulfonate oxygens and the hydronium hydrogens, as outlined in ref 26. For the hydrated proton and water parameters, the MSEVB3 model,31 which was parametrized for use with the SPC/ Fw water model,49 was used. The SCI-MS-EVB) algorithm44 was applied to perform the simulations of the multiple excess protons in these systems. The DREIDING force field45 was used for the sodium ion parameters. The suite of demanding reactive MD simulations reported in this paper consumed approximately 100 million CPU hours on the IBM Blue Gene/Q supercomputer “Mira” at the Argonne National Laboratory Leadership Computing Facility.

Figure 1. Initial structures for the PSFA bundle (a, c) and lamellar (b, d) morphologies. The top panels look along the z-axis while the bottom panels look along the x-axis. Backbone atoms are in red, sidechain atoms are in yellow, and water is not shown for clarity but was placed around the side chains for the simulations.

using the density of Nafion and the mass of the chains. Every 10 carbons had a side chain attached. The atomic positions of the side chains relative to the attaching backbone carbon were taken from a side chain in an equilibrated random morphology. Periodic boundary conditions mean that the two ends of an individual chain are adjacent to each other, so the positioning of the chains in the z-direction was randomized so that there would be no “joint” created in the lamellar sheet. In the bundle morphology, six chains were positioned in a spiral around a z-axis cylinder of radius 10 Å, such that the spiral completed one rotation from the bottom of the cylinder to the top, again placing the two ends of an individual chain adjacent to each other. A total of six more chains were placed in a spiral around a concentric cylinder of radius 15 Å which was rotated by 30° relative to the inner cylinder so the outer chains lay between the inner chains. Again, every 10 carbons had a side chain attached which pointed radially away from the center of the cylinders, the atomic positions of the side chains relative to the attaching backbone carbon were taken from a side chain in an equilibrated random morphology, and the starting position of the chains were randomized in the z-direction. Two different hydration levels of λ = 10 and 15 were simulated for each morphological system. All simulations were performed using the LAMMPS42 MD package, with our inhouse module used for the MS-EVB calculations.43 The systems were first relaxed using nonreactive constant NPT dynamics at 300 K and 1 atm for 6 ns. This was followed by 12 ns of nonreactive simulation in the constant NVT ensemble. The straight chain systems were equilibrated under a further 20 ns of nonreactive constant NVT ensemble due to their much greater size than the random morphology. Finally, 4 ns of nonreactive constant NVT simulation was performed and configurations sampled every 1 ns to use as the starting configurations for the four production reactive MD simulations. Equilibration of the systems using the reactive self-consistent iterative MS-EVB44 (SCI-MS-EVB) methodology was next

■

RESULTS AND DISCUSSION Structural Features of Equilibrated PFSA Morphologies. Shown in Figure 2 are the equilibrated morphologies for the bundle and lamellar PFSA morphologies. The structures built are stable, with no dissolution of the morphology into random arrangements. The one large difference that is seen is the flattening of the bundle morphology in the radial dimensions, which can be seen most clearly in Figure 2a. In fact, initial attempts to build bundle morphologies were done without using spiraled bundles. The flattening of these straight bundle structures resulted in the periodic images combining into structures very similar to the lamellar structures already built. Examples of these structures are shown in the Supporting Information. In order to explore a wider range of morphological structures, and inspired by the figures in the rod-bundle morphology study,6 the bundles were spiraled to approximate polymer entanglement. This flattening is interesting as it validates several theories postulated from SAXS studies. The initial rod-bundle study proposed a ribbon-like form for the rod-bundles.5 In addition, as outlined in the Introduction, 3178

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C

Table 1. Densities of the Equilibrated Morphologies in g/ cm3 at Different Hydration Levels bundle lamellae random

λ = 10

λ = 15

1.66 1.68 1.61

1.62 1.63 1.53

the best agreement with these experimental values. The trend in density also follows well the level of chain entanglement, with the random morphology showing the lowest densities and the lamellar model showing the highest densities. The chains in the lamellar structure are the most free as there is no physical entanglement of the polymer chains. The bundle morphology chains are simply twisted together, while the random morphology chains are extremely entangled. Clearly the less entangled polymer chains can pack together more efficiently to lower their density. Proton Diffusion Coefficients. Shown in Figure 4 are the self-diffusion coefficients for the hydrated excess proton center Figure 2. Equilibrated structures for the bundle (a, c) and lamellar (b, d) morphologies at λ = 15 (λ = 10 structures look similar). The top panels look along the z-axis while the bottom panels look along the xaxis. Backbone atoms are in red, side-chain atoms are in yellow, and water/hydronium atoms are in blue.

Kreuer and co-workers have already proposed on electrostatic grounds that the morphology should be locally flat.14 The rod bundle flattening has an interesting effect on the water channels in the bundle morphology also. Given the initial rod-like shape of the bundle, the water formed continuous channels in both the x- and y-directions. With the flattening of the rods into ribbons, the flatter part of the bundles can interact with its own periodic image, cutting off the water channel at that point. Since the polymer chains are spiraled, the flat part of the bundle corkscrews along the bundle length, cutting off the water channel in the x- and y-direction at various lengths along the bundle. This results in a tortuous yet locally flat shape to the hydrophilic domains in the bundle morphology, as seen in Figure 3. Density. As a simple check of the representability of these structures, we can examine the density of the equilibrated structures, as shown in Table 1. Experimental measurements show that the density of Nafion is 1.6−1.7 g/cm3 at these hydration levels.50 The two straight chain morphologies show

Figure 4. Self-diffusion constants for the hydrated proton CEC shown in units of 10−7 cm2/s plotted as a function of hydration level. Experimental data were obtained from ref 51 and included in panel a. Panel b shows simulation results alone for ease of comparison.

of excess charge (CEC) as defined in our previous publications.21−31,34 These diffusion coefficients were obtained from the long-time slope of the CEC mean-squared displacement. The diffusion coefficients for the three morphologies show different values. The random morphology diffusion coefficients match well with previous reactive MD simulations of random morphologies.25,26 The change to more regular morphologies increases the diffusion coefficients significantly. The bundle morphology has diffusion coefficients about 1.5 times as large and the lamellar morphology has diffusion coefficients about twice as large as those of the random morphology. When compared to the experimental diffusion coefficients obtained from conductivity experiments51 using the Nernst−Einstein equation, we see that the calculated diffusion coefficients are still off by a factor of 2−4. However, given that the MS-EVB3 model used here underestimates the hydrated proton self-diffusion coefficient in bulk water by a factor of 2−3 due to neglect of nuclear quantum effects and other possible factors, this mismatch is not too surprising. Moreover, the relationship of the proton self-diffusion to the experimental conductivity measurements is not straightforward (for example, assumptions in the latter of ideal solution behavior). However, the large relative difference in diffusion coefficients when we change the morphology is the most interesting feature here, and it is important to further elaborate on these differences in the next few subsections of the paper. Water Diffusion Coefficients. Given that the SPC/Fw water model used in this study matches very well the bulk water diffusion,49 we can be confident of our results on the diffusion

Figure 3. Water in the bundle morphology at λ = 15 with the polymer removed to highlight the channels formed: (a) looks in the x direction; (b) looks in the y direction. The original system is replicated twice on either side to show the connectivity of the water channels. 3179

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C coefficient of water. Since the pulsed-field gradient spin-echo NMR technique used to determine the water diffusion coefficient cannot differentiate the protons from hydronium and water molecules, the diffusion coefficients calculated in this study are derived from the long-time slope of the mean-squared displacement (MSD) of all water and hydronium molecules. Figure 5 shows these diffusion coefficients as a function of

order in that model. The trend in the height of the first peak in the RDF corresponds well with the diffusion coefficient; however, we showed in previous work the RDF between sulfur and the hydrated proton CEC is generally not a good predictor of transport properties. The question then remains how does changing the morphology of the system from a random morphology allow the excess protons to associate less with the sulfonate groups and dissociate faster.

Figure 5. Self-diffusion constants for all water and hydronium oxygens shown in units of 10−7 cm2/s plotted as a function of water volume fraction. Experimental data were obtained from ref 51.

Figure 6. RDFs between the center of excess charge (CEC) of the hydrated excess protons and the sulfur of the sulfonate groups compared among the three morphologies for (a) λ = 10 and (b) λ = 15.

water volume fraction for ease of comparison to experimental results from ref 51. To calculate the water volume, we calculated the Voronoi tessellation52 of all atoms in the system. The volume of water can be calculated as the sum of the volumes of the Voronoi cells of all water and hydronium atoms. The water volume fraction is the volume of water over the volume of the simulation box. The diffusion coefficients from the bundle and random morphologies show the best match to experiment. The lamellar morphology shows much too high diffusion coefficients compared to experiment, suggesting that this morphology is not representative of the true physical morphology of Nafion. It should be noted here that the comparison to experiment is complicated by the correct delay time to use.53 Using longer time delays (1 s compared to ∼10 ms) for the NMR experiment results in slower, limiting diffusion times. However, given the length of our simulations, we compare our simulation results to experiments with the shorter delay times. Since a sizable portion of the hydrated excess proton CEC diffusion is related to the water diffusion coefficient, we can surmise that a large proportion of the increase in lamellar CEC diffusion is caused by the larger water diffusion rate. The random and bundle morphologies show similar water diffusion rates, especially at λ = 15, so the larger proton CEC diffusion coefficient of the bundle morphology seen in Figure 4 is presumably due to greater Grotthuss mechanism activity in the bundle morphology. Radial Distribution Functions. In order to understand how the change in morphology affects overall proton transport, we first start by examining the local environment of the excess proton. The radial distribution function (RDF) between sulfur and the hydrated proton CEC tells us the arrangement of protons around the negatively charged sulfonate group. As we have seen previously, this RDF shows a broad peak between 3 and 6.5 Å, which can be further resolved into two peaks corresponding to the two hydration shells around the sulfonate group. As in our previous paper,26 the RDF drops off to background levels beyond 6.5 Å since the proton spends negligible time outside the second solvation shell of sulfonate groups. The lamellar model shows a slightly higher RDF than the other models beyond 6.5 Å due to the higher long-range

Surface Area. Given that the RDF analysis shows protons are within two solvation shells of sulfonates on the hydrophobic/hydrophilic interface at all times, it is important to understand the features of this interface. One way to accomplish this is to examine the surface area of the hydrophobic/hydrophilic interface. To estimate the surface area, we calculated the Voronoi tessellation52 of all atoms in the system. The surface area of an interface between a region A and a region B can then be calculated as the sum of the areas of the Voronoi faces between atoms in A and atoms in B. In addition, the volume of a region can be calculated as the sum of the volumes of the Voronoi cells of all atoms in the region. Given that Porod analysis to determine surface area from scattering experiments has a resolution of at most 10 nm−1, the side chain surface area will not be resolved. Thus, we determined the surface area of the backbone atoms with the water and hydronium atoms, though this may underestimate the surface area somewhat. Plotted in Figure 7 is the surface area per

Figure 7. (a) Surface area of hydrophobic/hydrophilic interface per sulfonate group as a function of water volume fraction. (b) Surface area of hydrophobic/hydrophilic interface over the volume of simulation box as a function of water volume fraction.

sulfonate group, which gives an idea of the spacing of the sulfonate groups along the interface. Given the possible underestimation of the surface area, it is difficult to definitively say which morphology shows best agreement with experiment; however, there is general agreement from all morphologies. Interestingly, the trend in surface area corresponds well with the proton and water diffusion coefficients, further corroborat3180

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C

two sulfonates. However, when we expand the analysis to include two solvation shells, we can see a more varied behavior. At λ = 10, we see that the lamellar and bundle morphologies show a greater probability for having protons being associated with four sulfonate groups. We have shown already in previous work that this degree of “sharing” was correlated with higher diffusion rates for the 3M membrane.26 Water Clusters. We can also examine the clustering of the water in the hydrophilic domains. In these simulations, however, since we are examining only high hydration level systems, the clustering of the water domains is extensive, with the vast majority of waters in one large cluster at all times. Figure 9 shows the cumulative probability distribution for the

ing the theory that the interface is the most important region for proton transport. Similar to the density of the systems, the surface area of the systems correlates well with the entanglement of the systems, with the entangled random system showing the highest surface areas. This suggests that the value of 40 Å2 per sulfonate group of the lamellar morphology is the most favorable surface area per sulfonate group, and the entanglement of the polymer chains in the bundle and random morphologies prevents them from reaching this low surface area. Association to Sulfonate Groups. A useful metric to determine the population of waters and hydrated protons in various parts of the hydrophilic pore is to count the number of sulfonate groups to which the water or hydrated proton is associated (Figure 8). We can define association as being within

Figure 9. Cumulative probability distribution of the size of the water clusters in each morphology. The cluster sizes are scaled by the maximum cluster size possible in each simulation to aid comparison between systems with different numbers of waters. This figure only shows greater than 0.95 cluster size for clarity, while the inset shows the full range. Figure 8. Number of sulfonate groups to which water or hydrated excess protons are associated. Association is defined as either being within a cutoff of 4.5 Å (a, c) or a cutoff of 6.5 Å (b, d). These cutoffs are chosen to select the first and second solvation shells, respectively. Λ = 10 are solid lines, and λ = 15 are dashed lines.

clustering of the water domains. Since the straight chain morphologies have 3 times more water than the random morphology, it was necessary to scale the results by the number of waters in each system (i.e., the maximum cluster size) to aid comparison between the systems. The results show that the lamellar morphology has uniformly good water clusters, with even λ = 10 simulations showing similar levels of water clustering to λ = 15. Bundle and random morphologies show comparable levels of water clustering at both hydration levels. Overall, the results correlate well with the water diffusion coefficients seen in Figure 4. Sodium Simulations. In order to further probe the diffusive behavior of ions in PFSAs, we performed simulations of the three morphologies with varying fractions of the excess protons replaced with sodium ions, as outlined in the Simulation Details section. Figure 10 shows the self-diffusion

one solvation shell or two solvation shells and so get a sense of what percentage of water and hydrated protons are residing in each of these populations. By examining waters first, we can see that as the hydration level increases the average number of sulfonate groups to which a water molecule is associated decreases, and the proportion of waters that are unassociated with any sulfonate groups increases. All morphologies exhibit similar behavior when we look at the first solvation shell, especially at λ = 15. When we look at the second solvation shell, however, we can see that at λ = 15 the lamellar morphology shows a difference with the other two morphologies. While 20% of waters are outside the second solvation shell of all sulfonate groups in the random and bundle morphologies, only 10% of waters are outside the second solvation shell of sulfonate groups in the lamellar morphology. This is not what one would expect from the higher water diffusion coefficients in the lamellar morphology seen in Figure 5, since water in the bulk like center of the pore should presumably diffuse faster than those at the interface. Examining the hydrated proton coordination, we can see that the behavior for the first solvation shell is again consistent between morphologies and for the proton is consistent between λ = 10 and λ = 15 as well. It shows that hydrated protons spend about one-third of time within one solvation shell of a sulfonate and about two-thirds of time outside the first solvation shell of all sulfonates, with a small probability of being shared between

Figure 10. Self-diffusion constants for the sodium ions shown in units of 10−7 cm2/s plotted as a function of hydration level. (a) Comparison with experimental measurements from ref 37. (b) Hydrated excess proton CEC diffusion constants from Figure 4 are shown as dashed lines and shifted slightly left for ease of comparison. 3181

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C

decreases while the self-diffusion coefficient of sodium increases. This effect is most pronounced in the lamellar and bundle morphologies, while the random morphology shows more mixed effects. In the random morphology, the sodium diffusion coefficient is clearly larger in the χp = 0.5 system than in the χp = 0, but not so for the χp = 0.25 system. In addition, the hydrated excess proton diffusion coefficient cannot be statistically differentiated between the different χp systems in the random morphology. The next sections will concentrate on elucidating the reasons for the slowdown in proton diffusion in conjunction with the speed-up in the sodium diffusion as the two ions are mixed. Water Structure. Figure 13 shows the RDF between the oxygens of all waters in the system compared between all

constants for the sodium ions in the system with 100% sodium ions, obtained from the long time slope of the MSDs. Similar to the water diffusion constants, we see a better match to experimental for the sodium diffusion constants in the bundle and random morphologies compared to the lamellar morphology. We can also see that the sodium diffusion coefficients are statistically identical to the excess proton diffusion coefficients (already seen in Figure 4) for all λ = 10 systems but lower than the excess proton diffusion coefficients at λ = 15. Reactive excess proton simulations show a large increase in diffusion coefficients compared to nonreactive simulations, so it is interesting that sodium ions are able to transport at similar rates to reactive hydroniums in our simulations. However, as mentioned previously, given that the MS-EVB3 model used here underestimates the hydrated proton self-diffusion coefficient in bulk water (likely due to some nuclear quantum effects), this does not contradict the experimental observation that membranes with protons show higher conductivity than membranes with sodium ions. In addition, when we look at the S−Na RDFs in Figure 11 we can see that the sodium ions are

Figure 13. RDFs between oxygens of water molecules compared among the four sodium mixtures for the bundle morphology at λ = 10.

sodium mixture simulations. For conciseness and since all systems show comparable behavior, we have arbitrarily chosen the bundle morphology at λ = 10 to display. All other systems are shown in the Supporting Information for completeness. The 100% H+ system (or χp = 1) shows the classic water−water RDF one would expect, with a large first solvation shell peak at about 3 Å and a second solvation shell peak at about 4 Å. The introduction of sodium ions to the system has very little effect on the first solvation shell peak but flattens the second solvation shell peak, indicating sodium is having a significant impact in the structure of water in the system. Hydrogen Bond Network Disruption. We can extend the water clustering methodology used in Figure 9 to examine hydrogen bonding between waters and hydrated protons in

Figure 11. RDFs between sodium ions and the sulfur of the sulfonate group compared among the three morphologies for (a) λ = 10 and (b) λ = 15.

much more attracted to the sulfonate groups than the hydrated protons, so one might assume that they would diffuse slower than the excess protons, again highlighting the importance of the passing between sulfonates diffusion mechanism. Sodium Mixture Simulations. Looking at the diffusion coefficients for the simulations with both sodium ions and excess protons in Figure 12, one can observe an interesting trend emerge. As the fraction of cations that are protons (χp) is decreased, the self-diffusion coefficient of the excess proton

Figure 12. Self-diffusion constants for all excess proton CEC and sodium shown in units of 10−7 cm2/s plotted as a function of χp, the fraction of ions that are protons, for λ = 10 in the top row and λ = 15 in the bottom row. 3182

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C each system, which provides more understanding of the shortrange disturbances in the water system. Instead of only using a distance cutoff between neighboring water or hydrated proton oxygens to define whether they are within the same cluster, we can use whether these oxygens are hydrogen bonded to each other (using the geometric definition of a hydrogen bond, where the O−H−O angle should be below 30° and the O−O distance should be below 3.5 Å) to define whether they are in the same cluster. This method produces a figure quite similar to Figure 9, as shown in Figure 14a, with all systems showing

Figure 15. Self-diffusion constants for all water and hydrated proton oxygens shown in units of 10−7 cm2/s plotted as a function of χp, the fraction of ions that are protons for (a) λ = 10 and (b) λ = 15.

and continuous components to get a sense of the contribution the Grotthuss mechanism to overall proton transport rCEC(t ) − rCEC(0) = ΔrCEC(t ) = Δrc(t ) + Δrd(t )

(1)

⟨|ΔrCEC(t )|2 ⟩ = ⟨|Δrc(t )|2 ⟩ + ⟨|Δrd(t )|2 ⟩ + 2⟨Δrc(t )Δrd(t )⟩ (2)

where Δrc(t) and Δrd(t) are understood here to be vectors and they are respectively the “continuous” and “discrete” components of the total displacement of the hydrated proton CEC, ΔrCEC(t). We have used the decomposition method outlined in detail in ref 25, where each 100 fs proton displacement is assigned to the discrete component if there has been a change in the identity of the most probable (“pivot”) hydronium during that 100 fs or is assigned to the continuous component otherwise. Figure 16 shows the diffusion

Figure 14. (a) Cumulative probability distribution of the size of the hydrogen bonded water clusters compared among the four sodium mixtures for the bundle morphology at λ = 10. The cluster sizes are scaled by the maximum cluster size possible in each simulation to aid comparison between systems with different numbers of waters. (b) Median hydrogen bonded water cluster size for each morphology as a function of χp, the fraction of ions that are protons. Solid lines are λ = 15, and dashed lines are λ = 10.

large-scale connectivity of the hydrogen bond network. There is a small increase in the size of clusters ∼4 waters in size when compared to distance based clustering, which can be seen in snapshots of the systems to be Eigen cations (H9O4+) which are not hydrogen bonded to the neighboring waters. As more sodium is added to the systems, however, the extent of the hydrogen bond network is decreased, with the cumulative probability lines shifting left, as seen in Figure 14a. In order to show this change more clearly across systems, we chose a cumulative probability of 0.50 to be our reference point, since this defines the median of the distribution, and determined the scaled cluster size at this point for all systems. Figure 14b shows this median scaled cluster size plotted for each system at all mixtures of sodium and hydrated protons. The median scaled cluster size decreases for all systems as sodium is added to the system, with λ = 10 systems more affected by the addition of sodium than the λ = 15 systems, but with similar levels of change between the different morphological models. As sodium is added to the system, both water structure and hydrogen bond structure decreases. This would suggest that vehicular diffusive transport should get faster in this less ordered structure which was seen in other simulations of ion mixtures.35 However, water diffusion coefficients are essentially unchanged at all values of χp as shown in Figure 15, and sodium ions diffuse faster at higher χp, as water structure and hydrogen bonding increases. It is possible, however, that the excess proton diffusion is being affected by this change in hydrogen bonding, as the Grotthuss mechanism depends on hydrogen bonding to occur. We can further investigate this by decomposing the motion of the hydrated excess CEC into vehicular and hopping components. Decomposition of the Proton Diffusion. As shown in previous work on reactive proton transport in PFSAs,22,25 we can decompose the displacement of the proton into discrete

Figure 16. Diffusion coefficients for the excess proton mean-squared displacement that has been decomposed into its discrete and continuous components at (a) λ = 10 and (b) λ = 15. The thick lines are the diffusion coefficients from the normal MSD, the dashed lines are the diffusion coefficients from the discrete MSD component, and the thin lines are the diffusion coefficients from the continuous component.

coefficients obtained from the long time slope of the decomposed continuous and discrete MSDs, ⟨|Δrc(t)|2⟩ and ⟨|Δrd(t)|2⟩. Both decomposed diffusion coefficients are much higher than the diffusion coefficients from the total displacement due to the effect of the correlation between the continuous and discrete displacements, ⟨Δrc(t)Δrd(t)⟩. Overall, however, the expected decrease in the discrete displacement in relation to the continuous displacement due to the decrease in hydrogen bonding does not occur. Instead, we see that the discrete and continuous components of the displacement change in step with each other. Radial Distribution Functions of the Ions. Previous work on fluoride−chloride ion co-ion effects in polyelectrolytes showed that some of the change in fluoride dynamics could be attributed to a change in the attraction of fluoride to the counterion with the addition of chloride.34,35 Figure 17 shows the radial distribution functions for sodium and the oxygen of the pivot hydronium with the sulfur of the sulfonate group compared between the different mixture simulations. Again, for 3183

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C

and decreasing ion interactions, we can now correlate a decrease in water structure to a decrease in ion diffusion. This is a surprising result when compared to previous work, which showed that while chloride ions were also associated with decreasing water structure, this disordered water structure was associated with faster water and ion mobility.

■

CONCLUSIONS In this study we have examined two prevailing models for the morphology of hydrated Nafion and compared them to the random morphology model generally used for atomistic simulations. Structurally, these two morphologies consisted of long straight chains of Nafion arranged in a rod-like bundle or a lamellar sheet. After equilibration, the rod-bundle morphology tended to flatten out into ribbon like bundles to form locally flat hydrophobic/hydrophilic interfaces and locally flat hydrophilic regions. When we compared properties obtained from experimental measurements to those calculated from our simulations, we found the bundle morphology showed the best agreement, with the lamellar model, as the extreme extension of the locally flat model, showing too low surface area and too fast water diffusion. The random model has too much polymer entanglement due to it not representing the persistence length of Nafion correctly, and thus its density is too low and water diffusion is too slow. In agreement with previous work,26 we see that transport along the hydrophobic/hydrophilic interface is the primary mechanism for proton transport. Proton diffusion coefficients are highly correlated with the surface area of the hydrophobic/ hydrophilic interface, with lower surface areas showing higher diffusion coefficients. Given that morphology in PFSA membranes can be changed by various processing methods such as extrusion and annealing, methods that can reduce the surface area of the hydrophobic/hydrophilic interface could increase the proton transport performance of the membrane. Mixtures of sodium ions with excess protons exhibited interesting results, where the addition of sodium tended to decrease water structure and hydrogen bonding, which caused less association of the cations with the sulfonate group and slowed the diffusion of the cations. This is in contrast to previous studies of mixtures of chloride ions with other ions, which showed that the decrease in water structure and hydrogen bonding caused by chloride ions was associated with an increase in the diffusion of ions. When we also take into account the lower diffusion coefficient of the hydrated excess proton in our simulations when compared to experiment, this unusual co-ion effect helps explain the decrease in ionic

Figure 17. Radial distribution functions between the sulfur of the sulfonate group and the (a) sodium ion and (b) oxygen of the hydrated proton for the bundle morphology at λ = 10 compared between different sodium mixture simulations.

conciseness and since all systems show comparable behavior, we have chosen the bundle morphology at λ = 10 to display. All other systems are shown in the Supporting Information for completeness. We can see in Figure 17b that the height of the RDF of the hydrated proton to the sulfonate group is decreased as the number of excess protons in the system is decreased, as one would expect. The much greater overall attraction of the sodium to the sulfonate group compared to the hydrated proton would mean the addition of sodium tends to drive out the hydrated protons from the sulfonate group. However, in Figure 17a, we can see that the height of the RDF of the sodium to the sulfonate group barely changes as the number of sodium ions in the system is decreased. We can see further evidence of this counterintuitive behavior in the RDFs between the cations in Figure 18. In Figure 18b we see again that the hydrated proton interaction decreases as the number of excess protons in the system decreases, as one would expect. However, Figure 18a shows the interaction between sodium ions increases as the number of sodium ions in the system decreases. Finally, Figure 18c shows the interaction of hydrated protons and sodium ions increases with the number of excess protons in the system (i.e., increasing χp). This helps to introduce a new way of interpreting these RDFs in Figures 17 and 18. All ion interactions become more structured with an increase in χp, which can be tied to the increase in water structure with χp described earlier. This increase in water structure makes disassociation of the ions from the sulfonates more difficult, increasing association with the sulfonates and hence increasing association with each other. This perspective forces a reconsideration of the results (Figure 12). Instead of the obvious observation that the increase in diffusion constant of sodium and the decrease in diffusion constant of the hydrated proton CEC upon mixing of the ions, we can instead see that there is a slowing of both ion diffusion constants with decreasing χp. Given that we have already associated decreasing χp with decreasing water structure

Figure 18. Radial distribution functions between (a) sodium ions, (b) hydrated proton oxygens, or (c) sodium ions and hydronium oxygens for the bundle morphology at λ = 10 compared between different sodium mixture simulations. 3184

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C

(12) Loppinet, B.; Gebel, G.; Williams, C. E. Small-Angle Scattering Study of Perfluorosulfonated Ionomer Solutions. J. Phys. Chem. B 1997, 101, 1884−1892. (13) Aldebert, P.; Dreyfus, B.; Pineri, M. Small-Angle NeutronScattering of Perfluorosulfonated Ionomers in Solution. Macromolecules 1986, 19, 2651−2653. (14) Kreuer, K. D.; Portale, G. A Critical Revision of the NanoMorphology of Proton Conducting Ionomers and Polyelectrolytes for Fuel Cell Applications. Adv. Funct. Mater. 2013, 23, 5390−5397. (15) 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. (16) Jorn, R.; Savage, J.; Voth, G. A. Proton Conduction in Exchange Membranes across Multiple Length Scales. Acc. Chem. Res. 2012, 45, 2002−2010. (17) Daly, K. B.; Panagiotopoulos, A. Z.; Debenedetti, P. G.; Benziger, J. B. Viscosity of Nafion Oligomers as a Function of Hydration and Counterion Type: A Molecular Dynamics Study. J. Phys. Chem. B 2014, 118, 13981−13991. (18) Lucid, J.; Meloni, S.; MacKernan, D.; Spohr, E.; Ciccotti, G. Probing the Structures of Hydrated Nafion in Different Morphologies Using Temperature-Accelerated Molecular Dynamics Simulations. J. Phys. Chem. C 2013, 117, 774−782. (19) Knox, C. K.; Voth, G. A. Probing Selected Morphological Models of Hydrated Nafion Using Large-Scale Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 3205−18. (20) Knight, C.; Voth, G. A. The Curious Case of the Hydrated Proton. Acc. Chem. Res. 2012, 45, 101−109. (21) Petersen, M. K.; Wang, F.; Blake, N. P.; Metiu, H.; Voth, G. A. Excess Proton Solvation and Delocalization in a Hydrophilic Pocket of the Proton Conducting Polymer Membrane Nafion. J. Phys. Chem. B 2005, 109, 3727−30. (22) Petersen, M. K.; Voth, G. A. Characterization of the Solvation and Transport of the Hydrated Proton in the Perfluorosulfonic Acid Membrane Nafion. J. Phys. Chem. B 2006, 110, 18594−600. (23) Feng, S.; Voth, G. A. Proton Solvation and Transport in Hydrated Nafion. J. Phys. Chem. B 2011, 115, 5903−12. (24) Feng, S. L.; Savage, J.; Voth, G. A. Effects of Polymer Morphology on Proton Solvation and Transport in Proton-Exchange Membranes. J. Phys. Chem. C 2012, 116, 19104−19116. (25) Tse, Y. L. S.; Herring, A. M.; Kim, K.; Voth, G. A. Molecular Dynamics Simulations of Proton Transport in 3M and Nafion Perfluorosulfonic Acid Membranes. J. Phys. Chem. C 2013, 117, 8079− 8091. (26) Savage, J.; Tse, Y. L. S.; Voth, G. A. Proton Transport Mechanism of Perfluorosulfonic Acid Membranes. J. Phys. Chem. C 2014, 118, 17436−17445. (27) Savage, J.; Voth, G. A. Persistent Sub-Diffusive Proton Transport in Perfluorosulfonic Acid Membranes. J. Phys. Chem. Lett. 2014, 5, 3037−3042. (28) Schmitt, U. W.; Voth, G. A. Multistate Empirical Valence Bond Model for Proton Transport in Water. J. Phys. Chem. B 1998, 102, 5547−5551. (29) Schmitt, U. W.; Voth, G. A. The Computer Simulation of Proton Transport in Water. J. Chem. Phys. 1999, 111, 9361−9381. (30) Day, T. J. F.; Soudackov, A. V.; Cuma, M.; Schmitt, U. W.; Voth, G. A. A Second Generation Multistate Empirical Valence Bond Model for Proton Transport in Aqueous Systems. J. Chem. Phys. 2002, 117, 5839−5849. (31) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. An Improved Multistate Empirical Valence Bond Model for Aqueous Proton Solvation and Transport. J. Phys. Chem. B 2008, 112, 467−82. (32) Agmon, N. The Grotthuss Mechanism. Chem. Phys. Lett. 1995, 244, 456−462. (33) Grotthuss, C. J. T. Sur la Decomposition de l’Eau et des Corps qu’Elle Tient en Dissolution a l’Aide de l’Électricité Galvanique. Ann. Chim. 1806, 58, 54−54.

conductivity seen in experimental measurements of sodium and proton mixtures in PFSA membranes.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b11168. Comparison of the equilibrated straight chain structures, RDFs between oxygens of water molecules, RDFs between sulfurs of the sulfonate group and the sodium ions, RDFs between sulfurs of the sulfonate group and the hydrated protons, RDFs between the sodium ions, and RDFs between the hydrated protons (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Phone 773-702-9092; Fax 773-702-9106; e-mail gavoth@ uchicago.edu (G.A.V.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the Department of Energy (DOE), Office of Basic Energy Sciences (BES), Division of Chemical Sciences, Geosciences, and Biosciences, through Award DE-SC0005418. An allocation of computer time was provided by the ASCR Leadership Computing Challenge (ALCC) Program at the Argonne Leadership Computing Facility (ALCF) at Argonne National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC02-06CH11357.

■

REFERENCES

(1) Hamrock, S. J.; Yandrasits, M. A. Proton Exchange Membranes for Fuel Cell Applications; Taylor & Francis: New York, 2006; pp 219− 244. (2) Mauritz, K. A.; Moore, R. B. State of Understanding of Nafion. Chem. Rev. 2004, 104, 4535−85. (3) Gierke, T. D.; Munn, G. E.; Wilson, F. C. The Morphology in Nafion Perfluorinated Membrane Products, as Determined by Wideand Small-Angle X-Ray Studies. J. Polym. Sci., Polym. Phys. Ed. 1981, 19, 1687−1704. (4) Schmidt-Rohr, K.; Chen, Q. Parallel Cylindrical Water Nanochannels in Nafion Fuel-Cell Membranes. Nat. Mater. 2008, 7, 75−83. (5) Rubatat, L.; Rollet, A. L.; Gebel, G.; Diat, O. Evidence of Elongated Polymeric Aggregates in Nafion. Macromolecules 2002, 35, 4050−4055. (6) Rubatat, L.; Gebel, G.; Diat, O. Fibrillar Structure of Nafion: Matching Fourier and Real Space Studies of Corresponding Films and Solutions. Macromolecules 2004, 37, 7772−7783. (7) Haubold, H. G.; Vad, T.; Jungbluth, H.; Hiller, P. Nano Structure of Nafion: A SAXS Study. Electrochim. Acta 2001, 46, 1559−1563. (8) Starkweather, H. W. Crystallinity in Perfluorosulfonic Acid Ionomers and Related Polymers. Macromolecules 1982, 15, 320−323. (9) Chu, B.; Wu, C.; Buck, W. Light-Scattering Characterization of Poly(Tetrafluoroethylene). 2. PTFE in Perfluorotetracosane Molecular-Weight Distribution and Solution Properties. Macromolecules 1989, 22, 831−837. (10) Rosi-Schwartz, B.; Mitchell, G. R. Extracting Force Fields for Disordered Polymeric Materials from Neutron Scattering Data. Polymer 1996, 37, 1857−1870. (11) Aldebert, P.; Dreyfus, B.; Gebel, G.; Nakamura, N.; Pineri, M.; Volino, F. Rod Like Micellar Structures in Perfluorinated Ionomer Solutions. J. Phys. (Paris) 1988, 49, 2101−2109. 3185

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186

Article

The Journal of Physical Chemistry C (34) Tse, Y. L. S.; Sarode, H. N.; Lindberg, G. E.; Witten, T. A.; Yang, Y.; Herring, A. M.; Voth, G. A. Chloride Enhances Fluoride Mobility in Anion Exchange Membrane/Polycationic Systems. J. Phys. Chem. C 2014, 118, 845−853. (35) Tse, Y. L.; Voth, G. A.; Witten, T. A. Ion Mixing, Hydration, and Transport in Aqueous Ionic Systems. J. Chem. Phys. 2015, 142, 184905. (36) Okada, T.; Möller-Holst, S.; Gorseth, O.; Kjelstrup, S. Transport and Equilibrium Properties of Nafion (R) Membranes with H+ and Na+ Ions. J. Electroanal. Chem. 1998, 442, 137−145. (37) Chaudhury, S.; Agarwal, C.; Pandey, A. K.; Goswami, A. SelfDiffusion of Ions in Nafion-117 Membrane Having Mixed Ionic Composition. J. Phys. Chem. B 2012, 116, 1605−11. (38) Okada, T.; Arimura, N.; Satou, H.; Yuasa, M.; Kikuchi, T. Membrane Transport Characteristics of Binary Cation Systems with Li + and Alkali Metal Cations in Perfluorosulfonated Ionomer. Electrochim. Acta 2005, 50, 3569−3575. (39) Saito, M.; Arimura, N.; Hayamizu, K.; Okada, T. Mechanisms of Ion and Water Transport in Perfluorosulfonated Ionomer Membranes for Fuel Cells. J. Phys. Chem. B 2004, 108, 16064−16070. (40) Daly, K. B.; Benziger, J. B.; Debenedetti, P. G.; Panagiotopoulos, A. Z. Molecular Dynamics Simulations of Water Sorption in a Perfluorosulfonic Acid Membrane. J. Phys. Chem. B 2013, 117, 12649− 12660. (41) Feng, S.; Savage, J.; Voth, G. A. Effects of Polymer Morphology on Proton Solvation and Transport in Proton-Exchange Membranes. J. Phys. Chem. C 2012, 116, 19104−19116. (42) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular-Dynamics. J. Comput. Phys. 1995, 117, 1−19. (43) Yamashita, T.; Peng, Y.; Knight, C.; Voth, G. A. Computationally Efficient Multiconfigurational Reactive Molecular Dynamics. J. Chem. Theory Comput. 2012, 8, 4863−4875. (44) Wang, F.; Voth, G. A. A Linear-Scaling Self-Consistent Generalization of the Multistate Empirical Valence Bond Method for Multiple Excess Protons in Aqueous Systems. J. Chem. Phys. 2005, 122, 144105. (45) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94, 8897−8909. (46) Zawodzinski, T. A.; Springer, T. E.; Davey, J.; Jestel, R.; Lopez, C.; Valerio, J.; Gottesfeld, S. A Comparative-Study of Water-Uptake by and Transport through Ionomeric Fuel-Cell Membranes. J. Electrochem. Soc. 1993, 140, 1981−1985. (47) Frisch, M. J.; et al. Gaussian 09, Revision B.01; Gaussian Inc.: Wallingford, CT, 2009. (48) Jang, S. S.; Molinero, V.; Cagin, T.; Goddard, W. A. NanophaseSegregation and Transport in Nafion 117 from Molecular Dynamics Simulations: Effect of Monomeric Sequence. J. Phys. Chem. B 2004, 108, 3149−3157. (49) Wu, Y.; Tepper, H. L.; Voth, G. A. Flexible Simple Point-Charge Water Model with Improved Liquid-State Properties. J. Chem. Phys. 2006, 124, 024503. (50) Morris, D. R.; Sun, X. Water-Sorption and Transport Properties of Nafion 117 H. J. Appl. Polym. Sci. 1993, 50, 1445−1452. (51) Kreuer, K. D.; Schuster, M.; Obliers, B.; Diat, O.; Traub, U.; Fuchs, A.; Klock, U.; Paddison, S. J.; Maier, J. Short-Side-Chain Proton Conducting Perfluorosulfonic Acid Ionomers: Why They Perform Better in Pem Fuel Cells. J. Power Sources 2008, 178, 499−509. (52) Rycroft, C. H. Voro++: A Three-Dimensional Voronoi Cell Library in C++. Chaos 2009, 19, 041111. (53) Zhao, Q.; Majsztrik, P.; Benziger, J. Diffusion and Interfacial Transport of Water in Nafion. J. Phys. Chem. B 2011, 115, 2717−2727.

3186

DOI: 10.1021/acs.jpcc.5b11168 J. Phys. Chem. C 2016, 120, 3176−3186