Onset of static and dynamic universality among molecular models of polymers.

A quantitatively accurate prediction of properties for entangled polymers is a long-standing challenge that must be addressed to enable efficient development of these materials. The complex nature of polymers is the fundamental origin of this challenge. Specifically, the chemistry, structure, and dynamics at the atomistic scale affect properties at the meso and macro scales. Therefore, quantitative predictions must start from atomistic molecular dynamics (AMD) simulations. Combined use of atomistic and coarse-grained (CG) models is a promising approach to estimate long-timescale behavior of entangled polymers. However, a systematic coarse-graining is still to be done for bridging the gap of length and time scales while retaining atomistic characteristics. Here we examine the gaps among models, using a generic mapping scheme based on power laws that are closely related to universality in polymer structure and dynamics. The scheme reveals the characteristic length and time for the onset of universality between the vastly different scales of an atomistic model of polyethylene and the bead-spring Kremer-Grest (KG) model. The mapping between CG model of polystyrene and the KG model demonstrates the fast onset of universality, and polymer dynamics up to the subsecond time scale are observed. Thus, quantitatively traceable timescales of polymer MD simulations can be significantly increased.

Introduction

Entangled polymers are widely used in many industrial applications. Despite many years of basic polymer research and industrial use, quantitatively accurate predictions of polymer properties is a long-standing problem that is of great importance for efficient improvement of those properties. The complex structures, dynamics, and physical properties of polymers can change dramatically over different time scales[1-3]. However, phenomena occurring over a wide range of timescales are closely related to each other; i.e., parameters at the atomistic scale can affect properties at the meso and macro scales. Therefore, predictions of polymer properties require atomistic molecular dynamics (AMD) simulations.

While recent advances in computation have enabled a wide range of AMD simulations of chemical and biological polymers[4-12], fully atomistic simulations of entangled polymer dynamics over long timescales are beyond the capability of current platforms. Thus, coarse-grained (CG) MD simulations have been used for long timescales. They can accelerate simulations by reducing the number of degrees of freedom and increasing the fundamental time scale. Accessible polymer characteristics depend on the coarse-graining level, which determines the smallest length scale to trace the polymer dynamics. Models with various CG levels have been developed for long timescale behaviors of polymers[13-18], and some of the models have great potential for more widespread use. However, one critical issue remains unresolved: namely, systematically determining the CG level to which polymers can be coarse-grained while still appropriately tracing static and dynamic polymer properties[18-22]. One path to attack this challenge is the direct comparison of polymer properties between atomistic and CG models. All CG models can be quantitatively compared with atomistic models through a generic mapping scheme based on experimentally established universality in the structures and dynamics of polymers having different chemistries[1,2,21,23-29]. Here, we examined the gap of length and time scales among models using the mapping scheme that assumed universality for three properties closely related to polymer dynamics that have been described in the Rouse theory[30]. These are the mean square end-to-end distance, the end-to-end relaxation time, and the diffusion coefficient. The gaps were estimated by mapping between the atomistic model of polyethylene (PE) and the bead-spring Kremer–Grest (KG) model[14,17], a linkage between the vastly different scales of AMD and high-level CGMD. The gaps were also estimated by mapping between the multiscale molecular model of polystyrene (PS)[31] and the KG model, which involves linkage among the three different scales AMD, middle-level CGMD, and high-level CGMD. PE and PS are common commercial polymers that have been extensively studied experimentally and theoretically, and several reliable molecular models are currently available. The KG model is the simplest very high-level CG model; it has been widely used for decades. With the exception of plain molecular geometries and excluded volume effects, chemical details are entirely omitted in the KG model. Despite its simplicity, it is highly useful and has been used in a wide range of studies on polymer nanocomposites[32,33], polymer welding[34], polymer brushes[35-37], poly-electrolyte gels[38], thermoresponsive polymers[39], ring polymers[40], polymer collapse[41], healing of polymer interfaces[42,43], and biopolymeric motions[44]. With using the generic mapping scheme, the spatial and temporal gaps among molecular models are systematically and accurately estimated. The gaps provide useful information for determination of appropriate CG level, showing a significant advance that will enable a quantitative solution for the challenging problems described above.

