Analytical Rescaling of Polymer Dynamics from Mesoscale Simulations

The development of a systematic approach to bridge time scales between different hierarchical levels of description is an important goal in many areas of the physics of complex systems.GellMann Furthermore, understanding the dynamics of polymeric liquids is relevant for many technological applications. Polymeric materials are processed in their liquid state, and the custom tailoring of new materials requires detailed predictions of their mechanical and dynamical properties to be made based on their chemical structure. To this end, molecular dynamics (MD) simulations shed light on the properties of complex systems; however, MD is limited by the precision of the calculations, which degrades with the number of computer iterations.shadowing Quantitative predictions of the dynamics of unentangled polymer liquids are possible through simulation runs that are computationally demanding because the polymer diffusion coefficient, $D$, scales with the degree of polymerization, $N$, as $D\propto N^{-1}$. Even more demanding are simulations of liquids of long, entangled, polymeric chains, where the diffusion coefficient scales as $D\propto N^{-2}$. For entangled polymer liquids it is therefore difficult to obtain well-equilibrated samples or to simulate the system for several relaxation cycles, which would improve the precision of the calculated time-correlation functions.

To reach the long time and length scales of interest, special strategies need to be employed. For example, some gain in computational time has been achieved by speeding up the equilibration process through an end-bridging Monte Carlo algorithm.Mavran Moreover, if only qualitative, and not quantitative, predictions are sought, it is possible to adopt simplified intra- and intermolecular potentials, which can speed up each simulation step,Hess ; Likhtman or purely phenomenological coarse-grained descriptions, which reproduce the expected dynamical scaling behavior.Likhtman ; entangleds However, quantitative computational methods that precisely relate chemical structure to the long-time properties of the polymeric liquid are still lacking.

There is a need for computationally efficient, predictive methods to simulate polymer liquids and complex fluids in the long-time regime. The strategy we are pursuing here is to develop coarse-graining methods and dynamical rescaling procedures which start from first principles theory.Ottinger ; YAPRL The goal is to obtain quantitative predictions of real polymer dynamics directly from properly rescaled fast mesoscale (MS) simulations of coarse-grained liquids.

In a coarse-grained description, the system is specified by a set of relevant mesoscopic variables while smaller length scale variables are omitted.MSMD This is done at the expense of entropy and dissipation (friction), which are underestimated by coarse-graining.Ottinger Because of the reduced molecular degrees of freedom and the simplified energy landscape, MS-MD simulations require less computational time than atomistic simulations of the same systems.YAPRL ; DePablo ; Klein ; M-Plathe However, because the effective energy landscape of the coarse-grained representation is artificially smooth, MS-MD simulations predict accelerated dynamics, which need to be rescaled to produce realistic values.

The enhanced diffusion in coarse-grained systems arises from the “soft” nature of the intermolecular potential. This is advantageous in enabling larger time steps to be used in integrating the equations of motion, and an efficient sampling of the energy landscape, which leads to good statistical averages of the structural properties on the large scale.Nielsen However, the measured dynamics of coarse-grained systems is too fast, and as of yet, there is no clear procedure to quantitatively re-scale these accelerated dynamics.

In an effort to develop quantitative rescaling methods, it is custom to build a numerical “calibration curve” obtained from direct comparison of MS-MD time correlation functions with atomistic simulations. The “calibration curves” are parametric, normal mode dependent, and specific to the system against which they are optimized, as well as to its thermodynamic conditions. Building these parametric curves in part defeats the purpose of the coarse-graining procedure as it requires running many atomistic simulations to optimize the fitting parameters. If the level of coarse-graining is low, i.e. if the average is performed over a small number of atoms, the correction to the dynamics is only perturbative and the parametric rescaling works well.Kremer This explains the success of united-atom (UA) simulations.MONDE ; Mavran However, because the gain in computational time increases with the level of coarse-graining, simulations of slightly coarse-grained systems,e.g. UA-MD, afford a limited gain in time and length scales.

The most gain in computational efficiency is achieved through large-scale coarse-graining, as the one adopted in this paper. This is most useful when studying bulk physical quantities, e.g. viscosity, or systems with large-scale fluctuations, e.g. approaching spinodal decomposition. macromol Once simulations of heavily coarse-grained systems are combined with local scale simulations in a multiscale procedure, they provide the complete description of the system at all lengthscales of interest.multis

