Addressing the Folding of Intermolecular Springs in Particle Simulations: Fixed Image Convention

Abstract

Mesoscopic simulations of long polymer chains and soft matter systems are conducted routinely in the literature in order to assess the long-lived relaxation processes manifested in these systems. Coarse-grained chains are, however, prone to unphysical intercrossing due to their inherent softness. This issue can be resolved by introducing long intermolecular bonds (the so-called slip-springs) which restore these topological constraints. The separation vector of intermolecular bonds can be determined by enforcing the commonly adopted minimum image convention (MIC). Because these bonds are soft and long (ca 3–20 nm), subjecting the samples to extreme deformations can lead to topology violations when enforcing the MIC. We propose the fixed image convention (FIC) for determining the separation vectors of overextended bonds, which is more stable than the MIC and applicable to extreme deformations. The FIC is simple to implement and, in general, more efficient than the MIC. Side-by-side comparisons between the MIC and FIC demonstrate that, when using the FIC, the topology remains intact even in situations with extreme particle displacement and nonaffine deformation. The accuracy of these conventions is the same when applying affine deformation. The article is accompanied by the corresponding code for implementing the FIC.

1. Introduction

Understanding the structure–property relations in polymers and soft matter systems is key to designing advanced materials with tailor-made properties for applications in biology [1], nanocomposites [2,3,4,5], rheology [6,7,8], and more. Polymeric materials exhibit a rich spectrum of relaxation mechanisms, the longest of them spanning the millisecond to second regime, depending on the chain molar mass. The latter phenomena can be investigated with mesoscale simulation techniques such as Brownian dynamics [9] and dissipative particle dynamics [10,11], which treat large collections of atoms as single coarse-grained sites, in that manner drastically decreasing the degrees of freedom and affording huge savings in computational time.

Coarse-grained polymer chains are, however, prone to unphysical intercrossing, thus violating the network topology and yielding artificially faster dynamics and lower viscosity. Much research has been devoted to restoring the proper viscoelastic properties of coarse-grained polymers by enforcing the topological constraints in an artificial manner, such as the TWENTANGLEMENT model [12], and the slip-link [13,14] and slip-spring (SS) [15,16,17] models, all of them resting on the assumptions of reptation theory [18,19,20]. In the slip-spring model, the entanglements are modeled with entropic springs which are generated/destroyed at the chain ends and are able to slide along the chain contours during the course of the simulation, modeling in this way the mechanism of reptation. The slip-spring model has been particularly successful for restoring the viscoelastic properties of linear [8,17,21,22], star [23,24], and cross-linked [25,26,27] polymers, in bulk [8,17] and interfacial [5,21,28] geometries under quiescent and nonequilibrium [8,16,29,30] conditions.

A common boundary condition for investigating homogeneous samples is the periodic boundary condition (PBC), where the domain is treated as infinite along the periodic directions [9]. This is achieved by replicating the central box along these directions; these replicas are called box images. It is convenient to allow the particles to diffuse outside the central box in order to compute their transport properties. In doing so, the relative distance between two particles i and j (which may lie in different box images) can be corrected by enforcing a suitable convention such as the minimum image convention (MIC).

According to the MIC, each particle i interacts with the image of particle j that is closest to it, in such a way that the absolute projections of their separation vector on the basis vectors of the simulation box do not exceed half the corresponding edge lengths of the box. The MIC is applicable to non-orthogonal boxes as well, after adjusting for the shear displacement of the periodic images [9,31]. Care must be exercised, however, in the case of particle pairs connected by soft springs; in a situation where the spring projection on one or more of the coordinate axes exceeds half the box edge length along that axis, the separation vector is redefined with respect to the (new) closest pair of images of the particles, resulting in a violation of the system topology.

This issue can become especially prevalent in simulations of coarse-grained polymers that have been subjected to extreme deformation, due to the inherent softness of the intra- and intermolecular effective bonds (springs) involved. For example, the equilibrium length of the intermolecular springs that are used to mimic the entanglements (slip-springs [15,16]) is ca 3–20 nm [32]. Subjecting the system to strong flows can overextend the aforementioned bonds significantly and lead to topology violations. In addition, the stability of the simulation is severely affected during extensional flow experiments [33,34], due to the decreased cross-section of the sample which leads to frequent instabilities and topology violations. Another drawback of the MIC is that it entails considerable computational cost due to the calculation of remainders [9], whereas choosing an optimal algorithm depends on several factors such as the compiler optimization and processor architecture [35].

