Event-chain Monte Carlo Algorithm Research Articles

Direction-sweep Markov chains

We discuss a non-reversible, lifted Markov-chain Monte Carlo (MCMC) algorithm for particle systems in which the direction of proposed displacements is changed deterministically. This algorithm sweeps through directions analogously to the popular MCMC sweep methods for particle or spin indices. Direction-sweep MCMC can be applied to a wide range of reversible or non-reversible Markov chains, such as the Metropolis algorithm or the event-chain Monte Carlo algorithm. For a single two-dimensional tethered hard-disk dipole, we consider direction-sweep MCMC in the limit where restricted equilibrium is reached among the accessible configurations for a fixed direction before incrementing it. We show rigorously that direction-sweep MCMC leaves the stationary probability distribution unchanged and that it profoundly modifies the Markov-chain trajectory. Long excursions, with persistent rotation in one direction, alternate with long sequences of rapid zigzags resulting in persistent rotation in the opposite direction in the limit of small direction increments. The mapping to a Langevin equation then yields the exact scaling of excursions while the zigzags are described through a non-linear differential equation that is solved exactly. We show that the direction-sweep algorithm can have shorter mixing times than the algorithms with random updates of directions. We point out possible applications of direction-sweep MCMC in polymer physics and in molecular simulation.

Multithreaded event-chain Monte Carlo with local times

We present a multithreaded event-chain Monte Carlo algorithm (ECMC) for hard spheres. Threads synchronize at infrequent breakpoints and otherwise scan for local horizon violations. Using a mapping onto absorbing Markov chains, we rigorously prove the correctness of a sequential-consistency implementation for small test suites. On x86 and ARM processors, a C++ (OpenMP) implementation that uses compare-and-swap primitives for data access achieves considerable speed-up with respect to single-threaded code. The generalized birthday problem suggests that for the number of threads scaling as the square root of the number of spheres, the horizon-violation probability remains small for a fixed simulation time. We provide C++ and Python open-source code that reproduces all our results. Program summaryProgram title: ParaSpheres.CPC Library link to program files:https://doi.org/10.17632/c3rjk5k3z9.1Developer’s repository link:https://github.com/jellyfysh/ParaSpheresLicensing provisions: GNU GPLv3.Programming languages: Python 3, C++, Fortran90.Nature of problem: Multithreaded event-chain Monte Carlo for hard spheres.Solution method: Event-driven irreversible Markov-chain Monte Carlo algorithm using local times.Additional comments: The collection of programs is complete with shell scripts that allow one to reproduce all data, and all the figures of the paper. Change of density and system size is straightforward. The manuscript is accompanied by a frozen copy of the GitHub repository that is made publicly available on GitHub (repository https://github.com/jellyfysh/ParaSpheres, commit hash e2aa5b9727fb080ebe65581586c0f6133efa495d).

Event-chain Monte Carlo with factor fields.

We study the dynamics of one-dimensional (1D) interacting particles simulated with the event-chain Monte Carlo algorithm (ECMC). We argue that previous versions of the algorithm suffer from a mismatch in the factor potential between different particle pairs (factors) and show that in 1D models, this mismatch is overcome by factor fields. ECMC with factor fields is motivated, in 1D, for the harmonic model, and validated for the Lennard-Jones model as well as for hard spheres. In 1D particle systems with short-range interactions, autocorrelation times generally scale with the second power of the system size for reversible Monte Carlo dynamics, and with its first power for regular ECMC and for molecular dynamics. We show, using simulations, that the autocorrelation time grows only with the square root of the system size for ECMC with factor fields. Mixing times, which bound the time to reach equilibrium from an arbitrary initial configuration, grow with the first power of the system size.

Mixing and perfect sampling in one-dimensional particle systems

We study the approach to equilibrium of the event-chain Monte Carlo (ECMC) algorithm for the one-dimensional hard-sphere model. Using the connection to the coupon-collector problem, we prove that a specific version of this local irreversible Markov chain realizes perfect sampling in single steps, whereas the reversible local Metropolis algorithm requires single steps for mixing. This confirms a special case of an earlier conjecture about scaling of mixing times of ECMC and of the lifted forward Metropolis algorithm, its discretized variant. We also prove that sequential ECMC (with swaps) realizes perfect sampling in single events. Numerical simulations indicate a cross-over towards mixing for the sequential forward swap Metropolis algorithm, that we introduce here. We point out open mathematical questions and possible applications of our findings to higher-dimensional models.

All-atom computations with irreversible Markov chains.

