Atomistic and Coarse-Grained Molecular Dynamics Simulation of a

Article pubs.acs.org/Macromolecules

Atomistic and Coarse-Grained Molecular Dynamics Simulation of a Cross-Linked Sulfonated Poly(1,3-cyclohexadiene)-Based Proton Exchange Membrane Qifei Wang,† Nethika S. Suraweera,† David J. Keffer,*,† Suxiang Deng,‡ and Jimmy Mays‡ †

Department of Chemical and Biomolecular Engineering, University of Tennessee, Knoxville, Tennessee 37996, United States Department of Chemistry, University of Tennessee, Knoxville, Tennessee 37996, United States

‡

S Supporting Information *

ABSTRACT: Atomistic and coarse-grained (CG) molecular dynamics (MD) simulations were conducted for a cross-linked and sulfonated poly(1,3-cyclohexadiene) (xsPCHD) hydrated membrane with λ(H2O/HSO3) = 10 and 20. From the atomistic level simulation results, nonbonded pair correlation functions (PCFs) of water−water, water−H3O+ ion, H3O+−H3O+, polymer−water, and polymer−H3O+ ion pairs were obtained and studied. The water self-diffusivity and H3O+ vehicular self-diffusivity were also obtained. Membrane structure was further studied at CG level using a multiscale modeling procedure. Nonbonded PCFs of polymer−polymer pairs were obtained from atomistic simulation of hydrated membrane with λ = 10 and 20. Two sets of CG nonbonded potentials were then parametrized to the PCFs using the iterative Boltzmann inversion (IBI) method. The CGMD simulations of xsPCHD chains using potentials from above method satisfactorily reproduced the polymer−polymer PCFs from atomistic MD simulation of hydrated membrane system at each hydration level. The transferability of above two set of CG potentials was further tested through CGMD simulation of hydrated membrane at an intermediate hydration level (λ = 15). Limited transferability was observed, presumably due to the use of an implicit solvent. Using an analytical theory, proton conductivity was calculated and compared with that from experimental measurement under similar conditions. Good agreement was obtained using inputs from both atomistic and CG simulation. This study provides a molecular level understanding of relationship between membrane structure and water and H3O+ ion transport in the xsPCHD membrane.

1. INTRODUCTION The unique phase separation behavior of cross-linked and sulfonated poly(1,3-cyclohexdiene) (xsPCHD) polymer in xsPCHD homopolymer, xsPCHD/PEG (poly(ethylene glycol)) blend, and xsPCHD−PEG block copolymer proton exchange membrane (PEM) generates higher proton conductivity than that of Nafion membrane at high temperatures.1 Such phase behavior and many other properties strongly depend on the conformation of the polymer in solution or bulk.2,3 Information on the molecular level (10−100 nm) membrane structures furthers the understanding of the mechanism of high proton conductivity in these systems. Different molecular modeling techniques have been used to study the structures of PEMs. Classical molecular simulation has proved to be useful in the study of chain conformation when the length scale is less that 10 nm roughly.4,5 Previously, nanoscale structures of hydrated perfluorosulfonic acid (PFSA) PEMs, including Nafion, have been intensively studied though molecular simulations.6−16 The pair correlation functions (PCFs) of water−water, water−hydronium, hydronium ion− hydronium ion, polymer−water, polymer−hydronium ion, and polymer−polymer pairs provided very useful information on the water and hydronium ion distribution in the membrane, which is crucial on generation of good conductivity.6,7 © XXXX American Chemical Society

Because of the computational limitations of classical molecular simulation, multiscale modeling is needed to study PEM structure at higher length scale. The coarse-grained (CG) based multiscaling modeling technique has been developed to study the structural and transport properties of long chain polymers at length scale that are larger than the classical molecular simulation can reach. In the CG procedure, normally a fully atomistic (or united atom (UA)) simulation of a short-chain (oligomer) system are conducted first. From the atomistic simulation, the CG structural distribution functions for different interactions modes assigned in the CG model are obtained. The CG potentials are then parametrized to these distribution functions. There are a variety of methods used to obtain the effective potentials. The iterative Boltzmann inversion (IBI) method17 is widely used because of its advantage to reproduce the structure of atomistic sampling compared to the existing methods, especially when complicated nonbonded structures are obtained from the atomistic simulation.18 The CG modeling work on polymer solutions, suspensions, or hydrated proton exchange membranes (PEMs) shows that the Received: February 23, 2012 Revised: July 18, 2012

A

dx.doi.org/10.1021/ma300383z | Macromolecules XXXX, XXX, XXX−XXX

Macromolecules

Article

which nonbonded pair correlation functions (PCFs) of water− water, water−hydronium ion, hydronium−hydronium ion, polymer−water, polymer−hydronium ion, and polymer−polymer pairs are obtained and studied. Second, the nonbonded PCFs for xsPCHD polymer from the atomistic simulations were used to generate CG potentials, using the IBI method. The CG potentials are implemented in CGMD simulations of the same system while water and H3O+ ion are treated implicitly. The nonbonded PCFs for xsPCHD polymer from CG level simulation are compared with those from the atomistic level simulation in order to validate the CG model. Third, the transferability of the derived CG potentials is studied at an intermediate hydration level (λ = 15). Finally, water and H3O+ ion transport in the membrane are studied using an analytical theory. Proton conductivity is calculated from the theory and compared with that of experimental measurement at different water contents.

solvent plays an important role in determining the polymer structure.19−23 However, there is no well-established method to obtain the CG effective potentials for solvent−polymer and solvent−solvent interactions to reproduce the structures from atomistic sampling when the solvent is treated explicitly in these systems. One may easily encounter convergence problems when apply the IBI to these inhomogeneous systems since interactions between beads of the same type may change with the local environment. Alternatively, a successful way to deal with solvent in these works is by treating the solvent explicitly at the atomistic scale while treating it implicitly at the CG scale. Since the polymer structure at the atomistic level is reproduced using the IBI method at CG level, the effect of the solvent on the polymers structure is still accounted for in the CG simulation. Using this method, Fischer et al.21 developed an implicit CG model for polyoxylethylene solution and studied the potential transferability to different chain length and mixture concentrations. The same method is adopted by Bedrov et al.22 in their CG modeling work on poly(ethylene oxide)−poly(propylene oxide)−poly(ethylene oxide) triblock copolymer in solution and Grest et al.23 in their work on nanoparticles in suspension. Depending on different environments (concentrations, CG models, densities, temperatures), CG potentials show different transferability.5,21,23,24 One has to test the transferability of derived CG potentials before applying them to different systems, such as PEMs with different hydration levels. The membrane structure dictates the water and charge transport in the membrane. Recently, an analytical theory has been developed to provide an understanding of water and charge transport in aqueous systems.25,26 The theory account for acidity, confinement, and connectivity, three structural factors that affect the water and charge transport in aqueous system, including PEMs. There are numerous ways that the disordered morphology of the hydrated polymer membrane could be characterized. Many of these different characterizations are not independent of each other. In previous work we sought to distill the characterization of the structure into as few independent measures as possible. We found three measures were sufficient: acidity, confinement, and connectivity. It is absolutely the case that the values for acidity, confinement, and connectivity are related to values of other characterizations such as backbone hydrophobicity and rigidity, and vice versa. In other words, the choice of traits is possibly not unique, and other models could perhaps be postulated based on other choices. By making the model a function of three independent structural traits, it is implicitly a function of other traits which are correlated with the selected three traits. Molecular simulations provide structural information for these three factors. The theory has successfully described self-diffusivity of water and charge across the full range of hydration for hydrated Nafion membrane. Theoretically, it can be applied to the study of water and charge transport in any hydrated PEMs, including xsPCHD membranes. Once the charge self-diffusivity is faithfully predicted, the proton conductivity can be calculated and further compared with that of experiments at different hydration levels. In this work, we studied the structures of hydrated xsPCHD membrane at the atomistic and CG level. Previously, a multiscale model for sulfonated cross-linked xsPCHD homopolymer was developed for the unhydrated state (a melt).27 The same CG procedure is used here. First, an atomistic MD simulation of the hydrated short-chain xsPCHD system at two different water contents (λ(H2O/HSO3) = 10 and 20) was conducted, from

