# Supplementary Materials for
Molecular simulation from modern statistics: Continuous-time, continuous-
space, exact

Philipp Höllmer, A. C. Maggs, Werner Krauth∗

∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
The PDF file includes:
Materials and Methods
Supplementary Text
Fig. S1

## Materials and Methods

In this section, we describe the SPC/Fw water model that we studied in this
paper. We further specify the simulation protocols for the different sampling
algorithms. Finally, we report on our analysis procedure of the time series of
the electric polarization.

### SPC/Fw Water Model

The atomistic flexible simple point-charge water model SPC/Fw defines a water
molecule by three charged interaction sites that represent the oxygen and
hydrogen atoms (_21_). We treat the canonical ensemble in a cubic box with
periodic boundary conditions, i.e., $N$ water molecules with $N_{a}=3N$ atoms
in a periodically repeated cubic box of side length $L$ at a given temperature
$T\sim$300\text{\,}\mathrm{K}$$. Within any water molecule
$i\in\\{1,\ldots,N\\}$, the intramolecular interactions of the SPC/Fw model
consist of two harmonic bond potentials $U_{\text{bond}}^{i,1}$ and
$U_{\text{bond}}^{i,2}$ that lead to fluctuations of the O-H bond lengths
around their equilibrium length. Likewise, a harmonic bending potential
$U^{i}_{\text{bend}}$ yields a fluctuation of the H-O-H opening angle around
an equilibrium value. The intermolecular interactions between two different
water molecules $i$ and $j$ consist of a Lennard-Jones potential
$U^{ij}_{\text{LJ}}$ between the two oxygen atoms, and a Coulomb potential
$U^{ij}_{\text{C}}$ between all nine pairs of charged atoms. The
intermolecular potentials explicitly include the interactions between all
periodic images of the two involved water molecules. The total potential $U$
of an all-atom configuration $\mathbf{x}$ is a sum of these factor potentials:

$U(\mathbf{x})=\sum_{i=1}^{N}\left[U_{\text{bond}}^{i,1}(\mathbf{x})+U_{\text{bond}}^{i,2}(\mathbf{x})+U_{\text{bend}}^{i}(\mathbf{x})\right]+\sum_{i=1}^{N}\sum_{j=1}^{i-1}\left[U_{\text{LJ}}^{ij}(\mathbf{x})+U_{\text{C}}^{ij}(\mathbf{x})\right].$
(1)

The bond and Lennard-Jones factor potentials depend on two atomic positions;
the bending and Coulomb factor potentials depend on three and six atomic
positions, respectively [see also Ref. (_28, Section V A_)]. The empirical
parameters of the different potentials of the SPC/Fw water model are as in
Ref. (_21_).

### Simulation Protocols—Molecular Dynamics

We use the feature release from February 8, 2023, of Lammps (_45, 8_) for the
molecular-dynamics simulations of the SPC/Fw water model in this paper. We use
a spherical cut-off for the Lennard-Jones potential. The Coulomb potential is
treated by a particle–particle particle–mesh solver (_12_). Unless explicitly
specified otherwise (as, e.g., in Fig. 4B), the solver uses a target accuracy
of $10^{-6}$. Coulomb interactions are partly treated in discretized
reciprocal space, whereby the grid size is chosen to meet the target accuracy
based on analytic error estimates obtained from a specific global charge
distribution (_40, 15, 41_). Thermostatting is achieved by integrating
Nosé–Hoover-chain equations of motion for the canonical ensemble (_46_) with a
time-reversible measure-preserving Verlet integrator (_47_). The time step is
$1\text{\,}\mathrm{fs}$, and the temperature is relaxed in a timespan of
roughly $300$ time steps. This simulation protocol for molecular dynamics in
Lammps is validated by comparing numerical results for $N=2$ SPC/Fw water
molecules of all sampling methods of this paper (see Fig. S1).

### Simulation Protocols—Metropolis Algorithm

We use version $2.07$ of the DL_MONTE software package (_48, 39, 49_) to
sample the canonical ensemble of the SPC/Fw water model with the reversible
Metropolis algorithm (_25_), and minimally amended the software to output the
electric polarization. We use a spherical cut-off for the Lennard-Jones
potential. The Coulomb interaction is treated by an Ewald summation (_50_)
with a target-accuracy tolerance of $10^{-6}$. The simulations with the event-
chain Monte Carlo algorithm in this paper move a single atom at a time. In
order to directly compare the dynamics of the reversible Metropolis algorithm
with the non-reversible event-chain Monte Carlo algorithm (see Fig. 3C), each
trial proposes a new position of a single random atom in a proposal cube
around its original position. The new position is accepted with the usual
Metropolis criterion based on the change of energy. The size of the proposal
cube is separately adapted for the hydrogen and oxygen atoms during the
simulation to obtain a target acceptance rate of $37\text{\,}\mathrm{\char
37\relax}$. This simulation protocol for the Metropolis algorithm in DL_MONTE
is validated by comparing numerical results for $N=2$ SPC/Fw water molecules
of all sampling methods of this paper (see Fig. S1).

### Simulation Protocols—Event-Chain Monte Carlo

We develop version $2.0$ of JeLLyFysh (_51, 32_) for the event-chain Monte
Carlo simulations of the SPC/Fw water model in this paper. In addition to the
straight variant implemented in versions $<2.0$, version 2.0 implements the
generalized Newtonian event-chain Monte Carlo variant (see Supplementary Text
for details). We use the factorization in eq. (1). In particular, we group all
the charge–charge Coulomb interactions between two water molecules into a
single factor. Theory then predicts an optimal logarithmic increase of the
number of events from these molecular Coulomb factors (_28_) (that is
numerically confirmed in Fig. 3B).

