Molecular Modeling of Proton Transport in the Short-Side-Chain

J. Phys. Chem. B 2008, 112, 2937-2949

2937

Molecular Modeling of Proton Transport in the Short-Side-Chain Perfluorosulfonic Acid Ionomer Iordan H. Hristov,† Stephen J. Paddison,*,‡ and Reginald Paul† Department of Chemistry, UniVersity of Calgary, 2500 UniVersity DriVe NW, Calgary, T2N 1N4 Canada, and Department of Chemical and Biomolecular Engineering, UniVersity of Tennessee, KnoxVille, Tennessee 37996 ReceiVed: NoVember 13, 2007; In Final Form: December 18, 2007

An explanation for the superior proton conductivity of low equivalent weight (EW) short-side-chain (SSC) perfluorosulfonic acid membranes is pursued through the determination of hydrated morphology and hydronium ion diffusion coefficients using classical molecular dynamics (MD) simulations. A unique force field set for the SSC ionomer was derived from torsion profiles determined from ab initio electronic structure calculations of an oligomeric fragment consisting of two side chains. MD simulations were performed on a system consisting of a single macromolecule of the polymer (EW of 580) with the general formula F3C-[CF(OCF2CF2SO3H)(CF2)7]40-CF3 at hydration levels corresponding to 3, 6, and 13 water molecules per sulfonic acid group. Examination of the hydrated morphology indicates the formation of hydrogen bond “bridges” between distant sulfonate groups without significant bending of the polytetrafluoroethylene backbone. Pair correlation functions of the system identify the presence of ion cages consisting of hydronium ions hydrogen-bonded to three sulfonate groups at the lowest water content. Such structures exhibit very low SsOH3+ separations, well below 4 Å and severely inhibit vehicular diffusion of the protons. The number of sulfonate groups in the first solvation shell of a given hydronium ion correlates well with the differences between Nafion and the SSC polymer (Hyflon). The calculated hydronium ion diffusion coefficients of 2.84 × 10-7, 1.36 × 10-6, and 3.47 × 10-6 cm2/s for water contents of 3, 6, and 13, respectively, show only good agreement to experimentally measured values at the lowest water content, underscoring the increasing contribution of proton shuttling or hopping at the higher hydration levels. At the highest water content, the vehicular diffusion accounts for only about 1/5 of the total proton transport similar to that observed in Nafion.

Introduction With the prospective of serving as a highly efficient and environmentally friendly power source in a broad range of applications including vehicular, stationary, and portable, the proton exchange membrane (PEM) fuel cell is the focus of research devoted to the engineering of catalysts based on nonnoble metals and an electrolyte that possesses chemical and mechanical stability accompanied with high proton conductivity (i.e., >0.1 S cm-1) under hot (i.e., >100 °C) and low humidity conditions. These requirements for the electrolyte notwithstanding the synthetic and characterization efforts of nearly two decades, have failed to discover a PEM superior to Nafion, a perfluorosulfonic acid (PFSA) ionomer with a polytetrafluoroethylene (PTFE) backbone that only possesses sufficient proton conductivity when fully hydrated. Clearly, the route to the successful development of a high performance fuel cell electrolyte will transpire through a molecular level understanding of how the chemical composition and degree of hydration of a PEM determine the transport of protons. Hence, not only is the presence of water critical to the function of the PEM fuel cell, but too much water in the membrane results in negative effects that impact both the design and efficient operation of the device. The fuel cell has to be either pressurized or a cooling apparatus employed to prevent * Author to whom correspondence should be addressed. E-mail: [email protected]. † University of Calgary. ‡ University of Tennessee.

evaporation of the water. Low operating temperatures (i.e., less than the condensation point of water at 1 atm) of the device necessitates the use of expensive Pt catalysts on the electrodes and purified hydrogen fuel to avoid CO poisoning. An electroosmotic drag of the water further aggravates the situation by drying out the anode and flooding the cathode. The hydration level determines not only the morphology of the polymer (i.e., the connected domains or pathways through which transport occurs) but also the number of available charge carriers by affecting the propensity of the sulfonic acid groups to dissociate. The efficient transport of these protons through a mechanism similar to structural diffusion of protons in bulk water (i.e., Grotthuss hopping or shuttling) is also dependent on the presence of water. How the hydrated morphology, itself a function of backbone chemistry and choice and density of the protogenic groups, determines the proton transport is not understood and therefore the goal of the present work. The utilization of classical molecular dynamics (MD) methods and simulations in the study of PEMs was recently reviewed,1 and hence the discussion here is limited to providing only a context for the present investigation. Vishnyakov and Neimark2 in an early MD study of Nafion of an equivalent weight (EW) equal to 1164 with potassium counterions observed isotropic clusters of up to 100 water molecules as the degree of hydration was increased but not a continuous hydrophilic phase (even at water contents equivalent to 17 wt %). They computed water and K+ diffusion coefficients with the Einstein relation but observed values for the former that were between 2 and 3 times

10.1021/jp7108434 CCC: $40.75 © 2008 American Chemical Society Published on Web 02/19/2008

2938 J. Phys. Chem. B, Vol. 112, No. 10, 2008 lower than the experimentally measured diffusivities. Spohr et al.3 also simulated mixtures of Nafion 1100 and water with classical MD combined with a simplified empirical valence bond (EVB) model. Their system consisted of two macromolecules of the polymer each possessing 40 monomer units and in the majority of the cases under an imposed initial system configuration of a slab. They observed a significant increase in the contribution due to Grotthuss shuttling of the protons when the water content was changed from 5 to 10 H2O/SO3- and did not observe any role of the polymer (i.e., dynamics of the backbone and/or side chains) in the transport of protons. Urata and coworkers in two separate investigations studied the swelling of Nafion membranes with neat water4 and aqueous methanol5 as a function of the weight percent of absorbed solvent utilizing a united atom approach in their MD simulations. In the former study, they were able to resolve static and dynamic properties of the water and side chains entirely consistent with earlier results obtained from calculations of the dielectric saturation of the water in PFSA pores6,7 and electronic structure calculations of the Nafion side chain.8,9 The latter investigation revealed that the presence of methanol not only affected gross morphological features including less spherical solvent aggregation but also the preferential location of CH3OH to the PTFE backbone (as opposed to association with the sulfonic acid groups). Jang et al.10 examined, through fully atomistic MD simulations, the differences between a random and a diblock copolymer of Nafion with an EW of 1150 consisting of 4 macromolecules each with a MW of 1.15 × 104 g/mol and water content of 15 H2O/SO3-. Although they did not observe any differences in the diffusion of hydronium ions they did witness differences in the distribution of the water within the polymeric matrix and water diffusion coefficients about 25% smaller in the random copolymer. Metiu and Voth, in two parallel and collaborative thrusts, utilized MD simulations to study the structure and dynamics of Nafion membranes with sodium counterions11,12 and the solvation and transport of hydrated protons in Nafion (acid form),13,14 respectively. Similar to previous classical MD simulations, their studies with Na+ cations (as opposed to protons) revealed that the hydrophilic phase was not connected at low (i.e., 5 wt %) hydration levels resembling “droplets” but formed more continuous “channels” as the water content was increased. They also observed the tendency of the Na+ ions to populate the center of the water phase at the high hydration in contrast to the electrostatic pining of the Na+ to several sulfonate groups at 5 wt %. They employed second generation (i.e., MSEVB2) and self-consistent multistate (SCI-MS-EVB) empirical valence bond methodologies to understand proton hopping in their investigations of the acid form of Nafion. The computed radial distribution functions (RDFs) showed the propensity for the formation of solvent separated Eigen-like cation/sulfonate pairs and a surprising anticorrelation between vehicular (i.e., as hydronium ions) and structure diffusion (proton shuttling). Dupius et al. very recently have examined the temperature and hydration effects including the dynamics of water molecules and hydronium ions of Nafion through extensive classical MD simulations.15-17 Their studies revealed the significant effects of temperature on the diffusion coefficients of both water and hydronium ions. There is only one reported computational study18 of the morphology and transport properties of a shortside-chain (SSC) (Hyflon) system. Despite some interesting parallels between Nafion, Hyflon, and the almost identical Aciplex membrane, the general applicability of the conclusions is not yet clear. A lack of averaging of the mean square displacements (MSD) of H3O+ and the entirely identical force

