Molecular Dynamics Simulations for the Description of Experimental

Jul 10, 2015 - Vaibhaw Kumar , C. Rebecca Locker , Pieter J. in 't Veld , and Gregory C. Rutledge. Macromolecules 2017 50 (3), 1206-1214. Abstract | F...
2 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/Macromolecules

Molecular Dynamics Simulations for the Description of Experimental Molecular Conformation, Melt Dynamics, and Phase Transitions in Polyethylene Javier Ramos,* Juan F. Vega, and Javier Martínez-Salazar Biophym, Departamento de Física Macromolecular, Instituto de Estructura de la Materia, IEM-CSIC, C/Serrano 113 bis, 28006 Madrid, Spain S Supporting Information *

ABSTRACT: Long molecular dynamics simulations of the melt dynamics, glass transition and nonisothermal crystallization of a C192 polyethylene model have been carried out. In this model, the molecules are sufficiently long to form entanglements in the melt and folds in the crystalline state. On the other hand, the molecules are short enough to enable the use of atomistic simulations on a large scale of time. Two force fields, widely used for polyethylene, are taken into account comparing the simulation results with a broad set of literature experimental data. Although both force fields are able to capture the general physics of the system, TraPPe-UA is in a better quantitative agreement with the experimental data. According with the simulation results some fundamental aspects of polyethylene physical parameters are discussed such as the characteristic ratio (Cn = 8.2 and 7.6 at 500 K, for TraPPe-UA and PYS force fields, respectively), the isothermal compressibility (α = 8.57 × 10−4 K−1), the static structure factor and the melt dynamics regimes corresponding to an entangled polymer. Furthermore, the simulated Tg (187.0 K) obtained for linear PE is in a very good agreement with the extrapolated Tg values (185−195 K) using the Gordon−Taylor equation. Finally, the simulation of the nonisothermal crystallization process supports the view of a mixed state of adjacent and nonadjacent re-entry model. The simulated two phase model reproduces very well the initial fold length expected for high supercoolings and the segregation of the system in ordered and disordered layers. The paper highlights the importance of combining simulation techniques with experimental data as a powerful means to explain the polymer physics.

1. INTRODUCTION Research activity in the field of macromolecular science requires the use of a variety of techniques. This is due to the intrinsic characteristics of polymeric materials, which reveal interesting physical phenomena in the length scale from subnanometer up to microns or, alternatively, in the time scale from picoseconds to years. Covering these wide length and time scales demands the combination of powerful experimental techniques and computational tools. At this respect major advances in computational power and methodologies can nowadays help to face classical problems, and to establish a shared knowledge between theorists and experimentalists. Melt dynamics, glass transition, and crystallization process represent three paradigms of polymer physics, of which important challenges remain to be solved. Among all synthetic polymeric materials, polyethylene (PE) represents a model system in the field. The advantage of this polymer is the simplicity of its chemical structure, and the extensive set of experimental data. A key component in molecular dynamics (MD) simulations is the description of forces among atoms and their movements relative to each other. This description is given by the so-called “force field” (FF), which can be assumed to be a set of energy © XXXX American Chemical Society

terms describing both intra- and intermolecular interactions among the atoms of the system. The FF primary consists of valence (stretch, bending and torsion terms) and nonbonded potential terms (van der Waals and Coulomb). An appropriate choice of the FF and its parameters is essential, in order to obtain realistic results from simulations. A large number of FFs for polymers are available in the literature, optimized for different systems and even for particular physical processes. In the case of PE two different FFs are being widely used. Paul, Yoon, and Smith proposed a FF, designated from now on as PYS, that was parametrized using experimental data and quantum calculations on short n-alkanes.1 PYS-FF has been shown to give a good description of physical properties and melt dynamics of PE melts.2−4 Rutledge’s group has recently used PYS-FF, with minor modifications.5−7 These authors have demonstrated that this version captures essential qualitative aspects of the initial stages of the crystallization of alkanes and short chains of linear PE. In particular, these authors have shown that homogeneous nucleation of PE models is feasible at Received: April 21, 2015 Revised: June 29, 2015

A

DOI: 10.1021/acs.macromol.5b00823 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules Table 1. Force Fields and Parameter Values for the Bonded and Non-Bonded Interactions Used in This Work Bond Potential: Ebond = kb(r − r0)2 PYSa ‑1

CHx−CHx

‑2

kb, kJ mol nm

TraPPE-UAb ‑1

ro, nm

146300

kb, kJ mol nm

‑2

0.154

PYSa ‑1

CHx−CH2−CHx

TraPPE-UAb

‑1

‑1

θo, deg

kθ, kJ mol rad 250.8

kθ, kJ mol rad

109 i Torsion Potential: Etorsion = ∑n=3 0 Ci cos (φ)

Co, kJ mol−1

C1, kJ mol−1

C2, kJ mol−1

‑1

114 TraPPE-UAb

C3, kJ mol−1

Co, kJ mol−1

C1, kJ mol−1

