Massively parallel simulations of diffusion in dense polymeric structures.
An original computational technique to generate close-to-equilibrium dense polymeric structures is proposed. Diffusion of small gases are studied on the equilibrated structures using massively parallel molecular dynamics simulations running on the Intel Teraflops (9216 Pentium Pro processors) and Intel Paragon (1840 processors). Compared to the current state-of-the-art equilibration methods this new technique appears to be faster by some orders of magnitude. The main advantage of the technique is that one can circumvent the bottlenecks in configuration space that inhibit relaxation in molecular dynamics simulations. The technique is based on the fact that tetravalent atoms (such as carbon and silicon) fit in the center of a regular tetrahedron and that regular tetrahedrons can be used to mesh the three-dimensional space. Thus, the problem of polymer equilibration described by continuous equations in molecular dynamics is reduced to a discrete problem where solutions are approximated by simple algorithms. Practical modeling applications include the construction of butyl rubber and ethylene-propylene-dimer-monomer (EPDM) models for oxygen and water diffusion calculations. Butyl and EPDM are used in O-ring systems and serve as sealing joints in many manufactured objects. Diffusion coefficients of small gases have been measured experimentally on both polymeric systems, and in general the diffusion coefficients in EPDM are an order of magnitude larger than in butyl. In order to better understand the diffusion phenomena, 10,000 atoms models were generated and equilibrated for butyl and EPDM. The models were submitted to a massively parallel molecular dynamics simulation to monitor the trajectories of the diffusing species. The massively parallel molecular dynamics code used in this paper achieves parallelism by a spatial-decomposition of the workload which enables it to run large problems in a scalable way where both memory cost and per-timestep execution speed scale linearly with the number of atoms being simulated. It runs efficiently on several parallel platforms, including the Intel Teraflops at Sandia. There are several diffusion modes observed depending if the diffusion is probed at short time scale (anomalous mode) or long time scale (normal mode). Ultimately, the diffusion coefficient that need to be compared with experimental data corresponds to the normal mode. The dynamics trajectories obtained with butyl and EPDM demonstrated that the normal mode was reached for diffusion within one nanosecond of simulation. In agreement with experimental evidences, the oxygen and water diffusion coefficients were found larger for EPDM than butyl.
Introduction
Many technological processes depend on the design of polymers with desired permeation characteristics of small molecules.
Examples include gas separation with polymeric membranes, food packaging, and encapsulant of electronic components in
polymers that act as barriers to atmospheric gases and moisture. To solve the polymer design problem successfully, one needs
to relate the chemical composition of the polymer to the diffusivities of the the penetrants molecules within it. One way of
establishing structure-property relationships is the development of empirical relationships (namely QSPRs - quantitative
structure property relationships). A more sophisticated approach is the application of simulation techniques that rely directly
upon fundamental molecular science. While QSPRs can generally be established from the the monomers configurations,
simulations involve the generation of configurations representing the whole polymer. From these configurations, structural,
thermodynamics, and transport properties are estimated. In principle, simulations can provide exact results for a given model
representation of the polymer/penetrant system. In practice, computer time consideration necessitates the introduction of
approximation, and particularly for polymeric systems, the use of high performance computing.
In this paper, we propose an original technique to generate close-to-equilibrium polymeric structures. We are presenting data
that suggest our technique is some orders of magnitude faster than the current state-of-the-art methods used to prepare and
equilibrate dense polymeric systems. The technique is used to generate initial model structures for polymers present in O-ring
systems. Oxygen and water diffusions though the O-ring are then probed using massively parallel molecular dynamics.
Results are compared with similar simulations and experimental data.
Methodology
As pointed out in the introduction, the ability to represent the molecular level structure and mobility of polymeric structures
is a prerequisite for simulating diffusion in them. One fundamental parameter characterizing molecular structures is the
potential energy. Polymer chains in a equilibrium melt or amorphous glass remains essentially at minimum energy
configurations. The potential energy is the sum of bond and bond angle distortion terms, tortional potentials, as well as
intermolecular and intramolecular repulsive and attractive interactions (i.e., van der Waals), and electrostatic (i.e.,
Coulombic) interactions. These energy terms are expressed as a function of the atoms coordinates constituting the polymer
and a set of parameters computed from experimental data or quantum mechanics calculations. The functional forms of the
energy terms and their associated parameters is called a forcefield. In the present study, we are making use of the
commercially available CHARMm forcefield [1], as well as forcefield parameters taken from Muller-Plathe et al. [2]
Molecular mechanics and molecular dynamics. Molecular mechanics is the procedure by which one locates local minima
of energy. Molecular mechanics simply consists of a minimization routine (conjugate gradient, steepest descent, etc.) that
finds the first energy minimum from a starting configuration. In the present paper, a parallel implementation of conjugate
gradient is utilized [3]. The main disadvantage of molecular mechanics is that thermal fluctuations are not explicitly taken
into account, and therefore cannot be used for diffusion calculations. Moreover, the procedure followed in generating
minimum energy configurations does not correspond to any physical process of polymer formation. Nonetheless, static
minimum energy structures provide satisfactory starting configurations for molecular dynamics simulations.
Molecular dynamics (MD) follows the temporal evolution of a microscopic model system through numerical integration of
the equations of motions for all the degrees of freedom. MD simulations can be performed in the microcanonical ensemble
(constant number of molecules, volume, and total energy, or NVE), the canonical ensemble (constant number of molecules,
volume, and temperature, or NVT), as well as the isothermal-isobaric ensemble (constant number of molecules, pressure, and
temperature, or NPT). The major advantage of MD is that it provides detailed information of short-time molecular motions.
Its limitation resides in computer time consideration. Hundreds of CPU hours on a vector supercomputer are required to
simulate a nanosecond of actual atomistic motions. However, computer time can be decreased by making use of massively
parallel processing. In the present work we are using a large-scale atomic/molecular massively parallel simulator (LAMMPS)
[4]. LAMMPS is a new parallel MD code suitable for modeling large molecular systems. LAMMPS has been written as part
of a CRADA (Cooperative Research and Development Agreement) between two DOE labs: Sandia and Lawrence Livermore,
and three industrial partners: Bristol-Myers Squibb, Cray Research, and Dupont. LAMMPS is capable of modeling a variety
of molecular systems such as bio-membranes, polymers, liquid-crystals, and zeolites. The code computes two kinds of forces:
(1) short-range forces such as those due to van der Waals interactions and molecular bond stretching, bending, and torsions,
and (2) long-range forces due to Coulombic effects. In the latter case, LAMMPS uses either Ewald or particle
particle/particle-mesh (PPPM) techniques to speed the calculation [5].
LAMMPS achieves parallelism by a spatial-decomposition of the workload which enables it to run large problems in a
scalable way where both memory cost and per-timestep execution speed scale linearly with the number of atoms being
simulated [6]. It runs efficiently on several parallel platforms, including the Intel Teraflops at Sandia (9216 processors), the
large Intel Paragon at Sandia (1840 processors) and large Cray T3D machines at Cray Research.
Generate initial structure. Creating an atomic-level model of a completely equilibrated dense polymer melt is a challenging
task, given current computational capabilities. van der Vegt et al. [7] discuss the various approaches to this problem. The
common techniques used in simulations of liquids which consists of melting an idealized structure generally takes tens of
picosecond of equilibration using MD. Unfortunately, the equilibration time for dense polymers is many orders of magnitude
larger than are feasible with MD [8]. Consequently, one needs other efficient ways of preparing initial polymer structures,
that resemble the equilibration structure. One of the most common methods has been to pack chains into the simulation box
at the experimental density, either by random placement followed by energy minimization or with Monte Carlo chain growth
techniques. However, these methods have tended to produce structures which are rather inhomogeneous and which lead to an
overestimation of solubility values for small permeants. In the present paper, we are proposing and comparing two alternative
construction methods: the compression box technique and the lattice construction technique.
Compression box technique. van der Vegt et al. [7] suggest that a more efficient way of producing near-equilibrium
structures is by starting with a dilute model polymer and compressing it slowly until the target experimental density is
reached. A set of chains is built with the correct dihedral angle distributions, at about 1/8 of the experimental density. Using
an NPT MD technique with a pressure ramp, the model system is compressed to the desired density over ~ 500 picoseconds.
During this compression, only the repulsive part of the nonbond (van der Waals) interaction is used, to avoid "clustering" of
the polymer. Then the structure is further equilibrated with the full nonbond interactions at the new density, for ~ 1000
picoseconds. van der Vegt et al. [7] observed that a model of poly(dimethylsiloxane) built with this "slow-compression"
technique had significantly less stress in the dihedral angles than a model built with a "random packing" method.
Furthermore, the slowly compressed model yielded small-permeant solubility results which were in much better agreement
with experiment.
Lattice construction technique. Many natural and synthetic polymers including all structures considered in this paper are
formed with hydrocarbon chains. More precisely, these polymers are composed of linear chains of tetravalent carbon atoms
to which are attached molecular groups (methyl, phenyl,. The lattice technique is based on the fact that any tetravalent
atom (such as carbon and silicon) fits in the center of a regular tetrahedron where its four bonds are perpendicular to the faces
of the tetrahedron (cf. Figure 1a). It is well known that regular tetrahedrons can be used to mesh the three-dimensional space.

Figure

1. Tetrahedron meshing. a) Carbon atom in center of tetrahedron. b) Two dimensional projection of meshed
space.
Using the previous observation, the experimental polymer volume is meshed and the polymer chains are constructed by
generating random walks following adjacent tetrahedrons (cf. Figure 1.b). One may notice from Figure 1 that the distance
between two adjacent tetrahedrons is equal to a carbon-carbon bond length, and the angle between three adjacent tetrahedrons
is equal to a carbon-carbon-carbon angle. Therefore, bond and angle energy terms are directly minimized by the lattice
construction procedure. Excluded volume interactions are treated by keeping track of the occupancy of each tetrahedron. As
the chain construction progresses the next atom position is always chosen in a non occupied tetrahedron. In the unprobable
event where all adjacent tetrahedrons are already occupied the construction routine back-tracks to the previous position. In
order to keep homogeneous density the first atom of the chain is chosen at random in the least occupied region of the volume
box. The goal of above excluded volume procedure is to keep the intermolecular and intramolecular repulsive forces to a
minimum value.
The lattice construction technique can be used with cubic cells instead of tetrahedrons. In such a case, each cell represents a
monomer, and the cell length, width, and height are those of the boundary box of the monomer. The advantage of the cubic
cells lattice is a reduction of the computational complexity since the objects manipulated are monomers instead of atoms (for
instance, for butyl, there are 12 atoms per monomer). The disadvantage of the cubic lattice is that bond lengths and bond
angles between monomers are no longer necessarily valid. Therefore, structures generated with cubic lattices must be energy
minimized prior using MD simulations in order to restore the correct bond lengths and bond angles.
An important parameter characterizing equilibrium melts or amorphous glasses is the mean square end-to-end distance of the
polymer chains. According to Flory "random coil hypothesis", for which there is ample experimental support [9], at
equilibrium the mean square end-to-end distance of any polymer chain in the bulk is related to the number n of skeletal bond
lengths l as , where is the characteristic ratio of the polymer. It is also known that the correct end-to-end
distance starting from any random initial configuration can be reached by MD using a simulation time proportional to O(n 3 )
[10]. Our lattice construction technique makes use of Flory's result while avoiding the computational complexity of MD
simulations. For each chain to be constructed the initial atom position is chosen at random (in the least occupied region) and
the final atom position is chosen within a distance equal to the value given by the Flory equation. Then, a path is constructed
between the two chosen positions. At each step of the construction the non-occupied adjacent tetrahedrons are ranked using
their respective distance to the final position. The tetrahedron or cubic cell corresponding to the shortest distance is chosen as
the next position. When the last position is reached, if the path length is greater than n, then the chain is deleted and two new
positions are chosen. Note that this situation is unlikely to occur . However, in the likely event where the path
length is smaller than n, additional atoms are added to the chain until the correct length is reached. As illustrated in Figure 2,
the addition of new atoms is carried out by deleting at random a bond along the chain, thus creating two non-bonded atoms
from which the chain is extended. Thus, the chain extension is grown using the same path construction as above with the
exception that the initial and final positions are not chosen at random but are the location of two non-bonded atoms.

