Reactive Molecular Dynamics Study of Proton Transport in Polymer

ARTICLE pubs.acs.org/JPCC

Reactive Molecular Dynamics Study of Proton Transport in Polymer Electrolyte Membranes Myvizhi Esai Selvan,† David J. Keffer,*,†,‡ and Shengting Cui† † ‡

Department of Chemical and Biomolecular Engineering, University of Tennessee, Knoxville, Tennessee 37996-2200, United States Department of Chemical and Biomolecular Engineering, Yonsei University, Seoul 120-749, Republic of Korea ABSTRACT: Dynamical properties of water and protons in Nafion with an equivalent weight of 1144 are studied using the recently developed reactive molecular dynamics (RMD) algorithm at various water contents. The structural diffusion of a proton along the aqueous domains is modeled via a mechanism similar to that observed in bulk aqueous systems. The algorithm implements reactivity in classical MD simulations by three steps: (i) satisfaction of the trigger, (ii) instantaneous reaction, and (iii) local equilibration. Two different schemes (Method 1 and Method 2) of execution of the algorithm are investigated, which differ in terms of the range of the local environment involved in the reaction. Both methods are parametrized to the experimental rate constant of proton transport in bulk water. The mean lifetime of the protons increased with the water content in Nafion. Water diffusivities of the reactive systems using the two schemes increased with the hydration level and were higher than the diffusion coefficients in the nonreactive system. The more detailed scheme in Method 2 which included the relaxation of the extended hydrogen bonding network around the proton-hopping reaction lowered the water diffusivity compared to that of Method 1 and also affected both the structural and vehicular components of diffusion. The total charge diffusion, vehicular component, and structural component increased with water content, and the two components of the charge diffusion were found to be negatively correlated. The negative correlation is due to preferential structural diffusion of the proton toward the sulfonate group, whereas vehicular diffusion tends to move the H3O+ ion away from the sulfonate group.

I. INTRODUCTION Proton exchange membrane fuel cells (PEMFCs) are compact and highly efficient power generators for stationary and transportation applications at low-temperature operating conditions with tremendous potential to reduce oil dependence and carbon emissions.1 The proton exchange membrane (PEM) is the polymer electrolyte employed in a PEMFC and in a hydrated state is responsible for the conduction of proton from anode to cathode. An ideal PEM is expected to have mechanical integrity, chemical stability, long-term durability, high proton conductivity at extreme operating conditions of low humidity or high temperature, and low cost of production. Currently, no existing membrane exhibits all of the above properties. The proton conduction critically depends on the PEM morphology which in turn is influenced by the water content2
5 and polymer chemistry.6
9 Toward this end, a molecular level understanding of the structure
property relationship is important for the improvement not only in performance of the existing PEMs but also in the design of novel PEMs with superior characteristics.10,11 Proton diffusion in bulk water occurs through a combination of both vehicular (conventional mass diffusion as a H3O+ ion) and structural diffusion (translocation of H+ from a water molecule to another through rearrangement of hydrogen bonds)12,13 mechanisms. Structural diffusion is responsible for the anomalous high r 2011 American Chemical Society

mobility of protons in aqueous solutions compared to other cations of the same size.14 The mechanism involves continual interconversion between the two limiting configurations, namely, Zundel and Eigen complexes, accompanied by breaking and forming of hydrogen bonds in the solvation shells around the excess charge.12,15 In a PEMFC, the proton diffuses from the anode to the cathode through the aqueous domain in a PEM. These aqueous clusters differ from that of bulk water because of their dimension in nanometers and the presence of an acidic moiety from the side chain of the PEM. There is evidence from theory and simulation suggesting that both proton transport mechanisms are active in the PEM16
18 with an increase in structural diffusion occurring at high hydration level.7 Experimentally, proton conductivity can be measured through electrochemical impedance spectroscopy,19 DC four-point probe electric conductivity,20 high-throughput screening,21 and nuclear magnetic resonance (NMR) technique.22,23 Diffusion coefficients of the charge can be measured through pulsedfield gradient NMR spectroscopy, but they are lower than the diffusivities calculated from ionic conductivity through the Nernst
Einstein equation due to the failure to capture Received: April 13, 2011 Revised: August 10, 2011 Published: August 10, 2011 18835

dx.doi.org/10.1021/jp203443c | J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C the Grotthuss mechanism.22,23 The experimental measurements do not provide direct information regarding the molecular-level interpretation of the proton transport mechanism. To understand the connection between the structure and transport, the contribution and behavioral change in the two components of proton diffusion within the highly acidic and confined aqueous domains of the membrane need to be studied. Simulation or molecular modeling of proton transport needs to capture the disparate time scales associated with vehicular diffusion (nanoscale) and structural diffusion (pico time scale). Modeling techniques like Car
Parinello ab initio molecular dynamics (CPAIMD),13,24 mixed quantum and classical mechanics technique (QM/MM),25,26 empirical valence bond (EVB) model27 and its family,28
30 the mixed MD/MC algorithm,31 and Q-HOP MD32 made significant contributions toward the understanding of proton transport in bulk water, and some of them have been extended to study proton transport in PEMs.33
37 Other theoretical methodologies like electronic structure calculations,38 ab initio molecular dynamics39,40 and meta-dynamics,41 and nonequilibrium statistical models42,43 have also been implemented. Molecular-level modeling of proton transport in perfluorosulfonic acid (PFSA) membranes has been reviewed by Elliott and Paddison.44 Classical MD simulations5,8,9 have been used to study proton dynamics, but they can only follow the vehicular diffusion since structural diffusion cannot be dictated by classical mechanics. Vehicular diffusion has been studied as a function of polymer chemistry (EW, monomeric sequence), hydration, and temperature. Vehicular diffusion increases with water content and temperature45 but is much lower than experimental values at high λ because of the absence of structural diffusion in the models. There is experimental46 and theoretical evidence43 that separation of sulfonic acid groups can affect the transport properties, and the vehicular component increased with the decrease in EW of the polymer at intermediate and high water contents.8 The influence of the nanophase segregation structure and water content on the transport properties of water molecules and protons in hydrated membranes has been analyzed by studying water-soluble dendrimers grafted onto different linear polymers using classical molecular dynamics simulation with the incorporation of Q-HOP to enable structural diffusion.33,36 It was concluded from the above simulations that an increase in nanophase segregation scale improved proton conductivity.33,36 Q-HOP MD has been used to study the activation energy and mean residence times for proton transfer reaction as a function of hydration and temperature with a single proton undergoing structural diffusion in Nafion. The effect of the water percolation network in Nafion on proton conductivity has also been investigated.37,47 A strong anticorrelation between the structural and vehicular component was found by the application of self-consistent multistate EVB, and the effect of the sulfonate ion on charge diffusion by acting as traps has been analyzed.35 This work uses a newly developed reactive molecular dynamics (RMD) algorithm to analyze proton transport in the PEM, Nafion. The RMD algorithm models the structural diffusion of a proton as a chemical reaction and implements the reactivity through a three-step procedure instead of the manipulation of interactive potential or empirical formulations. The algorithm has been previously applied to study structural diffusion of a proton in bulk water48 and extended to other systems which have the same transition states (TSs) (Zundel and Eigen cations) involved, including aqueous hydrochloric

ARTICLE

