Table Of ContentResolving Dynamic Properties of Polymers through Coarse-Grained Computational
Studies
K. Michael Salerno,1 Anupriya Agrawal,2,3 Dvora Perahia,3 and Gary S. Grest1
1Sandia National Laboratories, Albuquerque, NM, 87185
2Department of Mechanical Engineering and Materials Science, Washington University, St. Louis, MO 63130
3Department of Chemistry, Clemson University, Clemson, SC 29634
Coupled length and time scales determine the dynamic behavior of polymers and underlie their
unique viscoelastic properties. To resolve the long-time dynamics it is imperative to determine
6 which time and length scales must be correctly modeled. Here we probe the degree of coarse
1 graining required to simultaneously retain significant atomistic details and access large length and
0 time scales. The degree of coarse graining in turn sets the minimum length scale instrumental
2 in defining polymer properties and dynamics. Using linear polyethylene as a model system, we
probe how coarse graining scale affects the measured dynamics. Iterative Boltzmann inversion is
n
used to derive coarse-grained potentials with 2-6 methylenegroups per coarse-grained bead from a
a
fully atomistic melt simulation. We show that atomistic detail is critical to capturing large scale
J
dynamics. Using these models we simulate polyethylene melts for times over 500 µs to study the
3
viscoelastic properties of well-entangled polymer chains.
1
]
t Polymer properties depend on a wide range of cou-
f
o pled length and time scales, with unique viscoelastic
s properties stemming from interactions at the atomistic
.
t level. The need to probe polymers across time and
a !"%$ !"#$
m lengthscalestocapturepolymerbehaviormakesprobing
dynamics, and particularly computational modeling, in-
-
d herently challenging. With increasing molecular weight,
n polymer melts become highly entangled and the long-
o timediffusiveregimebecomescomputationallyinaccessi- !"’$ !"&$
c
bleusingatomisticsimulations. Inthesesystemsthedif-
[
fusivetimescaleincreaseswithpolymerizationnumberN
1 faster than N3, becoming greater than 1010 times larger
v
9 than the shortest time scales even for modest molecular FIG.1. Asingle C96H194 PE chain representedwith increas-
weightpolymers. Whileitisclearthatthelargestlengths ing degree of coarse graining λ = 2, 3, 4, and 6 methylene
8
2 scales of polymer dynamics are controlled by entangle- groupsperCGbead. Beaddiametercorrespondstothemin-
3 ments, the shortest time and length scales required to imum in thenonbonded interaction for each CG model.
0 resolve dynamic properties are not obvious. This knowl-
.
1 edge is critical for developing models that can trans-
grainedwhilestillappropriatelycapturingpolymerprop-
0 pose atomistic details into the long time scales needed
erties and dynamics [5]. The current study probes the
6
to model long, entangled polymer chains.
1 effects of the degree of coarse graining of polymers on
: One path to overcoming this computational challenge their dynamic and static properties.
v
is to coarse grain the polymer, reducing the number of With the vast efforts to coarse grain polymers, many
i
X degrees of freedom and increasing the fundamental time models have emerged with differences in the number of
r scale. The effectiveness of this process depends on re- atoms combined in each bead and the procedure for de-
a
taining the smallest length scale essential to capturing termining the interaction potentials. One of the most
the polymer dynamics. The process of coarse graining commoncoarse-grained(CG) models for polymers is the
amounts to combining groups of atoms into pseudoatom united atom (UA) model, which combines each CH
n
beads and determining the bead interaction potentials group into a pseudoatom. The UA interaction param-
[1, 2]. Simple models like the bead-spring model [3], eters are determined phenomenologically to reproduce
capture characteristicsdescribedby scalingtheories,but physicalproperties suchasdensities and criticaltemper-
disregardatomisticdetailsandcannotquantitativelyde- atures[6–9]. AnothermodelcommonlyusedistheMAR-
scribe properties like structure, local dynamics or den- TINImodel,whichutilizes thesameapproach,matching
sities. Immense efforts have been made to systemati- bulkdensitiesandcompressiblitiesofshortalkanechains
cally coarse grain polymers and bridge the gap of time at a larger scale of four CH groups per CG bead [10].
2
and length scales while retaining atomistic characteris- Moreadvancedmethodssuchasforcematching,iterative
tics [4]. One critical issue underlying the coarsegraining Boltzmanninversion,andoptimizedrelativeentropy[11–
process is the degree to which a polymer can be coarse 13]haverecentlybeendevelopedtoincorporateatomistic
2
detailintotheCGmodel. Withthesemethodsthereisan
openquestionas to the number ofatoms to representby
1
a single bead and the effect of this coarse-graining scale ol] CG2
m CG3
on the measured properties of the system. One critical al/ CG4
physicsquestionremainsunresolved: namelydefiningthe kc 0.5 CG4-M
shortestlengthscaleinapolymerthatisfundamentalto y [ CG6
g UA
the macroscopic dynamics and properties.[1, 3, 4, 6–13] er 0
n
Here, this issue is addressed through the development E
d
of CG models with increasing degree of coarse graining e
using iterative Boltzmann inversion. By examining how nd -0.5
o
well these CG models describe both the static and dy- nb
o
namic properties of a polymer melt, using polyethylene N -1
(PE)asamodelsystem,weprobethisoutstandingques- 0.4 0.5 0.6 0.7 0.8 0.9 1
tion. The backbone of PE consists of -CH2- methylene Radius [nm]
groups that provide a natural unit or scale for coarse
graining. ThoughthechemicalstructureofPEissimple,
it is a thermoplastic materialuseful in a largenumber of FIG.2. TabulatedpairpotentialbetweenCGbeads. TheUA
and CG4-M potentials are included for comparison.
applications, with tunable mechanical properties deter-
mined by the degree of branching.
Polyethylenechainshavepreviouslybeenstudiedusing
CG models with beads of λ = 3−48 methylene groups potential with modified dihedral coefficients that better
per bead [14–19]. These studies were able to capture reproduce the properties of long alkanes [24]. With this
the radius of gyration as a function of molecular weight modified potential the mean squared radius of gyration
and the pair correlationfunction between CG beads. As hRg2i and end to end distance hR2i match experimental
mostofthese studies useda largedegreeofcoarsegrain- values [25, 26] better than with original OPLS param-
ing(λ∼20)tostudydynamicalproperties,anextracon- eters [27]. For the CGλ models, hR2i for n=96 chains
straintwasneededtopreventchainscuttingthrougheach is within 20% of the atomistic value, while hR2i for the
other [20]. With this extra constraint,the mean squared MARTINI model is 50% too high. Static properties for
displacement(MSD),stressautocorrelationfunctionand different chain lengths are reported in the Supplement.
shear viscosity of linear and branched PE [20–22] have Tabulated CG angle and bond potentials were deter-
been studied for long, entangled chains. However, these minedbyBoltzmanninversionoftheatomisticbondand
studies did not account for or study the effects of the angle distributions in Fig. S1. Torsion terms were omit-
coarse-grainingdegree λ on dynamic properties. ted in all CG models, which may accountfor the shorter
Hereforthefirsttime,weelucidatetheeffectofcoarse- endtoenddistanceslistedinTableSIfortheCG2model.
graining degree on the ability to capture both the struc- Tabulatednonbondedpotentialsweredeterminedbyiter-
tureanddynamicsofPE.Weareabletocapturepolymer ative Boltzmann inversion[4]. The intermolecular radial
chain dynamics for lengths up to C H and time distribution function g(r) fromthe atomistic simulation,
1920 3842
scales of 400 µs using models that accurately represent shown in Fig. S2, was used as the target for iteration of
atomistic detail. Accessing large length and time scales the nonbonded potentials shown in Fig. 2. Also shown
allowsus to measurequantities like the plateaumodulus are 6-12 Lennard Jones pair potentials for the united
which depend on a hierarchy of length and time scales. atom (UA) model of Yoon et al. [6], and the MARTINI
Coarse-grained beads shown in Fig. 1 represent λ (CG4-M)model[10]. TheMARTINIparameterǫwasre-
methylenegroups. Westudyλ = 2,3,4and6andrefer duced from 0.8365 kcal/mol to 0.803 kcal/mol to match
to these models as CGλ. We find that for surprisingly the density ρ = 0.72 g/cm3 of atomistic simulations for
small λ the chains cross and diffuse rapidly, indicating n=96 chains. For each CG model a pressure correction
that CG features directly link to macroscopic polymer isappliedtomatchthedensityρ=0.72g/cm3 forn=96.
motion. With this result, we further probed the CG6 The similarity in length and energy scales between the
model polymer including non-crossing constraints, and CG4 and CG4-M models is evident in Fig. 2. For each
comparing with models with unconstrained dynamics. CGλ model all beads identical interactions,howeverend
The tabulated CG PE potentials were derived from beads have an extra hydrogen atom mass.
a single fully-atomistic simulation of a melt of C H The CG6 model has a surprisingly large equilibrium
96 194
PE chains at 500K. The study was then generalized to bond distance relative to the bead diameter. Therefore
melts of C H with n=96, 480 for the fully atom- a modified soft segmental repulsive potential [28] was
n 2n+2
istic model, andn=96,480,960and1920at500K using added between CG beads to inhibit chain crossing. We
the CG models. Atomistic simulations used a version of used a segmental bead diameter of 0.5 nm. This scheme
the OptimizedPotentialsforLiquidSimulations(OPLS) increases the pressure in our samples by about 80 atm
3
CG2 104
2m]110034 CGCCC4UGGG-MA346 !"’(& 2m] 103 )*& +,#-.& 0246810 time scale factor! !"’(& 1/& !"’!("%#&$%&
2|r(t)-r(0)| [n!"ii102 atomistic !"#$%& 2(t)-r(0)| [n"mcm 110012 1 2 CH32 / B4ead5 6 !"/’0%&
101 |r!c 100
100100 101 102 103 104 10-1100 101 102 103 104 105 106
Time [ns]
Time [ns] (scaled)
104 !"’(& FIG. 4. MSD of center of mass scaled with the same α scale
+)& factorasthefortheinner-CH2-groupsinFig.3(b). Thesolid
!"#$%& line has slope t1. Inset: The alpha scale factor for different
2m]103 !"’(%& coarse-grained models and theUA and Martini models.
n
2 ["
|r(t)-r(0)|!ii102 !")’*%& idceatnetrlmyifnaesttehretdhyannaimnicatsocmaliisntgicfasicmtourlaotfiotnhse[C29G–3m6]o.dTelos
101 we compare the MSD of the inner 24 methylene groups
+),#&
(4, 6, 8 or 12 beads) for CG models and the inner 24
carbonatomsforatomisticsimulationsforn=96and480
100 as shownin Fig. 3 (a). The mobility of the chains in the
100 101 102 103 104 105 106
CG models is larger than in atomistic simulations. By
Time [ns] (scaled)
scaling the time for each of the CG models we create a
single collapsed curve for each chain length for both the
FIG. 3. (a) The MSD of the inner 24 -CH2- groups of each
atomistic and CG data as shown in Fig.3 (b). Notably,
polymer chain at 500 K. (b) Same data as in (a), scaled by
α. The solid lines represent the scaling predictions t1 for the a single scaling factor α is required to collapse atom-
diffusiveregime and t1/4 for thereptation regime. istic and CG data for each model, independent of chain
length. As seen in Fig. 3b the MSD has reached the dif-
fusiveregimewhereMSD∼t1 evenforthe longestchain
at fixed density compared to simulations with no con- length n = 1920. Over intermediate time scales, the
straint. We do not re-derive the potential including the chains show the expected t1/4 scaling predicted by rep-
softsegmentalbead,howeverFig.S3showsthatthenon- tation theory [37]. These results demonstrate that one
crossingbead induces only smallchanges in g(r) relative cancapture longtime and lengthscales with CGmodels
to a model with no constraint. By eliminating the finest while accounting for atomistic details.
degreesoffreedom,CGmodelsallowasignificantlylarger The MSD of the center of mass was then measured to
time step than atomistic models. We use a time step δt test the scaling factor α. Figure 4 shows the MSD of the
= 20 fs for the CG6, CG4 and CG4-M models, 10 fs for chain center of mass for chain lengths n=96, 480, 960,
the CG3 model and 2 fs for the CG2 model, compared and 1920. These data have been scaled by the same α
to 1 fs for the atomistic model. as the monomer MSD, producing an excellent collapse.
Coarsegrainingreducesthe number of degreesoffree- The scale factor α as a function of CG model is shown
dom in a system, creating a smoother free-energy land- inset in Fig. 4 along with α for the MARTINI [10] and
scape compared with fully-atomistic simulations. This UA [6] models. Our CG potentials have a much larger
speedup can be addressed by including frictional and time scaling factor than the MARTINI and UA models,
stochastic forces [15]. This approach allows large CG similar to the time-scaling factor found previously for
scales, however rigorously correct dissipation requires a PE for a single λ [30]. Values of α are also comparable
sophisticated generalized Langevin kernel, and simplifi- to those found previously for polystyrene, modeled at
cations are usually employed[5]. It has been shownthat a similar coarse-graining level [34]. The UA model has
forsmallλCGmodelswithoutaddedfrictionorstochas- long been considered approximate to the fully-atomistic
tic forces can be employed, but CG dynamics are signif- simulationandindeedis≈40%fasterthanfully-atomistic
4
107 0.95
CG2
CG3
CG4
CG4-M 0.9
CG6 3m]
Pa] g/c
G(t) [106 sity [ 0.85
n
e CG2
D
0.8 CG3
CG4
CG4-M
CG6
105 0.75
10-1 100 101 102 103 104 105 200 250 300 350 400
Time [ns] (scaled) Temperature [K]
FIG. 5. ModulusG(t) for each of theCG polymer models at FIG. 6. Density versus temperature for the n=96 samples
500K.Filledandopensymbolsrepresentthen=96andn=480 cooled at 1.4K/ns.
chain length, respectively. Solid lines represent the n=1920
length, while n=960 chains are omitted for clarity.
by dividing the datasetsintwo andmeasuringthe varia-
tionintheplateauvalue. Similarplateaumodulusvalues
simulations. Interestingly, the time scaling factor is not werefound by Paddingand Briels[20] for n≤1000with
monotonic in CG level, with the CG2 and CG6 models λ=20 model with a non-crossing constraint.
exhibiting the largest speedup. The potential depths in The thermal expansion coefficient is another way to
Fig. 2, relate to the value of α, as described previously assess the validity of the CG models. The linear ther-
by Depa and Maranas [30]. mal expansion coefficients for the n=96 samples in the
The polymer entanglement mass M governs many temperature range from 495K to 480K are 5.0, 3.4, 3.5,
e
properties of the polymer melt and provides information 3.3, 3.4 ×10−4T−1 respectively for the CG2, CG3, CG4,
about chain mobility within the polymer mesh. Experi- CG4-M,andCG6models,comparedwith3.1×10−4T−1
mentally,M =ρRT/G0 isdeterminedfromtheplateau for the atomistic model. Hence, all the CG models with
e N
modulus G0 of the stress relaxation function G(t) [25]. λ≥3 agree with the atomistic thermal expansion.
N
Experimental values for polyethylene are 1.6-2.5 MPa, Although one can derive CG potentials using mul-
correspondingto Me of 1300-2000g/mol[25, 26, 38, 39]. tiple temperatures [40], our CG models are developed
The relaxationmodulus ineachofourCGmodelswas in the traditional way at a single temperature. Hence
measuredfor the four different chainlengths via equilib- there is uncertainty about the validity of these models
rium stress correlations using the Green-Kubo relation away from the chosen point [41]. Shown in Fig. 6 are
G(t) = (V/k T)hσ (t)σ (0)i where σ are the off- temperature-density data for the CG models from 400K
B αβ αβ αβ
diagonal components xy, xz, and yz of stress. Figure 5 to 200K. For λ ≤ 4, the curves show a pronounced den-
shows G0 for each of the CG models for n=96 and 480 sityincreasebetween250Kand330K,correspondingtoa
N
and for λ ≥ 3 for n = 1920. The times for each model semi-crystalline state. The CG6 model does not capture
havebeenscaledbythecorrespondingvalueofα.Though crystallinity and we expect that coarsermodels will not,
it shows similar behavior, the UA model is omitted be- either. Previous studies of bead-spring polymer models
cause the zero-pressure density is higher than the other indicatedthatacommensuratebondlengthandbeaddi-
models,makingcomparisondifficult. Thecurvescollapse ameterleadstocrystallization[42]. Themeltingtemper-
fortheshort-timet−1/2 regime,withlonger,moreentan- ature for n=96 is about 400K, so crystallization occurs
gled chains forming progressively more distinct plateau atalowertemperaturethanexpected,yetobservationof
regions. The plateau modulus is measured as the value anysemi-crystallinephaseisremarkable. Thissurprising
of the relaxation modulus in the plateau region, roughly featureindicatesthatalthoughourpotentialsarederived
between 20 and 600 ns. Using the longest chain length, at500Kthey maybe usefulawayfromthis temperature.
n = 1920, the plateau modulus G0 =2.2±0.3 MPa for Here we have shown that the smallest length scale
N
CG6,2.1±0.3MPaforCG4and2.1±0.7MPaforCG3, needed in the hierarchy of length scales to correctly
allwithintheexperimentalrange. ForCG4-MG0 issig- describe macroscopic behavior and properties is rather
N
nificantlyhigherG0 =3.7±0.4MPaanddoesnotagree small, between 4 and 6 monomers. The CG4 model of-
N
with experimental values. Uncertainties are measured fersaspeedupofmorethanfourordersofmagnitudeover
5
atomistic simulations, which includes contributions from (2007).
thetimestepδt,scalefactorα,and(3λ)2reductioninthe [11] S. Izvekov and G. A. Voth,J. Phys.Chem. B 109, 2469
numberofpairwiseinteractions. ThespeedupoftheCG4 (2005).
[12] M. S.Shell, J. Chem. Phys. 129, 144108 (2008).
model is about three times faster than the CG3 model
[13] V.Ru¨hleandC.Junghans,Macromol.Theor.Simul.20,
andis comparablethe CG6model. Howeverbecause the
472 (2011).
CG4 model does not require a crossing constraint, we
[14] H. Fukunaga, J.-i. Takimoto, and M. Doi, J. Chem.
prefer this model. With this realized speedup polymers Phys. 116, 8183 (2002).
as long as n = 1920 can now be simulated for over 500 [15] J.T.PaddingandW.J.Briels,J.Chem.Phys.115,2846
µs. Reaching this time scale allows probing some of the (2001).
most unique and intrinsic properties of polymers includ- [16] H.S.Ashbaugh,H.A.Patel,S.K.Kumar, andS.Garde,
ing the plateau modulus and intermediate t1/4 scaling J. Chem. Phys. 122, 104908 (2005).
[17] X. Guerrault, B. Rousseau, and J. Farago, J. Chem.
in the mean squared displacement. From the computa-
Phys. 121, 6538 (2004).
tional viewpoint, the CG models developed significantly
[18] L.-J. Chen, H.-J. Qian, Z.-Y. Lu, Z.-S. Li, and C.-C.
reduce the resources needed to study polymers for long Sun,J. Phys.Chem. B110, 24093 (2006).
times. Our results for the plateaumodulus anddiffusion [19] D. Curco´ and C. Alema´n, Chem. Phys. Lett. 436, 189
dynamicsshowthatwithoutaddingextraconstraintsthe (2007).
CG4 model captures the atomistic detail needed for cor- [20] J.T.PaddingandW.J.Briels,J.Chem.Phys.117,925
(2002).
rect dynamics from monomer to polymer scale.
[21] J. T. Padding and W. J. Briels, J. Chem. Phys. 118,
10276 (2003).
AAandDPacknowledgefinancialsupportfromGrant
[22] L.Liu,J.T.Padding,W.K.denOtter, andW.J.Briels,
No.DE-SC007908andanallotmentoftimeontheClem- J. Chem. Phys. 138, 244912 (2013).
son University Palmetto cluster. This research used re-
[23] See Supplemental Material [url], which includes Refs.
sourcesatthe NationalEnergyResearchScientific Com- [43–45].
puting Center, which is supported by the Office of Sci- [24] S. W. I. Siu, K. Pluhackova, and R. A. B¨ockmann, J.
ence of the U.S. Department of Energy under Contract Chem. Theory Comput.8, 1459 (2012).
DE-AC02-05CH11231. This work was supported by the [25] L. J. Fetters, D. J. Lohse, S. T. Milner, and W. W.
Graessley, Macromolecules 32, 6847 (1999).
Sandia Laboratory Directed Research and Development
[26] L. J. Fetters, D. J. Lohse, and W. W. Graessley, J.
Program. Research was carried out in part, at the Cen- Polym. Sci. Pol. Phys.37, 1023 (1999).
ter for IntegratedNanotechnologies,a U.S. Department
[27] W.L.Jorgensen,D.S.Maxwell, andJ.Tirado-Rives,J.
of Energy, Office of Basic Energy Sciences user facility. Am. Chem. Soc. 118, 11225 (1996).
Sandia National Laboratories is a multi-program labo- [28] T. W. Sirk, Y. R. Slizoberg, J. K. Brennan, M. Lisal,
ratory managed and operated by Sandia Corporation, a andJ.W.Andzelm,J.Chem.Phys.136,134903(2012).
wholly owned subsidiary of Lockheed Martin Corpora- [29] P. Depa, C. Chen, and J. K. Maranas, J. Chem. Phys.
134, 014903 (2011).
tion, for the U.S. Department of Energy’s National Nu-
[30] P. K. Depa and J. K. Maranas, J. Chem. Phys. 123,
clear Security Administration under contract DE-AC04-
094901 (2005).
94AL85000. [31] I.Y.LyubimovandM.G.Guenza,J.Chem.Phys.138,
12A546 (2013).
[32] I. Y. Lyubimov, J. McCarty, A. Clark, and M. G.
Guenza, J. Chem. Phys. 132, 224903 (2010).
[33] V. Harmandaris, Korea-Aust.Rheol. J. 26, 15 (2014).
[1] S.O.Nielsen,C.F.Lopez,G.Srinivas, andM.L.Klein, [34] V. A. Harmandaris and K. Kremer, Soft Matter 5, 3920
J. Phys. Condens. Mat. 16, R481 (2004). (2009).
[2] Y.Li,B.C.Abberton,M.Kr¨oger, andW.K.Liu,Poly- [35] D. Fritz, K. Koschke, V. A. Harmandaris, N. F. A.
mers 5, 751 (2013). vanderVegt, andK.Kremer,Phys.Chem.Chem.Phys.
[3] K. Kremer and G. S. Grest, J. Chem. Phys. 92, 5057 13, 10412 (2011).
(1990). [36] B. Hess, S.Leon, N.vanderVegt, and K.Kremer, Soft
[4] F. Mu¨ller-Plathe, Chem. Phys. Chem. 3, 754 (2002). Matter 2, 409 (2006).
[5] J.T.PaddingandW.J.Briels,JournalofPhysics: Con- [37] P. G. de Gennes, J. Chem. Phys.72, 4756 (1980).
densed Matter 23, 233101 (2011). [38] J. F. Vega, S. Rastogi, G. W. M. Peters, and H. E. H.
[6] W.Paul, D.Y.Yoon, andG.D.Smith,J.Chem.Phys. Meijer, J. Rheol. 48, 663 (2004).
103, 1702 (1995). [39] V. R. Raju, G. G. Smith, G. Marin, J. R. Knox, and
[7] S.K.Nath,F.A.Escobedo, andJ.J.dePablo,J.Chem. W. W. Graessley, J. Polym. Sci. Pol. Phys. 17, 1183
Phys.108, 9905 (1998). (1979).
[8] M.G.MartinandJ.I.Siepmann,J.Phys.Chem.B102, [40] T.C.Moore,C.R.Iacovella, andC.McCabe,J.Chem.
2569 (1998). Phys. 140, 224104 (2014).
[9] M.MondelloandG.S.Grest,J.Chem.Phys.103,7156 [41] P.Carbone, H.A.K.Varzaneh,X.Chen, andF.Mller-
(1995). Plathe, J. Chem. Phys.128, 064904 (2008).
[10] S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tiele- [42] R. S. Hoy and N. C. Karayiannis, Phys. Rev. E 88,
man, and A. H. de Vries, J. Phys. Chem. B 111, 7812 012601 (2013).
6
[43] M.Tuckerman,B.J.Berne, andG.J.Martyna,J.Chem. Chem. Phys.137, 174107 (2012).
Phys.97, 1990 (1992). [45] S. Plimpton, J. Comput. Phys.117, 1 (1995).
[44] R. E. Isele-Holder, W. Mitchell, and A. E. Ismail, J.
