Detailed Atomistic Simulation of the Segmental Dynamics and

2978

Macromolecules 2004, 37, 2978-2995

Detailed Atomistic Simulation of the Segmental Dynamics and Barrier Properties of Amorphous Poly(ethylene terephthalate) and Poly(ethylene isophthalate) Nikos Ch. Karayiannis,†,‡ Vlasis G. Mavrantzas,†,‡ and Doros N. Theodorou*,†,§ Institute of Chemical Engineering and High Temperature Chemical Processes (ICE/HT-FORTH) and Department of Chemical Engineering, University of Patras, Patras GR 26504, and Department of Materials Science and Engineering, School of Chemical Engineering, National Technical University of Athens, 9 Heroon Polytechniou Street, Zografou Campus, Athens GR 15780, Greece Received August 25, 2003; Revised Manuscript Received February 9, 2004

ABSTRACT: We present results from detailed atomistic simulations concerning the structural, conformational, dynamic, and barrier properties of the amorphous (glassy and melt) phases of two polyisomers, PET [poly(ethylene terephthalate)] and PEI [poly(ethylene isophthalate)]. First, well-relaxed atomistic configurations of the two polyesters are generated following a succession of equilibration stages consisting of energy minimizations, temperature annealings and coolings, and compressions-decompressions, as proposed by Hofmann et al. (Hofmann D.; et al. Macromol. Theory Simul. 2000, 9, 293). This equilibration cycle is significantly extended here by subjecting the resulting configurations to an additional molecular dynamics (MD) simulation at a high temperature for 2 ns, affording extra relaxation at all length scales. With the statistically uncorrelated configurations generated via these extended equilibration runs, isothermal-isobaric (NPT) at P ) 1 atm as well as canonical (NVT) (with the densities set at the corresponding experimentally measured values) MD simulations are performed in the melt state at T ) 450 and 600 K to probe differences in the structure and segmental dynamics of the two polyesters. The simulations reveal significant differences in the local relaxation dynamics of the phenyl rings between the two polyisomers: In PET, these rings are found to exhibit significantly higher mobility than in PEI; this is attributed to the way in which phenyls are connected to their adjacent ester groups. Structural and conformational properties, on the other hand, are predicted to be practically identical in the two polyesters. Finally, transition-state theory (Gusev, A. A.; Suter, U. W. J. Chem. Phys. 1993, 99, 2228) is employed to calculate the rate constants of diffusive jumps between sorption sites, and hence the lowconcentration self-diffusivity, of oxygen (O2) molecules in well-relaxed atomistic configurations of the two polymers, whose densities exactly match the experimentally measured values. Calculated O2 diffusivities are found to be in excellent qualitative and quantitative agreement with experiment: O2 diffusivity in PEI is predicted to be 1.8 times smaller than in PET, in agreement with the experimental finding that PEI is 2-2.5 times less permeable to O2 than PET in the amorphous state.

1. Introduction Poly(ethylene terephthalate) (PET) is extensively used nowadays as the base packaging resin in the bottle manufacturing and food packaging industries, due to its superior barrier properties in comparison with polyolefins, polycarbonates, polystyrene, and other polymers. PET and related alkylene or isomeric polyesters meet the rapidly growing commercial demands for the design of polymeric materials with specified barrier properties, driven by the trend to replace traditional materials, such as glass and metal, with thermoplastics. In general, plastic packaging materials are less brittle and lighter and offer a higher degree of design flexibility than glass and metal. On the other hand, no plastic material is perfectly impermeable: sorption and diffusion of small molecules are possible, to certain extents, in any polymer. Thus, many technologically important processes rely on the design of polymers with tailored characteristics of permeability and selectivity toward fluid molecules. To solve this design problem success* To whom correspondence should be addressed at the National Technical University of Athens. Phone: +30-210-772-3157. Fax: +30-210-772-3112. E-mail: [email protected]. † ICE/HT-FORTH. ‡ University of Patras. § National Technical University of Athens.

fully, one needs to elucidate the molecular-level mechanisms that govern small-molecule transport in a polymer material and how it is related to the chemical constitution, molecular architecture, and processing history of the material. Minimizing permeability (i.e., maximizing the barrier properties) without significantly increasing the production cost of a bottle-grade resin is intimately related to the commercial value of the contained product: In certain cases it can result in safer, more nutritious foods; in other cases it can substantially increase the maximal shelf life of the contained substances with favorable environmental consequences. End-product PET, processed under a very large spectrum of conditions, can exist in many different forms.1-5 It assumes an essentially amorphous, clear form when quenched from the melt. It undergoes thermal crystallization when heated above the glass transition temperature or cooled slowly from the molten state during various stages of container manufacturing; this is generally undesirable, because it reduces clarity. During the stretch-blowing step in bottle manufacturing by blow molding of injection-molded preforms, the preforms are subjected to rapid deformation under high strain rates, which results in the development of crystallinity by a strain-induced mechanism. Unlike thermal crystallization, strain-induced crystallization is

10.1021/ma0352577 CCC: $27.50 © 2004 American Chemical Society Published on Web 03/26/2004

Macromolecules, Vol. 37, No. 8, 2004

Figure 1. Schematic of the PET (left) and PEI (right) repeat units.

not accompanied by haze. PET also exhibits controllable permeability; its barrier properties are relatively low (permeable) in the amorphous state and become drastically higher (less permeable) as the degree of crystallinity increases.6,7 To allow the fabricator more control over thermal crystallization during preform injection or stretch-blow molding, PET bottle-grade resins are modified by isophthalic acid (IA).8,9 Typical IA substitution levels are on the order of 2-4%. IA incorporation widens the process window of PET, resulting in clearer, brighter bottles without significantly affecting the capability of the resin to undergo strain-induced crystallization. Alternatively, taking advantage of the flexibility of stepgrowth condensation polymerizations, the para-, meta-, and ortho-isomers of the constituent diacid can be used to form commercial products: PET, poly(ethylene isophthalate) or PEI, and poly(ethylene phthalate) or PEP, respectively. The repeat units of PET and PEI are compared in Figure 1. Recently, poly(ethylene-2,6naphthalene) or PEN, which is structurally similar to PET, has become competitive to PET in certain performance-driven markets on the basis of its superior strength, heat stability, and barrier properties.8,10-12 Differences in the physical properties of PET and PEN have generally been attributed to the assumed increased rigidity conferred upon PEN by its constituent naphthyl rings, which are clearly larger than the phenyl rings in PET. To overcome the higher cost associated with the production and use of PEN,13 industrial strategies nowadays favor more and more the use of random copolymers or blends based on PET,14,15 such as PETPEI (or PETI)2,7-9,13-16 and PET-PEN,8,11,13,17 with the second component (PEI or PEN) being present in relatively small concentrations. PET and PEI are distinguished solely by the geometric connectivity of ethylene glycol fragments to their phenyl rings (para and meta, respectively). The improved barrier properties of PETI and PET-PEN copolymers have been demonstrated by a large number of experimental investigations;2,7-9,13-16,18-22 in addition, it has been observed that the random incorporation of the comonomer leads to significant changes in the crystallization behavior of the copolyester. This happens because of the increased number of conformational defects along the polymer chain, which prevent copolymer crystallization (random copolyesters exhibit lower melting points and reduced thermal crystallization properties compared to the pure PET form). PETI packaging resins contain up to 2-5% PEI per mole, with an upper limit of about 40%. (Larger PEI concentrations lead to purely amorphous products with relatively poor mechanical properties).13,14,16,19-22 On the other hand, PEN concentration in commercially used PET-PEN copolymers does not exceed 2% per mole, due to cost considerations.13 Why PET-PEI copolymers (even at small PEI concentrations) exhibit improved barrier properties relative to pure PET is an issue that has not been elucidated so

Dynamics and Barrier Properties of PET and PEI 2979

far. Tonelli23,24 calculated conformational characteristics (values of unperturbed mean-square end-to-end distances and torsional potential energies) and correlated them with the inherent dynamic flexibilities (tendency of phenyl rings to flip) in PEP, PEI, and PET polyisomers. He stated that gas diffusivity in these polyesters is controlled by the different dynamic flexibilities characterizing them, rather than by the corresponding static free volumes. Shanks and Pavel25 reported O2 and CO2 diffusivities in amorphous PET samples and related alkylene and isomeric polyesters, and tried to correlate them with free volume data and the number of CH2 groups present in the alkylene part of the monomeric unit in each polymer. The general observation was that O2 and CO2 diffusivities increase with increasing number of CH2 groups; it was also suggested that gas diffusivities in polyesters with an even number of CH2 groups are higher than in those with an odd number. In an effort to shed light on the mechanisms governing the improved barrier properties of PEI relative to PET, we have followed here a hierarchical modeling approach which provides a quantitative description of gas sorption and transport through these two polyisomers, starting from the atomistic level. Keys to the success of our work are (a) the use of a highly detailed atomistic force field in the simulations, (b) the generation of realistic atomistic configurations, representative of the true polyester structure, on the computer, and (c) an analysis of the structural, dynamic, and barrier properties of the two polyesters through averaging over many, statistically uncorrelated, configurations, all of which have been preequilibrated at given temperature and pressure conditions. This paper is organized as follows: In section 2 we present a concise literature survey of simulation approaches for the prediction of gas solubilities and diffusivities in amorphous (molten or glassy) polymers and of past simulation efforts on PET and related materials. In section 3 we describe the methodology followed here to generate well-relaxed structures of the polyesters and study their properties. Results concerning the static properties, segmental dynamics (local relaxation), and barrier properties (O2 diffusivity) of PET and PEI are presented in section 4. In section 5 we summarize our findings and conclusions and discuss plans for future work. 2. Prediction of Gas Permeability in Amorphous Polymers Permeation is a phenomenon encompassing mechanisms that span a wide spectrum of time and length scales, exceeding by many orders of magnitude the corresponding scales that can be tracked by conventional atomistic molecular dynamics (MD) and Monte Carlo (MC) simulation techniques. Mathematically, defining the permeability of a polymer to a smallmolecule gas amounts to calculating the product of diffusivity and solubility of the gas in the polymer.26,27 Solubilities and diffusivities are governed by the penetrant size and interactions with the polymer as well as by the shape, size, connectivity, and time scales of thermal rearrangement of unoccupied space within the polymer. At high enough temperatures (well above the glass transition temperature of the polymer, Tg) the polymer matrix undergoes fast local thermal motions (segmental relaxation) that cause openings (often described as the accessible volume) among chains capable

2980 Karayiannis et al.

