Performance of Coarse Graining in Estimating Polymer Properties: Comparison with the Atomistic Model.

Combining atomistic and coarse-grained (CG) models is a promising approach for quantitative prediction of polymer properties. However, the gaps between the length and time scales of atomistic and CG models still need to be bridged. Here, the scale gaps of the atomistic model of polyethylene melts, the bead-spring Kremer-Grest model, and dissipative particle dynamics with the slip-spring model were investigated. A single set of spatial and temporal scaling factors was determined between the atomistic model and each CG model. The results of the CG models were rescaled using the set of scaling factors and compared with those of the atomistic model. For each polymer property, a threshold value indicating the onset of static or dynamic universality of polymers was obtained. The scaling factors also revealed the computational efficiency of each CG model with respect to the atomistic model. The performance of the CG models of polymers was systematically evaluated in terms of both the accuracy and computational efficiency.

1. Introduction

Because of scientific and industrial interest, entangled polymer materials have been studied for many years. However, accurate quantitative prediction of polymer properties is still a long-standing issue because of the slow relaxation nature of the polymer dynamics far from the nanoscale of the molecular architecture [1,2,3]. Atomistic molecular dynamics (MD) simulations trace the dynamics of polymer molecules using realistic length and time units in a straight forward manner, but the simulation time (up to a few microseconds) is not sufficient to cover polymer relaxation (up to seconds, minutes, hours, months, or years), even using current high-performance computers. Thus, coarse graining of molecular models and simulations was developed to cover larger timescales [4,5,6,7,8,9]. To accelerate simulations, in most cases, the number of degrees of freedom is reduced and a Langevin equation of motion is often used. Importantly, the realistic units of length, time, temperature, and so forth are lost during such coarse graining, and thus coarse-grained (CG) MD (CGMD) simulations only return qualitative values. These “scale gaps” make it difficult to evaluate the accuracy of the estimated polymer properties and determine how much MD simulations are accelerated using CG models. To solve this problem and quantitatively evaluate CG models, the results of CGMD simulations should be rescaled by using realistic units. A general way to rescale the results is application of the scaling law of polymer dynamics experimentally determined for polymer materials [1,2,10,11,12,13,14,15,16,17]. For simple straight-chain polymers, the scaling law can be determined regardless of the difference of the chemistry of the molecular architecture. Importantly, the above character of the law can be applied to different molecular models of polymers even though each model has a different concept, such as the presence/absence of chemistry, employment of the Newton/Langevin equation of motion, and so forth. This was investigated in several previous reports. Hamandaris et al. [18] compared the polymer dynamics of the atomistic polystyrene (PS) model with that of the CG PS model. The CG PS model was rescaled by a temporal scaling factor determined by comparison of the unentangled short chain dynamics. Masubuchi et al. [19] investigated the compatibility of the dynamical and static polymer properties among different models of the multi-chain slip-link model, multi-chain slip-spring model, and bead–spring model. Polymer diffusion of the three models coincided by using spatial and temporal scaling factors, showing the compatibility of the models. Takahashi et al. [20] compared the polymer dynamics of the atomistic PE model with that of the bead-spring model. They determined a set of spatial and temporal scaling factors to successfully rescale the results of the bead–spring model, but the scaling behavior around the onset of entanglement was different for each polymer property. This quantifies the applicable range of bridging between the models and suggests the possibility of further improvements for CG models. Nevertheless, few studies evaluated both the accuracy and computational efficiency of CG models rescaled by the above scheme.

In this work, we investigated the scale gaps of the atomistic model of polyethylene (PE) melts, bead–spring Kremer–Grest (KG) model [5,8], and dissipative particle dynamics (DPD) with slip-spring model [21]. A single set of spatial and temporal scaling factors was determined between the atomistic model and each CG model. The results of the CGMD simulations were rescaled using the set of scaling factors and compared with those of the atomistic model. Threshold values that indicate the onset of static and dynamic universality among the polymer models were obtained for almost all of the results of comparison of the polymer properties, showing the accuracy of the CG models. The computational efficiency of each CG model was quantitatively estimated with respect to the atomistic model. The performance of the CG models of polymers was systematically and quantitatively evaluated in terms of both the accuracy and computational efficiency.