acid solutions48 and water confined in carbon nanotubes (CNTs).49 To implement the RMD algorithm one must know a priori the chemical reactions of interest to be modeled. This can be taken as both an advantage and a disadvantage. The disadvantage is that we need to know all the possible reactions that occur, and a separate set of parameters is needed for each reaction. Structural diffusion of protons in a PEM can take place along either the surface of a hydrophobic/hydrophilic interface or the core of the aqueous domain which has a different environment.50 The advantage is that in a complex system like a PEM where confinement can dictate the aqueous environment and acidity can define the hydrophobic/hydrophilic interface understanding the individual contribution of reactions occurring in these regions will help us to relate structure to transport property. In the present work, we have assumed that structural diffusion of a proton away from the interface follows the same mechanism as used in the studies of bulk water, bulk HCl solutions, and water in CNTs. Currently, we have modeled only the structural diffusion of a proton in the PEM via the above mechanism. However, we have implemented two different methods of implementation of the RMD algorithm in the PEM system to understand the effect of the parameters involved in the algorithm on the transport property of the charge. The paper is organized as follows: The methodology, including a brief description of the RMD algorithm, and the details of the simulations are described in Section II. Results and discussion are presented in Section III, and the findings are summarized in Section IV.

II. SIMULATION METHODOLOGY A. Reactive Molecular Dynamics Algorithm. The RMD algorithm48,51 is a generalized algorithm that can represent any chemical reaction via three steps: (i) satisfaction of a set of geometric and energetic triggers, (ii) instantaneous reaction, and (iii) local equilibration. The algorithm takes input from both macroscopic description (stoichiometry, rate constant, k, activation energy, Ea, heat of reaction, ΔHr) and quantum mechanical (QM) description (transition state, reactant, and product structures) of a chemical reaction and reproduces a model that is far more detailed than a macroscopic model and less detailed than the QM model. The RMD algorithm is implemented at the end of each classical MD step. The model has an advantage of using the existing nonreactive potentials, thereby reducing the development time involved in the parametrization of reactive potentials. The methodology for a reaction can be extended from one environment (bulk water48) to the other (aqueous HCl48 solution and water confined in CNTs49) as long as they have the same TS structures. The information regarding the TS is embedded within the triggers. Therefore, the algorithm has the benefit of identifying its local environment to capture the effect at least qualitatively. The RMD algorithm models the structural diffusion of a proton in an aqueous system as a chemical reaction nH2 O

H3 Oþ þ H2 O rsf H2 O þ H3 Oþ

ð1Þ

where n is the additional number of water molecules present during the reaction. The first step of the algorithm is to check whether the reactants (H3O+ and H2O) satisfy a set of geometric and energetic triggers. The energetic trigger (Ea,f) is identified 18836

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C

ARTICLE

from the Ea required to overcome the energy barrier for the proton to be moved from one water molecule to the other. The geometric trigger checks whether the reactants are in a favorable configuration for the reaction to take place and is based on the reactant, product, and TS structures. The use of nonreactive potentials prevents the molecules from reaching the TS. These triggers are parametrized to reproduce the macroscopic property of k and Ea. Hence, the algorithm was first applied to proton transport in bulk water48 which has the most reliable k and Ea from both experiment52
54 and simulation.32,55 Instantaneous reaction takes place in the second step of the algorithm where the proton is explicitly transferred from one water molecule to the other by converting the covalent bond between the H* and O of reactant H3O+ into a hydrogen bond and vice versa with O of reactant H2O. The instantaneous reaction involves a change in identity of the molecules (potentials and charges), movement of proton, and breaking of bonds that cannot be captured by the classical MD simulation. Therefore, a reaction is always accompanied by the change in potential energy of the system and disturbance to the complex 3-D hydrogen bonding network. Local equilibration is the final step that relaxes selected molecules to satisfy the target ΔHr and to re-establish a reasonable hydrogen bonding network. There is a need to re-establish the hydrogen bonding network because the reaction results in products that are not in equilibrated configuration. For example, the O
O radial distribution function (RDF) for H2O
H2O is different than that for H3O+
H2O, and when we change the identities of the reactants during the reaction the RDFs are affected. Therefore, one of the purposes of the local equilibration is to establish a reasonable configuration for the products. The ΔHr for reaction depicted by eq 1 is zero since the products and reactants are identical. Since we know ΔHr for the structural diffusion of a proton and local structure from the RDF, a local equilibration scheme can be devised. For this purpose, we have identified an objective function, Fobj, that has an energetic and geometric component as given below Fobj ¼ w1 ΔgðrÞRMS þ w2 ΔU RMS

ð2Þ

The energetic component, ΔU , is the root-mean-square (rms) error of the energy difference of the system before (Ubefore) and after (Uafter) the reaction (since ΔHr = 0) sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   Uafter
Ubefore 2 ΔU RMS ¼ ð3Þ Ubefore rms