Straight event-chain Monte Carlo (used in Fig. 3C) periodically aligns the
unit velocity $|\mathbf{v}_{a}|=1$ of the active atom with the coordinate
axes. We change the velocity after a chain time of $\tau_{\text{chain}}=0.2N$
(_52_). In an event triggered by a factor $M$, the velocity is transferred to
another atom which belongs to $M$. For factors with more than two atoms, the
next active atom is chosen according to the ratio lifting scheme (_27, 28_).

Newtonian event-chain Monte Carlo assigns a velocity label to every atom, but
only moves a single one at any time. The velocities are initialized so that
the average speed of the active atom during the simulation is approximately
one: $\langle|\mathbf{v}_{a}|\rangle\approx 1$ (see Supplementary Text for
details). This is possible because the velocities have no kinematic meaning as
they have in molecular dynamics. The physical temperature only rescales the
event rates [see eq. (5) in the Supplementary Text]. The velocities are
resampled in large time intervals of $\tau_{\text{chain}}=$10\,000${N}$. This
is because Newtonian event-chain Monte Carlo rotates dipoles most efficiently
in the limit $\tau_{\text{chain}}\rightarrow\infty$, while a finite value of
$\tau_{\text{chain}}$ ensures irreducibility (_53, 31_). Although our
generalized Newtonian event-chain Monte Carlo can consider a general mass
matrix of the atoms in its events, we set the masses of all atoms to be equal.
Then, the events of factors $M$ of two atoms can be treated in a deterministic
way in which the two involved velocities and the active atom always change.
For factors $M$ of more than two atoms, we apply a stochastic scheme that
corresponds to a generalized ratio lifting scheme (see the Newtonian-pair and
Newtonian-general schemes in the Supplementary Text).

The events of the bond factors are computed exactly [see eq. (6) in the
Supplementary Text], while the bending factors are treated with a piecewise-
linear bounding potential followed by a Poisson thinning procedure [see (_32,
Section 4.4.5_)]. The cell-veto algorithm (_34_) bundles most molecular
Coulomb factors and relies on a cell-occupancy system that tracks the
molecular barycenters [we generally adapt the setup described in Ref. (_32,
Section 5.3.4_)]. Very close molecule pairs or surplus molecules in the cell-
occupancy system are separately treated with a piecewise-linear bounding
potential. The cell bounds for $10$ different Fibonacci vectors (see
Supplementary Text for details) are estimated by varying the position of an
atom in one cell, and the position of a dipole in the other cell. Here, the
dipole is aligned with the direction of the gradient of the charge–dipole
Coulomb interaction. The maximum event rate multiplied by an empirical
prefactor yields the cell bound. During the simulation, the cell bounds in the
Walker table are rescaled to the actual speed and charge of the active atom.
Events from the cell-veto algorithm are confirmed with a Poisson thinning
procedure that compares the actual event rate of the Coulomb interaction
between two molecules with the sampled cell bound. To compute the real event
rate, we compute the gradients of the Coulomb potential between the nine pairs
of charged atoms [see also eq. (5) in the Supplementary Text] with an Ewald
summation (_50_) that is tuned to machine precision without any assumption on
the global charge distribution. We also use the cell-veto algorithm to treat
most Lennard-Jones factors, although we cut off the interaction beyond the
closest images for consistency with the other sampling algorithms. We
carefully checked that an alternative spherical cut-off does not change the
results of this paper. The cell-veto algorithm for the Lennard-Jones
interaction relies on a cell-occupancy system that tracks the oxygens. Very
close or surplus oxygen pairs are treated directly [see eq. (6) in the
Supplementary Text]. The cell bounds for $10$ different Fibonacci vectors are
estimated by varying the oxygen positions evenly in the two cells, and again
including an empirical prefactor. The parameters of the cell-occupancy
systems, as the cell sizes and the number of nearby cells that are excluded
from the cell-veto algorithm, are tuned to minimize computer time. As a rule
of thumb, the largest possible cells that do not yield any surplus water
molecules, and two excluded layers are found to be a good choice.

These simulation protocols for straight and Newtonian event-chain Monte Carlo
in JeLLyFysh are validated by comparing numerical results for $N=2$ SPC/Fw
water molecules of all sampling methods of this paper (see Fig. S1).

### Simulation Protocols—Creation of Initial Configurations

We create initial configurations in the range from $N=64$ to $N=2744$ SPC/Fw
water molecules with hydrogen mass $m_{H}=$1.0079\text{\,}\mathrm{Da}$$ and
oxygen mass $m_{O}=$15.999\,400\,1\text{\,}\mathrm{Da}$$ at a density of
$\rho=$0.97\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$ using the software
package Playmol (_54_) (commit 67eb56c from 26 November 2019). Initial
configurations are equilibrated using a molecular-dynamics simulation of
Lammps in the isothermal–isobaric ensemble at pressure
$P=$1\text{\,}\mathrm{atm}$$. Similar to thermostatting, barostatting is
achieved by integrating the appropriate Nosé–Hoover-chain equations of motion
with a Verlet integrator (_46, 47_) with a time step of
$1\text{\,}\mathrm{fs}$. This results in initial configurations for the
simulations in the canonical ensemble with the different sampling algorithms
at slightly different densities for different $N$. We confirmed that these
small density variations do not influence the results of this paper.

### Analysis of the Electric Polarization

We sample the electric polarization $\mathbf{P}$ (or the total electric dipole
moment $\mathbf{P}=\sum_{i=1}^{N}\mathbf{p}_{i}$, where $\mathbf{p}_{i}$ is
the molecular dipole moment of the water molecule $i$) of the SPC/Fw water
system with the different sampling algorithms. The generated “time”-series
$\mathbf{P}(t)$ (where the “time” $t$ only has a a physical meaning in
molecular dynamics) yields the normalized autocorrelation function

