Knotting behaviour of polymer chains in the melt state for soft-core models with and without slip-springs

The hPF-MD approach and its applications to atomistic and coarse-grained (bio)macromolecular systems have been extensively presented in previous publications [2, 12, 13, 18–25]. Here, we briefly recall the main ideas. In hPF-MD, the intramolecular interactions (bond, angle...) are the same as the standard molecular dynamics (MD) simulations, while pairwise interactions between nonbonded particles in standard MD simulations are transformed into an interaction of a particle with an external potential depending on the density field. For a system composed of two different types of particles, the potential energy in the density field is

$$\begin{align}\hfill W\left[\rho \left(r\right)\right]& =\frac{1}{{\rho }_{0}}\int \mathrm{d}r\left(\frac{{k}_{\text{B}}T}{2}\sum\limits _{i,j}{\chi }_{ij}{\rho }_{i}\left(r\right){\rho }_{j}\left(r\right)\right.\hfill \\ \hfill & \quad +\left.\frac{1}{2\kappa }{\left(\sum\limits _{i}{\rho }_{i}\left(r\right)-{\rho }_{0}\right)}^{2}\right).\hfill \end{align} \tag{ 1 }$$

By applying the saddle point approximation, it is possible to obtain the mean-field external potential ${V}_{i}^{\text{field}}\left(r\right)$ acting on an individual particle (type i) at position r from the functional derivative of the potential energy W[ρ(r)] with respect to the local density:

$$\begin{align}\hfill {V}_{i}^{\text{field}}\left(r\right)& =\frac{\delta W\left[\rho \left(r\right)\right]}{\delta {\rho }_{i}\left(r\right)}=\frac{1}{{\rho }_{0}}\left({k}_{\text{B}}T\sum\limits _{j}{\chi }_{ij}\left(r\right){\rho }_{j}\left(r\right)\right.\hfill \\ \hfill & \quad +\left.\frac{1}{\kappa }\left(\sum\limits _{i}{\rho }_{i}\left(r\right)-{\rho }_{0}\right)\right).\hfill \end{align} \tag{ 2 }$$

Here, the Flory–Huggins parameter χ_ij represents the strength of the mean field interaction between particles of type i and j, ρ₀ is the average number density of the system, κ is the compressibility factor for the system, and ρ_i and ρ_j are number densities of particles of type i and j in the density field at position r, respectively. The forces acting on particles are computed by interpolating the gradients of the external potential on spatial grids. More details about the implementation of density fields and its force calculations in MD simulations can be found in references [2, 22, 26]. Replacing hard-core nonbonded pair interaction by the interaction with a field on a lattice turns MD into an O(N) algorithm, but also makes the nonbonded interactions necessarily soft-core. Polymer chains are thus able to cross each other and fail to show entangled dynamics such as reptation [27]. To remedy this shortcoming, we recently employed slip-springs in hPF-MD simulations of polymer melts, which successfully restored their entangled dynamics [7].

We study a model system of polystyrene melts, which has been introduced in previous publications [12, 13, 21]. The coarse-grained force field of polystyrene melts was originally developed by Qian and co-workers [28]. In this model, two different bead types (R and S) are defined to reproduce the tacticity of the polymer chain with the bead placed at the centre of mass of the repeating unit. Coarse-grained atactic polystyrene chains are created by randomly generating the sequence of R and S beads. Specifically, by using the iterative Boltzmann inversion method [29] to retain probability distribution functions of their all-atom counterparts, coarse-grained bond interactions, namely bonds, angles and dihedral angles, are derived. The coarse-grained nonbonded interactions are computed through the density-functional field, whose update time interval is Δt_field = 1 hPF-MD time step. In this density-functional field, the grid spacing is chosen as 0.55 nm, comparable to the average bond length between connected coarse-grained beads. The latter is l ≈ 0.52 nm for both standard MD and hPF simulations. The Flory–Huggins parameter χ is 0 for a homogeneous polymer melt. To probe the effect of different incompressibility conditions, two incompressibility parameters κ = 0.1 and 0.05 mol × kJ⁻¹ are employed. Their influence on the monomer diameter is derived from the potential of mean force (PMF) between monomers. Specifically, the PMF is computed by Boltzmann inversion of the intermolecular RDF, and the monomer diameter is identified as the length where it equals 0. This definition can be applied in the same way to both MD with pairwise interactions and to hPF. The resulting monomer diameters for MD and the two hPF models with κ = 0.05 and κ = 0.1 are 1.01 nm, 1.24 nm, and 1.27 nm, respectively. All details of the simulated systems such as the number of chains, chain length, box sizes and the total number of slip-springs in hPF simulations (if present at all) are listed in table 1.

