KDE Based Coarse–graining of Semicrystalline Systems with Correlated Three–body Intramolecular Interaction

Coarse–grained (CG) molecular dynamics (MD) simulations are commonly used to study polymer physics for their ability in probing molecular processes at substantially longer time and length scales than all–atom models.^{1, 2, 3, 4} While including more atoms in a bead reduces the computational cost of the model, the reduced resolution weakens its ability to accurately represent the structural properties of the underlying polymer. For example, Salero et al.⁵ found that polyethylene (PE) models with flexible chains did not crystallize at all when the coarse–grained beads represent more than four methylene groups. They attributed the phenomenon to the exceedingly large equilibrium bond–length compared with the equilibrium pair distance.^{6, 7} At the lowest degree of coarse–graining, united–atom model (UAM) that lumps hydrogen atoms with their bonded carbon atoms along the chain backbone are more capable at reproducing structural features. The UAM of polyethylene, initially developed by Paul et al.⁸ and subsequently modified by Waheed et al.,^{9, 10} yields a hexagonal crystal structure.¹¹ The UAM developed by Zubova et al. reproduced the expected orthorhombic crystal structure by adjusting the displacements of beads from the chain axes and the van der Waals radii.^{12, 13} However, given that united atom models have a relatively small reduction in the degrees of freedom compared to all–atom models, the range of time–scales the models can access remains highly limited. Therefore, many MD studies of polyethylene crystallization^{9, 14, 15, 16, 17, 18, 19} that employed UAMs used systems with chain lengths no longer than 1000 methyl units because the drastically reduced crystallization rate with longer chains necessitates unfeasible simulation time, although real polyethylene chains are much longer. Thus, developing CG models with efficient computation and accurate representation of structure can strengthen the capability of MD simulations in the investigations that require large systems and long time scales.

It has been shown that the accuracy of CG models in representing water,²⁰ RNA,²¹ and proteins^{22, 23} increase when considering multi–body interactions. Larini et al.²⁰ introduced three–body non–bonded interactions into their single–site water CG model by adding the relative angular orientation of the triplets of water molecules into the energy formulation; the addition greatly improved the model’s ability to reproduce the radial distribution function (RDF) and angular distribution function (ADF). Ejtehadi et al.²⁴ and Krishna et al.²⁵ showed that three–body interactions play an important role in accurately representing the protein folding mechanism and the structure of HIV-1 capsid, respectively. Also, Schaposchnikov et al.²⁶ showed that the contribution of three–body effect is 20% to 40% of ligand–capped nanocrystals. Multibody interactions can also improve the representability of polymer structures. Fukunaga et al.²⁷ presented a CG model that groups three CH₂ units into a CG bead and showed that the bond–lengths and torsion angles distributions strongly correlated with the bond–angle distributions at room temperature.

CG model interactions, e.g. non–bonded interactions, can be described either by a functional form or interpolated from tabulated values. Although, the tabulated potentials are more flexible in describing complex structures, they are also more prone to sampling noise. The noise–induced energetic fluctuations causes magnified oscillations when computing the forces, leading to a reduced stable time step. When training CG models with a structure–based method such as iterative Boltzmann inversion (IBI),²⁸ the structural distribution functions, e.g. RDF, are commonly constructed using histograms, which has well–known disadvantages such as limited smoothness with a small sample size and sensitivity to the bin size. Even though smoothing methods such as B-spline²⁷ and cubic spline²⁰ can be used to smooth the potential surface parametrized using histograms, McCabe et al.²⁹ showed that with the kernel density estimation (KDE), we can readily obtain analytically differentiable and smooth structural distribution functions regardless of sample size.

In this paper, we develop an approach for constructing CG models with three–body interactions to reproduce correlated bond–length and bond–angle distribution functions. The resulting potentials are generated on a two–dimensional grid and evaluated using bicubic interpolation. To generate accurate and smooth energy derivatives, the KDE is used to calculate the structural distribution functions that enables analytical differentiation of energy. We demonstrate the new method by developing coarse–grained models for PE, a model semicrystalline polymer. The models are named as CG–BA model. We also study how the KDE bandwidths affect the activation energy of diffusion, computational performance, and the crystallization kinetics of the resulting CG models.