Hristov et al. field parameters used for the three systems have also lead to the qualitatively wrong conclusion that the mobilities of H3O+ are practically the same in Hyflon and Nafion. Previously, we developed a kinetic model based on a nonequilibrium statistical mechanical approach19-21 for direct determination of the effects of geometrical parameters of the water- and ion-containing pores within the membrane on the diffusion of hydrated protons. This “top down approach” relating membrane morphology to measurable properties is general and may be applied to any hydrated PEM system. The only systemspecific parameters used in the model are the radius of the pore, the protrusion of the anionic groups into the pore, and their distribution within the pore. These studies revealed that a larger pore radius, shorter side-chains and more uniform side-chain distribution will lead to higher proton diffusion coefficients.20 However, with this top down approach the effects of experimentally controllable parameters like ionomer EW or molecular weight cannot be directly determined. An alternative, bottom up approach is to build copolymers with different EW distribution, molecular weight and blockiness (nonrandomness of monomer arrangement) and computationally model properties such as morphology and the diffusion of protons and water molecules. Jang et al., for example, in the work referred to earlier10 studied the effect of blockiness on the properties of Nafion membranes. Blocky structures with high degree of clustering of the same monomers exhibited larger hydrophilic domains and better proton conductivity. The correspondence between hydrophilic domain size and efficient proton transport is evident even when comparisons are made across different families of polymers. On the basis of a simple cubic lattice model, Kreuer22 derived morphologies of Nafion and sulfonated polyetherketone membranes from small-angle X-ray scattering (SAXS) experiments. This model suggests that Nafion (which has higher proton diffusivities at low to intermediate hydration levels) has a morphology characterized with wide channels and smaller separation between the sulfonic acid groups. These findings from explicit MD simulations and experiment reaffirm the conclusions drawn from our previous statistical mechanic treatment on the role that morphology plays on proton diffusion. When the objective is to screen a large number of polymer structures the two approaches may be combined. The bottom up approach can be terminated after the hydrated membranes have been equilibrated and before a lengthy production run is carried out. Using the available morphological criteria for pore radius, distribution, and protrusion of the side-chains, one can eliminate polymers that would lead to lower diffusion results and simulate the proton transport only for the most promising cases. Another way to combine the two approaches and expedite the simulations was presented by Spohr et al.23-25 in which the polymer system is divided into structureless pore walls, modeled by a simple volume exclusion potential and explicit side-chain fragments, treated with an all atom potential. Excluding the motion of a large part of the polymer has allowed simulation times of tens of nanoseconds revealing slow conformational changes in the side-chain dihedral angles. The present work bears of these findings and is discussed later. The morphology and properties of real PEMs are also dependent on the method used in the processing of the polymer. For example, water uptake behavior is different for extruded versus cast ionomers.26 Extruded polymers also show anisotropic mechanical properties,27 which suggests similar anisotropy will be observed in all other morphology dependent properties.28 The complexity and time scale of the membrane preparation (includ-

Short-Side-Chain Perfluorosulfonic Acid Ionomer

J. Phys. Chem. B, Vol. 112, No. 10, 2008 2939

Figure 2. Labels for the various atom types used in the SSC force field. Figure 1. Structure of the repeat unit of the SSC ionomer, where in our simulations n ) 1 and m ) 7.

ing hydration) makes it impossible to reliably simulate its effect on the structure. When there is experimental information available, such as SAXS spectra, one can verify10 if the simulated system reflects the SAXS characteristics of the real one. Furthermore, using fairly well-developed techniques29,30 of structure refinement it is possible to bring an initial morphology into accord with the measured SAXS spectrum. The nature of these hydrated ionomer systems poses an immense challenge for the prediction of their three dimensional (3D) structure, perhaps even eclipsing that of proteins. Determining the structure from X-ray measurements is not possible as only a small part of the polymer is crystalline. Homology modeling that aids in the simulation of the structure of one protein based on the X-ray structures of other evolutionary related proteins is not applicable with ionomers. There are no statistical measures like the accessible surface area of a protein that could help rule out implausible ionomer morphologies. Furthermore, the available force fields are either too general or only limited to Nafion. Hence, in the present work we first develop a unique force field for the SSC ionomer (described in Part I of Simulation Methods) and an algorithm for building the 3D structure of ionomers in Part II. Dihedral angles play a crucial role in determining the morphology of a polymer. While the bond lengths and angles of a specific group are restricted about a single equilibrium value, dihedral potentials are modeled with periodic functions allowing for a number of low-energy conformations. By specifying the sequence of dihedral angles one can uniquely determine the morphology of the polymer, an approach widely used in proteins.31 Recent ab initio modeling of the SSC ionomer has suggested that the conformation of the backbone has an important effect on chemical properties under low hydration conditions. Extensive electronic structure calculations of oligomeric fragments possessing two side chains revealed that the extent of the separation of the sulfonic acid groups along with the conformation of the backbone significantly affect the propensity of the acid groups to dissociate.32,33 Kinked conformations of the backbone give rise to closer proximity of the sulfonic acid groups, stronger binding of the water molecules, and enhanced proton dissociation at lower degrees of hydration. As these properties needed to be correctly reproduced in our SSC-specific force field, we take special care in our treatment of the torsion potentials (described below). The methodology we have developed is applied to studying the structural and transport properties of the SSC polymer with an EW of 580 whose generic structure is indicated in Figure 1.