6.56 16.93 3.59 −27.09 8.39 16.77 Nonbonded Lennard-Jones: ELJ = 4εij[(σij/rij)12 − ((σij/rij)6)]; εij = (εiεj;)1/2σij = 1/2)(σi + σj) PYSa CH2 CH3 a

θo, deg

258.6

PYSa CHx−CH2−CH2−CHx

ro, nm

188100c

0.153 Angle Potential: Ebond = kb(θ − θ0)2

C2, kJ mol−1

C3, kJ mol−1

1.13

−26.29

TraPPE-UAb

εi kJ·mol−1

σi, nm

εi, kJ·mol‑1

σi, nm

0.468 0.468

0.401 0.401

0.382 0.814

0.395 0.375

Data adopted from ref 1. bData adopted from refs 10−12. cData adopted from ref 54.

discuss our simulation results and compare them with existing experimental data of static and dynamic properties of the melt. Finally we deal with results obtained for the simulations of the glass transition and crystallization processes and compare them with those experimentally obtained.

deep supercooling using MD simulations. Furthermore, they have characterized the crystal nucleus as a cylinder shape in order to calculate the thickness of the crystal as a function of the nucleus size.7 Gee, Lacevic, and Fried modified the parameters of the PYS-FF with the purpose of enhancing computational efficiency for the study of the crystallization process in polymers.8,9 However, in this case a major modification was done to the torsional potential increasing the trans to gauche energy barrier as compared to the bond torsional potential of the original PYS-FF. In this way, the chain stiffness increases, enhancing computational efficiency due to an acceleration of the process. This procedure allowed these authors to observe a forced crystallization at time scales accessible via large-scale MD simulations. The second FF is the so-called TraPPe-UA (transferable potentials for phase equilibria, united atom version), developed by Siepmann’s group and parametrized against fluid-phase equilibrium data of linear and branched alkanes.10−12 MD simulations using the TraPPe-UA FF have been extensively used for the study of static and dynamic properties of PE melts, including the effect of molecular length and short chain branching.13−31 Yamamoto et al. have also used this FF for the study of “induced” crystallization processes in polypropylene (PP),32 but as far as we know, it has not been yet employed to simulate the homogeneous crystallization process of PE. The objective of this contribution is to examine fundamental physical processes in C192H386 (from now C192) polymer model using both PYS and TraPPe-UA FFs. We mainly focus our attention on the ability of each FF to simulate complex phenomena as melt dynamics, glass transition and crystallization process. There are various reasons for choosing C192 model: (i) the entanglement length of PE is about 60−160 carbons as determined by experiments and simulations26−30,33−38 (ii) C192 is long enough to display chain folded crystals;39−43 and (iii) there are experimental results on static properties and on melt dynamics of model PE samples with similar molecular weight44 and of the crystallization and melting behavior of similar n-alkanes and PE samples.39−43 The paper is organized as follows. We first present the computational methods and the simulation details. Then we

2. COMPUTATIONAL METHODS 2.1. Minimization and Molecular Dynamics Details. All minimizations and MD simulations were performed with GROMACS 4.6.x package45 using Metrocubo Graphical Processor Units (GPU) workstations46 and CESGA supercomputing facilities. A steepest descent algorithm was used for energy minimizations during 10000 steps. For MD simulations, the integration of the motion equations were performed through the velocity Verlet algorithm with a time step of Δt = 2 fs. The nonbonded interactions were calculated using the Verlet buffer cutoff scheme implemented in GROMACS 4.6 using a van der Waals cutoff of 1.4 nm. Simulations were performed on either the canonical (NVT) or the isothermal−isobaric (NPT) ensembles. The Nose−Hoover thermostat with a relaxation time (τT) of 0.5 ps was used to control the temperature.47 The pressure was set up to 1 atm in all NPT simulations and isotropically controlled with the Parrinello−Rahman barostat with a time constant (τP) for coupling of 5 ps.48 Periodic boundary conditions (PBC) were applied at all of the three directions to remove surface effect and to mimic bulk state. Unless otherwise stated, the structures were dumped every 20 ps for subsequent analysis using “in-house” codes. Molecular graphics have been generated using the VMD package.49 2.2. System Description and Preparation of Initial Melt Structures. The simulation boxes involve 75 chains of 192 CHx (where x = 2 for middle and 3 for terminal groups) united atom (UA) units. Thus, the hydrogen atoms were not explicitly taken into account in our model. The initial structures for MD were prepared as follows. The generation and packing of the 75 alkane chains were initially performed using a modified recoil growth algorithm with periodic boundary conditions.26 The initial melt density was set up to 0.740 g/ cm3 which is close to the experimental density at 473 K.50 Subsequently, the chains in the simulation box were minimized with a fixed box length. Further equilibration was performed along 10 ns of NVT-MD runs at 410, 450, 500, 550, and 600 K and then the box volume (density) was relaxed by 1 ns of NPT-MD runs. After the equilibration period, NPT MD simulations were extended to 600 ns (0.6 μs) at each temperature mentioned above using either PYS or B