$\Phi_{\mathbf{P}}(\tau)=\frac{\langle\mathbf{P}(t)\cdot\mathbf{P}(t+\tau)\rangle_{t}}{\langle\mathbf{P}(t)\cdot\mathbf{P}(t)\rangle_{t}},$
(2)

where we explicitly assume that $\langle\mathbf{P}\rangle=0$ for
$t\rightarrow\infty$. The first part of each trajectory in the canonical
ensemble is not considered in the equilibrium average. We observe an
exponential decay in all sampling algorithms and fit an exponential
$\sim\exp(-\tau/\tau_{\Phi})$ to extract the time constant $\tau_{\Phi}$ and
its standard error [see also Ref. (_21_)]. We naïvely parallelize the
Metropolis and event-chain Monte Carlo algorithms by running
$20\text{\,}50\text{\,}$ independent simulations (where the number of runs
grows with $N$). For every $\tau$, we then compute the median autocorrelation
function $\Phi_{\mathbf{P}}(\tau)$ and estimate its error with the bootstrap
method. These errors are then considered in the least-squares fit.

The units of the autocorrelation time $\tau_{\Phi}$ for the different sampling
algorithms are physical time for molecular dynamics, Monte-Carlo trials for
the Metropolis algorithm, and continuous Monte-Carlo time for the event-chain
Monte Carlo algorithms. As a first measure of efficiency of the different
dynamics, we compare the autocorrelation distances $d_{\Phi}$ (see Fig. 3C).
It gives the average cumulative distance moved by the atoms in each of the
sampling schemes. Since only a single atom moves in our implementation of the
event-chain Monte Carlo algorithm, it follows that
$d_{\Phi}^{\text{ECMC}}=\langle|\mathbf{v}_{a}|\rangle^{\text{ECMC}}\,\tau_{\Phi}^{\text{ECMC}}$,
where $\langle|\mathbf{v}_{a}|\rangle^{\text{ECMC}}$ is the average speed of
the active atom during the simulation. For molecular dynamics, we can compute
the average speed of the hydrogen and oxygen atoms $t\in\\{\ch{H},\ch{O}\\}$
with the Maxwell–Boltzmann distribution:

$\langle\left|\mathbf{v}_{t}\right|\rangle^{\text{MD}}=\sqrt{\frac{8k_{\mathrm{B}}T}{\pi
m_{t}}},$ (3)

which yields
$d_{\Phi}^{\text{MD}}=N(2\langle|\mathbf{v}_{\ch{H}}|\rangle^{\text{MD}}+\langle|\mathbf{v}_{\ch{O}}|\rangle^{\text{MD}})\,\tau_{\Phi}^{\text{MD}}$.
In the Metropolis algorithm, each trial samples a random displacement vector
$\mathbf{d}_{t}=[\operatorname{ran}(-\delta_{t},\delta_{t}),\operatorname{ran}(-\delta_{t},\delta_{t}),\operatorname{ran}(-\delta_{t},\delta_{t})]^{T}$
for a random atom, where $\delta_{t}$ is tuned independently for the oxygen
and hydrogen atoms $t\in\\{\ch{H},\ch{O}\\}$ to accept the displacement with
probability $p\approx 0.37$. We can compute the average length
$\langle|\mathbf{d}_{t}|\rangle^{\text{Met}}$ of the sampled displacement
vector as

$\langle|\mathbf{d}_{t}|\rangle^{\text{Met}}=\frac{1}{8\delta_{t}^{3}}\int_{-\delta_{t}}^{\delta_{t}}\mathrm{d}x\,\int_{-\delta_{t}}^{\delta_{t}}\mathrm{d}y\,\int_{-\delta_{t}}^{\delta_{t}}\mathrm{d}z\,\sqrt{x^{2}+y^{2}+z^{2}}\approx$0.961$\,\delta_{t}.$
(4)

Neglecting correlations between rejection probability and the displacement for
this first measure of efficiency, we get
$d_{\Phi}^{\text{Met}}=p\,(2\langle|\mathbf{d}_{\ch{H}}|\rangle^{\text{Met}}/3+\langle|\mathbf{d}_{\ch{O}}|\rangle^{\text{Met}}/3)\,\tau_{\Phi}^{\text{Met}}$.

Finally, the computer autocorrelation time of our generalized Newtonian event-
chain Monte Carlo algorithm (see Fig. 4B) is obtained by combining the
autocorrelation distance (see Fig. 3C), the event rates (see Fig. 3B), and the
computer time per event in JeLLyFysh (see Fig. 3A). Similarly, we measure the
computer time per time step of molecular dynamics in Lammps for different
target accuracies of its particle–particle particle–mesh solver (see Fig. 4A).
This is then combined with the autocorrelation time in the number of time
steps.

## Supplementary Text

In this supplementary text, we generalize Newtonian event-chain Monte Carlo to
smooth factor potentials depending on an arbitrary number of atoms. We further
prove that this algorithm satisfies the necessary global-balance condition.
Finally, we report on our generation of Fibonacci vectors on the unit sphere.

### Generalized Newtonian Event-Chain Monte Carlo

Newtonian event-chain Monte Carlo was originally formulated for the hard-
sphere model that contains peculiar stepwise-changing two-body factor
potentials (_29_). Among various event-chain Monte Carlo variants, it was
shown to escape faster from sparse hard-disk packings (_53_) and to produce
favorable rotation dynamics in tethered hard-disk dipoles (_31_). Its
superiority is confirmed for the rotation dynamics of SPC/Fw water molecules
in this paper, but requires its generalization to smooth factor potentials
that depend on an arbitrary number of atom positions.