TABLE 1: SSC Force Field Parameters for the Harmonic Stretching Potential, Ebonda atom 1

atom 2

r0 [nm]

k [kJ mol-1 nm-2]b

H1 H2 H3 O3H O4:c S6: C2S: C: C:O C1:O C1:

O1 O2 O3H S6 S6: C2S: C2O F: O3 C1: C1:

0.1000 0.0982 0.0974 0.1628 0.1451 0.1886 0.1555 0.1348 0.1392 0.1566 0.1564

209200 454364 292880 292880 585760 292880 292880 253240 292880 292880 179627

a Ebond ) 0.5k(θ - θ0).2 b Force constants k are taken from ref 10. A colon (:) in the atom name represents a wildcard (e.g., C: covers C1, C1T, C1O, C2O, C2S, C2SI, etc.).

c

The polymeric system is taken to be comprised of a single macromolecule consisting of forty monomer units and the properties investigated at three distinct hydration levels: λ ) 3, 6, and 13 (where λ has its conventional definition as the number of water molecules per sulfonic acid group). We limit the scope of our simulations in this initial study to the structural characterization of the system and an investigation of vehicular proton diffusion (i.e., H3O+). In our subsequent communications, we will address more thoroughly the water diffusion (including the electroosmotic drag) and the proton transport through proton hopping or shuttling. The latter part requires an entirely new EVB parametrization, because under hydration conditions the sulfonic acid group will need to be represented in the EVB model. This methodology is currently under development. Computational Methods I. SSC-Specific Force Field. The force field we present here is based on the generic DREIDING force field developed by Goddard et al.34 and supplemented by several authors investigating the Nafion system.10,13,14,17 Our adaptation of the DREIDING force field to the SSC polymer is based on extensive electronic structure calculations performed by Paddison and Elliott.35 Six different torsion profiles were computed in their studies of a two-side chain oligomeric fragment altogether comprising more than 200 geometries. The average bond lengths and angles were extracted from these structures, and their values were used as the equilibrium parameters r0 and θ0 in the force field potentials. The force constants of the stretching and bending motion were assumed to be the same as in Nafion (see Figure 2 and Tables 1 and 2). The atom charges were calculated through a Mulliken population analysis of the ab initio structures of SSC fragment,

2940 J. Phys. Chem. B, Vol. 112, No. 10, 2008

Hristov et al.

TABLE 2: SSC Force Field Parameters for the Harmonic Bending Potential, Eanglea atom 1

atom 2

atom 3

θ0 [deg]

k [kJ mol-1 rad-2]b

H1 H2 C:c F1O O3 O3 C: F: F2O C2S: C2O F2S O4: O4 C2S: C2S: C1O

O1 O2 C: C1O C1O C1O C: C: C2O C2O C2S: C2S: S6: S6 S6 S6: O3

H1 H2 C: C1: C1: F1O: F: F: O3 O3 S6: S6: O4: O3H O3H O4: C2O

109.47 113.40 113.64 106.43 112.23 105.60 108.62 109.02 111.82 106.96 114.52 106.76 123.71 108.07 100.36 107.09 126.40

502.08 330.64 471.45 472.04 470.70 470.70 472.04 470.70 470.70 470.70 490.91 478.15 509.36 509.36 456.31 456.31 470.70

a Eangle ) 0.5k(θ - θ0).2 b Force constants k are adapted from refs 10 and 23. c A colon (:) in the atom name represents a wildcard (e.g., C: covers C1, C1T, C1O, C2O, C2S, C2SI).

TABLE 3: SSC Force Field Parameters for the Nonbonding Interactions, ECoulomb and ELJa atom 1

q

[kJ mol-1]b

σ [nm]c

H1 O1 H2 O2 F1 C1O F1O C2S C2SI H3 F2O F2S O3H S6 S6I C1 O4 O4I O3

0.4100 -0.8200 0.4606 -0.3818 -0.2709 0.4462 -0.2741 0.3715 0.3234 0.3882 -0.2569 -0.2414 -0.5140 1.2412 1.0237 0.5497 -0.4512 -0.5876 -0.5392

0.0418400 0.7732032 0.0418400 0.7732032 0.2075264 0.3978984 0.2075264 0.3978984 0.3978984 0.0004184 0.2075264 0.2075264 0.4004088 1.4392956 1.4392956 0.3531296 0.4004088 0.4004088 0.4004088

0.08018 0.31655 0.08018 0.31655 0.30249 0.34730 0.30249 0.34730 0.34730 0.28464 0.30249 0.30249 0.30332 0.35903 0.35903 0.34599 0.30332 0.30332 0.30332

ECoulomb ) k q1q2/r; k ) 138.93547 nm kJ mol-1 and ELJ ) 4 x12[((σ1 + σ2)/2r)12 - ((σ1 + σ2)/2r)6]. b,c The Lennard-Jones parameters and σ are obtained from ref 10. a

while the Lennard-Jones parameters and σ were assumed to remain unchanged from Nafion (Table 3). The last component of the SSC force field is the torsion potential, and it is the corner stone of our force field because it ties together all the other potentials and ensures agreement with the ab initio torsion barriers. From Paddison and Elliott’s oligomer geometries, one can further extract any structural information like dihedral angles, including those whose profile was not explicitly studied. One such dihedral is the O4:sS6: sC2S:sC2O (where the colon indicates a wildcard ensuring that both the neutral and ionized side-chains are covered, see Figure 2). Obtaining a molecular mechanics (MM) representation of the profile of this dihedral angle (assuming that it does not change after proton dissociation) makes it possible to go beyond torsions of the neutral system modeled in this ab initio study. Another compelling reason for including this torsion potential in the MM energy balance is the accompanying exclusion of the 1-4 nonbonding interactions (i.e., between the

TABLE 4: SSC Force Field Parameters for the Torsion Potential, Ediha atom 1 atom 2 atom 3 atom 4 H3 O3H S6: C2S: C2O C1: O4:

O3H S6 C2S: C2O O3 C1: S6:

S6 C2S: C2O O3 C1O C1: C2S:

C2S:b C2O O3 C1O C1: C1: C2O

a [kJ mol-1] 7.4185 -9.5835 -34.0433 8.4275 38.2339EVB 5.6356 12.0085

b

c [deg]

1.3349 -21.037 1.1453 17.673 1.0306 -0.931 1.1671 14.264 0.9381 32.494 1.0258 -49.553 0.4842 -6.642

a Edih ) a cos(bφ - c). b A colon (:) in the atom name represents a wildcard (e.g., C: covers C1, C1T, C1O, C2O, C2S, C2SI, etc.).

