Massively Parallel Simulations of Diffusion in Dense Polymeric Structures

Jean-Loup Faulon
Computational Materials Science
Sandia National Laboratories
MS 1111
Albuquerque, New Mexico 87185-1111

J. David Hobbs
Department of Chemistry and Geochemistry
Montana Tech of the University of Montana
1300 West Park Street
Butte, Montana 59701-8997

David M. Ford
Assistant Professor
Department of Chemical Engineering
Texas A&M University
College Station, Texas 77843-3122

Robert T. Wilcox
Computational Materials Science
Sandia National Laboratories
MS 1111
Albuquerque, New Mexico 87185-1111

Abstract:

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.

Keywords:

 

Polymer, Gas Diffusion, Molecular Builder, Molecular Dynamics, Teraflop

 

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³) [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, d = 3 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³. 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:

- [ CH2 - C (CH3) 2 ] - , - [ CH2 - C(CH3) = CH2 ] -

isobutylene (98 %) isoprene (2%)

EPDM:

- [ CH2 - CH2 ] - [ CH2 - CH (CH3) ] - , [ CH2 = CH - CH2 - CH = CH (CH3) ]

ethylene-propylene (97%) 1,4 hexadiene (3%)

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 (full-atom) 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.  

PIB 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³, 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³) 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.

PIB Model B was built at the experimental density of 0.91 g/cm³ 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_BT 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 0.16 k_BT 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)

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³ 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 n = 1. The normal mode is reached for both systems within 1000 picosecond. The diffusion coefficient for model A is 24.09 10^-6 cm²/s and 25.25 10^-6 cm²/s for model B. Both diffusion coefficients are close to the value found by Muller-Plathe et al. [2] (30.1 10^-6 cm²/s) 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 y = x 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³, 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³. 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₂ and water (10^-6cm²/s)

 

Model system	Diffusion coef. (this work)	Diffusion coef. (experiment) [11]
O2 in rubber	0.285	0.081
H2O in rubber	0.159	not reported
O2 in EPDM	0.781	0.177
H2O 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 free 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.

References

[1] Brooks, B.R.; Bruccoleri, R. E.; Olafson, B.D.; States, D.J.; Swaminathan, S. and Karplus, M. "CHARMm: A program for macromolecular energy, minimization and dynamics calculations", J. Comp. Chem. 1983, 4, 187.

[2] Muller-Plathe, F. Rogers, SC., van Gunsteren, W. F. J. Chem. Phys. 1993, 98, 9895.

[3] Plantenga, T. Personal communication 1997.

[4] Plimpton, S. J., Hendrickson, B. A. J. Comput. Chem. 1996.

[5] Plimpton, S. J., Pollock R., Stevens M., "Particle--Mesh Ewald and rRESPA for Parallel Molecular Dynamics Simulations", in Proc. of Eighth SIAM Conf on Parallel Processing for Scientific Computing, Minneapolis, MN, March 1997.

[6] Plimpton, S. J. J. Comput. Physics 1995, 117, 1.

[7] van der Vegt, N., F. A., Briels, W. J., Wessling, M., Strathmann, H. J. Chem. Phys. 1996, 105, 8849.

[8] Volume relaxation times for polymer can be on the order of years (physical aging). Cf. Eisenberg, A. In Physical Properties of Polymers, J. E. Mark, et al. Eds., American Chemical Society, Washington, DC.

[9] Flory, P. J. Principles of Polymer Chemistry, Cornell University Press: Ithaca, New York, 1953.

[10] de Gennes, P. G. Scaling Concepts in Polymer Physics, Cornell University Press: Ithaca, New York, 1979.

[11] Pauly, S. in Polymer Handbook, 3rd ed., edited by J. Branrup and E. H. Immergut, Wiley: New York, 1989.

[12] Faulon, J. L., Mathews, J. P., Carlson, G. A., Hatcher, P. G. Energy & Fuels 1994, 8, 408.

[13] Takeuchi, H., Okazaki, K. J. Chem. Phys. 1990, 92, 5643. Takeuchi, ibid. 1990, 93, 2062. Takeuchi, ibid. 1990, 93, 4490.

[14] Muller-Plathe, F. J. Chem. Phys. 1991, 94, 3192.

[15] Sok, R. M., Berendsen, H. J. C., van Gunsteren, W. F. J. Chem. Phys. 1992, 96, 4699.