Efficient Equilibration of Hard Spheres with Newtonian Event Chains
Marco Klement Michael Engel




marco.klement@fau.de
	

michael.engel@fau.de

Institute for Multiscale Simulation, Friedrich-Alexander University Erlangen-Nürnberg, 91058 Erlangen, Germany




An important task in the simulation of hard spheres and other hard particles is structure prediction via equilibration.
Event-driven molecular dynamics is efficient because its Newtonian dynamics equilibrates fluctuations with the speed of sound.
Monte Carlo simulation is efficient if performed with correlated position updates in event chains.
Here, we combine the core concepts of molecular dynamics and event chains into a new algorithm involving Newtonian event chains.
Measurements of the diffusion coefficient, nucleation rate, and melting speed demonstrate that Newtonian event chains outperform other algorithms.
Newtonian event chains scale well to large systems and can be extended to anisotropic hard particles without approximations.


Efficient Equilibration of Hard Spheres with Newtonian Event Chains
    Michael Engel
    December 30, 2023
===================================================================









§ INTRODUCTION


Hard particles receive much attention as model particles in the statistical mechanical theory of colloidal self-assembly.<cit.>
Research is often tasked with finding the stable phase, which can be performed in simulation via equilibration from a disordered fluid.<cit.>
During equilibration, the system follows a trajectory in configuration space from a basin of low probability corresponding to a phase with higher free energy into a basin of high probability corresponding to a phase with lower free energy.
Configuration space can either be sampled deterministically with molecular dynamics or stochastically with Monte Carlo.
The equilibrium phase, reached after a sufficiently long simulation, is independent of the algorithm used.
However, the time to reach equilibrium depends crucially on details of the simulation trajectory.
Many ideas have been proposed to improve efficiency by decreasing the length of the trajectory connecting a starting configuration to the basin of lowest free energy.

A classic way to simulate particles is to solve Newton's equations of motion.
This is the approach of molecular dynamics.
In the case of hard particles, molecular dynamics involves the prediction of collision events.<cit.>
Event-driven molecular dynamics (EDMD) moves all particles fully collectively and thus equilibrates density fluctuations with the speed of sound.
EDMD is generally considered to be highly efficient.<cit.>
In contrast, local Monte Carlo (LMC) moves only one randomly chosen particle at a time.<cit.>
LMC neglects particle momentum, which means fluctuations propagate and relax much slower resulting in poor performance in comparison to EDMD.
Advantages of LMC are its simple implementation, good scaling to large systems, and easy generalization to anisotropic particles.
Equilibration via Monte Carlo that hopes to rival EDMD in simulation efficiency requires collective moves. <cit.>
A promising approach is updating particles in chains.<cit.>
Chain moves can be performed rejection-free in the form of event chain Monte Carlo.<cit.>
In an event chain, a randomly selected particle is displaced until it collides with another particle.
Next, the collision partner is displaced until it collides with yet another particle and the process iterated.
The chain of collision events terminates once a preselected chain length has been reached.<cit.>

Existing collective Monte Carlo algorithms neglect a core concept of molecular dynamics, the momentum of the particles in the system.
As a result, fluctuations decay slowly, which lowers efficiency.
Here, we propose combining the main advantage of EDMD, Newtonian dynamics, with the main advantage of event chain Monte Carlo, collective moves in chains.
We call the new algorithm Newtonian event chains (NEC).
We evaluate equilibration efficiency of NEC and compare it to existing algorithms.
While EDMD is most efficient per particle displacement, NEC comes close and is more efficient when analyzed in terms of CPU time.
The only free parameter is the duration of each event chain.
In the limit of short chains, NEC reduces to a variant of molecular dynamics with random permutation of particle updates.







§ ALGORITHMS AND METHODS




 §.§ Existing hard sphere simulation algorithms




Local Monte Carlo (LMC).
A local Monte Carlo trial move consists of selecting a random particle and displacing it to a new position randomly chosen within a sphere of radius L_step.
The trial move is accepted with probability given by the Metropolis criterion.
For hard particles the Metropolis criterion states that moves that generate one or more overlaps are always rejected (<ref>).
Trial moves that do not generate an overlap are always accepted.
LMC is simple to implement but not efficient.

