The simulation of polymer solutions often requires the development of methods that accurately include hydrodynamic interactions. Resolution on the atomistic scale is too computationally expensive to cover mesoscopic time and length scales on which the interesting polymer phenomena are observed. Therefore, coarse-graining methods have to be applied. In this work, the solvent is simulated using the well-established multi-particle collision dynamics scheme, and for the polymer, different coarse-graining methods are employed and compared against the monomer resolved Kremer–Grest model by their resulting diffusion coefficients. This research builds on previous work [Ruiz-Franco *et al.*, J. Chem. Phys. **151**, 074902 (2019)], in which star polymers and linear chains in a solvent were simulated and two different coarse-graining methods were developed, in order to increase computational efficiency. The present work extends this approach to ring polymers and seeks to refine one of the authors’ proposed model: the penetrable soft colloid model. It was found that both proposed models are not well suited to ring polymers; however, the introduction of a factor to the PSC model delivers satisfying results for the diffusion behavior by regulating the interaction intensity with the solvent.

## I. INTRODUCTION

In the last decades, computer simulations have become a powerful tool to understand physical principles of soft matter systems, such as colloidal suspensions and polymeric fluids. Conventional simulations of polymers at the atomistic level using full-atom molecular dynamics are limited in terms of length scales (typically nanometer) as well as time scales (typically microseconds). In order to overcome such limitations, several coarse-grained methods have been developed to capture the dynamics at mesoscopic scales while neglecting the behavior at the atomistic level. For example, for polymers, these strategies involve the substitution of interactions on the atomistic level by effective interaction potentials between coarse-grained segments of a polymer. A nowadays widely used approach approximates the polymer by a chain of spherical beads connected by springs.^{1} An important example is the well-known Kremer–Grest model^{2} that connects adjacent beads by nonlinear springs and bead-bead repulsion is modeled as a truncated Lennard-Jones potential, also known as the Weeks–Chandler–Anderson (WCA) potential.^{3} In these approaches, the solvent is often taken into account in an effective way. In the simplest case, it acts as a local heat bath that transmits uncorrelated stochastic forces onto the beads where the strength of these forces depends not only on temperature but on the bead size and the fluid viscosity.

However, such a simplified procedure neglects solvent-induced non-local effects that appear as hydrodynamic interactions between different beads. While equilibrium properties do not depend on capturing hydrodynamic interactions, for example, the dynamics of polymers under non-equilibrium conditions often involve fluid flow, such as in phoretic transport,^{4} polymers under shear,^{5} or in the presence of external forces.^{6} In order to keep the computational expenses of resolving hydrodynamic interactions reasonable, the hydrodynamic interactions can be treated by another level of coarse-graining, such as by Stokesian dynamics,^{7} the Lattice Boltzmann (LB) method,^{8} Dissipative Particle Dynamics (DPD),^{9} or Multiparticle collision dynamics (MPCD).^{10,11}

In particular, MPCD, which had first been introduced by Malevanets and Kapral more than 20 years ago,^{12} has since been widely used for the numerical simulation of complex fluids at the nano- and the micrometer scale. Examples include simulations of polymers,^{13–19} colloids,^{20–22} microswimmers,^{23–25} and biological cells.^{26} The MPCD method models the solvent by effective coarse-grained fluid particles that perform alternating streaming and collision steps, where momentum is conserved locally, and it thus approximates solutions of the Navier–Stokes equations and includes thermal fluctuations.^{11} Polymer dynamics is usually efficiently coupled to the solvent in the collision step.

Performing simulations of the bead-resolved polymer dynamics, with or without the presence of the explicit MPCD solvent, can become very time-consuming, in particular, when modeling large semidilute suspensions or polymer melts. Therefore, there has been great interest in finding even stronger levels of coarse-graining of polymers. In these approaches typically star polymers or dendrimers are treated as effective soft colloids,^{27} where the multitude of configurations is averaged out and individual polymers interact with each other only via an effective pair potential,^{28} instead of the monomer-resolved interactions.