We follow the initial introduction of non-reversible event-chain Monte Carlo
as a lifted continuous-time Markov process (_20_) [see also the formulation as
a piecewise-deterministic Markov process in Ref. (_55_)]. For Newtonian event-
chain Monte Carlo for $N_{a}$ three-dimensional atoms, the lifting framework
(_56, 57_) [see also Refs (_24_) and (_58, Appendix A_)] extends the physical
all-atom configuration $\mathbf{x}$ to $(\mathbf{x},\mathbf{v},i)$, with
auxiliary lifting variables that represent an all-atom velocity $\mathbf{v}$
and an activity label $i$, respectively. The physical sample space
$\mathbf{x}\in\Omega$ is thus extended to the lifted sample space
$(\mathbf{x},\mathbf{v},i)\in\widehat{\Omega}=\Omega\times\mathcal{V}^{3N_{a}}(E_{\text{kin}})\times\mathcal{N}$,
where $\mathcal{N}=\\{1\ldots N_{a}\\}$ the set of atom indices, and
$\mathcal{V}^{3N_{a}}(E_{\text{kin}})=\\{\mathbf{v}\in\mathbb{R}^{3N_{a}}:\mathbf{v}^{T}\cdot\mathbf{M}\cdot\mathbf{v}=2E_{\text{kin}}\\}$
with a positive-definite symmetric mass matrix $\mathbf{M}$ and total kinetic
energy $E_{\text{kin}}$. The Markov process targets the lifted stationary
distribution
$\hat{\pi}(\mathbf{x},\mathbf{v},i)=\pi(\mathbf{x})\times\mu_{\mathcal{V}}(\mathbf{v})\times\mu_{\mathcal{N}}(i)$
that separates into the factorized Boltzmann distribution
$\pi(\mathbf{x})\propto\exp[-\beta U(\mathbf{x})]=\prod_{M}\exp[-\beta
U_{M}(\mathbf{x})]$, and the uniform distributions
$\mu_{\mathcal{V}}(\mathbf{v})$ on $\mathcal{V}^{3N_{a}}$, and
$\mu_{\mathcal{N}}(i)$ on $\mathcal{N}$.

Given a lifted configuration $(\mathbf{x}(t_{0}),\mathbf{v}(t_{0}),i)$ at time
$t_{0}$, Newtonian event-chain Monte Carlo continuously moves the single
active atom $i$ starting from its position $\mathbf{x}_{i}(t_{0})$ with its
constant velocity $\mathbf{v}_{i}=\mathbf{v}_{i}(t_{0})$ up to an event at
time $t_{\text{ev}}>t_{0}$, which interrupts the motion:
$\mathbf{x}_{i}(t)=\mathbf{x}_{i}(t_{0})+\mathbf{v}_{i}\,(t-t_{0})$ for
$t_{0}<=t<t_{\text{ev}}$. The velocities $\mathbf{v}_{j}$ of all other atoms
$j\neq i$ in the all-atom velocity $\mathbf{v}$ are, at this point,
hypothetical, that is, mere labels. The time-dependent event rate
$\lambda(t)$, that is, the probability density to interrupt the piecewise-
deterministic motion of the active atom $i$, is given by a sum of factor event
rates $\lambda_{M}(t)$:

$\lambda(t)=\sum_{M}\lambda_{M}(t)=\sum_{M}\beta\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}\left(\mathbf{x}(t)\right)\cdot\mathbf{v}_{i}\right].$
(5)

Every factor $M$ can be considered as statistically independent. Each of them
stochastically generates a candidate event time $t_{\text{ev},M}$ in an
inhomogeneous Poisson process based on its factor event rate
$\lambda_{M}(t)\geq 0$ (that is only nonzero for factor potentials $U_{M}$
that are actually changed by the motion of the active atom):

$\operatorname{ran}_{M}(0,1)=\exp\left[-\int_{t_{0}}^{t_{\text{ev},M}}\lambda_{M}(t)\,\mathrm{d}t\right].$
(6)

Here, $\operatorname{ran}_{M}(0,1)$ is a uniformly distributed random number
between $0$ and $1$ that is drawn separately for each factor $M$. In eq. (6),
the cumulative increments of the factor potential $U_{M}$ under the motion of
the active atom with its velocity $\mathbf{v}_{i}$ since $t_{0}$ are equal to
an exponentially distributed random number with mean $1/\beta$ at the time
$t_{\text{ev},M}$. If the computation of $t_{\text{ev},M}$ from eq. (6) is
tedious or even impossible, a Poisson thinning procedure may yield candidate
event times (_36_).

At the event time $t_{\text{ev}}=\min_{M}t_{\text{ev},M}$, the motion of the
active atom $i$ is interrupted at the lifted configuration
$(\mathbf{x}=\mathbf{x}(t_{\text{ev}}),\mathbf{v},i)$. An event changes the
lifting variables and sets the initial lifted configuration
$(\mathbf{x},\mathbf{v}^{\prime},j)$ for the next piece in the piecewise-
deterministic Markov process. The event is realized by a unique event-factor
$M_{\text{ev}}=\operatorname*{arg\,min}_{M}t_{\text{ev},M}$. Generalized
Newtonian event-chain Monte Carlo proposes an update of the all-atom velocity
$\mathbf{v}\rightarrow\mathbf{v}^{\prime}$ at time $t_{\text{ev}}$ with a
force kick in the direction of the gradient of the event-factor potential:

$\mathbf{v}^{\prime}=\mathbf{v}-2\,\frac{\mathbf{v}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M_{\text{ev}}}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M_{\text{ev}}})^{T}\cdot\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M_{\text{ev}}}}\,\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M_{\text{ev}}}.$
(7)

This force kick leaves the kinetic energy
$2E_{\text{kin}}=\mathbf{v}^{T}\cdot\mathbf{M}\cdot\mathbf{v}$ invariant and
applying it twice yields the original $\mathbf{v}$ (_59_). The event-factor
potential only acts on a small subensemble $k\in M_{\text{ev}}$ of atoms, and
thus only the velocities $\mathbf{v}_{k}$ of these contributing atoms are
possibly modified. For the simple choice of an identity mass matrix
$\mathbf{M}=\mathbf{I}$, eq. (7) may be written as

