The Shinoda-DeVane-Klein (SDK) model is herein demonstrated to be a viable coarse-grain model for performing molecular simulations of polyethylene (PE), affording new opportunities to advance molecular-level, scientific understanding of PE materials and processes. Both structural and dynamical properties of entangled PE melts are captured by the SDK model, which also recovers important aspects of PE crystallization phenomenology. Importantly, the SDK model can be used to represent a variety of materials beyond PE and has a simple functional form, making it unique among coarse-grain PE models. This study expands the suite of tools for studying PE *in silico* and paves the way for future work probing PE and PE-based composites at the molecular level.

## I. INTRODUCTION

Polyethylene (PE) is a widely used material with an annual production in excess of 110 × 10^{6} tonnes, and with applications ranging from packaging to construction materials.^{1} Therefore, there is significant interest in understanding both PE processing^{2,3} and the structures that PE adopts upon cooling.^{4–11} As PE-based materials and processing become increasingly sophisticated (e.g., PE nanofibers^{12}), computer simulations have the potential to continue providing molecular-level insights into PE, shaping scientific understanding of complex materials, and informing technological advancements. Given that industrially utilized PE comprises long, entangled polymer chains, future work will likely require large-scale simulations of polymer mixtures at long time scales (e.g., milliseconds and beyond). In contrast, atomistic polymer simulations are generally limited to time and length scales on the order of 100 ns and 10 nm.^{13} United-atom (UA) models^{14–18} address such limitations by representing each heavy atom and its hydrogen atoms (i.e., each methylene unit in the case of PE) as a single interaction site; however, the time and length scales accessible to UA models are still fairly limited. Coarse-grain models that group multiple heavy atoms and corresponding hydrogen atoms into effective interaction sites (i.e., so-called beads) are needed to realize future transformative simulation-based studies.

To this end, a variety of coarse-grain molecular models^{19–25} have been developed for PE where groups of adjacent methylene units along the polymer chain backbone are represented by coarse-grain beads. Usage of a coarse-grain model reduces the number of particles required to represent a PE chain, enabling faster computations. Coarse-grain models also generally exploit simpler interactions compared to their atomistic counterparts, further enhancing the time and length scales that are amenable to molecular simulations. However, coarse-grain PE models have generally been developed for the sole purpose of simulating pure PE melts (e.g., see Refs. 20 and 22–26). It remains unclear which of the previously developed coarse-grain models are capable of properly capturing phenomenology related to PE crystallization and the semicrystalline solid state of PE. In an attempt to provide a route to studying PE crystallization with coarse-grain simulations, Moyassari *et al.*^{19} recently developed several new coarse-grain models for simulating neat PE at different temperatures; however, leveraging their models to study PE phase behavior requires switching models as temperature changes, potentially leading to artifacts. Moreover, their models, like other coarse-grain PE models, are only appropriate for pure PE systems and are not part of larger molecular modeling suites for mixed material systems. Leveraging existing coarse-grain PE models to conduct simulations of PE mixtures (e.g., PE composites) would require substantial model development and optimization. Consequently, it remains an open challenge to realize a coarse-grain model that is appropriate for representing PE under disparate conditions, and within the context of mixed systems, while also rendering long time and length scales accessible.

Instead of developing a new coarse-grain PE model, this study adopts an alternative strategy to address the aforementioned challenge. More specifically, this study confirms that the Shinoda-DeVane-Klein (SDK) model^{27}—an existing coarse-grain model—properly captures PE properties and behaviors. The SDK model can be used to construct coarse-grain versions of molecules containing aliphatic, aromatic, polar, and ionic groups and can be used to simulate mixtures of these molecules.^{27–34} The SDK approach has been previously used to study complex phenomena such as lipid-nanoparticle interactions,^{35,36} membrane fusion,^{37} and dendrimersome formation.^{38} The SDK model continues to be both extended and developed under the *SPICA* name (e.g., see Ref. 39). However, the SDK model’s aliphatic coarse-grain beads were optimized to reproduce the properties of short alkanes in the liquid state at 303.15 K for the purpose of simulating surfactant systems.^{27} It was unclear whether the SDK model would recover the characteristics of long-chain PE systems and processes. On the basis of disparate static properties and dynamical phenomena, it is herein demonstrated that the SDK model is a suitable coarse-grain model for probing entangled PE melts and PE crystallization phenomenology. Furthermore, time scales on the order of tens of microseconds were accessible in this study thanks to the SDK model. These time scales combined with the broader capabilities of the SDK model afford new opportunities for probing and understanding PE processes and materials at the molecular level.

## II. METHODS

For reference, a brief overview of the SDK model^{27} is provided in this section with a focus on the elements relevant to modeling PE. More details about the SDK model and its subsequent development can be found in Refs. 27–34 and 39. In this study, PE chains were represented using the SDK CM and CT bead types,^{27} which correspond to –(CH_{2})_{3}– and –CH_{2}CH_{2}CH_{3}, respectively. These bead types were connected to form PE chains via harmonic bond potentials,

where *U*_{bond} is the change in potential energy associated with bond stretching, *k*_{bond} is the bond spring constant, *l* is the instantaneous bond length, and *l*_{0} is the equilibrium bond length. The bond parameters relevant to PE are provided in Table I. Each group of three consecutive beads along a PE chain backbone also interacted via a harmonic bond potential,

where *U*_{angle} is the potential energy change associated with the deviation of the bond angle (*θ*) from its equilibrium value (*θ*_{0}) and *k*_{angle} is the angle spring constant. Bond angle parameters relevant to PE are provided in Table II.

Bond type . | k_{bond} (kcal/mol)/Å^{2}
. | l_{0} (Å)
. |
---|---|---|

CT-CM | 6.160 | 3.65 |

CM-CM | 6.160 | 3.64 |

Bond type . | k_{bond} (kcal/mol)/Å^{2}
. | l_{0} (Å)
. |
---|---|---|

CT-CM | 6.160 | 3.65 |

CM-CM | 6.160 | 3.64 |

Angle type . | k_{angle} (kcal/mol)/rad^{2}
. | θ_{0} (deg)
. |
---|---|---|

CT-CM-CM | 1.190 | 175.0 |

CM-CM-CM | 1.190 | 173.0 |

Angle type . | k_{angle} (kcal/mol)/rad^{2}
. | θ_{0} (deg)
. |
---|---|---|

CT-CM-CM | 1.190 | 175.0 |

CM-CM-CM | 1.190 | 173.0 |

Beads belonging to different chains and beads belonging to the same chain but separated by more than 2 bonds interacted via pairwise 9-6 Lennard-Jones potentials of the form,

where *r* is the distance between a pair of beads, while the values of *ε* and *σ* depend on the types of beads involved in the pair interaction (see Table III). The SDK representation of PE does not include electrostatics.

Bead types . | ε (kcal/mol)
. | σ (Å)
. |
---|---|---|

