Structure, Dimensions, and Entanglement Statistics of Long Linear

Dec 16, 2008 - Onset of Entanglements Revisited. Topological Analysis. C. Tzoumanekas , F. Lahmar , B. Rousseau and D. N. Theodorou. Macromolecules 0 ...
0 downloads 0 Views 698KB Size
Download PDF
442

J. Phys. Chem. B 2009, 113, 442–455

Structure, Dimensions, and Entanglement Statistics of Long Linear Polyethylene Chains Katerina Foteinopoulou, Nikos Ch. Karayiannis,* and Manuel Laso Institute for Optoelectronics and Microsystems (ISOM) and ETSII, UPM, Jose´ Gutie´rrez Abascal 2, E-28006 Madrid, Spain

Martin Kro¨ger Polymer Physics, Department of Materials, ETH Zu¨rich, Wolfgang-Pauli-Strasse 10, CH-8093 Zu¨rich, Switzerland ReceiVed: September 17, 2008; ReVised Manuscript ReceiVed: October 31, 2008

This work elucidates the effect of both temperature and molecular length on the conformational and structural properties as well as on the entanglement statistics of long amorphous, polydisperse, and molten linear polyethylene (PE). A large number of PE samples are modeled in atomistic detail, with average molecular lengths ranging from C24 up to C1000 over a wide range of temperatures in the interval of 300 e T e 600 K under constant pressure (P ) 1 atm). By employing enhanced chain-connectivity-altering moves, full-scale equilibration is achieved within modest computational time even for the longest molecules at ambient conditions. At a second stage, direct geometrical analysis is applied on all equilibrated polymer configurations providing the corresponding primitive paths and intermolecular entanglements. Simulation findings on the characteristic ratio, density, and atomic packing are in excellent agreement with available experimental data. The same holds for the calculated plateau modulus; simulation predicts 1.8 ( 0.1 MPa. Regarding the primitive path statistics, the average contour length and the number of entanglements are found to exhibit a simple exponentialtype dependency on temperature. For the polydisperse samples studied here, a superposition of Poissonians (often represented by a negative binomial) describes best the distribution of entanglements of the primitive paths. I. Introduction Polymer structure, dynamics, and rheology have been and still remain the focus of extensive experimental, theoretical, and modeling works.1-7 Recent advances in different fields of computer science enable us now to study, from first principles, the static properties and packing of macromolecular systems for length scales that range from the bond length up to the chain end-to-end distance and also the dynamical behavior of atomistic polymer melts for time scales that reach up to 10 µs. For the dynamical analysis of polymers at or beyond equilibrium, molecular dynamics (MD) is the obvious and most straightforward computational method4,6,8,9 since it is attractively simple in its core implementation and generally applicable over diverse chemical constitutions and molecular architectures. However, MD clearly lacks the ability to efficiently equilibrate entangled polymers of high molecular weight (MW) since chains exhibit slow reptation dynamics3 with the center-of-mass diffusivity, Dcm, scaling with chain length, N, as Dcm ∼ N-b (b ≈ 2.3).10,11 Even in a state-of-the-art and massively parallel implementation, MD requires computational times that exceed many months, or even years, to dynamically equilibrate melts consisting of chains with the correct statistics, free of finitesize effects, and lengths of up to C500 (where the number index denotes the number of carbon atoms along the backbone) for the chemically simplest synthetic polymer (polyethylene, PE).12 The sluggish motion and dynamics exhibited by long linear macromolecules is intimately linked to the presence of topological constraints (entanglements) whose existence and nontrivial * To whom correspondence [email protected].

should

be

addressed.

E-mail:

effect on slow polymer dynamics constitute a central piece of the reptation theory for melts and solutions of linear chains.1,3,13-15 Concerning static properties and equilibration of macromolecules, Monte Carlo (MC) methods are excellent competitors (and companions) to conventional MD for the robust configurational sampling of polymer models in either on-lattice or continuum simulations.4,6,16-18 One of the main ingredients for the success of MC in rapidly and nonphysically relaxing chemically simple macromolecular entities is its inherent ability to integrate moves (algorithms) to function at maximum efficiency for a specific polymer system. It is widely accepted that the long succession of uncorrelated, stochastically attempted conformational transitions (Markov chain)8 employed in these “smart” MC algorithms is able to explore very effectively the configurational space of dense macromolecular systems. Nowadays, bulk systems of macromolecules characterized by relatively simple chemical constitution like polyolefins or polydienes can be equilibrated at all length scales within modest computational time through the application of advanced, chainconnectivity-altering MC moves accompanied by a set of localized algorithms. Employing chain-connectivity-altering moves like the end-bridging (EB),19 the double-bridging (DB),20 and their intramolecular variants allowed rapid creation of wellequilibrated atomistic configurations of (i) linear systems, polydisperse21-23 or monodisperse20,23 PE melts with molecular lengths up to C1000, and (ii) nonlinear (branched) structures, H-shaped,24 triarm stars,25 and linear triarm blends26,27 of PE chains. In parallel, chain-connectivity-altering moves have been successfully applied on chemically more complex macromolecules like poly(ethylene oxide) (PEO) and poly(ethylene

10.1021/jp808287s CCC: $40.75  2009 American Chemical Society Published on Web 12/16/2008

Long Linear Polyethylene Chains terephthalate) (PET) either in atomistic representation (PEO)28 or in coarse-grained models (PET).29 Within a hierarchical modeling approach, the large MC-based trajectory of fully equilibrated and statistically uncorrelated atomistic configurations offers an excellent starting point for successive (i) MD simulations to study the local and global chain dynamics and relaxation spectra,12,30 (ii) transition-state theory (TST)-based simulations to predict the barrier properties to small penetrant molecules at conditions of infinite dilution,31 (iii) Gibbs ensemble MC simulations to study phase equilibria,32 (iv) systematic development of coarse-grained models (for example, based on the concept of super blobs33,34 or for simulations of dissipative particle dynamics (DPD) on soft-core spheres35), and (v) topological analysis of entanglement statistics and the corresponding identification of primitive paths based on geometrical minimization rules.36,37 Here, the accurate identification of the static and dynamic behavior of chain entanglements is of critical importance for the prediction of the dynamical and rheological properties of polymer melts and solutions at and beyond equilibrium. Numerous theoretical studies have proposed constitutive models to predict the rheological properties of entangled polymer melts.38 In particular, the predictions of the primitive path statistics can be connected with recently developed sliplink-based models.34,39-41 Since the original work of Everaers et al.,42 an increasing number of simulation studies has focused on the identification of chain entanglements and on the calculation of statistical properties that are related to primitive paths. These modeling efforts have been applied on several chemically different polymers modeled also at diverse levels of detail (from fully atomistic to lattice-based coarse-grained representations).29,36,37,42-55 Detailed reviews about the various methodological approaches employed for the extraction of entanglement statistics can be found in refs 56-59. Differences between the original annealing versus newer geometric approaches concerning efficiency and reliability have been directly tested and summarized in ref 52. Qualitative deviations are attributed to disentanglement that occurs during annealing. In this work, we employ a state-of-the-art geometric algorithm (Z1) which proceeds by transforming the physical picture of topological interchain constraints (as conceived by DoiEdwards)3 into a pure mathematical problem of identifying the shortest multiple disconnected path subject to geometrical constraints arising from the configuration of the corresponding atomistic system.36,43,52,56 A direct consequence of the specific mathematical formulation is that the Z1 algorithm provides as output the minimum average contour length, 〈Lpp〉 (where 〈 〉 denotes averaging over all primitive paths of a given configuration), an approximate but definition-independent quantity. In addition, by mapping the extracted interior nodes of each primitive path into topological kinks, the number of entanglements, Z, is returned along with the corresponding contour length, Lpp, for each primitive path. In contrast to the contour length, the number of kinks per PP is subject to the definition adopted; for Z1, entanglements can exist either as a lone kink assigned to a single chain (when the topological obstacle is a straight line) or as a pair of kinks each belonging to a different parent chain. This article presents results from a hierarchical modeling approach to study the effect of temperature on structure, dimensions, and entanglement statistics of atomistic polydisperse linear PE samples with molecular lengths ranging from C24 up to C1000. In the first stage, very long trajectories of statistically uncorrelated structures are generated through MC simulations

J. Phys. Chem. B, Vol. 113, No. 2, 2009 443 TABLE 1: Average Chain Length, Nav, Molecular Weight, MW, Polydispersity Index, I, and Total Number of Chains, Nch, for All Linear PE Systems Studied Here Nav

MW (g/mol)

I

Nch

Nav

MW (g/mol)

I

Nch

C24 C34 C48 C70 C78 C100 C122 C142

338 478 674 982 1094 1402 1710 1990

1.083 1.083 1.083 1.083 1.083 1.083 1.083 1.083

50 50 50 80 40 48 48 22

C174 C224 C270 C320 C400 C500 C1000

2438 3138 3782 4482 5602 7002 14002

1.083 1.083 1.053 1.053 1.053 1.053 1.053

32 32 32 32 16 16 8

employing the panoply of atomistic chain-connectivity-altering algorithms.19,20 In the second stage, direct topological analysis is applied on the MC-based trajectory (for all configurations over the whole range of molecular lengths and densities) using the Z1 algorithm43,56 to extract the corresponding primitive path statistics. The paper is organized as follows. Section II provides details about the modeled PE systems and the applied atomistic force field. Section III describes, briefly, the moves that form the MC scheme and presents results that document its computational efficiency in equilibrating the long- and short-range system conformations. Simulation data on the dependence of density, chain dimensions, and interatomic packing on temperature and chain length are shown in section IV. In section V, we report the primitive path statistics as a function of temperature and chain length. Finally, section VI, summarizes the main results of the present work and lists potential applications for future studies. II. Systems Studied: Force Field All simulated linear PE systems are polydisperse, characterized by a uniform distribution of molecular lengths in the closed interval [(1 - ∆)Nav, (1 + ∆)Nav], where Nav is the number average chain length and ∆ is the half-width of the uniform chain length distribution reduced by Nav. In the limit of long chains (N f ∞), the polydispersity index, I, is related to ∆ through I ) 1 + ∆2/3.21 In our simulations, for the whole temperature range, ∆ is set equal to 0.50 (I ) 1.083) for systems with Nav e C224 and to 0.40 (I ) 1.053) otherwise (Nav > C224). The average chain length, Nav, for the modeled PE samples along with the corresponding molecular weight, MW, polydispersity index, I, and the total number of chains, Nch, is reported in Table 1. All MC simulations are cast in the nNchPTµ* ensemble, where n is the total number of interacting sites, P is the pressure (P ) 1 atm), and µ* is the reduced chemical potential spectrum properly specified to generate a uniform distribution of chain lengths.19 Cubic cells are employed with periodic boundary conditions applied in all three dimensions. Original configurations at T ) 450 K were borrowed from NPT-MC trajectories of well-equilibrated, strictly monodisperse PE analogues simulated through the DB-IDR algorithm.20,23 Initial configurations for all other temperatures are further generated from progressive quenching (cooling and heating for lower and higher temperatures, respectively) of fully relaxed systems at T ) 450 K. After an initial equilibration period, which depends both on the applied temperature and on the size of the system, production runs are carried out with a record frequency of either 200 000 (T g 400 K) or 400 000 (T ) 300 or 350 K) MC steps. The total simulation time, reported in Table 2 for all systems, ranges from a minimum of 200 million MC steps for the C24 at T ) 600 K up to a maximum of 8 billion MC steps for the C1000 at T ) 300 K.