This paper presents a derivation of a first-principles approach to rescale the dynamics from MS-MD simulations of a coarse-grained polymer melt to the values of an atomistic description. The rescaling is analytical, general and transferible. The favorable comparison with atomistic simulations and experimental data in different thermodynamic conditions supports the validity of the proposed procedure. The theoretical basis of the rescaling rests on the fact that the accelerated dynamics is a result of the missing dissipation due to the eliminated degrees of freedom.

The paper is organized as following. In Section 2 we briefly review our coarse-grained model, while the coarse-grained potential and mesoscale simulations are described in Section 3. In the following section we present our atomistic model and the calculation of the energy due to the internal degrees of freedom. Rescaling of the friction coefficient is discussed in Section 5, followed by a comparison of the diffusion coefficients predicted from MS-MD simulations after rescaling, with united-atom simulation and experimental data. A brief discussion concludes the paper.

Our coarse-grained representation models polymers as interacting soft-colloidal particles with repulsive interaction of the order of the size of the macromolecule, defined by the radius-of-gyration, $R_{g}$.YAPRL The total distribution function is derived from the solution of a generalized Ornstein-Zernike equation, treating coarse-grained sites as auxiliary sites, and atomic sites as real sitesK2002 . In reciprocal space, the total distribution function is

The structure of the coarse-grained polymer liquid is described by the total distribution function,YAPRL approximated for polymer chains with $N\geq 30$ as

Each soft-colloidal particle interacts with other colloids through an effective potential of the range of the overall polymer dimension, $R_{g}$. While hard-sphere systems are best described by a Percus-Yevick closure, the Hyper-Netted Chain (HNC) closure works best for systems with soft potentials,mcquarrie including the mesoscopically coarse-grained polymer melts investigated here,YAPRL ; edjcp and polymer coils in dilute or semidilute solutions.Louis The HNC potential between a pair of coarse-grained units is derived by applying the closure $\beta v(r)=h(r)-ln[h(r)+1]-c(r)$, where the direct correlation function is defined by the Ornstein-Zernike relation in reciprocal space $c(k)=h(k)/(1+\rho_{s}h(k))$.SMPLQ

The potential between two spheres is calculated numerically, after Fourier transform, from the analytical expression, Eq.(1), by adopting the Debye form of the monomer intramolecular distribution, $\omega^{mm}(k)$. The Debye approximation has been shown to better represent simulation data than its Padé approximant.edjcp Tabulated HNC potentials are input to the mesoscale simulations.

MS-MD simulations of polymer liquids are performed in the microcanonical ensemble, where each molecule is represented as an interacting soft-colloidal particle. In the initialization step, all particles are placed on a lattice with periodic boundary conditions. Each site is given an initial velocity and subsequently, the system is evolved using a velocity Verlet integrator. Equilibrium is induced in the ensemble by rescaling the velocity at regular intervals until the desired average temperature is reached. At this stage, velocity rescaling is discontinued and trajectories are collected over a traversal of $\sim 8R_{g}$, while the temperature is monitored to assure that it fluctuates around the desired equilibrium value. Because of the form of Eq.(2), the simulation uses reduced quantities of distance, $R_{g}=1$, mass, $m=1$, and energy, $k_{B}T=1$.

A typical MS-MD simulation for our model is performed on a single-CPU workstation, and consists of $\sim 3000$ polymers, evolving for a duration of $\sim 4$ hours. The MS-MD simulation provides identical structural information to that of the analogous atomistic or united-atom (UA) simulations on length scales equal and larger than the polymer $R_{g}$. However, the MS-MD requires a much smaller computational power than the atomistic and UA-MD simulation, which is typically performed on a liquid of $\sim 500$ polymers, on a parallel supercomputer. The convenient requirements in computational power of the MS-MD simulation allows one to increase considerably the number of particles and the simulation box size without dramatically affecting the computational time, thus improving the precision on the large-scale data collected.