Event-driven molecular dynamics (EDMD).
In molecular dynamics each particle has position and velocity, which are updated by solving Newton's equations of motion.
Molecular dynamics for hard particles can only be performed event-driven.
This means the times of all collision events in the system have to be computed.
Collisions are sorted and evaluated in temporal order, which is slow and difficult to scale to large systems.
The pair of particles involved in the current collision event is moved up to their collision positions.
Next, new velocities are computed from the equations for elastic collisions.
EDMD is difficult to parallelize and scales poorly to large systems.
It also cannot be applied to anisotropic particles without an iterative solver, which lowers precision and performance, because calculating collision times for such particles involves nonlinear equations that usually have no analytic solution.<cit.>

Event chain Monte Carlo.
An event chain starts with a randomly selected particle i with radius r_i, located at position x⃗_i and a normalized displacement vector d⃗.
Particle i is displaced along d⃗ until it collides with a first other particle j, its collision partner.
The position of particle i is then updated to be the contact point with the collision partner, given by

    x⃗_i' = x⃗_i + d⃗(d⃗·x⃗_ij - √( (r_i+r_j)^2 - x_ij^2 + (d⃗·x⃗_ij)^2 ))

with x⃗_ij=x⃗_j - x⃗_i.
Next, it has to be decided with which particle and in which direction the event chain is continued.
We denote the new particle that is chosen to continue the event chain as i' and the new displacement direction with d⃗'.

Several variants for selecting i' and d⃗' have been proposed in the literature.<cit.>
Straight event chains (SEC) is the simplest choice.
The chain is continued with the collision partner and in the same direction,

    SEC:   i'=j,  d⃗' = d⃗.

In contrast, in reflected event chains (REC) the direction is reflected,

    REC:   i'=j,  d⃗' = d⃗_r = 2 x⃗_ijd⃗·x⃗_ij/x_ij^2 - d⃗.

Detailed balance is obeyed in REC because the equation for d⃗' is symmetric.
The event chain ends when it reaches a preselected length.
Once the sum of all displacements distances x⃗_i'-x⃗_i adds to above the chain length L_chain, the last particle is displaced only by the remaining distance.

We modify the event chain variants SEC and REC to mimic trajectories generated by Newtonian dynamics more closely.
For this purpose an event chain is sought that follows the path of largest momentum transfer.
This is achieved by distinguishing two cases.
Collisions with obtuse collision angles, d⃗·d⃗_r>0, are handled as before in REC.
Collisions with acute angles, d⃗·d⃗_r<0, are now handled as total reflections.
Total reflections continue the chain with the original particle and mirror the displacement vector.
We call the resulting algorithm obtuse reflected event chains (OREC).
It obeys

    OREC:  
    i'=j,  d⃗' = d⃗_r     if d⃗·d⃗_r>0,
     
    
    i'=i,  d⃗' = -d⃗_r    if d⃗·d⃗_r<0.

Again, the equation for d⃗' is symmetric, which ensures detailed balance.
Other variants of event chains handling collisions with acute and obtuse angles in different ways are possible.
We tested other variants but found that OREC is the most efficient variant among those.
For this reason we will focus on the well-known event chain variant SEC and its slight modification OREC in the following.





 §.§ Newtonian event chains


Newtonian dynamics is beneficial for efficient equilibration.
We therefore wish to augment event chains by assigning a velocity vector v⃗_i to each particle in the system.
Velocities are initialized according to the Maxwell-Boltzmann distribution.
As before, we start an event chain by selecting a particle i at random.
This particle is then displaced along its normalized velocity vector d⃗=v⃗_i/v_i up to the first collision with collision partner j.
The collision is handled simply as an elastic collision of two spheres with the same mass taking into account the velocity vectors of both particles.
From the collision equations we obtain the new velocities v⃗'_i and v⃗'_j.
Next, the chain proceeds with the collision partner as in SEC but now in direction of the updated velocity.