To determine appropriate baseline KDE bandwidth parameters, we first compare the representation and sampling errors of the target structural distributions sampled from atomistic simulations. The representation error measures the loss of detail incurred by increasing the KDE bandwidth and is computed as the $L_{2}$ norm of the difference between the distributions calculated at a particular bandwidth and the reference bandwidth. The sampling error measures the average deviation of the distributions sampled across different systems. The definitions of these errors for the radial distribution functions are:

where $g_{w_{r}}$ is the RDF calculated using bandwidth $w_{r}$, $g_{0}$ is the RDF calculated at the reference bandwidth $w_{r_{0}}=0.01$ Å, and $\sigma_{w_{r}}$ is the point–wise standard deviation of $g_{w_{r}}$ calculated from the 15 different target atomistic systems. The representation and sampling errors of the bond–angle probability distributions are calculated as:

where $H$ represents the bandwidth matrix, defined by $w_{l}$ and $w_{\theta}$, and $\sigma_{H}$ denotes the standard deviation of the probability distributions computed using a given bandwidth matrix. The bond–length and bond–angle bandwidths used for the reference bandwidth matrix are $w_{l_{0}}=0.002$ Å  and $w_{\theta_{0}}=0.0031~{}\mathrm{rad}$, respectively. As shown in Figure 3, increasing the kernel bandwidth parameters reduces the sampling errors but results in substantially larger increases in the representation errors. The baseline bandwidth parameters were chosen where the representation and sampling errors have equal magnitudes: $w_{r}^{*}=0.07$ Å  $w_{l}^{*}=0.016$ Å, and $w_{\theta}^{*}=0.021~{}\mathrm{rad}$.

Subsequently, we conducted IBI training varying different bandwidths individually in the target distributions and studied the effect of each bandwidth on the resulting potentials. For $w_{r}$ study, we had 4 trainings using different $w_{r}$ with constant $w_{l}$ and $w_{\theta}$, e.g., $(nw^{*}_{r},w^{*}_{l},w^{*}_{\theta})$ where n = 1, 2, 4, and 8. Studies of $w_{l}$ and $w_{\theta}$ were conducted similarly. Therefore, there were 10 training conducted in total since $(w^{*}_{r},w^{*}_{l},w^{*}_{\theta})$ is shared in all three studies.

It is important to note that we used the baseline bandwidth set $(w^{*}_{r},w^{*}_{l},w^{*}_{\theta})$ in the CG simulations of all the trainings while changing the bandwidths of target distributions. If the bandwidths used in training are the same as the ones in the target distributions, the smoothing of the CG and target distributions will negate each other and the CG models would be trained to reproduce the same target structure, i.e., all the converged potentials will produce almost the same distribution functions if smoothed with the same bandwidths. Also, training with large bandwidths will have uniqueness issues; i.e., a set of strongly smoothed target distributions can be produced by more than one set of potentials due to the over–smoothed distribution functions of CG systems.

The trained CG models successfully reproduced the target distributions in all cases. A comparison between the target and CG distribution function for the $(w^{*}_{r},w^{*}_{l},w^{*}_{\theta})$ case is shown in Figure 4 as an example. In addition, the pair potentials of different cases are presented in Figure 5(a). The bond–angle potential of the case $(w^{*}_{r},w^{*}_{l},w^{*}_{\theta})$ is presented in Figure 5(b) as an example. It can observed that there are four major potential wells corresponding to the 4 peaks shown in Figure 4(c), representing different local conformations of PE chains. Expectation projection of bond–angle potentials in $l$- and $\theta$-direction, as defined in the caption, are presented in Figure 5(c) and (d), respectively. It is shown that as larger bandwidths are used, the BA potential would have less prominent peaks and crests, indicating that the energy barriers between the conformations are smaller. It is not recommended to use $w_{l}>4w_{l}^{*}$ or $w_{\theta}>4w_{\theta}^{*}$ if accurate representation of local conformation is needed, since some gauche conformations would not be produced due to disappeared corresponding energy wells.