2. SIMULATION METHOD 2.1. Atomistic Simulation of Hydrated xsPCHD Membrane with λ = 10 and 20. The force field used for xsPCHD polymer in the melt,27 which invoked potentials from the literature for each of the components in the polymer,6,28−36 was largely used again here. The model is fully atomistic with the exception of hydrogen atoms bonded to carbons, which are modeled as CHx groups, as shown in Figure 1a. The atomistic

Figure 1. Atomistic (a) and CG (b) representations of xsPCHD polymer.

model includes bond stretching, bond bending, bond torsion, and intramolecular and intermolecular nonbonded interactions via both Lennard-Jones and Coulombic potentials. The only difference between the atomistic potential used in the simulation of the melt and that used in the simulation of the hydrated polymer is that the SO3H group in the melt is dissociated into SO3− and H3O+ ions in the hydrated system. The SPC/Fw model by Wu et al.37 was used for water, and a fully flexible model of H3O+ used in hydrated Nafion simulations by Cui et al.6 work was used. The number of molecules in the simulation systems are summarized in Table 1. The spherically truncated chargeneutralized method developed by Wolf et al.38was used to evaluate the electrostatic energy. We first simulated in the isobaric−isothermal (NpT) ensemble and implemented the Hamiltonian-based thermostat39,40 and barostat41 with controller frequencies set to 10−4 fs−1. The XI-RESPA NPT algorithm developed by Tuckerman et al.42 was used to integrate the equations of motion. The large time step was 1 fs, and the B


Macromolecules

Article

stretching, bond bending, bond torsion, and nonbonded interactions for both intramolecular (for beads separated by at least four bonds) and intermolecular pairs. The same CG nonbonded potential is used for both intramolecular and intermolecular interactions. The atom fragments that make up the A and B beads carry a negative charge. The corresponding positive charge is carried in the solvent, which is implicit in the CG simulations. One cannot have a Coulombic interaction with an implicit solvent. Therefore, the CG potentials for AA, AB, and BB modes include all interactions: repulsion, dispersion, and electrostatic. There are not explicit charges in the CG model. Probability distribution functions (PDFs) for stretching, bending, torsion, and nonbonded pair correlation functions (PCFs) between CG beads were generated from the atomistic simulations of the hydrated xsPCHD membrane. As included in the Supporting Information, a comparison of bonded PDFs from atomistic simulations of both melt and hydrated xsPCHD systems shows that there are no significant changes in the bonded PDFs. Consequently, the bonded CG potentials for the xsPCHD melt are directly used in the CG simulation of hydrated xsPCHD membrane. The nonbonded potential for intramolecular (for beads over four bonds) and intermolecular interactions is generated by the IBI method.17 In the IBI method, a tabulated potential is numerically determined by simulation iteration. The interaction potential is refined iteratively via

Table 1. Number of Molecules for Each Component, Densities, and Times in the Atomistic Simulation Systems λ

no. of polymers

no. of H2O

no. of H3O+

total no. of particles

density (g/cm3)

simulation time (ns)

10 20

40 40

6480 13680

720 720

35920 57520

1.14 1.05

10 4

small time step was 0.1 fs. The cutoff distance used was 15 Å for this level simulation. The state point was set at 1 atm and 353 K, corresponding to future fuel cell applications. Following the equilibration procedure described elsewhere,5,24 we gradually equilibrated to the correct density. Data production lasted for 10 ns for λ = 10 and 4 ns for λ = 20 system in NVT simulations. The equilibrium configuration is shown in Figure 2. To compare the

⎛ g (r ) ⎞ αβ , i ⎟ ϕαβ ,i + 1(r ) = ϕαβ ,i(r ) + kBT ln⎜⎜ ⎟ g ( r ) ⎝ αβ ⎠

(1)

with initial guess ϕαβ ,0(r ) = −kBT ln(gαβ (r ))

(2)

where gαβ(r) is the target PCF, which is from atomistic simulation. In CG level simulation, water and H3O+ ion are treated implicitly, as commonly adopted in the CG modeling work of polymer solutions.21,22 Because of the complexity of CG model, there are still 10 nonbonded interactions modes (AA, AB, AC, AD, BB, BC, BD, CC, CD, and DD) that need to be refined. First, two iterations were done on all interactions modes simultaneously, and then an additional 1−3 iterations were done on each mode individually. The number of iterations used for each mode depends on the complexity of the PCF. Generally, each mode took 3−5 iterations to reproduce the structure of atomistic sampling, which is chosen to stop the iteration cycle. A pressure correction is not applied since the simulations are conducted in the NVT ensemble and pressure is not reported. The whole process took roughly one month. It is computationally expensive to apply this approach to a complicated CG model. Future work will involve improving the CG model and potential generation approach to speed up the CG procedure. In the CGMD simulations, we use the same density as that used in the atomistic simulation. A larger system of 64 chains is used at this level of simulation. We again simulated in the NVT ensemble under the same temperature as the atomistic simulation. The cutoff distance used is 25 Å for the CG simulation. The time step of the CGMD simulation is twice that used in the atomistic MD simulation. Based on wall-clock time, the CGMD simulations are ∼400 times faster than the atomistic simulations.

Figure 2. Snapshots of equilibrium configurations from atomistic (a−c, e−g) and CGMD (d, h) simulations of xsPCHD membrane at T = 353 K, p = 1 atm with λ = 10 (a−d) and λ = 20 (e−h). In atomistic simulation, water and hydronium ion are treated explicitly (a, e). From (a) and (e), water and hydronium ion (b) and (f) and xsPCHD polymer (c) and (g) are taken, respectively. In CG simulation, water and hydronium ion are treated implicitly.