The NEC algorithm requires only a few more calculations than the SEC algorithm but otherwise no additional effort.
It reads

    NEC:     v⃗_i' = v⃗_i + x⃗_ijv⃗_ij·x⃗_ij/x_ij^2,
       v⃗_j' = v⃗_j - x⃗_ijv⃗_ij·x⃗_ij/x_ij^2,
        i'=j,  d⃗' = v⃗'_j/v'_j

with v⃗_ij=v⃗_j - v⃗_i.
The event chain ends when it reaches a preselected duration.
Once the sum of all displacement times x⃗_i'-x⃗_i / v_i adds to above the chain duration T_chain, the last particle is displaced only by the remaining time.

It is important to preselect a duration for event chains in NEC as parameter instead of preselect a length.
To see this, assume that all particle velocities are initially distributed according to Maxwell-Boltzmann statistics.
Because we perform collisions using Newtonian dynamics, the velocity distribution after a fixed elapsed time or, alternatively, a preselected event chain duration is not affected.
As demonstrated in Ref., the invariance of the Maxwell-Boltzmann velocity distribution guarantees local detailed balance, which in turn is sufficient to ensure a statistically valid Monte Carlo simulation.
Hence, NEC with preselected chain duration is a statistically correct algorithm.

In contrast, let us assume a preselected length was chosen for event chains in NEC.
Internal time of event chains that contain many slow particles then advances faster than internal time of event chains that contain many fast particles.
As a result, such event chains are slightly more likely to terminate with a fast particle than expected from Maxwell-Boltzmann.
This means the velocity distribution skews towards higher velocities.
Clearly, NEC with a preselected chain length violates the local balance condition.
We directly observe this effect.
Reference values for pressure and local order are reproduced correctly when the chain duration is fixed (<ref>) but not if the chain length is fixed.





 §.§ Measured quantities
 

There are many ways to evaluate equilibration efficiency, including autocorrelation times of thermodynamic quantities,<cit.> positional auto correlation such as the diffusion coefficient,<cit.> and evolution of specific quantities.<cit.>
Here we focus on measuring the diffusion coefficient, nucleation rate, and melting speed.
We find that these three quantities are roughly equivalent for evaluating the efficiency of hard sphere simulation algorithms.
Because Monte Carlo does not have an internal time, we analyze the advancement of the simulation by the number of displacements N_disp.
This number is the number of trial moves for LMC, the number of particle translations in all event chains, and the number of time steps multiplied with the number of particles N in the system for EDMD.
We interpret N_disp as a measure for the length of the trajectory in configuration space.

Diffusion coefficient.
The diffusion coefficient as a function of the number of displacements is calculated from the mean square displacement by the relationship

    D_disp = lim_N_disp→∞∑_i=1^N(Δx⃗_i - ⟨Δx⃗⟩)^2/6   N_disp,

where Δx⃗_i is the total displacement vector of particle i after N_disp particle displacements in the system and we subtract the mean displacement vector ⟨Δx⃗⟩ obtained by averaging over all particles.
Not all displacements require the same computational cost.
We therefore also determine the diffusion coefficient as a function of CPU time t_cpu,

    D_cpu = lim_t_cpu→∞∑_i=1^N (Δx⃗_i - ⟨Δx⃗⟩)^2/6 N  t_cpu.

While D_disp is an absolute measure of the effectiveness of displacements for generating diffusion, D_cpu is affected by practical factors such as the implementation of the simulation algorithm and the computer hardware.

Nucleation rate.
We model nucleation as a Bernoulli process.
The probability to find a new nucleus within time τ is p.
There are N_t=t/τ trials to create a nucleus within the simulation time t.
On average we find N_n=N_t p nuclei with a variance of Var(N_n) = N_t p (1-p).
For sufficiently small τ we can approximate 1-p ≈ 1, and thus Var(N_n) = N_t p = N_n.
The nucleation rate R = N_n / t has standard deviation σ = √(N_n)/t.
Because Monte Carlo does not have an internal time, we use the number of displacements N_disp or CPU time t_cpu instead.
Nucleation rate is measured by performing several independent simulations starting in a fluid phase.
When we detect a critical nucleus, we increment N_n and reset the simulation to the fluid phase.
We stop once we find the twentieth nucleus and record t_cpu and N_disp at the end of the sequence.