of accommodating the penetrant and their connectivity to undergo rapid redistribution. During these thermal motions of the surrounding chains, penetrant molecules are also carried along by hopping between accessible cavities. As a consequence, penetrant diffusion proceeds through a large succession of small, local moves (jumps) of the fluid molecule caused by rapid density fluctuations. The time scale between penetrant moves is set by the relaxation time of these density fluctuations on the length scale of the penetrant size. At low temperatures (e.g., below Tg), diffusion takes place through a more complicated mechanism.28 Glassy polymers are not in thermodynamic equilibrium. In addition, they are much less mobile than high-temperature melts; their configurations are trapped in regions around local potential energy minima, with transitions from one region to another being prohibited by highenergy barriers. The distribution of accessible cavities remains largely unaltered within the time scale of the diffusive process, undergoing small fluctuations.29 Consequently, a gas molecule inserted into a glassy polymer spends most of its time trapped in clusters of accessible volume which remain disjointed most of the time. Only infrequently does the molecule execute hops from one accessible volume cluster to another through channels opened momentarily by thermal fluctuations in regions of low density or enhanced molecular mobility.28,30 Such elementary diffusive jumps can be described by a characteristic, first-order rate constant.31-33 Diffusivity in glassy polymers is consequently orders of magnitude slower than in high-temperature melts; it is also strongly dependent on the number and connectivity of clusters of accessible volume and on the distribution of rate constants governing penetrant transitions.27,30-34 Diffusion of small molecules in polymer melts can be safely calculated by employing MD and invoking the Einstein relation, in the long-time, hydrodynamic or Fickian limit.30,35 In contrast, in a glassy polymer, diffusion is too slow (by about 2-3 orders of magnitude) to be tracked by brute-force MD. In this case, the distribution of rate constants characterizing the elementary diffusive jumps of a low-molar-mass substance through the polymer matrix can be calculated by applying transition-state theory (TST) on atomistically detailed models, evolving via a succession of infrequent (rare) events.36 For each transition, the “reaction trajectory” leading from a local energy minimum to another through a saddle point in configuration space is tracked, and the transition rate constant is evaluated. For polymeric glasses, TST was first introduced by Arizzi31 and Gusev et al.32,35 According to the Gusev-Suter TST method, a 3-dimensional (3-D) grid of high resolution (spacing smaller than 3 Å) is initially laid on a wellrelaxed polymer structure, covering the entire volume of the simulation cell. Next, a test particle, described as a united atom with dimensions representative of the gas penetrant, is inserted into the polymer matrix at each point on the grid, and the resulting nonbonded intermolecular potential energy, associated with interactions with all polymer atoms, is calculated. These calculations allow separation of the polymer matrix into clusters of accessible volume (low-energy regions, containing local minima of the energy) and regions of excluded volume, characterized by higher interaction energies. Each local minimum of the energy is associated with a sorption “microstate” in the configuration space of the penetrant. Borders between microstates are


defined as high-energy surfaces separating the local energy minima, using the grid construction. Penetrant diffusion is viewed as resulting from infrequent transitions between adjacent microstates. A first-order rate constant, kifj, is assigned to each transition from microstate i to microstate j.33 For low enough temperatures, for which the polymer matrix can be considered as frozen, the rate constant for the transition from microstate i to microstate j is given by31,32,35

kifj )

kBT Qij h Qi

(1)

where h is Plank’s constant, kB is the Boltzmann constant, and Qi and Qij are the partition functions of the molecule in microstate i and on the separating surface between microstates i and j, respectively. Applications of the TST as formulated by Gusev et al.31,35 to rigid (frozen) polymeric structures suggested that small thermal vibrations of the atoms comprising the polymer matrix should also be taken into account. Thus, in a modified formulation, Gusev et al.32,35 account also for the thermal motion of polymer atoms under the simplification/assumption that the matrix atoms execute independent “elastic” motions, each matrix atom being tethered to its average equilibrium position through a harmonic spring. This simplification restricts the study of transport properties to molecules of small size (He, O2, N2, CH4); it is not appropriate for larger molecules (CO2) that may force the polymer atoms in the vicinity of the penetrant to undergo substantial local relaxations to accommodate the guest molecule. In the modified version of the Gusev-Suter TST, every polymer atom i experiences small displacements, ∆ri ) ri - 〈ri〉, from its equilibrium position 〈ri〉. These displacements follow a Gaussian probability density function of the form32,35

(

W(∆r1,∆r2,...,∆rN) ∝ exp -

N

∑ i)1

3∆ri2

2〈∆ri2〉

)

(2)

where 〈∆ri2〉 is the mean square displacement from the equilibrium position, usually assigned a common (average) value, 3∆2, ∆2 being the “smearing factor”, for all atoms. With this form of W(∆r1,∆r2,...,∆rN), given the value of ∆, the jump rate constants can readily be calculated through eq 1.32,35 Greenfield and Theodorou calculated the distribution of jump rate constants and of the lengths traversed during the diffusion of CH4 molecules through atomistically detailed model configurations of atactic polypropylene using a multidimensional transition-state theory,33 which explicitly incorporates the degrees of freedom of the matrix in the transition path determination and rate constant calculation. Transition lengths between microstates were found to be very narrowly distributed, with an average value around 5 Å; in contrast, the corresponding distribution of rate constants was calculated to be very broad, covering up to 18 orders of magnitude.33 After coarse-graining from the microstate to a “macrostate” level, each macrostate comprising several microstates connected through saddle points of energy on the order of kBT or less, a reverse Monte Carlo strategy was adopted to generate large networks of macrostates that exhibit the spatial correlations, connectivity, and transition rate constant



distribution determined through the atomistic calculations. Penetrant motion in these networks was tracked with kinetic Monte Carlo simulation, and the diffusivity was extracted from the slope of the mean square displacement (msd) versus time at long times, when Fickian diffusion is established.31,37,38 Alternatively, diffusivity can be calculated theoretically39 from the rate constant distribution by applying the time-dependent effective medium approximation40,41 (EMA). The Gusev-Suter TST method has been used in recent years by Hofmann and co-workers,27,42-45 Fried and Ren,46 Lo´pez-Gonza´lez et al.,47 and Kucukpinar and Doruker48 to predict the permeability of small gas molecules in different polymeric membranes. In all cases the results obtained have been promising, except for CO2, whose anisotropic shape and large size cause nonnegligible segmental rearrangements in the polymer matrix that cannot be captured by the smearing factor ∆2 of Gusev-Suter.44 Key to the successful implementation of TST for the reliable prediction of barrier properties is the generation of well-equilibrated model configurations of the polymer matrix, exhibiting realistic packing of the constituent atoms and chains. For relatively simple polymers, such as polyethylene (PE), polypropylene (PP), polyisoprene (PI), and polybutadiene (PB) above Tg, this issue has been solved through the design and implementation of a set of very drastic chain-connectivity-altering Monte Carlo moves, such as end-bridging and double-bridging.49-52 For more complicated polymers, however, such as polyesters, polyimides, and poly(amide-imide)s, these algorithms are not yet available; consequently, one should resort to other ways of generating well-relaxed initial structures. As far as PET and its copolyesters are concerned, a careful literature survey shows that, despite the huge body of experimental work addressing a number of issues related to their structural, dynamic, and barrier properties, computer simulation studies of these polymers have been rare. One reason for this is the high glass transition temperatures Tg of pure, amorphous PET and PEI (around 80 and 65 °C, respectively), rendering both polyesters glassy at room temperature.2,7-9,13,20-22,53 Simulating glassy polymers reliably is much more difficult than simulating melts, because of the intrinsically nonequilibrium, history-dependent nature of the former and their very long relaxation times, which can clearly not be tracked by brute-force MD. The chemical complexity of PET and PEI, which for the moment rules out the use of advanced Monte Carlo equilibration methods49-52 even in the melt state, constitutes another difficulty. Simulation studies of structure and relaxation in amorphous PET with atomistic MD include those of Hedenqvist et al.54 and Boyd and Boyd.55 Diffusion and sorption calculations of (a) CH4 molecules in high-temperature PET and PEN melts have been reported by Bharadwaj and Boyd56 and (b) O2 and CO2 in PET and related alkylene and isomeric polyesters by Shanks and Pavel.25 We also mention the work of Cail et al.57 on the orientational behavior of PET chains by means of rotational isomeric state Metropolis Monte Carlo (RMMC) simulation,58 and of Tonelli on the conformational characteristics of PET, PEI, and PEP.23,24

inflexible and polar groups along the backbone, such as PET and PEI, sampling well-equilibrated melt configurations is still a major challenge. This challenge can in principle be overcome by adopting a coarse-graining strategy59 similar to that developed by Tscho¨p et al. for polycarbonates,60 Aoyagi et al.,61 and Mu¨ller-Plathe.62 According to these strategies59 the detailed atomistic structure of the rigid polymer is mapped onto a coarsegrained one involving much fewer degrees of freedom, wherein whole groups (e.g., the benzene rings) are represented as single effective interaction sites. Intramolecular interaction parameters (effective bond stretching, bending, and torsional potentials) for the coarsegrained model are determined from the corresponding distribution functions of unperturbed single-chain conformations sampled with atomistic MC. Intermolecular interactions in the coarse-grained model, usually represented as hard- or soft-sphere repulsions, are extracted from the intermolecular pair distribution functions obtained from atomistic MD of smallmolecular-weight analogues in the liquid state. A coarse-grained equilibration methodology is not available for aromatic polyesters. In view of this, an alternative pragmatic approach that can be followed is that proposed by Hofmann and co-workers27,42-45,63-66 for the simulation of the transport properties of small molecules in amorphous rubbery polysiloxane and stiff glassy polyimide systems. This methodology, based on commercial simulation packages that afford the use of validated, detailed atomistic force fields, is modified here and implemented for PET and PEI. Following the methodology of Hofmann et al.,27,42-45,63-66 the simulation of PET and PEI polyesters has been carried out by employing the Cerius2 (version 4.2) and Materials Studio (version 2.1) software packages from Accelrys Inc.67 The simulation procedure starts off by choosing the proper repeat units for PET and PEI (through the “Polymer Builder” menu of Cerius2) in an explicit atom (EA) description (with all atoms being considered as individual, spherical interacting sites) and growing them bond-by-bond to build a single PET (or PEI) chain. Here, a “parent” chain consisting of 80 repeat units is generated (containing a total of 1762 atoms, with molar mass 15362 g/mol). The initial packing cell, created by employing the “Amorphous Builder” module of Cerius2, has large enough dimensions that the resulting cubic simulation box is characterized by very low density, close to 0.1 g/cm3, to avoid unrealistic packings of atoms along the chain (e.g., polymer atoms spearing through phenyl rings, or two phenyl rings being concatenated).27 To create equilibrated structures, representative of the polyester system under the specified conditions of pressure (P) and temperature (T), an equilibration cycle is next applied, consisting of many, relatively short, NPT and NVT MD stages. According to this procedure, the initially built packing cell is subjected to a static structure optimization using a molecular mechanics (MM) algorithm68,69 until the maximum energy gradient with respect to all microscopic degrees of freedom becomes less than 10 kJ/(mol nm). This is followed by a long succession of 12 MD stages, consisting of NVT and NPT MD simulations (corresponding to compressions, decompressions, annealings, and coolings under various conditions of temperature and pressure), aimed at creating a final structure with realistic density, topology, atomic packing, and low-potential-energy charac-