atoms O4: and C2O). The Coulomb interaction between these two atoms, especially in the case of an ionized chain, is quite large, which makes fitting the ab initio data very difficult. It is for this reason we have seven torsion potentials as shown in Table 4. It should be noted that our force field does not include all possible dihedral angles in the polymer. One such group of dihedral angles that are absent are those involving fluorine atoms. The reason we limit the number of dihedral angle types in the force field is to minimize the computational work during the MD simulations. As our force field (including the torsion potential) completely replaces the potentials defined in the parent force field, any terms from the latter that do not carry over to the SSC force field should be considered zero. This implies, however, that such terms are not simply missing but that their effect has been absorbed into the parameters of the remaining terms. Our force field should be considered complete in the sense that we have obtained the best possible correspondence to the ab initio calculations given the number of fitting parameters in the torsion potential. We now describe exactly how the values of these parameters were obtained. The torsion energy of a oligomer fragment with atom coordinates Ri is defined as the difference of the molecule’s ab initio energy and the sum of the energies of the other bonding and nonbonding potentials:

Etors(Ri) ) EQM - Ebonds - Eangles - ECoulomb - ELJ (1) The index i runs over all of the 200 SSC oligomer structures, e.g., R1 represents the coordinates of the atoms of the SSC structure in the first data point of the first torsion profile. This assignment bestows an additional role to the torsion term in enabling it to absorb the error of the model, which arises from the assumption that an ab initio energy can be split into a sum of MM terms. Alternatively, the torsion energy in eq 1 can be viewed as the residual energy after all other classical terms have been accounted. The torsion energy is modeled as a sum of dihedral potentials according to 7

Efit(Ri) )

Nm

i am cos(bmφm,l - c m) ∑ ∑ m)1 l)1

(2)

where the outer sum is over the seven distinct types of dihedral angles and the inner sum is over all instances of the given type in the oligomer. We choose to optimize the set of parameters i,j {a,b,c} so that the difference ∆E fit-tors between any two oligomer structures i and j is zero:

E ifit - E jfit - (E itors - E jtors) ) i,j i,j ∆E i,j fit - ∆E tors ≡ ∆E fit-tors (3)

In the ideal case, where the ab initio torsion profiles match those i,j derived from our force field, the ∆E fit-tors must vanish. Hence,



Figure 3. Comparison of the ab initio (B3LYP/6-311G**) and classical torsional profiles including contributions from all MM terms: (a) around the C1O-O3 bond and (b) around the C2S-S6 bond.

i,j the absolute value (or the square) of ∆E fit-tors is a measure of the quality of the fit for oligomer structures i and j. This allows us to define a penalty function for the entire set as

Structures

P(a,b,c) )

∑ i>j

i,j (∆E fit-tors )2

(4)

As we do not limit i and j to belong to data points from the same torsion profile, our fit will ensure that the ab initio energy differences are evenly matched across all the profiles.

The form of our dihedral potential, Edih ) a cos(bφ - c), warrants further description. The constant term that is present in other definitions of a dihedral potential36 is redundant in our case owing to the fact that we are only fitting to energy differences between molecules possessing the same type and number of dihedral angles (i.e., different conformations of the same oligomer). Even if a constant term was included in the potential, it cannot improve the fit, as it will cancel in eq 4. Finally, we note that the absence of such a term has no effect on the forces.


Hristov et al.

Figure 4. A single cone dihedral angle φ (red) can be used to uniquely determine all cone atom positions (green) given a force field (for further computational details see ref 42).

When the penalty function is at a minimum its derivatives (hereafter referred to as “forces”) with respect to am, bm, and cm will be zero for all m ) {1, 2,... 7}. By deriving analytical expressions for these forces that act on the parameters, we can employ an efficient minimization scheme for P, for example, steepest descent (SD) minimization.37 The first derivative of the penalty function with respect to each of the unknown parameters is straight forwardly obtained

Fam ) -

∂P

Structures

)-2

∂am

∑ i>j

i,j ∆E fit-tors

∂∆E i,j fit

)

∂am

Structures

-2

∑ i>j

i,j ∆E fit-tors

(

∂E ifit

-

∂am

)

∂E jfit ∂am

Nm Structures

)-2

∑ ∑ l)1 i >j

i,j i ∆E fit-tors (cos(bmφm,l - c m) j - cm)) (5) cos(bmφm,l

Fbm ) -

∂P

)

∂bm

Nm Structures

2

∑ ∑ l)1 i>j

i,j i i ∆E fit-tors am(sin(bmφm,l - cm)φm,l j j - cm)φm,l ) (6) sin(bmφm,l

F cm ) -

∂P

)

∂cm Nm Structures

-2

∑ ∑ l)1 i>j

i,j i ∆E fit-tors am(sin(bmφm,l - c m) j - cm)) (7) sin(bmφm,l

Because an SD algorithm had already been developed in the context of our general MD code,38 it was expedient to use this algorithm for the present purposes. The actual implementation involved treating the 21 (i.e., 7 × 3) unknowns (am, bm, and cm) as the coordinates of seven dummy atoms in 3D space. The forces that act on these atoms were calculated by eqs 5-7 with the “energy” of the system being the penalty function P in eq 4. The lowest possible value for P along with the corresponding values of the dummy atom coordinates (i.e., the force field parameters) were obtained at the completion of the minimization run and the latter are collected in Table 4. The most difficult torsion profile to accurately match was that around the C1O-O3 bond,35 where the side chain attaches to the polymer backbone (see Figure 2). This profile has the highest rotational barrier of more than 38 kJ/mol and is highly

Figure 5. A spherical SSC PFSA polymer generated with the help of the restraining potential Erestrain ) k(r - r0)2 with k ) 1.0 × 5 kJ mol-1 nm-2 and r0 ) 1 nm. The restrain acts only on the backbone CF2 carbon atoms (C1 in our nomenclature). A second restrain with smaller r0 may be applied to the sulfur atoms to produce an inverted micelle morphology.

asymmetric. The latter suggests significant contributions arising from Coulomb and Lennard-Jones interactions resulting from the close proximity of the atoms and occurs when the O-C bond eclipses the C-F bond. A comparison of the ab initio results with our MM force field energies is shown in Figure 3a and indicates very good agreement in most regions of the energy surface. A second torsion profile for rotation about the C2SS6 bond is displayed in Figure 3b and likewise shows good correspondence between the electronic structure calculations and the MM results. II. Building a 3D Polymer Structure. The molecular weight distribution of the Nafion polymer has recently been reported in the literature39 showing a peak in the 105-106 Da range. Given a typical Nafion membrane with an EW of 1100 such a molecular weight translates into chain lengths varying between N ) 90-900 monomer units. As data for the SSC ionomer is not available, information obtained from Nafion can be used as a guideline. Hence, an SSC polymer of chain length N will have (n + m)N backbone dihedral angles. If it is assumed that each dihedral angle is found in only three low-energy conformations (e.g., trans, + gauche and - gauche) then the number of low-energy configurations of the shortest polymer works out to be 10344 (see ref 40). This enormous configurational phase space cannot be traversed by any computational means or even by the real polymer in the time that the universe has existed.41 Consequently, the polymer can explore only a small part of the entire phase space volume, limited to some neighborhood of the initial configuration. Hence, when a polymer builder does not emulate the real system’s conformational choices during construction, no subsequent propagation techniques could guarantee that the model system will be in the same part of the phase space as the real polymer. The current available simulation methodology for generating the 3D structure of a polymer places the emphasis on the “deep” equilibration of the system, that is, techniques like annealing are usually employed to accelerate the location of an equilibrium structure. However, unlike proteins in their native state, the hydrated polymers are characterized as nonequilibrium, “living systems”.23 In this case, the benefits of annealing may be overshadowed by throwing the simulation system into regions of phase space that are not accessible to the real polymer. This