2. Methodology

2.1. Bead–Spring KG Model

In this work, the bead–spring model for polymer chains developed by Kremer and Grest [8] was used as the coarse-grained model without chemical parameters. The model consists of simple beads whose non-bonded interactions are described by the following repulsive LJ potential: where is the unit of energy, is the unit of length, and r is the distance between two beads. The bonds between beads are described by the following finite extensible nonlinear elastic (FENE) potential: where is the spring constant and is the maximum extension of the spring. For polymer motion in the KG model, the following Langevin-type equation is used: where is the mass of the KG bead, are the coordinates of the bead, is an arbitrary friction constant, is a random force, is Boltzmann’s constant, and is the temperature. , , and T were set to 1, and was set to 0.5. All of the simulations using the KG model were performed using the ensemble. The monomer density was . The velocity Verlet integrator [22] was used with three-dimensional periodic boundary conditions and a time step of , where is the time unit of the KG model. The number of KG beads () varied from 20 to 250. For –80, equilibrium simulations with a total of time steps were performed for three independent initial structures. For –250, equilibrium simulations with a total of time steps were performed for six independent initial structures. The total number of beads in the system () was fixed at about 25,000. Each KG MD (KGMD) simulation system was equilibrated prior to data acquisition.

2.2. DPD with the Slip-Spring Model

We used the DPD method [6,23,24] to reproduce the behavior of the polymer melt. The DPD method was proven to be an effective mesoscopic simulation tool to study fluid events occurring at millisecond timescales and micrometer length scales via by tracking the motion of coarse-grained particles (composed of a group of atoms or molecules). Many researchers studied the thermodynamic and morphological behavior of polymers using the DPD method, for example, self-assembly of dendritic copolymers [25,26,27], micelle fusion or vesicle formation [28,29], and nanocomposite systems [30,31,32,33].

The fundamental equation in the DPD method is Newton’s equation of motion. Newton’s equation of motion for particle i is given by where is the velocity vector, is the conservative force, is the pairwise random force, and is the dissipative force. The conservative force is softly repulsive: where and . Here, is a parameter to determine the magnitude of the repulsive force between particles i and j, and is the cutoff distance. The random force and dissipative force are given by and respectively, where , is the noise parameter, is the friction parameter, and is a random number based on the Gaussian distribution. Here, and are r-dependent weight functions:

The temperature is controlled by a combination of dissipative and random forces. The noise parameter and friction parameter are connected by the fluctuation–dissipation theorem:

Usually, reduced units are used for reporting DPD results. The DPD unit of length is the cutoff radius (), the unit of mass is the DPD bead mass (), and the unit of energy is . Hence, = = = 1.

In this work, the number of DPD beads was varied from 5 to 100. Each particle and its nearest neighbor in the polymer are connected by a harmonic spring described by with a spring constant of k = 2 [21].

All of the simulations were performed with the ensemble. The particle density and temperature were set to 3.0 and 1.0 , respectively. Periodic boundary conditions were applied in all three dimensions. In the DPD simulation, the thermostat is achieved by pairwise random and dissipative forces. These forces are coupled by the fluctuation-dissipation theorem [24]. The noise amplitude and friction coefficient are included in the dissipative and random forces, respectively. They satisfy = , which reproduces the canonical ensemble. and were set to 3.0 and 4.5, respectively.

The total number of beads in the system () is 3000 to 24,000 depending on the chain length. In our simulations, the interaction parameters in Equation (1) are described by = 25 .

In the DPD method, a soft potential is also used as the interaction between particles. However, it is known to be unable to reproduce entanglement effects [34]. To capture the entanglement behavior of the polymer, slip-springs [21,35,36] were introduced with the DPD polymer chains. The details of the procedure are given elsewhere [21,36]. The number of slip-springs is 300 to 2400, namely a slip-spring is placed for every ten beads.

