Event-Chain Monte-Carlo Simulations of Dense Soft Matter Systems

We discuss the rejection-free event-chain Monte-Carlo algorithm and several applications to dense soft matter systems. Event-chain Monte-Carlo is an alternative to standard local Markov-chain Monte-Carlo schemes, which are based on detailed balance, for example the well-known Metropolis-Hastings algorithm. Event-chain Monte-Carlo is a Markov chain Monte-Carlo scheme that uses so-called lifting moves to achieve global balance without rejections (maximal global balance). It has been originally developed for hard sphere systems but is applicable to many soft matter systems and particularly suited for dense soft matter systems with hard core interactions, where it gives significant performance gains compared to a local Monte-Carlo simulation. The algorithm can be generalized to deal with soft interactions and with three-particle interactions, as they naturally arise, for example, in bead-spring models of polymers with bending rigidity. We present results for polymer melts, where the event-chain algorithm can be used for an efficient initialization. We then move on to large systems of semiflexible polymers that form bundles by attractive interactions and can serve as model systems for actin filaments in the cytoskeleton. The event chain algorithm shows that these systems form networks of bundles which coarsen similar to a foam. Finally, we present results on liquid crystal systems, where the event-chain algorithm can equilibrate large systems containing additional colloidal disks very efficiently, which reveals the parallel chaining of disks.


Event-Chain Monte-Carlo Algorithm
Since its first application to a hard disk system [57], Monte-Carlo (MC) simulations have been applied to virtually all types of models in statistical physics, both on-lattice and off-lattice.MC samples microstates {a, b, . ..} of a thermodynamic ensemble statistically according to their Boltzmann weight π a = exp (−E a /k B T ) (in the following, we use thermal units, k B T ≡ 1).One advantage of MC schemes over, for example molecular dynamics simulations, is that it only requires knowledge of microstate energies E a rather than interaction forces.In its simplest form, the Metropolis-Hastings algorithm [29,57], a MC simulation is easily implemented for any system by offering moves from a microstate a to a microstate b of the system and accepting or rejecting these moves according to the Metropolis rule for the Boltzmann distribution, i.e., based on the energy difference between states ∆E a→b = E b − E a .The standard Metropolis rule is defined by the acceptance probability The Metropolis rule (1) satisfies detailed balance for the Boltzmann distribution π a = exp (−E a ) if moves are also offered with a symmetric trial probability p trial (a → b) = p trial (b → a) (which is typically fulfilled trivially for standard local MC moves).A characteristic of the Metropolis rule (1) is also that moves that are energetically downhill, ∆E a→b < 0, are always accepted (p Metr (a → b) = 1).We point out that detailed balance (for symmetric trial probabilities) and this maximal acceptance rate for energetically downhill moves determine the Metropolis rule for the Boltzmann distribution uniquely.The Metropolis acceptance rate for moves that are energetically uphill is exponentially small, however, resulting in frequent rejections.
Typical local MC moves, such as single spin flips in spin systems or single particle moves in off-lattice systems of interacting particles, are often motivated by the actual dynamics of the system.Sampling with local moves can become slow, however, in certain physically relevant situations, most notably, close to a critical point, where large correlated regions exist, or in dense systems, where acceptable moves become rare.
Cluster algorithms deviate from the local Metropolis MC rule and construct MC moves of large non-local clusters, ideally, in a way that the MC move of the cluster is performed rejection-free.For lattice spin systems, the Swendsen-Wang [83] and Wolff [94] algorithms are the most important cluster algorithms with enormous performance gains close to criticality, where they reduce the dynamical exponent governing the critical slowing down in comparison to local MC simulations.These cluster algorithms still fulfill detailed balance but based on a non-trivial trial probability that derives from the cluster construction rules.
Figure 1: The upper scheme illustrates the billiard-type construction of an EC with a total displacement length and with three participating spheres.The lower picture shows the example of a typical EC in a dense hard-sphere system: the coherent movement of long line-like clusters in each ECMC move accelerates the sampling by nearly two orders of magnitude.
The event-chain (EC) MC algorithm introduced by Krauth et al. [5] has been developed to decrease the autocorrelation time in hard disk or sphere systems.It performs rejection-free displacements of several spheres in a single MC move (see Fig. 1) and can, therefore, also be classified as a cluster algorithm.The basic idea is to perform a displacement in a billiard-type fashion, transferring the displacement to the next disk upon a collision, which is called an event and gives the algorithm its name (see Fig. 1).As opposed to cluster algorithms such as the Swendsen-Wang or Wolff algorithms, ECMC algorithms do no longer fulfill detailed balance but only global balance.
In the following, we will give an introduction into the ECMC algorithm starting with the example of simple hard sphere systems.We then discuss how the algorithm realizes the general concept of lifting moves, which is employed in order to achieve global balance.Finally, we review the adaptation and generalization of the ECMC algorithm from athermal systems with hard interactions to systems with soft potentials specified by an interaction energy.Many soft matter systems can be constructed from hard spheres as basic building blocks, eventually with additional interactions, for example spring-like connections to form a polymer chain.We discuss in some detail several applications of the EC algorithm in dense soft matter systems, namely polymer systems, liquid crystal systems as well as mixed systems such as liquid crystal colloids.

EC Algorithm For Hard Spheres
We consider the conceptually simple system of hard (impenetrable) disks (of diameter σ) in two dimensions in order to introduce the ECMC algorithm.Hard disks are the epitome of a system that is easily described and quickly implemented in a (naive) simulation [57], but exhibits a non-trivial phase transition that has been debated for a long time [81].The two-dimensional melting phase transition is believed to be a two-step Kosterlitz-Thouless melting via an intermediate hexatic phase [26,47,95], which is, however, difficult to confirm unambiguously in simulations [4,55].
In a hard disk system, an EC move is constructed by first selecting a random starting particle, an EC displacement direction, and a total displacement length .Both direction and total displacement length are adjustable parameters of the algorithm (see Fig. 1); the sampling of EC directions of directions (uniform or non-uniform) will determine whether detailed or global balance is fulfilled, while the displacement length can be adjusted for performance.We start the EC by moving the chosen particle in the EC direction, which is only possible until it touches another particle after a certain displacement r, which will be typically much smaller than the intended displacement in a dense system.Then, instead of declining such a move as in local MC schemes, the displacement r of the first particle is subtracted from the initial total displacement length and the remaining displacement − r is carried over to the hit particle, which we attempt to displace next by the remaining distance − r.If also the displacement direction is carried over to the hit particle the EC is called straight; this is the case we will focus on in the following and which is illustrated in Fig. 1.This procedure of lifting the EC to the next particle upon collision continues until no displacement length is left.Then an EC has been constructed that moved as a line-like cluster without the possibility of rejection, and we start over by choosing a new starting particle and a new random starting direction to construct the next EC with the same total displacement length .
The total displacement length is an adjustable algorithm parameter, which determines the average number of disks n EC that are moved in an EC move (see Fig. 1).This number is given by n EC ∼ /λ 0 , where λ 0 denotes the mean free path; in general, the efficiency of the EC algorithm depends crucially on n EC .For smaller , the efficiency decreases and approaches traditional local MC (corresponding to the case < λ 0 ).For very large , ECs comprise large parts of the system and give rise to motion similar to a collective translation of disks, which is also inefficient.
We note that the EC algorithm will be very well suited to simulate dense systems but does not allow to simulate maximally dense, i.e., jammed systems since an EC can not displace any particle if all particles touch each other.In jammed hard disk systems, the mean free path λ 0 approaches zero such that the number of participating disks n EC diverges, and the EC does no longer terminate or essentially comprises the entire system resulting in collective translation.This will happen for n EC ∼ /λ 0 > N .In a dense hard disk system of area L 2 and close to the close-packing volume fraction η hcp = π/2 √ 3, we have λ 0 ≈ σ(η hcp − η)/2η hcp resulting in a more concrete condition (η hcp − η)L 2 /σ for the EC algorithm to remain effective close to close-packing.
We will show in the next section that the EC algorithm ensures global balance, regardless of how the EC displacement direction is chosen.If the EC displacement directions are chosen randomly the straight EC algorithm also satisfies detailed balance over many EC events.The most performant version of the straight EC algorithm, however, the so-called xy-version or irreversible straight EC algorithm, breaks detailed balance and uses only displacements along coordinate axes and only in positive direction.Breaking detailed balance but preserving global balance leads to performance gains.Ergodicity must also be ensured, which is done via periodic boundary conditions for the irreversible EC algorithm; moreover, several computations can be simplified for the xy EC algorithm.In Ref. 16, the irreversible EC simulations are roughly 70 times faster than a simulation employing local MC moves.We can obtain similar performance gains in dense soft matter systems.

