# Hard-disk dipoles and non-reversible Markov chains

Philipp Höllmer<EMAIL_ADDRESS>Bethe Center for Theoretical
Physics, University of Bonn, Nussallee 12, 53115 Bonn, Germany  A. C. Maggs
<EMAIL_ADDRESS>CNRS UMR7083, ESPCI Paris, Université PSL, 10 rue
Vauquelin, 75005 Paris, France  Werner Krauth<EMAIL_ADDRESS>Laboratoire de Physique de l’Ecole normale supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université de Paris, Paris, France

###### Abstract

We benchmark event-chain Monte Carlo (ECMC) algorithms for tethered hard-disk
dipoles in two dimensions in view of application of ECMC to water models in
molecular simulation. We characterize the rotation dynamics of dipoles through
the integrated autocorrelation times of the polarization. The non-reversible
straight, reflective, forward, and Newtonian ECMC algorithms are all event-
driven, and they differ only in their update rules at event times. They
realize considerable speedups with respect to the local reversible Metropolis
algorithm. We also find significant speed differences among the ECMC variants.
Newtonian ECMC appears particularly well-suited for overcoming the dynamical
arrest that has plagued straight ECMC for three-dimensional dipolar models
with Coulomb interactions.

##  I Introduction

Markov-chain Monte Carlo [1] (MCMC) is a computational method for sampling
high-dimensional probability distributions $\pi$ including the Boltzmann
distribution of statistical physics. In MCMC, an initial sample $x$, at a time
$t=0$, is drawn from a distribution $\pi^{\\{0\\}}$ rather than from the
“equilibrium” distribution $\pi$. The sample $x$, at time $t$, moves to the
sample $y$, at time $t+1$, with a conditional probability $P_{xy}$, element of
a time-independent transition matrix $P$. The probability distribution
$\pi^{\\{t\\}}$ then moves to $\pi^{\\{t+1\\}}=\pi^{\\{t\\}}P$. The
intermediate probability distributions are usually not known explicitly, but
they evolve towards $\pi$ for $t\to\infty$. The transition matrix $P$ must
satisfy a balance condition that expresses the stationarity of $\pi$.
Reversible transition matrices satisfy the detailed-balance condition
$\pi_{x}P_{xy}=\pi_{y}P_{yx}$ for all $x$ and $y$. Detailed balance is
equivalent to the statement that the equilibrium flow
$\mathcal{F}_{xy}=\pi_{x}P_{xy}$ from $x$ to $y$ equals the reverse flow from
$y$ to $x$. Non-reversible MCMC algorithms satisfy a weaker global-balance
condition $\pi_{x}=\sum_{y}\pi_{y}P_{yx}$ for all $x$. Under conditions of
irreducibility and aperiodicity, [2] they converge to the equilibrium
distribution $\pi$ with non-zero net flows
$\mathcal{F}_{xy}-\mathcal{F}_{yx}$. Samples are then drawn from the
equilibrium distribution $\pi$ for $t\to\infty$ by a non-equilibrium random
process with non-vanishing net flows.

In past decades, research and applications have focused almost exclusively on
reversible Markov chains (and some close relatives, such as sequential schemes
[3, 4]). Reversible Markov chains are straightforward to set up for arbitrary
probability distributions $\pi$, and are easy to conceptualize, in particular
because of the real-valued eigenvalue spectrum of their transition matrices.
Popular reversible Markov chains include the Metropolis- [1] or heat-bath [5,
6, 7] algorithms. Interest in non-reversible Markov chains has increased in
recent years, ever since it was understood [8, 9] that such algorithms often
approach equilibrium faster than reversible Markov chains. Systematic non-
reversible Markov-chain schemes are becoming available. [10, 11]

Event-chain Monte Carlo [11, 12, 13] (ECMC) is a family of non-reversible
local Markov chains that have been applied to, e.g., particle [14] and spin
systems, [15, 16, 17] polymers [18, 19] and field-theoretical models. [20]
ECMC algorithms are defined in continuous time, and they can be interpreted as
piecewise deterministic Markov processes. [21] They have found applications in
Bayesian statistics. [22, 23] A number of features set ECMC algorithms apart
from reversible Markov chains. First, they may reach the “equilibrium”
distribution $\pi$ on fast ballistic time scales rather than on the diffusive
time scales that are typically associated with local reversibility.
Considerable speedups were demonstrated in analytically solvable test cases
[8, 9, 24] and recovered in applications of practical relevance. [14, 19, 13,
20] Second, ECMC, for a given model (probability distribution $\pi$), covers
closely related local MCMC algorithms that can behave quite differently.
Algorithms adopt a wide range of factor sets for the breakup of the Boltzmann
distribution. On a finer scale, even for a given set of factors, different
variants may have widely different behaviors, as we will discuss in the
present paper. Third, ECMC need not evaluate the ratio $\pi_{y}/\pi_{x}$ to
move from a sample $x$ to a sample $y$ (in other words, need not evaluate the
change in energy), because the factors are statistically independent. [12, 25]

The sampling of the Boltzmann distribution without evaluating the energy is
implemented in a general-purpose open-source ECMC application [26] for
classical molecular simulation in Coulomb systems, where the evaluation of the
energy represents for other methods the computational bottleneck. In the test
case of three-dimensional water with the simple-point-charge-flexible-water
(SPC/Fw) potential [27] and the “periodic” variant of straight ECMC (see
Section III.2 for a definition), the simulated water molecules however
resisted rotation and the polarization remained dynamically arrested for long
times. [28] The present paper studies many variants of ECMC for a simpler
model of polar molecules, and suggests that its periodic variant is ill-suited
to these systems.