water−water distribution, atomistic simulation of bulk water was also conducted. The same water model and simulation method was used. Data production lasted for 5 ns. 2.2. CGMD Simulations of Hydrated xsPCHD Membrane with λ = 10 and 20. As shown in Table 1 and Figure 1b, the CG model for xsPCHD proposed in previous work27 is used here although the hydrogen on SO3H is transferred to H3O+. There are three CG beads on xsPCHD chains (bead A, B, and C) and two cross-linkers (D) at the two ends. The bead is placed at the center-of-mass of the atoms in the corresponding fragment. The mapping scheme is the same as that adopted in the previous work. In this model, the CG interaction modes include bond C


Macromolecules

Article

3. RESULTS AND DISCUSSION In this section, we discuss (1) the results of structures from atomistic MD simulation of hydrated xsPCHD membrane with λ = 10 and 20, (2) comparison of xsPCHD polymer structures from atomistic simulation of hydrated membrane with different water contents, (3) generation of CG potentials for hydrated xsPCHD membrane with λ = 10 and 20, including the test of the transferability of the CG potentials, and (4) study of water and charge transport in hydrated membranes with different water contents. In Figure 2, we show equilibrated snapshots from both the atomistic and the CGMD simulations for hydrated membrane with λ = 10 and 20. In Figures 2a and 2e, all molecules are shown to make it clear that we are simulating a multiphase system in which nanoscale segregation of the hydrophobic (polymer backbone) and hydrophilic (aqueous) phases is present. In Figures 2b and 2f, all polymer chains are rendered invisible to better indicate the shape of the aqueous domain. In Figures 2c and 2g, all water molecules and H3O+ ions are rendered invisible to better show the shape of the polymer domain. In Figures 2d and 2h, the xsPCHD structure at the CG level is shown. Interactive structures are available on an archived Web site.43 One of the features that is clearly visible in these snapshots is the segregation at the nanoscale into hydrophobic (polymer) and hydrophilic (aqueous) domains. The sulfonate ions are located at the interface between the hydrophobic and hydrophilic domains. This sort of nanophase segregation has been observed for Nafion by many researchers using molecular simulation.7−9,12 Characterization of the volume of the aqueous domain, the interfacial surface area, and the connectivity of the aqueous domain depicted qualitatively in these snapshots has been shown to have direct impact on the mobility of water and protons in the membrane.26 Experimental studies of xsPCHD have shown that water channels are visibly connected at a 10 nm scale while the polymer domain is at 10−100 nm length scale.1 Our simulations have a dimension of 10−15 nm. Therefore, they are capable of capturing the structure of the aqueous domain but not of the polymer domain. Work with Nafion has shown that MD simulation of these smaller systems can still yield quantitive agreement with experimental measurements of water and charge self-diffusivities.26 We would not expect simulations of this size to capture any crystallinity effects (as exist in Nafion) in the polymer phase, which if present span longer length scales. 3.1. Aqueous Phase Structure from Atomistic Simulation. In Figure 3, we show the water−water, H3O+−H3O+, and water−H3O+ PCFs of the CG beads obtained from atomistic simulation. These distribution functions are based on the analysis of configurations from atomistic MD simulation. Although the aqueous phase will be treated implicitly at the CG level, these PCFs provide useful information on water and H3O+ distributions in the membrane. The following PCFs are based on positions of the center-of-mass of water molecules and hydronium ions. In Figure 3a, the water−water PCF for the hydrated membranes is compared to that of bulk water. There is significant difference in the water−water structure, including an enhanced first peak and the presence of a small second at about 4 Å. These same features are present in the water−water PCF in hydrated Nafion at λ = 9,7 in which similar nanoscale segregation occurs. It is worth noting that, following standard convention, the PCFs are normalized by the bulk density.

Figure 3. Pair correlation functions (PCFs) for CG beads involved in (a) water−water, (b) H3O+−H3O+, and (c) water−H3O+ interactions (black line: λ = 10; red line: λ = 20; water−water PCF of bulk water is also shown in (a)).

In a homogeneous system, like bulk water, there is no ambiguity associated with the bulk density. In an inhomogeneous system, such as the hydrated xsPCHD membrane, we have normalized by a density that includes the total system volume (since that is unambiguously available). Consequently, some of the discrepancy between these PCFs is due to the definition of the PCF. The PCF here is based on center-of-mass position of a single water molecule. The features on the water−water PCF will change if a different CG model of water is used.44 We note that aqueous D


Macromolecules

Article

Figure 4. PCFs for CG beads involved in polymer−water interactions (a, c) and polymer−H3O+ interactions (b, d). (a, b) λ = 10; (c, d) λ = 20.

structures, as shown in Table 1, we expect the correlation with water to be similar. The presence of the hydroxide group on the B bead lowers the correlation. We see the weakest correlation with particles of type D, which contain the cross-linking groups and are consequently surrounded by more polymer than water. In Table S.1 (Supporting Information), the hydration number of A, B, C, and D beads is calculated from polymer−water PCFs within the distance of 5.2 Å. The distance 5.2 Å was chosen to correspond with the end of the first peak. The values of the hydration numbers shows the exact same trend as that observed from Figure 4a. There are on average between 5 and 6 water molecules around each A and B beads within a radius of 5.2 Å. There are 3.26 water molecules within the same radius of C beads and 2.0 water around D beads, reflecting both their nonionic character. In Figure 4b, we show the PCFs for polymer−H3O+ ion pairs. For all bead types, the peak heights are greater for the hydronium ion than they are for water molecules in Figure 4a. This indicates that H3O+ ions are more strongly bound to polymer than are the water molecules, due to the stronger Coulombic attraction between the hydronium ion and the sulfonate anion. In Figure 4b, there is an apparent difference between the PCFs of A and B. The major peak around 5 Å is split into two peaks in the B PCF. We believe this is due to the existence of a hydroxyl group on bead B. The orientations of some H3O+ ions are likely altered due to hydrogen bonding with the hydroxyl group. In Table S.2 (Supporting Information), the coordination number of A, B, C, and D beads with hydronium ions is calculated from the PCFs within the distance of 5.4 Å. We see coordination numbers for A and beads close to unity for low hydration. The coordination number of B bead is slightly larger than that of A bead, which is again probably due to the existence of hydrogen bonding

phase is treated as implicit at CG level. Therefore, the aqueous phase structures are only studied based on atomistic simulation results. In Figure 3b, the H3O+−H3O+ PCF is shown. There is very little structure in the PCF, commensurate with the fact that two cations in solution will not aggregate, due to Coulombic repulsion. What little structure there is can be attributed to the fact that the anions are fixed in close proximity to each other due to their interactions with the sulfonate anions attached to the polymer backbone. Since there may be association between the cation and anion, some weak cation−cation association may be evident. In Figure 3c, the water−H3O+ ion PCF is shown. We observe two peaks in this PCF. The first, larger peak occurs at about 2.60 Å, and the second peak occurs at a distance of 4.95 Å. The same feature has been found in the atomistic simulation work of Nafion at different water contents,6,7 in which the PCF between oxygen of water and oxygen of H3O+ is used to represent water− H3O+ ion PCF. The first peak is much higher that the second one because of the tighter binding of water molecule with H3O+ ion at short distance. This feature is also in agreement with above atomistic modeling work of hydrated Nafion. In Figure 4, we examine the PCFs obtained from atomistic simulations between the center-of-mass of various polymer fragments corresponding to CG beads with water molecules and hydronium ions. Figure 4a shows the PCFs between the beads on the polymer and a water molecule at the lower level of hydration, λ = 10. We see the strong correlations with beads of type A and B because these beads contain SO3− ions, which interact strongly with water molecules. Both A and B beads show a first peak around 5.0 Å, a second peak around 7.0 Å, and a third peak around 14.0 Å. Since beads A and B have very similar chemical E