All simulations are performed using GPU-accelerated large-scale molecular simulation toolkit (GALAMOST [26]) package with our version of the slip-spring model. The initial configurations are obtained from the reference MD simulations after their density is converged at a temperature of T = 500 K. The equilibration of the hPF systems is achieved by utilizing the fast equilibration procedure reported previously [12]. Specifically, we perform a pure hPF simulation to pre-equilibrate the system up to 4 × 10⁸ steps, with the grid spacing slowly decreasing from l_grid ≈ 1 nm to l_grid ≈ 0.55 nm. The mean squared internal distance is tracked throughout the equilibration step to ensure a complete equilibration [30]. Data are then collected from production simulations continued from these equilibrated systems with l_grid ≈ 0.55 nm. Both hPF simulations with and without slip-springs are run in an NVT ensemble using a Langevin thermostat with a frictional coefficient of 226 g × mol⁻¹ × ps⁻¹.

Slip-springs are modelled by employing a soft Lennard-Jones potential with its parameters mapped to the corresponding hard-core MD systems. The slip-spring bond length follows a Gaussian distribution, whose mean value is close to half of the tube diameter of polystyrene melts $\langle {l}_{\text{ss}}\rangle =\frac{1}{2}{d}_{\mathrm{T}}\approx 3.75\enspace \mathrm{n}\mathrm{m}$ [31]. It should be noted here that the slip-spring bond length in hPF models is longer than that employed in DPD models compared to the respective bead sizes. Slip-spring motion is governed by a Monte Carlo (MC) algorithm. In our example CG-PS systems, the MC motion of slip-springs is activated every Δt_MC = 100 hPF time steps (3 ps). This is close to the relaxation time of a single monomer, which we estimated as the transition time of the monomer's mean squared displacement g₁(t) from ballistic (g₁(t) ∼ t²) to Rouse motion (g₁(t) ∼ t^0.5). The mobility of slip-springs is controlled by a hopping frequency v_hop, which is a control parameter determined by matching corresponding experimental or MD simulation data. In the present model, comparison with MD diffusion coefficients yielded v_hop = 1.25 × 10⁻⁴ ps⁻¹. We note that this implementation of a hopping-frequency-controlled migration is slightly different from our original slip-spring model [5, 6], where MC blocks have a fixed amount of time steps. Nonetheless, both methods successfully restore entanglement dynamics. If a slip-spring reaches a chain end, it attempts a relocation move. Slip-springs are only allowed to relocate from one chain end to another: in a relocation attempt, a new slip-spring is created by connecting a randomly chosen chain end to any other bead within a distance closer than 0.95 × d_T ≈ 7.125 nm. A Metropolis MC trial between the initial and the new slip-spring position determines whether the relocation is successful; the 'loosing' slip-spring is discarded. More details about the implementation of slip-springs in the coarse-grained polystyrene (CG-PS) model can be found in reference [25].

DPD [1, 32, 33] is a mesoscopic simulation method that describes particle interactions in terms of purely pairwise conservative, dissipative, and random contributions. The coarse-grained, soft-core conservative force allows for large integration steps at low computational costs, which makes DPD a popular technique for simulations of soft matter [34]. The coupled dissipative and random forces account for the system's friction and act as a thermostat. Here, we follow the methodology of Groot and Warren [1], where the conservative force between two beads i and j is linear and defined by a cutoff radius r_c and a repulsion parameter a_ij (equation (3)).

$$\begin{equation}{ \overrightarrow {F}}_{ij}^{\mathrm{C}}={a}_{ij}\left(1-\frac{{r}_{ij}}{{r}_{\mathrm{c}}}\right){ \overrightarrow {e}}_{ij}.\end{equation} \tag{ 3 }$$