444 J. Phys. Chem. B, Vol. 113, No. 2, 2009

Foteinopoulou et al.

TABLE 2: Total Time (In Millions of Steps) for the MC Simulations on Linear PE Systems As a Function of Temperature, Ta T (K) Nav C24 C34 C48 C70 C78 C100 C122 C142 C174 C224 C270 C320 C400 C500 C1000

300

350

400

450

500

550

600

2800 (9.00) 2800 (14.1) 2800 (19.4) 5200 (41.5) 4800 (20.8) 4800 (29.8) 4600 (30.9) 6000 (16.1) 6400 (25.1) 6800 (29.7) 7200 (52.0) 7200 (50.7) 7200 (23.6) 7600 (32.4) 8000 (18.7)

2400 (3.00) 2400 (4.14) 2400 (5.20) 3200 (9.99) 3200 (5.67) 3600 (7.08) 4000 (7.50) 4000 (3.87) 4400 (5.64) 4400 (5.83) 4400 (11.0) 4400 (10.5) 4800 (5.42) 5200 (6.30) 5200 (3.30)

1000 (1.20) 1000 (1.63) 1000 (1.97) 1400 (3.58) 1200 (1.95) 1200 (2.46) 1200 (2.60) 1600 (1.50) 1800 (1.88) 1800 (1.94) 1800 (3.25) 1800 (3.28) 1800 (1.60) 1800 (1.77) 2000 (0.993)

800 (0.690) 800 (0.885) 800 (1.04) 1000 (1.75) 1000 (1.04) 1000 (1.19) 1000 (1.27) 1200 (0.825) 1600 (0.937) 1600 (0.978) 1600 (1.57) 1600 (1.49) 1600 (0.851) 1600 (0.990) 1800 (0.507)

600 (0.457) 600 (0.559) 1000 (0.623) 800 (0.993) 800 (0.589) 800 (0.704) 800 (0.757) 1000 (0.473) 1000 (0.568) 1000 (0.584) 1000 (0.815) 1200 (0.827) 1000 (0.497) 1000 (0.497) 1000 (0.376)

400 (0.376) 400 (0.414) 400 (0.472) 600 (0.676) 600 (0.455) 600 (0.515) 600 (0.525) 800 (0.376) 800 (0.431) 800 (0.433) 1000 (0.560) 1000 (0.677) 1000 (0.423) 1000 (0.403) 1000 (0.332)

200 (0.342) 200 (0.327) 200 (0.391) 400 (0.517) 400 (0.402) 400 (0.423) 400 (0.414) 600 (0.340) 600 (0.357) 600 (0.373) 600 (0.438) 600 (0.448) 800 (0.339) 600 (0.353) 600 (0.330)

a Numbers in parentheses denote the corresponding equilibration time, τdec (in millions of MC steps), of the autocorrelation function of the end-to-end unit vector, 〈u(t)u(0)〉.

The force field employed in the course of the present work is identical to the one used for the simulation of strictly monodisperse linear PE analogues at T ) 450 K.36 It adopts the united-atom (UA) representation where carbon atoms along with their bonded hydrogens are lumped into single interacting sites. Methylene (CH2) and methyl (CH3) groups are treated as equivalent sites with respect to all forms of bonded or nonbonded potentials. Bond lengths are kept fixed at an equilibrium value of 1.54 Å. A simple harmonic potential governs the bending angles with an equilibrium value set at θ0 ) 114°.60 The distribution of torsional (dihedral) angles is controlled by the nine term sum-of-cosines formula proposed by Toxvaerd.61 Finally, all nonbonded interactions are described by a typical 12-6 Lennard-Jones potential with the parameters defined in the TraPPE model.62 The 1-4 intramolecular interactions (between sites separated by three bonds) affecting dihedral conformations are solely controlled by the torsion angle potential. More details about the energetic functions and the corresponding parameters employed here can be found in refs 23 and 36. III. Monte Carlo Scheme Both the end-bridging (EB)19,21 and the double-bridging (DB)20,23 chain-connectivity-altering moves along with their intramolecular analogues [intramolecular-end-bridging (IEB)18 and intramolecular-double-rebridging (IDR)20,23] are incorporated in the MC scheme since all simulated PE systems are

Figure 1. Percentage acceptance rates (in logarithmic scale) as a function of applied temperature, T, for all moves participating in the general MC scheme. The reported values were obtained by setting the following parameters for the moves: (i) flip: maximum change in the angle of δφmax ) (10° (see ref 63), (ii) Con.Rot.: intramolecular rebridging version with δφmax ) (10° (see ref 19), and (iii) Vol.Fluc.: maximum isotropic change in the box length of δl ) (0.2Å. Experimental values for the glass transition temperature, Tg, of PE fall in the range of 150-200 K.

characterized by a finite degree of polydispersity. These moves are accompanied by a set of localized algorithms applied either on a single site [end-mer rotation,8 internal libration (flip),63 and reptation64,65] or on a succession of sites along the chain [concerted rotation, (ConRot)66]. Preliminary simulations revealed that the intermolecular variant of the reptation algorithm,67 while preserving the same acceptance rate, was more efficient than the conventional pattern. On the basis of this finding, conventional reptation64,65 was replaced by the intermolecular analogue in all MC schemes.18 In addition, it was found that the chain-connectivity-altering moves that entail the construction of two trimer bridges (i.e., DB and IDR) exhibited very low acceptance rates for temperatures equal to and below T ) 400 K. Consequently, in this temperature range, DB and IDR are removed from the MC suite, and their portions are equally redistributed between the remaining single-bridge analogues (EB and IEB). For the general scheme, that is, the one that is used for the whole temperature range, the following mix of moves was employed: rotation (3%), flip (6%), intermolecular reptation (23%), concerted rotation (20%), endbridging (33%), and intermolecular end-bridging (13%), where numbers in parentheses denote percentage attempt probabilities. Since all simulations are cast in the semigrand statistical ensemble (nNchPTµ*), the remaining portion (2%) is allocated to (isotropic) volume fluctuations to sample over the density of the system. In order to facilitate comparisons, we use the same attempt probabilities for each participating move in all simulations regardless of the simulated (average) chain length and the applied temperature. Figure 1 presents the acceptance probabilities (in logarithmic scale) for all moves participating in the general MC scheme. It is evident that a temperature drop of even 300 K, accompanied by an increase of the total mass density of the system, affects only marginally the acceptance rates and performance of the standard localized MC moves (flip, rotation, and reptation). In contrast, the chain-connectivityaltering moves (EB and IEB) involving reconnection of two sites through a trimer bridge are significantly affected; their acceptance rates drop by a factor of ∼40. Even worse, in the

Long Linear Polyethylene Chains

Figure 2. Evolution of the orientational autocorrelation function of the end-to-end unit vector, 〈u(t) · u(0)〉, as a function of MC steps (in millions) for the linear polydisperse C1000 PE system in the whole temperature range studied here.