In a generic Markov chain, the sample space $\Omega$ is distinct from the
moves $\Delta$. The latter are then part of a move set and lead from samples
$x\in\Omega$ to samples $y\in\Omega$. In the local reversible Metropolis
algorithm for particle systems, for example, $\Delta=(\delta,i)$ consists in a
small random displacement $\delta$ of a random particle $i$. ECMC, in
contrast, belongs to the class of lifted Markov chains, [8, 9] in which parts
of the move set are integrated into a “lifted” sample space. In the above
example, the proposed move $\Delta^{\prime}=(\delta^{\prime},i^{\prime})$ at
time $t+1$, after a move $\Delta$, at time $t$, is then no longer
independently sampled, but it rather depends on the lifted configuration, and
in particular on $\Delta$. The probability distribution lives itself in an
extended space. Nevertheless, in the cases that we consider here, the
stationary probability distribution in the lifted space factorizes into the
original distribution $\pi$ and a distribution for the lifting variables,
while the transition matrix non-trivially couples the two sectors.

In this paper we caricature the SPC/Fw water model as tethered hard-disk
dipoles. This replaces a three-dimensional model by a two-dimensional one, and
lumps the vibrational, bending, Lennard-Jones, and Coulomb interactions into a
hard-disk potential. The simplified dipole model retains the polar nature of
molecules, and it can be simulated orders of magnitude faster than the system
of flexible Coulomb dipoles. It also greatly simplifies the implementation of
the different algorithms. This allows us to scan straight, [11, 29]
reflective, [11] forward, [30] and Newtonian [31] ECMC with thousands of
parameter sets. We benchmark the ECMC algorithms against the local Metropolis
algorithms with different displacement sets. We characterize the speed of
algorithms _via_ the autocorrelation of the polarization, the mode that with
ECMC equilibrates slowly in the SPC/Fw water model. [28]

For dipole-model parameters that roughly correspond to those of SPC/Fw water,
we find a $50$-fold speedup for the fastest non-reversible ECMC algorithm with
respect to the optimized local reversible Metropolis algorithm. We also find a
more than order-of-magnitude speed difference among the different ECMC
variants. The different versions of straight ECMC are clearly the slowest ones
at low density, while they are comparable with forward and reflective ECMC at
high densities. Newtonian ECMC is by far the fastest. We suggest that this is
due to its absence of intrinsic parameters.

The content of this paper is as follows: In Section II, we define the two-
dimensional tethered hard-disk dipole model and motivate its parameters with
respect to physical systems of the SPC/Fw model in three dimensions. We also
discuss the polarization, which tracks the ability of dipoles to rotate. In
Section III, we discuss the reversible and non-reversible MCMC algorithms that
we have implemented and illustrate their behavior for the case of a single
dipole. We also discuss subtle aspects concerning the irreducibility of the
Markov chains. In Section IV, we benchmark the different variants of ECMC
against the local reversible Metropolis algorithm for a range of densities and
system sizes, and we evidence the superiority of the non-reversible ECMC
algorithms. In Section V, we project how our conclusions can be extended to
the more complex water models.

##  II Hard-disk dipole model

We consider tethered dipoles in a two-dimensional square box of sides $L$ with
periodic boundary conditions. Each dipole consists of two hard disks of radius
$\sigma$, with a flat inner-dipole interaction that constrains its extension
$\rho$ (the separation of the disk centers) between the contact distance
$r=2\sigma$ and the tether length $R$ (see Fig. 1a). Dipole configurations
without overlapping disks and with all dipole extensions between $r$ and $R$
all have the same Boltzmann weight, whereas all other configurations have zero
statistical weight. As the minimum dipole extension at contact equals
$2\sigma$ (see Fig. 1a), each dipole configuration is also a configuration of
hard disks. Any dipole system is parameterized by the number of dipoles $N$,
the tether ratio $\eta=R/r$, and the hard-disk density
$D=2N\pi\sigma^{2}/L^{2}$.

Figure 1: Tethered hard-disk dipole model. (a): Dipoles with two hard disks of
radius $\sigma$, and minimum extension $r=2\sigma$ and tether length $R$. (b):
Probability distribution of the oxygen–hydrogen (OH) and hydrogen–hydrogen
(HH) bond lengths for a single SPC/Fw water molecule at room temperature. (c):
Five dipoles with polarization $\mathbf{p}_{i}$. (d): Total polarization
$\mathbf{p}$ of (c), together with its probability distribution.

In a single three-dimensional water molecule at room temperature (as for
example described in the SPC/Fw model), the oxygen–hydrogen distance
fluctuates by $\sim$2.3\text{\,}\mathrm{\char 37\relax}$$ and the
hydrogen–hydrogen distance by $\sim$3.4\text{\,}\mathrm{\char 37\relax}$$ (see
Fig. 1b). In the hard-disk dipole model, we thus adopt a tether ratio
$\eta=1.1$, for which the $\sim$2.75\text{\,}\mathrm{\char 37\relax}$$
fluctuations of the extension are quite similar. The tether length $R$ is thus
only $10\%$ larger than the contact distance between two disks, so that two
dipoles cannot lock into a crossed state that would be difficult to
disentangle in two dimensions. Our results for the dynamics of two-dimensional
models may well extend to three dimensions.

We measure dynamical information through the integrated autocorrelation
functions of components of the total polarization $\mathbf{p}$:

$\mathbf{p}=\sum_{i=1}^{N}\mathbf{p}_{i},$ (1)