We propose an alternative convention for determining the separation vectors of molecular bonds, which we call the fixed image convention (FIC). In FIC, the MIC is enforced only once during the formation of the bond and the minimum separation vector is stored in memory, in terms of a shift vector. In subsequent iterations, the separation vector is calculated simply by adding a shift vector to the absolute separation vector. The FIC is applicable to boxes with either constant or varying size (i.e., simulations in the isothermal-isobaric ensemble). Note, however, that subjecting the system to affine deformation requires deforming the shift vectors as well. Given that the distance between the images is fixed, the system topology remains intact even for very large deformations. In general, this operation is much faster than the usual MIC, especially when working with non-orthogonal boxes.

The article is structured as follows: Section 2.1 illustrates the coordinate system used in this work, and Section 2.2 reports the formulation of the various conventions for determining the separation vector. Section 3.1 compares the accuracy of the MIC and the FIC under conditions of extreme displacement and affine/non-affine deformation. Section 3.2 assesses the stability of the conventions under random displacement and strong shear deformation. Section 3.3 discusses the efficiency of the MIC and the FIC. Finally, Section 4 summarizes the main findings of the work. The code for conducting the comparisons between the MIC and FIC is provided via a GitHub repository [36].

2. Methods

2.1. Coordinate System and Constraints

The edge vectors of a parallelepiped simulation box can be described according to Equation (1):

$$\begin{array}{l}{l}_1=(L_x,0,0)\\ l_2=({\Delta }_{xy},L_y,0)\\ l_3=({\Delta }_{xz},{\Delta }_{yz},L_z)\end{array}$$

(1)

with Δ_xy, Δ_xz, and Δ_yz indicating the displacement of the y/z/z face along the x/x/y axis of an originally orthogonal box with dimensions L_x, L_y, and L_z. In addition, the shearing components can be constrained according to Equation (2):

$$-\frac{L_{\alpha }}{2}\le {\Delta }_{\alpha \beta }\le \frac{L_{\alpha }}{2}$$

(2)

Adding or subtracting multiples of the box edge vectors does not affect the magnitude of the separation vectors when enforcing suitable minimum image conventions [9]. The aforementioned restrictions impose no loss of generality; any set of arbitrary non-collinear edge vectors can be rotated, inverted, and translated in order to conform to these constraints. For convenience, we will refer to the configuration of a box with the vector:

$$\mathit{L}=\left(L_x,L_y,L_z,{\Delta }_{xy},{\Delta }_{xz},{\Delta }_{yz}\right)$$

(3)

2.2. Conventions of Separation Vectors

2.2.1. Absolute Distance (AD)

Let r_i, r_j denote the absolute position of two particles i, j. By absolute, it is meant that the particles may lie anywhere in the Cartesian domain regardless of the bounds of the box. Their absolute separation vector is defined according to Equation (4).

$${\mathit{r}}_{ij}^{\mathrm{AD}}={\mathit{r}}_j–{\mathit{r}}_i$$

(4)

2.2.2. Minimum Image Convention (MIC)

In orthogonal boxes, the calculation of the minimum separation vector can be performed according to the MIC as follows:

$${\mathit{F}}_{\mathrm{orth}}^{\mathrm{MIC}}\left({\mathit{r}}_{ij}^{\mathrm{AD}},\mathit{L}\right)=\left[\begin{array}{c}{x}_{ij}^{\mathrm{AD}}-L_x\mathrm{nint}\left(x_{ij}^{\mathrm{AD}}/L_x\right)\\ y_{ij}^{\mathrm{AD}}-L_y\mathrm{nint}\left(y_{ij}^{\mathrm{AD}}/L_y\right)\\ z_{ij}^{\mathrm{AD}}-L_z\mathrm{nint}\left(z_{ij}^{\mathrm{AD}}/L_z\right)\end{array}\right]$$

(5)