case of the DB and IDR moves, where two trimer bridges are formed, acceptance rates drop by more than 2 orders of magnitude, rendering both moves practically inefficient at low temperatures. Even with an acceptance probability of approximately 0.02% at standard conditions, EB remains very powerful to guarantee long-range equilibration within modest CPU time. For comparison, the more complicated DB and IDR algorithms (involving the displacement of twice as many atoms) at T ) 450 K do not exceed a single successful transition within 80 000 trials. Still, both moves are very efficient and vastly superior to conventional techniques20 for the rapid equilibration of the computationally more demanding, strictly monodisperse PE analogues of the same MW (C1000). Moving to even lower temperatures (not shown in the present manuscript) around and below the glass transition temperature (Tg), the acceptance rates of the localized moves, along with those of the chainconnectivity-altering algorithms, get significantly reduced. A solution to this performance issue of particular use for dense polymer packings can be provided by casting all localized moves to either an “adaptive bias”68 or “min-map bias”69 pattern. The computational efficiency of various MC schemes incorporating chain-connectivity-altering moves in equilibrating, at all length scales, atomistic PE melts of both linear and branched chain architectures has been well established previously.12,18,20,21,23,24 MC schemes have also been employed efficiently in a large temperature range; for example, simulations have been reported in the temperature range of 250-600 K for C78 and C1000 systems to study the PE barrier properties31 and in the range of 150-600 K for bulk and grafted C156 systems to determine the glass transition temperature of PE.70 In the present work, we examine the performance of the employed MC scheme both with respect to the applied temperature (from 600 down to 300 K) and to the employed chain lengths, emphasizing the longchain systems. The long-range equilibration time of the PE chains can be quantified by the evolution of the orientational autocorrelation function of the end-to-end unit vector, 〈u(t) · u(0)〉, where 〈 〉 denotes averaging over all chains and recorded configurations. As the chain system progressively loses the traces of its initial conformation, 〈u(t) · u(0)〉 decays from the original value of unity to zero. Accordingly, in a 〈u(t) · u(0)〉 versus t plot (where t is measured in MC steps), the faster the 〈u(t) · u(0)〉 curve drops to zero the faster the system equilibrates. Figure 2 presents the corresponding 〈u(t) · u(0)〉 versus t curves

J. Phys. Chem. B, Vol. 113, No. 2, 2009 445

Figure 3. Evolution of the orientational autocorrelation function of the end-to-end unit vector, 〈u(t) · u(0)〉, as a function of CPU time at T ) 300 K for linear polydisperse PE systems characterized by various average molecular lengths. The 〈u(t) · u(0)〉 versus t curves are scaled, where necessary, to correspond to the same total number of sites (8000) for all simulated systems. Also shown are corresponding 〈u(t) · u(0)〉 curves obtained from NPT MD simulations executed in parallel over 16 processors.

for the C1000 PE system as obtained for the whole temperature range studied here. As clearly seen, the long-range equilibration of the high-MW (and definitely entangled) C1000 PE system is reached very fast. For example, at T ) 450 K, less than 5 million MC steps are required (which is translated into about 3-4 h in real time on a typical modern processor). Not surprisingly, the equilibration time increases as the temperature is reduced, and approximately 100 million MC steps are needed for full longrange decorrelation at T ) 300 K; still, the required CPU time remains very reasonable (about 3-4 days in real time). For comparison, although state-of-the-art and massively parallel implementations of MD codes running on modern supercomputers allow us to access simulation times on the order of many microseconds, the long-range equilibration of an atomistic C1000 system at standard conditions remains well beyond the capabilities of any MD simulation. Figure 3 shows the decay of 〈u(t) · u(0)〉 as a function of CPU time for selected chain lengths at T ) 300 K. All MC-based CPU times are properly scaled to correspond to simulation runs on an Intel Pentium Xeon at 1.8 GHz bearing the same total number of interacting sites. Also shown in the same plot are the corresponding relaxation curves for the C78, C142, and C500 systems as obtained from NPT MD simulations starting from exactly the same initial configuration as the MC-based ones. All MD simulations were executed through the LAMMPS software package (version 2001)71 in parallel, with each run spreading over a total of 16 processors. The chain length dependence of the long-range equilibration, exhibited here from the MC simulations based on chain-connectivity-altering algorithms, is qualitatively identical to the unique behavior reported for poly-21,23 and monodisperse20 PE melts at higher temperatures; simulated under the same conditions, longer chains equilibrate faster than the shorter ones with respect to their longrange dimensions (an explanation can be found in the appendix D of ref 21). By comparing the 〈u(t) · u(0)〉-versus-t curves from MC and MD simulations at standard conditions, it is immediately apparent that MC is significantly faster than MD for the long-range equilibration even for the short C78 chains. In addition, it is shown that as the MW of the chains increases, so

446 J. Phys. Chem. B, Vol. 113, No. 2, 2009

Foteinopoulou et al. TABLE 3: Optimal Parameters a0(T) and G∞(T) of the Hyperbolic Expression (Equation 1) Used to Fit the Simulation Data for the Dependence of Mass Density, G(N,T), on Molecular Length, N, at a Given Temperature, T (Figure 4) along with Corresponding Values for the Correlation Coefficient, rfit, for Each Hyperbolic Fit temperature (K)

R0(T)

F∞(T) g/cm3

rfit

300 350 400 450 500 550 600

2.552 2.760 3.152 3.576 4.076 4.746 5.615

0.8642 0.8350 0.8089 0.7826 0.7567 0.7326 0.7093

0.9984 0.9977 0.9974 0.9976 0.9952 0.9933 0.9899

dependence on chain length is accurately described by a classical hyperbolic expression of the form Figure 4. Mass density as a function of average chain length (C24 e N e C1000) from MC simulations for the whole temperature range studied here (300 e T e 600 K). Also shown are lines as obtained from fittings with hyperbolic expressions at each temperature.

does the performance gap between MC and MD. For example, for the C500 (T ) 300 K), MC is by at least four orders of magnitude more efficient than MD in equilibrating chain dimensions. The decorrelation time, τdec, of the orientational autocorrelation function 〈u(t) · u(0)〉 is calculated from the corresponding relaxation curves (like the ones shown in Figures 2 and 3) as the integral τdec ) ∫〈u(t) · u(0)〉dt. The whole set of the calculated values of τdec (in MC steps) is listed in parentheses in Table 2 as obtained for each simulated system without any scaling associated with the total number of interacting sites. These data can be used to directly compare the total simulation time and the corresponding τdec (also reported in Table 2), providing an accurate estimation of the number of the statistically uncorrelated configurations that can be used for the calculation of chain dimensions as well as for the extraction of the corresponding primitive paths. IV. Results from Atomistic Simulations IV.1. Torsion Angle Distribution. Results for the torsion angle distribution (not shown here) agree with those reported in previous works.70 It is important to note that as the temperature drops, (i) the majority of torsion angles assumes conformations that lie in the vicinities of either the trans or gauche( regimes as there exist no angles with values at around (60° due to high-energy barriers and (ii) a gradual, but still significant, increase of trans conformations [-40o e φ e 40o] is observed as the trans population increases from approximately 54 at T ) 600 K to 70% at standard conditions (T ) 300 K). As expected (further discussed in section IV.3), the aforementioned increase of the trans population of dihedral angles is the main factor for the considerable expansion of chain dimensions of PE with decreasing temperature. IV.2. Mass Density. Since all present simulations are cast in the isobaric semigrand statistical ensemble, the volume fluctuations of the simulation cell in the equilibrated part of the MC trajectory allow for the accurate calculation of the mass density, F(N,T), of each modeled linear PE system. Figure 4 presents the mass density as a function of (average) chain length for all temperatures. For any given temperature, the density

F(N, T) )

F∞(T) a0(T) 1+ N

(1)

where F∞(T) is the density in the limit of infinite molecular weight at a given temperature, T. As can be seen in Figure 4, all simulation-based curves agree remarkably well with the hyperbolic expression, eq 1. Corresponding F∞(T) and a0(T) parameters for each simulated temperature are reported in Table 3 along with correlation coefficients, rfit, of the hyperbolic fits. The simulated infinite MW density at standard conditions F∞(T ) 300 K) ) 0.864 g/cm3 is in very good agreement (the relative error being approximately 1%) with experimental values reported in the literature. Specifically, the corresponding experimental values are (i) 0.854-0.855 g/cm3 by Mattozzi et al.72 obtained from extrapolations of available experimental data,73,74 (ii) 0.855 g/cm3 by Flaconneche at al.75 after ref 76, and (iii) 0.860-0.863 g/cm3 for PEB-x samples (where x denotes the number of ethyls branched per 100 backbone carbons) by Fetters et al.77 Very good agreement is also observed at T ) 450 K; for the C24 system, the simulation finding of F ) 0.679 g/cm3 lies within less than 3% from the reported experimental value of 0.698 g/cm3.78 At infinite chain length, the absolute error between our density value of 0.783 g/cm3 and the experimental one78 of 0.766 g/cm3 is as low as 2%. Also, our density calculations for the intermediate C142 or C174 systems are consistent with those reported by Alexiadis et al.70 for C156. Figure 5 illustrates the effect of applied temperature on mass density for PE systems characterized by various chain lengths as well as for the simulation-based extrapolations at infinite chain length, F∞(T). Perfect fittings (rfit e -0.9995) for all systems studied imply a linear dependence (in the form of reduction) of the density on temperature. Furthermore, by calculating the slopes of the simulation data shown in Figure 5, we estimate the isobaric thermal expansion coefficient, a, defined by

a)

( )

j 1 ∂V j ∂T V

P

(2)

j is the specific volume (reciprocal of mass density). where V From the slopes of the density curves shown in Figure 5, the isobaric thermal expansion coefficient at T ) 400 K is found to be equal to a ) 7.2 and 6.9 (10-4 K-1) for C1000 and for N f ∞, respectively, in very good agreement with the experimental

Long Linear Polyethylene Chains

J. Phys. Chem. B, Vol. 113, No. 2, 2009 447

Figure 5. Mass density as a function of temperature for linear polydisperse PE systems characterized by various average chain lengths. Also shown is the corresponding curve of the simulation-based estimation of the mass density at infinite chain length, F∞(T) (see also Table 3).

data range of 6.96-7.38 (10-4 K-1) reported by Han et al.79 for various metallocene-produced, ethylene-based copolymers at T ) 394 K. In an analogous approach, isothermal compressibility, β, is calculated via

β)-

( )

j 1 ∂V j ∂P V

(3)

T