3. Sampling Uncorrelated, Relaxed Structures of PET and PEI Melts and Glasses As mentioned in the previous section, for polymers with complex chemical constitutions bearing bulky



Table 1. Successive Steps in the Original Hofmann et al. Equilibration Cycle (Ref 27), Used To Generate Relaxed Amorphous, Glassy Polyester Structuresa equilibration step

simulation conditions

duration (ps)

1 2 3 4 5 6 7 8 9 10 11 12

NVT, 600 K NVT, 300 K NPT, 1000 bar, 300 K NVT, 600 K NVT, 300 K NPT, 30000 bar, 300 K NVT, 600 K NVT, 300 K NΡT, 5000 bar, 300 K NVT, 600 K NVT, 300 K NPT, 1 bar, 300 K

50 50 50 50 100 50 50 100 50 50 100 1000

a In the course of this work, as described in section 3, additional steps were introduced to generate better-relaxed configurations.

teristics. The successive stages of the equilibration cycle and the conditions imposed at each simulation step are listed in Table 1. An advantage of the Hofmann et al. method is its speed when taking into account the relaxation times of dense, glassy structures. It suffers, however, from two major drawbacks: (a) In the end, only one equilibrated glassy polymer structure is created (as prescribed by the final NPT MD simulation stage), with a density near the experimental one. Repeating stages 5-12 does not lead to different configurations, since the final glassy structures created are more or less “trapped” always in the vicinity of the same minimum-energy configuration. (b) The short duration of all intermediate MD stages does not guarantee that the configuration finally generated is fully equilibrated at all length scales; in particular, the long-chain (end-to-end vector) characteristics are most probably not well equilibrated. Both disadvantages can be overcome by generating a large number of additional uncorrelated polymer configurations, in the spirit of ref 68. To this end, the entire procedure is repeated from the very beginning, and properties are estimated by averaging over all final configurations. Alternatively, one can carry out a long MD simulation run at very high temperature, where equilibration is much easier to achieve, followed by successive coolings of carefully selected uncorrelated melt configurations down to the glassy state (see below). Although, in both cases, the procedure followed for generating polymer glasses does not correspond to well-defined histories of vitrification, arithmetic averaging can lead to realistic predictions for the properties of interest, on the basis of the assumption that the generated configurations obey a probability distribution comparable to that encountered in the real-life glassy polymer.68 All MD simulation runs prescribed in Table 1 were conducted using the Nose´-Hoover thermostat70,71 (with the value of the parameter Q used to scale the fictitious mass of the Nose´-Hoover temperature set equal to 1) and the Andersen barostat72 (with the value of the cell mass parameter set equal to 20 amu). Isotropic deformations of the simulation box were used in the NPT ensemble, and the integration time step was set equal to dt ) 1 fs. Also, in all simulations, the PCFF force field (polymer-consistent force field)73,74 was chosen to describe all energetic interactions. To validate the computations, some of the intermediate equilibration stages were repeated with the more accurate COMPASS

force field (condensed-phase optimized molecular potentials for atomistic simulation studies).75,76 Predicted density and radius-of-gyration values with both force fields were very similar. COMPASS is widely considered as the most accurate and advanced general-purpose force field available, especially when compared to the simplified united-atom models77-79 often used for chemically simpler chainlike systems.80 Its main drawback is the large CPU time required to calculate all terms in the expression for the potential energy, many of which involve complicated cross terms. For the systems studied here (constructed from single 80-mer long PET or PEI parent chains), completing the equilibration cycle described in Table 1 with PCFF took approximately 25 days of real time, as compared to 5-6 months if the cycle had been completed with COMPASS. On the basis of these estimates and given that no significant differences were detected in the topology and molecular packing of the PET or PEI structures generated with the two force fields, PCFF was the force field of choice. All Cerius2 simulations were performed on a single-processor Origin R10000 SGI workstation at 180 MHz; the corresponding Materials Studio simulations were executed on a singleprocessor Intel Pentium4 at 1.8 GHz. As briefly mentioned above, a drawback of the Hofmann et al.27,42-45,63-66 method is the short (0.050.1 ns) duration of all intermediate MD steps in the equilibration cycle (apart from the concluding NPT step). To come around this problem and enable sampling of as many uncorrelated PET and PEI configurations as possible, the equilibration cycle was significantly extended in our work as follows. (1) The volume of the final configuration resulting from the successive MD simulations described in Table 1 is gradually expanded in the course of NPT MD simulations (through the successive application of negative pressures, the level and duration of which are chosen by trial-and-error runs) so that the system density drops to 0.7 g/cm3. A static structure optimization follows, until the maximum energy gradient becomes less than 5 kJ/(mol nm). (2) The resulting configuration is subjected to a 2 ns long NVT MD simulation at T ) 2000 K, with system configurations being recorded every 0.5 ns over the last 1.5 ns of the run. (3) Each recorded configuration is subjected to successive NPT MD simulations with the temperature set at T ) 600 K and the pressure being gradually adjusted, for (a) the density to match the experimentally measured one at the same temperature and (b) the stress tensor components to attain values corresponding to the level of the applied pressure. The latter requirement was dictated by the observation that, when the set pressure is abruptly changed from a very high level (e.g., 1 GPa) to a significantly lower level (e.g., 0.1 MPa), the stress tensor components in the system do not always attain consistent values; instead, significant levels of residual stresses may remain in many of the offdiagonal terms. To get rid of those, a milder compression procedure has to be followed, consisting of successive applications of smaller step changes down to the desired value (0.1 MPa). In the present work, the optimal compression rate was decided with trial and error: the system was found to respond smoothly to the applied pressure changes only when successive pressure levels differed by less than 1 order of magnitude each time, resulting in five to six intermediate compressions to


bring the initial pressure of 1 GPa down to the atmospheric pressure. Such a simulation strategy was seen to be necessary in all simulations carried out in this work. These NPT runs at 600 K have been the second most CPU-time-consuming step, after the long NVT run at 2 ns. The duration of each intermediate NPT run was exactly 100 ps (which was found to be long enough for the density to reach asymptotic values), except for the last run, for which the simulation lasted twice as long, to ensure that the average value of all off-diagonal stress components was identically zero. Since no experimental data were available for the density of amorphous PEI at T ) 600 K, the target density in that case was assumed to be roughly equal to that of PET. (4) Each configuration obtained is subjected to static energy minimization, followed by an NVT MD simulation at T ) 600 K for 100 ps. (5) Each resulting configuration is subjected to a 3 ns long NVT MD run at T ) 600 K. By averaging over these runs, all properties of interest at this temperature are calculated. For each polyester, the original Hofmann equilibration cycle (HEC) was applied to four different initial configurations. However, only one of these was subjected to the additional equilibration stages introduced in this work, from which four uncorrelated samples were obtained at T ) 2000 K. (Had all four HEC structures been subjected to the extended cycle, 16 ()4 × 4) uncorrelated samples would have been produced at T ) 2000 K). The choice to produce only four was dictated by CPU time limitations, as 1 ns of MD was seen to take about 25 days of real time to complete. Each one of the four uncorrelated samples was subjected next to the NVT MD simulation at T ) 600 K, from which three more statistically different structures were obtained. This means that, in the end, the total number of available samples was 12. All these could have been used in the subsequent TST calculations; here, however, only 10 of these samples were utilized, as they were found enough to reliably calculate configurational averages for the properties of interest. Note that one has the freedom to generate additional uncorrelated structures by extending the NVT MD at T ) 2000 K to times longer that 3 ns; the same holds true for the NVT MD run at T ) 600 K. The extended cycle allowed us to accumulate reliable data on the structural, conformational, and dynamic properties of PET and PEI at T ) 600 K, and compare them in detail. The same cycle was employed for studying differences in the dynamic and barrier properties of the two polyesters at two lower temperatures (T ) 450 K and T ) 300 K), where again model densities were forced to match experimental values. NVT MD trajectories 3 ns long at 450 and 300 K were obtained successively from the final T ) 600 K configurations by repeating the last three steps of the scheme described above. The low-density (0.7 g/cm3) and high-temperature (2000 K) levels prescribed in the extended equilibration cycle guarantee system equilibration at all length scales within the simulation time of 2 ns (which corresponds to about 50 days of real time), as the system undergoes significant global and local rearrangements that lead to decorrelation from the configuration of the initial structure. Exploratory NPT MD test simulation runs at a lower temperature (T ) 1000 K), and at P ) 1 atm, proved insufficient to provide complete system de-


Figure 2. Time evolution of the mean square chain end-toend distance 〈R2〉 of a PEI model system during the 2 ns long NVT MD simulation (T ) 2000 K, F ) 0.70 g/cm3). Circles denote frames from the trajectory taken for cooling to lower temperatures.

Figure 3. Time evolution of the calculated pressure of PET and PEI systems during the 3 ns long NVT MD simulation at T ) 600 K.

correlation within the 2 ns period: on one hand, even at 1000 K, the system density is still quite high, so that molecular mobility remains sluggish; on the other hand, despite the use of the simpler (as compared to COMPASS) PCFF force field, potential energy calculations are still too expensive to allow monitoring of system evolution for times longer than 2 ns within reasonable CPU time. It is only at very high temperature (close to 2000 K) and low density (0.7 g/cm3) that polymer atoms in the generated PET and PEI polyester structures exhibit reasonably high mobilities, traceable by bruteforce MD in reasonable CPU time. That PET and PEI melts are equilibrated within the 2 ns long MD simulation at T ) 2000 K and F ) 0.7 g/cm3 is supported by Figure 2, which displays the instantaneous values of the mean square end-to-end distance, 〈R2〉, of the single, 80mer long PEI parent chain. 〈R2〉 is seen to exhibit significant fluctuations ((5000 Å2) in the course of the MD run, indicative of ergodic and robust sampling by the MD method. The solid symbols in the figure denote the frames along the trajectory from which structures were chosen to be cooled to the temperatures of interest (T ) 600 K initially, and then down to 450 and 300 K) for analysis of their structural, dynamic, and barrier properties. As stated above, the densities of the initial PET and PEI structures that were subjected to the 3 ns long NVT MD runs at T ) 600 and 450 K were set



Table 2. Simulation Results for (a) the Density of Purely Amorphous PET and PEI Glass and Melt Systems and (b) O2 Diffusivity through Thema polyester