with nint denoting the nearest integer rounding function, and ${\mathit{F}}_{\mathrm{orth}}^{\mathrm{MIC}}$ a vector-valued function which returns the minimum image of a vector ${\mathit{r}}_{ij}^{\mathrm{AD}}$ subject to the dimensions of an orthogonal periodic box, L.

In triclinic boxes, the minimum image vector is determined by applying the Lees–Edwards boundary condition [31], which requires a prerequisite step. Initially, one has to correct for the displacement of the periodic images [9]:

$${\mathit{r}}_{ij}^{\mathrm{cor}}={\mathit{F}}_{}^{\mathrm{cor}}\left({\mathit{r}}_{ij}^{\mathrm{AD}},\mathit{L}\right)=\left[\begin{array}{c}{x}_{ij}^{\mathrm{AD}}-{\Delta }_{xy}\mathrm{nint}\left(y_{ij}^{\mathrm{AD}}/L_y\right)-{\Delta }_{xz}\mathrm{nint}\left(z_{ij}^{\mathrm{AD}}/L_z\right)\\ y_{ij}^{\mathrm{AD}}-{\Delta }_{yz}\mathrm{nint}\left(z_{ij}^{\mathrm{AD}}/L_z\right)\\ z_{ij}^{\mathrm{AD}}\end{array}\right]$$

(6)

Subsequently, the minimum separation vector is retrieved by inputting ${\mathit{r}}_{ij}^{\mathrm{cor}}$ to Equation (5):

$${\mathit{r}}_{ij}^{\mathrm{MIC}}={\mathit{F}}_{}^{\mathrm{MIC}}\left({\mathit{r}}_{ij}^{\mathrm{AD}},\mathit{L}\right)={\mathit{F}}_{\mathrm{orth}}^{\mathrm{MIC}}\left({\mathit{r}}_{ij}^{\mathrm{cor}},\mathit{L}\right)$$

(7)

2.2.3. Fixed Image Convention (FIC)

The FIC entails storing the difference between the minimum and the absolute separation vectors during the formation of the bond, in terms of the shift vector:

$${\mathit{r}}_{ij}^{\mathrm{shift}}={\mathit{F}}^{\mathrm{MIC}}\left({\mathit{r}}_{ij}^{\mathrm{AD}},\mathit{L}\right)-{\mathit{r}}_{ij}^{\mathrm{AD}}$$

(8)

In the subsequent iterations, the separation vector is calculated simply by adding the absolute separation vector and the (constant) shift vector:

$${\mathit{r}}_{ij}^{\mathrm{FIC}}={\mathit{r}}_{ij}^{\mathrm{AD}}+{\mathit{r}}_{ij}^{\mathrm{shift}}$$

(9)

Deforming the box entails some additional operations. Let F be a linear transformer describing the mapping between a reference ($v^{\mathrm{ref}}$) and a transformed (v) vector:

$$v=\mathit{F}\cdot v^{\mathrm{ref}}$$

(10)

with “ref” denoting the vectorial quantity before the deformation; details regarding the analytical calculation of F are reported in the Appendix A. Applying F to a parallelepiped transforms its edge vectors according to Equation (11).

$$\left[l_1 l_2 l_3\right]=\left[\mathit{F}\cdot l_1^{\mathrm{ref}} \mathit{F}\cdot l_2^{\mathrm{ref}} \mathit{F}\cdot l_3^{\mathrm{ref}}\right]$$

(11)

Depending on the application [37], the coordinates of the constituent particles may $\left({\mathit{r}}_i^{}=\mathit{F}\cdot {\mathit{r}}_i^{\mathrm{ref}}\right)$ or may not $\left({\mathit{r}}_i^{}={\mathit{r}}_i^{\mathrm{ref}}\right)$ be transformed according to $\mathit{F}$. According to the type of applied deformation, we have the following cases:

• Affine deformation: the transformation is applied to both the box and the constituent particles. In this case, the shift vectors are transformed accordingly:

  $${\mathit{r}}_{ij}^{\mathrm{shift}}=\mathit{F}\cdot {\mathit{r}}_{ij}^{\mathrm{shift},\mathrm{ref}}$$

  (12)