Methodology

For quantitative comparison among polymer models, the length and time scales of CG models should be rescaled by a generic way. Because of the Gaussian nature of long polymers, the length unit is usually arbitrary when the linkage to atomistic models is attempted for the static properties of the global polymer structure. Therefore, we assume universality for the onset of entanglement. The onset specifies the characteristic molecular weight at which the power law exponents describing the relationships between dynamical parameters and the molecular weight change significantly. If this universality is assumed, the length units of two different polymer models can be linked to each other. Once the length unit is fixed, the unit of time can be determined. In CG models, it is determined by comparing dynamical measurements with those obtained by the corresponding atomistic model. For example, the mean square displacement (MSD) of chain centers has been used[22,31,45-51], based on reptation theory that predicts the inflection point in MSD at the characteristic time of entanglement. This strategy is useful because the characteristic time of entanglement does not depend on the molecular weight. Despite the successful use of the linkage strategy, the universality of static and dynamic polymer properties is not fully attained for molecular models used in polymer simulations. Gaussian statistics assumed in Rouse models are not observed in atomistic models unless the molecular weight is sufficiently high. The non-Gaussian nature affects the dynamics, which deviate from predictions of the Rouse model. While these deviations are often concealed in the power law expressions[52,53], this issue should be carefully considered when linking different models.

To obtain a reasonable linkage between the atomistic and CG models, we performed the following four steps. (i) The power law relations between the molecular weight M and the mean square end-to-end distance 〈R 2〉, the end-to-end relaxation time τ R, and the diffusion coefficient D, were evaluated for the atomistic and CG models (see Supplement Material for the simulation conditions). (ii) To determine the mass scaling factor, the critical molecular weight M c was estimated from the change in τ R – M and D – M power law exponents that indicate the onset of entanglement. The entanglement molecular weight, M e, was also estimated by primitive pass analysis (PPA)[54,55] (see Supplement Material for the results of PPA). (iii) The spatial and temporal scaling factors were calculated using a set of 〈R 2〉, τ R, and D for each M, under the conditions of M ≤ M c. (iv) The validity of the scaling factors was evaluated by comparison of 〈R 2〉–M, τ R–M, and D – M power law relations obtained from AMD with those obtained from CGMD, and rescaled by the mass, length, and time scaling factors. Based on these results, a simple and reasonable mapping scheme was evaluated.

Results and Discussion

Evaluation of mapping scheme

Table 1 lists the representative parameters used for the estimation of scaling factors. The variables m, σ, and τ are the units of mass, length, and time for KGMD, respectively. The mass scaling factor can be accurately determined from the M or M C results (for more details of the M e or M C values, see Supplement Material for the results of PPA). The spatial and temporal scaling factors can be determined by using any two of the following property combinations: 〈R 2〉A/〈R 2〉K, τ R,/τ R,, and DA/DK (suffixes “AMD” and “KGMD” denote the type of simulation). Therefore, three sets of two factors can be obtained from each condition of M. For the spatial scaling factor, the largest differences between the sets were 19%, 11%, and 2.2% for the smallest M, M , and M c values, respectively. For the temporal scaling factor, the largest differences between the sets were 43%, 23%, and 4.5% for the smallest M, M e, and M c values, respectively. Thus, the scaling factors converge at M = M c, irrespective of the particular combination of these three properties. The results also imply that the AMD and KGMD simulations will be in quantitative agreement for the three power laws mentioned above. Therefore, the following equations can be used for estimating the spatial and temporal scaling factors, where 〈R 2〉c and τ c are 〈R 2〉 and τ R at M = M c, respectively. Note that Eqs (1) and (2) are avaliable for any two different MD simulations. From Table 1, it is clear that D A/D K should not be used to estimate the scaling factors because the standard errors or deviations were higher for D than for 〈R 2〉 and τ R.

Table 1

Representative parameters for estimation of scaling factors.