T (Κ)

F(simulation) (g/cm3)

F(exptl) (g/cm3)

D(TST) (10-8 cm2/s)

D(exptl) (10-8 cm2/s)

PET PET PEI PEI

300 600 300 600

1.340 ( 0.02 1.135 (imposed) 1.332 ( 0.04 1.135 (imposed)

1.336,b 1.333d 1.135 1.346,b 1.356d

5.30 ( 2.5 120 ( 50 7.1 ( 4.0 165 ( 80

0.95c (0.65)b 0.27b

a Also shown are the corresponding experimentally measured data (where available). The results have been obtained by averaging over four uncorrelated structures produced through the equilibration procedure described in Table 1. b Reference 8. c Reference 87. d I. Dairanieh (BP). Personal communication.

(through successive steps of NPT MD runs) equal to experimentally known values. Typical profiles of the corresponding instantaneous pressure in the course of a 3 ns long run are shown in Figure 3; they have been obtained from two different simulations with a PET and a PEI sample, respectively. Despite its disadvantages, the extended equilibration cycle proposed above, relying solely on the application of brute-force MD and energy minimization, is capable of generating uncorrelated PET and PEI configurations which remain stable with simulation time and exhibit reasonable properties. To further reduce the statistical error, all results presented in the next section concerning the structural and dynamic properties of PET and PEI have been obtained by averaging over the four independent initial configurations chosen from the NVT MD run at T ) 2000 K (those shown by the solid symbols in Figure 2), which then were cooled to T ) 600 K (and later to T ) 450 K and T ) 300 K) and subjected to the 3 ns long NVT MD run. The same holds also for the barrier properties of the two polyesters. At the higher temperatures studied (T ) 450 and 600 K), averaging over these four different samples was seen to be enough to capture the corresponding properties (and how they differ between the two polymers) with the desired accuracy. At room temperature (T ) 300 K), however, a larger number of samples (about 10) was needed to be averaged over for the accurate prediction of O2 diffusivity in the polyester matrices. We summarize this section by pointing out that the extended equilibration strategy implemented in this work, involving in total 19 different steps (and quite expensive in CPU time), has allowed us to (a) sample well-relaxed polyester melt structures at temperatures as low as T ) 600 and 450 K and analyze their segmental dynamics in detail and (b) generate, by cooling these structures to lower temperature (T ) 300 K), a large number of statistically uncorrelated glassy configurations at the experimental density and calculate their O2 barrier properties. 4. Results 4.1. Static Properties of Amorphous PET and PEI. 4.1.1. Density. As noted in section 3, the PET and PEI model systems employed in this work consisted of a single parent macromolecular chain 80 repeat units long. Application of the Hofmann et al. equilibration cycle (presented in Table 1), starting off at four different initial configurations for each polyester (PET and PEI), led to four different final glassy structures (at T ) 300 K and P ) 1 atm) with average densities as shown in Table 2. Also reported in the same table are the corresponding experimental densities of the purely amorphous PET and PEI at the same conditions of temperature and pressure. The data of Table 2 show that simulation and experimental density values for

both polyesters are in good agreement, the maximum relative error being less than 2%. This should be considered as an excellent comparison, given that density predictions even for simpler polymers, such as PE, PP, PI, and PB, are rarely better than within 2% of the experimental data. Our predicted values, however, do not follow the trend observed in the experimental measurements, namely, that amorphous PEI is slightly denser than amorphous PET by about 0.5-1%. This difference, being within the statistical error of the simulations, is too small to be captured by our modeling approach; it will, however, prove critical to the correct prediction of the transport properties of the two systems (see below). Typical atomistic snapshots of the amorphous PET and PEI cells (at the end of the Hofmann et al. equilibration cycle) subject to periodic boundary conditions and of the corresponding parent chains fully unwrapped in space are shown in Figures 4 and 5, respectively. Such glassy structures constituted the starting points for the subsequent extended equilibration cycle simulations described in section 3, whose outcome was the generation of fully relaxed PET and PEI melt samples at two different temperatures (T ) 600 and 450 K), as well as of glassy structures at T ) 300 K. In all these samples, densities exactly match the corresponding experimental values (equal to, e.g., 1.135 g/cm3 for PET at T ) 600 K). Results concerning the conformational, structural, dynamic, and barrier properties of PET and PEI at these temperatures are reported and discussed in detail in the following paragraphs. 4.1.2. Conformational Properties. Chain dimensions in PET and PEI are quantified through the calculation of the equilibrium mean square chain endto-end distance 〈R2〉. This is used to define the characteristic ratio C∞ in the limit of infinitely long chains according to

C∞ )

〈R2〉

∑i

(3)

Nili2

where Ni and li are the number and length of bonds of kind i, respectively, along the main chain backbone. Experimentally measured C∞ values range from 3.5 to 5.8 for PET54,81 and from 4.3 to 7.8 for PEI.81 Tonelli suggested that chain dimensions in PEI should be smaller than in PET (C∞ ) 4.13 as compared to 4.67).24 Figure 6 displays the time evolution of the mean square chain end-to-end distance 〈R2〉 for the four different samples of PET and PEI subjected to the NVT MD simulations of the extended equilibration cycle at T ) 600 K. The curves present significant differences in both polymers and are seen to remain practically unaltered



Figure 5. Same as Figure 4 but with the PET (top) and PEI (bottom) parent chains fully unwrapped in space.

Figure 4. Typical snapshots of amorphous cells of singleparent-chain, 80-mer long PET (top) and PEI (bottom) systems subject to periodic boundary conditions, obtained by the Hofmann et al. equilibration cycle. Carbon, oxygen, and hydrogen atoms are displayed in blue, red, and black, respectively (T ) 300 K, P ) 1 atm).

in time, indicative of the inability of MD to equilibrate the long-range characteristics of the two polyesters within the available simulation time (3 ns, corresponding to 75 days of real time) at this temperature. The two polyester structures are so dense that 3 ns is too short a time for the PET and PEI chains to diffuse and assume dimensions characteristic of a fully equilibrated bulk system. Only through extensive averaging over many, independently generated, configurations can one attempt to obtain a better estimate of 〈R2〉. In the present case, such averaging leads to characteristic ratios equal to 4.0 ( 1.3 and 4.3 ( 1.3 for PET and PEI, respectively. Here, we point the readers’ attention to the small size of the systems simulated (each constructed from just one parent chain); this increases the statistical uncertainty of the calculated C∞ values, but also leaves open the possibility of system size effects, given that the box dimension (27.5 Å) was in all cases comparable to the mean chain radius of gyration (about 25 Å). An estimate of the chain stiffness in the two polyesters can also be provided by sampling continuous unperturbed chains through, for example, the RMMC algorithm.57,58

Figure 6. Time evolution of the mean square chain end-toend distance 〈R2〉 of the four different initial configurations of PET and PEI, during the 3 ns long NVT MD simulation at T ) 600 Κ.

C∞ values calculated with this method were again seen to be broadly distributed, but their mean was not significantly different from the values reported above from the bulk simulations. More insight into the conformational properties of PET and PEI systems can be gained by analyzing the distributions of the different kinds of torsion angles along the chain. Torsion angles were labeled here according to the Hedenqvist et al.54 convention, as listed in Table 3. For PET, the labeling is also schematically



Figure 7. Notation used to represent different atoms in the PET repeat unit (see Table 3).

Figure 10. Same as Figure 8, but for the distribution of the ΟD-CD-O-C torsion angle.

Figure 8. Distribution of the CP-CP-CP-CD torsion angle in PET and PEI melts. The results have been obtained by averaging over all four samples subjected to NVT MD simulation at T ) 450 and 600 K. The coplanar conformation with CP-CP and CP-CD bonds antiparallel corresponds to 0°.

Figure 11. Same as Figure 8, but for the distribution of the CD-O-C-C torsion angle.

Figure 9. Same as Figure 8, but for the distribution of the CP-CP-CD-OD torsion angle. The trans conformation corresponds to 0°. Table 3. Symbols Used To Represent Atoms of Different Kinds in PET and PEI Repeat Units (According to Ref 54) atom identity

symbol

atom identity

symbol

phenyl carbon carbonyl carbon (sp2) glycol carbon (sp3)

CP CD C

hydrogen carbonyl oxygen glycol oxygen

H ΟD O

explained in Figure 7; it is entirely analogous for PEI. Torsion angle distributions (obtained again as configurational averages over available different material samples) for the two polyesters at T ) 600 and 450 K are shown in Figures 8-12. In all cases, the convention is adopted that the trans conformation corresponds to the zero angle. Figure 8 presents the distribution of the CP-CP-CP-CD torsion angle, Figure 9 that of the CP-

Figure 12. Same as Figure 8, but for the distribution of the O-C-C-Ο torsion angle.

CP-CD-OD angle, Figure 10 that of the OD-CD-O-C angle, Figure 11 that of the CD-O-C-C angle, and Figure 12 that of the O-C-C-O angle. According to the torsion angle distribution functions shown in Figures 8-12, the static conformations of PET and PEI are practically identical: Although individual torsional angle distributions, as calculated from a single MD simulation (see below), are not smooth enough to analyze and compare, the curves obtained after averaging over four or more different trajectories are identical between the two polyisomers. This holds for all torsion



Figure 13. Distribution of the Ο-C-C-Ο torsion angle in PEI. The results have been obtained by averaging separately over each of the four samples subjected to NVT MD simulation at T ) 600 K. The trans conformation corresponds to 0°.

Figure 14. Total radial distribution function g(r) as obtained from our NVT MD simulation runs with glassy PET and PEI systems (T ) 300 K, P ) 1 atm). All atom types were used in accumulating g(r).