Figure

2. Increase of polymer chain length. a) Bond i-j is deleted. b) Chain extension is carried out by growing a new
chain from position i and j.
Note that the lattice construction procedure allows one to construct cross-linked polymers since the basic operation is to
generate a chain between two given positions. Hence, cross-linking chains are constructed by generating chains between pairs
of branch-points.
Once the polymeric structures have been constructed, minimization followed by NVT MD simulations are used to further
equilibrate the structures. While minimization is carried out using the aforementioned conjugate gradient algorithm, MD is
performed using the massively parallel LAMMPS code.
Diffusion calculations. Once polymeric structures have been generated and equilibrated using the compression box or lattice
construction technique, diffusion calculations can be carried out. Microscopically, the diffusion coefficient can be calculated
from the motion of the diffusing particles, provided they have been traced long enough so that they perform random walks.
There are several formulations to derive coefficients of diffusion. With MD simulations the most often used expression is the
so called Einstein relation
where D is the diffusion coefficient, is the number of spatial dimension, t is the time, and r(t) is the position vector of a
particle at time t. The angle brackets denote the ensemble average, which in MD simulations is realized by averaging over all
particles of the diffusing species. The calculation of diffusion coefficients rests on the fact that for sufficiently long times the
mean-square displacement of a diffusing particle increases linearly with time. There are however, cases in which the
mean-square displacement is not linearly proportional to time, but obeys a different power law
where n is lower than 1 (normal diffusion). This process is called anomalous diffusion. It is caused by some mechanism
which forces the particles onto a path that is not a random walk. These can be obstacles in the way of the diffusant and thus
inhibit random motion, such as to force the diffusant to remain inside a cavity. Anomalous diffusion persists only on short
time scales. At long enough time scales (and hence length scales) the trajectory of the diffusant becomes randomized and a
change to normal diffusion occurs.
Results and discussion
The main goal of our study is to measure oxygen and water diffusion coefficient in two polymeric materials, butyl rubber and
ethylene-propylene-dimer-monomer (EPDM). Butyl and EPDM are used in O-ring systems and serve as sealing joints in
many manufactured objects. Butyl rubber is a copolymer of isobutene and isoprene (Figure 3). EPDM is a terpolymer of
ethylene, propylene, and 1,4 hexadiene (Figure 3). The experimental density for both butyl and EPDM is 0.91 g/cm 3 . Small
gases diffusion coefficients have been measured experimentally on both polymeric systems, and in
general the diffusion coefficients in EPDM are an order of magnitude larger than in butyl. Presently this
decrepency is poorly understood.
Butyl rubber:
EPDM:

