# Exchanging Replicas with Unequal Cost, Infinitely and Permanently

- Sander RoetSander RoetDepartment of Chemistry, Norwegian University of Science and Technology (NTNU), N-7491Trondheim, NorwayMore by Sander Roet
- ,
- Daniel T. ZhangDaniel T. ZhangDepartment of Chemistry, Norwegian University of Science and Technology (NTNU), N-7491Trondheim, NorwayMore by Daniel T. Zhang
- , and
- Titus S. van Erp
*****Titus S. van ErpDepartment of Chemistry, Norwegian University of Science and Technology (NTNU), N-7491Trondheim, Norway*****Email: [email protected]More by Titus S. van Erp

## Abstract

We developed a replica exchange method that is effectively parallelizable even if the computational cost of the Monte Carlo moves in the parallel replicas are considerably different, for instance, because the replicas run on different types of processor units or because of the algorithmic complexity. To prove detailed-balance, we make a paradigm shift from the common conceptual viewpoint in which the set of parallel replicas represents a high-dimensional superstate, to an ensemble-based criterion in which the other ensembles represent an environment that might or might not participate in the Monte Carlo move. In addition, based on a recent algorithm for computing permanents, we effectively increase the exchange rate to infinite without the steep factorial scaling as a function of the number of replicas. We illustrate the effectiveness of this replica exchange methodology by combining it with a quantitative path sampling method, replica exchange transition interface sampling (RETIS), in which the costs for a Monte Carlo move can vary enormously as paths in a RETIS algorithm do not have the same length and the average path lengths tend to vary considerably for the different path ensembles that run in parallel. This combination, coined ∞RETIS, was tested on three model systems.

This publication is licensed under

### License Summary*

You are free to share (copy and redistribute) this article in any medium or format and to adapt (remix, transform, and build upon) the material for any purpose, even commercially within the parameters below:

Creative Commons (CC): This is a Creative Commons license.

Attribution (BY): Credit must be given to the creator.

*Disclaimer

This summary highlights only some of the key features and terms of the actual license. It is not a license and has no legal value. Carefully review the actual license before using these materials.

### License Summary*

You are free to share (copy and redistribute) this article in any medium or format and to adapt (remix, transform, and build upon) the material for any purpose, even commercially within the parameters below:

Creative Commons (CC): This is a Creative Commons license.

Attribution (BY): Credit must be given to the creator.

*Disclaimer

This summary highlights only some of the key features and terms of the actual license. It is not a license and has no legal value. Carefully review the actual license before using these materials.

### License Summary*

You are free to share (copy and redistribute) this article in any medium or format and to adapt (remix, transform, and build upon) the material for any purpose, even commercially within the parameters below:

Creative Commons (CC): This is a Creative Commons license.

Attribution (BY): Credit must be given to the creator.

*Disclaimer

This summary highlights only some of the key features and terms of the actual license. It is not a license and has no legal value. Carefully review the actual license before using these materials.

## 1. Introduction

*s*

^{o}to a new state

*s*

^{n}is exactly balanced by the number of transitions from the new to the old state: ρ(

*s*

^{(o)})π(

*s*

^{(o)}→

*s*

^{(n)}) = ρ(

*s*

^{(n)}) π(

*s*

^{(n)}→

*s*

^{(o)}), where ρ(·) is the state space equilibrium probability density and π(·) are the probabilities to make a transition between the two states given the set of possible MC moves. Further, the transition is split into a generation and an acceptance/rejection step such that π(

*s*→

*s*′) =

*P*

_{gen}(

*s*→

*s*′)

*P*

_{acc}(

*s*→

*s*′). In the case that the sampled state space is the configuration space of a molecular system at constant temperature,

*P*

_{gen}might relate to moving a randomly picked particle in a random direction over a small random distance, and ρ(

*s*) is proportional to the Boltzmann weight e

^{–βE(s)}, with β = 1/

*k*

_{B}

*T*the inverse temperature,

*k*

_{B}the Boltzmann constant, and

*E*(

*s*) the state’s energy.

*P*

_{gen}(

*s*→

*s*′) =

*P*

_{gen}(

*s*′ →

*s*) as in the less generic Metropolis scheme. At each MC step, the new state is either accepted or rejected based on the probability above. In case of a rejection, the old state is maintained and resampled. This scheme obeys detailed-balance. In addition, if the set of MC moves are ergodic, equilibrium sampling is guaranteed. When ergodic sampling, even if mathematically obeyed, is slowed down by a rough (free) energy landscape, replica exchange MC becomes useful.

*S*= (

*s*

_{1},

*s*

_{2}, ···,

*s*

_{N}) representing the system in a set of

*N*independent “parallel universes”. The Metropolis scheme applied to the superstate yields

_{i}(·) is the specific probability density of ensemble

*i*. For example, the move that attempts to swap the first two states, implying

*S*

^{o}= (

*s*

_{1},

*s*

_{2}, ···,

*s*

_{N}) and

*S*

^{n}= (

*s*

_{2},

*s*

_{1}, ···,

*s*

_{N}), will be accepted with a probability

## 2. Methods

### 2.1. Finite Swapping

1. | Define | ||||

2. | Assign | ||||

3. | If the worker is finished with its MC move in ensemble | ||||

4. | Take a uniform random number ν between 0 and 1. If ν > | ||||

5. | Among the free ensembles, pick a random pair ( | ||||

6. | Try to swap the states of ensembles | ||||

7. | Select one of the free ensembles at random and assign the worker to that ensemble for performing a new standard move. Go to step 3. |

*S*

^{(o)}→

*S*

^{(n)}, is the same as from new to old,

*S*

^{(n)}→

*S*

^{(o)}. Instead, by writing $S=({s}_{1},\mathcal{E})$, from ensemble 1’s perspective, we have that the number of $({s}_{1}^{(o)},{\mathcal{E}}^{(o)})\to ({s}_{1}^{(n)},{}^{a}\mathcal{E}^{(n)})$ transitions should be equal to the number of $({s}_{1}^{(n)},{\mathcal{E}}^{(o)})\to ({s}_{1}^{(o)},{}^{a}\mathcal{E}^{(n)})$ transitions when the standard move is applied where ${}^{a}\mathcal{E}^{(n)}$ refers to

*any*new environment. In the SI, we show a similar detailed-balance condition for the replica exchange moves. At step 6 we sample only ensemble

*i*and

*j*or, alternatively, all free ensembles get a sample update. This would mean resampling the existing state of those not involved in a swap (“null move”). This makes the approach more similar to the superstate sampling albeit using only free ensembles, as described in the SI. The null move does not reduce the statistical uncertainty, but we mention it here as it makes it easier to explain the infinite swapping approach. But for the detailed-balance conditions to be valid it is imperative that occupied ensembles are not sampled.

*K*is less than the number of ensembles

*N*. The case

*K*=

*N*is valid but would reduce the number of replica exchange moves to zero as only one ensemble is free at the maximum. Reducing the

*K*/

*N*ratio will generally imply a higher acceptance in the replica exchange moves as we can expect a higher number of free ensembles whose distributions have significant overlap. What gives the optimum number of workers is therefore a nontrivial question that we will further explore in Section 4. However, for case

*K*< N we can maximize the effect of the replica exchange moves by taking the

*P*

_{RE}parameter as high as possible. In fact, we can simulate the effect of the limit

*P*

_{RE}→ 1 without having to do an infinite number of replica exchange moves explicitly. This leads to an infinite swapping (13) version of our algorithm.

### 2.2. Infinite Swapping

*P*

_{RE}= 1 – δ, we will loop through steps 4–6 for many iterations (

*n*

_{it}= ∑

_{n = 0}

^{∞}

*n*(1 – δ)

^{n}δ = 1/δ in the limit δ → 0) before getting to step 7. When δ vanishes and

*n*

_{it}becomes infinitely large, we expect that all possible swaps will be executed an infinite number of times. Since the swaps obey detailed-balance between unoccupied ensembles, these will essentially sample the distribution of eq 3 (for the subset

*S** of unoccupied ensembles). Hence, when the loop is exited, each possible permutation σ ∈

*S** has been sampled

*n*

_{it}× ρ(σ)/∑

_{σ}ρ(σ) times. By lumping all the times that the same permutation was sampled and normalizing by division with

*n*

_{it}, we simply sample all of the possible permutations in one go using fractional weights that sum up to 1. This is then the only sampling step, as the single update in step 3 can be skipped due to its negligible 1/

*n*

_{it}weight.

*e*

_{1},

*e*

_{2},

*e*

_{3},

*e*

_{4}) containing four states (

*s*

_{1},

*s*

_{2},

*s*

_{3},

*s*

_{4}). Which state is in which ensemble after this step is irrelevant. We can now define a weight-matrix

*W*

*W*

_{ij}∝ ρ

_{j}(

*s*

_{i}). Essential to our approach is the computation of the permanent of the

*W*-matrix, perm(

*W*), and that of the

*W*{

*ij*}-matrices in which the row

*i*and column

*j*are removed.

*W*) = ∑

_{j=1}

^{4}

*W*

_{1j}perm(

*W*{1

*j*}). As the permanent of the 1 × 1 matrix is obviously equal to the single matrix value, the permanent of arbitrary dimension could in principle be solved recursively using this relation. Based on the permanents of

*W*, we will construct a probability matrix

*P*

*P*

_{ij}is the chance to find state

*s*

_{i}in ensemble

*e*

_{j}. As for each permutation each state is in one ensemble and each ensemble contains one state, the

*P*-matrix is bistochastic: both the columns and the rows sum up to 1. If we consider

*S*

_{ij}

^{*}the set of permutations in which state

*s*

_{i}is in

*e*

_{j}, we can write

*P*

_{ij}= ∑

_{σ∈Sij*}ρ(σ)/∑

_{σ′∈S}*ρ(σ′). We can, however, also use the permanent representation in which

*n*×

*n*matrices, will not work for permanents as only some but not all row and column operations have the same effect to a permanent as to a determinant. One can for instance swap rows and columns without changing the permanent. Multiplying a row by a nonzero scalar multiplies the permanent by the same scalar. Hence, this will not affect the P-matrix based on eq 5. Unlike the determinant, adding or subtracting to a row a scalar multiple of another row, an essential part of the Gaussian elimination method, does change the permanent. This makes the permanent computation of a large matrix excessively more expensive than the computation of a determinant. Yet, recent algorithms based on the Balasubramanian-Bax-Franklin-Glynn (BBFG) formula (17−20) scale as $\mathcal{O}({2}^{n})$. This means that the computation of the full

*P*-matrix scales as $\mathcal{O}({2}^{n}\times {n}^{2})$, which seems still steep but is nevertheless a dramatic improvement compared to factorial scaling.

*N*= 7 in the factorial approach, while we reach

*N*= 12 in the BBFG method using a mid-to-high-end laptop (DELL XPS 15 with an Intel Core i7-8750H). If matrix size of

*N*= 20 is the target, the BBFG method can perform a full

*P*-matrix determination in ∼711 s, while it would take ∼15.3 × 10

^{6}years in the factorial approach. The BBFG method is the fastest completely general solution for the problem of computing a

*P*-matrix from any

*W*-matrix. For several algorithms, the

*W*-matrix has special characteristics that can be exploited to further increase efficiency. For instance, if by shuffling the rows and columns the

*W*-matrix can be made into a block form, where squared blocks at the diagonal have only zeros at their right and upper side, the permanent is equal to the product of the block’s permanents. For instance, if