Balance Conditions And Lifting For Hard Spheres
The standard Metropolis algorithm (1) employs detailed balance to ensure stationarity of the probability distribution of states by pairwise balancing of probability fluxes between microstates in the condition (2).Since circular probability fluxes also fulfill the stationarity requirement for the probability distribution, detailed balance is not a necessary condition and, often, performance can be gained by finding an algorithm that fulfills the weaker condition of global balance, b  ,b,c) and the lifting variable is indicated by color, i.e., the current active particle is indicated by a colored arrow.The EC displacement direction is to the right.Proposed infinitesimal displacements are indicated by opaque spheres.In case of a physical rejection at MC step τ 1 the pivot particle is changed and rejection is avoided.Periodic boundaries ensure ergodicity.One way to read this scheme, is that the lower row is the displacement-reversed version (except trivial translations) of the upper row (obtained be reversing all displacements dw → −dw), resulting in a central symmetry, especially (c,red) = (a,green).
(where the sums are over all microstates), which is simply the continuity equation for a stationary probability π a , i.e., the Boltzmann distribution in our case.While the last equality in (3) is a trivial consequence of normalization, the second equality is the actual global balance condition.All quantities in eq. ( 3) are transition probabilities per MC step; we call the algorithmic time measured in MC steps τ in the following.The algorithmic time τ measured in MC steps must not be identified with an actual physical time t as for any MC simulation.
Global balance algorithms should become particularly performant if maximal global balance is achieved, which means that probability backflow is forbidden leading to the additional constraint In particular, this implies a rejection-free algorithm, i.e., J a→a = 0. Then we can also take a hydrodynamic point of view and see probability density as liquid mass density, unidirectional transition probabilities as liquid velocities, and the probability currents as liquid current densities.Then, maximal global balance corresponds to stationarity of the mass distribution in the presence of fluid flow, i.e., a divergence-free liquid current density.
In Fig. 2, we show how maximal global balance is achieved in EC algorithms for a hard sphere system.In order to investigate balance conditions, it is often convenient to decompose the EC into infinitesimal moves dw, which add up to the total EC displacement in the end [27,59,60].To ensure maximal global balance lifting transitions are introduced, which change an additional lifting variable that enlarges the state space of the system [14].For hard sphere systems, the additional lifting variable is simply the number i of the active/moving particle.This will be the case throughout this article.The microstates are then characterized not only by the positions of all particles but also by specifying the active particle.
For infinitesimal EC moves, only two spheres can collide at once, and it is sufficient to consider a simplified setting of only two hard spheres (green and red) as shown in the cartoon in Fig. 2, where proposed infinitesimal displacements are indicated by opaque spheres.The physical states are labeled with letters (a,b,c), and the lifting variable is indicated by color (red or green particle is active).Physical transitions change the physical states by changing particle positions; this happens by moving the active sphere in the EC displacement direction (i.e., to the right from (a,green) at τ 0 to (b,green) at τ 1 in Fig. 2).If a physical transition leads to an overlap, instead of rejecting the move as in local MC, we perform a lifting transition and change the lifting variable with a lifting probability that exactly equals the rejection probability of the local MC rule in the EC algorithm.In Fig. 2), this means that the active particle is changed to the red sphere if a collision is proposed at τ 1 .Then the red active particle continues to move (to the right from (b,red) to (c,red).
In this way, rejections are avoided, and the entire rejection probability flow is redirected to a lifting move probability flow, which ensures maximal global balance if all physical configurations are equally probable, i.e., if the Boltzmann distribution for hard spheres holds: • The total physical inflow J phys a,green→b,green to (b,green) equals the rejected physical flow J phys b,green→coll (that would lead to collisions on the upper right picture) because it involves moving the same sphere by the same distance dw and because physical states a and b are equally probable.
• By construction of the EC algorithm, all collisions give rise to lifting such that the lifting flow J lift b,green→b,red equals the rejected physical flow J phys b,green→coll .This leads to J phys a,green→b,green = J lift b,green→b,red , so far.
• Apart from a translation, the states a and c and the forbidden overlapping configurations (upper right and lower left pictures in Fig. 2) are identical (please note the periodic boundary conditions).
• Therefore, the missing physical inflow J phys coll→b,red (from the forbidden overlapping configuration on the lower left picture) into state (b,red) equals the rejected physical flow J phys b,green→coll and, thus, the lifting flow J lift b,green→b,red .
• Again, the missing physical inflow J phys coll→b,red into state (b,red) equals the physical outflow J phys b,red→c,red from (b,red) because it involves moving the same sphere by the same distance dw and because physical states b and c are equally probable.
All in all, we obtain J phys a,green→b,green = J lift b,green→b,red = J phys b,red→c,red , which are exactly the maximal global balance conditions for states (b,green) and (b,red) proving maximal global balance for any EC move on a state space that is extended by the lifting variable denoting the active particle.Global balance holds, regardless of how the EC displacement direction or displacement length is chosen.If the reverse EC displacement directions are offered with the same probability (for example by drawing EC directions completely uniformly in all directions) the resulting EC algorithm also satisfies detailed balance over the course of many EC moves.
For pair interactions between hard spheres the lifting variable is the particle number.We can also introduce hard walls to the system, which represent additional steric one-particle interactions.For such one-particle interactions, the lifting variable will not be the particle number (as there are no interaction partners, the active particle must stay the same) but the EC direction itself.It can be shown that classical reflection of the EC direction upon collision with the hard wall will ensure maximal global balance.Analogously to the forbidden overlapping configurations Fig. 2, both the rejected outflow in a collision with a wall and the missing inflow from forbidden configurations are exactly equal to the lifting flow from reflected EC directions.

Soft Interaction Energies
In many applications other than pure hard core systems, we have to deal with systems containing hard core repulsion alongside with other soft interaction energies E = E(r 1 , ..., r N ) that can depend on all particle positions r n in a canonical ensemble.These are not handled by the hard sphere EC scheme described so far.One naive strategy to handle an additional interaction is the following: to ensure ergodicity, an EC is limited to a total displacement , i.e., the sum of all particle displacements, after which a random new particle is chosen as active and a new isotropic direction is set.Interpreting a whole EC as a single MC move, an additional interaction will cause an energy change ∆E when a hard sphere EC (which ignores the additional soft interaction) is executed.
To properly sample the Boltzmann distribution in the presence of the additional interaction, the EC as a whole is then accepted or rejected by a standard Metropolis filter (1).
This strategy of handling additional interaction energies by re-introduction of Metropolis sampling and, thus, rejections will surely decrease the efficiency of the simulation.All advantages of the EC algorithm are effectively lost.There is, however, a way to include arbitrary N -particle interaction energies into the EC framework.Within the lifting formalism it is possible to extend the algorithm from hard spheres to arbitrary pair potentials and continuous spin models by employing factorized Metropolis filters [59,64,67]; this concept can be further generalized to N -particle interactions with N > 2 according to Ref. 27.In principle, these expansions of the EC algorithm presented in the following can be reduced to our simple example of a two-sphere collision from Fig. 2 by generalizing the notion of a collision [67].Effectively, we assign a virtual hard sphere radius to each interaction, which determines which configurations are energetically forbidden and trigger a generalized collision.This can be done by introducing the so-called rejection distance, which is further explained in Fig. 4.
In a pure hard sphere simulation, a sphere is continuously displaced until it hits another sphere which triggers a lifting event (or the total displacement reaches ).We can resolve all events into unique binary collisions or rejections, which then leads to a simple lifting probability flux during each event (see Fig. 2).When a particle is displaced in the presence of continuous interaction energies, it is not clear which interaction causes the rejection, since the rejection in the standard Metropolis filter is based on the sum of all changes in all interaction energies (pair or many-particle interactions).In fact, the rejection is caused by all interactions at once.By using an infinitesimal factorized Metropolis filter one can handle each interaction independently.
Each particle i has a set S i of interactions with other particles; each interaction energy E ∈ S i depends on the position r i of particle i and on the positions of other particles in the set M E of interaction partners participating in E. Pairwise interactions, such as elastic springs E = E spring (ij) in a bead-spring polymer or a van-der-Waals interaction between particles, are interactions with one other particle j (M Espring(ij) = {j}; all particles j = i can be interaction partners); three-particle interactions such as a bending energy in a polymer are interactions E = E bend (jik) with two other particles j, k that form two neighboring bonds with i and thus define a bond angle (M E bend (ijk) = {j, k}; all pairs jk with j = i, k = i can be interaction partners).
In the following, we consider again infinitesimal displacements dw of particle i in the EC direction.Such an infinitesimal move changes the configuration from b to a and leads to the total energy change dE i,b→a , which is the sum of all changes in interactions in S i , dE i,b→a = E∈Si dE i,b→a .
We construct a global balance algorithm starting from a Metropolis filter for the acceptance of physical moves based on the energy dE i,b→a .For soft interactions we switch, however, from the standard Metropolis filter (1) to a factorized one [59,64,67], where the Boltzmann weights for the sum dE i,b→a = E∈Si dE i,b→a are factorized and the Metropolis acceptance rule is applied to each factor separately.This is equivalent to switching the sum and maximum operation on the energy changes regarding each interaction, It will be important that the factorized Metropolis filter still has the detailed balance property (2).The probability of rejecting the next infinitesimal step 1 − p fact is proportional to the sum of each infinitesimal energy change of each interaction.If one interaction causes a rejection, the whole move is rejected.Since we employ infinitesimal steps, only one interaction can reject at once.In analogy to the hard sphere case, these rejection events can be seen as a generalized collision events.In other words, the factorized Metropolis filter uniquely assigns a rejection to one specific interaction at the expense of an increased rejection rate, because energy increase from one interaction cannot be compensated by an energy decrease of another interaction.
The factorized Metropolis filter can now be used to construct a maximal global balance algorithm on an enlarged state space by redirecting all physical rejection events into lifting events.This is done analogously to the hard sphere case outlined above.If a rejection should take place according to the factorized Metropolis filter, we perform a compensating lifting move instead.We lift to one of the particles that participate in the interaction that caused the rejection.For a pair interaction this already determines the particle to lift to uniquely, for many-body interactions an additional rule is needed to select the particle to lift to properly [27].Therefore, we first discuss the simpler case of pair interactions and show how we can realize maximal global balance.

Pair Interactions
For states in the enlarged state space, we use the notation a i for a physical state a with an active particle i.In the following we try to move the active particle i by an infinitesimal step dw along the EC direction d resulting in a physical move from a i into b i .For pairwise interactions E(ij) (i.e., S i = {E(ij)|j = i}) the transition rate is given by the factorized Metropolis filter This transition rate is decreased by rejections if a move has dE ai→bi (ij) > 0, i.e., is energetically uphill.If p phys (a i → b i ) > 0 because configuration b i can be reached by a move dw of particle i from a i , it follows that p phys (b i → a i ) = 0 because the EC direction is fixed, and we cannot reach a i from b i by moving i into the same direction dw (we would need to move into the opposite direction −dw).Therefore, physical transitions are unidirectional.
Because only one pair interaction E(ij) can reject at once in the factorized Metropolis filter, it is sufficient to consider the two interacting particles i and j in deriving the lifting probabilities ) from global balance, see Fig. 3.We want to apply the global balance condition (3) to incoming physical and lifting flows to state b i .Using detailed balance of the factorized Metropolis filter we obtain for the physical inflow from a i to b i lifting variable Lifting moves from b j to b i are triggered by a move dw of particle j into configuration b j from a different configuration c j , which is obtained from b j by displacing particle j by −dw (c j is analogous to the forbidden configuration on the lower left in Fig. 2).