Polymer chains are modelled using a standard bead-spring model, where DPD beads are connected by a weak Hookean spring in addition to their nonbonded interactions. As a computationally inexpensive answer to the physical shortcomings of soft-core conservative interactions, slip-springs have been employed in DPD simulations of polymer melts [5, 6, 8] and solutions [35, 36] by different groups. In this study, we follow the work of Langeloth et al [5]: initially, a fixed number of slip-springs is distributed between pairs of nonbonded beads, following a similar distance criterion as in the slip-spring hPF model. During the simulation, slip-springs are allowed to migrate along the chains governed by a Metropolis MC criterion. If slip-springs reach the end of a chain, a relocation move to a different, randomly chosen chain end is attempted as described in the previous subsection. Slip-springs are frozen after a number of MC migration attempts and act as fixed bonds in the next set of DPD steps. By performing alternating blocks of MC and DPD steps, entanglement dynamics such as reptation and constraint release are restored. For further details, we refer to the initial work in reference [5].

In addition to unrestricted melts, we present simulations of systems confined by an array of nanotubes (figure 1). Nanotubes are modelled by a purely repulsive potential and placed on a hexagonal grid in the xy-plane. They are infinite in z-direction. The nanotubes are positioned in the centre and on the corners of the simulation box, so that one periodic image contains two tubes. They interact with DPD beads via the conservative potential: in a modified version of equation (3), a_ij and r_ij are replaced by the nanotubes' repulsion parameter a_NT,i and the distance between a bead and their 'surface' r_NT,i . A nanotube's surface is defined by its radius R_NT, which is also used to evaluate its effective excluded volume ${V}_{\text{NT}}=\pi {R}_{\text{NT}}^{2}{l}_{z}$, where l_z is the box size in z-direction. The nanotubes' repulsion parameter a_NT,i is chosen to be twice as repulsive as the regular a_ij . Notably, interactions between beads and nanotubes are still soft-core, which theoretically allows beads to penetrate the tubes. In this case, r_NT,i becomes formally negative. However, we have not observed polymers crossing into the nanotubes, so the repulsion appears to be large enough to keep them out.

Zoom In Zoom Out Reset image size

DPD simulation results are presented in reduced units: time, distance, mass, and energy are given in units of t_DPD, r_c, m_DPD, and k_B T, respectively. If need be, our simulations can be mapped onto a polystyrene model [36]. All DPD parameters are taken from reference [1] for a density of $\rho =3\enspace {r}_{\mathrm{c}}^{-3}$. The repulsion parameters for bead-bead and bead-nanotube interactions are ${a}_{ij}=25\enspace {k}_{\mathrm{B}}T{\times}{r}_{\mathrm{c}}^{-1}$ and ${a}_{\text{NT},i}=50\enspace {k}_{\mathrm{B}}T{\times}{r}_{\mathrm{c}}^{-1}$. Bonds and slip-springs have a force constant of ${k}_{\mathrm{B}}={k}_{\text{SS}}=2\enspace {k}_{\mathrm{B}}T{\times}{r}_{\mathrm{c}}^{-2}$. The integration step is Δt = 0.06 t_DPD. The monomer radius is 1 r_c by definition, and an average bond length of l ≈ 1.21 r_c emerges from the sum of bonded and nonbonded interactions.

Our simulations feature chains of length N = 25, 50, 75, 100, and 150. For the unconfined melts, 1668, 834, 556, 417, and 278 of these chains are simulated in a cubic ${\left(24\enspace {r}_{\mathrm{c}}\right)}^{3}$ box, respectively. Chains under confinement are investigated for different nanotube surface-to-surface distances d_NT. On a hexagonal grid, only two of the parameters d_NT, R_NT, l_x , and l_y can be chosen independently. Here, we fix the nanotube radius at R_NT = 10 r_c and the x-dimension as l_x = 38, 42, 48 and 56 r_c with l_x > l_y . Consequently, ${l}_{y}={l}_{x}/\sqrt{3}=21.939,24.249,27.713,32.332\enspace {r}_{\mathrm{c}}$ and d_NT = l_y − 2R_NT = 1.9, 4.2, 7.7, 12.3 are determined. The z-length is l_z = 21 r_c for all systems. To match a density of $3\enspace {r}_{\mathrm{c}}^{-3}$ in the accessible volume V_box − 2V_NT, the total number of beads in the d_NT = 1.9, 4.2, 7.7 and 12.3 systems is set to 12 900, 25 000, 45 000, and 74 500, respectively. Exceptions are the $\left({d}_{\text{NT}}=4.2,N=\left[75,150\right]\right)$ and $\left({d}_{\text{NT}}=12.3,N=\left[75,150\right]\right)$ systems where the total number of beads is 24 900 and 74 550, respectively, to allow an integer number of chains. Simulations are performed with and without slip-springs for every setup. If present, the number of slip-springs is always 10% of the total number of beads. Details of the slip-spring method can be found in reference [5]. Here, we use the same block sequence lengths of 500 DPD steps and 500 MC migration attempts, a pattern chosen to match the spatial correlation of slip-springs in both blocks [5]. The unconfined melts are equilibrated for 5 × 10⁵ (N = [25, 50, 75]) and 1.5 × 10⁶ (N = [100, 150]) time steps, which is longer than any chains' longest relaxation time. Production runs are performed for 3 × 10⁶ time steps. All confined melts are equilibrated for 10⁶ (N = [25, 50, 75]) and 2 × 10⁶ (N = 100, 150) time steps, while data is extracted from 10⁶ (N = [25, 50, 75]) and 5 × 10⁶ (N = [100, 150]) time steps of production. If not denoted otherwise, errors are the standard deviation of the mean of all chains.