Figure 6. Polymer morphology snapshots at the end of the production run (2 ns) for the three levels of hydration: (a) λ ) 3, (b) λ ) 6, and (c) λ ) 13. Each system shown represents a 3 Å deep cross section of a 3 × 3 supercell so that the continuity can be observed in the x and y directions. The three supercells are drawn to scale. Hydronium ion oxygen atoms (blue) and sulfur atoms (yellow) are emphasized. The other atoms seen are oxygen (red), hydrogen (white), carbon (light blue), and fluorine (green).

issue is even more pertinent in the cases where the annealing takes place with significant amounts of water already in the system. The approach that we have adopted aims at building and folding the growing polymer chain as it may occur during synthesis based on the following rules: 1. A newly added fragment to the growing end of the polymer (i.e., CF2, CF3, CFO, SO3, groups or H) will have time to find its lowest energy conformation before the next fragment addition takes place. 2. The part of the polymer already built is much larger than the new fragment, hence the former will not have time to relax its structure upon fragment addition. 3. The entire polymer will relax its structure in the case of steric clashes but only to the point to avoid the particular clash and allow for a new addition to take place. The polymer is represented as a subset of all the dihedral angles, such that there is only one dihedral angle (cone dihedral) around a given bond as shown in Figure 4.42 Now we require that the cone dihedral angle φ is chosen so that the cone energy Econe(φ) is minimized. The latter is defined as the energy of all interactions, bonding and nonbonding, of the cone atoms and the constructed polymer. Because the bond lengths and angles of the cone atoms are set at their respective equilibrium values by construction, there are no stretching and bending contributions to Econe. The only bonding contribution is due to the dihedral angle potentials (e.g., one-fourth of the C1T-C1C1O-C1 dihedral energy is attributed as a bonding energy of the C1 atom, and because our force field does not include the C1T-C1-C1O-O3 dihedral potential or any dihedral potentials that include fluorine atoms, the C1 dihedral energy is the only bonding contribution to Econe). The nonbonding component of Econe includes the usual Coulomb and Lennard-Jones energy terms between the cone atoms and the rest of the already constructed polymer. When the system is periodic, the nonbonding interactions are evaluated with the periodic boundary conditions taken into account. There are also interactions between the cone atoms in the central cell and their periodic images. One can use different approaches to minimize Econe(φ). Our algorithm generates the cone atom positions and energies for

Figure 7. Hydrogen bond chain from the snapshot in Figure 6a (λ ) 3) using the minimum image convention. The sulfonate oxygens have been left out. The numbers on the sulfur atoms represent the index of the side chain (i.e., 1-40). Coloring scheme as indicated in Figure 6.

all values of φ between 0-360° with a step of 5°. The lowest energy value will determine the conformation of the cone dihedral angle and the positions of the cone atoms. According to our second rule, the structure of the constructed polymer is kept frozen upon the addition of new fragments. However, when the lowest Econe is above some threshold value (e.g., when a steric clash has occurred) the only way to relax the structure is to optimize the entire polymer. The third rule requires that this optimization is only a short one (i.e., about 15-30 steps) as it is reasonable to expect that a new addition may take place as soon as Econe has been lowered. We control how often the polymer is optimized beyond the end fragment rotations by varying the energy threshold value. As this optimization alleviates the somewhat unphysical assumption in the second rule, we select a threshold energy of about 150 kJ/mol per cone atom, so that an optimization is triggered at relatively short intervals.


Hristov et al.

Figure 8. Polymer chain of Figure 6a (λ ) 3) without the minimum image convention. Water and hydronium ions are left out for clarity. The sulfur atoms from the hydrogen bond bridges in Figure 7 are emphasized.

Figure 9. S-S pair correlation plot (between the sulfur atoms of the ionized sulfonic acid group) for the three hydration levels.

A distinct advantage of our polymer builder is that the desired polymer can be built directly into the periodic simulation box. Even though a polymer chain may span many simulation cells it always experiences the interactions of the periodic images and folds accordingly at every fragment addition step. This is in stark contrast to the technique that is currently employed where a polymer chain is built in an infinite simulation box first and then put into a periodic box. The latter technique complicates the folding of the polymer enormously due to the sheer size of the configurational phase space. One solution that tries to alleviate the problem is to use multiple copies of one

much shorter polymer chain, which can be equilibrated a lot faster. However, the resulting increase in polymer chain mobility in this segmentation approach may affect the transport properties of both water and hydronium ions. Another unique feature of our polymer builder is that the morphology of the polymer can be designed a priori. We can explore the properties of realistic polymer systems with ideal morphologies that could not otherwise be obtained (e.g., cylindrical, inverted micelle, etc). This is achieved by introducing atom-specific restrain potentials, an example of which is shown in Figure 5 with a harmonic potential. The energy



Figure 10. An ion cage exhibiting very short S-S distances. The rest of the side chain atoms in the polymer have been omitted for clarity. The two lower hydronium ions are also fully coordinated with three sulfonate groups (only two are shown).

Figure 11. S-O pair correlation plot (between the ionized sulfonic acid sulfur atom and the hydronium ion oxygen atom) at each of the three distinct hydration levels.

contribution from these constraints is included in the expression for Econe(φ), and the restrain forces taken into account during optimization. III. Simulation Details. All atom simulations of a single polymer system (n ) 1, m ) 7, EW of 580 g/mol, see Figure 1) were carried out at three different hydration levels: λ ) 3, 6, and 13 water molecules per sulfonic acid group. The simulation procedure begins with building a single 40 unit random SSC polymer in a 3D periodic box, as described in the previous section. The size of the simulation box is determined by the desired density of the hydrated polymer systems: 1.67 g/cm3, a value taken from Nafion.10 The nonbonding interactions are calculated with the isotropic periodic sum method43 with the cut off of the Coulomb and Lennard-Jones interactions at