$\mathbf{v}_{k}^{\prime}=\mathbf{v}_{k}-2\,\frac{\sum_{j\in
M_{\text{ev}}}\mathbf{v}_{j}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M_{\text{ev}}}}{\sum_{j\in
M_{\text{ev}}}\left|\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M_{\text{ev}}}\right|^{2}}\,\boldsymbol{\nabla}_{\mathbf{x}_{k}}U_{M_{\text{ev}}}.$
(8)

In Newtonian event-chain Monte Carlo, only a single active atom moves with its
velocity at any time. At an event, it must choose the next active atom $k\in
M_{\text{ev}}$ from the event factor. It also chooses to either change the
velocity $\mathbf{v}\rightarrow\mathbf{v}^{\prime}$ of all contributing atoms
according to eq. (7), or to keep the all-atom velocity constant
$\mathbf{v}\rightarrow\mathbf{v}$. Generalized Newtonian event-chain Monte
Carlo offers two schemes to treat events. The _Newtonian-general_ scheme
applies to general mass matrices $\mathbf{M}$ and factors that depend on an
arbitrary number of atoms. The _Newtonian-pair_ scheme only applies to
distance-dependent pair potentials and identity mass matrix
$\mathbf{M}=\mathbf{I}$.

In the Newtonian-general scheme, the probability for each possible choice
$(k,\mathbf{v}^{\prime})$ or $(k,\mathbf{v})$ of the lifting variables is
given by
$\max[0,-\boldsymbol{\nabla}_{\mathbf{x}_{k}}U_{M_{\text{ev}}}\cdot\mathbf{v}_{k}^{\prime}]/C$
or
$\max[0,-\boldsymbol{\nabla}_{\mathbf{x}_{k}}U_{M_{\text{ev}}}\cdot\mathbf{v}_{k}]/C$,
respectively. Here, $C$ is a common normalization factor. If the force kick is
not applied and $\mathbf{v}$ stays constant, the active atom always changes
because the factor event rate was positive at the time of the event [see eqs
(5) and (6)]. Likewise, if the event-factor potential only depends on a single
atom (as for external potentials), the force kick is always applied.

In the Newtonian-pair scheme for distance-dependent pair potentials between
two atoms $a$ and $b$,
$U_{M_{\text{ev}}}=U_{M_{\text{ev}}}(|\mathbf{x}_{a}-\mathbf{x}_{b}|)$ and an
identity mass matrix $\mathbf{M}=\mathbf{I}$, we alternatively exploit the
inherent translational invariance to treat the event deterministically. We
always apply the force kick in eq. (8) to modify both $\mathbf{v}_{a}$ and
$\mathbf{v}_{b}$, and always change the active atom within the event factor.
In this case, eq. (8) is equivalent to an elastic Newtonian collision between
two particles of equal mass, and this alternative scheme recovers the original
formulation of Newtonian event-chain Monte Carlo for the hard-sphere model
(_29_).

Samples of all-atom configurations $\mathbf{x}$ that can be used to compute
observables are taken at periodic time intervals (_29_). Furthermore,
Newtonian event-chain Monte Carlo requires resamplings of the all-atom
velocity $\mathbf{v}$ and the active atom $i$ in periodic time intervals
$\tau_{\text{chain}}$ to become irreducible (_31_). We sample $\mathbf{v}$ and
$i$ from their respective stationary distributions
$\mu_{\mathcal{V}}(\mathbf{v})$ and $\mu_{\mathcal{N}}(i)$. For $\mathbf{v}$,
this is equivalent to sampling a multivariate Gaussian distribution with mean
$\boldsymbol{\mu}=\mathbf{0}$ and covariance matrix
$\boldsymbol{\Sigma}=\mathbf{M}^{-1}$, followed by corrections of the total
velocity $\mathbf{1}^{T}\cdot\mathbf{v}\rightarrow 0$ and
$\mathbf{v}^{T}\cdot\mathbf{M}\cdot\mathbf{v}\rightarrow 2E_{\text{kin}}$. The
activity label $i$ is sampled uniformly from the the set
$\mathcal{N}=\\{1\ldots N_{a}\\}$.

### Proof of Correctness of Generalized Newtonian Event-Chain Monte Carlo

In order to prove the correctness of the Newtonian event-chain Monte Carlo
algorithm for general smooth interactions, we must show that the global
balance condition is satisfied. The total probability flow into any lifted
configuration $(\mathbf{x},\mathbf{v},i)\in\widehat{\Omega}$ consists of a
physical flow $\mathcal{F}^{\text{phys}}(\mathbf{x},\mathbf{v},i)$ and a
lifting flow $\mathcal{F}^{\text{lift}}(\mathbf{x},\mathbf{v},i)$ and must
equal its statistical weight $\hat{\pi}(\mathbf{x},\mathbf{v},i)$:

$\hat{\pi}(\mathbf{x},\mathbf{v},i)=\mathcal{F}^{\text{phys}}(\mathbf{x},\mathbf{v},i)+\mathcal{F}^{\text{lift}}(\mathbf{x},\mathbf{v},i).$
(9)

The physical flow into $(\mathbf{x},\mathbf{v},i)$ stems from the continuous
movement of the active atom $i$ that was not interrupted by an event with the
event rate given in eq. (5). With
$\mathbf{x}^{\prime}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{i}-\mathbf{v}_{i}\mathrm{d}t,\ldots,\mathbf{x}_{N_{a}})$,
we get