Figure

3. Butyl rubber and EPDM structural formulas.
The results outlined in this section were obtained using the lattice construction code and the energy minimization program both running on a SGI
R10000 platform. MD simulations were carried out using the MD LAMMPS code running on both the 9216 processors Teraflops Intel and 1840
processors Intel Paragon. All butyl MD runs were carried out on the Paragon using 216 nodes, EPDM and polyisobutylene MD runs were
performed on the Teraflops using 125 nodes. We observed a speed up factor of 7 between the Teraflops and the Paragon for the same number of
nodes.
Polymer construction. In prior calculations of diffusion coefficients we probed the performances of the lattice construction code for dense polymeric
systems. It has been suggested by van der Vegt et al. [7] that polymer building techniques which create models directly at the experimental density
may introduce significant nonequilibrium strain in the angle bending and dihedral torsion degrees of freedom. Furthermore, the time required to
anneal this strain may be beyond MD time scales, thus prohibiting the generation of true equilibrium polymer structures. van der Vegt et al. [7]
showed that model polymers built by slow compression from a low density initial state exhibit less strain in the angle and dihedral degrees of
freedom, and they suggest that these structures may be more representative of an equilibrated polymer. Since the lattice construction technique we
introduce here is novel, and it creates models directly at experimental densities, we thought it appropriate to compare it to the "slow-compression"
building method. The polymer used in this study was poly(isobutylene), or PIB, this polymer is enssentially an non cross-linked butyl rubber. The
forcefield parameters were taken from Muller-Plathe et al. [2] The Lorenz-Berthelot mixing rules were used for nonbond interactions
between different atomic species, and all studies were performed at 300K.
Model A was built with the slow-compression method described in the Methodology section. First, a set of polymer chains was built at a density
of 0.11 g/cm 3 , which is approximately 1/8 of the experimental density. There were 21 chains with
lengths distributed randomly between 80 and 120 monomer units. The full-atom model consisted
of 25104 atoms, 25076 bonds, 50040 angles, and 74594 dihedrals. To anneal out nonbond (van der
Waals) overlaps in the initial configuration, an MD simulation with a small timestep was run for a
few picoseconds (ps). We then performed two MD stages as suggested by van der Vegt et al. [7]
The first stage was a compression using only the repulsive part of the nonbond potential (i.e., the
Lennard-Jones interaction between species i and j was truncated at ij ) but full bonded
interactions (bonds, angles, dihedrals). The model was compressed from its initial density to the
experimental density, 0.91 g/cc, over 525 picoseconds, using constant-temperature MD with a
pressure ramp. The Nose-Hoover method was used to control both temperature and pressure,
with time constants of 0.1 and 0.5 ps, respectively. The second stage was initiated from the last
configuration of the first stage, with the attractive nonbond interactions turned on by increasing
the Lennard-Jones cutoffs to 74.5 nm. A constant pressure MD simulation was performed such
that the experimental density (0.91 g/cm 3 ) was maintained. This run lasted for about 120
picoseconds. About 20 picoseconds into the second stage, the polymer seemed to be equilibrated;
changes in the various components of potential energy were not significant or systematic.
Model B was built at the experimental density of 0.91 g/cm 3 using the lattice construction
technique described earlier. It was created to be close in size to Sample A, although it ended up
2.8% smaller. It contained 21 chains with lengths distributed randomly between 80 and 120
monomer units. It was a full-atom model consisting of 24408 atoms, 24387 bonds, 48690 angles,
and 72658 dihedrals. Prior to running dynamic simulations, this model was minimized to reduce
nonbond overlaps using the procedure described in the Methodology section. A MD run at
constant volume and temperature was then performed for 380 picoseconds. The Nose-Hoover
method with a time constant of 0.1 ps was used to maintain temperature at 300K. Equilibration
was reached within the first 10 picoseconds.
Our first comparison of the two polymer models (i.e, PIB A and B) is given in Table I, which
shows the various contributions to potential energy in the systems. The averages and standard
deviations were taken over the last 10 ps of the dynamics for Model A and for the last 80
picoseconds of the dynamics for Model B. All values represent total amounts for each
contribution; since the two models were approximately the same size, we did not normalize the
values. Interestingly, we observe differences between the two models in the van der Waals, angle,
and dihedral contributions; in each case the model created with the lattice contruction technique
(Model B) has a higher energy. The total potential energy difference between the models is 8650
kJ/mol, or 0.14 k B T per atom. The largest differences are found in the angle and dihedral energies,
accounting for about 80% of the total difference. Both of these observations are consistent with those
of van der Vegt et al. [7] for poly(dimethylsiloxane) polymer models. They observed a total difference
of per atom between models made by direct packing and slow compression methods, with
about 70% of that due to angle, dihedral, 1-4 nonbond, and 1-5 nonbond contributions.