MC time
with π bi = π bj = π b /N by symmetry in a N -particle system.Using the global balance condition (3) for flows to state b i , and dividing by π bi we obtain the lifting probability as : Event-driven approach to determine the rejection distance for a pair interaction between i and j.
Condition (13) has to be solved for the rejection distance w ij at which an event/lifting occurs.The black curve is the interaction energy along the EC direction.Only on red segments with increasing energy a rejection can occur.Among all pair interactions, the rejection will be assigned to the interaction with the smallest rejection distance.
The last equality holds for a translationally invariant pair interaction, where dE cj →bj (ij) = −dE ai→bi (ij).In the other lifting direction, this implies which leads to a rejection-free algorithm because the rejection probability of the physical moves ( 6) is exactly redirected into a lifting probability.We also conclude that , also lifting transitions are unidirectional proving maximal global balance.
Global balance is established for each move a i → b i independently This generalizes the picture from Fig. 2 to soft interactions.In fact, for hard spheres [dE ai→bi (ij)] + = 0 before a collision and [dE ai→bi (ij)] + = ∞ at collision such that lifting to the colliding particle j happens with probability one upon collision, exactly as in Fig. 2.
For an efficient implementation, an event-based approach is chosen [59,67], in which we determine the distance to the next rejection.This rejection distance corresponds to the distance to the next collision for hard spheres.
The probability p(w ij ) to reject and lift from particle i to particle j because of the interaction E(ij) and after moving a distance in the interval [w ij , w ij + dw] is given by the probability to not encounter a rejection/lifting up to w ij and then lift with probability p lift E(ij) (w ij )dw = [dE ai→bi (ij, w ij )] + , see eq. (11).As a result, we obtain a Poisson-type distribution for the rejection distance w ij of particle i caused by interactions with particle j, which depends on all uphill energy differences caused by particle j if particle i is moved, see Fig. 4. Interpreting −∂ w E(ij, w) as force onto particle i exerted by particle j, the rejection length distribution depends on the line integral over all opposing forces.
Transforming the probability distribution (12), we find that the distribution of the usable (i.e., uphill) energy + is a simple exponential.It follows that the lifting distance w ij caused by particle j can be determined with the proper distribution by drawing a probability u ij = p(∆U ij ) from a uniform distribution in [0, 1], which translates into an exponentially distributed ∆U ij = − ln u ij , and by solving the equation for w ij .This is illustrated in Fig. 4. Because all interaction energies E(ij) with different partners j are independent, rejection probabilities from the factorized Metropolis filter are additive by construction, and we have a Poissonian rejection length distribution, we can determine a rejection distance for all possible interaction energies, and the shortest distance triggers the next lifting event.This is analogous to the first-reaction method in a Gillespie algorithm [23].In this sense, the rejection distance is a virtual collision radius, which can be assigned to each interaction, and the shortest collision radius triggers a generalized collision.This procedure based on (13) gives a simple and fast event-driven ECMC algorithm for arbitrary pair potentials.