angles, both in the vicinity of the phenyl ring and ester group and in the alkyl part of the polyesters. Some asymmetries showing up in the pattern of the distribution function for the glycol group (O-C-C-O) at T ) 450 K for PET (see Figure 12) are attributed to incomplete averaging, similarly to what has been reported in the literature for glassy (rigid) polymeric systems,54,82 especially at low temperatures. For the PET and PEI systems, these asymmetries are caused by (a) the short duration of the MD simulations (total run time for each sample equal to 3 ns), (b) the rather small population of distinct torsion angles in each simulation box, and (c) the large residence times in the trans, gauche+, and gauche- wells of the potential energy function (this is related to the low molecular mobility in these dense systems, resulting in a limited number of conformational transitions). The application of the very complex PCFF force field intensifies the phenomenon. To further illustrate that configurational averaging alleviates many of the problems related with poor (and sometimes nonergodic) sampling, Figure 13 shows the individual curves for the distribution of the O-C-C-O torsion angle in PEI as obtained from the four independent NVT MD runs at T ) 600 K. It is seen that the population of gauche- angles is significantly overestimated in the samples labeled as 1 and 2, and underestimated in the samples labeled as 3 and 4. The “mean” distribution, however (shown by the thin broken line in Figure 12), that arises after averaging over all samples is more symmetric and smoother. As far as the distribution of individual torsion angles is concerned, inspection of Figures 8-12 shows the following: (a) The CP-CP-CP-CD torsion angle for both polyisomers is peaked around 0° (with a fluctuation window of (20°), indicative of the coplanarity of the carbonyl carbon and the phenyl ring; in fact, the additional curves of Figure 9 suggest that the entire carbonyl group is coplanar with the phenyl ring. (b) Glycol C-O bonds assume trans and gauche conformations with almost equal probabilities (see the distribution of CD-O-C-C torsion angles in Figure 11); however, the positions of the gauche+ and gauche- maxima are seen to deviate somewhat from the conventional values of +120° and -120°. This is attributed to excluded-volume interactions between the neighboring atoms. (c) The distribution of the torsion angle associ-

ated with the C-C bond in the glycol group (Figure 12) exhibits peaks of almost the same height at the trans and gauche populations. (d) Figure 10 also demonstrates that the OD-CD-O-C torsion angle prefers the cis conformation (-180° and +180°). Although this seems rather strange at first (the cis conformation brings bonded atoms closer than the trans one, causes stronger repulsive interactions, and should be avoided), it can be explained by noticing that the orientation of the CD-O bond is greatly affected by the orientation of bonds in the CP-CD-O-C sequence, which assumes largely the trans conformation. Thus, it is the sum of the torsional energies associated with both torsion angles (OD-CD-O-C and CP-CD-O-C) that should be minimized, and this causes the CD-O bond in the ODCD-O-C sequence to be mainly in the cis conformation. 4.1.3. Structural Properties. Information about the structural features of PET and PEI at the atomic level is also provided by the (total) radial distribution function (RDF), g(r), whose Fourier transform determines the X-ray diffraction pattern of the sample. Calculated RDFs for the glassy PET and PEI systems at T ) 300 K and P ) 1 atm are shown in Figure 14. All atoms, regardless of type, were included in the calculation for this figure. Peaks at small distances reflect mainly intramolecular arrangements: The first peak at around 1.1 Å corresponds to the bond distance between carbon and hydrogen atoms (Cp-H, C-H); it is also broad enough to cover the bond distance of ∼1.2 Å of oxygen and carbon atoms in the carbonyl group (CD-OD). The second intramolecular peak comes from the CD-O (∼1.37 Å), CP-CP (∼1.41 Å), C-O (∼1.43 Å), and C-C (∼1.53 Å) bond lengths. The third and subsequent intramolecular peaks are attributed to distances between atoms two bonds apart, such as between (a) hydrogen atoms in Η-C-H sequences (∼1.77 Å), (b) hydrogen and carbon atoms in Η-CP-CP sequences (∼2.16 Å), and (c) carbon atoms in CP-CP-CP sequences (∼2.44 Å). In general, all peaks up to 4 Å are of intramolecular origin; subsequent peaks reflect mainly relative positions between intermolecular neighbors and phenyl-phenyl intrachain correlations,83 very similar in the two polymers, except from the distances at ∼5-6 Å. Differences there are caused by differences in the sequence (para versus meta) of the repeat units in the two polyisomers.


Figure 15. PET and PEI static structure factors S(k) as calculated from our NVT MD simulation runs at T ) 600 K (melt state) and 300 K (glassy state).

By Fourier transforming the partial radial distribution functions g(r), each weighted with the product of scattering factors for the corresponding pair of atoms, one can construct the static structure factor S(k) of PET and PEI. The patterns obtained from our simulation data for the two polymers at two different temperatures, T ) 600 K (melt state) and T ) 300 K (glassy state), are shown in Figure 15. Simulated patterns are seen to follow quite well experimentally measured X-ray diffraction curves:84 In particular, the three peaks with the highest intensities in the experimental intensity versus k plots (see, e.g., p 125 in ref 84) at 1.5, 3, and 5.6 Å-1 are well reproduced in Figure 15. The first peak reflects interpacking and intrapacking (phenyl-phenyl correlations),83 while the next two are of intramolecular origin. The agreement between simulated and experimentally measured X-ray diffraction patterns holds also for the “diffuse” part of the pattern observed at higher wavevectors, arising from the shorter characteristic intramolecular lengths. 4.2. Segmental Dynamics of PET and PEI Melts. The study of the torsional angle distributions, the total radial distribution function, and the static structure factor presented in section 4.1 proves that PET and PEI exhibit very similar, almost identical, static (structural, conformational, and packing) properties. Consequently, the static distribution of accessible volume in the two polymers should not be expected to be significantly different, which cannot explain why PET and PEI have so different barrier properties. To analyze this, we turn our attention in this section to the study of dynamical properties in the two polymers, especially those related to their local or segmental dynamics. All results reported have again been calculated as configurational averages over the four independent, statistically uncorrelated, configurations subjected to 3 ns long NVT MD simulations at T ) 600 and 450 K. Unfortunately, due to the problem of long relaxation times, no results are presented concerning the long-length-scale (chain) dynamics of the two polyesters, since it is governed by time scales significantly longer that the 3 ns tracked in the present MD studies. Local relaxation in the two PET and PEI systems can be probed by calculating the decay of the time autocorrelation functions for the different torsional angles along the chain; these functions express the rate with which the corresponding angles lose the memory of their


Figure 16. Time decay of the orientational autocorrelation function Μu(t) of the unit vector normal to the phenyl ring in PET and PEI, as calculated from 3 ns NVT MD simulation runs at T ) 450 and 600 K.

Figure 17. Same as Figure 16, but for the unit vector normal to the ester group.

initial conformation (i.e., how fast they decorrelate in time). According to Takeuchi and Roe,85,86 the autocorrelation function Pφ(t) of torsion angle φ at time t is defined through

Pφ(t) )

〈cos φ(t) cos φ(0)〉 - 〈cos φ(0)〉2 〈cos2 φ(0)〉 - 〈cos φ(0)〉2

(4)

where angular brackets denote averaging over all different torsion angles of the same kind along the chain, over all time origins throughout the simulation trajectory, over all four different MD trajectories. Alternatively, local mobility can be expressed through the orientational autocorrelation functions of the unit vectors normal to the plane of coplanar atoms or groups, such as the phenyl rings; similarly to the torsional time autocorrelation functions, these functions also provide a measure of how quickly the corresponding vectors “forget” their initial orientation and decorrelate. The orientational autocorrelation function Mu(t) of any normal unit vector u at time t is given by30

Mu(t) ) 〈u(t)‚u(0)〉

(5)

Figures 16 and 17 present the time evolution of the orientational autocorrelation functions Mu(t) of the unit vectors normal to the phenyl ring and the ester group, respectively, in PET and PEI at T ) 600 and 450 K.



Figure 18. Decay of the time autocorrelation function, Pφ(t), of the CP-CP-CD-OD torsion angle in PET and PEI, as obtained from 3 ns NVT MD simulations at T ) 450 and 600 Κ.

Figure 19. Same as Figure 18, but for the CD-O-C-C torsion angle.

The figures show that in PET the phenyl ring presents enhanced mobility, similar to that of the adjacent Od CsO ester group; in contrast, in PEI, the phenyl ring cannot follow the motion of the neighboring ester group, and consequently, its mobility is significantly reduced. This important difference in the motion of the phenyl ring in the two polyisomers is attributed to the way in which the phenyl ring is “hooked” to the two ester groups on its left and right: In PET (see Figure 1, left) the symmetric way in which the two ester groups are attached to the phenyl ring allows the latter to rotate freely (and to flip) around its axis of symmetry without displacing the former. In contrast, in PEI (see Figure 1, right) the connectivity of the phenyl ring and its two side ester groups is such that no symmetry axis exists: For the phenyl ring to flip, the adjacent carbonyl atoms must undergo significant displacements, which is not favored energetically. Although more pronounced at the higher temperature, the enhanced mobility exhibited by the phenyl ring in PET relative to PEI is also observable at T ) 450 K. The consequences of such a difference for the segmental dynamics of the two polyesters are very important. As suggested by Tonelli,24 it can significantly affect the rate at which the free volume is rearranged in the neighborhood of the phenyl ring. Tonelli uses the term “dynamic flexibility” to describe this: The facile ring motion in PET results in fast redistributions of the free volume locally. This reduces the free energy barrier that must be overcome by the penetrant gas molecules; thus, PET is characterized by a higher gas diffusivity than PEI. Here we should also mention that, in rigid systems such as PET and PEI, transitions between different torsion angle conformations appear to be highly synergetic: isolated conformational jumps can only be completed by overcoming high-energy barriers caused by excluded-volume interactions (overlaps with nearestneighbor atoms not participating in the transition) associated with dense atom packing. Such a synergy becomes more and more necessary at lower temperatures, where thermal fluctuations fade away. On the basis of these considerations, the reduced mobility exhibited by the phenyl ring in PEI, especially at T ) 600 K, should induce reduced mobilities also in the adjacent neighboring groups. This is indeed verified in Figures 18-20, where the CD-O-C-C sequence of atoms is seen to be more mobile than the CP-CP-CD-

Figure 20. Same as Figure 18, but for the O-C-C-O torsion angle.

OD and Ο-C-C-O sequences, the latter exhibiting less frequent transitions between different conformations. Even for the glycol group (Ο-C-C-O), which might have been expected to exhibit enhanced mobility similar to that of an alkyl group independently of the meta or para arrangement, the mobility recorded is smaller in PEI than in PET, which should again be attributed to a reduced dynamic flexibility under the influence of the neighboring rigid phenyl rings. Local mobility differences in the vicinity of the phenyl ring in PET and PEI can be analyzed in detail by looking at “time series” of randomly selected CP-CP-CD-OD torsion angles along the PET and PEI chains. For T ) 600 K, such representative trajectories are recorded in Figure 21. The figure shows that the torsion angle CPCP-CD-OD in both polymers undergoes not only large torsional librations but also full 180° “flips” between its preferred states. However, the number of such flips is significantly higher in PET than in PEI. For example, the torsion angle denoted by the solid curve in Figure 21b (belonging to PEI) undergoes only three changes by roughly 180° at times t = 0.37, 1.1, and 2.25 ns within the 3 ns period examined. The first of these changes is reversed practically immediately through a fast recrossing event. On the contrary, within the same time period (3 ns), the corresponding PET angle whose trajectory is depicted by the solid line in Figure 21a undergoes more than 10 actual transitions. The higher mobility in PET can also be evidenced by the larger amplitude



