#  Enhanced path sampling using subtrajectory Monte Carlo moves

Daniel T. Zhang  Norwegian University of Science and Technology, Department of
Chemistry, NO-7491 Trondheim, Norway  Enrico Riccardi  Department of
Informatics, UiO, Gaustadalléen 23B, 0373 Oslo, Norway  Titus S. van Erp
Norwegian University of Science and Technology, Department of Chemistry,
NO-7491 Trondheim, Norway

(20 September 2022)

###### Abstract

Path sampling allows the study of rare events like chemical reactions,
nucleation and protein folding via a Monte Carlo (MC) exploration in path
space. Instead of configuration points, this method samples short molecular
dynamics (MD) trajectories with specific start- and end-conditions. As in
configuration MC, its efficiency highly depends on the types of MC moves.
Since the last two decades, the central MC move for path sampling has been the
so-called shooting move in which a perturbed phase point of the old path is
propagated backward and forward in time to generate a new path. Recently, we
proposed the subtrajectory moves, stone-skipping (SS) and web-throwing (WT),
that are demonstrably more efficient. However, the one-step crossing
requirement makes them somewhat more difficult to implement in combination
with external MD programs or when the order parameter determination is
expensive. In this article, we present strategies to address the issue. The
most generic solution is a new member of subtrajectory moves, wire fencing
(WF), that is less thrifty than the SS, but more versatile. This makes it
easier to link path sampling codes with external MD packages and provides a
practical solution for cases where the calculation of the order parameter is
expensive or not a simple function of geometry. We demonstrate the WF move in
a double well Langevin model, a thin film breaking transition based on
classical force fields, and a smaller ruthenium redox reaction at the ab
initio level in which the order parameter explicitly depends on the electron
density.

##  I Introduction

Rare event simulation techniques aim to sample events that require an
exceedingly long CPU/wall time to be simulated with standard molecular
dynamics (MD). In classical full atom simulations of protein folding, for
example, the longest reported Lindorff-Larsen  _et al._ (2011) MD runs
generated by the special-purpose molecular dynamics Anton 1 supercomputer are
around 1 ms, allowing the study of fast-folding proteins. The most recently
released Anton 3 supercomputer is even able to generate 100 $\mu$s/day for a
million atom system Shaw  _et al._ (2021). Despite this a remarkable speed, it
is still not fast enough to study the folding of all proteins. For instance,
the tryptophan synthase $\beta_{2}$ subunit has an experimentally measured
Goldberg  _et al._ (1990) folding rate of $k_{f}=0.001$ s-1. Hence, the
protein needs on average 1000 seconds to fold. The Anton 3 computer would thus
need 27,379 wall time years to generate one single transition. For ab initio
MD (AIMD) the situation is even worse as the quantum mechanical force
evaluation is orders of magnitude slower than computing the gradient of a
classical force field potential. In addition, no special purpose AIMD
computers exist today.

Yet, rare event simulations allow the calculation of rate constants and the
study reaction mechanisms orders of magnitude faster than MD, oftentimes
without sacrificing any molecular-level resolution. Peters (2017) (Replica
exchange) transition interface sampling (RE)TIS van Erp  _et al._ (2003); van
Erp (2007) is such a method that exploits the idea of transition path sampling
(TPS) Dellago  _et al._ (1998a) to focus the CPU time on the actual barrier
crossing event via a Monte Carlo (MC) sampling of MD paths.

RETIS and TIS employ a series of path sampling simulations, each sampling a
different path ensemble. The path ensembles differ with respect to a minimal
progress requirement, i.e. the number of interfaces (defined by fixed values
of the reaction coordinate/order parameter) that have to be crossed.Cabriolu
_et al._ (2017) Combining the results of all path ensembles allows the
computation of rate constants and other properties with an exponentially
reduced CPU time compared to a single MD simulation.

For instance, a classical simulation study on methane hydrate formation Arjun
and Bolhuis (2020) using TIS and RETIS reports on a crystallisation rate of
$10^{-17}$ nuclei per second per simulation volume. In other words, in a
system as small as those used in atomistic simulations, the process for
forming a single critical nucleus takes physically 3 years. Naturally, the
hypothetical wall time for reaching this with MD is astronomical for any
supercomputer. Likewise, RETIS simulations Moqadam  _et al._ (2018) reproduced
rate constant of water dissociation at the AIMD level in reasonable agreement
with experiments suggesting it happens once per 11 hours for each water
molecule. Eigen and de Maeyer (1958); Natzle and Moore (1985) As it required
30 minutes to produce 1 ps MD time in the 32 water molecules system, a naive
straightforward AIMD approach would need 0.7 billion centuries wall time to
generate a single dissociation.

Despite being orders of magnitude times faster than plain MD, simulations like
the above are still computationally expensive and can require months to years
to obtain satisfactory statistical accuracy. A further increase in efficiency
is therefore desirable. There are essentially three approaches to achieve
this: i) reducing the cost of the MC moves, ii) reducing the number of
required trajectories, and iii) parallelization of the algorithms. Partial
path sampling (PPTIS) Moroni  _et al._ (2003) and milestoning Faradjian and
Elber (2004) can be viewed as realisations of i) by sampling more restrictive
path ensembles with a reduced average path length. Unfortunately, this
introduces additional approximations. Strategies ii) and iii), on the other
hand, allow for a speed-up while still producing exact results, identical to
those from hypothetical unattainably long MD simulations. In fact, RETIS
successfully employs strategy ii) by complementing the shooting moves with
replica exchange moves between path ensembles. RETIS is thus more CPU
efficient compared to TIS. However, regarding strategy iii), TIS has the
advantage that path ensembles can be run in parallel completely independently,
while replica exchange moves require the progress of the sampling in the path
ensembles to be synchronized such that processing units do not have to wait
for each other. As a result, RETIS might not always outperform TIS based on
wall time, which is the reason why previously mentioned hydrate formation
study was partly based on TIS. Arjun and Bolhuis (2020) The recently
introduced $\infty$RETIS algorithm Roet  _et al._ (2022) is expected to solve
this issue for future studies based on a fundamentally new replica exchange
technique for cost-unbalanced replicas.

In fact, $\infty$RETIS implicitly applies the cost-free replica exchange moves
an infinite number of times after each shooting move. Still, replica exchange
moves alone are not ergodic and should, therefore, only be used in combination
with another MC move like shooting. Dellago  _et al._ (1998b) To further push
strategy ii), the principle MC move should be changed in order to reduce both
the rejection rate and the resemblance between accepted paths. This is exactly
what subtrajectory moves aim to establish. These MC moves resemble PPTIS
Moroni  _et al._ (2003) and milestoning Faradjian and Elber (2004) in the
sense that they evolve via series of shorter paths (subtrajectories/subpaths),
but differently to those methods, these subpaths are just intermediates
between sampled paths that are extended to their full lengths. Sampled paths,
therefore, have no configuration point in common with the previous path and
the statistical inefficiency is typically reduced by a factor equal to the
number of intermediate subtrajectories. So while the creation of a full new
path becomes more expensive, this is more than offset by the fact that far
fewer trajectories are needed to achieve a certain statistical accuracy. In
addition, the approach can be combined with a high-acceptance protocol, which
minimizes the number of rejections. As a result, most path ensembles obtain a
nearly 100% acceptance.Riccardi  _et al._ (2017)

The two moves presented in Ref. Riccardi  _et al._ , 2017, stone skipping (SS)
and web throwing (WT), however, have one element, the one-step crossing
condition, which can hinder the practical implementation with external MD
programs or when the calculation of the order parameter is computationally
expensive. In SS and WT, the subtrajectories are launched from a configuration
point of a previous (sub) path that is just before or after the path
ensemble’s interface. At this configuration point, velocities are generated
such that the interface is crossed again within a single time step. The
velocity randomization and one-step crossing test is reiterated several times
until the condition is fulfilled. The procedure is based on the idea that
generation of new random velocities followed by a one-step crossing test is
relatively cheap compared to generating MD steps, especially if the test can
be performed without new force calculations. This might not always be the
case. Present path sampling codes Swenson  _et al._ (2019a, b); Lervik  _et
al._ (2017); Riccardi  _et al._ (2020) use external MD codes for performing
the MD steps. PyRETIS version 2 has for instance couplings to Gromacs,Abraham
_et al._ (2015a) Lammps,Plimpton (1995) openmmEastman  _et al._ (2013), and
CP2KHutter  _et al._ (2014). In order to reduce the number of stop/restart
calls to these programs, a “time step” in the RETIS program is often several
(10-1000) MD steps by the external MD engine. This complicates the one-step
crossing condition as it actually involves not one, but several steps which is
costly and not easy to predict without actually performing these steps.
Another issue arises when the calculation of the order parameter is expensive
such as the those used in nucleation studies. Stukowski (2012); Winczewski
_et al._ (2016)

In this article, we discuss several approaches to tackle this issue. The most
generically applicable solution is a new member of the subtrajectory family
called wire fencing (WF). The approach is slightly more wasteful with respect
to the number of MD steps compared to SS, but very versatile and does not
require any code modifications of the external engines. We illustrate the WF
move on three model systems, a simple 1D double well potential, a Gromacs thin
film breakage application, and a CP2K study on ruthenium redox reactions.

##  II subtrajectory moves

The schematic main idea of the three subtrajectory moves is shown in Fig. 1.

Figure 1: Cartoon representation of the three subtrajectory moves: stone
skipping, web throwing, and wire fencing. The old path is shown in blue. Four
subtrajectories are shown in orange. The final new path consists of the fourth
subtrajectory and its extensions colored in green.

These are the stone skipping (SS), web throwing (WT), and the new wire fencing
(WF) move. The commonality is that an arbitrary number of partial trajectories
(subtrajectories/subpaths) are generated before the completion of a new full
trajectory. The subtrajectories obey different start- and end-conditions and
are, due to this, considerably shorter than full trajectories. The
subtrajectories are not part of the sampling, but just intermediate steps
between one full trajectory to another. The $[i^{+}]$ path ensemble that is
being sampled in Fig. 1 consists of paths starting at $\lambda_{A}$, crossing
$\lambda_{i}$ at least once, and ending at either $\lambda_{A}$ or
$\lambda_{B}$. In Fig. 1, the old full trajectory is colored blue. In the
example, the new trajectory is generated via 4 subtrajectories. The first
subtrajectory is obtained from a shooting move from the old trajectory. Then,
a next subtrajectory is generated from the previous one until the number of
predetermined subtrajectories (4 in this case, colored in orange) is reached.
The final subtrajectory is extended backward and forward in time until
reaching a stable state. The new full trajectory comprises the last
subtrajectory and the extensions colored in green. The difference between the
three moves lies in the way the shooting of subpaths is executed.

The SS move resembles a flat stone that collides with the water’s surface
after a skilful throw. The move starts by selecting randomly any of the
crossing points of the old path with $\lambda_{i}$, generates new velocities
that also establish a crossing, and then proceeds until $\lambda_{B}$ is
crossed or $\lambda_{i}$ is crossed again. The process is then repeated by
selecting the subpath’s last crossing with $\lambda_{i}$ for shooting off the
next subpath. Finally, the last subpath is extended and possibly accepted or
rejected.Riccardi  _et al._ (2017)

The WT move has been named after a gesture of the famous Marvel character
swinging between skyscrapers. Here, an additional interface needs to be
defined, _the surface of unlikely return_ (SOUR), at the state $A$ side of the
$\lambda_{i}$ interface. If this interface, $\lambda_{\rm sour}$, is crossed
towards the direction of state $A$, it is assumed to be highly unlikely that
the MD trajectory will end up in state $B$ rather than $A$ (defined by the
last interface, $\lambda_{B}$, and the first interface, $\lambda_{A}$,
respectively). The first subpath is then shot from a random crossing point
with either $\lambda_{i}$ or $\lambda_{\rm sour}$ at a path segment of the old
path that connects these two interfaces. After the velocities of the system’s
atoms are re-set, like in the SS move, the subpath is continued till
$\lambda_{\rm sour}$ or $\lambda_{i}$ is crossed, but is only kept if the
subpath connects $\lambda_{\rm sour}$ and $\lambda_{i}$ again like the segment
of the old path. If not, the subpath is rejected and a new crossing point is
taken randomly from the same segment. If both $\lambda_{\rm sour}$ and
$\lambda_{i}$ are crossed, the subpath replaces the segment. The process is
repeated until the selected number of subpaths, accepted or rejected, has been
completed. The final accepted subpath is extended in both time-directions to
make a full new path. Note that a rejection of a subpath does not imply a
rejection of the MC move itself, but just redirects the process of achieving a
new path from an old path. The time-direction is chosen such that from
$\lambda_{\rm sour}$ the trajectory is propagated backward in time and from
$\lambda_{i}$ forward in time. Due to the placement of $\lambda_{\rm sour}$,
it is nearly guaranteed that the backward extension reaches state $A$. As
$\lambda_{i}$ is also crossed, it is ensured that the path is valid for the
$[i^{+}]$ ensemble, thought it might still be rejected due to a final
acceptance/rejection step, as required by detailed balance Metropolis  _et
al._ (1953).

The WF move, further discussed in Sec. VII, differs with the other moves by
its location of the shooting points. In the WF move, these might be any point
with a corresponding value of the reaction coordinate that is larger than
$\lambda_{i}$ and lower than $\lambda_{B}$ (or $\lambda_{\rm cap}$ if a so-
called _cap interface_ is set, see Sec. VII). From this point, no specific
requirements are needed for the velocities so that they are most conveniently
generated from a Maxwell-Boltzmann distribution for the temperature of
interest. From the new phase point, MD steps are generated forward and
backward in time until $\lambda_{B}$ (or $\lambda_{\rm cap}$) or $\lambda_{i}$
is crossed. The subpath is accepted unless it reaches $\lambda_{B}$ (or
$\lambda_{\rm cap}$) in both time directions. In that case it would be
rejected and the next shot is taken again from latest accepted subpath or the
previous segment of the old path if no accepted subpaths yet exist. After
finishing the number of desired subpaths, the last accepted one is extended to
the stable states, like in SS and WT. While the WF move is slightly more
wasteful with respect to the MD moves compared to SS, the velocity generation
is much simpler which can have both practical and fundamental advantages
compared to SS and WT. These are further discussed in Sec. VI. The name of the
WF move is derived from the visual resemblance between the set of full paths
and subpaths and the top of a wire fence.

The subtrajectory moves go against strategy i) as these MC moves require more
MD steps than just the number of MD steps for generating a new path. These
moves are nevertheless more efficient because they utilize strategy ii): the
statistical inefficiency of the sampling is reduced and, therefore, fewer
trajectories are needed to achieve a desired statistical error. Like with the
standard shooting move, a final acceptance/rejection step should ensure that
the correct statistical distribution of paths is sampled. However, due to the
complexity of the subtrajectory move, the design and mathematical validation
of the acceptance rule is substantially more complex and is derived from the
so-called superdetailed balance Frenkel and Smit (2002) principle.