The structural component, Δg(r)rms, is the rms error of the current distance (rij) between two atoms, i and j, and the ) from the RDF and can be described as expected distance (rtarget ij vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u u 1 j ¼ Npairs rij
rijtarget 2 ΔgðrÞRMS ¼ t ð4Þ target Npairs j ¼ 1 rij

∑

where Npairs refers to the number of pairs of atoms considered. The satisfaction of energy or structural component is distributed by the weighting factors, w1 and w2. The multivariate nonlinear function described by eq 4 is optimized using the Polak
Ribiere conjugate gradient method.56 A selected set of molecules are relaxed based on Fobj, while the rest of the system is held rigid. The list of molecules relaxed is discussed later. During local equilibration, the MD time clock is stopped, and once completed the MD time step resumes.

There are a number of parameters and variables that need to be identified in the algorithm. The geometric triggers involved in the first step of the algorithm and the following variables in Fobj must be defined: (1) selection of molecules to be relaxed, (2) identification of the structural points (Npairs), and (3) specification of w1 and w2. Therefore, it is possible to approach the execution of the RMD algorithm in numerous ways. One way to confirm the selection of a scheme for the structural diffusion of a proton is to compare the water diffusivities in a reactive and nonreactive system (where protons are not allowed to hop) or pure water (TIP3P model). Experimentally, the selfdiffusion coefficient of water is not affected by protons at dilute concentrations.57 Hence the addition of the RMD algorithm should not alter the water diffusivities in a reactive system from that of a nonreactive system, at least under dilute conditions. The structural diffusion of a proton in bulk water, aqueous HCl solution, and CNTs was modeled with six geometric triggers and a local equilibration scheme that relaxed the four hydrating molecules around the reactants and the H+ transferred with Npairs = 11, w1 = 1, and w2 = 10. A high weight to the energy component was chosen to compensate for the low contribution of the energy component to Fobj compared to the structural component. The specific details of the triggers and local equilibration are already referred to in our previous studies conducted.48,49 The six geometric triggers are chosen based on the equilibrium structures of H3O+, H2O,58 and the two protonated complexes—Zundel59 and Eigen60 cations that play a major role in charge diffusion.12,15 The six geometric triggers are as follows: (1) rOO,Zundel, checks the O(H3O+)
O(H2O) separation since the reactants must be close enough to form a Zundel ion; (2) rOH,Equilb, checks whether the covalent bond distance between the O(H3O+) and the proton to be transferred, H*, exceeds the equilibrium bond distance; (3) θOHO, checks whether the angle formed by the O(H3O+), H*, and O(H2O) is nearly linear, similar to that observed in the Zundel ion; (4) θHOH, checks whether the angles formed by the H*, the O of the H2O, and each of the two H of the H2O is near the equilibrium H
O
H bond angles of the H3O+ ion since the lone pair of electrons in the water should point toward the H*; (5) rOO,Eigen, checks whether each of the two nonreactive H atoms of the reactant H3O+ ion forms a hydrogen bond with a H2O molecule since the reactant H3O+ ion forms an Eigen ion; and (6) rOO,hyd, checks whether each of the two H atoms of the reactant H2O molecule forms a hydrogen bond with a H2O molecule since an Eigen cation is formed around the product H3O+ ion. These triggers are parametrized to reproduce the macroscopic property of k and Ea. The restructuring of hydrogen bonds in the solvation shell around the excess charge with the transfer of H+ has been observed by ab initio studies.15,61 The RMD algorithm aims at capturing the proton transport description to an extent though not as detailed like a quantum mechanical description. Therefore, we decided to relax the four hydrating molecules during the local equilibration. Other parameters in Fobj were chosen by trial and error with preliminary short runs until a reasonable water diffusivity was obtained. The sensitivity and adaptability of the algorithm to different environments has been previously tested and proved.48,49 The pH dependence on the structural diffusion of a proton was investigated by implementing the algorithm to aqueous HCl systems of different concentrations (0.22
0.83 M).48 The algorithm qualitatively captured the decrease in the total diffusion of charge with an increase in acidity of the HCl solution. It also 18837

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C

ARTICLE

Table 1. Description of the Methods Analyzed with the RMD Algorithma method

N

Npairs

w1

w2

1

39

12

1

10

2

72

26

1

10

a

N represents the degree of freedom of the Fobj; Npairs is the number of the structural points used to characterize the final configuration of products and hydrating molecules; w1 and w2 are parameters in eq 2.

Figure 1. Typical structure of the reactants surrounded by six hydrating molecules in the bulk water system. Atom and molecule labels serve for identification purposes in Table 2. O of H3O+, green; O of H2O, red; H, white.

demonstrated that the primary reason for this reduction was due to a decrease in structural diffusion. Proton diffusion in water confined in CNTs of varying radii (5.42
10.85 Å) was studied to analyze the effect of confinement on structural diffusion.49 The axial structural component of the total charge diffusion was found to decrease with the increase in confinement. Thus, with the aid of the geometric and energetic triggers, the model is able to capture the essential features of the structural diffusion of proton irrespective of its surroundings and physical conditions that can affect the energy and solvation structure of the proton. Having analyzed the individual effects of confinement and acidity, the RMD algorithm is applied in this work to study proton transport in PEMs, where both nanoscale confinement and acidity are present. One can categorize the structural diffusion of a proton in a PEM into the following three reactions based on the reactants, hydration level, and hydrating molecules nH2 O

H3 Oþ þ H2 O rsf H2 O þ H3 Oþ H3 O þ þ H 2 O

mH2 O=ðn
mÞSO3


rsf

nH2 O

H2 O þ H3 Oþ

SO3 H þ H2 O rsf SO3
þ H3 Oþ

ð5aÞ ð5bÞ ð5cÞ

The reaction given by eq 5a is similar to that of bulk water (eq 1) where it is hydrated only by the surrounding water molecules. Equations 5b and 5c may occur along the interface of hydrophobic and hydrophilic regions where the nonreactive H atoms of the reactants can be hydrated by the O(SO3
) in addition to water molecules. Only eq 5a will have TS structures of Zundel and Eigen cations similar to that found in bulk water, whereas eq 5b will have TS structures of “Zundel-like” and “Eigen-like” complexes, in which one or more of the hydrating water molecules are replaced by a sulfonate group. Since the algorithm has been parametrized for the structural diffusion of proton in bulk water only the reaction having the same transition state is modeled in this work. A Separate set of triggers and parameter values are required for modeling the other two reactions. The last two reactions are likely more prevalent at low water contents. Therefore, in this work, we limit the simulation to water contents λ g 6H2O/SO3H. The chosen triggers and local equilibration method in the RMD algorithm seemed sufficient in a simple system of bulk water at infinitely dilute concentration and proved efficient in

aqueous HCl and water confined in CNTs. However, when moved to complex system of PEM (Nafion) which has both confinement and acidity the current choice of the scheme exhibited some unexpected behaviors that are discussed below. Specifically, the water diffusivities in the reactive Nafion system are not a precise match to that of the nonreactive system. Improving the local equilibration is the next logical step for investigating the sensitivity of water diffusivity to the presence of structural diffusion in PEMs. Therefore, we have introduced a second scheme of execution of the RMD algorithm in this study. In this new scheme, Npairs has been increased to include more details of the solvation shell structure, and a larger number of molecules are relaxed to equilibrate the structure better. The previous scheme applied in bulk water and other systems is referred to henceforth as Method 1, while the new scheme is referred to as Method 2. The high concentration of H3O+ ions in confined aqueous regions of the PEM sometime places the product H3O+ ion next to another H3O+ ion after the reaction. Therefore, to have an environment in PEM similar to that of bulk water, hydration by six molecules is checked instead of four hydrating molecules from the initial assessment. This addition results in the reacting hydronium ion and water molecule having to participate in four hydrogen bonds. This can also be conceived as inclusion of two more geometric triggers (rOO,h3o4 checks whether the oxygen of the reactant H3O+ ion is close enough to the oxygen of the fourth water molecule and rOH,h2o4 checks whether the oxygen of the reactant H2O molecule is close enough to the hydrogen of a water molecule other than two H2O molecules hydrating the H atoms of reactant H2O molecule). A typical structure of the reactants surrounded by six hydrating molecules in bulk water is shown in Figure 1. Each molecule and atom in Figure 1 is labeled for future reference. Table 1 lists the method number, the total degree of freedom, N (number of atoms relaxed  dimensionality), Npairs, w1, and w2. In Method 2, the six hydrating molecules and the two reactants except O1 are relaxed based on Fobj, which causes N in the new scheme to increase from 39 to 72. A list of the 26 Npairs used in Method 2 is given in Table 2. The rtarget values are obtained from ij the RDF between the atoms i and j from classical MD simulation in the previous paper.48 Method 2 concentrates not only on the first solvation shell but also on the second solvation shell. We have included a new scheme with a major modification in the local equilibration method to point out that we understand the reason behind the quantitative disagreement between the simulation and experimental/estimated results and to prove the existence of ways to improve the results discussed in Section III. From the general point of view of RMD as a coarse-graining procedure, there is always a balance between the gain in computational efficiency and the loss of model fidelity through the reduction of degrees of freedom. It is frequently difficult a priori to judge the best choice in resolution. From this 18838

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C

ARTICLE

Table 2. Values of rtarget Used in Equation 4 Obtained from ij RDF between Atoms i and ja

a

Npairs

i
j

rtarget (Å) ij

Npairs

i
j

rtarget (Å) ij

1

O3
O2

2.80

14

O6
O2

3.75

2 3

O4
O2 O5
O1

2.80 2.55

15 16

H4
O2 H3
O2

2.85 2.85

4

O6
O1

2.55

17

H4
O5

2.85

5

O1
O2

2.55

18

H5
O5

2.85

6

O3
H1

1.85

19

H3
O6

2.85

7

O4
H2

1.85

20

H5
O6

2.85

8

O5
H3

1.60

21

O1
O8

3.30

9

O6
H4

1.60

22

O2
O7

2.80

10 11

O3
O1 O4
O1

4.90 4.90

23 24

H5
O3 H5
O4

4.10 4.10

12

H5
O2

1.60

25

H5
O7

4.10

13

O5
O2

3.75

26

O2
H6

1.90

All the distances are measured in Å. The labels appear in Figure 1.

perspective, Methods 1 and 2 vary the degree of detail in the RMD model, with Method 2 containing greater detail. We acknowledge that Method 2 is not the only way to improve the results. Since a detailed understanding of the effect of parameters of the RMD algorithm onto the system property is required to proceed in the right direction, we have done a comparison of the results obtained from Methods 1 and 2. Generally, implementation of the new local equilibration methods involves the following steps: (1) the parameters in Fobj are defined; (2) the method is applied to bulk water, and the triggers are reparameterized to match the experimental k and Ea,f; (3) the water diffusivity in bulk water is confirmed not to be affected by the inclusion of the algorithm; and (4) the method is applied to PEM. This work is based on the basic assumption that the structural diffusion of a proton along the center of aqueous channels where the water behaves more bulklike62 is similar to the mechanism found in bulk water. We stress that this RMD approach is a coarse-grained approach, with the intent to capture observed transport phenomena, without resorting to the quantum mechanical solution of electron distribution. From this perspective, the reaction describing proton transport in eq 1 is equally valid in bulk solution and in a PEM. Since the triggers are parametrized to match k and Ea for bulk water and they were able to identify the physical conditions like temperature,48 confinement,49 and presence of an acidic group,48 we believe at least to an extent the algorithm can capture the environment in a hydrated membrane. Partial validation of this RMD approach is available. This RMD approach has been used to model proton transport in highly concentrated bulk hydrochloric acid solutions and in carbon nanotubes. Based in part on these simulations, an analytical model for charge transport in Nafion membranes was developed that generates quantitative agreement in the reduced self-diffusivity of charge compared to experimental measurements.47 Therefore, there is some demonstration that this RMD approach is applicable to systems other than just that to which it was applied, namely, bulk water. B. Computational Details. NVT simulations are performed with Nafion at 300 K at four hydration levels (λ = H2O/SO3H) of 6, 9, 15, and 22. In contrast to previous work, a lower hydration

level than λ = 6 is not chosen in this study since we are interested only in the structural diffusion of a proton following the mechanism similar to that observed in the bulk aqueous system, and finding protons following that mechanism might be difficult at λ < 6. In the simulations, we have used the same interaction potentials for the polymer electrolyte, water, and hydronium molecules as used earlier in the simulations of the bulk hydrated membrane.8 Even the initial configurations for the current simulations are obtained from the production runs of the previous work. The equilibration details of the structure are not discussed here. The hydrated membrane simulations from our previous work8 are referred to henceforth as nonreactive systems since they were modeled using classical MD simulation and proton transport involved only the vehicular diffusion. A Nafion polymer with an equivalent weight (EW) of 1144 each consisting of 15 monomers is considered for the study. All the systems at λ = 6, 9, 15, and 22 have 40 Nafion ionomers, the corresponding 600 H3O+ ions, and 3000, 4800, 8400, and 12 600 H2O molecules, respectively. The densities are obtained from an empirical function derived from experimental results.63 The justification of this model and the successful comparison with similar simulations performed have been discussed previously.3,8 The CFX in the backbone is represented as an united atom with a single Lennard-Jones (LJ) parameter,11,64 and the backbone does not carry any charge. The side chain of Nafion is modeled with parameters defined by Vishnyakov and Neimark.65,66 The detailed potential model of Nafion and sources for the force constants are reported in our earlier work.8 Water is modeled using the TIP3P model58 with a flexible OH bond,67 while the H3O+ ion uses the same model with modified charges on O and H which is similar to that used by Urata et al.2 The bonded interactions include stretching, bending, and bond torsion. Coulombic and LJ interactions contribute to nonbonded interactions. The two time scale r-RESPA68 approach with a time step of 2.0 fs for the large time step and 0.2 fs for the intramolecular interactions is used for the integration of the equations of motion. The temperature is maintained by the Nose
Hoover thermostat69,70 with a slightly higher frequency of 0.01 fs
1 to remove the excess heat that was not removed by the local equilibration after the reaction. LJ potentials are cut off at 10 Å. The site
site reaction field method is incorporated to account for electrostatic interactions,71 and long-range LJ interactions are taken into account. The parameters for Method 1 of the RMD algorithm are the same as that of bulk water, aqueous HCl, and water confined in CNTs with a minor modification in the trigger values. We have the same six geometric triggers, one energetic trigger, and objective function for the local equilibration. However, we have slightly modified the definition of a reaction in this work. Earlier, a reaction was taken into account for k calculation when the H+ rattling between two molecules ended up with either of the parent molecule or another H2O molecule at the end of rattling. In the current work, we define a reaction taking place only when the proton ends up with a molecule other than the parent molecule at the end of rattling. In short, only if a H+ undergoes an odd number of rattles it is considered as a reaction. Due to this change, triggers were reparameterized (rOO,Zundel, 2.35
2.75 Å; rOH,Equilb, g0.965 Å; θOHO, g146
; θHOH, 65
143
; rOO,Eigen, 2.35
2.93 Å; rOO,hyd, 2.50
3.27 Å; Ea,f, g
23.719 kcal/mol) to match the experimental rate constant. At the end, we got a better 18839

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C

ARTICLE

Figure 2. Fitting of the rate constant to experimental value (ref 53) for the new set of triggers in Method 1 and Method 2. Occurrence of the reaction is taken into account only when the proton ends in a water molecule other than the parent water molecule.

agreement of the total charge diffusivity with the experimental value because of the better contribution of the structural component compared to the old definition of the reaction. Method 2 of the RMD algorithm also uses the new definition of reaction. The energetic trigger and eight geometric triggers were parametrized to match the experimental k (rOO,Zundel, 2.35
2.74 Å; rOH,Equilb, g0.965 Å; θOHO, g146
; θHOH, 65
143
; rOO,Eigen, 2.35
2.93 Å; rOO,hyd, 2.50
3.29 Å; rOO,h3o4, e3.65 Å; rOH,h2o4, e2.6 Å; Ea,f, g
23.417 kcal/mol). For the parametrization process, bulk water simulations at infinite dilute concentration (650 H2O molecules containing one H3O+ ion) were run. The potential models and simulation methodology are the same as in the previous work.48 The properties were calculated from 96 independent runs with 1 ns long production runs. RMD simulations of Nafion were run with both Method 1 and Method 2 schemes. The total length of the simulations at λ = 6, 9, 15, and 22 are 0.86, 0.46, 0.40, and 0.27 ns for Method 1 and 1, 1, 1, and 0.5 ns for Method 2. It was sufficient to run these simulations for a relatively short time since we are interested only in the dynamics of H2O and the H3O+ ion. Moreover, a newly developed confined random walk (CRW) diffusion model72 can be used to calculate the diffusivity in a long time limit from short MD or RMD runs.

III. RESULTS AND DISCUSSION This section mainly concentrates on transport properties of water and hydronium ions in Nafion since the structural properties like radial distribution function (RDF), density profile, hydration histograms, and cluster size distributions have been already discussed in detail by the earlier work of classical MD simulations.3,8 The results from RMD bulk water simulations using Method 1 and 2 schemes are also discussed in this section to validate the procedure before application in the PEM. The results of Nafion RMD simulation obtained from both methods are compared, and the advantage and disadvantage of each method is pointed out. A. Bulk Water Simulation. The parametrization process of the triggers is discussed in great detail in the earlier study,48 and

Figure 3. Charge diffusivity and its two components (vehicular and structural) as a function of temperature compared to the reference values. Solid symbols represent reference values, and empty symbols correspond to simulation values. Details of the reference values can be found in the text. (a) Method 1. (b) Method 2.

hence they are not reported again in this work. However, the fitted k and measured transport properties of H2O and H3O+ ions are discussed here to show that both the schemes are parametrized and can be applied to define structural diffusion of the proton in Nafion following the mechanism similar to that of bulk water. The geometric and energetic triggers were parametrized so that the rate constant from simulation matched the experimental value. Figure 2 compares the fitted k for Methods 1 and 2 with the experimental measurements.53 We find that the fitted values are in good agreement with the experimental values with an average relative error of 5%. Total charge diffusivity, structural component, and vehicular component along with their corresponding reference values are plotted for Method 1 and Method 2 in Figure 3(a) and (b), respectively. The total charge diffusivity is compared to that of experimental value;14,73 structural diffusivity is compared to the value estimated from Einstein relation using the average lifetime of a proton;12,74 and vehicular diffusivity is compared to the difference between the experimental value and structural component (under the assumption they are uncorrelated). As can be seen from Figure 3(a), the current definition of reaction (where only the protons that end up in a water molecule other than the parent water molecule are considered to have undergone a reaction) and new trigger values 18840

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C

Figure 4. Water diffusivities in the bulk water system at infinite dilute concentration as a function of temperature using the local equilibration schemes listed in Table 1. Pure water simulation using the TIP3P model is used as reference [ref 48].

of Method 1 have reproduced the transport property of the charge and its two components with better match to the reference values than the old definition of the reaction.48 Using Method 2 for the execution of the RMD algorithm has slightly increased the vehicular and structural component of the charge diffusion which has resulted in an overall increase in the total charge diffusivity when compared to the experimental value. Figure 4 compares the water diffusivity obtained by both methods to the experimental bulk water diffusivty,75,76 pure bulk water diffusivity using the TIP3P model,48 and water diffusion coefficient of a nonreactive system.48 As discussed in Section II.A, the method to test our choice of local equilibration and triggers is to confirm the water diffusivity is not affected by the inclusion of the RMD algorithm. The TIP3P water model is known to give diffusion coefficients higher than experimental value and is used as a benchmark to compare the diffusivities from nonreactive and reactive systems. Method 2 gives much higher diffusivities than Method 1, implying that the local relaxation in Method 2 is not as efficient as Method 1 in bulk water. However, as can be seen from later discussions, Method 2 gives lower H2O diffusivity in Nafion than Method 1, and that is the reason we have pursued Method 2 for a complete analysis. B. Reaction Statistics. The reactivity of the structural diffusion can also be measured in terms of lifetime or mean residence time of proton, τP. The τP for bulk water at 300 K is around 1.7 ps.53 We have measured the average τP in Nafion from the direct count of reactions divided by the simulation time, and they are plotted in Figure 5. In both systems, the τP is found to decrease with increase in the water content due to the increase in number of reactions with the increase in λ. The number of H3O+ ions satisfying the given set of triggers increases with hydration level. Even at the highest water content τP is nearly 10 times higher than that of bulk water. Though the triggers in both the methods were parametrized to fit the experimental rate constant (i.e., τP is the same for both methods in bulk water), τP from Method 1 is consistently higher than Method 2 in Nafion because the latter is more restrictive in the aqueous nanochannels of PEM with the two additional triggers. Devanathan et al. applied Q-HOP MD to study proton transport in a system of four Nafion polymers each

ARTICLE

Figure 5. Average τP of the protons in Nafion for Method 1 and 2 reactive systems measured from the total number of reactions taken place as a function of water content.

Figure 6. Histogram of the fraction satisfying the geometric and energetic triggers in Method 1 at various hydration levels in Nafion. Each fraction is calculated based on the satisfaction of the previous trigger. Description of the triggers can be found in Section II.A.

with 10 side chains. They allowed structural diffusion of a single proton, while the rest of them were allowed to diffuse only through vehicular mechanism. At λ = 5, 10, and 15 they measured the mean τP of the proton to be 250, 6.6, and 2.9 ps, respectively, from the time correlation of population operator. The reason for τP in the RMD system being consistently higher than Q-HOP MD systems is likely due to the differences in the coarse-graining of the proton-shuttling reaction. The reactivity in an RMD algorithm is set by the triggers. The concept of having geometric and energetic triggers satisfied for the occurrence of a reaction in the RMD algorithm helps us to identify the triggers that promote the decrease in τP with λ. The fraction satisfying each trigger discussed in Section II.A is shown in Figure 6 for all the water contents of the Method 1 RMD system of PEM. The two methods are not compared because the number of triggers involved varies. Each fraction is calculated based on the satisfaction of earlier triggers, and that is why only six triggers are presented in the plot. The hierarchy of the trigger is based on the manner they were tested in the code and was 18841

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C

ARTICLE

chosen for computational efficiency. From Figure 6, we observe that the energetic trigger provides the most severe limitation in preventing the H3O+ ions from reacting, and hydration of the reactant water molecule is also affected by the change in λ. In other words, the energetic trigger is the most sensitive to change in water content, reflecting changes in both acidity and confinement. An analogous result was observed for proton transport in HCl solution48 and water confined in CNTs.49 C. Water Diffusivity. Diffusion of protons and water in the aqueous domains of a PEM undergoes nonlinear time-dependent confined diffusion initially, and then they follow Einstein’s linear relation in the infinite time limit. The water diffusion coefficient in a nonreactive system can be unambiguously calculated using the Einstein expression D ¼ lim

τf∞

ƽrðt þ τÞ
rðtÞ2 æ 2dτ

ð6Þ

where d is the dimensionality of the system; τ is the observation time; and Æ[r(t + τ)
r(t)]2æ is the mean square displacement (MSD). In a reactive system, since the molecules undergo change in identity whenever the reaction occurs, the diffusivity is measured in the time segments in which they did not react. The production run of the nonreactive system by Liu et al.8 was run 2 ns long, which is sufficient for the structural data analysis but was found insufficient for the diffusion coefficient calculation. The exponent, m, in MSD  τm, falls in the subdiffusive range of 0.69
0.85 for λ = 3
22. Recently, we have developed and implemented a CRW simulation approach72 that can extend the short-time MSD data for diffusion in the presence of nanoscale confinement obtained from MD simulation to statistically reliable MSDs for a longer time limit. Therefore, instead of running the simulation longer than 2 ns for the sole purpose of MSD data collection, a simple confined random walk (CRW) simulation can be integrated with classical MD simulation of a short time scale to get diffusivities in the long time limit where m is 1
MD/CRW algorithm. The computational expense added to the classical MD simulation by the RMD algorithm can also be overcome by the CRW approach. Currently, we have applied the CRW simulation to obtain the diffusivities in the long time limit from the MD or RMD simulation MSD. The MD/CRW algorithm applied to the nonreactive system did not affect the diffusivities significantly since m was in the range of 0.69
0.85.72 Though the reactive system RMD simulations were not run as long as the nonreactive system, they had m in the range of 0.60
0.95 due to the increased diffusion rate. Water diffusivity in the nonreactive system8 and reactive systems using Method 1 and Method 2 schemes and experimentally measured17 values is shown in Figure 7 as a function of water content. We find the water diffusivities in all four systems to increase with the hydration level. The TIP3P water model58 inherently has a higher diffusivity value than the experimentally measured bulk water value.77,78 Therefore, water diffusivity in the nonreactive system is higher than the experimental value at all λ but qualitatively follows the same trend. Implementation of the RMD algorithm is not expected to affect the water diffusivity in a reactive system from that of a nonreactive system. However, in the reactive system using Method 1 scheme, the water diffusivity is found to be nearly 4.15 times higher than the nonreactive system at the highest water content. This increase can be attributed to the improper local equilibration in Method 1 which seems to be insufficient in a complex system like Nafion

Figure 7. Comparison of diffusion coefficients of water from experimental (ref 17), nonreactive (ref 8), and reactive systems at different hydration levels. The water diffusivity in the nonreactive system is expected to be equal to that of the reactive system.

unlike the simple system of bulk water for which triggers were parametrized. The more detailed local equilibration in Method 2 has lowered the water diffusivity in a reactive system compared to Method 1 but at the cost of increasing water diffusivity in bulk water as shown in Figure 4. Therefore, it is possible to lower the water diffusivity in a reactive system to match the results from a nonreactive system by improving the local equilibration scheme, but a compromise has to be made between the bulk water and PEM simulations. The increase in energy of the system due to the instantaneous reaction was decreased by 67% on average by both the local equilibration methods, but since Method 2 included a higher number of molecules in relaxation it might have had better influence on the hydrogen bonding network. Further improvement can be made in the local equilibration with this basic understanding of the scheme. D. Charge Diffusion. One of the features of the RMD algorithm is that one can be specific about the chemical reactions that take place in a system. Only the structural diffusion of a proton having the same TS as that of bulk water is included in this work. So the total charge diffusion measured is a combination of vehicular component and structural diffusion contributed only by the reaction given by eq 5a. The total charge diffusion can be measured using the Einstein relation in eq 6 by following the trajectory of the center of mass position of the H3O+ ion. To decouple the individual contributions of the vehicular and structural component to the total charge diffusion, the total displacement of the H3O+ ion in each step can be written as35 ΔB r tot ¼ Δ B r veh þ Δ B r struct

ð7Þ

where Δr Bveh is the displacement due to vehicular diffusion or in other words corresponds to the displacement due to classical MD simulation, and Δr Bstruct is the structural component displacement. Whenever there is no reaction in a MD step, Δr Bveh becomes the sole contributor to Δr Btot. Bstruct = 0 and Δr is approxiWhen a reaction takes place in a MD step, Δr Bstruct mately equal to 2.55
2.65 Å which is the distance between + the reacting O(H3O )
O(H2O), and Δr Btot is the sum of both Δr Bstruct. Therefore, the trajectory of Br veh is Bveh and Δr 18842

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C

ARTICLE

Figure 8. Charge diffusivity from experimental conductivity (refs 17 and 79), classical H3O+ ion diffusivity in a nonreactive system (ref 8), total charge diffusivity, Dtot, in a reactive system, and its (i) vehicular component, Dveh, and (ii) structural component, Dstruct, as a function of hydration level. (a) Method 1 execution of the RMD algorithm. (b) Method 2 execution of the RMD algorithm.

contributions. The discussion on Dcorr is provided later in the paper. Figures 8(a) and (b) show the experimental charge diffusivity,17,79 H3O+ ion diffusivity from the nonreactive system of Liu et al.,8 and Dtot from the reactive system using Method 1 and Method 2 RMD schemes, respectively, along with Dveh and Dstruct as a function of water content. Experimentally the charge diffusivity is found to increase with water content, and the same trend is followed by all the diffusivities shown in the plots. Addition of the RMD algorithm was found to have not much altered the vehicular diffusion of H3O+ ions in bulk water as shown in Figures 3(a) and (b). The vehicular component of the charge diffusion was almost equal to the H3O+ ion diffusivities from the nonreactive system even for water confined in CNT systems.49 However, in the application of RMD algorithm using Methods 1 and 2 to PEMs, Dveh are higher than the diffusivities of H3O+ ions in the nonreactive system for the same reason of finding an increase in the water diffusivities of a reactive system from that of a nonreactive system which is due the insufficient relaxation of the molecules during local equilibration. Corresponding to the increase in a number of reactions with the hydration level, Dstruct also increases as seen from Figures 8(a) and (b). In both Methods 1 and 2, Dtot is lower than the sum of Dveh and Dstruct at all λ implying they are negatively correlated (Dcorr is negative), and at the lowest water content Dtot is even lower than Dveh. In Figure 8(a), Dtot is higher than the experimentally observed values because the Dveh is too high. The change in the local equilibration method in Method 2 has not only affected the Dveh and water diffusivity but also lowered Dstruct. Using Method 2 implementation of the RMD algorithm has lowered Dtot and Dveh when compared to Method 1. In Figure 8(b), Dtot is found to be close to Dveh at all water contents, and at higher λ, Dtot is also closer to experimental values. At lower water contents, as shown in Figure 8(b), Dtot is less than the experimental value, and this might be due to modeling of only the reactions that are similar to those observed in bulk water. E. Correlation. The correlation between the structural and vehicular components of the total charge diffusion can be obtained from the Dcorr term in eq 9c. Dcorr is calculated from the relation Dcorr ¼ Dtot
ðDveh þ Dstruct Þ

continuous, andBr struct is discontinuous. By substituting eq 7 in eq 6, total charge diffusivity, Dtot, can be written as Dtot ¼ lim

τf∞

ÆΔ B r 2veh æ þ ÆΔ B r 2struct æ þ 2ÆΔ B r veh Δ B r struct æ 2dτ

ð8Þ The above definition can be decomposed as Dveh  lim

τf∞

Dstruct  lim

ÆΔ B r 2veh æ 2dτ

τf∞

Dcorr  lim

τf∞

ÆΔ B r 2struct æ 2dτ

2ÆΔ B r veh Δ B r struct æ 2dτ

ð9aÞ

ð9bÞ ð9cÞ

where Dveh and Dstruct are vehicular and structural diffusion coefficients, respectively. Dcorr is the correlation term that denotes the coupling between the structural and vehicular

ð10Þ

If Dcorr = 0, then both the components of the total charge diffusion are uncorrelated. Proton transport in the previously studied system like bulk water,48 aqueous HCl, and water confined in CNTs49 exhibited uncorrelation between the two components of Dtot. MS-EVB2 has also observed negligible negative correlation of the two components in bulk water.35 Figure 9 plots the Dcorr/Dtot as a function of hydration level for the two schemes of execution of algorithm. The plot shows a decrease in Dcorr/Dtot with increase in λ for Method 1, implying the two components of the charge diffusion are anticorrelated at low water contents, and at high water contents, as they approach bulk water behavior, their anticorrelation diminishes. In Method 2, Dcorr/Dtot is almost constant for λ = 6
22 with both the components of charge diffusion anticorrelated. This might suggest that with Method 2 of local equilibration the studied water content range is not sufficient to observe a decrease in the anticorrelation. Petersen et al.35 observed a strong anticorrelation between the two components in Nafion at λ = 15 and quantified the time scale of these correlated motions as 425 ps 18843

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C

ARTICLE

Figure 9. Dcorr/Dtot as a function of water content for reactive systems using Method 1 and 2 RMD schemes. If Dcorr = 0, the structural and vehicular components are uncorrelated. Negative values indicate both components are anticorrelated.

using the distinct portion of the van Hove correlation function. An analytical model47 developed to model the charge diffusion through the aqueous domain of Nafion quantitatively reproduced the experimental results. The model had the three factors
acidity, confinement, and connectivity, incorporated to characterize charge diffusion along with other assumptions like the two components were uncorrelated. Inclusion of the correlation factor between the two components of charge diffusion could have lowered charge diffusion coefficients obtained from the model compared to the experimental value. The relative magnitude of the anticorrelation remains significant even at low values of λ because not only does the structural component decrease with a decrease in hydration but also do the vehicular component and the total value. Furthermore, although the structural component has a much smaller magnitude than the vehicular component, from a statistical point of view, it is not surprising that the correlation term has a magnitude that falls between the two components. Understanding the anticorrelation between the two components cannot be done with direct decoupling. Hence, we formulated two functions to study the correlation between the O(H3O+) and an almost stationary point in the system like S(SO3
) for the structural (fstruct) and vehicular components (fveh). After every reaction, we identify the S*(SO3
) that is closest to the product O(H3O+). fstruct is calculated as fstruct ¼ ð B r OðH3 Oþ Þ
B r SðSO3
Þ Þ
ð B r OðH2 OÞ
B r SðSO3
Þ Þ pdt pdt

Figure 10. fstruct as a function of hydration level for the two reactive systems of Nafion. fstruct measures the difference in position of the product H3O+ ion and H2O molecule relative to the S of the closest sulfonic group. If fstruct < 0, then the H3O+ ion is moving toward SO3
, and if fstruct > 0, then the H3O+ ion is moving away from the SO3
group.

Figure 11. fveh as a function of hydration level. If fveh < 0, then the H3O+ ion is moving toward SO3
, and if fveh > 0, then the H3O+ ion is moving away from the SO3
group.

numerous reactions. Since the two components are anticorrelated if the structural diffusion tends to bring the H3O+ ion closer to the sulfonic group, then the vehicular diffusion has to move the H3O+ ion away from the sulfonic group. The above phenomenon can be tested by measuring fveh as a function of time, where

ð11Þ

fveh ðtÞ ¼ XðtÞ
Xð0Þ, ½XðtÞ ¼ j B r OðH3 Oþ Þ
B r SðSO3
Þ jt  pdt

where Br is the position vector of the corresponding molecule. When fstruct is less than zero, it means the reaction is bringing the H3O+ ion closer to the sulfonic group, and when fstruct is greater than zero, it implies the product H3O+ ion is moving away from the sulfonic group. fstruct averaged over all the reactions is plotted as a function of λ in Figure 10 for the two methods. We find the fstruct to approach zero from a negative value, implying the tendency of the reaction to bring the H3O+ ion closer to the hydrophobic/hydrophilic interface reduces or the cluster size is big enough to avoid the effect of sulfonate ions. It is to be kept in mind that the numbers in the plot represent the overall effect of

ð12Þ In other words, fveh tracks the trajectory of the H3O+ ion relative to S*(SO3
). If fveh(t) is negative, then the H3O+ ion is moving toward the sulfonic group, and if fveh(t) is positive, then the H3O+ ion is moving away from the sulfonic group. Figure 11 shows the fveh as a function of time at different λ for the reactive systems using Method 1 and 2 schemes. fveh is averaged over the trajectories of all H3O+ ions after each reaction and is calculated for a time period of 1000 fs or shorter if they undergo another reaction within the specified time frame. From the plot, we can 18844

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C see that as λ increases the magnitude of fveh also increases and is always positive. Combining the results shown in Figures 10 and 11, we can confirm that the structural diffusion brings the H3O+ ion closer to the sulfonic group, while vehicular diffusion moves it in the opposite direction. A similar conclusion was deduced by Petersen et al.35 from the distinct portion of the van Hove space
time correlation for the hydronium
sulfonate ion pair where they observed the fluctuating bond topology resetting the hydronium ion back to a position relative to the sulfonate ion, preventing it from diffusing away. Though the decrease in anticorrelation found in Method 1 or constant negative correlation in Method 2 cannot be described quantitatively by the above formulations, they at least explain the reason behind the negative coupling of Dveh and Dstruct.

IV. CONCLUSIONS Proton diffusion in Nafion occurs through a combination of both vehicular and structural diffusion mechanisms. In this work, the RMD algorithm is applied to model the structural diffusion of protons in Nafion (EW = 1144) that share the same transition state as in bulk water at λ = 6, 9, 15, and 22 at 300 K. The algorithm partitions the structural diffusion of protons into three steps: (i) trigger satisfaction, (ii) instantaneous reaction, and (iii) local equilibration. The algorithm was implemented in Nafion using two different approaches (Method 1 and Method 2). Method 1 increased the water diffusivity in a reactive Nafion system compared to a nonreactive system. Therefore, to lower the water diffusivity in a reactive Nafion system, Method 2 was introduced. Method 2 had a much more detailed local equilibration scheme and increased number of geometric triggers than that of Method 1. Both the methods had their triggers in step 1 of the algorithm parametrized to match the simulation rate constant in bulk water to that of experimental value. This work mainly focuses on dynamics of water, proton transport, and the correlated motion described by the structural and vehicular component of the charge diffusion. The lifetime of the proton in Nafion decreased with the increase in water content. The triggers that slowed the reaction at lower water contents were the energetic trigger and geometric trigger that checked the hydration of reactant water molecule by two water molecules. Water diffusivities from both methods increased with the hydration level but were higher than the values from the nonreactive system. However, the water diffusivities in Nafion from Method 2 were lower than Method 1 due to the improvement in local relaxation of the molecules in the last step of the RMD algorithm, but Method 2 caused an increase in bulk water diffusivity compared to Method 1, implying there is uncertainty in choosing the right method since a compromise has to be made between the results obtained from bulk water and Nafion. The contribution and correlation between the structural and vehicular component to the total charge diffusion were analyzed. The two components were found to be negatively correlated with their anticorrelation decreasing with the increase in λ in Method 1 and almost constant negative correlation with the change in λ in Method 2. Two simple functions were formulated using the position of S(SO3
) as a reference point to describe the anticorrelation. From the functions, it was proved that in both methods as the vehicular diffusion tried to move the H3O+ ion away from the sulfonate group the structural diffusion moved the H3O+ ion toward the sulfonate group. The total charge diffusivity obtained from both

ARTICLE

the methods was found to increase with the increase in water content. We understand that the dynamic properties of water and charge are strongly dependent on the RMD scheme of implementation, and these analyses would help us formulate a better scheme.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: dkeff[email protected].

’ ACKNOWLEDGMENT This research project is supported by the U.S. Department of Energy’s (DOE) Office of Basic Energy Sciences program (grant number DE-FG02-05ER15723). This work benefitted from the resources of the Center for Computational Sciences at Oak Ridge National Laboratory, which is supported by the Office of Science of the DOE under contract DE-AC05-00OR22725. ’ REFERENCES (1) Costamagna, P.; Srinivasan, S. J. Power Sources 2001, 102, 253. (2) Urata, S.; Irisawa, J.; Takada, A.; Shinoda, W.; Tsuzuki, S.; Mikami, M. J. Phys. Chem. B 2005, 109, 4269. (3) Cui, S. T.; Liu, J. W.; Esai Selvan, M.; Keffer, D. J.; Edwards, B. J.; Steele, W. V. J. Phys. Chem. B 2007, 111, 2208. (4) Fujimura, M.; Hashimoto, T.; Kawai, H. Macromolecules 1981, 14, 1309. (5) Devanathan, R.; Venkatnathan, A.; Dupuis, M. J. Phys. Chem. B 2007, 111, 8069. (6) Hristov, I. H.; Paddison, S. J.; Paul, R. J. Phys. Chem. B 2008, 112, 2937. (7) Kreuer, K. D. J. Membr. Sci. 2001, 185, 29. (8) Liu, J. W.; Suraweera, N.; Keffer, D. J.; Cui, S. T.; Paddison, S. J. J. Phys. Chem. C 2010, 114, 11279. (9) Jang, S. S.; Molinero, V.; Cagin, T.; Goddard, W. A. J. Phys. Chem. B 2004, 108, 3149. (10) Kreuer, K. D. ChemPhysChem 2002, 3, 771. (11) Li, Q. F.; He, R. H.; Jensen, J. O.; Bjerrum, N. J. Chem. Mater. 2003, 15, 4896. (12) Agmon, N. Chem. Phys. Lett. 1995, 244, 456. (13) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1995, 103, 150. (14) Robison, R. A.; Stokes, R. H. Electrolyte Solutions; Butterworths: London, 1959. (15) Markovitch, O.; Chen, H.; Izvekov, S.; Paesani, F.; Voth, G. A.; Agmon, N. J. Phys. Chem. B 2008, 112, 9456. (16) Zawodzinski, T. A.; Springer, T. E.; Davey, J.; Jestel, R.; Lopez, C.; Valerio, J.; Gottesfeld, S. J. Electrochem. Soc. 1993, 140, 1981. (17) Kreuer, K. D. Solid State Ionics 1997, 97, 1. (18) Petersen, M. K.; Wang, F.; Blake, N. P.; Metiu, H.; Voth, G. A. J. Phys. Chem. B 2005, 109, 3727. (19) Lee, C. H.; Park, H. B.; Lee, Y. M.; Lee, R. D. Ind. Eng. Chem. Res. 2005, 44, 7617. (20) Slade, S.; Campbell, S. A.; Ralph, T. R.; Walsh, F. C. J. Electrochem. Soc. 2002, 149, A1556. (21) Ivanovskaya, A.; Fan, J.; Wudl, F.; Stucky, G. D. J. Membr. Sci. 2009, 330, 326. (22) Zawodzinski, T. A.; Neeman, M.; Sillerud, L. O.; Gottesfeld, S. J. Phys. Chem. 1991, 95, 6040. (23) Ochi, S.; Kamishima, O.; Mizusaki, J.; Kawamura, J. Solid State Ionics 2009, 180, 580. (24) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471. (25) Billeter, S. R.; van Gunsteren, W. F. J. Phys. Chem. A 1998, 102, 4669. 18845

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846

The Journal of Physical Chemistry C (26) Intharathep, P.; Tongraar, A.; Sagarik, K. J. Comput. Chem. 2006, 27, 1723. (27) Warshel, A.; Weiss, R. M. J. Am. Chem. Soc. 1980, 102, 6218. (28) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361. (29) Wu, Y. J.; Chen, H. N.; Wang, F.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2008, 112, 467. (30) Chen, H. N.; Voth, G. A.; Agmon, N. J. Phys. Chem. B 2010, 114, 333. (31) Schmidt, R. G.; Brickmann, J. Ber. Bunsen-Ges. Phys. Chem. Chem. Phys. 1997, 101, 1816. (32) Lill, M. A.; Helms, V. J. Chem. Phys. 2001, 115, 7993. (33) Jang, S. S.; Lin, S. T.; Cagin, T.; Molinero, V.; Goddard, W. A. J. Phys. Chem. B 2005, 109, 10154. (34) Seeliger, D.; Hartnig, C.; Spohr, E. Electrochim. Acta 2005, 50, 4234. (35) Petersen, M. K.; Voth, G. A. J. Phys. Chem. B 2006, 110, 18594. (36) Jang, S. S.; Goddard, W. A. J. Phys. Chem. C 2007, 111, 2759. (37) Devanathan, R.; Venkatnathan, A.; Rousseau, R.; Dupuis, M.; Frigato, T.; Gu, W.; Helms, V. J. Phys. Chem. B 2010, 114, 13681. (38) Paddison, S. J. J. New Mater. Electrochem. Syst. 2001, 4, 197. (39) Roudgar, A.; Narasimachary, S. P.; Eikerling, M. Chem. Phys. Lett. 2008, 457, 337. (40) Choe, Y. K.; Tsuchida, E.; Ikeshoji, T.; Yamakawa, S.; Hyodo, S. Phys. Chem. Chem. Phys. 2009, 11, 3892. (41) Glezakou, V. A.; Dupuis, M.; Mundy, C. J. Phys. Chem. Chem. Phys. 2007, 9, 5752. (42) Eikerling, M.; Kornyshev, A. A.; Kuznetsov, A. M.; Ulstrup, J.; Walbran, S. J. Phys. Chem. B 2001, 105, 3646. (43) Paddison, S. J.; Paul, R.; Zawodzinski, T. A. J. Chem. Phys. 2001, 115, 7753. (44) Elliott, J. A.; Paddison, S. J. Phys. Chem. Chem. Phys. 2007, 9, 2602. (45) Venkatnathan, A.; Devanathan, R.; Dupuis, M. J. Phys. Chem. B 2007, 111, 7234. (46) Xu, H.; Kunz, H. R.; Bonville, L. J.; Fenton, J. M. J. Electrochem. Soc. 2007, 154, B271. (47) Esai Selvan, M.; Calvo-Munoz, E.; Keffer, D. J. J. Phys. Chem. B 2011, 115, 3052. (48) Esai Selvan, M.; Keffer, D. J.; Cui, S. T.; Paddison, S. J. J. Phys. Chem. C 2010, 114, 11965. (49) Esai Selvan, M.; Keffer, D. J.; Cui, S.; Paddison, S. J. Mol. Simul. 2010, 36, 568. (50) Perrin, J. C.; Lyonnard, S.; Volino, F. J. Phys. Chem. C 2007, 111, 3393. (51) Jiang, B. W.; Selvan, M. E.; Keffer, D. J.; Edwards, B. J. J. Phys. Chem. B 2009, 113, 13670. (52) Loewenstein, A.; Szoke, A. J. Am. Chem. Soc. 1962, 84, 1151. (53) Luz, Z.; Meiboom, S. J. Am. Chem. Soc. 1964, 86, 4768. (54) Glick, R. E.; Tewari, K. C. J. Chem. Phys. 1966, 44, 546. (55) Lapid, H.; Agmon, N.; Petersen, M. K.; Voth, G. A. J. Chem. Phys. 2005, 122. (56) Press, W. H.; Vetterling, W. T.; Teukolsky, S. A.; Flannery, B. P. Numerical Recipes in Fortran 77: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, 1992; Vol. 1. (57) McCall, D. W.; Douglass, D. C. J. Phys. Chem. 1965, 69, 2001. (58) Jorgensen, W. L.; J. Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (59) Zundel, G.; Metzger, H. Z. Phys. Chem. 1968, 58, 225. (60) Wicke, E.; Eigen, M.; Ackermann, T. Z. Phys. Chem. 1954, 1, 340. (61) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parinello, M. J. Phys. Chem. 1995, 99, 5749. (62) Lu, Z. J.; Polizos, G.; Macdonald, D. D.; Manias, E. J. Electrochem. Soc. 2008, 155, B163. (63) Weber, A. Z.; Newman, J. Chem. Rev. 2004, 104, 4679. (64) Cui, S. T.; Siepmann, J. I.; Cochran, H. D.; Cummings, P. T. Fluid Phase Equilib. 1998, 146, 51. (65) Vishnyakov, A.; Neimark, A. V. J. Phys. Chem. B 2001, 105, 9586.

ARTICLE

(66) Vishnyakov, A.; Neimark, A. V. J. Phys. Chem. B 2001, 105, 7830. (67) Neria, E.; Fischer, S.; Karplus, M. J. Chem. Phys. 1996, 105, 1902. (68) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990. (69) Nose, S. Mol. Phys. 1984, 52, 255. (70) Nose, S. J. Chem. Phys. 1984, 81, 511. (71) Hummer, G.; Soumpasis, D. M.; Neumann, M. Mol. Phys. 1992, 77, 769. (72) Calvo-Munoz, E.; Esai Selvan, M.; Xiong, R.; Ojha, M.; Keffer, D. J.; Nicoloson, D. M.; Egami, T. Phys. Rev. E 2010, 83, 011120. (73) Cornish, B. D.; Speedy, R. J. J. Phys. Chem. 1984, 88, 1888. (74) Meiboom, S. J. Chem. Phys. 1961, 34, 375. (75) Holz, M.; Heil, S. R.; Sacco, A. Phys. Chem. Chem. Phys. 2000, 2, 4740. (76) Easteal, A. J.; Price, W. E.; Woolf, L. A. J. Chem. Soc., Faraday Trans. I 1989, 85, 1091. (77) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2001, 114, 363. (78) Markovitch, O.; Agmon, N. J. Chem. Phys. 2008, 129, 084505. (79) Kreuer, K. D.; Paddison, S. J.; Spohr, E.; Schuster, M. Chem. Rev. 2004, 104, 4637.

18846

dx.doi.org/10.1021/jp203443c |J. Phys. Chem. C 2011, 115, 18835–18846