Macromolecules

Article

Figure 5. Comparison of PCFs from atomistic simulation of xsPCHD as dry melt (red line) and hydrated with λ = 10 (black line) and λ = 20 (blue line) conditions. Significant changes are observed.

Figure 4c shows the PCFs between the center of mass of atom fragments in the atomistic simulation that will correspond to the CG beads on the polymer and water molecules at the higher level

between hydroxyl group in B bead and hydronium ion cluster, which causes the hydronium ion to be locally more stable bound to a B bead. F


Macromolecules

Article

of hydration, λ = 20. We see all the peak positions of polymer− water PCF at λ = 10 are reproduced at λ = 20. The magnitude of the peaks of polymer−water PCF at λ = 20 at are lower that that at λ = 10. Analogous behavior for the decrease in the height of the first and second peak of the SO3−−water PCF as a function of water content has also been observed for Nafion and SSC PFSA irrespective of EW.6−8 In other words, the water molecules are less tightly bound to the anions with increasing water content. One cannot exclusively analyze the hydration of the polymer beads as a function of water content based on the PCF alone because the PCF is normalized by the system density, which varies with water content. Therefore, hydration numbers, integrals of the PCF up to a given distance (5.2 Å in this case), are presented in Table S.1. We observe a slight increase in the hydration number of each bead with increasing water content. In Figure 4d, we show the PCFs for polymer−H3O+ ion pairs at λ = 20. Similar to that of polymer−water PCFs, all peaks of the polymer−H3O+ ion PCFs at lower hydration level appear at the higher hydration level with reduced magnitudes. From Table S.2, we see the coordination number of each type of bead with hydronium ions at λ = 20 is almost half of that at λ = 10 and less than unity since the hydronium ions are moving further from the sulfonate group into the aqueous phase. 3.2. Polymer Phase Structure from Atomistic Simulation. To understand the effect of hydration on polymer structure, we compare the structures of xsPCHD in the hydrated membrane with different hydration levels. The bonded stretching, bending, and torsion distributions of xsPCHD melt and hydrated membrane are shown as Figures S.1, S.2, and S.3 in the Supporting Information. No apparent change on the structures is observed. The CG nonbonded PCFs are shown in Figure 5. In all cases, these distribution functions were generated from atomistic simulations of short chains. In Figure 5, we compare the PCFs for each of the 10 polymer− polymer nonbonded modes for the melt (λ = 0) and for the hydrated membrane at λ = 10 and λ = 20. As a reminder, the A and B beads in the melt contain the entire sulfonic acid group, whereas the A and B beads in the hydrated membrane contain only the sulfonate anion. For those modes that involve ion−ion interactions in the hydrated state (AA, AB, and BB) there are significant differences in the PCFs. The first peak in the melt is reduced significantly in magnitude and pushed to slightly larger distances in the hydrated system due to electrostatic repulsion. Some modes are relatively unchanged by the introduction of water, including the AC, AD, BC, and CC modes. However, the structure around the D bead, the cross-linking group, has changed significantly, resulting in very different PCFs for BD, CD, and DD. We also note that the nonbonded interaction modes are correlated with each other. A change in the distribution of one mode may cause changes in the distributions of the other modes. 3.3. Comparison of Structures and Potentials from CG Simulation. One aim of this research is to derive a CG force field for xsPCHD in a hydrated membrane. Previously, we derived a CG model for xsPCHD in the melt using the IBI method.45 First and foremost, we cannot use the CG nonbonded potential from the melt for the hydrated state due to the fact that A and B are now charged beads in the hydrated state. These beads are fundamentally different. Second, the effective potentials obtained from the IBI method are dependent on thermodynamic state (density and temperature).17 Its transferability to temperature, chain length, or density needs to be tested.24 Several studies have shown in general that CG potentials

are transferable with respect to chain length.18,24,27,46,47 The previous study on xsPCHD melts demonstrated that the PCFs did not have any significant chain length dependence, thus establishing its transferability with respect to chain length.27 In this work, we therefore did not explicitly simulate longer chain lengths. We simulated the xsPCHD hydrated membrane at the same temperature and pressure as that of melt xsPCHD simulation. The density of the xsPCHD melt and hydrated membrane are of course different. In the previous section, we compared the nonbonded PCFs for xsPCHD in the melt and in the hydrated membrane. The nonbonded xsPCHD PCFs show visible changes with different water contents (λ). In the work of aqueous poly(oxyethylene) solutions, Fischer et al.21 show that the nonbonded PCFs have similar trends but apparent differences at different water contents, which depend very much on the mixture and CG model used. In this case, one has to derive CG potentials from systems at varying water contents and validate the transferability. Based on this work, the CG potentials may be different for xsPCHD membrane with different water contents, and it is quite necessary to develop different sets of CG potentials based on these PCFs and test their transferability. To do so, we made the aqueous phase implicit at the CG level as mentioned above. Thus, 10 polymer−polymer nonbonded modes were iteratively refined to reproduce the PCFs from atomistic simulation. Effective potentials from the last iteration are shown in Figure 6. These effective potentials are quite different from the potentials used in atomistic simulation, since the effective potentials are generated from the PCFs. For comparison, also shown are the CG potentials from the xsPCHD melt. The first obvious observation is that there is significant change in the CG potentials between the hydrated and melt states. The most distinctive difference occur for AA, AB, and BB, modes in which the beads have become charged in the hydrated state. In these cases the depth of the attractive well has been greatly reduced, and the position of the minima has been pushed to larger distances. For interactions between nonionic beads, such as CC, the well depth has increased in the hydrated state. This can be explained by remembering that the water is treated implicitly. In an atomistic simulation, we observe phase segregation into hydrophobic and hydrophilic regions because the water− water interaction is much stronger than any other interaction. However, in the CG simulation where water is not explicitly present, the mechanism for phase segregation is the attraction of hydrophobic beads, like C. This observation is also observed for the strictly uncharged CD and DD modes. For modes involving one charged bead and one uncharged bead, no uniform difference between the potentials of the hydrated and melt states is observed. To further understand the features of these potentials, we can revisit the PCFs from which the potentials were extracted. In Figure 5, the first peaks of the CC, CD, and DD (uncharged− uncharged) modes are sharper than that of the other modes for the hydrated system. Correspondingly, the potentials of CC, CD, and DD show deep first minima in Figure 6. As noted above, the PCFs of the charged−charged modes (AA, AB, and BB) show much smaller and more distant first peaks, which corresponds to the observed changes in the CG potential. The PCFs of the charged−uncharged modes (AC, AD, BC, and BD) modes are relatively unchanged by the introduction of water, with the exception of BD. However, the CG potentials may be significantly changed even in the absence of strong changes in the PCF, e.g., AD. G