• Nonaffine deformation: the transformation is applied only to the box, and the shift vectors remain unaffected:

  $${\mathit{r}}_{ij}^{\mathrm{shift}}={\mathit{r}}_{ij}^{\mathrm{shift},\mathrm{ref}}$$

  (13)

Finally, in order to restart a simulation with the FIC, the shift vectors should be exported to the hard drive alongside the other properties of the simulation (box dimensions, coordinates, etc.) and be read during the restart procedure.

3. Results

3.1. Comparisons between MIC and FIC

We will first assess the accuracy of the MIC and FIC for the six representative scenarios (scn) that are depicted in

Table 1. Parameters of scenarios A–F. The columns under the operation header denote the components of the maximum applied strain tensor (Equation (14)) and whether the transformation has been applied to the box and/or the particle positions. The columns under the reference/deformed header correspond to the dimensions of the box and the magnitude of the separation vectors before/after the deformation *.

scn	Operation				Reference					Deformed
	εyy	εxy	box	coord	${\mathit{L}}_{\mathit{y}}^{\mathbf{r}\mathbf{e}\mathbf{f}}$	${\mathbf{\Delta }}_{\mathit{x}\mathit{y}}^{\mathbf{r}\mathbf{e}\mathbf{f}}$	${\mathit{r}}_{\mathit{i}\mathit{j}}^{\mathbf{A}\mathbf{D},\mathbf{r}\mathbf{e}\mathbf{f}}$	${\mathit{r}}_{\mathit{i}\mathit{j}}^{\mathbf{M}\mathbf{I}\mathbf{C},\mathbf{r}\mathbf{e}\mathbf{f}}$	${\mathit{r}}_{\mathit{i}\mathit{j}}^{\mathbf{F}\mathbf{I}\mathbf{C},\mathbf{r}\mathbf{e}\mathbf{f}}$	${\mathit{L}}_{\mathit{y}}^{}$	${\mathbf{\Delta }}_{\mathit{x}\mathit{y}}^{}$	${\mathit{r}}_{\mathit{i}\mathit{j}}^{\mathbf{A}\mathbf{D}}$	${\mathit{r}}_{\mathit{i}\mathit{j}}^{\mathbf{M}\mathbf{I}\mathbf{C}}$	${\mathit{r}}_{\mathit{i}\mathit{j}}^{\mathbf{F}\mathbf{I}\mathbf{C}}$
A	1	0	F	T	10	0	12.5	2.5	2.5	10	0	25	5	15
B	1	0	F	T	10	4	12.5	4.71	4.71	10	4	25	5.38	15.52
C	1	0	T	F	10	0	12.5	2.5	2.5	20	0	12.5	7.5	2.5
D	1	0	T	T	10	0	12.5	2.5	2.5	20	0	25	5	5
E	0	1	T	F	10	0	12.5	2.5	2.5	10	5	12.5	2.5	2.5
F	0	1	T	T	10	0	12.5	2.5	2.5	10	5	17.67	3.53	3.53

* In all cases: ε_xx = ε_yx = 0; $L_x^{\mathrm{ref}}=L_y^{\mathrm{ref}}=L_x=10$; ${\mathit{r}}_i^{\mathrm{ref}}={\mathit{r}}_i^{}=\left(0, 0\right)$, ${\mathit{r}}_j^{\mathrm{ref}}={\mathit{r}}_{ij}^{\mathrm{AD},\mathrm{ref}}=\left(0,12.5\right)$.

. In each case, we apply a linear transformer to the box and/or the constituent particles which is related to the (nonsymmetrical) engineering strain tensor as follows:

$$\mathit{\epsilon }\doteq \frac{\partial }{\partial v^{\mathrm{ref}}}\left(v-v^{\mathrm{ref}}\right)=\mathit{F}-I$$

(14)

Depending on whether the deformation has been applied to the box and/or the particle coordinates, we have the following cases: displacement (box, coord) = (F, T), nonaffine deformation (box, coord) = (T, F), and affine deformation (box, coord) = (T, T), with box and coord denoting Booleans.

In all cases, the particle i resides at the origin of the box; hence, it remains unaffected. The particle j has been placed in the first image of the box along the y-axis, mimicking in this manner the formation of a bond between two chains which reside in different periodic images. For simplicity, the calculations have been conducted in 2D; calculations in 3D yield similar qualitative findings.

