Introduction
In today’s society, industries have extremely high demands that force them to mass produce to meet these needs. In the manufacturing process, industries mostly opt for the use of plastic due to its wide variety of applications, durability and cost. However, the depletion of oil resources along with increasing environmental pollution have resulted in a global concern that has led to the search for another more sustainable approach for the environment1.
Poly (ε-caprolactone) (PCL) is a widely recognized and extensively researched biodegradable polymer, highly valued for its ability to degrade under physiological conditions2. Derived from petroleum, this material is semicrystalline and hydrophobic, with key properties such as crystallinity and biocompatibility strongly dependent on its molecular weight3. PCL offers a unique combination of characteristics, including high solubility, a low melting point and glass transition temperature, and high compatibility in mixtures, making it a material with significant potential for applications such as tissue engineering2 and additive manufacturing4,5.
A fundamental understanding of PCL’s properties at the molecular level is crucial for optimizing its performance in these advanced applications. While atomistic molecular dynamics (MD) simulations are powerful tools for investigating the dynamic evolution of molecular systems and condensed matter, they are often limited by accessible length and time scales. Simulating complex polymer-based materials, especially high molecular-weight entangled chains over industrially relevant scales, presents a significant challenge due to the computational cost and time constraints.
To overcome these limitations, multiscale modeling approaches are increasingly employed6–8. Specifically, systematic coarse-grained (CG) models are developed to reduce the resolution of the system by lumping groups of chemically bonded atoms into single beads. This allows for the simulation of larger systems over extended time periods, significantly improving computational efficiency while still capturing essential material properties7. The development of such CG models often follows a “bottom-up” approach, where information from more detailed, fine-grained (atomistic) simulations is systematically used to derive the parameters for the reduced-resolution models. This methodology typically involves two key steps for building a CG model: i) defining the CG mapping: This establishes the correspondence between the atomistic and the reduced resolutions, determining how atoms are grouped into CG beads; ii) defining the CG effective potentials: This involves deriving the interaction potentials for the CG beads based on the statistics obtained from the atomistic system. Techniques such as iterative Boltzmann inversion (IBI)9 are commonly utilized to match the local structural distributions of the CG model to those of the atomistic model. This ensures that the CG model accurately reproduces the structural and, with appropriate time scaling, the dynamical properties of the underlying atomistic system.
This systematic linking of simulation methodologies across different scales is essential for quantitative prediction of polymer behavior over broad spatiotemporal scales. For PCL, all-atom molecular dynamics simulations, often employing force fields like the Optimized Potentials for Liquid Simulations All-Atom (OPLS-AA) force field10–12, have been used to study fundamental properties such as melt density, transition temperatures and mechanical characteristics13–15. Yungerman et al., utilized the OPLS-AA and Class II Polymer Consistent Force Field (PCFF) to investigate the properties of telechelic PCL with diacrylates as reactive functionalities15. Their results indicated that both PCFF and OPLS-AA were capable of describing satisfactorily the geometry of individual PCL chains in dilute and melt systems, and both force fields provided properties such as density and Young’s modulus comparable to the experimental data. The articles by Di Pasquale et al. and Lavino et al. primarily employed the OPLS-AA force field in their atomistic investigation of PCL chains consisting of up to 30 repeating units within various solvent environments, specifically water-acetone mixtures16–18. Ezquerro et al. chose OPLS-AA force field over AMBER force field to model polymer melt due to its superior performance in predicting PCL’s density14. The research specifically investigated neat PCL polymer matrices of relatively short chains containing 20 monomers, which were then blended into nanocomposite systems containing either ungrafted silica nanoparticles or silica nanoparticles grafted with polyethylene glycol. Sohrabian et al. also investigated a nanocoposite system, but employed COMPASS (Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies) force-field to model the PCL matrix19. Recently, Nie et al. introduced a priori temperature rescaling method (Hybrid Iterative Boltzmann Inversion - HIBI) for coarse-graining linear PCL chains in melt with the molecular weight of 5700 g/mol6. The method aims to compensate for the loss of entropy by scaling up the temperature a priori and successfully reproduces structural distributions, heat capacity, and Young’s modulus, however, its prediction for viscosity falls short compared to the all-atom systems6.
The current work has as its primary objective to develop a reliable and relatively simple systematic coarse-grained model for linear PCL with the capability of being extended to experimentally relevant time and length scales. We based the model on all-atom simulations of PCL which employed OPLS-AA force field adapted for studying long chains. The atomistic data were first compared to the literature data and then the polymer chains were coarse-grained at the monomeric level. The IBI method was used to capture the structural features of the studied melts and the obtained potentials were employed to simulate CG systems.
A fundamental novelty of this work resided in the investigation of PCL chains across a wide range of molecular weights, specifically from 10 up to 125 monomers (equivalent to 1143 to 14270 g/mol). This expanded range, covering the molecular weight regime from unentangled (Rouse-like) to mildly-entangled systems, provides a more robust platform for studying the influence of chain length effects on fundamental structural and dynamical properties of PCL. Study of the chain length effects is particularly relevant for industrial applications, such as additive manufacturing (3D printing), where properties like rheological behavior (e.g., viscosity) are affected by the degree of crystallinity and thus by the molecular weight of the polymer melt. The ultimate goal is to facilitate the development of the computational approaches, which would speed up the process of tuning the properties of PCL-based materials and by doing so reducing the ecological and economical impact of the property optimization of these new generation of materials.
Methods and simulation details
In our study we employed all-atom molecular dynamics simulations to study linear chains of poly (ε-caprolactone) (PCL, see the structure in Figure 1(a)) of multiple molecular weights. The simulated systems are listed in Table 1. In order to describe the interatomic forces in the systems, two types of force fields were combined, the original OPLS-AA10 and the modified L-OPLS11. More specifically, L-OPLS was used as the main source of information, since this force field is based on the parameters of OPLS-AA but was optimized for long hydrocarbon chains, as well as for alcohols and esters11. This optimization is in line with the objective of this work, namely, to provide a reliable and relatively simple model of PCL that could be extended to experimentally relevant time and length scales as well as to relevant molecular weights. In addition, it provides an improvement with respect to the previously published data, which relied on OPLS-AA for description of oligomeric PCL. The full set of parameters used here together with the corresponding source can be found elsewhere20. All the simulations were performed by open-access simulation package Gromacs21.
Figure 1.
(a) Chemical structure of PCL. (b) Snapshot of the PCL chain of 30 monomers in melt. One chain is highlighted while the surrounding molecules are transparent for a better visibility. (c) Schematic illustration of the coarse-graining procedure for PCL10. Each orange sphere represents one CG unit placed in the center of mass of a monomer.
Table 1. Details about the systems under study.
Mw refers to the molecular weight of the chain, Nc denotes the number of chains in the system, N is the number of monomers per chain, ρA refers to the density of the atomistic system and ρCG to the density of the coarse-grained system.
Label of
the system | N | Mw [g/mol] | Nc | ρA [kg/m3] | ρCG [kg/m3] |
|---|
| PCL10 | 10 | 1143 | 70 | 928 ± 7 | 930±30 |
| PCL30 | 30 | 3426 | 70 | 928 ± 4 | 940±20 |
| PCL50 | 50 | 5709 | 70 | 928 ± 3 | 950±10 |
| PCL100 | 100 | 11416 | 70 | 928 ± 2 | 939±9 |
| PCL125 | 125 | 14270 | 70 | 928 ± 2 | 939±9 |
Due to the lack of data on PCL modelled by L-OPLS force field, we designed the methodology as follows. In the first stage, we simulated the PCL chains in water. This step in our methodology serves as a comparative analysis with prior research. A key consideration was the lack of available studies that investigate a range of PCL’s molecular weights in a molten state. Consequently, we opted to use solvent simulations to verify our model. The radius of gyration Rg of these systems was then confronted with the published data, to corroborate the behaviour of our systems. The detailed description of the preparation of these systems, the full set of data and the results used for the comparison of Rg with the literature data are publicly available20. Briefly, after the preparation of a single chain molecule with the desired molecular weight (see Table 1), the chain was placed in vacuum to achieve a collapsed configuration. In continuation, so-obtained molecule was simulated in water at 300 K. The results obtained from these simulations were highly satisfactory.
The initial state of systems in melt was composed of the partially collapsed chains (see the description in ref. 20 for more information about the preparation of these chains). 70 polymer chains were randomly placed in a simulation box, big enough to avoid overlaps. The density was then adjusted by performing a short simulation at 50 bar. This preparation process is similar to the process published in ref. 22 and is used to remove gaps between chains and obtain a homogeneous density. The systems were then equilibrated in a similar way as published in ref. 23. Briefly, the systems were first equilibrated at 600 K, then they were cooled down to 500 K and then equilibrated at this temperature. The length of the equilibration run at each temperature depended on the chain length and the radius of gyration was monitored as an indicator of the equilibrated structure. Before the production run, a run of 100 ns was performed applying the same conditions as the production run. Namely, the PME method was used for the calculation of the electrostatic interactions, the velocity rescaling with a stochastic term was applied to maintain the temperature at 500 K and the Parrinello-Rahman barostat kept the pressure at 1 bar. The time step was 2 fs and the bonds containing hydrogens were constrained by the LINCS algorithm24. The atomistic production run was 400 ns long for chain lengths of N = 10 − 100 and 600 ns long for PCL125. All runs, the equilibration runs and the productions runs, are publicly available25. The coarse-grained simulations were performed at constant temperature of 500 K and constant pressure of 1 bar, using Nosé-Hoover thermostat and the Parrinello-Rahman barostat. The time step was 5 fs and the bonds were constrained by the LINCS algorithm24. The coarse-grained simulations are also publicly available26.
Results and discussion
Polymer chains in melt typically exhibit dense and homogeneous packing, where monomers are isotropically surrounded (see also Figure 1(b)) and long-range intermolecular forces are effectively screened. This screening behavior in the melt state, where chains “screen” each other from long-range forces, causes the chain to statistically behave like an ideal chain without the excluded volume effect. This leads to a particular theoretical scaling of fundamental structural properties with the number of monomers in the chain and since in the validation step we focused on confronting the results from the atomistic model with systems in solution, in what follows we focus on confirming this theoretical assumption in melt for both atomistic and coarse-grained model derived systematically from the atomistic configurations.
In the first place, to characterize the polymer’s local structure, monomer distribution functions were calculated directly from the atomistic simulation data. These distributions, encompassing bond lengths, angles, dihedrals, and monomeric radial distribution functions, served as critical target data for the subsequent coarse-graining procedure. In this procedure, the atoms were grouped into “super-atoms” or beads. More specifically, in our case each bead represented one monomer, i.e., the internal beads of the chain consisted of 18 atoms and the terminal ones of 19 and 20, due to the presence of the carboxyl and hydroxyl groups at the extremes (see also Figure 1(a)). The centre of the bead was then placed in the centre of mass of each monomer. One example of such a coarse-graining procedure is shown in Figure 1(c).
The monomer distribution functions measured from the atomistic simulations are shown in Figure 2 together with the sketch of the measured internal geometrical parameters. Considering the calculation of the non-bonded radial distribution function, note that the short-range intramolecular contributions were excluded, since those short-range interactions are intended to be described by the bonded potentials. More specifically, the exclusion was defined by the number of bonds separating two beads (i.e., CG units) in a chain. In our case, only the intramolecular distances among beads placed more than 3 bonds away from each other are included in the radial distribution function.
Figure 2.
(a) Schematic illustration of the bonded interactions and the parameters plotted in the monomeric distributions (b–e) calculated from the atomistic simulations. Each bead represents one single monomer. (b) Distribution of the bond length between monomers. (c) Distribution of the angles between monomers. (d) Distribution of the dihedrals between monomers. (e) Monomeric radial distribution function for the studied systems.
Despite covering a wide range of molecular weights, no significant differences in these distributions across studied PCL chains were found. This observed consistency in local structural features, regardless of the overall chain length, implies that the immediate chemical environment and inherent conformational flexibility of the PCL monomer are largely independent of the macroscopic chain size in the melt. This independence significantly simplifies the development of coarse-grained models, as a potential derived from shorter oligomers can reliably be applied to longer chains for these specific local properties, forming a robust basis for bottom-up polymer modeling.
In continuation, we applied Iterative Boltzmann Inversion (IBI), as one of the most powerful methods used to derive effective mesoscale potentials from atomistic simulations9. We followed closely the method published elsewhere8, therefore we will only summarize briefly the main concepts. The coarse-grained (CG) potentials were generated in a two-stage process, beginning with an Initial Direct Boltzmann Inversion to convert the target atomistic distributions into initial estimates of the potentials. These potentials were then iteratively refined using the IBI method. The refinement was typically structured sequentially: bonded, angular, dihedral, and non-bonded potentials were tuned one after the other in cycles, in this given order. Before the pressure correction, all the contributions were also tuned at the same time. The convergence of the Iterative Boltzmann Inversion is quantified by defining a norm as the integral of the absolute difference between the target and the current step's unnormalized distributions. To ensure consistent scaling across all components, each distribution is first normalized such that the target distribution's integral equals unity. The overall convergence norm is then calculated as the mean of the norms from all distinct contributing potentials. IBI iterations stops when this mean norm drops below a set tolerance, or when the improvement between successive steps falls below a specified small percentage. In our IBI procedure, the maximum number of iterations for each potential contribution (bonded, angular etc.) was 20. Each CG simulation within an iteration was run for up to 1500000 steps.
Figure 3 illustrates the results of the IBI procedure in comparison to the target atomistic data, which were chosen to be the data from the longest simulated polymers, PCL125. The notable agreement observed between the IBI curves and the target data across all distributions confirms that such a simple and versatile method is capable of accurately reproducing the structural characteristics derived from atomistic simulations throughout the entire studied molecular weight range.
Figure 3. Results of the IBI procedure.
Each bead represents one single monomer. (a) Distribution of the bond length between monomers. (b) Distribution of the angles between monomers. (c) Distribution of the dihedrals between monomers. (d) Monomeric radial distribution function for the studied systems. The dashed line represents the data from the atomistic simulation of PCL125.
Once ensuring that the sampled distributions achieved good agreement with the target data obtained from atomistic simulations, the derived potentials were subsequently utilized to simulate the coarse-grained (CG) systems. It is noteworthy that the reduction in the system's degrees of freedom achieved through the CG methodology directly results in a significant enhancement of computational efficiency. Specifically, a comparative performance analysis of the atomistic and CG simulations for the longest chain, PCL125, demonstrated a speed-up factor exceeding 100. This performance metric was established using a computational cluster configured with four Intel(R) Xeon CPUs E5606 @ 2.13 GHz.
In what follows we present structural analysis of both atomistic and coarse-grained systems to demonstrate that our simple methodology renders consistent results at atomistic as well as at the mesoscale.
One important measure of the internal molecular structure is the plot of internal distances, Rn2/n, as a function of the separation n along the chain. When measured among pairs of atoms along the atomistic backbone, this quantity is directly related to the characteristic ratio Cn (or C∞), which serves as a crucial measure of a polymer chain’s intrinsic flexibility or extension in its unperturbed state. The Cn calculated directly from the atomistic simulations is shown in Figure 4. The characteristic ratio Cn is defined as where is the mean squared distance between atoms, n represents the separation (of the atoms) along the chain backbone, and bA is the average bond length in the atomistic simulations. In Figure 4 the maximum separation equals to nmax = 7N − 1, where N is the number of monomers in the chain. The value of Cn changes with n. As n gets very large, Cn approaches a maximum value, which we call C∞. The value of C∞ visually estimated from the plateau at long n in Figure 4 is C∞ ≈ 5. This estimated value is in good agreement with the previously published experimental results (C∞ = 3.9−627–32) and with the theoretical estimation from the rotational isomeric state calculations (C∞ = 4.1 − 4.531).
Figure 4. Characteristic ratio Cn as a function of the distance between atoms along the backbone of the atomistic chains.
In Figure 5(a) the internal distances, Rn2/n are shown as a function of the separation n along the chain for studied PCL systems. Note that in this case, in order to be able to compare the atomistic and CG systems, n denotes the separation of the monomers (i.e., the difference in their indices) and not of the atoms. Giving the connection between the internal distances and the characteristic ratio, the differences in Rn2/n in Figure 5(a) reflect the differences in chain stiffness. As depicted in the figure, the value of Rn2/n generally increases with increasing separation n. For longer PCL chains, particularly those with a higher number of monomers (e.g., PCL100 and PCL125), the function demonstrates a tendency to reach a plateau at larger n values. For n equal to the number of monomers N, Rn2 represents the squared value of the average chain span, i.e., the end-to-end distance Re2. The visual agreement between the atomistic and coarse-grained results in Figure 5(a) is particularly good for PCL10, and PCL30, which highlights the consistent representation of these internal distances across different simulation resolutions. For longer chains, the CG model seems to overestimate the chain stiffness, leading to higher values of the plateau at high values of n. This deviation of larger-scale dimensions is attributed to the fact, that while IBI performs very well at reproducing local structural correlations, it often struggles to accurately capture larger-scale, long-range properties of polymer chains. The phenomenon is not inherent to IBI and was reported in various studies dealing with the coarse graining8,32.
Figure 5.
(a) Monomer internal distances as a function of the separation n between monomers, for the atomistic (shaded area) and CG systems (lines). The largest error bars for the CG systems are also plotted. (b) Radius of gyration as a function of the number of monomers, N. The inset shows the ratio of the end-to-end distance, Re, to the radius of gyration, Rg. The dashed line in both cases represent the theoretical prediction, namely Rg ∼ N0.5 and Re/Rg = √6. The error bars were calculated by the standard block method, using 4 blocks of simulation data.
Figure 5(b) presents the radius of gyration Rg as a function of the number of monomers N for both the atomistic and coarse-grained representations of PCL. Our analysis of the atomistic data for Rg reveals a clear upward trend with increasing number of monomers, consistent with theoretical scaling predictions for polymers in the melt, namely Rg ∼ N0.5. This finding is important as it confirms that our atomistic models accurately capture the unperturbed conformational state of PCL in a dense environment.
When examining the coarse-grained results for Rg, we observe that these values for chains with N > 30 are generally above those obtained from the atomistic simulations for the same chain lengths. This outcome is in line with slightly more stretched character of these CG chains in comparison to their atomistic analogues observed in Figure 5(a). Despite the quantitative differences, the CG Rg values seem to maintain a similar scaling trend with N as their atomistic counterparts, indicating that the coarse-graining procedure successfully preserves the fundamental chain-length dependence of structural dimensions.
The inset of Figure 5 provides further insight into the conformational characteristics by illustrating the ratio of the end-to-end distance and the radius of gyration, Re/Rg. For ideal chains, this ratio is theoretically Re/Rg = √6. Our atomistic PCL chains in melt show Re/Rg ratio that approaches this theoretical ideal value as chain length increases, more specifically, the theoretical value is reached for N > 30. The CG models also reflect this trend, although minor deviations from the atomistic data are present.
Figure 6(a) presents the asphericity distribution functions for the studied chains. The asphericity quantifies the shape of the molecules and is measured in simulation by calculating the geometrical inertia tensor:
Figure 6.
(a) Asphericity distribution for atomistic (points) and coarse-grained (lines) PCL systems. The dashed black line represents the theoretical expression reported in ref. 36 for unentangled poly(ethylene oxide). (b) Form factor divided by the number of particles in the molecule as a function of the q vector for atomistic (points) and coarse-grained (lines) PCL systems. The solid line indicates the ideal scaling in the intermediate range of q values.
where α,β denote the Cartesian coordinates and Np represents the total number of particles in the chains, i.e., the number of atoms in the atomistic systems and Np = N for the CG systems. The eigenvalues of Tαβ, labelled as λ1, λ2, λ3, are subsequently used to estimate the asphericity parameter from the following relation:
The more spherical the molecule, the closer its asphericity value is to 0. The distributions in Figure 6(a) for both atomistic and CG systems are identical within the obtained accuracy, with a sole exception of PCL125 system, which appears slightly more spherical at atomistic scale. The collected data have some features of the bimodal distributions expected by the theoretical model of a chain in θ conditions33. The bimodal character originates from the two distinct molecular conformations present in θ state: a collapsed conformation with lower asphericity and a coil conformation. The distributions are in a perfect agreement with the data reported previously for other common unentangled synthetic polymers in melt (see the dashed line in Figure 6(a) and the Supplementary Information in ref. 34 which reports the unpublished data for refs. 8,35 and data from ref. 36) and for poly(lactic acid) in melt34, which confirms that the used model reproduces the expected behaviour of polymers in dense environments. It is important to note that, in contrast to poly(lactic acid) (PLA)34, a related sustainable polymer, we found no evidence that the asphericity parameter depends on chain length, even with the shortest PCL10 chain.
To gain further insight into the internal molecular packing, the single-molecule form factor, W(q), was calculated. For a robust comparison between the coarse-grained and atomistic chains, W(q) was subsequently normalized by the number of particles Np per chain. The form factor probes the internal molecular packing at different characteristic length scales, represented by 1/q, where q is a wave vector. Focusing on the length scales larger than the bond length (i.e., ⟨bA⟩ ≈ 0.148 nm for atomistic and ⟨bCG⟩ ≈ 0.66 nm for CG systems) and lower than the radius of gyration, 1/b ≥ q ≥ 1/Rg, data in Figure 6(b) show a perfect overlap between the CG and atomistic systems for the given chain length. In addition, within this q regime, a fractal regime W(q) ∼ q−1/ν is anticipated and confirmed. The determined exponent ν = 0.5 is typical for ideal Gaussian chains. In line with previous findings, this observation reinforces that the chains are following their predicted behavior. It is interesting to note, that a sudden drop of W(q) at q ≈ 3 nm−1 and a consequent maximum at q ≈ 11 nm−1 is a consequence of a spherical nature of the CG units (see e.g., ref. 37 for W(q) of a compact sphere).
As also described and demonstrated above, IBI method excels at reproducing the structural distributions between coarse-grained super-atoms, matching those obtained from atomistic simulations. However, CG models often exhibit faster dynamics compared to their atomistic counterparts, due to the implicit nature of friction and the averaged degrees of freedom. To accurately reproduce the dynamical properties, it is necessary to rescale the dynamics of the CG systems, often by applying a time scaling factor that compensates for the lower monomeric friction coefficient in the coarse-grained model38–40. In order to obtain this factor, we calculate two dynamical properties: the characteristic time related to the rotational movement of the chain and mean squared displacement of the center of mass.
The orientational autocorrelation function of the end-to-end vector for a chain provides a hint about the terminal dynamics of the chain, i.e., about the chain relaxation at long times. This function, denoted as C(t), quantifies the temporal correlation of the chain’s end-to-end vector orientation and is defined as:
where is the end-to-end vector of the polymer chain at time t+t0, and is the end-to-end vector at the time origin t0. The angle brackets ⟨·⟩ denote an ensemble average.
The data on C(t) collected from the atomistic and CG simulations are show in Figure 7(a). Note that the autocorrelation functions for the CG chains are considerable slower than those of their atomistic analogues and that the time axis in the CG data was rescaled (see below).
Figure 7.
(a) Correlation function C(t) of the end-to-end vector as a function of time for all simulated systems. (b) Mean squared displacement of the center of mass of the chains (in nm2) as a function of time. Lines correspond to the CG systems and points to the atomistic ones. The time t for the CG data in (a,b) was multiplied by the shift factor tCG to match the atomistic data. The thick solid and dashed black lines in (b) are guides for the eye and illustrate the power laws indicated above them.
Analysis of the decay of C(t) allows for the extraction of relaxation times that govern the overall reorientation of the chain. While theoretical models like the Rouse model describe the dynamics of unentangled, short polymers, predicting an exponential decay for its Rouse modes, real polymer melts often exhibit more complex features, such as stretched exponentials, due to factors like intramolecular bond correlations, excluded volume effects, and chain stiffness. Therefore, we fitted the data with the following function:
where A ≤ 1 is a prefactor, τ is the characteristic time (either τA or τCG) and β is the stretch exponent. The results obtained from this analysis for the atomistic samples are summarized in Table 2. The same procedure was repeated for the CG samples, obtaining τCG and βCG. Comparing the characteristic times in Table 2, the faster dynamics in CG systems is corroborated. Comparing the stretch exponents and the data overlap in Figure 7, it is clear that the C(t) functions for CG chains exhibit behavior more akin to simple exponential decay, as it can be anticipated for models with reduced degree of details. Their systematically higher βCG values, relative to the βA values of comparable molecular weight chains, lead to noticeable deviations in the C(t) functions at early stages.
Table 2. Results of the analysis of the dynamical properties.
τ and β refer to the characteristic time and stretch exponent obtained from the fit of the end-to-end correlation function. Parameters labelled with A were obtained from the atomistic simulations and those labelled with CG from the CG simulations. The error bars were calculated by the standard block method using 4 blocks of simulation data. tCG is the shift factor obtained as the ratio tCG = τA/τCG. Units labelled with * are units in CG simulations, which do not correspond to the physical units.
Label of
the system | τA [ns] | βA | τCG [ns]* | βCG | tCG |
|---|
| PCL10 | 1.42±0.04 | 0.73±0.04 | 0.12±0.01 | 0.89±0.06 | 12±1 |
| PCL30 | 13.5±0.7 | 0.65±0.04 | 1.1±0.3 | 0.67±0.05 | 12±3 |
| PCL50 | 42±5 | 0.62±0.05 | 2.7±0.4 | 0.6±0.1 | 16±3 |
| PCL100 | 350±110 | 0.48±0.06 | 12±3 | 0.55±0.08 | 29±14 |
| PCL125 | 640±80 | 0.43±0.02 | 20±3.5 | 0.6±0.1 | 32±7 |
The scaling factor tCG, essential to map the CG data and to account for differences in the friction in both types of simulations, was calculated as tCG = τA/τCG. Note that after applying this scaling factor and rescaling the time axis in the C(t) time dependence, there is an apparent visual mismatch in CG and atomistic data for PCL100 and PCL125. We attribute this to the fact that by rescaling the data by tCG, i.e., by a ratio of characteristic times obtained from the stretched exponential function, the overlap should happen primarily around the time scales where C(t) ≈ 1/e, thus around the characteristic time τA. Unfortunately, data from the atomistic simulations in this given time window are noisy and the rescaled CG data are expected to show a good match with the atomistic systems at the time scales not reached in our simulations.
To describe the translational dynamics, we calculated the mean squared displacement of the center of mass, g1(t):
where RCOM,i(t0) is the position of the center of mass of the ith chain at time t0. The data for the atomistic systems and the rescaled data for the CG systems are shown in Figure 7(b). At intermediate time scales, before reaching the diffusive regime, g1 exhibits a sub-diffusive behavior very similar to the one reported for multiple unentangled systems, g1 ∼ t0.8541. Detailed analysis of the origin of this sub-diffusive motion is beyond the scope of this work, however, as described in ref. 41, it indicates that unentangled polymers do not behave as ideal chains in a uniform medium. Instead, their dynamics are intrinsically linked to interactions with their polymeric environment, acting as a “fixed-point” effect on their dynamic properties41. Concerning the agreement between the atomistic and rescaled coarse-grained g1(t)s, as the chain length increases the time window of the perfect overlap shrinks, and particularly, for PCL100 and PCL125, the data only match at the long time scales. This is in line with our previous observation about the expected agreement at terminal times scales, i.e., in the diffusive regimes, which are in the case of the longest chains beyond the simulation time window.
In order to provide an alternative estimation of tCG, we calculated the scaling factor from the g1 data. We chose a characteristic time for the center-of-mass diffusion, t1 as the time at which the chain diffuses its own size, i.e., at the time when g1(t1) = Rg2. Consequently, we estimated tCG as the ratio of these characteristic times in the atomistic and in the CG simulations for the given chain length, i.e., tMSDCG = t1,A/t1,CG. Since only the chains with N < 100 reached the diffusive regime, the estimated values of tMSDCG are 15, 15 and 20 for PCL10, PCL30 and PCL50. Assuming that the value of the scaling factor should reach a constant value at high molecular weights8, we used tMSDCG = 20 for PCL100 and PCL125 to be able to plot the full set of data in Figure 8. Figure 8 shows that while this alternative estimation leads to a better apparent agreement in g1(t) for PCL10, PCL30 and PCL50 and at intermediate time scales also for PCL100 and PCL125, it overestimates the friction in rotational dynamics in CG systems with low molecular weight and N < 100 (see the C(t) functions in Figure 8(a)).
Figure 8.
(a) Correlation function C(t) of the end-to-end vector as a function of time for all simulated systems. (b) Mean squared displacement of the center of mass of the chains (in nm2) as a function of time. Lines correspond to the CG systems and points to the atomistic ones. The time t for the CG data in (a,b) was multiplied by the shift factor tMSDCG estimated from the onset of to the diffusive regime in the mean squared displacement. The thick solid and dashed black lines in (b) are guides for the eye and illustrate the power laws indicated above them.
Since our primary goal was to ensure the transferability of the CG potentials to much longer polymer chains, we rely on the fact that the dynamic scaling factor is expected to converge to a chain-length independent stable value once the chains are long enough; this asymptotic value can then be applied in future CG models for longer systems. Summarizing the observed behavior, the anticipated convergence is clearly supported by the scaling factor derived from the correlation function analysis (tCG), as the values obtained for PCL100 and PCL125 are statistically indistinguishable within their estimated error bars. Consequently, the scaling factors determined for these longest simulated chains serve as the recommended, converged value for future models. In contrast, the scaling factor estimated from the MSD, t MSD CG, was only reliably calculated for shorter chains (N < 100), which prevents us from confidently extrapolating its stable value to greater chain lengths.
To provide a quantitative analysis of the dynamical properties, we calculated the terminal relaxation time using the data reported in Table 2 from the following relation:
where Γ denotes the gamma function. The data are plotted as a function of the number of monomers in Figure 9. Note that the terminal relaxation time obtained in this way can be related to the terminal relaxation in rheological measurements. In addition, the Rouse model predicts the scaling of τend ∼ N2 for the unentangled linear polymers42 while as the entangled regime is reached, a scaling of τend ∼ N3.4 is anticipated43. Looking at the illustrated power laws in Figure 9, it seems that the crossover between the Rouse regime and the entanglement regime occurs at chain lengths around 50 monomers. The critical molecular weight Mc, which marks the onset to the entanglement regime, is usually a multiple of the entanglement molecular weight, Me, with Mc/Me being between 2 and 744. Having in mind the limited number of points in our study, the estimation of Me, or the number of monomers per entanglement strand Ne is not directly possible from the plot in Figure 9.
Figure 9. Terminal relaxation time τend as a function of the number of monomers N.
The black lines are guides for the eye and indicate the scaling regimes for the Rouse model (solid line) and for the entangled regime (dashed line). The grey arrow points out to the value of 2Ne estimated from the packing length.
As the estimation of crossover in Figure 9 is only visual, we proceed with a more rigorous way related to the structural property of a system called packing length, p. Packing length relates the size of the polymer, namely its end-to-end distance Re, to the density and the monomer mass Mw:
where NA denotes the Avogadro’s constant. Considering that in our systems the density is independent of the molecular weight (see Table 1) and the scaling of Re with N and thus with Mw is ideal, the packing length does not exhibit chain length effects. Therefore, we averaged the obtained p over the studied chain lengths, which resulted in ⟨p⟩ = 0.27 ± 0.01 nm. This value is close to the value of various synthetic polymers, such as polybutadiene (p = 0.21 nm) and polyisoprene (p = 0.31 nm)44. In addition, is in between the values reported for experimental PCL systems, namely p = 0.166 nm in ref. 27 and p = 0.315 nm in ref. 45. It is interesting to note that the calculation of the packing length from experiments is usually done through the relation GN ∝ p−3, where GN is the entanglement modulus27. Therefore, even a small difference among published GN from the rheological experiments may cause a huge difference in the final value of p. Having estimated the value of p from simulations, we can calculate Me such as:
where nt is a dimensionless number with a value of approximately 2144. The resulted value equals to Me = 4850 g/mol, which results in Ne ≈ 43. Assuming for simplicity that Mc = 2Me, the onset to the entanglement regime would lie around N ≈ 86 (see the arrow in Figure 9). The experimental values for Me of PCL range between 1600 g/mol to 3900 g/mol27,46–48, with the latest reported value of 8626 g/mol45. During the quantitative comparison, one should keep in mind that the experimental values were estimated from the rheological measurements and thus effects such as polydispersity, crystallization and temperature effects should be considered. Giving that our estimation of Me relies solely on static properties, the agreement between the reported values and the one obtained in this work is satisfactory. This agreement is corroborated by the visual transition observed in Figure 9.
Conclusions
This study successfully developed and validated a systematic coarse-graining approach for poly(ε-caprolactone) (PCL) in the melt state, demonstrating its capability to quantitatively capture key structural and dynamical properties across a wide range of molecular weights.
Our atomistic simulations, combining the L-OPLS and OPLS-AA force field for a better affinity with the long chains, provided robust and consistent local structural features (bond lengths, angles, dihedrals, and radial distribution functions) that were largely independent of chain length. This consistency significantly simplified the coarse-graining process, allowing for the derivation of effective mesoscale potentials from shorter oligomers that could be reliably applied to longer chains. The Iterative Boltzmann Inversion (IBI) method proved highly effective in reproducing the structural characteristics from atomistic simulations across the entire studied molecular weight range, confirming its versatility and accuracy for local correlations.
Analysis of structural properties at both atomistic and coarse-grained levels revealed good agreement for internal distances, particularly for shorter chains (PCL10, PCL30). For longer chains, the coarse-grained model showed a tendency to overestimate chain stiffness and larger-scale dimensions, a known limitation of IBI in capturing long-range properties. Despite these deviations in chain stretching, the coarse-grained model successfully maintained the fundamental chain-length dependence of structural dimensions, as evidenced by the consistent scaling trend of the radius of gyration, Rg, with the number of monomers, N. Both atomistic and coarse-grained Rg values exhibited scaling consistent with theoretical predictions for ideal Gaussian chains in the melt (Rg ∼ N0.5), indicating that the chains adopt an unperturbed configuration. Furthermore, the ratio of end-to-end distance to the radius of gyration (Re/Rg) for atomistic PCL chains approached the theoretical ideal value of √6 for chains with N > 30, a trend also reflected in the coarse-grained models. The analysis of conformational shapes, including asphericity distributions and form factors, also demonstrated good agreement between the atomistic and coarse-grained representations, further validating the structural accuracy of the coarse-grained model.
Regarding the dynamical properties, we presented results for orientational and translational dynamics of the entire chains. The coarse-grained model was appropriately time-scaled to account for the difference in the chain friction. Both dynamical properties, relaxation of the chain end-to-end vector and the mean squared displacement of the center of mass were used to extract the scaling factors. Since the dynamic scaling factor converges to a chain-length independent stable value for long chains, the factors derived from the longest chains simulated (PCL100 and PCL125) provide the recommended, converged value for use in future long-chain models. In addition, the dependence of the terminal relaxation time of the chains on the chain length was analysed, with the aim to shed light on the transition from the unentangled to entangled behavior. Using the packing length calculated from the atomistic simulations, we estimated the entanglement molecular weight to be around 4850 g/mol, which compares reasonable well with the reported experimental values estimated from the rheological spectra and with the visual inspection of the N dependence of the terminal relaxation times from the atomistic simulations.
In summary, this work provides a reliable and relatively simple systematic coarse-grained model for linear PCL, offering a valuable tool for extending computational studies to experimentally relevant time and length scales. While the coarse-grained model accurately reproduces local structural correlations and general scaling trends, further refinement may be beneficial for precise quantitative prediction of large-scale, long-range properties in very long chains. This methodology represents a significant step towards facilitating the development of computational approaches for tuning the properties of PCL-based materials, contributing to more sustainable material development. In addition, it may serve as a first step towards quantitative predictions of rheological properties of PCL-based systems with well-controlled composition and architecture. Taking advantage of the well-defined samples in simulations, the confrontation of the simulation results with the experimentally available data and theoretical models may lead to better understanding of the structure-properties-performance relationship, without having to deal with complex features present in experimental PCL samples such as self-nucleation46.