$\begin{split}\mathcal{F}^{\text{phys}}(\mathbf{x},\mathbf{v},i)&=\hat{\pi}(\mathbf{x}^{\prime},\mathbf{v},i)\left\\{1-\beta\sum_{M}\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}^{\prime}_{i}}U_{M}\cdot\mathbf{v}_{i}\right]\right\\}\\\
&=\hat{\pi}(\mathbf{x},\mathbf{v},i)\left\\{1-\beta\sum_{M}\max\left[0,-\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}\cdot\mathbf{v}_{i}\right]\right\\}\\\
&\eqqcolon\hat{\pi}(\mathbf{x},\mathbf{v},i)+\sum_{M}\mathcal{F}^{\text{phys}}_{M}(\mathbf{x},\mathbf{v},i),\end{split}$
(10)

where the second line uses the detailed-balance property of the factorized
Metropolis filter that yields the event rates in eq. (5) and allows to treat
all factors as statistically independent (_20_). The lifting flow into
$(\mathbf{x},\mathbf{v},i)$ stems from interrupted motions of lifted
configurations $(\mathbf{x},\mathbf{v}^{\prime},j)$ with different lifting
variables $\mathbf{v}^{\prime}$ and $j$ but the same configuration
$\mathbf{x}$:

$\begin{split}\mathcal{F}^{\text{lift}}(\mathbf{x},\mathbf{v},i)&=\beta\sum_{M}\sum_{j\in\mathcal{N}}\int_{\mathcal{V}^{3N_{a}}(E_{\text{kin}})}\mathrm{d}^{3N_{a}}\mathbf{v}^{\prime}\,\hat{\pi}(\mathbf{x},\mathbf{v}^{\prime},j)\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}\cdot\mathbf{v}^{\prime}_{j}\right]p^{M}_{(\mathbf{v}^{\prime},j),(\mathbf{v},i)}\\\
&\eqqcolon\sum_{M}\mathcal{F}^{\text{lift}}_{M}(\mathbf{x},\mathbf{v},i).\end{split}$
(11)

Here, $p_{(\mathbf{v}^{\prime},j),(\mathbf{v},i)}$ is the probability to
change the lifting variables $(\mathbf{v}^{\prime},j)$ to $(\mathbf{v},i)$.
For every factor $M$, the lifting flow
$\mathcal{F}^{\text{lift}}_{M}(\mathbf{x},\mathbf{v},i)$ in eq. (11) exactly
cancels the physical flow
$\mathcal{F}^{\text{phys}}_{M}(\mathbf{x},\mathbf{v},i)$ in eq. (10) so that
the global-balance condition in eq. (9) is satisfied. A trivial (inefficient)
solution would be to simply invert the velocity of the active atom $i$ in an
event while keeping it active:
$p^{M}_{(\mathbf{v}^{\prime},j),(\mathbf{v},i)}=\delta_{ij}\,\delta^{(3)}(\mathbf{v}_{i}-\mathbf{v}^{\prime}_{i})$.

We presented two schemes for the update of the lifting variables in an event
of generalized Newtonian event-chain Monte Carlo. The first stochastic
Newtonian-general scheme applies for a general number of atoms on the event
factor. The second deterministic Newtonian-pair scheme exploits the
translational invariance of distance-dependent pair potentials. For the
Newtonian-general scheme, we get

$\begin{split}p^{M}_{(\mathbf{v}^{\prime},j),(\mathbf{v},i)}=&\frac{\max\left[0,-\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}\cdot\mathbf{v}_{i}\right]}{C}\bigg{[}\\\
&\quad\delta^{(3N_{a})}\left(\mathbf{v}-\mathbf{v}^{\prime}+2\,\frac{(\mathbf{v}^{\prime})^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\,\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}\right)\\\
&\quad+\delta^{(3N_{a})}\left(\mathbf{v}-\mathbf{v}^{\prime}\right)\bigg{]}\\\
=&\frac{\max\left[0,-\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}\cdot\mathbf{v}_{i}\right]}{C}\bigg{[}\\\
&\quad\delta^{(3N_{a})}\left(\mathbf{v}^{\prime}-\mathbf{v}+2\,\frac{\mathbf{v}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\,\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}\right)\\\
&\quad+\delta^{(3N_{a})}\left(\mathbf{v}^{\prime}-\mathbf{v}\right)\bigg{]},\end{split}$
(12)

where the first term applies the force kick, while the second does not. From
the necessary condition

$1=\sum_{k\in\mathcal{N}}\int_{\mathcal{V}^{3N_{a}}(E_{\text{kin}})}\mathrm{d}^{3N_{a}}\mathbf{v}\,p^{M}_{(\mathbf{v}^{\prime},j),(\mathbf{v},k)},$
(13)

it follows that

$\begin{split}C&=\sum_{k\in\mathcal{N}}\bigg{\\{}\max\left[0,-\boldsymbol{\nabla}_{\mathbf{x}_{k}}U_{M}\cdot\left(\mathbf{v}^{\prime}_{k}-2\frac{(\mathbf{v}^{\prime})^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\left(\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}\right)_{k}\right)\right]\\\
&\quad\quad\quad\,\,\,\,+\max\left[0,-\boldsymbol{\nabla}_{\mathbf{x}_{k}}U_{M}\cdot\mathbf{v}_{k}^{\prime}\right]\bigg{\\}}\\\
&=\sum_{k\in\mathcal{N}}\bigg{\\{}\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{k}}U_{M}\cdot\left(\mathbf{v}^{\prime}_{k}-2\frac{(\mathbf{v}^{\prime})^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\left(\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}\right)_{k}\right)\right]\\\
&\quad\quad\quad\,\,\,\,+\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{k}}U_{M}\cdot\mathbf{v}_{k}^{\prime}\right]\bigg{\\}}.\end{split}$
(14)

Here, the second rewriting follows from expressing the trivial equality

$\mathbf{v}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}+\left(\mathbf{v}-2\frac{\mathbf{v}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}\right)^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}=0$
(15)