*W*

_{14}=

*W*

_{24}=

*W*

_{34}= 0, we have two blocks, 3 × 3 and 1 × 1. If

*W*

_{13}=

*W*

_{14}=

*W*

_{23}=

*W*

_{24}= 0, we can identify two blocks of 2 × 2, etc. Identification of blocks can hugely decrease the computation of a large permanent. Another speed-up can be made if all rows in the

*W*-matrix are a sequence of ones followed by all zeros, or can be made into that form after the previously mentioned column and row operations. This makes an order $\mathcal{O}({n}^{2})$ approach possible. We will further discuss this in Section 3.1.

3 | If the worker is finished with its MC move in a specific ensemble, the new state is accepted or rejected (but not yet sampled) according to eq 1. The ensemble is free. | ||||

4 | Determine the | ||||

5 | Pick randomly one of the free ensembles | ||||

6 | Pick one of the available states ( | ||||

7 | The worker is assigned to do a new standard move in ensemble |

## 3. Application: ∞RETIS

_{0}, λ

_{1}, ···, λ

_{M}) with λ

_{k+1}> λ

_{k}, that can be viewed as nonintersecting hypersurfaces in phase space having a fixed value of the reaction coordinate. A MC move that generates a trial path not fulfilling the path ensemble’s criteria is always rejected. RETIS defines different path ensembles based on the direction of the paths and the interface that has to be crossed, but all paths start by crossing λ

_{0}(near the reactant state/state

*A*) and they end by either crossing λ

_{0}again or reaching the last interface λ

_{M}(near the product state/state

*B*). There is one special path ensemble, called [0

^{–}], that explores the left side of λ

_{0}, the reactant well, while all other path ensembles, called [

*k*

^{+}] with

*k*= 0, 1, ···

*M*– 1, start by moving to the right from λ

_{0}reaching at least λ

_{k}.

_{0}in the positive direction reaches λ

_{M}without recrossing λ

_{0}. It provides the rate of the process when multiplied with the flux through λ

_{0}(obtained from the path lengths in [0

^{–}] and [0

^{+}] (22)) and is usually an extremely small number. The chance that any of the sampled paths in the [0

^{+}] path ensemble crosses λ

_{M}is generally negligible, but a decent fraction of those (∼0.1–0.5) will cross λ

_{1}and some even λ

_{2}. Likewise, paths in the [

*k*

^{+}],

*k*> 0, path ensembles have a much higher chance to cross λ

_{k+1}than a [0

^{+}]-path as they already cross λ

_{k}. This leads to the calculation of

*M*local conditional crossing probabilities, the chance to cross λ

_{k+1}given λ

_{k}was crossed for

*k*= 0,1, ···,

*M*– 1, whose product gives an exact expression for the overall crossing probability with an exponentially reduced CPU cost compared to MD.

^{–}] and [0

^{+}] ensembles are performed by exchanging the end and start points of these paths that are then continued by MD at the opposite site of the λ

_{0}interface.

### 3.1. *W*-Matrix in RETIS

*M*+ 1 interfaces, λ

_{0}, λ

_{1}, ···, λ

_{M}, there are also

*N*=

*M*+ 1 ensembles, [0

^{–}], [0

^{+}], [1

^{+}], ···, [(

*M*– 1)

^{+}]. For

*K*workers, the size of the

*W*-matrix is, hence, either (

*N*–

*K*+ 1) × (

*N*–

*K*+ 1) or (

*N*–

*K*) × (

*N*–

*K*) as swappings are executed when 1 of the

*K*workers is free, while the remaining

*K*– 1 workers occupy path ensembles that are locked and do not participate in the swap. The smallest matrix occurs when one worker is occupying both [0

^{–}] and [0

^{+}] during the point exchange move, as described in the SI.

*L*+ 1,

*X*= (

*x*

_{0},

*x*

_{1}, ···,

*x*

_{L}), the plain path probability density ρ(

*X*) is given by the probability of the initial phase point times the dynamical transition probabilities to go from one phase point to the next: ρ(

*X*) = ρ(

*x*

_{0})ϕ(

*x*

_{0}→

*x*

_{1})ϕ(

*x*

_{1}→

*x*

_{2}) ··· ϕ(

*x*

_{L–1}→

*x*

_{L}). Here, the transition probabilities depend on the type of dynamics (deterministic, Langevin, Nosé-Hoover dynamics, etc). The weight of a path within a specific path ensemble ρ

_{j}(

*X*) can be expressed as the plain path density times the indicator function

**1**

_{ej}and possibly an additional weight function

*w*

_{j}(

*X*): ρ

_{j}(

*X*) = ρ(

*X*) ×

**1**

_{ej}(

*X*) ×

*w*

_{j}(

*X*). The indicator function equals 1 if path

*X*belongs to ensemble

*e*

_{j}. Otherwise, it is 0. The additional weight function

*w*

_{j}(

*X*) is part of the high-acceptance protocol that is used in combination with the more recent path generation MC moves such as stone skipping (27) and wire fencing. (28) Using these “high-acceptance weights”, nearly all of the CPU-intensive moves can be accepted as they are tuned to cancel the

*P*

_{gen}-terms in the Metropolis–Hastings scheme, eq 1, and the effect of the nonphysical weights is undone in the analysis by weighting each sampled path with the inverse of

*w*

_{j}(

*X*).

*s*

_{i}=

*X*) is difficult to compute, determining

**1**

_{j}(

*s*

_{i}) and

*w*

_{j}(

*s*

_{i}) is trivial. It is therefore a fortunate coincidence that we can replace

*W*

_{ij}= ρ

_{j}(

*s*

_{i}) with

*P*-matrix does not change if we divide or multiply a row by the same number, as mentioned in Section 2. Except for [0

^{–}], all path ensembles have the same start and end conditions and only differ with respect to the interface crossing condition. A path that crosses interface λ

_{k}automatically crosses all lower interfaces λ

_{l<k}. Reversely, if the path does not cross λ

_{k}, it will not cross any of the higher interfaces λ

_{l>k}. This implies that if the columns of

*W*

_{ij}are ordered such that the first column (

*e*

_{1}) is the first available ensemble from the sequence ([0

^{–}], [0

^{+}], [1

^{+}], ···, [(

*M*– 1)

^{+}]), the second column (

*e*

_{2}) is the second available ensemble, and so on, most rows will end with a series of zeros.

*W*-matrix can be brought into a block form such that the permanent can be computed faster based on smaller matrices. In particular, if [0

^{–}] is part of the free ensembles, it will always form a 1 × 1 block as there is always one and no more than one available path that fits in this ensemble.

*w*

_{j}(

*X*) = 1 and each row in the

*W*-matrix (after separating the [0

^{–}] ensemble if it is part of the free ensembles) is a sequence of ones followed by all zeros. The

*W*-matrix can hence be represented by an array (

*n*

_{1},

*n*

_{2},

*n*

_{3}, ···

*n*

_{n}), where each integer

*n*

_{i}indicates the number of ones in row

*i*. As we show in the SI, the permanent of such a

*W*-matrix is simply the product of (

*n*

_{i}+ 1 –

*i*): perm(

*W*) = ∏

_{i}(

*n*

_{i}+ 1 –

*i*). Further, the

*P*-matrix can be constructed from the following order $\mathcal{O}({n}^{2})$ method.

*W*-matrix such that

*n*

_{1}≤

*n*

_{2}≤ ··· ≤

*n*

_{n}. We then fill in the

*P*-matrix from top to bottom for each row using

*P*-matrices from a large

*W*-matrix, up to several thousands, within a second of CPU-time. The above method applies whenever the rows of the

*W*-matrix can be transformed into sequence of ones followed by all zeros. Besides RETIS without high acceptance, this would apply to other MC methods like subset sampling (29) or umbrella sampling (30) with semi-infinite rectangular windows.

## 4. Results and Discussion

*k*

^{+}] has a 0.1 probability to cross λ

_{k+1}and has an average path length of approximately

*k*/10 seconds (see Section 2). The worker is paused for a number of seconds equal to the path length before it can participate in replica exchange moves to mimic the time it would take to perform all of the necessary MD steps. While this artificial simulation allows us to investigate the potential strength of the method to tackle extremely rare events, it cannot reveal the effect of correlations between accepted paths when fast exploration of the reaction coordinate’s orthogonal directions is crucial. To analyze this effect, we also ran a two-dimensional (2D) membrane permeation system with two slightly asymmetric channels. (31) Finally, to study our algorithm with a more generic

*W*-matrix that needs to be solved via the BBFG formula, we also ran a set of underdamped Langevin simulations of a particle in a double-well potential (32) using the recent wire fencing algorithm with the high-acceptance protocol. (28) All simulation results were performed using five independent runs of 12 h. Errors were based on the standard deviations from these five simulations, except for the MSVS process, where a more reliable statistical error was desired for the comparison with analytical results. Here, block errors were determined on each of the five simulations based on the running average of the overall crossing probability. The block errors were finally combined to obtain the statistical error in the average of the five simulations.

### 4.1. Memoryless Single Variable Stochastic (MSVS) Process

^{–50}with the more accurate estimates for the simulations having a large number of workers. Also, the true value is within one standard deviation of the reported averages for 70% of the data points, as is expected from the standard Gaussian confidence intervals. Figure 1a shows the scaling of the MD time (solid lines) and number of MC moves (dashed lines) of the MSVS simulations (orange) compared to linear scaling (black) and the expected scaling for standard replica exchange (REPEX) in which ensembles are updated in cohort (purple).

MSVS | two-channel system | double well with wire fencing | |||||
---|---|---|---|---|---|---|---|

K | P_{cross}/10^{–50} | K | P_{cross}/ 10^{–5} | perm./10^{–6} | K | P_{cross}/10^{–7} | rate/10^{–7} |

1 | 0.61 ± 0.33 | 1 | 1.52 ± 0.17 | 1.28 ± 0.14 | 1 | 5.91 ± 0.18 | 2.59 ± 0.07 |

5 | 1.47 ± 1.04 | 2 | 1.63 ± 0.24 | 1.37 ± 0.20 | 2 | 5.70 ± 0.13 | 2.51 ± 0.06 |

10 | 0.86 ± 0.51 | 3 | 1.52 ± 0.07 | 1.28 ± 0.06 | 3 | 5.57 ± 0.19 | 2.45 ± 0.08 |

15 | 0.68 ± 0.08 | 4 | 1.42 ± 0.10 | 1.19 ± 0.08 | 4 | 5.20 ± 0.30 | 2.34 ± 0.12 |

20 | 1.02 ± 0.13 | 5 | 1.40 ± 0.12 | 1.18 ± 0.10 | 5 | 5.05 ± 0.41 | 2.23 ± 0.18 |

25 | 1.02 ± 0.17 | 6 | 1.54 ± 0.06 | 1.30 ± 0.05 | 6 | 5.49 ± 0.29 | 2.42 ± 0.13 |

30 | 1.26 ± 0.24 | 7 | 1.48 ± 0.08 | 1.24 ± 0.07 | 7 | 4.99 ± 0.39 | 2.21 ± 0.17 |

35 | 1.05 ± 0.15 | 8 | 1.46 ± 0.08 | 1.23 ± 0.06 | 8 | 4.88 ± 0.43 | 2.15 ± 0.19 |

40 | 1.05 ± 0.14 | 9 | 1.42 ± 0.10 | 1.20 ± 0.08 | |||