In order to describe a knot in a chain mathematically, the chain is required to be closed. The Alexander polynomials [37] are then used as topological invariants to characterize the knots. In this work, we utilize the Kymoknot software [38] (version 1.0), which enables us to analyse knots in open and closed polymers. For the closure of open chains, the minimally-interfering closure scheme [38, 39] is applied. Within this closure scheme, two distinct ways of connecting the chain ends are compared—they are either connected by direct bridging, or via the closest points of the convex hull of the chain portion. The end-to-end distance of the chain is thus compared with the sum over the distances between the chain ends and their closest points on the convex hull. If the former is smaller, simple bridging is applied, otherwise the chain is closed via the convex hull. This closure scheme has shown to be robust and computationally efficient and leads to the least amount of additional entanglements [39].

The smallest region of the chain that has the same topology as the entire chain after closure is considered as the knotted region. Various search schemes can be used for determining knots in a polymer chain which may return different knotted regions [38]. We use a bottom-up approach, which starts from very short, unknotted portions of the chain. These portions are gradually increased until the physical knot of the chain portion equals the knotting type of the whole chain. The remainder of the chain is unknotted. The knot size is then defined as the number of beads within the knotted region.

The kinetics of the constant transition between the knotted and the unknotted state can be analysed as well. For a single chain i, we introduce a state function h_i (t), which equals 1 if chain i contains at least one knot and is 0 otherwise. We then define a correlation function C(t) for the knotting dynamics using this state function:

$$\begin{equation}C\left(t\right)=\frac{1}{M}\sum\limits _{i}\frac{\langle \left({h}_{i}\left(t\right)-\langle h\rangle \right)\left({h}_{i}\left(0\right)-\langle h\rangle \right)\rangle }{{\langle \left({h}_{i}\left(t\right)-\langle h\rangle \right)\rangle }^{2}}.\end{equation} \tag{ 4 }$$

Here, we average over all M chains in the system. The relaxation time of a knot τ_K is obtained as the integral of the correlation function C(t)

$$\begin{equation}{\tau }_{\text{K}}={\int }_{0}^{\infty }C\left(t\right)\mathrm{d}t.\end{equation} \tag{ 5 }$$

In this section, we discuss the differences in the knotting behaviour of our hard-core standard MD model and the soft-core hPF model derived from it. In the hPF model, the softness of nonbonded interactions can be controlled by tuning the compressibility parameter κ. Specifically, a smaller value of κ gives a less compressible system, in which less particle overlapping occurs. The chain dimensions of all polystyrene models are shown in figure 2. Among these models, the chain structures are found to be consistent, indicating that the soft-core and hard-core natures do not necessarily alter polymer structural properties on the chain scale. Both the ensemble-averaged squared end-to-end distance $\langle {R}_{\text{ete}}^{2}\rangle$ and radius of gyration $\langle {R}_{\text{g}}^{2}\rangle$ show the expected linear scaling with the number of monomer N (dashed lines in figure 2) [40]. It should be noted that the compressibility parameter κ must have an impact on the structural behaviour of polymer chains; here, the structural deviation for the employed range of κ is within 3%.

Zoom In Zoom Out Reset image size