Macromolecules

Article

Figure 6. Comparison of CG potentials from IBI method for xsPCHD under melt (red line) and hydrated with λ = 10 (black line) λ = 20 (blue line) conditions. Significant changes are observed.

validation of the bonded stretching, bending and torsion potentials are provided in previous work.27 In general, the agreement between atomistic and CG scales is very good.

In order to validate the CG potential, one should compare structures obtained from the atomistic and CG simulations of the same system. That validation has been performed here. The H


Macromolecules

Article

Figure 7. Comparison of PCFs from atomistic (red line) and CG simulation (black line) using effective potential from the IBI method for hydrated membrane with λ = 10.

modes for each case. All peaks present in the atomistic simulation are also present in the CG simulation. For the PCFs in which multiple peaks are present, the relative magnitude of the peaks observed in the atomistic simulations is reproduced by the CG simulations in all cases except for the

We also validate the CG model through comparison of the nonbonded PCFs. In Figures 7 and 8, a comparison of PCFs from the atomistic and CG simulations of the hydrated membrane is provided for hydration levels of λ = 10 and 20, respectively. In general, good agreement is obtained for all I


Macromolecules

Article

Figure 8. Comparison of PCFs from atomistic (red line) and CG simulation (black line) using effective potential from the IBI method for hydrated membrane with λ = 20.

indicates the convergence of the IBI method, though guaranteed by theory, is difficult to achieve in practice, especially when a complex CG model is used. However, we

CC mode. There is admittedly some quantitative discrepancy for a few of the modes (AC, CC, and CD). Further iterations provide little improvement on matching of these modes. This J


Macromolecules

Article

Figure 9. Comparison of PCFs from CGMD simulations of xsPCHD under λ = 15 using CG potentials obtained from atomistic simulations of λ = 20 (red line) and λ = 10 (black line) conditions.

believe that this level of agreement is sufficient to allow the CG simulations of hydrated sxPCHD to capture the essential structures of the atomistic simulations.

As an important aspect, the transferability of CG potentials from the IBI method needs to be carefully checked. Previous studies have shown that the CG potentials are very concentration K


Macromolecules

Article

Figure 10. Comparison of PCFs from CGMD simulations of xsPCHD large (red line) and small (black line) systems using CG potentials obtained from atomistic simulations of λ = 10.

λ = 20 since the two sets of potentials obtained at those hydration levels have significant differences, as shown in Figure 6. In this work, we also tested a more limited transferability of the so

dependent,48 and limited concentration transferability is observed21 based on the comparison of the PCFs. One can observe that the CG potential is not transferable from λ = 10 to L


Macromolecules

Article

Figure 11. Snapshots of configurations from CGMD simulations and image from AFM measurement for xsPCHD membrane: (a) from CGMD simulation under λ = 15 condition, using CG potentials from λ = 10 condition; (b) from CGMD simulation under λ = 15 condition, using CG potentials from λ = 20 condition; (c) from CGMD simulations of large dimension system under λ = 10 condition, using CG potentials from λ = 10 condition; (d) AFM image from experiments.

derived two sets of the CG potentials through the use of both potentials in CGMD simulations on a system with an intermediate hydration level of λ = 15, where the aqueous phase is still treated as implicit. The CGMD is run in a NVT ensemble with equilibrium density estimated through an empirical relation.7 The comparison of the PCFs from CGMD simulation using the two sets of the potentials is shown in Figure 9. The potentials generate all the peak positions for all the 10 nonbonded modes. However, we observed great deviations on the magnitudes of the peaks. Further examination of the configure snapshots shows that very different system configurations are generated from the two set of the potentials (Figure 11a,b). We believe that here the derived CG potentials are also strongly concentration dependent. One reason could be we applied a higher CG level here (by mapping 6 or 7 united atoms as a CG bead and using implicit solvent). Good concentration transferability of CG potentials has been reported by Fischer et al.21 in their modeling of aqueous poly(oxyethylene) solutions, where a finer mapping scheme and the including of local degree of freedoms are assumed to be the reason. The goal of developing the CG model is to apply it to study the membrane structure at larger length scales. Previously, we have shown that the model for xsPCHD is able to predict the structure of longer chain systems for the anhydrate membrane. We chose the CG potential set based on λ = 10, but enlarged the

dimension of the system to be 8 times of the original system in volume by putting more xsPCHD polymer in the simulation box. The box length reached 18 nm. The density, temperature, and pressure are the same as that of the original small system at λ = 10. The CGMD simulation lasted for 20 ns* (unscaled CG time). We first compared the PCFs from the large system to that of small system, as shown in Figure 10. Since nothing but the length scale changes, significant changes on the PCFs are not expected. No significant changes on the 10 nonbonded PCFs are observed when the length scale increases. The PCFs from CGMD simulations of the large and small systems show very similar trends. This indicates the CG potentials can be used to study the structures of the larger system. At this larger length scale, the AFM image has shown that the polymer domain (bright region) and aqueous domain (dark region) are segregated, as shown in Figure 11d. The aqueous channels are roughly around 10 nm in diameter. The CGMD simulation on the large system captures the features on the membrane morphology. Figure 11c shows the snapshot (equilibrium configuration) from the CGMD simulation on xsPCHD membrane at length scale of 18 nm. Clearly, we observed the segregation between the polymer and aqueous domains (empty space in the box since the aqueous phase is treated as implicit at CG level). On the basis of the observations on AFM images, we are most probably able to observe one roughly cylindrical water channel at the simulation length scale M


Macromolecules

Article