Structure analysis.
We quantify the ordering process during the fluid-to-solid transformation with the global bond orientational order parameter Q_6 and the local bond orientational order parameter q_6.<cit.>
The SANN algorithm<cit.> is used to identify nearest neighbors.
Nuclei are detected with the common cluster algorithm.<cit.>
We determine melting speed by measuring the necessary CPU time to reach Q_6=0.1.





 §.§ Implementation and setup




We use the hard particle Monte Carlo (HPMC) package<cit.> of HOOMD-blue<cit.> without changes as implementation for LMC.
All variants of event chain algorithms are implemented without parallelization into a private fork of HPMC.
Implementations employ an AABB tree as neighbor list.
EDMD is our own implementation using a priority queue with constant complexity as described in Refs. .
All simulations are performed on Intel Xeon E3-1240 v5 CPUs with 3.5 GHz.

We compare the computational costs for displacing particles in <ref>.
LMC is the simplest algorithm and displaces particles the fastest.
The three event chain algorithms displace particles roughly in the same time, which is slower by a factor of 2.5 to 3.5 compared to LMC, depending on whether the system size is N=2^14=16 384 or N=2^17=131 072.
EDMD is the algorithm that requires most computational effort.
It displaces particles a factor of between 5 and 14 slower than LMC.
The performance of our implementations is comparable to that of Ref. , which has been obtained, however, using Intel Xeon E5-2680 CPUs with 2.8 GHz, i.e. a similar microarchitecture with 20% slower clock rate than the CPUs used in our tests.
If computational cost is corrected via division by CPU clock rate, which certainly is a strong simplification but nevertheless provides a rough comparison, then our implementation is 6% faster for LMC, 37% slower for SEC, and 44% faster for EDMD.
An explanation for the slower performance of SEC is discussed at the end of <ref>.

Systems of N spheres with diameter d are initialized in either a random fluid starting configuration or a face-centered cubic (FCC) solid.
Simulations are performed in the isochoric ensemble.
Pressure is calculated from the virial expression as in previous works.<cit.>
Except for the parameter study optimizing Monte Carlo parameters in the next section, the step size of LMC is tuned for maximal speed to about 15% to 20% trial move acceptance probability.
Chains with chain length of approximately one particle diameter or equivalent chain duration are used for all event chain variants except when noted otherwise.







 §.§ Validation

To validate our implementation of the three event chain variants, we measure dimensionless pressure and the local bond orientational order parameter in a fluid of N=131 072 particles at 49% volume fraction and compare both to the reference values<cit.> β P d^3 = 11.3894(1) and q_6=0.38678.
All our algorithms give results in agreement to these numbers within their respective standard errors (<ref>).
Pressure averages to β P d^3 = 11.3907(14) for SEC, β P d^3 = 11.3906(13) for OREC, and β P d^3 = 11.3886(18) for NEC.
The equilibrium value for the local bond orientational order parameter is q_6 = 0.38678(3) for SEC, q_6 = 0.38679(3) for OREC, and q_6 = 0.38680(4) for NEC.
Our standard errors in these measurements are obtained by block averaging.
Pressure is averaged over four independent runs.







§ EFFICIENCY ANALYSIS




 §.§ Parameter dependence
		



The efficiency of Monte Carlo algorithms is affected by the choice of the Monte Carlo parameters step size, chain length, and chain duration.
<ref> shows the two efficiency measures diffusion coefficient and nucleation rate recorded as a function of these parameters.
A narrow peak for LMC, half an order of magnitude wide, is found for both efficiency measures near L_step=3L_MFP.
Similar behavior was observed for hard discs in two dimensions.<cit.>
Maximum efficiency occurs near trial move acceptance probability α=20%.
As is often the case, rejected moves are faster to execute than accepted moves due to early returns from the overlap search and the fact that rejected moves do not require an update of the particle position.
The data for SEC and OREC initially follows the curve for LMC and then increases further beyond the LMC peak.
A broad plateau appears for chain length in the range 1<L_chain/d<100 in the case of SEC and OREC and for equivalent chain duration in the case of NEC.
OREC performs slightly better and NEC much better than LMC for all chain lengths or chain durations.

