Comparing selected morphological models of hydrated nafion using large scale molecular dynamics simulations

Experimental elucidation of the nanoscale structure of hydrated Nafion, the most popular polymer electrolyte or proton exchange membrane (PEM) to date, and its influence on macroscopic proton conductance is particularly challenging. While it is generally agreed that hydrated Nafion is organized into distinct hydrophilic domains or clusters within a hydrophobic matrix, the geometry and length scale of these domains continues to be debated. For example, at least half a dozen different domain shapes, ranging from spheres to cylinders, have been proposed based on experimental SAXS and SANS studies. Since the characteristic length scale of these domains is believed to be ~2 to 5 nm, very large molecular dynamics (MD) simulations are needed to accurately probe the structure and morphology of these domains, especially their connectivity and percolation phenomena at varying water content.

COMPARING SELECTED MORPHOLOGICAL MODELS OF HYDRATED NAFION USING LARGE SCALE MOLECULAR DYNAMICS SIMULATIONS by Craig K. Knox A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Chemical Engineering The University of Utah December 2010 Copyright © Craig K. Knox 2010 All Rights Reserved The Unive rs i ty of Utah Graduat e School STATEMENT OF DISSERTATION APPROVAL The dissertation of Craig K. Knox has been approved by the following supervisory committee members: Gregory A. Voth , Chair 2/5/2010 Date Approved Jules J. Magda , Member 8/20/2009 Date Approved Valeria Molinero , Member 8/20/2009 Date Approved Terry A. Ring , Member 8/20/2009 Date Approved Grant D. Smith , Member 8/20/2009 Date Approved and by JoAnn S. Lighty , Chair of the Department of Chemical Engineering and by Charles A. Wight, Dean of The Graduate School. ABSTRACT Experimental elucidation of the nanoscale structure of hydrated Nafion, the most popular polymer electrolyte or proton exchange membrane (PEM) to date, and its influence on macroscopic proton conductance is particularly challenging. While it is generally agreed that hydrated Nafion is organized into distinct hydrophilic domains or clusters within a hydrophobic matrix, the geometry and length scale of these domains continues to be debated. For example, at least half a dozen different domain shapes, ranging from spheres to cylinders, have been proposed based on experimental SAXS and SANS studies. Since the characteristic length scale of these domains is believed to be ~2 to 5 nm, very large molecular dynamics (MD) simulations are needed to accurately probe the structure and morphology of these domains, especially their connectivity and percolation phenomena at varying water content. Using classical, all-atom MD with explicit hydronium ions, simulations have been performed to study the first-ever hydrated Nafion systems that are large enough (~2 million atoms in a ~30 nm cell) to directly observe several hydrophilic domains at the molecular level. These systems consisted of six of the most significant and relevant morphological models of Nafion to-date: (1) the cluster-channel model of Gierke, (2) the parallel cylinder model of Schmidt-Rohr, (3) the local-order model of Dreyfus, (4) the lamellar model of Litt, (5) the rod network model of Kreuer, and (6) a ‘random' model, iv commonly used in previous simulations, that does not directly assume any particular geometry, distribution, or morphology. These simulations revealed fast intercluster bridge formation and network percolation in all of the models. Sulfonates were found inside these bridges and played a significant role in percolation. Sulfonates also strongly aggregated around and inside clusters. Cluster surfaces were analyzed to study the hydrophilic-hydrophobic interface. Interfacial area and cluster volume significantly increased during the simulations, suggesting the need for morphological model refinement and improvement. Radial distribution functions and structure factors were calculated. All nonrandom models exhibited the characteristic experimental scattering peak, underscoring the insensitivity of this measurement to hydrophilic domain structure and highlighting the need for future work to clearly distinguish morphological models of Nafion. I dedicate this thesis to my darling Lissa, who inspires so much goodness and encourages so wholeheartedly during even the most distressing times; my precious children, Anna, Beka, Emma, and Gabe, who are more exciting and fascinating than any scientific discovery; my wonderful parents, Clifford and De'ann, who give so much of themselves and cheer so lovingly; and my family and friends, who help so generously and love so freely. I love all of you so very much! Thanks for your support these many years! "Let us realize that the privilege to work is a gift, that the power to work is a blessing, that love of work is success" - Pres. David O. McKay TABLE OF CONTENTS ABSTRACT ....................................................................................................................... iii LIST OF FIGURES ........................................................................................................... ix ACKNOWLEDGMENTS .................................................................................................. x 1 INTRODUCTION ........................................................................................................ 1 1.1 General Background .............................................................................................. 2 1.2 Literature Review .................................................................................................. 4 1.3 Motivation and Justification .................................................................................. 9 1.4 Objectives ............................................................................................................ 12 2 CHOICE OF MORPHOLOGICAL MODELS .......................................................... 16 2.1 Background .......................................................................................................... 16 2.2 Random Model .................................................................................................... 18 2.3 Cylinder Model .................................................................................................... 18 2.4 Sphere-Rod Model ............................................................................................... 19 2.5 Sphere Model ....................................................................................................... 21 2.6 Slab Model ........................................................................................................... 21 2.7 Rod Model ........................................................................................................... 22 3 SYSTEM BUILDING AND PREPARATION .......................................................... 24 3.1 Background .......................................................................................................... 24 3.2 Growing Chain Algorithm ................................................................................... 25 3.3 Carving Out Water Shapes .................................................................................. 28 3.4 Filling Random Voids .......................................................................................... 30 3.5 Hydronium Placement ......................................................................................... 31 3.6 Polymer Details .................................................................................................... 32 3.7 System Size .......................................................................................................... 34 3.8 Performance ......................................................................................................... 34 4 SIMULATION DETAILS .......................................................................................... 36 4.1 Background .......................................................................................................... 36 viii 4.2 Forcefield ............................................................................................................. 36 4.3 Annealing ............................................................................................................. 37 4.4 Production Runs ................................................................................................... 40 4.5 Performance ......................................................................................................... 40 4.6 Validation ............................................................................................................. 41 5 CLUSTERS AND BRIDGES ..................................................................................... 43 5.1 Background .......................................................................................................... 43 5.2 Rapid Bridge Formation ...................................................................................... 45 5.3 Percolation ........................................................................................................... 47 5.4 Sulfonate Participation and Penetration ............................................................... 53 6 SULFONATES ........................................................................................................... 56 6.1 Background .......................................................................................................... 56 6.2 Aggregation ......................................................................................................... 57 6.3 Mutual Spacing and Proximity ............................................................................ 58 6.4 Ion Pairing ............................................................................................................ 61 7 INTERFACE............................................................................................................... 63 7.1 Background .......................................................................................................... 63 7.2 Deformation and Structure ................................................................................... 67 7.3 Charge Saturation and Coverage ......................................................................... 71 8 IONOMER PEAK ...................................................................................................... 78 8.1 Background .......................................................................................................... 78 8.2 Basis of Morphological Models ........................................................................... 80 8.3 Insensitivity to Structural Differences ................................................................. 81 9 CONCLUSIONS......................................................................................................... 85 9.1 Summary .............................................................................................................. 85 9.2 Significance ......................................................................................................... 89 9.3 Original Contributions ......................................................................................... 92 9.4 Future Work ......................................................................................................... 93 REFERENCES ................................................................................................................. 96 LIST OF FIGURES 1. Morphological models of Nafion .................................................................................. 17 2. Flow chart of chain growth algorithm .......................................................................... 26 3. Example nucleation and chain growth .......................................................................... 28 4. Diagrams of water cutout and random void filling procedures .................................... 29 5. Bridges, sulfonate penetration, and cluster changes ..................................................... 46 6. Representative bridge between clusters ........................................................................ 47 7. Cluster count ................................................................................................................. 48 8. Size of biggest cluster ................................................................................................... 50 9. Rapid sulfonate aggregation around clusters ................................................................ 57 10. Radial distribution function and structure factor plots ............................................... 59 11. Water iso-density surfaces of random model .............................................................. 68 12. Water density profile of random model ...................................................................... 68 13. Top and bottom plot total cluster volume ................................................................... 69 14. Interfacial area ............................................................................................................ 71 15. Top and bottom sulfonate surface (interface) coverage .............................................. 74 ACKNOWLEDGMENTS I sincerely thank Prof. Gregory Voth for this great opportunity and challenging project. I gratefully acknowledge helpful discussions with Prof. Voth and Dr. Matt Petersen (Sandia) as well as with many others: Dr. Jan Andzelm (Army Research Laboratory), Dr. Bryan Pivovar (National Renewable Energy Laboratory), Prof. Valeria Molinero, Dr. Ly Le (Molinero group), Prof. Seung Soon Jang (Georgia Tech), Prof. Grant Smith, Prof. Oleg Borodin (Smith group), Dr. Gregory Mills (UC Santa Barbara), Dr. Craig Gittleman (GM), Dr. Shyam Kocha (Nissan), Prof. Dean Wheeler (BYU), James Throckmorton (Drexel), Sam Liston (CHPC), and many members of the Voth group (Prof. Gary Ayton, Dr. Philip Blood, Prof. Michele Ceotto, Shulu Feng, Dr. Eric Heatwole, Dr. Brian Hopkins, Prof. Satoru Iuchi, Dr. Sergey Izvekov, Dr. Sven Jakobtorweihen, Dr. Wei Jiang, Dr. Vinod Krishna, Jim Lai, Dr. Luca Larini, Hui Li, Dr. Pu Liu, Dr. Lanyuan Lu, Dr. Edward Lyman, Dr. C. Mark Maupin, Prof. Will Noid, Prof. Francesco Paesani, Dr. Sterling Paramore, Prof. W. James Pfaendtner, Marissa Saunders, Prof. Qiang Shi, Richard Swenson, Prof. Ian Thorpe, Prof. Arun Venkatnathan, Prof. Yanting Wang, Dr. Zunjing Wang, Prof. Kim Wong, Dr. Yujie Wu, Jianqing Xu, Dr. Takefumi Yamashita, Prof. Tianying Yan, Dr. Yong Zhang). I also gratefully acknowledge funding from the Department of Energy (DOE) and the Army Research Office Multidisciplinary University Research Initiative (ARO-MURI) as well as generous computer allocations from the Department of Defense (DoD). 1 INTRODUCTION This doctoral dissertation research focuses on advancing the understanding of the structure of hydrated Nafion, which is an amorphous ionomer polyelectrolyte membrane consisting of copolymer chains of nonpolar (N) tetrafluoroethylene and polar (P) perfluorosulfonic vinyl ether comonomers and which is currently the most popular commercial proton exchange membrane (PEM) for fuel cell applications, by using classical atomistic molecular dynamics (MD) simulations and modeling to compare properties of some of the most promising and relevant morphological models, which have been largely based on small-angle X-ray scattering (SAXS) and small-angle neutron scattering (SANS) experimental studies. The aim for this model comparison is to help guide future experiments and simulations in the ongoing effort to probe the true molecular-level hydrophilic domain structure of Nafion, with the hope that this understanding may help advance future PEM fuel cell design and performance, especially with regard to proton conductance and fuel crossover through the membrane. This dissertation will first discuss some relevant literature studies as well as motivation and objectives for performing this research. Then, the expected significance of this research will be highlighted and the important details of the approach will be described, including the model building and simulations. Justifications and clarifications will be given where necessary. The major results will then be presented and discussed. 2 And finally, conclusions will be made, followed by suggestions for future work and a summary of original contributions. 1.1 Nafion is both a random co-polymer and a polyelectrolyte consisting of nonpolar (N) tetrafluoroethylene (TFE) monomers, [-CF2-CF2-], and polar (P) perfluorosulfonic vinyl ether (PSVE) monomers, [-CF2-*CF(O-CF2-*CF(CF3)-O-CF2-CF2-SO3H)-], where *C indicates a carbon atom that is a chiral center. The N monomer compromises the skeletal backbone of Nafion chains and imparts hydrophobicity. The P monomer, which contains a ‘pendant' sidechain with a sulfonic acid headgroup at the end, imparts hydrophilicity. Hydrated Nafion is well-known for its nanophase segregation behavior (hydrophilic reverse micelle-like clusters embedded in a hydrophobic matrix), resulting from its amphiphilic nature. To date, it is the most commonly used polymer electrolyte or proton exchange membrane (PEM) for fuel cell applications. General Background PEM fuel cell research is currently aimed at developing novel membranes with higher proton conductance at higher operating temperatures and lower hydrations than Nafion in order to facilitate the use of cheaper electrocatalysts (e.g., Ni instead of Pt) and to decrease catalyst poisoning. The design of these alternative membranes may greatly benefit from detailed knowledge of the mechanism for proton transport in Nafion at the molecular level, which is poorly understood. Before proton transport in Nafion can be understood, detailed knowledge of the molecular level structure and morphology will be needed. The structure of the hydrophilic domains in hydrated Nafion has been strongly debated for over 40 years and also remains poorly understood. 3 Proton transport is an important design consideration in PEM fuel cells because it strongly affects the current load (and hence the power) that the fuel cell may sustain during operation. Higher proton transport may result in higher current and hence higher power output. Fuel cells function similarly to batteries, except that fuel cells are open systems (fuel is input and product is output) while batteries are closed systems (batch operation with no input or output streams). For the fuel cell to operate, protons, typically generated from hydrogen gas at the anode, must transport through the PEM (electrolyte) to reach the cathode, where they typically combine with oxygen gas from air to form water. Since the PEM is an electrical insulator, electrons produced at the anode must flow out of the fuel cell to perform electrical work (power electrical devices) and then are consumed at the cathode. Proton transport through the membrane is strongly coupled with the external flow of electrons and vice versa. Molecular dynamics (MD) simulation is becoming an increasingly powerful computational method for studying nanoscale structure and dynamics (e.g., transport phenomena) in complex condensed-phase systems. Given a forcefield or potential that describes particle (e.g., atom and molecule) interactions, classical MD calculates the trajectory (motion) of each particle in a system over time by numerically integrating Newton's equations of motion: 2 2 1 2 ( , ,..., ) dt f r r r m d ri i n i      = (1) 4 where i f  is the total force acting on particle i, which is a function of the positions ( j r ) of all other particles (i ≠ j) in the system (containing n particles) and which contains the forcefield information (e.g., bonding, Van der Waals, and Coulombic terms, etc.); m is mass; t is time; arrows denote vectors; and the expression 2 2 dt d ri  represents the acceleration of particle i. This system of ordinary differential equations can be readily solved using standard finite difference numerical methods (e.g., velocity Verlet, leap frog, or Gear predictor-corrector integration, etc.), which have been discussed in detail elsewhere [1]. The trajectory data typically contain particle positions, velocities, forces, and energies at each step in time. Given these data, statistical mechanical relationships allow the estimation of essentially any structural, dynamic, or thermophysical property or behavior of any fluid or solid material, in principle. In practice, however, time and length-scale sampling issues often limit the accuracy and precision of this methodology. Because most of the essential physics and chemistry is contained within particle interactions, the parameterization and development of forcefields may quantitatively and even qualitatively affect the results. Also, reactions and other complex phenomena often require simulations based on higher levels of theory (e.g., ab initio or quantum MD). 1.2 The literature is replete with studies on or related to Nafion and similar PEMs. This field is vast and continues to grow rapidly. For brevity, this literature review will only focus on the major contributions and most notable works in the field. However, many other important and excellent studies exist in the literature. Literature Review 5 Mauritz and Moore [2] compiled a detailed review covering the experimental, mostly SAXS/SANS, and simulation/modeling progress towards elucidating the structure of hydrated Nafion. This review spans the major contributions to this field from the 1970s to the present and describes each of the significant morphological models of Nafion and the experimental scattering data from which they have been derived. This review also discusses the "morphological complexity of Nafion membranes and the limited structural information that can be extracted from the relatively diffuse scattering profiles" as well as the "complicating" issues, namely "the random chemical structure of the copolymer, the complexity of co-organized crystalline and ionic [hydrophilic] domains, vast morphological variations with solvent swelling, the relatively low degree of crystallinity [i.e., amorphous nature], and the diffuse, heterogeneous nature of the morphology that leads to a wide range of domain dimensions." And finally, the authors state, "The simple fact remains that this polymer yields very little characteristic detail in the dimensions probed by these [small-angle scattering and wide-angle X-ray diffraction] methods. Thus, the quest for a universally accepted model continues with a spirited debate in the current literature." Schmidt-Rohr and Chen [3] recently compared the predicted scattering profiles of several popular morphological models of Nafion in the literature using a novel continuum-level modeling algorithm [4], which involves a numerical Fourier transformation of the scattering density of a 3D lattice to calculate the scattering intensity, and then used this study to propose a new model: the parallel cylinder model, which model resulted in a predicted scattering profile that appears to best fit the "ionomer peak", the small-angle upturn ("matrix knee") primarily due to the crystalline domains, 6 and the overall power-law behavior of an experimental scattering profile. These authors also used qualitative arguments to explain how this new model may explain some of the observed transport behavior of real Nafion. However, the authors did not calculate or study any other properties or behavior in this comparison. Furthermore, some of the models (e.g., the "cluster-channel" or spheres with connecting rods model of Gierke, the "slab" model, and the "worm" model) also closely fit the experimental scattering curve and behavior, suggesting that scattering alone may not be able to adequately resolve the ongoing structural debate. In a pioneering study, Jang et al. [5] explored the effect of monomeric sequence on the segregation of nanoscale hydrophilic and hydrophobic domains of hydrated Nafion using atomistic MD simulation and a random Monte Carlo build approach (randomly growing polymer chains in an empty box and then placing water in the leftover void spaces). In this study, two extreme sequences were used: blocky and dispersed, both at a water content of λ = 15. To date, this is the only attempt to understand this effect. Although the simulation cell sizes (~8 nm cube length, or ~37,000 atoms) were too small to probe multiple hydrophilic domains, the results of this study suggest that blocky sequences lead to more nanophase segregation than dispersed sequences. Nanophase segregation was observed in both cases. Since the predicted density of the blocky and dispersed systems were within ~5% and 10%, respectively, of the experimental value (~1.75 g/cm3 at 300 K), this study also helped validate the use of the modified DREIDING forcefield of Jang et al. [6, 7] with the F3C water and explicit hydronium model [8] in hydrated Nafion MD simulations. Also, in both cases, all sulfonate head groups were in the water-phase, consistent with experiment. 7 Seeliger et al. [9] performed two kinds of MD simulations of hydrated Nafion: (1) a two-state empirical valence bond (EVB) model, to treat proton hopping, with a flexible representation of the polymer chains (DREIDING [7]; 0.25 fs time step) and (2) a rigid SPC/E water [10] with a rigid explicit hydronium model as well as a rigid-bond representation of the polymer chains (2.5 fs time step). In addition, two kinds of small (~4 nm cube length, or ~5,000 atoms) model systems were used: random and slab. No details of the random building process were provided. The slab-like boxes were constructed by first generating a slab-like void (using restraints) in an equilibrated Nafion melt (nonhydrated polymer chains) and then placing the solvent molecules (water and hydronium) in the void space. During this process, the sulfonate groups were restrained to lie on the surface of the slab. Snapshots of three of the rigid systems after 30 ns suggest the formation of large water clusters "connected by narrow bridges" and "filamentous [clusters] with irregularly shaped narrow cylindrical pores" at high (λ = 10) and low (λ = 5) water contents, respectively. The slab system "[evolved] into a more filamentous structure, although… remnants of the original structure [were] still visible," even after 30 ns, suggesting the very long simulation time required to equilibrate such a glassy system. Ether oxygen-ether oxygen and ether oxygen-sulfur separation distance distributions in the Nafion sidechains from the flexible model suggest significant torsional conformational flexibility with a small preference for trans over gauche for the former and conformational rigidity with a strong preference for trans (extended) for the latter. A difference in directional (parallel and perpendicular directions to initial slab surface) proton diffusion using the flexible model (with EVB) indicates anisotropy and/or incomplete equilibration. 8 Wescott et al. [11] performed novel mesoscale simulations of hydrated Nafion using a self-consistent mean field theory approach. Classical MD simulations were used to parameterize a coarse-grained (CG) representation of the polymer "beads" (each bead represented one monomer) that was then used to propagate a density distribution on a cubic lattice (64x64x64 nodes) by numerically solving Langevin equations. A wide range (λ = 2, 4, 8, 16) of water contents was studied at room temperature for 5,000 time steps (~150 μs) each. The initial lattice configuration of monomer beads was homogeneous. Lower water contents led to isolated, nearly spherical domains while higher water contents led to deformed elliptical and "barbell" shapes due to domain merging. Domain shapes and sizes were "generally consistent with some experimental observations." The authors concluded that, because of the grid resolution (~0.9 nm) used, the proposed subnanometer-sized bridging structures responsible for the low percolation threshold in Nafion "remain elusive" and "cannot be resolved by mesoscale simulations." Although most of the research in the Voth group has been and continues to be focused on biophysical MD simulations (e.g., studying proteins, lipid bilayers, BAR domains, Actin filaments, ribosomes, etc.), a few members of the group have recently studied some materials science related systems, including water interfaces, ionic liquids, mesoporous silica PEMs, Nafion, and Nafion-like phosphonic acid based PEMs. Much of this research has employed the group's own multistate empirical valence bond (MS-EVB) model [12-15] to accurately treat proton hopping in water inside these systems. This method has been extended by a former postdoc in the group, Dr. Feng Wang, to handle multiple excess protons using a self-consistent iterative (SCI) approach known as 9 SCI-MS-EVB [16], allowing detailed studies of proton transport in small (~4 nm cube length, or ~5,000 atoms with 40 excess protons) Nafion systems to be performed. Dr. Matt Petersen, a former grad student and postdoc in the group, performed these pioneering Nafion studies [17, 18]. These studies focused on the excess proton and sulfonate interactions within a single hydrophilic domain and found that the excess proton prefers the solvent separated ion pair position with respect to the sulfonate (i.e., one water molecule in between the sulfonate and the proton complex), proton hopping toward sulfonates is anticorrelated with vehicular (Brownian) diffusion away from them, the future position of an excess proton is strongly correlated to the past position of a sulfonate with a time correlation of ~0.5 ns, and the proton complex is amphiphilic (due to its lone pair of electrons) and thus prefers hydrophilic/hydrophobic interfaces, suggesting the possibility of interesting bridging structures between domains in Nafion and partly explaining its low percolation threshold. Current and future fuel cell research in the group aims to continue Petersen's work by exploring more experimentally accessible quantities such as infrared (IR) spectra and applying the SCI-MS-EVB methodology to more model systems of Nafion, perhaps small subsystems of my proposed model systems. Also, future research may attempt to interface the SCI-MS-EVB method with coarse-graining (CG) or a mesoscale field theory approach to study proton transport in Nafion at larger length-scales. 1.3 After nearly three decades of research, experimental elucidation of the nanoscale structure of hydrated Nafion and its influence on macroscopic proton conductance remains largely unsolved and extremely challenging. While it is generally agreed that Motivation and Justification 10 hydrated Nafion is organized into distinct hydrophilic domains within a surrounding hydrophobic matrix (i.e., reverse micelle-like structures), the geometry and length scale of these domains continues to be debated. For example, at least half a dozen different domain shapes, ranging from spheres to cylinders to slabs, have been proposed based primarily on experimental SAXS/SANS morphological model studies. Likewise, different domain spatial orientations, size distributions, and swelling behavior with changing water content and temperature have also been proposed. These models have been mostly based on the position and shape of the characteristic "ionomer" peak in the experimental scattering spectra, along with some qualitative arguments based on predicted proton and water transport behavior in the models. Spectroscopy and microscopy experimental methods (e.g., NMR, Mossbauer spectroscopy, FTIR, TEM, AFM, etc.) have also been performed to try to learn more about the structure of Nafion, but they have been limited by resolution. Although these studies have complimented the scattering experiments, they have not made significant advances to this endeavor because of the complicating issues mentioned earlier in the Literature Review. In short, experimental studies have been unable to clearly determine the structure of Nafion. Without further research, experimental and/or simulation based, it is not clear which of these morphological models best represent real hydrated Nafion, especially in a complex fuel cell environment. Clearly, SAXS/SANS data alone cannot sufficiently determine this; complementary experimental and/or simulation approaches must be combined with the existing scattering data to solve this puzzle. But which experiments/simulations should be chosen? To help answer this question, this study has made a detailed structural comparison of some of the most promising and relevant 11 morphological models to find ways in which experiments/simulations may clearly distinguish between these models, thus guiding the future identification of the most representative model. Only simulations/theory can make such a model comparison since experiments cannot test ‘hypothetical' structures. Furthermore, the complexity of the problem is beyond current analytical theoretical treatments, so numerical calculations (i.e., simulation) are the only viable solution. Among the many simulation techniques currently available (e.g., quantum, classical/atomistic, coarse-grained, continuum), only atomistic and coarse-grained (CG) MD simulations currently have the necessary balance of computational performance, accuracy, and resolution required for this molecular-level comparison. Quantum simulations would not be currently feasible because of the huge computational cost. Continuum simulations do not have the necessary resolution to accurately probe the relevant nanoscale phenomena without first resorting to some kind of atomistic/CG simulation interface to a field theory framework (e.g., SPAM). Between the two most promising methods that may work, this study selected atomistic MD because a CG study would most likely require preliminary atomistic simulations, regardless, to parameterize and validate an accurate CG model before performing the actual study. Fortunately, the time-consuming task of parameterizing and validating an atomistic model has already been performed [5-8]. The atomistic MD simulations of this work can be used to develop CG models for future simulation studies. Atomistic MD has higher accuracy and resolution than CG, although these benefits come with higher computational cost. However, the computational demands of these atomistic MD 12 simulations are currently feasible, and the necessary computational resources were readily available for this project. 1.4 The objectives of this research are: Objectives 1) To help guide the search for the true hydrophilic domain structure in hydrated Nafion a) To determine the significant structural properties and behavior of each chosen morphological model b) To find those properties that may help distinguish between models in order to identify the most representative model c) To propose future experiments/simulations that may help clearly identify the true hydrophilic domain structure or the best model 2) To search for possible ‘bridging' mechanisms or structures that may facilitate proton transfer between hydrophilic domains in the various models 3) To understand how Nafion sidechains interact with the hydrophilic domains The significance of this research follows: 1) It may help refine the morphological models to better represent hydrated Nafion by providing molecular-level insights that are currently inaccessible to experiment 2) It may help ‘bridge the gap' between experiments and simulations in this field 3) It may help guide experiments/simulations to determine the true hydrophilic domain structure 4) It may lead to the development of novel PEMs with superior performance and cost 13 5) Knowledge of this structure is essential for the understanding of proton transport in Nafion, which is essential for improving PEM fuel cell operation and performance The objectives of this research project are (1) to help guide the search for the true hydrophilic domain structure in hydrated Nafion, (2) to search for possible ‘bridging' mechanisms or structures that may facilitate proton transfer between hydrophilic domains in the various models, and (3) to understand how Nafion sidechains interact with the hydrophilic domains. The first objective, guiding the structural search, is the principal focus and central theme of this project. Although this is a very ambitious goal overall, this study adopts a pragmatic approach that attempts to utilize tools readily available to answer important questions, which are closely related to this goal and which can be answered at this time with these tools. The purpose of this approach is to advance the understanding of the structure of Nafion and come closer towards choosing the best model, not to make this choice or to directly probe the true structure of real Nafion at this time, which indeed seems unattainable without more information and experimental/simulation advances. In this context, this study may be viewed as foundational, whereupon future studies may be built and further developed. Determining the significant structural properties and behavior of each chosen morphological model accomplishes this objective. Little is known about how each model would actually behave or the similarities/differences that would be observed between them, aside from predicted scattering profiles and, in some cases, limited (mostly qualitative) transport properties. To address this deficiency, ‘idealized' systems that 14 resemble the hypothetical models as closely as possible have been simulated and the resulting structural properties from these simulations have been quantitatively determined. This information may be used to find those properties that may help distinguish between models. By finding those properties that are sensitive to the choice of model, future experimental and/or simulation studies may be proposed to measure them in real Nafion. These measurements can then be compared with the model calculations of this work to clearly determine the best model and the true hydrophilic domain structure. Although the other two objectives, looking for a bridging mechanism and understanding the sidechain interactions, are secondary, they may significantly advance this field and have been pursued in this work. Fortunately, the same simulations above can be used to study both of these issues as well. Searching for a bridging mechanism is mostly concerned with examining, both visually and quantitatively, the connecting water structures that may form in between domains over time in each model. Understanding how the sulfonic acid group (sulfonate) at the end of each Nafion sidechain interacts with the hydrophilic domains can be achieved by quantitatively studying the domain surfaces of each model where the sulfonates should preferentially aggregate. In summary, six relevant morphological models of hydrated Nafion have been selected, built, and simulated using classical atomistic MD to probe the significant structural properties and behavior of each model, especially distinguishable characteristics that may help aid future studies. These properties have been analyzed and compared and recommendations for future work have been made that may lead to finding ways of determining the best model or probing the true molecular-level structure of 15 Nafion. The most time-consuming part of this study has been building and simulating these models, which will be discussed below. 2 CHOICE OF MORPHOLOGICAL MODELS 2.1 The morphological models of Nafion are based primarily on SAXS and SANS scattering data. The location and shape of the signature ‘ionomer' peak in the spectra has been the principal tool for model development. This peak is believed to be due to intercluster correlations of the hydrophilic domains. Background These models propose structural details of the hydrophilic clusters, such as geometry and shape, size distribution, and spacing as well as how these details may change with water content. Some of the models attempt to explain observed transport properties, although detailed molecular simulations have not been performed. Models were created by assuming form factors based on proposed geometry and then by fitting resulting structure factors to experimentally measured scattering intensities. Scattering spectra suggest general cluster sizes and spacings of ~30 to 50 Å. Because of this intercluster spacing, large systems are needed to probe multiple clusters at once. Previous simulation studies have not been large enough to accomplish this at the molecular level. Each model system (Figure 1) has been simulated at a water content of λ ~ 7.4, where λ is the ratio of the number of water molecules per sulfonate; in addition to this water content, the random model system has also been simulated at λ = 15 for 17 Figure 1. Morphological models of Nafion shown as schematics. The model names used in this work are labeled with their respective abbreviations in parentheses. Dark represents water; polymer is not shown for clarity. comparison to other studies. These two water contents are equivalent to ~10 and 20 wt% (~20 and 40 vol%), typical operating conditions in real PEM fuel cells. Also, many of the detailed geometrical parameters of the selected morphological models have only been provided in the literature at these water contents. Each model chosen for this study is considered to be a plausible representation of the structure of real Nafion. The set of models has been selected to span the widely diverse range of nearly all proposed model shapes and geometries in the literature. The six selected morphological model systems of hydrated Nafion will now be individually described. 18 2.2 A ‘random' model with randomly placed water molecules in a box of Nafion polymer chains (basically a model of random water shapes and sizes in polymer) [3, 5] has been studied. This model has been used before in the literature to try to predict a priori the true structure of Nafion without making any assumptions [5]. Due to practical time-scale limitations, such approaches have not succeeded in determining the true structure of Nafion but have succeeded in better understanding intracluster phenomena. Although this model does not directly assume any particular cluster shapes or sizes, it does indirectly assume a random distribution of shapes and sizes. This generally leads to a random distribution of intercluster spacings. Such uncorrelated or weakly correlated behavior lacks a strong scattering peak. Hence, this model is not expected to have the ionomer peak in the scattering spectra due to the short timescales accessible to simulation. Much longer timescales, in theory, would change this model into an accurate depiction of true Nafion, but such timescales are not computationally feasible. Random Model 2.3 The parallel ‘cylinder' model of Schmidt-Rohr is the newest morphological model of Nafion to-date. It consists of an array of aligned rigid tubes of varying diameter and spacing filled with water and surrounded by polymer [3]. No connecting bridge structures have been proposed for this model to explain the observed percolation threshold in Nafion. It was proposed to best fit the ionomer peak and the power law behavior at lower and higher scattering angle, although other models in that study also reasonably fit those features. The model was also proposed to better explain the low Cylinder Model 19 freezing behavior of water in Nafion, although no studies have been performed to quantify or validate this argument. The mean center-to-center separation distance between cylinders and the mean cylinder radius were proposed to be ~38 Å and ~12 Å, respectively. These same values were also used in this study to closely mimic the proposed hypothetical model. In accordance with Schmidt-Rohr's cylinder design, this study used nonoverlapping hard cylinders with a thin polymer shell (~7 Å) coating the outside of each one so that the cylinders were not allowed to touch one another. Cylinders were randomly placed based on these criteria and using the same replacement and shrinking radius approach of Schmidt-Rohr as well as periodic boundary conditions. The axis of each cylinder was placed parallel to the z-direction to maintain perfect alignment. Cylinder radii were chosen from the same ‘slanted' distribution reported by Schmidt-Rohr. Cylinders (30) were placed in a 300 Å x 300 Å x 300 Å box using this approach, which resulted in a total cylinder volume of ~17.1% of the total volume of the box, corresponding to the target water content. 2.4 The cluster-channel (‘sphere-rod') model of Gierke is among the oldest and most well known morphological models of Nafion. It consists of water spheres with connecting water rods between them, all surrounded by polymer [2, 3, 19]. Typically, the spheres are arranged on a lattice and are approximated with a monodisperse size distribution. However, in this work, the spheres have been randomly placed (nonoverlapping) and given random sizes from a broad, slanted distribution, similar to the approach used for the cylinder model, in order to better mimic the intuitive disorder of Sphere-Rod Model 20 amorphous systems and to more closely fit the observed scattering spectra, especially the ionomer peak. Since this model includes proposed connecting rods of water in between some of the spheres, it may be considered a type of network model that attempts to explain percolation phenomena. The mean center-to-center separation distance between spheres and the mean sphere and connecting rod radii were ~47 Å, ~18 Å, and ~5 Å, respectively. The mean separation distance corresponds to the location of the ionomer peak in the scattering spectrum. The size of spheres and rods affects the water content. Varying these values would result in different water contents and/or different ionomer peak locations (different correlations between clusters). This study used nonoverlapping hard spheres with a thin polymer shell (~5.5 Å) coating the outside of each one so that the spheres were not allowed to touch one another. Spheres were randomly placed based on these criteria and using the same replacement and shrinking radius approach of the cylinder model as well as periodic boundary conditions. Spheres within the average separation distance (~47 Å) from one another were linked together with rods. Sphere radii were chosen from a similar ‘slanted' distribution to that from the cylinder model, which basically results from the increasing difficulty of placing additional spheres as more and more spheres reside in the box and which is overcome by occasionally shrinking the newly placed spheres a little to compensate for this effect. Rod radii were chosen from a normal distribution with a 1 Å standard deviation. Spheres (190) with connecting rods (72) between neighboring pairs were placed in a 300 Å x 300 Å x 300 Å box using this approach, which resulted in a total sphere volume of ~16.6% of the total volume of the box (the total rod volume was ~1%), corresponding to the target water content. 21 2.5 The local-order (or hard ‘sphere') model of Dreyfus consists of randomly placed water spheres (without connecting rods) surrounded by polymer [2]. Geometrically, this model is similar to the sphere-rod model (cluster-channel model of Gierke), except it does not assume connecting bridges between spheres, and thus it is not a network model. Sphere Model The mean center-to-center separation distance between spheres and the mean sphere radius were the same as those of the sphere-rod model. The same building approach (randomly placed nonoverlapping hard spheres with a thin polymer shell and with occasional shrinking sphere radius) was used as mentioned above, except that no connecting rods were used for this model. Spheres (190) were placed in a 300 Å x 300 Å x 300 Å box using this approach, which resulted in a total sphere volume of ~16.9% of the total volume of the box, which closely matches that of the sphere-rod model. On average, the spheres used in this model were slightly larger than the spheres of the sphere-rod model to closely match water content while maintaining the same number of spheres. Due to the random nature of the cluster placement algorithm, the location and distribution of spheres are both completely uncorrelated between the two models. 2.6 The lamellar (‘slab') model of Litt consists of alternating parallel slabs or slices, each one filled with water or polymer and sandwiched in between two slabs of the other [3]. Since it lacks bridges between slabs, it is not a network model. It is a model containing elongated structures in 2D, whereas the cylinder model contains 1D elongated structures. The other models do not contain elongated clusters. Slab Model 22 The mean center-to-center separation distance between water slabs and the mean water slab thickness were ~47 Å and ~8.5 Å, respectively. This corresponds to a mean polymer slab thickness of ~38.5 Å (polymer thickness = water separation - water thickness). Nonoverlapping slabs were used as in the other models, but the random placement approach was more straightforward because of the 1D nature of the slab locations (placed along the z-direction), so shrinking slabs and overlap checks were not needed, making the building algorithm simpler. Slab thicknesses were taken from a normal distribution with a standard deviation of 0.9 Å and 3.5 Å for the water and polymer slabs, respectively. These values correspond to ~10% standard deviation in thickness of either slab type. Slabs (6 water and 6 polymer) were placed in a 300 Å x 300 Å x 300 Å box using this approach, which resulted in a total slab volume of ~17% of the total volume of the box. 2.7 The network (‘rod') model of Kreuer consists of a randomly connected and entangled network or mesh of flexible water tubes or worms in polymer [3]. Geometrically, this model is similar to the sphere-rod model with rod diameters that are larger than the sphere diameters. Likewise, it is a network model and includes bridges (rods) by design. Rod Model Because of the geometric similarity between the rod and sphere-rod models, the same building algorithm was used for both. Of course, the building parameters were significantly different. For the rod model (think rod-sphere model), small spheres or nodes were used to make rod connections, resulting in a rounded-end rod representation instead of a flat-end one. The mean radius of the spheres was the same as that of the 23 rods, ~6.3 Å, which was ~25% larger than the connecting rod size in the sphere-rod model. The spheres or nodes in the rod model were ~1/3 the size of those in the sphere-rod model. The mean center-to-center separation distance between spheres in the rod model was the same as that in the sphere-rod model (~47 Å). The rod cutoff distance between spheres was extended to 65 Å, compared to the 47 Å cutoff of the sphere-rod model, to create a more connected network of rods and to increase the rod count an order of magnitude (to closely match the desired water content). Sphere and rod radii were both chosen from normal distributions with standard deviations of 1 Å and 2 Å, respectively, resulting in more mono-disperse spheres and a broader distribution of rod sizes. The distribution of sphere sizes in the rod model was not as "slanted" as that of the sphere-rod model because of the lower potential for sphere overlap (and hence less frequent need for periodically shrinking spheres) in the former. Spheres (198) with connecting rods (830) were placed in a 300 Å x 300 Å x 300 Å box, resulting in a fully connected (percolating) 3D network of rods. The average number of connecting rods per sphere was more than four (remember that rods are shared between pairs of spheres). 3 SYSTEM BUILDING AND PREPARATION 3.1 Building and preparing systems for MD simulation is often the most time-consuming task in terms of human effort and labor. One of the most difficult steps of this task generally involves carefully placing the atoms in a box to generate an initial configuration of coordinates. Many biological systems have known crystal structures that are readily available from online databases (e.g., PDB), making the generation of the initial configuration easier. Similarly, some materials science systems also have or are related to known crystal structures. However, amorphous systems, which are common in materials science, and systems with incomplete or unknown crystal structures, which are common in biology, require more work to generate initial configurations. Many methods exist for ‘guessing' the coordinates of atoms in macromolecules in these situations: the rotational isomeric state (RIS) method, which randomly chooses known dihedral states (e.g., trans, gauche+, gauche-) from a known distribution (weighed according to sum of Boltzmann factors, exp[ -βUstate ], where β = kbT and kb is the Boltzmann constant, T is absolute temperature, and Ustate is the potential energy of a specific dihedral state) and places atoms and/or grows chains bond-by-bond accordingly; melting a similar crystal structure; simulating the polymerization of a box of free (unreacted) monomers; randomly or systematically (using chemical intuition) arranging atoms in a box (off-lattice or on a lattice), which can also be used for chain growth; and other chain growing Background 25 variations on the RIS method. Additionally, these methods can be combined and extended by ab initio quantum calculations, minimization, molecular mechanics, MD simulation, and coarse-graining and other field theory approaches. In the case of coarse-graining, a coarse-grained (CG) representation of an atomistic system, which consists of grouping together parts of molecules to reduce the number of pair interactions, may be built and then randomly or systematically back-mapped into atomistic detail. 3.2 This study has adopted and modified the Monte Carlo (MC) chain growing algorithm of Theodorou and Suter [20, 21], which is related to the RIS method mentioned above. A Fortran 90/95 code was independently designed and developed specifically for this project. The random MC building algorithm has been used before in the literature to build small (~40,000 atoms) random-like Nafion systems at condensed-phase density for subsequent MD simulation [5]. The principal advantages of using this algorithm are its speed and its precise control of monomer sequence and chain length, all of which have been useful for this work. Growing Chain Algorithm The implementation of the algorithm (Figure 2) used in this work attempts to grow chains by placing atoms sequentially in a box under periodic boundary conditions, which is a trick for simulating bulk systems without edge effects, by choosing random dihedrals (torsion angles) and using other bonded information (bond lengths, bond angles, and coordinates of bonding partners) to determine the 3D Cartesian coordinates (x, y, z) of each grown atom inside the box, checking for possible overlap with the other atoms that have already been placed, and then trying again or continuing onward, accordingly. If the current atom overlaps another already placed atom, according to the 26 Figure 2. Flow chart of chain growth algorithm. van der Waal (VDW) radii of the atoms and the distance (nearest image convention) between them, then the placement attempt is rejected and another attempt is made; otherwise, the atom is placed in the box and placement attempts for the next atom are made. If many rejections occur for a particular atom then placement attempts for the previous atom start over as needed, and so on. In other words, forward and backward progress can result at times, depending on atom overlap. Reducing the system density and/or the VDW atomic radii can reduce the likelihood of atom overlap, thus reducing the required number of iterations and the amount of time needed to build a system. This algorithm is analogous to performing traditional MC using the standard Metropolis 27 acceptance/rejection criteria based on a hard-wall potential (U = ∞ if r < σ, U = 0 otherwise; where U is the potential energy, r is the distance between atoms, and σ is the location of the hard-wall, typically near the VDW radii). The amorphous builder code written for this study has modified this algorithm to build chains in parallel instead of one after the other in order to better represent real chain polymerization. This code can grow polymer chains around existing water shapes or it can place water in the empty voids afterwards, allowing for a wide range of model system designs, which will be discussed below. It also uses different VDW radii for polymer and solvent atoms for better atom overlap control, especially when building model systems that require the chains to grow around already placed water. This code also allows parts of chains to have fixed (known) Cartesian coordinates. The rest of the chain can be randomly grown given these fixed coordinates. This may be useful for implementing coarse-graining and back-mapping to atomistic detail. Also, this may be used for predicting the structures of macromolecules with incomplete crystallographic data. Multiple polymer chains are first ‘nucleated' by randomly placing them in a box and then chain growth commences on all of the chains in parallel (Figure 3). While atom placement attempts are being made, including going back and redoing previous atom placements because of too many unsuccessful attempts further down the chain, attention remains fixed on the current chain until the original progress point (the furthest atom down that chain for which placement attempts have been made, which is also the point of backward progression) has been successfully passed, at which time attention is then turned to the next chain. After each successful atom addition to a chain (during forward 28 Figure 3. Example nucleation and chain growth. From left to right, chain seeds are nucleated in a box and then chains are grown, bond-by-bond. Periodic boundary conditions are not shown to demonstrate the length of the chains. progress), chain growth continues with the next chain, and so on until all of the chains have been completely grown. This method of alternating attention on each chain in turn evenly balances the increasing difficulty of chain growth as the chains get longer and better mimics real polymerization. 3.3 All of the model systems, with the exception of the random model, have been built by first placing water into an empty box and then growing polymer around the water using the MC algorithm and code mentioned above ( Carving Out Water Shapes Figure 4). The polymer-water separation gap can be carefully controlled by adjusting the VDW radii of the water and polymer. The water was placed by first carving out representative model shapes from a large box of water and then removing the excess water (outside of these shapes), leaving only the water inside the shapes. This allows for the construction of tailor-made water 29 Figure 4. Diagrams of (top) water cutout and (bottom) random void filling procedures. cutout shapes grow chains water empty grow chains around shapes fill voids with water 30 clusters that conform to any morphological model. Sizes, shapes, locations, and spacings of these clusters can be controlled and adjusted as needed. For some of the models with fewer numbers of clusters (e.g., slab and cylinder), reaching the target water content required a trial-and-error iteration approach because of the statistics of smaller samples (individual clusters). A Python code was designed and developed to randomly carve out different geometries from a water box before using the MC chain growth code. To closely mimic the proposed morphological models, the carver code can randomly pick sizes and spacings from distributions, allow/disallow cluster overlap (e.g., the cylinder model was designed and built with nonoverlapping cylinders with a thin polymer shell surrounding each cylinder, as originally proposed in the literature), and connect neighboring clusters with rods (e.g., the sphere-rod and rod models used this feature). The water clusters in the cylinder, sphere-rod, sphere, rod, and slab models were built using this code, followed by chain growth. 3.4 The ‘random' model was built using the reverse approach: first growing Nafion polymer chains in an empty box, using the MC algorithm discussed above, and then placing water into the randomly dispersed void spaces in between polymer ( Filling Random Voids Figure 4). The end result of this is similar to the regular approach of cutting out random shapes, sizes, and spacings of water clusters and then growing chains around those clusters, but the procedures are simpler for the reverse approach because it avoids the complexity of directly defining and creating random objects. 31 A Fortran 90/95 code was designed and developed for the purpose of placing the water into the random void spaces in between polymer. This code superimposes a box of water over the box containing the grown polymer chains and deletes those waters that overlap with polymer and adds the nonoverlapping waters into the polymer box, thus filling only the random void spaces with water. To reach specific water contents, the polymer-water separation gap can be tuned and/or waters can be randomly deleted from already filled voids. Also, the density can be adjusted during the chain growth stage beforehand to create bigger and more numerous voids, and multiple polymer boxes can be grown as desired. 3.5 After placing water and growing chains, in all model systems, each sulfonate's (SO3 -) nearest water molecule (H2O) was mutated into an explicit hydronium ion (H3O+) for charge neutralization. At the water contents of interest to this study, all sulfonates in the models occur in the de-protonated form, as suggested by pKa calculations [22] and FT-IR measurements [23] in Nafion. A Fortran 90/95 code was designed and developed to quickly find the closest water to each sulfonate, and to convert that water into a hydronium ion. Given the strong interaction between sulfonates and hydroniums and the fully de-protonated state of the sulfonates at these water contents, this approximation seems reasonable. Hydronium Placement Notice that the sulfonates were not constrained to be on the hydrophilic cluster surfaces even though they are believed to primarily aggregate at these interfaces. By not constraining the sulfonates, their aggregation behavior can be predicted and accurately studied. This will be discussed in more detail in further sections. 32 3.6 The equivalent weight (EW) of the Nafion chains in this work was 1,100 g/mol, which corresponds to ~6.5 N/P monomers and which is the most commonly used EW in PEM fuel cells. EW is defined as the molecular weight of the polymer per sulfonate. It is a measure of sulfonate content, by which the average spacing between sulfonates may be roughly inferred (assuming a particular monomer sequence and chain configuration). Lower EW (higher sulfonate content) membranes are becoming more popular for PEM fuel cells because of their higher proton conductance, despite their poor mechanical properties and the difficulty of processing them, all of which have been improved in recent years through the use of novel nanocomposites and fillers. Polymer Details The results of Jang et al. [5] suggest an influence from Nafion monomer sequences in Nafion on nanophase segregation and hydrophilic domain morphology, highlighting the need for realistic monomer sequencing in simulation studies. However, no other simulation study has ever utilized realistic monomer sequences. Since Nafion is a random co-polymer, each grown chain in this study has been given a unique, random monomer sequence (e.g., …NNNPNNNNNNNNNNNPPNNNNNNN…) based on the monomer reactivity ratios (N 8.0, P 0.08) and steady-state monomer concentrations (N 44 mol%, P 56 mol%) in a typical Nafion polymerization reactor [24]. The degree of randomness (DR) [5] obtained was ~1 (‘dispersed' sequences). Similarly, each chain in these model systems has been given a unique, random length based on an approx. lognormal experimental molecular weight (MW) distribution (mean = 5.0 or ~105 g/mol, std. deviation = 0.3) [25]. All previous simulation studies have used unrealistically short Nafion chains (e.g., 80 mers/chain, even though the 33 experimental MW suggests an average of at least ~800 mers/chain or ~100 P's per chain) for faster polymer equilibration and rearrangement as a way of trying to overcome the pitfalls of very limited simulation times. However, the effect of chain length on hydrophilic domain morphology is unknown. Intuitively, shorter chains may overemphasize end-of-chain interactions and underemphasize chain entanglement and intrachain interactions and behavior, thus possibly affecting the resulting structure and dynamics. Also, since formal equilibration of systems containing even the shortest of these simulated Nafion chains has never been demonstrated in the literature, nor currently seems computationally feasible with traditional approaches, using short chains does not seem to have any proven practical advantages over using full-length chains except to help shrink simulation system size, which will be discussed below. Further studies that compare large systems of short and long chains are needed to address this issue. Since the sidechain monomer (P) of Nafion has two chiral centers and the possibility of head/tail flipping, eight (2x2x2=8) different isomers of this monomer have been or will be incorporated into the model system chains using equally weighted probabilities, leading to atactic chains. In the absence of experimental data, the assumption of equal probabilities seems reasonable given the amorphous nature of Nafion and its random monomer sequences. Although the chirality and head/tail flipping of P monomers are not expected to significantly influence these simulation results, these minor details have been included to make this study as realistic as possible and did not require much effort or work. 34 CF3 terminated chain ends approximate the initiators used for polymerization of Nafion, which, after reaction, become the chain ends of real Nafion. This assumption is negligible for long chains. 3.7 Since the characteristic length scale of a single hydrophilic domain in Nafion is believed to be ~2 to 5 nm, which is approximately equal to the box size of a typical ~5,000 atom PEM system, very large (~20+ nm box length) MD systems are needed to accurately probe the structure and morphology of several of these nanoscale domains at once, especially their connectivity and percolation phenomena at varying water content. Using classical, all-atom MD with explicit hydronium ions, very large hydrated Nafion systems (~2 million atoms in a cube with a side length of ~30 nm) have been simulated to directly observe several hydrophilic domains at the molecular-level. This corresponds to 200 polymer chains and ~21,000 counter ion pairs. Each model system explicitly includes and accounts for every atom in every molecule, whether Nafion, water, or hydronium. Smaller systems would not be able to achieve the objectives of this study. System Size 3.8 The MC chain grower code mentioned above uses a unique and original grid-hashing technique that was independently developed for this project to significantly speed up the atom overlap check, thus making the entire build process very fast for even multimillion atom systems (e.g., using this grid-hashing technique, the code takes ~2 hours to build a ~2 million atom Nafion system at a density of ~1.6 g/cm3 with a VDW scale factor of 40% for polymer atoms and 100% for water atoms while running on one Performance 35 Intel Pentium 4 processor; for comparison, this would take ~2 weeks without the grid-hashing technique). The random void filler and hydronium mutater codes mentioned above also use this same grid library for fast distance and neighbor searches. 4 SIMULATION DETAILS 4.1 MD simulations have been performed on one representative configuration of each model system to further explore structural properties, in particular how they change with time from their initial ‘hypothetical' model state. Background 4.2 The all-atom modified DREIDING forcefield of Jang et al. [5-7] and the flexible F3C water and explicit hydronium ion model [8] have been used. The modified DREIDING and F3C forcefields have been chosen for a close comparison with previous simulation studies, for simplicity, and to balance computational expense with needed accuracy. In particular, a polarizable forcefield alternative would be too computationally expensive while rigid water and hydronium model [10, 26] alternatives would not drastically reduce the computational expense. Using a flexible water model also allows clear validation with previous simulation studies and sets the groundwork for future studies that may use dissociable proton models (e.g., MS-EVB), which generally require flexible water models. Forcefield In this case, a classical explicit hydronium model (nondissociable) is justified because this study probes relatively large-scale (nanoscale) water structure without focusing on small-scale or localized excess proton behavior (e.g., the local environment 37 surrounding a couple of sulfonates). In this regard, the objectives of this work are very different than those of most modelers in the field. Also, this research does not attempt to study proton conductance directly but rather focuses on the hydrophilic (especially water) domain structure itself, from which the effect on proton transfer may be inferred. For example, the vehicular (Brownian) component of hydronium diffusion may be calculated using this approach, but the important proton hopping (Grotthuss mechanism whereby protons move via the concerted formation/breakage of covalent bonds in surrounding water molecules) component of this diffusion may not be directly calculated, nor by extension the total diffusion. The latter would require a reactive forcefield model, e.g., the multistate empirical valence bond (MS-EVB) model, to accurately handle the atomic interactions in a changing bond topology. Although this approach has been used before in the literature and by some members of the Voth research group to study small (~5,000 atoms with 40 excess protons) Nafion systems [9, 17, 18], it is not currently computationally feasible to simulate systems as large as those of this work, each of which would require ~21,000 excess protons! Fortunately, however, the approach of this study should result in accurate global nanoscale structure since recent sulfonate-oxygen/ hydronium-oxygen radial distribution function (RDF) findings that compare results from an explicit hydronium model and the MS-EVB model show close agreement at a distance of 5 Å and beyond [18], suggesting that the MS-EVB model, and hence proton hopping, does not greatly influence long-range structure in Nafion. 4.3 Many novel and exotic performance enhancing and equilibration techniques exist for polymer simulation such as multiple timestepping (calculating the expensive long- Annealing 38 range electrostatic interactions less frequently than short-range interactions), short-range electrostatic approximations (mapping long-range electrostatics onto short-range potentials), coarse-grain and field theory approaches, resolution exchange (coarse-graining and back-mapping), double-bridging hybrid MC/MD (changing chain topologies during simulations to accelerate chain equilibration) and MC ‘cluster' (group of atoms) moves, phantom chains then slow push-off (gradually ramping up nonbonded interactions over time to allow chains to initially slide past one another without excluded volume effects), and 4th dimensional ‘ghost' particle methods (adding another dimension of space to reduce excluded volume effects and enhance chain motion). This study has adopted a simple and efficient simulated annealing technique that heats up the polymer to enhance chain motion while holding the water and hydronium fixed in space to preserve the initial model structures. Each representative model system was annealed with NVT (constant volume and temperature) MD at 600 K (above the glass transition temperature of Nafion) for at least ~0.5 ns before starting ‘production' NPT (constant temperature and pressure) MD simulations. This way, the polymer chains relax and rearrange, especially the sulfonates of the sidechains, without disrupting the assumed solvent model structures. However, this will not formally equilibrate the chain backbones, which are long, highly entangled, and glassy and hence would require prohibitively long simulation times to completely equilibrate. Fortunately, this global equilibration should not be necessary because the purpose of these simulations is to study the behavior of assumed model structures as opposed to predicting a priori the true structure(s). Furthermore, the MC building algorithm should yield satisfactory random backbone configurations for an amorphous polymer such as 39 Nafion. And finally, this work represents state-of-the-art research in this field. This issue is beyond the realm of current computational capabilities using well-established, traditional MD methods. Using more advanced and sophisticated equilibration tricks mentioned above (e.g., ghost particles, 4th dimensional equations of motion, double-bridging MC, etc.) would be beyond the scope of this project and would require extensive coding and more research, and none of these have ever been applied in the study of a polymer as complex as Nafion or of systems as large as those of this study. Whether or not the backbones do significantly influence the sulfonate and water morphology is unknown. The hypothesis of this work is that they do not because of the amorphous nature of Nafion, the mobility of the sulfonates, the high density of sulfonates, and the strong sulfonate-solvent interactions. Hopefully future studies will be able to explore this assumption further to test the validity of this hypothesis. For now, local equilibration and relaxation of the mobile sidechains in an essentially fixed ‘scaffolding' of very slow backbones should suffice and has already been observed during short simulated annealing runs (significant sulfonate rearrangement and polymer-water interface coverage were observed, in agreement with most morphological models). The two λ = 7.4 random models came from two λ = 15 random models after ~20 ns NPT simulation and after random water deletion. One of the random models had been annealed and the other had not. The annealing consisted of heating the system from 300 K to 600 K over 0.1 ns with a Berendsen barostat (1 atm, 4.57x10-5 bar-1 compressibility, and 100 fs relaxation time) and temperature reassignment (from Boltzmann distribution) every 100 timesteps, running hot (600 K) for 1 ns at constant volume (no barostat) with same frequency of temperature reassignment, and cooling from 600 K to 300 K over 0.5 40 ns using the same barostat and temperature reassignment as for the heating stage. During the heating, annealing, and cooling stages, no restraints or constraints (no fixed atoms) were used. This differs from the annealing protocol of the other models in that they were built at the lower water content (λ = 7.4) and were annealed in the following manner: running hot at 600 K for 0.5 ns using velocity rescaling every 100 timesteps and while holding the waters and hydronium ions fixed in space in order to preserve the assumed morphology. The cylinder model was also annealed for 6 ns (~10x longer) in order to assess the adequacy of the 0.5 ns annealing time. 4.4 After simulated annealing, NPT MD simulations at 300 K and 1 atm for at least 5 ns were performed on each model system for comparison to experimental studies as well as other simulations. These simulations were the production runs from which structural properties were calculated and analyzed. In the case of the cylinder and sphere-rod models, more than 6x longer simulations were performed (~34 ns) to assess equilibration and convergence. Production Runs 4.5 NAMD was chosen as the MD code for all of the simulations in this study because of its excellent parallel scalability and performance. All simulations employed classical atomistic MD using NAMD [27] with the particle mesh Ewald (PME) method and a 1 Å mesh spacing for long-range electrostatic interactions, a 12 Å short-range cutoff with a smooth switching function from 10 to 12 Å, a Langevin thermostat (300 K temperature, 0.5 ps-1 damping coefficient) or velocity rescaling (300 K temperature) and Performance 41 a modified Nose-Hoover Langevin piston (1 atm pressure, 200 fs oscillation period, 100 fs damping decay timescale) for temperature and pressure control, respectively, and a time step of 1.0 fs for production runs. Although MD simulations of the large systems of this study were computationally expensive (1 ns of simulation time for each system required ~24,000 CPU*hours, or the equivalent of running on 1,000 processors for 1 day), the computational resources were available for this project. 4.6 The total system density of the λ ~ 7.4 morphological models ranged from ~1.68 to 1.74 g/cm3, which is within ~5 to 10% of the experimental value (1.86 g/cm3 at λ = 7) [28]. The slightly higher value of λ and the lower initial density (~1.6 g/cm3 used to facilitate MC chain growth) of this study may partially account for the lower observed density compared to the experimental value. The forcefield and details of the morphological models may also contribute. Overall, such close agreement to the experimental value helps validate this work. Similarly, close agreement to the observed densities in previous simulation studies was also found. For example, Devanathan et al. reported a density of 1.78 g/cm3 at λ = 7 [29] using a random model with the same forcefield as that of this study. Close agreement to previous simulation and experimental work was also found for the higher water content random models. The density of the λ = 15 random models of this study was ~1.63 g/cm3 after ~20 ns of NPT (thermostat and barostat) simulation. Annealing did not seem to significantly affect these results. For comparison, the simulated density of Jang et al. at the same water content and using the same model (random), forcefield, and temperature was ~1.60 g/cm3 for the DR = 1.1 case Validation 42 (dispersed sequences). The corresponding experimental value is ~1.75 g/cm3 at λ = 15 [28, 30]. 5 CLUSTERS AND BRIDGES 5.1 Experimental elucidation of the nanoscale structure of the proposed hydrophilic reverse-micelle clusters in hydrated Nafion remains challenging and largely unresolved. In particular, the geometry and distribution of these domains continues to be debated, and the connecting structures in between domains are not understood. These bridging connections have never been directly observed at the molecular level, either in experiments or simulations, but have been hypothesized and discussed in the literature [2, 19]. Some of the proposed morphological models have assumed and even included such channels or bridges in their design (e.g., the sphere-rod model proposed ~10 Å diameter rods of comparable length between spherical micelles and the rod model proposed somewhat wider and longer rods). For those models that have not explicitly included such connecting structural elements (e.g., the cylinder, slab, sphere, and random models), these observed bridges may help explain how clusters may connect and percolate in such models. Percolation is believed to be vital for proton transfer between domains and through the membrane. In the case of those models that already incorporate such features, additional connecting bridges quickly formed, leading to a more connected network of clusters. These additional connections may play a vital role in the low percolation threshold of Nafion and help promote proton transport at low hydration. Background 44 To analyze the hydrophilic clusters, an algorithm and corresponding C++ computer code were independently developed to group together molecules according to a distance criterion and cluster definition. The distance criterion approximates hydrogen-bond lengths in order to group all of the hydrophilic molecules within a hydrogen-bond network together into a ‘cluster'. Molecules are recursively grouped together into clusters by starting at an unassigned molecule and searching for its neighbors that are within a surrounding sphere (the distance criterion is the radius of this sphere), and then searching for their neighbors and so on until no more neighbors are found (until recursion ends). Each neighbor within this distance cutoff is assigned to the same cluster as the starting molecule. Then, another unassigned molecule starts the next cluster and recursive neighbor search and this process repeats until all molecules have been assigned to clusters. Since this algorithm does not limit cluster size, all of the molecules may belong to only one very large cluster, or each molecule may belong to its own one-molecule cluster, or anything in between these two extreme cases may occur. Hydrogen-bond angles and other typical hydrogen-bonding criteria were not considered in this algorithm to simplify the analysis. These additional criteria are not expected to significantly alter the clustering analysis because of the focus on slower nanoscale structural details rather than faster highly localized ones. The expensive neighbor search of this algorithm benefits from the same fast grid-hashing technique mentioned in the description of the chain growing algorithm above. Two cluster definitions were used to classify ‘hydrophilic' species: (1) ‘W', which represents the oxygen atoms of water and hydronium, and (2) ‘W+S', which represents the sulfur and oxygen atoms of sulfonate, in addition to the species in the W definition. Hydrogen atoms were not explicitly used in 45 this definition in order to decrease computational expense. This approximation is consistent with not using more detailed hydrogen-bond criteria. This clustering code reads in the Cartesian coordinates (x,y,z) of the atoms in the cluster definition as well as the distance cutoff and periodic boundary conditions, groups the atoms into clusters, and then outputs the total number of clusters, the number of atoms in each cluster, and the coordinates of all of the atoms that belong to each cluster. This code requires less than one minute running on a single core CPU of an Intel Core2 Duo processor for the large Nafion systems of this study. 5.2 Rapid formation of intercluster connecting ‘bridges' was observed in all of the models during equilibration. The cylinder and slab models are clear examples ( Rapid Bridge Formation Figure 5) because they can easily be viewed in 2D, whereas the other models are better visualized in 3D. However, similar bridging phenomena occurred in the other models as well. Visual formation of the connecting bridges was observed within ~2 ns for all of the models. These bridges typically range from narrow water-wire to wider tube and mesh-like structures and span the distance between clusters. The bridges primarily contain water but also contain significant amounts of sulfonates and hydronium ions. The presence of sulfonates inside the bridges is surprising and does not seem to have been discussed in the literature. Sulfonates may actively participate and direct in bridge formation and the transport of ions through the bridges. This may contrast with the possibly less active role of surrounding the bridge walls previously proposed [2, 19]. Figure 6 shows a representative bridge from the cylinder model, containing sulfonates and hydronium ions in addition to waters. Careful inspection of Figure 5 also reveals 46 Figure 5. Bridges, sulfonate penetration, and cluster changes. Rapid inter-cluster connecting bridge formation, sulfonate penetration into clusters, and water/hydronium cluster surface-increasing deformation during postanneal equilibration (no restraints) @ 300 K and 1 atm, as seen from 2D face-on snapshot looking down (bird's eye view) through 3D cylinder (left) and slab (right) models. The initial and final (after 34 ns and 5 ns for the cylinder and slab models, respectively) hydronium ion locations are shown as white and red dots, respectively. The final locations of sulfonates and waters are shown in yellow and blue, respectively. The rest of the polymer is not shown for clarity. sulfonate and hydronium participation inside bridges, as seen from the yellow (sulfonate) and red (hydronium) dots within the blue (water) connecting structures between white circles, which represent the initial cylinder surface coated with hydroniums. Future work needs to be done to directly quantify the number and size distribution (e.g., diameter and length) of the bridges. Such an analysis would require a definition of bridges to distinguish them from the cluster definition. This bridge definition may involve coordination (lower water coordination may indicate bridge), distance (shortest distance through part of a cluster may indicate bridge), motion (less translational/rotational water motion may indicate bridge), cluster ID information (a change in cluster identity may indicate bridge), or some combination of these techniques. 47 Figure 6. Representative bridge between clusters. In cylinder model after 34 ns. A perspective view was used to enhance 3D depth perception, and a ball-and-stick (CPK) atom and bond representation was used. Inset contains the same bridge without the exaggerated depth perception (orthographic view). Oxygen, hydrogen, and sulfur atoms are shown in red, white, and yellow, respectively, except hydronium oxygen atoms are shown in gold to clearly distinguish them from those of water. Blue dots represent the hydrophilic cluster surfaces. The rest of the polymer is not shown for clarity. 5.3 Figure 7 Percolation shows the general decrease in total number of clusters over time for each model and each cluster definition. This number decreases as several clusters merge together into one large cluster of clusters via bridge formation, indicating enhanced network connectivity between hydrophilic domains. Even though a percolating network develops in each model, cluster count remains high due to the large population of very 48 Figure 7. Cluster count. Top left and top right plot total # of clusters vs. time (ns) for the W and W+S cluster definitions, respectively, for a cluster cutoff radius of 3 Å. Bottom left and bottom right plot the same for a cluster cutoff radius of 5 Å. The models are labeled c (cylinder), sr (sphere-rod), s (sphere), r (rod), sl (slab), rn (random, no-anneal), and rp (random, postanneal), as illustrated in Figure 1. small groups of lone molecules. By definition, the clustering algorithm assigns each molecule to a ‘cluster' regardless of size, so even a single isolated water molecule will be counted as a cluster. Therefore, a final cluster count of unity is extremely unlikely and not preferred. Model differences and the impact of the cluster definitions on total cluster count suggest the vital role of sulfonates within the clusters and bridges. For a 5 Å hydrogen-bond cutoff and for the sphere-rod and cylinder models, including sulfonates in the definition of a cluster led to a ~4x decrease in the observed number of clusters. Without 49 including sulfonates in the definition, many more clusters are counted due to the fragmented structure of water and hydronium near ‘hidden' sulfonates. The cluster count of the random model systems experienced the biggest change (over one order of magnitude!) with cluster definition (highest and lowest cluster count for W and W+S cluster definitions, respectively), suggesting more penetration of sulfonates into clusters and bridges in those models, perhaps due to a more dispersed and closer spaced domain structure. The slab model may have the least penetration of sulfonates as suggested by it having the largest cluster count of all models according to the W+S definition. Sulfonates may have more difficulty penetrating deeply into the slabs due to more confined cluster geometry. This may also be due to the small sample size (one) of the model and the noticeable variation in initial sulfonate density between slabs (discussed below) due to the Monte Carlo chain nucleation and growth algorithms and the small number of initial slabs (six). The rod model initially started with a cluster count of one (W cluster definition) because of its built-in percolation design. The peak in cluster count (W definition) followed by a gradual decrease shows the perturbation imposed by annealing and the system response afterward. This may be due to charge saturation and overly strong aggregation and deep cluster penetration of sulfonates. In the case of the W+S definition, the monotonic decrease is due to the high initial number of sulfonates far from the clusters during annealing followed by a fast drop in this number as these sulfonates begin to join clusters. Cluster count was also calculated for a more restrictive hydrogen-bond cutoff distance (3 Å). Going from the looser criterion (5 Å) to the more restrictive one resulted in approximately an order of magnitude increase in cluster count for both cluster 50 definitions. Also, the more restrictive cutoff distance was much more sensitive to dynamics (e.g., bond vibrations, angle bends, and translational and rotational motion) than the looser one, as demonstrated by the larger fluctuations. As seen in Figure 8, the size of the biggest cluster rapidly increased to over 50% and 90% in all of the models within 5 ns for the W and W+S cluster definitions, respectively, and for the longer hydrogen-bond cutoff. For the shorter hydrogen-bond Figure 8. Size of biggest cluster. Top left and top right plot the size of the largest cluster (expressed as percentage of total # of molecules in cluster definition) vs. time (ns) for the W and W+S cluster definitions, respectively, for a cluster cutoff radius of 3 Å. Bottom left and bottom right plot the same for a cluster cutoff radius of 5 Å. The models are labeled c (cylinder), sr (sphere-rod), s (sphere), r (rod), sl (slab), rn (random, no-anneal), and rp (random, postanneal), as illustrated in Figure 1. 51 cutoff, percolation was not observed according to the W definition but was observed for the W+S definition, suggesting the important role of sulfonates in the bridges. The size of the biggest cluster is expressed in terms of the proportion of molecules within the largest cluster compared to the number in the entire box, where the actual molecules counted depend on the specific cluster definition. For example, a size of 50% according to the W cluster definition means half of the waters and hydroniums in the system are found in this cluster. A cluster that large is a strong indicator of percolation, meaning that single cluster is really a cluster of clusters composed of many connected smaller clusters. Under the W+S definition, the cylinder model had much smaller fluctuations than the sphere-rod model in the biggest cluster size plot for both hydrogen-bond cutoff distances, suggesting shorter-lived bridges in the sphere-rod model and more stable ones in the cylinder and others. For the W definition and the more restrictive cutoff, the rod model actually experienced a decrease in biggest cluster size due to sulfonates displacing waters in bridges. The cylinder model in this same plot experienced big fluctuations, suggesting short-lived water bridges (without sulfonates) between clusters. For the longer cutoff, the initial size of the biggest cluster in the rod and random models was over 90% using either cluster definition. By design, the rod model started with a 100% biggest cluster size according to the W definition, meaning all of the waters and hydroniums in the system were connected together into a single cluster or percolating network of rods. Although a small fraction (~2%) of these waters and hydroniums left this perfectly connected network and formed tiny, isolated ‘islands', the initial network remained connected and stable over time, suggesting strong self-assembly driving forces 52 that maintain rod connection. The formation of such islands is reasonable and expected because of the complex nanophase environment. The lifetime of these islands are unknown. Since both random model systems came from higher hydration systems run for over 20 ns, they also started with a percolating network and maintained percolation over time. The size of the biggest cluster in the other models (cylinder, sphere-rod, sphere, and slab) started small (< ~20%) and rapidly grew to include essentially all of the initial clusters. The remaining percentage consists of many small clusters or lone molecules. This rapid cluster growth clearly demonstrates bridge formation. As bridges form, the clusters link and merge together, leading to a percolating network of clusters. This is an indirect yet conclusive measure of percolation and cluster bridging. It demonstrates that all of the models strongly prefer a percolating network of clusters at this water content, in agreement with percolation threshold experimental studies [31]. This merging of clusters via bridge formation and percolation may promote proton transport through the membrane. Upon comparing the two cluster definitions, a significant difference in the final size of the biggest cluster and the rate of this cluster's growth was observed. By including sulfonates in the cluster definition, a larger and hence better connected network formed in a shorter time interval, indicating faster percolation and the formation of additional bridge connections directly involving sulfonates and suggesting the important role of sulfonates in guiding percolation. Applying the W+S definition also led to noticeably narrower deviations between models. This may suggest a universal role of sulfonates in the percolation behavior of Nafion, independent of morphological model 53 specifics. Overall, the relatively tight convergence of all of the models in terms of cluster count and biggest cluster size is surprising. Cluster geometry and distribution do not seem to significantly impact these measures. At first glance, the difference in the size of the fluctuations in cluster count and biggest cluster size with cluster definition seems inconsistent. Without considering the large difference in cluster count scale for the two definitions (Figure 7), the W definition apparently exhibited smaller fluctuations than W+S. However, the W definition also experienced larger fluctuations in biggest cluster size. This discrepancy is resolved by considering the difference in scale in cluster count, which is almost an order of magnitude. Using a percentage change basis, the fluctuations in cluster count were comparable for both definitions (+/- ~10%). Using an absolute scale, however, the fluctuations were actually much larger in the W definition (+/- ~100) compared to the W+S definition (+/- ~25). The larger fluctuations in both cluster count and biggest cluster size for the W definition may result from the motion of sulfonates within bridges and clusters, which should disrupt the hydrogen-bond network of water and hydronium. 5.4 Figure 6 Sulfonate Participation and Penetration clearly shows active participation of sulfonates in intercluster bridges within the cylinder model. Similar sulfonate-containing bridge structures were observed in the other models as well. Closer visual inspection reveals that sulfonates are part of the narrow hydrogen-bond network in these channels and agglomerate together with only a hydronium ion and a few waters between them, suggesting relatively low local hydrations (λ ~ 3 to 5) within the bridges even though the bigger clusters are at bulk hydration (λ = 7.4). Further bridge analysis may better quantify this hydration difference 54 between clusters and bridges and explore the distribution of hydration in bridges. Lower hydrations coupled with the presence of sulfonates in the bridges may suggest higher local viscosity and lower diffusion within bridges. The rapid formation of intercluster bridges and the participation of sulfonates within these bridges suggest the strong driving force for percolation and the strong hydrophilicity of the sulfonates at this water content. Future studies aimed at understanding how different polymer architectures and acid headgroups may affect these driving forces may be needed. In addition to sulfonate participation within the bridges, significant sulfonate penetration into the clusters was observed. Sulfonates were observed to burrow throughout the clusters, at varying depths (as seen from the yellow dots within the white outlines of the clusters in Figure 5). Little is known about the location of the sulfonates with respect to the clusters from experimental studies. Most morphological models have proposed a sort of surface coating film of sulfonates surrounding the cluster surfaces, but essentially nothing is known about how many, if any, sulfonates are inside the clusters and what the corresponding sulfonate spatial distribution may be within the clusters. The presence of sulfonates within the clusters may strongly affect proton transport. Petersen et al. [17] suggested strong proton ‘sink' behavior of the sulfonates, as observed from anticorrelated proton hopping (Grotthuss mechanism) vs. vehicular (Brownian) proton diffusion as well as strong time and spatial correlations from the van Hove sulfonate-proton correlation function, which is a radial distribution function over time as well as over distance. With respect to the cluster penetration behavior of the sulfonates in this work, the proton sink behavior of the sulfonates just mentioned may reduce proton 55 transport within the clusters, owing to the more dispersed presence of harder-to-avoid potential proton sinks throughout the hydrophilic domains. Further work aimed at investigating how sidechain length and architecture may alter this cluster penetration, and how this may in turn affect proton transport and the percolation threshold, may be insightful for future PEM design. Further studies involving detailed SCI-MS-EVB [16, 17] proton transport simulations through and near these bridges should be carried out to better understand how these bridges affect proton transport, especially at lower hydrations near the percolation threshold. Such simulations may be compared with proton transport studies from Donnan dialysis [32] and electro-osmotic drag experiments [33, 34]. Intuitively, these bridges may largely determine or affect the rate-limiting step of proton conductance through the membrane. 6 SULFONATES 6.1 Radial distribution functions (RDFs) were calculated according to the following expression: Background gab (r) = pab (r) pid ab (r) (2) where gab(r) is the RDF between species a and b, pab(r) is the actual (real) probability of finding any b-particle a spherical distance r away from any particle of a, and pid ab(r) = ρbVs(r) represents the corresponding probability in an ideal gas, where ρb is the total density of b-particles (Nb/V, where Nb and V are the total number of b-particles in the system and V is the total system volume) and Vs(r) is the volume of a spherical shell (histogram bin) located at r. pab(r) is measured by counting the number of b-particles within spherical shells centered on a-particles and then performing an ensemble average of the histogram counts, normalized by the number of timestep samples and the number of a-particles in the system. The number of radial distribution function bins (spherical shells) was 1500, corresponding to a bin size of ~0.1 Å given half of the typical box lengths (~300/2 Å = 150 Å) of this study. Due to periodic boundary conditions (PBCs), RDFs cannot be calculated beyond distances of half the box size without resorting to 57 extrapolation techniques because the maximum pair separation distance (and hence the maximum r) within a system under PBCs is one half of the smallest box dimension. 6.2 Rapid sulfonate aggregation around the water/hydronium clusters was observed during simulated annealing in all of the nonrandom models, in agreement with the proposed morphological models and as a consequence of the strong electrostatic interactions between the sulfonates and hydronium ions. As a clear example, Aggregation Figure 9 shows this aggregation in the cylinder and slab models. These models have been chosen Figure 9. Rapid sulfonate aggregation around clusters. As seen from 2D face-on snapshot looking down through 3D cylinder (left) and slab (right) models during simulated annealing of the polymer (while holding the water/hydronium fixed) @ 600 K. The initial and final (after ~0.5 ns) sulfonate locations are shown as red and yellow dots, respectively. The hydronium ions are shown in white and coat the cluster surfaces. The water and the rest of the polymer are not shown for clarity. 58 as the example because of their simple 2D nature. For the other models, this aggregation phenomenon is not visualized as easily in only 2D, although it still occurs in them. Initially, the sulfonates were approximately uniformly distributed throughout the spaces between clusters because no constraints were placed on them during MC chain growth in order to predict rather than assume their spatial distribution. Although most (~85% according to the sulfonate interfacial coverage) of the sulfonates moved near (within ~5 Å), on the surface of, or inside the clusters within a very short time (~0.5 ns), a small fraction remained far away. More than 10x longer annealing (~6 ns) did not seem to significantly alter this, perhaps demonstrating very high charge density saturation of the sulfonates, which prevents further agglomeration. Due to their mutual repulsion, the sulfonates seem incapable of moving closer together without a higher dielectric medium to increase charge screening. This charge density saturation may suggest a preference for clusters with more surface area that will allow the sulfonates to spread out more, thus decreasing the charge density, while closely interacting with the hydronium ions and water. 6.3 Figure 10 Mutual Spacing and Proximity shows RDF and structure factor (SF) plots for several combinations of species and will be used repeatedly in the remaining sections of this dissertation. The S-S and S-Ow RDF and corresponding SF plots, respectively, where S represents the sulfur atom in sulfonates and Ow represents the oxygen atom in water and hydronium are shown in Figure 10. S-S may be the most significant RDF in terms of studying proton transport because it reveals the mean sulfonate-sulfonate distance (based on the location of the 1st peak), which is believed to strongly influence how rapidly protons can hop 59 Figure 10. Radial distribution function and structure factor plots. Ow-Ow RDF (top left) with inset, Ow-Ow SF (top right) with error bars for the cylinder and sphere-rod models, S-S RDF (middle left), S-S SF (middle right), S-Ow RDF (bottom left), and S-Oh RDF (bottom right). The models are labeled ‘c' (cylinder), ‘sr' (sphere-rod), ‘s' (sphere), ‘r' (rod), ‘sl' (slab), ‘rn' (random, no-anneal), and ‘rp' (random, postanneal) as illustrated in Figure 1. 60 between sulfonates [35], according to the Grotthuss mechanism. Protons need to have a short distance from one sulfonate to the next for rapid transport. Even though charge repulsion and saturation spread the sulfonates further away from each other, their strong attraction to protons and water brings them closer together. This is a delicate balance of competing forces. Short sulfonate-sulfonate distances are believed to enhance proton conductance. Indeed, lower equivalent weight (less molecular weight per sulfonate or higher sulfonate content) Nafion membranes demonstrate higher proton conductance [36]. The mean sulfonate-sulfonate distance from S-S RDF in this work was ~5 to 8 Å, in close agreement with previous simulation studies [5, 37, 38]. The slab model had a high peak (1st) at r ~ 5-8 Å, a broad low peak (2nd) at r ~ 10-20 Å, and a deep dip from r ~ 25-40 Å. This model had the highest peak in S-S SF probably because of the deep dip. The other nonrandom models (cylinder, sphere-rod, sphere, and rod) exhibited similar features as the slab model but all had higher 1st peaks, lower 2nd peaks, and weaker dips (shifted to closer distances). In contrast, the random models had the highest and sharpest 1st peak at r ~ 5 Å, followed by a deep dip at r ~ 10-20 Å in place of the usual 2nd peak. Since the RDF dip was shifted to shorter distance and less broadened, the usual peak in S-S SF was shifted to higher k (~0.26 Å-1) and significantly broadened. The behavior of the random models in this work was in close agreement with the results of Venkatnathan et al. [38] in terms of 1st peak height and the short location of the deep dip. The dip location of the latter was shifted to even shorter distances (~7 to 10 Å). This may be due to a slightly lower hydration (λ = 6) and a much smaller box size (~40 Å) used in that work. 61 In all models, the S-Ow RDF (Figure 10) had a sharp peak at r ~ 4 Å, a small peak at r ~ 7 Å, and a broad dip at longer distance. The two peaks reveal the 1st and 2nd solvation shells of water and ions surrounding each sulfonate. The peak locations and shapes resemble those of previous studies [29, 38]. The height of the 1st peak in all of the models was much higher than that of the S-S RDFs, demonstrating the strong S-Ow attraction and the repulsive S-S interaction. The behavior of the nonrandom models was very similar. The random models once again had a shifted dip towards shorter distances. They also had lower 1st and 2nd peaks. The corresponding features in S-Ow SF (not shown) reveal a negative region (due to neglecting the imaginary part) at short k (< 0.07 Å-1) that suggests destructive interference, followed by a weak ionomer peak at k ~ 0.13 Å-1, in close agreement with the experimental ionomer peak location from the total structure factor. The random models had a deeper interference region and, as before, a flattened peak shifted to higher k (shorter distance). The slab model had the sharpest and highest peak once again, possibly indicating stronger nanophase segregation and sulfonate correlation. 6.4 In all models, the S-Oh RDF ( Ion Pairing Figure 10) had a very high and sharp peak at r ~ 4 Å, corresponding to the contact ion-pair position. The solvent-separated ion-pair position (r ~ 6 Å) was weakly present. This peak may be better resolved using Os-Ow RDFs, where Os represent oxygen atoms in sulfonates, but these RDFs have not been included here because of the well-known limitations of using explicit hydronium ions. Dissociable models (e.g., MS-EVB) would be needed to accurately probe the short length-scale behavior of excess protons and to accurately distinguish between contact ion-pair and 62 solvent-separated ion-pair preferences. Such short-range structural details have not been focused on in this work in order to concentrate attention on probing the larger cluster structures. For more detailed short-range structure, using the MS-EVB methodology, please refer to the work of Petersen et al. [17, 18]. That work found a strong tendency for the solvent-separated ion-pair position (2nd solvation shell) with an Eigen-type structure of the excess proton (H9O4 +, a hydronium-like ion solvated by three water molecules). The overall features of the S-Oh RDF were similar to those of the S-Ow RDF, with the exception of relatively smaller peaks in the case of S-Ow due to the very strong attraction of negatively charged sulfonates and positively charged hydronium ions. Surprisingly, S-Oh SF (not shown) more closely resembled S-S SF rather than S-Ow, despite the close visual comparison with S-Ow RDF. This may be due to more subtle quantitative differences (e.g., peak heights) between the RDFs. Similarly, the S-Oh SFs exhibited more model differences than visually observed in the corresponding S-Oh RDFs, which had fairly comparable curves between all models, although differences in peak heights between models were apparent nonetheless. 7 INTERFACE 7.1 A numerical algorithm for calculating solvent accessible surface area (SASA) and volume using a grid based approach adapted from Faulon et al. [39] was developed for this study. The algorithm proceeds as follows: Background Step 1) Find atomic, excluded, 'possible' included or interstitial, and external points (cells in grid); external points are the remaining cells after this step. Each atom is superimposed over the nearby cells in the grid. The distance from this atom to each nearby cell is then used to determine the state of that cell. In order of highest-to-lowest precedence, the cell states are: ‘atomic' (inside an atom) r <= rvdw where, r = center-to-center distance from an atom to a cell rvdw is the van der Waals (VDW) radius of an atom ‘excluded' (less than a probe radius away from an atom) rvdw < r <= rexc where, rexc = rvdw + rprobe 64 rprobe is the solvent probe radius (typically 1.5 Å, comparable to a water molecule radius) ‘possible' included or interstitial (within some maximum void size) rexc < r <= rext where, rext = rvdw + rmaxvoid rmaxvoid is an arbitrary maximum void radius, which is typically 2*rprobe, meaning the largest acceptable void or interstitial must be less than twice the size of the solvent probe ‘external' (outside of this range) rext < r Step 2) Recursively find included points by starting at external points and spreading through the grid of possible points until blocked by excluded or atomic cells. Included cells are external cells, but they are labeled separately from the external cells already found for convenience. Included cells must touch at least one face (alternatively, this restriction can be reduced to only having to touch at least one edge or one corner) of at least one other included or external cell, basically forming a chain of cells to the outside world (where the probe lives). All 'possible' cells will either become included (external) cells or interstitial cells. Interstitial cells are the remaining cells after this step and do not touch included or external cells; they are completely surrounded by atomic or excluded cells, 65 basically forming the void space or interior (the inside world that the probe cannot reach). From this step forward, only possible cells will change; no other cells (atomic, excluded, or external) will change. Recursion is done in stages to prevent running out of stack memory space when the recursion depth goes too deep. Step 3) Find surface points (faces), which lie on the faces between solvent accessible (included or external) and solvent excluded (atomic, excluded, or interstitial) cells. Step 4) Calculate solvent accessible surface area as: SASA = n*a where, SASA is the solvent accessible surface area n is the number of surface points (faces) a is the area of a face, which is the same for all faces due to uniform cell sizes in the grid Step 5) If desired, calculate the fraction of this surface covered by a particular species [5] by counting the number of ‘covered' surface points (faces) within an arbitrary cutoff distance of these atoms and dividing by n, the total number of surface points found in Step 3. The fraction of atoms covered by surface points can be similarly calculated by counting the number of atoms within the cutoff distance of the points. The number of cells in the grid can be increased to decrease cell size (typically 0.5 Å) and thus increase resolution and accuracy subject to the amount of available 66 computer memory and time. The expensive distance search of this algorithm benefits from the same fast grid-hashing technique mentioned in the description of the chain growing algorithm above. The same two cluster definitions mentioned previously (W and W+S) were used to determine the cluster boundaries (the hydrophilic-hydrophobic interface). The corresponding C++ code reads in the Cartesian coordinates (x,y,z) and VDW radii of the atoms in the cluster definition as well as the probe radius, periodic boundary conditions, grid cell size, and (optionally) the coordinates and radii of the covering atoms for the surface coverage calculation in Step 5. The code outputs the total cluster volume (vol%), the total solvent accessible surface or interfacial area (Å2), and (optionally) the total surface coverage (% of surface covered by atoms and % of atoms covered by surface). Using a grid cell size of 1 and 0.5 Å, respectively, the code requires ~1 and ~5 minutes on a single core CPU of an Intel Core2 Duo processor and ~0.5 and ~2 GB memory for the large Nafion systems of this study. This difference in resolution typically changes SASA less than ~3%. Commonly used surface area codes, such as MSMS [40], Surf [41], and g_sas from GROMACS [42], were not able to successfully perform this calculation for such large systems due to the difficulty of handling multiple clusters (more than one surface), over prediction of surface area due to incorrect interstitial regions (mistaken surfaces were visually observed in void spaces where the probe should not be able to reach), and memory or data structure errors (segmentation faults) using default compiler, source code, and options. Using Surf, the order of atoms affected the results, leading to segmentation faults in many cases. In this work, the probe radius was 1.5 Å and the grid cell size was 0.5 Å. 67 7.2 Deformation of the clusters, due to rearrangement of sulfonates, hydroniums, and waters was observed near the beginning of the equilibration, which did not involve any restrained or fixed atoms, after annealing ( Deformation and Structure Figure 5). Initially, clusters were assumed to be rigid and smooth in all nonrandom models. However, the cluster surfaces (interfaces) became rougher and more diffuse during equilibration. In the case of the cylinder model, for example, the rigid and straight cylinders turned into slightly bending tubular structures with wider and narrower sections along the tubes and a more tortuous path. The elongated fibrillar structure of the cylinders remained fairly stable during the short timescales (~34 ns) explored in this work. For other models, similar cluster transformations were observed. This is likely a result of bridge formation and sulfonate penetration into the clusters as well as sidechain movement and water and hydronium rearrangement around the sulfonates. This phenomenon contributed to higher SASA and may be related to alleviation of high charge density of sulfonates and hydroniums. Figure 11 shows water iso-density surface changes in the high water content (λ = 15) random model (annealed) during equilibration. This illustrates the deformation and movement observed in the random models and their complex structure. Figure 12 shows a representative water density slice from the same system. Notice the bulk-like water density (~1 g/cm3) in some of the clusters in red. This slice also suggests percolating bridge pathways between islands as seen from the connections between hot-colored regions. 68 Figure 11. Water iso-density surfaces of random model. On left, H2O/H3O+ iso-density surfaces shown as red solid blobs with sulfur atoms in yellow, averaged over ~10 ns at 300 K, 1 atm, and λ = 15 using 1 nm grid spacing and an iso-value of ~10% of the maximum density; on right, same red solid blobs compared to pink translucent blobs that represent the initial iso-density surfaces. Figure 12. Water density profile of random model. H2O/H3O+ density (g/cm3) at a representative z-slice of the box and for a single timestep (after ~10 ns). x and y values represent positions (nm) along these respective dimensions in the box. Hot and cold colors indicate regions of high and low water density, respectively. 69 In all of the models, the clusters increase in volume ~40% (There is a ~8 vol% increase) to give the sulfonates and hydroniums more room to spread out, thus alleviating the charge density saturation (Figure 13). The spread between models is modest (~2 vol%), suggesting a universal trend. The W+S cluster definition exhibits a smaller increase (~5 vol% or ~20% change), suggesting a substantial inclusion of sulfonates within or on the surface of clusters during annealing. This may be understood by Figure 13. Top and bottom plot total cluster volume (expressed as a percentage of the total box volume) vs. time (ns) after annealing for two different cluster definitions: W, representing water and hydronium oxygen atoms, and W+S, representing sulfur atoms and sulfonate oxygen atoms in addition to water and hydronium, respectively. The models are labeled ‘c' (cylinder), ‘sr' (sphere-rod), ‘s' (sphere), ‘r' (rod), ‘sl' (slab), ‘rn' (random, no-anneal), and ‘rp' (random, postanneal) as illustrated in Figure 1. 70 considering that the sulfonates themselves take up volume, yet the difference between cluster definitions initially (at the end of annealing) was much greater than that at the end of the production runs, indicating that a large proportion of sulfonates were already part of clusters at the beginning of the postannealing production runs, in agreement with the initially high sulfonate coverage. The random model systems had higher volumes in part because they came from higher water content systems. Cluster density increases ~35% going from the W (~0.61 g/cm3) to the W+S (~0.83 g/cm3) cluster definitions because of the significant presence of sulfonates in the clusters. The low density within clusters suggests the presence of sidechains inside clusters. Intuitively, mechanical stability suggests that cluster densities should be comparable to the bulk system density (~1.7 g/cm3). By accounting for the mass of the sidechains in the cluster density calculations, the difference between cluster densities and bulk density is resolved. By including one ether linkage of the sidechain (the one nearest the sulfonate), in addition to water, hydronium, and sulfonate, cluster density increases from ~0.83 to 1.3 g/cm3. Accounting for the entire sidechain (both ether linkages) increases cluster density further to ~1.9 g/cm3, in close agreement with the bulk density, suggesting the penetration of most of the sidechain into clusters, in addition to the sulfonate end. The presence of the ether linkages of the sidechain may affect proton transport within clusters. Due to the recently observed amphiphilic nature of excess protons in aqueous environments [43], these more hydrophobic moieties may actually enhance proton transport by increasing the hydrophilic-hydrophobic interfacial area, giving protons more preferred paths to travel. 71 7.3 Figure 14 Charge Saturation and Coverage shows the specific solvent accessible surface area (SASA) or hydrophilic-hydrophobic interfacial area of each model over time. These surfaces surround the clusters and represent the interface between the nanophases. Each nonrandom model rapidly increased ~100% (doubled) in SASA for the W cluster definition, indicating that none of the proposed models correctly account for interfacial Figure 14. Interfacial area. Top and bottom plot specific solvent-accessible surface area (SASA) or interfacial area (m2/g) vs. time (ns) for the W and W+S cluster surface definitions, respectively. The models are labeled ‘c' (cylinder), ‘sr' (sphere-rod), ‘s' (sphere), ‘r' (rod), ‘sl' (slab), ‘rn' (random, no-anneal), and ‘rp' (random, postanneal) as illustrated in Figure 1. 72 area and suggesting the need for model refinement. The rapid rise in interfacial area may alleviate high charge density of sulfonates and excess protons. By increasing cluster surface area, sulfonates and protons gain more space to interact, preventing unfavorable charge repulsion between like species. Surprisingly, SASA only decreased ~5% with the inclusion of sulfonates into the definition of clusters (W+S). The W+S cluster definition had been anticipated to result in much lower SASA because sulfonates were ‘fragmenting' the water structure into groups of small clusters, thus increasing SASA in terms of the W definition. However, the smaller than anticipated decrease in SASA may suggest cluster fragmenting behavior due to the presence of the ether linkages of the sidechains in the clusters, which were not included in either cluster definition but which has been mentioned above in terms of cluster density observations. As expected from its small clusters, the rod model had the highest SASA of the nonrandom models. Surprisingly, the sphere-rod and sphere models had lower SASA than the slab model and even the cylinder model, both of which have larger clusters with intuitively less initial area. In reality, the slab model initially had more interfacial area than the sphere-rod and sphere models as well as the cylinder model, probably because of initially rough edges. The observed model differences may be due to increased surface ‘roughening', whereby a more diffuse and undulating surface forms to spread out sulfonates and accommodate more charge. The random model systems exhibited the highest area and experienced an area decrease, probably resulting from the higher water content systems from which they were initially derived. Nevertheless, the random models may be viewed as upper bounds to SASA; the nonrandom models may be viewed as lower bounds. Good convergence 73 and equilibration of this property are observed during the course of the simulations as evidenced by small fluctuations and the lack of drift. The initial spread of SASA in the nonrandom models for the W cluster definition was ~350 to 600 m2/g. The final spread was narrower: ~950 to 1,100 m2/g. Combined with the higher limit from the random models, all of the models exhibited close agreement in SASA. The large magnitude of SASA (~1,000 m2/g) apparently has never before been observed or discussed in the literature. To illustrate its significance, this value is comparable to that of zeolites [44], which are well known for their extremely high surface areas routinely used in the fields of catalysis and materials science. Further analogy suggests that ~5 g of hydrated Nafion would have an interface large enough to cover an entire football field! Such a large interfacial area in such a small amount of material underscores the strong nanophase segregation behavior of this ionomer and the role this interface has in shaping structure and morphology. Intuition suggests a significant effect on dynamics as well. Clearly, morphological model refinement and improvement are needed to account for such high interfacial area of hydrophilic domains. Figure 15 shows sulfonate coverage on the interface over time for both cluster definitions and for a sulfonate ‘sphere of influence' of 5 Å [5]. This spherical cutoff approximates the influence each sulfonate exerts on its local environment in terms of an excluded body effect. Sulfonate interfacial coverage decreased for both definitions, although the W definition resulted in a rapid initial increase before the gradual decrease. The resulting peak (~65 to 80%) in the W definition occurred during annealing. This indicates initial saturation of sulfonates on the interfaces before equilibration, as discussed previously with regard to charge density saturation. Before annealing, the 74 Figure 15. Top and bottom sulfonate surface (interface) coverage, expressed as a percentage of SASA within the ‘sphere of influence' of each sulfonate, vs. time (ns) for the W and W+S cluster surface definitions, respectively. The sulfonate radius of influence was 5 Å. The models are labeled ‘c' (cylinder), ‘sr' (sphere-rod), ‘s' (sphere), ‘r' (rod), ‘sl' (slab), ‘rn' (random, no-anneal), and ‘rp' (random, postanneal) as illustrated in Figure 1. coverage was low (< 25%) because the initial model building algorithm did not bias sulfonate placement, resulting in fairly uniform ‘dispersed' sulfonate distributions between clusters and relatively few sulfonates near clusters. During annealing, sulfonates rapidly over-aggregated on the clusters due to strong electrostatic interactions between sulfonates and hydronium ions, which initially coated the cluster surfaces by design and which were held stationary during annealing. After annealing, the models responded by increasing SASA and cluster volume, resulting in lower sulfonate coverage (~65%). 75 Within the first 1-2 ns of equilibration, all sulfonates were quantitatively observed to be within 5 Å of the interface in all of the models, as measured by the 100% coverage of sulfonates by the cluster surfaces. This indicates complete aggregation and solvation of sulfonates, in agreement with the morphological models and previous simulations [2, 5, 9]. This also demonstrates the rapid motion of sulfonates and especially water and hydronium despite the glassy backbone structures, in validation of the building, annealing, and simulation methodologies used in this work. The difference in sulfonate coverage between cluster definitions was small (<10%) because of the small accompanying change in SASA. In contrast to the low initial value and peak of the W definition, the W+S definition started high (~75 to 90%) and monotonically decreased. The high initial value was a result of including sulfonates as part of the definition of clusters and then measuring the coverage of those sulfonates on the clusters (essentially measuring how much they cover themselves as well as water and hydronium). The rod model had the lowest surface coverage and smallest peak using the W definition because it had the highest SASA of the nonrandom models. Also, it probably had more overlapping ‘covered' regions of area, which would result in less overall covered area given the same total area and number of sulfonates, compared to the random models, as seen from its lower surface coverage under the W+S definition despite its lower SASA. The random models had lower surface coverage and no observable peak because of their initially high SASA. The general trend of models with higher SASA having lower surface coverage held true, except for the case of the rod and random models just mentioned. This may indicate that either overlap of covered regions was not 76 significant or that it was universally consistent regardless of model. Based on obs