Figure 21. Time series of two, randomly chosen CP-CP-CDOD torsion angles along the chain, in (a, top) PET and (b, bottom) PEI, in the course of the 3 ns long NVT MD run at T ) 600 K. The trans conformation corresponds to 0°. Periodic boundary conditions apply across the top and bottom of the figure.

Figure 22. Same as Figure 21 but for the dot product u(t)‚ u(t-∆t), where u is the unit vector normal to the plane of the phenyl ring at time t and ∆t ) 10 ps, for the same two monomer units examined in Figure 21 in (a, top) PET and (b, bottom) PEI.

of librations executed by angles trapped in a trans conformation over the time scale of observation in the two polymers (dotted curves in parts a and b of Figure 21 for PET and PEI, respectively). For example, the PET angle traced with the dotted line in Figure 21a undergoes librations roughly between +50° and -50° as compared to (20° librations executed by the PEI angle traced with the dotted line in Figure 21b. That the CP-CP-CD-OD torsion angle undergoes 180° flips does not necessarily imply that these are accompanied by corresponding full 180° rotations of the phenyl rings in space. To analyze this, in Figure 22 we present trajectories of the quantity u(t)‚u(t-∆t), where u is the unit vector normal to the plane of the phenyl ring at time t and ∆t ) 10 ps, for the same PET and PEI repeat units whose CP-CP-CD-OD torsion angles were recorded and discussed above (Figure 21). A drastic reorientation of the phenyl ring within ∆t will appear as a sharp downward peak in the plots of Figure 22; a flip of the phenyl by 180° within ∆t will result in a peak reaching down to -1. The results shown in Figure 22 indicate that only in PET do entire rotations of the phenyl group occur. PEI displays no such phenyl ring rotations in space (Figure 22b). Consequently, local segmental dynamics is richer in PET, including both enhanced torsional librations and flips accompanied by full rotations of the phenyl rings, whereas in PEI local segmental dynamics is restricted, involving smaller torsional librations and a limited number of CP-CP-

CD-OD torsion angle changes by ca. 180°, which tend to be reversed within a short time. Similar conclusions were drawn by looking at time series of representative CP-CP-CD-OD torsion angles in PET and PEI at T ) 450 K. At this lower temperature, however, only torsional librations were recorded for both polymers; no 180° flips were ever monitored for the torsion angles randomly selected. This can be seen in Figure 23, which presents the trajectories of two CPCP-CD-OD torsion angles, one belonging to PET (solid line) and another belonging to PEI (dotted line). The absence of flips in the course of the MD run is immediately recognized in this figure. Given the synergy between CP-CP-CD-OD torsion angle flips and adjacent phenyl ring rotations, absent in the course of the 450 K run are also phenyl ring rotations by 180° (results not shown here) for both PET and PEI; this explains why the respective 〈u(t)‚u(0)〉 and Pφ(t) time autocorrelation functions of the two polymers at T ) 450 K, presented in Figures 16 and 18, come close to each other. The reader can further observe from Figure 23 that, at T ) 450 K, the amplitude of the CP-CP-CD-OD torsion angle librations in PET is not very much higher than in PEI. This observation is also consistent with the information obtained from the decay of the 〈u(t)‚u(0)〉 time autocorrelation functions of the unit vectors normal to the ester group in the two polymers, reported in Figure 17.


Figure 23. Time series of a single, randomly chosen CP-CPCD-OD torsion angle along the chain in PET (solid line) and PEI (dotted line), in the course of the 3 ns long NVT MD run at T ) 450 K. The trans conformation corresponds to 0°.

4.3. Oxygen Diffusivity in PET and PEI Glasses. On the basis of reported experimental data, oxygen (O2) diffusivity in amorphous PEI is 2-2.5 times smaller than in amorphous PET; to a smaller degree, this holds also for O2 solubility.8 As explained in section 2, O2 diffusivity in PET and PEI is strongly dependent on the accessible volume (the part of the total free volume that is large enough to accommodate the penetrant molecule), its distribution and connectivity (between accessible cavities) in the polymer matrix, and the frequency with which fast transition paths connecting accessible clusters are formed. Channel formation depends critically on the mobility of the polymer matrix, i.e., on the amplitude of vibrations and torsional librations and on the frequency of transitions between different conformational states. Consequently, O2 diffusivity in PET and PEI should strongly depend on the static (density, accessible volume and its connectivity) and dynamic (modes of motion, distribution of transition rate constants) properties of the system polyester/penetrant. In the melt state, where both polyesters are expected to be of similar density (e.g., F ) 1.135 g cm-3 at T ) 600 K), any differences in O2 diffusivity between the two polyisomers should arise mainly from differences in the segmental mobility. On the basis of the results presented in the previous section (at T ) 600 K), PEI exhibits a significantly lower mobility than PET; accordingly, O2 should diffuse more slowly in PEI than in PET. Since the initial version of Gusev et al.32 TST assumes the polymer matrix to be frozen, accommodating differences in PET and PEI segmental dynamics in the calculation of O2 diffusivity requires that we resort to the improved TST version of Gusev-Suter. In the improved TST version,32,35 thermal fluctuations of the atoms in the polymer matrix are included through the use of the smearing factor ∆2 (see, e.g., eq 2). A strict calculation of this factor could in principle correctly accommodate the significant differences recorded here in the segmental dynamics between the two polyesters. However, the way this factor is proposed to be calculated by Gusev and Suter (by running short-time MD simulations for each polymer and calculating the msd of all atoms in the polymer for a time equal to the mean residence time in a microstate) poses questions on our ability to capture the faster and more extensive phenyl


motions in PET relative to PEI, since these motions may not be accurately reflected in the corresponding msd’s. Even if we attribute different smearing factors to different kinds of atoms in the polymer, we may not fully account for the differences in the phenyl ring motion between PET and PEI, given that smearing factors are based on average msd values. Despite these possible limitations, and lacking a more robust methodology, such as multidimensional TST,29,33 we have adopted the approximate self-consistent scheme method of Gusev and Suter in this work. To calculate oxygen diffusivities in the two polyisomers, the Gsnet and Gsdiff macros of Accelrys, Inc. were used. Both programs are based on the original code developed by Gusev and Suter. The Gsnet macro calculates the free energy surface of the system; it identifies all sorption sites for the given penetrant (O2) and their connectivity, and it determines the intersite rate constants governing the diffusive jumps. The Gsdiff macro performs a kinetic Monte Carlo simulation of a large number of phantom walkers (representing the penetrant) “hopping” between neighboring sorption sites with transition probabilities based on the rate constants calculated before by Gsnet. By running short MD simulations (on the order of 100 ps) with the two polymers, we first calculated the mean square displacement averaged over all atoms in the polymer matrix, for both PET and PEI. Then, the GSnet macro was used to calculate the smearing factors, using as input information the van der Waals radii of the penetrant molecule; the latter was described as a united atom whose short-range interactions with the matrix are captured through a potential of the form

[ ( ) ( )]

UvdW ) 2

σ rij

9

-3

σ rij

6

(6)

with σ and being the PCFF force field parameters for nonbonded interactions between the united-atom oxygen and the polymer atoms. For O2, values of σ ) 3.460 Å and ) 0.2344 kcal mol-1 were ascribed to the unitedatom penetrant molecule. The cutoff distance for the nonbonded interactions was equal to 10 Å. (Test runs with larger values revealed no dependence of diffusivity on cutoff radius.) The smearing factor was calculated through a self-consistent (SC) scheme, involving information about the msd of the polymer atoms from their equilibrium position as obtained from the short MD calculations. The two factors of the self-consistent scheme adopted in the original TST approach by GusevSuter and implemented in the Gsnet and Gsdiff macros are the smearing factor and the most probable residence time (of the penetrant in the sorption sites). Initially, an arbitrary value is set for the smearing factor, and the distribution of the residence times is computed. The most probable residence time (τ) is identified, and the new value of the smearing factor is deduced by the msd data (obtained from short MD runs) for time equal to τ. Convergence occurs when the relative error between two successive values of the smearing factor is less than a critical tolerance. The initial estimate for the smearing factor was 0.404 Å for both polymer matrixes. Convergence was assumed whenever calculated smearing factors between two successive SC cycles differed by less than 2.5%; this required a minimum of three and a maximum of five cycles. With this scheme, the following estimates for the mean smearing factors, applied to the coordinates of all atoms, were obtained: ∆ ) 0.432 Å


Macromolecules, Vol. 37, No. 8, 2004 Table 4. Simulation Results for O2 Diffusivity in PET and PEI Glasses (T ) 300 K)a F (imposed in the simulations) polyester (g/cm3) PET PEI

1.335 1.356

F (exptl) (g/cm3) 1.336b (1.333d) 1.356d (1.346b)

