Atomistic simulation and molecular dynamics of model ... - CiteSeerX

Atomistic simulation and molecular dynamics of model systems for perÑuorinated ionomer membranes James A. Elliott,¤ Simon Hanna,*a Alice M. S. Elliottb and Graham E. Cooleyb a H. H. W ills Physics L aboratory, University of Bristol, T yndall Avenue, Bristol, UK BS8 1T L . E-mail : S.Hanna=bristol.ac.uk b National Power plc, Harwell International Business Park, Harwell, Didcot, Oxfordshire, UK OX11 0QA Received 30th June 1999, Accepted 24th August 1999

An atomistic model for perÑuorinated ionomer membranes (PIMs), in particular NaÐon materials, is presented and used in conjunction with NVT molecular dynamics simulations to investigate the dynamic and conÐgurational properties of these polymers. It is found that the electrostatic term in the force Ðeld is responsible for the formation of an apparently phase separated morphology which is selectively conductive, favouring the passage of cations. SpeciÐcally, the mobility of H O` ions is found to be D3.2 times greater 3 than that of OH~ ions, under the application of an external electric Ðeld. This phenomenon is shown to be consistent with a jump di†usion model of ion transport in PIMs. There is also evidence for the existence of water in two distinct environments in the simulations : both tightly bound to ion exchange groups, and more loosely associated with the Ñuorocarbon matrix.

Introduction PerÑuorinated ionomer membranes Ionomers are polymeric materials which have the capacity to form intermolecular ionic bonds. The presence of even small quantities of ionic groups (under 10 mol%) in materials such as polyethylene and polystyrene gives rise to dramatically elevated glass transition temperatures, anomalous viscoelastic properties, prodigious water uptake and order of magnitude increases in conductivity.1,2 It is widely believed that the interesting properties of these materials derive from the microscopic phase separation of ionic material from the amorphous matrix, but the precise manifestation of this remains in question at present. However, it is known that the electrochemical properties of ionomers are highly sensitive to the water content of the polymer.3,4 PerÑuorinated ionomer membranes (PIMs), consisting of a polytetraÑuoroethylene (PTFE) backbone with sulfonic acid groups substituted at intervals along the chain, are of particular commercial interest due to their relatively high chemical and mechanical stability. The membrane studied in this paper is manufactured by E. I. Du Pont de Nemours and Company, under the trade name NaÐon. The structural formula of NaÐon is shown below in Fig. 1. The index m is usually equal to unity, in which case the value of n determines the ratio of polar to non-polar material in the membrane. It is conventional to refer to the equivalent weight (EW) of a membrane, which is deÐned as the mass of dry polymer (in grams) which contains one mole of sulfonate ion exchange groups.

ion transport properties at a fundamental level. Also, it was considered desirable to be able to compare the detailed structural information available from an equilibrated molecular model with existing experimental data from nuclear magnetic resonance (NMR)5h7 and infrared (IR)8,9 spectroscopy studies. The agglomeration of absorbed water in the models can be visualised directly, and used to aid the interpretation of small angle X-ray scattering (SAXS) data from PIMs10,11. Finally, the behaviour of the simulations in the absence of electrostatic forces, and in the presence of an external electric Ðeld, can help to elucidate the role played by charged interactions in the formation of a phase separated morphology. In the past, modelling of such highly complex systems as PIMs has been difficult due to the excessive computational demands of simulating thousands of atoms interacting over a sufficiently long time to witness the emergence of interesting structural features. However, thanks to improvements in computer speeds, there have been a number of recently published studies on related polyelectrolyte systems (e.g., see refs. 12 and 13). In the current paper, we consider the aggregated structure formed by 26 anions in an aqueous environment. The result is a realistic, detailed morphological paradigm for PIMs which is consistent with existing phenomenological models. A review of these can be found in ref. 14.

Purpose of molecular modelling The motivation for constructing an atomistic model for PIMs was to try to improve the understanding of their macroscopic ¤ Present address : Department of Materials Science and Metallurgy, University of Cambridge, Cambridge, UK CB2 3QZ.

Fig. 1 The structural formula of NaÐon. The index m is equal to unity, and therefore the value of n determines the ratio of polar to non-polar material in the polymer (i.e., the equivalent weight).

Phys. Chem. Chem. Phys., 1999, 1, 4855È4863 This journal is ( The Owner Societies 1999

4855

Method A molecular model of perÑuorinated ionomer membranes The models were built using up to four distinct molecular fragments, denoted ““ sulfonate ÏÏ, ““ water ÏÏ, ““ hydroxonium ÏÏ and ““ hydroxyl ÏÏ. It was decided to represent the membrane using a suspension of such fragments rather than a fully connected polymer chain in order to produce higher molecular mobilities, so that a model could undergo signiÐcant structural reorganisation over the course of a simulation. Also, the initial conÐguration of a polymer chain can exert a powerful inÑuence over the region of phase space explored by the molecular dynamics. To nullify this e†ect would have required averaging results from very many simulations, generated using di†erent starting conÐgurations. However, the use of molecular fragments means that the constraining inÑuence of chain-bending forces on ionic aggregation is less signiÐcant than in the real polymer. A sulfonate fragment is shown in Fig. 2. The outline on the right-hand side of the diagram shows the partial charges assigned to each atom using methods which are described in the following discussion. The sulfonate fragment represents a portion of the pendant side groups attached to the Ñuorocarbon backbone. The side group was severed from the chain to increase its mobility but it still possesses a distinct ionic sense, with the Ñuorinated top end being relatively non-polar and the oxygenated bottom end being relatively polar. The omission of the Ñuorocarbon links between ion exchange groups will result in an e†ective increase in the ratio of polar to non-polar material in the model (i.e., a lower equivalent weight) compared to the real polymer. The sulfonic acid group was represented as completely ionised, enabling the dissociated hydrogen ion to move freely between di†erent sulfonate fragments in the polymer, with the net charge of the combination being zero. The requirement for mobile hydrogen ions required careful thought in order to achieve realistic results. Lone protons proved very destructive to the simulations due to their extremely high charge to mass ratio and small size. Each proton was therefore associated with a water molecule, to form a hydroxonium ion, as shown in Fig. 3c. In this way, the cations were free to di†use between sulfonate fragments in a natural way, as opposed to being bound to a particular site. The hydroxonium ion was built from a water molecule (shown in Fig. 3b) with a net charge of ]e, and equal charges on each hydrogen. In some simulations hydroxyl ions (shown in Fig. 3a) were also used, which possess a net charge of [e.

Fig. 2 Diagram of a sulfonate fragment featured in the molecular models of PIMs. The atomic shading scheme is : medium grey for carbon or Ñuorine, light grey for sulfur and dark grey for oxygen. The outline on the right-hand side of the Ðgure shows the partial charges assigned to each atom, giving a net charge of [ e.

4856

Phys. Chem. Chem. Phys., 1999, 1, 4855È4863

Fig. 3 Diagram of the “ water-like Ï fragments featured in the molecular models : (a) hydroxyl ion, (b) water molecule and (c) hydroxonium ion. The atomic shading scheme is : dark grey for oxygen and light grey for hydrogen. The outlines show the partial charges assigned to each atom, giving a net charge of [e, zero and ] e on fragments (a), (b) and (c), respectively.

It should be emphasised that we have ignored the possibility of charge transfer between ions. It is not possible to take account of this e†ect in a classical molecular dynamics simulation, and ab initio calculations are currently out of the question for such large systems. Thus, only a lower bound on the rate of di†usion of protons in the membrane can be obtained. The complete models were assembled from their constituent parts using the Cerius2 software package.15 The standard procedure was to Ðll space with a given number of sulfonate fragments in random orientations, neutralise these with the addition of hydroxonium ions and then Ðnally add the required amount of water. In the cases where hydroxyl ions were used, neutralisation was achieved by replacing a corresponding number of water molecules by hydroxonium ions. The models can be characterised by just a single parameter, which is the ratio of water to sulfonate fragments. This can either be viewed as varying the equivalent weight at constant water content, or vice versa. The calculation of the partial charges of the molecular fragments was mainly handled by the molecular orbital package (MOPAC) using the AM1 Hamiltonian.16 A semi-empirical method was chosen over ab initio techniques due to the size of the molecular fragments under consideration. It was decided to neglect polarisability because accounting for polarisation e†ects would add little to the simulations, since the main driving force for structural reorganisation in PIMs is thought to be the minimisation of unfavourable electrostatic interactions. Although MOPAC was used to calculate partial charges for the sulfonate and hydroxyl fragments, it was decided to use a standard model for water. The TIP3P model has been well tested, and found to reproduce the thermodynamic and structural properties of liquid water extremely well.17,18 It was found that this model is more polar than MOPAC would predict. Hence, for the sake of internal consistency, hydroxonium ions were constructed from the TIP3P water molecule. There will inevitably be a slight polarity mismatch in the models between MOPAC and TIP3P fragments, but this is outweighed by the advantages of using a sound model for water. A modiÐed version of the Dreiding II force Ðeld19 was used for carrying out simulations of PIMs, in the absence of a suitable speciÐc force Ðeld, as the conÐgurational and dynamical properties were of primary interest.

The production of a speciÐc force Ðeld for PIMs is beyond the scope of this paper. However, an attempt to derive a speciÐc force Ðeld for polytetraÑuoroethylene (PTFE) was recently published by Holt et al.20 based on semi-empirical molecular orbital calculations on perÑuorohexadecane. The parameters which they obtained from MOPAC AM1 energies indicated higher barriers to torsional motion and stronger intermolecular attractions than in alkane systems. Reasonable intramolecular geometries and intermolecular packing arrangements for PTFE were obtained with these parameter sets. Although AM1 partial charges were used in the work described in the present paper, it was felt that torsional e†ects would not be a signiÐcant factor in the models, due to the lack of connectivity between molecular fragments and their small size. In addition, the inclusion of PTFE-like torsional behaviour would have required substantial modiÐcation to the van der Waals parameters of Dreiding II, causing serious repercussions throughout the rest of the force Ðeld. For these reasons, the force-Ðeld of Holt et al. was not used. The Hamiltonian for an atom in the Dreiding II force Ðeld is : H \ KE ] V ] V (1) bo nb where V and V are the bonded and non-bonded potential bo nb energy terms, respectively. These two terms can be further decomposed into : V \V ]V ]V (2) bo b h Õ V \V ]V ]V (3) nb vdW e ext where V is the bond stretching energy, V is the bond bending b h energy, V is the dihedral energy, V is the van der Waals Õ vdW energy, V is the electrostatic energy and V is the energy of e ext charged atoms in an external electric Ðeld. The non-bonded terms are only calculated for those atoms which are separated by more than three bonds, and a 6 Ó cut-o† was used for the short ranged van der Waals forces. The bonded terms were represented by the standard harmonic and cos(3/) potentials, with k \ 700 kcal mol~1 Ó~2, b k \ 100 kcal mol~1 degree~2 and k \ 1 kcal mol~1. The h Õ bond angle force constants of the SO group in the sulfonate 3 fragments were increased in order to represent the e†ects of electron delocalisation in this area, and to prevent the oxygen atoms from approaching each other too closely due to their large negative charges. The actual electrostatic interaction is not evaluated explicitly as the oxygens are 1È3 bonded neighbours. The modiÐed force constants are given in Table 1. The electrostatic interactions in the simulations were computed using an Ewald summation technique.18,21 In addition, a uniform external electric Ðeld was applied in some cases. The Ðeld strength used in the simulations (B107 V m~1) was considerably larger than a realistic value for the Ðeld inside a fuel cell (B104 V m~1). This increase was necessary in order to produce measurable changes in di†usion over the short timescale of the simulation.

The simulation time step was 1 fs, and simulations with larger time steps showed signs of instability. The results were converted to Cerius2 trajectory Ðles using custom written routines. These trajectories were then analysed using the Cerius2 package.

Results The benchmark model consisted of 26 sulfonate fragments, 26 hydroxonium ions and 74 water molecules (872 atoms, 6.3% water by mass). This model corresponds to an 1100 EW membrane under ambient humidity conditions, according to mass swelling isotherm data. The size of the simulation box was adjusted to give the correct overall density for the massweighted combination of ionomer and water. On examination of the molecular dynamics trajectories for this, and other, systems, the main feature that was apparent was the formation of a phase-separated morphology, which manifested itself as a clustering of the hydrated sulfonate groups. Phase-separation was only observed when the electrostatic term in the force Ðeld was used ; without electrostatics, the molecular species remained mixed. A fuller discussion of the morphology and its inÑuence on ionic di†usion is given at the end of the section. Molecular energetics The total bond-stretch, angle and torsion energies of the models, together with the non-bonded vdW and electrostatic potentials, were continually monitored during the NVT dynamics simulations. The three bond-energy terms, together with the vdW energy, are shown in Fig. 4 over 1 ns for a typical simulation at 300 K. The electrostatic energy is plotted using a di†erent scale in Fig. 5. It is immediately clear that the electrostatic energy is the dominant term in the force Ðeld Hamiltonian, eqn. (1), possessing a magnitude which is nearly 75% of the total energy after equilibration. This equilibration was manifest by a decrease in magnitude of the electrostatic energy by approximately 15% in the Ðrst 0.4 ns. After this period, there were random Ñuctuations with a standard deviation of 0.6% due to energy exchange with the thermostat. Fluctuations in the temperature of the system were found to be approximately 2% at 300 K. Di†usion The di†usion of the atomic species was examined by calculating a self-di†usion coefficient from their mean-squared displacement (MSD) during the course of a simulation. When no external Ðeld was applied, the self-di†usion coefficient D, hereafter just ““ di†usion coefficient ÏÏ, was obtained from the MSD

Simulation techniques The MD simulations were carried out using the DL–POLY software (version 2.10) from Daresbury Laboratory,22 using a Berendsen-type thermostat23 with a relaxation time of 0.4 ps.

Table 1 Bond angle force constants used in the original and modiÐed Dreiding II force Ðelds Bond angle

Dreiding II /kcal mol~1 degree~2

ModiÐed Dreiding II /kcal mol~1 degree~2

OÈSÈO OÈSÈC

100 100

200 150

Fig. 4 Graph of the bond energies and the non-bonded vdW energy as function of time for the benchmark model during a 1 ns NVT dynamics simulation at 300 K.


4857

Fig. 5 Graph of the electrostatic energy as a function of time for the benchmark model during a 1 ns NVT dynamics simulation at 300 K.

via the Einstein relation.18 However, in non-equilibrium simulations (i.e., with an external electric Ðeld), the di†usion coefficients were calculated from the steady-state particle Ñuxes.18,24 In each case, the di†usion coefficients were calculated by averaging the MSDs over all molecules of a particular species. The ensemble averaged MSDs were plotted during the calculation of the di†usion coefficients, and were found to be well-behaved once steady-state conditions had been achieved. The results from the benchmark simulation discussed in the previous section are shown in Table 2. The error in the di†usion coefficients was difficult to quantify because it is systematically dependent on the model size and total simulation time. As a guide, the values given were reproducible to two signiÐcant Ðgures. It can be seen that the di†usion coefficients of the oxygen atoms in both the hydroxonium ions (carrying a positive charge ]e) and water molecules (neutral) are approximately equal. This should be compared with the di†usion coefficient of the sulfur atoms, a component of the sulfonate fragments only, which is lower by a factor of 3.6. The di†erence in mobility is not surprising, and can be attributed simply to the smaller size of the water-like moieties. The addition of 10 hydroxyl ions to the benchmark simulation produced a more interesting situation. Due to the negative charge ([e) carried by the hydroxyls, it was necessary to replace 10 water molecules by hydroxoniums in order to ensure overall charge neutrality. The simulations were identical in all other respects, and the resulting di†usion coefficients are listed in Table 3.

Table 2 Di†usion coefficients of atomic species in the benchmark NVT simulation at 300 K for 1 ns Species

Fragment (charge)

D/Ó2 ps~1 atom~1

Sulfur Oxygen Oxygen

Sulfonate ([e) Hydroxonium (]e) Water (O)

5.30 ] 10~5 1.92 ] 10~4 1.91 ] 10~4

As in the benchmark simulation, the mobilities of the hydroxonium and water fragments were comparable in magnitude, and considerably greater than that of the sulfonates. SigniÐcantly, the di†usion coefficient of the hydroxonium ions was greater than that of the hydroxyls by a factor of 1.64. This was somewhat surprising, as the hydroxyl ions are considerably smaller than the hydroxoniums. It seems that the motion of hydroxyl ions was restricted by virtue of their negative charge. In order to test this hypothesis, the equilibrium structure of the models was perturbed, Ðrst, by turning o† the electrostatic interactions, and second, by applying a uniform external electric Ðeld. By comparing the results obtained from simulations with and without electrostatics, it can be proved that the difference in mobility between the hydroxonium and hydroxyl ions is principally due to charged interactions. The di†usion coefficients from simulations including hydroxyl ions, both with and without electrostatic forces, are shown in Table 4 for comparison. The simulations were identical in all other respects. It is immediately clear that the mobility of all species was higher with electrostatics switched o†. On closer inspection, it can be seen that the mobility of fragments without electrostatics varied according to their size, with the ratio of the hydroxonium to hydroxyl di†usion coefficients being 0.61. With the electrostatics switched on, this ratio rose to 1.64, i.e., the hydroxyls went from being more mobile (by virtue of their smaller size) to less mobile (by virtue of their negative charge) than the hydroxonium ions. Interestingly, even the motion of the neutral water molecules was restricted to a greater extent than that of the hydroxoniums when electrostatic interactions were present. In the absence of electrostatics, the simulation does not form a phase-separated morphology. It seems likely, therefore, that there are two mechanisms by which the electrostatic interactions a†ect the di†usion behaviour. First, the presence of electrostatics causes the formation of a phase-separated morphology, and secondly, the action of the charged interactions limits the ability of the ions to move within these morphological constraints. This proposition will be examined in more detail below. The di†usion results constitute clear evidence for selective conduction, i.e., the preferential conveyance of cations over anions through the model. It seems that the mechanism for this ionic selectivity operates at the atomic level, and has been captured by the relatively simple molecular models presented here. In the process, it has been shown that electrostatic interactions are required for the operation of selectivity. The reasons why the motion of negatively charged ions is Table 3 Di†usion coefficients in an NVT simulation containing hydroxyl fragments at 300 K for 1 ns Species

Fragment (charge)

D/Ó2 ps~1 atom~1

Sulfur Oxygen Oxygen Oxygen

Sulfonate ([e) Hydroxonium (]e) Hydroxyl ([e) Water (O)

4.93 ] 10~5 1.00 ] 10~4 6.07 ] 10~5 1.14 ] 10~4

Table 4 Comparison of di†usion coefficients for NVT simulations with and without electrostatic forces, but identical in all other respects, at 300 K for 1 ns D/Ó2 ps~1 atom~1 Species

Fragment (charge)

No electrostatics

Electrostatics

Mobility decrease factor


Sulfonate ([e) Hydroxonium (]e) Hydroxyl ([e) Water (O)

3.05 ] 10~4 2.46 ] 10~3 4.01 ] 10~3 3.74 ] 10~3

4.93 ] 10~5 1.00 ] 10~4 6.07 ] 10~5 1.14 ] 10~4

6.2 24.6 66.1 32.8

4858


restricted will be considered in more detail when discussing a model for ion di†usion in PIMs. The selective conductivity of the molecular models was enhanced in the presence of an external electric Ðeld. Table 5 shows the di†usion coefficients of species in simulations including hydroxyl ions, both with and without an external electric Ðeld. The Ðeld strength used was 4.3 ] 107 V m~1. As would be expected, the external Ðeld caused a rise in the di†usion coefficients of all the species. However, the increase in mobility of the hydroxonium ions was nearly three times that of the hydroxyl ions or water molecules. In addition, the ratio of hydroxonium di†usion to hydroxyl di†usion rose from 1.64, with no Ðeld, to 3.2 when the Ðeld was switched on. It seems that selective conductivity was exaggerated by the presence of an external electric Ðeld, irrespective of any structural disturbances which may have been caused. The origin of selective conductivity in the molecular models can be understood by considering the motion of single charged fragments in the simulations. Fig. 6 shows the MSD of single hydroxonium and hydroxyl ions averaged over 1 ns under the inÑuence of an external electric Ðeld during an NVT simulation at 300 K. The motion of the ions is characterised by long periods of stagnancy, punctuated by relatively abrupt movements. If one of these staccato movements is deÐned as a ““ jump ÏÏ when it is bounded by a plateau on either side, then it is clear that the hydroxonium jumped four times during the simulation, whereas the hydroxyl jumped only once. It is this di†erence between the frequencies of jumping which gives rise to the dramatic di†erence in mobility of the charged fragments when averaged over the entire model. It would appear that the sulfonate fragments play a dominant role in the conduction of ions through the model by providing an electrostatic barrier to the passage of negatively charged fragments. This can be demonstrated by observing the motion of individual fragments when the electrostatic forces are switched o†, as shown in Fig. 7. The punctuated equilibrium behaviour seen in Fig. 6 is replaced by a type of motion which is almost Brownian.

Fig. 6 The MSD of oxygen atoms in single hydroxonium and hydroxyl ions averaged over 1 ns during an NVT simulation at 300 K under the inÑuence of an external electric Ðeld.

Fig. 7 The MSD of oxygen atoms in single hydroxonium and hydroxyl ions during an NVT simulation at 300 K in the absence of electrostatic forces.

An example of the jumping motion of the hydroxonium fragments between sulfonate groups is shown explicitly in Fig. 8. A series of 15 frames are displayed, each separated in time by 0.5 ps, taken from a 1 ns NVT simulation at 300 K with an externally applied electric Ðeld. A hydroxonium fragment, which is initially associated with the sulfonate fragment on the lower right of the frame, moves abruptly to another sulfonate on the upper left of the frame in a time of order 1 ps. Fig. 6 demonstrates that the hydroxyl fragments also show jump diffusion. In this case it is clear that the jumps will not take place between sulfonate groups. In fact, it appears that the jumps occur when the hydroxyls di†use past constrictions in the channels in the phase-separated morphology. A possible alternative explanation for the observed ordering of di†usion coefficients should also be considered. That is, it is possible that the motion of the hydroxyl fragments is restricted by their hydration spheres (which will tend to be larger than those surrounding the hydroxoniums). We investigated this possibility by performing similar simulations to those described above, but omitting the sulfonate fragments. These are e†ectively simulations of ionised water. They show that, without an external Ðeld, the mobilities of the hydroxyl and hydroxonium groups is approximately equal, whereas, when an external Ðeld is applied, the mobility of the hydroxyl fragments is greater than that of the hydroxoniums by a factor of 3. This is the opposite to the ordering found in simulations involving sulfonates, and in line with expectations based on ion size. Thus it is clear, that the anomalous di†usion of the ionic species in the simulations of PIMs, is a consequence of the presence of the sulfonate ions. The jump di†usion behaviour corresponds closely with that predicted by the cluster-network model of Hsu and Gierke25 for NaÐon. They postulated the existence of clusters, within which the movement of ions was e†ectively unrestricted, connected by narrow tubes which presented a signiÐcant electrostatic barrier to di†usion (of anions in particular). Although the model was based on continuum concepts, it seems that the underlying physical intuition has been vindicated by the

Table 5 Comparison of di†usion coefficients in NVT simulations with and without an external electric Ðeld (4.3 ] 107 V m~1), but identical in all other respects, at 300 K for 1 ns D/Ó2 ps~1 atom~1 Species

Fragment (charge)

No Ðeld

External Ðeld

Mobility increase factor


Sulfonate ([e) Hydroxonium (]e) Hydroxyl ([e) Water (0)

4.93 ] 10~5 1.00 ] 10~4 6.07 ] 10~5 1.14 ] 10~4

1.75 ] 10~4 1.63 ] 10~3 3.01 ] 10~4 8.25 ] 10~4

3.5 16.2 5.0 7.2


4859

Fig. 9 Comparison of the cumulative RDF between sulfonate fragments over a 1 ns NVT simulation at 300 K, with and without the use of electrostatic forces.

Fig. 8 A hydroxonium fragment (in black) jumping from one sulfonate fragment (in grey) to another (also in grey) during a 1 ns NVT simulation at 300 K with an externally applied electric Ðeld. The Ðeld was applied in the vertical direction, pointing from the bottom of the page to the top. The sulfur atoms in the sulfonate groups are labelled in the Ðrst and last frames, and the times quoted are relative to the Ðrst frame.

results of atomistic modelling. However, the molecular structure of the percolation network is more subtle than the phenomenological parody of the cluster-network model. The conÐgurational properties of the simulations will now be described in detail.

The RDFs between sulfur atoms in the sulfonate fragments and oxygen atoms in the hydroxonium and hydroxyl ions provided important information about the state of aggregation of ionic material in the simulations. The sulfonateÈ hydroxonium RDF for simulations with and without electrostatics, shown in Fig. 10, will be discussed Ðrst. There are two sharp peaks at 3.3 and 4.2 Ó, which correspond to two distinct hydroxoniumÈsulfonate bound states. These correspond to the ““ bound water ÏÏ molecules detected in IR spectroscopy studies,8,9 which are predominantly hydrogen-bonded to sulfonate groups. Further out, beyond a hydroxonium depleted region, there are a series of di†use peaks between 6 and 10 Ó. These may be interpreted as ““ free water ÏÏ, which is more loosely associated with the matrix. It is clear from the heights of the sharp peaks that there is an enhanced probability of Ðnding a hydroxonium ion bound to a sulfonate group, i.e., in a hydroxoniumÈsulfonate boundstate. Moving on to the sulfonateÈhydroxyl RDFs, shown in Fig. 11, it can be seen that the occupancy of hydroxylÈsulfonate bound states is much lower than that of the hydroxoniumÈ sulfonate states. This would seem to be at variance with the hypothesis proposed in the previous section that hydroxonium fragments are more mobile than hydroxyls. However, it is their greater affinity for sulfonates which enables the hydroxoniums to pass more easily through the model than the hydroxyls. This is because the conductive channels formed by the polar and non-polar material are naturally bounded by a layer of sulfonate groups. Being less able to bind to these, the

Radial distribution functions The normalised radial distribution function (RDF) G (r) is AB deÐned here as the spherically averaged distribution of interatomic vector lengths between two species A and B, totalling N and N , in a unit cell of volume v : A B N (r)v AB G (r) \ (4) AB 4pr2drN N A B where N (r) is the unnormalised RDF. The cut-o† radius was AB set to the linear dimensions of the unit cell in order to exclude self-terms. The cumulative RDFs between sulfonate fragments for simulations with and without electrostatics are shown in Fig. 9. It is clear that the presence of electrostatic interactions gives rise to a preferred spacing of approximately 5.4 Ó between sulfur atoms in adjacent fragments. This e†ect is due to the repulsive forces between the negatively charged sulfonate groups, which are a common feature of phenomenological models for ion clustering. 4860


Fig. 10 Comparison of the cumulative hydroxoniumÈsulfonate RDF over a 1 ns NVT simulation at 300 K, with and without the use of electrostatic forces.

Fig. 11 Comparison of the cumulative hydroxylÈsulfonate RDF over a 1 ns NVT simulation at 300 K, with and without the use of electrostatic forces.

Fig. 13 Comparison of the cumulative RDF between water molecules over a 1 ns NVT simulation at 300 K, with and without the use of electrostatic forces.

hydroxyls are restricted to the amorphous ““ free water ÏÏ phase, thus considerably restricting their passage through the tortuous conduction pathways of the membrane. Having examined the associations between the sulfonate fragments and the charged water-like fragments, the state of aggregation of the ions themselves will now be considered. The RDFs between hydroxoniums in simulations with and without electrostatics are shown in Fig. 12. The reduction of the nearest neighbour peak of the charged simulation in the RDF compared to the uncharged simulation suggests that the hydroxoniums are more likely to be present as ““ bound water ÏÏ in the former. In fact, the peak at 4.2 Ó, together with its shoulder, may be attributed to hydroxoniums which are bound to sulfonate fragments. The ““ bound water ÏÏ peak at 4.2 Ó can be clearly identiÐed if Fig. 12 is compared to the RDF between neutral water fragments, shown in Fig. 13. The similarity of the hydroxonium and water RDFs in the uncharged simulations, contrasted with the di†erence between these two when electrostatic inter-

actions were present, is further evidence for a division between the bound and free water phases. Having shown that electrostatic interactions play a crucial part in the establishment of a realistic ionomer morphology, the e†ect of varying water content on the molecular models will now be described. The intersulfonate RDF is shown in Fig. 14 as a function of water content. The water contents of the models were parametrised from a mass swelling isotherm, and the total number of water-like fragments for each model is given in Table 6. It should be noted that the models with varying water content contained no hydroxyl ions, and a constant number of hydroxonium ions. The size of the simulation box was adjusted in each case to give the correct overall density for the mass-weighted combination of ionomer and water. Fig. 14 shows that as the water content of the models was increased, the mean separation of the sulfonate groups also increased. This was accompanied by the appearance of a nextnearest neighbour peak at just over 10 Ó, indicating a rise in the quality of the spatial ordering.

Fig. 12 Comparison of the cumulative RDF between hydroxonium ions over a 1 ns NVT simulation at 300 K, with and without the use of electrostatic forces.

Fig. 14 Comparison of the cumulative RDF between sulfonate fragments over a 1 ns NVT simulation at 300 K, as a function of water content.

Table 6 Water contents of the molecular models, with their corresponding physical condition, parametrised empirically using mass swelling isotherm data Physical condition

Percentage water by mass

Number of water like fragments

Number of water molecules

Density /g cm~3

Box size (edge length)/A

Dry Ambient 100% RH

1.6 6.3 15.9

26 100 253

0 74 227

2.082 2.031 1.925

20.0 21.0 22.2


4861

The principal intersulfonate spacings for the 100% relative humidity (RH) model, at approximately 6.0, 10.5, and 12.5 Ó are close to the Ðrst, second and third nearest neighbour spacings in a face-centred cubic (f.c.c.) packing arrangement. If we assume this conÐguration of the sulfonate groups is applicable to all three models, then two results emerge. First, the nearest neighbour spacing is not simply related to the dimensions of the simulation box, dispelling any thought that the observed ordering may be a consequence of the small simulation box size. Second, the fraction of the simulation box occupied by the hydrated sulfonate fragments increases with water content, from 21% for the ““ dry ÏÏ system to 35% for the 100% RH system. This result again suggests a preference for the water to be associated with the ionic groups rather than the Ñuorocarbon matrix, as the ionomer swells. It also reinforces the notion that the sulfonate fragments are aggregated to some degree, and that the increase in intersulfonate spacings with water content is not simply due to the increase in simulation box size. To examine further the e†ect of simulation box size, two test simulations were performed, one with, and one without an external Ðeld, using models with eight times the volume of the majority of the models considered here. These simulations conÐrmed the formation of a phase-separated structure in the presence of electrostatic interactions, with a feature size comparable with the smaller simulations. In particular, the Ðne structure in the RDFs above 6 Ó was preserved, as well as the general trend in di†usion coefficients, conÐrming that the trends displayed by the smaller models were representative, and not due to Ðnite size e†ects. There were few other trends which varied systematically with water content, although the results were limited by the size of the models and therefore the number of intersulfonate spacings which could be sampled. Having described the internal structure of the models using RDFs, a more direct way of visualising the conducting channels which have been described in this discussion will now be presented. Free volumes The state of aggregation of water-like fragments in the molecular models was visualised directly by replacing these fragments by a 3-dimensional Connolly surface. A Connolly surface is formed by rolling a spherical probe of given radius over the van der Waals (vdW) surface of the model, i.e., it is the locus of points formed by the intersection of the probe with the vdW surface of the nearest atom.26,27 After the dele-

tion of all water-like fragments from an equilibrated model, the construction of such a surface allows the direct visualisation of conduction pathways formed by the ionic phase. Colouring the surface by its proximity to the remaining atoms permits the species which border these pathways to be identiÐed. In addition, it is possible to calculate the amount of free volume as a function of the water content of the models. The Connolly surfaces from models which were empirically parametrised to ambient and 100% RH using gravimetric swelling data will be examined. The models were equilibrated at 300 K for 1 ns, the water-like fragments removed and the remaining space Ðlled by a Connolly surface using a probe radius of 1.4 Ó. The resulting images are shown in Fig. 15, using a colour scheme that is related to the proximity of the surface to the remaining atoms (which are hidden for the purpose of clarity). Green areas of the surface correspond to areas closest to Ñuorine atoms, red areas to oxygen atoms and yellow areas to sulfur atoms. The conductive channels formed by the water-like fragments can be clearly seen in Fig. 15. At ambient humidity, Fig. 15a, the channels are narrow and relatively discontinuous. Isolated regions, or areas with a predominantly green perimeter, correspond to bottlenecks for ionic di†usion. At 100% RH, Fig. 15b, the channels are virtually continuous and visibly swollen compared to Fig. 15a. The red and green areas can clearly be delineated. In the cluster-network model, red areas would correspond to clusters, and green areas to connective tubes. However, from a molecular perspective, these two intermingle seamlessly to create an intricate network of channels with varying resistivity. The amount of free volume was found to increase linearly as a function of water content. However, in the case of real PIMs, evidence for a non-linear swelling isotherm at high water contents has been obtained from small angle X-ray scattering.10,11,28,29 This would seem to indicate that there is some sort of higher order ionic agglomeration occurring above the level of clustering. It is clear that the molecular models discussed were not sufficiently large to reproduce this feature of ionomer morphology. However, as indicated above, experiments with larger box sizes showed that phase separated structures were still formed, and therefore that those observed in smaller models were not due to Ðnite size e†ects.

Concluding remarks The results from the MD simulations show that it is possible to capture some of the salient dynamical and conÐgurational

Fig. 15 Three-dimensional Connolly surfaces produced from models containing : (a) 6.3% and (b) 15.9% water by mass. These were generated by removing the water-like fragments after NVT equilibration at 300 K for 1 ns. Although the remaining sulfonate fragments have been hidden for clarity, the surface is coloured by its proximity to their atomic components. Green areas are closest to Ñuorine atoms, red areas to oxygen atoms and yellow areas to sulfur atoms.

4862


properties of PIMs in an atomistic model of a tractable size. Despite the use of a suspension of molecular fragments, rather than a fully connected polymer chain, the result is an apparently phase separated morphology which displays the correct type of electrochemical transport behaviour. The internal structure of the simulations is closely analogous to the phenomenological cluster-network of Hsu and Gierke.25 The dimensions of the models (typically 21 Ó cubed) are too small to enable a quantitative comparison with the various small-angle X-ray scattering models which have been proposed for PIMs. These normally involve ionic clusters with an average separation, or diameter, of 35È55 Ó.11 In fact, the small box size, together with the absence of a fully connected polymer chain, is no doubt responsible for the smaller size of the aggregates in the simulations described here. The structure of water in the models is in agreement with results from IR spectroscopy. In particular, there is evidence for a ““ bound water ÏÏ phase which is strongly associated with the sulfonic acid groups, and a more loosely attached ““ free water ÏÏ phase. Selective conductivity of ions through the model has been demonstrated by applying an external electric Ðeld to the simulation. Cation mobilities have been measured which are approximately 3.2 times greater than anion mobilities, for ions of similar size and mass (i.e., H O` and OH~). The di†erence 3 in mobility is attributed to di†erent conduction mechanisms for the two types of ions : the cations move by jump di†usion from one sulfonate group to the next along water-Ðlled connecting channels, whereas the motion of the anions is e†ectively blocked by the presence of the same sulfonate groups. We believe that the work presented here is a signiÐcant Ðrst step in understanding the properties of PIMs at a molecular level, and that the technique could usefully be applied to other ionomeric materials in the future.