(18 nm). In Figure 11c, we observed such a cylindrical channel. In addition, the cylindrical volume in the middle of the snapshot shows a diameter of about half of the length of the simulation box, about 9 nm. This agrees with the experimental observation on the diameter of the water channels (around 10 nm). We shall show below that the surface area and free volume ratios of the small and large CGMD simulation at λ = 10 are consistent. 3.4. Water and Charge Transport in Hydrated Membranes with Different Water Contents. For crosslinked xsPCHD polymer, the chain relaxation is very slow. With finite computational time, the atomistic simulation can not provide reliable information on chain relaxation and dynamics. We observed a faster dynamics in the CG simulation. However, we cannot compare the dynamics with that of atomistic simulation due to above reason. Therefore, the dynamics of the xsPCHD polymer in the hydrated membranes are not presented. While the atomistic simulations are far too short to capture the relaxation times of even the very short chain polymers simulated here, they are capable of capturing relaxation of the solvent, which is much faster. Evidence for this comes from the ability of combined MD and confined random walk simulations generating water and charge self-diffusivities in Nafion found to be in quantitative agreement with experiment.26 Typically, the self-diffusivity is extracted via the Einstein relation from the slope of the mean-square displacement (MSD) versus observation time in the infinite time limit. The MSD of water and the vehicular component of the MSD of charge are plotted in Figure 12 as a function of observation time. (Note that in these simulations we do not include a proton-hopping mechanism and thus are limited to measuring only one component of the overall charge self-diffusivity.) In the long-time limit of the Einstein relation, there should be a linear relationship between the MSD and the observation time. In the log−log plot, we show the slope is indeed unity at long time, indicating that we have simulated sufficiently long to obtain reliable self-diffusivities of the water molecules and hydronium ions for systems with λ = 10 and 20. Invoking the Einstein relation then yields self-diffusivities of 1.18 × 10−9 m2/s for water and 2.57 × 10−10 m2/s for the vehicular component of the hydronium ion for the λ = 10 system, 2.58 × 10−9 m2/s for water, and 8.8 × 10−10 m2/s for the vehicular component of the hydronium ion for λ = 20 system. Using the same method, a value of 9.24 × 10−9 m2/s is obtained for bulk water at the same temperature and pressure. The self-diffusivity of water in xsPCHD membrane with λ = 10 is 12.8% of that in bulk water. For comparison, the self-diffusivity of water in Nafion membrane with λ = 11.83 is 9.7% of that in bulk water at 300 K.6 In both cases, the water diffusion coefficient is significantly smaller than in bulk water, due to increased acidity, increased confinement, and decreased connectivity in the membrane.26 In addition to MD simulations, confined random-walk (CRW) theory can also be used to study the diffusion in the presence of nanoscale confinement, like water in proton exchange membrane.25,26 The details of the CRW theory are provided in the Supporting Information. The theory contains two parameters describing confinement: a cage size and cage-to-cage hopping probably. The theory fit to MSDs of very short MD simulations can reproduce the MSDs of much longer MD simulations, where Einstein relation can be used to extract the self-diffusivity. Furthermore, the CRW analysis generates cluster sizes and cluster-to-cluster hopping probabilities. In Figure 13, we present the MSDs from both the theory and MD simulations for λ = 10 and 20 systems. For each of the systems, we showed the MSD data from two different MD simulations. The first set of the

Figure 12. Mean-square displacement (MSD) for water (a) and hydronium ion (b) as a function of observation time for hydrated membrane with λ = 10 and 20 (black lines: MSDs from atomistic MD simulations; red lines: MSD from curve fittings). MD simulations were run long enough to reach the long time limit. The diffusion coefficient (D) for each component is reported in Tables 2 and 3.

Table 2. CG Beads Representations in the Hydrated xsPCHD Membranea CG bead

molecular fragments

A B C D

(C6H9SO3−) (C6H10OHSO3−) (C6H10) (C6H10SCl)

a

The same symbols A and B in bulk CG xsPCHD model are still used while H atom in SO3H is dissociated in these beads.

Table 3. Properties from the Atomistic MD and CRW Simulations Applied to the Diffusion of Water in Hydrated xsPCHD Membrane with λ = 10 and 20a CRW theory λ

Rcage (Å)

pcage

D0 (10−9) (m2/s)

D (10−9) (m2/s)

MD simulation D (10−9) (m2/s)

10 20

15.8 26

1.41 × 10−2 3.8 × 10−2

3.29 4.23

1.21 ± 0.25 2.79 ± 0.40

1.18 ± 0.20 2.83 ± 0.12

a

In both cases, MD simulations were run long enough to reach the long time limit (Figure 12a).

simulations was only run for short observation time (0.2 ns), while in the second set, simulation was run long enough to reach N


Macromolecules

Article

these hydrated membranes, water and charge transport occur in the nanoscale aqueous domain. The model is parametrized to water data and used to predict proton conductivity. We parametrized the model as follows. The diffusivity through open windows, Do, is based upon the expression D(c , SA) = D(c = 0, SA = 0) exp( −k1c) exp(−k 2SA) (6)

where c is the charge concentration, defined as charges in aqueous volume, and SA is the interfacial surface area. The calculation of the surface area and aqueous volume adopted the method used by Liu et al.7 work on hydrated Nafion membrane. The constants D (c = 0, SA = 0), k1, and k2 are taken from Esai Selvan et al.'s work26 and are based upon experimental and simulated measurements of bulk and confined solutions. The diffusivity through blocked windows, Db, is given by eq 5, and the probability of cluster-to-cluster hops, pcage, is obtained from the CRW theory. The fraction of blocked pores, pEMA, is calculated to match the observed water self-diffusivities from the atomistic simulations, reported above. The same value of pEMA, which characterized the percolative structure of the membrane and which is obtained on the basis of water mobility, is used to predict proton conductivity. In the prediction of conductivities, both the vehicular and structural components of the charge self-diffusivities are described by the EMA model. The values of the structural and vehicular components of the charge self-diffusivity through open windows, Do, are determined from eq 6 with the appropriate sets of parameters. The values of pcage and pEMA determined above for water are used. The effective self-diffusivity of charge is generated via the EMA, eq 3. The conductivity is related to the selfdiffusivity via50

Figure 13. Mean-square displacement (MSD) for water as a function of observation time for hydrated membrane with λ = 10 and 20 (black dashed lines: MSDs from short time MD simulations; black solid lines: MSDs from long time MD simulations; red line: MSDs from CRW simulations). The diffusion coefficient (D) for each component is reported in Tables 2 and 3.

the infinite-time limit of the Einstein relation. Clearly, we observe that the MSDs from the theory match that from both the shorttime and long-time MD simulations. As shown in Table 3, the water self-diffusivity calculated based on MSDs from long-time MD simulation is also faithfully reproduced by the theory at each of the water contents. As water contents of PEM increases, both cage size and cage-to-cage hopping probability increase. This same trend has been reported on the modeling work of Nafion membrane.26 Compared to water diffusivity, charge diffusivity is more difficult to study. The charge diffusivity has both the structural (proton hopping) and vehicular components (cross term ignored),26 while the classical MD simulation can only calculate the vehicular component of the total charge diffusivity. An analytical theory has also been developed recently for charge transport in aqueous systems. The theory accounts for three factors: connectivity, acidity, and confinement that affect charge transport. The theory has successfully predicted the water and charge diffusivities for Nafion membrane with different water contents. Here we applied the theory to hydrated xsPCHD membrane with different water contents. A brief description of this theory is provided below. We use the effective medium approximation (EMA) from percolation theory to account for the effect of connectivity on self-diffusivity.49 Based on a bimodal distribution of diffusivities through blocked (b) and open (o) connections, the effective self-diffusivity from EMA can be analytically evaluated to yield ⎛ Deff 1 = ⎜⎜A + Do 2⎝

A2 +

z 2

⎞ 4f ⎟ − 1 ⎟⎠

A = 1 − pEMA + fpEMA −

D b = pcage Do

σ=

D |z|2 F 2c RT

(7)

In eq 7, all parameters are well-defined except for charge diffusivity, which is obtained from the above theory. After clarifying all the variables, we conducted the calculation and reported some interesting results in Tables 4 and 5. Table 4 Table 4. Interfacial Surface Area and Aqueous Volume Calculated for xsPCHD Membrane with Different Water Contents λ

description

surf. area (Å2/water)

aq vol (Å3/water)

10 10 10 15 15 20

atomistic CG-small CG-large CG CG atomistic