D D (exptl) (10-9 cm2/s) (10-9 cm2/s) 9.6 ( 6.2 5.4 ( 3.7

9.5c (6.5)b 2.7b

a The densities of the corresponding structures involved in the TST calculations have been imposed to match experimentally known data (reported in Table 2). b Reference 8. c Reference 87. d I. Dairanieh (BP). Personal communication.

Figure 24. Results for the mean square displacement (〈r2〉) of oxygen molecules as a function of time in well-relaxed PET and PEI structures at T ) 300 K, as calculated from the selfconsistent theory approach of Gusev and Suter. Curves, displayed in logarithmic coordinates, are averages over 2000 uncorrelated kinetic Monte Carlo simulation paths in a single structure. The dimensions of the simulation boxes are 26.74 and 26.60 Å for PET and PEI, respectively.

for PET and 0.426 Å for PEI. These values were used for the identification of all possible oxygen sorption sites (microstates) in PET and PEI and for the calculation of the rate constants governing penetrant transitions between all pairs of microstates. Oxygen diffusivity was predicted by running kinetic Monte Carlo simulations to monitor the mean square displacement of about 2000 penetrant walkers on the network of microstates over a time period of 10-4 s. For this purpose the Gsdiff macro was invoked. Figure 24 presents the msd of oxygen molecules in well-relaxed PET and PEI structures as computed by TST. The figure shows that, in both polymers, oxygen motion on time scales on the order of hundreds of nanoseconds is not diffusive (Fickian); instead, a region of “anomalous diffusion” is observed. The crossover from “anomalous” to normal diffusion takes place when molecules have been displaced by a length roughly equal to the size of the simulation box. This observation is in agreement with the findings32 of the initial TST approaches for He and Ar dynamics in PC and He and O2 in PIB, and the work of Greenfield and Theodorou,37 which indicated that the crossover to normal diffusion is a simulation artifact imposed by the periodic boundary conditions. The true duration of the anomalous regime can be captured only through repetition of the TST calculations (and the preceding equilibration procedure) in progressively larger boxes, or through the reconstruction of a large network of sorption sites conforming to the connectivity and rate constant distributions extracted from detailed atomistic simulations.37 However, as shown in the work of Karayiannis et al.,39 the long-time diffusivity in a system of given connectivity is a function only of the distribution of rate constants governing the jumps between sorption sites. Short-range spatial correlations may extend the duration of anomalous diffusion, but their effect on the longtime diffusivity is imperceptible. As long as the distributions of the geometric and topological features of the network of sorption sites, and also the intersite rate constants, are correctly realized and calculated by the atomistic simulations and the TST approach, respectively, the estimated diffusivity is not distorted by the

accelerated crossover into the Fickian regime brought about by the periodic boundary conditions. The sensitivity of penetrant diffusivity on the value of the smearing factor used in the calculations can be estimated by inspecting “instantaneous” estimates in the course of the self-consistency (SC) scheme. For example, employing the typical Gusev-Suter TST approach on a randomly picked, well-relaxed (through the extended equilibration cycle) PET sample yielded the following values for the smearing factor (through Gsnet) and the corresponding diffusivity coefficient (through Gsdiff): (1) initial guess, ∆0 ) 0.404 Å, D0 ) 6.005 × 10-9 cm2/s, (2) first iteration, ∆1 ) 0.444 Å (relative error 9.9%), D1 ) 13.05 × 10-9 cm2/s (relative error 117%), (3) second iteration, ∆2 ) 0.425 Å (relative error 4.3%), D2 ) 7.80 × 10-9 cm2/s (relative error 40%), (4) third iteration, ∆3 ) 0.433 Å (relative error 1.9%, convergence achieved), D3 ) 9.90 × 10-9 cm2/s (relative error 27%). It is seen that, at the beginning of the SC iterations, values of the smearing factors between successive iterations differ considerably, sometimes even more than by 10% (while the corresponding error in the diffusivity values is 1 order of magnitude higher; see also Table 4). However, as the number of iterations increases, the relative error between such SC iterations decreases monotonically, and at the end it becomes comparable to the value reported for the error of the mean over different structures. This error (at convergence) is also equal to the difference between the mean values of the smearing factors for PET and PEI: ∆ ) 0.432 ( 0.008 Å for PET and ∆ ) 0.426 ( 0.005 Å for PEI. Reduction of the relative error to smaller levels (below 1%) was not possible for most samples studied with the Gusev-Suter TST approach. This is an extra reason sampling over a large enough number of statistically uncorrelated, well-equilibrated structures for the two copolyesters is needed. In the results presented in this work, we have averaged over 10 available PET and 10 PEI structures, to produce reliable estimates of penetrant diffusivity in the two polymers with the Gusev-Suter TST method. Results from such TST calculations are summarized and compared to experimental data in Tables 2 and 4. The data in Table 2 refer to diffusivity values calculated by applying TST on PET and PEI structures equilibrated with the Hofmann et al. cycle, presented in Table 1. In these structures the density was not adjusted; it was calculated as a byproduct of the simulation runs (see relevant description in section 3). In contrast, the data presented in Table 4 refer to diffusivity values calculated by applying TST on PET and PEI structures equilibrated with the new extended equilibration cycle discussed in section 3, in which the density was always adjusted to match experimentally measured values. As shown in Table 2, although pre-



dicted and experimental density values are quite close (the error being below 2% for both systems; see section 4.1), the qualitative comparison of diffusivities based on the predicted densities with the Hofmann et al. scheme is not satisfactory: As a result of the subtle but very important discrepancy in the densities (which is apparently amplified by the computations of accessible volumes and their distribution throughout the material samples), predicted O2 diffusivities are off by a factor of 6 in the case of PET and by 1 order of magnitude in the case of PEI compared to experimental data.87,88 As if this were not enough, they suggest that O2 diffuses faster in PEΙ than in PEΤ, which is obviously incorrect. Erroneous density values are considered as the primary reason for the discrepancy between predicted and measured O2 diffusivities in PET and PEI listed in Table 2. A second possible source of error is incomplete equilibration afforded by the Hofmann et al. cycle. Thus, in Table 4, a second set of results are reported for O2 diffusivity in PET and PEI glasses (T ) 300 K), obtained by applying TST on structures fully relaxed with the extended equilibration cycle reported in section 3. As already pointed out, in these structures the density was always imposed to match the experimentally known values for amorphous PET and PEI. This second set of predictions for O2 diffusivity in PET and PEI glasses are clearly in excellent agreement with the measured values: (a) Predicted and experimental values are remarkably close to each other. (b) Predicted values for O2 diffusivity in PET are systematically higher than in PEI, suggesting that PEI is a better barrier for O2. The ratio of O2 diffusivities in PEI and PET is predicted to be 1.8, remarkably close to the measured values87,88 of 2-2.5. This excellent qualitative and quantitative agreement should be attributed to (a) the exhaustive equilibration of PET and PEI structures at the temperatures for which results have been presented, thanks to the new extended equilibration cycle introduced for the first time here, (b) the use of an accurate atomistic force field (PCFF), and (c) carrying out the comparison between predicted and measured values at exactly the same density levels. Despite the simple way in which the higher segmental mobility in PET relative to PEI enters the Gusev-Suter TST (through the calculation of a single parameter, the smearing factor), the theory is apparently able to discriminate correctly between the two polymers. Whether the Gusev-Suter approach remains equally successful for larger penetrants, such as CO2, in these and other polyesters will be the subject of future investigations. Despite the heuristic, computationally intensive methodology adopted for the generation of melt and glassy samples, the results of this work are very encouraging, because it is the first time that (a) the different barrier properties of PET and PEI are quantitatively explained in terms of the different local mobilities, primarily of phenyl rings, in the two polymers and (b) O2 diffusivity in PET and PEI glasses is calculated in excellent agreement with experiment. Past simulation efforts overestimated O2 diffusivity in PET by almost 3 orders of magnitude.25

temperatures of interest were generated through the successful implementation of an extended equilibration cycle, involving 19 different steps (energy minimizations and NVT-NPT MD relaxations) in total. Our simulation data proved that the two polyisomers are characterized by (practically) identical static properties: Torsion angle distributions, total radial distribution functions, and chain characteristic ratios were calculated to be similar in the two polymers. However, our results revealed a significant difference in the segmental dynamics between the two polymers, which is related to phenyl ring rearrangements: In PET, the phenyl ring was seen to undergo facile motions around the axis of rotation formed by the bonds connecting it to the rest of the chain, leading to fast decorrelation of the torsional angles about the latter bonds. Due to a cooperative effect, the faster local dynamics of the phenyl migrates throughout the PET repeat unit and is evidenced even in the glycol group. In PEI, on the other hand, where the phenyl ring possesses no natural in-plane axis of rotation as a result of its para configuration, phenyl ring motions are more sluggish. The rate at which corresponding torsion angles decorrelate was calculated to be 2-3 times slower. Differences in the local phenyl ring dynamics between PET and PEI were evidenced most clearly at the higher temperatures simulated (e.g., T ) 600 K). In the second part of this work, the Gusev-Suter TST method, employing a self-consistently calculated smearing parameter to account for motions of the polymer matrix atoms, was employed to predict O2 diffusivity in the two polyesters. Two sets of data were presented: The first referred to O2 diffusivities calculated by applying TST on structures preequilibrated with the Hofmann et al. cycle,27 characterized by the densities obtained at 1 atm through this preequilibration procedure; these densities differed from the experimental ones by 1-2%. The second set of predictions concerned diffusivities calculated by applying TST on structures preequilibrated with the extended equilibration cycle introduced for the first time in this work; the densities of these structures were set to exactly match experimental values. Only in the latter case we were able to obtain consistent diffusivity predictions for O2 in PET and PEI. The predicted ratio of O2 diffusivities in PEI and PET was 1.8; this compares very favorably with the experimentally measured values that lie between 2 and 2.5.87,88 One should note that the difference in diffusivity cannot be attributed solely to the difference in mass density between PET and PEI. Computer specimens generated by the new, extended equilibration scheme do not yield correct diffusivities, or the correct ratio thereof, if the difference in matrix dynamics is disregarded (i.e., if the same smearing factor is used in the self-consistent Gusev-Suter calculation). This, combined with the observed structural similarity between the two polymers, leads to the conclusion that the difference in segmental dynamics is most probably the major factor behind the difference in diffusivities. The methodology presented here is very general and can be applied to other PET-based copolymers, blends, and polymeric membranes. Current efforts focus on (a) modifying the TST approach to account for larger, nonspherical penetrant molecules, such as CO2, and (b) developing a robust simulation approach59 for equilibrating complex polymer systems through an imple-

5. Conclusions. Future Plans We have presented results from detailed, atomisticlevel simulations concerning the structural, conformational, dynamic, and barrier properties of two polyisomers, PET and PEI. Well-relaxed configurations at the