as a sum over $k\in\mathcal{N}$ and from grouping the positive and negative
terms separately. The normalization factor $C$ in
$p^{M}_{(\mathbf{v}^{\prime},j),(\mathbf{v},i)}$ is the same for every
possible next active atom $i$. Because of the $\delta^{(3N_{a})}$-functions in
eq. (12), the normalization factor $C$ can also be written in terms of the
velocity $\mathbf{v}$ after the force kick by replacing
$\mathbf{v}^{\prime}\rightarrow\mathbf{v}$ in eq. (14). This yields

$\begin{split}\mathcal{F}^{\text{lift}}_{M}(\mathbf{x},\mathbf{v},i)&=\beta\,\hat{\pi}(\mathbf{x},\mathbf{v},i)\max\left[0,-\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}\cdot\mathbf{v}_{i}\right]\frac{1}{C}\sum_{j\in\mathcal{N}}\bigg{\\{}\\\
&\hphantom{=}\quad\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}\cdot\left(\mathbf{v}_{j}-2\frac{\mathbf{v}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\left(\mathbf{M}^{-1}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}\right)_{j}\right)\right]\\\
&\hphantom{=}\quad+\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}\cdot\mathbf{v}_{j}\right]\bigg{\\}}\\\
&=\beta\,\hat{\pi}(\mathbf{x},\mathbf{v},i)\max\left[0,-\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}\cdot\mathbf{v}_{i}\right]\\\
&=-\mathcal{F}^{\text{phys}}_{M}(\mathbf{x},\mathbf{v},i),\end{split}$ (16)

where we used that the velocity $\mathbf{v}$ and activity label $i$ are
uniformly distributed over their sample spaces
$\mathcal{V}^{3N_{a}}(E_{\text{kin}})$ and $\mathcal{N}$, respectively. Eq.
(16) concludes the proof that the presented Newtonian-general scheme for the
update of the lifting variables in an event by a general factor potential
satisfies the global-balance condition of eq. (9). The reflective (_23_) and
forward (_30_) event-chain Monte Carlo variants can likewise be generalized to
factor potentials depending on an arbitrary number of atoms.

For the Newtonian-pair scheme for distance-dependent pair potentials and
identity mass matrix $\mathbf{M}=\mathbf{I}$, we get

$\begin{split}p^{M}_{(\mathbf{v}^{\prime},j),(\mathbf{v},i)}&=\begin{cases}\delta^{(3N_{a})}\left(\mathbf{v}-\mathbf{v}^{\prime}+2\,\frac{(\mathbf{v}^{\prime})^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\,\boldsymbol{\nabla}_{\mathbf{x}}U_{M}\right)&\text{if
$i,j\in M$ and $i\neq j$,}\\\ 0&\text{otherwise,}\end{cases}\\\
&=\begin{cases}\delta^{(3N_{a})}\left(\mathbf{v}^{\prime}-\mathbf{v}+2\,\frac{\mathbf{v}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\,\boldsymbol{\nabla}_{\mathbf{x}}U_{M}\right)&\text{if
$i,j\in M$ and $i\neq j$,}\\\ 0&\text{otherwise.}\end{cases}\end{split}$ (17)

With $i,j\in M$ and $i\neq j$, this yields

$\begin{split}\mathcal{F}^{\text{lift}}_{M}(\mathbf{x},\mathbf{v},i)&=\beta\,\hat{\pi}(\mathbf{x},\mathbf{v},i)\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}\cdot\left(\mathbf{v}_{j}-2\frac{\mathbf{v}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}{(\boldsymbol{\nabla}_{\mathbf{x}}U_{M})^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}}U_{M}}\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}\right)\right]\\\
&=\beta\,\hat{\pi}(\mathbf{x},\mathbf{v},i)\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}\cdot\mathbf{v}_{j}-2\frac{\mathbf{v}_{i}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}+\mathbf{v}_{j}^{T}\cdot\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}}{\left|\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}\right|^{2}+\left|\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}\right|^{2}}\left|\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}\right|^{2}\right]\\\
&=\beta\,\hat{\pi}(\mathbf{x},\mathbf{v},i)\max\left[0,-\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}\cdot\mathbf{v}_{i}\right]\\\
&=-\mathcal{F}^{\text{phys}}_{M}(\mathbf{x},\mathbf{v},i).\end{split}$ (18)

Here, we again used the uniform distributions of $\mathbf{v}$ and $i$ and the
translational invariance
$\boldsymbol{\nabla}_{\mathbf{x}_{i}}U_{M}=-\boldsymbol{\nabla}_{\mathbf{x}_{j}}U_{M}$.
Thus, also the Newtonian-pair scheme satisfies the global-balance condition of
eq. (9). We note that it always applies the force kick.

### Fibonacci Vectors

The cell-veto algorithm was previously only used for the straight variant of
event-chain Monte Carlo (_34, 28, 32_). Its finite set of possible velocities
$\mathcal{D}$ of the single active atom $a$ allows to set up a Walker table
for every velocity $\mathbf{v}_{a}\in\mathcal{D}$. For instance, for factor
pair potentials between the active atom $a$ and another atom $b$ at positions
$\mathbf{x}_{a}$ and $\mathbf{x}_{b}$,
$U_{M}(\mathbf{x})=U_{M}(\mathbf{x}_{a},\mathbf{x}_{b})$, every cell bound
$q_{M}^{\text{cell}}(\mathcal{C}_{a},\mathcal{C}_{b},\mathbf{v}_{a})$ for the
pair of cells $\mathcal{C}_{a}$ and $\mathcal{C}_{b}$ in the Walker table may
be written as

$q_{M}^{\text{cell}}(\mathcal{C}_{a},\mathcal{C}_{b},\mathbf{v}_{a})=\max_{\mathbf{x}_{a}\in\mathcal{C}_{a},\mathbf{x}_{b}\in\mathcal{C}_{b}}\beta\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{a}}U_{M}(\mathbf{x}_{a},\mathbf{x}_{b})\cdot\mathbf{v}_{a}\right].$
(19)