where $\mathbf{p}_{i}$ is the oriented separation vector between the two disks
in the $i$th dipole, possibly corrected for periodic boundary conditions (see
Fig. 1c and d).

For $N=1$, analytic results are available for the dynamics of some of the
algorithms discussed here. [29] The center-of-mass motion of the single dipole
then factors out, and the sample space $\Omega$ simply corresponds to the
possible values of the polarization vector, in other words to the two-
dimensional ring with inner radius $r$ and outer radius $R$ (see Ref. 29). The
single-dipole polarization $\mathbf{p}$ is uniformly distributed in this ring.

##  III Reversible and non-reversible MCMC

In the present section, we introduce MCMC algorithms that will be used in
Section IV. We also illustrate them for a single dipole ($N=1$) and review
some of their fundamental properties. We only consider algorithms where at
each moment a single disk moves, as this can also be realized within ECMC in
the presence of long-range Coulomb interactions. [25, 26] Dipole
rotations—that relax the polarization—are thus pieced together from subsequent
displacements of single disks. This has already proven efficient in ECMC for
polymer models and for elongated hard needles. [18, 19] Explicit dipole
rotations appear difficult to put into place in this context.

The local MCMC moves discussed here feature small single-disk displacements at
each time step (for the Metropolis algorithm) or continuous moves in time (for
the ECMC algorithms). Again, this is motivated by the future applications to
more complicated systems such as water, for which non-local moves including
cluster updates [32, 33] appear out of reach. Although we thus restrict our
attention to a subclass of MCMC algorithms with very similar design
principles, we will show in Section IV that their properties differ
substantially.

In the ECMC algorithms, a single active disk moves with a constant velocity in
a given direction. The straight-line trajectories stop at events which
correspond to hard-disk collisions or to the moment when the dipole extension
reaches the tether length. In addition, resamplings take place at predefined
MCMC times, which are separated by a chain time $\tau_{\text{chain}}$.
Resamplings also interrupt the straight-line trajectory, and for example draw
a new velocity and a new active disk.

###  III.1 Reversible local Metropolis algorithms

The local Metropolis algorithm, [1] for the dipole model, selects a random
disk $i$ at each time step. For the cross-shaped displacement set, the
proposed move for $i$ is sampled uniformly between $-\delta$ and $\delta$ with
range $\delta$, randomly in $x$\- or in $y$-direction. For the square-shaped
displacement set, the proposed move is uniformly sampled from a square of
sides $2\delta$ centered at zero, with both the $x$ and the $y$ components of
the displacement sampled independently between $-\delta$ and $\delta$. In both
cases, the move is accepted if it violates no constraints. If the move is
rejected, the configurations at times $t$ and $t+1$ are the same. The
Metropolis algorithm, with the given displacement sets, satisfies the
detailed-balance condition. For $N=1$, the polarization $\mathbf{p}$ performs
a random walk in the sample space $\Omega$. The algorithm can be considered
irreducible for all $N$, though some blocked configurations [34] may pose
problems. [35]

###  III.2 Straight ECMC

In straight ECMC, [11] the velocity is of constant absolute value, and it
changes only at resamplings. We analyze several variants for updating the
angle $\phi$ that describes the velocity in two dimensions. In the “periodic”
variant, the angle alternates between $\phi=0$ and $\phi=\pi/2$. In arbitrary
dimension, the velocity cycles through the unit vectors aligned with the
positive coordinate axes. The periodic variant, used in most previous works,
is implemented in the general-purpose ECMC application. [26] In the present
paper, we also analyze the “random” variant where $\phi$ is uniformly sampled
in the interval $[0,2\pi)$ and, finally, the “sequential” variant,[29] where
$\phi$ is incremented by a constant $\Delta\phi$ at each resampling. The
polarization trajectory between resamplings of straight ECMC is a straight
line. In the periodic variant, subsequent straight-line trajectories are at
right angles. They are at an angle $\Delta\phi$ in the sequential variant (see
Fig. 2a for an example in $N=1$). The described variants have very different
time evolutions. The polarization trajectories of the sequential variant, for
example, persistently rotate for $N=1$, with a total rotation angle that
diverges as $\Delta\phi$ goes to zero. [29]

Figure 2: Sample trajectories for a single tethered dipole. The polarization
$\mathbf{p}$ is the oriented vector between the two disk centers. (a):
Sequential variant of straight ECMC. [29] (b): Reflective ECMC. [11] (c):
Forward ECMC. [30] (d): Newtonian ECMC. [31]

At a resampling of straight ECMC at time $t_{\text{res}}$ (usually a multiple
of $\tau_{\text{chain}}$), the angle $\phi$ is updated, but the active disk
must also be chosen randomly in order to satisfy global balance. To motivate
this, we show for $N=1$ that the two disks $1$ and $2$ must be active with the
same probability at time $t_{\text{res}}$. Up until the next resampling, all
polarizations $\mathbf{p}(t)$ lie on a segment of the ring $\Omega$, and each
of the two directions of motion on that segment corresponds to one of the two
disks being active (see Fig. 3a). Global balance requires the flows into
$\mathbf{p}(t)$ to be independent of its position. It can be reached from two
polarizations $\mathbf{p}(t_{\text{res}})$ and
$\mathbf{p}^{\prime}(t_{\text{res}})$. Depending on the value of
$\mathbf{p}(t)$, the corresponding polarizations $\mathbf{p}(t_{\text{res}})$
and $\mathbf{p}^{\prime}(t_{\text{res}})$ either both correspond to disk $1$
being active at $t_{\text{res}}$, or both to disk $2$, or one to $1$ and one
to $2$ (see Fig. 3a). Disks $1$ and $2$ must thus be equally likely
($p_{1}=p_{2}$) to be active right after the resampling, a “random”-mode
condition that is precisely implemented by the resampling of the active disk
at $t_{\text{res}}$.

