Theoretical Reconstruction of Realistic Dynamics of Highly Coarse-Grained cis-1,4-Polybutadiene Melts
The theory to reconstruct the atomistic-level chain diffusion from the accelerated dynamics that is measured in mesoscale simulations of the coarse-grained system, is applied here to the dynamics of cis-1,4-Polybutadiene melts where each chain is described as a soft interacting colloidal particle. The rescaling formalism accounts for the corrections in the dynamics due to the change in entropy and the change in friction that are a consequence of the coarse-graining procedure. By including these two corrections the dynamics is rescaled to reproduce the realistic dynamics of the system described at the atomistic level. The rescaled diffusion coefficient obtained from mesoscale simulations of coarse-grained cis-1,4-Polybutadiene melts shows good agreement with data from united atom simulations performed by Tsolou et al. The derived monomer friction coefficient is used as an input to the theory for cooperative dynamics that describes the internal dynamics of a polymer moving in a transient regions of slow cooperative motion in a liquid of macromolecules. Theoretically predicted time correlation functions show good agreement with simulations in the whole range of length and time scales in which data are available. The theory provides, from data of mesoscale simulations of soft spheres, the correct atomistic-level dynamics, having as solo input static quantities.
Computer simulations of macromolecular liquids described at the atomistic level are extremely useful because they bridge information between the microscopic molecular structure of the polymeric system, at the given thermodynamic conditions, and its macroscopic properties, such as viscosity, diffusion, and dynamical-mechanical response, which are observed experimentally.Binder (); Frenkel () It is unfortunate that, despite the progress already occurred in the computational hardware and software, the study of the dynamics of polymeric liquids by atomistic computer simulations is still limited by the impossibility of simulating polymer dynamics in the wide range of time and length scales relevant for these systems. It is known that long simulation trajectories deteriorate with time, following a law defined by the Lyapunov exponent characteristic of the system, and the resulting long-time dynamics is affected by errors.shadowing (); Tisley ()
Atomistic simulations of polymeric liquids are limited either in the length and number of macromolecules, or in the maximum timescale that can be reached. The long timescale is a regime of great interest for these systems, because of the relevance of dynamical-mechanical properties of macromolecules in that regime for their industrial applications. Because the limitation concerns only the region of long-time, large-scale properties, it is possible to shift the focus of a simulation to this region of interest by the reduction of the simulated degrees of freedom through the averaging of local scale properties, i.e. adopting a so-called “coarse-graining” procedure.
While simulations of coarse-grained systems afford to represent larger length- and longer time-scales than atomistic simulations, they also are marred by the shortcoming that the measured dynamics is unphysically fast. This is the advantage that allows one to perform simulations in the long-time regime, but it is also a problem because the dynamics in the simulation becomes too fast and measured diffusion coefficients are too large. Depending on the level of coarse-graining the dynamics can become orders of magnitude faster than the atomistic dynamics.ShortPaper (); LongPaper ()
Being aware that the dynamics has to be rescaled to recover the correct timescale,Kroger (); Kroger1 (); Voth () one resorts to the most conventional method, which is to build a calibration curve calculated by superimposing the long-time dynamics of the coarse-grained and atomistic simulations.Padding (); Klein () The hope is for this calibration curve to be transferable to other systems, or to be identical, and so applicable, for thermodynamic conditions close to, but different than the phase point for which the curve was optimized in the first place. Otherwise, if an atomistic simulation has to be ran every time we need to rescale a coarse-grained simulation there is no advantage in running simulations of coarse-grained (CG) systems. The rescaling factor of the conventional calibration curve is, however, purely a numerical correction. There is no physical motivation for one to assume that this correction has to be identical for any system or thermodynamic condition different from the ones in which it was optimized. Given that the effective potentials between coarse-grained units are free energies, and as such they are parameter dependent, it is very unlikely that numerically optimized corrections will be identical for different systems. However, if the calibration curve is applied for conditions of temperature, density, and degree of polymerization, very close in the phase diagram to the ones in which the calibration curve was originally optimized, there is hope that the error is very contained. Alternatively, in the case that very limited levels of coarse-graining are applied, e.g. combination of a few atoms to obtain a new coarse-grained unit, the error incurred in using the calibration curve in different conditions, can be small enough to be within the numerical error of the calculation. This is in fact the reason for the observed good agreement in the rescaled dynamics of systems with contained level of coarse-graining.polymer (); pb (); Vagelis ()
The goal of this paper is orthogonal to the more conventional calculation of an optimized calibration curve for systems with limited level of coarse-graining. In fact our goal is to investigate the physical motivation that leads to the accelerated dynamics in first place and to provide, through our analysis, a theoretical approach to calculate, from first principle equations, the correction factor that needs to be applied to the fast dynamics measured in coarse-grained simulations to recover the slower and more realistic dynamics of the atomistic description. We study the most extreme level of molecular coarse-graining, where one polymer chain is represented as a soft sphere, because i) the description is analytical allowing the formal study of the rescaling problem, ii) it provides a solid test of the reconstruction procedure because possible errors would be maximized, iii) this level of coarse-graining is associated with the largest computational gain, so it is important.
Recently we have proposed an original procedure to rescale the dynamics measured in simulations of coarse-grained polymer melts.ShortPaper (); LongPaper () The procedure we proposed does not require to run atomistic simulations to calibrate the dynamics of the coarse-grained simulation. Once one parameter, specifically the diameter of the hard-core monomer potential, is fixed, the rescaling factors are fully determined. The advantage of this method is clear, given that the coarse-grained simulation can reach larger-scale longer-time properties than the atomistic simulation. Because no atomistic simulations are needed, the measured dynamics can be directly rescaled to obtain the dynamics of the real system we wanted to study.
The method has been formally derived, and then applied to systems of polyethylene (PE) melts. The choice of PE as a test system was motivated by the wealth of experimental and simulation data available in the literature. We have shown that the proposed theory is able to predict diffusion coefficients in agreement with both atomistic simulations and experiments for systems in both the unentangled and entangled regimes, for a range of temperatures and densities. For entangled system the approach of dynamical reconstruction includes an extra loop in the calculation of the friction coefficient, which accounts for the fact that in real systems described at the atomistic level, the dynamics is slowed down by the presence of entanglements. We hypothesized for the entanglements to relax in time following the same dynamical mechanism of single chain interdiffusion. In this paper entangled chain dynamics is not investigated as most of the samples in the atomistic simulations are either unentangled or in the region of crossover from the unentangled to the entangled regimes.
One of the questions that still needed to be answered was if our method was able to generate high quality predictions also for macromolecular liquids with a more complex monomeric structure than the simple polyethylene. A first study of this issue is presented in this manuscript which investigates the agreement between theoretically predicted diffusion coefficients for cis-1,4-Polybutadiene (PB) chains and data from molecular dynamics simulations performed by Tsolou, Mavrantzas and coworkers.Tsolou (); Mavrantzas1 (); Mavrantzas2 () The approach is completely general and can be applied to polybutadiene chains with different tacticity. For these systems the difference in local chain semiflexibility, which corresponds to different overall chain dimension, i.e. a different radius-of-gyration () and statistical segment length (), ultimately leads to a different diffusion coefficients for polybutadienes with different tacticity. In this paper we limit our calculations to cis-1,4-Polybutadiene samples because those are the samples investigated by atomistic simulations and provide a effective test of our procedure for dynamical reconstruction.
Our rescaling method considers two contributions emerging from the consequence of applying the coarse-graining process and derives the theoretical corrections that need to be applied to the dynamics of the coarse-grained description to recover the modalities of the original atomistic system. To calculate the analytical correction we adopt for the atomistic descriptions a simple bead-and-spring description of the chain where each PB polymer of monomeric units is represented as a collection of beads connected by semi flexible springs of length , and the overall chain dimension is defined by the radius-of-gyration of the molecule, . This model is a Rouse approach modified to include chain semiflexibility and has been shown to represent well the diffusion of polymer melts in the unentangled regime.DoiEdw () From the comparison of this model and the soft sphere representation of the chain, it is possible to evaluate the two contributions that are responsible for the accelerated dynamics, and which have to be included to correct the measured dynamics and bring it to the realistic values of the atomistic (here bead-and-spring) description.
The first correction emerges from the consideration that due to the process of averaging the intramolecular degrees of freedom, the coarse-grained system experiences a change in its entropy, because a number of local states are neglected. While the atomistic system devotes time to explore local configurations, the coarse-grained system doesn’t. To recover the correct dynamics those states need to be included a posteriori in the form of an entropic contribution to the free energy that rescales the time measured in the CG simulations.
The second correction comes from the change of the “shape” of the molecule once it is coarse-grained. Macromolecules described at the atomistic level become represented by chains of effective atoms, or as bead-and-springs, chains of soft blobs, or even each molecule as a soft sphere. In all these multiple forms the surface of each molecule available to the surrounding molecules, i.e. the “solvent”, is different. The hydrodynamic radius, , of the coarse-grained unit in each description is different and so is its friction coefficient, , defined by Stokes’ law as , where is the viscosity of the medium. By solving the memory function for the friction coefficient in each description we are able to calculate the correction factor that has to be applied to each friction coefficient to recover the description at a different level of coarse-graining.
In this paper our rescaling method is first briefly summarized and then applied to study the dynamics of cis-1,4-PB and compared with simulation data by Tsolou et.al.Tsolou () who also investigated the effect of temperature and pressure on the simulations.Mavrantzas1 (); Mavrantzas2 () The global dynamics of the polymer is represented by the diffusion coefficient and the decay of the self-correlation for the end-to-end molecular vector. These dynamical quantities are reconstructed using the method by Lyubimov at al.,ShortPaper (); LongPaper () from the information contained in the trajectories of the mesocale simulations of the coarse-grained polymer melts and show good agreement with the data from united atom simulations. For the internal dynamics the theory for cooperative motionMarina (); MGuenza1 (); Richter () is used to predict monomer dynamics and the dynamical structure factors, which compare well with the simulations. If atomistic simulations are not available, the theory of dynamical reconstruction is able to predict the correct dynamics from the rescaling of the mesoscale simulations having as an input the polymer radius-of-gyration, which can be calculated from simple structural models for polymers. The accuracy of adopting a freely rotating chain model for the calculation of the polymer radius of gyration is discussed in the last section. A brief summary concludes the paper.
Ii Rescaling of the free-energy and rescaling of the friction
In this section we briefly present an overview of the main steps in our procedure. A complete and detailed presentation of the same has been already publishedShortPaper (); LongPaper () and will not be repeated here. We start from the consideration that the procedure of coarse-graining corresponds to eliminating some internal degrees of freedom, by combining groups of atoms into one effective CG unit. Specifically, our rescaling theory has been developed for an extended level of coarse-graining where intramolecular coordinates are fully averaged out and the polymer is represented as an isotropic sphere centered on the center-of-mass (com) of the macromolecule. There are some advantages in the choice of this description. First, the CG representation is fully isotropic, and even if it is known that the shape of a polymer is closer to an ellipsoid than to a sphere,DoiEdw () the total correlation function of the polymeric liquid, i.e. the structure of the polymeric liquid, statistically averaged over all the possible configurations is well reproduced by this model; second the formalism is analytical, which allows us to derive formally many of the physical quantities of interest for both the static and the dynamic properties of the coarse-grained system as a function of the atomistic description; third all the consequences of the coarse-graining procedure are enhanced with this extreme level of coarse-graining and are easier to study in such a description. If the procedure adopted were not precise, the error would be maximized because of the high level of coarse-gaining; the fact that the theory, which is predictive, is found to be in quantitative agreement with both simulated and experimental data is encouraging.ShortPaper (); LongPaper ()
One disadvantage of having such a coarse-grained representation of the molecule is that no information is collected from the mesoscale simulation on the internal dynamics of the polymer. In this way, the soft sphere representation only provides information on the diffusion coefficient of the molecule. However, we show in this paper that from the knowledge of the diffusion coefficient and the related monomer friction coefficient for unentangled chains it is possible to recover correctly the internal dynamics of the chain by applying our theory for the cooperative dynamics of a group of interacting chains, which calculates the single chain dynamics for a macromolecule whose dynamics is coupled by the presence of intermolecular forces with the other interpenetrating macromolecules in the liquid.
In every CG model, the averaging of the intramolecular degrees of freedom leads to a speeding up of the dynamics. When a polymer is represented as a soft sphere the dynamics in the mesoscopic simulation is orders of magnitude faster than in the atomistic description. One of the reasons for this acceleration is the fact that while a macromolecule needs to explore a large number of internal chain configurations for each position of the center-of-mass, the soft sphere instead is free to move in the three-dimensional space undergoing Brownian motion.entropy1 (); entropy2 () To take into account the missing contribution of the time necessary to explore the internal degrees of freedom, the time measured in the mesoscale simulation of the soft sphere is multiplied by a rescaling factor that is calculated as the ratio of the energy due to the internal degrees of freedom in the two representations, i.e an atomistic and a soft sphere representations. For the atomistic representation we adopt an analytical bead-and-spring model, as previously described.ShortPaper (); LongPaper () This model has been shown to provide a quantitative description of the dynamics for polymeric liquids and for proteins in theta solvent.Richter (); BJ ()
The first correction term is calculated from the definition of the internal energy for a simple bead-and-spring model and for the soft sphere CG model. The mesoscale (MS) molecular dynamic (MD) simulations are performed using a potential expressed in dimensionless units, which is a common procedure. The time from the simulations has to be properly normalized so that the dimensionless simulation time from MS simulations, , once dimensionalized and rescaled reads as
where is the molecular mass of one chain and is the radius-of-gyration of the molecule, which is an input quantity in our approach. The internal degree-of-freedom averaged out in the coarse-grained description are accounted for by introducing the factor where is the number of beads.
The second rescaling is due to the change in the shape of the molecule represented in the two levels of coarse-graining, and so the effective friction in the two representations. In the long time the com mean-square displacement is related to the molecular diffusion coefficient as
where is the diffusion coefficient of the system, which has to be calculated. If we use as a starting point the diffusion coefficient measured in the mesoscale simulation , the mean square displacement in the rescaled formalism is give by
where the problem reduces to the calculation of the friction coefficient in the monomer and soft-sphere representations. The friction coefficients are calculated from the solution of the memory function in the two representations (bead-and-spring and soft-sphere).
The formalism presented so far, adopts a Rouse-like description of the single chain diffusion for the atomistic level representation of the system. At the coarse-grained level the chain is represented as a soft-sphere. Eq.3 shows how the diffusion coefficient measured in a mesoscale simulation of a liquid of soft-spheres, , needs to be rescaled to give the mean-square displacement of the real chain, , as a function of the rescaled time, . The contributions that are still unknown are the friction coefficients in the atomistic and coarse-grained descriptions, which can be solved analytically starting from the definition of the memory function.Zwanzig ()
For a bead-and-spring model the monomer (bead) friction is given by:
where the memory function is defined asSchweizer ()
where , is the monomer density and , are the unit vectors characterizing the direction of the forces acting on monomers and , respectively. is the projected dynamic structure factor, which includes both intra- and inter-molecular contributions, i.e. incoherent and coherent scattering, and is approximated by its non-projected form as , which is an acceptable approximation when the Langevin equation is expressed as a function of the slow variables.LongPaper ()
In Eq.(5) is the monomer pair distribution function of the molecular liquid, and is the force due to the effective potential between two monomers, obtained from the solution of the Ornstein-Zernike (OZ) equation by applying the Percus Yevick closure. The monomer potential, which in the atomistic-level simulation is a Lennard-Jones potential, is approximated here by an effective hard-sphere with an effective diameter, , to mimic the properties of the L-J potential. All the physical quantities that appear in the equation are known and the analytical expression for the monomer friction is LongPaper ()
This equation contains quantities that are well defined once the system of interest is selected. For example is the pair distribution function at contact, is the density fluctuation correlation length with , and already defined in the text, and is defined as . The only physical quantity that needs to be determined is the hard-sphere diameter, . This is defined once the Lennard-Jones potential is mapped onto an effective repulsive hard-sphere system as described later in this section of the paper.
In the soft colloidal representation, the friction coefficient is analytically calculated from the definition of the memory function as
where is the chain density. All the other physical quantities that appear in the integral, for example the pair distribution function, , and its related force, , are now defined in relation to the description of the polymeric liquid as a liquid of soft spheres, which are point particles interacting through a soft repulsive potential. Each of these quantities is analytical and calculated from the solution of the OZ equation with the hypernetted chain closure approximation.YAPRL ()
The regime of interest in our calculations is the diffusive limit where, in reciprocal space, the wave vector . Here the dynamic structure factor can be approximated as
with the total correlation function of the soft-sphere representation is given in the limit of long chains () by the approximated expressionYAPRL ()
The resulting equation for soft particle friction is LongPaper ()
where the value of each physical parameter (temperature and density) and molecular parameters (degree of polymerization and radius-of-gyration) that enter this equation is defined once the system to simulate is selected. These parameters are identical to the ones that are input to the equation of the monomer friction, Eq.(6), given that the two equations are representations for different levels of coarse-graining of the same system.
When the two equations for the monomer and the soft colloidal friction coefficients are introduced in Eq.(3) the diffusion coefficient for a single chain in a liquid is recovered as
where we have assumed that the diffusion coefficients that appear in Eq.6 and Eq.10 are identical, as the long-time relaxation of the dynamic structure factor in both CG descriptions is guided by molecular diffusion.
To summarize, as the first step the thermodynamic and molecular structure, i.e. the radius of gyration or equivalently the persistence length, of a polymer melt has to be defined. Then, MS MD simulations are performed for a liquid of point particles interacting through a soft pair potential, , as described in several of our previous papers, where
The mean-square-displacement of the soft spheres is measured from the MS MD trajectories as a function of time (dimensionless), and the resulting diffusion coefficient is rescaled following Eq.11 to obtain the reconstructed diffusion coefficient.
For the solution of the rescaling factors in Eq.11 the only parameter that has to be defined is the value of the effective hard sphere diameter, , in Eq.6. In our formalism the elementary interaction between monomers belonging to different chains, which in the UA MD is a Lennard-Jones potential, is approximated by a hard-sphere interaction with an effective diameter, , to allow for the analytical solution of the monomer memory function. Under fixed thermodynamic conditions, the value of should depend only on the local monomer structure, and be independent of the degree of polymerization, because the monomer interaction potential is a local property of the chain. Therefore in our model, once the value of is chosen, this value is kept fixed for all samples with different molecular weights. Note that for convenience we use the term monomer to identify the group, where or 3.
To evaluate the value of we analyze the dimensionless quantity , which is identical to if the chain obeys strictly Rouse dynamics in the long-time limit. To fulfill this requirement the chain has to have much smaller than the entanglement degree-of-polymerization and larger than the value of for which chains start to obey Gaussian statistics, i.e. the central limit theorem applies.
Figure 1 a) displays the dimensionless friction coefficient, , from Eq.(6) as a function of the hard sphere diameter, , for PB samples of increasing degree of polymerization, . The parameter values (, , , ) are chosen to match those of UA MD simulations in ref. Tsolou () as reported in Table 1. Horizontal lines represent the value for each sample, which is the Rouse value of the physical quantity in the diffusive regime, and reproduces the scaling behavior of the diffusion for unentangled, short-chains. For comparison Figure 1 b) displays the analogous plot for the PE samples, which is reproduced from refs.ShortPaper () and LongPaper ().
The change of with observed in Figure 1 is a consequence of two different effects: for short chains, end effects are important, while for long chains the crossover to entangled dynamics starts to be important. In this way, it could be argued that in our method the actual fitting parameter is not , but the length of the chosen sample for which is determined. Clearly when entangled chains are forced to behave according to the Rouse expression the hard sphere diameter needs to be artificially decreased in order to compensate the overlooked increase in the monomer friction, which is associated with entanglement effects growing with increasing . Whereas the range of lengths in the unentangled region can be rather wide, it is always possible to choose the chain length knowing experimentally measured or theoretically estimated value of ,entengl_Fetters () and by optimizing our choice of 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 information is the value of and the value of for a sample of short, unentangled, chains ().LongPaper ()
In our previous study of PE melts we fixed the hard sphere diameter following the same procedure described here, and we selected a sample with to calculate . The obtained value of Å (see Figure 1b), was larger than the carbon-carbon bond length Å, but smaller than the Lennard-Jones parameter Å typically used in UA MD simulations, which seems to be a reasonable choice of the parameter. Compared to PE, PB chains are more flexible and so they are less entangled. The estimated entanglement degree-of-polymerization, , is calculated from the generalized formula for the plateau value in the shear relaxation modulus,entengl_Fetters () which gives backbone carbon atoms for PB, compared to for PE. Notice that the value of from shear measurements is different from the value obtained from the analysis of scattering experiments which gives, for example, for PE samples .
It is important to notice that even if there is some freedom in choosing the reference sample, the deviations of from the chosen value is not large. However, small deviations of the value of can lead to different values of the diffusion coefficient. Figure 2 displays a comparison of the predicted diffusion coefficient as a function of when different values are chosen for the hard-sphere parameter . It is possible to see that for values of obtained by enforcing Rouse dynamics for three different chain lengths that obey diffusive Rouse dynamics in the long-time regime, the predicted diffusion coefficients are all in reasonable agreement with the values of the UA MD simulation.
Notwithstanding the fact that for PB there is a wide range of degrees-of-polymerization, , for which the hard-sphere diameter, , could be optimized, the best agreement between theory and simulations is obtained for values of the parameter optimized for short chains, . This can be explained by considering that in PB samples the crossover region from unentangled to entangled dynamics is extended (see for example Figure 3, where deviations from the scaling occur already for ). Samples with approaching are not following the unentangled scaling behavior, which is necessary condition for the optimization of using Rouse-like diffusion.
Considering the broader range of unentangeled PB chains, we select for the calculation of the hard sphere diameter the PB sample with . The obtained value of Å (which is smaller than the value for PE of Å given the flexibility of PB) is consistent with the fact that in PB chains half of the carbon bonds are shorter double bond and Å .
Iii Theoretical predictions of center-of-mass diffusion
The solution of Eq.(3) gives the diffusion coefficient for the center-of-mass of a polymer described at the atomistic level, having as an input the diffusion coefficient calculated from the mesoscale simulation of the CG polymer liquid. The purpose of this procedure is to have ultimately mesoscale simulations that, once properly rescaled, can provide directly the diffusion coefficient, without need of running atomistic simulations. To test our procedure we first run the mesoscale simulations, then we rescale the diffusion coefficient, and finally we compare the predicted values of against atomistic simulations and/or experiments. Thermodynamic and molecular parameters entering our equations have to be consistent with the parameters of the system against which we compare our approach.
We study the dynamics of cis-1,4-Polybutadiene (PB) melts of increasing degree of polymerization and compare our results with the data from the simulations performed by Tsolou et al.Tsolou () The input parameters to our theory are displayed in Table 1, including the degree-of-polymerization, , and the monomer density . Table 1 also includes the mean-square-radius-of-gyration calculated from the united atom simulations . From this radius-of-gyration the value of the semiflexibility parameter, , is derived, as the value that corresponds to chains with the desired overall dimension, i.e. . In this way the parameter is not an independent parameter but is determined by the value. The parameter does not enter the equation, but is used in the calculations reported in the last section of this paper. It is important to notice that the values for could be taken from experimental data, and that, in principle, no atomistic simulations of the system under study are necessary for our approach.
The UA MD simulations were performed in the ensemble, at the constant pressure, , and temperature, , reported in the table.Tsolou () Note that the values of for systems PB160, PB200, PB240 were corrected from 257Å to 275Å, from 304Å to 340Å and from 401Å to 410Å correspondingly, after consultation with the authors of ref.Tsolou (). The correction was motivated by the inconsistency observed with the data reported in reference Tsolou () when the calculation of is performed using a standard model, such as the Freely-Rotating-Chain, described later in in this paper.
|T=413K, P=1 atm|
In the mesoscale simulations of a polymer melt, each chain is represented as a point particle interacting through a soft-core potential derived from the solution of the Ornstein Zernike equation applying the Hyper Netted Chain closure. MS MD simulations were implemented in the microcanonical ensemble on a cubic box with periodic boundary conditions. We used reduced units such that all the units of length were scaled by () and energies were scaled by . More details of our simulation procedure have been reported in previous papers.ShortPaper (); jaypaper (); multis ()
Table 2 reports the diffusion coefficient directly calculated from the MS MD in the soft sphere representation, , and the diffusion coefficient calculated using the dynamical reconstruction procedure, . Finally, for comparison, the Table shows the diffusion coefficient measured in the united atom simulation, . The diffusion 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 directly from the atomistic simulation. It is important to notice that, once the parameter , which is characteristic of the polymer considered, is determined the diffusion coefficient is calculated without any input from the dynamics of the atomistic simulations, so the procedure is predictive.
The samples here are relatively short and these calculations are performed for unentangled and slightly entangled systems. For strongly entangled systems, we proposed in a previous paper a perturbative version of our approach that accounts for the fact that both the tagged chain and the surrounding chains, that provide the entanglements, relax following the same dynamics.LongPaper () In this manuscript the PB samples are in the unentangled and slightly entangled regimes, which are well represented by the theory without the one-loop perturbation.
On Figure 3 the diffusion coefficient is presented as a function of degree of polymerization . Filled symbols represent UA MD simulations from ref. Tsolou () and open symbols represent MS MD simulations rescaled as in Eq.(3). In analogy with the figure in Ref.Tsolou () from which the UA MD data are taken, we report three scaling regimes in terms of power law dependence of .For , there is a faster than Rouse regime with with , which is attributed to the free-volume effects due to chain ends which is significant only for very short chains FreeVolume1 (); FreeVolume2 (); FreeVolume3 (). In the intermediate regime for a Rouse-like behavior can be observed where the scaling exponent is close to 1. For longer chains with the value of is typical of the crossover to reptation dynamics.
The data in this plot are not closely following the reported scaling exponents because as the degree of polymerization increases the density of the system changes (Fig.4). The different regimes are even less prominent in the MS MD simulations, where we observe a smoother transition between regimes than in UA MD data. In our calculations features like free volume effects due to chain ends (causing faster than Rouse decay of the diffusion coefficient for short chains) enter only indirectly as a consequence of the interplay of the input parameters, i.e. , , . The increasing of the density with enters our analytical equation both directly, through the liquid density, and indirectly through the effective radius of gyration of the polymer, given that the compactness of the chain changes with the density (see also Section V of this paper).
The cis-1,4-PB chains are more flexible in comparison to PE chains, therefore we should expect larger entanglement length (for PE from scattering ). Figure 3 shows that crossover from Rouse to reptation regime starts at . The largest system shown on Figure 3 is for which is still a very weakly entangled system (1 or 2 entanglements per chain). Overall the agreement between simulated and reconstructed diffusion coefficients is quite good.
Iv Internal dynamics
The coarse-grained model adopted in this paper simplifies the macromolecular structure by reducing the description of the molecule to an effective site that is centered 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 long-time diffusive behavior of the center-of-mass of the molecule. However from the knowledge of the diffusion coefficient the monomer friction coefficient can be obtained, and the latter can be used as an input to the Langevin equation that describes the internal dynamics of the polymer chain at any length scale of interest.
Specifically we follow our approach for the Cooperative Dynamics of a group of macromolecules, or Cooperative Dynamics Generalized Langevin Equation (CDGLE).Manychain () Conventional theories of polymer dynamics, such as the Rouse theory for unentangled polymer dynamics and the “reptation” model for entangled dynamics, predict diffusive motion at short time following the ballistic regime; this is however in disagreement with simulations and experiments, which show a sub-diffusive regime before the com starts to follow Brownian motion. Both theories are mean-field approaches of single chain dynamics in an uniform bath.DoiEdw () Our approach focuses on the dynamics of a group of polymer chains in a melt and relates the anomalous subdiffusive behavior to the presence of cooperative motion of a group of polymer chains in a dynamically heterogeneous liquid, as observed in simulations data of polymer melts and experiments.het1 (); het2 (); het3 (); het4 () Theoretical predictions of the com and monomer mean-square displacements are shown to be in quantitative agreement with computer simulations of unentangledMGuenza1 (); MGuenza2 () and slightly entangled polymer melts,pallavi () and with scattering experiments of Neutron Spin Echo.Richter ()
The physical picture underlying the theory builds on the fact that in polymer melts the dynamics appears heterogeneous with the motion of a tagged chain correlated to the dynamics of a group of chains comprised inside the range of the intermolecular potential. The latter has a range of the order of the correlation hole, i.e. of the radius-of-gyration of the polymer, as it emerges from the solution of the OZ equation. The number of chains comprised in a volume of the order of is given by , where , the statistical segment length, and the other quantities have been defined earlier on. The number of interpenetrating chains, , increases with increasing density, degree of polymerization, and the stiffness of the polymer.
In the Cooperative Dynamic approach the dynamics is described by a set of coupled equations of motion (eom). Each equation is expressed in the space coordinates of the monomer , belonging to molecule and positioned at , and contains a balance of different forces acting on the monomer. These include the viscous force, , the intramolecular force ,which contains the structural matrix , the time-dependent intermolecular mean-field force , and the random interactions with the surrounding liquid, given by the random force .
Here is the matrix of intramolecular connectivity, which reduces to the Rouse matrix when infinitely long and completely flexible macromolecules are considered, and is defined as
Here is the equilibrium averaged bond correlation matrix
and is a structural matrix, with all the elements equal zero except for , with , and with . The matrix is a function of the local semiflexibility parameter , which is related to the persistence length of the polymer. For our samples the values of are reported in Table 1 and are calculated from the molecule radius of gyration, which can be obtained from experiments or from atomistic simulations.LongPaper ()
Through Eq.(16) the eom includes a complete microscopic description of the structure and local flexibility 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 undergoing slow cooperative dynamics are coupled by the presence of intra- and inter-molecular interaction potentials. The intrinsic, chemical dependent semiflexibility of the macromolecule enters explicitly through the description of the matrix.
The intermolecular potential is time dependent, as it is a function of the relative position coordinates of the centers-of-mass of a pair of molecules. As the two molecules move with respect to each other the intermolecular force changes. From the initial ensemble of dynamically correlated chains, molecules diffuse in time and finally move outside the range of the potential, , in a timescale of .
The set of coupled equations is solved after transformation into normal modes of motion and numerically using a self-consistent procedure that calculate the effective potential at any given distance. Once and are defined, the first optimized against the simulation data and the second from our rescaling procedure, the set of equations has no adjustable parameters.
By applying consideration of symmetry the set of coupled eoms in normal coordinates can be simplified and reduced to two sets of uncoupled equations in the relative and collective normal-mode coordinates . We assume that the eoms for internal modes () do not contain intermolecular forces. This approximation is justified on the basis that polymer local dynamics is affected in a similar way by inter- and intra-molecular excluded volume interactions, which in polymer liquids tend to compensate each other. Intermolecular forces, which enter the dynamics through the eom for the first () normal mode, still perturb the dynamics on the local scale through the linear combination of the normal modes.
In the following, we present an overview of some of the quantities analyzed in the original simulationsTsolou () and their comparison with the predictions of the Cooperative Dynamics theory having as an input parameter the reconstructed friction coefficient calculated from the soft-sphere simulations and the rescaling procedure as described in the first part of this paper.
Figure 5 presents the end-to-end vector, time decorrelation function for the sample of PB. The Cooperative Dynamic theory with the reconstructed monomer friction coefficient is directly compared against UA MD data. This function represents the rotational relaxation of the chain and cannot be measured from the soft colloid simulation directly. The sample has chains shorter than the entanglement length and the Cooperative Dynamics theory reproduces the atomistic simulation data rather well.
Figure 6 shows the center of mass mean square displacement (com MSD) for three samples of cis-1,4-Polybutadiene melts with chains of increasing degree of polymerization, and . The long time diffusive dynamics is calculated from the diffusion coefficients obtained from the rescaling procedure. For time shorter than the longest relaxation time, , the mean-square displacement of the com shows a subdiffusive behavior that cannot be reproduced by the Rouse approach, i.e. the conventional theory of unentangled chain dynamics. In the cooperative dynamics theory the subdiffusive dynamics is a consequence of the cooperative motion of the interpenetrating chains inside the correlation hole region, coupled by the effective intermolecular potential between the coms of the macromolecules. The theory, with the monomer friction from rescaled MS MD simulations as an input, is compared against the UA MD simulation data. The diffusion coefficient for PB240 is underestimated and for PB320 is slightly overestimated by the theory in comparison with the simulations, which is also can be seen in Table 2 and Figure 3. These deviations are within the error related to . Dot-dash lines in the inset of Figure 6 represent results calculated with the upper and the lower values of taken from the data of the atomistic simulations, as reported in Table 1 for the PB240 sample. The number of correlated chains is optimized to reproduce the subdiffusive behavior giving for PB240, PB320 and PB400 respectively, with a trend that qualitatively agree with the predicted scaling behavior. The error in the radius-of-gyration does not allow for a more precise calculation of the parameter .
Starting from the theory of cooperative dynamics, the single chain dynamic structure factor can be calculated as
Eq.17 describes the dynamics of a ”tagged” polymer chain, which moves in a group of interpenetrating chains. The intermolecular interaction between chains affects the dynamics resulting in the subdiffusive motion of the single-chain center-of-mass, as depicted in Figure 6, and in the subdiffusive decay of the low dynamics of .
The average time-dependent distance between two monomers in the tagged chain depends on the presence of the other chains dynamically correlated, and it is expressed as a function of relative and collective mode coordinates. Once and are determined, the theory predicts the full decay of the dynamic structure factor, at the different experimental wave vectors. Figures 7 and 8 display the normalized dynamic structure factor obtained from the theory and compared with UA MD simulations for Å. For unentangled system PB96 shown on Figure 7 the agreement is excellent for all values of scattering vector and time range. The Rouse model (not shown in this figure, see ref. Tsolou ()) fails to describe these data with satisfactory agreement, while the cooperative dynamics theory is found to be in quantitative agreement for samples at different molecular weights, over the whole time range, and for different wave vectors.
On Figure 8 the normalized dynamic structure factor is shown for the weakly entangled system PB400. The agreement for the smallest scattering vector which represents global dynamics is excellent at all times, while for intermediate and large values of the agreement between theory and simulations is slightly less accurate. The theoretical decay, which is slightly faster than simulations at long times, suffers from the lack of properly describing the entanglement effects. A more recent version of the cooperative dynamics theory includes chain uncrossability, i.e. entanglements, and shows the correct slowing down of the relaxation in the long-time regime for entangled samples.marinainprogress ()
Figure 9 displays the monomer mean square displacement for cis-1,4-Polybutadiene sample with . Theoretical prediction is compared against UA MD data from ref. Tsolou (). The agreement in both subdiffusive and diffusive regions is very good. The theory accounts for the local semiflexibility and chain connectivity, together with the cooperative motion of the chains due to the presence of intermolecular interactions.
V Theoretical predictions using the freely rotating chain model
The purpose of the rescaling procedure is to predict diffusion coefficients from the mesocale simulations properly rescaled when the atomistic simulations or the experiments are not available. One key input quantity to the procedure is the molecular radius of gyration. In the previous sections we adopted the values of determined from the atomistic simulations. When there are no data available for the radius of gyration, neither from atomistic simulations nor from experiments, it is convenient to use some simple model, like the freely rotating chain model to provide , as it captures the semiflexibility of the polymer chain. In this section we present a critical discussion of adopting a freely rotating chain model to describe the PB chain, using a constant flexibility parameter .
For the cis-1,4-polybutadiene chain with carbon atoms along its backbone there are double bonds with Å, single bonds with Å and single bonds with Å. The average bond length over all bonds is calculated as
which in the limit gives the value of Å.
The radius of gyration for a freely rotating chain (FRC) model with semiflexibility parameter is given as
Figure 10 shows the density dependence of the semiflexibility parameter in the UA MD simulations.
Fitting simultaneously all samples to data gives an average value of the semiflexibility parameter for all densities, . Figure 11 displays the radius of gyration squared over degree of polymerization as a function of from united atom simulations and from the freely-rotating-chain model. The curvature at small is due to chain-end effects. Samples with small present the largest deviation between the two sets of data, while the constant flexibility hypothesis works best for the large samples. The deviation is due to the fact that the united atom simulations are performed in the NPT ensemble where increasing corresponds to an increase of the liquid density, as shown in Figure 4. The hypothesis of a constant semiflexibility parameter works well for the long chains for which coarse-graining methods are most useful.
With this in mind, we report in Figure 12 the diffusion coefficient calculated from the mesoscale simulations, properly rescaled, where now the input radius-of-gyration is taken from the freely rotating chain model. As it is expected the best agreement between UA simulations and theoretical predicted diffusion is for the samples with larger . The disagreement at small is related to the fact that the density is changing with increasing in that region of the plot, and the radius of gyration should be properly corrected for this difference. In conclusion, if an exact value of the radius of gyration as a function of thermodynamic conditions and degree of polymerization is known, the described procedure should be able to provide an accurate estimate of the dynamics starting from the mesoscale simulations of the coarse-grained system, which are computationally very efficient. The figure also reports a number of samples for which neither experiments nor simulations have been performed (, , , and ). The most relevant conclusion of this part of our study is the fact that it is possible to make predictions for data that are not reported in the UA MD simulations, obtaining information for new systems for which neither simulations nor experiments are available.
In a couple of recent papers we have presented an original approach to rescale the dynamics measured in mesoscale simulations of coarse-grained systems and recover the realistic dynamics as measured in atomistic simulations performed in identical thermodynamic conditions. The relevance of this method is in its ability of predicting diffusion coefficient, with good precision, directly from the simulation of a coarse-grained system, without the need of performing an atomistic simulation. Simply put, the method does not need a calibration curve. While the previous papers focused on the coarse-grained dynamics of polyethylene samples, this is the first paper in which the method is applied to a polymer with a different monomeric structure.
In this paper we use an “extreme” level of coarse-graining as the whole chain is represented as a soft colloidal particle. The choice of this representation is motivated by two reasons. The first is that a possible error in the procedure would be clearly made evident in this level of coarse-graining, and we believe that the procedure would produce an even better agreement if the macromolecule is coarse-grained at a lower level, e.g. as a collection of soft particles, or blobs.anthony (); anthony1 () The second is that for this representation all the physical quantities that enter our approach are analytically solved and the rescaling procedure is formally derived.
The paper presents a comparison between theoretical predictions of coarse-grained simulations properly rescaled following our first-principles procedure, and data from united-atom simulations for cis-1,4-Polybutadiene melts. The simulations were performed by Tsolou et al. in a ensemble.Tsolou () The rescaled diffusion coefficient shows good agreement with the diffusion measured in atomistic simulations, which supports the validity of the proposed procedure.
As the coarse-graining model represents the whole polymeric chain as a soft colloidal particle, the only dynamical information that can be collected from the rescaled mesoscale simulation is the diffusion coefficient of the center-of-mass of the polymer. From this, however, the friction coefficient of the monomer is derived. This paper shows how the monomer friction coefficient is used as an input to traditional Langevin equations for the dynamics of polymeric liquids. Specifically we use here the cooperative dynamics model for unentangled chains, given that the samples in the UA simulations are in the regime of unentangled and slightly entangled dynamics. From the theory of cooperative dynamics, which represents the cooperative motion of a group of interacting polymers in a dynamically heterogeneous liquid, all time correlation functions of interest can be calculated. Specifically we present here the comparison with the quantities reported in the original paper by Tsolou et al.Tsolou () and show that the theory, with the rescaled friction coefficient, is able to reproduce quantitatively the dynamics observed in the end-to-end time decorrelation function, and in the dynamic structure factor for unentangled and slightly entangled chains, and from the global to the local intramolecular dynamics. For entangled chains the local dynamics, represented by the high regime in the dynamic structure factor and in the monomer mean-square-displacement, shows a slight departure from the theory as the cooperative dynamics approach adopted here does not yet include the effect of entanglements.marinainprogress ()
We acknowledge support from the National Science Foundation through grant DMR-0804145. Computational resources were provided by Trestles through the XSEDE project supported by NSF. We are grateful to V. G. Mavrantzas for the informative discussion of his paper.
- (1) K. Binder Ed. Monte Carlo and Molecular Dynamic Simulations in Polymer Science (Oxford University Press, New York, 1995).
- (2) D. Frenkel, and B. Smit, Understanding Molecular Simulation (Academic, New York, 2000).
- (3) T. Sauer, C. Grebogi, and J. A. Yorke Phys. Rev. Lett. 79, 59 (1997).
- (4) M. P. Allen, and D. J. Tildesley Computer Simulation of Liquids (Oxford Science Publications, Oxford, 1992).
- (5) I. Y. Lyubimov, J. McCarty, A. Clark, and M. G. Guenza, J. Chem. Phys. 132, 224903 (2010).
- (6) I. Y. Lyubimov, M.G. Guenza Phys. Rev. E 84, 031801 (2011).
- (7) P. Ilg, H. C. Öttinger, and M. Kröger Phys. Rev. E 79, 011802 (2009).
- (8) P. Ilg,and M. Kröger J. Rheol. 55, 69 (2011).
- (9) Izvekov and Voth, J. Chem. Phys. 125, 151101 (2006).
- (10) J. T. Padding, W. J. Briels J. Phys. Cond. Matt. 23 233201 (2011).
- (11) S. O. Nielsen, C. F. Lopez, G. Srinivas, M. L. Klein J. Phys. Condens. Mat. 16, R481 (2004).
- (12) Y. Li, S. Tang, B. C. Abberton, M. KrÃ¶ger, C. Burkhart, B. Jiang, G. J. Papakonstantopoulos, M. Poldneff, W. K. Liu Polymer 53, 5935 (2012).
- (13) J. T. Padding and W. J. Briels J. Chem. Phys. 117, 925 (2002).
- (14) V. A. Harmandaris and K. Kremer Soft Matter 5, 3920 (2009).
- (15) G. Tsolou, V. G. Mavrantzas, D. Theodorou Macromol. 38, 1478 (2005).
- (16) G. Tsolou, V.A. Harmandaris, V.G. Mavrantzas Macromol. Theory Simul. 15, 381 (2006).
- (17) P.S. Stephanou, C. Baig, G. Tsolou, V.G. Mavrantzas, and M. Kroger J.Chem.Phys 132, 124904 (2010).
- (18) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics (Clarendon Press, Oxford, 1986).
- (19) M. G. Guenza J. Phys: Cond. Matt. 20, 033101 (2008).
- (20) Marina Guenza, Phys. Rev. Lett. 88, 025901 (2002).
- (21) M. Zamponi, A. Wischnewski, M. Monkenbusch, L. Willner, D. Richter, P. Falus, B. Farago and M. G. Guenza J. Phys. Chem. 112, 16220 (2008).
- (22) E. Caballero-Manrique, J. K. Brey, W. A. Deutschman, F. W. Dahlquist, and M. G. Guenza, Bioph. J. 93, 4128 (2007).
- (23) Y. Rosenfeld Phys. Rev. A 15, 2545 (1977).
- (24) J. A. Armstrong, C. Chakravarty, P. Ballone J. Chem. Phys. 136, 124503 (2012).
- (25) R. Zwanzig Nonequilibrium Statistical Mechanics (Oxford University Press, New York, 2001).
- (26) K. S. Schweizer, J. Chem. Phys. 91, 5802 (1989).
- (27) G. Yatsenko, E. J. Sambriski, M. A. Nemirovskaya, and M. Guenza, Phys. Rev. Lett. 93, 257803 (2004).
- (28) L.J. Fetters, D.J. Lohse, C.A. Garcia-Franco, P. Brant, and D. Richter Macromol. 35, 10096 (2002).
- (29) J. McCarty, I. Y. Lyubimov, M. G. Guenza Macromol. 43, 3964 (2010).
- (30) J. McCarty, I. Y. Lyubimov, and M. G. Guenza J. Phys. Chem. B 113, 11876 (2009).
- (31) F. Bueche, Physical Properties of Polymers (Interscience, New York, 1962).
- (32) V.A. Harmandaris, M. Doxastakis, V.G. Mavrantzas, D.N. Theodorou J.Chem. Phys. 116, 436 (2002).
- (33) E. von Meerwall, S. Beckman, J. Jang, and W.L. Mattice J.Chem. Phys. 108, 4299 (1998).
- (34) M. G. Guenza, J. Chem. Phys. 110,7574 (1999).
- (35) B. Frick, G. Dasseh, A. Cailliaux, and C. Alba-Simionesco Chem. Phys. 292, 311 (2003).
- (36) A. Arbe, J. Colmenero, B. Farago, M. Monkenbusch, U. Buchenau, and D. Richter Chem. Phys. 292, 295 (2003).
- (37) J. Colmenero, F. Alvarez, and A. Arbe, Phys. Rev. E 65, 041804 (2002).
- (38) G. Tsolou, V. A. Harmandaris, and V. Mavrantzas J. Non-Newtonian Fluid Mech. 152, 184 (2008).
- (39) Marina Guenza, Macromol. 35, 2714 (2002).
- (40) P. Debnath, and M. G. Guenza, Phil. Mag. 88, 33 (2008).
- (41) M. G. Guenza (in preparation).
- (42) A. J. Clark, and M. G. Guenza, J. Chem. Phys. 132, 044902 (2010).
- (43) A. J. Clark, J. McCarty, I. Y. Lyubimov, M. G. Guenza Phys. Rev. Lett. 109, 168301 (2012).