1.5 nm. All optimizations performed in this work (unless otherwise stated) consist of first carrying out a SD minimization,37 followed by an adopted basis Newton-Raphson minimization.44 After the polymer’s structure was optimized, it was introduced into a simulation box filled with flexible three-center (F3C)45 water molecules. The interaction energy between each water molecule and the polymer was then calculated, and all but the required 40λ water molecules removed, beginning with the H2O molecules possessing the highest energy (i.e., those that overlap with the polymer). This procedure results in three hydrated polymer systems with 120, 240, and 520 water molecules and a total of 1768, 2128, and 2968 atoms, respectively. The water phase is first optimized (SD only) keeping the polymer frozen, followed by optimization of the


Hristov et al.

Figure 12. Fraction of hydronium ions with a given number of sulfonate neighbors. A sulfonate group is considered a neighbor if its sulfur atom is within 4.3 Å from the oxygen atom of a hydronium ion. The numbers have been averaged from the production run trajectory (2 ns).

whole system (polymer and water). All sulfonic acid protons are then transferred to their second closest neighbor water molecule producing solvent separated hydronium ions and the sulfonate-terminated side chains. Following proton dissociation, the atoms in the side chain are renamed according to our scheme (see Figure 2), and the corresponding new potentials are used henceforth. The system is again optimized in stages: first the hydronium ions (SD only, keeping water and polymer frozen), then water and hydronium ions (SD only, keeping the polymer frozen), and finally the whole system. The optimized hydrated polymer systems are then equilibrated at 300 K for up to 450 ps using a canonical (NVT) simulation. The NVT run is stopped when the predicted change in Epot is less than 1% per ns. All data is collected from a 2 ns microcanonical (NVE) production run with data points saved every 200 steps. The time step in all MD simulations is 1 fs. The H3O+ diffusion coefficients were calculated from an average of five initial configurations R(0), taken from the same trajectory about 25 ps apart, using the Einstein relation: 〈(R(t) - R(0))2〉 ) 6Dt. A second averaging is done over the diffusion coefficients of the 40 hydronium ions. The average temperature from the NVE production run is 315 K. Results and Discussion Polymer Morphology. Snapshots of typical morphologies at the three different water contents (λ ) 3, 6 and 13) are shown in Figure 6. Cursory examination of Figure 6 shows the increasing distinctions in the distribution of the water and hydronium ions through the polymer as the water content is increased; an observation which is consistent with studies of Devanathan et al. in Nafion membranes over a similar hydration range.16 The water appears to only form isolated clusters or domains with the sulfonate groups and hydronium ions at the lowest hydration level (λ ) 3). However, connectivity between these domains begins to emerge at λ ) 6 and might be considered as the onset of channel formation. There is also the emergence of a separation of the hydronium ions from the sulfonate groups by several water molecules. At the highest

hydration level (λ ) 13), these channels appear to permeate the morphology of the system. The sulfonate groups and hydronium ions are found at much greater average separations from each other. However, contact ion pairs between one sulfonate group and two hydronium ions are still visible in this snapshot (Figure 6c). Careful examination of the central section of Figure 6a reveals the presence of a hydrogen bond chain that spans about onethird of the supercell. This particular chain is composed of the following species: (going clockwise) H2O-H3O-SO3-H3OH2O-SO3-H3O-H2O-SO3-H3O-SO3-H3O-H2O, which is shown in Figure 7. The side chain index numbers (on the sulfur atoms) indicate that in this particular snapshot the sulfonate groups are not from neighboring side chains but are two side chains apart in the case of the 1-4 hydronium/water bridge and five side chains apart in the case of the 31-37 hydronium ion bridge. One should note that the distances between the sulfur atoms shown in Figure 7 are measured using the minimum image convention. These distances are also used in the paircorrelation plots presented later. However, because the polymer spans several simulation cells the separation of the same sulfonate groups in the actual polymer is quite different. The single polymeric fragment of Figure 6a is shown in Figure 8 without the minimum image convention and with both the water molecules and hydronium ions left out for clarity. This rendition of the macromolecule reveals much greater separation between the sulfur atoms (e.g., 22.7 Å for 1-4 and 25.6 Å for 31-37). Furthermore, it is evident that the polymer backbone consists of a few almost perfectly straight segments. This preference in the conformation of the backbone with the CF2 groups forming a helical pitch is similar to that witnessed in the QM calculations of Paddison and Elliott for a two unit oligomer32 and a three unit oligomer.46 This result in our MM simulations is similar despite the polymer being placed in a periodic simulation box. The periodicity of the model probably allows for the formation of hydrogen bond “bridges” while avoiding any significant bending of the backbone. Similarly, the real polymer may



Figure 13. Initial configuration average of the mean square displacement of hydronium ions in the 2 ns production run.

organize itself into straight bands with different threads held together by intermolecular hydrogen bond bridges. Further quantitative information concerning molecular interactions in the hydrated SSC systems is obtained in examining the RDF. The interactions of the sulfonate groups are quantified through the sulfur-sulfur g(R) and are plotted in Figure 9 for all three hydration levels. These pair correlation plots are derived from trajectory snapshots taken every 200 steps during the entire production run. One striking feature in this data is the intense peak for λ ) 3 at approximately 4.5 Å indicating the vast majority of the sulfonate groups are in very close proximity to one another. This separation between the sulfur atoms is significantly smaller than the one found in the case of only a single bridging hydronium ion (see Figure 7) and is actually very close to the separation between a sulfur atom and a hydronium ion in a contact ion pair (see Figure 11 and the discussion below). It is also somewhat smaller than observed in the simulations of Devanathan et al.16 with 1134 EW Nafion where at a similar degree of hydration they observed a broader peak with a shoulder at 4.6 Å and a maximum at about 5.2 Å. Scanning the trajectories in this part of the plot revealed a close packing structure that is shown in a wire-frame rendition of the atoms in Figure 10 (similar color convention as before). The hydronium ions in these structures are found in ion cages with all three hydrogen atoms hydrogen-bonded to the oxygen atoms in three different sulfonate groups. There is no water to mitigate the strong ionic interactions in such clusters and this results in very small separations of the sulfur atoms. Because the hydronium ions in these cages are immobilized to a large extent, it is insightful to know what fraction of the hydronium ions exist in such a state and how this depends on the degree of hydration; further analysis is presented later. The intensity of the 4.5 Å peak greatly diminishes for λ ) 6 broadening into two peaks at around 5.0 and 7.5 Å with a probability of about 1. This may be viewed as being entirely determined by the number density of the sulfur atoms and suggests there is little interaction between the sulfonate groups at this water content. The peak completely disappears at the highest hydration level with the majority of the SO3- groups seemingly very well separated by over 10 Å, the average