Each run was performed for time steps with a time step , where is the time unit for DPD with the slip-spring model, to reproduce the equilibrium polymer melt. The reported data are the average over the last time steps.

2.3. Atomistic Model of PE

Most accurate classical MD simulations are performed using all-atom (AA) molecular models. However, it is often beyond the capability of current computers to trace the slow polymer motion of the AA model up to a meaningful timescale [18]. Thereforem united-atom (UA) molecular models were developed and used as quantitative atomistic models. In UA models, a group of atoms, such as , is merged into one UA particle whose non-bonded interactions are described by Lennard-Jones (LJ) potentials. Such UA models have three main benefits in terms of computational efficiency: (i) the degrees of freedom can be reduced, (ii) the time step for molecular motion can be extended, and (iii) computationally expensive Coulomb interactions can be omitted for UA particles that do not have an overall charge. Because of these benefits, UA models “accelerate” the timescale of AA molecular motion. A previous study showed that the dynamics of AA polystyrene (PS) systems is accelerated by about 120 times by using the UA PS model [18]. Importantly, this acceleration ratio between AA and UA models cannot be evaluated before direct comparison of the dynamics of the two simulations. This indicates that the UA model has a scale gap problem, like CG models. However, in terms of only computational efficiency, UA models are used for quantitative molecular simulations. In fact, UA models are widely used for MD simulations of polymers. For UA models of the common polymer PE, the anisotropic UA (AUA) [37], optimized potentials for liquid simulations (OPLS)-UA [38], and transferable potentials for phase equilibria force field (TraPPE)-UA [39] models are widely accepted [40,41,42,43,44,45,46]. MD simulations using UA PE models were recently performed for a wide range of polymer nanocomposites [47,48], polymer interfaces [49,50,51], ring polymers [52], nucleation of polymer droplets [53], the Fermi–Pasta–Ulam problem in realistic systems [54,55], and better understanding of the macroscopic mechanical properties [56,57,58,59,60].

In this work, the TraPPE-UA [39] model was used for the atomistic model of PE melts. The TraPPE-UA PE model consists of two different types of UAs ( and ), whose non-bonded interactions are described by LJ 12–6 potentials. The bond lengths between the UA particles are kept rigid using the LINCS algorithm [61]. The bond angle bending interactions are described using harmonic potentials. Rotations along the bonds in the aliphatic backbone are described using standard torsional potentials. These dihedral potentials are also used for the 1–4 non-bonded interactions. Atomistic MD (AMD) simulations using the above UA model were performed for PE melts with molecular weight (M) ranging from to . Please note that the critical molecular weight of the model was estimated to from the scaling law and primitive path analysis [20]. Therefore both unentangled and entangled motions of PE chains can be observed in this work. The MD package GROMACS [62] was used for the simulations. To obtain well-equilibrated initial atomistic structures, over long-time MD simulations with intermittent pressure increasing and temperature decreasing processes were performed with the constant particle number, pressure, and temperature (NPT) ensemble. The equilibrated systems were confirmed by comparison with the results of a previous study [46]. The product runs using the equilibrated initial structures were performed with the constant particle number, volume, and temperature (NVT) ensemble using the Nosé–Hoover thermostat [63,64]. The density and temperature were set to and , respectively. Non-bonded interactions were cut off beyond . The Verlet leapfrog integrator [65] was used with three-dimensional periodic boundary conditions and a time step of . For –, equilibrium simulations with a total of time steps (=) were performed for three independent initial structures. For –, equilibrium simulations with a total of time steps (=) were performed for six independent initial structures. For –, equilibrium simulations with a total of time steps (=) were performed for six independent initial structures.

For easy understanding, here we state the two clear differences between AMD using TraPPE-UA and KGMD. The one is a force field. In TraPPE-UA, a bond, angle, torsion, and non-bonded repulsive and diffusive forces are quantitatively considered based on quantum and molecular mechanics. In contrast, KGMD considers only qualitative FENE-bond and non-bonded repulsive forces. The other is an equation of motion. In AMD, Newton’s equation of motion is simply but strictly used. In KGMD, Langevin equation is modified for polymer dynamics, and the arbitrary parameter is employed for easy expression of friction between polymer chains.