25.25 18.60 18.05 15.42a 15.64a 12.67

42.61 39.07 39.22 41.98a 41.85a 38.44

a

Results are scaled by a factor of 1.358 based on the ratio of the surface area from atomistic and CG simulation of the λ = 10 system.

(3)

f + pEMA − fpEMA z 2

−1

shows the calculated surface area and aqueous volume for xsPCHD membrane with different water contents from both atomistic (λ = 10 and 20) and CG simulations (λ = 10 and 15). The interfacial surface area and aqueous volume has been reported as monotonically decreasing with water content.7 We observed this trend from data shown in Table 4. Note here a scaling factor of 1.358 based on the ratio of the surface area from atomistic and CG simulation of λ = 10 system is used to scale the

(4) (5)

where f = Db/Do, Db is the diffusivity through a blocked window, Do is the diffusivity through an open window, and z is the coordination number of each cluster to the other clusters, chosen to be 6 here. The fraction of the blocked windows is pEMA. In O


Macromolecules

Article

4. CONCLUSIONS Atomistic and coarse-grained (CG) molecular dynamics (MD) simulations were conducted for a cross-linked and sulfonated poly(1,3-cyclohexadiene) (xsPCHD) hydrated membrane with λ(H2O/HSO3) = 10 and 20. From the atomistic level simulation results, nonbonded pair correlation functions (PCFs) of water− water, water−H3O+ ion, H3O+−H3O+, polymer−water, and polymer−H3O+ ion pairs are obtained and studied. The water self-diffusivity and H3O+ vehicular self-diffusivity are also obtained. Membrane structure was further studied at the CG level using a multiscale modeling procedure. Nonbonded PCFs of polymer−polymer pairs are obtained from atomistic simulation of hydrated membrane with λ = 10 and 20. Two sets of CG nonbonded potentials are then parametrized to the PCFs using the iterative Boltzmann inversion (IBI) method. The CGMD simulations of xsPCHD chains using potentials from above method satisfactorily reproduce the polymer−polymer PCFs from atomistic MD simulation of hydrated membrane system at each hydration level. The transferability of above two set of CG potentials is further tested through CGMD simulation of hydrated membrane at an intermediate hydration level (λ = 15). Poor transferability is observed, based on the comparison of the nonbonded PCFs. This is probably due to the relatively high coarse-graining level applied here. Some important degrees of freedom are neglected in the model. Specifically, the solvent, water, is treated implicitly. The CG potentials, which are pairwise between explicit CG particles, cannot compensate for changes in the relative amount of implicit solvent. This limits ability of the CG potentials to predict the structures of the membrane at different hydration levels. The CG potentials were useful in allowing the simulation of larger systems at the same hydration levels at which the potentials were generated. These larger systems give us clearer picture of the morphology of the aqueous domain. Using an analytical theory, proton conductivity is calculated and compared with that from experimental measurement under similar conditions. Good agreement is obtained using both the atomistic and CG simulations. The study provides a molecular level understanding of water and H3O+ ion distribution and transport in the membrane. The study on effect of different molecular structures of xsPCHD and xsPCHD based copolymer and blends on the conductivity is underway.

Table 5. Properties from the Analytical Percolation Theory and Experiments for Diffusion of Charge in Hydrated xsPCHD Membrane with λ = 10 and 20a theory −9

a

λ

pEMA

Dv (10 ) (m2/s)

Ds (10−9) (m2/s)

Deff (10−9) (m2/s)

σ (mS/cm)

expt σ (mS/cm)

10 20

0.4810 0.2412

0.257 0.88

0.43 1.89

0.21 2.14

25.5 147

18.13 125.8

Experimental data are taken from ref 1.

surface areas and aqueous volumes from CG simulations of λ = 15 system to the right scale. The use of scaling factor is required because the CG procedure reduces interfacial surface area. We did a comparison on the surface area for λ = 10 system, as shown in Table 4. The reduction is about 26% for surface area. The calculation of interfacial surface area and aqueous volume for λ = 15 system is based on configurations obtained from CGMD simulations using the two sets of potentials from λ = 10 and λ = 20. We also note that the small and large CGMD simulations at λ = 10 have the same surface area and free volume ratios, indicating that the more complete image of the morphology from the larger system in Figure 11d is structurally consistent with the image from the smaller system in Figure 2d. All components of the charge diffusivity in xsPCHD are shown in Table 5. We observed that the total charge diffusivity increases with increasing hydration for three reasons. First, the acidity (defined as excess protons per unit volume of aqueous domain) is lower with increasing hydration. Second, the confinement (defined as interfacial surface area per water molecule) is lower with increasing hydration. Both acidity and hydration in general tend to reduce the mobility of water and charge. Third, the connectivity of the aqueous domain is better in the more hydrated membrane, as indicated by the fact that the fraction of blocked windows, pEMA, has decreased by a factor of 2.0 Taken together, these three factors increase the self-diffusivity of charge by a factor of 11 from λ = 10 to λ = 20. In Table 5, we compared the proton conductivity from theory with that from the experiment at two water contents. The calculated conductivity is very close to that of experiment measurement for both the two cases, with average error of 28.8%. Given the approximate nature of the theory, the agreement is relatively good. Roughly the same level of accuracy was obtained by Esai Selvan et al.26 in the application of the theory to proton transport in Nafion. We also applied the same procedure to λ = 15 system, where the structural information are obtained from CGMD simulations using the two sets of potentials. The application the theory to the structure from CGMD using potential derived from λ = 20 system generates a conductivity of 77.4 mS/cm, which falls in the range of 25.5 mS/cm (for λ = 10) and 147 mS/cm (for λ = 20). The application of the theory to the structure generated from the CGMD simulations (surface area and aqueous volume given in Table 4.) for the λ = 10 system generates a conductivity of 27 and 24 mS/cm for the small and large CG systems, respectively, slightly higher than that (25.5 mS/cm) using structures from atomistic simulation. The major reason for this similarity is that the three factors that affect charge transport in PEMacidity, confinement, and connectivityare adequately reproduced by the CG model. Thus, we have successfully demonstrated the relationship between the structure of the hydrated xsPCHD membrane and used the structural information to predict the conductivity, which is used as a standard to characterize a good PEM.

■

ASSOCIATED CONTENT

S Supporting Information *