45 | 0.93 ± 0.09 | 10 | 1.44 ± 0.08 | 1.21 ± 0.07 | |||

50 | 1.00 ± 0.07 | 11 | 1.41 ± 0.09 | 1.19 ± 0.08 | |||

12 | 1.30 ± 0.15 | 1.09 ± 0.12 | |||||

Literature/Theoretical Result | |||||||

1.23 ± 0.16^{b} | 1.06 ± 0.14^{b} | 2.79 ± 0.70^{d} | |||||

5.84 ± 0.13^{e} | 2.58 ± 0.06^{e} | ||||||

1.00^{a} | 1.61^{c} | 1.37^{c} | 5.83^{c} | 2.58^{c} |

_{i}via: ϵ

^{2}= ∑

_{i}ϵ

_{i}

^{2}. It is henceforth clear that additional computational power should not aim to lower the error in a few path ensembles that were already low compared to other path ensembles. We therefore measure the effectiveness of the additional workers by calculating computational efficiencies. The efficiency of a specific computational method is here defined as the inverse computer time, CPU- or wall-time, to obtain an overall relative error equal to 1: ϵ = 1.

^{–50}. The somewhat fluctuating behavior of these curves should hence be viewed as statistical noise as the confidence interval of these efficiencies depends on the statistical error of this error. Despite that, clear trends can be observed in which the CPU-time efficiency is more or less constant within statistical fluctuations, while the wall-time efficiency shows an upward trend. If we neglect the effect of replica exchange moves on the efficiency, we can relate these numerical results with theoretical ones (22,33) for any possible division of a fixed total CPU-time over the different ensembles. A common sense approach would be to aim for the same error ϵ

_{i}in each ensemble (which implies doing the same number of MC moves per ensemble) or to divide the total CPU-time evenly over the ensembles. These two strategies correspond to the case

*K*= 1 or standard RETIS and

*K*=

*N*or standard TIS, respectively. Ref (33) showed that these two strategies provide the same efficiency, and in the SI, we derive that this leads to a wall-time efficiency as a function of the number of workers (

*K*) equal to

*K*/56250, which is denoted by the continuous purple line in Figure 1d (the hypothetical wall-time efficiency for a parallel simulation that uses the CPU hours equally efficient as TIS and serial RETIS, see eq 50 in the SI). The optimum division, however, would give a slightly better wall-time efficiency equal to

*K*/50,000 which is the continuous black line in this figure. Also shown in Figure 1d are the expected theoretical efficiencies based on the numerical distribution of MC moves in each ensemble. This hybrid numerical/theoretical result is shown by the small purple dots and the dashed purple line. This shows that ∞RETIS, at least for a system in which the path length grows linearly with the ensemble’s rank, naturally provides a division of the computational resources that is even better than TIS (

*K*=

*N*) or RETIS (

*K*= 1). Yet, due to statistical inaccuracies this is only evident for the

*K*= 15 case if we base our analysis on the numerical block errors.

*K*=

*N*, which is essentially equivalent of running independent TIS simulations (i.e., without doing any replica exchange moves). We do not expect this to apply to more complex systems where the replica exchange move is a proven weapon for efficient sampling.

### 4.2. Two-Channel Simulations

_{M–1}, that are in the lower barrier channel. While from the average fraction it still looks like the simulations sampled both channels for any number of workers, 4 out of the 5 simulations in the

*K*=

*N*= 12 case solely visited one of the two channels. This is in agreement with previous TIS results. (31) The

*K*= 11 case already provides a dramatic improvement, but is still expected to be suboptimal due to the relatively low frequency of replica exchange moves compared to

*K*< 11. From this 2D system, it would indicate that having

*K*≈

*N*/2 is a safe bet for optimum efficiency.

### 4.3. Double-Well 1D Barrier Using Wire Fencing

*W*-matrix. The same bottleneck can be seen in Figure 1f where both efficiencies plummet with more than two workers. The reported scaling deficiency is of little significance for actual molecular systems where the creation of a full path takes minutes to hours rather than subseconds.

## 5. Conclusions

## Supporting Information

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpca.2c06004.

Details of the simulations, proof that the replica exchange method with cost unbalanced replicas conserves the equilibrium distribution at the individual ensemble level; $\mathcal{O}({n}^{2})$ algorithm for computing the

*P*-matrix from a*W*-matrix for the case that the*W*-matrix consists of rows having a series with ones, followed by zeros; derivations of the theoretical results on the crossing probabilities, rate constant and permeability via Kramers’ theory that are shown in Table 1; discussion on the computational efficiencies with the derivations for the most optimal efficiencies; some additional simulation results on the relative transition probabilities through the lower and higher barrier channel for the two-channel system; our current code with which the model system calculations were done and the project started for advanced simulations (PDF)

## Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

## Acknowledgments

The authors thank Raffaela Cabriolu and Anders Lervik for their critical comments. They acknowledge funding from the Research Council of Norway through FRINATEK Project No. 275506.

## References

This article references 33 other publications.