separation of the tertiary carbons (i.e., -CFO-) when the backbone carbon atoms are in an anti conformation.32 Another type of contact ion pair between the sulfonate groups and the hydronium ions are those where the ionic cages are more labile, either due to the incorporation of water or to missing hydronium ions. This is the case of the hydrogen bond bridges in Figure 7. The observed S-S separation in this type of bridges gives rise to the broad peak near 7.5 Å that is visible in the pair-correlation plot for both λ ) 3 and λ ) 6. Not surprisingly, these peaks of the labile ionic cages also disappear in the highest hydration case. The ion cages release some of the trapped hydronium ions with an increase in hydration, similar to what was seen in the S-S pair correlation plot. Another effect is the increase in the spatial separation between the hydronium ion and the sulfonate groups once the cages break down. This can be seen from the S-O pair correlation plot in Figure 11. The plot for λ ) 3 exhibits a sharp peak for the contact ion pair just below 4 Å. A shoulder on the left side of the peak is also visible that can be attributed to the fully “sulfonated” hydronium ions in the ion cages. The typical S-O separation there is between 3.58 and 4.0 Å. Another broad peak for the more hydrated systems is visible at about 6 Å that corresponds to the solvent-separated species sulfonate group/water/hydronium ion. The intensity of the first peak goes down and becomes more Gaussian as the ion cages disappear with the increase in λ. Simultaneously, it becomes more likely to find the hydronium ion further away from the sulfonate groups. Now we turn our attention to further quantification of these observations. The mobility of the hydronium ions is of interest and is directly linked to the resilience of these ion cages. As we have seen above, any defects in the ion cage composition lead to much greater separation of the hydronium ions and the sulfonate groups making the hydronium ions more effective, or at least more mobile, as charge carriers in the hydrated membrane. In Figure 12, we show a histogram of the number of sulfonate groups present in the first solvation shell (4.3 Å) of a hydronium ion. When the S-O separation between the two species is within this cutoff the sulfonate group is considered

2948 J. Phys. Chem. B, Vol. 112, No. 10, 2008 a neighbor of the hydronium ion. The histogram is drawn in such a way so that a direct comparison can be made with Nafion.16 For the SSC polymer at λ ) 3, about one-fourth of the hydronium ions are trapped in ion cages with three sulfonate neighbors. Roughly, equal numbers of hydronium ions are found with one sulfonate group, while slightly more have two neighbors. This trend in the histogram indicates that as the hydration level increases there is a dramatic drop in the number of fully “sulfonated” hydronium ions (i.e., with 3 SO3neighbors) and a significant increase in the number of free hydronium ions (i.e., with zero neighbors). It is intriguing that the SSC curves for λ ) 3, 6, and 13 resemble the Nafion plots of Devanathan et al. but for a higher hydration level, namely, λ ) 5, 9, and 20. This suggests that the SSC ionomer may have a greater number of free hydronium ions than Nafion at the same water content. Another feature that emerges is that the correspondence between the two PFSA membranes changes with the hydration level. While at the lowest water content the improvement is ∆λ ) 2 (i.e., 3f5), at the highest λ the improvement is already ∆λ ) 7 (i.e., 13f20). These simulation results can be used to explain the differences found in the experimental conductivity plots for the two polymer systems.47 In particular, it was found that low EW SSC ionomers have superior conductivity to Nafion 1100 at the same hydration level and with increasing the water content this superiority becomes even more pronounced. Hydronium Ion Diffusion. The average of the MSD of the hydronium ions for the three water contents is shown in Figure 13. We point out the importance of the initial configuration averaging of the MSD to obtain accurate results. The slopes of each of the MSD lines directly gives 6DH3O+. We computed values for the diffusion coefficient of the hydronium ion of 2.84 × 10-7, 1.36 × 10-6, and 3.47 × 10-6 cm2/s for water contents corresponding to λ ) 3, 6, and 13, respectively. The agreement with the experimental results of Kreuer et al.48 is nearly quantitative at the lowest water content attesting to the fact that the dominating transport mechanism under minimal hydration is vehicular. The diffusion coefficients at λ ) 6 and λ ) 13 are much lower than those obtained from the conductivity measurements indicating the importance of a component in the mobility of hydrated protons due to structural diffusion. Because our force field is based on the one used by Jang et al.10 to study the hydronium ion diffusion in Nafion, it is remarkable to see how well the two polymer systems are differentiated in the simulations. The only hydration level investigated for the Nafion system was λ ) 15 where one finds an average hydronium diffusion coefficient about 2.5 times smaller than the one we found here for λ ) 13. This is in qualitative agreement with the conductivity difference between a SSC polymer of EW 800 and Nafion 1100.47,48 As the conductivity of the polymer samples is proportional to the total diffusion coefficient (and not just the vehicular component), it is important to know if the vehicular contribution changes significantly between the two polymer systems. The computed vehicular diffusion coefficient in the case of Nafion corresponds to 19.5% of the total measured diffusion, whereas in the case of the SSC polymer the computed one is about 20.4%. As the hydration level of the Nafion system is a bit higher it is reasonable to expect a higher contribution from the structural diffusion and consequently a smaller vehicular component. We can conclude that the two polymer systems exhibit the same partition between their vehicular and hopping modes of proton