2.4. Rescaling of the CG Models

Regarding the unentangled/entangled nature of polymer molecules, some scaling laws of the polymer properties have a critical molecular weight that indicates the change of the polymer dynamics. For simple explanation of the scaling laws, the Rouse model [66] and reptation model [1] were suggested for unentangled motion () and entangled motion (), respectively. Thus, the scaling factor to rescale the CGMD simulation results can be determined by comparison of the unentangled or entangled polymer dynamics based on the Rouse or reptation model. Determination from unentangled motion is beneficial in terms of the computational efficiency. Generally, relaxation of unentangled motion is much faster than that of entangled motion. Thus, the scaling factor can be determined with less effort. Hermandaris et al. [18] compared the unentangled polymer dynamics of the atomistic model with that of the CG model by using the mean-squared displacements (MSDs) of the short chain motions, and they determined the temporal scaling factor based on the Rouse model. They observed hundreds of nanoseconds of polymer dynamics by quantified CGMD simulations of PS melts. Nevertheless, determination of the scaling factor based on the Rouse model should be carefully performed. The scaling law of the Rouse model is supported by the Gaussian nature of the motion of the segment with sufficiently high molecular weight (i.e., the long chain). However, several atomistic and CG models do not show Gaussian nature at [67,68]. The non-Gaussian nature becomes the cause of deviation of the scaling law from that of the Rouse model. Moreover, there is a concern that a lot of accurate CG models have non-Gaussian nature in the unentangled state. Saleno et al. [69] compared the polymer dynamics of the atomistic PE model with that of the coarse-grained PE model for different degrees of coarse graining. They found that a small degree of coarse graining (4–6 monomers per coarse-grained particle) is required to correctly describe the macroscopic behavior of the polymer properties. The small degree of coarse graining indicates non-Gaussian nature of the coarse-grained particle. Therefore, the deviation of the scaling law for unentangled motion should be carefully examined for each molecular model.

Another strategy is to use the scaling law of entangled motion. This is computationally expensive, but it is a promising method because some scaling laws are robust with respect to the degree of coarse graining for . In other words, the entangled motion based on the reptation model surpasses the non-Gaussian nature of the segment that consists of the short chain. In this work, we determined the spatial and temporal scaling factors focusing on the entangled motion. Takahashi et al. [20] focused on the critical point of the scaling laws of the self-diffusion coefficient and end-to-end relaxation time with respect to the molecular weight. They determined a single set of scaling factors that attain the best accuracy for rescaling the KG model to the units and dimensions of the UA model of PE. However, their scheme was optimized for comparison of the above two models. To determine the scaling factors in a more general and robust way, we performed the following six steps: (i) for each molecular model, compute the scaling laws with respect to the molecular weight for the self-diffusion coefficient D–M, end-to-end relaxation time –M, mean-squared end-to-end vector –M, and mean-squared radius of gyration –M (some other kinds of scaling law are basically addable), (ii) determine the critical molecular weight of each molecular model from the critical point of D–M or –M scaling laws, and determine a scaling factor of the mass, (see Table 1 for actual values), (iii) define the undetermined scaling factors for length and time, (iv) rescale the scaling laws of the CG models using the undetermined scaling factors, (v) calculate the relative root-mean-squared deviations (rRMSDs) between the scaling laws of the atomistic model and those of the rescaled CG model, and (vi) determine the set of spatial and temporal scaling factors that give the smallest sum of the rRMSDs. The definition of rRMSD is as follows: where is the rRMSD of the scaling law for property X, p and q are undetermined scaling factors, n is the number of data points at , and is the deviation of property X from the results of AMD.

Table 1