CM CM | 0.4200 | 4.5060 |

CM CT | 0.4440 | 4.5455 |

CT CT | 0.4690 | 4.5850 |

Bead types . | ε (kcal/mol)
. | σ (Å)
. |
---|---|---|

CM CM | 0.4200 | 4.5060 |

CM CT | 0.4440 | 4.5455 |

CT CT | 0.4690 | 4.5850 |

The SDK model and PE were studied using Molecular Dynamics (MD) simulations. Three PE chain lengths were considered in this study: *n-*C_{480}H_{962}, *n-*C_{720}H_{1442}, and *n-*C_{960}H_{1922}. For each chain length, a low-density system composed of 400 chains was constructed and underwent an energy minimization. Each system was then heated from 10 K to 500 K under isochoric conditions over the course of 2 × 10^{6} time steps (effective time: ∼41 ns—see discussion later in this section for more details). Each system was pressurized by subjecting it to a NPT simulation at 500 K that was 10 × 10^{6} time steps long (∼204 ns) with a barostat set point of 1.0 atm. Each system was then simulated at 1 atm and 500 K under NPT conditions. This second NPT simulation was 100 × 10^{6} time steps long (∼2 *μ*s) for the *n-*C_{480}H_{962} system, 200 × 10^{6} (∼4 *μ*s) for the *n-*C_{720}H_{1442} system, and 400 × 10^{6} (∼8 *μ*s) for the *n-*C_{960}H_{1922} system. The systems were further equilibrated under NVT conditions at 500 K. These NVT equilibration simulations were 200 × 10^{6} time steps long (∼4 *μ*s) for the *n-*C_{480}H_{962} and *n-*C_{720}H_{1442} systems, and 400 × 10^{6} (∼8 *μ*s) for the *n-*C_{960}H_{1922} system. Long equilibration simulations were performed since the systems were expected to correspond to entangled polymer melts based on experimental results for PE (e.g., see Ref. 40). Variable simulation lengths were used since polymer melts with longer chains have slower dynamics.

As discussed in Sec. III B, *post hoc* analysis of the production simulations described below indicated that the above protocol yielded equilibrated PE melts. Moreover, the equilibration simulations are long with respect to the time scales associated with polymer chain rearrangement in this study’s systems as evidenced by comparing the simulation lengths reported here with the disengagement times discussed in Sec. III C. There are various algorithms for rapidly equilibrating polymer melts (e.g., see Ref. 41). However, such approaches require users to provide information about the target system (e.g., bond angles, bond lengths, and radial distributions as in the case of Ref. 41), and then the system is constructed according to the provided information. Given that the goal of this study was to test the ability of the SDK model to reproduce the properties of PE and given that the SDK model was parameterized^{27} using data at 303.15 K, equilibration strategies like those in Ref. 41 were not leveraged for the current study; future work could consider such approaches.

Following the equilibration simulations, a production NVT simulation at 500 K was performed for each system. The *n-*C_{480}H_{962} production simulation was 800 × 10^{6} time steps (∼16 *μ*s), the *n-*C_{720}H_{1442} simulation was 1.2 × 10^{9} time steps (∼24 *μ*s), and the *n*-C_{960}H_{1922} simulation was 2.0 × 10^{9} time steps (∼41 *μ*s). In general, long simulations were required to sample the long time scales associated with polymer dynamics and to sample the stress relaxation modulus discussed in Sec. III C.

For the production simulations, system configurations were saved every 100 000 time steps, as were individual chain properties (e.g., radius of gyration). To probe diffusive behaviors at the level of the individual SDK beads and on time scales of less than 100 000 time steps, average mean-squared displacements were evaluated on-the-fly and saved every 1000 time steps during the simulations. For longer time scales and for other diffusive behaviors (e.g., center-of-mass diffusion), mean-squared displacements were extracted from the saved configurations *post hoc*. The elements of the system’s symmetric pressure tensor were saved every time step. Additional system properties were saved every 10 000 time steps. Additional data analysis details are provided in Sec. III.

For the melting point estimate discussed in Sec. III D, a perfect crystal consisting of 3600 *n*-C_{171}H_{344} chains was constructed according to the experimental crystal structure for long-chain hydrocarbons.^{42} The crystal was heated from 1 K to 350 K over the course of 30 × 10^{6} time steps, while the crystal was allowed to relax anisotropically at 1 atm. The relaxed crystal was then heated from 350 K to 500 K at different heating rates (i.e., over the course of simulations ranging from 10 × 10^{6} to 30 × 10^{6} time steps in length). Three replicates were performed for each heating rate, and the simulations were performed under constant pressure conditions at 1 atm while allowing the systems to relax anisotropically. System configurations and properties were saved every 100 000 and 1000 time steps, respectively.

The crystallization simulation discussed in Sec. III D was conducted by quenching a configuration from the *n*-C_{720}H_{1442} production simulation. The configuration was first quenched to 300 K and simulated for 100 × 10^{6} time steps (∼2 *μ*s) under NPT conditions at 1 atm. Crystallization was not observed during this simulation. The *n*-C_{720}H_{1442} system was then quenched to 280 K and simulated under NPT conditions for 200 × 10^{6} time steps (∼4 *μ*s) during which time crystal nucleation occurred.

Dynamical processes occur more rapidly in simulations leveraging coarse-grain models than those using atomistic models as is discussed in Refs. 23, 43, and 44 and references therein. Therefore, the times associated with the SDK simulations should be scaled to provide an indication of anticipated time scales for experimental PE systems. This is why simulation lengths have been primarily reported in terms of time steps up to this point; the parenthetical times correspond to times scaled according to the 4.07 factor reported below. Pearson and co-workers previously performed a detailed analysis of self-diffusion coefficients at 448 K for linear PE samples with molecular masses ranging from 600 to 120 000 g/mol and provided an empirical relationship between self-diffusion coefficients and molecular masses.^{40} Therefore, an additional set of simulations was performed at 448 K to obtain a scaling factor for melt phase SDK PE chains. The *n*-C_{480}H_{962} system was quenched to 448 K and simulated for 40 × 10^{6} time steps under NPT conditions at 1 atm. A NVT production simulation was then performed at 448 K for 600 × 10^{6} time steps, and system configurations were saved every 100 000 time steps. The time-dependent average mean squared displacement for the center-of-mass positions of the polymer chains was determined using the configurations from the production simulation. For time windows exceeding 592 × 10^{6} time steps, the relative standard error associated with the mean squared displacement estimates exceeded 0.5%, and sampling issues started to arise; recall that the 448 K simulation was 600 × 10^{6} time steps in total length. A linear fit of the time evolution of the average center-of-mass mean squared displacement for the time window of 500 × 10^{6} to 592 × 10^{6} time steps yielded a diffusion coefficient of 1.765 × 10^{−7} ± 1 × 10^{−10} cm^{2}/s. Comparing to the value of 4.34 × 10^{−8} cm^{2}/s from the experimental empirical relationship of Pearson and co-workers,^{40} the SDK model yields PE chains that diffuse approximately 4.07-fold faster than their experimental counterparts. The 5 fs time step used in this study’s simulations corresponds to ∼20.4 fs with respect to experimental systems. Consequently and consistent with previous work,^{45} all times reported in this study are scaled times (i.e., include the 4.07 scaling factor), unless otherwise stated.