DOI: 10.1021/acs.macromol.5b00823 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules TraPPe-UA FF. The pressure in all NPT simulations was set up to 1 atm. 2.3. Cooling-Rate Simulations. The glass transition and crystallization processes of polymer chains can be studied by cooling an amorphous system at a constant rate. If the cooling rate is sufficiently low, the early stages of nucleation and crystallization can be studied.51,52 On the contrary, at high cooling rates a glass state is produced.53 In this report, simulations of several cooling rates from an amorphous system in the range of 2000 to 0.05 K ns−1 were performed. The cooling processes were performed from an initial temperature of 600 K to a final temperature of 50 K using a piecewise linear function, unless otherwise indicated. 2.4. Force Fields. Table 1 collects the potential parameters of the two FFs used in this study. The PYS-FF was optimized to reproduce the heat of vaporization and the pressure of short n-alkanes where CH2 and CH3 beads behave identically.1 However, the authors reported that they were unable to reproduce the chain-length dependence of the experimental PVT behavior. In spite of these deficiencies, this FF has been widely used to simulate melt properties and crystallization process of long alkanes and polyolefins.2−4 The other potential used in this study is TraPPE-UA, developed by the Siepmann group.10−12 This FF was parametrized using a broad set of fluid-phase equilibria data for n-alkanes, with a special stress on the transferability of the given interaction sites. Remarkably, this potential makes differences between CH2 and CH3 interaction sites. The TraPPE-UA FF has shown to be applicable to a number of polymers giving a good description of their behavior in the melt.13−31

FFs. From the values in Figure 1 it is possible to obtain the thermal expansion coefficient: α=

⎛ δ[lnV ] ⎞ 1 ⎛⎜ δV ⎞⎟ =⎜ ⎟ ⎝ ⎠ V δT P ⎝ δT ⎠ P

(1)

The results are listed in Table 2 in comparison with those found in the literature both from experiments and computer Table 2. Volumetric Thermal Expansion Coefficients Obtained from MD Simulations Using the Different Force Fields Compared with Those Obtained from the Literature force field a

PYS TraPPe-UAa MD PYS C150b MC C172c MC C224c exptl C142d exptl HDPEe exptl HDPEf exptl amorphous PEg

α × 104 K‑1 h 5.95 8.57 5.00 7.18 7.09 8.03 7.70 7.92 8.80

(−25.9%) (+6.7%) (−37.7%) (−10.6%) (−11.7%)

a

