Theoretical Reconstruction of Realistic Dynamics of Highly Coarse-Grained cis-1,4-Polybutadiene Melts I. Y. Lyubimov, M. G. Guenza1 Department of Chemistry and Institute of Theoretical Science, University of Oregon, Eugene, Oregon 97403 (Dated: 11 December 2013) The theory to reconstructthe atomistic-levelchain diffusion from the accelerateddynamics that is measured inmesoscalesimulationsofthecoarse-grainedsystem,isappliedheretothedynamicsofcis-1,4-Polybutadiene meltswhereeachchainisdescribedasasoftinteractingcolloidalparticle. Therescalingformalismaccountsfor thecorrectionsinthedynamicsduetothechangeinentropyandthechangeinfrictionthatareaconsequence 3 ofthecoarse-grainingprocedure. Byincludingthesetwocorrectionsthedynamicsisrescaledtoreproducethe 1 0 realistic dynamics of the system described at the atomistic level. The rescaled diffusion coefficient obtained 2 from mesoscale simulations of coarse-grained cis-1,4-Polybutadiene melts shows good agreement with data from united atom simulations performedby Tsolouet al. The derivedmonomer friction coefficient is used as n an input to the theory for cooperative dynamics that describes the internal dynamics of a polymer moving a J in a transientregions of slow cooperative motion in a liquid of macromolecules. Theoretically predicted time 2 correlation functions show good agreement with simulations in the whole range of length and time scales in whichdataareavailable. The theoryprovides,fromdata ofmesoscalesimulationsofsoftspheres,the correct ] atomistic-level dynamics, having as solo input static quantities. t f o s I. INTRODUCTION lations in the long-time regime, but it is also a problem . t becausethedynamicsinthesimulationbecomestoofast a m Computer simulations of macromolecular liquids de- and measured diffusion coefficients are too large. De- pending onthe levelofcoarse-grainingthe dynamics can - scribed at the atomistic level are extremely useful d become orders of magnitude faster than the atomistic because they bridge information between the micro- n dynamics.5,6 scopic molecular structure of the polymeric system, at o c the given thermodynamic conditions, and its macro- Being aware that the dynamics has to be rescaled [ scopic properties, such as viscosity, diffusion, and to recover the correct timescale,7–9 one resorts to the dynamical-mechanical response, which are observed mostconventionalmethod,whichistobuildacalibration 1 experimentally.1,2 It is unfortunate that, despite the curvecalculatedbysuperimposingthe long-timedynam- v 6 progressalreadyoccurredinthecomputationalhardware ics of the coarse-grained and atomistic simulations.10,11 6 andsoftware,thestudyofthedynamicsofpolymericliq- The hope is for this calibration curve to be transferable 1 uidsbyatomisticcomputersimulationsisstilllimitedby to other systems, or to be identical, and so applicable, 0 the impossibility of simulating polymer dynamics in the forthermodynamicconditionscloseto,butdifferentthan 1. wide range of time and length scales relevant for these the phasepointforwhichthe curvewasoptimizedinthe 0 systems. It is known that long simulation trajectories first place. Otherwise, if an atomistic simulation has to 3 deterioratewithtime,followingalawdefinedbytheLya- be ran every time we need to rescale a coarse-grained 1 punov exponent characteristicof the system, and the re- simulation there is no advantage in running simulations v: sulting long-time dynamics is affected by errors.3,4 of coarse-grained (CG) systems. The rescaling factor of i Atomistic simulations of polymeric liquids are limited the conventional calibration curve is, however, purely a X eitherinthe lengthandnumberofmacromolecules,orin numericalcorrection. Thereisnophysicalmotivationfor ar the maximum timescale that can be reached. The long one to assumethat this correctionhas to be identicalfor timescale is a regime of great interest for these systems, any system or thermodynamic condition different from because of the relevance of dynamical-mechanical prop- the ones in which it was optimized. Given that the ef- erties of macromolecules in that regime for their indus- fective potentials between coarse-grained units are free trial applications. Because the limitation concerns only energies,andassuchtheyareparameterdependent,itis the region of long-time, large-scale properties, it is pos- veryunlikely that numericallyoptimized correctionswill sible to shift the focus of a simulation to this region of beidenticalfordifferentsystems. However,ifthecalibra- interestbythereductionofthesimulateddegreesoffree- tion curve is applied for conditions of temperature, den- dom through the averaging of local scale properties, i.e. sity,anddegreeofpolymerization,verycloseinthephase adopting a so-called “coarse-graining”procedure. diagram to the ones in which the calibration curve was While simulations of coarse-grained systems afford originally optimized, there is hope that the error is very to represent larger length- and longer time-scales than contained. Alternatively, in the case that very limited atomisticsimulations,theyalsoaremarredbytheshort- levels of coarse-graining are applied, e.g. combination comingthatthe measureddynamics isunphysicallyfast. of a few atoms to obtain a new coarse-grained unit, the This is the advantage that allows one to perform simu- error incurred in using the calibration curve in different 2 conditions, canbe smallenoughto be within the numer- complex monomeric structure than the simple polyethy- icalerrorofthe calculation. Thisisinfactthe reasonfor lene. A first study of this issue is presented in this theobservedgoodagreementintherescaleddynamicsof manuscript which investigates the agreement between systems with contained level of coarse-graining.12–14 theoretically predicted diffusion coefficients for cis-1,4- The goal of this paper is orthogonal to the more con- Polybutadiene (PB) chains and data from molecular dy- ventional calculation of an optimized calibration curve namics simulations performed by Tsolou, Mavrantzas for systems with limited level of coarse-graining. In fact and coworkers.15–17 The approach is completely general our goal is to investigate the physical motivation that and can be applied to polybutadiene chains with differ- leads to the accelerated dynamics in first place and to ent tacticity. For these systems the difference in local provide, through our analysis, a theoretical approach to chainsemiflexibility,whichcorrespondstodifferentover- calculate, from first principle equations, the correction all chain dimension, i.e. a different radius-of-gyration factorthatneedstobeappliedtothefastdynamicsmea- (Rg) and statistical segment length (l), ultimately leads sured in coarse-grainedsimulations to recoverthe slower to a different diffusion coefficients for polybutadienes andmorerealisticdynamicsofthe atomisticdescription. with different tacticity. In this paper we limit our cal- We study the most extreme level of molecular coarse- culationstocis-1,4-Polybutadienesamplesbecausethose graining, where one polymer chain is represented as a arethesamplesinvestigatedbyatomisticsimulationsand soft sphere, because i) the description is analytical al- provideaeffectivetestofourprocedurefordynamicalre- lowing the formal study of the rescaling problem, ii) it construction. provides a solid test of the reconstruction procedure be- Our rescaling method considers two contributions cause possible errors would be maximized, iii) this level emerging from the consequence of applying the coarse- of coarse-graining is associated with the largest compu- graining process and derives the theoretical corrections tational gain, so it is important. that need to be applied to the dynamics of the coarse- Recently we have proposed an original procedure to graineddescription to recoverthe modalities of the orig- rescale the dynamics measured in simulations of coarse- inal atomistic system. To calculate the analytical cor- grained polymer melts.5,6 The procedure we proposed rection we adopt for the atomistic descriptions a simple doesnotrequiretorunatomisticsimulationstocalibrate bead-and-spring description of the chain where each PB the dynamics of the coarse-grained simulation. Once polymer of N monomeric units is represented as a col- oneparameter,specifically the diameter ofthe hard-core lection of N beads connected by semi flexible springs of monomerpotential,isfixed,therescalingfactorsarefully length l = 6/NRg, and the overall chain dimension determined. Theadvantageofthismethodisclear,given is defined by the radius-of-gyration of the molecule, Rg. p that the coarse-grainedsimulation can reachlarger-scale ThismodelisaRouseapproachmodifiedtoincludechain longer-timepropertiesthantheatomisticsimulation. Be- semiflexibility and has been shown to represent well the cause no atomistic simulationsare needed, the measured diffusion of polymer melts in the unentangled regime.18 dynamicscanbedirectlyrescaledtoobtainthedynamics From the comparison of this model and the soft sphere of the real system we wanted to study. representation of the chain, it is possible to evaluate the The method has been formally derived, and then ap- twocontributionsthatareresponsiblefortheaccelerated plied to systems of polyethylene (PE) melts. The choice dynamics, and which have to be included to correct the of PE as a test system was motivated by the wealth of measureddynamicsandbringit tothe realisticvaluesof experimental and simulation data available in the lit- the atomistic (here bead-and-spring) description. erature. We have shown that the proposed theory is The first correction emerges from the consideration able to predict diffusion coefficients in agreement with that due to the process of averaging the intramolecu- both atomistic simulations and experiments for systems lardegreesoffreedom,the coarse-grainedsystemexperi- in both the unentangled and entangled regimes, for a ences a change in its entropy, because a number of local range of temperatures and densities. For entangled sys- statesareneglected. While the atomistic systemdevotes tem the approach of dynamical reconstruction includes time to explore local configurations, the coarse-grained anextraloopinthecalculationofthefrictioncoefficient, system doesn’t. To recover the correct dynamics those whichaccountsforthefactthatinrealsystemsdescribed states need to be included a posteriori in the form of an at the atomistic level, the dynamics is slowed down by entropiccontributionto the free energythat rescalesthe the presence of entanglements. We hypothesized for the time measured in the CG simulations. entanglements to relax in time following the same dy- The second correction comes from the change of the namicalmechanismof single chaininterdiffusion. In this “shape”ofthemoleculeonceitiscoarse-grained. Macro- paper entangled chain dynamics is not investigated as molecules described at the atomistic level become rep- most of the samples in the atomistic simulations are ei- resented by chains of effective atoms, or as bead-and- ther unentangled or in the region of crossover from the springs, chains of soft blobs, or even each molecule as a unentangled to the entangled regimes. softsphere. Inallthesemultipleformsthesurfaceofeach One of the questions that still needed to be answered molecule available to the surrounding molecules, i.e. the was if our method was able to generate high quality “solvent”, is different. The hydrodynamic radius, r , hydr predictions also for macromolecular liquids with a more of the coarse-grained unit in each description is differ- 3 ent and so is its friction coefficient, ζ, defined by Stokes’ a function of the atomistic description;third all the con- law as ζ = 6πη r , where η is the viscosity of the sequences of the coarse-grainingprocedure are enhanced s hydr s medium. Bysolvingthememoryfunctionforthefriction with this extreme level of coarse-graining and are easier coefficientineachdescriptionweareabletocalculatethe to study in such a description. If the procedure adopted correction factor that has to be applied to each friction were not precise, the error would be maximized because coefficient to recover the description at a different level of the high level of coarse-gaining;the fact that the the- of coarse-graining. ory, which is predictive, is found to be in quantitative In this paper our rescalingmethod is firstbriefly sum- agreementwithbothsimulatedandexperimentaldatais marized and then applied to study the dynamics of cis- encouraging.5,6 1,4-PB and compared with simulation data by Tsolou One disadvantageofhaving sucha coarse-grainedrep- et.al.15 who also investigated the effect of temperature resentation of the molecule is that no information is col- and pressure on the simulations.16,17 The global dynam- lected from the mesoscale simulation on the internal dy- ics of the polymer is represented by the diffusion coeffi- namics of the polymer. In this way, the soft sphere cientand the decay ofthe self-correlationforthe end-to- representation only provides information on the diffu- endmolecularvector. Thesedynamicalquantitiesarere- sion coefficient of the molecule. However, we show in constructedusingthemethodbyLyubimovatal.,5,6from this paper that from the knowledge of the diffusion co- theinformationcontainedinthetrajectoriesofthemeso- efficient and the related monomer friction coefficient for calesimulationsofthe coarse-grainedpolymermelts and unentangled chains it is possible to recover correctly the show good agreement with the data from united atom internal dynamics of the chain by applying our theory simulations. For the internal dynamics the theory for for the cooperative dynamics of a group of interacting cooperative motion19–21 is used to predict monomer dy- chains, which calculates the single chain dynamics for a namics and the dynamical structure factors, which com- macromolecule whose dynamics is coupled by the pres- pare well with the simulations. If atomistic simulations ence of intermolecular forces with the other interpene- arenotavailable,the theoryofdynamicalreconstruction trating macromolecules in the liquid. is able to predict the correct dynamics from the rescal- IneveryCGmodel,theaveragingoftheintramolecular ing of the mesoscale simulations having as an input the degrees of freedom leads to a speeding up of the dynam- polymerradius-of-gyration,whichcanbecalculatedfrom ics. When a polymer is represented as a soft sphere the simple structural models for polymers. The accuracy of dynamics in the mesoscopic simulation is orders of mag- adoptingafreelyrotatingchainmodelforthecalculation nitude faster than in the atomistic description. One of of the polymer radius of gyrationis discussed in the last the reasons for this acceleration is the fact that while a section. A brief summary concludes the paper. macromolecule needs to explore a large number of inter- nalchainconfigurationsforeachpositionofthecenter-of- mass,thesoftsphereinsteadisfreetomoveinthethree- dimensional space undergoing Brownian motion.23,24 To II. RESCALING OF THE FREE-ENERGY AND take into account the missing contribution of the time RESCALING OF THE FRICTION necessary to explore the internal degrees of freedom, the time measured in the mesoscale simulation of the soft In this section we briefly present an overview of the sphere is multiplied by a rescaling factor that is calcu- main steps in our procedure. A complete and detailed lated as the ratio of the energy due to the internal de- presentation of the same has been already published5,6 grees of freedom in the two representations,i.e an atom- and will not be repeated here. We start from the con- istic and a soft sphererepresentations. For the atomistic sideration that the procedure of coarse-graining corre- representation we adopt an analytical bead-and-spring sponds to eliminating some internal degrees of freedom, model, as previously described.5,6 This model has been bycombininggroupsofatomsintooneeffectiveCGunit. shown to provide a quantitative description of the dy- Specifically, our rescaling theory has been developed for namics for polymeric liquids and for proteins in theta an extended level of coarse-graining where intramolecu- solvent.21,22 larcoordinatesarefully averagedoutandthe polymer is The first correction term is calculated from the def- representedasanisotropicspherecenteredonthecenter- inition of the internal energy for a simple bead-and- of-mass(com)ofthemacromolecule. Therearesomead- spring model and for the soft sphere CG model. The vantages in the choice of this description. First, the CG mesoscale(MS)moleculardynamic(MD)simulationsare representation is fully isotropic, and even if it is known performed using a potential expressed in dimensionless that the shape of a polymer is closerto anellipsoidthan units, which is a common procedure. The time from the to a sphere,18 the total correlation function of the poly- simulationshastobeproperlynormalizedsothatthedi- meric liquid, i.e. the structure of the polymeric liquid, mensionlesssimulationtimefromMSsimulations,t˜,once statistically averagedover all the possible configurations dimensionalized and rescaled reads as is well reproduced by this model; second the formalism is analytical, which allows us to derive formally many of the physicalquantities ofinterestforboth the static and m 3 t=t˜R N , (1) the dynamic properties of the coarse-grained system as g k T 2 r B 4 where m is the molecular mass of one chain and R is The contributions that are still unknown are the friction g the radius-of-gyration of the molecule, which is an in- coefficients in the atomistic and coarse-grained descrip- put quantity in our approach. The internal degree-of- tions, which can be solvedanalytically starting from the freedom averaged out in the coarse-grained description definition of the memory function.25 areaccountedforby introducingthe 3N factorwhere N Forabead-and-springmodelthemonomer(bead)fric- 2 is the number of beads. tion is given by: The secondrescalingis due to the changeinthe shape of the molecule represented in the two levels of coarse- N ∞ graining, and so the effective friction in the two repre- 1 ζ = dτΓ (τ) . (4) sentations. In the long time the com mean-square dis- m ∼ N a,b a,b=1Z0 placementisrelatedtothemoleculardiffusioncoefficient X as where the memory function is defined as26 h∆R2(t)i=6Dt , (2) Γa,b(t)∼= β3ρ dr dr′g(r)g(r′)F(r)F(r′)ˆr·ˆr′ Z Z where D =k T/(Nζ ) is the diffusion coefficient of the B m dRSQ (R;t)SQ(r r′+R;t), (5) system,whichhastobecalculated. Ifweuseasastarting × a,b | − | point the diffusion coefficient measured in the mesoscale Z simulationDMS =k T/ζ ,themeansquaredisplace- where β = 1/k T, ρ is the monomer density and ˆr, ˆr′ B soft B ment in the rescaled formalism is give by are the unit vectors characterizing the direction of the forcesactingonmonomersaandb,respectively. SQ(r,t) istheprojecteddynamicstructurefactor,whichincludes ∆R2(t) =6DMS ζsoft t , (3) both intra- and inter-molecular contributions, i.e. inco- h i Nζ herent and coherent scattering, and is approximated by m its non-projected form as SQ(r,t) S(r,t), which is an ≈ where the problem reduces to the calculation of the fric- acceptableapproximationwhentheLangevinequationis tioncoefficientinthemonomerandsoft-sphererepresen- expressed as a function of the slow variables.6 tations. The friction coefficients are calculated from the In Eq.(5) g(r) is the monomer pair distribution func- solution of the memory function in the two representa- tion of the molecular liquid, and F(r) is the force due tions (bead-and-spring and soft-sphere). to the effective potential between two monomers, ob- The formalism presented so far, adopts a Rouse-like tained from the solution of the Ornstein-Zernike (OZ) description of the single chain diffusion for the atomistic equation by applying the Percus Yevick closure. The levelrepresentationofthe system. Atthe coarse-grained monomer potential, which in the atomistic-level simula- levelthechainisrepresentedasasoft-sphere. Eq.3shows tion is a Lennard-Jones potential, is approximated here howthediffusioncoefficientmeasuredinamesoscalesim- by an effective hard-sphere with an effective diameter, ulation of a liquid of soft-spheres, DMS, needs to be d, to mimic the properties of the L-J potential. All the rescaledtogivethemean-squaredisplacementofthereal physicalquantitiesthatappearintheequationareknown chain, ∆R2(t) , as a function of the rescaled time, t. andtheanalyticalexpressionforthemonomerfrictionis6 h i 2 2 1 d d ζ (Dβ)−1ρg2(d) πN2d2R 15√2+40 +12√2 m g ≈ 3 12 " Rg (cid:18)Rg(cid:19) # 1 ξ 2 +ρπNh 12√2ξ7+12d4R3 1 2 ρ 03√2(Rg2−2ξρ2)2" ρ g − (cid:18)Rg(cid:19) ! 2 4 2 5 ξ ξ ξ ξ +4√2d3R4 5 14 ρ +2 ρ +3d2R5 5 14 ρ 4√2 ρ g − (cid:18)Rg(cid:19) (cid:18)Rg(cid:19) ! g − (cid:18)Rg(cid:19) − (cid:18)Rg(cid:19) ! 2 12√2e−ξ2ρdξ7 1+ d +ρ2πh2 1 40d3R6+15√2d2R7 − ρ(cid:18) ξρ(cid:19) # 012(Rg2−2ξρ2)3" g g ξ 2 24√2d4R3ξ2 144d3R4ξ2+6√2d4R5 2 9 ρ − g ρ − g ρ g − (cid:18) d (cid:19) ! 5 3 2 3 2 d d d d +12R2ξ7 4 7 +9 8ξ9 4 9 +15 g ρ (cid:18)ξρ(cid:19) − (cid:18)ξρ(cid:19) !− ρ (cid:18)ξρ(cid:19) − (cid:18)ξρ(cid:19) ! −e−ξ2ρd12ξρ4(d+ξρ) Rg2(d+3ξρ)(2d+3ξρ)−2ξρ2(2d2+5ξρd+5ξρ2) #! . (6) (cid:0) (cid:1) This equation contains quantities that are well defined 507 3 1183 679√3 + ρ h + ρ2 h2 , (10) once the system of interest is selected. For example ×512"r2 507 ch 0 1024 ch 0# g(d) is the pair distribution function at contact, ξ = ρ Rg/(√2+2πRg3ρ/N) is the density fluctuation correla- wherethevalueofeachphysicalparameter(temperature tionlengthwithρ,Rg andN alreadydefinedinthetext, and density) and molecular parameters (degree of poly- andh0 isdefinedash0 =h(k =0)=−(1−2ξρ2/Rg2)/ρch. merization and radius-of-gyration) that enter this equa- The only physical quantity that needs to be determined tion is defined once the system to simulate is selected. is the hard-sphere diameter, d. This is defined once the Theseparametersareidenticaltotheonesthatareinput Lennard-Jones potential is mapped onto an effective re- to the equation of the monomer friction, Eq.(6), given pulsive hard-spheresystemas describedlaterin this sec- that the two equations are representations for different tion of the paper. levels of coarse-grainingof the same system. In the soft colloidal representation, the friction coeffi- Whenthe twoequationsforthe monomerandthe soft cient is analytically calculated from the definition of the colloidalfrictioncoefficients areintroducedinEq.(3) the memory function as diffusioncoefficientfora singlechainina liquidis recov- ∞ ered as ζ =(β/3)ρ dt dr dr′g(r)g(r′)F(r)F(r′)ˆr ˆr′ soft ∼ chZ0 Z Z · D =DMSζsoft/(NζmRg 3mN/(2kBT)), (11) dRS(R;t)S(r r′+R;t), (7) p where we have assumed that the diffusion coefficients × | − | Z that appear in Eq.6and Eq.10are identical, as the long- where ρch = ρ/N is the chain density. All the other time relaxation of the dynamic structure factor in both physicalquantitiesthatappearinthe integral,forexam- CG descriptions is guided by molecular diffusion. ple the pair distribution function, g(r), and its related To summarize, as the first step the thermodynamic force, F(r), are now defined in relation to the descrip- and molecular structure, i.e. the radius of gyration or tion of the polymeric liquid as a liquid of soft spheres, equivalentlythepersistencelength,ofapolymermelthas which are point particles interacting through a soft re- to be defined. Then, MS MD simulations are performed pulsive potential. Each of these quantities is analytical for a liquid of point particles interacting through a soft andcalculatedfromthesolutionoftheOZequationwith pair potential, βv (r), as described in several of our soft the hypernetted chain closure approximation.27 previous papers, where The regime of interest in our calculations is the dif- fusive limit where, in reciprocal space, the wave vector βvsoft(r)=hsoft(r) ln[1+hsoft(r)] csoft(r) , (12) − − k 1/R . Here the dynamic structure factor can be g ≪ with approximated as c (k)=h (k)[1+ρ h (k)]−1 . (13) soft soft ch soft S(k;t) S(k) e−Dtk2 =(1+ρ h (k)) e−Dtk2, (8) ≈ ch soft Themean-square-displacementofthesoftspheresismea- suredfromthe MSMD trajectoriesasa function oftime withthetotalcorrelationfunctionofthesoft-sphererep- (dimensionless), and the resulting diffusion coefficient is resentation is given in the limit of long chains (N 30) ≥ rescaled following Eq.11 to obtain the reconstructed dif- by the approximated expression27 fusion coefficient. For the solution of the rescaling factors in Eq.11 the 39 3 ξ ξ 9r2 −3r2 only parameter that has to be defined is the value of hsoft(r)≈−16 πRρ 1+√2Rρ 1− 26R2 e 4R2g . the effective hard sphere diameter, d, in Eq.6. In our r g (cid:18) g(cid:19)(cid:18) g(cid:19) formalismtheelementaryinteractionbetweenmonomers (9) belonging to different chains, which in the UA MD is The resulting equation for soft particle friction is6 a Lennard-Jones potential, is approximated by a hard- sphere interaction with an effective diameter, d, to allow 2 forthe analyticalsolutionofthemonomermemoryfunc- √2ξ ζ =4√π(Dβ)−1ρ R ξ2 1+ ρ tion. Under fixed thermodynamic conditions, the value soft ∼ ch g ρ Rg ! ofdshoulddepend onlyonthe localmonomerstructure, 6 it could be argued that in our method the actual fitting parameter is not d, but the length of the chosen sam- ple for which d is determined. Clearly when entangled chains are forced to behave according to the Rouse ex- pressionthehardspherediameterneedsto beartificially decreasedinordertocompensatetheoverlookedincrease in the monomer friction, which is associated with entan- glement effects growing with increasing N. Whereas the range of lengths in the unentangled regioncan be rather wide, it is always possible to choose the chain length knowing experimentally measured or theoretically esti- mated value of N ,28 and by optimizing our choice of e d by testing the predictions of the theory against data from simulations of short chains. However comparison with simulations is not necessary, as the only needed in- formation is the value of N and the value of R for a e g sample of short, unentangled, chains (N <N ).6 e FIG.1. Dimensionlessmonomerfrictioncoefficientasafunc- In our previous study of PE melts we fixed the hard tion of the hard sphere diameter, based on Eq.(6). a) cis- sphere diameter following the same procedure described 1,4-Polybutadiene samples with N= 32, 56, 128, 320 (solid, here, and we selected a sample with N = 44 to cal- dashed, dot-dashed, and dotted lines correspondingly) and parametersasreportedinTableIandref.15. b)Polyethylene culate d. The obtained value of d = 2.1˚A (see Fig- samples with N=30, 44, 96, 270 (solid, dashed, dot-dashed, ure 1b), was larger than the carbon-carbon bond length anddottedlinescorrespondingly)andparametersasreported lb = 1.54˚A, but smaller than the Lennard-Jones param- inrefs.5and6. Horizontallinesrepresent1/N values,following eter σ 3.9˚A typically used in UA MD simulations, ≃ thediffusioncoefficientforunentangledchains,Dβζm =1/N. which seems to be a reasonable choice of the parame- ter. Compared to PE, PB chains are more flexible and so they are less entangled. The estimated entanglement and be independent of the degree of polymerization, be- degree-of-polymerization,N ,iscalculatedfromthegen- e cause the monomer interaction potential is a local prop- eralizedformula for the plateau value in the shear relax- ertyofthe chain. Thereforeinourmodel,once the value ation modulus,28 which gives N 330 backbone e(shear) ≃ ofdischosen,thisvalueiskeptfixedforallsampleswith carbonatomsforPB,comparedtoN 55forPE. e(shear) ≃ different molecular weights. Note that for convenience Notice that the value of N from shear measurements is e we use the term monomer to identify the CH group, different from the value obtained from the analysis of x where x=1, 2 or 3. scattering experiments which gives, for example, for PE To evaluate the value of d we analyze the dimension- samples Ne 130. ≃ less quantity Dζ /(k T), which is identical to N−1 if It is important to notice that even if there is some m B the chain obeys strictly Rouse dynamics in the long- freedomin choosing the reference sample, the deviations time limit. To fulfill this requirement the chain has to of d from the chosen value is not large. However, small have N much smaller than the entanglement degree-of- deviations of the value of d can lead to different values polymerization N and larger than the value of N = 30 of the diffusion coefficient. Figure 2 displays a compari- e forwhichchainsstarttoobeyGaussianstatistics,i.e. the son of the predicted diffusion coefficient as a function of central limit theorem applies. N when different values are chosen for the hard-sphere Figure 1 a) displays the dimensionless friction coeffi- parameter d. It is possible to see that for values of d cient,Dβζ ,fromEq.(6)asafunctionofthehardsphere obtainedbyenforcingRousedynamicsforthreedifferent m diameter, d, for PBsamples ofincreasingdegreeofpoly- chain lengths that obey diffusive Rouse dynamics in the merization,N. Theparametervalues(ρ ,T,N,R )are long-time regime, the predicted diffusion coefficients are m g chosentomatchthoseofUAMDsimulationsinref.15 as all in reasonable agreement with the values of the UA reported in Table I. Horizontal lines represent the 1/N MD simulation. value for each sample, which is the Rouse value of the Notwithstanding the fact that for PB there is a wide physical quantity Dβζm = 1/N in the diffusive regime, range of degrees-of-polymerization, 30 < N < 330, for and reproduces the scaling behavior of the diffusion for which the hard-sphere diameter, d, could be optimized, unentangled, short-chains. For comparison Figure 1 b) thebestagreementbetweentheoryandsimulationsisob- displays the analogous plot for the PE samples, which is tained for values of the parameter d optimized for short reproduced from refs.5 and6. chains, N 50. This can be explained by considering ≈ The change of d with N observedin Figure 1 is a con- that in PB samples the crossover region from unentan- sequenceoftwodifferenteffects: forshortchains,endef- gled to entangled dynamics is extended (see for example fectsareimportant,whileforlongchainsthecrossoverto Figure3,wheredeviationsfromtheN−1scalingoccural- entangled dynamics starts to be important. In this way, readyforN 100). SampleswithN approachingN are e ≈ 7 pare our results with the data from the simulations per- formed by Tsolou et al.15 The input parameters to our theory are displayed in Table I, including the degree-of- polymerization, N, and the monomer density ρ. Table I also includes the mean-square-radius-of-gyration calcu- lated from the united atom simulations R2 UA. From g thisradius-of-gyrationthevalueofthesemiflexibilitypa- (cid:10) (cid:11) rameter, g, is derived, as the value that corresponds to chains with the desired overall dimension, i.e. R . In g this way the parameter g is not an independent param- eter but is determined by the R value. The parameter g g does not enter the equation, but is used in the cal- culations reported in the last section of this paper. It is importanttonoticethatthevaluesforR couldbetaken g fromexperimentaldata,and that, in principle, no atom- isticsimulationsofthe systemunderstudy arenecessary for our approach. FIG. 2. Diffusion coefficient predicted from the rescaled MS The UA MD simulations were performed in the NPT MD, when different values of the hard-sphere diameter are ensemble, at the constant pressure, P, and tempera- chosen as an input to the reconstruction procedure. Assum- ture, T, reported in the table.15 Note that the values ing different values of d, calculated by enforcing Rouse dif- fusive behavior for N=56 (d=1.4672), 64 (d=1.4051), or 80 of Rg for systems PB160, PB200, PB240 were corrected (d=1.374) (circles, diamond, triangles correspondingly) leads from 257˚A2 to 275˚A2, from 304˚A2 to 340˚A2 and from to diffusion coefficients in good agreement with the UA MD 401˚A2 to 410˚A2 correspondingly,after consultation with simulation data (filled squares). the authors of ref.15. The correction was motivated by the inconsistency observed with the data reported in reference15 when the calculation of R is performed us- g not following the unentangled scaling behavior, which is ingastandardmodel,suchastheFreely-Rotating-Chain, necessaryconditionfortheoptimizationofdusingRouse- described later in in this paper. like diffusion. Considering the broader range of unentangeled PB chains, we select for the calculation of the hard sphere TABLE I. Simulation parameters for 1,4-cis PB chains of in- creasing lengths. diameter the PB sample with N = 56. The obtained value of d 1.47˚A (which is smaller than the value for PEof2.1˚A≃giventhe flexibility ofPB)is consistentwith System ρ[sites/˚A3] R2 UA[˚A2] g g the fact that in PB chains half of the carbon bonds are PB 32 0.0352375 45 5 0.6237 shorter double bond and l =1.34˚A. (cid:10) (cid:11)± db PB 48 0.0363767 70 5 0.6243 ± PB 56 0.0367159 85 7 0.6339 ± PB 64 0.0369744 95 10 0.6246 III. THEORETICAL PREDICTIONS OF ± CENTER-OF-MASS DIFFUSION PB 80 0.0373425 125 10 0.6374 ± PB 96 0.0375921 152 16 0.6394 ± The solution of Eq.(3) gives the diffusion coefficient PB 112 0.0377723 184 15 0.6490 ± for the center-of-mass of a polymer described at the PB 128 0.0379087 215 18 0.6544 ± atomistic level, having as an input the diffusion coef- PB 140 0.0379910 234 18 0.6523 ± ficient calculated from the mesoscale simulation of the PB 160 0.0381012 275 25 0.6595 ± CG polymer liquid. The purpose of this procedure is to PB 200 0.0382567 340 20 0.6551 have ultimately mesoscale simulations that, once prop- ± PB 240 0.0383610 410 30 0.6557 erly rescaled, can provide directly the diffusion coeffi- ± PB 320 0.0384923 576 30 0.6695 cient, without need of running atomistic simulations. To ± PB 400 0.0385714 678 30 0.6519 test our procedure we first run the mesoscale simula- ± tions, then we rescale the diffusion coefficient, and fi- T=413K, P=1 atm nallywecomparethepredictedvaluesofDagainstatom- istic simulations and/or experiments. Thermodynamic In the mesoscale simulations of a polymer melt, each andmolecularparametersenteringourequationshaveto chain is represented as a point particle interacting be consistent with the parameters of the system against throughasoft-corepotentialderivedfromthesolutionof which we compare our approach. theOrnsteinZernikeequationapplyingtheHyperNetted We study the dynamics of cis-1,4-Polybutadiene(PB) Chain closure. MS MD simulations were implemented melts of increasing degree of polymerization and com- in the microcanonical (NVE) ensemble on a cubic box 8 with periodic boundary conditions. We used reduced units such that all the units of length were scaled by ∗ R (r = r/R ) and energies were scaled by k T. More g g B details of our simulation procedure have been reported in previous papers.5,29,30 Table II reports the diffusion coefficient directly cal- culated from the MS MD in the soft sphere represen- tation, DMS, and the diffusion coefficient calculated us- t ing the dynamical reconstruction procedure, D. Finally, for comparison, the Table shows the diffusion coefficient measuredintheunitedatomsimulation,DUA. Thediffu- sion coefficient measured in the mesoscale simulations is several orders of magnitude faster than in the atomistic simulations. However, once it is rescaled the diffusion coefficient becomes very similar to the one obtained di- rectly from the atomistic simulation. It is important to notice that, once the parameter d, which is characteris- ticofthepolymerconsidered,isdeterminedthediffusion coefficient is calculated without any input from the dy- FIG. 3. Center of mass self-diffusion coefficient as a func- tionofdegreeofpolymerization,N,forcis-1,4-Polybutadiene namics of the atomistic simulations, so the procedure is meltswithparametersdefinedinTableI.Diffusioncoefficients predictive. reconstructedfromMSMDbyapplyingourprocedure(open The samples here are relatively short and these cal- symbol) are compared against UA MD data (filled symbol) culations are performed for unentangled and slightly en- from reference15. In analogy with thefigure from thesource, tangledsystems. Forstronglyentangledsystems,wepro- three scaling regimes in terms of power law dependence of posed in a previous paper a perturbative version of our D Nb are shown as dashed (b > 1), solid (b 1), and ∝ ≈ approachthataccountsforthe factthatboththetagged dot-dot-dashed(b 2) lines. ≈ chain and the surrounding chains, that provide the en- tanglements,relaxfollowingthe same dynamics.6 Inthis manuscript the PB samples are in the unentangled and bols representMS MD simulations rescaledas in Eq.(3). slightlyentangledregimes,whicharewellrepresentedby In analogy with the figure in Ref.15 from which the UA the theory without the one-loop perturbation. MD data are taken, we report three scaling regimes in termsofpowerlawdependence ofD N−b.ForN <80, ∝ there is a faster than Rouse regime with D with b > 1, TABLE II. Diffusion coefficient reconstructed from MS MD whichisattributedtothefree-volumeeffectsduetochain simulation compared against UAMD simulations. ends which is significant only for very short chains31–33. Intheintermediateregimefor80<N <160aRouse-like N DMS[˚A2/ns] D[˚A2/ns] DUA[˚A2/ns] behaviorcanbeobservedwherethescalingexponentbis t 32 3875 184.5 152.2 close to 1. For longer chains with N > 200 the value of b 2.1 is typical of the crossoverto reptation dynamics. 48 3400 68.1 69.4 ≈ The data in this plot are not closely following the re- 56 3425 51.9 46.4 ported scaling exponents because as the degree of poly- 64 3275 37.3 40.1 merization increases the density of the system changes 80 3174 25.1 24.9 (Fig.4). The different regimesare evenless prominentin 96 3114 18.0 23.2 the MS MD simulations, where we observe a smoother 112 3235 15.1 16.7 transition between regimes than in UA MD data. In 128 3283 12.6 15.5 our calculations features like free volume effects due to 140 3221 10.5 12.0 chain ends (causing faster than Rouse decay of the dif- 160 2963 9.1 11.0 fusion coefficient for short chains) enter only indirectly 200 3245 6.0 6.6 as a consequence of the interplay of the input parame- ters, i.e. ρ, R , N. The increasingof the density with N 240 3012 4.5 5.4 g enters our analyticalequationboth directly, throughthe 320 3386 3.2 2.6 liquiddensity,andindirectlythroughtheeffectiveradius 400 3016 1.8 1.6 of gyration of the polymer, given that the compactness d=1.4672˚A(PB56) ofthe chainchangeswith the density (see alsoSectionV of this paper). On Figure 3 the diffusion coefficient is presented as a The cis-1,4-PBchains are more flexible in comparison function of degree of polymerization N. Filled symbols to PE chains, therefore we should expect larger entan- representUA MDsimulations fromref.15 andopensym- glement length N (for PE from scattering N 130). e e ≃ 9 erativemotionof a groupof polymer chains in a dynam- ically heterogeneous liquid, as observed in simulations data of polymer melts and experiments.35–38 Theoreti- cal predictions of the com and monomer mean-square displacements are shown to be in quantitative agree- mentwith computer simulationsofunentangled20,39 and slightly entangled polymer melts,40 and with scattering experiments of Neutron Spin Echo.21 The physical picture underlying the theory builds on thefactthatinpolymermeltsthedynamicsappearshet- erogeneous with the motion of a tagged chain correlated to the dynamics of a group of n chains comprised in- sidetherangeoftheintermolecularpotential. Thelatter has a range of the order of the correlation hole, i.e. of theradius-of-gyrationofthe polymer,asitemergesfrom the solution of the OZ equation. The number of chains comprised in a volume of the order of R3 is given by g n ρN1/2l3,wherel,thestatisticalsegmentlength,and FIG.4. Densityasafunctionofdegreeofpolymerization,N, ≈ the other quantities have been defined earlier on. The for cis-1,4-Polybutadiene samples reported in Table I. number of interpenetrating chains, n, increases with in- creasing density, degree of polymerization, and the stiff- ness of the polymer. Figure 3 shows that crossover from Rouse to reptation In the Cooperative Dynamic approach the dynam- regime starts at N =200. The largest system shown on ics is described by a set of coupled equations of mo- Figure 3 is for N = 400 which is still a very weakly en- tion (eom). Each equation is expressed in the space tangledsystem(1or2entanglementsperchain). Overall coordinates of the monomer a, belonging to molecule the agreement between simulated and reconstructed dif- i and positioned at r(i)(t), and contains a balance of fusion coefficients is quite good. a different forces acting on the monomer. These in- cludethe viscousforce,ζdra(i)(t), theintramolecularforce dt IV. INTERNAL DYNAMICS −ks Nb=1Aa,brb(i)(t),which contains the structural ma- trix A, the time-dependent intermolecular mean-field forcePβ−1 ∂ ln n g(r(l)(t),r(k)(t)) , and the ran- The coarse-grained model adopted in this paper sim- ∂ra(i)(t) k<l plifies the macromolecular structure by reducing the de- dom interactions w(cid:2)iQth the surrounding li(cid:3)quid, given by scription of the molecule to an effective site that is cen- the random force F(i)(t). a tered on the center-of-mass of the molecule. For this reason, the only dynamics that can be predicted by this CG description, through the rescaling procedure, is the dr(i)(t) N ζ a = k A r(i)(t) long-time diffusive behavior of the center-of-mass of the eff dt − s a,b b molecule. However from the knowledge of the diffusion Xb=1 coefficient the monomer friction coefficient can be ob- 1 ∂ n + ln g(r(l)(t),r(k)(t)) tained, and the latter can be used as an input to the Langevin equation that describes the internal dynamics β∂ra(i)(t) "kY<l # of the polymer chain at any length scale of interest. +F(i)(t) . (14) a SpecificallywefollowourapproachfortheCooperative Dynamics of a group of macromolecules, or Cooperative Here A is the matrix of intramolecular connectivity, Dynamics Generalized Langevin Equation (CDGLE).34 which reduces to the Rouse matrix when infinitely long Conventional theories of polymer dynamics, such as the and completely flexible macromolecules are considered, Rousetheoryforunentangledpolymerdynamicsandthe and is defined as “reptation”modelfor entangleddynamics, predictdiffu- N sive motion at short time following the ballistic regime; A = M U−1M . (15) this is howeverin disagreementwith simulationsand ex- i,j k,i k,p p,j periments, which show a sub-diffusive regime before the kX,p=2 com starts to follow Brownian motion. Both theories HereUisthe equilibriumaveragedbondcorrelationma- aremean-fieldapproachesofsingle chaindynamics inan trix uniform bath.18 Our approach focuses on the dynamics of a group of polymer chains in a melt and relates the <lk lp > U = · , (16) anomaloussubdiffusivebehaviortothepresenceofcoop- k,p l l k p 10 andMisastructuralmatrix,withalltheelementsequal zero except M = 1/Nfor i = 1,...,N, M = 1 with 1,i i,i i = 2,...,N, and Mi−1,i = 1 with i = 2,...,N. The U − matrixisafunctionofthelocalsemiflexibilityparameter g =< l l > /(l l ), which is related to the persis- i i+1 i i+1 · tence length of the polymer. For our samples the values of g are reported in Table I and are calculated from the molecule radius of gyration,which can be obtained from experiments or from atomistic simulations.6 Through Eq.(16) the eom includes a complete micro- scopic description of the structure and localflexibility of the specific molecule under investigation, containing all the relevant parameters that define the intramolecular mean-force potential. Equations of motion for different beads belonging to the same chain or to two chains un- dergoing slow cooperative dynamics are coupled by the presence of intra- and inter-molecular interaction poten- FIG.5. End-to-endvectortimedecorrelationfunctionforcis- tials. Theintrinsic,chemicaldependentsemiflexibilityof 1,4-Polybutadiene with N=112. Data from UA MD simula- the macromolecule enters explicitly through the descrip- tions (symbols) are compared with predictions of the theory tion of the A matrix. for Cooperative Dynamics (solid line) where the number of The intermolecular potential is time dependent, as correlated chains is set to n = 7 and the monomer friction it is a function of the relative position coordinates of coefficient is reconstructed from MS MD simulations, using the centers-of-mass of a pair of molecules. As the two theprocedure described in this paper. molecules move with respect to each other the inter- molecular force changes. From the initial ensemble of n dynamically correlated chains, molecules diffuse in time erativeDynamictheorywiththereconstructedmonomer and finally move outside the range of the potential, Rg, friction coefficient is directly compared against UA MD in a timescale of τdecorr =Rg2/D. data. This function represents the rotational relaxation The set of coupled equations is solved after transfor- of the chain and cannot be measured from the soft col- mation into normal modes of motion and numerically loid simulationdirectly. The N =112 sample has chains using a self-consistent procedure that calculate the ef- shorter than the entanglement length and the Coopera- fective potential at any given distance. Once n and tive Dynamics theory reproduces the atomistic simula- ζm =kBT/(ND) are defined, the firstoptimized against tion data rather well. the simulation data and the second from our rescaling Figure 6 shows the center of mass mean square dis- procedure,the setofequationshasnoadjustableparam- placement (com MSD) for three samples of cis-1,4- eters. Polybutadiene melts with chains of increasing degree of By applyingconsiderationofsymmetry the setofcou- polymerization, N = 240,320 and 400. The long time pled eoms in normal coordinates can be simplified and diffusive dynamics is calculated from the diffusion coef- reducedtotwosetsofN uncoupledequationsintherela- ficients obtained from the rescaling procedure. For time tiveandcollectivenormal-modecoordinates. Weassume shorter than the longest relaxation time, τ R2/D, the that the eoms for internal modes (p>0) do not contain mean-square displacement of the com show≈s a sgubdiffu- intermolecular forces. This approximationis justified on sive behavior that cannot be reproduced by the Rouse the basis that polymer local dynamics is affected in a approach, i.e. the conventional theory of unentangled similar way by inter- and intra-molecular excluded vol- chaindynamics. In the cooperative dynamics theory the ume interactions, which in polymer liquids tend to com- subdiffusivedynamicsisaconsequenceofthecooperative pensate each other. Intermolecular forces, which enter motion of the interpenetrating chains inside the correla- the dynamics through the eom for the first (p = 0) nor- tion hole region, coupled by the effective intermolecular mal mode, still perturb the dynamics on the local scale potential between the coms of the macromolecules. The through the linear combination of the normal modes. theory,withthe monomerfrictionfromrescaledMS MD In the following, we present an overview of some of simulationsasaninput,iscomparedagainsttheUAMD the quantities analyzedin the originalsimulations15 and simulation data. The diffusion coefficient for PB240 is their comparison with the predictions of the Coopera- underestimated and for PB320 is slightly overestimated tive Dynamics theory having as an input parameter the by the theory in comparisonwith the simulations,which reconstructedfrictioncoefficientcalculatedfromthesoft- is also can be seen in Table II and Figure 3. These de- sphere simulations and the rescaling procedure as de- viations are within the error related to R . Dot-dash g scribed in the first part of this paper. lines in the inset of Figure 6 represent results calculated Figure 5 presents the end-to-end vector, time decorre- with the upper and the lower values of R taken from g lationfunctionfortheN =112sampleofPB.TheCoop- the data of the atomistic simulations, as reported in Ta-