Table

I. Potential energy contributions in PIB models (kcal/mol)
PIB Model van der Waals bonds angles dihedrals
Since the concern of this paper is the calculation of diffusion coefficients of small molecules in
polymers, we are interested in how the differences between Models A and B will affect such
calculations. We studied the diffusion of helium (He) through each of the polymer models, using the
Lennard-Jones parameters given by Muller-Plathe et al. [2] for He. In the case of Sample A, a
polymer configuration with a density of exactly 0.91 g/cm 3 was taken from the end of the second-stage runs. Then 100 He
atoms were added to the configuration. In the case of sample B, 100 He atoms were added to the final configuration. The overlaps of the new He atoms
with the polymer were relaxed using constant volume and temperature MD with a small time step. Diffusion MD runs were then performed with constant
volume and temperature (300K) for 360 picoseconds for model A and 770 picoseconds for model B. Figure 4 visualizes the mean-square displacement of
He molecules in samples A and B. The straight lines are drawn with slopes of exactly unity and one can see that at long times eq. 2 is valid for 1. The
normal mode is reached for both systems within 1000 picosecond. The diffusion coefficient for model A is
model B. Both diffusion coefficients are close to the value found by Muller-Plathe et al. this is due to the fact we are using the
same forcefield. Most importantly the diffusion coefficients differences between systems A and B are not significant. Hence, the lattice construction
technique appears to generate equilibrated polymer structures that have the same behavior that those created by the slow compression method in so far
as diffusion calculations are concerned.

