An Anisotropic Effective Model for the Simulation of Semiflexible Ring Polymers.

We derive and introduce anisotropic effective pair potentials to coarse-grain solutions of semiflexible ring polymers of various lengths. The system has been recently investigated by means of full monomer-resolved computer simulations, revealing a host of unusual features and structure formation, which, however, cannot be captured by a rotationally averaged effective pair potential between the rings' centers of mass [Bernabei M.; Soft Matter2013, 9, 1287]. Our new coarse-graining strategy is to picture each ring as a soft, penetrable disk. We demonstrate that for the short- and intermediate-length rings the new model is quite capable of capturing the physics in a quantitative fashion, whereas for the largest rings, which resemble flexible ones, it fails at high densities. Our work opens the way for the physical justification of general, anisotropic penetrable interaction potentials.

Introduction

By the simple process of joining the ends of a linear polymer chain, one obtains a ring polymer (RP).[1] While the architecture of ring polymers is very simple, they differ in many interesting ways from their linear counterparts and are the subject of active research in physics, biology, chemistry, and even pure mathematics. One interesting consequence is that for a dynamics that disallows strand crossing, there are different classes of configurations of a RP, which can never transform into each other. These are referred to as the topology classes or knot types of a RP. Knot theory is a fascinating and active branch of mathematics with many open problems concerning the enumeration and classification of knots.[2] A RP that has the topology of a circle is called an unknotted RP.

Unlike RPs, linear polymer chains can strictly speaking never be knotted, as every configuration of a linear polymer chain can be continuously transformed to a straight line, without the need for strand crossings. While this is a fundamental difference between linear polymer chains and RPs, it is nevertheless possible to extend the concept of knots to physical knots on linear polymer chains.[3] There are many works dealing with the properties of these physical knots on linear polymer chains,[4−11] in particular due to their relevance in biophysics, where they are for instance found in DNA[12,13] and can have significant effects on key processes.[14−16]

The topological constraint of a ring polymer has important consequences for its physical behavior. It took the work of many authors[17−24] to establish that the diameter of gyration Dg of an isolated ideal RP with fixed topology scales as ⟨Dg2⟩ ∼ N2ν, where ν ≈ 0.588 is the Flory exponent, which also describes the scaling behavior of the radius of gyration of self-avoiding linear chains.[25] This is true for all knot types, even for an ideal unknotted RP (i.e., without monomer excluded volume interactions, just keeping the topological constraints). An ideal linear polymer chain, on the other hand, remains more compact and exhibits the scaling law of a nonself-avoiding random walk ⟨Dg2⟩ ∼ N. Another important difference lies in the effective potential between the centers of mass of RPs. While the effective potential vanishes between infinitely thin linear polymer chains and also between an infinitely thin linear polymer chain and a RP, there remains a nonzero repulsive contribution between cyclic polymers with fixed topology.[26,27] Here one usually speaks of a topological potential. Furthermore, it was shown that the effective potential between two moderately sized RPs increases with the knot complexity of the rings.[27]

Also for concentrated systems, the topology of polymer chains plays an important role. The scaling of linear polymer chains in the melt is the one of a nonself-avoiding random walk ⟨Dg2⟩ ∼ N. Simulations of dense systems of RPs,[28,29] on the other hand, showed that while short chains also exhibit a Gaussian scaling behavior ⟨Dg2⟩ ∼ N, long chains are compact and thus scale as ⟨Dg2⟩ ∼ N2/3. In between there is a broad crossover region, where a ⟨Dg2⟩ ∼ N4/5 scaling provides a good description of the data. For the dynamics, it is expected that concatenations of ring polymers can have a significant effect, as they are permanent, in contrast to the entanglement of linear polymer chains. However, this implies that those concatenations are there in the first place, i.e., from the very synthesis of the sample on. Even in the absence of concatenations, there are important differences in the dynamics of RPs in the melt with respect to their linear counterparts. For instance, recent experiments[30] and simulations[31] revealed a power-law stress relaxation instead of the rubbery plateau found for linear chains.

For the large intermediate density domain between dilute solutions and melts, there are relatively few theoretical results despite the practical relevance of this regime for instance in the field of biophysics, where the topological interactions between chromatin loops plays a crucial role in the creation of chromosome territories.[26,32−34] A fruitful and modern approach for the economic description and simulation of macromolecules in this regime is the method of coarse-graining. The idea behind this method is to bridge the time and length scales in the system by describing the macromolecules via an effective model with a reduced set of suitably chosen effective degrees of freedom (dof). The microscopic information on the monomer-resolved model is underlying the effective model, as it determines the form of the effective potential, which describes the interaction between the macromolecules. The advantage of this method is not only that every time step in a simulation requires less computational effort due to the simplified representation but also that one can often choose a much larger time step in a simulation of the coarse grained model, as the dof that remain in the coarse-grained model change much slower in time than their counterparts in the monomer-resolved model.[35,36]