j on P can be deducted from a series where the dependence of V of NVT simulations at varied pressures. Since all MC simulations reported here are cast in the semigrand ensemble at constant pressure, we can alternatively calculate the isothermal compressibility from volume fluctuations of any given NPT MC simulation as8

β)

1 〈V2〉 - 〈V〉2 kBT 〈V〉

) ( ∂U T( ∂P∂T ) ∂V ) T

U

-P)

for all chain lengths, as documented in section III, allows the robust sampling of the chain dimensions as quantified by the mean square end-to-end distance, 〈R2〉, and the mean square radius of gyration, 〈Rg2〉. In the limit of large N, our macromolecular chains behave as Gaussian coils with 〈R2〉 ) 6〈Rg2〉 and 〈R2〉 ∼ N2V with Flory exponent V ) 0.5 to high precision. Figure 6 shows the dependence of 〈R2〉 (the nonindexed brackets denote averaging over all chains within a given sample) on N over the whole range of simulated systems. For C1000, for the ratio 〈R2〉/ MW, we obtain 1.38 and 1.47 Å2 mol g-1 at T ) 450 and 400 K, respectively, both close to the reported literature value of 1.25 Å2 mol g-1 at T ) 413 K.77 In addition, the square root of the mean square radius of gyration of the C1000 at T ) 400 K is found equal to 〈Rg2〉1/2 ) 58.6 Å, which agrees exceptionally well with the experimental value of 56.9 Å (T ) 413 K) proposed by Horton et al.81 Having established the linear dependence of chain dimensions on molecular length, we have access to the characteristic ratio, C(N), as

(4) C(N) )

where kB is Boltzmann’s constant. By calculating the volume fluctuations of the simulation cell for the C1000 system at T ) 400 K, we obtain β ) 1.92 ((0.40) × 10-3 MPa-1, in good agreement with the available experimental data range of 0.8-0.9 (×10-3 MPa-1) for various metallocene-produced, ethylenebased copolymers at T ) 394 K.79 The solubility parameter, δ, is defined as80

δ)

Figure 6. Mean square end-to-end distance, 〈R2〉, as a function of average chain length N as obtained from MC simulations for the whole range of applied temperatures (300 e T e 600 K).

T( βa ) - P

(5)

using the isobaric thermal expansion coefficient, a, and the isothermal compressibility, β, readily calculated from simulation data through eqs 2 and 3, respectively. For the linear C1000 at T ) 400 K, on the basis of the values of a and β reported earlier, we find δ ) 12.2 MPa1/2 to be compared with δ ≈ 18 MPa1/2 as obtained from experiments on ethylene-based copolymers at T ) 394 K.79 IV.3. Chain Dimensions. The rapid long-range equilibration afforded by the advanced chain-connectivity-altering algorithms

〈R2〉 (N - 1)l2

(6)

where l is the CH2-CHx (x ) 2, 3) bond length, being fixed at 1.54 Å in all simulations. On the basis of the extracted 〈R2〉 values, the calculated characteristic ratios are shown in Figure 7 as a function of chain length for the whole range of applied temperatures. As seen for all temperatures, after an initial monotonic increase of C(N) with increasing length for longer chains (N > 200), the characteristic ratio reaches an almost constant value, which increases only marginally with N. If we extrapolate the simulation-based results to the limit of infinite length, or equivalently to 1/(N - 1) f 0 (see Figure 14 and eq 9 of ref 23 for details), we can compute for each temperature the characteristic ratio at infinite chain length, C∞(T). Alternatively, one could use the expression proposed by Mattice et al.82 for the correlation between C(N), N, and C∞, solving for C∞ since C(N) is readily available from our simulations

C(N) 1 1 )1+ C∞ C∞ N - 1

[ ( )] dC(N) 1 d N-1

(7) 1/(N-1)f0

448 J. Phys. Chem. B, Vol. 113, No. 2, 2009

Foteinopoulou et al.

Figure 7. Characteristic ratio, C(N), as a function of molecular length, N, as obtained from MC simulations for the whole range of applied temperatures 300 e T e 600 K. Also shown is the corresponding curve from best fit using eq 7 at T ) 300 K.

TABLE 4: Characteristic Ratio at Infinite Chain Length, C∞, and Persistence Length, ξ, for Linear PE at Various Temperatures As Extracted from Present MC Simulations temperature 300 K 350 K 400 K 450 K 500 K 550 K 600 K a

C∞ C∞b C∞c ξ (Å)

10.98 10.78 10.80 9.225

9.415 9.393 9.402 8.020

8.738 8.742 8.741 7.498

8.255 8.220 8.231 7.126

7.886 7.867 7.868 6.842

7.653 7.606 7.616 6.663

7.361 7.323 7.335 6.438

a Using the extrapolation method of ref 23. b Using eq 7 after ref 82 (where the slope is calculated over the whole range of N). c From linear fits of the 〈R2〉 versus (N - 1) simulation data (see, for example, Figure 6).

Equation 7 appears to be valid even for relatively small N. Extrapolated values of C∞ using both approaches are given in Table 4. Additionally, the values for C∞ from a direct estimate from the slope of the 〈R2〉 versus (N - 1) curves at each temperature (i.e., Figure 6) are given in Table 4. According to ref 2, the characteristic ratio at infinite length of PE in the temperature range of 411-415 K, as deduced experimentally from intrinsic viscosities at various solvents, is approximately 6.6-6.8, a range of values which lays very far from our simulation findings of C∞(T ) 400 K) ) 8.738 and C∞(T ) 450 K) ) 8.255. Horton et al.81 estimate a characteristic ratio (T ) 413 K) of around 7.4, which is quite closer. In addition, very recent light scattering measurements on linear PE in diphenyl at T ) 400 K by Muraoka et al.83 determine characteristic ratios from 8.7 to 10.5, which clearly deviate form the established value of 6.7 and agree very well with our present findings. According to Flory, in the limit of N f ∞, the persistence length, ξ(T), at a given temperature can be calculated as2

ξ(T) )

[(

〈R2〉 l 2 (N - 1)l

)

Nf∞

]

l + 1 ) [C∞(T) + 1] (8) 2

Simulation-based values of the persistence length for the whole temperature range studied here are reported in Table 4. Figure 8 presents the natural logarithm of the mean square end-to-end distance, ln〈R2〉, as a function of temperature for the C24, C48, C100, and C1000 systems. Temperature coefficients of the linear PE can be obtained from the ratio d ln〈R2〉/dT, which

Figure 8. Natural logarithm of the mean square end-to-end distance, ln〈R2〉, as a function of temperature for the C24, C48, C100, and C1000 linear PE systems. The best fit on the C1000 system data in the temperature range of 400 e T e 500 K is also shown (solid line). Experimental values for the glass transition temperature, Tg, of PE fall in the range of 150-200 K.

is readily available as the slope of the simulation data of Figure 8. If we focus our attention to the temperature range of 400 e T e 500 K, linear regression yields a slope of d ln〈R2〉/dT ) -1.07 × 10-3 K-1 for the C1000 system, in excellent agreement with the corresponding temperature coefficient of -1.0((0.1) × 10-3 K-1 from experimental studies at 413 e T e 463 K.2,84 Notice that the trend exhibited by PE regarding the effect of temperature (chain expansion with decreasing temperature) is in contrast with the majority of polymers where chains tend to collapse with temperature drop. Not surprisingly, this specific behavior of chain dimensions has further implications concerning similar “anomalous” trends for some characteristic lengths related to the underlying primitive paths (PP). IV.4. Interatomic Packing. Local packing is quantified and analyzed in terms of the intermolecular pair distribution function, ginter(r), the intramolecular pair density function, w(r), and the total radial pair distribution function, g(r). These distributions are interrelated through21

g(r) ) ginter(r) +

n w(r) V

(9)

Figure 9 presents the dependence of the intermolecular pair distribution function ginter(r) on the radial distance, r, for the C1000 system at T ) 300, 450, and 600 K. Also shown are ginter(r) data for the C500 at T ) 300 K, which coincide with the corresponding curve of the C1000, suggesting that the MW (for large N) has no appreciable effect on local packing as quantified by ginter(r). As the temperature drops (and the density increases), a small shift is observed for the ginter(r) curves toward smaller distances, indicating that intermolecular neighbors are brought closer together because of dense packing. This space-filling behavior of the PE constituent atoms with increasing density is reminiscent of the one observed for coarser models of freely jointed hard-sphere chains at high densities.85 At small radial distances (r < 5Å), the “correlation hole” effect86 is apparent as the surrounding cloud of intramolecular neighbors (especially the bonded ones) around a reference atom prohibits the close proximity of intermolecular neighbors. It is only at ap-

Long Linear Polyethylene Chains

Figure 9. Intermolecular pair distribution function, ginter(r), as a function of the radial distance, r, for the C1000 system at T ) 300, 450, and 600 K. Also shown are corresponding simulation data for the C500 system at T ) 300 K.

J. Phys. Chem. B, Vol. 113, No. 2, 2009 449

Figure 11. Total radial pair distribution function, g(r), as a function of radial distance for the C1000 (T ) 300, 450, and 600 K) and C500 (T ) 300 K) systems as obtained from present MC simulations. Inset: Static structure factor, S(k), as a function of wavenumber, k, of the linear C1000 PE system at T ) 300, 450, and 600 K.

chain length (for large N). Again, at lower temperatures, the characteristic maxima and minima of the total radial pair distribution function become more pronounced, indicating a clear tendency of pairs of sites to assume characteristic distances (and conformations). The calculation of g(r) from simulations is very useful since it is correlated to the static structure factor, S(k), through a 1-D Fourier transform8

S(k) ) 1 +