While the aforementioned approaches using effective potentials are usually applied in equilibrium conditions, where fluid flow and hydrodynamic interactions can be neglected, recently Ruiz-Franco *et al.* investigated how linear and star polymers can be coarse grained as effective colloids and how they can be coupled to an MPCD fluid,^{29} allowing for incorporating hydrodynamic effects in the system. The authors introduced two different methods for describing linear polymers and star polymers as soft colloids in an MPCD fluid, which they validate against the fully monomer-resolved model (MRM) for polymers. The first method models the polymer as a rigidly connected distribution of monomers, which captures the real polymer density distribution obtained from the MRM. This method, also termed effective monomer model (EMM), had been proven to work well for linear polymers. The second method considers the polymer as a single penetrable soft colloid (PSC), where the solvent particles are allowed to penetrate the colloid with a certain probability, which is again related to the monomer density obtained from the MRM. This methods works particularly well for star polymers.^{29}

In this work, we extend the works from Ruiz-Franco *et al.*^{29} to ring polymers. Recently, there has been a great interest in understanding the dynamics and structural properties of ring polymer solutions,^{30–32} in particular, because of their importance in biological systems, such as chromatin organization.^{33} In previous studies, monomer-resolved models of ring polymers had been coupled to an MPCD fluid,^{34–36} and in particular the importance of hydrodynamic flows for ring polymers out of equilibrium had been demonstrated.^{36} We first introduce the methods, such as MPCD, the MRM, and the coarse-grained models developed in Ref. 29 (see Sec. II). We then introduce two modifications of the PSC model that we apply, in addition to the other model, to ring polymers. In Sec. III, we show our results on comparing the diffusive motion obtained of the effective descriptions of ring polymers for the different coarse-grained soft colloid models. We provide summary and conclusion in Sec. IV.

## II. METHODS

In this section, we introduce several methods to simulate a polymer in an MPCD solvent. Each method uses a different level of coarse-graining.

### A. Multi-particle collision dynamics (MPCD)

*n*point-like solvent particles each with mass

*m*are initialized. Their initial positions are random, and their velocities follow a Boltzmann distribution with temperature

*T*. The dynamics consist of two alternating steps: streaming and collision.

^{11,12,37}During the streaming step, the particles move ballistically for the duration of one time step

*h*with

*i*at time

*t*and $v\u20d7i(t)$ represents its velocity. For the collision step, the simulation box is divided into cubic collision cells of side length

*a*. Particles that share a collision cell interact with each other. The collision exchange of momentum is realized by rotating the particles’ relative velocity vectors $v\u20d7i\u2212v\u20d7cm$ with respect to the center-of-mass velocity of the collision cell $v\u20d7cm$, which is defined as

*n*

_{c}being the number of solvent particles inside the collision cell. The velocities after collision can then be calculated using a standard universal rotation matrix $D\u0302(R\u20d7,\alpha )$, which rotates the solvent particle relative velocities about a randomly chosen unit vector $R\u20d7$, i.e., the rotation axis, by a fixed angle

*α*,

*a*of the collision cell remains fixed during the entire simulation, their center points have to be shifted by a uniformly picked random vector with components ∼

*U*(−

*a*/2,

*a*/2) before each collision step. This is done as to not violate Galilean invariance. Galilean invariance and the local conservation of linear momentum are necessary to reproduce hydrodynamic interactions correctly.

^{11}

*ν*=

*ν*

_{kin}+

*ν*

_{coll}of the solvent can be tuned by the choice of parameters for collision cell side length

*a*, average number of solvent particles per collision cell $nc$, rotation angle

*α*, and time step size

*h*.

^{38}It can be separated into a kinetic part

*ν*

_{kin}, in which momentum is transported mainly through particle displacements, and a collisional part

*ν*

_{coll}, in which momentum is transported mainly through particle collisions

^{11}

^{,}

*k*

_{B}represents the Boltzmann constant and $nc$ refers to the mean number of particles per cell. The temperature

*T*is controlled using cell-level Maxwell–Boltzmann scaling.

^{37}The reference units for length, mass, energy, and time are

*a*,

*m*,

*k*

_{B}

*T*, and $\tau =ma2/kBT$, respectively, and they are generally set to

*a*= 1,

*m*= 1, and

*k*

_{B}

*T*= 1, following

*τ*= 1. We chose the remaining parameters as follows:

*h*= 0.1

*τ*and

*α*= 3

*π*/5. Furthermore, the cubic simulation box side length

*l*= 46

*a*and solvent particle number

*n*= 10

^{6}, which results in $nc=n/l3=10.3$. These are all within the usual range for an MPCD solvent,