M conditions	AMD (standard error)				KGMD (standard error)
	M [g/mol]	〈R 2〉 [nm 2]	τ R [ps]	D [10−3 nm 2/ps]	M [m]	〈R 2〉 [σ 2]	τ R [τ]	D [10−3 σ 2/τ]
Smallest	282.5	2.61962 (0.00087)	41.19 (0.22)	4.50 (0.17)	20	29.409 (0.020)	414.8 (4.5)	3.51 (0.26)
≈Me	703.4	8.5151 (0.0098)	266.2 (2.1)	1.473 (0.060)	50	80.42 (0.19)	3108 (49)	0.967 (0.078)
M c	983.9	12.515 (0.021)	552.1 (5.2)	0.914 (0.039)	70	114.72 (0.39)	6910 (130)	0.640 (0.055)

The above scaling factors should be evaluated by comparing 〈R 2〉 – M, τ R – M, and D – M power laws from AMD simulations with those obtained from the rescaled KGMD data. Figure 1(a) plots the 〈R 2〉 – M power law. The rescaled KGMD results begin to coincide with those from AMD at M = M e, and are almost equal to the AMD values at M ≥ M c. This indicates that the spatial scaling factor is reasonable. In contrast, a discrepancy between AMD and KGMD is observed at M < M e. This implies that the spatial scaling factors calculated at M < M e are inadequate for quantitative static mapping. Figure 1(b) plots the τ R – M power law. Similar to the 〈R 2〉 results, those from rescaled KGMD start to coincide with the AMD data at M = M e and are almost identical to AMD values at M ≥ M c. This indicates that the temporal scaling factor is reasonable. In contrast, a disagreement between AMD and KGMD was seen at M < M e. This implies that the temporal scaling factors calculated at M < M e are inadequate for quantitative dynamic mapping. For the D – M power law, the rescaled KGMD results are almost identical to the AMD data (see Fig. S1). Despite the discrepancies between the two types of simulations at M < M e for 〈R 2〉 and τ R (even at M < M e), the results of KGMD were close to the AMD results. This implies that the use of D may not lead to an accurate estimation of the scaling factors.

Figure 1

Comparison of power laws between AMD of PE and rescaled KGMD. (a) 〈R 2〉 – M power law. (b) τ R – M power law.

These results suggest a simple and reasonable mapping scheme consisting of three steps. (i) Estimate M c for atomistic and CG models from τ R – M (and D – M) power law. (ii) Determine the spatial and temporal scaling factors using Eqs (1) and (2). (iii) Evaluate the accuracy of the factors by comparing the 〈R 2〉 – M, τ R – M, and D – M power laws obtained from AMD with those obtained from rescaled CGMD. Note that the accuracy of the above scheme relies on the accuracy of the power laws. The uncertainty of power laws reported in previous studies[56-58] implies that a lot of data for long-time trajectories are required for the accurate estimation of the laws. Therefore we carried out long-time MD simulations of atomistic PE (up to 1 microseconds) and bead-spring KG (up to 1 billion time steps) models for a number of initial structures (see Supplement Material for the simulation conditions).

Mapping between atomistic PE and KG models