Figure 10. Same as in Figure 9 but for the intramolecular pair density function, w(r). The increasing area under the curve with increasing T at larger r arises from differences in the volume with temperature.

proximately r ≈ 5.1 Å that the intensity of ginter(r) reaches the first peak, which, at ambient conditions (T ) 300 K), corresponds also to a global maximum. Figure 10 shows the dependence of the intramolecular pair density function, w(r), on the radial distance for the C1000 at T ) 300, 450 and 600 K. The (global) maximum of the first peak at r ) 1.54 Å corresponds to the characteristic distance of the (constant) bond length. The second (broader) peak corresponds to fluctuations of the characteristic 1-3 distance (between atoms two bonds apart) and can be trivially calculated as 2l sin(θ0/2) ≈ 2.58 Å. Successive local maxima at characteristic radii are attributed primarily to different conformations of the 1-4 torsion angles, with a preference toward the trans population. Not surprisingly, as the temperature drops, the relative intensities of the peaks (maxima) and depths become more pronounced, and at standard conditions, local ordering (as indicated by the presence of characteristic peaks) seems to exist for distances up to 9 Å. The total radial pair distribution function, g(r), is shown in Figure 11 for the C1000 system at T ) 300, 450, and 600 K (and for C500 at T ) 300 K). Since g(r) is the accumulation of the intra- and intermolecular contribution, it remains unaffected by

4πFN k

∫0∞ [g(r) - 1]r sin(kr)dr

(10)

where FN denotes the number density of the system (FN ) n/V). S(k) is readily available through X-ray scattering experiments, which allow a direct comparison with corresponding atomistic simulations. In the inset of Figure 11, the S(k) curves are presented as obtained from our MC simulations on the C1000 system at T ) 300, 450, and 600 K. Past simulations of strictly monodisperse PE samples at T ) 450 K reported an excellent agreement between simulation findings and experimental data for the whole range of wave numbers.20 Our present S(k) curve for C1000 at T ) 450 K coincides with the corresponding one from monodisperse PE model of the same MW, with the only exception of a small discrepancy in the intensity of the first peak (which reflects intermolecular interactions); the present value of around 2.1 is higher than the older one of approximately 1.8.20 This discrepancy can be solely attributed to the small difference (around 0.7%) in the value of the mass density as obtained from the equilibrated part of the MC trajectory. This is expected since the total simulation time of the present MC trajectories exceeds by almost 1 order of magnitude the ones reported in ref 20, ensuring better predictions from ensembleaveraging an extended set of configurations. V. Entanglement Statistics V.1. Average Values. We now proceed to the second stage of the hierarchical approach to extract the primitive paths (PP) of all configurations. To this end, we implement the geometrical algorithm as embedded in the “Z1” code, which proceeds by successively reducing the number of nodes and the contour length of the shortest path (entanglement) network until convergence. Details about the algorithm and its application on

450 J. Phys. Chem. B, Vol. 113, No. 2, 2009

Foteinopoulou et al. TABLE 5: Bond Length of the Shortened Path, bpp(T), Effective Tube Diameter, app(T), Entanglement Molecular Weight under the Assumption of a Random Coil of Equidistant Step Length, Ne,coil(T), and Plateau Modulus GN0 from Extrapolations to Parent Chains of Infinite Molecular Length and at Various Temperatures temperature (K) bpp (T)a (Å) app (T)b (Å) Ne,coil (T)c GN0 d (MPa) 300 350 400 450 500 550 600

Figure 12. Average contour length of the primitive paths, 〈Lpp〉, as a function of molecular length for the whole temperature range studied here (300 e T e 600 K). Also shown are lines from best fits at T ) 300, 400, 500, and 600 K for N g 224. The slope of the lines are used to construct bpp(T) (values reported in Table 5).

computer-generated models of diverse chemical detail can be found in refs 36, 43, 52, and 56. The computational code can be applied to configurations online at www.complexfluids.ethz.ch/cgi-bin/Z1. Alternatively, and also for atomistic or coarse-grained systems, one could use the “CReTA” algorithm developed by Tzoumanekas and Theodorou based on performing random aligning string moves on repeat units that are treated as hard spheres.29,37,53,57,87 It is well established that both topological algorithms, when applied on exactly the same system configurations, provide almost identical results (subject to the statistical error of the length minimization procedure) about the average contour length, 〈Lpp〉, the principal quantity of interest characterizing the PPs, and that both geometric methods outperform annealing approaches used in past works. Here, 〈Lpp〉 is a characteristic measure for the primitive paths, which is not subjected to ambiguities related to different definitions.36,37,52 Given the computational efficiency of the Z1 algorithm,88 all existing configurations of the MC-generated trajectories are subjected to topological analysis. Figure 12 presents the average contour length of the primitive paths, 〈Lpp〉, as a function of average molecular length of the parent atomistic PE chains for the whole range of temperatures. In all cases studied here, 〈Lpp〉 is seen to depend on N as

〈Lpp〉 ) bpp(T)(N - 1)Vpp

(11)

where the definition of the bond length of the PP, originally introduced by Everaers et al.,42 bpp(N,T) ≡ Lpp(N,T)/(N - 1), has been used. The quantity bpp(T) in eq 11 is the bond length of the shortest path in the limit of N f ∞, which, at a given temperature, is assumed to be constant over all distances between kinks. As can be clearly seen by the perfect agreement between best linear fits and the simulation data of Figure 12 for large N (N > 224), 〈Lpp〉 scales linearly with N(Vpp ) 1). The values for the bond length of the PP, bpp(T), are obtained from the slope of the Lpp(N,T) versus N data for N > 224. The corresponding values of the bond length, bpp(T), at all temperatures are reported in Table 5. Not surprisingly, bpp(T) is seen to increase with decreasing temperature in a trend reminiscent of the behavior of the characteristic ratio (and therefore also the Kuhn length) at infinite molecular weight (see Table 4).

0.5839 0.5405 0.5014 0.4678 0.4397 0.4127 0.3873

43.33 40.73 40.89 41.15 43.86 43.19 44.14

74.21 75.36 81.55 87.96e 99.75 104.65 113.97

1.66 1.84 1.88 1.90 1.80 1.83 1.77

a Obtained from the slope of the 〈Lpp〉 versus (N - 1) curves for N > 224 (see Figure 12). b Calculated as app(T) ) C∞(T)l2/bpp(T). c Calculated from Ne,coil ) app /bpp. d Obtained through GN0 ) 4F∞RT/ 5Me,coil, Me,coil being the entanglement molecular weight resulting from the reported Ne,coil values (fourth column). e Experimental value of Ne ≈ 82 (T ) 443 K) reported by Fetters et al.90

Figure 13. Natural logarithm of the average contour length of the primitive paths, ln〈Lpp〉, as a function of temperature, T, for parent chains with average molecular lengths of C48, C100, C320, C400, C500, and C1000. Also shown are lines as obtained from best fits for all chain lengths over all temperatures. Experimental values for the glass transition temperature, Tg, of PE fall in the range of 150-200 K.

Figure 13 shows the dependence of the natural logarithm of the average contour length of the primitive paths, ln〈Lpp〉, on temperature for atomistic PE systems of various average chain lengths. For all simulated systems, independent of molecular length, ln〈Lpp〉 scales linearly with T for the whole temperature range studied here, implying that 〈Lpp〉(T) ) L0 exp(cLppT), where L0 and cLpp are coefficients that are readily available from best linear fits. The cLpp parameter is analogous of the thermal coefficient but now refers to the dimensions (contour length) of the primitive paths. It is found to be equal to -1.22 and -1.32 (×10-3 K-1) for the C500 and C1000 systems, respectively, indicating a dependence of 〈Lpp〉 on T, which is independent of chain length and similar to, although slightly higher than, the corresponding one of 〈R2〉 on T (see section IV.2). We should note that the present finding cannot be immediately generalized to other systems as PE exhibits a negative temperature coefficient for chain dimensions, in contrast to the majority of polymers. Simulations for polymers of different chemical constitution are required to conclude about the universality of this finding, that is, a strong correlation between the thermal coefficients of 〈Lpp〉 and 〈R2〉.

Long Linear Polyethylene Chains

J. Phys. Chem. B, Vol. 113, No. 2, 2009 451

The effective tube diameter (or PP Kuhn segment37), under the assumption that primitive paths behave as Gaussian coils (i.e., in the large molecular weight regime), is defined by

app(N, T) ≡

〈R2(N, T)〉 〈Lpp(N, T)〉

(12)

In the limit of infinite chain length (N f ∞), simple calculations lead to app(T) ) C∞(T)l2/bpp(T); values of app(T) based on the corresponding C∞(T) and bpp(T) results are reported in Table 5. The value at T ) 450 K, calculated here for polydisperse PE melts at infinite chain length (app(450 K) ) 41.2 Å), is in perfect agreement with app ) 40 Å as extracted for monodisperse PE samples (for lengths up to C1000).36 In addition, our finding of app(500 K) ) 43.9 Å is an excellent match with the neutronspin echo measurements of Wischnewski et al.,89 suggesting that app(NSE, 509 K) ≈ 45 Å. The dependence of bpp and app on the molecular length of the parent atomistic chains at T ) 300, 450, and 600 K is shown in Figure 14. Values for app are shown for N g 70 since the definition (eq 12) does not make (physical) sense for smaller chains and is related to the tube diameter only in the large molecular weight (polymeric) regime. On the basis of our data, the decrease of bpp with increasing N is well captured by a modified stretched exponential decay

bpp(N, T) ) b∞pp(T) + b1(T)exp(-Nγ(T))

(13)