The method of coarse-graining is well-established and has for instance found successful application for polymer chains,[37−39] star polymers,[40−43] star-shaped polyelectrolytes,[44,45] dendrimers,[46−48] and block copolymers.[49−51] The identification of the relevant degrees of freedom is an essential part in the design of an effective model. One often uses isotropic effective models, where the macromolecules consisting of many individual monomers are reduced to their center of mass. For the semiflexible ring polymers such a model has already been investigated in ref (52). Clustering was observed in monomer-resolved simulations of semiflexible ring polymers as well as in the corresponding isotropic effective model. However, it was also shown that the monomer-resolved system shows anisotropic features that cannot be accounted for in the isotropic effective model. Also, the correlation functions stemming from the isotropic effective model are markedly different from the microscopically derived ones. Anisotropy is particularly strong for rings with high bending stiffness or few monomers, as they have a strong tendency to orient with respect to other rings in their proximity. This motivates us to introduce an anisotropic effective model for the description of semiflexible ring polymers in this article. In this model, we will define the effective particles as soft disk-like molecules which are described not only by their center of mass but also by the direction in which their faces are oriented. An anisotropic effective model was already used successfully for the description of hard disk-like macromolecules,[53] but to the best of our knowledge this approach was up to now never applied to penetrable macromolecules, where the centers of mass of the macroparticles can coincide. Penetrable particles are particularly interesting as they allow for clustering, which often leads to a rich phase behavior. For instance, point-particles interacting with a certain class of ultrasoft potentials form so-called cluster crystals.[54−56] Unlike in an ordinary crystal multiple soft particles can sit on top of each other at the same lattice site in a cluster crystal. Another peculiar feature of this state of matter is that by compressing it one only changes the occupation number of particles per lattice site, while the lattice constant remains invariant. Monomer-resolved simulations of semiflexible ring polymers, on the other hand, show the formation of the cluster glass phase,[57,58] which is an arrested state that also contains some of the features found in the cluster crystal phase. In both cases the overall structure of the system is frozen, while individual particles can hop between the lattice sites of the cluster crystal or the stacks found in the cluster glass. Elongated dendrimers, which unlike hard rod-like particles exhibit local antinematic order,[200] are another interesting example for a system of penetrable particles, which behaves distinctively different to its solid counterpart. By creating an anisotropic effective model for the semiflexible ring polymers, we aim at a model that is still computationally cheap and that improves the description obtained by the isotropic model, especially for the case of high densities. In addition, the analysis of the interactions in the anisotropic effective potential allows us to get a better understanding for the interaction between the anisotropic, penetrable nanoparticles in systems of semiflexible RPs.

The remainder of this article is structured as follows. In section we first present the Hamiltonian of the monomer-resolved model which we use for the description of semiflexible RPs and then introduce an anisotropic effective model for such a system. In section A of the Appendix we give more details about the derivation of the effective interactions for such a model. We carried out molecular dynamics simulations of the monomer-resolved model and Monte Carlo simulations of the effective models; details concerning these simulations are given in section . In section we present the anisotropic effective potential and discuss its features. Results of Monte Carlo simulations with this potential, which show that the inclusion of anisotropy in the effective model can significantly improve the agreement with the monomer-resolved model, are presented in section , whereas in section we briefly discuss the effects of truncation of the expansion of the potential on the quality of the results. Conclusions are given in section . In the Appendix, we explain the expansion of the anisotropic pair-correlation function of a system of two RPs, which contains all the information for calculating the effective potential, as a sum of suitably chosen basis functions.

Anisotropic Effective Model

Monomer-Resolved Model for Semiflexible Ring Polymers

The derivation of the anisotropic effective model is based on a microscopic model of semiflexible ring polymers, each consisting of N monomers. They are described with the bead–spring model by Kremer and Grest[59] and an additional rigidity term. Thus, any two monomers interact via the truncated and shifted Lennard-Jones potentialThis potential is purely repulsive, accounting then for monomer excluded volume interactions. Bonded monomers also interact through a finitely extensible nonlinear elastic potential (FENE)Rigidity is introduced via the bending potentialwhere θ is the angle between two consecutive bond vectors. [Note that this is not the Kratky–Porod model (linear in the cosine). We expect the same qualitative results for Kratky–Porod rings with the same N’s and persistence lengths as the model simulated here.] The potential Vbend vanishes for θ = 0, when the polymer chain does not bend at the respective angle. We choose ϵ = kBT, k = 30kBT/σ2, R0 = 1.5σ, and κ = 30kBT, where kB is the Boltzmann constant and T the temperature. These are precisely the parameters employed in the simulation study of ref (52). The corresponding dynamics does not allow for chain crossings and thus topology is preserved.

In ref (58) the characteristic ratio[60] of this polymer model was estimated by carrying out simulations of isolated linear chains. Excluded volume interactions were switched off except for mutually connected monomers in order to obtain long-range Gaussian statistics. C∞ was obtained by analyzing the long-s limit of the ratio ⟨R2(s)⟩/s⟨b2⟩, where R(s) is the distance between two monomers i, j with s = |i – j|, and b is the bond length (⟨b2⟩ = 0.94). The authors reported a value of C∞ ∼ 15, which is typical for stiff polymers.[60] We can give an estimate for the persistence length of the model by mapping it to the freely rotating chain model using the relation cos θ = (C∞ – 1)/(C∞ + 1), where θ is the bending angle of the freely rotating chain model.[60] The persistence length is then obtained as sb = −b/ln(cos θ) ∼ 7.3. We carried out simulations of ring polymers with N = 20, 50, and 100 monomers, which have the contour to persistence length ratio N/s ∼ 2.7, 6.7, and 13.3, respectively.