**1**Metropolis, N.; Rosenbluth, A.; Rosenbluth, M.; Teller, A.; Teller, E. Equation of state calculations by fast computing machines.*J. Chem. Phys.*1953,*21*, 1087– 1092, DOI: 10.1063/1.1699114Google Scholar1https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaG3sXltlKhsw%253D%253D&md5=5dac586da7183d2d81bf35fc83e8c715Equation-of-state calculations by fast computing machinesMetropolis, Nicholas; Rosenbluth, Arianna W.; Rosenbluth, Marshall N.; Teller, Augusta H.; Teller, EdwardJournal of Chemical Physics (1953), 21 (), 1087-92CODEN: JCPSA6; ISSN:0021-9606.A general method, suitable for fast computing machines, is described for investigating such properties as equations of state for substances consisting of interacting individual mols. The method consists of a modified Monte Carlo integration over configuration space. Results for the 2-dimensional rigid-sphere system were obtained on the Los Alamos MANIAC. These results were compared to the free-vol. equation of state and to a 4-term virial coeff. expansion.**2**Hastings, W. K. Monte-Carlo Sampling methods using Markov chains and their Applications.*Biometrika*1970,*57*, 97– 109, DOI: 10.1093/biomet/57.1.97Google ScholarThere is no corresponding record for this reference.**3**Swendsen, R. H.; Wang, J. S. Replica Monte-Carlo simulation of spin-glasses.*Phys. Rev. Lett.*1986,*57*, 2607– 2609, DOI: 10.1103/PhysRevLett.57.2607Google Scholar3https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BC2sfotFaisQ%253D%253D&md5=6f748450b1063d67006b4e0ff4273fddReplica Monte Carlo simulation of spin glassesSwendsen; WangPhysical review letters (1986), 57 (21), 2607-2609 ISSN:.There is no expanded citation for this reference.**4**Marinari, E.; Parisi, G. Simulated tempering - a new Monte-Carlo scheme.*Europhys. Lett.*1992,*19*, 451– 458, DOI: 10.1209/0295-5075/19/6/002Google Scholar4https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK38XltFyks7w%253D&md5=1490d83eab16a3e141e7525ec20208f7Simulated tempering: a new Monte-Carlo schemeMarinari, E.; Parisi, G.Europhysics Letters (1992), 19 (6), 451-8CODEN: EULEEJ; ISSN:0295-5075.A new global optimization method (Simulated Tempering)is proposed for simulating effectively a system with a rough free-energy landscape (i.e., many coexisting states) at nonzero finite temp. This method is related to simulated annealing; but in the present method, the temp. becomes a dynamic variable, and the system is always kept at equil. Application of the method to the Random-Field Ising Model showed a dramatic improvement over conventional Metropolis and cluster models. The conditions under which the method has optimal performances are discussed.**5**Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding.*Chem. Phys. Lett.*1999,*314*, 141– 151, DOI: 10.1016/S0009-2614(99)01123-9Google Scholar5https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK1MXotVKrsLc%253D&md5=0fec0ff81ca7806c1e1ac29e5f50ce19Replica-exchange molecular dynamics method for protein foldingSugita, Y.; Okamoto, Y.Chemical Physics Letters (1999), 314 (1,2), 141-151CODEN: CHPLBC; ISSN:0009-2614. (Elsevier Science B.V.)We have developed a formulation for mol. dynamics algorithm for the replica-exchange method. The effectiveness of the method for the protein-folding problem is tested with the penta-peptide Met-enkephalin. The method can overcome the multiple-min. problem by exchanging non-interacting replicas of the system at several temps. From only one simulation run, one can obtain probability distributions in canonical ensemble for a wide temp. range using multiple-histogram re-weighting techniques, which allows the calcn. of any thermodn. quantity as a function of temp. in that range.**6**Siepmann, J. I.; Frenkel, D. Configurational bias Monte-Carlo - A new sampling scheme for flexible chains.*Mol. Phys.*1992,*75*, 59– 70, DOI: 10.1080/00268979200100061Google Scholar6https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK38XosVGmsA%253D%253D&md5=38387a1aee004db147b9669742844fb9Configurational bias Monte Carlo: a new sampling scheme for flexible chainsSiepmann, Joern Ilja; Frenkel, DaanMolecular Physics (1992), 75 (1), 59-70CODEN: MOPHAM; ISSN:0026-8976.An approach which allows efficient numerical simulation of systems contg. flexible polymers is proposed. A new type of Monte Carlo move is introduced which makes it possible to carry out large-scale conformational changes of the chain mol. in a single trial move. The scheme is based on the self-avoiding random-walk algorithm of Rosenbluth and Rosenbluth (1955). Results of calcns. of mean-square end-to-end lengths for single chains on a 2-dimensional square lattice are compared with those from other simulations.**7**Vlugt, T. J. H.; Krishna, R.; Smit, B. Molecular simulations of adsorption isotherms for linear and branched alkanes and their mixtures in silicalite.*J. Phys. Chem. B*1999,*103*, 1102– 1118, DOI: 10.1021/jp982736cGoogle Scholar7https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK1MXmsFyktQ%253D%253D&md5=f8a60998b144732ec078a212b458498bMolecular Simulations of Adsorption Isotherms for Linear and Branched Alkanes and Their Mixtures in SilicaliteVlugt, T. J. H.; Krishna, R.; Smit, B.Journal of Physical Chemistry B (1999), 103 (7), 1102-1118CODEN: JPCBFK; ISSN:1089-5647. (American Chemical Society)The configurational-bias Monte Carlo (CBMC) technique has been used for computing the adsorption isotherms for linear and branched 2-methylalkanes on silicalite. The carbon nos. of the alkanes ranged from four to nine. For branched alkanes inflection behavior was obsd. for all carbon nos. studied. The inflection was found to occur at a loading of four mols. per unit cell. Below this loading the branched alkanes are seen to be located predominantly at the intersections of the straight and zigzag channels. To obtain loadings higher than four, the branched alkane must seek residence in the channel interiors which is energetically more demanding and therefore requires disproportionately higher pressures; this leads to the inflection behavior. Linear alkanes with six and more carbon atoms also were found to exhibit inflection behavior. Hexane and heptane show inflection due to commensurate "freezing"; the length of these mols. is commensurate with the length of the zigzag channels. This leads to a higher packing efficiency than for other linear alkanes. Available exptl. data from the literature are used to confirm the accuracy of the predictions of the CBMC simulations. Furthermore, the temp. dependency of the isotherms are also properly modeled. For purposes of fitting the isotherms we found that the dual-site Langmuir model provides an excellent description of the simulated isotherms for linear and branched alkanes. In this model we distinguish between two sites with differing ease of adsorption: site A, representing the intersections between the straight and zigzag channels, and site B, representing the channel interiors. CBMC simulations of isotherms of 50-50 binary mixts. of C5, C6, and C7 hydrocarbon isomers show some remarkable and hitherto unreported features. The loading of the branched isomer in all three binary mixts. reaches a max. when the total mixt. loading corresponds to four mols. per unit cell. Higher loadings are obtained by "squeezing out" of the branched alkane from the silicalite and replacing these with the linear alkane. This "squeezing out" effect is found to be entropic in nature; the linear alkanes have a higher packing efficiency and higher loadings are more easily achieved by replacing the branched alkanes with the linear alkanes. The mixt. isotherms can be predicted quite accurately by applying the appropriate mixt. rules to the dual-site Langmuir model. This model allows the mixt. isotherm to be predicted purely on the basis of the parameters describing the isotherms of the pure components. The sorption selectivity exhibited by silicalite for the linear alkane in preference to the branched alkane in mixts. of C5, C6, and C7 hydrocarbon isomers, provides a potential for the development of a novel sepn. technique based on entropy-driven sorption selectivity.**8**Frenkel, D.; Smit, B.*Understanding Molecular Simulations from Algorithms to Applications*; Academic Press: San Diego, California, U.S.A., 2002.Google ScholarThere is no corresponding record for this reference.**9**Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. Transition path sampling and the calculation of rate constants.*J. Chem. Phys.*1998,*108*, 1964, DOI: 10.1063/1.475562Google Scholar9https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK1cXkvFChsQ%253D%253D&md5=4bb2f7b8316dbb4526be81bd24083814Transition path sampling and the calculation of rate constantsDellago, Christoph; Bolhuis, Peter G.; Csajka, Felix S.; Chandler, DavidJournal of Chemical Physics (1998), 108 (5), 1964-1977CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)We have developed a method to study transition pathways for rare events in complex systems. The method can be used to det. rate consts. for transitions between stable states by turning the calcn. of reactive flux correlation functions into the computation of an isomorphic reversible work. In contrast to previous dynamical approaches, the method relies neither on prior knowledge nor on explicit specification of transition states. Rather, it provides an importance sampling from which transition states can be characterized statistically. A simple model is analyzed to illustrate the methodol.**10**Swendsen, R. H.; Wang, J.-S. Nonuniversal critical dynamics in Monte Carlo simulations.*Phys. Rev. Lett.*1987,*58*, 86– 88, DOI: 10.1103/PhysRevLett.58.86Google Scholar10https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BC2sfotF2kug%253D%253D&md5=b4ee869c3c93309088026cfe7fb888c9Nonuniversal critical dynamics in Monte Carlo simulationsSwendsen; WangPhysical review letters (1987), 58 (2), 86-88 ISSN:.There is no expanded citation for this reference.**11**Peters, E. A. J. F.; de With, G. Rejection-free Monte Carlo sampling for general potentials.*Phys. Rev. E*2012,*85*, 026703 DOI: 10.1103/PhysRevE.85.026703Google Scholar11https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC38XktFOitrg%253D&md5=ba0ac4a010e7296311bd304ed480bcc6Rejection-free Monte Carlo sampling for general potentialsPeters, E. A. J. F.; de With, G.Physical Review E: Statistical, Nonlinear, and Soft Matter Physics (2012), 85 (2-2), 026703/1-026703/5CODEN: PRESCM; ISSN:1539-3755. (American Physical Society)A Monte Carlo method to sample the classical configurational canonical ensemble is introduced. In contrast to the Metropolis algorithm, where trial moves can be rejected, in this approach collisions take place. The implementation is event-driven; i.e., at scheduled times the collisions occur. A unique feature of the new method is that smooth potentials (instead of only step-wise changing ones) can be used. In addn. to an event-driven approach, where all particles move simultaneously, we introduce a straight event-chain implementation. As proof of principle, a system of Lennard-Jones particles is simulated.**12**Michel, M.; Kapfer, S. C.; Krauth, W. Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps.*J. Chem. Phys.*2014,*140*, 054116 DOI: 10.1063/1.4863991Google Scholar12https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2cXjslOisbw%253D&md5=6befc6841f295121adb328b99dfdec66Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal stepsMichel, Manon; Kapfer, Sebastian C.; Krauth, WernerJournal of Chemical Physics (2014), 140 (5), 054116/1-054116/8CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)In this article, we present an event-driven algorithm that generalizes the recent hard-sphere event-chain Monte Carlo method without introducing discretizations in time or in space. A factorization of the Metropolis filter and the concept of infinitesimal Monte Carlo moves are used to design a rejection-free Markov-chain Monte Carlo algorithm for particle systems with arbitrary pairwise interactions. The algorithm breaks detailed balance, but satisfies maximal global balance and performs better than the classic, local Metropolis algorithm in large systems. The new algorithm generates a continuum of samples of the stationary probability d. This allows us to compute the pressure and stress tensor as a byproduct of the simulation without any addnl. computations. (c) 2014 American Institute of Physics.**13**Plattner, N.; Doll, J. D.; Dupuis, P.; Wang, H.; Liu, Y.; Gubernatis, J. E. An infinite swapping approach to the rare-event sampling problem.*J. Chem. Phys.*2011,*135*, 134111, DOI: 10.1063/1.3643325Google Scholar13https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC3MXht1OhtrrN&md5=7914460986203590dd157c918ea97bf2An infinite swapping approach to the rare-event sampling problemPlattner, Nuria; Doll, J. D.; Dupuis, Paul; Wang, Hui; Liu, Yufei; Gubernatis, J. E.Journal of Chemical Physics (2011), 135 (13), 134111/1-134111/14CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)We describe a new approach to the rare-event Monte Carlo sampling problem. This technique utilizes a symmetrization strategy to create probability distributions that are more highly connected and, thus, more easily sampled than their original, potentially sparse counterparts. After discussing the formal outline of the approach and devising techniques for its practical implementation, we illustrate the utility of the technique with a series of numerical applications to Lennard-Jones clusters of varying complexity and rare-event character. (c) 2011 American Institute of Physics.**14**Plattner, N.; Doll, J. D.; Meuwly, M. Overcoming the Rare Event Sampling Problem in Biological Systems with Infinite Swapping.*J. Chem. Theory Comput.*2013,*9*, 4215– 4224, DOI: 10.1021/ct400355gGoogle Scholar14https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC3sXht1WltrrE&md5=339b80b1f094873fab385cc94e1e80a9Overcoming the Rare Event Sampling Problem in Biological Systems with Infinite SwappingPlattner, Nuria; Doll, J. D.; Meuwly, MarkusJournal of Chemical Theory and Computation (2013), 9 (9), 4215-4224CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)Infinite swapping (INS) is a recently developed method to address the rare event sampling problem. For INS, an expanded computational ensemble composed of a no. of replicas at different temps. is used, similar to the widely used parallel tempering (PT) method. While the basic concept of PT is to sample various replicas of the system at different temps. and exchange information between the replicas occasionally, INS uses the symmetrized distribution of configurations in temp. space, which corresponds to the infinite swapping limit of PT. The effect of this symmetrization and the enhanced information exchange between replicas is evaluated for three different biol. systems representing different sampling problems in biol.: (1) blocked alanine dipeptide, which is a small system and therefore optimal to evaluate sampling efficiency quant., (2) Villin headpiece, which is used as a test case for the protein folding process, and (3) neuroglobin, which is used to evaluate the effects of enhanced information exchange between replicas for sampling the substate space of a folded protein. For these three test systems, PINS is compared to PT, and it is found that in all cases the sampling with PINS is substantially more efficient.**15**Yu, T.-Q.; Lu, J.; Abrams, C. F.; Vanden-Eijnden, E. Multiscale implementation of infinite-swap replica exchange molecular dynamics.*Proc. Natl. Acad. Sci. U.S.A*2016,*113*, 11744– 11749, DOI: 10.1073/pnas.1605089113Google Scholar15https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC28Xhs1SisbbJ&md5=366a8c35808828f9dfef3229251419deMultiscale implementation of infinite-swap replica exchange molecular dynamicsYu, Tang-Qing; Lu, Jianfeng; Abrams, Cameron F.; Vanden-Eijnden, EricProceedings of the National Academy of Sciences of the United States of America (2016), 113 (42), 11744-11749CODEN: PNASA6; ISSN:0027-8424. (National Academy of Sciences)Replica exchange mol. dynamics (REMD) is a popular method to accelerate conformational sampling of complex mol. systems. The idea is to run several replicas of the system in parallel at different temps. that are swapped periodically. These swaps are typically attempted every few MD steps and accepted or rejected according to a Metropolis-Hastings criterion. This guarantees that the joint distribution of the composite system of replicas is the normalized sum of the symmetrized product of the canonical distributions of these replicas at the different temps. Here the authors propose a different implementation of REMD in which (1) the swaps obey a continuous-time Markov jump process implemented via Gillespie's stochastic simulation algorithm (SSA), which also samples exactly the aforementioned joint distribution and has the advantage of being rejection free, and (2) this REMD-SSA is combined with the heterogeneous multiscale method to accelerate the rate of the swaps and reach the so-called infinite-swap limit that is known to optimize sampling efficiency. The method is easy to implement and can be trivially parallelized. Here the authors illustrate its accuracy and efficiency on the examples of alanine dipeptide in vacuum and C-terminal β-hairpin of protein G in explicit solvent. In this latter example, the authors' results indicate that the landscape of the protein is a triple funnel with two folded structures and one misfolded structure that are stabilized by H-bonds.**16**Lu, J.; Vanden-Eijnden, E. Methodological and Computational Aspects of Parallel Tempering Methods in the Infinite Swapping Limit.*J. Stat. Phys.*2019,*174*, 715– 733, DOI: 10.1007/s10955-018-2210-yGoogle ScholarThere is no corresponding record for this reference.**17**Balasubramanian, K. Combinatorics and Diagonals of Matrices. Ph.D. thesis, Loyola College, Madras, India, 1980.Google ScholarThere is no corresponding record for this reference.**18**Bax, E. Finite-difference Algorithms for Counting Problems. Ph.D. thesis, California Institute of Technology: Pasadena, United States, 1998.Google ScholarThere is no corresponding record for this reference.**19**Bax, E.; Franklin, J.*A finite-Difference Sieve to Compute the Permanent*; CalTech-CS-TR1996; pp 96– 04.Google ScholarThere is no corresponding record for this reference.**20**Glynn, D. G. The permanent of a square matrix.*Eur. J. Comb.*2010,*31*, 1887– 1891, DOI: 10.1016/j.ejc.2010.01.010Google ScholarThere is no corresponding record for this reference.**21**van Erp, T. S. Reaction rate calculation by parallel path swapping.*Phys. Rev. Lett.*2007,*98*, 268301 DOI: 10.1103/PhysRevLett.98.268301Google Scholar21https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD2sXntF2hsLk%253D&md5=ae152f551a7cdf3287a3bd9c454fdbd7Reaction Rate Calculation by Parallel Path Swappingvan Erp, Titus S.Physical Review Letters (2007), 98 (26), 268301/1-268301/4CODEN: PRLTAO; ISSN:0031-9007. (American Physical Society)The efficiency of path sampling simulations can be improved considerably using the approach of path swapping. For this purpose, we devise a new algorithmic procedure based on the transition interface sampling technique. In the same spirit of parallel tempering, paths between different ensembles are swapped, but the role of temp. is here played by the interface position. We test the method on the denaturation transition of DNA using the Peyrard-Bishop-Dauxois model. We find that the new algorithm gives a redn. of the computational cost by a factor of 20.**22**Cabriolu, R.; Refsnes, K. M. S.; Bolhuis, P. G.; van Erp, T. S. Foundations and latest advances in replica exchange transition interface sampling.*J. Chem. Phys.*2017,*147*, 152722, DOI: 10.1063/1.4989844Google Scholar22https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2sXhs1WhsL3F&md5=b0c93a1fee607319883f352b2ba6f83aFoundations and latest advances in replica exchange transition interface samplingCabriolu, Raffaela; Skjelbred Refsnes, Kristin M.; Bolhuis, Peter G.; van Erp, Titus S.Journal of Chemical Physics (2017), 147 (15), 152722/1-152722/17CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)A review. Nearly 20 years ago, transition path sampling (TPS) emerged as an alternative method to free energy based approaches for the study of rare events such as nucleation, protein folding, chem. reactions, and phase transitions. TPS effectively performs Monte Carlo simulations with relatively short mol. dynamics trajectories, with the advantage of not having to alter the actual potential energy surface nor the underlying phys. dynamics. Although the TPS approach also introduced a methodol. to compute reaction rates, this approach was for a long time considered theor. attractive, providing the exact same results as extensively long mol. dynamics simulations, but still expensive for most relevant applications. With the increase of computer power and improvements in the algorithmic methodol., quant. path sampling is finding applications in more and more areas of research. In particular, the transition interface sampling (TIS) and the replica exchange TIS (RETIS) algorithms have, in turn, improved the efficiency of quant. path sampling significantly, while maintaining the exact nature of the approach. Also, open-source software packages are making these methods, for which implementation is not straightforward, now available for a wider group of users. In addn., a blooming development takes place regarding both applications and algorithmic refinements. Therefore, it is timely to explore the wide panorama of the new developments in this field. This is the aim of this article, which focuses on the most efficient exact path sampling approach, RETIS, as well as its recent applications, extensions, and variations. (c) 2017 American Institute of Physics.**23**Dellago, C.; Bolhuis, P. G.; Chandler, D. Efficient transition path sampling: Application to Lennard-Jones cluster rearrangements.*J. Chem. Phys.*1998,*108*, 9236– 9245, DOI: 10.1063/1.476378Google Scholar23https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK1cXjtFSitb8%253D&md5=0d02dabdf6c5811d0bf11847d8949c7eEfficient transition path sampling: Application to Lennard-Jones cluster rearrangementsDellago, Christoph; Bolhuis, Peter G.; Chandler, DavidJournal of Chemical Physics (1998), 108 (22), 9236-9245CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)We develop an efficient Monte-Carlo algorithm to sample an ensemble of stochastic transition paths between stable states. In our description, paths are represented by chains of states linked by Markovian transition probabilities. Rate consts. and mechanisms characterizing the transition may be detd. from the path ensemble. We have previously devised several algorithms for sampling the path ensemble. For these algorithms, the numerical effort scales with the square of the path length. In the new simulation scheme, the required computation scales linearly with the length of the transition path. This improved efficiency allows the calcn. of rate consts. in complex mol. systems. As an example, we study rearrangement processes in a cluster consisting of seven Lennard-Jones particles in two dimensions. Using a quenching technique we are able to identify the relevant transition mechanisms and to locate the related transition states. We furthermore calc. transition rate consts. for various isomerization processes.**24**van Erp, T. S.; Moroni, D.; Bolhuis, P. G. A Novel Path Sampling Method for the Sampling of Rate Constants.*J. Chem. Phys.*2003,*118*, 7762– 7774, DOI: 10.1063/1.1562614Google Scholar24https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD3sXjtVyitrk%253D&md5=20f1b212a17a30c02a8c7d05c34cd499A novel path sampling method for the calculation of rate constantsvan Erp, Titus S.; Moroni, Daniele; Bolhuis, Peter G.Journal of Chemical Physics (2003), 118 (17), 7762-7774CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)We derive a novel efficient scheme to measure the rate const. of transitions between stable states sepd. by high free energy barriers in a complex environment within the framework of transition path sampling. The method is based on directly and simultaneously measuring the fluxes through many phase space interfaces and increases the efficiency with at least a factor of 2 with respect to existing transition path sampling rate const. algorithms. The new algorithm is illustrated on the isomerization of a diat. mol. immersed in a simple fluid.**25**Swenson, D. W. H.; Prinz, J.-H.; Noe, F.; Chodera, J. D.; Bolhuis, P. G. OpenPathSampling: A Python Framework for Path Sampling Simulations. 2. Building and Customizing Path Ensembles and Sample Schemes.*J. Chem. Theory Comput.*2019,*15*, 837– 856, DOI: 10.1021/acs.jctc.8b00627Google Scholar25https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BB3cvjtFGmug%253D%253D&md5=931c9a476f458250b7f3fdbc39a3ce0cOpenPathSampling: A Python Framework for Path Sampling Simulations. 2. Building and Customizing Path Ensembles and Sample SchemesSwenson David W H; Bolhuis Peter G; Swenson David W H; Prinz Jan-Hendrik; Chodera John D; Prinz Jan-Hendrik; Noe FrankJournal of chemical theory and computation (2019), 15 (2), 837-856 ISSN:.The OpenPathSampling (OPS) package provides an easy-to-use framework to apply transition path sampling methodologies to complex molecular systems with a minimum of effort. Yet, the extensibility of OPS allows for the exploration of new path sampling algorithms by building on a variety of basic operations. In a companion paper [ Swenson et al. J. Chem. Theory Comput. 2018 , 10.1021/acs.jctc.8b00626 ] we introduced the basic concepts and the structure of the OPS package, and how it can be employed to perform standard transition path sampling and (replica exchange) transition interface sampling. In this paper, we elaborate on two theoretical developments that went into the design of OPS. The first development relates to the construction of path ensembles, the what is being sampled. We introduce a novel set-based notation for the path ensemble, which provides an alternative paradigm for constructing path ensembles and allows building arbitrarily complex path ensembles from fundamental ones. The second fundamental development is the structure for the customization of Monte Carlo procedures; how path ensembles are being sampled. We describe in detail the OPS objects that implement this approach to customization, the MoveScheme and the PathMover, and provide tools to create and manipulate these objects. We illustrate both the path ensemble building and sampling scheme customization with several examples. OPS thus facilitates both standard path sampling application in complex systems as well as the development of new path sampling methodology, beyond the default.**26**Riccardi, E.; Lervik, A.; Roet, S.; Aarøen, O.; van Erp, T. S. PyRETIS 2: An improbability drive for rare events.*J. Comput. Chem.*2020,*41*, 370– 377, DOI: 10.1002/jcc.26112Google Scholar26https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC1MXitF2mtL%252FP&md5=0432297a909865f505f1685eba1dd412PyRETIS 2: An improbability drive for rare eventsRiccardi, Enrico; Lervik, Anders; Roet, Sander; Aaroen, Ola; van Erp, Titus S.Journal of Computational Chemistry (2020), 41 (4), 370-377CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)The algorithmic development in the field of path sampling has made tremendous progress in recent years. Although the original transition path sampling method was mostly used as a qual. tool to sample reaction paths, the more recent family of interface-based path sampling methods has paved the way for more quant. rate calcn. studies. Of the exact methods, the replica exchange transition interface sampling (RETIS) method is the most efficient, but rather difficult to implement. This has been the main motivation to develop the open-source Python-based computer library PyRETIS that was released in 2017. PyRETIS is designed to be easily interfaced with any mol. dynamics (MD) package using either classical or ab initio MD. In this study, we report on the principles and the software enhancements that are now included in PyRETIS 2, as well as the recent developments on the user interface, improvements of the efficiency via the implementation of new shooting moves, easier initialization procedures, anal. methods, and supported interfaced software. © 2019 The Authors. Journal of Computational Chem. published by Wiley Periodicals, Inc.**27**Riccardi, E.; Dahlen, O.; van Erp, T. S. Fast decorrelating Monte Carlo moves for efficient path sampling.*J. Phys. Chem. Lett.*2017,*8*, 4456– 4460, DOI: 10.1021/acs.jpclett.7b01617Google Scholar27https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2sXhsVWnurnL&md5=989aed1feea51562657c27d8e788d047Fast Decorrelating Monte Carlo Moves for Efficient Path SamplingRiccardi, Enrico; Dahlen, Oda; van Erp, Titus S.Journal of Physical Chemistry Letters (2017), 8 (18), 4456-4460CODEN: JPCLCD; ISSN:1948-7185. (American Chemical Society)Many relevant processes in chem., physics, and biol. are rare events from a computational perspective as they take place beyond the accessible time scale of mol. dynamics (MD). Examples are chem. reactions, nucleation, and conformational changes of biomols. Path sampling is an approach to break this time scale limit via a Monte Carlo (MC) sampling of MD trajectories. Still, many trajectories are needed for accurately predicting rate consts. To improve the speed of convergence, the authors propose two new MC moves, stone skipping and web throwing. In these moves, trajectories are constructed via a sequence of subpaths obeying superdetailed balance. By a reweighting procedure, almost all paths can be accepted. Whereas the generation of a single trajectory becomes more expensive, the reduced correlation results in a significant speedup. For a study on DNA denaturation, the increase was found to be a factor 12.**28**Zhang, D. T.; Riccardi, E.; van Erp, T. S. Enhanced path sampling using subtrajectory Monte Carlo moves. e-Print archive. https://doi.org/10.48550/arXiv.2210.07026.Google ScholarThere is no corresponding record for this reference.**29**Au, S.-K.; Beck, J. L. Estimation of small failure probabilities in high dimensions by subset simulation.*Probabilistic Eng. Mech.*2001,*16*, 263– 277, DOI: 10.1016/S0266-8920(01)00019-4Google ScholarThere is no corresponding record for this reference.**30**Torrie, G.; Valleau, J. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella samping.*J. Comput. Phys.*1977,*23*, 187, DOI: 10.1016/0021-9991(77)90121-8Google ScholarThere is no corresponding record for this reference.**31**Ghysels, A.; Roet, S.; Davoudi, S.; van Erp, T. S. Exact non-Markovian permeability from rare event simulations.*Phys. Rev. Res.*2021,*3*, 033068 DOI: 10.1103/PhysRevResearch.3.033068Google Scholar31https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3MXit1Oru7nL&md5=4f9289705e98613305b7eaa2168878fcExact non-Markovian permeability from rare event simulationsGhysels, An; Roet, Sander; Davoudi, Samaneh; van Erp, Titus S.Physical Review Research (2021), 3 (3), 033068CODEN: PRRHAI; ISSN:2643-1564. (American Physical Society)Permeation of compds. through membranes is important in biol. and engineering processes, e.g., drug delivery through lipid bilayers, anesthetics, or chem. reactor design. Simulations at the at. scale can provide insight in the diffusive pathways and they give ests. of the membrane permeability based on counting membrane transitions or on the inhomogeneous soly.-diffusivity model described by the Smoluchowski equation. For many permeants, permeation through a membrane is too slow to gather sufficient statistics with conventional mol. dynamics simulations, i.e., permeation is a rare event. Recent attempts to improve the description of the dynamics of such rare permeation events have been based on milestoning, which allows the study of processes at timescales beyond those achievable by straightforward mol. dynamics. The approach is not relying on an overdamped description, but, still, it uses a Markovian approxn. which is only valid for small permeants that are not disruptive to the membrane structure. To overcome this fundamental limitation, we show here how replica exchange transition interface sampling (RETIS) can effectively be used on this problem by deriving an effective set of equations that relate the outcome of RETIS simulations and the permeability coeff. In addn., we introduce two new path Monte Carlo (MC) moves specifically for permeation dynamics, that are used in combination with the ordinary path generating moves, which considerably increase the efficiency. The advantage of our method is that it gives exact results, identical to brute force mol. dynamics, but orders of magnitude faster.**32**Van Erp, T. S. Dynamical Rare Event Simulation Techniques for Equilibrium and Nonequilibrium Systems. In*Advances in Chemical Physics*; John Willey & Sons Inc., 2012; Vol. 151, p 27.Google ScholarThere is no corresponding record for this reference.**33**van Erp, T. S. Efficiency analysis of reaction rate calculation methods using analytical models I: The two-dimensional sharp barrier.*J. Chem. Phys.*2006,*125*, 174106, DOI: 10.1063/1.2363996Google Scholar33https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD28XhtF2htrjN&md5=c87b391368e0fba95be6a9c8e687cb85Efficiency analysis of reaction rate calculation methods using analytical models I: The two-dimensional sharp barriervan Erp, Titus S.Journal of Chemical Physics (2006), 125 (17), 174106/1-174106/20CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)The efficiency is analyzed of different methods for the calcn. of reaction rates in the case of a simple two-dimensional anal. benchmark system. Two classes of methods are considered: the first is based on the free energy calcn. along a reaction coordinate and the calcn. of the transmission coeff., the second on the sampling of dynamical pathways. Scaling rules are given for how this efficiency depends on barrier height and width, and simple optimization rules are handed out for the method-specific parameters. It is shown that the path sampling methods, using the transition interface sampling technique, become exceedingly more efficient than the others when the reaction coordinate is not the optimal one.