^{14,39}and they result in a kinematic viscosity (with unit $\nu 0=a2kBT/m$) of

### B. Monomer-resolved model (MRM) for polymers

*N*monomers each with mass

*M*. The monomer–monomer interactions are simulated using the purely repulsive Lennard-Jones potential

*U*

_{LJ}(

*R*) for excluded-volume interactions as well as the finitely extensible nonlinear elastic (FENE) potential

*U*

_{FENE}(

*R*) for bonds.

^{40}This combination of potentials to simulate polymers is a well-established model developed by Grest and Kremer.

^{2}The purely repulsive Lennard-Jones potential

^{3}

^{,}

*θ*(

*x*) denotes the heavyside step function, and

*R*is the absolute distance between two monomers. We set the energy and length scale to

*ɛ*=

*k*

_{B}

*T*and

*σ*=

*a*, respectively. The FENE potential

*K*= 30

*ɛ*/

*σ*

^{2}for the spring constant and

*R*

_{max}= 1.5

*σ*for the maximum allowed distance between two neighboring monomers.

^{2}The linear combination of the two potentials, Eqs. (6) and (7), then leads to an equilibrium monomer–monomer distance of $\u2248\sigma $. In the streaming step the solute monomers do not interact with the fluid particles. A molecular dynamics of the polymer beads is performed by integrating Newton’s equations of motion using Eqs. (6) and (7) numerically with the velocity Verlet algorithm

^{41}at discrete time steps

*h*

_{MD}.

^{11,13}Analogous to the solvent, monomers are treated as point-like particles for the collision step with according mass

*M*and velocity $V\u20d7i$. The center-of-mass velocity of each cell with

*n*

_{c}solvent particles and

*N*

_{c}monomers is then calculated as