##  III Superdetailed balance

The term superdetailed balance was first introduced within the context of
configurational bias MC (CBMC), Siepmann and Frenkel (1992); Vlugt  _et al._
(1999); Frenkel and Smit (2002) which is an effective method to study the
adsorption of polymers in nanoporous materials such as zeolites. In this
algorithm, polymers are removed and then regrown atom by atom such that any
overlap between the polymer and the zeolite’s walls and other polymers is
avoided. In this growth process, several attempted branch formations are
tested and potentially rejected. Therefore, a specific final accepted
configuration could, in principle, be obtained from the old configuration via
an infinite number of ways (construction paths). As a result, the Metropolis-
Hastings Hastings (1970) rule for deriving acceptance probabilities becomes
impractical as it requires the knowledge on the generation probabilities of
all these branches, accepted and rejected, that need to be summed up. This
issue is overcome in CBMC using the superdetailed balance principle, which can
be formulated in terms of a construction path $\chi$ and its inverse
$\overline{\chi}$.Riccardi  _et al._ (2017) That is, we not only require
detailed balance between any possible old state and new state, but we require
this for any specific route that connects these two states:

$\displaystyle P_{\rm acc}={\rm min}\left[1,\frac{P({\rm path}^{(n)})P_{\rm
gen}({\rm path}^{(n)}\rightarrow{\rm path}^{(o)}\textrm{ via
}\overline{\chi})}{P({\rm path}^{(o)})P_{\rm gen}({\rm
path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{ via }\chi)}\right]$ (1)

where $P_{\rm gen}({\rm path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{ via
}{\chi})$ is the generation probability to to generate the new state (path in
our case) from the old state via construction path $\chi$ and $P_{\rm
gen}({\rm path}^{(n)}\rightarrow{\rm path}^{(o)}\textrm{ via
}\overline{\chi})$ is the generation probability to to generate the old state
from the new state via the reverse construction path $\overline{\chi}$.

In subtrajectory moves the “construction” path does not only describe the MD
extensions of the final path, but also the sequence of subtrajectories
including the failed ones. For SS and WT, the unsuccessful velocity
generations, that do not obey the one-step crossing condition, should also be
considered as part of the construction path $\chi$. In other words, $\chi$
consists of several steps and the generation probability “via $\chi$” is given
by the product of generation probabilities of each step.

For each construction path $\chi$ there should exist an unique reverse
construction path $\overline{\chi}$. Roughly said, when $\chi$ represents a
sequence of algorithmic steps, $\overline{\chi}$ will typically consist of the
reverse steps in reverse order. However, some groups of consecutive steps
might actually happen in the same order. In fact, there is no unique way to
define “a reverse”, but for a given definition there will be a one-to-one
relation between any possible $\chi$ and its reverse $\overline{\chi}$ and
with that, valid acceptance/rejection rules can be derived based on the
superdetailed balance, Eq. 1.

Yet, the definition of the reverse should be chosen such that the acceptance
probability is computable and not negligibly small in the majority of cases.
Therefore, the mathematical definition for the inverse is taken such that the
probabilities of most of the algorithmic steps in the expressions for $P_{\rm
gen}({\rm path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{ via }{\chi})$ and
$P_{\rm gen}({\rm path}^{(n)}\rightarrow{\rm path}^{(o)}\textrm{ via
}\overline{\chi})$ will cancel.

For instance, if we represent the construction path as vector containing the
different steps in chronological order, $\chi$ could look like

$\displaystyle\chi=[s^{0},t^{1},t^{2},s^{3},s^{4},t^{5},s^{6}]$ (2)

which shows that there were 6 subtrajectories generated of which there were 3
failed trials $t^{1},t^{2}$ and $t^{5}$. The initial step involves cutting out
the very first subtrajectory $s^{0}$ from the old path, while the final step
implies not only the generation of the last subtrajectory $s^{6}$ but also its
extension to a full trajectory. The reverse construction path in this case is
conveniently defined as

$\displaystyle\overline{\chi}=[s^{6},s^{4},t^{5},s^{3},s^{0},t^{1},t^{2}]$ (3)

So the order of the steps is not completely reversed, but the reverse order
takes place on groups of consecutive steps, a group being a successful
subtrajectory with all its failed trials that follow. The reason for this
inverse is that Eq. 2 shows that trial trajectory $t^{5}$ can be generated
starting from $s^{4}$, but this is not necessarily the case from $s^{6}$.
Reversely, as $s^{6}$ was generated from $s^{4}$, they share a common
configuration point, which makes it possible to to generate $s^{4}$ from
$s^{6}$. There is, however, no reason whatsoever that $s^{6}$ and $t^{5}$
share a common configuration point. Hence, if we would consider the reverse to
be $\overline{\chi}=[s^{6},t^{5},s^{4},\ldots]$, $P_{\rm gen}({\rm
path}^{(n)}\rightarrow{\rm path}^{(o)}\textrm{ via }\overline{\chi})$ would
most likely be zero as $\overline{\chi}$ itself cannot be generated. In
contrast, the inverse based on the grouped reordering, Eq. 3, contains
generation probabilities like the probability to generate $t^{5}$ given
$s^{4}$ which appear both in $\chi$ and $\overline{\chi}$. Therefore, all the
generation probabilities of failed trajectories cancel in Eq. 1. Likewise, all
failed velocity generations in SS and WT that do not obey the one-step
crossing condition cancel out for the same reason as is shown in Ref. Riccardi
_et al._ , 2017.

Excluding all the failed steps that will cancel in Eq. 1, we can write for
$P_{\rm gen}({\rm path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{ via
}{\chi})$:

$\displaystyle P_{\rm gen}({\rm path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{
via }{\chi})\propto P_{\rm sel}(s^{0}|{\rm path}^{(o)})\times$ $\displaystyle
P_{\rm sel}(r^{0,3}|s^{0})P_{\rm gen}(v^{0,3})P_{\rm MD}(s^{3}|x^{0,3})\times$
$\displaystyle P_{\rm sel}(r^{3,4}|s^{3})P_{\rm gen}(v^{3,4})P_{\rm
MD}(s^{4}|x^{3,4})\times$ $\displaystyle P_{\rm sel}(r^{4,6}|s^{4})P_{\rm
gen}(v^{4,6})P_{\rm MD}(s^{6}|x^{4,6})\times$ $\displaystyle P_{\rm sel}({\rm
td})P_{\rm MD}({\rm path}^{(n)}|s^{6})$ (4)

Here, $P_{\rm sel}(s^{0}|{\rm path}^{(o)})$ is the probability for selecting
$s^{0}$ from the old path ${\rm path}^{(o)}$ and $P_{\rm sel}(r^{0,3}|s^{0})$
is the selection probability of choosing point $r^{0,3}$ from the subpath
$s^{0}$ as the shooting point. Since $r^{0,3}$ is a shooting point to go from
$s^{0}$ to $s^{3}$, it is a configuration point that $s^{0}$ and $s^{3}$ have
in common. $P_{\rm gen}(v^{0,3})$ is the probability for generating the
velocities $v^{0,3}$ that are the velocities of $s^{3}$ at the corresponding
configuration point $r^{0,3}$. $P_{\rm MD}(s^{3}|x^{0,3})$ is the chance that
starting from phase point $x^{0,3}=(r^{0,3},v^{0,3})$, the MD integrator
produces subpath $s^{3}$ by integrating the equations of motion forward and
backward in time. The MD integrator can be based on actual Newtonian MD,
Langevin, Brownian, etc. Likewise, $P_{\rm MD}({\rm path}^{(n)}|s^{6})$ is the
chance that the new path ${\rm path}^{(n)}$ is produced by extending the final
subpath $s^{6}$. Finally, $P_{\rm sel}({\rm td})$ is the selection probability
for direction of time along the new path. Note that the direction of time
along the subpaths is irrelevant in WT and WF. In SS, subpaths do have a sort
of direction as the next shooting always takes place at the last $\lambda_{i}$
crossing.Riccardi  _et al._ (2017)

For the reverse construction path, Eq. 3, we can write

$\displaystyle P_{\rm gen}({\rm path}^{(n)}\rightarrow{\rm path}^{(o)}\textrm{
via }\overline{\chi})\propto P_{\rm sel}(s^{6}|{\rm path}^{(n)})\times$
$\displaystyle P_{\rm sel}(r^{6,4}|s^{6})P_{\rm gen}(v^{6,4})P_{\rm
MD}(s^{4}|x^{6,4})\times$ $\displaystyle P_{\rm sel}(r^{4,3}|s^{4})P_{\rm
gen}(v^{4,3})P_{\rm MD}(s^{3}|x^{4,3})\times$ $\displaystyle P_{\rm
sel}(r^{3,0}|s^{3})P_{\rm gen}(v^{3,0})P_{\rm MD}(s^{0}|x^{3,0})\times$
$\displaystyle P_{\rm sel}({\rm td})P_{\rm MD}({\rm path}^{(o)}|s^{0})$ (5)

Now, it becomes apparent that most terms will cancel out in Eq. 1 when we take
the ratio between Eq. 5 and Eq. 4. First of all, the time-direction is chosen
with a 50% probability such that $P_{\rm sel}({\rm td})=0.5$. Then, we can use
the fact that a path probability can be written in terms of a phase point
probability times the MD generation probability

$\displaystyle P({\rm path})=\rho(x)P_{\rm MD}({\rm path}|x)$ (6)

where $\rho(x)$ is the phase space equilibrium density for any phase point $x$
that is part of the path. Riccardi  _et al._ (2017) For a phase point
$x=(r,v)$ this can be split into

$\displaystyle\rho(x)=\rho_{r}(r)\rho_{v}(v)$ (7)

where $\rho_{r}$ and $\rho_{v}$ are, respectively, the configuration
(Boltzmann) distribution and the velocity Maxwell-Boltzmann distribution
(possibly subjected to bond- and angle constraints if applicable). Further, as
generating new velocities in Eqs. 4 and 5 is based on the velocity
distribution, $P_{\rm gen}(v)=\rho_{v}(v)$, we can substitute all $P_{\rm MD}$
terms in Eqs. 4 and 5, e. g.:

$\displaystyle P_{\rm gen}(v^{4,6})P_{\rm MD}(s^{6}|x^{4,6})$
$\displaystyle=\frac{P_{\rm gen}(v^{4,6})P(s^{6})}{\rho(x^{4,6})}$
$\displaystyle=\frac{\rho_{v}(v^{4,6})P(s^{6})}{\rho(x^{4,6})}=\frac{P(s^{6})}{\rho_{r}(r^{4,6})}$
$\displaystyle P(s^{6})P_{\rm MD}({\rm path}^{(n)}|s^{6})$
$\displaystyle=P({\rm path}^{(n)})$ (8)

Applying these operations to Eqs. 4 and 5, we get

$\displaystyle P_{\rm gen}({\rm path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{
via }{\chi})\propto P_{\rm sel}(s^{0}|{\rm path}^{(o)})P_{\rm sel}({\rm
td})\times$ $\displaystyle P_{\rm sel}(r^{0,3}|s^{0})P_{\rm
sel}(r^{3,4}|s^{3})P_{\rm sel}(r^{4,6}|s^{4})\times$ $\displaystyle
P(s^{3})P(s^{4})P({\rm
path}^{(n)})/[\rho_{r}(r^{0,3})\rho_{r}(r^{3,4})\rho_{r}(r^{4,6})]$
$\displaystyle P_{\rm gen}({\rm path}^{(n)}\rightarrow{\rm path}^{(o)}\textrm{
via }\overline{\chi})\propto P_{\rm sel}(s^{6}|{\rm path}^{(n)})P_{\rm
sel}({\rm td})\times$ $\displaystyle P_{\rm sel}(r^{6,4}|s^{6})P_{\rm
sel}(r^{4,3}|s^{4})P_{\rm sel}(r^{3,0}|s^{3})\times$ $\displaystyle
P(s^{4})P(s^{3})P({\rm
path}^{(o)})/[\rho_{r}(r^{6,4})\rho_{r}(r^{4,3})\rho_{r}(r^{3,0})]$ (9)

In the ratio of these two equations more terms will cancel out as
$r^{\alpha,\beta}=r^{\beta,\alpha}$. Further, $P_{\rm sel}({\rm td})=0.5$ as
stated before. In all subtrajectory moves, $P_{\rm sel}(r|s^{\alpha})$ is
either a fixed number (SS and WT) or it depends on $s^{\alpha}$, but not on
$r$ (WF). In SS, the shooting point is selected from the last crossing with
$\lambda_{i}$ and therefore $P_{\rm sel}(r|s^{\alpha})=2$ (the phase point
just before or after $\lambda_{i}$). In WT, it is randomly chosen from a
crossing with either $\lambda_{i}$ or $\lambda_{\rm sour}$ and therefore
$P_{\rm sel}(r|s^{\alpha})=4$. With stochastic dynamics one can also opt to
choose only the inner points Riccardi  _et al._ (2017) such that $P_{\rm
sel}(r|s^{\alpha})=2$. In WF any point of the subpath that lies between
$\lambda_{i}$ and $\lambda_{B}$ (or $\lambda_{\rm cap}$) can be chosen. In all
these cases the $P_{\rm sel}(r|s^{\alpha})$ terms with identical $s^{\alpha}$
cancel out in the ratio. That means that the only terms that remain depend on
the first and last subpath ($s^{0}$ and $s^{6}$), or on the full paths (${\rm
path}^{(o)}$ and ${\rm path}^{(n)}$):

$\displaystyle\frac{P_{\rm gen}({\rm path}^{(n)}\rightarrow{\rm
path}^{(o)}\textrm{ via }\overline{\chi})}{P_{\rm gen}({\rm
path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{ via }{\chi})}$ $\displaystyle=$
(10) $\displaystyle\frac{P_{\rm sel}(s^{6}|{\rm path}^{(n)})P_{\rm
sel}(r^{6,4}|s^{6})P({\rm path}^{(o)})}{P_{\rm sel}(s^{0}|{\rm
path}^{(o)})P_{\rm sel}(r^{0,3}|s^{0})P({\rm path}^{(n)})}$ $\displaystyle=$
$\displaystyle\frac{P_{\rm sel}(r^{6,4}|{\rm path}^{(n)})P({\rm
path}^{(o)})}{P_{\rm sel}(r^{0,3}|{\rm path}^{(o)})P({\rm path}^{(n)})}$
$\displaystyle=\frac{P({\rm path}^{(o)})/M^{(n)}}{P({\rm
path}^{(n)})/M^{(o)}}$

where in the third expression we contracted the selection probabilities
involving the two-steps (first selecting $s^{0}$ or $s^{6}$, then selecting
$r^{0,3}$ or $r^{6,4}$) to the chance of selecting the very first successful
crossing point from the existing full path. Finally, the latter was replaced
by $1/M^{(n)}$ and $1/M^{(o)}$ where $M^{(n)}$ and $M^{(o)}$ are the number of
different equal probable possibilities to select a shooting point for
generating a subtrajectory from the new and old full path, respectively.

If we substitute Eq. 10 into Eq. 1, we obtain a rather simple expression for
the acceptance:

$\displaystyle P_{\rm acc}={\rm min}\left[1,\frac{M^{(o)}}{M^{(n)}}\right]$
(11)

In SS, $M^{(o)}$ and $M^{(n)}$ are simply proportional to the number of
crossing points of the old and new path with $\lambda_{i}$, while for WT these
are proportional to the number segments that can be cut out of these
trajectories that connect $\lambda_{\rm sour}$ and $\lambda_{i}$. Riccardi
_et al._ (2017) In WF, these relate to the number of points between
$\lambda_{i}$ and $\lambda_{B}$. If a so-called _cap_ -interface is defined,
$M^{(o)}$ and $M^{(n)}$ relate to the number of points between $\lambda_{i}$
and $\lambda_{\rm cap}$ excluding any points lying on a segment $\lambda_{\rm
cap}\rightarrow\lambda_{\rm cap}$ without crossing $\lambda_{i}$.

Eq. 11 can also be combined with an early rejection scheme as was introduced
in Ref. van Erp  _et al._ , 2003. In the standard approach one would complete
the MC move, compute the acceptance probability, Eq. 11, take a uniform random
number $\alpha$ between 0 and 1, and then accept if $\alpha<P_{\rm acc}$ and
reject otherwise. In the early rejection scheme, the random number $\alpha$ is
taken first and the move is rejected as soon as $M^{(n)}>M^{(o)}/\alpha$. In
normal shooting, this provides a considerable speed up since long paths have a
high chance to get rejected. Using the early rejection scheme a lot of
unnecessary MD steps can be avoided as these path can be stopped whenever they
exceed the predetermined maximum length. Yet, for the subtrajectory moves the
_high-acceptance_ scheme is preferable as we discuss in Sec V. In the next
section we show why the subtrajectory moves allow us to sample fewer
trajectories than with standard shooting via a reduction of the statistical
inefficiency.

##  IV Statistical inefficiency

The principal property that is computed in the $[i^{+}]$ ensemble is the local
crossing probability $P_{A}(\lambda_{i+1}|\lambda_{i})$. This is the history
dependent conditional probability that the system, given it crosses
$\lambda_{A}$ and then crosses $\lambda_{i}$, crosses $\lambda_{i+1}$ before
$\lambda_{A}$. In the _post hoc_ analysis, this local crossing probability is
simply the fraction of sampled path in the $[i^{+}]$ ensemble that happen to
cross $\lambda_{i+1}$ in addition to $\lambda_{i}$. Once these are accurately
enough determined, the global crossing probability
$P_{A}(\lambda_{B}|\lambda_{A})$ is obtained from van Erp  _et al._ (2003);
Cabriolu  _et al._ (2017):

$\displaystyle
P_{A}(\lambda_{B}|\lambda_{A})=\prod_{i=0}^{n-1}P_{A}(\lambda_{i+1}|\lambda_{i})$
(12)

where $\lambda_{0}=\lambda_{A}$ and $\lambda_{n}=\lambda_{B}$. The above
expression is exact since the local crossing probabilities include the full
history dependence ($\lambda_{A}\rightarrow\lambda_{i}$) in their condition.
van Erp (2012) An alternative approximate expression for the global crossing
probability is used in partial path TIS Moroni  _et al._ (2003) in which the
amount of spatial memory is reduced though not set to zero, as in milestoning
Faradjian and Elber (2004). The global crossing probability gives the rate of
the transition when multiplied with $f_{A}$, the conditional flux through
$\lambda_{A}$.

In TIS, the flux is calculated by straightforward MD where the system is
prepared in state $A$ and then the number of crossings with $\lambda_{A}$ per
time unit is computed. If a spontaneous transition to state $B$ takes place,
which is unlikely for a rare event, the simulation is paused, reinitiated in
state $A$ and then continued. RETIS computes the flux term differently as it
does not use a single continuous MD simulation. Instead, the $[0^{-}]$ path
ensemble is introduced to explore the $A$ state, and the flux is derived from
the average path lengths in $[0^{-}]$ and $[0^{+}]$.van Erp (2007) In addition
to rate constants, the overall crossing probability can also be used to
compute permeability coefficients Ghysels  _et al._ (2021) and activation
energies.Dellago and Bolhuis (2004); van Erp and Bolhuis (2005)

Considering the $j$-th path in the simulation for path ensemble $[i^{+}]$, the
main output of sample $j$ (the generated path) that is relevant for the
computation of the crossing probability is simply the observation of whether
it crosses $\lambda_{i+1}$ or not. We can describe this by a characteristic
function $h_{j}$ which equals 1 if $\lambda_{i+1}$ is crossed and 0 otherwise.
The simulation estimate of the local crossing probability, $p(m)$, after $m$
MC moves is then expressed as

$\displaystyle p(m)=\frac{1}{m}\sum_{j=0}^{m-1}h_{j}\approx
P_{A}(\lambda_{i+1}|\lambda_{i})$ (13)

where the index counter starts from zero for mathematical convenience.

For finite $m$, the value of $p(m)$ will not be exact and the absolute error,
$\epsilon_{a}$, is defined as the the standard deviation of the mean
$\sigma_{p(m)}$. This is essentially the standard deviation in possible $p(m)$
results if the simulation experiment would be carrier out multiple times.
Mathematically we can write this as

$\displaystyle\epsilon_{a}=\sigma_{p(m)}=\sqrt{\left\langle(p(m)-p)^{2}\right\rangle}$
(14)

where $p=p(\infty)=P_{A}(\lambda_{i+1}|\lambda_{i})$ and the brackets
$\left\langle\cdot\right\rangle$ refer to the perfect ensemble sampling
average. This can be viewed as the hypothetical average that is obtained after
repeating the simulation an infinite number of times starting with initial
conditions that are randomly drawn form a perfect statistical equilibrium
distribution. In other words, we have $\left\langle
p(1)\right\rangle=\left\langle p(m)\right\rangle=p$. Further, since detailed
balance MC moves conserve the equilibrium distribution, Frenkel and Smit
(2002) the absolute value of the index $j$ is irrelevant and $\left\langle
h_{0}\right\rangle=\left\langle h_{1}\right\rangle=\left\langle
h_{j}\right\rangle=p$ and $\left\langle h_{j}h_{k}\right\rangle=\left\langle
h_{0}h_{k-j}\right\rangle$ for any $j,k$. Using this, one can show that van
Erp (2006):

$\displaystyle\sigma_{p(m)}^{2}$
$\displaystyle=\frac{\sigma_{p(1)}^{2}}{m}{\mathcal{N}},\quad{\mathcal{N}}=\left[1+2n_{c}\right]$
(15)

where ${\mathcal{N}}$ is called the statistical inefficiency and $n_{c}$ is
the correlation number which is the integral of the correlation function
$C(j)$:

$\displaystyle n_{c}=\sum_{j=1}^{\infty}C(j),\quad
C(j)=\frac{\left\langle(h_{0}-p)(h_{j}-p)\right\rangle}{\left\langle(h_{0}-p)^{2}\right\rangle}$
(16)

As the output $h_{j}$ of a single sample is either 1 with a probability $p$ or
0 with a probability $(1-p)$, the sample standard deviation $\sigma_{p(1)}$
can be simplified

$\displaystyle\sigma_{p(1)}^{2}$
$\displaystyle=\left\langle(p(1)-p)^{2}\right\rangle=\left\langle(h_{0}-p)^{2}\right\rangle$
$\displaystyle=p(1-p)^{2}+(1-p)(0-p)^{2}=p(1-p)$ (17)

Via Eqs. 15, 16, and 17, we can write for the relative error:

$\displaystyle\epsilon_{r}=\frac{\epsilon_{a}}{p}=\sqrt{\frac{1-p}{p}\frac{\mathcal{N}}{m}}$
(18)

Eq. 18 shows that for a fixed number of MC moves $m$, the larger the local
crossing probability $p=P_{A}(\lambda_{i+1}|\lambda_{i})$, the lower the
relative error. Hence, the result in simulation $[i^{+}]$ converges faster
when the difference between $\lambda_{i}$ and $\lambda_{i+1}$ is small, but
this will obviously increase the number of path ensembles needed. Analytical
results on model systems suggest that the optimum placement of interfaces in
TIS is achieved when $p\approx 0.2$ for all ensembles.van Erp (2006) In RETIS
the optimum is expected to be slightly higher as this would lead to more
successful swaps. Likewise, the optimum is also slightly higher if the
weighted histogram analysis method (WHAM) Ferrenberg and Swendsen (1989) is
used instead of single-point matching to determine the total crossing
probability. In this approach, the crossing statistics of path ensemble
$[i^{+}]$ is not limited to the fraction of paths crossing $\lambda_{i+1}$,
but also the fractions for crossing $\lambda_{i+2}$, $\lambda_{i+3}$, etc. are
used to get a slightly more accurate estimate of Eq. 12.van Erp  _et al._
(2016); Rogal  _et al._ (2010)

If the sampling between successive MC moves is completely uncorrelated, we
have that
$\left\langle(h_{0}-p)(h_{j}-p)\right\rangle=\left\langle(h_{0}-p)\right\rangle\cdot\left\langle(h_{j}-p)\right\rangle=0\cdot
0=0.$ This would imply that $C(j)=n_{c}=0$ and ${\mathcal{N}}=1$. In this
case, if $p=0.2$, there are about $m=400$ trajectories required to obtain an
$\epsilon_{r}=10\%$ error. For ${\mathcal{N}}>1$, one would need
$m/{\mathcal{N}}=m_{u}=400$ to get the same error. Here, $m_{u}$ is called the
number of effectively uncorrelated samples.

In general, $C(j)\neq 0$ except for the limit $j\rightarrow\infty$ as
correlation decays. If a MC move is rejected at step $j$, then the previous
sample is kept and recounted such that sample $j$ is identical to sample
$j-1$. Hence, if there are $j$ consecutive rejections, sample $j$ is identical
to sample 0 causing correlation over multiple steps. Even if the $j$-th step
is accepted, it tends to have some similarity with the previous sample.
Therefore, there is an high probability for $h_{j}=h_{j-1}$, even if the
samples are not identical. The correlations lead to a sampling output
$(h_{0},h_{1},h_{2},\ldots)$ with long rows of consecutive zeros and
consecutive ones.

To illustrate this effect with a mathematical example: suppose that the MC
move has a probability $\pi_{R}$ to remain unchanged such that $h_{j}=h_{j-1}$
and a probability $\pi_{M}=1-\pi_{R}$ to actually make a move that potentially
(but not necessarily) changes the output: the new sample yields $h_{j}=1$ with
a probability $p$ and $h_{j}=0$ with a probability $(1-p)$. As shown in the
appendix, for this mathematical model the statistical inefficiency equals:

$\displaystyle{\mathcal{N}}=\frac{2-\pi_{M}}{\pi_{M}}$ (19)

This shows that for a typical MC acceptance probability of 50%, the effect of
rejections alone causes the statistical inefficiency to be equal to 3. The
situation is usually worse in complex systems and also more difficult to
identify than merely by the presence of rows of consecutive ones or zeros. For
instance inter- and intramolecular changes of reactants could temporarily
boost or reduce the probability of a transition. The same kind of fluctuations
in the temporary transition probability can be caused by the local solvent
structure and the position and orientation of catalytic molecules. These
describe degrees of freedom that are orthogonal to the reaction coordinate.

We can examine this by a slightly more complex model where we assume that
there are two phases $\alpha$ and $\beta$, described by the orthogonal degrees
of freedom, which occur with probabilities $P_{\alpha}$ and
$P_{\beta}=1-P_{\alpha}$. Let $p_{\alpha}$ and $p_{\beta}$ be the
corresponding local crossing probabilities along the reaction coordinate for
these phases such that: $p=P_{\alpha}p_{\alpha}+P_{\beta}p_{\beta}$. Analogous
to the above, let $\pi_{\rho}$ be the chance to not update the phase, and
$\pi_{\mu}=1-\pi_{\rho}$ be the chance to freshly choose between phase
$\alpha$ or $\beta$ with respective probabilities $P_{\alpha}$ and
$P_{\beta}$. As shown in the appendix, in this case the statistical
inefficiency equals:

$\displaystyle{\mathcal{N}}=\frac{2K_{s}-\pi_{\mu}(2K_{s}-1)}{\pi_{\mu}}$ (20)

where $K_{s}$ is a system parameter that does not depend on the type of MC
move:

$\displaystyle
K_{s}=\frac{P_{\alpha}P_{\beta}(p_{\alpha}-p_{\beta})^{2}}{p(1-p)}=\frac{(p-p_{\alpha})(p_{\beta}-p)}{p(1-p)}$
(21)

Note that $K_{s}=0$ whenever $p_{\alpha}=p_{\beta}$, which gives
${\mathcal{N}}=1$. This would be the case if all TIS interfaces are placed at
isocommittor surfaces, which partly supports the hypothesis of Ref. Vanden-
Eijnden  _et al._ , 2008 that stated that path sampling simulations are most
efficient if the reaction coordinate $\lambda$ equals the committor. However,
although this surely minimizes the statistical inefficiencies, the mean path
lengths in the path ensembles also depend on the choice of the reaction
coordinate $\lambda$. If this is included in the analysis, the hypothesis is
at least not generally true. van Erp (2012)

Now assume that not all generated paths are saved and analyzed, but instead
only every $N_{s}$-th path is kept. While this will cause a reduction in the
number of samples from $m$ to $m/N_{s}$, it does not necessarily reduce the
number of uncorrelated samples $m_{u}$ as the statistical inefficiency between
saved samples is also reduced. In particular, the “remain” probability between
saved samples changes from $\pi_{\rho}$ to $\pi_{\rho}^{N_{s}}$ and,
therefore, the “move” probability changes from $\pi_{\mu}$ to
$1-\pi_{\rho}^{N_{s}}=1-(1-\pi_{\mu})^{N_{s}}$. The statistical inefficiency
between saved samples is henceforth:

$\displaystyle{\mathcal{N}}(N_{s})=\frac{2K_{s}-(1-(1-\pi_{\mu})^{N_{s}})(2K_{s}-1)}{1-(1-\pi_{\mu})^{N_{s}}}$
(22)

Eq. 22 shows that the statistical efficiency indeed goes down with increasing
$N_{s}$ up to an asymptote equal to 1. Taking the power series up to first
order in $\pi_{\mu}$, we see that the initial downfall is inversely linear:

$\displaystyle{\mathcal{N}}(N_{s})\approx\frac{2K_{s}-N_{s}\pi_{\mu}(2K_{s}-1)}{N_{s}\pi_{\mu}}\approx\frac{{\mathcal{N}}(1)}{N_{s}}$
(23)

where we assumed $N_{s}\pi_{\mu}\ll 1$. As a result, saving every $N_{s}$-th
path instead of all paths will not affect much the post-simulation analysis in
terms of accuracy. The reduction in the number of data points from $m$ to
$m/N_{s}$ is compensated by a lower statistical inefficiency such that the
number of uncorrelated samples $m_{u}$ remains nearly unchanged. While this
allows for obvious data storage savings, reducing both the memory and time for
writing to disk, it also paves the way to reduce MD steps as is shown in Fig.
2.

Figure 2: Illustration of wasted MD steps in shooting and WF. Panel a) shows
six consecutive paths being generated by the shooting move where only the
solid golden paths, with index 0 and 5, are being saved. Panel b) gives an
equivalent scenario in the WF algorithm showing that considerable fewer MD
steps are needed to obtain the same paths 0 and 5 via $N_{s}=5$
subtrajectories. Still, WF is not as thrifty as SS and WT since only parts of
the subtrajectories, shown in panel c), actually contribute to the sampling
progress to get from path 0 to 5. The additional steps in panel b) are
seemingly “wasted” but still needed for the superdetailed balance relation. SS
and WT do not generated wasted MD steps, but rely on a one-step crossing
condition as discussed in the main text.

The figure illustrates a hypothetical MC sequence in path sampling of six
consecutive paths, labeled 0 to 5, where the shooting point has an order
parameter larger than $\lambda_{i}$. If only every fifth path is saved, only
path 0 and 5 are considered as in Fig. 2-a). Although the intermediate paths
contribute for their decorrelation, it is clear that many MD steps can be
omitted, as exploited by the subtrajectory moves. Fig. 2-b) shows a scenario
where the same final path is being generated with a set of hypothetical WF
subtrajectories resembling the top scenario. Instead of five full
trajectories, only four short subtrajectories and one full trajectory are
needed to establish a new full path (path 5) from the old one (path 0). Based
on this principle alone, the relative efficiency gain $\eta$ of subtrajectory
moves compared to standard shooting is expected to be

$\displaystyle\eta(N_{s})=\frac{N_{s}L_{p}}{L_{p}+(N_{s}-1)L_{s}}$ (24)

where $L_{p}$ and $L_{s}$ are, respectively, the average length of a full path
and a subpath. Still, if we purely focus on the MD steps that are required to
allow for the progression from path 0 to path 5, even fewer MD steps are
needed as shown in Fig. 2-c). Yet, the “extra” (wasted) MD steps in panel b)
are required for the superdetailed balance as discussed in Sec. III. Wasted MD
steps are avoided in SS and WT where the shooting always happens at an
interface (see Figs. 1-a) and b)). The price to be paid for this is the
additional complication with regard to the one-step crossing condition (see
Sec. VI). But even with a slightly higher MD waste, the WF move requires
considerably less MD steps than standard shooting.