Coarse-grain temporal scaling factors are dependent on system conditions (e.g., temperature and pressure). The 4.07 scaling factor for this study’s simulations is based on a coarse-grain SDK *n*-C_{480}H_{962} melt at 448 K and a density of 0.76 g/cm^{3}. In turn, this study’s scaling factor might not be appropriate for PE systems at substantially different system conditions. However, work probing PE crystal nucleation with the SDK model suggests that it is reasonable to use the 4.07 scaling factor in the context of supercooled PE melts and PE crystallization.^{45}

One final simulation was performed to obtain the thermal expansion data discussed in Sec. III A. The final configuration from the 448 K NPT simulation of the *n*-C_{480}H_{962} melt was used to launch a cooling simulation. Over the course of 19 200 000 time steps (∼391 ns), the system was cooled from 448 K to 400 K, while the properties of the system were monitored (e.g., system volume). System properties were saved every 1000 time steps.

The following methodological details apply to all MD simulations in this study unless otherwise stated. All MD simulations were performed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS),^{46} in conjunction with periodic boundary conditions and a 15.0 Å cutoff for Lennard-Jones interactions. Particle positions and momenta (as well as simulation cell sizes in the case of NPT simulations) were evolved using the approach of Shinoda and co-workers,^{47} which involves using Nosé-Hoover^{48,49} chain thermostats.^{50} Temperature was controlled using four-member chain thermostats with coupling constants of 1 ps (unscaled time). For the NPT simulations, four-member chain barostats were used with coupling constants of 10 ps (unscaled time). All NPT simulations used isotropic barostats except for the *n*-C_{171}H_{344} simulations and the 280 K crystallization simulation; these simulations used anisotropic barostating. Anisotropic barostating was achieved by allowing the orthogonal directions of the simulation cells to relax independently while maintaining the same set point pressure for each direction. In contrast to the other simulations in this study, the simulations heating the *n*-C_{171}H_{344} crystal from 350 K to 500 K and the simulation cooling the *n*-C_{480}H_{962} melt from 448 K to 400 K all leveraged Berendsen^{51} thermostats (time constant: 1 ps unscaled time) and Berendsen barostats (time constant: 5 ps unscaled time and modulus: 23 000 atm), while particle equations of motion were integrated according to the velocity-Verlet approach.^{52} The thermal ramping was achieved through variation in the set points of the thermostats as internally implemented in LAMMPS.^{46}

## III. RESULTS AND DISCUSSION

The thermodynamic, structural, and dynamical properties of this study’s PE melt systems demonstrate that the SDK model affords a suitable coarse-grain representation of PE. Moreover, preliminary results reveal that the SDK model can be used to probe PE crystallization.

### A. Thermodynamic properties

The focus of this subsection is the thermodynamic properties of PE melts. The properties of the PE crystalline phase when using the SDK model will be discussed in Sec. III D. This study’s SDK PE melt at 448 K (M = 6,735 g/mol) exhibited a density of 0.758 g/cm^{3}, which agrees well with densities from previous experimental and computational studies on PE melts at comparable conditions (see Table IV). The thermal expansion coefficient (*α*_{P}) of SDK PE melts was estimated using the relationship,^{53}

and the volume-temperature data obtained from cooling the *n*-C_{480}H_{962} melt from 448 K to 400 K. The thermal expansion coefficient for SDK melts between 400 K and 448 K is approximately 6.54 × 10^{−4} K^{−1}. While this value is low compared to previous experimental results^{54,55} as provided in Table V, it is well within the range of thermal expansion coefficients from previous computational studies with other models (also see Table V). It is worth noting that the previous models^{14–16,19,56} also generally have difficulty in reproducing experimental thermal expansion coefficients, potentially highlighting a common challenge when constructing molecular models.

Source . | Density (g/cm^{3})
. | Molecular mass (g/mol) . | Melt conditions . |
---|---|---|---|

This work | 0.758 | 6 735 | 448 K, 1 atm |

Moyassari et al.^{19} (own CG model) | 0.783 | 17 957 | 450 K, 1 atm |

Sliozberg et al.^{41} (Paul et al. UA model^{16}) | 0.77/0.82 | 14 029 | 500 K/400 K, 1 atm |

Ramos et al.^{57} (Paul et al. UA model^{16}) | 0.795 | 2 695 | 500 K |

Ramos et al.^{57} (TraPPE-UA model^{56}) | 0.723/0.754 | 2 695 | 500 K/450 K |

Ramos et al.^{58} (TraPPE-UA model^{56}) | ∼0.772 | 14 015 (no branching) | 450 K, 1 atm |

Foteinopoulou et al.^{59} (Karayiannis et al. UA model^{14,15}) | ∼0.775 | Long chain limit | 450 K, 1 atm |

Dodd and Theodorou^{18} (own UA model) | ∼0.769 | 1 067 g/mol | 450 K, 0.1 MPa |

Foteinopoulou et al.^{60} (Karayiannis et al. UA model^{14,15}) | 0.783 | Long chain limit | 450 K, 1 atm |

Fukunaga et al.^{25} (own CG model) | ∼0.81/∼0.83 | Long chain limit | 500 K/450 K, 1 atm |

Ramos et al.^{61} | 0.756 | 2 300 000 and 3 600 000 | 463 K |

Pearson et al.^{40} | 0.761 | 3 600 | 448 K |

Pearson et al.^{40} | 0.766 | 12 400 | 448 K |

Source . | Density (g/cm^{3})
. | Molecular mass (g/mol) . | Melt conditions . |
---|---|---|---|

This work | 0.758 | 6 735 | 448 K, 1 atm |

Moyassari et al.^{19} (own CG model) | 0.783 | 17 957 | 450 K, 1 atm |

Sliozberg et al.^{41} (Paul et al. UA model^{16}) | 0.77/0.82 | 14 029 | 500 K/400 K, 1 atm |

Ramos et al.^{57} (Paul et al. UA model^{16}) | 0.795 | 2 695 | 500 K |

Ramos et al.^{57} (TraPPE-UA model^{56}) | 0.723/0.754 | 2 695 | 500 K/450 K |

Ramos et al.^{58} (TraPPE-UA model^{56}) | ∼0.772 | 14 015 (no branching) | 450 K, 1 atm |

Foteinopoulou et al.^{59} (Karayiannis et al. UA model^{14,15}) | ∼0.775 | Long chain limit | 450 K, 1 atm |

Dodd and Theodorou^{18} (own UA model) | ∼0.769 | 1 067 g/mol | 450 K, 0.1 MPa |