mentation of the advanced set of chain-connectivityaltering MC algorithms49-52 for coarse-grained polymer models. Acknowledgment. This work was supported by BPAmoco, Naperville, IL. We express our sincere gratitude to Dr. I. Dairanieh (BP-Amoco) for his continuous help and guidance throughout the whole project, Dr. N. Zacharopoulos, K. Doulas, and N. Vergadou for their help with the use of the Cerius 2 software, Drs. J. Golab and S. Sakellarides of BP-Amoco for many fruitful and stimulating discussions, and Prof. R. Snurr and the members of his research group for their hospitality during N.Ch.K.’s stay at Northwestern University, Evanston, IL. Interactions with Dr. Dieter Hofmann, Dr. Matthias Heuchel, and Mr. Martin Siegert in the context of the PERMOD project, funded by the European Commission (Project Number G5RD-CT-200000200) are deeply appreciated. References and Notes (1) Zhao, J.; Wang, J. J.; Li, C. X.; Fan, Q. R. Macromolecules 2002, 35, 3097. (2) Karayannidis, G. P.; Sideridou, I. D.; Zamboulis D. N.; Bikiaris, D. N.; Sakalis, A. J. J. Appl. Polym. Sci. 2000, 78, 200. (3) Llana, P. G.; Boyce, M. C. Polymer 1999, 40, 6729. (4) Huang, J. M.; Chu, P. P.; Chang, F. C. Polymer 2000, 41, 1741. (5) Yazdanian, M.; Ward, I. M.; Brody, H. Polymer 1985, 26, 1779. (6) Sekelik, D. J.; Stepanov, E. V.; Nazarenko, S.; Schiraldi, D.; Hiltner, A.; Baer, E. J. Polym. Sci., Polym. Phys. 1999, 37, 847. (7) Cho, C. K.; Kim, J. D.; Cho, K.; Park, C. E.; Lee, S. W.; Ree, M. J. Adhes. Sci. Technol. 2000, 14, 1131. (8) Light, R. R.; Seymour, R. W. Polym. Eng. Sci. 1982, 22, 857. (9) Lee, S. W.; Ree, M.; Park, C. E.; Jung, Y. K.; Park, C. S.; Jin, Y. S.; Bae, D. C. Polymer 1999, 40, 7137. (10) Hu, Y. S.; Liu, R. Y. F.; Zhang, L. Q.; Rogunova, M.; Schiraldi, D. A.; Nazarenko, S.; Hiltner, A.; Baer, E. Macromolecules 2002, 35, 7326. (11) Kit, K. M.; Schultz, J. M.; Gohil, R. M. Polym. Eng. Sci. 1995, 35, 680. (12) Brolly, J. B.; Bower, D. I.; Ward, I. M. J. Polym. Sci., Polym. Phys. 1996, 34, 769. (13) McDowell, C. C.; Partin, J. M.; Freeman, B. D.; McNeely, G. W. J. Membr. Sci. 1999, 163, 39. (14) Sakellarides, S. L. ANTEC 1997, 822. (15) Schiavone, R. J. ANTEC 2001, 2313. (16) Sakellarides, S. L. ANTEC 1996, 938. (17) Bravard, S. P.; Boyd R. H. Macromolecules 2003, 36, 741. (18) Tatsumi, T.; Ito, E.; Hayakawa, R. J. Polym. Sci., Polym. Phys. 1992, 30, 701. (19) Li, B.; Yu, J.; Lee, S.; Ree, M. Eur. Polym. J. 1999, 35, 1607. (20) Li, B.; Yu, J.; Lee, S.; Ree, M. Polymer 1999, 40, 5371. (21) Wu, T. M.; Chang, C. C.; Yu, T. L. J. Polym. Sci., Polym. Phys. 2000, 38, 2515. (22) Karayannidis, G. P.; Bikiaris, D. N.; Papageorgiou, G. Z.; Pastras, S. V. J. Appl. Polym. Sci. 2002, 86, 1931. (23) Tonelli, A. E. Polymer 2002, 43, 6069. (24) Tonelli, A. E. J. Polym. Sci., Polym. Phys. 2002, 40, 1254. (25) Shanks, R.; Pavel, D. Mol. Simul. 2002, 28, 939. (26) Paul, D. R.; Koros, W. J. J. Polym. Sci., Polym. Phys. 1976, 14, 675. Petropoulos, J. H. J. Membr. Sci. 1990, 53, 229. (27) Hofmann, D.; Fritz, L.; Ulbrich, J.; Schepers, C.; Bo¨hning, M. Macromol. Theory Simul. 2000, 9, 293. (28) Takeuchi, H. J. Chem. Phys. 1990, 93, 2062. (29) Greenfield, M. L.; Theodorou, D. N. Macromolecules 1993, 26, 5461. (30) Theodorou, D. N. In Diffusion in Polymers; Neogi P., Ed.; Marcel Dekker: New York, 1996; p 67. (31) Arizzi, S. Diffusion of small molecules in polymeric glasses: a modelling approach. Ph.D. Thesis, Massachusetts Institute of Technology, Boston, 1990. (32) Gusev, A. A.; Arizzi, S.; Suter, U. W. J. Chem. Phys. 1993, 99, 2221. Gusev, A. A.; Suter, U. W. J. Chem. Phys. 1993, 99, 2228.

Macromolecules, Vol. 37, No. 8, 2004 (33) Greenfield, M. L.; Theodorou, D. N. Macromolecules 1998, 31, 7068. (34) Gray-Weale, A. A.; Henchman, R. H.; Gilbert, R. G.; Greenfield, M. L.; Theodorou, D. N. Macromolecules 1997, 30, 7296. (35) Gusev, A. A.; Mu¨ller-Plathe, F.; van Gunsteren, W. F.; Suter, U. W. Adv. Polym. Sci. 1994, 116, 207. (36) Ha¨nggi, P.; Talkner, P.; Borkovec, M. Rev. Mod. Phys. 1990, 62, 251. (37) Greenfield, M. L.; Theodorou, D. N. Macromolecules 2001, 34, 8541. (38) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1991, 95, 8866. (39) Karayiannis, N. C.; Mavrantzas, V. G.; Theodorou, D. N. Chem. Eng. Sci. 2001, 56, 2789. (40) Kirkpatrick, S. Phys. Rev. Lett. 1971, 27, 1722. (41) Argyrakis, P.; Milchev, A.; Pereyra, V.; Kehr, K. W. Phys. Rev. E 1995, 52, 3623. (42) Hofmann, D.; Fritz, L.; Ulbrich, J.; Paul, D. Comput. Theor. Polym. Sci. 2000, 10, 419. (43) Hofmann, D.; Fritz, L.; Ulbrich, J.; Paul, D. Polymer 1997, 38, 6145. (44) Heuchel, M.; Hofmann, D. Desalination 2002, 144, 67. (45) Hofmann, D.; Heuchel, M.; Yampolskii, Y.; Khotimskii, V.; Shantarovich, V. Macromolecules 2002, 35, 2129. (46) Fried, J. R.; Ren, P. Comput. Theor. Polym. Sci. 2000, 10, 447. (47) Lo´pez-Gonza´lez, M.; Saiz, E.; Guzma´n, J.; Riande, E. J. Chem. Phys. 2001, 115, 6728. (48) Kucukpinar, E.; Doruker, P. Polymer 2003, 44, 3607. (49) Pant, P. V. K.; Theodorou, D. N. Macromolecules 1995, 28, 7224. (50) Mavrantzas, V. G.; Boone, T. D.; Zervopoulou, E.; Theodorou, D. N. Macromolecules 1999, 32, 5072. (51) Karayiannis, N. C.; Mavrantzas, V. G.; Theodorou, D. N. Phys. Rev. Lett. 2002, 88, 105503. Karayiannis, N. C.; Giannousaki, A. E.; Mavrantzas, V. G.; Theodorou, D. N. J. Chem. Phys. 2002, 117, 5465. (52) Karayiannis, N. C.; Giannousaki, A. E.; Mavrantzas, V. G. J. Chem. Phys. 2003, 118, 2451. (53) Coburn, J. C.; Boyd, R. H. Macromolecules 1986, 19, 2238. (54) Hedenqvist, M. S.; Bharadwaj, R.; Boyd, R. H. Macromolecules 1998, 31, 1556. (55) Boyd, S. U.; Boyd, R. H. Macromolecules 2001, 34, 7219. (56) Bharadwaj, R. K.; Boyd, R. H. Polymer 1999, 40, 4229. (57) Cail, J. I.; Stepto, R. F. T.; Taylor, D. J. R.; Jones, R. A.; Ward, I. M. Phys. Chem. Chem. Phys. 2000, 2, 4361. (58) Honeycutt, J. D. Comput. Theor. Polym. Sci. 1998, 8, 1. (59) Zacharopoulos, N.; Vergadou, N.; Theodorou, D. N. To be published. (60) Tscho¨p, W.; Kremer, K.; Batoulis, J.; Burger, T.; Hahn, O. Acta Polym. 1998, 49, 61. Tschop, W.; Kremer, K.; Hahn, O.; Batoulis, J.; Burger, T. Acta Polym. 1998, 49, 75. (61) Aoyagi, T.; Sawa, F.; Shoji, T.; Fukunaga, H.; Takimoto, J.; Doi, M. Comput. Phys. Commun. 2002, 145, 267. (62) Mu¨ller-Plathe, F. Chem. Phys. Chem. 2002, 3, 754. (63) Hofmann, D.; Ulbrich, J.; Fritsch, D.; Paul, D. Polymer 1996, 37, 4773. (64) Fritz, L.; Hofmann, D. Polymer 1997, 38, 1035. (65) Fritz, L.; Hofmann, D. Polymer 1998, 39, 2531. (66) Hofmann, D.; Fritz, L.; Paul, D. J. Membr. Sci. 1998, 144, 145. (67) www.accelrys.com/cerius2/; www.accelrys.com/mstudio/. (68) Theodorou, D. N.; Suter, U. W. Macromolecules 1985, 18, 1467. (69) Theodorou, D. N.; Suter, U. W. Macromolecules 1986, 19, 139. (70) Nose´, S. Mol. Phys. 1984, 52, 255. (71) Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (72) Andersen, H. C. J. Comput. Phys. 1983, 52, 24. (73) Sun, H.; Mumby, S. J.; Maple, J. R.; Hagler, A. T. J. Am. Chem. Soc. 1994, 116, 2978. (74) Sun, H. Macromolecules 1995, 28, 701. (75) Sun, H.; Ren, P.; Fried, J. R. Comput. Theor. Polym. Sci. 1998, 8, 229. (76) Sun, H. J. Phys. Chem. B 1998, 102, 7338. (77) Smith, G. D.; Boyd, R. H. Macromolecules 1990, 23, 1527. (78) Vacatello, M.; Avitabile, G.; Corradini, P.; Tuzi, A. J. Chem. Phys. 1980, 73, 548. Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638. Paul, W.; Yoon D. Y.; Smith G. D. J. Chem. Phys. 1995, 103, 1702. Martin, M. G.; Siepmann J. I. J. Phys. Chem. B 1998, 102, 2569.



(79) Nath, S. K.; Khare, R. J. J. Chem. Phys. 2001, 115, 10837. Bourasseau, E.; Haboudou, M.; Boutin, A.; Fuchs, A. H.; Ungerer, P. J. Chem. Phys. 2003, 118, 3020. (80) Tashiro, K.; Nimmanpipug, P.; Rangsiman, O. J. Phys. Chem. B 2002, 106, 12884; Heuchel, M.; Hofmann, D.; Pullumbi, P. Macromolecules 2004, 37, 201. Mavrantza, I.-E.; Prentzas, D.; Mavrantzas, V. G.; Galiotis, C. J. Chem. Phys. 2001, 115, 3937. (81) Youk, J. H.; Jo, W. H.; Yoo, D. I. Polym. Bull. 1997, 39, 257. (82) Jin, Y.; Boyd, R. H. J. Chem. Phys. 1998, 108, 9912. (83) Schubash, H. R.; Nagy, E.; Heise, B. Colloid Polym. Sci. 1981, 259, 789.

(84) Gupta, M. R.; Yeh, G. S. Y. J. Macromol. Sci. Phys. 1978, 15, 119. (85) Takeuchi, H.; Roe, R.-J. J. Chem. Phys. 1991, 94, 7446. (86) Takeuchi, H.; Roe, R.-J. J. Chem. Phys. 1991, 94, 7458. (87) Weinkauf, D. H.; Paul, D. R. J. Polym. Sci., Polym. Phys. 1992, 30, 837. (88) Michaels A. S.; Vieth, W. R.; Barrie, J. A. J. Appl. Phys. 1963, 34, 13.

MA0352577

Detailed Atomistic Simulation of the Segmental Dynamics and

Recommend Documents