Eq. 24 levels off to a constant $L_{p}/L_{s}$ for increasing $N_{s}$.
Likewise, Eqs. 23 and 22 show that the trend
${\mathcal{N}}(N_{s})={\mathcal{N}}(1)/N_{s}$ is not sustained for increasing
$N_{s}$ as ${\mathcal{N}}$ ultimately levels off to 1. It is henceforth
assumed that while efficiency initially increases quite rapidly as function of
$N_{s}$, it can not surpass $L_{p}/L_{s}$ and ultimately even decreases when
${\mathcal{N}}(N_{s})$ levels off. Clearly, for the $[0^{-}]$ and $[0^{+}]$
ensemble where $L_{p}=L_{s}$ no gain is expected and one could set $N_{s}=1$
if data storage latency would not be an issue. Therefore, as a rule off thumb,
$N_{s}$ can be set approximately equal to $L_{p}/L_{s}$ such that for
$L_{p}>L_{s}$ the cost of the MC move is less than doubled, while Eq. 24
reaches more than 50% of its anyways unattainable maximum of $L_{p}/L_{s}$.

Although the essence of the above analysis is correct, there is however a
caveat: rejections leave a much heavier mark on the subtrajectory move than on
standard shooting. If, for instance, the extension of the fifth and last
subpath in Fig. 2-b) is rejected, it would imply a complete reset to the
latest accepted full path (path 0) since subpath 4 is not a valid trajectory
and extending subpath 4 after the rejection would violate detailed balance. As
a result, all MD steps of subpaths 1 to 5 are trashed as the next move starts
from path 0 again. Instead, the MC chain will only fall back to path 4
(assuming path 4 was accepted) in standard shooting. It is therefore clear
that rejections in the subtrajectory move approach should be avoided even more
than in the shooting method. This can be achieved with the high-acceptance
procedure that is discussed in the next section.