Hristov et al. transport at least in the case of λ between 13 and 15 water molecules per sulfonic acid group. Conclusions By carefully developing a SSC PFSA-specific force field and a polymer builder that allows us to simulate a single long polymer chain, we have obtained new insights into the morphology of the PEM systems and their transport properties. The polymer backbone is able to assume straight configurations while at the same time forming strong hydrogen bond bridges. Ion cages with fully “sulfonated” hydronium ions are typical of the low hydration level and explain the previously unresolved peaks (in Nafion) in the S-S and S-O pair correlation plots. The vehicular proton diffusion coefficients obtained here are in excellent agreement with the experimental values at the lowest hydration level indicating that under these conditions the proton shuttling does not play a major role. At the highest hydration level, the vehicular component accounts for only one-fifth of the total diffusion, exactly as observed in Nafion. The trends seen in a structural parameter representing the number of sulfonate groups around a hydronium ion allow us to suggest a plausible explanation for the differences between the conductivity of Nafion and the SSC PFSA ionomer. Acknowledgment. I.H.H. thanks the Alberta Ingenuity Fund for a scholarship and NSERC for additional support, S.J.P. acknowledges the ARO for partial support under contract number W911NF-07-1-0085, and R.P. thanks NSERC for support. References and Notes (1) Elliott, J. A.; Paddison, S. J. Phys. Chem. Chem. Phys. 2007, 8, 2602. (2) Vishnyakov, A.; Neimark, A. V. J. Phys. Chem. B 2001, 105, 9586. (3) Seeliger, D.; Hartnig, C.; Spohr, E. Electrochim. Acta 2005, 50, 4234. (4) Urata, S.; Irisawa, J.; Takada, A.; Shinoda, W.; Tsuzuki, S.; Mikami, M. J. Phys. Chem. B 2005, 109, 4269. (5) Urata, S.; Irisawa, J.; Takada, A.; Shinoda, W.; Tsuzuki, S.; Mikami, M. J. Phys. Chem. B 2005, 109, 17274. (6) Paul, R.; Paddison, S. J. J. Chem. Phys. 2001, 115, 7762. (7) Paul, R.; Paddison, S. J. J. Phys. Chem. B 2004, 108, 13231. (8) Paddison, S. J.; Zawodzinski, T. A., Jr. Solid State Ionics 1998, 115, 333. (9) Paddison, S. J.; Pratt, L. R.; Zawodzinski, T. A., Jr. J. New Mater. Electrochem. Sys. 1999, 2, 183. (10) Jang, S. S.; Molinero, V.; Cagin, T.; Goddard, W. A. J. Phys. Chem. B 2004, 108, 3149. (11) Blake, N. P.; Petersen, M. K.; Voth, G. A.; Metiu, H. J. Phys. Chem. B 2005, 109, 24244. (12) Blake, N. P.; Mills, G.; Metiu, H. J. Phys. Chem. B 2007, 111, 2490. (13) Petersen, M. K.; Wang, F.; Blake, N. P.; Metiu, H.; Voth, G. A. J. Phys. Chem. B 2005, 109, 24244. (14) Petersen, M. K.; Voth, G. A. J. Phys. Chem. B 2006, 110, 18594. (15) Venkatnathan, A.; Devanathan, R.; Dupius, M. J. Phys. Chem. B 2007, 111, 7234. (16) Devanathan, R.; Venkatnathan, A.; Dupius, M. J. Phys. Chem. B 2007, 111, 8069. (17) Devanathan, R.; Venkatnathan, A.; Dupius, M. J. Phys. Chem. B 2007, 111, 13006. (18) Brandell, D.; Karo, J.; Liivat, A.; Thomas, J. O. J. Mol. Model. 2007, 13, 1039. (19) Paddison, S. J.; Paul, R.; Zawodzinksi, T. A. J. Electrochem. Soc. 2000, 147, 617. (20) Paddison, S. J.; Paul, R.; Zawodzinksi, T. A. J. Chem. Phys. 2001, 115, 7753. (21) Paul, R.; Paddison, S. J. J. Chem. Phys. 2005, 123, 224704. (22) Kreuer, K. D. J. Membr. Sci. 2001, 185, 29. (23) Spohr, E.; Commer, P.; Kornyshev, A. A. J. Phys. Chem. B 2002, 106, 10560. (24) Spohr, E. Mol. Simul. 2004, 30, 107.

Short-Side-Chain Perfluorosulfonic Acid Ionomer (25) Commer, P.; Hartnig, C.; Seeliger, D.; Spohr, E. Mol. Simul. 2004, 30, 755. (26) Yamamoto, Y.; Ferrari, M. C.; Baschetti, M. G.; De Angelis, M. G.; Sarti, G. C. Desalination 2006, 200, 636. (27) Arcella, V.; Ghielmi, A.; Tommasi, G. Ann. N.Y. Acad. Sci. 2003, 984, 226. (28) Mauritz, K. A.; Moore, R. B. Chem. ReV. 2004, 104, 4535. (29) Svergun, D. I. Biophys. J. 1999, 76, 2879 and 77, 2896. (30) Grishaev, A.; Wu, J.; Trewhella, J.; Bax, A. J. Am. Chem. Soc. 2005, 127, 16621. (31) Lesk, A. M. Introduction to Protein Science - Architecture, Function and Genomics; Oxford University Press: New York, 2004. (32) Paddison, S. J.; Elliott, J. A. J. Phys. Chem. A 2005, 109, 7583. (33) Paddison, S. J.; Elliott, J. A. Solid State Ionics 2006, 177, 2385. (34) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897. (35) Paddison, S. J.; Elliott, J. A. Phys. Chem. Chem. Phys. 2006, 8, 2193. (36) see for example http://en.wikipedia.org/wiki/AMBER#Functional_ form (accessed December 17, 2007). (37) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical recipes in Fortran 77: the art of scientific computing, 2nd edition; Cambridge University Press: New York, 1992; Vol 1. (38) http://www.pyweave.com (accessed December 17, 2007). (39) Curtin, D. E.; Lousenberg, R. D.; Henry, T. J.; Tangeman, P. C.; Tisack, M. E. J. Power Sources 2004, 131, 41. (40) We can see from Figure 1 that for a polar side-chain fragment there is one bond in the backbone around which we can define a dihedral angle. Similarly, for a non-polar fragment there is also one bond. Hence, the total number of dihedral angles in the backbone (involving only carbon atoms)

J. Phys. Chem. B, Vol. 112, No. 10, 2008 2949 is (n + m)N. If each of these dihedral angles has c possible conformations the total number of configurations for the polymer would be c(n + m)N. Taking c ) 3, n ) 1, m ) 7, and N ) 90, we get approximately 10344. (41) In protein folding, this is known as the Levinthal paradox; see, for example, Huang, K. Lectures on Statistical Physics and Protein Folding; World Scientific Publishing: River Edge, NJ, 2005. (42) The procedure that we have used to determine the cone atom positions is the following: The atom that causes the backbone to grow has the highest priority of all cone atoms, and for clarification it is referred to as the “alpha” atom. All other cone atoms are considered as branches having equal but lower priority than the alpha atom. We designate them as “beta” atoms. In Figure 4, the alpha atom is C1 and the beta atoms are F1O and O3. For a given dihedral angle φ, the position of the alpha atom is first ascertained. Its coordinates are required to satisfy the equilibrium conditions for both the C1O-C1 distance and the C1-C1O-C1 angle. Once the alpha atom has been fixed, the beta atom positions are determined by rotating the alpha atom through 120° to positions C1′ and C1′′, respectively. We then scale the C1O-C1′ bond to the equilibrium distance of C1O-F1O and scale the C1O-C1′′ bond to the equilibrium distance of C1O-O3. (43) Brooks, B. R.; Wu, X. J. Chem. Phys. 2005, 122, 044107. (44) Chu, J.-W.; Trout, B. L.; Brooks, B. R. J. Chem. Phys. 2003, 119, 12708. (45) Levitt, M.; Hirshberg, M.; Sharon, R.; Laidig, K. E.; Daggett, V. J. Phys. Chem. B 1997, 101, 5051. (46) Paddison, S. J.; Elliott, J. A. Solid State Ionics 2007, 178, 561. (47) Kreuer, K. D.; Paddison, S. J.; Spohr, E.; Schuster, M. Chem. ReV. 2004, 104, 4637. (48) Kreuer, K. D.; Schuster, M.; Obliers, B.; Diat, O.; Traub, U.; Fuchs, A.; Klock, U.; Paddison, S. J.; Maier, J. J. Power Sources 2008, doi: 10.1016/j.jpowsour.2007.11.011.

Molecular Modeling of Proton Transport in the Short-Side-Chain

Recommend Documents