Figure 3: Global balance and irreducibility for sequential ECMC. (a): Flow
into polarization $\mathbf{p}(t)$ from $\mathbf{p}(t_{\text{res}})$ and
$\mathbf{p}^{\prime}(t_{\text{res}})$ on a segment of $\Omega$ as in Fig. 2a.
Events take place at the boundaries of the segment, and the active disk
changes. (b): Single-dipole distribution $\pi(\mathbf{p}/\sigma)$ for
sequential ECMC with small $\Delta\phi$. The snake-mode distribution is
incorrect, and the probability $p_{1}$ that disk $1$ is active right after a
resampling is smaller than $p_{2}$. (c): Same as (b) but for different initial
configuration and chain time. The snake-mode distribution is incorrect
although both disks are equally likely to be active after a resampling
($p_{1}=p_{2}$).

To illustrate the relevance of the random-mode requirement, we test for $N=1$
the case when the active disk is maintained at a resampling while the angle
$\phi$ is incremented in the sequential variant of straight ECMC (“snake”
mode). Clearly, snake mode is incorrect, as it samples $\Omega$ non-uniformly
($p_{1}\neq p_{2}$, see Fig. 3b). It seems to be incorrect even when both
disks happen to be initially active with practically the same probability
($p_{1}=p_{2}$, see Fig. 3c). In both cases, the snake-mode trajectory is
repetitive, and lack of irreducibility in the lifted sample space
$\Omega\times\\{\phi\\}$ induces a non-uniform distribution in $\Omega$. For
$N>1$ at very low density, sequential ECMC in snake mode might not be
irreducible, because a single dipole deterministically rotates on a closed
trajectory. However, no difference between snake mode and random mode is
detected for $N>1$ at the higher densities considered in Section IV.

###  III.3 Reflective ECMC

Reflective ECMC [11] differs from the straight ECMC only in its handling of
events. At an event, the velocity of the incoming active disk is reflected by
the line connecting the active and the target disks, and it becomes the
velocity of the target disk. The absolute value of the velocity is thus
preserved. For $N=1$, the polarization follows straight lines which are
symmetrically reflected off the inner and outer ring boundaries, with a single
reflection angle $\alpha$ (see Fig. 2b).

Reflective ECMC satisfies global balance for any number of dipoles. For $N=1$,
and for an initial velocity that is nearly perpendicular to the separation
between the two disks, the polarization may follow a “whispering-gallery-mode”
trajectory which never visits small separations $\rho\gtrsim r$ (compare with
Ref. 23, Fig. 3), showing that resamplings of the active disk and of its
velocity angle are required for irreducibility. For all initial velocities at
$N=1$, even for those that visit all of sample space $\Omega$, the probability
distribution in the lifted sample space does not separate into independent
distributions, and the sampled stationary distribution in $\Omega$ is non-
uniform, and thus incorrect. With periodic resamplings, the stationary
distribution is uniform in $\Omega$. For $N>1$ dipoles, in a large enough box,
reflective ECMC fails to be irreducible without resamplings, as the dynamics
is then deterministic and a single dipole may effectively rotate in a
stationary closed trajectory. Nevertheless, at high enough density, we have
not observed any difference in static observables with or without resamplings.
It thus appears safe to omit them in practical simulations. As we will see in
Section IV, resamplings do not improve the decorrelation. Without resamplings,
reflective ECMC is free of intrinsic parameters, and somewhat easier to set
up. Computational overhead is avoided.

In previous work on hard disks [11] (rather than hard-disk dipoles),
reflective ECMC was benchmarked against straight ECMC and found to be
considerably slower. In contrast, in the dipole model studied here, reflective
ECMC will prove more powerful than straight ECMC, even when the latter runs
with optimized intrinsic parameters.

###  III.4 Forward ECMC

Forward ECMC [30] resembles reflective ECMC in that, after an event, the
target-disk velocities of both algorithms are located in the same quadrant of
the coordinate system with axes parallel and orthogonal to the line connecting
the disks at the event. All velocities are of unit absolute value, but forward
ECMC incorporates a random element into each event. More precisely, the
component orthogonal to the line connecting the disks at contact is uniformly
sampled between $0$ and $1$ (while reflecting the orthogonal orientation). Its
parallel component is determined so that the velocity is of unit norm. For a
single dipole, the polarization describes straight lines that are reflected
off the inner and outer ring boundaries. However, the outgoing reflection
angle is non-deterministic (see Fig. 2c). Forward ECMC is irreducible even for
$N=1$, and it requires no resamplings. Also, as we will show in Section IV,
resamplings of the velocity and the active disk in intervals of the chain time
do not speed up the algorithm.

###  III.5 Newtonian ECMC