The parameters of the polymer and for the MD should be in accordance with the solvent and MPCD specific values. For this work, a time step of *h*_{MD} = 0.002*τ* was selected. MD steps and MPCD steps have to be alternated correctly such that the time is synchronized. Hence, the MD algorithm has to be performed *h*/*h*_{MD} = 50 times for each MPCD step. Also, the monomer mass should be roughly equal to the total solvent mass per cell;^{38} hence, we set *M* = 10m. This ensures that both the embedded object and the solvent influence each other equally. In this work, a ring polymer with *N* = 100 monomers in a cubic box using periodic boundary conditions was chosen for all simulations.

We initialize the ring polymer as a regular polygon in 2D. The results of the MD-MPCD simulations are averaged over 400 independent runs, where each run consists of 5 × 10^{5} MD steps of pre-equilibration without solvent, by which a typical equilibrium conformation for the polymer is produced, and 5 × 10^{4} MPCD steps both for equilibration and for the actual simulation. Note: Only the MRM needs time for equilibration. The subsequent coarse-grained models do not need equilibration due to their inherent coarse-grained nature.

^{36}

^{,}

*S*

_{i,α}(

*t*) =

*R*

_{i,α}(

*t*) −

*R*

_{cm,α}(

*t*) is the

*α*-coordinate of the relative position vector of monomer

*i*with respect to the polymer center-of-mass $R\u20d7cm$. The eigenvalues of the individual instantaneous gyration tensors are sorted ({Λ

_{1}(

*t*) > Λ

_{2}(

*t*) > Λ

_{3}(

*t*)}) and then averaged over time and ensemble. Of interest are especially the square-roots of the sorted, averaged eigenvalues ${\Lambda 1(t),\Lambda 2(t),\Lambda 3(t)}={\lambda 1,\lambda 2,\lambda 3}$, which represent the three principal radii of gyration of the polymer along the corresponding three principal axes.

^{42}We further use them to calculate the radius of gyration

*R*

_{G}via

*N*= 100 in an MPCD-solvent are

### C. Effective monomer model (EMM)

The effective monomer model (EMM) is a coarse-grained model for polymers, which was introduced by Ruiz-Franco *et al.*^{29} The principle behind this model is to keep the polymer as monomer resolved for the interaction with the solvent, but remove internal interactions between the monomers. By integrating out the fluctuation degrees of freedom, the EMM describes a rigid but soft object and can be categorized between a polymer and a colloid. Therefore, we will refer to the coarse-grained polymers also using the term “colloid.”

*a*=

*σ*. From simulations using the monomer resolved model (Sec. II B), we obtain the monomer distribution relative to the center-of-mass as a radial density function

*ρ*

_{mon}(

*r*), with

*r*being the absolute distance to the center-of-mass. The radial monomer density fulfills the normalization condition

^{43}

*N*

_{mon}(

*r*) is the cumulative particle number from the center-of-mass at

*r*′ = 0 up to

*r*′ =

*r*. With that, we also define a cutoff radius

*R*

_{coll}, on which we impose the condition

*ρ*

_{mon}(

*R*

_{coll})

*σ*

^{3}< 10

^{−3}, to limit the polymer’s extension. Then, accordingly,

*N*

_{mon}(

*R*

_{coll}) =

*N*is the total monomer number.

To continue with initializing the effective polymer, a uniform random number $R\u223cU(0,1)$ is drawn for each of the cells and compared to the value of the volume density *ρ*_{mon}(*r*_{cell})*σ*^{3} at the distance *r*_{cell} of the cell’s center to the simulation box center. If $R\u2264\rho mon(rcell)\sigma 3$, we insert a monomer at the center of the respective cell. At the end, we have a number of monomers *N*_{eff} equal to the number of successful insertions with expectation value $Neff=N$. Figure 1 shows an example configuration generated through this process.

Instead of undergoing molecular dynamics following the intermolecular potentials as described for the MRM (Sec. II B), the colloid moves as a rigid body and keeps its monomer configuration fixed over the course of the simulation. Consequently, the EMM neglects the angular momentum of the polymer. The colloid displacement, and with it all effective monomers, is ballistic in the same way as for the solvent particles following Eq. (1) and uses its center-of-mass velocity $V\u20d7cm$ for the ballistic movement.

Ruiz-Franco *et al.*^{29} suggested that after calculating the new center-of-mass velocity, the effective monomer velocities should follow such that each monomer has the same velocity, which is the center-of-mass velocity $(V\u20d7j(t+h)=V\u20d7cm(t+h))$, before the next collision step. If we observe the colloid temperature, however, we find that this balancing makes the colloid cooler compared to the thermostated solvent. This is because the momentum exchange occurs through the effective monomers with mass *M*_{mon} as opposed to the entire colloid with mass *M* = *N*_{eff}*M*_{mon}. Through matching the velocities of the monomers with the center-of-mass velocity, however, they are no longer Maxwell–Boltzmann distributed and their temperature is on average comparably lower. Even after a collision with the thermostated solvent particles, the temperature has not adjusted sufficiently fast before the next averaging of the velocities occurs. For this reason, after a collision, the individual monomer velocities are kept for the next collision in this work. The average velocity is only used for the colloid displacement. Figure 2 shows the difference in temperature for both methods.

The monomer configuration is not redrawn during the simulation process, but colloid dynamics are evaluated by averaging over a sample size of 300 fixed rigid body configurations instead, with each run consisting of 5 × 10^{4} MPCD steps.

### D. Penetrable soft colloid model (PSC)

The penetrable soft colloid model (PSC) was also introduced by Ruiz-Franco *et al.*^{29} In this coarse-grained model, the monomer resolution is removed, and instead the polymer is depicted as one single spherical colloid that is to a certain extent penetrable for solvent molecules. Movement of the solvent particles inside the colloid leads to interaction between the two, where the colloid itself acts as one rigid particle with no internal interaction.

*I*of the colloid in addition to its total mass

*M*=

*NM*

_{mon}. For a radially symmetric object with axis of rotation through the center-of-mass, the moment of inertia tensor can be simplified to a scalar. In this case, it is calculated as

*ρ*

_{mon}(

*r*) and cutoff radius

*R*

_{coll}as defined in Sec. II C.

*r*denotes the distance of the solvent particle relative to the colloid center $R\u20d7$. Trial moves for which both the initial and the final position lie outside the colloid are always accepted. If either one of the solvent particle’s initial relative position $r1=|r\u20d71|=|r\u20d7i(t)\u2212R\u20d7(t+h)|$ or trial end relative position $r2=|r\u20d72|=|r\u20d7i(t+h)\u2212R\u20d7(t+h)|$ are inside the colloid (

*r*

_{1}≤

*R*

_{coll}or

*r*

_{2}≤

*R*

_{coll}), the acceptance probability is calculated as follows:

*ρ*

_{sol}(

*r*), which is the local volume fraction that is populated by the solvent after the monomer excluded volume is subtracted. Since monomers are assumed to be spherical the solvent density is calculated via

^{21}

*ν*

_{n}that is perpendicular to the colloid surface and the velocity component

*ν*

_{t}that is tangential to it. Linear and angular momentum conservation remain upheld. The new absolute velocity of the solvent particle then calculates as

*ρ*

_{sol}(

*r*) can be achieved in the end. The tangential unit vector $e\u0302t$ is chosen randomly under the condition that it is orthogonal to the normal unit vector: $e\u0302t\u22a5e\u0302n$.

*i*, the velocity and angular velocity of the PSC also have to be updated,

^{4}MPCD steps.

It turns out that neither of the two models developed by Ruiz-Franco *et al.*^{29} deliver satisfying diffusion coefficients compared to the MRM for a ring polymer (see Table I in Sec. III B for results). While the EMM works well for linear chains and the PSC is able to capture star polymers, ring polymers require the development of a different model. Therefore, we introduce two new coarse-graining models to test against the MRM for ring polymers: the PSC model with ellipsoidal shape (EPSC) and the PSC model with increased interaction factor (PSC-f).

Model . | n_{sim}
. | D_{l}/D_{0}
. | R_{H}/a
. | D/D_{0}
. |
---|---|---|---|---|

MRM | 400 | 0.001 28 ± 0.000 06 | 4.29 | 0.001 73 |

EMM | 300 | 0.001 15 ± 0.000 06 | 4.64 | 0.001 60 |

PSC | 400 | 0.007 02 ± 0.000 32 | 0.99 | 0.007 48 |

EPSC | 300 | 0.006 30 ± 0.000 30 | 1.10 | 0.006 56 |

PSC-f (th.)a | 400 | 0.001 91 ± 0.000 09 | 3.14 | 0.002 36 |

PSC-f (emp.)b | 400 | 0.001 32 ± 0.000 06 | 4.19 | 0.001 77 |

Model . | n_{sim}
. | D_{l}/D_{0}
. | R_{H}/a
. | D/D_{0}
. |
---|---|---|---|---|

MRM | 400 | 0.001 28 ± 0.000 06 | 4.29 | 0.001 73 |

EMM | 300 | 0.001 15 ± 0.000 06 | 4.64 | 0.001 60 |

PSC | 400 | 0.007 02 ± 0.000 32 | 0.99 | 0.007 48 |

EPSC | 300 | 0.006 30 ± 0.000 30 | 1.10 | 0.006 56 |

PSC-f (th.)a | 400 | 0.001 91 ± 0.000 09 | 3.14 | 0.002 36 |

PSC-f (emp.)b | 400 | 0.001 32 ± 0.000 06 | 4.19 | 0.001 77 |

^{a}

PSC-f whose factor was calculated from theoretical considerations.

^{b}

PSC-f whose factor was empirically approximated.

### E. PSC model with ellipsoidal shape (EPSC)

Within the standard PSC model we have approximated the ring polymer to a spherical shape. While this may be reasonable for star polymers or polymers in a collapsed state due to bad solvent conditions, for ring polymers or linear chain polymers in good solvent conditions, an ellipsoidal shape poses a much more appropriate choice. Considering that the gyration tensor obtained from the MRM suggests an ellipsoid with an aspect ratio of *λ*_{1}/*λ*_{3} ≈ 3 for our ring polymer, this can hardly be considered an approximately spherical shape. The EPSC is a refined version of the PSC discussed in Sec. II D, which does not assume spherical shapes but allows the polymer to be approximated by an ellipsoid with any set of three length parameters {*a*, *b*, *c*} that represent its extension along the three principal axes. It is known that the diffusion coefficient of ellipsoids is smaller than that of spheres with the same volume.^{44} This indicates that an ellipsoidal colloid could close the gap between discrepancies in the measured diffusion coefficients between the MRM and the standard PSC model.

#### 1. Parameterizing the EPSC

Because of the anisotropy of the colloid, its orientation is tracked over time. Hence, we define a rotation matrix $Q\u0302(t)$ that we use to track the orientation of the colloid. $Q\u0302(t)$ is updated every step using its angular velocity $\Omega \u20d7(t)$.

*t*in the frame of reference of the colloid (denoted by primed variables), we have $R\u20d7\u2032(t)=0$, $V\u20d7\u2032(t)=0$, $\Omega \u20d7\u2032(t)=0$, and $Q\u0302\u2032(t)=1$. In the frame of reference of the colloid, we parameterize the EPSC with the parameters {

*a*,

*b*,

*c*} and resulting volume

*V*

_{E},

*r*′,

*θ*′,

*ϕ*′) instead of Cartesian coordinates (

*x*,

*y*,

*z*),

*M*=

*NM*

_{mon}. For the monomer density function, we use the radially symmetric density function evaluated previously, but transform it according to the stretched spherical coordinates. That means, we transform

*r*→

*r*′ and then normalize again the monomer density

*ρ*

_{mon}(

*r*′),

*V*′ =

*abc r*′

^{2}sin

*θ*d

*r*′d

*θ*d

*ϕ*and

*N*=

*N*

_{mon}(

*R*

_{coll}). We can hereby calculate all elements to calculate the moment of inertia tensor for the ellipsoid $I\u0302$ and can write it in diagonal form

^{45}

*λ*

_{1},

*λ*

_{2},

*λ*

_{3}}, obtained from the MRM simulation, are equal to the semi-axis lengths of the ellipsoid of inertia.

^{46}In order to match a physical ellipsoid, they can simply be scaled with a factor. Using Eq. (31) then allows us to determine {

*a*,

*b*,

*c*},

#### 2. EPSC movement and coupling to the solvent

In principal, the EPSC is simulated in the same way as the spherical PSC. Each time step consists of three steps: EPSC ballistic movement, solvent particles ballistic movement with colloid interaction, and the multi-particle collision between the solvent particles.

*α*(

*t*) and the rotation axis $R\u20d7(t)$, both of which can be calculated from the angular velocity $\Omega \u20d7(t)$ following

*x*

_{s}, $ys$,

*z*

_{s}) define the scattering point coordinates in the colloid frame of reference. The new colloid angular velocity is calculated using the inverse moment of inertia tensor $I\u0302\u22121$ such that

The third step, which is the multi-particle collision among solvent particles, is performed the same way as if there was no solute present. The EPSC dynamics is evaluated by averaging over a sample size of 300 independent runs of 5 × 10^{4} MPCD steps each.

### F. PSC model with increased interaction factor (PSC-f)

If the theoretical diffusion coefficient is known, or it can be acquired by a simulation following the MRM algorithm (as described in Sec. II B), it is possible to adjust the acceptance probability *p*_{acc}(*r*_{1} → *r*_{2}) of the PSC algorithm in such a way that the desired diffusion coefficient is enforced. The algorithm of the PSC-f is based on the standard PSC algorithm. The only difference is in the calculation of the Monte Carlo acceptance probability *p*_{acc}(*r*_{1} → *r*_{2}) and consequences arising from that.

^{47}

*ρ*(

*r*

_{i}) describes the probability density to find state

*r*

_{i}. If we introduce a factor 0 <

*f*≤ 1, detailed balance would be upheld for

*f*≤ 1 the ensemble will converge to its equilibrium state. The factor

*f*only lowers the rate of this convergence. The PSC algorithm does not use the usual Metropolis algorithm, but the acceptance probability is conditional. Tailored to the PSC algorithm, it results in

*f*. One way to realize this is by distinguishing between four cases (as opposed to the two cases in the original PSC model) in order to determine the orientation of the normal unit vector $e\u0302n$:

*r*_{2}≤*r*_{1}and $\rho sol(r2)\rho sol(r1)\u22651\u21d2e\u0302n$ points outwards,*r*_{2}≤*r*_{1}and $\rho sol(r2)\rho sol(r1)<1\u21d2e\u0302n$ points outwards,*r*_{1}<*r*_{2}and $\rho sol(r1)\rho sol(r2)\u22651\u21d2e\u0302n$ points inwards,*r*_{1}<*r*_{2}and $\rho sol(r1)\rho sol(r2)<1\u21d2e\u0302n$ points inwards, with probability: $(1\u2212f)1\u2212f\rho sol(r1)\rho sol(r2)\u22121$.

*f*from theoretical considerations, which is shown in the Appendix. However, it turns out that the resulting diffusion constant does not resemble the desired diffusion constant sufficiently well (see Table I). Therefore, the factor has to be empirically estimated, keeping in mind that decreasing the factor increases the interaction rate and therefore decreases the diffusion coefficient. For the PSC-f, we tuned

*f*to fit the desired outcome regarding the diffusion coefficient and obtained a value of

^{4}MPCD steps each.

## III. RESULTS FOR RING POLYMERS

In this section, we present results for the monomer and solvent density profiles and calculate the diffusion coefficient of the ring polymer within the different models.

### A. Monomer and solvent density profiles

*ρ*

_{mon}(

*r*) is assessed through simulation using the MRM. In order to obtain a smooth analytical function, the discrete data are approximated using an empirically created fit function. Figure 5 shows the simulation data for a ring polymer organized into a histogram and their empirical fit, namely,

*A*,

*b*,

*d*,

*s*. The fit parameters for the ring polymer with

*N*= 100 result in

*A*= 2.246 × 10

^{−4}

*σ*

^{−8},

*b*= 66.16

*σ*

^{3},

*d*= 365.2

*σ*

^{5}, and

*s*= 2.128

*σ*.

From the monomer density function and using Eq. (16), we obtain the theoretical solvent density profile. Considering time for equilibration, the actual radial distribution of the solvent particles is compared with the theoretical distribution for all soft colloid models. The solvent density function in the original PSC and the PSC-f can be seen in Fig. 6. By introducing the Monte Carlo acceptance probability for solvent movement in the PSC models, the colloid creates an excluded volume, which enables a more realistic solvent density distribution. This is a feature that is not inherent to the MRM or the EMM because, in the standard MPCD model, there is no excluded volume for solvent particles.^{1} As shown in Fig. 6, the desired solvent density distribution is achieved with the PSC. Introducing an increased interaction factor in the PSC-f keeps this feature if the acceptance rules are properly executed.

### B. Diffusion constant

*D*can be estimated from the mean squared-displacement of the center-of-mass over time (MSD = $\Delta r\u20d7cm2(t)$) according to

^{48}Consequently, for a box with periodic boundaries, the diffusion coefficient is dependent on the system size. For this reason, the diffusion coefficient obtained by the simulation is only the finite-size diffusion coefficient

*D*

_{l}of system size

*l*. The infinite-system-size diffusion coefficient

*D*can be obtained from the finite-size diffusion coefficient

*D*

_{l}by solving for the hydrodynamic radius

*R*

_{H}in

^{49}

*l*and the dynamic viscosity of the solvent

*η*

_{sol}=

*ϱν*. The infinite-system-size diffusion coefficient can subsequently be calculated using the Stokes–Einstein relation

*l*= 46

*a*and a solvent dynamic viscosity

*η*

_{sol}/

*η*

_{0}= 7.143 $(\eta 0=kBTm/a2)$, the results regarding polymer diffusion for all employed models are summarized in Table I.

When we calculate the ratio between radius of gyration and hydrodynamic radius, we can further validate the accuracy of our methods. Literature suggests a ratio of *R*_{G}/*R*_{H} = 1.2533 for Gaussian ring polymers.^{50} Using the MRM, we find a ratio of *R*_{G}/*R*_{H} = 1.25, which is in excellent agreement with the theoretical value.

## IV. SUMMARY AND CONCLUSIONS

Neither the EMM nor the PSC model from Ruiz-Franco *et al.*^{29} are able to adequately reproduce the long-time diffusion constant of the ring polymer. Adapting the PSC model by enabling not only spherically shaped colloids but general ellipsoids (EPSC) does not significantly improve the diffusion behavior for ring polymers either, despite the fact that they have rather large aspect ratios *a*/*c* ≈ 3. On the other hand, introducing a factor to the PSC’s Monte Carlo acceptance probability (PSC-f), thus freely increasing interaction, provides the desired results and slows down the diffusion coefficient such that it matches the MRM. An advantage of this method is that it offers a lot of freedom since the factor can be chosen anywhere between 0 < *f* ≤ 1. As a result, the diffusion coefficient can be decreased arbitrarily within range while maintaining the desired solvent density profile. A disadvantage of this method is that the target diffusion coefficient has to be known, and the factor has to be approximated empirically.