Scaling factors of the UA model of PE melts, KG model, and DPD with slip-spring model. “Carbon num.” becomes the “dimension” of segment number.  1 e.u. = 6.9 × 10 J. 1 e.u./nm = 6.9 × 10 J/m = 6.9 MPa.

Property	Dimension	Symbol for	UA	(Unit)	KG	(Unit)	DPD	(Unit)
		Scaling Factor
—	Length	l^*	1	nm/nm	0.33	nm/σ	0.57	nm/r_c
—	Time	t^*	1	ps/ps	0.080	ps/τ_KG	0.028	ps/τ_DPD
—	Mass	m^*	1	g/mol/(g/mol)	14	g/mol/m_KG	32	g/mol/m_DPD
—	Energy	—	1	J/J	6.9×10^−21	J/(k_BT)_KG	6.9×10^−21	J/(k_BT)_DPD
—	(Carbon num.)  ^a	N^*	1	(carbons/carbons)	1.0	(carbons/KG segs.)	2.3	(carbons/DPD segs.)
k_BT	Energy	(k_BT)^*	1	e.u./e.u.  ^b	1.0	e.u./e.u.	1.0	e.u./e.u.
Segment density	(Carbon num.)/Length^3	ρ^*	1	(carbons)nm^−3(carbons)nm^−3	0.85	(carbons)nm^−3(carbons)nm^−3	1.3	(carbons)nm^−3(carbons)nm^−3
Critical segment number	(Carbon num.)	Nc*	1	(carbons/carbons)	1.0	(carbons/carbons)	1.0	(carbons/carbons)
Strand density	Length^−3	ν^*=ρ^*/Nc*	1	nm^−3/nm^−3	0.85	nm^−3/nm^−3	1.3	nm^−3/nm^−3
G^d	Energy/Length^3	(See Equation (22))	1	e.u.·nm^−3e.u.·nm^−3  ^c	0.036	e.u.·nm^−3e.u.·nm^−3	0.036	e.u.·nm^−3e.u.·nm^−3
G^u	Energy/Length^3	G^u	1	e.u.·nm^−3e.u.·nm^−3	1.0	e.u.·nm^−3e.u.·nm^−3	0.19	e.u.·nm^−3e.u.·nm^−3
A	—	A^*	1	—/—	1.2	—/—	0.14	—/—

3. Results and Discussion

3.1. Scaling Laws

First of all, the D–M, –M, –M, and –M scaling laws are computed. Takahashi et al. [20] shows that the critical molecular weight can be determined from the change in D–M and –M power law exponents that indicate the onset of entanglement, and the determined value corresponds to the critical molecular weight estimated from primitive path analysis [70,71,72]. Therefore, we employ their strategy as a simple and robust way to determine the critical molecular weight. The –M and –M scaling laws are shown in Figure 1a,b, respectively. For –M, the rescaled KGMD results are almost equal to the AMD results for , whereas the rescaled DPD results slightly deviate from the AMD results for . However, for –M, the rescaled KGMD and DPD results are almost equal to the AMD results for . The above results indicate that the determined spatial scaling factor has adequate accuracy. The –M scaling law is shown in Figure 1c. The rescaled KGMD and DPD results are almost equal to the AMD results for . This indicates that the temporal scaling factor is reasonable. The D–M scaling law is shown in Figure 1d. The rescaled KGMD and DPD results are almost equal to the AMD results for . Please note that both the spatial and temporal scaling factors were required for rescaling D. The results indicate that the set of scaling factors is reasonable. The scaling factors determined from our scheme are given in Table 1.

Figure 1

Comparison of the scaling laws for the UA PE model and rescaled CG models. (a) –M scaling law. (b) –M scaling law. (c) –M scaling law. (d) D–M scaling law.

3.2. Static Properties

In the following sections, scaling factors determined in the previous section (see Table 1) were used for rescaling the KGMD and DPD results to real units.