Additional tables and figures. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: dkeff[email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Chemical Sciences, Geosciences, and Biosciences Division of the Office of Basic Energy Sciences of the U.S. Department of Energy under Contract DEFG02-05ER15723. This research project used resources of the National Institute for Computational Sciences (NICS) supported by NSF under Agreement OCI 07-11134.5. S.D. and J.W.M. acknowledge support by the NSF-funded TN-SCORE program, NSF EPS-1004083, under Thrust 2. P


Macromolecules

■

Article

(36) Vishnyakov, A.; Neimark, A. V. J. Phys. Chem. B 2001, 105, 7830− 7834. (37) Wu, Y. J.; Tepper, H. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 024503. (38) Wolf, D.; Keblinski, P.; Phillpot, S. R.; Eggebrecht, J. J. Chem. Phys. 1999, 110, 8254−8282. (39) Nose, S. J. Chem. Phys. 1984, 81, 511−519. (40) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697. (41) Keffer, D. J.; Baig, C.; Adhangale, P.; Edwards, B. J. Mol. Simul. 2006, 32, 345−356. (42) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1993, 99, 2278−2279. (43) Keffer, D. J. Molecular Simulation Structures from the Computational Materials Research Group at the University of Tennessee Knoxville, TN, 2010. (44) Hadley, K. R.; McCabe, C. J. Phys. Chem. B 2010, 114, 4590− 4599. (45) Muller-Plathe, F.; Reith, D.; Putz, M. J. Comput. Chem. 2003, 24, 1624−1636. (46) Kremer, K.; Harmandaris, V. A.; Adhikari, N. P.; van der Vegt, N. F. A. Macromolecules 2006, 39, 6708−6719. (47) Harmandaris, V. A.; Reith, D.; Van der Vegt, N. F. A.; Kremer, K. Macromol. Chem. Phys. 2007, 208, 2109−2120. (48) Silbermann, J. R.; Klapp, S. H. L.; Schoen, M.; Chennamsetty, N.; Bock, H.; Gubbins, K. E. J. Chem. Phys. 2006, 124, 074105. (49) Keffer, D.; Davis, H. T.; McCormick, A. V. J. Phys. Chem. 1996, 100, 638−645. (50) Bockris, J. O. M.; Reddy, A. K. N. Modern Electrochemistry; Plenum Press: New York, 1970; Vol. I.

REFERENCES

(1) Mays, J. W.; Deng, S.; Mauritz, K. A.; Hassan, M. K.; Gido, S. P. Polyelectrolytes comprising sulfonated polydienes and poly(alkylene oxides) and related methods. US20090306295A1, 2009. (2) Yun, S. I.; Terao, K.; Hong, K. L.; Melnichenko, Y. B.; Wignall, G. D.; Britt, P. F.; Mays, J. W. Macromolecules 2006, 39, 897−899. (3) Yun, S. I.; Hong, K.; Melnichenko, Y. B.; Wignall, G. D.; Mays, J. W. Polym. Prepr. (Am. Chem. Soc., Div. Polym. Chem.) 2006, 47, 620−621. (4) Kremer, K. Macromol. Chem. Phys. 2003, 204, 257−264. (5) Wang, Q. F.; Keffer, D. J.; Nicholson, D. M.; Thomas, J. B. Macromolecules 2010, 43, 10722−10734. (6) Cui, S. T.; Liu, J. W.; Selvan, M. E.; Keffer, D. J.; Edwards, B. J.; Steele, W. V. J. Phys. Chem. B 2007, 111, 2208−2218. (7) Liu, J. W.; Suraweera, N.; Keffer, D. J.; Cui, S. T.; Paddison, S. J. J. Phys. Chem. C 2010, 114, 11279−11292. (8) Cui, S. T.; Liu, J. W.; Selvan, M. E.; Paddison, S. J.; Keffer, D. J.; Edwards, B. J. J. Phys. Chem. B 2008, 112, 13273−13284. (9) Devanathan, R.; Venkatnathan, A.; Dupuis, M. J. Phys. Chem. B 2007, 111, 8069−8079. (10) Glezakou, V. A.; Dupuis, M.; Mundy, C. J. Phys. Chem. Chem. Phys. 2007, 9, 5752−5760. (11) Idupulapati, N.; Devanathan, R.; Dupuis, M. J. Phys. Chem. B 2011, 115, 2959−2969. (12) Jang, S. S.; Molinero, V.; Cagin, T.; Goddard, W. A. Solid State Ionics 2004, 175, 805−808. (13) Karo, J.; Aabloo, A.; Thomas, J. O.; Brandell, D. J. Phys. Chem. B 2010, 114, 6056−6064. (14) Urata, S.; Irisawa, J.; Takada, A.; Shinoda, W.; Tsuzuki, S.; Mikami, M. J. Phys. Chem. B 2005, 109, 4269−4278. (15) Urata, S.; Irisawa, J.; Takada, A.; Shinoda, W.; Tsuzuki, S.; Mikami, M. J. Phys. Chem. B 2005, 109, 17274−17280. (16) Brandell, D.; Karo, J.; Liivat, A.; Thomas, J. O. J. Mol. Model. 2007, 13, 1039−1046. (17) Reith, D.; Putz, M.; Muller-Plathe, F. J. Comput. Chem. 2003, 24, 1624−1636. (18) Wang, Q. F.; Keffer, D. J.; Nicholson, D. M. J. Chem. Phys. 2011, 135, 214903. (19) Malek, K.; Eikerling, M.; Wang, Q. P.; Navessin, T. C.; Liu, Z. S. J. Phys. Chem. C 2007, 111, 13627−13634. (20) Malek, K.; Eikerling, M.; Wang, Q. P.; Liu, Z. S.; Otsuka, S.; Akizuki, K.; Abe, M. J. Chem. Phys. 2008, 129, 204702. (21) Fischer, J.; Paschek, D.; Geiger, A.; Sadowski, G. J. Phys. Chem. B 2008, 112, 13561−13571. (22) Bedrov, D.; Ayyagari, C.; Smith, G. D. J. Chem. Theory Comput. 2006, 2, 598−606. (23) Grest, G. S.; Wang, Q. F.; in’t Veld, P.; Keffer, D. J. J. Chem. Phys. 2011, 134, 144902. (24) Kamio, K.; Moorthi, K.; Theodorou, D. N. Macromolecules 2007, 40, 710−722. (25) Calvo-Munoz, E. M.; Selvan, M. E.; Xiong, R. C.; Ojha, M.; Keffer, D. J.; Nicholson, D. M.; Egami, T. Phys. Rev. E 2011, 83, 011120. (26) Selvan, M. E.; Keffer, D. J.; Calva-Munoz, E. J. Phys. Chem. B 2011, 115, 3052−3061. (27) Wang, Q. F.; Keffer, D. J.; Deng, S. X.; Mays, J. W. Polymer 2012, http://dx.doi.org/10.1016/j.polymer.2012.02.005. (28) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638−6646. (29) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. J. Am. Chem. Soc. 1996, 118, 11225−11236. (30) Jorgensen, W. L.; Ulmschneider, J. P.; Tirado-Rives, J. J. Phys. Chem. B 2004, 108, 16264−16270. (31) Jang, S. S.; Molinero, V.; Cagin, T.; Goddard, W. A. J. Phys. Chem. B 2004, 108, 3149−3157. (32) Dacquin, J. P.; Cross, H. E.; Brown, D. R.; Duren, T.; Williams, J. J.; Lee, A. F.; Wilson, K. Green Chem. 2010, 12, 1383−1391. (33) Abu-Sharkh, B. F. Comput. Theor. Polym. Sci. 2001, 11, 29−34. (34) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230−252. (35) Smith, J. C.; Karplus, M. J. Am. Chem. Soc. 1992, 114, 801−812. Q


Atomistic and Coarse-Grained Molecular Dynamics Simulation of a

Recommend Documents