Newtonian ECMC [31] mimics event-driven molecular dynamics. [36] All disks
have velocity labels in addition to their positions. The label of the active
disk indicates the time derivative of its position, that is, its displacement
in time. All other disks are stationary. Events (including those where the
maximum dipole separation is reached) take place as in molecular dynamics,
with all velocity labels treated on an equal footing, and with equal masses
for all disks. The identities of the active and the target disks are
interchanged in the event. During the simulation, the absolute value of the
active-disk velocity varies. The equilibrium distribution of all velocity
labels is uniform on the $2N$-dimensional unit sphere, and the root-mean-
square velocity $v_{\text{rms}}$ is conserved. The event rate per unit
distance of the Newtonian ECMC equals the one of the other variants, but the
event rate per unit time is smaller by $\sqrt{\pi}/2$, because of the
difference between $v_{\text{rms}}$ and the mean absolute value for the two-
dimensional Gaussian distribution of velocities. All the velocity labels are
used as lifting variables. At a resampling, these labels must be sampled from
the exact equilibrium distribution (the rescaled Maxwell–Boltzmann
distribution for all disks). Without resampling, the choice of the initial
velocity labels is arbitrary.

For $N=1$, the polarization trajectory of Newtonian ECMC reflects off the
inner and outer ring boundaries but considers both velocity labels at each
event (see Fig. 2d). Without resamplings, Newtonian ECMC breaks irreducibility
for this single dipole, and it samples an incorrect stationary probability
distribution. For $N\geq 2$, we find no influence of observable means on the
resampling rate, so that the irreducibility problem is probably again due to
the high symmetry of the ring geometry in Section 2. In the larger-$N$ dipole
systems in Section IV, we furthermore notice that the algorithm becomes faster
in the limit of infinite resampling time.

##  IV Autocorrelations for $N$ dipoles

In the present section, we characterize the local MCMC algorithms of Section
III via the autocorrelation dynamics of the polarization. All correlation
times $\tau_{\text{int}}$ refer to averaged integrated autocorrelation times
of single components of the polarization. The unit of time corresponds to a
single trial move for the Metropolis algorithm or to the mean event time for
ECMC. In ECMC, the computational effort is $\mathcal{O}(1)$ per event, as it
is per move in the Metropolis algorithm. For this reason, no effort is made to
compare actual CPU times, which would depend on the implementation.
Statistical errors are estimated from the difference between autocorrelation
times of the $x$ and $y$ components of the polarization which, by symmetry,
must be the same. For ECMC, we also take into account the uncertainty in the
mean event time. Single MCMC-run times, for each parameter set, are at least
three orders of magnitude longer than $\tau_{\text{int}}$.

###  IV.1 Intrinsic parameters of Metropolis and straight ECMC

The correlation time of the local Metropolis algorithm depends on its
intrinsic parameters, namely the range $\delta$ and the displacement set (see
Fig. 4a). The acceptance rate decreases with increasing $\delta$. The
correlation time can thus be expressed as a function of the acceptance rate.
As in many comparable models, the “one-half” rule [37] is roughly respected,
and the Metropolis algorithm converges best for a rejection rate on the order
of $50\%$. For $N=81$ dipoles at density $D=0.70$, we observe a broad optimum
between $20\text{\,}\mathrm{\char 37\relax}40\text{\,}\mathrm{\char 37\relax}$
for the square-shaped displacement set and an optimum close to $50\%$ for the
cross-shaped displacement set (see Fig. 4b).

Figure 4: Correlation time $\tau_{\text{int}}$ of the Metropolis algorithm
with cross- and square-shaped displacement sets for $N=81$ dipoles at density
$D=$0.70$$. (a): $\tau_{\text{int}}$ vs. the range $\delta$. (b):
$\tau_{\text{int}}$ vs. acceptance rate. The minimum of $\tau_{\text{int}}$
provides the reference in Fig. 5, that the ECMC algorithms are benchmarked
against.

All ECMC variants allow for resamplings so that their correlation times depend
on the intrinsic parameter $\tau_{\text{chain}}$. The performance of the
sequential variant of straight ECMC also depends on the angle increment
$\Delta\phi$. Without resamplings, straight ECMC is not irreducible and its
correlation time is infinite. The optimum $\tau_{\text{int}}$ is thus at a
finite $\tau_{\text{chain}}$, to be obtained by fine-tuning (see upper curves
in Fig. 5). In its sequential variant, the cases $\Delta\phi=\pi/2$ and
$\Delta\phi=\pi/3$ make the velocity cycle through a few values only. For
these cases, we notice that $\tau_{\text{int}}$ has two local minima, both of
which are rather large. In the absence of a heuristic for the choice of
intrinsic parameters, straight ECMC requires explicit fine-tuning, which
increases its computational complexity.

As discussed in Section III, reflective, forward and Newtonian ECMC appear
irreducible for $N>1$ even without resampling and are fastest in this limit
$\tau_{\text{chain}}\to\infty$ (see lower curves in Fig. 5). This was also
observed, e.g., for reflective ECMC in two-dimensional hard-disk systems, [11]
and for Newtonian ECMC in three-dimensional hard spheres. [31] For moderate
values of $\tau_{\text{chain}}$, Newtonian ECMC is comparable to other
variants. However, for much larger $\tau_{\text{chain}}$, its correlation
times again decrease strongly (see inset of Fig. 5). For $N=81$ dipoles at
density $D=0.70$, the variants that require no fine-tuning perform best.
Newtonian ECMC without resampling accelerates with respect to the local
Metropolis algorithm by a surprising factor $50$ at density $0.70$, that seems
to further increase at even higher density (see Section IV.2).

Figure 5: Correlation times of ECMC variants for different intrinsic
parameters (chain time $\tau_{\text{chain}}$, angle increment $\Delta\phi$)
for $N=81$ dipoles at density $D=$0.70$$. Optimal values of
$\tau_{\text{int}}$ are highlighted on the right $y$-axis. For reflective,
forward, and Newtonian ECMC, they agree with those for
$\tau_{\text{chain}}\to\infty$. The optimal performance of the Metropolis
algorithm is indicated as the benchmark reference (see Fig. 4). The inset
illustrates the evolution of $\tau_{\text{int}}$ for large
$\tau_{\text{chain}}$ for Newtonian ECMC.