Next, we examine the knotting probability and knot size in both soft-core and hard-core models using the knotting analysis presented in subsection 2.3. As shown in figure 3(a), the knotting probability P_K of hPF models is profoundly higher than the MD models over the whole range from N = 100 to N = 750 monomers. Moreover, the knotting probability P_K of hPF models continuously increases with the polymer chain length N. A linear fit yields slopes of k ≈ 3.9 × 10⁻⁴ and k ≈ 2.3 × 10⁻⁴ for the hPF models with κ = 0.1 and κ = 0.05, respectively. This indicates that the knotting probability in hPF models with κ = 0.1 is almost 50% higher than that with κ = 0.05. It is noted that the linearity of the knotting probability with the chain lengths in hard-core MD simulations is not as evidently seen as that in the soft-core models. This probably comes from the rather small low statistical value of P_K in MD simulations. Nevertheless, the knotting probability in hard-core MD models is sufficiently low in comparison to the soft-core models, implying that the degree of overlapping plays an important role. We further investigate the knot size where we focus on the dominating simple trefoil (3₁) knot. As illustrated in figure 3(b), the average size detected in the MD model is significantly larger than that in hPF models. Linear fits reveal that approximately 44% of all monomers of any chain are involved in one trefoil 3₁ knot in the MD model. In contrast, only ∼5% and ∼8% of the monomers per chain form a 3₁ knot in the hPF models with κ = 0.1 and κ = 0.05, respectively. With the similar characteristic ratio of the polystyrene melt in hPF and MD models, the overestimated knotting probability and underestimated knot size are expected to be intimately related to the soft-core interactions, implying that the overlapping of particles facilitates the localization and occurrence of the formation of a trefoil 3₁ knot. Recently, Meyer et al [16] studied the knotting properties of generic polymer models with various characteristic ratios using MD simulations and compared their results with the random walk (RW) polymer model. The RW model is found to overestimate the knotting probability and underestimate the knot size in comparison with the simulation model with hard-core repulsive interactions, consistent with our observations in hPF and MD models.

Zoom In Zoom Out Reset image size

In this section, we study the influence of slip-springs, which are initially designed to mimic the chain entanglements for dynamical properties, on the knotting behaviour of soft-core hPF and DPD models. The typical reptation behaviour of polymer melts was recovered by the introduction of slip-springs to these soft-core models in prior studies [5, 7]. The effect of slip-springs on the conformational properties of polymer chains (e.g. the radius of gyration) in the melt state was found to be negligible. How slip-springs, as artificial ad hoc interactions, affect the self-entanglements (knottedness) of polymer chains is still unknown.

3.2.1. Hybrid-particle-field simulations

We analyse the knotting behaviour in terms of the probability P_K to carry at least one knot in hPF models of polystyrene melts with and without slip-springs at various chain lengths. As displayed in figure 4(a), P_K of the hPF chains with slip-springs is in good agreement with that without slip-springs at all investigated chain lengths. Figure 4(b) shows the average trefoil 3₁ knot size for chain lengths between N = 100 and 750 in hPF models with and without slip-springs. Consistent with previous observations for hPF, we find the 3₁ knot size in hPF simulations with slip-springs to increases with increasing chain length. No significant difference is observed between systems with and without slip-springs. Both of these analyses suggest that slip-springs indeed have a negligible effects on the self-entanglements (knottedness) of polymer chains in hPF simulations within the tested parameter space.

Zoom In Zoom Out Reset image size

3.2.2. Dissipative-particle dynamics simulations of polymer melts

DPD simulations of polymer melts with and without slip-springs are investigated for their probability P_K to carry at least one knot per chain, as well as the size of the simple 3₁ trefoil knot. The results are shown as a function of the number of beads N in figures 5(a) and (b), respectively. They are comparable with those of the hPF model in terms of slip-springs, which show no significant impact on the knotting probability or the knot size. The knotting probability as well as the average knot size of the 3₁ trefoil knots both increase roughly linearly with the chain length.

Zoom In Zoom Out Reset image size

3.2.3. Dissipative-particle dynamics simulations of confined polymer melts