Anisotropic Effective Model

In earlier work,[52] Bernabei et al. carried out extensive monomer-resolved simulations of stiff ring polymers to obtain the structure of concentrated solutions of the same. In an attempt to coarse-grain the system in the simplest possible way, they also derived and employed an isotropic effective potential for their effective description, reducing thereby stiff RP’s into point-like effective particles, namely their centers of mass. At this level of approximation, the effective particles possess no other internal (spin-like) degrees of freedom, and thus the effective interaction is isotropic. The effective potential between these macroparticles was defined by calculating the pair correlation function giso(r) in an infinitely dilute system and usingto define the effective interaction potential between these point particles, where β = 1/(kBT). In the infinitely dilute case, the distribution of the centers of mass of the ring polymers in equilibrium is identical to the distribution of the point particles in the effective model. At higher densities, however, it turns out that multiparticle terms in the effective potential are necessary to obtain the correct equilibrium distribution of the centers of mass in the effective model. One of the reasons due to which multiparticle interactions become important is that two ring polymers that are sufficiently close and stiff will prefer to align parallel to each other. A third ring polymer interacting with those two will not see them as two independent rings but as a system of two rings that are correlated. When using the potential Veff(r) calculated in eq , one assumes that the free energy penalty of one ring with respect to a second is independent of the presence of another polymer in the vicinity of the second. If the ring polymers preferentially align parallel, this assumption is clearly violated and one has to correct the effective potential by introducing multiparticle terms. Bernabei et al. showed that already at moderate densities one can encounter strong correlations of the orientations of the semiflexible ring polymers, in particular if the chains contain only few monomers (e.g., N = 20).[52] Therefore, a more accurate coarse-graining which takes into account the ring anisotropy is called for.