N -Particle Interactions
In many applications, also three-particle interactions occur.This happens, in particular, for extended objects, such as rods or polymers, which can be described by bead-spring models.One prominent example are semiflexible polymer chains with a bending energy.Because the local bending angle involves three neighboring beads in a discrete model, the bending energy is a three-particle intra-polymer interaction in terms of bead positions.
In Ref. 27, the EC algorithm has been generalized to three-and many-particle interactions, thus broadening the range of applicability of rejection-free EC algorithms considerably.For N -particle interactions E, there are N − 1 interacting particles to which the EC can lift to avoid rejections.
In order to redirect the rejection probability [dE ai→bi ] + from physical moves into lifting moves and, thus, a rejection-free algorithm, lifting to the set M E of N − 1 interaction partners should be done with the analogous probability as for pairs, see eq. (11).We need, in addition, a set of conditional lifting probabilities λ ij to assure maximal global balance.The λ ij specify the probabilities to lift to one of the N − 1 interaction partners j ∈ M E ( j∈M E λ ij = 1); the total probability to lift to j becomes We first consider the important case of a three-particle interaction (between particles i, j, k), see also Fig. 5. Analogously to eq. ( 10), the global balance condition (3) applied to flows to state b i requires a backward lifting probability Inserting (15) this results in the condition lifting variable

MC time
Figure 5: Conditional lifting in a system with three-particle interaction.The conditional lifting probabilities can be derived from the balance condition.Analogously to Fig. 2 (i.e., using eq.( 5) and dividing by π b ) and decomposing the lifting probabilities according to eq. ( 15) and choosing λ red→green ≡ λ green , we get ).Note that this can be extended to match eq.(20).
As for pairs, lifting from particle i to j will only take place (λ ij > 0) if dE ai→bi > 0 giving rise to a rejection of a i → b i , and if the corresponding dE cj →bj < 0 because of eq. ( 17).Translational invariance of the interaction implies dE ai→bi + dE c k →b k + dE cj →bj = 0. Writing out eq. ( 17) analogously also for lifting to j and k, we can determine all conditional lifting probabilities uniquely [27], where Θ(x) is the Heaviside step function.
For arbitrary N , the global balance condition (17) generalizes to For N > 3, the choice of conditional lifting probabilities is not unique [27].They can be fixed by the additional requirement λ ji ≡ λ i , that they are independent of the particle j we are lifting from as in (18).This leads to the simple general result Regardless of the number of interacting particles N in an interaction in S i , the first step of calculating the rejection distance does not differ and is given by (13).For N = 2, the lifting destination when a rejection occurs is unambiguous, hence λ ij = 1.For N ≥ 3, the probability to lift from the active particle i to one of the other interacting particles j ∈ M E is also proportional to λ ij ∝ [−dE j ] + (see eq. ( 20)), which can be interpreted as the force which is exerted on the particle j by the interaction.Note that the conditional lifting probability vanishes when the particle, one would lift to, will increase the energy when moving along the current EC direction.We want to highlight that due to symmetry the sum of all forces must vanish j dE j = 0 and i can only trigger an event when dE i > 0 holds, which implies that for at least one other particle dE j < 0 must hold.In short, one most probably lifts to the particle, which decreases the interaction energy the fastest.

One-Particle Interactions
Finally, soft external one-particle potentials E(i) frequently occur in soft matter systems, for example, as confining potentials, gravity, or external electric fields for charged particles.Such one-particle potentials are handled analogously to hard walls.Because there are no interaction partners the lifting variable is the EC direction d.
The lifting probability for lifting the EC direction to a reflected direction d during an attempted move a i → b i .The reflection can be performed with respect to the equipotential surface of E(i), i.e., by lifting to d = (1 − 2P )d, where P is the projection operator onto the ∇ i E(i)-direction.This has been shown to satisfy global balance [67].

ECMC Algorithm
In summary, an ECMC simulation, which we use for soft matter systems, can be structured as follows.Starting from a random particle i in a random EC direction d and with a certain EC total displacement , the rejection distance w ij for each interaction E ∈ S i is calculated from eq. ( 13) and as shown in Fig. 4. The shortest distance determines which interaction caused the rejection.If the interaction consists of more than two particles, the conditional lifting probabilities λ ij are calculated from eq. ( 20) and as illustrated in Fig. 5. Then we lift to the corresponding interacting particle j after the corresponding rejection distance w ij .For one-particle interactions causing the rejection we lift the EC direction by reflection.
1. choose a random EC direction d, a total displacement and an active particle i , 2. calculate the minimal rejection distance r M E = min E ∈Si w iM E for each E according to eq. ( 13), 3. if rejection is triggered by a three-or more particle interaction: select j ∈ M E by applying conditional lifting probabilities from eqs. (18) or (20) (see Fig. 5), 4. move particle i: r i = r i + rd, 5. lift to rejecting particle i → j (or to new reflected EC direction d if rejection is triggered by a one-particle interaction), 6. let = − r, 7. if = 0 goto 1 else goto 2

Additional ECMC Simulation Features
We want to highlight some useful properties of an ECMC simulation.First, using the factorized Metropolis filter each interaction is independent, making it easy to implement the algorithm in a highly modular manner.
Furthermore, the rejection distance distribution and lifting statistics is intimately linked to the pressure in the system.Without loss of generality, let the EC direction be in x-direction, then the pressure can be expressed by [59] where . . .denotes an ensemble average over multiple ECs with length .Rewriting x final − x initial = + events (x j − x i ), where x i,j are the x-coordinates of particles i and j, one sees that the excess length per EC length determines interaction specific contributions to the pressure of an ideal gas, i.e., attractive interactions with a negative excess length decrease the pressure, and repulsive interactions increase the pressure.This procedure does not need any additional computations.For hard spheres, the pressure is measured far more efficient than with the usual detour via the contact theorem.
Every configuration visited is weighted properly even if a hard sphere collision just happened.This can be exploited by special MC moves.This is not unconditionally true in a molecular dynamics simulation where, for instance, the statistics is skewed if measurements are taken on particle collisions.Inclusion of special many-particle MC moves can be beneficial for particular systems.For polymer melt simulations, we introduce a bead-swapping move [36], which switches positions of two hard spheres upon collision obeying a standard metropolis filter for the spring energies involved, as explained later.
Many simulation techniques based on single particle moves can benefit from massive parallelization if the simulation system can be efficiently divided into independent pieces.For ECMC algorithms, massive parallelization will not be optimal because of the cluster nature of EC moves [37].In Ref. 37, a parallelization scheme for the EC algorithm was introduced, and extensive tests for correctness and efficiency were performed for the hard sphere system in two dimensions.For optimal parameters the algorithm achieves a speed up by a factor of roughly the number of cores used compared to a sequential EC simulation.For parallelization we use a spatial partitioning approach into simulation cells.We analyzed the performance gains for the parallel EC algorithm and find the criterion ∼ Lλ 0 / √ nσ for an optimal degree of parallelization, where n is the number of parallel threads, L the system size, and η the occupied area fraction.If the number n of parallel threads is chosen too large (massive parallelization), this choice for will drop below the optimal window < 10λ 0 .Because massive parallelization is not optimal, the parallel EC algorithm will be best suited for commonly available multicore CPUs with shared memory.This ECMC parallelization scheme can also be applied to other soft matter systems such as polymer melts [36].(A) Scheme for the rejection distance for hard sphere chains.If two bonded spheres do not collide, the active sphere can move freely until the minimal distance between both is reached.A usable energy ∆E is drawn (see Fig. 4) which determines where the rejection takes place.(B) Scheme of a swap move.In case of a hard sphere collision, the energy difference of each involved springs is calculated when the beads would switch positions.The swap move is accepted with a standard Metropolis filter, otherwise the EC is continued normally.The active particle is marked by a red circle.

Applications
We now briefly present some soft matter systems where we successfully employed ECMC algorithms, namely a dense polymer melt, a network-forming system of polymers with a short range attraction, hard needle models of liquid crystal and liquid crystal colloids from a mixture of needles and hard spheres.In all of these systems the ECMC algorithm speeds up simulations considerably at high particle densities as compared to standard local MC simulations.The simulations for these soft matter applications have been performed within the highly modular polyeventchain framework, which is particularly suited for soft matter applications and will be presented in detail elsewhere.Jellyfish provides a similar EC framework with a focus on all atom simulations [17,30].