The reason why the PSC model as proposed by Ruiz-Franco *et al.*^{29} does not work well for ring polymers lies most likely in the fact that it relies on the solvent density gradient to control the probability of interaction. Compared to that, interaction in the monomer-resolved model is better resembling the dependence on absolute solvent density. The reason why this model works for star polymers but not for linear chains and rings is probably that star polymers have a steeper monomer density function, and hence a higher solvent density gradient overall. As a result they possess higher interaction probabilities and consequently a slower diffusion.

For future work, it would be interesting to investigate the different coarse-grained models using varying degrees of polymerization *N*. Furthermore, testing the several coarse-grained models on many-polymer-systems would be the natural next step in this analysis.

## ACKNOWLEDGMENTS

We would like to acknowledge the Vienna Scientific Cluster (VSC) for providing us with the computational resources that were necessary to perform our simulations.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Lisa Sappl**: Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). **Christos N. Likos**: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Project administration (lead); Resources (lead); Supervision (equal); Validation (supporting); Writing – review & editing (lead). **Andreas Zöttl**: Formal analysis (equal); Investigation (equal); Methodology (supporting); Supervision (equal); Writing – original draft (supporting); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: THEORETICAL DERIVATION OF THE MONTE CARLO FACTOR

*D*

_{l}is inversely proportional to the solvent’s viscosity [Eq. (43)]. For the solvents that are of interest to these kind of simulations, usually the collisional viscosity dominates over the kinetic viscosity