We apply the irreversible event-chain Monte Carlo (ECMC) algorithm to the simulation of dense all-atom systems with long-range Coulomb interactions. ECMC is event-driven and exactly samples the Boltzmann distribution. It neither uses time-step approximations nor spatial cutoffs on the range of the interaction potentials. Most importantly, it need not evaluate the total Coulomb potential and thus circumvents the major computational bottleneck of traditional approaches. It only requires the derivatives of the two-particle Coulomb potential, for which we discuss mutually consistent choices. ECMC breaks up the total interaction potential into factors. For particle systems made up of neutral dipolar molecules, we demonstrate the superior performance of dipole-dipole factors that do not decompose the Coulomb potential beyond the two-molecule level. We demonstrate that these long-range factors can nevertheless lead to local lifting schemes, where subsequently moved particles are mostly close to each other. For the simple point-charge water model with flexible molecules (SPC/Fw), which combines the long-ranged intermolecular Coulomb potential with hydrogen-oxygen bond-length vibrations, a flexible hydrogen-oxygen-hydrogen bond angle, and Lennard-Jones oxygen-oxygen potentials, we break up the potential into factors containing between two and six particles. For this all-atom liquid-water model, we demonstrate that the computational complexity of ECMC scales very well with the system size. This is achieved in a pure particle-particle framework, without the interpolating mesh required for the efficient implementation of other modern Coulomb algorithms. Finally, we discuss prospects and challenges for ECMC and outline several future applications.

Irreversible Markov chains in spin models: Topological excitations

We analyze the convergence of the irreversible event-chain Monte Carlo algorithm for continuous spin models in the presence of topological excitations. In the two-dimensional XY model, we show that the local nature of the Markov-chain dynamics leads to slow decay of vortex-antivortex correlations while spin waves decorrelate very quickly. Using a Fréchet description of the maximum vortex-antivortex distance, we quantify the contributions of topological excitations to the equilibrium correlations, and show that they vary from a dynamical critical exponent at the critical temperature to in the limit of zero temperature. We confirm the event-chain algorithm's fast relaxation (corresponding to z = 0) of spin waves in the harmonic approximation to the XY model. Mixing times (describing the approach towards equilibrium from the least favorable initial state) however remain much larger than equilibrium correlation times at low temperatures. We also describe the respective influence of topological monopole-antimonopole excitations and of spin waves on the event-chain dynamics in the three-dimensional Heisenberg model.

Phase Transitions for Cuboc Orders in Stacked Kagome Heisenberg Systems

Using the event-chain Monte Carlo (MC) algorithm, we investigate phase transitions of the stacked Kagome Heisenberg systems with classical vector spins including up to the 3rd-nearest-neighbor couplings. In particular, we focus on two types of non-coplanar spin orders ---cuboc1 and cuboc2 orders, both of which have twelve-sublattice structures accompanying the translational-symmetry breaking. We perform event-chain MC simulations up to L=72, where L represents the linear dimension of the stacked Kagome lattice, and then find that the cuboc1 transition shows 2nd-order-transition behaviors with a tendency to a weak-1st-order transition up to L=72, while the cuboc2 transition is basically described by the 1st-order transition. We then discuss the above transitions in connection with the effective Landau-Ginzburg-Wilson theory with the O(3)xO(3) symmetry.

Event-chain Monte Carlo algorithms for three- and many-particle interactions

We generalize the rejection-free event-chain Monte Carlo algorithm from many-particle systems with pairwise interactions to systems with arbitrary three- or many-particle interactions. We introduce generalized lifting probabilities between particles and obtain a general set of equations for lifting probabilities, the solution of which guarantees maximal global balance. We validate the resulting three-particle event-chain Monte Carlo algorithms on three different systems by comparison with conventional local Monte Carlo simulations: i) a test system of three particles with a three-particle interaction that depends on the enclosed triangle area; ii) a hard-needle system in two dimensions, where needle interactions constitute three-particle interactions of the needle end points; iii) a semiflexible polymer chain with a bending energy, which constitutes a three-particle interaction of neighboring chain beads. The examples demonstrate that the generalization to many-particle interactions broadens the applicability of event-chain algorithms considerably.

Cell-veto Monte Carlo algorithm for long-range systems.

We present a rigorous efficient event-chain Monte Carlo algorithm for long-range interacting particle systems. Using a cell-veto scheme within the factorized Metropolis algorithm, we compute each single-particle move with a fixed number of operations. For slowly decaying potentials such as Coulomb interactions, screening line charges allow us to take into account periodic boundary conditions. We discuss the performance of the cell-veto Monte Carlo algorithm for general inverse-power-law potentials, and illustrate how it provides a new outlook on one of the prominent bottlenecks in large-scale atomistic Monte Carlo simulations.

Event-chain Monte Carlo algorithm for continuous spin systems and its application

The event-chain Monte Carlo (ECMC) algorithm is described for hard-sphere systems and general potential systems including interacting particle system and continuous spin systems. Using the ECMC algorithm, large-scale equilibrium Monte Carlo simulations are performed for a three-dimensional chiral helimagnetic model under a magnetic field. It is found that critical behavior of a phase transition changes with increasing the magnetic field.