where γ(T), b1(T), and b∞pp(T) are a priori temperature-dependent fitting parameters. In more detail, b∞pp(T) corresponds to the bond length of the shortened path for chains of infinite length (N f ∞), and γ(T) corresponds to the “rate” at which the bpp(N,T) ∞ (T). In the versus N curve approaches the plateau value bpp limiting case of γ(T) ) 1, eq 13 reduces to a first-order exponential decay expression, while the closer the value of γ(T) is to zero, the longer the chain length required to reach asymptotic behavior. Finally, the exponential prefactor b1(T) corresponds to the amplitude of the decrease from the high values of small N to the asymptotic b∞pp(T) at N f ∞. All optimal parameters after application of eq 13 on the bpp versus N simulation data (similar to the ones shown in Figure 14) are listed in Table 6. Close inspection of the values of bpp from infinite N extrapolations, either from application of eq 13 on the bpp(N,T) versus N curves or from linear fitting of eq 11 on the 〈Lpp〉(N,T) versus (N - 1) ones, suggests that, within the statistical uncertainty, both approaches lead to very comparable ∞ (T) (values listed in Table 6) and bpp(T) (listed in values for bpp Table 5), respectively. Furthermore, γ(T) (≈0.30) remains almost unaffected by the applied temperature, indicating a decay ratio of bpp(N,T), which approaches a constant with increasing N. Similarly to bpp, the dependence of app on the average chain length can be fitted accurately by a modified stretched exponential (growth) formula

app(N, T) ) a∞pp(T) + a1(T)exp(-Nδ(T))

(14)

with a∞pp(T), a1(T), and δ(T) being the corresponding parameters. The physical interpretation of the fitting parameters is identical to those for bpp, with the only difference that a1 should adopt negative values since it describes an increase with chain length. The optimal values for the effective tube diameter, as obtained

Figure 14. Bond length, bpp, and effective tube diameter, app, of the shortened path as a function of chain length of the parent atomistic PE chains at T ) 300, 450, and 600 K. Graph shown in double-x format. Fittings for both bpp and app with the modified stretched exponential decay functions are also presented.

from present simulations, are reported in Table 6. The extrapo∞ , agree very well with lated values at infinite chain length, app the corresponding ones shown in Table 5. Within the statistical error, a1 and δ remain constant, independent of the applied temperature. Furthermore, the low values of δ suggest that the growth ratio with N is a much slower “process” than the typical first-order exponential growth. While the bond length of the shortened path (bpp(N,T)) increases monotonically with decreasing temperature, the effective tube diameter (app(N,T)) remains quite unaffected by temperature, as demonstrated in Figure 15. More precisely, for the long-chain systems (N g C320), app(N,T) shrinks by about 2-3 Å in the interval 500 e T e 600 K, remains almost constant at 350 e T e 450 K, and then, at ambient conditions (T ) 300 K), increases by about 4-6 Å. This behavior is understood by analyzing the temperature behavior of 〈R2〉 and 〈Lpp〉 reported above; while for a given N, d ln〈Lpp〉/dT remains constant (see Figure 13), d ln〈R2〉/dT assumes different characteristic values depending on the applied range of temperatures (see Figure 8). Again, we should state that this behavior is characteristic of PE, as one would usually expect app(N,T) to decrease appreciably with T when 〈R2〉 and 〈Lpp〉 shrink and expand, respectively. Figure 16 presents the average number of entanglements per primitive path, 〈Z〉, as a function of molecular length of the parent chain. At ambient conditions, 〈Z〉 is found to be equal to 17.0 and 33.2 for the C500 and C1000 systems, respectively. Interestingly, at the same temperature (T ) 300 K), PE oligomers with an average length as low as C24 bear, on average, one entanglement per PP. The molecular weight between entanglements is defined through Ne,coil(N,T) ) app(N,T)/bp42 p(N,T) by invoking the assumption of a random coil of equidistant step length. The calculated values for Ne,coil(T) at infinite chain length as obtained from the corresponding values of app(T) and bpp(T) are also reported in Table 5. Our value at T ) 450 K of Ne,coil(450 K) ≈ 88 is in very good agreement with the experimental data of 82 (T ) 443 K) by Fetters et al.90 Moreover, the molecular weight between entanglements,N∞e,coil, ∞ ∞ can also be calculated using the app and bpp values calculated by fittings of the app(N,T) and bpp(N,T) data to the stretched exponential growth functions (eqs 13 and 14). Results are

452 J. Phys. Chem. B, Vol. 113, No. 2, 2009

Foteinopoulou et al.

TABLE 6: Optimal Parameters of the Stretched Exponential Formulas for bpp(N,T) (Equation 13) and for app(N,T) (Equation 14) from Fittings on the Dependence of the Primitive Path Bond Length, bpp(N,T), and of the Effective Tube Diameter, app(N,T), on the Chain Length at Various Temperaturesa temperature (K)

b∞pp (Å)

b1 (Å)

γ

a∞pp (Å)

a1 (Å)

δ

∞ a Ne,coil

GN0,∞b

300 350 400 450 500 550 600

0.6052 0.5620 0.5245 0.4930 0.4602 0.4356 0.4205

3.540 3.602 3.796 3.986 3.909 4.018 4.599

0.3059 0.3039 0.3051 0.3070 0.2993 0.2994 0.3123

45.26 40.02 40.50 40.89 41.70 43.04 43.34

-254.2 -233.6 -222.3 -218.0 -220.3 -221.9 -237.6

0.2551 0.2726 0.2616 0.2548 0.2495 0.2417 0.2460

74.78 71.21 77.22 82.94c 90.61 98.81 103.19

1.65 1.95 1.99 2.02 1.98 1.94 1.96

a ∞ Using the a∞pp and b∞pp, predictions for the interentanglement distance Ne,coil and the plateau modulus GN0,∞ are also provided. a Calculated as ∞ ∞ ∞ ∞ Ne,coil ) a∞pp/b∞pp. b Calculated as GN0,∞ ) 4F∞RT/5Me,coil , Me,coil being the entanglement molecular weight resulting from the reported Ne,coil values (eighth column). c Experimental value of Ne ≈ 82 (T ) 443 K) reported by Fetters et al.90

Figure 15. Effective tube diameter, app, as a function of temperature for atomistic parent chains characterized by various average chain lengths. Experimental values for the glass transition temperature, Tg, of PE fall in the range of 150-200 K.

Figure 16. Average number of kinks per primitive path, 〈Z〉, as a function of molecular length of the parent atomistic PE chains at all temperatures. Also shown are lines as obtained from best fits on the simulation data for N g C224.

reported in Table 6. Using this methodology, we find that N∞e,coil(450 K) ) 82.94, which is very close to the experimental value. It is worthwhile mentioning that the presence of kinks, in an average sense, for small systems does not necessarily mean that the system is also “rheologically entangled” as it may exhibit

very small disentanglement times. In the present study, entanglements have a topological definition, but in order to be related to the original reptation picture (or its several modifications), information for the dynamics is also required. There is no one-to-one relationship between the entanglement and the rheologically relevant critical molecular weight at which the shear viscosity changes its scaling behavior. However, Everaers et al.42 established a connection between the plateau modulus and degree of entanglement. We calculate the plateau modulus as G0N ) 4FRT/5Me,coil14,91 at infinite N using the values of density F∞ given in Table 3 and subsequently calculate Me,coil through the infinite N values of Table 5 or the values from fittings with stretched exponential functions of Table 6. Our findings of GN0 ∼ 1.80 MPa or GN0,∞ ∼ 1.93 MPa (300 e T e 600 K) are close together, and they are both clearly supported by experimental data, which report values for G0N in the range of 1.75-2.0 MPa.91 The dependence of 〈Z〉 on T is entirely described by the exponential model with 〈Z〉(T) ) Z0 exp(cZT), where the Z0 and cZ (thermal coefficient for Z) are in agreement with the corresponding trend of the average contour length of the PPs. Z0 and cZ parameters are calculated from linear fits on the simulation data of 〈Z〉 as a function of T on a linear/naturallogarithmic scale. In particular, the Z thermal coefficient (cZ) appears to be unaffected by the molecular length of the atomistic chains. For example, for the C48 system, cZ ) -1.87 (×10-3 Å-1) coincides, within the statistical error, with the corresponding values of -1.82 and -1.89 (×10-3 Å-1) extracted for the C500 and C1000 systems, respectively. Alternatively, one can define the number of entanglements Zcoil using the random coil assumption as 〈Zcoil〉 ) 〈N/(N - 1)〉(〈Lpp2〉/〈R2〉 - 1).3 The dependence of the ratio between 〈Z〉 and 〈Zcoil〉 on molecular length is similar to the bpp(N,T) dependence depicted in Figure 14, approaching a plateau value in the high molecular length regime, N g 224. The plateau values depend slightly on temperature (〈Z〉 ≈ 2.8〈Zcoil〉 and 〈Z〉 ≈ 2.5〈Zcoil〉 at T ) 300 and 600 K, respectively). These proportionality coefficients quantify the important fact that the shortest path is far from being a random walk with constant step length. Their existence shows that orientational correlations are partially, but not completely, lost at entanglement points (as determined via Z1). V.2. Distributions. Since the Z1 algorithm provides all relevant information about the contour length and the number of kinks (internal nodes) for each individual primitive path, we can calculate not only the average values of the PP statistics but also the corresponding distributions and identify their dependence on N and T. In all cases shown, distributions are obtained over all primitive paths (irrespective of the molecular length of the parent atomistic chain) and over all (equilibrated) MC configurations.