In principle, a Walker table must be precomputed for every velocity
$\mathbf{v}_{a}\in\mathcal{D}$ and for every possible cell $\mathcal{C}_{a}$
of the active atom. Symmetries, such as a translational invariance of the
factor potential, may heavily reduce the necessary number of Walker tables
(_28_). During the straight event-chain Monte Carlo simulation, the relevant
Walker table is determined by the cell of the currently active atom and its
velocity. The generalization of eq. (19) to more complex factors is
straightforward. For example, in order to treat the Coulomb interaction
between two water molecules in the SPC/Fw water model, the position and
orientation of the water molecule containing the active atom may be varied in
$\mathcal{C}_{a}$, while the position and orientation of the other molecule is
varied in $\mathcal{C}_{b}$. From a practical point of view, the cell bounds
$q_{M}^{\text{cell}}(\mathcal{C}_{a},\mathcal{C}_{b})$ are usually not
computed exactly but rather approximated, e.g., by considering a finite set of
positions $\mathbf{x}_{a}\in\mathcal{C}_{a}$ and
$\mathbf{x}_{b}\in\mathcal{C}_{b}$, and including a corrective multiplicative
prefactor (_32_). As long as the approximated cell bound satisfies
$\tilde{q}^{\text{cell}}_{M}(\mathcal{C}_{a},\mathcal{C}_{b},\mathbf{v}_{a})\geq
q^{\text{cell}}_{M}(\mathcal{C}_{a},\mathcal{C}_{b},\mathbf{v}_{a})$, a
Poisson thinning procedure (_36_) corrects _any_ overestimate. The quality of
$\tilde{q}^{\text{cell}}_{M}(\mathcal{C}_{a},\mathcal{C}_{b},\mathbf{v}_{a})$,
however, does influence the performance because higher cell bounds in the
cell-veto algorithm yield more events per unit distance that have to be
confirmed by computing the actual event rate.

The inherent discretization of continuous space in eq. (19) can be translated
to a continuous velocity space, as it appears, e.g., in Newtonian event-chain
Monte Carlo. Consider a finite number of unit vectors $\hat{\mathbf{v}}_{d}$
on the two-dimensional unit sphere, and let $\mathcal{V}_{d}$ be the
associated Voronoi cells under some distance function. We can then formally
compute the cell bounds for every Voronoi cell as

$q_{M}^{\text{cell}}(\mathcal{C}_{a},\mathcal{C}_{b},\mathcal{V}_{a})=\max_{\mathbf{x}_{a}\in\mathcal{C}_{a},\mathbf{x}_{b}\in\mathcal{C}_{b},\mathbf{v}_{a}\in\mathcal{V}_{a}}\beta\max\left[0,\boldsymbol{\nabla}_{\mathbf{x}_{a}}U_{M}(\mathbf{x}_{a},\mathbf{x}_{b})\cdot\mathbf{v}_{a}\right].$
(20)

The results are used to precompute Walker tables for every Voronoi cell
$\mathcal{V}_{a}$ and for every cell $\mathcal{C}_{a}$, where symmetries may
again heavily reduce the number of actually necessary tables. During the
Newtonian event-chain Monte Carlo simulation, the relevant Walker table is
determined by finding the cell of the currently active atom and the Voronoi
cell of its normalized velocity $\mathbf{v}_{a}/|\mathbf{v}_{a}|$. We can then
correct all cell bounds in the relevant Walker table to the actual speed of
the active atom by multiplying them with $|\mathbf{v}_{a}|$.

In this paper, as proof of concept, we map a generalized Fibonacci lattice
onto a two-dimensional unit sphere by the Lambert cylindrical equal-area
projection to generate $D$ unit vectors $\hat{\mathbf{v}}_{d}^{\text{fib}}$:

$\hat{\mathbf{v}}_{d}^{\text{fib}}=\left(\cos\phi_{d}\sin\theta_{d},\sin\phi_{d}\sin\theta_{d},\cos\theta_{d}\right)^{T},$
(21)

where $0\leq d<D$ and

$\begin{split}\phi_{d}&=\frac{2\pi d}{\varphi},\\\
\theta_{d}&=\arccos\left(1-\frac{2(d+\varepsilon)}{D-1+2\varepsilon}\right),\end{split}$
(22)

with the golden ratio $\varphi=(1+\sqrt{5})/2$ (_60, 61, 62_). The empirical
choice of the parameter $\varepsilon=$0.36$$ optimizes the average nearest-
neighbor distance of the Fibonacci vectors $\hat{\mathbf{v}}_{d}$. As a
distance function to construct the Voronoi cells $\mathcal{V}_{d}$, we use the
quick-to-evaluate cosine distance. Eqs (21) and (22) efficiently generate an
arbitrary number of Fibonacci vectors with a nearly uniform distribution.
However, there is no efficient inverse mapping from a general velocity
$\mathbf{v}_{a}$ to the closest Fibonacci vector $\hat{\mathbf{v}}_{d}$. This
is not a problem in this paper because we choose $D=10$ small. Then,
$\hat{\mathbf{v}}_{d}$ can be found by brute force. If considerably larger
values of $D$ become necessary, other point configurations on the two-
dimensional unit sphere may be considered (_62, 63_). For the uniform SPC/Fw
water systems in this paper, the Walker tables do not strongly depend on the
respective Fibonacci vector. The proper discretization of velocity space will,
however, impact the performance of Newtonian event-chain Monte Carlo in
nonuniform systems.

Figure S1: Validation for the SPC/Fw water model. Cumulative distribution
function of the absolute value of the electric polarization of $N=2$ SPC/Fw
water molecules in a cubic box of side length $L=$20\text{\,}\mathrm{\AA}$$.
The simulation protocols of the different sampling methods are described in
the Materials and Methods.