Foteinopoulou et al.^{60} (Karayiannis et al. UA model^{14,15}) | 0.783 | Long chain limit | 450 K, 1 atm |

Fukunaga et al.^{25} (own CG model) | ∼0.81/∼0.83 | Long chain limit | 500 K/450 K, 1 atm |

Ramos et al.^{61} | 0.756 | 2 300 000 and 3 600 000 | 463 K |

Pearson et al.^{40} | 0.761 | 3 600 | 448 K |

Pearson et al.^{40} | 0.766 | 12 400 | 448 K |

Source . | α_{P} (10^{−4} K^{−1})
. |
---|---|

This work | 6.54 |

Moyassari et al.^{19} (own CG model) | 6.21 |

Ramos et al.^{57} (UA model of Paul et al.^{16}) | 5.95 |

Ramos et al.^{57} (TraPPE-UA model^{56}) | 8.75 |

Foteinopoulou et al.^{60} (Karayiannis et al. UA model^{14,15}) | 6.9 |

Zhao et al.^{54} (high-density PE, 423 K) | 7.60 |

Zhao et al.^{54} (low-density PE, 423 K) | 7.93 |

Han et al.^{55} (high-density PE, 449 K) | 7.92 |

Han et al.^{55} (branched polyolefins, 394 K) | 6.95–7.35 |

Source . | α_{P} (10^{−4} K^{−1})
. |
---|---|

This work | 6.54 |

Moyassari et al.^{19} (own CG model) | 6.21 |

Ramos et al.^{57} (UA model of Paul et al.^{16}) | 5.95 |

Ramos et al.^{57} (TraPPE-UA model^{56}) | 8.75 |

Foteinopoulou et al.^{60} (Karayiannis et al. UA model^{14,15}) | 6.9 |

Zhao et al.^{54} (high-density PE, 423 K) | 7.60 |

Zhao et al.^{54} (low-density PE, 423 K) | 7.93 |

Han et al.^{55} (high-density PE, 449 K) | 7.92 |

Han et al.^{55} (branched polyolefins, 394 K) | 6.95–7.35 |

### B. Structural properties

This study’s PE melts also exhibit molecular-level structural properties comparable to previous experimental and computational studies as is demonstrated in this subsection. It is worth noting that it is not anticipated that the SDK model will capture the structural details of PE melts on length scales comparable to the size of the coarse-grain unit for the SDK model (i.e., three methylene units). Rather, this study focuses on higher level structural properties, such as characteristic ratio and mean squared radius of gyration.

Mean squared radius of gyration ($\u27e8Rg2\u27e9$) provides an indication of how compact polymer chains are within a polymer melt and can be calculated according to the following equations:

In Eq. (5), the sum is taken over the beads composing a polymer chain with the interior quadratic term being the squared distance between the position of the i-th bead of the chain (*r*_{i}) and position of the chain’s center of mass (*r*_{COM}). The latter is given by Eq. (6) where *m*_{i} is the molar mass of the i-th bead and *M* is the molar mass of the chain. The angular brackets in Eq. (5) indicate ensemble averaging, and this convention is followed throughout this study unless otherwise stated. The SDK model yields polymer melts with $\u27e8Rg2\u27e9$ values comparable to previous computational^{24,57,60–63} and experimental^{40,64} studies as is shown in Fig. 1(a). In particular, the $\u27e8Rg2\u27e9$ values for this study’s PE melts align well with the $\u27e8Rg2\u27e9$—molecular mass relationship established by Takahashi and co-workers.^{65}

A polymer chain’s characteristic ratio (*C*_{∞}) indicates the chain’s stiffness,^{69} and experimental work has previously determined characteristic ratios for PE samples.^{64,70,71} Computational studies^{15,41,67,68,72} have estimated characteristic ratios for polymer species by evaluating the *C*(*n*) function in Eq. (7) (or variants of this function) and then taking its infinite limit,

Here, *R*(*n*)^{2} represents the squared internal distance between a monomer and its n-th neighbor along the polymer chain backbone, while *l* is the separation between adjacent monomers (i.e., n = 1). It is possible to construct a coarse-grain analog of Eq. (7) by (1) replacing *R*(*n*)^{2} with the squared internal distance between a coarse-grain bead and its n-th coarse-grain neighbor along the polymer chain backbone and (2) replacing *l* with the separation between adjacent coarse-grain beads (i.e., *n*_{CG} = 1), thus arriving at

Here, the CG subscripts indicate that the quantities are for a coarse-grain representation. The *C*_{CG}(*n*_{CG}) function for this study’s *n*-C_{960}H_{1922} melt is provided in Fig. 1(b), and it displays expected asymptotic behavior. The infinite limit of Eq. (8) (*C*_{∞CG}) corresponds to the coarse-grain analog of *C*_{∞}. However, *C*_{∞CG} and *C*_{∞} are not equivalent, and the direct usage of *C*_{∞CG} can lead to an incorrect estimate of the characteristic ratio of the atomistic polymer species that a coarse-grain model represents as is now demonstrated.

As previously detailed by Auhl *et al.*,^{67} one can write ⟨*R*(*n*)^{2}⟩ according to

where *R*(*n*)^{2}, *n*, and *l* have the same meaning as in Eq. (7), while *θ* is the bond angle. In the limit of large *n*, the second term in the parentheses approaches zero implying that the first term corresponds to *C*_{∞},^{67} thus allowing one to write the following equation:

Given that *n*_{CG} = *n*/*λ*, where *λ* is the monomer scaling factor (i.e., the number of methylene units per coarse-grain bead in the case of SDK PE chains) and $\u27e8R(n)2\u27e9=\u27e8R(nCG)2\u27e9$ for a perfect coarse-grain model, Eq. (10) can be re-expressed as

In turn, Eq. (11) can be rewritten as follows:

Therefore, *C*_{∞CG} should be scaled by $\u27e8lCG\u27e92/(\lambda \u27e8l\u27e92)$ to obtain the characteristic ratio of the atomistic species that a CG model represents. The distinction between *C*_{∞CG} and *C*_{∞} will be particularly important when one is combining models of differing length scales (e.g., as part of a hierarchical backmapping strategy^{73}).

Fitting the *C*_{CG}(*n*_{CG}) function for this study’s *n*-C_{960}H_{1922} melt according to $CCG(nCG)=C\u221eCG+a1/nCG+a2/nCG2+a3/nCG3$ (*R*^{2} > 0.99) yielded a *C*_{∞CG} value of 4.761 ± 0.001. *C*_{∞} values have been previously estimated using such a fitting procedure.^{15,41} As shown in Table VI, 4.761 is much lower than *C*_{∞} values from experiments^{64,70,71} and previous atomistic simulations.^{15,41,60,62} This is consistent with the fact that the SDK model has less rigid coarse-grain bond angles compared to atomistic PE (e.g., the CHARMM forcefield^{74} that was used to parameterize the SDK model). For the SDK PE melts, $\u27e8lCG\u27e92/(\lambda \u27e8l\u27e92)=1.866$, so the SDK PE model corresponds to an atomistic polymer with *C*_{∞} = 8.884 ± 0.002, which is slightly elevated, but still reasonable for PE (see Table VI).

Source . | C_{∞}
. | Temperature (K) . |
---|---|---|

This work (C_{∞CG}) | 4.761 | 500 |

This work (C_{∞}) | 8.884 | 500 |

Sgouros et al.^{21} (own mesoscopic model, C_{2080}H_{4162}) | 8.93 | 450 |

Sliozberg et al.^{41} (Paul et al.UA model^{16}) | 7.14 | 500 |

Ramos et al.^{57} (Paul et al.UA model^{16}) | 7.61 | 500 |

Ramos et al.^{57} (TraPPE-UA model^{56}) | 8.18 | 500 |

Carbone et al.^{62} (Smit et al.UA model^{17}) | 8.27 | 450 |

Foteinopoulou et al.^{60} (Karayiannis et al.UA model^{14,15}) | 7.868 | 500 |

Ramos et al.^{61} (TraPPE-UA model^{56}) | 8.34 | 450 |

Ramos et al.^{61} (Karayiannis et al.UA model^{14,15}) | 8.06 | 450 |

Ramos et al.^{58} (TraPPE-UA model,^{56} no branching) | 8.4 | 450 |

Karayiannis et al.^{15} (own UA model) | 8.27 | 450 |

Harmandaris et al.^{63} (own model, building off Ref. 18) | 8.35 | 450 |

Fetters et al.^{71} | 8.3 | 298 |

Fetters et al.^{70} | 7.3 | 413 |

Horton et al.^{64} [using Eq. (11) in Ref. 64] | 7.4 | 413 |

Source . | C_{∞}
. | Temperature (K) . |
---|---|---|

This work (C_{∞CG}) | 4.761 | 500 |

This work (C_{∞}) | 8.884 | 500 |

Sgouros et al.^{21} (own mesoscopic model, C_{2080}H_{4162}) | 8.93 | 450 |

Sliozberg et al.^{41} (Paul et al.UA model^{16}) | 7.14 | 500 |

Ramos et al.^{57} (Paul et al.UA model^{16}) | 7.61 | 500 |

Ramos et al.^{57} (TraPPE-UA model^{56}) | 8.18 | 500 |

Carbone et al.^{62} (Smit et al.UA model^{17}) | 8.27 | 450 |

Foteinopoulou et al.^{60} (Karayiannis et al.UA model^{14,15}) | 7.868 | 500 |

Ramos et al.^{61} (TraPPE-UA model^{56}) | 8.34 | 450 |

Ramos et al.^{61} (Karayiannis et al.UA model^{14,15}) | 8.06 | 450 |

Ramos et al.^{58} (TraPPE-UA model,^{56} no branching) | 8.4 | 450 |

Karayiannis et al.^{15} (own UA model) | 8.27 | 450 |

Harmandaris et al.^{63} (own model, building off Ref. 18) | 8.35 | 450 |

Fetters et al.^{71} | 8.3 | 298 |

Fetters et al.^{70} | 7.3 | 413 |

Horton et al.^{64} [using Eq. (11) in Ref. 64] | 7.4 | 413 |

An alternative metric for assessing the stiffness of a polymer chain is its persistence length. The bond-bond correlation function [*C*_{bb}(*s*)] can be related to a chain’s persistence length according to^{75}

where the numerator corresponds to the dot product of the i-th monomer’s bond vector (*b*_{i} = *r*_{i+1} − *r*_{i}) with the bond vector of the *i* + *s* monomer (*b*_{i+s} = *r*_{i+s+1} − *r*_{i+s}) and *l*_{p} is the persistence length of the polymer chain. An exponential fit of the bond-bond correlation function for the *n*-C_{960}H_{1922} melt [see Fig. 1(c)] yields a persistence length of 0.825 ± 0.001 nm, which is slightly elevated for PE but still within the range of previously observed values for PE systems (see Table VII).

Source . | Persistence length (nm) . | Temperature (K) . |
---|---|---|

This work | 0.825 ± 0.001 | 500 |

Salerno and Bernstein (own A1 CG model)^{20} | 0.835 | 500 |

Yamamoto^{76} (earlier Yamamoto model^{77}) | 0.50 | 405 |

Foteinopoulou et al.^{60} (Karayiannis et al. UA model^{14,15}) | 0.682 | 500 |

Ramos et al.^{61} (TraPPE-UA model^{56}) | 0.922 | 450 |

Ramos et al.^{61}(Karayiannis et al. UA model^{14,15}) | 0.954 | 450 |

Ramos et al.^{58} (TraPPE-UA model,^{56} no branching) | 0.91 | 450 |

Ramachandran et al.^{78} | 0.65–0.90 | 398.15 |

Source . | Persistence length (nm) . | Temperature (K) . |
---|---|---|

This work | 0.825 ± 0.001 | 500 |

Salerno and Bernstein (own A1 CG model)^{20} | 0.835 | 500 |

Yamamoto^{76} (earlier Yamamoto model^{77}) | 0.50 | 405 |

Foteinopoulou et al.^{60} (Karayiannis et al. UA model^{14,15}) | 0.682 | 500 |

Ramos et al.^{61} (TraPPE-UA model^{56}) | 0.922 | 450 |

Ramos et al.^{61}(Karayiannis et al. UA model^{14,15}) | 0.954 | 450 |

Ramos et al.^{58} (TraPPE-UA model,^{56} no branching) | 0.91 | 450 |

Ramachandran et al.^{78} | 0.65–0.90 | 398.15 |

Characteristic ratio and persistence length values suggest that SDK PE chains have slightly elevated rigidity. For applications where chain rigidity is especially important in and of itself (e.g., hierarchical backmapping^{73}), future work might consider modifying the SDK parameters to reduce the rigidity of SDK PE chains, though such reparameterizations will likely be challenging given the substantial work that initially went into developing the SDK model. Moreover, the SDK model does recover PE properties that are dependent on chain rigidity as will be demonstrated in Secs. III C and III D (e.g., entanglement mass as shown in Sec. III C), and SDK PE chains crystallize via chain-folded structures to form semicrystalline lamellae that are of comparable thickness to previous experimental work (see Sec III D). Attempts to reduce the rigidity of SDK PE chains could negate these other agreements and would necessarily influence the broader properties of the model. Therefore, despite the fact that SDK PE chains exhibit slightly increased rigidity, it is still reasonable to study PE with the original SDK model.

For the SDK model, the PE chain lengths considered in this study are structurally consistent with Gaussian chains. For Gaussian chains,^{69,79} $\u27e8Ree2\u27e9/(6\u27e8Rg2\u27e9)=1$, with $Ree2$ being the mean squared end-to-end distance of the polymer chains. The evaluation of the left-hand side of the equality has been previously used to assess whether PE systems corresponded to Gaussian chains with values close to one being interpreted as Gaussian behavior.^{22} This study’s systems exhibit $\u27e8Ree2\u27e9/(6\u27e8Rg2\u27e9)$ values in the range of 1.01–1.02, which is consistent with previous work on PE chains of comparable length.^{22} The SDK model enables the simulation of large-scale PE melts with Gaussian behavior and structural properties that align well with previous work on PE.

### C. Entangled dynamics

The SDK model recovers expected dynamical behaviors for entangled PE melts. As summarized by Doi and Edwards^{69} as well as by Kremer and Grest,^{79} the time evolution of the average mean squared displacement (MSD) of polymer monomers [*g*_{1}(*t*)] transitions through several temporal scaling regimes for entangled polymer melts. Calculation of *g*_{1}(*t*) according to Eq. (14) reveals that the SDK PE melts exhibit such behavior [see Fig. 2(a)]. Additional MSD functions can also be defined for a polymer chain’s monomers with respect to the chain’s center of mass [see Eq. (15)] and for polymer chain’s center of mass [see Eq. (16)]. The behavior of these additional MSD functions also indicates that this study’s PE systems correspond to entangled polymer melts. For example, as shown in Fig. 2(b), the center-of-mass MSD function [*g*_{3}(*t*)] exhibits appropriate temporal scaling regimes while *g*_{2}(*t*) crosses *g*_{3}(*t*) at long times,

The slow, complex dynamics of polymer melts are highlighted in Fig. 2 by the microsecond time scales on which MSD temporal scaling is less than *t*^{1} (i.e., subdiffusive). The microsecond time scales associated with polymer dynamics in this study’s PE melts are exemplified by the slow decay of their end-to-end vector correlation functions [see Fig. 3(a)]. As de Gennes demonstrated,^{80} the end-to-end vector correlation function decays exponentially according to a polymer’s disengagement time (*τ*_{d}), the longest time scale associated with a polymer chain’s conformational rearrangement, and the time scale on which a polymer chain escapes the tube-like region of space to which it is initially confined by surrounding obstacles (e.g., other polymer chains). Exponential fits of the curves in Fig. 3(a) (*R*^{2} > 0.99) yielded disengagement times of 0.407 ± 0.001 *μ*s, 1.591 ± 0.003 *μ*s, and 4.262 ± 0.005 *μ*s for the *n*-C_{480}H_{962}, *n*-C_{720}H_{1442}, and *n*-C_{960}H_{1922} systems, respectively. These *τ*_{d} values demonstrate *M*^{3.4} scaling as shown in Fig. 3(b). As highlighted by Doi and Edwards’ review of polymer theories,^{69} theoretical treatments of polymer dynamics generally conclude that *τ*_{d} values are proportional to zero shear viscosities (*η*_{0}), and hence comparable scaling behavior is anticipated for both *τ*_{d} and *η*_{0} with respect to molecular mass. Consistent with these conclusions, the observed *τ*_{d} scaling factor of 3.4 agrees with previous experimental^{81,82} and computational^{23} results showing that *η*_{0} scales according to *M*^{3.4} for entangled long-chain polymer melts. Although some studies have observed slightly different scaling behaviors (e.g., *M*^{3.64} in the case of Ref. 40), it is broadly accepted that polymer melt viscosities empirically scale according to *M*^{3.4} for long polymer chains.^{69} Therefore, SDK PE chains exhibit known phenomenological behavior and relationships for entangled polymer melts.

The SDK model also recovers entanglement properties specific to PE systems, such as the PE entanglement mass (*M*_{e}). Consistent with previous *in silico* studies,^{22,23} the stress relaxation modulus [*G*(*t*)] for each PE system was calculated according to the Green-Kubo relation in the following equation:

where the ensemble average in the numerator is the average of the un-normalized autocorrelation functions for the off-diagonal elements of the system’s pressure tensor (*σ*) for *αβ* = *xy*, *xz*, and *yz*. For an entangled polymer system, *G*(*t*) exhibits a plateau regime at increasing molecular masses (see Fig. 4). The pressure corresponding to the plateau (*G*_{plateau}) is connected to the entanglement mass of a polymer species according to the following equation:^{23,69}

Here, *M*_{e} is the entanglement mass and *R* is the universal gas constant. Based on this relationship and a *G*_{plateau} value of ∼2.3 MPa for this study’s *n*-C_{920}H_{1922} melt, the SDK model yields PE chains with an entanglement mass of ∼1300 g/mol. There is substantial variation in both the computational and experimental estimates of *M*_{e} for PE (see Table VIII). Computational estimates range from ∼500–1700 g/mol, while experimental estimates range from 828 to 2400 g/mol. The SDK model yields a PE entanglement mass that is well within these ranges and toward their center.

Source . | M_{e} (g/mol)
. |
---|---|

This work | ∼1300 |

Moyassari et al.^{19} (own coarse-grain model) | ∼492 |

Salerno et al.^{23} (own coarse-grain models) | 1360–1425 |

Salerno et al.^{23} (Martini coarse-grain model^{83}) | 1360–1425 |

Sliozberg et al.^{41} (Paul et al. UA model^{16}) | 1390 |

Sgouros et al.^{21} (own mesoscopic model) | 766 |

Foteinopoulou et al.^{60} (Karayiannis et al. UA model^{14,15} at 500 K) | ∼1270–1400 |

Ramos et al.^{61} (Karayiannis et al. UA model^{14,15}) | 790 |

Ramos et al.^{58} (TraPPE-UA model,^{56} no branching) | 710 |

Foteinopoulou et al.^{59} (Karayiannis et al. UA model^{14,15}) | 855 |

Padding and Briels^{24} (coarse-grain blob model^{66}) | 1700 |

Padding and Briels^{24} (coarse-grain blob model,^{66} alternate) | 960 |

Litvinov et al.^{84} | 1760 |

Ramos et al.^{61} | 1140 |

Aguilar et al.^{85} | 2400 |

Fetters et al.^{86} | 1150 |

Fetters et al.^{70} | 828 |

Raju et al.^{87} | 1850 |

Source . | M_{e} (g/mol)
. |
---|---|

This work | ∼1300 |

Moyassari et al.^{19} (own coarse-grain model) | ∼492 |

Salerno et al.^{23} (own coarse-grain models) | 1360–1425 |

Salerno et al.^{23} (Martini coarse-grain model^{83}) | 1360–1425 |

Sliozberg et al.^{41} (Paul et al. UA model^{16}) | 1390 |

Sgouros et al.^{21} (own mesoscopic model) | 766 |

Foteinopoulou et al.^{60} (Karayiannis et al. UA model^{14,15} at 500 K) | ∼1270–1400 |

Ramos et al.^{61} (Karayiannis et al. UA model^{14,15}) | 790 |

Ramos et al.^{58} (TraPPE-UA model,^{56} no branching) | 710 |

Foteinopoulou et al.^{59} (Karayiannis et al. UA model^{14,15}) | 855 |

Padding and Briels^{24} (coarse-grain blob model^{66}) | 1700 |

Padding and Briels^{24} (coarse-grain blob model,^{66} alternate) | 960 |

Litvinov et al.^{84} | 1760 |

Ramos et al.^{61} | 1140 |

Aguilar et al.^{85} | 2400 |

Fetters et al.^{86} | 1150 |

Fetters et al.^{70} | 828 |

Raju et al.^{87} | 1850 |

### D. Melting and crystallization

The SDK model captures important properties and behaviors related to PE solid-liquid transitions, specifically melting and crystallization. To estimate the melting point (*T*_{m}) of SDK PE chains, *n*-C_{171}H_{344} crystals composed of extended chains were heated from 350 K to 500 K at different heating rates ranging from ∼0.25 to 0.74 K/ns. The disordering of each crystal’s constituent chains upon heating was assessed using the nematic order parameter described by Eppenga and Frenkel.^{88} More specifically, the second rank order parameter tensor (⟨*Q*⟩) was evaluated for each simulation configuration according to

where the sum is taken over the *N* coarse-grain beads in the configuration, *u*_{i}*u*_{i} is the dyadic product of the unit backbone orientation vector at bead *i* (*u*_{i}) with itself, and *I* is a second-rank identity tensor. For each nonterminal bead, *u*_{i} was estimated by normalizing the vector connecting the bead’s two intramolecular neighboring beads. For terminal beads, normalized bond vectors were used instead. The nematic order parameter (S) corresponds to the largest eigenvalue of < *Q* > such that a value of 1 corresponds to perfect alignment. Isotropic environments have S values approaching 0.^{88}

For each simulation, the melting of the *n*-C_{171}H_{344} crystal coincided with sharp decreases in nematic order and system density, and a sharp increase in average molecular potential energy [see Fig. 5(a)–5(b)]. While all three transitions are essentially coincident, they are not equivalent. Molecular potential energy is largely proportional to density as evidenced by their near linear relationship during the heating process [Fig. 5(c)]. The deviations from the linear fit, particularly, in the vicinity of the phase transition around ∼0.8 g/cm^{3} in Fig. 5(c), suggest that potential energy and density are not completely coupled, although they are strongly correlated. In contrast, nematic order and density do not exhibit a linear relationship during the heating simulations [Fig. 5(c)], indicating a complex relationship between orientational ordering and density. These results are consistent with previous studies on polymer crystallization, which distinguish densification and the development of orientational ordering.^{7,8,10,45,89} Moreover, the distinction between densification and orientational ordering is observed during the early-stages of crystallization in systems beyond polymers (e.g., see Ref. 90 and references therein).

The melting temperature of each heating simulation was estimated using its nematic order-temperature profile [i.e., S-T curves like those shown in Figs. 5(a)–5(b)]. The nematic order curves were used for two reasons. First, previous studies^{6,9,45,76,91–93} have extensively studied polymer crystallization using the *P*_{2} order parameter, which is analogous to S but suited to probing local environments. Second, the S-T curves could be fit (*R*^{2} > 0.99) to sigmoid functions of the form,

with *y*_{0}, *x*_{0}, *a*, and *b* serving as adjustable parameters. In turn, the temperature corresponding to the inflection point of the sigmoid fit for each simulation was taken to be the melting temperature of that simulation. The extracted melting temperatures exhibited a linear relationship with respect to heating rate [see Fig. 5(d)], and extrapolation to the 0 K/ns regime yielded a melting point estimate of 448.5 ± 0.5 K for SDK *n*-C_{171}H_{344} chains. The aforementioned strategy was used to estimate *T*_{m} since direct phase coexistence simulations at low supercoolings proved to be too time consuming; experimental PE crystal growth rates can be on the order of nanometers per second, if not smaller, at small supercoolings.^{94,95}

As noted by Waheed *et al.*,^{96} estimating melting points *in silico* by superheating crystals leads to an overestimation of *T*_{m} because a nucleation event is required to initiate the melting process. In particular, Waheed *et al.*^{96} obtained a *T*_{m} estimate for *n*-C_{20}H_{42} that was 35 K too high when using a superheating approach in conjunction with the atomistic alkane model of Paul *et al.*^{16} Subsequent work with the same atomistic model obtained a melting point for *n*-C_{50}H_{102} that agreed with experimental results,^{97} highlighting that the overestimation of *T*_{m} for *n*-C_{20}H_{42} likely arose from the melting point determination strategy leveraged by Waheed *et al.* rather than their atomistic model. The melting point estimate extracted for the SDK representation of *n*-C_{171}H_{344} is higher than the expected ∼400 K value based on experimental results.^{4,40} However, the discrepancy in the *n*-C_{171}H_{344} melting point corresponds to a relative error of ∼13%, which is comparable to the 11% relative error that Waheed *et al.* encountered for *n*-C_{20}H_{42}. We anticipate that the exact *in silico T*_{m} values for PE chains are reasonably close to experimental values when using the SDK model.

Since the pioneering work of Keller and co-workers that demonstrated that PE crystallizes via lamellae composed of chain folded structures (e.g., see Ref. 98), there has been much interest in chain folding during PE crystallization and PE lamellae (e.g., Refs. 4, 7–9, and 99). Therefore, to test the crystallization behavior of SDK PE chains, a *n*-C_{720}H_{1442} melt was quenched to 280 K. Local PE crystallinity was assessed using the *P*_{2} order parameter; this choice is also consistent with previous studies probing PE crystallization.^{6,9,45,76,91–93} More specifically and consistent with our previous work on PE,^{45} the *P*_{2} values for the system’s constituent beads were evaluated according to *P*_{2}(*i*) = ⟨(3cos^{2}Θ_{i,j} − 1)/2⟩, where Θ_{i,j} is the angle between the polymer backbones at the i-th and j-th beads, and the angular brackets indicate a local average over all beads within 6.35 Å of the i-th bead excluding itself. Beads were assigned to be crystalline if their *P*_{2} values were greater than or equal to 0.85. Starting from an initial metastable melt at 280 K with less than 0.25% of the beads having *P*_{2} ≥ 0.85, the *n*-C_{720}H_{1442} system underwent crystallization. SDK PE chains crystallize via the formation of lamellae [see Fig. 6(a)] that are composed of folded chains, as illustrated by Fig. 6(b).

To estimate the thickness of the lamella in Fig. 6(a), profiles of the crystalline bead population and average *P*_{2} values were calculated across the lamella along the lamella normal [Fig. 6(c)]. The positive and negative sides of each profile were fit to sigmoid functions of the form in Eq. (20) (*R*^{2} > 0.99), and the lamella thickness was taken to be the distance between the inflection points of the two sides of the profile. Accordingly, the crystalline bead population and average *P*_{2} profiles yielded lamella thicknesses of 7.75 ± 0.05 nm and 8.57 ± 0.05 nm. As such, the lamella shown in Fig. 6(a) is ∼8–9 nm thick, which agrees well with experiments.^{98,99} For example, Keller and O’Connor^{98} obtained a lamellar thickness of 9.2 nm by rapidly quenching a PE-xylene solution to 283 K, while the revised supercooling-initial fold length relationship of Barham *et al.*^{99} predicts an initial fold length of ∼7.6 nm during crystallization at 280 K. The SDK model thus yields PE crystals with appropriate lamellar morphology.

An important feature of a forcefield is the crystal phase that it yields upon crystallization. In general, previous models for probing PE crystallization, such as the modified UA model of Paul *et al.*,^{16,93,96,97,100} favor the formation of a hexagonal crystal structure rather than the orthorhombic phase—the crystal phase commonly associated with PE.^{42} An exception to this is the recent model of Zubova *et al.*,^{101,102} which was specifically designed with the orthorhombic phase of polyethylene in mind; the broader applicability and properties of their model remain to be determined (e.g., with respect to high-temperature melts and molecular-level properties such as $Rg2$). To assess the PE crystal structure favored by the SDK model, a short simulation was launched from the configuration corresponding to Fig. 6(a) in order to generate an average configuration and thus localize crystalline beads at their lattice positions. The simulation was 40 000 time steps in length and was identical to the production simulation that yielded the lamella-containing starting configuration, except for using a slightly shorter time step (i.e., 10.18 fs scaled time). Particle positions were averaged over the last 20 000 time steps of the simulation with a sampling frequency of 200 time steps. The average positions of crystalline beads composing the lamella exhibit a hexagonal structure as shown in Fig. 7(a). The *a* lattice parameter for this structure is approximately 5.0 Å, as evidenced by the second peak in the radial distribution function for the lamellar beads [see Fig. 7(b)], and as confirmed by visual analysis. It is important to note that PE adopts hexagonal crystal packing under certain conditions,^{103,104} and alkane chains can also exhibit pseudohexagonal phases, such as upon heating and just prior to melting.^{105–107} The extracted *a* value for the SDK hexagonal structure agrees well with the experimental value of approximately 5 Å for ultradrawn hexagonal PE at high temperatures.^{103} As such, even though the SDK model does not exhibit the more common orthorhombic phase of PE, it does still crystallize via a known PE crystal structure. Moreover, previous work^{10,108} suggests that hexagonal PE structuring might play an intermediary role during the formation of the orthorhombic phase. Consequently, it is reasonable to use the SDK model to probe PE crystallization.

Density is another important characteristic of a crystal phase. To estimate the density of the SDK PE hexagonal crystalline phase, a Voronoi analysis was performed on the configuration corresponding to the lamella pictured in Fig. 6(a), and Voronoi volumes were extracted for the crystalline beads composing the lamella. Based on these local volumes, the SDK PE crystalline phase has an estimated density of 0.93 g/cm^{3}. For comparison, Hsieh and Hu^{109} estimated a density of 1.02 g/cm^{3} for the crystalline portions of ultrahigh molecular weight PE fibers at room temperature. As such, the SDK model underestimates the PE crystal density by approximately 9%. Based on the experimental work of Tsubakihara *et al.*,^{103} there is an increase of approximately 10%–20% in molecular cross-sectional area perpendicular to the PE chain axis when ultradrawn PE transitions from the orthorhombic phase to the hexagonal phase at high temperatures. In turn, the ∼9% density discrepancy for the SDK PE crystal phase is consistent with SDK PE chains crystallizing via a hexagonal structure rather than the orthorhombic structure of PE. Nevertheless, the SDK model can be used to model both entangled PE melts and PE crystallization, which is surprising given that the SDK model was designed to reproduce small-molecule properties at 303.15 K. In contrast, it remains to be determined which, if any, of the other coarse-grain models are appropriate for probing PE crystallization, although there are some indications that it may be possible to leverage these models to study crystallization. For example, Salerno *et al.*^{23} observed sharp density increases upon cooling coarse-grain *n*-C_{96}H_{194} systems, while Fukunaga *et al.*^{25} observed the formation of lamellar crystalline structures when cooling a melt system composed of a single C_{600}H_{1202} chain to 300 K with their model. However, the crystallization details of these other models have yet to be determined, such as favored crystal structures and PE melting points. Moyassari *et al.*^{19} did leverage their coarse-grain models to probe crystallization and did establish that their low-temperature model yields a hexagonal structure; however, the broader utilization of their models might be problematic. For example, leveraging their models requires switching between models with temperature variations, and their models exhibit a low entanglement mass (see Table VIII). By providing crystallization and crystal phase information for the SDK model, a coarse-grain model that reasonably captures multiple facets of PE, this study has laid a foundation that the chemical physics community can leverage to probe PE crystallization with coarse-grain simulations and thus realize insights beyond the time and length scales of atomistic simulations. Moving forward, the scientific community can leverage the SDK model—a comparatively simple coarse-grain model—to study PE over extended time and length scales to realize physically accurate, bottom-up understandings of PE crystallization, PE materials, and PE-based composites.

## IV. CONCLUSIONS

The SDK model of PE was shown to reasonably capture the thermodynamic properties of PE melt and crystalline phases, molecular-level PE structure (e.g., the molecular mass dependence of $\u27e8Rg2\u27e9$ values), and dynamical characteristics of PE (e.g., PE entanglement mass). Furthermore, the estimated melting point obtained in this study suggests that SDK PE chains exhibit a melting point in the vicinity of experimental melting points for linear PE chains. As such, the SDK model is a suitable coarse-grain model for probing PE, particularly long-chain entangled melts and crystallization processes. Furthermore, SDK representations have been developed for many species beyond PE, and “off-the-shelf ” interaction parameter sets are already available for combining these species with the SDK bead types for PE, thus affording opportunities to explore PE in environments and processes beyond melts and melt-based crystallization (e.g., potentially PE crystallization from solution phases or in the presence of additives). The SDK model also renders longer time and length scales accessible compared to atomistic models, and SDK-based simulations can be conducted using LAMMPS (an open source MD code).^{46} The SDK model opens up new opportunities for exploring PE and PE-composites at the molecular level.

## ACKNOWLEDGMENTS

The authors thank the U.S. Army Research Laboratory for its support (Contract Nos. W911NF-16-2-0189 and W911NF-18-S-0269). This study was conducted using Temple University’s high-performance computing (HPC) resources and was thus also partially supported by the National Science Foundation through major research instrumentation Grant No. 1625061. This work was completed in part while K. Wm. Hall was an International Research Fellow of the Japan Society for the Promotion of Science (JSPS); K. Wm. Hall thanks JSPS for their support. M. L. Klein thanks H. R. H. Sheikh Saud for his support via a Sheikh Saqr Research Fellowship.