We finally study the knottedness of DPD chains confined by a regular array of nanotubes for different distances d_NT between two nanotubes and thus different degrees of confinement. Generally, confinement increases the polymer's compactness. An example plot of the squared radius of gyration $\langle {R}_{\text{g}}^{2}\rangle$ as a function of the chain length N is given as figure A1 in the supplementary information (https://stacks.iop.org/JPCM/33/244001/mmedia) for the strongest confinement (d_NT = 1.9) and the unconfined systems. $\langle {R}_{\text{g}}^{2}\rangle$ shrinks by roughly 35–40% for d_NT = 1.9, however, the scaling of $\langle {R}_{\text{g}}^{2}\rangle$ with N − 1 is barely affected. Further findings concerning the static and dynamic polymer properties will be reported elsewhere. The probability P_K of chains of length N to carry at least one knot is given in figure 6(a). Consistent with the observations made for unconfined melts, simulations with and without slip-springs are barely distinguishable. As the distance between the nanotubes d_NT decreases, the confined chains assume more compact structures, which is associated with a higher knottedness. For a chain length of N = 150, P_K is 18% in the pure melt, while it increases to 39% for a strongly confined system with a nanotube distance of d_NT = 1.9. For short chains (N ⩽ 50), the radius of gyration is of similar size as the range of the interstice region between three nanotubes. We do, however, still observe a distinct effect: for N = 25, we see an increase in P_K from 0.8% (unconfined) to 2.0% (d_NT = 1.9), which is comparable to the increase detected for longer chains. We briefly study the knottedness of a chain as a function of its confinement: a plot of the knotting probability P_K against the squared radius of gyration $\langle {R}_{\text{g}}^{2}\rangle$ for various chain lengths and degrees of confinement can be found in figure A2 (supplementary information). For a fixed confinement, P_K and $\langle {R}_{\text{g}}^{2}\rangle$ show the expected linear relation, as they are both linked by the chain length (${P}_{\text{K}}\propto N\propto \langle {R}_{\text{g}}^{2}\rangle$). For a fixed chain length N, on the other hand, the knotting probability increases with decreasing $\langle {R}_{\text{g}}^{2}\rangle$. However, we did not find a trivial relation between them, e.g. by a power law. This is likely due to the alignment of chains parallel to the nanotubes, and the therefore anisotropic nature of the confinement.

Zoom In Zoom Out Reset image size

The average 3₁ knot size for all systems containing nanotubes (figure 6(b)) reveals a rather weak confinement effect. The influence of slip-springs is negligible, and the slight deviations for longer chains (N ⩾ 100) are likely of statistical nature. Interestingly, the decreasing nanotube distance appears to lead to an increase in the knot size, where a maximum is reached for d_NT = 4.4. By analysing the spatial components of the radius of gyration, we find that chains are elongated parallel to the nanotubes. Therefore, generally, chains tend to form larger knots when brought into confinement. However, the size of the interstice region between three nanotubes decreases with decreasing nanotube distance. The size of the radius of gyration of longer chains (N ⩾ 50) thus competes with the diameter of this interstice region, especially for short nanotube surface-to-surface distances. In the case of d_NT = 1.9, a larger fraction of chains are spread over multiple interstitial regions. In this case, knots might be formed in parts of the chain that are confined within the interstitial region, leading to smaller knot sizes.