Initialization And Simulation Of Polymer Melts
The simulation of polymer melts by Molecular Dynamics (MD) or MC simulations is a challenging problem, in particular, for long chains at high densities, where polymers in the melt exhibit slow reptation and entanglement dynamics.For chain molecules of length N the disentanglement time increases ∝ N 3 [15], which makes equilibration of long chain molecules in a melt a numerically challenging problem if only local displacement moves of polymer segments are employed as in a typical off-lattice MC simulation with fixed bond lengths [3,9,28,41] or fluctuating bond lengths [73], or in MD simulations [1,49,69], or for highly confined and therefore extremely packed melts [72].In MD simulations, reptation dynamics has been successfully identified in the entangled regime [1,31,49,69].In MC simulation, reptation dynamics has been first identified in lattice models [48] or fluctuating bond lattice models [66] with Rouse-like local bond moves.
In MC simulations, non-local or collective MC moves can be introduced, for example, chain-topology changing double-bridging moves [40], which speed up equilibration; dynamic properties, however, are no longer realistic if such collective MC moves are employed.Likewise, reptation can also be introduced as explicit additional local MC reptation move [28,88] (slithering snake moves) to obtain faster equilibration of a polymer melt but reptation dynamics is no longer dynamically realistic then.
In order to apply the above scheme of ECMC simulations to polymer melts, flexible polymers are modeled by a bead-spring model consisting of hard spheres, which are connected by harmonic springs [36].These are the only additional soft pair interactions in the ECMC scheme for a flexible polymer melt.We choose the bond rest length to equal the hard-sphere diameter σ and a sufficiently large bond stiffness such that bond-crossing becomes inhibited by the impenetrable hard spheres and polymers remain entangled.
ECMC simulations can speed up MC simulations of polymers such that simulation speeds become comparable to MD simulation speeds and such that reptation dynamics can be clearly observed [36].In each EC move, all beads that would collide successively during a short time interval in a MD simulation are displaced at once.This gives rise to a ECMC dynamics which is effectively very similar to the realistic MD dynamics and statements about polymer dynamics are still possible.In Ref. 36 we compared the ECMC algorithm to MD simulations implemented in LAMMPS regarding the equilibration of polymer melts.We utilized the time, that the positional fluctuation of the most inner bead of each polymer needs to reach the diffusive regime, as gauge to introduce an ad-hoc interpretation of time in the MC simulation (see Table I and Fig. 5 in Ref. 36).This showed that LAMMPS equilibrates more efficiently for moderate chain lengths.For long polymers (N = 500) and when using an additional swap-move, ECMC simulations become equally performant, suggesting clear advantages for melts with even longer polymers.
The additional swap-move chanes the topology of entanglements locally.In contrast to the double-bridging move, which changes bonds, the swap move changes topology by changing bead positions.For this purpose we modify the EC move so that the EC does not directly transfer to the next bead upon hard sphere contact but, instead, a swap of the two touching spheres is proposed, see Fig. 6.Such an additional swap move is EC specific and allows for a local change of entanglements.An analogous swap move in a standard MC algorithm needs to select pairs of beads such that the swap move has a reasonable acceptance rate (the particles have to be reasonably close).Moreover, the selection rule has to satisfy detailed balance (for example, simply proposing the nearest neighbor for swapping will lead to a violation of detailed balance).Therefore, there is no straightforward analogue of the EC swap move in a standard MC simulation with local moves.
Here, we want to focus on another important aspect of EC moves in polymer melt simulations.EC moves can also be employed to generate initial configurations that are already representative of equilibrium configurations.This keeps the equilibration or warm-up phase of any subsequent simulation scheme (MC or MD) short.An early strategy to generate initial conditions is to slowly compress a very dilute equilibrated melt, which can take a long time for very dense target systems and becomes more difficult when chain lengths increase [49].Another method that has been proposed first distributes the chains without interaction (phantom chains) and tries to resolve overlaps in a way which distorts the chain statistics as least as possible [49].Over the years more sophisticated ways of resolving overlaps have been proposed, for instance a slowly increase of the strongly repulsive excluded volume potential [1].This procedure is called push-off and the rate at which the potential is increased and the potential range is extended determines to a large extent how much the original statistics of the chains are disturbed.The method can also be applied to more realistic polymer models [78].Recent work by Moreira et al. [61] shows how the procedure of Auhl et al. can be applied more efficiently and how the equilibration time can be shortened roughly by a factor of 6.An alternative strategy that performs even better for very long chains (again by a factor of ∼ 5) has been introduced by Zhang et al. [96].They employ a hierarchical ) captures how well the polymers are equilibrated and has to be as close to the equilibrium result as possible because this observable equilibrates very slowly.For equilibration a chain needs to move its own size, which is slow due to reptational dynamics.Because of strong fluctuations the warm-up times have not been measured precisely, but are of the order of a couple of hours.(D) The tangent correlation shows the correct non-ideal behaviour of chains in the melt [92,93].The configuration after rattling resembles the equilibrium quite well.The warm-up takes only a couple of minutes.
approach, which uses sequential backmapping from a coarse-grained representation in order to re-introduce molecular details.
We present an EC-based scheme that is similar to the push-off scheme and takes advantage of the fact that the EC algorithm is well suited to resolve overlaps.To do this, ECs are started in random directions on any sphere that overlaps with at least one other sphere until all overlaps of this sphere are resolved.This procedure ignores spheres that overlap with the active sphere, so that the algorithm itself does not need to be changed.This procedure corresponds to locally rattling the polymer system, and can create enough space around the overlapping bead to insert it.
For relatively dense melts, Flory [18,19] puts forward the hypothesis that the volume exclusion effect (chain swelling, R 2 ∼ N 2ν with ν ≥ 1/2) is exactly compensated by the pressure generated by surrounding chains, so that a chain shows ideal behaviour again.More recent results [92,93] show that this hypothesis is only a good approximation.In particular, the tangent-tangent correlation of long (N > 500) polymers in a melt falls algebraically, which deviates significantly from exponential decorrelation happening in the single polymer, i.e., free case.
This motivates the following scheme to generate pre-equilibrated polymer melts.Since excluded volume interactions between parts of a given chain are more or less screened as Flory suggests, we place phantom/ideal chains into the system in configurations that resemble non-reversal configurations by excluding a certain angular region for backward steps in generating the initial phantom chain configurations as described in Ref. 1.This captures the structure of the chains on all length scales quite well, but comes with the downside that the asymptotic behavior of the end-to-end distance (the asymptotic Flory ratio R 2 /(n b 2 )) must be known beforehand to tune the excluded angular region accordingly.Luckily, the internal structure of the chains depends not strongly on the chain length and the asymptotic behavior can be determined in a faster system of shorter chains (here, we carried out one long run.Theoretical and numerical values for this ratio can be found in Ref. 21 and references therein).
The resulting phantom chain system is then equilibrated for a short time, so that the phantom chains will relax towards ideal chain configurations, i.e., the phantom chains contract starting from the initial non-reversal configurations.Then the hard sphere diameter is increased from σ = 0 (phantom chains) to the target σ at once1 , which gives rise to overlaps.Then we use the EC rattling moves to resolve overlaps in the system.Switching on the hard sphere diameter will swell the chains and, therefore, have the opposite effect to the pre-equilibration.Both steps cancel each other well, resulting in high quality initial conditions.The quality can be assessed by comparison of equilibrium values of observables, such as R 2 , pair distribution g(r), or bead distributions, with ensemble values given by the initial condition generation as shown in Fig. 7.
For this procedure it is vital that the spring constant is quite large.Starting ECs on one particle until overlaps are resolved obviously breaks balance.Presumably, setting the spring constant to a relatively high value (in Fig. 7 k = 1000k B T /σ 2 , where the rest length is identical to the hard sphere diameter b 0 = σ) leads to less distortions along the chain.
In general, it is hard to quantitatively compare time scales and thus performance of the EC-based initialization scheme with other schemes from the literature.Hardware development over the years (Moore's law, number of cores, etc.), underlying polymer models (volume exclusion: hard spheres vs. Lennard-Jones, bonding: harmonic springs vs. FENE vs. tangent hard spheres vs. united atom force fields), criteria of deciding when equilibrium is reached, or simply untested parameter dependencies (regarding system size, chain length, spring constants) can have non-trivial influence, even on an otherwise robust benchmark.Despite these problems to compare quantitatively, the proposed EC-based method seems comparable to other state-of-the-art methods.This has to be checked further in future work.

Self-Assembly Of Filament Bundle Networks
The cytoskeleton is a mechanically important structure which consists of three classes of interacting filamentous proteins (microtubules, filamentous F-actin, intermediate filaments) in animal cells.The different constituents reflect the multi-functionality of the cytoskeleton and the different requirements with respect to force and time scales.F-actin structures in the cell cortex are most important for cell mechanics and cell motility [6,32].From the polymer physics point of view, F-actin is a semiflexible polymer since typical contours lengths are of the same order as its persistence length.Therefore, bending rigidity and thermal fluctuations are relevant to understand F-actin structures.Within the cell, F-actin forms locally very dense sub-cellular structures like bundles [2], networks and even networks of bundles to fulfill their biological functionality [32].In vivo, bundles and also networks are held together by crosslinking proteins [54].
Also minimal cytoskeletal in vitro systems [54], some of which are confined to droplets [33,34] or in microfluidic chambers [12,13,82], show formation of bundles and also networks of bundles under the influence of attractive interactions.In vitro, such attractive interactions can not only be induced by crosslinkers but also by counterions and depletion interactions.
The simulation of the cytoskeleton of a cell is a long-standing problem as its physical properties are often governed by the interplay of many long polymers, which are deformable by thermal fluctuations against their bending rigidity and subject to crosslinker-mediated attractive interactions [44].Theoretical work on crosslinker-mediated bundling of semiflexible polymers [35,42,44,45] and the related problem of counterion-mediated binding of semiflexible polyelectrolytes [7] show a discontinuous bundling transition above a threshold concentration of crosslinkers or counterions.The nature of the crosslinkers can influence the mechanical properties of polymer bundles [77].
Simulations of few filaments clearly confirm bundling in a single transition [44].The emerging bundles of semiflexible polymers (Fig. 8) are typically rather densely packed which causes difficulties in equilibrating bundled structures in traditional MC simulations employing local moves.In Ref. 44, MC simulations showed evidence for kinetically arrested states with segregated sub-bundles; in Ref. 80, kinetically arrested bundle networks have been observed for rather small semiflexible polyelectrolyte systems.Further numerical progress requires a faster equilibration of dense bundle structures.Simulations of larger systems consisting of cytoskeletal filaments and crosslinkers [8,10,20] show different competing phases from isotropic networks to bundles and other aggregates, for example, with bond-orientational order.In Ref. 20, it was also pointed out that the dynamics of bundling and polymerization in combination with entanglement of polymers strongly influence the resulting structure.This raises the question to what extent bundled structures can reach a genuine equilibrium under realistic conditions.
Further numerical progress requires a faster equilibration of dense bundle structures, for which the ECMC algorithm is ideally suited [37].In comparison to polymer melts, two additional interactions have to be included in this system: (i) Polymers are semiflexible such that we have to include a bending rigidity of each bead-spring polymer, which is an example of a three-particle interaction between three neighboring hard sphere beads along a polymer.(ii) There is a short-range attraction, which mimics a crosslinker-or counterion-mediated attraction between polymers.Therefore, all beads also interact pairwise with a short-ranged attractive square well potential of strength g and range d.We choose a potential strength g well above the critical value for bundling [35,44], and the range of the attractive square well potential is d = 0.4σ, i.e., comparable to the filament radius (bead radius).Both additional interaction energies are incorporated in the ECMC polymer simulation.Again, it is important to note that we expect a qualitatively realistic MC dynamics from ECs, which resembles collisions in MD simulations.The performance can be further increases by parallelization as discussed above [37].
We have applied this ECMC algorithm to a system consisting of many2 interacting semiflexible polymers, such as actin filaments, in a flat simulation box in three dimensions.Periodic boundary conditions are used only for the extended direction of the simulation box whereas the spatially short direction is non-periodic.Under the influence of the short-range attraction, we clearly observe formation of densely packed bundles.We also observe that, starting from an isotropic melt-like situation, it is not a single bundle that forms but the system evolves into a structure consisting of a network of bundles.
The simulation does not reach a truly equilibrated stationary state, but the network of bundles keeps evolving by coarsening processes as shown in Fig. 8A.This clearly demonstrates that the assembly dynamics plays an important role in structure formation in these systems, as similarly suggested in Ref. 20.As mentioned above, the ECMC dynamics is effectively very similar to the realistic MD dynamics.Therefore, the observed coarsening dynamics is not an artefact of the MC simulation but a generic property of this multi-filament system.On the time scales of our simulations (up to 10 7 MC sweeps) the bundle system keeps evolving by coarsening and no genuine equilibrium state is reached.
Our preliminary results suggest that the bundle networks form structures similar to foams [65], which also continue to coarsen over time, see Fig. 8.The dynamics of the coarsening process exhibits further analogies.It can be observed that cells in the bundle network, which enclose a comparatively large area, tend to increase their area at the expense of cells with smaller enclosed area.This behaviour is known from foams as a consequence of the von Neumann law for the area growth rate of cells [24] and of the correlation between relative bubble volume and the number of sides [25].However, the von Neumann law assumes an exchange of the enclosed medium by diffusion through the liquid interfaces and cannot be applied directly to the bundle networks considered here.Other empirical laws found in the context of foams like Aboav-Weaire law [89] for the mean number M (n) of sides of cells surrounding an n-sided cell can also be applied very well to the emerging bundle networks (see Fig. 8D).The Aboav-Weaire law is of particular interest since it provides an alternative approach to determine the variance µ 2 of the mean number of edges per cell (denoted as P (n) in Fig. 8C).Comparison of the value of µ 2 computed directly from the data to the values obtained from fit parameters of P (n) and the Aboav-Weaire law yields a reasonable agreement.
While the structure formation in dry foams is driven by the minimization of interfacial energies and diffusion [90], it appears that the bundle networks coarsen by a zipping mechanism [43,50] which leads -similar to the decrease in interfacial energy in foams -to a decrease of the polymer binding energy (i.e., to more adhesive contacts between polymers).At a fixed amount of polymers this gives rise to a minimization of bundle length, which plays an analogous role to the minimization of surface area in a classical three-dimensional foam.The foam-like bundle networks resemble the structures observed in different in vitro experiments [12,33], where also polygonal cell structures of the self-assembled networks have been observed.

Liquid Crystals And Colloidal Suspensions
Composite soft matter systems such as colloidal mixtures containing different colloidal particles, often of different size, represent challenging systems for simulations because their physics are typically governed by effective interactions which arise if microscopic degrees of freedom of one species are integrated out.Effective interactions are essential to characterize stability and potential self-assembly into crystalline phases but hard to access in a microscopic particle-based simulation.The process of integrating out microscopic degrees of freedom corresponds to the numerical evaluation of a potential of mean force between the colloidal species of interest.In order to measure the potential of mean force for one colloidal species accurately, all degrees of freedom must be properly equilibrated.
As an example of such a colloidal mixture we consider a two-dimensional system of needles and colloidal disks.This serves as a simple model system for liquid crystal (LC) colloids, which are colloidal particles suspended in a liquid crystal.Such LC colloids exhibit anisotropic effective interactions between colloidal particles if the LC is in an ordered, e.g., nematic phase [53,79].A nematic LC phase forms a rich variety of defect-structures around a spherical inclusion such as Saturn-ring disclination rings or a satellite hedgehog for normal anchoring and boojums for planar anchoring [51,75,87].In a dense nematic LC, the effective interactions can be governed both by depletion interactions [52] or by long-range elastic interactions mediated by director field distortions in the nematic hard needles.Because hard needles tend to align tangentially at a hard wall, we expect a quadrupolar elastic interaction, which is characteristic for planar anchoring at the colloidal disk [62,68,71,74,86] but also generic in two dimensions [85].This interplay has been studied in Ref. 63 employing an ECMC simulation.
The colloidal mixture of hard disks suspended in a nematic host of hard needles is particularly challenging as the hard needle system must be fairly dense to establish a nematic phase.While particle-based simulations exist for dilute rods in the isotropic phase [46,76], the regime of a nematic host is fairly unexplored up to now and simulations resorted to coarse-graining approaches.[68] So far, only single inclusions [70] or confining geometries [22] have been investigated by particle-based simulations.In order to measure the potential of mean force for larger colloidal disks accurately, all degrees of freedom of the needles but also the relatively slow degrees of freedom of the colloidal disks must be properly equilibrated.This is efficiently achieved by the ECMC algorithm.
The idea for the ECMC simulation of a needle system is to represent needles by their two endpoints, where only one endpoint is moving at a time [27]; the endpoints are connected by an infinitely thin, hard line.For efficient sampling, the distance of the two endpoints, i.e. the length l of the hard needle can fluctuate around its characteristic length l 0 in order to allow for independent motion of both endpoints in the MC simulation; the needle length l is restricted by a tethering potential V n (l) with V n (l) = 0 for l/l 0 ∈ [0.9, 1.1] and infinite V n (l) else, such that l 0 is the effective needle length.This constitutes an additional pair potential in the ECMC framework.
In two dimensions, the needle-needle interaction simplifies effectively to a collision of an endpoint with another needle.The remaining MC move distance is lifted to one of the endpoints of this needle.Therefore, we have a fluid of endpoints with an effective three-particle interaction (two endpoints of a passive needle and the active endpoint).For needles the probability λ ij to which endpoint is lifted (see eq. ( 20) is simply proportional to the the distance to the other endpoint, i.e., it is lifted with higher probability to the closer endpoint.
In a composite colloidal problem the ECMC algorithm also faces the problem of lifting between different species.This can be performed exactly according to the same rules as lifting between a single species.In the presence of additional disks, MC displacement of needle endpoints is also lifted to disks if a needle collides with the disk and vice versa as in Fig. 9A.The example of a needle system also shows that collision detection is often the computational bottleneck of ECMC simulations.Here, we use a sophisticated neighbor list system to achieve high simulation speeds also in the nematic phase.Each particle is confined to a container, which triggers an event when the particle leaves it.Then the neighbor list is updated, which ensures that the neighbor lists are always valid.Particles are added to the neighbor lists of the other particle and vice versa when their containers overlap.This way, for different particles different container shapes can be chosen.For the needles, a very narrow rectangle can be used, which limits the computational effort for calculating the distances to the next collision and makes the simulation significantly more efficient in the nematic phase.In particular, the anisotropy of the needles can be assigned particularly well without sacrificing any flexibility for the bookkeeping of the spheres.Furthermore, we optimize the updating of lists by putting them onto a collision grid.
In Fig. 9B, we show a detailed benchmarking for a two-dimensional pure needle system.With increasing EC length , the sampling efficiency is increased, where an EC length comparable to the mean free path can be seen as the limit of a local MC simulation.In Fig. 9B, we relax an isotropic initial configuration towards nematic equilibrium in a quite dense needle system at density ρ n = 10/l 2 0 .The resulting relaxation of the nematic order parameter is exponentially, where we use the decay constant of the exponential as efficiency measure.For EC lengths approaching the system size, no further improvements occur as naively expected.
In three dimensions, hard needles introduce no excluded volume and therefore lack any phase transitions (we have checked with ECMC simulations not shown here that, independent of the density, the nematic order parameter vanishes for hard needles in equilibrium).In three dimensions hard spherocylinders are a natural extension of needles that introduces a non-vanishing excluded volume.As an additional proof of concept of our ECMC implementation with three-particle lifting, we investigate the series of liquid crystal phase transitions by monitoring the pressure as a function of increasing volume fraction, see Fig. 10.We see a sequence of transitions from isotropic to nematic ordering, followed by a transition to smectic order, and ultimately, crystallization in complete agreement with literature results [56].
For the colloidal mixture of hard disks suspended in a nematic host of hard needles, the ECMC simulation reveals a surprising and robust tendency for parallel chaining of disks along the director axis which seems to contradict the chaining in a 45 • angle with respect to the director axis as predicted by quadrupolar elastic interactions in two dimensions.[86] This is a result of a dominant short-range depletion interaction, which strongly favors parallel chaining [63].The effective disk-disk interaction can be obtained as a potential of mean force by simulations of systems containing only two colloidal disks and is in good agreement with analytical calculations, where we add the depletion interaction mediated by needles on short scales and the elastic quadrupolar interaction mediated by a nematic needle LC which has weak planar anchoring at the colloidal disks and exhibits elastic anisotropy.
Here, we demonstrate the validity of this results on a larger scale for systems containing many colloidal disks.Figures 9C and D show the equilibration process of 160 and 640 disks with diameter σ = 1 in needle densities of ρ n = 20/l 2 0 and 10/l 2 0 , respectively, where l 0 is the equilibrium length of a needle.We clearly see that parallel chaining always corresponds to the equilibrium state of the composite system.

Discussion
We presented several applications that demonstrate how ECMC is an effective and fast simulation technique that is particularly suited for dense soft matter systems.Polymer melts, bundles of semiflexible polymers, and the composite system of a liquid crystal colloid are computationally challenging problems, where the ECMC techniques give an improved performance.
We started with an introduction to the algorithm, where we showed how arbitrary interactions between hard spheres can be easily included into the ECMC framework such that a multitude of soft matter systems can be modeled.On this algorithmic side, long-range interactions have not been covered explicitly here but can also be included in an effective manner [39].Future developments will also address variations of the basic EC moves presented here in order to sample even more efficiently.On a generic collision, one interaction rejects the move and lifting occurs which means the current EC direction has a component parallel to the energy gradient of this specific interaction.One can resample the perpendicular parts of the EC direction in order to avoid reshuffling the EC direction from time to time, i.e., a finite EC length.This procedure has been recently developed and is called forward EC [58].In specialized settings, the scheme seems very promising with regard to efficient sampling and should also be tested in the context of soft matter systems in the future.
On the application side in the soft matter field, we presented new results for the initialization and pre-equilibration of polymer melts, where a novel rattling algorithm based on ECs can be employed to efficiently remove overlaps in initial configurations.For large systems of semiflexible polymers with a short-range attraction, we demonstrated that ECMC simulations are capable to follow the formation of large networks of bundles and their subsequent foam-like coarsening behavior, see Fig. 8.For a two-dimensional hard disk and needle model for LC colloids, the structure formation by parallel chaining of disks could be followed for large systems over large time scales as shown in Fig. 9.Here we find that increasing the EC length, i.e., increasing the cluster size in our EC algorithm gives rise to considerable performance gains.
In summary, we addressed compact quasi-zero-dimensional hard sphere particles, one-dimensional hard rods or hard needles, and one-dimensional polymers as constituents in soft matter ECMC simulations.A natural extension for future work is the development of ECMC algorithms for two-dimensional triangulated surfaces with local elastic properties, which are impenetrable for other simulation objects, e.g., hard spheres or polymers.This will allow us to study fluctuating triangulated elastic membranes and, in a compound polymer-membrane system, polymer networks confined to elastic capsules.It was shown that the boundary conditions significantly contribute to the network structure [12].The elastic properties of a triangulated surface can be described by a set of elastic energies such as the TRBS model [11], where elastic constants like Young's modulus and Poisson ratio are adjustable and also bending stiffness can be included [84].
The actual implementation will be a very challenging problem.The interaction of a triangle with a point or a hard sphere is an effective 4-particle interaction.The rejection distance is efficiently calculable by the Möller-Trumbore intersection algorithm.To implement non-intersecting surfaces, interactions between two triangles have to be introduced, which can lead to a 6-particle interaction.Therefore, membrane simulations would also represent a further test of the general framework for N -particle interactions in the EC algorithm [27].The conditional lifting probabilities will look similar to the case of hard needles.
Based on the extensive experiences with EC based simulations presented here, we think this simulation technique is suitable to tackle large (biological) systems.In order to verify and evaluate the EC algorithm for triangulated surfaces, we can investigate the crumpling transition [38].In Ref. 91, the necessary biophysical conditions for the formation of tubular membrane protrusions have been investigated and, it has been suggested that actin filaments polymerizing against a soft membrane are sufficient to form protrusions, which seems a suitable next system to explore the possibilties of EC algorithms for dense soft matter systems.

p
Metr (a → b) = min 1, π b π a = min(1, exp(−∆E a→b )) = exp(−[∆E a→b ] + ) (1) with [x] + ≡ max(0, x) for a move a → b.Together with the trial probability p trial (a → b) that a move a → b is proposed we obtain the transition rate as the product p(a → b) = p trial (a → b)p Metr (a → b).The transition rates p(a → b) give rise to probability currents J a→b = π a p(a → b) from state a to b. Transition rates p(a → b) are given per MC time; likewise probability currents J a→b have units probability per MC time, which must not be identified with an actual physical time.The Metropolis rule is designed to fulfill detailed balance of probability currents between any states a and b, i.e., J a→b = π a p(a → b) = π b p(b → a) = J b→a .

Figure 2 :
Figure2: Scheme of the lifting formalism for hard spheres.The physical states are labeled with letters (a,b,c) and the lifting variable is indicated by color, i.e., the current active particle is indicated by a colored arrow.The EC displacement direction is to the right.Proposed infinitesimal displacements are indicated by opaque spheres.In case of a physical rejection at MC step τ 1 the pivot particle is changed and rejection is avoided.Periodic boundaries ensure ergodicity.One way to read this scheme, is that the lower row is the displacement-reversed version (except trivial translations) of the upper row (obtained be reversing all displacements dw → −dw), resulting in a central symmetry, especially (c,red) = (a,green).

Figure 3 :
Figure 3: Scheme of the lifting formalism for a continuous pair interaction (notation see Fig. 2).Using eq.(5) in the global balance condition (3) for incoming and outgoing currents for state b and dividing by π b yields p lift (red → green) − [dE b→c ] + = p lift (green → red) − [dE b→a ], where dE b→c = −dE b→a ≡ dE because of translational invariance.[dE]+ = 0 leads to p lift (red → green) = 0, yielding p lift (green → red) = [−dE] + and vice versa, which assures maximal global balance.Detailed balance for physical transitions (green) ensures time reversibility, so that a reversed EC would redo the original one, while an EC itself satisfies maximal global balance.For clarity we completed the lower row to show translational symmetry.

Figure 6 :
Figure 6: The background shows a polymer melt of flexible chains, which are presented as volumetric splines.(A)Scheme for the rejection distance for hard sphere chains.If two bonded spheres do not collide, the active sphere can move freely until the minimal distance between both is reached.A usable energy ∆E is drawn (see Fig.4) which determines where the rejection takes place.(B) Scheme of a swap move.In case of a hard sphere collision, the energy difference of each involved springs is calculated when the beads would switch positions.The swap move is accepted with a standard Metropolis filter, otherwise the EC is continued normally.The active particle is marked by a red circle.

Figure 7 :
Figure 7:(A) Wall time per particle to generate rattled configurations of variable polymer number M with a fixed length of N = 200.The wall time scales linearly in the numbers of particles, i.e., curves for different system sizes L × L are identical and diverge when the packing fraction η = M N (π/6)(σ 3 /L 3 ) approaches the random closed packed limit η rcp ≈ 0.64 (indicated as vertical line).Preparation of even larger systems (M × N = 10 6 ) for moderate η (≈ 0.4) take only ≈ 10 3 s.Preparation of pure hard sphere systems (N = 1) are roughly twice as fast.For the following we use M = 1000 chains à N = 100 beads at η ≈ 0.45 (ρ = 0.85/σ 3 , i.e., L ≈ 49σ).All data is generated on a single core of a i7-8650U.(B) The pair correlation g(r) is sensitive to the local structure of the hard spheres and resembles the equilibrium measurement quite well.The warm-up takes only a couple of minutes.(C) The internal bead-to-bead distance R 2 /(n b 2 ) captures how well the polymers are equilibrated and has to be as close to the equilibrium result as possible because this observable equilibrates very slowly.For equilibration a chain needs to move its own size, which is slow due to reptational dynamics.Because of strong fluctuations the warm-up times have not been measured precisely, but are of the order of a couple of hours.(D) The tangent correlation shows the correct non-ideal behaviour of chains in the melt[92,93].The configuration after rattling resembles the equilibrium quite well.The warm-up takes only a couple of minutes.

Figure 8 :
Figure 8: (A) A series of snapshots (after 0, 10 4 , 10 5 , 10 6 sweeps; in a sweep, each sphere is moved once on average); demonstrating coarsening.The snapshots are taken from a simulation of M = 1000 polymers à N = 100 beads in a 400 σ × 400 σ × 10 σ simulation box.The simulation performed with the following parameters: bending constant κ = 30 k B T σ, spring constant k = 100 k B T /σ 2 and rest length b 0 = σ, potential strength g = 0.7 k B T and range d = 0.4 σ.A random color is assigned to each polymer for clarity.(B) Example of a simple graph representation (white balls and bonds) as overlay over a network snapshot generated by an automated procedure to identify vertices (white balls) and edges (white bonds) of the bundle network.The color of cells codes for the number of enclosing edges.All data in Figs.C and D is measured from such graphs.For clarity, we use a comparatively small system.(C) Distribution P (n) of the number of edges n per cell for different simulation times.The curves are obtained by fitting log-normal distributions to the data.(D) Aboav-Weaire law M (n) for different simulation times, which measures the mean edge counts of neighboring cells.Cells with a high number of edges tend to have neighboring cells with fewer edges and vice versa.Error bars represent one standard deviation.The Aboav-Weaire law, the fits of P (n) in C as well as direct computation of the variance µ 2 of the number of edges per cell yield similar values µ 2 ∼ 5.2 for the simulation parameters given in A.

Figure 9 :
Figure9: (A) List of all geometrical cases where rejections/lifting between needles and disks occur (first two depict needle-needle, the third needle-disk, and the last two ones sphere-needle collisions).If a needle is not hit one of its endpoints (first and last picture) conditional lifting probabilities λ ij for a three-particle interaction are applied.(B) Efficiency gains with increasing EC length for a two-dimensional needle system in a volume of L × L at density ρ n = 10/l 2 0 .Starting at one random configuration an ensemble of simulations with different EC lengths are conducted.We show the difference between the nematic order parameter over time and the equilibrium value yielding an exponential.The decay rate of the order parameter is used as efficiency measure, where a simulation with = 5l 0 is 2.8 times more efficient than a simulation with = 0.1, where the latter took even 1.7 times the wall time, yielding a speed-up of ≈ 5 between = 0.1 and = 5 = L/2.(C,D) Evolution of a mixture of hard needles of length 0 = 1 and hard disks with a diameter σ = l 0 indicating an effective interaction between disks, which favors parallel chaining.(C) A system of size 20 0 × 20 0 containing needles at density ρ n = 20/l 2 0 and 160 disks of diameter σ = l 0 , (D) a larger system of 40 0 × 40 0 containing needles at density ρ n = 10/l 2 0 and 640 disks of diameter σ = l 0 .

Figure 10 :
Figure 10: Phase diagram of hard spherocylinders (with diameter D and length L = 5D) in three dimensions.We see four phases (I -isotropic, N -nematic, S -smectic, C -crystalline), in agreement with previous results by McGrother et al. [56] We added six snapshots at different volume fractions η to illustrate the different phases and phase transitions (also indicated by the blue regions).