3.1.1. Particle Displacement in Rigid Boxes

First, we will examine the effect of (over)extending an intermolecular bond connecting particles which lie in different periodic images of orthogonal boxes as a function of displacement.

Figure 1a,b illustrate results from scenario A (see Table 1) where the particle j has been displaced with respect to particle i along the y-axis. In particular, Figure 1a illustrates snapshots of the particle positions and the separation vectors ${\mathit{r}}_{ij}^{\mathrm{MIC}}$ and ${\mathit{r}}_{ij}^{\mathrm{FIC}}$ connecting particle i with the image of particle j, hereafter referred to as $j^{\prime }$. Figure 1b depicts the magnitude of the separation vectors (${\mathit{r}}_{ij}^{\mathrm{MIC}}$, ${\mathit{r}}_{ij}^{\mathrm{FIC}}$, and ${\mathit{r}}_{ij}^{\mathrm{AD}}$) as a function of the displacement in box units, i.e., Δ_y/L_y = ε_yy. For more details, please refer to the caption of Figure 1.

For small displacements ($\left|{\Delta }_y\right|<0.5L_y$), the separation vectors according to the MIC and FIC are the same. For larger displacements ($\left|{\Delta }_y\right|\ge 0.5L_y$), the MIC becomes unstable since ${\mathit{r}}_{ij}^{\mathrm{MIC}}$ changes sign; i.e., the image of the particle in Figure 1a (MIC) hops to the −y image of the box. Instances such as this should be avoided during the course of practical simulations, since they violate the system topology. The FIC, on the other hand, yields the desired behavior, i.e., the separation vector increases proportionally with displacement.

Figure 1c,d depict the same evaluations but for a box that has been sheared along the x-axis (scenario B in Table 1). The behavior of the MIC and FIC is qualitatively the same as in the previous case; ${\mathit{r}}_{ij}^{\mathrm{MIC}}$ swaps signs and the image of the particle hops from the −x to the −y box image when $\left|{\Delta }_y\right|\ge 0.5L_y$, whereas ${\mathit{r}}_{ij}^{\mathrm{FIC}}$ increases steadily with the displacement.

3.1.2. Affine and Nonaffine Deformation (Varying Box Size)

In the following, we will examine the accuracy of the MIC and FIC during affine and nonaffine deformation experiments.

Figure 2 illustrates results from elongation experiments according to the parameters of scenarios C (panels a, b) and D (panels c, d) in Table 1; panels a, c illustrate snapshots of the system, and panels b, d show the magnitude of the separation vectors during the deformation according to the AD, MIC, and FIC.

According to Figure 2a,b, enforcing the FIC under nonaffine deformation conditions leaves the separation vector unchanged. Given that the relative positions of the particles are unchanged, it makes sense that the separation vector should remain the same. In contrast, enforcing the MIC results in some interesting behavior: ${\mathit{r}}_{ij}^{\mathrm{MIC}}$ decreases with increasing strain until L_y = r_j, after which it increases indefinitely. Regardless, however, deforming the box while keeping the coordinates of the particles unchanged is not a well-defined process; thus, depending on the application, one may choose to apply either the MIC or the FIC.

Figure 2c,d present results from the same experiment, but with the difference that the particle coordinates are transformed along with the box (scn D in Table 1). In this situation, the MIC and FIC yield identical results. The projections of the separation vector change proportionally with the box dimensions; hence, they cannot become larger than half of the edge vectors of the box and swap box images.

Figure 3 illustrates results from the shearing deformations according to scenarios E (panels a, b) and F (panels c, d) in Table 1. Note that at ε_xy = 0.5, Δ_xy flips sign according to the constraint in Equation (2). Similarly, with the nonaffine elongation experiments, enforcing the FIC leaves the separation vector unchanged, whereas the MIC results in complicated behavior. Under affine shear deformation conditions, the MIC and FIC result in identical behavior (Figure 3c,d).

3.2. Stability

In this section, we will assess the performance of the MIC and FIC against bead-spring molecular mechanics models under conditions of extreme deformation.