Figure

4. Trajectories of helium particles versus time in PIB models A & B (cf. text). The straight line represents the normal diffusion regime.
We further probed the performances of the lattice construction technique with PIB systems of increasing sizes. As shown in Figure 5, the computational
complexity appears to scale linearly with the number of atoms. This is a substantial gain compare to other MD techniques in which the number of steps
of the simulation must be at least n 3 , where n is the number of atoms of the polymer chains (cf. Methodology section). As a consequence, we were able to
generate PIB structures up to 1,380,000 atoms. To the best of our knowledge, these models are the largest non-crystalline bonded atomic systems ever
generated and in general are several orders of magnitude larger than the current models used in polymer science. It is also important to note that the
structures generated by the lattice construction program are equilibrated using, at most, ten picoseconds MD simulations, rather than several hundred
picoseconds with other techniques.

Figure

5. Computational running time of the lattice construction code versus number of atoms and number of monomers.
Oxygen and water diffusion calculations. The lattice construction code was used to generate initial structures for butyl and EPDM. Both structures
contained approximately 10,000 atoms and were generated in a cubic box of 4.5 nm size, which lead to a density of 0.91 g/cm 3 . In both systems, oxygen
molecules were added up to 3% of the total weight, while water molecules were added up to 9% weight. Although 3% weight is the experimental value for
both oxygen and water in butyl and EPDM, a higher value for water was chosen to obtain a better statistical average when calculating diffusion
coefficients. Once the diffusing species were added, the two resulting structures were energy minimized and equilibrated with MD simulations. The
simulation time used for equilibration did not exceed 10 picoseconds.
In order to compute diffusion coefficients using eq. 1, MD was run up to a 1000 picoseconds simulation. This large simulation time was chosen with the
hope of being able to reach the normal diffusion mode (cf. eq. 2 and discussion below). Figures 6 and 7 visualize the mean-square displacement of
oxygen and water molecules in butyl and EPDM as a function of time. With the exception of water in EPDM the normal mode was reached in all cases.
The diffusion coefficients D were evaluated from the positions of these lines.
Figure 6a. Trajectories of oxygen particles versus time in butyl rubber. 6b. Trajectories of oxygen particles versus time in EPDM. The straight line
represents the normal diffusion regime.
Figure 7a. Trajectories of water particles versus time in butyl rubber. 7b. Trajectories of water particles versus time in EPDM.
The non convergence for diffusion calculations of water in EPDM may be attributed to the fact that several diffusing molecules can eventually fit into
EPDM cavities. Indeed, when visualizing the solvated EPDM model it was observed that some cavities contained several water molecules. Because of the
electrostatic attractions between water molecules, one can hypothesize that when several penetrants are present in one cavity, it is energetically favorable
for the penetrants to remain inside the cavity rather than jumping to neighboring cavities. Since 9% weight for water is overestimated compared to the
corresponding experimental value, a new EPDM model was constructed containing only 3% water. With this new model, the normal mode was reached
in less than 1000 picoseconds. The diffusion coefficient values are listed in Table II.