With using the present mapping scheme, the spacial and temporal gaps of representative static and dynamic properties can be accurately estimated. Figure 2(a) plots the static structure factor, (S(q), for M ≥ M e, where q is a spatial frequency that is equal to 2π/r (r is an intra- or intermolecular distance). The rescaled KGMD results begin to coincide with the AMD data at a threshold value q t = 1.78 rad/nm, indicating that the onset of static universality between AMD and KGMD occurs at a threshold length r t = 3.5 nm, which is equal to the square root of 〈R 2〉c. This implies that the onset of static universality is closely related to M c. Note that the two simulations can be reasonably linked, even though a true plateau is not reached at the highest M of the present study. Figure 2(b) plots the power law relation between M and the radius of gyration, . The rescaled KGMD results begin to coincide with the AMD values at M = M c. In contrast, the discrepancy between AMD and KGMD is observed at M < M c because the onset of static universality occurs at M = M c, in agreement with the S(q) results. The results for the radial distribution function, g(r) (see Fig. S2), are consistent with those discussed above. Thus, we have demonstrated that our mapping scheme reveals the characteristic length for the onset of universality between the vastly different scales of AMD and KGMD, for the four static properties 〈R 2〉, , S(q), and g(r). The spatial gap between AMD and KGMD is considerably large, raising the discrepancy at M < M c.

Figure 2

Results of mapping between AMD of PE and KGMD. (a) Static structure factors for M ≥ M e. (b) power law. (c) Relaxation modulus for M ≥ M e. Also plotted are [τ s,G(τ s)], where τ s, is a time at which a shoulder of G(t) occurs. (d) MSD of central monomers for M ≥ M e. Also plotted are [τ s,g 1(τ s)], that correspond to the time at which the shoulder of G(t) occurs.

Figure 2(c) plots the relaxation modulus, G(t), for M ≥ M e, where t is time. We also plotted [τ s,G(τ s)], where τ s is the time when the shoulder of G(t) occurs. For M = M e, M c, and 3509 g/mol, and τ s = 149, 316, and 12400 ps, respectively. It was confirmed that τ s is roughly proportional to τ R and the terminal relaxation time τ d (see Fig. S3). The rescaled KGMD results coincide with the AMD values at t ≥ τ s. The coincidence of the AMD and KGMD results at t < τ s is partially observed. Figure 2(d) plots the MSD of central monomers, g 1(t), for M ≥ M e. Also plotted are [τ s,g 1(τ s)], for times corresponding to the shoulder of G(t). The rescaled KGMD results coincide with AMD values at t ≥ τ s, in contrast to the large discrepancy between AMD and KGMD data at t < τ s. This indicates that the onset of dynamic universality for AMD and KGMD occurs at t = τ s, and that a delay occurs in the onset of dynamic universality depending on M, because τ s (which is roughly proportional to τ R and τ d) increases exponentially with increasing M. This delay arises from differences in chemical structure between the two simulation methods. In KGMD, the chemical details are absent, whereas in AMD, they are included explicitly. This leads to differences in the details of the modeled entanglement. Effects of these differences completely disappear at t ≥ τ s. For the large M (=3509 g/mol), the coincidence of the AMD and KGMD results at t < τ s is partially observed. This indicate that the delay in the onset of dynamic universality has a tendency to saturate with increasing M. The results of the autocorrelation function for the end-to-end vector C(t) were consistent with those for τ R (see Fig. S4). Thus, the mapping scheme reveals the characteristic time for the onset of universality between the vastly different scales of AMD and KGMD, for the five dynamic properties τ R, D, C(t), G(t), and g 1(t). The temporal gap between AMD and KGMD is considerably large, raising the discrepancy at t < τ s.

Mapping between multiscale PS and KG models

Harmandaris et al. performed multiscale MD (MSMD) simulations of PS melts[31]. The simulations were a combination of AMD and middle-level CGMD that are linked with a specialized mapping scheme. Data for 〈R 2〉, , D, and g 1(t) indicate that the rescaled CGMD provides quantitatively accurate results when compared to AMD and experimental data. Here, our mapping scheme was applied to estimate the gap between MSMD data (middle-level CGMD data rescaled by AMD data of PS) and that of KGMD. From the D – M power law of MSMD, M c, = 25000 g/mol and 〈R 2〉c, = 82.0 nm 2 (the suffix “MSMD” denotes MSMD simulation). The value τ c, = 2.59 × 108 ps is estimated from the convergence of the scaling factors discussed above: τ c,/τ c, = (〈R 2〉c,/〈R 2〉c,)(D c,/D c,)−1. The spatial and temporal scaling factors for this mapping were determined from the above values, KGMD values, and Eqs (1) and (2). Figure 3(a) and (b) plot the 〈R 2〉 – M and D – M power laws, respectively. The rescaled KGMD results are almost identical to the MSMD data, confirming the validity of the spatial and temporal scaling factors. Figure 3(c) plots the power law. The rescaled KGMD results are almost equal to those of MSMD, which indicates that the static mapping is reasonable. There was no discrepancy between two models, unlike the mapping between AMD and KGMD. Thus, the spatial gap between the middle-level CG and the KG model is smaller than that between the atomistic and KG models. Figure 3(d) plots the MSD of central monomers at M = 50000 g/mol. The interaction points per chain of atomistic, CG, and KG models are 7680 atoms (or 3840 united atoms), 960 particles, and 140 beads, respectively. This implies the efficient multiscale mapping among the three molecular models. The rescaled KGMD results are almost equal to those of MSMD, which indicates that the dynamic mapping is reasonable. Furthermore, the results for g 1(t) are a quantitative estimation of polymer dynamics up to a subsecond time scale (0.135 sec). This highlights the great potential of multiscale mapping for polymer molecular dynamics with realistic time scales. A direct comparison between MD results and experimental data becomes much easier for polymer dynamics, and enables high-throughput screening in polymer development. There was no discrepancy between two models, unlike the mapping between AMD and KGMD. The onset of universality in g 1(t) is much faster (t ≪ τ c) than that for mapping between AMD and KGMD. Thus, the temporal gap between the middle-level CG and the KG model is smaller than that between the atomistic and KG models.

Figure 3

Results of mapping between MSMD of PS and KGMD. (a) 〈R 2〉 – M power law. (b) D – M power law. (c) power law. (d) MSD of central monomers at M = 50000 g/mol.

Again we focus on the characteristic of the spatial and temporal gaps. The large gaps are observed between AMD and KGMD despite the number of united-atoms per PE chain is almost equal to that of beads per KG chain. In contrast, The gaps between middle-level CGMD and KGMD is small despite the number of CG particles per PS chain is about 7 times larger than that of beads per KG chain. This means that the difference of chemical details (e.g., intra- and intermolecular interactions) is more effective for the gaps than the difference of the plain geometries of polymer models. Significantly, the gaps become small when linking molecular models having not less than a certain CG level. Estimating this level is critical for the development of CG models, because it provides the threshold for the effective coarse-graining achieving a good balance between model accuracy and computational cost. Here the present mapping scheme can reasonably determine the spatial and temporal gaps for representative polymer properties. The quantified gaps are useful as the index for searching the appropriate CG level. Therefore, optimization techniques of coarse-graining for polymer molecular models will be systematically extended by using the present scheme.

Conclusions

The application of universal polymer behavior to the linkage among molecular models at various length and time scales was examined in terms of a mapping scheme focused on power laws that also show universality. The spatial and temporal gaps were evaluated by mapping between the atomistic model of PE and the coarse-grained KG model, a linkage between the vastly different scales of AMD and high-level CGMD. The scheme reveals the characteristic length and time for the onset of universality between AMD and KGMD for the representative static and dynamic properties. The large gaps raise the discrepancies of the properties at short length and time scales, showing the clear limit for bridging the two simulations while retaining atomistic characteristics. To overcome this limit, low-level (i.e., elaborate) CG models that introduce the effects of short length and time scales should be developed[18,20,22]. The scheme was then use to evaluate the gaps between the multiscale molecular models of PS and the KG model, which involves a linkage among the three different scales of AMD, middle-level CGMD, and high-level CGMD. The small gaps between middle- and high-level CGMD attain the fast onset of universality without the discrepancies of the properties at any length and time scales. Furthermore, it allows a linkage between the two models that quantitatively estimates polymer dynamics up to the subsecond timescale. This demonstrates the great potential to that the combination of various scales of molecular models permits multiscale mapping for polymer MD simulations. Here, quantitatively traceable MD time ranges can be significantly scaled up, which introduces the possibility of direct linkage between MD simulations and macroscale approaches[59,60]. The present mapping scheme is simple, and the scaling factors can be estimated with reasonable accuracy by using existing MD simulations and computational resources. Therefore, it is anticipated that it will be applicable to the sequential multiscale mapping among AMD, low-, middle-, and high-level CGMD for quantitative estimation of polymer properties at realistic length and time scales while maintaining a reasonable computational cost.