We would expect the diffusion coefficient of SEC to match the behavior reported in Ref.  and decrease after the point where a significant part of all particles in the system are translated in the same direction.
In contrast, OREC should not drop even for long chains because the displacement direction keeps getting updated.
We observe in <ref> that SEC and OREC both decrease in lockstep towards high L_chain.
This means our simulations did not yet reach the point where a sufficient number of particles in the system are moved.
The decrease of D_cpu must be explained solely due to details of our implementation.
In fact, HOOMD keeps track of particles neighborhoods with an AABB tree data structure. <cit.>
In the current implementation, the AABB tree is only updated but not continuously rebalanced during an event chain.<cit.>
This means the number of neighbor particle candidates gradually increases as the event chain progresses, which slows down simulations of long event chains.
Efficiency is affected less by overly long chains in NEC.
We do not attempt to improve the use of the AABB tree and select Monte Carlo parameters from the beginning of the plateau.





 §.§ Volume fraction dependence




Event chain Monte Carlo improves simulation efficiency near the fluid-to-solid transition.<cit.>
Indeed, we observe a performance increase of SEC compared to LMC by a factor of 5 in terms of the number of displacements and 2 in terms of CPU time (<ref>).
The advantage of SEC over LMC keeps growing towards higher volume fraction because particle dynamics becomes more collective.
Still, we do not see quite such a large increase near the critical point as reported in two dimensions.<cit.>
Apparently, collective motion is less important in three dimensions than in two dimensions.

Across all volume fractions, OREC outperforms SEC by a factor of 1.6 but still never reaches the speed of EDMD.
The performance gap between EDMD and all purely stochastic Monte Carlo algorithms grows with volume fraction demonstrating that the Newtonian dynamics of EDMD is superior to stochastic dynamics.
Indeed, the inclusion of Newtonian dynamics into event chain Monte Carlo in the form of NEC outperforms all other algorithms for all densities.





 §.§ System size dependence




The efficiency of algorithms is compared across system sizes in <ref>.
We choose LMC with step size optimized for maximum diffusion as reference.
The trial move acceptance probability of optimized (opt.) LMC is 15% to 20%.
In comparison, LMC with acceptance probability 50% is slower by at least a factor of 2.
As expected, D_disp is essentially independent of system size.
This makes sense because it only depends on local structure.
In contrast, D_cpu is affected by system size, most strongly for EDMD.

Because HOOMD-blue and HPMC are highly optimized codes, LMC and event chain Monte Carlo scale well to large numbers of particles.
Our own implementation of EDMD is most efficient for small system size.
Efficiency quickly decreases once system size reaches the CPU cache limit.
We believe our implementation of EDMD can be improved with more effort.
Still, it remains inherently difficult to achieve good scaling to large system sizes for EDMD because collisions can happen anywhere at all times, which makes it practically impossible to utilize CPU cache well.
In contrast, successive collisions in event chains occur in close proximity of another.
For this reason the performance advantage of NEC increases with system size.





 §.§ Melting the FCC crystal
 



We investigate equilibration by melting in <ref>.
The hard sphere system is initialized in a FCC crystal structure, equilibrated at 55% volume fraction, and then expanded to 49%.
All algorithms show a similar melting behavior as quantified by the decrease of the global bond orientational order parameter Q_6 (<ref>).
NEC melts within the shortest CPU time, next EDMD, then OREC, and finally SEC.
This sequence is in agreement with measurements of the diffusion coefficient (<ref>).
We observe in <ref> that melting speed is approximately inversely proportional to the diffusion coefficient, which suggests that melting is dominated by diffusion processes.

We remark on an unexpected behavior that we observed in early stages of our tests.
<cit.> propose restricting displacements in SEC in the +x, +y, and +z directions only.
If we implement such a restriction and align the FCC crystal along the coordinate axes then an artificial speed-up of SEC is observed.
Symmetry axis-aligned event chains melt the crystal faster than event chains with displacement directions not aligned with crystallographic axes and faster than event chains with randomly chosen displacement directions.
To avoid this unphysical dependence of melting on the orientation of the crystal and for better comparison of event chain variants we always choose displacement directions fully randomly in our implementation of SEC.
The use of random displacement directions slightly increases the computational cost for translating particles, which is observed in <ref>.