The correction in Helmhotz free energy, which accounts for the discarded internal degrees of freedom, is calculated starting from the atomistic representation of the liquid, where each chain is described as a collection of beads connected by springs defined by an effective intramolecular quadratic potential $U(r)=3k_{B}T/(2l^{2})\sum_{i,j=1}^{N}\textbf{A}_{i,j}\textbf{r}_{i}\cdot\textbf{r}_{j}$. Here $\bf{A}$ is the connectivity matrix, which represents the structure and local flexibility of the polymer, $\textbf{r}_{i}$ the position of unit $i$ in a chain of $N$ beads, and $\textbf{l}_{i}=\textbf{r}_{i+1}-\textbf{r}_{i}$ the bond vector connecting two adjacent beads.jpcmreview For polyethylene melts, each polymer is represented as a freely-rotating-chain (FRC), finite in size, with semiflexibility parameter $g=\left<(\textbf{l}_{i}\cdot\textbf{l}_{i+1})/(l_{i}l_{i+1})\right>=0.785$ and fixed bond length $l=1.54$Å. This bead-and-spring model of the FRC has been shown to represent correctly the dynamics of polyethylene melts as measured in united-atom simulationsmanychain and in experiments.Richter A FRC model was also successfully adopted to model the dynamics of polymer melts with different molecular architectures,jpcm and even proteins,jpcmreview ; protein once the proper semiflexibility parameters are selected.

To map the MS-MD simulation onto a real system, the reduced unit of time needs to be properly rescaled. Because energy is dissipated in internal degrees of freedom in the atomistic representation, the contribution due to the internal vibrational modes is included in the coarse-grained representation by rescaling the time in the MS-MD simulation, $\tilde{t}$, by the amount of internal free energy dispersed in vibrational modes in the atomistic description as $t=\tilde{t}R_{g}\sqrt{3mN/(2k_{B}T)}$, with the particle mass, $m$, and size $R_{g}$. By rescaling the internal free energy, we accounts for the change in entropy in the coarse-grained description. In the next section we derive the friction rescaling.

To account for the change in dissipation caused by coarse-graining, we start from the diffusion coefficient measured in the MS-MD simulation, $D^{MS}_{t}$, and we derive the center-of-mass (cm) diffusion coefficient of the polymer, $D_{cm}$, through the rescaling of the friction as

The correction factor is calculated from the ratio of the cm friction in the soft colloid representation $\zeta_{s}$ and the cm friction in the atomistic representation $N\zeta_{m}$, with $\zeta_{m}$ the monomer friction. Each friction coefficient is evaluated by solving the memory function in the corresponding Generalized Langevin Equation. These memory function definitions result from the straightforward application of Mori-Zwanzig projection operators to the Liouville equation, where either the monomer (atomistic description) or the center-of-mass (soft colloid description) of the polymer are assumed to be the “relevant slow variables.”jpcmreview ; SMPLQ

In the soft colloid representation the friction coefficient is defined as

In Eq.(4), the projected dynamics has been substituted with the real (unprojected) dynamics, which is a valid approximation when the Langevin equation is expressed as a function of slow variables for the diffusive regime.SMPLQ In the long-time regime, the relaxation of the liquid is dominated by the polymer center-of-mass diffusion, $D$, here represented by the cm of the colloidal particle. In Fourier space the dynamic structure factor of the liquid reads $S(k;t)\approx S(k)exp({-k^{2}Dt})$. Evaluating the integral with use of Eq.(2) gives

Eq.(5) expresses the friction coefficient of a soft colloidal particle.

In the atomistic description, the monomer friction coefficient is defined by the memory function as

To evaluate Eq.(6), we approximate the potential as an effective hard-core potential with a diameter $d$ to be defined.SMPLQ In hard-core fluids, $g(r)F(r)=g(d)\beta^{-1}\delta(r-d)$. Working in reciprocal space, the integrals in Eq. (6) can be performed analytically to give an expression for the dynamical quantity $D\beta\zeta_{m}$ depending on two length scales: $R_{g}$, and $d$. The result is lengthy, and it is not reported here.note The values of $R_{g}$ used are the ones reported in Tables 1 and 2.

The monomer hard-core diameter, $d=2.1$Å, is identical for all the samples, and is obtained by reproducing the scaling with $N$ of an unentangled melt, $D\beta\zeta_{m}=1/N$, for the PE 44 sample. This sample is chosen because its degree of polymerization is smaller than the entanglement one, $N_{e}=130$, and it is large enough to ensure Gaussian chain statistics. Once $d$ is defined, it is not changed for any other system considered, either unentangled or entangled. This is the only parameter that has to be fixed in our approach. Theoretically predicted values for unentangled systems recover the correct scaling behavior as $ND\beta\zeta_{m}\approx 1$. For entangled systems, we solve Eq.(6) by including a one-loop perturbation of the diffusion coefficient. For these entangled systems $ND\beta\zeta_{m}\propto N^{-1}$, in agreement with the known scaling behavior.