## Cited By

This article is cited by 2 publications.

- Wouter Vervust, Daniel T. Zhang, Titus S. van Erp, An Ghysels. Path sampling with memory reduction and replica exchange to reach long permeation timescales. Biophysical Journal
**2023**,*1858*https://doi.org/10.1016/j.bpj.2023.02.021 - Daniel T. Zhang, Enrico Riccardi, Titus S. van Erp. Enhanced path sampling using subtrajectory Monte Carlo moves. The Journal of Chemical Physics
**2023**,*158*(2) , 024113. https://doi.org/10.1063/5.0127249

## References

ARTICLE SECTIONSThis article references 33 other publications.

**1**Metropolis, N.; Rosenbluth, A.; Rosenbluth, M.; Teller, A.; Teller, E. Equation of state calculations by fast computing machines.*J. Chem. Phys.*1953,*21*, 1087– 1092, DOI: 10.1063/1.16991141https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaG3sXltlKhsw%253D%253D&md5=5dac586da7183d2d81bf35fc83e8c715Equation-of-state calculations by fast computing machinesMetropolis, Nicholas; Rosenbluth, Arianna W.; Rosenbluth, Marshall N.; Teller, Augusta H.; Teller, EdwardJournal of Chemical Physics (1953), 21 (), 1087-92CODEN: JCPSA6; ISSN:0021-9606.A general method, suitable for fast computing machines, is described for investigating such properties as equations of state for substances consisting of interacting individual mols. The method consists of a modified Monte Carlo integration over configuration space. Results for the 2-dimensional rigid-sphere system were obtained on the Los Alamos MANIAC. These results were compared to the free-vol. equation of state and to a 4-term virial coeff. expansion.**2**Hastings, W. K. Monte-Carlo Sampling methods using Markov chains and their Applications.*Biometrika*1970,*57*, 97– 109, DOI: 10.1093/biomet/57.1.97There is no corresponding record for this reference.**3**Swendsen, R. H.; Wang, J. S. Replica Monte-Carlo simulation of spin-glasses.*Phys. Rev. Lett.*1986,*57*, 2607– 2609, DOI: 10.1103/PhysRevLett.57.26073https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BC2sfotFaisQ%253D%253D&md5=6f748450b1063d67006b4e0ff4273fddReplica Monte Carlo simulation of spin glassesSwendsen; WangPhysical review letters (1986), 57 (21), 2607-2609 ISSN:.There is no expanded citation for this reference.**4**Marinari, E.; Parisi, G. Simulated tempering - a new Monte-Carlo scheme.*Europhys. Lett.*1992,*19*, 451– 458, DOI: 10.1209/0295-5075/19/6/0024https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK38XltFyks7w%253D&md5=1490d83eab16a3e141e7525ec20208f7Simulated tempering: a new Monte-Carlo schemeMarinari, E.; Parisi, G.Europhysics Letters (1992), 19 (6), 451-8CODEN: EULEEJ; ISSN:0295-5075.A new global optimization method (Simulated Tempering)is proposed for simulating effectively a system with a rough free-energy landscape (i.e., many coexisting states) at nonzero finite temp. This method is related to simulated annealing; but in the present method, the temp. becomes a dynamic variable, and the system is always kept at equil. Application of the method to the Random-Field Ising Model showed a dramatic improvement over conventional Metropolis and cluster models. The conditions under which the method has optimal performances are discussed.**5**Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding.*Chem. Phys. Lett.*1999,*314*, 141– 151, DOI: 10.1016/S0009-2614(99)01123-95https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK1MXotVKrsLc%253D&md5=0fec0ff81ca7806c1e1ac29e5f50ce19Replica-exchange molecular dynamics method for protein foldingSugita, Y.; Okamoto, Y.Chemical Physics Letters (1999), 314 (1,2), 141-151CODEN: CHPLBC; ISSN:0009-2614. (Elsevier Science B.V.)We have developed a formulation for mol. dynamics algorithm for the replica-exchange method. The effectiveness of the method for the protein-folding problem is tested with the penta-peptide Met-enkephalin. The method can overcome the multiple-min. problem by exchanging non-interacting replicas of the system at several temps. From only one simulation run, one can obtain probability distributions in canonical ensemble for a wide temp. range using multiple-histogram re-weighting techniques, which allows the calcn. of any thermodn. quantity as a function of temp. in that range.**6**Siepmann, J. I.; Frenkel, D. Configurational bias Monte-Carlo - A new sampling scheme for flexible chains.*Mol. Phys.*1992,*75*, 59– 70, DOI: 10.1080/002689792001000616https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK38XosVGmsA%253D%253D&md5=38387a1aee004db147b9669742844fb9Configurational bias Monte Carlo: a new sampling scheme for flexible chainsSiepmann, Joern Ilja; Frenkel, DaanMolecular Physics (1992), 75 (1), 59-70CODEN: MOPHAM; ISSN:0026-8976.An approach which allows efficient numerical simulation of systems contg. flexible polymers is proposed. A new type of Monte Carlo move is introduced which makes it possible to carry out large-scale conformational changes of the chain mol. in a single trial move. The scheme is based on the self-avoiding random-walk algorithm of Rosenbluth and Rosenbluth (1955). Results of calcns. of mean-square end-to-end lengths for single chains on a 2-dimensional square lattice are compared with those from other simulations.**7**Vlugt, T. J. H.; Krishna, R.; Smit, B. Molecular simulations of adsorption isotherms for linear and branched alkanes and their mixtures in silicalite.*J. Phys. Chem. B*1999,*103*, 1102– 1118, DOI: 10.1021/jp982736c7https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK1MXmsFyktQ%253D%253D&md5=f8a60998b144732ec078a212b458498bMolecular Simulations of Adsorption Isotherms for Linear and Branched Alkanes and Their Mixtures in SilicaliteVlugt, T. J. H.; Krishna, R.; Smit, B.Journal of Physical Chemistry B (1999), 103 (7), 1102-1118CODEN: JPCBFK; ISSN:1089-5647. (American Chemical Society)The configurational-bias Monte Carlo (CBMC) technique has been used for computing the adsorption isotherms for linear and branched 2-methylalkanes on silicalite. The carbon nos. of the alkanes ranged from four to nine. For branched alkanes inflection behavior was obsd. for all carbon nos. studied. The inflection was found to occur at a loading of four mols. per unit cell. Below this loading the branched alkanes are seen to be located predominantly at the intersections of the straight and zigzag channels. To obtain loadings higher than four, the branched alkane must seek residence in the channel interiors which is energetically more demanding and therefore requires disproportionately higher pressures; this leads to the inflection behavior. Linear alkanes with six and more carbon atoms also were found to exhibit inflection behavior. Hexane and heptane show inflection due to commensurate "freezing"; the length of these mols. is commensurate with the length of the zigzag channels. This leads to a higher packing efficiency than for other linear alkanes. Available exptl. data from the literature are used to confirm the accuracy of the predictions of the CBMC simulations. Furthermore, the temp. dependency of the isotherms are also properly modeled. For purposes of fitting the isotherms we found that the dual-site Langmuir model provides an excellent description of the simulated isotherms for linear and branched alkanes. In this model we distinguish between two sites with differing ease of adsorption: site A, representing the intersections between the straight and zigzag channels, and site B, representing the channel interiors. CBMC simulations of isotherms of 50-50 binary mixts. of C5, C6, and C7 hydrocarbon isomers show some remarkable and hitherto unreported features. The loading of the branched isomer in all three binary mixts. reaches a max. when the total mixt. loading corresponds to four mols. per unit cell. Higher loadings are obtained by "squeezing out" of the branched alkane from the silicalite and replacing these with the linear alkane. This "squeezing out" effect is found to be entropic in nature; the linear alkanes have a higher packing efficiency and higher loadings are more easily achieved by replacing the branched alkanes with the linear alkanes. The mixt. isotherms can be predicted quite accurately by applying the appropriate mixt. rules to the dual-site Langmuir model. This model allows the mixt. isotherm to be predicted purely on the basis of the parameters describing the isotherms of the pure components. The sorption selectivity exhibited by silicalite for the linear alkane in preference to the branched alkane in mixts. of C5, C6, and C7 hydrocarbon isomers, provides a potential for the development of a novel sepn. technique based on entropy-driven sorption selectivity.**8**Frenkel, D.; Smit, B.*Understanding Molecular Simulations from Algorithms to Applications*; Academic Press: San Diego, California, U.S.A., 2002.There is no corresponding record for this reference.**9**Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. Transition path sampling and the calculation of rate constants.*J. Chem. Phys.*1998,*108*, 1964, DOI: 10.1063/1.4755629https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK1cXkvFChsQ%253D%253D&md5=4bb2f7b8316dbb4526be81bd24083814Transition path sampling and the calculation of rate constantsDellago, Christoph; Bolhuis, Peter G.; Csajka, Felix S.; Chandler, DavidJournal of Chemical Physics (1998), 108 (5), 1964-1977CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)We have developed a method to study transition pathways for rare events in complex systems. The method can be used to det. rate consts. for transitions between stable states by turning the calcn. of reactive flux correlation functions into the computation of an isomorphic reversible work. In contrast to previous dynamical approaches, the method relies neither on prior knowledge nor on explicit specification of transition states. Rather, it provides an importance sampling from which transition states can be characterized statistically. A simple model is analyzed to illustrate the methodol.**10**Swendsen, R. H.; Wang, J.-S. Nonuniversal critical dynamics in Monte Carlo simulations.*Phys. Rev. Lett.*1987,*58*, 86– 88, DOI: 10.1103/PhysRevLett.58.8610https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BC2sfotF2kug%253D%253D&md5=b4ee869c3c93309088026cfe7fb888c9Nonuniversal critical dynamics in Monte Carlo simulationsSwendsen; WangPhysical review letters (1987), 58 (2), 86-88 ISSN:.There is no expanded citation for this reference.**11**Peters, E. A. J. F.; de With, G. Rejection-free Monte Carlo sampling for general potentials.*Phys. Rev. E*2012,*85*, 026703 DOI: 10.1103/PhysRevE.85.02670311https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC38XktFOitrg%253D&md5=ba0ac4a010e7296311bd304ed480bcc6Rejection-free Monte Carlo sampling for general potentialsPeters, E. A. J. F.; de With, G.Physical Review E: Statistical, Nonlinear, and Soft Matter Physics (2012), 85 (2-2), 026703/1-026703/5CODEN: PRESCM; ISSN:1539-3755. (American Physical Society)A Monte Carlo method to sample the classical configurational canonical ensemble is introduced. In contrast to the Metropolis algorithm, where trial moves can be rejected, in this approach collisions take place. The implementation is event-driven; i.e., at scheduled times the collisions occur. A unique feature of the new method is that smooth potentials (instead of only step-wise changing ones) can be used. In addn. to an event-driven approach, where all particles move simultaneously, we introduce a straight event-chain implementation. As proof of principle, a system of Lennard-Jones particles is simulated.**12**Michel, M.; Kapfer, S. C.; Krauth, W. Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps.*J. Chem. Phys.*2014,*140*, 054116 DOI: 10.1063/1.486399112https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2cXjslOisbw%253D&md5=6befc6841f295121adb328b99dfdec66Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal stepsMichel, Manon; Kapfer, Sebastian C.; Krauth, WernerJournal of Chemical Physics (2014), 140 (5), 054116/1-054116/8CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)In this article, we present an event-driven algorithm that generalizes the recent hard-sphere event-chain Monte Carlo method without introducing discretizations in time or in space. A factorization of the Metropolis filter and the concept of infinitesimal Monte Carlo moves are used to design a rejection-free Markov-chain Monte Carlo algorithm for particle systems with arbitrary pairwise interactions. The algorithm breaks detailed balance, but satisfies maximal global balance and performs better than the classic, local Metropolis algorithm in large systems. The new algorithm generates a continuum of samples of the stationary probability d. This allows us to compute the pressure and stress tensor as a byproduct of the simulation without any addnl. computations. (c) 2014 American Institute of Physics.**13**Plattner, N.; Doll, J. D.; Dupuis, P.; Wang, H.; Liu, Y.; Gubernatis, J. E. An infinite swapping approach to the rare-event sampling problem.*J. Chem. Phys.*2011,*135*, 134111, DOI: 10.1063/1.364332513https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC3MXht1OhtrrN&md5=7914460986203590dd157c918ea97bf2An infinite swapping approach to the rare-event sampling problemPlattner, Nuria; Doll, J. D.; Dupuis, Paul; Wang, Hui; Liu, Yufei; Gubernatis, J. E.Journal of Chemical Physics (2011), 135 (13), 134111/1-134111/14CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)We describe a new approach to the rare-event Monte Carlo sampling problem. This technique utilizes a symmetrization strategy to create probability distributions that are more highly connected and, thus, more easily sampled than their original, potentially sparse counterparts. After discussing the formal outline of the approach and devising techniques for its practical implementation, we illustrate the utility of the technique with a series of numerical applications to Lennard-Jones clusters of varying complexity and rare-event character. (c) 2011 American Institute of Physics.**14**Plattner, N.; Doll, J. D.; Meuwly, M. Overcoming the Rare Event Sampling Problem in Biological Systems with Infinite Swapping.*J. Chem. Theory Comput.*2013,*9*, 4215– 4224, DOI: 10.1021/ct400355g14https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC3sXht1WltrrE&md5=339b80b1f094873fab385cc94e1e80a9Overcoming the Rare Event Sampling Problem in Biological Systems with Infinite SwappingPlattner, Nuria; Doll, J. D.; Meuwly, MarkusJournal of Chemical Theory and Computation (2013), 9 (9), 4215-4224CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)Infinite swapping (INS) is a recently developed method to address the rare event sampling problem. For INS, an expanded computational ensemble composed of a no. of replicas at different temps. is used, similar to the widely used parallel tempering (PT) method. While the basic concept of PT is to sample various replicas of the system at different temps. and exchange information between the replicas occasionally, INS uses the symmetrized distribution of configurations in temp. space, which corresponds to the infinite swapping limit of PT. The effect of this symmetrization and the enhanced information exchange between replicas is evaluated for three different biol. systems representing different sampling problems in biol.: (1) blocked alanine dipeptide, which is a small system and therefore optimal to evaluate sampling efficiency quant., (2) Villin headpiece, which is used as a test case for the protein folding process, and (3) neuroglobin, which is used to evaluate the effects of enhanced information exchange between replicas for sampling the substate space of a folded protein. For these three test systems, PINS is compared to PT, and it is found that in all cases the sampling with PINS is substantially more efficient.**15**Yu, T.-Q.; Lu, J.; Abrams, C. F.; Vanden-Eijnden, E. Multiscale implementation of infinite-swap replica exchange molecular dynamics.*Proc. Natl. Acad. Sci. U.S.A*2016,*113*, 11744– 11749, DOI: 10.1073/pnas.160508911315https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC28Xhs1SisbbJ&md5=366a8c35808828f9dfef3229251419deMultiscale implementation of infinite-swap replica exchange molecular dynamicsYu, Tang-Qing; Lu, Jianfeng; Abrams, Cameron F.; Vanden-Eijnden, EricProceedings of the National Academy of Sciences of the United States of America (2016), 113 (42), 11744-11749CODEN: PNASA6; ISSN:0027-8424. (National Academy of Sciences)Replica exchange mol. dynamics (REMD) is a popular method to accelerate conformational sampling of complex mol. systems. The idea is to run several replicas of the system in parallel at different temps. that are swapped periodically. These swaps are typically attempted every few MD steps and accepted or rejected according to a Metropolis-Hastings criterion. This guarantees that the joint distribution of the composite system of replicas is the normalized sum of the symmetrized product of the canonical distributions of these replicas at the different temps. Here the authors propose a different implementation of REMD in which (1) the swaps obey a continuous-time Markov jump process implemented via Gillespie's stochastic simulation algorithm (SSA), which also samples exactly the aforementioned joint distribution and has the advantage of being rejection free, and (2) this REMD-SSA is combined with the heterogeneous multiscale method to accelerate the rate of the swaps and reach the so-called infinite-swap limit that is known to optimize sampling efficiency. The method is easy to implement and can be trivially parallelized. Here the authors illustrate its accuracy and efficiency on the examples of alanine dipeptide in vacuum and C-terminal β-hairpin of protein G in explicit solvent. In this latter example, the authors' results indicate that the landscape of the protein is a triple funnel with two folded structures and one misfolded structure that are stabilized by H-bonds.**16**Lu, J.; Vanden-Eijnden, E. Methodological and Computational Aspects of Parallel Tempering Methods in the Infinite Swapping Limit.*J. Stat. Phys.*2019,*174*, 715– 733, DOI: 10.1007/s10955-018-2210-yThere is no corresponding record for this reference.**17**Balasubramanian, K. Combinatorics and Diagonals of Matrices. Ph.D. thesis, Loyola College, Madras, India, 1980.There is no corresponding record for this reference.**18**Bax, E. Finite-difference Algorithms for Counting Problems. Ph.D. thesis, California Institute of Technology: Pasadena, United States, 1998.There is no corresponding record for this reference.**19**Bax, E.; Franklin, J.*A finite-Difference Sieve to Compute the Permanent*; CalTech-CS-TR1996; pp 96– 04.There is no corresponding record for this reference.**20**Glynn, D. G. The permanent of a square matrix.*Eur. J. Comb.*2010,*31*, 1887– 1891, DOI: 10.1016/j.ejc.2010.01.010There is no corresponding record for this reference.**21**van Erp, T. S. Reaction rate calculation by parallel path swapping.*Phys. Rev. Lett.*2007,*98*, 268301 DOI: 10.1103/PhysRevLett.98.26830121https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD2sXntF2hsLk%253D&md5=ae152f551a7cdf3287a3bd9c454fdbd7Reaction Rate Calculation by Parallel Path Swappingvan Erp, Titus S.Physical Review Letters (2007), 98 (26), 268301/1-268301/4CODEN: PRLTAO; ISSN:0031-9007. (American Physical Society)The efficiency of path sampling simulations can be improved considerably using the approach of path swapping. For this purpose, we devise a new algorithmic procedure based on the transition interface sampling technique. In the same spirit of parallel tempering, paths between different ensembles are swapped, but the role of temp. is here played by the interface position. We test the method on the denaturation transition of DNA using the Peyrard-Bishop-Dauxois model. We find that the new algorithm gives a redn. of the computational cost by a factor of 20.**22**Cabriolu, R.; Refsnes, K. M. S.; Bolhuis, P. G.; van Erp, T. S. Foundations and latest advances in replica exchange transition interface sampling.*J. Chem. Phys.*2017,*147*, 152722, DOI: 10.1063/1.498984422https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2sXhs1WhsL3F&md5=b0c93a1fee607319883f352b2ba6f83aFoundations and latest advances in replica exchange transition interface samplingCabriolu, Raffaela; Skjelbred Refsnes, Kristin M.; Bolhuis, Peter G.; van Erp, Titus S.Journal of Chemical Physics (2017), 147 (15), 152722/1-152722/17CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)A review. Nearly 20 years ago, transition path sampling (TPS) emerged as an alternative method to free energy based approaches for the study of rare events such as nucleation, protein folding, chem. reactions, and phase transitions. TPS effectively performs Monte Carlo simulations with relatively short mol. dynamics trajectories, with the advantage of not having to alter the actual potential energy surface nor the underlying phys. dynamics. Although the TPS approach also introduced a methodol. to compute reaction rates, this approach was for a long time considered theor. attractive, providing the exact same results as extensively long mol. dynamics simulations, but still expensive for most relevant applications. With the increase of computer power and improvements in the algorithmic methodol., quant. path sampling is finding applications in more and more areas of research. In particular, the transition interface sampling (TIS) and the replica exchange TIS (RETIS) algorithms have, in turn, improved the efficiency of quant. path sampling significantly, while maintaining the exact nature of the approach. Also, open-source software packages are making these methods, for which implementation is not straightforward, now available for a wider group of users. In addn., a blooming development takes place regarding both applications and algorithmic refinements. Therefore, it is timely to explore the wide panorama of the new developments in this field. This is the aim of this article, which focuses on the most efficient exact path sampling approach, RETIS, as well as its recent applications, extensions, and variations. (c) 2017 American Institute of Physics.**23**Dellago, C.; Bolhuis, P. G.; Chandler, D. Efficient transition path sampling: Application to Lennard-Jones cluster rearrangements.*J. Chem. Phys.*1998,*108*, 9236– 9245, DOI: 10.1063/1.47637823https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADyaK1cXjtFSitb8%253D&md5=0d02dabdf6c5811d0bf11847d8949c7eEfficient transition path sampling: Application to Lennard-Jones cluster rearrangementsDellago, Christoph; Bolhuis, Peter G.; Chandler, DavidJournal of Chemical Physics (1998), 108 (22), 9236-9245CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)We develop an efficient Monte-Carlo algorithm to sample an ensemble of stochastic transition paths between stable states. In our description, paths are represented by chains of states linked by Markovian transition probabilities. Rate consts. and mechanisms characterizing the transition may be detd. from the path ensemble. We have previously devised several algorithms for sampling the path ensemble. For these algorithms, the numerical effort scales with the square of the path length. In the new simulation scheme, the required computation scales linearly with the length of the transition path. This improved efficiency allows the calcn. of rate consts. in complex mol. systems. As an example, we study rearrangement processes in a cluster consisting of seven Lennard-Jones particles in two dimensions. Using a quenching technique we are able to identify the relevant transition mechanisms and to locate the related transition states. We furthermore calc. transition rate consts. for various isomerization processes.**24**van Erp, T. S.; Moroni, D.; Bolhuis, P. G. A Novel Path Sampling Method for the Sampling of Rate Constants.*J. Chem. Phys.*2003,*118*, 7762– 7774, DOI: 10.1063/1.156261424https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD3sXjtVyitrk%253D&md5=20f1b212a17a30c02a8c7d05c34cd499A novel path sampling method for the calculation of rate constantsvan Erp, Titus S.; Moroni, Daniele; Bolhuis, Peter G.Journal of Chemical Physics (2003), 118 (17), 7762-7774CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)We derive a novel efficient scheme to measure the rate const. of transitions between stable states sepd. by high free energy barriers in a complex environment within the framework of transition path sampling. The method is based on directly and simultaneously measuring the fluxes through many phase space interfaces and increases the efficiency with at least a factor of 2 with respect to existing transition path sampling rate const. algorithms. The new algorithm is illustrated on the isomerization of a diat. mol. immersed in a simple fluid.**25**Swenson, D. W. H.; Prinz, J.-H.; Noe, F.; Chodera, J. D.; Bolhuis, P. G. OpenPathSampling: A Python Framework for Path Sampling Simulations. 2. Building and Customizing Path Ensembles and Sample Schemes.*J. Chem. Theory Comput.*2019,*15*, 837– 856, DOI: 10.1021/acs.jctc.8b0062725https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A280%3ADC%252BB3cvjtFGmug%253D%253D&md5=931c9a476f458250b7f3fdbc39a3ce0cOpenPathSampling: A Python Framework for Path Sampling Simulations. 2. Building and Customizing Path Ensembles and Sample SchemesSwenson David W H; Bolhuis Peter G; Swenson David W H; Prinz Jan-Hendrik; Chodera John D; Prinz Jan-Hendrik; Noe FrankJournal of chemical theory and computation (2019), 15 (2), 837-856 ISSN:.The OpenPathSampling (OPS) package provides an easy-to-use framework to apply transition path sampling methodologies to complex molecular systems with a minimum of effort. Yet, the extensibility of OPS allows for the exploration of new path sampling algorithms by building on a variety of basic operations. In a companion paper [ Swenson et al. J. Chem. Theory Comput. 2018 , 10.1021/acs.jctc.8b00626 ] we introduced the basic concepts and the structure of the OPS package, and how it can be employed to perform standard transition path sampling and (replica exchange) transition interface sampling. In this paper, we elaborate on two theoretical developments that went into the design of OPS. The first development relates to the construction of path ensembles, the what is being sampled. We introduce a novel set-based notation for the path ensemble, which provides an alternative paradigm for constructing path ensembles and allows building arbitrarily complex path ensembles from fundamental ones. The second fundamental development is the structure for the customization of Monte Carlo procedures; how path ensembles are being sampled. We describe in detail the OPS objects that implement this approach to customization, the MoveScheme and the PathMover, and provide tools to create and manipulate these objects. We illustrate both the path ensemble building and sampling scheme customization with several examples. OPS thus facilitates both standard path sampling application in complex systems as well as the development of new path sampling methodology, beyond the default.**26**Riccardi, E.; Lervik, A.; Roet, S.; Aarøen, O.; van Erp, T. S. PyRETIS 2: An improbability drive for rare events.*J. Comput. Chem.*2020,*41*, 370– 377, DOI: 10.1002/jcc.2611226https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC1MXitF2mtL%252FP&md5=0432297a909865f505f1685eba1dd412PyRETIS 2: An improbability drive for rare eventsRiccardi, Enrico; Lervik, Anders; Roet, Sander; Aaroen, Ola; van Erp, Titus S.Journal of Computational Chemistry (2020), 41 (4), 370-377CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)The algorithmic development in the field of path sampling has made tremendous progress in recent years. Although the original transition path sampling method was mostly used as a qual. tool to sample reaction paths, the more recent family of interface-based path sampling methods has paved the way for more quant. rate calcn. studies. Of the exact methods, the replica exchange transition interface sampling (RETIS) method is the most efficient, but rather difficult to implement. This has been the main motivation to develop the open-source Python-based computer library PyRETIS that was released in 2017. PyRETIS is designed to be easily interfaced with any mol. dynamics (MD) package using either classical or ab initio MD. In this study, we report on the principles and the software enhancements that are now included in PyRETIS 2, as well as the recent developments on the user interface, improvements of the efficiency via the implementation of new shooting moves, easier initialization procedures, anal. methods, and supported interfaced software. © 2019 The Authors. Journal of Computational Chem. published by Wiley Periodicals, Inc.**27**Riccardi, E.; Dahlen, O.; van Erp, T. S. Fast decorrelating Monte Carlo moves for efficient path sampling.*J. Phys. Chem. Lett.*2017,*8*, 4456– 4460, DOI: 10.1021/acs.jpclett.7b0161727https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BC2sXhsVWnurnL&md5=989aed1feea51562657c27d8e788d047Fast Decorrelating Monte Carlo Moves for Efficient Path SamplingRiccardi, Enrico; Dahlen, Oda; van Erp, Titus S.Journal of Physical Chemistry Letters (2017), 8 (18), 4456-4460CODEN: JPCLCD; ISSN:1948-7185. (American Chemical Society)Many relevant processes in chem., physics, and biol. are rare events from a computational perspective as they take place beyond the accessible time scale of mol. dynamics (MD). Examples are chem. reactions, nucleation, and conformational changes of biomols. Path sampling is an approach to break this time scale limit via a Monte Carlo (MC) sampling of MD trajectories. Still, many trajectories are needed for accurately predicting rate consts. To improve the speed of convergence, the authors propose two new MC moves, stone skipping and web throwing. In these moves, trajectories are constructed via a sequence of subpaths obeying superdetailed balance. By a reweighting procedure, almost all paths can be accepted. Whereas the generation of a single trajectory becomes more expensive, the reduced correlation results in a significant speedup. For a study on DNA denaturation, the increase was found to be a factor 12.**28**Zhang, D. T.; Riccardi, E.; van Erp, T. S. Enhanced path sampling using subtrajectory Monte Carlo moves. e-Print archive. https://doi.org/10.48550/arXiv.2210.07026.There is no corresponding record for this reference.**29**Au, S.-K.; Beck, J. L. Estimation of small failure probabilities in high dimensions by subset simulation.*Probabilistic Eng. Mech.*2001,*16*, 263– 277, DOI: 10.1016/S0266-8920(01)00019-4There is no corresponding record for this reference.**30**Torrie, G.; Valleau, J. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella samping.*J. Comput. Phys.*1977,*23*, 187, DOI: 10.1016/0021-9991(77)90121-8There is no corresponding record for this reference.**31**Ghysels, A.; Roet, S.; Davoudi, S.; van Erp, T. S. Exact non-Markovian permeability from rare event simulations.*Phys. Rev. Res.*2021,*3*, 033068 DOI: 10.1103/PhysRevResearch.3.03306831https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BB3MXit1Oru7nL&md5=4f9289705e98613305b7eaa2168878fcExact non-Markovian permeability from rare event simulationsGhysels, An; Roet, Sander; Davoudi, Samaneh; van Erp, Titus S.Physical Review Research (2021), 3 (3), 033068CODEN: PRRHAI; ISSN:2643-1564. (American Physical Society)Permeation of compds. through membranes is important in biol. and engineering processes, e.g., drug delivery through lipid bilayers, anesthetics, or chem. reactor design. Simulations at the at. scale can provide insight in the diffusive pathways and they give ests. of the membrane permeability based on counting membrane transitions or on the inhomogeneous soly.-diffusivity model described by the Smoluchowski equation. For many permeants, permeation through a membrane is too slow to gather sufficient statistics with conventional mol. dynamics simulations, i.e., permeation is a rare event. Recent attempts to improve the description of the dynamics of such rare permeation events have been based on milestoning, which allows the study of processes at timescales beyond those achievable by straightforward mol. dynamics. The approach is not relying on an overdamped description, but, still, it uses a Markovian approxn. which is only valid for small permeants that are not disruptive to the membrane structure. To overcome this fundamental limitation, we show here how replica exchange transition interface sampling (RETIS) can effectively be used on this problem by deriving an effective set of equations that relate the outcome of RETIS simulations and the permeability coeff. In addn., we introduce two new path Monte Carlo (MC) moves specifically for permeation dynamics, that are used in combination with the ordinary path generating moves, which considerably increase the efficiency. The advantage of our method is that it gives exact results, identical to brute force mol. dynamics, but orders of magnitude faster.**32**Van Erp, T. S. Dynamical Rare Event Simulation Techniques for Equilibrium and Nonequilibrium Systems. In*Advances in Chemical Physics*; John Willey & Sons Inc., 2012; Vol. 151, p 27.There is no corresponding record for this reference.**33**van Erp, T. S. Efficiency analysis of reaction rate calculation methods using analytical models I: The two-dimensional sharp barrier.*J. Chem. Phys.*2006,*125*, 174106, DOI: 10.1063/1.236399633https://chemport.cas.org/services/resolver?origin=ACS&resolution=options&coi=1%3ACAS%3A528%3ADC%252BD28XhtF2htrjN&md5=c87b391368e0fba95be6a9c8e687cb85Efficiency analysis of reaction rate calculation methods using analytical models I: The two-dimensional sharp barriervan Erp, Titus S.Journal of Chemical Physics (2006), 125 (17), 174106/1-174106/20CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)The efficiency is analyzed of different methods for the calcn. of reaction rates in the case of a simple two-dimensional anal. benchmark system. Two classes of methods are considered: the first is based on the free energy calcn. along a reaction coordinate and the calcn. of the transmission coeff., the second on the sampling of dynamical pathways. Scaling rules are given for how this efficiency depends on barrier height and width, and simple optimization rules are handed out for the method-specific parameters. It is shown that the path sampling methods, using the transition interface sampling technique, become exceedingly more efficient than the others when the reaction coordinate is not the optimal one.

## Supporting Information

## Supporting Information

ARTICLE SECTIONSThe Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpca.2c06004.

Details of the simulations, proof that the replica exchange method with cost unbalanced replicas conserves the equilibrium distribution at the individual ensemble level; $O\left({n}^{2}\right)$ algorithm for computing the

*P*-matrix from a*W*-matrix for the case that the*W*-matrix consists of rows having a series with ones, followed by zeros; derivations of the theoretical results on the crossing probabilities, rate constant and permeability via Kramers’ theory that are shown in Table 1; discussion on the computational efficiencies with the derivations for the most optimal efficiencies; some additional simulation results on the relative transition probabilities through the lower and higher barrier channel for the two-channel system; our current code with which the model system calculations were done and the project started for advanced simulations (PDF)

## Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.