§ CONCLUDING REMARKS




A central result of this paper is the finding that equilibration trajectories are shorter if particle dynamics more closely mimics Newtonian dynamics.
By bringing the geometry of event chains closer to trajectories generated by molecular dynamics, first in the variant OREC and then in the form of NEC, simulations are sped up significantly.
We believe conservation of momentum in NEC leads to a memory effect that synchronizes the motion of particles in close proximity and allows the system to equilibrate more collectively.

We can understand the effectiveness of particle displacements for advancing the system in configuration space by analyzing mean square displacement in <ref>.
For short times, i.e. for a small number of displacements, N_disp≲ 1, ballistic motion in EDMD is responsible for slope 2 scaling while Monte Carlo algorithms have slope ≤ 1.
All algorithms show signs of a slowdown at intermediate times, 10<N_disp<100, indicative of caging.
EDMD and NEC gradually increase their advantage over other algorithms and escape caging after the fewest number of displacements.
For long times, N_disp≫ 1000, all algorithms behave indistinguishably and evolve diffusively.
Because displacements in NEC are faster to compute, NEC outperforms EDMD when analyzed in terms of CPU time.

In summary, NEC combines aspects of Newtonian dynamics with event chain Monte Carlo.
NEC is an improvement over past event chain variants and requires little additional implementation effort.
It is currently the algorithms with best CPU performance available for hard spheres.
As all event chain algorithms, NEC scales better than EDMD to large systems.
Future work can generalize NEC to anisotropic hard particles and implement parallelization.






We thank Praveen Bommineni and Sebastian Kapfer for helpful discussions.
This work has been funded by Cluster of Excellence Engineering of Advanced Materials grant EXC 315/2 of the German Research Foundation (DFG).
Support by the Central Institute for Scientific Computing (ZISC), the Interdisciplinary Center for Functional Particle Systems (FPS), and computational resources and support provided by the Erlangen Regional Computing Center (RRZE) are gratefully acknowledged.











35
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Frenkel(2015)]Frenkel2015
  author author D. Frenkel, 10.1038/nmat4178 journal
  journal Nat. Mater. volume 14,  pages 9 (year 2015)NoStop
[Manoharan(2015)]Manoharan2015
  author author V. N. Manoharan, 10.1126/science.1253751 journal journal Science volume
  349, pages 1253751 (year 2015)NoStop
[Boles, Engel, and Talapin(2016)]Boles2016
  author author M. A. Boles, author M. Engel,  and  author D. V. Talapin, 10.1021/acs.chemrev.6b00196 journal journal Chem. Rev. volume 116, pages
  11220 (year 2016)NoStop
[Agarwal and Escobedo(2011)]Agarwal2011
  author author U. Agarwal and author F. A. Escobedo, 10.1038/nmat2959 journal
  journal Nat. Mater. volume 10,  pages 230 (year 2011)NoStop
[Damasceno, Engel, and  Glotzer(2012)]Damasceno2012a
  author author P. F. Damasceno, author M. Engel,   and author S. C. Glotzer,  10.1126/science.1220869 journal journal Science volume 337, pages
  453 (year 2012)NoStop
[Ni et al.(2012)Ni,
  Gantapara, de Graaf, van
  Roij, and Dijkstra]Ni2012
  author author R. Ni, author A. P. Gantapara,
  author J. de Graaf, author R. van Roij,  and author M. Dijkstra, 10.1039/c2sm25813g journal journal Soft Matter  volume 8, pages 8826 (year
  2012)NoStop
[Alder and Wainwright(1959)]Alder1959
  author author B. J. Alder and author T. E. Wainwright, 10.1063/1.1730376 journal
  journal J. Chem. Phys. volume 31,  pages 459 (year 1959)NoStop
[Isobe(1999)]Isobe1999
  author author M. Isobe, 10.1142/S0129183199001042 journal
  journal Int. J. Mod. Phys. C volume
  10, pages 1281 (year 1999)NoStop