The scale gap can be approximately evaluated from the scaling factors, but its details should be carefully determined by comparison of the physical properties. To reveal the details of the spatial scale gap, some static properties of the MD simulations were compared. The radial distribution function is useful to reveal the local structure of polymer melts: where is the number density of particles and is the number of particles (i.e., carbon atoms/beads) in the region between r and when particle j is on . The term can be defined for the total, intermolecular, and intramolecular contributions. describes the distance dependence of the probability for beads/carbon atoms to meet. for all particles, intermolecular particles, and intramolecular particles is shown in Figure 2a–c, respectively. For all of the results, the peaks of the three models do not coincide at nm, but the plateaus/curves at nm do coincide. For DPD, violation of the excluded volume effect is observed, showing that the DPD particle can pass other DPD particles even after implementing slip-springs.

Figure 2

Comparison of the radial distribution function (RDFs) for the UA PE model and rescaled CG models. (a) RDF for all particles. (b) RDF for intermolecular particles. (c) RDF for intramolecular particles.

The static structure factor of an individual chain is useful to determine the small and large internal structures of polymer melts: where q is the spatial frequency equal to . The fractal scattering of is expected to be equal to () and independent of the chain length. The static structure factors for the three models are shown in Figure 3. The curves coincide for rad/nm, but the fractal scattering of the CG models is clearly different from that of the atomistic model for rad/nm. For the CG models, the fractal scattering for 1.0 rad/nm 5.0 rad/nm is almost the same as that for an ideal chain. In contrast, the fractal scattering of the atomistic model for 2.0 rad/nm 10 rad/nm is clearly different from the ideal value. This indicates that the expected cancellation of the dispersion forces for polymer melts [73] is not entirely satisfied for the atomistic model. The fractal scattering of the atomistic model for 2.0 rad/nm 10 rad/nm is estimated to be . Please note that the fractal scattering of the atomistic model becomes almost ideal for 1.0 rad/nm 2.0 rad/nm, indicating the onset of cancellation of the dispersion forces at q=2.0 rad/nm (i.e., nm). In other words, of the atomistic model and CG models do not coincide at length scales smaller than 3.1 nm owing to the effect of dispersion forces. Even though the LJ interactions were cut off beyond 1.2 nm, internal structures up to 3.1 nm were subject to the influence of dispersion forces.

Figure 3

Comparison of the static structure factors for the UA PE model and rescaled CG models.

For the static properties, the largest threshold value of length is 3.1 nm. This is approximately equal to at , indicating mismatch of the static structures between the CG models and the atomistic model within the distance between the strands of the polymer network.

3.3. Dynamic Properties

To clarify the details of the temporal scale gap, some dynamic properties of the MD simulations were compared. The MSD of the chain center shows the relaxation motion of mass transfer: where is the coordinates of the chain center. From the Rouse and reptation models, the scaling-law sequence for the MSD is expected to be approximately where is the segment number corresponding to , N is the segment number per chain, is a specific short time, and and are the Rouse relaxation times that correspond to and N, respectively. The curves for the three models are shown in Figure 4a. For (=980 g/mol), the curves coincide for ps. For g/mol, the curves coincide for 10,000 ps. These threshold values of 550 ps and 10,000 ps are proportional to , showing the delay in the onset of dynamic universality for the polymer melts. The matching between the KG and atomistic models is slightly better than that between the DPD and atomistic models. However, the shape of the curve for the DPD model is close to that of the KG model. The threshold values of Equation (17) are unclear from the results themselves. Thus, we plotted against time (Figure 4b). The scaling law for is approximately satisfied at about 10,000 ps, irrespective of the differences of the models. In contrast, the scaling laws of the CG models for deviate from those of the atomistic model. This indicates the difficulty in determining accurate scaling factors from the MSDs of short chain motions. The shape of the curve for the DPD model is close to that of the KG model, indicating the similarity of the short chain motions of the two models. Please note that D can be determined from because is expected for .

Figure 4

Comparison of the MSD of the chain center for the UA PE model and rescaled CG models. (a) . (b) .

The time-correlation function of the end-to-end vector shows the relaxation motion of the longest chains:

Logarithmic plots of for the three models are shown in Figure 5a. For , the curves coincide for ps. For g/mol, the curves coincide for ps. However, the DPD results at g/mol are slightly smaller than the others for ps. This behavior corresponds to the results of . For DPD, the value is slightly larger than the others, i.e., mass transfer is slightly faster.

Figure 5

Comparison of the time-correlation functions of the end-to-end vector for the UA PE model and rescaled CG models.

Please note that can be determined from because it is expected that .

The relaxation modulus shows the viscoelastic behavior of polymer melts: where are the off-diagonal stress components , , and . For the rheological property, it is generally believed that the scaling factor should be determined from the rheological response. Indeed, mismatch of the scaling factor was demonstrated between direct determination from comparison of the rheological response and the use of the unit of (energy per volume). For convenience, we defined two different types of G scaling factors: the G scaling factor from direct comparison and the G scaling factor after using the unit . Thus, and satisfy the following relation: where is the scaling factor for . The and values are given in Table 1. rescaled by the scaling factor of energy per volume () based on the unit of the relaxation modulus is shown in Figure 6a. The KGMD results are almost the same as the AMD results, despite not using . This is because of the similarity of the polymer networks. For rubber elasticity theory, the shear modulus of polymer networks can be determined by where and are the number density of the network strand and segment, respectively, A is a constant dependent on the network functionality and fluctuations imposed to the network nodes, and is the segment number between entanglements. Please note that is equal to , which is referred to as the “Kuhn length of the primitive path” and can be determined not only by the scaling law but also by primitive path analysis [70,71,72]. From Equations (20) and (21) and Table 1, can be rewritten as

Figure 6

Comparison of for the UA PE model and rescaled CG models. (a) rescaled by the scaling factor of energy per volume () based on the unit of the relaxation modulus. (b) using the scaling factor determined by Equation (22).

Thus, the similarity of the polymer networks can be examined using the A and values. For KGMD, the value is 0.85 times that of AMD, indicating that the strand densities of the two model are approximately the same (see Table 1). From Equation (22) and , the value was determined to be 1.2, showing that the polymer network of the KG model is very similar to that of the atomistic model. In contrast, the results of DPD greatly deviate from those of AMD. For DPD, the value is 1.3 times that of AMD, indicating that the strand densities of the two models are relatively close. However, between AMD and DPD is 0.19. Thus, the value is 0.14, showing that the polymer network of the DPD model is considerably different from that of the atomistic model. Nevertheless, using Equation (22), the results successfully rescaled for DPD. using the scaling factor determined by Equation (22) is shown in Figure 6b. All of the results coincide for ps.

For the dynamic properties, the largest threshold value of time is proportional to , showing the delay in the onset of dynamic universality for the polymer melts. This differs from the onset of static universality in terms of whether the observed threshold value continues to depend on the chain length. If the threshold value of time continues to depend on the chain length, it becomes difficult to effectively switch the use of AMD and CGMD with respect to the timescale. The delay in the onset of dynamic universality may disappear for relatively large molecular weight, but this condition is beyond the capability of current computers.

3.4. Computational Efficiency

The computational cost of MD simulations can generally be estimated by the number of particles in the simulation system and the number of times of numerical integration . Recent advances in calculation methods have made it possible to reduce the computational cost of interaction calculations by up to while maintaining the accuracy, where O is a function indicating the order. Thus, the total computational cost of a MD simulation can be determined by . When considering numerical integration up to the same timescale for systems with the same particle scale, the relations of the total time step and number of particles between CGMD and AMD can be expressed as