Table

II. Diffusion coefficients for O 2 and water (10 -6 cm 2 /s)
Model system Diffusion coef. (this work) Diffusion coef. (experiment) [11]
O 2 in rubber 0.285 0.081
O in rubber 0.159 not reported
O 2 in EPDM 0.781 0.177
O in EPDM 1.921 not reported
The diffusion coefficients for EPDM are larger than for butyl as expected from experiments. A tentative explanation was provided when visualizing the
volumes of the equilibrated polymers. Free volumes are void spaces in model structures that are accessible by penetrant molecules. The free
volumes were computed using a program developed by one of the authors of this paper. [12] The determination of the free volume within a given model
system is conducted as follows. The radii of all the polymer atoms are augmented by the length equal to the penetrant radius, and the unoccupied volume
of the resulting model system is then calculated. For both butyl and EPDM it was observed that the free volumes did not percolate for penetrant having a
diameter greater 0.30 nm (diameter of atomic oxygen is 0.35 nm). In other words, all cavities having an entrance size greater than 0.30 nm were not
connected. This observation is consistent with the hoping mechanism picture proposed by some authors. [13][14][15] Penetrant molecules spend
relatively long times in cavities before performing infrequent jumps between adjacent cavities. The same authors have performed visual inspection of
polymer models in the vicinity of penetrants, and it appeared that jumping events occur after channels between neighboring cavities are formed. Once
the channels are formed, the penetrants slip through it without much effort. If such a picture is true, since EPDM has higher diffusion coefficients than
butyl, EPDM should have a greater number of channels than butyl. We probed the free volume distribution for butyl and EPDM versus the cavity
entrance size. As shown in Figure 8, butyl and EPDM have a different free volume distribution. EPDM contains significantly more free volumes having
an entrance size smaller than 0.30 nm than butyl. It is important to note that these small free volumes provide links between larger cavities since the free
volume network percolates for entrance sizes smaller than 0.30 nm. Consequently, the curves presented in Figure 8 are consistent with the hoping
mechanism picture, and according to this picture, EPDM has higher diffusion coefficients because it comprises more channels between cavities.