Figure 4a (left) depicts a periodic face-centered cubic (FCC) lattice made by 5 × 5 × 5, four-particle unit cells with lattice constant, a_cell = 10 Å. Each particle interacts with its 12 nearest neighbors (n_bond = 3000 total interactions) via a harmonic potential of the form:

$$E_{\mathrm{b}}\left(\mathit{r}\right)=k_{\mathrm{b}}{\left(\mathit{r}-{\mathit{r}}_{\mathrm{b}}\right)}^2$$

(15)

where k_b = 10⁻⁵ kJ/(molŲ), and ${\mathit{r}}_{\mathrm{b}}=a_{\mathrm{cell}}/\sqrt{2}$ is the first neighbor distance. The particles at the periodic boundaries interact with their nearest images subject to the MIC or the FIC; thus, the energy of this reference state is zero.

The particle coordinates are randomized and the configuration is subjected to an energy minimization via a conjugate gradient algorithm [38,39] that has been implemented in the EMSiPoN code [8,17,25,28] (see Appendix B). The evolution of the normalized energy during the course of the minimization is illustrated in Figure 4b. The FIC yields a higher energy at the start of the optimization (Figure 4b, step = 0) because it describes the topology of the overextended bonds properly. The initial energy is lower when enforcing the MIC because the topology is violated since the separation vectors are redefined with respect to the (new) nearest image particles. Due to these topology violations, the MIC performs rather poorly, and the system cannot return to its initial state. In contrast, the FIC keeps the topology intact, and the system is able to return to its initial state after the minimization.

Figure 5a (panel ε_xy = 0) depicts a periodic 5 × 5 square lattice with a_lat = 2 Å. Each particle interacts with its four closest neighbors (n_bond = 50) according to Equation (15) with k_b = 1 kJ/mol and r_b = a_lat. The lattice is gradually sheared along the x-axis for ε_xy up to 5; the evolution of the normalized energy during this operation is shown in Figure 5b.

For shear deformations up to ε_xy = 2.5, the energy increases parabolically, regardless of the image convention. After this threshold, the ends of the vector ${\mathit{r}}_{ij}^{\mathrm{MIC}}$ attach to the image particles in the +x box image; see green bonds in Figure 5a. As a result, the energy starts to decrease with increasing strain, whereas at ε_xy = 5, the system returns to its initial (undeformed) state. When enforcing the FIC, on the other hand, the energy increases indefinitely with increasing strain, indicating that the extreme shear displacement is described properly.

3.3. Efficiency

The implementation of the MIC entails increased computational cost due to the calculation of the remainders in Equation (5); the computational time increases further when working with nonorthogonal boxes (see Equation (6)) due to the correction for shear displacement. Several algorithms have been proposed in the literature for enforcing the MIC [9,35]. Choosing an optimal algorithm depends on several factors such as the processor architecture and the compiler optimizations [35].

The FIC is generally cheaper since it requires adding a shift vector to the absolute separation vector. The shift vectors are computed only once during the formation of the bond and are stored in memory. The additional burden on system memory scales as O(n); it is negligible considering the modern computer architectures. It is worth mentioning that, in cases where the system is subjected to affine deformation, the shift vectors must be transformed according to Equation (12), thus increasing the cost of the operation slightly.

4. Conclusions

We have developed the fixed image convention (FIC) for determining the separation vectors of the intra- and intermolecular (effective) bonds in particle simulations. The accuracy of the FIC has been assessed in terms of conducting comparisons with the commonly adopted minimum image convention (MIC). The experiments for particle displacement and nonaffine deformation demonstrate that, when the projections of the separation vectors exceed half of the box projections, the minimum image separation vector swaps periodic images, and in this manner, the topology of the system is violated. The FIC, on the other hand, keeps the topology of the system intact, even at extreme deformations. Subjecting the system to affine deformation (i.e., applying the linear transformer to both the box vectors and particle coordinates) yields the same (stable) behavior, regardless of the convention. Given that in the FIC the separation vectors between the images are fixed, the system is able to return to its initial (undeformed) state regardless of how it is deformed. The computational cost for enforcing the FIC is generally lower than that of the MIC. The FIC is generic and applicable to several types of particle-based simulations and bonding schemes.