[Donev, Torquato, and Stillinger(2005)]Donev2005a
  author author A. Donev, author S. Torquato,   and author F. H. Stillinger,  10.1016/j.jcp.2004.08.014 journal journal J. Comput. Phys. volume 202, pages 737 (year 2005)NoStop
[Engel et al.(2013)Engel,
  Anderson, Glotzer, Isobe,
  Bernard, and Krauth]Engel2013
  author author M. Engel, author J. A. Anderson,
  author S. C. Glotzer, author M. Isobe, author
  E. P. Bernard,  and author
  W. Krauth, 10.1103/PhysRevE.87.042134 journal journal Phys.
  Rev. E volume 87, pages 042134
  (year 2013)NoStop
[Isobe(2016)]Isobe2016
  author author M. Isobe, 10.1080/08927022.2016.1139106 journal journal Mol. Simul. volume
  42, pages 1317 (year 2016)NoStop
[Metropolis et al.(1953)Metropolis, Rosenbluth, Rosenbluth,
  Teller, and Teller]Metropolis1953
  author author N. Metropolis, author A. W. Rosenbluth, author M. N. Rosenbluth, author A. H. Teller,  and author E. Teller, 10.1063/1.1699114 journal
  journal J. Chem. Phys. volume 21,  pages 1087 (year 1953)NoStop
[Dress and Krauth(1995)]Dress1995
  author author C. Dress and author W. Krauth,  10.1088/0305-4470/28/23/001 journal journal J. Phys. A. Math. Gen. volume 28,  pages L597 (year 1995)NoStop
[Liu and Luijten(2004)]Liu2004
  author author J. Liu and author E. Luijten,  10.1103/PhysRevLett.92.035504 journal
  journal Phys. Rev. Lett. volume 92,  pages 035504 (year 2004)NoStop
[Whitelam and Geissler(2007)]Whitelam2007
  author author S. Whitelam and author P. L. Geissler, 10.1063/1.2790421 journal
  journal J. Chem. Phys. volume 127,  pages 154101 (year 2007)NoStop
[Jaster(1999)]Jaster1999
  author author A. Jaster, 10.1016/S0378-4371(98)00337-9 journal journal Physica A volume
  264, pages 134 (year 1999)NoStop
[Bernard, Krauth, and Wilson(2009)]Bernard2009
  author author E. P. Bernard, author W. Krauth,   and author D. B. Wilson,  10.1103/PhysRevE.80.056704 journal journal Phys. Rev. E volume 80, pages 056704 (year 2009)NoStop
[Peters and de With(2012)]Peters2012
  author author E. A. J. F.  Peters and author G. de With, 10.1103/PhysRevE.85.026703
  journal journal Phys. Rev. E volume 85, pages 026703 (year
  2012)NoStop
[Michel, Kapfer, and Krauth(2014)]Michel2014
  author author M. Michel, author S. C. Kapfer,
   and author W. Krauth, 10.1063/1.4863991 journal journal J.
  Chem. Phys. volume 140, pages 054116
  (year 2014)NoStop
[Kampmann, Boltz, and Kierfeld(2015)]Kampmann2015
  author author T. A. Kampmann, author H.-H. Boltz,
   and author J. Kierfeld,  10.1016/j.jcp.2014.10.059 journal journal J. Comput. Phys. volume 281, pages 864 (year 2015)NoStop
[Isobe and Krauth(2015)]Isobe2015
  author author M. Isobe and author W. Krauth,  10.1063/1.4929529 journal journal J. Chem. Phys. volume 143, pages 084509 (year 2015)NoStop
[Harland et al.(2017)Harland, Michel, Kampmann, and  Kierfeld]Harland2017
  author author J. Harland, author M. Michel,
  author T. A. Kampmann,  and  author J. Kierfeld, 10.1209/0295-5075/117/30001 journal journal Europhys. Lett. volume 117, pages 30001 (year 2017)NoStop
[Weigel(2018)]Weigel2018
  author author R. F. B.  Weigel, title "Equilibration of
  Orientational Order in Hard Disks via Arcuate Event-Chain Monte Carlo",  @noop  type Master's thesis, school
  Friedrich-Alexander-Universität Erlangen-Nürnberg (year 2018)NoStop