##  V High-acceptance procedure

As discussed in the previous section, a rejection in the subtrajectory moves
implies a large amount of wasted MD steps. An early rejection scheme, as the
one used in TIS and RETIS with standard shooting (see Sec. III), is also not
so helpful as a rejection cannot be made until the generation of the last
subtrajectory has been initiated. It is, therefore, preferable to combine the
subtrajectory moves with the _high-acceptance_ scheme. Riccardi  _et al._
(2017) The approach uses the following two tricks. First, if the final
subtrajectory has a backward extension ending in state $B$, the MC move is not
directly rejected. Instead, the extension forward in time is completed and, if
it ends in state $A$, the path is time-reversed providing an $A\rightarrow B$
path. The consequence is that the time-direction selection probability $P_{\rm
sel}({\rm td})$ in Eq. 4 is no longer 0.5 for all paths as an $A\rightarrow B$
path can be generated in two ways: either by choosing the correct time
direction immediately, or in reverse. This implies an extra factor two in the
generation probabilities $P_{\rm gen}$, in Eqs. 1 and 10, of the $A\rightarrow
B$ paths compared to $A\rightarrow A$ paths. We henceforth write:

$\displaystyle\frac{P_{\rm gen}({\rm path}^{(n)}\rightarrow{\rm
path}^{(o)}\textrm{ via }\overline{\chi})}{P_{\rm gen}({\rm
path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{ via }{\chi})}$
$\displaystyle=\frac{P({\rm path}^{(o)})q({\rm path}^{(o)})/M^{(n)}}{P({\rm
path}^{(n)})q({\rm path}^{(n)})/M^{(o)}}$ (25)

where

$\displaystyle q({\rm path})=\begin{cases}1&\text{ if }{\rm
path}\in\\{A\rightarrow A\\}\\\ 2&\text{ if }{\rm path}\in\\{A\rightarrow
B\\}\end{cases}$ (26)

The second trick is to slightly change the sampling distribution. Instead of
sampling the correct physical path distribution, $P(\rm path)$, restrained to
the path ensemble’s requirements, an alternative path distribution
$\tilde{P}(\rm path)$ is sampled. From Eqs. 1 and 25, the acceptance
probability thus becomes

$\displaystyle P_{\rm acc}={\rm min}\left[1,\frac{\tilde{P}({\rm
path}^{(n)})P({\rm path}^{(o)})q({\rm path}^{(o)})M^{(o)}}{\tilde{P}({\rm
path}^{(o)})P({\rm path}^{(n)})q({\rm path}^{(n)})M^{(n)}}\right]$ (27)

and to maximize the acceptance, we choose the sampling distribution in
ensemble $[i^{+}]$ as

$\displaystyle\tilde{P}({\rm path})$ $\displaystyle=w_{i}({\rm path}){\bf
1}_{[i^{+}]}({\rm path})\textrm{ with }$ $\displaystyle w_{i}({\rm path})$
$\displaystyle=q({\rm path})M_{\lambda_{i}}({\rm path})$ (28)

where ${\bf 1}_{C}(x)$ is the indicator function that equals 1 if $x$ is part
of set $C$ and 0 otherwise. A subscript $\lambda_{i}$ was added to the last
term $M$, as the number of equal probable possibilities for a first shooting,
generally depends on the interface $\lambda_{i}$. Substituting Eq. 28 in Eq.
27 implies that with high-acceptance

$\displaystyle P_{\rm acc}={\bf 1}_{[i^{+}]}({\rm path}^{(n)})$ (29)

In other words, the new path will always be accepted unless the MC move led to
a path not obeying the ensemble’s definition: starting at $\lambda_{A}$,
ending at $\lambda_{A}$ or $\lambda_{B}$, and having at least one crossing
with $\lambda_{i}$. By construction, the crossing of $\lambda_{i}$ is always
achieved in the subtrajectory moves if the starting condition at $\lambda_{A}$
is met. Hence, the only necessary rejection is when the extension of the final
successful subtrajectory ends at $\lambda_{B}$ in both time-directions.

If no successful subtrajectories were generated after $N_{s}$ attempts,
$s^{0}$ could be extended. However, since this would regenerate the old
trajectory in deterministic dynamics and otherwise a trajectory that is highly
correlated with the old one, it is preferable to reject the move. Other
potential reasons for a rejections could be due to non-convergence of the
atomistic forces in AIMD level calculations. Another potential issue is jumpy
order parameters, Haji-Akbari (2018) such that $M_{\lambda_{i}}$ can be zero
even if the path is actually valid. This issue is further discussed in Sec.
VII.

Exact natural averages can still be obtained by weighting each sample $j$ with
the inverse of of $w_{i}(j)$. For instance, the estimated local crossing
probability, previously defined by Eq. 13, can now be expressed as

$\displaystyle
p(m)=\frac{\sum_{j=0}^{m-1}w_{i}(j)^{-1}h_{j}}{\sum_{j=0}^{m-1}w_{i}(j)^{-1}}\approx
P_{A}(\lambda_{i+1}|\lambda_{i})$ (30)

The effect of the weighting implies that different samples have different
contribution. If a sample $j^{\prime}$ has a much lower than average
$w_{i}^{-1}$ factor, the sample could essentially be removed from Eq. 30
without significantly affecting the estimate $p(m)$. Yet, thanks to this
sample not being rejected, sample $j^{\prime}-1$ is more different than
$j^{\prime}+1$ than it would be in case that $j^{\prime}$ was rejected. This
shows the power of the high-acceptance approach.

The improved acceptance in the subtrajectory move will slightly reduce the
acceptance in the replica exchange move. For instance, if a path $j$ from
ensemble $[i^{+}]$ will be exchanged with a path $k$ from ensemble
$[(i+1)^{+}]$, the acceptance becomes Riccardi  _et al._ (2017):

$\displaystyle P_{\rm acc}={\bf 1}_{[(i+1)^{+}]}(j)\times{\rm
min}\left[1,\frac{w_{i}(k)w_{(i+1)}(j)}{w_{i}(j)w_{(i+1)}(k)}\right]$ (31)

Without high-acceptance, the factor in Eq. 31 after the multiplication sign
equals 1. This means that whenever $j$, the path originating from $[i^{+}]$ is
valid for $[(i+1)^{+}]$, the swap will be accepted. Note that any path in
$[(i+1)^{+}]$ is also valid in $[i^{+}]$. This lower acceptance is not
dramatic since replica exchange moves do not require any MD steps. Therefore,
replica exchange moves have negligible CPU cost. The only exception is the
$[0^{-}]\leftrightarrow[0^{+}]$ swap in which two new paths are generated.
Without high-acceptance, this move is always accepted. For SS and WT the
acceptance remains 100%, but this is not the case for WF. We can solve this
problem for WF in RETIS by sampling the $[0^{-}]$ and $[0^{+}]$ ensembles with
the standard shooting method without high acceptance. Due to this $w_{0^{+}}$
and $w_{0^{-}}$ equal 1 irrespective to the the paths and swapping between
these two ensembles will always be accepted. The absence of high-acceptance is
partly compensated by early rejection (see Sec. III). Moreover, in these
ensembles there is no difference between the average path length of a subpath
and a full path, making the subtrajectory moves anyways not so effective for
these ensembles.

The high-acceptance protocol eliminates the more serious drawbacks of
rejections in the subtrajectory moves compared to shooting. In the next
section we discuss how the one-step crossing condition can be met.

##  VI One-step crossing condition

As discussed above, SS and WT are very thrifty algorithms with respect to the
number of generated MD steps. Yet, the one-step crossing condition puts a
challenge to the implementation. One can eliminate the one-step crossing
condition via the new but less thrifty WF algorithm that is further discussed
in Sec. VII. In this section we discuss a few algorithmic solutions to
overcome the one-step crossing condition in SS and WT. These two approaches
assume that one time step in (RE)TIS is effectively also one MD step.

The one-step crossing can be achieved in different ways. The most
straightforward way is to generate velocities from a Maxwell-Boltzmann
distribution, execute an MD step, calculate the new order parameter, and if
the crossing is established, then the two frames comprising the crossing are
extended at the side above $\lambda_{i}$ to create a new subpath. The problem
with this approach is that, after each velocity generation, an MD step, and
therefore a force calculation is required. Especially if $\lambda_{i}$ is at a
steep slope of the potential energy surface, the two trajectory frames forming
the crossing of a given interface might be rather far apart in
$\lambda$-space. In such cases, if one of the two frames is located in the
very proximity of the interface, it might be extremely unlikely to re-generate
a new one step crossing from the configuration furthest to the $\lambda_{i}$
interface given a random approach to generate velocities.

There are essentially two strategies to reduce the cost for fulfilling the
one-step criterion: i) generate atom velocities from a Maxwell-Boltzmann
distribution and predict the next step’s order parameter without performing an
actual MD step, and ii) generate velocities in a way such that the crossing is
likely achieved after very few attempts. Strategy i) assumes that generating
new velocities is rather computationally inexpensive and the expense of the
one-step crossing condition is mostly provided by the force calculation. This
is the case for AIMD level simulations as these typically consist of just a
few (hundreds of) atoms, while requiring a high CPU demand for the force
calculation. In large classical MD systems with a significant number of atoms,
the velocity generation might actually be equally expensive as a force
calculation. In that case, strategy ii) might be preferable.

###  VI.1 Prediction strategy

The velocity-Verlet Frenkel and Smit (2002) MD integrator propagates a phase
point $x(t)=(r(t),v(t))$ deterministically to a next phase point $x(t+\Delta
t)=(r(t+\Delta t),v(t+\Delta t))$. The integrator is most conveniently
expressed via “intermediate velocities” at $t+\Delta t/2$:

$\displaystyle v(t+\Delta t/2)$ $\displaystyle=v(t)+f(t)\Delta t/(2m)$
$\displaystyle r(t+\Delta t)$ $\displaystyle=r(t)+v(t+\Delta t/2)\Delta t$
$\displaystyle v(t+\Delta t)$ $\displaystyle=v(t+\Delta t/2)+f(t+\Delta
t)\Delta t/(2m)$ (32)

where $m$ is the mass and $f$ are the forces. We used a simplified notation
here, but one should realize that for an $N$ particle system both $r,v$ and
$f$ are 3N-dimensional vectors and $m$ is actual a $3N\times 3N$ diagonal mass
matrix. Further, the forces are determined from the positions: $f(t)=f(r(t))$.

Eq. 32 suggests that one MD step requires two force evaluations, but this is
not the case when the steps of Eq. 32 are called repeatedly in a loop. After
the force calculation at the third step, required to determine $v(t+\Delta
t)$, the forces are stored such that these can be used at the first step of
the next cycle. With the same reasoning, if the forces are known already at
time $t$ from its previous step, a new force evaluation is only needed to
determine $v(t+\Delta t)$, but not $r(t+\Delta t)$. This means that if the
order parameter only depends on geometry, $\lambda=\lambda(r)$, its value at
$t+\Delta t$ can also be determined without the need of doing an actual force
calculation.

When testing the one-step crossing for the selected configuration with
randomized velocities, a new (single step) MD trajectory is started with no
information available from previous MD step. However, the selected
configuration is also part of the previous subpath, so the corresponding
forces could have been known, in principle. When not available, the forces can
be reobtained from the trajectory data without further electronic structure
calculation in AIMD or from the gradient of force field potential in classical
MD. In particular, let $x_{1}=(r_{1},v_{1})$ and $x_{2}=(r_{2},v_{2})$ be two
consecutive phase points of the latest subpath that define a crossing. This
means that $x_{2}$ follows from $x_{1}$ through a single MD step and both
points are at opposite sides of the interface. Therefore, both points are
viable points for shooting off the next subpath. By inverting Eq. 32 we can
derive

$\displaystyle f_{1}=\frac{2\,m\,\left(r_{2}-r_{1}-\Delta
t\,{v}_{1}\right)}{\Delta t^{2}},\quad
f_{2}=\frac{2\,m\,\left(r_{1}-r_{2}+\Delta t\,{v}_{2}\right)}{\Delta t^{2}}$
(33)

So Eq. 33 directly provides the forces on the two potential shooting points by
reading the trajectory data from the subpath. Given that one of these two
points is selected as a shooting point and new randomized velocities are
generated, the coordinates after one MD step can be determined without any
additional force calculation but using just the first two steps of Eq. 32.
Hence, the value of the order parameter after one step can be asserted.

If the prediction suggests that a crossing might be achieved, the MD step is
completed and then the next subtrajectory is generated. If the velocities do
not lead to a crossing, a new velocity randomization is attempted until the
crossing condition is met. As in SS the shooting point selection has to be
maintained, the computation of Eq 33 only needs to be done once for the
generation of each subpath. Naturally, if the MD step integrator is more
complex than velocity-Verlet (due to thermostats, barostats, constraints,
stochasticity), then the prediction becomes more difficult. The method also
works best if a MD step is computationally expensive while regenerating
velocities is relatively cheap. This method is therefore more suitable for
simulations with AIMD level. The approach has been implemented in the PyRETIS
software, and it can be directly used with the CP2K Hutter  _et al._ (2014)
external MD engine. Note that the use of the plain velocity-Verlet MD
integrator is rather common in path sampling since the generation of paths is
already thermostated via the shooting move that allows a change of energy,
while the individual paths have NVE dynamics.

###  VI.2 Alternative velocity generation

The mathematically simple form of Eq. 10 is due to the many terms conveniently
canceling out. For instance, the terms in Eqs. 9, $\rho_{r}(r^{0,3})$,
$\rho_{r}(r^{4,3})$, and $\rho_{r}(r^{4,6})$ in $P_{\rm gen}({\rm
path}^{(o)}\rightarrow{\rm path}^{(n)}\textrm{ via }{\chi})$ cancel out with,
respectively, $\rho_{r}(r^{3,0})$, $\rho_{r}(r^{4,3})$, and
$\rho_{r}(r^{6,4})$ in $P_{\rm gen}({\rm path}^{(n)}\rightarrow{\rm
path}^{(o)}\textrm{ via }\overline{\chi})$ due to the fact that
$r^{\alpha,\beta}=r^{\beta,\alpha}$. However, whereas consecutive (accepted)
subtrajectories share a common configuration point, they do not necessarily
share of common phase point as $v^{\alpha,\beta}\neq v^{\beta,\alpha}$. Here,
$v^{\alpha,\beta}$ refers to the velocities of $s^{\beta}$ at the
configuration point $r^{\alpha,\beta}$, and $v^{\beta,\alpha}$ refers to the
velocities of $s^{\alpha}$ at an identical configuration point. These
velocities have typically not the same orientation nor amplitude. Luckily, the
$\rho_{v}(v^{\alpha,\beta})$ terms still cancel out via $P_{\rm
gen}(v^{\alpha,\beta})=\rho_{v}(v^{\alpha,\beta})$ and
$\rho(x^{\alpha,\beta})=\rho_{r}(r^{\alpha,\beta})\rho_{v}(v^{\alpha,\beta})$
in Eq. 8.

Now, suppose that in a $N$ particle system not all $3N$ velocity components
are regenerated from a Maxwell-Boltzmann distribution, but some velocities
components are kept and some others are inverted (multiplied with -1). These
two velocity groups do not cancel out in Eq. 8 as they are not part of $P_{\rm
gen}$ which implies that the final results changes from
$P(s^{6})/\rho_{r}(r^{4,6})$ to
$P(s^{6})/[\rho_{r}(r^{4,6})\rho_{v}(u^{4,6})]$ where $u^{4,6}$ are the
velocity components that are either unchanged or inverted. Since the
equilibrium velocity distribution is symmetric $\rho_{v}(v)=\rho_{v}(-v)$, and
$u^{4,6}$ is identical to $u^{6,4}$ except for some components having
different sign, all the $\rho_{v}(u^{\alpha,\beta})$ terms cancel in the
ratio, Eq. 10, just like the $\rho_{r}(r^{\alpha,\beta})$ terms.

This allows for different strategies. For instance, if the dynamics is
stochastic, all velocities can simply be inverted. This option was used for WT
in Ref. Riccardi  _et al._ , 2017. Inverting the velocities of specific atoms
or molecules whose coordinates determine the order parameter could also be
effective. The other velocities could be either kept unchanged, randomnized or
a combination. For instance, in protein folding simulations inverting the
velocities of all protein atoms while leaving the velocities of the solvent
molecules (partly) unchanged would make the sampling less diffusive.
Reinspection of Eq. 32 shows that the coordinates of the atoms with the
inverted velocities are mapped exactly back after 1 MD step to the previous
coordinates regardless of the velocities of the other atoms. As a result, the
one-step crossing condition is automatically fulfilled.

This approach requires, however, a single MD step resolution at the interface
crossing. In large molecular systems, it is not desirable to save trajectory
coordinates every MD step as it could overwhelm hard disk capacity and will
result in a loss of effective CPU efficiency due to time that is spend for
writing to disk. An adaptive scheme could be adopted when the frequency of
order parameter determination and the data retention is intensified whenever
the system approaches an interface. Since trajectories can later be swapped in
a replica exchange move, this adaptive approach would have to be carried out
for all interfaces or, at least, in the proximity of neighboring interfaces.
The latter choice might still lead to path ensembles receiving a trajectory
missing the right resolution at the relevant interface. That part of the
trajectory would then have to be reintegrated by MD. While all these issues
can be solved in theory, it puts quite some challenges to the implementation.
Moreover, if the integrator is not deterministic, but involves a thermostat or
barostat, the one-step crossing might still not be guaranteed. Several
velocity generation steps might still needed. These challenges lead us to
derive the WF move that straightforwardly can be implemented in present path
sampling codes like OpenPathSampling Swenson  _et al._ (2019a, b) and PyRETIS
Lervik  _et al._ (2017); Riccardi  _et al._ (2020) with, potentially, any MD
engine.

##  VII Wire Fencing

Compared to the SS and WT moves, the shooting point selection of the WF move
is constructed to avoid the one-step crossing issue altogether. Instead of
restricting the shooting point to sets of crossing points at an interface, WF
allows any phase point between the path ensemble’s specific ensemble
interface, $\lambda_{i}$, and interface $\lambda_{\mathrm{B}}$ to be picked.
To increase the efficiency of the WF move in systems with asymmetric free
energy barriers (See Fig. 3), the selection range and the boundaries of the
subtrajectories can be changed by replacing $\lambda_{\mathrm{B}}$ with an
user-defined _cap_ -interface, $\lambda_{\mathrm{cap}}$ with
$\lambda_{i}<\lambda_{\mathrm{cap}}\leq\lambda_{\mathrm{B}}$ value.

Figure 3: Illustration of an asymmetric barrier where the placement of a cap-
interface, $\lambda_{\rm cap}$, in WF can avoid the generation of long
subtrajectories and too many shooting points being in the basin of attraction
of state $B$.

The presence of a relative flat downhill region after the barrier’s maximum
and before a stable product state implies that transition paths can become
very long. If accepted, the paths will have a large fraction of points at the
right side of the free energy barrier from which shooting has a very high
chance to generate a failed $\lambda_{B}\rightarrow\lambda_{B}$ trajectory.
This problem was also addressed by the spring-shooting method.Brotzakis and
Bolhuis (2016)

In an AIMD level simulation of aqueous silicate condensation, Moqadam  _et
al._ (2017) this issue was solved by defining $\lambda_{B}$ in the RETIS
algorithm at the position of $\lambda_{\rm cap}$ in the figure. After the
simulation was completed, all paths reaching $\lambda_{B}$ were extended in a
straight forward MD simulation. The introduction of the $\lambda_{\rm cap}$
interface make these post-simulation MD extensions redundant.

We will first outline the WF algorithm without a cap-interface (or
$\lambda_{\mathrm{cap}}=\lambda_{B}$) using the high-acceptance protocol. The
introduction of the $\lambda_{\mathrm{cap}}$ only requires a few modifications
that we discuss afterwards.

  1. 1.

From the old path, count the number of frames $M_{\lambda_{i}}^{(o)}$ between
$\lambda_{i}$ and $\lambda_{B}$. If $M_{\lambda_{i}}^{(o)}=0$ we immediately
reject the full MC move. Otherwise continue with the next step.

  2. 2.

Subdivide the $M_{\lambda_{i}}^{(o)}$ points into groups where each group are
the points lying on a segment connecting $\lambda_{i}$ with $\lambda_{B}$ or a
segment connecting $\lambda_{i}$ with itself.

  3. 3.

Select one segment as $s^{0}$ based on a weighted random selection such that
each segment has a chance to be selected proportional to the number of points
it has.

  4. 4.

Set two counter $n_{s}$ and $n_{a}$ equal to zero: $n_{s}=n_{a}=0$. Then start
the following loop: step 5 til 12.

  5. 5.

Select at random one of the configuration points of the last subpath,
$s^{n_{s}}$, as the new shooting point.

  6. 6.

Generate random velocities from a Maxwell-Boltzmann distribution.

  7. 7.

Starting from the configuration point with the new velocities, apply the MD
integrator to go backward and forward in time until $\lambda_{i}$ or
$\lambda_{\mathrm{B}}$ is crossed.

  8. 8.

Increase the $n_{s}$ counter by one: $n_{s}=n_{s}+1$.

  9. 9.

If both time-directions crosses $\lambda_{\mathrm{B}}$, the trial subpath is
rejected. In that case, the previous successful subpath is kept,
$s^{n_{s}}=s^{n_{s}-1}$. Go to step 12. Otherwise, continue with next step.

  10. 10.

Increase the $n_{a}$ counter by one: $n_{a}=n_{a}+1$.

  11. 11.

Accept the trial subpath such that it becomes $s^{n_{s}}$.

  12. 12.

If $n_{s}<N_{s}$, go to step 5. Otherwise, continue with next step.

  13. 13.

If no accepted subpaths have been generated, $n_{a}=0$, stop and reject the
move. Otherwise, continue with the next step.

  14. 14.

Extend the last subpath $s^{N_{s}}$ in both time-directions with MD until
$\lambda_{\mathrm{A}}$ or $\lambda_{\mathrm{B}}$ is hit. If the path ends at
$\lambda_{\mathrm{B}}$ at both time-directions, the whole MC move is rejected.
Otherwise, continue to the next step

  15. 15.

If the path is $\lambda_{B}\rightarrow\lambda_{A}$, reverse the time-direction
of the path.

  16. 16.

Now a new full path has successfully been established. Let $q^{(n)}$ be 2 if
it is a $\lambda_{A}\rightarrow\lambda_{B}$ path. Otherwise, it is 1. Let
$M_{\lambda_{i}}^{(n)}$ be the number of frames between $\lambda_{i}$ and
$\lambda_{B}$. The weight-factor of the path is
$w^{(n)}=q^{(n)}M_{\lambda_{i}}^{(n)}$ that is needed for computing proper
path ensemble averages, Eq. 30, and for a possible swap move via Eq. 31.

Figure 4: Illustration of the $s^{0}$ selection from the old path. Selectable
shooting points are shown in green, end-points by open black circles, and all
other points in red. a) shows the “jumpy orderparameter” case that leads to an
immediate rejection as no selectable points are present. b) and c) show the
selectable points without and with cap-interface, respectively.

The scenario of the potential rejection at step 1, is shown if Fig. 4-a) which
can occur due to a jumpy character of the order parameter. Haji-Akbari (2018)
A typical example is nucleation where the time steps in path sampling is
usually chosen to consist of many MD steps Moroni  _et al._ (2005) for the
reason that computing order parameters for nucleation is rather costly. As a
result, occasionally the order parameter, defined by the size of the largest
cluster, can make sudden jumps such that more than one interface is crossed in
a single RETIS time step.

The path shown in Fig. 4-a) is a valid path in $[i^{+}]$ such that ${\bf
1}_{[i^{+}]}=1$, but $w_{i}=0$ since $M_{\lambda_{i}}=0$. If only WF moves are
applied, such a path has zero probability to be generated. Yet, its
contribution in Eq. 30 to the average, if hypothetically sampled, would be
$w_{i}^{-1}=\infty$ and, therefore, the sampling average becomes ill-defined.
This can be solved by not allowing $w=0$ weights:

$\displaystyle w_{i}({\rm path})$ $\displaystyle={\rm min}\left[1,q({\rm
path})M_{\lambda_{i}}({\rm path})\right]$ (34)

Introducing this small modification of Eq. 28 solves the “division by zero”
problem and has further no impact of the implementation nor on the robustness
of the algorithm. The existence of jumpy trajectories implies that the WF move
is no longer ergodic. A path like the one in Fig. 4-a) can never be made from
a WF move and, vice versa, it can not be destroyed by the WF move if it is fed
as initial path to the algorithm. However, the full sampling remains ergodic
due to the replica exchange moves.

Step 2 is further illustrated in Fig. 4-b). We can identify two groups of
selectable shooting points (in green), one group of seven point lying on a
$\lambda_{i}\rightarrow\lambda_{i}$ segment and one group of nine points on a
$\lambda_{i}\rightarrow\lambda_{B}$ segment. So these segments are selected as
$s^{0}$ with a 7/16 and 9/16 probability, respectively. In the next step, the
points of the selected segment have an equal probability to be selected for
the first shooting.

Despite that all the green points have the same 1/16 probability to be
selected for shooting off the first subpath, the two-step selection process is
needed to fix $s^{0}$. With a single step selection, it could be possible to
first obtain a failed trial path $t^{1}$ that starts from a point at the first
group, followed by a successful subtrajectory that is launched from a point of
the second group. This will break the superdetailed balance as it would not be
possible to generate $t^{1}$ from $s^{0}$ in the reverse path (see the example
construction paths in Eqs. 2 and 3).

The introduction of the cap-interface changes the initial $s^{0}$ selection as
is shown in Fig. 4-c) where, for the same path as panel b), there are now
three groups of two points that can be chosen. Note that not all the points
between $\lambda_{i}$ and $\lambda_{\rm cap}$ are selectable as the points on
a $\lambda_{\rm cap}\rightarrow\lambda_{\rm cap}$ segment should be excluded.
The algorithm is further identical as described above with $\lambda_{\rm cap}$
instead of $\lambda_{B}$ in the main loop (step 5 til 12). Outside the main
loop (step 13-15), $\lambda_{B}$ is not replaced by $\lambda_{\rm cap}$ since
the final extension always shall reach the $A$ or $B$ states. In the final
step (16), $M_{\lambda_{i}}^{(n)}$ is replaced the number of frames between
$\lambda_{i}$ and $\lambda_{\rm cap}$ excluding those on $\lambda_{\rm
cap}\rightarrow\lambda_{\rm cap}$ segments.

##  VIII Numerical results

We tested the WF algorithm on three model systems: a simple one-dimensional
system for which we can perform full RETIS simulations with high convergence,
and two challenging complex systems based on classical MD and AIMD, where our
analysis is more qualitative based on a single path ensemble simulation. The
one-dimensional system describes a single particle in a double-well potential
that is moving following the underdamped Langevin equation as previously
described in Ref. van Erp, 2012. The purpose of these simulations is to show
numerically that the WF method leads indeed to exact results. In addition, due
to the high degree of convergence that can be reached, we also draw some
conclusions on the efficiency compared to standard shooting. However, it
should be taken into account that a larger boost factor is expected for more
complex high-dimensional systems.

The other two systems are part of ongoing projects on which we plan to report
extensively at later publications. The classical MD system describes the thin
film breakage in oil-water mixtures based on the studies Ref. Aarøen  _et al._
, 2022, 2021; Riccardi and Tichelkamp, 2019; Riccardi  _et al._ , 2014. The
system size of this simulation is over 100,000 atoms making the one-step
crossing impracticable as it requires a stop/restart at every MD step.
Instead, in our single path ensemble simulation the coordinates were recorded
every 50 MD steps. The AIMD system describes the electron transfer between
ruthenium ions in a redox reaction taking place in liquid water. To determine
the relative position of the moving electron, the Kohn-Sham orbitals are
projected on maximally localized Wannier Functions Marzari  _et al._ (2012)
whose centers can viewed as “electron positions”. This implies that in order
to compute the order parameter from a configuration point, a full electronic
structure calculation is required. A cheap prediction scheme as described in
Sec. VI is therefore not suitable. For both systems, we show the usefulness of
the cap-interface in practical simulations.

###  VIII.1 Double well 1D barrier

Despite the model’s simplicity, several popular rare event simulation methods,
like forward flux sampling (FFS) Allen  _et al._ (2009); Escobedo  _et al._
(2009) and other splitting based methods,Booth and Hendricks (1984); Melnik-
Melnikov and Dekhtyaruk (2000); Villenaltamirano and Villenaltamirano (1991)
have shown that they can easily fall into a kind of sampling trap when applied
to this system yielding a too low rate and non time-symmetric transition
paths.van Erp (2012)

The double well barrier system consists of a one-dimensional particle moving
in the following potential van Erp (2012)

$V(z)=z^{4}-2z^{2}$ (35)

with underdamped Langevin dynamics. In reduced units, the Boltzmann constant
and mass are set to unity, $k_{B}=m=1$, while the temperature and friction
coefficient are set equal to $T=0.07$ and $\gamma=0.3$. The equations of
motion are propagated using an MD time step equal to $\mathrm{d}t=0.025$. In a
straightforward MD run, the particle will mostly oscillate within one of the
potential minima at $z=-1$ and $z=1$, but also (very) infrequently cross the
transition state at z = 0. During the oscillatory movement, the total energy
of the particle will fluctuate by the random force of the Langevin dynamics.
As a result, the system is effectively two-dimensional in phase space where
the velocity can be considered as an orthogonal degree of freedom. The reason
that FFS and other splitting type methods underestimate the crossing rate is
due to an insufficient sampling of the tail in the velocity distribution.van
Erp (2012) Path sampling methods like RETIS that are based on both forward and
backward in time propagation do not have this issue.

We defined eight RETIS interfaces: $\lambda_{A}=\lambda_{0}=-0.99$,
$\lambda_{1}=-0.8$, $\lambda_{2}=-0.7$, $\lambda_{3}=-0.6$,
$\lambda_{4}=-0.5$, $\lambda_{5}=-0.4$, $\lambda_{6}=-0.3$, and
$\lambda_{B}=\lambda_{7}=1.0$, and ran four RETIS simulations using the
PyRETIS code Lervik  _et al._ (2017); Riccardi  _et al._ (2020) consisting of
200,000 cycles. In all simulations (Shooting, WF*, WF, WF-cap), each path
ensemble either employs only shooting or only WF as main MC move in addition
to replica exchange moves. In simulation “Shooting” all path ensembles employ
the shooting move. In the other simulations the WF move is used for most path
ensembles. However, simulation WF* uses normal shooting in the $[0^{-}]$
ensemble, while simulations WF and WF-cap use the shooting move in both the
$[0^{-}]$ and $[0^{+}]$ ensemble as was suggested in Sec. III. The WF-cap
simulation uses a cap-interface at $\lambda_{\rm cap}$ = 0.1 At each cycle,
all path ensembles are updated with an ensemble move (shooting or WF) or with
replica exchange moves with a 50%-50% probability. In case that a replica
exchange move is selected, another 50%-50% probability determines whether the
$[0^{-}]\leftrightarrow[0^{+}]$, $[1^{+}]\leftrightarrow[2^{+}]$, $\ldots$,
$[5^{+}]\leftrightarrow[6^{+}]$ swaps will be attempted or the
$[0^{+}]\leftrightarrow[1^{+}]$, $[2^{+}]\leftrightarrow[3^{+}]$, $\ldots$,
$[4^{+}]\leftrightarrow[5^{+}]$ swaps. In the latter case, the $[0^{-}]$ and
$[6^{+}]$ ensembles simply duplicate the previous path (null move). In the WF
simulations, the number of subpaths was arbitrarily set equal to $N_{s}=6$ for
all path ensembles.

The results are shown in Fig. 5 and in table 1 where they are compared with
Kramers’ theory Hanggi  _et al._ (1990) which, for this system, can be
considered as a nearly exact reference. Fig. 5 shows that the WF based
simulations rapidly converge close to the Kramers’ value of the rate
confirming the exactness of the superdetailed balance relations and the
correct implementation in the PyRETIS code. The results based on shooting are
further off, but have a significant lower computational cost per RETIS cycle
(see table 1).

$0$$50$$100$$150$$200$$0$$1$$2$$3$RETIS
cycles/$10^{3}$Rate/$10^{-7}$KramersWF*WFWF-capShooting Figure 5: Total
running average of computed rate as function of RETIS cycles.

Based on the relative errors from the block averaging analysis and the cost
per cycle, we can compute the CPU efficiency time for each method, shown in
the last column. Based on these numbers, we can see that the WF*, WF and WF-
cap simulations are 2.4, 2.5 and 2.7 times more efficient than the simulation
in which all path ensembles use the standard shooting move as their main MC
move. Note that an improvement of more than a factor 2 is rather remarkable
given the low dimensionality of the system.

Table 1: Simulation data for the double well 1D barrier system. The cost column describes the total number of calculated MD steps. The errors are based on block averaging using single standard deviations. The final column shows the CPU efficiency times van Erp (2006) corresponding to the number of required MD steps for obtaining a relative error equal to 1. Simulation “Shooting” uses the standard shooting move as main MC move in all path ensembles. The other simulations use the WF move in all ensembles except for $[0^{-}]$ (WF*, WF, WF-cap) and $[0^{+}]$ (WF, WF-cap). WF-cap uses a cap-interface at $\lambda_{\rm cap}=0.1.$ Simulation | Rate$/10^{-7}$ |  $\epsilon_{r}$ (%) | Cost/$10^{7}$ | Cost$\cdot\epsilon_{r}^{2}$/$10^{11}$
---|---|---|---|---
Shooting | 2.30 | 6.46 | 5.32 | 222.0
WF* | 2.58 | 2.19 | 19.56 | 93.9
WF | 2.69 | 2.28 | 16.98 | 88.3
WF-cap | 2.54 | 2.29 | 15.72 | 82.4
Kramers | 2.58 |  |  |

In table 2, we further examine the acceptance probabilities of the different
moves. It is apparent that in all simulations the main MC move has a nearly
100% acceptance in the path ensembles where the WF move is employed thanks to
the high-acceptance protocol. The acceptance is marginally lower at the last
path ensembles $[5^{+}]$ and $[6^{+}]$ from which there is a higher
probability to generate $\lambda_{B}\rightarrow\lambda_{B}$ paths. The
shooting move has a lower acceptance, but has the advantage that all swapping
moves with the $[0^{-}]$ ensemble are accepted if shooting is the main move in
both $[0^{-}]$ and $[0^{+}]$. Since $[0^{-}]$ can only swap with $[0^{+}]$,
these are the computational expensive $[0^{-}]\leftrightarrow[0^{+}]$ swaps.
The other swapping moves are inexpensive as they do not require any MD steps.

Therefore, an anticipated lower acceptance for these swapping moves in the WF
simulations would not be dramatic. However, even this is not always the case.
At first sight this appears counter-intuitive. Given a pair of paths in two
neighboring ensembles, the standard swap should always have an acceptance
probability that is equal or higher than the acceptance based on Eq. 31.
However, this effect can be canceled by the path distributions not being the
same. Since the altered path distribution in the high-acceptance scheme, Eq.
28, overrepresents paths with many points between $\lambda_{i}$ and
$\lambda_{B}$ or $\lambda_{\rm cap}$, the $[i^{+}]$ path ensemble is likely to
contain a higher fraction of paths crossing $\lambda_{i+1}$. From the data of
table 2, this seems indeed the case in the majority of path ensembles.

Table 2: Acceptance ratio (%). Simulation “Shooting” uses the standard shooting move as main MC move in all path ensembles. The other simulations use the WF move in all ensembles except for $[0^{-}]$ (WF*, WF, WF-cap) and $[0^{+}]$ (WF, WF-cap).  ens. | Shooting | WF* | WF | WF-cap
---|---|---|---|---
main | swap | main | swap | main | swap | main | swap
$[0^{-}]$ | 84.6 | 100.0 | 84.5 | 83.9 | 84.3 | 100.0 | 84.3 | 100.0
$[0^{+}]$ | 84.2 | 57.8 | 100.0 | 49.0 | 84.0 | 55.6 | 84.0 | 55.8
$[1^{+}]$ | 48.8 | 15.5 | 100.0 | 17.9 | 100.0 | 16.3 | 100.0 | 16.8
$[2^{+}]$ | 37.8 | 13.4 | 100.0 | 19.7 | 100.0 | 19.9 | 100.0 | 20.2
$[3^{+}]$ | 32.2 | 11.5 | 100.0 | 18.0 | 100.0 | 18.5 | 100.0 | 18.4
$[4^{+}]$ | 30.1 | 12.2 | 100.0 | 20.3 | 100.0 | 20.6 | 100.0 | 20.4
$[5^{+}]$ | 30.0 | 14.7 | 99.9 | 28.3 | 99.8 | 28.2 | 100.0 | 26.6
$[6^{+}]$ | 29.1 | 16.7 | 99.2 | 34.2 | 99.2 | 33.9 | 100.0 | 30.8

###  VIII.2 Thin Film Breakage

A system of 1100 dodecane molecules layered on a slab of 23936 water molecules
is studied in the NPT ensemble via full atom TIS simulations using the GROMACS
2020.1 simulation package Abraham  _et al._ (2015b) as external engine. The
dodecane molecules are simulated according to the OPLS-AA force field
Jorgensen  _et al._ (1996) and the water molecules with the TIP4p/2005 model
Abascal and Vega (2005). The preparation of the initial equilibrated system is
explained in detail by Ref. Aarøen  _et al._ , 2022. The temperature is set to
300 K and is controlled with a velocity rescaling method Bussi  _et al._
(2007) employing a coupling time of 0.1 ps. Pressure is controlled by the
Berendsen barostat and its normal component is maintained constant at 1 bar,
with a time constant of 1.0 ps and compressibility coefficient of $4.7\cdot
10^{-5}\mathrm{bar}^{-1}$. The velocity-Verlet algorithm is used to solve the
Newton equations of motion with a timestep of 0.002 ps. Periodic boundary
conditions are applied in all directions, with the z direction being
perpendicular to the 2D film. The box size is set to equal box size of
$15\times 15\times 5.1983$ nm.

The order parameter of the system is calculated through discretizing the
system into $85\times 85$ tiles in the $x$ and $y$ direction such that the
order parameter value becomes the number of empty dodecane tiles that also
have empty neighbors in the x and y direction. Such a definition provides a
way to measure the presence of low-density regions, in addition to any
breakage or “hole” formation that occurs within a trajectory. The sensitivity
of the order parameter is determined by the specified discretizing size. In
our case, the order parameter values fluctuated between 0 and 5 during an
equilibrium run at $T=300$ K. Based on this, we set $\lambda_{A}=5$. We
further defined $\lambda_{B}=100$ as preliminary analysis showed that from
this point on the hole tends to grow further with negligible chance to close
again.