This work. bC150 (MD simulations by Yi et al. in ref 7. cC172 and C224 (MC simulations by Foteinopolou et al. in ref 23. dC142 (experimental value in ref 50). eHDPE (experimental value in ref 55). fHDPE (experimental value in ref 56). gAmorphous PE in ref 57. hRelative error to ref d, in brackets.

3. RESULTS AND DISCUSSION 3.1. Physical Properties, Conformational Features, and Melt Dynamics. As a check of the performance of the two FFs used in this work, equilibrium physical properties, conformational features and dynamics of the polymer melt are compared with experimental results. Density and Thermal Expansion Coefficient. Density values were obtained from long NPT simulations at different temperatures. Figure 1 shows the values of the density, ρ,

simulations. In general, the simulations using PYS FF and other FFs from the literature give rise to an underestimation of the temperature dependence of density values. Special good agreement is found for the thermal expansion coefficient α obtained with the TraPPe-UA FF, showing a difference with respect to the experimental data of around 7%. Coil Dimensions and Their Temperature Coefficients. Table 3 presents the results of the mean-squared end-to-end distance and the radius of gyration to molecular weight ratios, ⟨R02⟩/M and ⟨Rg2⟩/M, of the equilibrated systems at different temperatures obtained from the two FFs. The ratio between ⟨R02⟩ and ⟨Rg2⟩ is 6.35 ± 0.04 for TraPPe-UA and 6.32 ± 0.05 for PYS. The agreement with Flory’s random coil hypothesis is then noticeable in both cases, a result that has been obtained earlier for well-equilibrated PE melts.58−60 It is noteworthy the slightly larger characteristic ratios (Cn = ⟨R02⟩/nl2) obtained in the case of TraPPe-UA when compared to those obtained from the PYS case. The higher computed values obtained with TraPPE-UA are in agreement with those obtained from the end-bridging MC method using the FF developed by Mavrantzas and Theodorou (C∞ ∼ 9.13)61 and with the values recently obtained using the TraPPe-UA for PE (C∞ ∼ 8.1− 8.5).22,23,26−30 These values deserve a brief historical revision. Indeed, the characteristic ratio of linear PE has been deeply discussed along decades. One can find values of C∞ ∼ 5.7−6.4 obtained from indirect measurements of the radius of gyration by means of solution viscosity measurements under Θconditions.62 Even higher values of C∞ ∼ 8.7−10.5 were obtained from light scattering measurements in Θ-conditions.63 Boothroyd et al. reported experimental data for ⟨Rg2⟩ obtained by small angle neutron scattering (SANS) of two molten PE samples in the temperature range of 380 to 480 K.64 At 450 K, they obtained values for C∞ of 6.7 and 8.1 for both samples. Our computed values obtained at 410 K (Cn = 9.3 and 8.4, for TraPPe-UA and PYS, respectively) are still larger than those

Figure 1. Comparison of density−temperature curves of C192 n-alkane using the two FFs and experimental data from ref.50 Error bars indicate standard error of the mean from independent runs. The standard error of the values obtained is less than 1%.

obtained for equilibrated melts of C192 using the two FFs. The simulated results are compared with PVT experiments obtained for C142 and PE (Mn = 28.0 kg·mol−1) in the same temperature range.50 As it can be clearly seen, the simulations performed with the Trappe-UA FF match with the experimental results. However, the simulations performed with PYS FF overestimate the values of ρ around a 10%. In addition, the two FFs give a different temperature dependence of ρ. We have obtained the volumetric thermal expansion coefficients from the MD simulations for the two C

DOI: 10.1021/acs.macromol.5b00823 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Table 3. Physical Properties and Coil Dimensions Obtained from the MD Simulations Using the Two Force Fields Explored in This Work at Different Temperatures T, K

ρ, g·cm‑3

⟨R02⟩/M, nm2·mol·kg‑1

410 450 500 550 600

0.837 0.818 0.795 0.771 0.748

14.20 12.93 12.85 12.61 11.95

410 450 500 550 600

0.779 0.754 0.723 0.691 0.658

15.70 15.00 13.88 13.34 12.70

Cna

⟨Rg2⟩/M, nm2·mol·kg‑1

κ × 103, K‑1

2.22 2.06 2.04 1.99 1.90

−0.74

2.46 2.34 2.20 2.10 2.01

−1.10

PYS

a

8.40 7.66 7.61 7.47 7.08 TRAPPE-UA 9.27 8.89 8.18 7.89 7.51

Cn(R0) = ⟨R02⟩m0/Ml02

Figure 2. (a) Intermolecular pair distribution function, ginter(r), as a function of the radial distance, r, for the C192 system at T = 450 K for the two FFs studied. (b) Same as in Figure 2a but for the intramolecular pair density function, gintra(r).

obtained in molten PE using SANS by Horton et al.65 and Fetters et al.66 (C∞ = 7.4 ± 0.8) at 413 K. The computed values at 500 K (Cn = 8.2 and 7.6 for TraPPe-UA and PYS, respectively) are still slightly larger than the more recent results obtained by Zamponi et al. (Cn = 7.2 ± 0.2) at 509 K.44 It is worth mentioning here that the SANS experiments were performed on hydrogenated polybutadiene samples. The polybutadiene parent polymer is obtained by anionic polymerization, and contains a minimum amount of 1.9 vinyl groups per 100 carbon atoms. Then, the equivalent hydrogenated polymer contains 1.9 ethyl branches every 100 carbon atoms.67 It has been reported that this amount of short chain branches along the backbones can affect the global size of the chains.26−30 The slightly higher values of the characteristic ratio obtained in TraPPe-UA model, compared to the PYS case, can be attributed to the different torsional potential employed.13−17 However, a calculation of the backbone dihedral angles shows that the probability of encountering a dihedral angle in trans conformation in the sample is quite similar in TraPPe-UA and PYS models in the temperature range explored (for instance, 66.4% and 65.3%, respectively at 450 K, an intermediate temperature within the experimental range of 413−509 K). Then, the only difference in backbone dihedral angles would not explain the differences found between TraPPe-UA and PYS FFs. Let us recall that the local chain packing is generally used to characterize the organization of monomers in the melt. In Figure 2a the intermolecular pair distribution functions

obtained from the two FFs at 450 K are plotted. The typical broad peaks at around 5, 10, and 15 Å attributed to intermolecular packing are shown, being the pattern similar in the two FFs. However, the height and position of the nearest-neighbor peak slightly depend on the FF used, indicating a stronger packing of chains for PYS model than for TraPPe-UA. This obviously accounts for the higher values of macroscopic density obtained using the PYS-FF as discussed above. We have also detected a clear difference between the two FFs concerning the intramolecular radial distribution function, gintra(r), i.e., the bonded restraints and packing among segments of the same chain in the melt. As it can be observed in Figure 2b the intramolecular distances at 450 K for the second (1−3 distance), third (1−4 gauche distance) and fourth (1−4 trans distance) peaks are higher for the TraPPe-UA FF as compared to PYS. Consequently, the chain dimensions obtained for the former are higher. Another interesting feature of coil dimensions is their temperature dependence, given by κ = d[ln⟨Rg2⟩]/dT. This property has not been usually tested by means of computer MD simulations. The role of κ on the packing length and entanglement features, such as the plateau modulus or the entanglement length, is significant. This is particularly important for cases where the melt rheology is studied over a wide temperature range. Experimental κ values have been obtained by means of θ-condition solution measurements62,68,69 and melt state experiments, including thermoelastic measurements on networks70,71 and SANS.64 It has been reported that κ varies for PE between −1.06 and −1.25 × 10−3 D

DOI: 10.1021/acs.macromol.5b00823 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 3. (a) Natural logarithm of the mean-squared radius of gyration, ln⟨Rg2⟩, as a function of temperature for the C192 system studied with the different FFs: (blue ●) PYS and (■) TraPPe-UA. The solid lines represent the best linear fits to the data. Error bars mean the standard error of the mean. The inlet shows the temperature dependence of the trans population. (b) Mean-squared radius of gyration, ⟨Rg2⟩, as a function of the trans population fraction. The trans (φ = 180°) population has been calculated by the integration of the normalized torsion distribution in the range of 120−240°.

K−1. We have extracted the values of κ from the slope of the curve ln⟨Rg2⟩ versus T, which proves to a good approximation the linearity over the temperature range explored (see Figure 3a). A very good agreement is found between the simulated and the experimental values in the case of TraPPe-UA FF. Fotainepolou et al. studied this specific behavior of κ in PE chains using the TraPPe-UA FF and Monte Carlo (MC) simulations, obtaining an identical result.22,23 The values given by the simulations performed by PYS are quite different, similarly to those found by Mondello et al. for short nalkanes.72,73 It has been suggested that an increase of the trans population of dihedral angles is the main factor for the considerable expansion of chain dimensions of PE with decreasing temperature.22,23 In the inlet to Figure 3a, we have compared the temperature dependence of the trans population obtained from the simulations. The differences in population are very small, although the trends of the curves are parallel to those observed in the temperature coefficient, κ. This result indicates that not only the dihedral angle population plays a role in the characteristic coil sizes of the systems. In fact, there is not a unique correlation between ⟨Rg2⟩ and the trans population for both FFs, as it can be easily demonstrated by plotting both magnitudes (see Figure 3b). This result demonstrates that not only the valence terms of the potential have an effect on the dimensions of polymeric chains, and points toward the importance of nonbonded chain interactions in these features. Static Structure Factor of the Melt. In order to analyze structural properties of the chains, the static structure factor has been calculated by means of the expression:

Figure 4. Static structure factor, S(q) − 1, as a function of wavenumber, q, of the C192 system as obtained from the simulations with the two force fields. The symbols are experimental data for PE obtained from the literature at 430 K.74

experiments, as it is able to describe with high accuracy both the positions and the intensity of the characteristic peaks. The differences found in the static structure factor obtained from the different simulations concerning the locations and intensity on the peaks, are in agreement with those found also in density and intermolecular distances discussed above. Melt Dynamics and Dynamic Structure Factor. The inner segments mean-squared, g1(t), and the center of mass, g3(t), displacements can be directly obtained from the MD simulations as g1(t ) =

N

S(q ⃗ ) =

1 ⎯r )|2 ⟩ ⟨|∑ exp(iq ⃗ ·→ n N n=1

1 2k + 1

N /2 + k

∑

⎯r (t ) − → ⎯r (0)]2 ⟩, ⟨[→ n n

n = N /2 − k

(3)

where k = 20

(2)

In Figure 4 we show S(q) − 1 at 450 K obtained from TraPPe-UA and PYS models, together with the experimental pattern obtained from the literature for a PE sample at 430 K.74 Both FFs reproduce the typical features of an amorphous PE. The nearest neighbor shell causes an “amorphous halo” around the maximum qmax of S(q) − 1. Oscillations around zero appear as q increases. An excellent agreement between the simulated and the experimental patterns is observed both at low (corresponding to the intermolecular correlations) and high (representing the intramolecular correlations) values of q. The results given by TraPPe-UA FF fits very well with the

g 3(t ) =

1 N

N

∑ ⟨[ rcm⃗ ,n(t ) − rcm⃗ ,n(0)]2 ⟩ n=1

(4)

From the results obtained for g1(t) and g3(t) we have tested the time−temperature superposition principle in segmental motion within the temperature range covered by the simulations.75 First, we have applied a vertical factor, bT, to account for the chain dimensions temperature difference, by using the corresponding simulated κ values (see Table 3). Second, the curves at the different temperatures have been horizontally shifted to a reference temperature (in this case we E

DOI: 10.1021/acs.macromol.5b00823 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules have selected the common temperature of TR = 500 K). The results obtained in the superposition for the two FFs can be observed in the Supporting Information, Figure S.1. An Arrhenius temperature dependence of the horizontal shift factor has been obtained for the two cases. From the slope, it is possible to extract an activation energy for the diffusion process for the two FFs, which nicely agrees with the experimental value of linear n-alkanes with similar molecular weight and linear PE (Ea = 24.2−24.7 kJ·mol−1).76−80 The shapes of the g1(t) and g3 (t) curves are very similar in the two cases but they are shifted each other due to (i) the slightly different conformation of the chains (vertical shift), and (ii) the different chain dynamics (horizontal shift). It is possible to normalize the dynamics using characteristic properties of the chain; for example, the mean squared radius of gyration, ⟨Rg2⟩, at the corresponding common temperature of TR = 500 K, and the chain Rouse time, τR, defined as g3(τR) = ⟨Rg2⟩. The values of both ⟨Rg2⟩ and τR used for the normalization are listed in Tables 3 and 4 for each model and the results can be observed in Figure 5.

the reptation theory. A deviation from the Rouse dynamics (t1/2) within the intermediate time scale is obtained (dashed line in Figure 5). The corresponding exponents, quantifying the scaling of g1(t) with t, are close to 1/2 and 1/4. The upturn takes place at the characteristic entanglement time, τe, which should be independent of molecular length. The value obtained is τe = 2.3 ns at 450 K using TraPPe-UA FF, which is in perfect agreement with our previous works for higher molecular weight26−30 and recent results from the literature [(τe = 1.6− 2.9 ns at 450 K].18−21,81 A higher value is obtained in the case of PYS, accounting for a slower segmental melt dynamics. A second upturn in g1(t) to a t1/2 dependence takes place at higher times, associated with the whole chain Rouse regime at τR = τeZ2, where Z is the number of entanglements, which in our case is around 3. To obtain the number of entanglements in our system it is necessary to estimate the step length of the primitive path. For this purpose, we have adopted the procedure described by Likhtman and McLeish,82 widely used in the literature.26−30 The reader is referred to these works for detailed information. Finally, in the original tube theory at a characteristic time τd = 3τeZ3 = 3τRZ the chain completely leaves its original tube and the dynamics are determined by diffusion, with a time dependence t1. The values of all characteristic times obtained are listed in Table 4. Regarding the behavior of g3(t) in the Rouse region, it can be observed (see Figure 5) that the power law tz is in disagreement with the ideal exponent predicted by the Rouse model for free diffusion (z = 1). Our calculations show a characteristic time exponent of z = 0.77 (irrespective of the FF used), in agreement with the experimental results obtained by Zamponi et al.44 This subdiffusive displacement of the center of mass has already been recognized by computer simulations in both unentangled and entangled polymeric liquids.83−89 Additional intrachain and intermolecular effects have been proposed as the origin for the Rouse model breakdown of the time evolution of g3(t). It has been shown that long-range interactions on length scales shorter than Rg lead to this failure. This is due to segmental correlation hole effect, which induces repulsive interactions when bringing chain segments together and causes density fluctuations in the systems.90−92 A few theoretical approaches, based on generalized Langevin equations, have been developed in order to explain the observed anomalies within the correlation hole effect framework. Guenza proposed that the inclusion of the local intermolecular forces improves the description of the global dynamics in the short-time domain, and causes a nonuniversal scaling with a power law exponent z decreasing as the chain length increases.93−95 A recent approach states that the correlation hole effect alone is insufficient to explain the observed subdiffusive behavior. Farago et al.96 have considered the transient collective flow effect of viscoelastic hydrodynamic interactions, which are not screened even beyond the static monomer correlation length. Figure 5 demonstrates that the two FFs effectively describe the melt dynamics of C192, which shows the typical features of a slightly entangled system. Unfortunately, there are no experimental data to directly compare with g1(t). Alternatively, a comparison can be made using results from the normalized single-chain intermediate dynamic structure factor recently obtained from neutron spin echo (NSE) at 509 K for a set of monodisperse PEs.44 These experimental results include a sample with a molecular weight similar to the one used in our simulations (Mw = 2.80 kg·mol−1). The polymers used in that study were obtained by anionic polymerization of butadiene

Table 4. Dynamic Features Obtained from the MD Simulations

a

T, K

τe, ns

τR, ns

410 450 500 550 600

7.3 3.7 2.1 1.2 0.73

58.9 29.5 16.6 9.5 5.9

410 450 500 550 600

4.2 2.3 1.2 0.66 0.43

34.1 18.8 9.6 5.3 3.5

τd, ns

D ×107, cm2 s‑1

PYS 1350 675.7 380.2 217.6 135.1 TraPPe-UA 781.1 430.6 219.9 121.4 80.2

Ea, kJ·mol‑1

1.7 2.8 5.3 8.8 13.4

24.2 (22.3)a

2.8 4.5 8.6 15.6 22.9

24.7 (23.2)a

Obtained from the self-diffusion coefficient temperature dependence.

Figure 5. Normalized mean−squared displacements of the center of mass, g3(t) (solid symbols) and of the inner chain segments, g1(t) (open symbols) for the (red ●) PYS and (blue ■) TraPPe-UA force fields. Dashed line indicates the expected g1(t) for a nonentangled Rouse chain. The vertical dotted lines represent the reduced characteristic times (τe/τR and τd/τR) predicted by the reptation model (see text).

From this figure, we can conclude that the two FFs essentially capture the same dynamics of the polymeric system. The characteristic signature of an entangled polymeric system is observed, given by the four characteristic regimes predicted by F

DOI: 10.1021/acs.macromol.5b00823 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 6. Comparison of experimental44 (symbols) and MD results (lines) of the dynamic structure factor for C192 obtained from simulations with (a) PYS and (b) TraPPE-UA for q = 1.15 nm−1 (squares), 0.96 nm−1 (circles), and 0.77 nm−1 (triangles).

with subsequent hydrogenation (and deuteration), and are known to contain at least ∼2.0% of ethyl branches. Although we have recently probed by combination of simulations and experiments that such a low amount of branches causes differences in the melt dynamics of PE,26−30 we have tried a direct comparison between our simulations and the experiments performed by Zamponi et al. The normalized dynamic structure factor can be obtained directly from simulations using the following expression: N ⎯r (t ) − → ⎯r (0))]⟩ ∑n , m = 1 ⟨exp[iq ⃗ · (→ S(q ⃗ , t ) n m = N → ⎯ → ⎯ S(q ⃗ , 0) ∑n , m = 1 ⟨exp[iq ⃗ · ( rn(0) − rm(0))]⟩ (5)

From experiments and simulations performed at sufficiently low momentum transfers q, the structure factor is determined by the diffusion of the center of mass of the chains, ⟨ΔRcm2(t)⟩, with internal modes being negligible:

Figure 6 illustrates the intermediate dynamic structure factor for various values of the scattering vector q (0.77, 0.96, and 1.15 nm−1) as obtained directly from the MD simulations for the C192 (2.69 kg·mol −1) system in comparison with the experimental results for a monodisperse PE sample (Mn = 2.80 kg·mol−1) at 509 K.44 It is clear from this comparison that TraPPe-UA FF (Figure 7b) is in good agreement with the experimental data (Figure 6a).

The values of D obtained in the two simulated systems are listed in Table 4. Our calculation with TraPPE-UA FF shows that D = 4.5 × 10−7 cm2·s−1 for C192 at 450 K. This calculated value agrees quite well with the experimental data taken from different sources at the same temperature.76−80 PYS-FF gives rise to a lower value of D, as expected due to the slower dynamics. By the other hand, the values of D obtained by the simulations at different temperatures are usually fitted to the Arrhenius type equation ln D = ln D∞ − Ea,D/T, where Ea,D is the activation energy for diffusion. We have also fitted our results of D to this equation in order to obtain Ea,D values (listed in Table 4 in brackets). The results are in agreement with those obtained previously from the more rigorous time− temperature superposition approach, as in this latter case, the effect of temperature in chain dimension (i.e., the vertical shift) is taken into account. It is evident from the results shown and discussed in this section that the comparison of simulations and experimental data including density, molecular size, static properties, melt dynamics and temperature dependences obtained for C192 have shown an excellent agreement in the case of MD simulations performed with the TraPPe-UA FF. Previous studies have already reported results in this direction, but generally focusing their attention in static features of both unentangled and entangled linear PE, or in dynamic properties of shorter PE chains.2−4,13−17 Then, the objective of the next section is to show the ability of the FFs studied to simulate more complex phenomena as glass transition and crystallization. 3.2. The Glass Transition Temperature. Upon cooling, a polymer freezes below the glass transition temperature, Tg, in a glassy state. In the temperature interval of cooling the polymer melt evolves from the liquid to the solid state, and the transition can be defined as the intersection of straight line extrapolations from the glassy and liquid states of the volume-

S(q , t ) ⎡ q ⎤ 2 = exp⎢ − ⟨ΔRcm (t )⟩⎥ ⎣ ⎦ S(q , 0) 6

(6)

Then an estimation of the center of mass mean-square displacement is possible, which in the infinite time limit is directly related to the self-diffusion coefficient, D: D = lim

t →∞

Figure 7. Apparent glass transition temperature versus cooling rate for the atomistic model C192 obtained from MD simulations using (blue ●) PYS and (■) TraPPe-UA force fields. The inlet shows the master curve obtained by only shifting the data horizontally using the extrapolated glass transition, Tg0 (see text). The lines represent the fit to eq 8.

We can anticipate a better description of the experiments using the TraPPe-UA FF by considering the inclusion of a minimum amount of short chain branches (∼2.0% of ethyl branches) in our C192 molecular model, as expected in a hydrogenated polybutadiene (results not published yet). G

2 ⟨ΔRcm (t )⟩ 6t

(7)

DOI: 10.1021/acs.macromol.5b00823 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

interactions. All of these features are actually depicted by the specific FF used in the simulation. Indeed, as both flexibility and free volume increase (higher degree of conformational freedom) and chain−chain interaction decreases (less efficient chain packing), Tg decreases. The difference in the population of gauche and trans conformational states for these FFs are not very high, as seen in preceding section. However, a substantial difference in the intramolecular distances between the TraPPeUA and PYS FFs (see Figure 2b) could give rise to an extra free volume in the former, which might account for the lower value of Tg0. 3.3. The Crystallization Process. To study the crystallization process of the entangled C192 system, we have performed MD simulations at very low Γ values. In Figure 8,

temperature curve. The equilibrated C192 melt system was cooled down at different rates, Γ = 1 to 2000 K·ns−1 (high enough to prevent polymer crystallization). Then, the specific volume as a function of temperature was calculated for the two FFs. A change of the slope of the curves indicates the apparent Tg, which is known to be dependent on the cooling rate.97 Moreover, the high cooling rates used in the simulations generally lead to a higher estimate of Tg as compared to experiments. To account for this effect, at least 10 different cooling rates were used within the range indicated above. The minimum cooling rate in our simulations is 1 order of magnitude lower than the one used in the most recent simulations of PE.5−7 The plot of the apparent Tg versus cooling rate is shown in Figure 7. The results in this figure reveal that Tg decreases nonlinearly with the logarithm of Γ. This behavior has been already reported both experimentally97 and by simulations.98−101 In general, it has been shown that the dependency is a nonlinear equation, adopting the form: Tg = T g0 +

A ln(BΓ)

(8)

This expression can be extracted from the assumption of a Vogel−Fulcher dependence of the relaxation times. The application of eq 8 to calculated data provides the value of Tg0. Similar extrapolation procedures have been also applied in other computer simulations.98−101 The resulting Tg0 values show similar dependency with the experimental data, in spite of the big difference in cooling rate employed, typically 107 - 1010 K·min−1 in simulations and 10 K·min−1 in experiments. The application of eq 8 to the results in Figure 7 gives estimated values of Tg0 = 187.0 and 214.1 K for TraPPe-UA and PYS FFs, respectively. Remarkably, the Tg0 value obtained for PYS is very close to that obtained using a different extrapolation method with the same FF by Yi et al. (215.22 and 223.04 K for both C150 and C1000 systems, respectively).7 The values of A and B in eq 8 are nearly the same in the two cases, A = 368.7 K and B = 0.023 ns·K−1. In fact it is possible to obtain a master curve from the two simulation sets by simply plotting Tg − Tg0 versus Γ, as it can be observed in the inlet to Figure 7. It is well-known that the literature is plenty of contradictory experimental Tg values in the case of linear PE. Illers underlined that an exact identification of Tg of linear amorphous PE from experiments is quite difficult due to the characteristic high rate of crystallization which prevents any freezing of the pure amorphous state.102 Because of this fact, it has been usual to obtain the Tg of PE by extrapolation. The extrapolation of Tg values in ethylene copolymers using the Gordon−Taylor equation gives rise to a value of 185−195 K. This is in contrast with other values of Tg for amorphous PE reported in the literature. In one case the Tg has been identified with the γ transition (120−150 K), whereas in other case is identified as the β transition of semicrystalline linear PE (240−260 K).103−105 The extrapolated value obtained from the TraPPeUA system matches the expected value for an amorphous PE free of constraints provoked by the presence of crystals, i.e. that obtained by Gordon−Taylor equation. Previous estimates of Tg in PE by computer simulations resulted in values located between 200 and 285 K. This disparity of values could be attributed to the different FF and/or to the extrapolation methods used.106−108 It should be noted at this point that Tg in polymers depends on factors as chain stiffness, chain bulkiness (free volume) and the strength of polymer chain−chain

Figure 8. Specific volume (v) and specific heat per mol of CH2 (Cp) versus temperature (T) at different cooling rates: (a) PYS and (b) TraPPe-UA. The crystallization temperature peaks in the Cp−T curve during cooling are marked by arrows. The stripped square indicates the temperature range in which the crystallization has been observed at experimental conditions.

we display the evolution of specific volume and specific heat at different Γ from 0.05 to 1 K·ns−1. A sudden drop of specific volume and a single peak of specific heat clearly indicate a phase transition at a crystallization temperature, Tc, and a semicrystalline polymer is formed when temperature is below the corresponding Tc (indicated by the arrows). A continuous shift of the Tc peak toward lower temperatures and a broadening of the transition with increasing Γ are clearly observed. Both PYS (Tc ∼ 300 K) and TraPPe-UA (Tc ∼ 325 K) FFs give values of Tc close to those experimentally obtained by ultrafast calorimetry (Tc ∼ 370 K) for long n-alkanes (C122, C162, C390) at Γ of 103−104 K·s−1. Please note that these experimental Γ values are 4 orders of magnitude lower than those used in the simulations (see Table 5).109 It is also very interesting the similarity obtained between our simulations and the recent crystallization experiments in linear PE using ultrafast calorimetry performed by Toda et al.110 Both H

DOI: 10.1021/acs.macromol.5b00823 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules Table 5. Degree of Crystallinity, α, Crystallization Temperature Peak, Tc, and Crystal Thickness, lg* PYS Γ, K ns 1 0.5 0.1 0.05

‑1

TraPPe-UA

α

Tc, K

lg*, nm

α

Tc, K

lg*, nm

0.06 − 0.12 0.25

247 − 281 297