Following the method of Li et al. ³⁰, the influence of the bandwidth parameters on the stable integration time step was determined by finding the largest time step for which the energy drift over a 1–ps NVE simulation was less than 1% of the initial kinetic energy of the system. Table 1 compares the performance of CG models of different bandwidth cases. It can be observed that the potentials trained in larger bandwidth cases can tolerate larger critical time steps and hence have better computational performance. It is verified that the highest–frequency mode is the bond vibration, since the critical time step significantly increases when the $l$-direction stiffness is reduced by using target distributions with larger $w_{l}$. From the perspective of time step–per–second, our CG–BA model has a performance of 29.78 time step/s with an 8000–monomer system and one MPI task, which is about 91% of the performance of the HB model. The slowdown is due to more computationally expensive bicubic interpolation compared to the simple linear interpolation used in the HB model and the cost induced by extra MPI gathering and communication of bond information.

To understand the effect of bandwidths on the crystalline phase chain diffusion behavior of semicrystalline polymers, we calculated the activation energy of the potentials trained with different targets. We constructed a crystal system in a periodic simulation cell with x, y, and z dimensions of 5.8, 5.3, and 20.5 nm, respectively. A total of 154 chains were placed parallel to the z-axis and are bonded across the periodic boundaries, representing an infinite lamellar thickness. The system were simulated in the anisotropic NPT ensemble for 50 ns at 1 atm, and various temperatures, namely 240, 270, 300, 330, and 370 K. Mean–squared displacements were measured to quantify diffusion. The diffusion coefficients were calculated with $D=\frac{g(t)}{2t}|_{t=50\mathrm{ns}}$ and the activation energy $E_{a}$ was calculated with the Arrhenius relationship. The results are shown in Figure 6.

Figure 6(a), (b), and (c) show the relationship between the activation energy and the $w_{r}$, $w_{l}$, and $w_{\theta}$, respectively. While increasing $w_{r}$ also tends to reduce the activation energy of crystalline diffusion, $w_{l}$ and $w_{\theta}$ demonstrate a stronger correlation with $E_{a}$, i.e., potentials from larger bandwidth trainings have a reduced activation energy of crystalline diffusion. We attribute the lower activation energy to the reduced energy barriers in $U(l,\theta)$ caused by stronger smoothing from larger bandwidths. It was shown that the defect movements are the agents of the chain–direction slip (or chain diffusion) in the crystal phases and account for the $\alpha_{c}$-relaxation.³⁹ This indicates that the state transitions of torsion angles in local chain segments, which form conformational defects and enable defect movements, would also facilitate diffusion. The torsion angle transitions in the atomistic representation are implicitly represented by the state transitions of bond–angle triplets in our CG mapping as shown in Figure 1, meaning that energy barrier between the local conformation states in $U(l,\theta)$ strongly affects the activation energy of crystalline diffusion.

We then identified the energy barrier $E_{b}$ of the BA potentials. More specifically, for a given BA potential, we first found all the local minima (states) in $U(l,\theta)$ and identified the trans state $T$. Then, we identified the minimum energy pathway from $T$ to all other states with the algorithm developed by Fu et al.⁴⁰ and calculated the maximum energy on the paths, which were regarded as the energy barriers of the paths. Eventually, the average energy barriers of all the paths are used as the energy barrier $E_{b}$ of the BA potential.

The relationship between the activation energy $E_{a}$ and the energy barrier $E_{b}$ is shown in Figure 6(d), where the data points follow a straight line. It verifies our hypothesis that the energy barrier in $U(l,\theta)$ correlates strongly with the activation energy $E_{a}$; and the reduction in the activation energy observed in crystalline diffusion $E_{a}$ is accounted by the decrease in $E_{b}$ due to stronger smoothing. Because the energy landscape is smoothened by the coarse–graining, the activation energy is lowered when compared with the experimental values of 5 to 22 kcal/mol³⁹, and 23 to 30 kcal/mol⁴¹, which will facilitate changes in local conformations. However, the controllability on $E_{b}$ provided by the bandwidth parameters offers a useful tool to adjust $E_{a}$ of CG models for different applications.