Thus, the ratio of the total computational cost of CGMD to that of AMD can be estimated by where is the time step for numerical integration of the equation of motion. For KGMD, the total computational cost is 4.2 times that of AMD, indicating that the computational efficiency of the KG model is lower than that of the UA model. This is partly because of that the number of beads (or site) in single chain of KG model is same as that of UA model. Ordinary, CG models for polymer chains supposes that one segment of CG models contains several repeat units. In this work, however, “repeat units per segment” is simply determined from mass (or carbon number) scaling factor shown in Table 1. Therefore “one unit per segment” is basically possible, and an inversion of computational efficiency is also possible eventuality. The inversion of computational efficiency between UA and KG model in this work can be thought of as an overestimating of arbitrary parameter in KGMD. The large makes dynamics slow because the friction between chains becomes large [74]. It can be assumed that the of AMD using TraPPE-UA is considerably smaller than that of KGMD. This fact indicates an existence of real case that the optimization of basic parameter of KG model becomes necessary for efficient coarse graining. For DPD, the total computational cost is 0.52 times that of AMD. Thus, coarse graining of the UA model using the DPD model has the benefit of speeding up the MD simulations by about 2 times. The DPD with slip-spring model can be considered to be an option for coarse graining of the UA PE model.

In Table 2 we summarize the overall performance of UA, KG, and DPD models to simulate PE melts.

Table 2

Comparison of overall performance of molecular models for PE melts simulations. “A”, “B”, “C”, and “D” mean “very good”, “good”, “adequate”, and “no good” performance, respectively.  “very good” after using Equation (22).

		UA	KG	DPD
Static properties	〈R^2〉	A	B	C
	〈RG2〉	A	C	C
	g(r)	A	B	B
	S(q)	A	B	B
Dynamic properties	D	A	A	A
	τ_R	A	A	B
	g_1(t)	A	B	B
	C(t)	A	B	C
	G(t)	A	A	A ^a
Computational cost		C	D	B

4. Conclusions

We investigated the scale gaps among the atomistic model of PE melts, bead–spring KG model, and DPD with slip-spring model. First, the single set of spatial and temporal scaling factors between the atomistic model and each CG model was determined. The method to determine the scaling factors was systematic and robust because some representative scaling laws of straight-chain polymer melts that can be determined from experimental/computational results were used. The rescaled scaling laws of the CG models coincide with that of the atomistic model for , showing that the entangled motion is more robust than the non-Gaussian nature of the segment which consists of the short chain. The results of the CGMD simulations were rescaled using the well-estimated scaling factors and compared with those of the atomistic model. Threshold values that indicate the onset of the static and dynamic universality among the polymer models were obtained for almost all of the results of comparison of the polymer properties. For the static properties, the largest threshold value of the length is 3.1 nm. This value is approximately equal to at , indicating the mismatch of the static structures between the CG models and the atomistic model within the distance between the strands of the polymer network. For the dynamic properties, the largest threshold value of time is proportional to , showing the delay in the onset of the dynamic universality for the polymer melts. This is in contrast to the onset of static universality in terms of whether the observed threshold value continues to depend on the chain length. If the threshold value of time continues to depend on the chain length, it becomes difficult to effectively switch the use of AMD and CG MD with respect to the timescale. The delay in the onset of dynamic universality may disappear for relatively large molecular weight, but this condition is beyond the capability of current computers. From the results, and thus the A value, the similarity of the polymer networks among the molecular models was quantitatively evaluated. The KG model has a similar polymer network to the atomistic model, whereas the polymer network of the DPD model considerably different from that of the atomistic model. The computational efficiency of each CG model was quantitatively estimated with respect to the atomistic model. The total computational cost of KGMD is 4.2 times that of AMD, whereas that of DPD is 0.52 times that of AMD. The KG model is clearly a poor choice for coarse graining of the UA PE model, whereas the DPD with slip-spring model is a possible choice for coarse graining, and it has the benefit of speeding up the MD simulations by about two times.

The performance of the CG models of polymers was systematically and quantitatively evaluated in terms of both the accuracy and computational efficiency. Using our scheme, the following three points can be quantitatively evaluated for CG models: (i) the “scale gap” of the static and dynamic properties from that of the atomistic model, (ii) the similarity of the polymer network to that of the atomistic model, and (iii) the computational efficiency relative to the atomistic model. The above three points are crucial elements for further development of better CG models and/or coarse-graining schemes. In future studies, the performance of other CG models will be extensively evaluated using our scheme.