Long Linear Polyethylene Chains The dependence of the distribution for a variety of PP statistics on chain length has been reported in ref 36 for atomistic monodisperse linear PE samples at T ) 450 K. Since we are able to confirm these results, we focus on the effect of temperature solely. Figure 17 presents the probability distributions of (a) Lpp(N,T), (b) app(N,T), (c) bpp(N,T), and (d) Z(N,T) of the primitive paths obtained from the atomistic C1000 system at T ) 300, 400, 500, and 600 K. Regarding Lpp (Figure 17a), the temperature drop not only leads to a general shift toward larger values (due to the presence of more topological constraints) but also significantly broadens the distribution of the contour length. For example, at T ) 300 K, the lower and higher parts of the Lpp distribution are characterized by values of around 300 and 900 Å, respectively. In sharp contrast, the distribution of the effective tube diameter, app, remains, within the statistical error, practically unaffected by any change in temperature, as clearly seen in Figure 17b. In addition, as a result of the shape of the app distribution, the most probable value for app (≈13-16 Å) remains much smaller than its mean value (40-44 Å). The app distribution can be accurately captured by a beta function 0 , w1, w2, w3, ab, ac), defined as (with parameters app

J. Phys. Chem. B, Vol. 113, No. 2, 2009 453

[ ( [ (

P(app) ) a0pp + ab 1 +

)

]

w2-1 w2 + w3 - 2 (app - ac) w1w2 - w1 w2 + w3 - 2 1(app - ac) w1w3 - w1

)

]

w3-1

The shape of the measured p(bpp) remains practically unaltered upon varying temperature since the distribution is only shifted to higher values with decreasing temperature (with a constant shifting factor for all of its parts). As seen in Figure 17c, the form of p(bpp) is perfectly captured by Gaussian distributions at all temperatures. The probability distribution for the number of entanglements per chain Z, for the polydisperse C1000 system, is obviously broader than the one for monodisperse samples of the same molecular weight (results presented in ref 36), with their variance being now larger than the 〈Z〉. For the monodisperse samples, it was shown that the distribution of Z agrees with the shape of the Poisson distribution P(Z) ) Pµ(n) ≡ µne-µ/ n!, with µ ) 〈Z〉 - 1 and n ) Z - 136 (more generally µ ) 〈Z〉 - k, n ) Z - k, k g 0). Here, we find that the probability distribution for Z for the simulated polydisperse samples (C1000) is well fitted by a superposition of Poisson distributions, itself often represented by a negative binomial distribution (Figure

Figure 17. Probability distribution functions of (a) the contour length, Lpp, (b) the effective tube diameter, app, (c) the bond length, bpp, and (d) the number of kinks, Z, of primitive paths from parent atomistic chains of a linear C1000 PE system at various temperatures. Also shown are curves from fits with the Gaussian expression for bpp (c) and a negative binomial for Z (d).

454 J. Phys. Chem. B, Vol. 113, No. 2, 2009 17d). For large values of 〈Z〉, superpositioned Poissonians and Gaussians become similar; however, a superposition of Poisson distributions describes better our results which exhibit nonvanishing skewness. The probability distributions of the number of entanglements of smaller polydisperse systems with 〈Z〉 < 10, with variance of the same order with 〈Z〉, may be again fitted well by a classical Poisson distribution (results not shown). VI. Conclusions and Outlook We have employed a hierarchical modeling approach to study the packing, dimensions, structure, and entanglement statistics of long, linear polydisperse polyethylene (PE) samples ranging from C24 up to C1000 over a wide variety of applied temperatures. In its first stage, we have equilibrated the atomistic configurations through application of an extended MC scheme based on powerful chain-connectivity-altering moves that are far superior to any available relaxation method. We have studied the dependence of the density, characteristic ratio, and static structure factor on chain length and temperature. Direct comparison with available experimental data reveals very good to excellent qualitative and quantitative agreement for all calculated quantities. In the second stage, a direct topological analysis has been applied to all equilibrated MC-based configurations in order to extract the configuration of the shortest disconnected path (entanglement network) and related quantities. We have found that its contour length, the number of entanglements, and the entanglement molecular weight of the primitive paths (PP) exhibit a simple exponential-type of dependence on temperature for all simulated systems, while the effective tube diameter remains practically unaltered upon temperature variation. Furthermore, present simulation findings provide insight on the effect of temperature on the shape and size of the distribution for each PP-related quantity. Computer-generated atomistic configurations of the PE systems, simulated in the present work, will serve as excellent starting points for successive modeling studies about the dynamics of the atomistic samples or of chemically simpler entities (blobs)33,34 that can be obtained from the atomistic counterparts through a systematic coarse-graining procedure. We are currently employing the present hierarchical approach (MC simulations and topological analysis) to handle random packings of long, freely jointed chains of tangent hard spheres for volume fractions ranging from dilute ones up to the close vicinity of the maximally random jammed (MRJ) state and to investigate the connection between entanglements (intermolecular constraints) and knots (intramolecular constraints) regarding their scaling dependence on packing density.68,85 Acknowledgment. We are indebted to Prof. Doros Theodorou and Dr. Christos Tzoumanekas (NTUA, Greece) for fruitful discussions about the available topological algorithms used to extract primitive paths. We are deeply grateful for financial support of the present project by the EC through Contracts G5RD-CT-2002-00720 (PMILS) and NMP3-CT2005-016375 (NSF-EC collaboration project MNIBS) and by CICYT through Contract MAT2005-25569-E. Generous allocations of computational time on the “Magerit” (CeSViMa, UPM, Spain) and “Hreidar/Gonzales” (ETH Zurich, Switzerland) supercomputing clusters and additional computational time on the “Palu” supercomputer of CSCS (Switzerland) are deeply appreciated. References and Notes (1) de Gennes, P.-G. Scaling Concepts in Polymer Physics; Cornell University: Ithaca, NY, 1979.

Foteinopoulou et al. (2) Flory, P. J. Statistical Mechanics of Chain Molecules; Hanser Verlag: Mu¨nchen, Germany, 1989. (3) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Clarendon: Oxford, U.K., 1986. (4) Binder, K., Ed. Monte Carlo and Molecular Dynamics Simulations in Polymer Science; Clarendon: Oxford, U.K., 1995. (5) Larson, R. G. The Structure and Rheology of Complex Fluids; Oxford University Press: New York, 1999. (6) Kotelyanskii, M.; Theodorou, D. N., Eds. Simulation Methods for Polymers; Marcel Dekker: New York, 2004. (7) Dealy, J. M.; Larson R. G. Structure and Rheology of Molten Polymers; Hanser: Munich, Germany, 2006. (8) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987. (9) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic: London, 2001. (10) Lodge, T. P. Phys. ReV. Lett. 1999, 83, 3218. (11) Wang, S.-Q. J. Polym. Sci., Polym. Phys. 2003, 41, 1589. (12) Karayiannis, N. C.; Mavrantzas, V. G. Macromolecules 2005, 38, 8583. (13) McLeish, T. C. B. AdV. Phys. 2002, 51, 1379. (14) Larson, R. G.; Sridhar, T.; Leal, L. G.; McKinley, G. H.; Likhtman, A. E.; McLeish, T. C. B. J. Rheol. 2003, 47, 809. (15) Milner, S. T. Macromolecules 2005, 38, 4929–4939. (16) Theodorou, D. N. In Bridging Time Scales: Molecular Simulations for the Next Decade; Nielaba, P., Mareschal, M., Ciccotti, G., Eds.; Springer: Berlin, Germany, 2002. (17) Theodorou, D. N. In Computer Simulations in Condensed Matter: From Materials to Chemical Biology; Ferrario, M., Ciccotti, G., Binder, K., Eds.; Springer: Berlin, Germany, 2006; Vol. 2. (18) Karayiannis, N. C.; Mavrantzas, V. G. In Multiscale Modelling of Polymer Properties (Computer Aided Chemical Engineering 22); Laso, M., Perpe`te, E., Eds.; Elsevier: Amsterdam, The Netherlands, 2006; pp 3167. (19) Pant, P. V. K.; Theodorou, D. N. Macromolecules 1995, 28, 7224. (20) Karayiannis, N. C.; Mavrantzas, V. G.; Theodorou, D. N. Phys. ReV. Lett. 2002, 88, 105503. (21) Mavrantzas, V. G.; Boone, T. D.; Zervopoulou, E.; Theodorou, D. N. Macromolecules 1999, 32, 5072. (22) Uhlherr, A.; Leak, S. J.; Adam, N. E.; Nyberg, P. E.; Doxastakis, M.; Mavrantzas, V. G.; Theodorou, D. N. Comput. Phys. Commun. 2002, 144, 1. (23) Karayiannis, N. C.; Giannousaki, A. E.; Mavrantzas, V. G.; Theodorou, D. N. J. Chem. Phys. 2002, 117, 5465. (24) Karayiannis, N. C.; Giannousaki, A. E.; Mavrantzas, V. G. J. Chem. Phys. 2003, 118, 2451. (25) Peristeras, L. D.; Economou, I. G.; Theodorou, D. N. Macromolecules 2005, 38, 386. (26) Peristeras, L. D.; Rissanou, A. N.; Economou, I. G.; Theodorou, D. N. Macromolecules 2007, 40, 2904. (27) Rissanou, A. N.; Peristeras, L. D.; Economou, I. G. Polymer 2007, 48, 3883. (28) Wick, C. D.; Theodorou, D. N. Macromolecules 2004, 37, 7026. (29) Kamio, K.; Moorthi, K.; Theodorou, D. N. Macromolecules 2007, 40, 710. (30) Harmandaris, V. A.; Mavrantzas, V. G.; Theodorou, D. N.; Kro¨ger, ¨ ttinger, H. C.; Vlassopoulos, D. Macromolecules 2003, M.; Ramirez, J.; O 36, 1376. (31) (a) Gestoso, P.; Karayiannis, N. C. In Multiscale Modelling of Polymer Properties (Computer Aided Chemical Engineering 22); Laso, M., Perpe`te, E., Eds.; Elsevier: Amsterdam, The Netherlands, 2006; pp 201240. (b) Gestoso, P.; Karayiannis, N. C. J. Phys. Chem. B 2008, 112, 5646. (32) Wick, C. D.; Siepmann, J. I.; Theodorou, D. N. J. Am. Chem. Soc. 2005, 127, 12338. (33) (a) Padding, J. T.; Briels, W. J. J. Chem. Phys. 2002, 117, 925. (b) Kindt, P.; Briels, W. J. J. Chem. Phys. 2007, 127, 134901. (34) (a) Curco´, D.; Alema´n, C. Chem. Phys. Lett. 2007, 436, 189. (b) Curco´, D.; Alema´n, C. J. Comput. Chem. 2007, 28, 1929. (35) Guerrault, X.; Rousseau, B.; Farago, J. J. Chem. Phys. 2004, 121, 6538. (36) Foteinopoulou, K.; Karayiannis, N. C.; Mavrantzas, V. G.; Kro¨ger, M. Macromolecules 2006, 39, 4207. (37) Tzoumanekas, C.; Theodorou, D. N. Macromolecules 2006, 39, 4592. (38) (a) Marrucci, G. J. Non-Newtonian Fluid Mech. 1996, 62, 279. (b) Mead, D. W.; Larson, R. G.; Doi, M. Macromolecules 1998, 31, 7895. (c) Milner, S. T.; McLeish, T. C. B. J. Rheol. 2001, 45, 539. (d) Likhtman, A. E.; McLeish, T. C. B. Macromolecules 2002, 35, 6332. ¨ ttinger, H. C. J. Rheol. 2000, 44, 1293, and (39) Fang, J.; Kro¨ger, M.; O references therein. (40) (a) Schieber, J. D.; Nair, D. M.; Kitkrailard, T. J. Rheol. 2007, 51, 1111. (b) Schieber, J. D. J. Chem. Phys. 2003, 118, 5162.