###  IV.2 Density and size dependence of speedups

For both versions of the Metropolis algorithms, $\tau_{\text{int}}$ naturally
increases at higher densities and climbs steeply for $D\gtrsim 0.70$, around
the liquid–hexatic phase-transition density of the hard-disk system without
the tether constraint (see Fig. 6a, again for $N=81$). For the Metropolis
algorithm, the square-shaped displacement set is somewhat preferable over the
cross-shaped displacement set for small densities, but the two become equally
fast at large densities (see Fig. 6b). This illustrates that, while step-size
control is a key feature in reversible MCMC, the details of a local reversible
algorithm do not really depend on the displacement sets. For simplicity, we
chose the optimal acceptance rates for $D=0.70$ for these runs at other
densities (see Fig. 4).

Figure 6: Performance of the Metropolis algorithm for $N=81$ dipoles
(optimized intrinsic parameters). (a): Correlation time $\tau_{\text{int}}$
vs. density $D$ for the cross- and square-shaped displacement sets. (b):
Speedup for the square-shaped displacement set compared to the cross-shaped
displacement set.

For the correlation times of the straight ECMC variants, we perform systematic
scans as in Section IV.1, for each density at $N=81$. For forward, reflective,
and Newtonian ECMC, we use $\tau_{\text{chain}}=\infty$ without any fine-
tuning. For all ECMC variants, $\tau_{\text{int}}$ increases with density (see
Fig. 7a). The periodic variant, which is fastest for two-dimensional hard
disks, is the slowest variant for hard-disk dipoles. The random and sequential
variants are somewhat faster. Reflective and forward ECMC resemble each other
in performance at all $D$. While they are considerably faster than the
straight variants at low $D$, this gap vanishes at $D=$0.72$$. Newtonian ECMC
is the fastest algorithm for all densities and, as mentioned, it requires no
fine-tuning. All ECMC algorithms outperform the Metropolis algorithm by a
considerable margin (see Fig. 7b). With increasing density, the speedup of
reflective and of forward ECMC drops from roughly $30$ to $15$, while the
speedup of the straight ECMC variants increases with density until they become
comparable. These trends were similarly observed near the liquid–hexatic
transition of two-dimensional hard disks for periodic and reflective ECMC [11]
(where the periodic variant, however, was considerably faster than reflective
ECMC at high densities). The speedup realized by Newtonian ECMC improves
slightly from low to high densities, as was likewise found for three-
dimensional hard spheres. [31] We observe the highest speedup of roughly a
factor of $60$ at the highest density $D=$0.72$$ that we studied.

Figure 7: Performance of ECMC variants and of the Metropolis algorithm with a
square-shaped displacement set (where applicable: with optimized intrinsic
parameters). (a): Correlation time $\tau_{\text{int}}$ for $N=81$ dipoles vs.
density $D$. (b): Speedup of ECMC with respect to the Metropolis algorithm.
(c): Correlation time $\tau_{\text{int}}$ at density $D=$0.70$$ vs. number of
dipoles $N$.

We finally study the dependence of the correlation time $\tau_{\text{int}}$ of
the number of dipoles $N$ at the density $D=$0.70$$. For the reflective,
forward and Newtonian ECMC variants, we simply set
$\tau_{\text{chain}}=\infty$, given our findings at $N=81$. For periodic,
random and sequential ECMC, we infer the optimal $\tau_{\text{chain}}$ and
$\Delta\phi$ from the $N=81$ case (performing occasional sweeps through
intrinsic parameters as cross-checks). For the Metropolis algorithm, we used
the optimal acceptance rate of $N=81$. We observe that $\tau_{\text{int}}$
increases linearly with $N$ for all algorithms (see Fig. 7c). This implies
that the observed speedups of the non-reversible algorithms stay constant.

##  V Conclusions and outlook

In this paper, we have systematically studied local MCMC algorithms for
complex molecules that we caricatured as hard-disk dipoles with parameters
inspired by the SPC/Fw water system. Moving from three to two spatial
dimensions and from a liquid of charged molecules to a model of hard-disk
dipoles has greatly simplified the algorithm development and its
implementation. This allowed us to scan thousands of parameters sets for
different ECMC variants and to benchmark them against the reversible local
Metropolis algorithm. At this exploratory stage, this would have required
prohibitive computing resources for the original three-dimensional model. [28]
The broad spread of the performance of different ECMC versions came as a
surprise. We found that algorithms with intrinsic parameters perform
considerably less well than others and demonstrated a 60-fold speedup of
Newtonian ECMC with respect to the Metropolis algorithm. Our study has also
brought out subtleties of ECMC, that are hidden in simple liquids. For the
case of a single dipole, we discovered strict requirements for resampling in
order to ensure irreducibility of straight, reflective and Newtonian ECMC.
However, these requirements did not seem to play a role for denser systems of
dipoles. Event-based randomness as in forward ECMC (or as in all ECMC
simulations of soft interactions) strictly ensures irreducibility.