We lastly compare the observed knottedness of our DPD and hPF chains with results reported for different limiting models. We focus on the knotting probabilities of systems without slip-springs. As discussed above, the same observations apply for the slip-spring systems. Recently, Meyer et al [16] investigated the knotting probabilities in self-avoiding model polymer chains with an adjustable bending potential and compared them to those of a RW. We follow their approach to compare the knotting probability of our models on the chains' Kuhn scale: by utilizing the ratio N^⋆ = N/C_∞, we normalize both the DPD and hPF chain lengths N to the number of Kuhn segments per chain. Here, C_∞ is the characteristic ratio and derived from the bond length l between two beads as ${C}_{\infty }=\langle {R}_{\text{ete}}^{2}\rangle /\left({l}^{2}{\times}N\right)$. The results are shown in figure 7 along with the data reported in reference [16] for their fully flexible chain and the RW. The fully flexible chain consists of purely repulsive Lennard-Jones beads and a harmonic bond known to reproduce the same chain structures as finitely extensible (FENE) bonds. This, along with a reported characteristic ratio of C_∞ = 2.1, makes it very similar to the standard Kremer–Grest (KG) model (C_∞,KG = 1.7) [41, 42]. We thus refer to it as a KG-like chain. Coarse-grained MD chains, hPF chains with κ = 0.1 and 0.05, and DPD chains have characteristic ratios of C_∞ = 1.55, C_∞ = 1.49, C_∞ = 1.45 and C_∞ = 1.45, respectively, which were derived from their bond lengths of 0.52 nm (hPF) and 1.21 r_c (DPD). We find the knotting probability of our hard-core MD simulations to be in the range of, but lower than, that of KG-like chains. The difference probably arises from the different number densities of these two models, as higher densities in polymer systems are known to favour higher knotting probabilities [16]. However, a comparison of both models' number densities is not straightforward and depends on the spatial quantity utilized for mapping. The knotting probability increases for the κ = 0.05 and κ = 0.1 hPF models and the DPD model, respectively, with those of the DPD chain being somewhat lower than for the RW. We note that the Kuhn-scale-reduced DPD model should be taken with a grain of salt, as the mapping of DPD chains usually happens on even coarser scales; however, it serves well for the comparison between models. The P_K of the less compressible hPF model (κ = 0.05) lies between those of the more flexible hPF (κ = 0.1) and DPD chains and the hard-core MD model. These results are generally consistent with the findings of Meyer et al: comparing their fully-flexible KG-like chain and the RW, they identify the knotting behaviour to be closely connected to the models' bond–bond correlation functions: since the formation of a knot requires a negative correlation on the respective length scale, the bond–bond correlations as well as their fluctuations determine at which scale and how likely knots are formed [16]. For hard-core models, knots are thus rather rare. In the soft-core limit, RWs allow many and strongly localized knots. Meyer et al conclude that 'real' polymer chain melts, such as the KG-like chain, are rather poorly described by RWs on a local scale [16]. It becomes apparent from figure 7 that soft-core models, especially our DPD system, are much closer to the RW on all scales. Consequently, their knotting probability is much closer to a RW than to a KG-like chain. However, weak local repulsion still exists, and the full knotting probability of the RW is not yet reached. The deviations from a local random-walk behaviour become stronger with increasing excluded-volume effects when going from DPD to the κ = 0.1 and κ = 0.05 hPF models and the hard-core MD model, which all have bending potentials. This causes a decrease of the knotting probability. We note that, if need be, the local structure of hard-core MD or KG-like chains can at least in parts be restored to the hPF models by tuning the compressibility parameter κ.

Zoom In Zoom Out Reset image size

In addition to the static properties related to knotting, we investigate the dynamics of the knotted-to-unknotted transition in the hPF model. Figure 8 displays the autocorrelation function (ACF) of this transition in hPF simulations of polystyrene melts without (a) and with slip-springs (b). In both models, the ACFs exhibit no chain-length dependence. Moreover, we find that the ACFs decay to zero in about ∼10 ps. In order to quantitatively characterize the relaxation time of knot formation dynamics, we fit the ACF by a stretched exponential function

$$\begin{equation}C\left(t\right)=\mathrm{exp}\left[-{\left(\frac{t}{{\tau }^{{\ast}}}\right)}^{\beta }\right].\end{equation} \tag{ 6 }$$

One example is shown in figure 9 for hPF simulations (κ = 0.1) of N = 100 and N = 750 with (solid lines) and without slip-springs (dashed lines). The relaxation time τ_knot is then calculated as

$$\begin{equation}{\tau }_{\text{knot}}={\int }_{0}^{\infty }C\left(t\right)\mathrm{d}t=\frac{{\tau }^{{\ast}}}{\beta }{\Gamma}\left(\frac{1}{\beta }\right),\end{equation} \tag{ 7 }$$

where Γ represents the complete gamma function. The calculated relaxation times are summarized in table 2. The relaxation time τ_knot evidently does not depend on the chain length. The knot relaxation happens on a time scale comparable to the time which a monomer needs to diffuse its own diameter. The same results are observed for our DPD models (not shown). Indeed, these results of fast knotting dynamics are not surprising. In the soft-core hPF model, the nonbonded interactions are so soft that the bonds can cross each other. As illustrated in the inset of figure 8(b), the knotted-to-unknotted transition can easily take place at the monomer length and time scale. This transition is essentially insensitive to the chain length and effectively dictated by the monomer relaxation time. A similarly fast relaxation of knotted structures in soft-core models has qualitatively been observed by Meyer et al, who also reported a change in knotting probabilities on drastically shorter time scales than the structural chain relaxation. Moreover, slip-springs in our systems are expected to alter the polymer dynamics only above the entanglement length (N_e) and time (τ_e) scales. This is consistent with the observation in table 2 that the knot formation dynamic is genuinely unaffected by the addition of slip-springs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size