*ν*

_{kin}≪

*ν*

_{coll}, and therefore, the kinetic viscosity may be neglected and the viscosity is approximately inversely proportional to the time step

*h*. This leads to an inverse relationship between the diffusion coefficient and the expected Monte Carlo rejection probability

*p*

_{rej}(

*r*

_{1}→

*r*

_{2}) = 1 −

*p*

_{acc}(

*r*

_{1}→

*r*

_{2}),

*D*

_{PSC}, then the theoretically desired value

*D*

_{theo}can be obtained using the MRM. The two differ by a factor

*f*

_{D}, which can be expressed in terms of the expected Monte Carlo rejection probabilities as

*f*used in the Monte Carlo scheme directly changes the acceptance probability

*f*(

*r*

_{1},

*r*

_{2}) for

*r*

_{1}and

*r*

_{2}in

*r*

_{1}are distributed following

*ρ*

_{sol}(

*r*) as provided by the Monte Carlo scheme. The trial end position $r2=|r\u20d71+hv\u20d7|$ is dependent on the initial position and its probability density is also a function of

*r*

_{1}. The variable that connects

*r*

_{1}and

*r*

_{2}, which is $\Delta r=|hv\u20d7|$, on the other hand, is independent of

*r*

_{1}and follows only the Maxwell–Boltzmann distribution. The joint probability density