We believe that our observations will carry over from the caricature two-
dimensional dipoles to the three-dimensional SPC/Fw systems and related
models. There are two reasons why we have not yet studied this system
directly. One is, as mentioned, the scale of the computing requirements for a
full-fledged three-dimensional scan. The other is that further algorithm
development is still required. The cell-veto algorithm [38] (which reduces the
computational complexity of the Coulomb problem to $\mathcal{O}(N\log N)$ per
sweep of events, without evaluating the Coulomb energy) is in the present
version of our open-source project only implemented for the periodic variant
of straight ECMC. However, we expect it to generalize to all ECMC variants for
flexible water molecules with an explicit Coulomb interaction. The benchmark
against Metropolis MCMC will then be much more favorable, as the change of the
Coulomb energy can there only be computed in $\mathcal{O}(N^{3/2})$ per sweep
of moves. [39] The key question will be whether the dynamic arrest of straight
ECMC for the three-dimensional water system can be overcome as in the two-
dimensional hard-disk dipole model studied here. The extreme dependence of the
performance of non-reversible MCMC on details of the algorithm was also
evidenced in the escape times from a tightly confined initial configuration,
[35] which might be relevant for overcoming dynamic arrest.

Finally, we point out that any tethered-dipole configuration is also a valid
hard-disk configuration, with an added upper limit on the distance between
pairs of hard disks. For the tethered dipole system, the global hexatic
orientational order parameter $\psi_{6}$ (see Refs. 40, 14), at densities
$D\gtrsim 0.70$, decorrelates faster than for the hard-disk system which there
already enters the mixed hexatic–fluid phase. This observed effect might allow
for an improved sampling algorithm for hard disks. It is also unclear why
straight ECMC performs worse for dipoles than the reflective, forward, and
Newtonian variants, while for simple hard disks straight ECMC is clearly the
most efficient.

###### Acknowledgements.

P. H. acknowledges support from the Studienstiftung des deutschen Volkes and
from Institut Philippe Meyer. W. K. acknowledges support from the Alexander
von Humboldt Foundation.

## Additional remarks

The figures in this paper were prepared with the Matplotlib [41] and Visual
Molecular Dynamics [42] (VMD) software packages. Integrated autocorrelation
times were computed by the Python emcee package. [43]