Acknowledgements The authors are indebted to National Power plc for full Ðnancial support of this work, which forms part of their ongoing research into fuel cell materials.

References 1 2

M. F. Hoover, J. Polym. Sci. : Polym. Symp., 1974, 45, 1. A. Eisenberg, J. Polym. Sci. : Polym. Symp., 1974, 45, 99.

3 H. L. Yeager and A. Steck, J. Electrochem. Soc., 1981, 128, 1880. 4 H. L. Yeager, B. OÏdell and Z. Twardowski, J. Electrochem. Soc., 1982, 129, 85. 5 N. G. Boyle, V. J. McBrierty and D. C. Douglass, Macromolecules, 1983, 16, 75. 6 N. G. Boyle, V. J. McBrierty and A. Eisenberg, Macromolecules, 1983, 16, 80. 7 S. Schlick, G. Gebel, M. Pineri and F. Volino, Macromolecules, 1991, 24, 3517. 8 S. R. Lowry and K. A. Mauritz, J. Am. Chem. Soc., 1980, 102, 4665. 9 M. Falk, Can. J. Chem., 1980, 58, 1495. 10 E. J. Roche, M. Pineri, R. Duplessix and A. M. Levelut, J. Polym. Sci. : Polym. Phys. Ed., 1981, 19, 1. 11 T. D. Gierke, G. E. Munn and F. C. Wilson, J. Polym. Sci. : Polym. Phys. Ed., 1981, 19, 1687. 12 S. Neyertz and D. Brown, Electrochim. Acta, 1998, 43, 1343. 13 J. Ennari, M. Elomaa and F. Sundholm, Polymer, 1999, 40, 5035. 14 M. A. F. Robertson and H. L. Yeager, in Ionomers, ed. M. R. Tant, K. A. Mauritz and G. L. Wilkes, Chapman and Hall, London, 1997, ch. 7. 15 Cerius2 modelling environment version 3.5, 1995, Molecular Simulations Inc., 240/250 The Quorum, Barnwell Road, Cambridge, CB5 8HE, UK. 16 J. J. P. Stewart, J. Comput.-aided Mol. Des., 1990, 4, 1. 17 W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926. 18 M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, Oxford, 1987. 19 S. L. Mayo, B. D. Olafson and W. A. Goddard, J. Phys. Chem., 1990, 94, 8897. 20 D. B. Holt, B. L. Farmer, K. S. Macturk and R. K. Eby, Polymer, 1996, 37, 1847. 21 A. R. Leach, Molecular modelling : principles and applications, Addison Wesley, Longman, Harlow, Essex, 1996. 22 DL–POLY is a package of molecular simulation routines written by W. Smith and T. R. Forester, copyright The Council for the Central Laboratory of the Research Councils, Daresbury Laboratory at Daresbury, Nr. Warrington, 1996. 23 H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Dinola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684. 24 F. Muller-Plathe, Acta Polym., 1994, 45, 259. 25 W. Y. Hsu and T. D. Gierke, J. Membrane Sci., 1983, 13, 307. 26 M. L. Connolly, Science, 1983, 221, 709. 27 M. L. Connolly, J. Appl. Crystallogr., 1983, 16, 548. 28 M. Fujimura, T. Hashimoto and H. Kawai, Macromolecules, 1981, 14, 5, 1309. 29 J. A. Elliott, PhD Thesis, University of Bristol, 1998.

Paper 9/05267D


4863