Phase transitions and ordering structures of a model of a chiral helimagnet in three dimensions

Phase transitions in a classical Heisenberg spin model of a chiral helimagnet with the Dzyaloshinskii--Moriya (DM) interaction in three dimensions are numerically studied. By using the event-chain Monte Carlo algorithm recently developed for particle and continuous spin systems, we perform equilibrium Monte Carlo simulations for large systems up to about $10^6$ spins. Without magnetic fields, the system undergoes a continuous phase transition with critical exponents of the three-dimensional \textit{XY} model, and a uniaxial periodic helical structure emerges in the low temperature region. In the presence of a magnetic field perpendicular to the axis of the helical structure, it is found that there exists a critical point on the temperature and magnetic-field phase diagram and that above the critical point the system exhibits a phase transition with strong divergence of the specific heat and the uniform magnetic susceptibility.

Event-chain algorithm for the Heisenberg model: Evidence for z≃1 dynamic scaling.

We apply the event-chain Monte Carlo algorithm to the three-dimensional ferromagnetic Heisenberg model. The algorithm is rejection-free and also realizes an irreversible Markov chain that satisfies global balance. The autocorrelation functions of the magnetic susceptibility and the energy indicate a dynamical critical exponent z≈1 at the critical temperature, while that of the magnetization does not measure the performance of the algorithm. We show that the event-chain Monte Carlo algorithm substantially reduces the dynamical critical exponent from the conventional value of z≃2.

Event-chain Monte Carlo for classical continuous spin models

We apply the event-chain Monte Carlo algorithm to classical continuum spin models on a lattice and clarify the condition for its validity. In the two-dimensional XY model, it outperforms the local Monte Carlo algorithm by two orders of magnitude, although it remains slower than the Wolff cluster algorithm. In the three-dimensional XY spin glass model at low temperature, the event-chain algorithm is far superior to the other algorithms.

Monte Carlo simulation of dense polymer melts using event chain algorithms.

We propose an efficient Monte Carlo algorithm for the off-lattice simulation of dense hard sphere polymer melts using cluster moves, called event chains, which allow for a rejection-free treatment of the excluded volume. Event chains also allow for an efficient preparation of initial configurations in polymer melts. We parallelize the event chain Monte Carlo algorithm to further increase simulation speeds and suggest additional local topology-changing moves ("swap" moves) to accelerate equilibration. By comparison with other Monte Carlo and molecular dynamics simulations, we verify that the event chain algorithm reproduces the correct equilibrium behavior of polymer chains in the melt. By comparing intrapolymer diffusion time scales, we show that event chain Monte Carlo algorithms can achieve simulation speeds comparable to optimized molecular dynamics simulations. The event chain Monte Carlo algorithm exhibits Rouse dynamics on short time scales. In the absence of swap moves, we find reptation dynamics on intermediate time scales for long chains.

Two-dimensional melting: from liquid-hexatic coexistence to continuous transitions.

The phase diagram of two-dimensional continuous particle systems is studied using the event-chain Monte Carlo algorithm. For soft disks with repulsive power-law interactions ∝r^{-n} with n≳6, the recently established hard-disk melting scenario (n→∞) holds: a first-order liquid-hexatic and a continuous hexatic-solid transition are identified. Close to n=6, the coexisting liquid exhibits very long orientational correlations, and positional correlations in the hexatic are extremely short. For n≲6, the liquid-hexatic transition is continuous, with correlations consistent with the Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) scenario. To illustrate the generality of these results, we demonstrate that Yukawa particles likewise may follow either the KTHNY or the hard-disk melting scenario, depending on the Debye-Hückel screening length as well as on the temperature.

Hard-disk equation of state: First-order liquid-hexatic transition in two dimensions with three simulation methods

We report large-scale computer simulations of the hard-disk system at high densities in the region of the melting transition. Our simulations reproduce the equation of state, previously obtained using the event-chain Monte Carlo algorithm, with a massively parallel implementation of the local Monte Carlo method and with event-driven molecular dynamics. We analyze the relative performance of these simulation methods to sample configuration space and approach equilibrium. Our results confirm the first-order nature of the melting phase transition in hard disks. Phase coexistence is visualized for individual configurations via the orientational order parameter field. The analysis of positional order confirms the existence of the hexatic phase.

Addendum to “Event-chain Monte Carlo algorithms for hard-sphere systems”

We extend the event-chain Monte Carlo algorithm from hard-sphere interactions to general potentials. This event-driven Monte Carlo algorithm is nonlocal and rejection free and allows for the breaking of detailed balance. The algorithm uses a discretized potential, but its running speed is asymptotically independent of the discretization. We apply the algorithm to two-dimensional soft spheres and discuss its possible implementation directly in the continuum limit.