To test our approach we compare the rescaled dynamics, predicted from mesoscale simulations, with experiments and UA-MD simulations. Each sample that we investigate is in well-defined thermodynamic conditions of density and temperature, and has a specific radius of gyration. Those quantities enter as an input to our mesoscale simulation, and also in the expressions for the rescaling of the energy and friction coefficient (see Tables 1 and 2).

Comparison with simulation data are limited here to United-Atom simulations, but our theory is general and comparison could be made with atomistic simulations as well. The UA-MD simulations reported in this paper cover a regime from unentangled,MONDE to slightly entangled dynamics (two entanglements per chain).Mavran We use as an input of our approach the radius of gyration, as measured in each simulation. These values are very close to the theoretical $R_{g}$ values calculated using a FRC model with semiflexibility parameter $g=0.785$.

The experimental samples considered in this paperRichter cover a region at the crossover from unentangled to entangled dynamics comparable to the one in UA-MD simulations. However, the values of the radius-of-gyration for those samples are not known, because only the degree of polymerization is reported. For those samples we assume as input values of $R_{g}$ those calculated using a FRC approach.

The predicted cm diffusion coefficient of a polymer chain, $D_{cm}$ is compared in Figure 1 against the data from simulations,Mavran ; MONDE and from experiments.Richter We also show, as a guide to the eye, lines with the scaling behavior of unentangled and entangled systems. The agreement between calculated and measured diffusion coefficients is good over a range of the degree-of-polymerization, covering unentangled as well as entangled polymer melts. Because each simulation and experimental value is taken in slightly different thermodynamic conditions, the data points do not perfectly align along the lines of the scaling exponents in the figure. However we observe a good agreement between predicted theoretical values and measured ones in simulations or experiments.

To test further the validity of our procedure, we calculate the decay of the rotational time-correlation function for the molecular end-to-end vector with input parameters for polyethylene and the rescaled monomer friction coefficient $\zeta=k_{B}T/(ND_{cm})$, and compare it against UA-MD simulations (see Figure 2). Predicted and measured decays are in excellent agreement for the unentangled samples, suggesting that the proposed procedure holds for different normal modes of motion.

Mesocale simulations of coarse-grained systems are becoming increasingly important, as they are computationally efficient and allow for the study of systems on larger length and time scales than their atomistic counterparts. However, while structural properties on large length scales are well described by MS-MD simulations, the dynamics is unrealistically fast due to the simplified free energy landscape. In this paper we have presented an analytical, first-principles, approach to scale the dynamics measured in mesoscale simulations down to realistic atomistic values. The rescaling procedure takes into account the averaged internal degrees of freedom and enhanced dissipation due to the coarse-grainining procedure. The agreement of predicted long-time dynamics with data from simulations and experiments is quantitative. Whereas previous efforts at dynamical rescaling have used numerical calibration curves, which are specific of the system under study, our approach is analytical and thus general and transferable: it is readily applicable to systems with different thermodynamic parameters and to polymer chains of increasing degree of polymerization crossing from the unentangled to the entangled regime. The development of a general scheme to rescale the dynamics from MS-MD simulations promises to be useful in multiscale modeling techniques and fast equilibration methods employed in computer simulations of complex fluids.

In summary, the development of schemes to rescale the dynamics from MS-MD simulations will certainly be beneficial in the advancement of multiscale modeling techniques of complex fluids. Equilibrium and non equilibrium simulations of polymer melts with different architectures should be natural implementations of the approach for future work. It should also be possible to extend this model to rescale systems represented at an intermediate level of coarse-graining as collections of soft colloidal beads, so that the internal dynamics of entangled chains can be simulated.

We acknowledge support from the National Science Foundation. Simulation trajectories for the UA-MD simulations were kindly provided by Gary G. Grest and Vlasis G. Mavrantzas. We thank Glenn T. Evans for the careful reading of the manuscript and helpful suggestions.