Figure

8. Free volume distribution in butyl rubber and EPDM.
Conclusion
We have proposed a new technique to generate close-to-equilibrium dense polymeric structures. This technique leads to diffusion results that are similar
to results obtained by previously published methods, while being several order of magnitude faster. We have shown the new technique to scale linearly
with the number of atoms, and have used it to generate dense polymer models comprising up to 1,380,000 atoms. To the best of our knowledge these
models are the largest non-crystalline bonded atomic systems ever generated and are many times larger than the current models used in polymer science.
We have used the new technique to construct 10,000 atoms models of butyl rubber and EPDM for the purpose of simulating the diffusion of oxygen and
water molecules. The simulations were carried out using LAMMPS molecular dynamics code, and were run on Sandia's Intel Teraflop and Sandia's
Intel Paragon. In agreement with experimental results the diffusion coefficients in EPDM were found an order of magnitude larger than in butyl. A
tentative explanation of the differences between the coefficients was advanced when comparing the free volume distributions of the two polymer models:
diffusion is facilitated in EPDM due to a larger number of channels between cavities than in butyl.

Acknowledgments

We are pleased to acknowledge the funding provided by the Accelerated Strategic Computing Initiative (ASCI) of the U.S. Department of Energy, Sandia
National Laboratories under contract DE-AC04-76DP00789.



--R

"CHARMm: A program for macromolecular energy, minimization and dynamics calculations"

Personal communication

"Particle-Mesh Ewald and rRESPA for Parallel Molecular Dynamics Simulations"



Principles of Polymer Chemistry
Scaling Concepts in Polymer Physics
in Polymer Handbook




--TR
gas diffusion;teraflop;molecular builder;molecular dynamics;polymer