Long Linear Polyethylene Chains (41) (a) Masubushi, Y.; Ianniruberto, G.; Greco, F.; Marrucci, G. J. NonCryst. Solids 2006, 352, 5001. (b) Masubuchi, Y.; Ianniruberto, G.; Greco, F.; Marrucci, G. J. Non-Newtonian Fluid Mech. 2008, 149, 87. (c) Greco, F. Eur. Phys. J. E 2008, 25, 175. (42) Everaers, R.; Sukumaran, S. K.; Grest, G. S.; Svaneborg, C.; Sivasubramanian, A.; Kremer, K. Science 2004, 303, 823. (43) Kro¨ger, M. Comput. Phys. Commun. 2005, 168, 209. (44) Kremer, K.; Sukumaran, S. K.; Everaers, R.; Grest, G. S. Comput. Phys. Commun. 2005, 169, 75. (45) Sukumaran, S. K.; Grest, G. S.; Kremer, K.; Everaers, R. J. Polym. Sci., Part B: Polym. Phys. 2005, 43, 917. (46) Shanbhag, S.; Larson, R. G. Phys. ReV. Lett. 2005, 94, 076001. (47) Zhou, Q.; Larson, R. G. Macromolecules 2005, 38, 5761. (48) Leo´n, S.; van der Vegt, N.; Delle Site, L.; Kremer, K. Macromolecules 2005, 38, 8078. (49) Hoy, R. S.; Robbins, M. O. Phys. ReV. E. 2005, 72, 061802. (50) Hoy, R. S.; Robbins, M. O. J. Polym. Sci., Part B: Polym. Phys. 2006, 44, 3487. (51) Shanbhag, S.; Larson, R. G. Macromolecules 2006, 39, 2413. (52) Shanbhag, S.; Kro¨ger, M. Macromolecules 2007, 40, 2897. (53) Spyriouni, T.; Tzoumanekas, C.; Theodorou, D.; Mu¨ller-Plathe, F.; Milano, G. Macromolecules 2007, 40, 3876. (54) Vladkov, M.; Barrat, J.-L. Macromolecules 2007, 40, 3797. (55) Kim, J. M.; Keffer, D. J.; Kro¨ger, M.; Edwards, B. J. J. NonNewtonian Fluid Mech. 2008, 152, 168. (56) Kro¨ger, M. Models for Polymeric and Anisotropic Liquids; Lecture Notes in Physics 675; Springer: Berlin Germany, 2005. (57) Tzoumanekas, C.; Theodorou, D. N. Curr. Opin. Solid State Mater. Sci. 2006, 10, 61. (58) Hoy, R. S.; Grest, G. S. Macromolecules 2007, 40, 8389. (59) Larson, R. G.; Zhou, Q.; Shanbhag, S.; Park, S. J. AIChE J. 2007, 53, 542. (60) Van der Ploeg, P.; Berendsen, J. J. Chem. Phys. 1982, 76, 3271. (61) Toxvaerd, S. J. Chem. Phys. 1997, 107, 5197. (62) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569. (63) Mavrantzas, V. G.; Theodorou, D. N. Macromolecules 1999, 31, 6310. (64) Bishop, M.; Ceperley, D.; Frisch, H. L.; Kalos, M. J. Chem. Phys. 1980, 72, 3228. (65) Vacatello, M.; Avitabile, G.; Corradini, P.; Tuzi, A. J. Chem. Phys. 1980, 73, 548. (66) Dodd, L. R.; Boone, T. D.; Theodorou, D. N. Macromolecules 1999, 32, 5072. (67) In the intermolecular reptation, a randomly selected end-mer is excised from a randomly selected chain and is reconstructed as an end-mer in a randomly selected new chain (different than the original one). By construction, it is applicable only in polydisperse chain systems. (68) Karayiannis, N. C.; Laso, M. Macromolecules 2008, 41, 1537. (69) Laso, M.; Karayiannis, N. C.; Mu¨ller, M. J. Chem. Phys. 2006, 125, 164108.

J. Phys. Chem. B, Vol. 113, No. 2, 2009 455 (70) Alexiadis, O.; Mavrantzas, V. G.; Khare, R.; Beckers, J.; Baljon, A. R. C. Macromolecules 2008, 41, 987. (71) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (72) Mattozzi, A.; Hedenqvist, M. S.; Gedde, U. W. Polymer 2007, 48, 5174. (73) Allen, G.; Gee, G.; Wilson, G. J. Polymer 1960, 1, 456. (74) Olabisi, O.; Simha, R. Macromolecules 1975, 8, 206. (75) Flaconneche, B.; Martin, J.; Klopffer, M. H. Oil Gas Sci. Technol. 2001, 56, 261. (76) Runt, J. P. Encyclopedia of Polymers Science and Engineering; Wiley: New York, 1985; pp 482-519. (77) Fetters, L. J.; Lohse, D. J.; Richter, D.; Witten, T. A.; Zirkel, A. Macromolecules 1994, 27, 4639. (78) Pearson, D. S.; Ver Strate, G.; von Meerwall, E.; Schilling, F. C. Macromolecules 1987, 20, 1133. (79) Han, S. J.; Lohse, D. J.; Condo, P. D.; Sperling, L. H. J. Polym. Sci., Polym. Phys. 1999, 37, 2835. (80) Allen, G.; Gee, G.; Wilson, G. J. Polymer 1960, 1, 456. (81) Horton, J. C.; Squires, G. L.; Boothroyd, A. T.; Fetters, L. J.; Rennie, A. R.; Glinka, C. J.; Robinson, R. A. Macromolecules 1989, 22, 681. (82) Mattice, W. L.; Helfer, C. A.; Sokolov, A. P. Macromolecules 2003, 36, 9924. (83) Muraoka, Y.; Kamide, K.; Suzuki, H. Br. Polym. J. 2007, 15, 107. (84) Ciferri, A.; Hoeve, C. A. J.; Flory, P. J. J. Am. Chem. Soc. 1961, 83, 1015. (85) (a) Karayiannis, N. C.; Laso, M. Phys. ReV. Lett. 2008, 100, 050602. (b) Laso, M.; Karayiannis, N. C. J. Chem. Phys. 2008, 128, 174901. (c) Foteinopoulou, K.; Karayiannis, N. C.; Laso, M.; Kro¨ger, M.; Mansfield, M. L. Phys. ReV. Lett., in press. (d) Laso, M.; Karayiannis, N. C.; Foteinopoulou, K.; Kro¨ger, M.; Mansfield, M. L. Soft Matter, submitted for publication. (86) Schweizer, K. S.; Curro, J. G. Phys. ReV. Lett. 1987, 58, 246. (87) Theodorou, D. N. Chem. Eng. Sci. 2007, 62, 5697. (88) We should note that the new version of the Z1 code does not necessitate the use of supercells in the case of systems polluted by size effects when chain dimensions exceed the box size. In the present work, for the three larger systems (C400, C500, and C1000), chain dimensions exceed significantly the size of the simulation cell, but there is no effect on any static or volumetric property as revealed by simulations with a much larger simulation box size. (89) Wischewski, A.; Monkenbush, M.; Willner, L.; Richter, D.; Likhtman, A. C.; McLeish, T. C. B.; Farago, B. Phys. ReV. Lett. 2002, 88, 058301. (90) Fetters, L. J.; Lohse, D. J.; Milner, S. T.; Graessley, W. W. Macromolecules 1999, 32, 6847. (91) Liu, C.; He, J.; van Ruymbeke, E.; Keunings, R.; Bailly, C. Polymer 2006, 47, 4461.

JP808287S