[Manousiouthakis and Deem(1999)]Manousiouthakis1999
  author author V. I. Manousiouthakis and author M. W.  Deem, 10.1063/1.477973 journal journal J. Chem. Phys. volume
  110, pages 2753 (year 1999)NoStop
[Swendsen and Wang(1987)]swendsen1987nonuniversal
  author author R. H. Swendsen and author J.-S. Wang, 10.1103/PhysRevLett.58.86 journal
  journal Phys. Rev. Lett. volume 58,  pages 86 (year 1987)NoStop
[R?žička and Allen(2014)]Ruzicka2014
  author author Š. Růžička and author M. P.  Allen, 10.1103/PhysRevE.90.033302
  journal journal Phys. Rev. E volume 90, pages 033302 (year
  2014)NoStop
[Steinhardt, Nelson, and  Ronchetti(1983)]Steinhardt1983
  author author P. J. Steinhardt, author D. R. Nelson,  and author M. Ronchetti, 10.1103/PhysRevB.28.784 journal journal Phys. Rev. B volume
  28, pages 784 (year 1983)NoStop
[van Meel et al.(2012)van
  Meel, Filion, Valeriani, and Frenkel]Meel2012
  author author J. A. van Meel, author L. Filion,
  author C. Valeriani,  and  author D. Frenkel, 10.1063/1.4729313 journal journal J.
  Chem. Phys. volume 136, pages 234107
  (year 2012)NoStop
[Filion et al.(2010)Filion,
  Hermes, Ni, and Dijkstra]Filion2010a
  author author L. Filion, author M. Hermes,
  author R. Ni,  and author M. Dijkstra, 10.1063/1.3506838 journal journal J. Chem.
  Phys. volume 133, pages 244115
  (year 2010)NoStop
[Anderson, Eric Irrgang, and  Glotzer(2016)]Anderson2015a
  author author J. A. Anderson, author M. Eric
  Irrgang,  and author S. C. Glotzer, 10.1016/j.cpc.2016.02.024 journal journal Comput. Phys. Commun. volume 204, pages 21 (year 2016)NoStop
[Anderson, Lorenz, and  Travesset(2008)]Anderson2008a
  author author J. A. Anderson, author C. D. Lorenz,  and author A. Travesset, 10.1016/j.jcp.2008.01.047 journal journal J. Comput. Phys. volume 227, pages 5342 (year
  2008)NoStop
[Glaser et al.(2015)Glaser,
  Nguyen, Anderson, Lui,
  Spiga, Millan, Morse, and  Glotzer]Glaser2015
  author author J. Glaser, author T. D. Nguyen,
  author J. A. Anderson, author P. Lui, author
  F. Spiga, author J. A. Millan, author D. C.  Morse,  and author S. C.  Glotzer, 10.1016/j.cpc.2015.02.028
  journal journal Comput. Phys. Commun. volume 192, pages 97 (year
  2015)NoStop
[Wang et al.(2018)Wang,
  Mbah, Przybilla, Apeleo
  Zubiri, Spiecker, Engel, and  Vogel]Wang2018
  author author J. Wang, author C. F. Mbah,
  author T. Przybilla, author B. Apeleo Zubiri, author E. Spiecker, author
  M. Engel,  and author
  N. Vogel, 10.1038/s41467-018-07600-4 journal journal Nat.
  Commun. volume 9, pages 5259
  (year 2018)NoStop
[Bommineni et al.(2018)Bommineni, Varela-Rosales, Klement, and  Engel]Bommineni2018
  author author P. K. Bommineni, author N. R. Varela-Rosales, author M. Klement,  and author M. Engel, http://arxiv.org/abs/1811.00061   (year
  2018), http://arxiv.org/abs/1811.00061 arXiv:1811.00061
  NoStop
[Howard et al.(2016)Howard,
  Anderson, Nikoubashman, Glotzer, and Panagiotopoulos]Howard2016
  author author M. P. Howard, author J. A. Anderson, author A. Nikoubashman, author S. C. Glotzer,  and author A. Z. Panagiotopoulos, 10.1016/j.cpc.2016.02.003 journal journal Comput. Phys. Commun. volume 203, pages 45 (year 2016)NoStop