The easiest way to incorporate the correlations of the orientations of rings is to introduce them via additional degrees of freedom (dof) in the effective description. This is precisely what we do in this article. For this purpose, we need to first come up with a suitable definition for the orientation of a RP. To this end, we make use of the gyration tensorwhere rα( (α = x, y, z, Cartesian components) denotes the position of the ith monomer with respect to the center of mass of the ring to which this monomer belongs. The eigenvectors of Sαβ are the principal axes of an ellipsoid that approximates the shape of the macromolecule: If a RP is flat, which means that all its monomers lie in one plane, the ellipsoid has one zero eigenvalue with a corresponding eigenvector that is perpendicular to that plane. Also in the more general case, where the monomers do not all lie in the same plane, we define the normalized eigenvector corresponding to the smallest eigenvalue of Sαβ as the direction vector d of the RP. Note that d and −d are equivalent for reasons of symmetry. The ring polymers in the anisotropic effective model we propose are described via the position vectors of their centers of mass, R(, and their direction vector, d(. Henceforth, we describe the stiff rings as soft circular disks and ignore differences in the other two eigenvectors of the gyration tensor Sαβ. This choice is motivated by the limit of infinite bending stiffness, where the rings assume flat and precisely circular conformations. Our model therefore amounts to the minimal anisotropic extension of the spherically symmetric effective interaction between the centers of mass. We emphasize, however, that there is no a priori guarantee that this will be an improvement over the isotropic model at finite densities and in particular at high concentrations: this depends on the degree in which the RP’s at high concentration maintain their anisotropic shape and properties encoded in the high-dilution limit in which the anisotropic pair potential is derived. Accordingly, the introduction of such a potential is not a straightforward part of a systematic strategy of introducing more and more detail into the effective description of the system.

In order to determine the anisotropic effective potential Veff, we carried out monomer-resolved simulations of two ring polymers inside a large simulation box. The effective pair potential is then defined such that it exactly reproduces the correlation functions of the effective degrees of freedom in this infinitely dilute, monomer-resolved simulation. In the effective model, two ring polymers are described by a total of 10 degrees of freedom—three for each center of mass and two for each direction vector of each ring polymer. However, due to translation, rotation, and mirror symmetry, the distinct configurations (those that cannot be related by symmetry transformations) of a system with two effective particles are reduced and can be specified by four parameters only.

A convenient choice for these variables is illustrated in Figure and reads as follows:where r ≡ R(2) – R(1) is the connection vector between the centers of mass of the two rings, r̂ ≡ r/r the unit vector in the direction of r, and d⊥( the component of the director d( perpendicular to r; 0 ≤ φ ≤ π denotes the angle between vectors d⊥(1) and d⊥(2). By selecting the appropriate sign of d(, we can always choose cos θ to lie in the interval [0, 1].

Figure 1

Illustration of the effective variables r = |r|, θ, and φ with which relative configurations of ring polymers are described.

Let us define the ideal case as the system where the effective particles do not interact, and thus every orientation and position of the effective particles occur with equal probability, independently of the configuration of the other effective particle. With the effective coordinates defined in eq the probability density in the ideal system, Pid(r, cos θ1, cos θ2, φ) is proportional to r2 and constant in both cos θ and φ. This simple behavior of Pid(r, cos θ1, cos θ2, φ), makes eq a particularly convenient choice of the effective coordinates. In the simulations with two ring polymers, we obtain a probability density P(r, cos θ1, cos θ2, φ) that is different from the ideal distribution Pid. We define as a generalized version of the radial distribution function the anisotropic pair correlation function asThus, the quantity g(r, cos θ1, cos θ2, φ) describes the factor by which configurations in the effective anisotropic model have to be enhanced or suppressed with respect to the ideal case, in order to obtain a distribution for the effective dof that is identical to the distribution obtained in a monomer-resolved simulation in the infinitely dilute case. As in the isotropic case (4) the relation to the associated effective potential readsFrom the anisotropic effective potential, we can deduce the isotropic pair-correlation function viaThe associated isotropic effective potential is then given by (4).

The effective potential between two identical ring polymers remains invariant if we swap the orientations of their respective director with respect to the connection vector, i.e., if we swap the values of the polar angles θ1 and θ2 of the two rings:This symmetry is violated for the effective potential between bidisperse ring polymers, e.g., for ring polymers with a different number of monomers. However, apart from this symmetry, there would be no differences in the procedure of calculating the effective potential between different types of ring polymers.

Simulation Details

Derivation of the Effective Interaction

To determine the effective potential, we carry out constant NVT molecular dynamics simulations with two ring polymers. We simulate rings of N = 20, 50, and 100 monomers. For these simulations we use the LAMMPS simulation package.[61] The polymer rings are placed in a simulation box, which is large enough to prevent multiple interactions via the periodic boundary conditions. The temperature in the simulation is maintained by the use of Langevin dynamics. The corresponding equations of motion read as[62]Here, r is the position of the ith monomer, m its mass, and F(t) is the deterministic force acting on it, which includes the microscopic forces originating from potentials (1)–(3) and the force originating from a biasing potential. The bias potential Vbias = (r – r)2k/2 introduces a harmonic spring with spring constant k between the centers of mass of the ring polymers. The spring is relaxed for r = r. We carried out simulations for different values of r, starting at r = 0 and increasing it up to some maximum value r in steps of σ/2. For N = 20, 50, and 100, r was chosen as 10σ, 20σ, and 30σ, respectively. These values for r are much bigger than the infinite-dilution diameters of gyration (Dg0 = 5.9σ, 13σ, and 21.5σ for N = 20, 50, and 100, respectively). For k we chose the values 2.5ϵ/σ2 and ϵ/σ2 for all ring sizes, and for N = 20 we also carried out simulations with k = 5ϵ/σ2. The quantity η(t) is a random force, with ⟨η(t)⟩ = 0, which is related to the friction coefficient γ by the fluctuation dissipation relation ⟨ηα(t)ηβ(t′)⟩ = 2γmkBTδδαβδ(t – t′), α and β denoting Cartesian components. Our unit of time is set by t0 = (mσ2/ϵ)1/2, and the friction coefficient γ is chosen as 1/t0. We integrate the equations of motion with a time step of Δt = 10–3t0 and use 2 × 108 timesteps for equilibrating the system and collect data during another 2 × 109 timesteps.

We sample histograms P((Q), where Q refers to a bin in the 4D space of the effective coordinates. P((Q) gives the probability for a state in the jth simulation to have effective coordinates in bin Q. The histograms have 128 bins in r and 16 bins in cos θ1, cos θ2, and φ direction. As discussed in section A of the Appendix, we use the self-consistent histogram method by Ferrenberg and Swendsen[63,64] to combine the P((Q) histograms for simulations with identical ring sizes but different biasing potentials to arrive at an estimate for P(Q) in the unbiased system.

Many-Body Effective Fluid

Using the anisotropic effective potential, we carry out standard metropolis Monte Carlo (MC) simulations for the anisotropic effective model. The values of the anisotropic effective potential have been calculated on a discrete grid in the (r, cos θ1, cos θ2, φ) space, and we use linear interpolation to estimate the values of exp(−βVeff) in between the grid points. Having in mind a comparison with both the monomer-resolved simulation results and the isotropic effective potential of ref (52), we choose the same number of particles and effective densities that were used in those simulations. For rings with N = 20, 50, and 100 monomers we simulate systems of n = 2400, 1600, and 1200 rings, respectively, varying in each case accordingly the cubic box size L as to achieve the desired density ρ = n/L3. As the effective potential Veff is bounded, a random distribution of the particles in the simulation box can be used as initial condition. We have implemented two types of MC moves: the first one translates a randomly chosen particle in a random direction, and the second randomly rotates the particle’s director by some angle. The distance by which the particles are displaced and the angle by which they are rotated is randomly selected in an interval starting at 0 and going up to some maximum value. For both moves, this maximum value of the interval is chosen such that the acceptance ratio is approximately 15%. We use 9 × 106 MC moves to equilibrate the system. During this equilibration period the individual soft particles diffuse to several times their own diameter. Afterward, during 15 × 106 MC moves equilibrium configurations are generated. We store the configurations every 20 × 103 moves and use them to compute the physical observables that are presented in section .

The gain in computational efficiency for the simulations in going from a monomer-resolved to a coarse-grained simulation is considerable. The relevant quantity to consider here would be the velocity through phase space. This can conveniently be characterized by means of the mean-square displacement of the rings per unit of CPU time. If one ignores the detailed implementation aspects, the CPU time spent on a single sweep over all monomers in the former and a run over all effective ring particles in the latter, which strictly speaking depend on both the number of monomers N per ring and the overall density of rings, are of similar magnitude. However, the diffusion per sweep in the coarse-grained simulation is significantly larger than that for the monomer resolved simulation; i.e., for the case of N = 50 and ρ* = 20 this results in a factor of approximately 104. The reason for this dramatic improvement is twofold. First of all, the translation/rotation of an effective ring corresponds to a much more time-consuming collective movement of the constituents. The second even more important contribution arises from the steric interaction that are present in the monomer resolved simulations and prevent the unphysical crossing of chain segments. In the coarse-grained simulations such a restriction is absent; i.e., the effective rings are penetrable and can move apparently through each other. On this level of description this effect is not an unphysical process but should be interpreted as a shortcut connecting initial and final configurations that are connected by a much more time-consuming and physically realizable pathway of folding and collective monomer movements.

Anisotropic Effective Potential

We commence by recalculating the isotropic, i.e., angularly average, effective pair potential Veffiso(r) between the stiff rings as a way of comparison with the previously derived results in ref (52). Results are summarized in Figure , reproducing indeed the previously derived ones.[52] The potential for r = 0 is finite, as the rings are allowed to overlap. It also features a local minimum there, whereas its maximum is located at r ≈ 0.25Dg0 for all ring types investigated. Here, Dg0 is the average diameter of gyration of a free ring polymer. The height of the potential barrier at small distances of r decreases by increasing the number of monomers N on the ring polymers.

Figure 2

Effective center-of-mass pair potential in the isotropic effective model for rings with different numbers of monomers N. The center-of-mass separation r is scaled with Dg0, the average diameter of gyration of a free ring polymer. The solid line shows the angularly averaged effective pair potential, while the dashed lines are results of ref (52).

Going over to the anisotropic effective model, we proceed with computing the aforementioned pair correlation function g(r, cos θ1, cos θ2, φ) for a system of two ring polymers. As the latter depends on four effective coordinates and is therefore difficult to visualize, we first introduce a reduced pair correlation function G(r, d(1)·d(2)), which expresses the relative joint probability density of observing the two ring polymers at a distance r and with the directors mutually oriented at the value given by their scalar product, d(1)·d(2), over the same quantity for noninteractive rings. Moreover, we introduce a reduced version of this function, G(r, d(1)·d(2)), by dividing over its isotropic counterpart, i.e.Results are summarized in Figure . One can see that the angular distribution of the directors changes significantly for r ≈ 0.25Dg0, which is approximately the position of the maximum of Veffiso(r). For r < 0.25Dg0 the angle between the directors is biased toward π/2, while for 0.25Dg0 < r < Dg0 they prefer to align parallel with respect to each other. The position of the maximum of Veffiso(r) coincides approximately with the distance r where interpenetrated configurations of the rings become subdominant and where they are more likely to align parallel to each other. The transition between these two domains is particularly steep for N = 20 and becomes smoother for rings with a larger number of monomers. When the rings interpenetrate each other, the distribution of angles between the directors is rather wide, while it gets narrow after the transition where the bias toward parallel alignments of the rings is very strong in particular for the smallest rings with N = 20. It is readily visible from Figure that anisotropy is particularly important for smaller rings. For r > Dg0, the distribution of the angle between the directors becomes flat, as the rings are then well separated and hence do not interact. Note that by definition, eq , the quantity G(r, d(1)·d(2)), is a normalized probability distribution for fixed r, and in Figure it is therefore meaningless to compare the plotted function at different r values.

Figure 3

Infinite-dilution limit of the quantity G(r, d(1)·d(2)), which quantifies the distribution of the scalar product between directors for different values of r. We visualize this distribution for the ring sizes N = 20, N = 50, and N = 100.

The relative orientation between the vectors r, d(1), and d(2) is of course not completely determined by the scalar product d(1)·d(2); the function G(r, d(1)·d(2)) contains less information than the full correlation function g(r, cos θ1, cos θ2, φ). In particular, when the directors are parallel, i.e., d(1)·d(2) = 1, the angle between the connection vector r and the directors d(1) and d(2) is still arbitrary. We denote a configuration with d(1)||d(2)||r as →→ and a configuration with d(1)||d(2)⊥r as ↑↑. From the reduced pair correlation function G(r, d(1)·d(2)) alone, we cannot say which of these two configurations is more probable, as the scalar product d(1)·d(2) is identical to 1 in both cases. Using the full anisotropic pair correlation function g(r, cos θ1, cos θ2, φ), we can compare the corresponding effective potentials:For the →→ case the value of the φ coordinate is immaterial. However, due to the finite bin size of the grid on which we have calculated g(r, cos θ1, cos θ2, φ), the choice of φ makes a small difference, even for V→→(r). We compute V→→(r) from the average of g(r, cos θ1 = 1, cos θ2 = 1, φ) in φ.

In Figure we see that V↑↑(r) increases significantly when r approaches Dg0, while V→→ stays close to 0 until much smaller distances r. We can understand this results if we imagine the rings as disks with diameter Dg0. In the ↑↑ configuration, the rings lie in the same plane and will therefore start to overlap as soon as r ≤ Dg0. Since the rings are not perfect circles and their shape fluctuates, they can feel each other also for distances r which are slightly larger than Dg0. In the →→ configuration two disks overlap only if the distance between their centers of mass is smaller than their thickness. These results tell us that the peak in the reduced pair correlation function g(r, d(1)·d(2)) for r ≈ 0.25Dg0 and d(1)·d(2) ≈ 1 is mostly due to →→ like configurations. However, as soon as the rings can overlap in →→ type configurations the effective potential increases very fast for smaller r, and we come to a regime where other configurations of the directors are more favorable. As a comparison, we also consider a configuration with d(1)⊥d(2)||r, which we denote by ↑→. The corresponding effective potential is given byAs in the →→ case, the value of the φ coordinate is irrelevant for calculating V↑→(r), which we compute from the arithmetic mean of g(r, cos θ1 = 0, cos θ2 = 1, φ) in φ.

Figure 4

Effective potential for three different fixed configurations of the directors and the connecting vector. As a comparison, we also plot the pair potential in the isotropic effective model. (a) N = 20; (b) N = 50; (c) N = 100. The effective potentials are shown only for r values for which we have relatively good statistics. We also show a sketch of the ↑↑, →→, and ↑→ configurations.

For small distances r one ring interpenetrates the other in microscopic configurations of type ↑→. While V↑→(r) starts to increase at larger r values than V→→(r), the increase is slower and converges to a constant for r → 0. This is intuitive to understand since it requires only a finite amount of bending energy to deform two rings such that one can fit into the other. The required bending energy is smaller if the rings are larger. The ↑→ becomes dominant over the →→ configuration at an r value below the threshold r ≈ 0.25Dg0. This is also the r value at which we find the transition in the reduced pair correlation function g(r, d(1)·d(2)) between a regime where configurations with parallel directors, as in →→ , are preferred, to a regime where they are suppressed and other configurations like ↑→ become dominant.

In the Appendix we explain how one can expand the angular part of g(r, cos θ1, cos θ2, φ) into a series of suitably chosen basis functions f(cos θ1, cos θ2, φ). The expression for this expansion is given in eq , and the corresponding coefficients c(r) can be determined by calculating particular ensemble averages as shown in eq . We plot these coefficients c(r) for l1, l2 ≤ 2 in Figure . From the fast change of c(r) for r ≈ 0.25Dg0 one can once more see the transition between two regimes for r in which the distribution of the directors of the ring polymers is very different. We can again see that this transition is smoother for larger rings. The magnitude of coefficients c(r)/c0,0,0(r) with (l1, l2) ≠ (0, 0) tells us about the significance of the corresponding anisotropy in g(r, cos θ1, cos θ2, φ). Anisotropy is more important for smaller rings and becomes more pronounced after the transition at r ≈ 0.25Dg0, where the rings prefer parallel configurations.

Figure 5

First coefficients in the expansion of g(r, cos θ1, cos θ2, φ) divided by the coefficient c0,0,0(r). (a) N = 20, (b) N = 50, and (c) N = 100.

Monte Carlo Simulations of the Anisotropic Effective Model

We carried out Monte Carlo simulations of systems of effective particles described by the anisotropic effective model for different ring sizes N and various densities. We define the reduced density in our simulation as ρ* ≡ nDg03/L3, where n is the number of rings in the sample. In order to assess the quality of the anisotropic effective model, we compare our results to results of full monomer-resolved simulations from ref (52) and the results of simulations using the isotropic effective model. As we can see in Figure for all choices of the number of monomers N the effective models are in good agreement with the full monomer-resolved simulations at low densities ρ*. This is an important consistency check for the effective models, in which the interactions have been chosen such that the distribution of the effective degrees of freedom agrees with their distribution in the full monomer-resolved simulations, in the limit of small densities.

Figure 6

Pair correlation function g(r) at low reduced densities ρ* for a simulation of many ring polymers in the full monomer-resolved simulation (symbols), the anisotropic effective model (solid line), and the isotropic effective model (dashed line).

In Figure , we present results for the smallest rings with N = 20 monomers at higher densities. There is a dramatic improvement of the accuracy as one compares the isotropic with the anisotropic model. While the former fails for ρ* > 2 the anisotropic effective model works up to ρ* ≅ 5 and even gives a semiquantitatively correct description of the system at ρ* = 5.97. At the highest densities, we see the development of a peak in g(r) at r ≅ 0.3Dg0. This peak in the pair correlation function is associated with the emergence of stacks of parallel rings, and its position describes the typical distance of rings in these stacks.[52] Interestingly, the isotropic effective potential has a maximum for r ≈ 0.25Dg0, which is close to the typical distance of the rings in the stacks, and one could wonder why the rings prefer to align at a distance which seems to have a very high free energy penalty according to the isotropic effective potential. The answer to this apparent paradox lies in the strongly peaked nature of the anisotropic effective interaction, which we could observe in Figures , 4, and 5. While the average configuration of the angular degrees of freedom at distances r ≈ 0.3Dg0 has a high free energy penalty, a certain class of configurations, where the directors of the rings are almost parallel, is much more favorable. Obviously, stacking can not be observed in the isotropic effective model, where particles possess no directional degrees of freedom.

Figure 7

Pair correlation function g(r), at high densities, for a simulation of many ring polymers with N = 20 monomers in the full monomer-resolved simulation (symbols), the anisotropic effective model (solid line), and the isotropic effective model (dashed line).

As a further characteristic of the short-range coordination of the rings, we consider the average number ⟨n⟩ of neighbors within a distance r from the center of mass of a randomly chosen ring. This is expressed asFor the rings with N = 20 monomers we present results for ⟨n⟩ in Figure . Once more, the good agreement between the full monomer-resolved and the anisotropic effective model, even at the highest densities investigated, is confirmed: small differences appear only for 0.25Dg0 ≤ r ≤ 0.4Dg0. Evidently, ⟨n⟩ does not contain more information than the g(r) plot in Figure , but it nevertheless clarifies the meaning of the disagreement between the g(r) curves in the monomer-resolved and the anisotropic effective model. At the highest densities in the full simulation, the centers of mass move a bit closer to each other than they do in the anisotropic effective simulation. This manifests itself as a shift of the peaks in the g(r) curves. The difference in the height of the peaks is partly a consequence of the shift, since a peak in g(r) has to be higher at smaller distances if it amounts to the same amount of average neighbors as a peak at a larger distance r. The fact that the ⟨n⟩ curves for the anisotropic effective and the full simulation in Figure agree for r ≥ 0.4Dg0 shows us that the peaks in the g(r) curves indeed correspond to the same amount of average nearby particles that are simply accumulated at slightly different distances.

Figure 8

⟨n⟩ in a simulation of many ring polymers with N = 20 monomers in the full monomer-resolved simulation (symbols), the anisotropic effective model (solid line), and the isotropic effective model (dashed line).

We proceed now with the longer rings, N = 50. As can be seen by the pair correlation curves in Figure a, also in this case the inclusion of anisotropy improves the agreement with the monomer-resolved simulations significantly for densities ρ* from 2.3 to 10.2. In Figure b we present g(r) for higher ρ*. In the full monomer-resolved simulations we can see a peak emerging in the pair correlation function g(r) at r = 0 on increasing the density. As described in refs (52 and 58), the monomer-resolved system forms stacks of quasi-parallel oblate rings that are fully penetrated by bundles of elongated rings. In this phase, the deformation of the penetrating rings is particularly strong. The effective description, on the other hand, breaks down if the internal configurations of the rings in the monomer-resolved system differ significantly to the internal configurations in the system with only two ring polymers. Therefore, the anisotropic effective model should not be expected to be a quantitative description at the high densities in which this phase is formed. Agreement with the monomer-resolved model here is less satisfactory but the improvement over the isotropic model is still spectacular.

Figure 9

Pair correlation function g(r) for a simulation of many ring polymers with N = 50 monomers in the full monomer-resolved simulation (symbols), the anisotropic effective model (solid line), and the isotropic effective model (dashed line). The two plots show different reduced density ρ* ranges.

For the rings with N = 100 we find that anisotropy does not play a key role any more, at least not for the full model at the investigated densities. This had to be expected, as we could already see in Figure and 5 that anisotropy is less pronounced for larger ring sizes. As we saw in Figure c, both the isotropic and the anisotropic model give good results for g(r) up to ρ* ≈ 2.5. In Figure , we see that for higher densities the inclusion of anisotropy does not yield results that are in better agreement with the full monomer-resolved simulations. The results in the isotropic effective model even seem to be in better agreement with the full model, which is attributed to multiparticle interactions that can change the configurations of the large and therefore more deformable rings significantly. The already small correlation between the directors, which is present in the dilute case, might therefore be even smaller at high densities. In the anisotropic effective model, we then overestimate the angular correlations between the directors and arrive at results that can be slightly worse than those of the isotropic model. Interestingly, at ρ* = 20.0, which is the highest density investigated, the anisotropic model appears to crystallize. At this density we see the emergence of columns that are closed over the periodic boundary conditions and organize in a hexagonal 2D lattice structure.

Figure 10

Pair correlation function g(r) for a simulation of many ring polymers with N = 100 monomers in the full monomer-resolved simulation (symbols), the anisotropic effective model (solid line), and the isotropic effective model (dashed line). The two plots show different reduced density ρ* ranges.

Finally, let us focus exclusively on orientational correlations. We define P(d(1)·d(2)) as the probability density distribution for the scalar products between the directors d(1) and d(2) of two ring polymers which are a distance r < 0.6Dg0 away from each other. In Figure , we present results for P(d(1)·d(2)) for simulations in the monomer resolved and in the anisotropic effective model. If the directors were uncorrelated, P(d(1)·d(2)) would be equal to 1. For low densities ρ* we obtain good agreement for all ring sizes investigated. Since we only look at the directional correlation of close by ring polymers, P(d(1)·d(2)) can show strong anisotropic features even for ρ* → 0. As expected, the anisotropy in P(d(1)·d(2)) is stronger for smaller rings. When the density is increased, the distribution always shifts toward parallel configurations in the effective model. This happens because less volume per ring is available for higher ρ* and by aligning parallel the rings occupy less space. Typically one observes the same trend in the monomer-resolved simulation; only for the N = 100 rings we find more parallel rings for ρ* = 2.5 than for ρ* = 17.0. In contrast to the effective model the rings in the monomer-resolved simulation can deform, and their interaction with other rings can therefore be more isotropic at higher densities ρ*. This explains why for the large rings with N = 100 monomers, which deform more easily than the smaller rings, the correlation between the directors is much weaker than in the effective model and can even decrease with density. For N = 50 one can see that the number of orthogonal rings in the monomer-resolved model at high densities is significantly larger than in the effective simulation. As described in refs (52 and 58) for N = 50 and ρ* ≥ 12.8 one observes that oblate rings are interpenetrated by elongated prolate rings. Since the directors of the oblate and the interpenetrating prolate rings can be orthogonal to each other, one observes perpendicular directors for N = 50 even at the highest densities investigated. In the anisotropic effective model, on the other hand, this interpenetration is disfavored, and we observe almost no orthogonal close-by rings at ρ* = 20 for N = 50.

Figure 11

P(d(1)·d(2)) is the probability density to find the scalar product d(1)·d(2) between the directors of two close by rings (r < 0.6Dg0). Here we show P(d(1)·d(2)) for a simulation of many ring polymers in the full monomer-resolved simulation (symbols) and the anisotropic effective model (solid line). (a) N = 20, (b) N = 50, and (c) N = 100.

Truncation of the Expansion of the Anisotropic Potential

Instead of working with a fully tabulated effective potential on a four-dimensional grid, it can be advantageous to use the analytical expansion on basis functions presented in the Appendix. Such expansions are truncated after some term, and here we shortly discuss the quality of such truncations for the problem at hand. To test the quality of the expansion of g(r, cos θ1, cos θ2, φ) we also carried out Monte Carlo simulations, where we used the effective potential associated with the expanded correlation function as the pair interaction between our effective particles. We took the 14 coefficients c(r) for which l1, l2 ≤ 4 into account and truncated the rest of the expansion. While g(r, cos θ1, cos θ2, φ) can never be negative, the truncated expanded version of g can accidentally become smaller than zero. Wherever this happens g = 0 and Veff = ∞ are used in the simulation. In Figure the pair correlation function g(r) obtained in this simulation is shown in comparison with the g(r) function, which we computed previously employing the full anisotropic effective potential. For both N = 20 and N = 50 we obtain reasonable results with the truncated effective interaction, given by only 14 expansion coefficients. For the full effective interaction, which we store on a 4D grid, we save 163 = 4096 entries for each value of r (see section ). At intermediate densities the results obtained with the expanded effective interaction are a significant improvement with respect to the isotropic effective model. However, one has to be aware that for N = 20 the coefficients of higher order modes can still be quite high, especially for r between 0.2Dg0 and 0.7Dg0. In Figure we see that the coefficient for the mode with (l1, l2, m) = (0, 2, 0) can be larger than the coefficient of the isotropic expansion mode. The reason for the high contribution of higher order modes for N = 20 is of course the strongly peaked nature of g(r, cos θ1, cos θ2, φ) for these small rings, which we can also observe in Figure . The convergence of g(r, cos θ1, cos θ2, φ) is poor for the N = 20 rings due to the strong anisotropy of their effective interaction. However, our results show that the expansion modes up to l1, l2 ≤ 4 already capture the main features of the effective interaction. For N = 50 the degree of anisotropy is weaker, and therefore the convergence of the expansion of g(r, cos θ1, cos θ2, φ) is better.

Figure 12

Pair correlation function g(r) for a simulation of many ring polymers in the anisotropic effective model. For the dashed line we expanded the pair-correlation function g before computing the associated effective pair potential. For the expansion we took the 14 coefficients for which l1, l2 ≤ 4 into account. The solid line shows results of a simulation with the unexpanded effective pair interaction. (a) N = 20 and (b) N = 50.

Conclusions

We have introduced a minimal anisotropic model to coarse-grain ring polymers with a finite bending rigidity as soft, penetrable disks. For the shortest (N = 20) and the intermediate (N = 50) sized rings, this model represents a dramatic improvement over the isotropic coarse-graining, in which the relative orientations between the rings are all integrated upon and a radially symmetric interaction results instead. The approach is capable of distinguishing between the relative orientations at infinite dilution, and it carries this distinction also to highly concentrated systems, where it reproduces well the salient features of the structure as seen in the full monomer-resolved simulations. Whereas this is valid more for N = 20 and N = 50, which have a contour length to persistence length ratio of N/s ∼ 2.7 and 6.7, respectively, some important features, such as the penetration of elongated rings in columns formed by oblate rings (found for N = 50), are suppressed or even lost in the effective description, as genuine many-body effects come into play. For the largest rings, N = 100, for which we obtain N/s ∼ 13.3, the contour length is much larger than the persistence length, and they thus resemble more flexible objects. In this case the anisotropic potential at high concentrations fails to describe the structural correlations. This indeed reflects the fact that such rings undergo, at high concentrations, conformational changes (shrinking, interpenetration) that are quite distinct from the assumptions that go into the anisotropic, soft disk model, rendering it thereby very inaccurate. We therefore expect that our anisotropic model yields quantitative results over a broad density range, for systems of polymer rings with a contour to persistence length ratio of N/s ≲ 10.

Our work provides, thus, an accurate and efficient general scheme for the coarse-graining of semiflexible ring polymers, as it allows for a very dramatic reduction of their degrees of freedom, while at the same time introducing a realistic class of systems for which anisotropic generalizations of the ultrasoft, penetrable effective interactions are physically meaningful. Future work will focus on the investigations of the structural and phase behavior of mixtures of stiff rings and of the dynamics of the structure formation in the same.