To understand how the bandwidth parameters can be used to accelerate crystallization kinetics to allow studies to be conducted at larger length and time scales, we simulated the crystallization process with 5 randomly generated systems using the trained CG potentials. For each system, 280 chains with 285 CG beads representing a molecular weight of 7980 g/mol were randomly placed in a cubic simulation box with side lengths of 10.6 nm. The chain length was selected so that polymer chains can crystallize into large crystallites and chains can pass through a crystallite multiple times and form loops. We then ran simulations in the anisotropic NPT ensemble at 300 K, 1 atm for 100 ns, during which the systems crystallized.

The crystallization of semicrystalline polymer consists of two processes. In primary crystallization, which is commonly modeled with the Avrami equation⁴², crystallites grow radially from nuclei and form spherulites. This process majorly consists of the growth of lamellae in the direction perpendicular to the stems (transverse direction). In secondary crystallization, lamellae thicken by reorganizing the polymer chains in adjacent amorphous phases in a manner similar to chain slips. We identified the characteristics of crystallites with the following procedures. After dividing the simulation boxes into bins of 7 Å, we calculated the $p_{2}$ parameter of the beads with a cutoff radius of 14 Å, and then identified the crystallites by finding the contiguous clusters of bins with average $p_{2}$ parameter greater than 0.4. To measure the secondary crystallization rate, we then calculated the average stem length in crystallites. To quantify the crystallite growth in the transverse direction, cross–sectional area of crystallites in the transverse plane were also estimated. The results are reported in Figure 7.

The crystallite growth follows a typical polymer crystallization process. The transverse growth starts and saturates quickly, providing the initial growth of crystallites, combined with stem growth that is responsible for the second–stage crystallite growth. Such a trend can also be observed from the rendered images in Figure 7(d). In image (i) and (ii), the system evolved from two nuclei to multiple small crystallites, meaning multiple nucleation sites formed which indicates sufficiently large domain size. The primary crystallization dominates the process from image (ii) to image (iii) which illustrates the system state only 2 ns later. Image (iv) demonstrates a fully crystallized state and after which no significant change was observed. The half–time of crystallization, defined as the time to reach 50% of relative crystallinity, is about 10 ns; when comparing to the experimental results by Zhuravlev et al.⁴⁴, we can roughly estimate that the acceleration is at most 100 times. More accurate estimation is impossible due to missing data for isothermal crytallization of polyethylene at 300 K. The CG model with accelerated crystallization kinetics is by no means ideal to reproduce realistic time scale, but does provide a viable and efficient approach to study the crystallization processes.

Multiple models are available to describe crystallite growth. The Avrami equation⁴² is commonly used to model crystallization kinetics, defined as:

where $\chi_{c}$ and $\chi_{c,\infty}$ denote the crystallinity and final crystallinity, $k$ is the crystallization rate, and $n$ is the Avrami exponent. The Avrami model assumes isotropic growth hence can not account for the secondary crystallization process. Among the works^{45, 46, 47, 48, 49, 50, 51, 52} that attempt to take the secondary crystallization into consideration, the Hay model^{53, 54} and Velisaris–Seferis model⁴³ were shown to generate the best fit according to Kelly and Jenkins.⁵⁵ We therefore selected the Velisaris–Seferis model because it converges to a finite value at infinite time. It is defined as:

where $k_{1}$ and $k_{2}$ signify the primary and secondary crystallization rates, respectively, $\chi_{c}$ denotes the crystallinity, $\chi_{c,\infty}$ is the final crystallinity, $w_{1}$ and $w_{2}$ are the weights of primary and secondary crystallization and sum to one, $n_{1}$ and $n_{2}$ are the Avrami exponents for primary and secondary crystallization. We selected $n_{1}=2$ and $n_{2}=1$ since the primary crystallization and the secondary crystallization are volume nucleation with two–dimensional and one–dimensional growth, respectively. We fitted the data in Figure 7(a) with the Velisaris–Seferis model to analyze the crystal growth. Naturally, the data shown in Figure 7(b) and (c) are fitted with the Avrami model. To distinguish the fitted parameters, the transverse growth rate is denoted by $k_{A}$ and the stem growth rate is denoted by $k_{l}$, which has the same physical meaning as $k_{1}$ and $k_{2}$, respectively. All the fittings matched with the data points well as shown in Figure 7 and the parameters are reported in Figure 8.

It can be observed that the trends in the row of $k_{l}$ and $k_{A}$ agree well with the row of $k_{2}$ and $k_{1}$ (e.g., the trend in Figure 8(a) is similar to the one in Figure 8(d)), respectively, which validates the fittings and post–processing since the parameters with the same physical implications are affected similarly. When considering the time to reach a certain relative crystsallinity level, the primary crystallization is generally 2 to 5 times faster than the secondary crystallization, consistent with the conclusion by Banks et al.⁵⁶ that secondary crystallization is less than one order of magnitude slower than the primary crystallization. The parameter $w_{r}$ reduces the secondary crystallization rate when $w_{r}$ is large, due to the deeper pair potential well in the ($8w_{r}$, $w_{l}$, $w_{\theta}$) case posing higher energy barrier to chain slips. $w_{l}$ has little effect on the crystallization kinetics in the investigated range. Increasing $w_{\theta}$ greatly increases $k_{A}$ and $k_{1}$; the increase of $k_{A}$ and $k_{1}$ is attributed to a weaker constraint at $\theta=\pi$ resulting from larger $w_{\theta}$, which facilitates the folding of the chains and accelerates the crystallite growth in the transverse direction.

In this work, we developed an extension to the IBI method that considers correlated multi–variable intramolecular interactions. This method also incorporates the KDE representation in distributions and energy formulations, which provides controllability to the smoothness of the resulting potential functions and the characteristics of the resulting models.

We trained 10 CG models of PE with different kernel bandwidths using this method. The models successfully reproduced the correlated distributions of bond–length and bond–angle of PE with an $L_{2}$ norm error of about 1%. The conformation states of bond–angle triplets, which represents the underlying atomistic conformations, were reproduced with appropriate bandwidths. We found that the activation energy of crystalline chain diffusion of the CG model is determined by the energy barrier $E_{\mathrm{b}}$ between the conformational states. By adjusting the bandwidths of the target distributions, we could tune the magnitude of the peaks and crests in the $l-\theta$ energy landscape, giving us controllability towards the system dynamics. Increasing the bandwidths will accelerate the dynamics by smoothing out the peaks and crests at the cost of weakened ability to accurately reproduce the local conformations; it will also reduce the stiffness of the BA potential and the maximum frequency of the model, leading to increased critical time step and computational performance. By adjusting the bandwidth, one could focus on accurately reproducing the local conformations, or efficiently cover long time scale in simulations

In the study of isothermal crystallization of supercooled PE, the models successfully crystallized the systems and reproduced the primary and secondary crystallization kinetics. The model trained with larger $w_{r}$ would slow down the secondary crystallization rates because deep pair potential well creates resistance for chain slip, whereas $w_{l}$ have little effect on the crystallization kinetics within the investigated range. Additionally, increasing $w_{\theta}$ greatly increases the primary crystallization rate, indicating the flexibility of chains has a significant impact on the crystallite growth in the transverse direction. The CG model presented in this work shows controllable characteristics in diffusion dynamics and crystallization kinetics. We look forward to further testing the capability of the model in simulating mechanical properties.

The authors acknowledge support from the National Science Foundation (NSF) under award Division of Civil, Mechanical, and Manufacturing Innovation (CMMI) 1653830. The authors acknowledge Research Computing at Arizona State University for providing HPC resources that have contributed to the research results reported within this paper.