*f*(

*r*

_{1}, Δ

*r*) =

*f*

_{1}(

*r*

_{1}) ·

*f*

_{Δ}(Δ

*r*) is therefore separable and simply a product of the individual probability densities

*r*, they are oriented such that

*θ*

_{Δ}denotes the angle between the vectors $r\u20d71$ and $\Delta r\u20d7$ (see Fig. 8). The rejection probability

*p*

_{rej}(

*r*

_{1}→

*r*

_{2}) uses the quantities

*r*

_{1}and

*r*

_{2}. From geometrical considerations, it is possible to derive the transformation (Δ

*r*,

*θ*

_{Δ}, $\varphi \Delta $) → (

*r*

_{2},

*θ*

_{2},

*ϕ*

_{2}) with

*J*,

*θ*

_{2}denotes the angle between the vectors $r\u20d71$ and $r\u20d72$. The Jacobi determinant

*J*helps to transform the volume element. The joint probability density for

*r*

_{1}and

*r*

_{2}can now be calculated by transforming and subsequent integration over the other variables,

*ρ*

_{sol}(

*r*), the expression in Eq. (A5) can be calculated numerically, although it is wise to rearrange Eq. (A11) for that by joining the hyperbolic sine and exponential function, as otherwise the hyperbolic sine increases too quickly. The factor we obtained using this theoretical approach was

*f*= 0.997 38).

## REFERENCES

*Computer Simulation of Liquids*

*The Lattice Boltzmann Method*

*Advanced Computer Simulation Approaches for Soft Matter Sciences III*

*Hierarchical Methods for Dynamics in Complex Molecular Systems*

*Elements of the Random Walk: An Introduction for Advanced Students and Researchers*

*Theoretische Physik 1*

*Introductory Lectures on Knot Theory*

*Advances in Microfluidics*