## References

  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of State Calculations by Fast Computing Machines,” J. Chem. Phys. 21, 1087–1092 (1953).
  * Levin, Peres, and Wilmer [2008] D. A. Levin, Y. Peres, and E. L. Wilmer,  _Markov Chains and Mixing Times_ (American Mathematical Society, 2008).
  * Berg [2004] B. A. Berg,  _Markov Chain Monte Carlo simulations and their statistical analysis: with web-based Fortran code_ (World Scientific, 2004).
  * Ren, O’Keeffe, and Orkoulas [2007] R. Ren, C. J. O’Keeffe, and G. Orkoulas, “Sequential Metropolis Algorithms for Fluid Simulations,” Int. J. Thermophys. 28, 520–535 (2007).
  * Glauber [1963] R. J. Glauber, “Time-Dependent Statistics of the Ising Model,” J. Math. Phys. 4, 294–307 (1963).
  * Creutz [1980] M. Creutz, “Monte Carlo study of quantized SU(2) gauge theory,” Phys. Rev. D 21, 2308–2315 (1980).
  * Geman and Geman [1984] S. Geman and D. Geman, “Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images,” IEEE Trans. Pattern Anal. Mach. Intell. PAMI-6, 721–741 (1984).
  * Diaconis, Holmes, and Neal [2000] P. Diaconis, S. Holmes, and R. M. Neal, “Analysis of a nonreversible Markov chain sampler,” Ann. Appl. Probab. 10, 726–752 (2000).
  * Chen, Lovász, and Pak [1999] F. Chen, L. Lovász, and I. Pak, “Lifting Markov Chains to Speed up Mixing,” Proceedings of the 17th Annual ACM Symposium on Theory of Computing , 275 (1999).
  * Turitsyn, Chertkov, and Vucelja [2011] K. S. Turitsyn, M. Chertkov, and M. Vucelja, “Irreversible Monte Carlo algorithms for efficient sampling,” Physica D 240, 410–414 (2011).
  * Bernard, Krauth, and Wilson [2009] E. P. Bernard, W. Krauth, and D. B. Wilson, “Event-chain Monte Carlo algorithms for hard-sphere systems,” Phys. Rev. E 80, 056704 (2009).
  * Michel, Kapfer, and Krauth [2014] M. Michel, S. C. Kapfer, and W. Krauth, “Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps,” J. Chem. Phys. 140, 054116 (2014).
  * Krauth [2021] W. Krauth, “Event-Chain Monte Carlo: Foundations, Applications, and Prospects,” Front. Phys. 9, 229 (2021).
  * Bernard and Krauth [2011] E. P. Bernard and W. Krauth, “Two-Step Melting in Two Dimensions: First-Order Liquid-Hexatic Transition,” Phys. Rev. Lett. 107, 155704 (2011).
  * Michel, Mayer, and Krauth [2015] M. Michel, J. Mayer, and W. Krauth, “Event-chain Monte Carlo for classical continuous spin models,” EPL 112, 20003 (2015).
  * Nishikawa  _et al._ [2015] Y. Nishikawa, M. Michel, W. Krauth, and K. Hukushima, “Event-chain algorithm for the Heisenberg model: Evidence for $z\simeq 1$ dynamic scaling,” Phys. Rev. E 92, 063306 (2015).
  * Kimura and Higuchi [2017] K. Kimura and S. Higuchi, “Anomalous diffusion analysis of the lifting events in the event-chain Monte Carlo for the classical XY models,” EPL 120, 30003 (2017).
  * Müller, Kampmann, and Kierfeld [2020] D. Müller, T. A. Kampmann, and J. Kierfeld, “Chaining of hard disks in nematic needles: particle-based simulation of colloidal interactions in liquid crystals,” Sci. Rep. 10, 12718 (2020).
  * Kampmann  _et al._ [2021] T. A. Kampmann, D. Müller, L. P. Weise, C. F. Vorsmann, and J. Kierfeld, “Event-Chain Monte-Carlo Simulations of Dense Soft Matter Systems,” Front. Phys. 9, 96 (2021).
  * Hasenbusch and Schaefer [2018] M. Hasenbusch and S. Schaefer, “Testing the event-chain algorithm in asymptotically free models,” Phys. Rev. D 98, 054502 (2018).
  * Davis [1984] M. H. A. Davis, “Piecewise-Deterministic Markov Processes: A General Class of Non-Diffusion Stochastic Models,” J. R. Stat. Soc. Series B Stat. Methodol. 46, 353–376 (1984).
  * Bierkens, Fearnhead, and Roberts [2019] J. Bierkens, P. Fearnhead, and G. Roberts, “The Zig-Zag process and super-efficient sampling for Bayesian analysis of big data,” Ann. Stat. 47, 1288–1320 (2019).
  * Bouchard-Côté, Vollmer, and Doucet [2018] A. Bouchard-Côté, S. J. Vollmer, and A. Doucet, “The Bouncy Particle Sampler: A Nonreversible Rejection-Free Markov Chain Monte Carlo Method,” J. Am. Stat. Assoc. 113, 855–867 (2018).
  * Lei and Krauth [2018] Z. Lei and W. Krauth, “Mixing and perfect sampling in one-dimensional particle systems,” EPL 124, 20003 (2018).
  * Faulkner  _et al._ [2018] M. F. Faulkner, L. Qin, A. C. Maggs, and W. Krauth, “All-atom computations with irreversible Markov chains,” J. Chem. Phys. 149, 064113 (2018).
  * Höllmer  _et al._ [2020] P. Höllmer, L. Qin, M. F. Faulkner, A. Maggs, and W. Krauth, “JeLLyFysh-Version1.0 — a Python application for all-atom event-chain Monte Carlo,” Comput. Phys. Commun. 253, 107168 (2020).
  * Wu, Tepper, and Voth [2006] Y. Wu, H. L. Tepper, and G. A. Voth, “Flexible simple point-charge water model with improved liquid-state properties,” J. Chem. Phys. 124, 024503 (2006).
  * Qin [2020] L. Qin,  _Application of irreversible Monte Carlo in realistic long-range systems_ , Phd thesis, Université Paris Sciences et Lettres (2020).
  * Qin, Höllmer, and Krauth [2020] L. Qin, P. Höllmer, and W. Krauth, “Fast sequential Markov chains,” (2020), arXiv:2007.15615 [cond-mat.stat-mech] .
  * Michel, Durmus, and Sénécal [2020] M. Michel, A. Durmus, and S. Sénécal, “Forward Event-Chain Monte Carlo: Fast Sampling by Randomness Control in Irreversible Markov Chains,” J. Comput. Graph. Stat. 29, 689–702 (2020).
  * Klement and Engel [2019] M. Klement and M. Engel, “Efficient equilibration of hard spheres with Newtonian event chains,” J. Chem. Phys. 150, 174108 (2019).
  * Dress and Krauth [1995] C. Dress and W. Krauth, “Cluster algorithm for hard spheres and related systems,” J. Phys. A: Math. Gen. 28, L597–L601 (1995).
  * Liu and Luijten [2004] J. Liu and E. Luijten, “Rejection-Free Geometric Cluster Algorithm for Complex Fluids,” Phys. Rev. Lett. 92, 035504 (2004).
  * Börözky [1964] K. Börözky, “ Über stabile Kreis- und Kugelsysteme,” Ann. Univ. Sci. Budapest. Eötvös Sect. Math. 7, 79–82 (1964).
  * Höllmer  _et al._ [2021] P. Höllmer, N. Noirault, B. Li, A. C. Maggs, and W. Krauth, “Sparse hard-disk packings and local markov chains,” (2021), arXiv:2109.13343 [cond-mat.stat-mech] .
  * Alder and Wainwright [1957] B. J. Alder and T. E. Wainwright, “Phase Transition for a Hard Sphere System,” J. Chem. Phys. 27, 1208–1209 (1957).
  * Krauth [2006] W. Krauth,  _Statistical Mechanics: Algorithms and Computations_ (Oxford University Press, 2006).
  * Kapfer and Krauth [2016] S. C. Kapfer and W. Krauth, “Cell-veto Monte Carlo algorithm for long-range systems,” Phys. Rev. E 94, 031302 (2016).
  * Berthoumieux and Maggs [2010] H. Berthoumieux and A. C. Maggs, “Collective dispersion forces in the fluid state,” EPL 91, 56006 (2010).
  * Nelson and Halperin [1979] D. R. Nelson and B. I. Halperin, “Dislocation-mediated melting in two dimensions,” Phys. Rev. B 19, 2457–2484 (1979).
  * Hunter [2007] J. D. Hunter, “Matplotlib: A 2d graphics environment,” Comput. Sci. Eng. 9, 90–95 (2007).
  * Humphrey, Dalke, and Schulten [1996] W. Humphrey, A. Dalke, and K. Schulten, “VMD – Visual Molecular Dynamics,” J. Mol. Graph. 14, 33–38 (1996).
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, “emcee: The MCMC Hammer,” Publ. Astron. Soc. Pac. 125, 306 (2013).