To obtain an initial reactive trajectory, we ran an equilibrium run at $T=375$
K until the thin-film broke down. For our single path ensemble analysis we
further defined $\lambda_{\mathrm{i}}=10.0$ as the interface that has to be
crossed. In addition, we set the cap-interface $\lambda_{\mathrm{cap}}=15.0$.
We then created 1000 trajectories using standard shooting and WF with
$N_{s}=10$. Three exemplary trajectories from the WF simulation are shown in
Fig. 6-a).

$0$$20$$40$$60$$80$$100$$120$$0$$20$$40$$60$$80$$100$Time [ps]Order
parameter$\lambda_{\mathrm{A/B}}$$\lambda_{\mathrm{cap}}$$\lambda_{\mathrm{i}}$a)

$0$$50$$100$$150$$200$$250$$300$$-1$$-0.5$$0$$0.5$$1$Time [fs]Order
parameter$\lambda_{\mathrm{A/B}}$$\lambda_{\mathrm{cap}}$$\lambda_{\mathrm{i}}$b)

Figure 6:  Exemplary trajectories from the WF algorithm in the $[i^{+}]$ path
ensemble showing the progress of the order parameter versus time. The stable
state interfaces $\lambda_{A}$, $\lambda_{B}$, the cap-interface $\lambda_{\rm
cap}$, and the ensemble interface $\lambda_{i}$ are shown as well. The two
different panels represent the a) classical MD level simulation of the thin
film breakage and b) the AIMD level simulations of the ruthenium self-exchange
reaction.

From the sample size of 1000 MC moves, the acceptance in WF was equal to 73.4%
and 35.0% for standard shooting. The limited sample size prohibit accurate CPU
efficiency analysis, but a qualitative assertion of the sampling effectivity
can be obtained by viewing the simulated path lengths as function of the MC
step.

$0$$100$$200$$\mathrm{Shooting}$$0$$200$$400$$600$$800$$1{,}000$$0$$100$$200$Shooting/WF
movesPath length [ps]WFa)

$0$$200$$400$$600$$800$$\mathrm{Shooting}$$0$$20$$40$$60$$80$$100$$0$$200$$400$$600$$800$Shooting/WF
movesPath length [fs]WFb)

Figure 7: Path length vs MC move for WF and standard shooting for a) classical
MD system of thin film breakage and b) AIMD system of ruthenium self-exchange
reaction.

Fig. 7-a) shows that the WF sampling has much more frequent transitions
between long and short paths whereas shooting is mostly stuck in the short
path domain. Once the shooting move manages to produce a long path, the path
remains in the MC chain due to a long series of rejections (e.g. around step
500 where the same path length remains for a number of steps due to
rejections). This indicates that the shooting move is struggling to properly
sample path space. Even if the acceptance is not extremely low for the short
paths, it fails to make regular switches to the longer paths. Moreover, if a
long path is generated, the subsequent moves are likely rejected such that
other longer paths are not likely found.

###  VIII.3 Ruthenium-Ruthenium Self-Exchange Reaction

We studied self-exchange reaction between two ruthenium ions in aqueous
solution described by the following chemical reaction

$\mathrm{Ru^{2+}}+\mathrm{Ru^{3+}}\rightarrow\mathrm{Ru^{3+}}+\mathrm{Ru^{2+}}$
(36)

The simulation system consisted of two ruthenium ions, 63 $\mathrm{H_{2}O}$
molecules and one OH- ion. The dynamics were propagated using NVE velocity-
Verlet and the CP2K Hutter  _et al._ (2014) simulation package. The effect of
temperature was introduced via the randomization of velocities from a Maxwell-
Boltzmann distribution at a temperature of 300 K. We used a time step of 0.5
fs and periodic boundary conditions were applied to a cubic box with an edge
length of 12.4138 Å. Further simulation details on functional and basis sets
are explained in Ref. Tiwari and Ensing, 2016.

To monitor the reaction progress, the electron transfer has been “followed” by
transforming the occupied Kohn-Sham orbitals Kohn and Sham (1965) into
maximally localized Wannier functions (MLWF) Marzari  _et al._ (2012) and
computing the distance between the center of these localized functions (X)
describing the moving electron to each of the ruthenium ions. The order
parameter of the system is then defined as

$\lambda=\frac{\left(d_{\mathrm{Ru}-X}-d_{\mathrm{Ru}^{\prime}-X}\right)}{d_{\mathrm{Ru}-\mathrm{Ru}^{\prime}}}$
(37)

where $d_{\mathrm{Ru}-X}$ is the distance between X and the initial ruthenium
electron donor, $d_{\mathrm{Ru}^{\prime}-X}$ is the distance between X and the
initial ruthenium electron acceptor and $d_{\mathrm{Ru}-\mathrm{Ru}^{\prime}}$
is the distance between the two ruthenium ions in the system. In this
formulation, $\lambda_{\mathrm{}}=-1$ and $\lambda_{\mathrm{}}=+1$ define the
reactant state and product state, respectively.
$\mathrm{Ru^{2+}}$/$\mathrm{Ru^{3+}}$ have 5/6 d-electrons and
$\mathrm{H_{2}O}$/$\mathrm{OH^{-}}$ have 8 valence electrons. This means there
are a total of 523 MLWFs in the system. The order parameter, Eq. 37, requires
the location X of the transferring electron, which is one of the centers of
these 523 MLWFs. To identify which is X, each Wannier center is being linked
to either a ruthenium or oxygen atom that is closest. Then, if one ruthenium
ion has 6 associated MLWFs, X is set to be the one that is the farthest away
from this ruthenium ion. In the case that both ruthenium ions have 5
associated MLWFs, one of the oxygens has an excess MLWF (9 instead of 8), and
the center that is farthest away from this oxygen is set as X.

To qualitatively compare standard shooting and WF for this system we run two
single path ensemble simulations representing $[i^{+}]$ with
$\lambda_{i}=-0.736$, $\lambda_{A}=-0.99$, and $\lambda_{B}=+0.99$. The value
for $\lambda_{i}=-0.736$ was chosen from preliminary runs where we aimed for a
20% probability that a path ends up at state $B$. In the WF simulation, an
additional $\lambda_{\rm cap}=+0.95$ was set to avoid
$\lambda_{B}\rightarrow\lambda_{B}$ rejections due to the selection of
shooting points lying within the basin of attraction of state $B$. Here, we
only applied a rather modest number of subtrajectories $N_{s}=2$. Higher
performances might be obtained with a larger number of subpaths. Exemplary
trajectories of the WF simulation are show in Fig 6-b).

Due to the relatively low value of $N_{s}$, the subpath contribution to the
total WF computational cost is only 15%. The acceptance probability increased
from the shooting move’s 48% to WF’s 96%. Similarly to the classical MD
system, the WF simulation seems to show a better sample exploration when we
look at the path length as function of the MC step (Fig. 7-b)). The standard
shooting algorithm seems not to be able to produce any paths larger than 300
fs. The WF algorithm, however, started with a short initial path but was able
to quickly move up to to the 600 fs range and making regular transitions
between the shorter and longer paths. So also here, the sampling quality of
the WF algorithm appears substantially superior to the one of standard
shooting.

##  IX Concluding Remarks

We reviewed the recently developed subtrajectory moves stone skipping (SS) and
web throwing (WT) and added a new member to this group: wire fencing (WF).
These moves are more efficient than the standard shooting move which has been
the main MC move for path sampling simulations during the last two decades.
The subtrajectory moves proceed from a complete old path to a complete new
path via a series of intermediate short paths (subpaths/subtrajectories).
While this increases the average cost of a MC step, the correlations between
paths are substantially reduced leading to a lower statistical inefficiency.
The use of shorter paths resembles approximate path sampling methods like
PPTIS or milestoning. However, the subtrajectory moves are still exact like
standard shooting as they are based on mathematically rigorous superdetailed
balance relations. The approach is preferably combined with a high-acceptance
protocol in which the sampling distribution of the paths is adjusted in order
to maximize the acceptance of newly generated trajectories. The effect of the
biased distribution is undone in the post-simulation analysis using
appropriate reweighting. The SS and WT move, however, require a one-step
crossing condition which complicates their implementation and we discussed
several solutions for this issue. The new WF does not rely on the one-step
crossing condition and is therefore the most practical solution to the
aforementioned problem even if it is slightly more wasteful than SS and WT.
The WF move is in particular useful when the path sampling code uses an
external MD engine and/or when the computation of the order parameter is
costly. We showed the exactness and the efficiency gain of the WF approach in
a RETIS simulation where the transition rate of an underdamped Langevin
particle in a double well potential has been computed and compared with the
analytical Kramers’ expression. Thereafter, we showed qualitatively how the WF
move performs in a classical MD system describing the thin film breaking
process and in an AIMD level system describing an electron transfer process
between ruthenium ions in aqueous solution. In both cases, the WF move seems
to allow a faster sampling through path space than standard shooting, which
was concluded from the rapid switches that WF made between the shorter and
longer paths.

###### Acknowledgements.

We acknowledge funding from the Research Council of Norway (toppforsk project
Theolight, grant no. 275506) and computational resources from NOTUR (project
nn9254k).

## Author Declarations

The authors have no conflicts to disclose.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

The algorithmic developments have been implemented in the current pyretis2.dev
version (the current release is PyRETIS-2 Riccardi  _et al._ (2020)) and will
be included in the forthcoming main release (PyRETIS-3). The code is already
available at https://gitlab.com/pyretis following the FAIR principle for
scientific data and software and data. Riccardi  _et al._ (2019); Lamprecht
_et al._ (2020)

##  Appendix A Analytical expressions for the statistical inefficiency in
model systems

Sec. IV introduces a model where at each MC move $j$ there is a chance of
$\pi_{R}$ that the state of the system remains essentially unchanged and a
chance of $\pi_{M}=1-\pi_{R}$ to “throw a dice”. The latter implies that at
step $j$ the output value of $h_{j}$ equals 1 with a probability $p$ and 0
with a probability $1-p$. Let us consider the conditional probability that
$h_{j}=0$ given that $h_{0}=0$: $P(h_{j}=0|h_{0}=0)$. We can distinguish two
scenarios. Scenario 1 relates to the case that all $j$ moves implied a
“remain” and therefore $h_{j}=h_{0}=0$. Scenario 2 is related to the situation
that at least once the dice was thrown. In this scenario we have that $h_{j}$
is either 1 or 0 with respective probabilities $p$ and $1-p$. The probability
of having scenario 1 equals $\pi_{R}^{j}$ and that of scenario 2 equals
$1-\pi_{R}^{j}$. Therefore,

$\displaystyle P(h_{j}=0|h_{0}=0)$
$\displaystyle=\pi_{R}^{j}+(1-\pi_{R}^{j})(1-p)$
$\displaystyle=(1-p)+p\pi_{R}^{j}$ (38)

Likewise, we can derive all the other conditional probabilities:

$\displaystyle P(h_{j}=1|h_{0}=0)$
$\displaystyle=(1-\pi_{R}^{j})p=p-p\pi_{R}^{j}$ $\displaystyle
P(h_{j}=0|h_{0}=1)$
$\displaystyle=(1-\pi_{R}^{j})(1-p)=(1-p)+(1-p)\pi_{R}^{j}$ $\displaystyle
P(h_{j}=1|h_{0}=1)$
$\displaystyle=\pi_{R}^{j}+(1-\pi_{R}^{j})p=p+(1-p)\pi_{R}^{j}$ (39)

Let us call $p_{kl}=P(h_{j}=k\wedge h_{0}=l)=P(h_{j}=k|h_{0}=l)P(h_{0}=l)$.
From Eqs. 38 and 39 we can derive:

$\displaystyle p_{00}$ $\displaystyle=(1-p)^{2}+p(1-p)\pi_{R}^{j}$
$\displaystyle p_{10}$ $\displaystyle=p(1-p)-p(1-p)\pi_{R}^{j}=p_{01}$
$\displaystyle p_{11}$ $\displaystyle=p^{2}+p(1-p)\pi_{R}^{j}$ (40)

and from this we can compute

$\displaystyle\left\langle(h_{0}-p)(h_{j}-p)\right\rangle$
$\displaystyle=p_{00}p^{2}-p_{10}(1-p)p$ $\displaystyle-
p_{01}(1-p)p+p_{11}(1-p)^{2}$ (41)

In the above expression, all the $\pi_{R}$-independent terms cancel. This is
expected since we know the result equals 0 if $\pi_{R}=0$. The remaining
$\pi_{R}$-dependent terms sum up to

$\displaystyle p(1-p)\pi_{R}^{j}\left[p^{2}+2p(1-p)+(1-p)^{2}\right]$
$\displaystyle=p(1-p)\pi_{R}^{j}\left[p+(1-p)\right]^{2}=p(1-p)\pi_{R}^{j}$
(42)

From Eqs. 16, 17, and 42 we derive that

$\displaystyle C(j)=\pi_{R}^{j}\Rightarrow n_{c}=\frac{\pi_{R}}{1-\pi_{R}}$
(43)

and via Eq. 15:

$\displaystyle{\mathcal{N}}=1+2\frac{\pi_{R}}{1-\pi_{R}}=\frac{1+\pi_{R}}{1-\pi_{R}}$
(44)

As $\pi_{M}=1-\pi_{R}$, 44 is equivalent to Eq. 19 of Sec. IV.

In the second model we assume $\pi_{R}=0$, but there are two phases
$x=\alpha,\beta$ that have have, respectively, probabilities $P_{\alpha}$ and
$P_{\beta}$ and local crossing probabilities $p_{\alpha}$ and $p_{\beta}$. Let
$\pi_{\rho}=1-\pi_{\mu}$ be the probability that the MC maintains the previous
phase. The inverse probability $\pi_{\mu}$ implies throwing the dice to
determine the phase $x$ such that the selection probability for $x$
corresponds to $P_{\alpha}$ and $P_{\beta}=1-P_{\alpha}$. After the phase $x$
is set, $h_{j}$ will be set to 1 or 0 with respective probabilities $p_{x}$
and $(1-p_{x})$. Given that the phase of the first sample equals $x_{0}=x$,
the chance that the $j$-th sample has the same or opposite phase equals,
respectively, $\pi_{\rho}^{j}+(1-\pi_{\rho}^{j})P_{x}$ and
$(1-\pi_{\rho}^{j})(1-P_{x})$. This leads to following conditional
probabilities akin Eqs. 38 and 39:

$\displaystyle P(h_{j}=0|x_{0}=x)$
$\displaystyle=\pi_{\rho}^{j}(1-p_{x})+(1-\pi_{\rho}^{j})(1-p)$
$\displaystyle=\pi_{\rho}^{j}(p-p_{x})+(1-p)$
$\displaystyle=\pi_{\rho}^{j}P_{y}(p_{y}-p_{x})+(1-p)$ $\displaystyle
P(h_{j}=1|x_{0}=x)$ $\displaystyle=\pi_{\rho}^{j}p_{x}+(1-\pi_{\rho}^{j})p$
$\displaystyle=\pi_{\rho}^{j}(p_{x}-p)+p$
$\displaystyle=\pi_{\rho}^{j}P_{y}(p_{x}-p_{y})+p$ (45)

where $y\in(\alpha,\beta)$ and $y\neq x$. Hence, analogous to Eqs. 40

$\displaystyle p_{k0}$
$\displaystyle=\sum_{x=\alpha,\beta}P_{x}(1-p_{x})P(h_{j}=k|x_{0}=x)$
$\displaystyle p_{k1}$
$\displaystyle=\sum_{x=\alpha,\beta}P_{x}p_{x}P(h_{j}=k|x_{0}=x)$

which leads to

$\displaystyle p_{00}$
$\displaystyle=(1-p)^{2}+\pi_{\rho}^{j}\sum_{x}P_{x}P_{y}(1-p_{x})(p_{y}-p_{x})$
$\displaystyle=(1-p)^{2}+\pi_{\rho}^{j}P_{\alpha}P_{\beta}(p_{\alpha}-p_{\beta})^{2}$
$\displaystyle p_{10}$
$\displaystyle=p(1-p)+\pi_{\rho}^{j}\sum_{x}P_{x}P_{y}(1-p_{x})(p_{x}-p_{y})$
$\displaystyle=p(1-p)-\pi_{\rho}^{j}P_{\alpha}P_{\beta}(p_{\alpha}-p_{\beta})^{2}=p_{01}$
$\displaystyle p_{11}$
$\displaystyle=p^{2}+\pi_{\rho}^{j}\sum_{x}P_{x}P_{y}p_{x}(p_{x}-p_{y})$
$\displaystyle=p^{2}+\pi_{\rho}^{j}P_{\alpha}P_{\beta}(p_{\alpha}-p_{\beta})^{2}$
(47)

Analogous to Eqs. 41 and 42 we find that

$\displaystyle\left\langle(h_{0}-p)(h_{j}-p)\right\rangle$
$\displaystyle=\pi_{\rho}^{j}P_{\alpha}P_{\beta}(p_{\alpha}-p_{\beta})^{2}$
(48)

and like Eq. 43:

$\displaystyle C(j)$
$\displaystyle=\frac{P_{\alpha}P_{\beta}(p_{\alpha}-p_{\beta})^{2}}{p(1-p)}\pi_{\rho}^{j}=K_{s}\pi_{\rho}^{j}$
$\displaystyle\Rightarrow n_{c}$
$\displaystyle=K_{s}\frac{\pi_{\rho}}{1-\pi_{\rho}}=K_{s}\frac{1-\pi_{\mu}}{\pi_{\mu}}$
(49)

where we used $\pi_{\mu}=1-\pi_{\rho}$ and Eq. 21. Substition of Eq. 49 in Eq.
15 leads to Eq. 20.

## References

  * Lindorff-Larsen  _et al._ (2011) K. Lindorff-Larsen, S. Piana, R. O. Dror, and D. E. Shaw, Science 334, 517 (2011).
  * Shaw  _et al._ (2021) D. E. Shaw, P. J. Adams, A. Azaria, J. A. Bank, B. Batson, A. Bell, M. Bergdorf, J. Bhatt, J. A. Butts, T. Correia, R. M. Dirks, R. O. Dror, M. P. Eastwood, B. Edwards, A. Even, P. Feldmann, M. Fenn, C. H. Fenton, A. Forte, J. Gagliardo, G. Gill, M. Gorlatova, B. Greskamp, J. Grossman, J. Gullingsrud, A. Harper, W. Hasenplaugh, M. Heily, B. C. Heshmat, J. Hunt, D. J. Ierardi, L. Iserovich, B. L. Jackson, N. P. Johnson, M. M. Kirk, J. L. Klepeis, J. S. Kuskin, K. M. Mackenzie, R. J. Mader, R. McGowen, A. McLaughlin, M. A. Moraes, M. H. Nasr, L. J. Nociolo, L. O’Donnell, A. Parker, J. L. Peticolas, G. Pocina, C. Predescu, T. Quan, J. K. Salmon, C. Schwink, K. S. Shim, N. Siddique, J. Spengler, T. Szalay, R. Tabladillo, R. Tartler, A. G. Taube, M. Theobald, B. Towles, W. Vick, S. C. Wang, M. Wazlowski, M. J. Weingarten, J. M. Williams, and K. A. Yuh, in  _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_, SC ’21 (Association for Computing Machinery, New York, NY, USA, 2021).
  * Goldberg  _et al._ (1990) M. E. Goldberg, G. V. Semisotnov, B. Friguet, K. Kuwajima, O. B. Ptitsyn, and S. Sugai, FEBS Lett. 263, 51 (1990).
  * Peters (2017) B. Peters,  _Reaction rate theory and rare events_ (Elsevier, Amsterdam, Netherlands, 2017).
  * van Erp  _et al._ (2003) T. S. van Erp, D. Moroni, and P. G. Bolhuis, J. Chem. Phys. 118, 7762 (2003).
  * van Erp (2007) T. van Erp, Phys. Rev. Lett. 98, 268301 (2007).
  * Dellago  _et al._ (1998a) C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, J. Chem. Phys. 108, 1964 (1998a).
  * Cabriolu  _et al._ (2017) R. Cabriolu, K. M. S. Refsnes, P. G. Bolhuis, and T. S. van Erp, J. Chem. Phys. 147, 152722 (2017).
  * Arjun and Bolhuis (2020) A. Arjun and P. G. Bolhuis, J. Phys. Chem. B 124, 8099 (2020).
  * Moqadam  _et al._ (2018) M. Moqadam, A. Lervik, E. Riccardi, V. Venkatraman, B. K. Alsberg, and T. S. van Erp, Proc. Natl. Acad. Sci. USA 115, E4569 (2018).
  * Eigen and de Maeyer (1958) M. Eigen and L. de Maeyer, Proc. R. Soc. Lond. A Math. Phys. Sci. 247, 505 (1958).
  * Natzle and Moore (1985) W. C. Natzle and C. B. Moore, J. Phys. Chem. 89, 2605 (1985).
  * Moroni  _et al._ (2003) D. Moroni, P. Bolhuis, and T. van Erp, J. Chem. Phys. 120, 1044 (2003).
  * Faradjian and Elber (2004) A. K. Faradjian and R. Elber, J. Chem. Phys. 120, 10880 (2004).
  * Roet  _et al._ (2022) S. Roet, D. T. Zhang, and T. S. van Erp, “Exchanging replicas with unequal cost, infinitely and permanently. Preprint at https://arxiv.org/abs/2205.12663,”  (2022).
  * Dellago  _et al._ (1998b) C. Dellago, P. G. Bolhuis, and D. Chandler, The Journal of Chemical Physics 108, 9236 (1998b), https://doi.org/10.1063/1.476378 .
  * Riccardi  _et al._ (2017) E. Riccardi, O. Dahlen, and T. S. van Erp, J. Phys. Chem. Lett. 8, 4456 (2017).
  * Swenson  _et al._ (2019a) D. W. H. Swenson, J.-H. Prinz, F. Noe, J. D. Chodera, and P. G. Bolhuis, J. Chem. Theory Comput. 15, 813 (2019a).
  * Swenson  _et al._ (2019b) D. W. H. Swenson, J.-H. Prinz, F. Noe, J. D. Chodera, and P. G. Bolhuis, J. Chem. Theory Comput. 15, 837 (2019b).
  * Lervik  _et al._ (2017) A. Lervik, E. Riccardi, and T. S. van Erp, J. Comput. Chem. 38, 2439 (2017).
  * Riccardi  _et al._ (2020) E. Riccardi, A. Lervik, S. Roet, O. Aaroen, and T. S. van Erp, J. Comput. Chem. 41, 370 (2020).
  * Abraham  _et al._ (2015a) M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, SoftwareX 1-2, 19 (2015a).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Eastman  _et al._ (2013) P. Eastman, M. S. Friedrichs, J. D. Chodera, R. J. Radmer, C. M. Bruns, J. P. Ku, K. A. Beauchamp, T. J. Lane, L.-P. Wang, D. Shukla, T. Tye, M. Houston, T. Stich, C. Klein, M. R. Shirts, and V. S. Pande, J. Chem. Theory Comput. 9, 461 (2013).
  * Hutter  _et al._ (2014) J. Hutter, M. Iannuzzi, F. Schiffmann, and J. VandeVondele, WileyWIREs Comput Mol Sci 4, 15 (2014).
  * Stukowski (2012) A. Stukowski, Model. Simul. Mater. Sci. Eng. 20, 045021 (2012).
  * Winczewski  _et al._ (2016) S. Winczewski, J. Dziedzic, and J. Rybicki, Comput. Phys. Commun. 198, 128 (2016).
  * Metropolis  _et al._ (1953) N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
  * Frenkel and Smit (2002) D. Frenkel and B. Smit,  _Understanding molecular simulations from algorithms to applications_ (Academic press, San Diego, California, U.S.A., 2002).
  * Siepmann and Frenkel (1992) J. I. Siepmann and D. Frenkel, Mol. Phys. 75, 59 (1992).
  * Vlugt  _et al._ (1999) T. Vlugt, R. Krishna, and B. Smit, J. Phys. Chem. B 103, 1102 (1999).
  * Hastings (1970) W. Hastings, Biometrika 57, 97 (1970).
  * van Erp (2012) T. van Erp, Adv. Chem. Phys. 151, 27 (2012).
  * Ghysels  _et al._ (2021) A. Ghysels, S. Roet, S. Davoudi, and T. S. van Erp, Phys. Rev. Research 3, 033068 (2021).
  * Dellago and Bolhuis (2004) C. Dellago and P. G. Bolhuis, Mol. Simu. 30, 795 (2004).
  * van Erp and Bolhuis (2005) T. van Erp and P. Bolhuis, J. Comput. Phys. 205, 157 (2005).
  * van Erp (2006) T. S. van Erp, J. Chem. Phys. 125, 174106 (2006).
  * Ferrenberg and Swendsen (1989) A. Ferrenberg and R. Swendsen, Phys. Rev. Lett. 63, 1195 (1989).
  * van Erp  _et al._ (2016) T. S. van Erp, M. Moqadam, E. Riccardi, and A. Lervik, J. Chem. Theory Comput. 12, 5398 (2016).
  * Rogal  _et al._ (2010) J. Rogal, W. Lechner, J. Juraszek, B. Ensing, and P. G. Bolhuis, J. Comp. Phys. 133, 174109 (2010).
  * Vanden-Eijnden  _et al._ (2008) E. Vanden-Eijnden, M. Venturoli, G. Ciccotti, and R. Elber, J. Comp. Phys. 129, 174102 (2008).
  * Haji-Akbari (2018) A. Haji-Akbari, J. Chem. Phys. 149, 072303 (2018).
  * Brotzakis and Bolhuis (2016) Z. F. Brotzakis and P. G. Bolhuis, J. Chem. Phys. 145, 164112 (2016).
  * Moqadam  _et al._ (2017) M. Moqadam, E. Riccardi, T. T. Trinh, A. Lervik, and T. S. van Erp, Phys. Chem. Chem. Phys. 19, 13361 (2017).
  * Moroni  _et al._ (2005) D. Moroni, P. R. ten Wolde, and P. G. Bolhuis, Phys. Rev. Lett. 94, 235703 (2005).
  * Aarøen  _et al._ (2022) O. Aarøen, E. Riccardi, T. S. v. Erp, and M. Sletmoen, Colloids and Surfaces A: Physicochemical and Engineering Aspects 632, 127808 (2022).
  * Aarøen  _et al._ (2021) O. Aarøen, E. Riccardi, and M. Sletmoen, RSC advances 11, 8730 (2021).
  * Riccardi and Tichelkamp (2019) E. Riccardi and T. Tichelkamp, Colloids and Surfaces A: Physicochemical and Engineering Aspects 573, 246 (2019).
  * Riccardi  _et al._ (2014) E. Riccardi, K. Kovalchuk, A. Y. Mehandzhiyski, and B. A. Grimes, Journal of dispersion science and technology 35, 1018 (2014).
  * Marzari  _et al._ (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Allen  _et al._ (2009) R. J. Allen, C. Valeriani, and P. R. ten Wolde, J. Phys.-Condes. Matter 21, 463102 (2009).
  * Escobedo  _et al._ (2009) F. A. Escobedo, E. E. Borrero, and J. C. Araque, J. Phys.-Condes. Matter 21, 333101 (2009).
  * Booth and Hendricks (1984) T. E. Booth and J. S. Hendricks, Nucl. Techno.-Fus. 5, 90 (1984).
  * Melnik-Melnikov and Dekhtyaruk (2000) P. Melnik-Melnikov and E. Dekhtyaruk, Probab. Eng. Eng. Mech. 15, 125 (2000).
  * Villenaltamirano and Villenaltamirano (1991) M. Villenaltamirano and J. Villenaltamirano, in  _Queueing, Performance and Control in Atm_ , North-Holland Studies in Telecommunication, Vol. 15, edited by J. W. Cohen and C. D. Pack (Elsevier Science Publ B V, Amsterdam, 1991) pp. 71–76, 13th International Teletraffic Congress ( ITC-13 ), Copenhagen, Denmark, Jun 19-26, 1991.
  * Hanggi  _et al._ (1990) P. Hanggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
  * Abraham  _et al._ (2015b) M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, SoftwareX 1, 19 (2015b).
  * Jorgensen  _et al._ (1996) W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, Journal of the American Chemical Society 118, 11225 (1996).
  * Abascal and Vega (2005) J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).
  * Bussi  _et al._ (2007) G. Bussi, D. Donadio, and M. Parrinello, J. Chem. Phys. 126, 014101 (2007).
  * Tiwari and Ensing (2016) A. Tiwari and B. Ensing, Faraday Discussions 195, 291 (2016).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Riccardi  _et al._ (2019) E. Riccardi, S. Pantano, and R. Potestio, Interface Focus 9, 20190005 (2019).
  * Lamprecht  _et al._ (2020) A.-L. Lamprecht, L. Garcia, M. Kuzak, C. Martinez, R. Arcila, E. Martin Del Pico, V. Dominguez Del Angel, S. Van De Sandt, J. Ison, P. A. Martinez, _et al._ , Data Science 3, 37 (2020).

TOC