An Always-Accepting Algorithm for Transition Path Sampling

TL;DR

This work introduces two complementary transition path sampling enhancements, ARA-TPS and AAA-TPS, to dramatically improve the efficiency of sampling rare reactive trajectories in overdamped systems. ARA-TPS guarantees that every proposed path is reactive, while AAA-TPS removes acceptance/rejection waste by reweighting trajectories a posteriori, preserving the correct transition path ensemble via $ ilde{P}_{ m AB}[X]= ilde{P}[X]$ weighted by $ ilde{P}_{ m AB}[X]$ and corresponding weights $ ilde{P}[X']/ ilde{P}[X]$. Through two-dimensional model tests and a CO$_2$ clathrate-hydrate nucleation study, the authors demonstrate significant gains in decorrelation and path-space exploration, particularly when using Gaussian shooting-point selection and reweighting. The methods are simple to implement, compatible with existing TPS workflows, and enable adaptive switching between strategies, providing a practical route to accelerate exploration of complex transition pathways in chemical and biomolecular systems. The results indicate substantial improvements in sampling efficiency and channel-switching dynamics, making previously inaccessible pathways tractable under industrially relevant conditions.

Abstract

We present a one-way shooting algorithm for transition path sampling that accepts every proposed trajectory yet samples the correct transition path ensemble for systems with overdamped stochastic dynamics. The method is based on two key elements: a procedure to propose trajectories that are always reactive, and a reweighting scheme that corrects for the bias introduced by always accepting the proposed paths. This approach significantly improves the efficiency of transition path sampling by eliminating the cost associated with generating trajectories that are then rejected. We demonstrate the algorithm by investigating the formation of CO$_2$ clathrate hydrates along different reaction mechanisms, showing that the increased efficiency allows proper sampling of the formation of crystalline hydrates at temperatures and pressures that are difficult to access with conventional algorithms.

Figures (10)

• Figure 1: Schematic representation of the always-reactive path generation algorithm. The old path is divided at $x_s$ and a new trajectory $X'_{\text{new}}$ is initiated from there (dashed lines). If $X'_{\text{new}}$ ends up in $\mathrm{B}$ (case (a)) a new reactive trajectory $X'(L')$ is generated by reindexing $X'_{\text{new}}$ to $X'_{\text{fw}}$ and by adding it to the backward segment ($X_\text{bw}$ in blue). If $X'_{\text{new}}$ ends up in $\mathrm{A}$ (case (b)), a new reactive trajectory $X'(L')$ is generated by reversing and relabeling $X'_{\text{new}}$, now denoted $X'_{\text{bw}}$, and appending it to the forward segment of the reindexed existing path ($X_\text{fw}$ in red).
• Figure 2: Results of the numerical simulations for the two-dimensional models. Standard double well: (A) Potential Energy Surface (PES) and stable state definitions. (B) $p(x|\text{TP})$ obtained via the normalized histogram of the configurations on transition paths. (C) Sum of bin-by-bin absolute difference between $p(x|\text{TP})$ obtained with the methods investigated and the reference transition path ensemble as a function of the number of force evaluations (main) and the number of TPS trials (inset). Bi-stable double well: (D) and (E) are the same as (A) and (B) for this PES. (F) Autocorrelation function of transition channel: to each transition path in the ensemble $X(L)$ a value of $1$ or $-1$ is assigned depending on whether the average $x^1$ on the trajectory is positive or negative; the autocorrelation function of this quantity is plotted as a function of the number of force evaluations (main) and of the number of TPS trials (inset). (G) Average fractional overlap between successive paths as a function of the number of force evaluations (main) and the number TPS trials (inset). For each trajectory in the ensemble, the fractional overlap with successive trajectories is computed until no configurations are shared. The vertical dotted lines represent how many force evaluations (main) or TPS trials (inset) are needed, on average, to obtain a completely new trajectory. Shaded areas represent the 10th and 90th percentiles of the distribution of path overlaps.
• Figure 3: Summary of numerical simulation results for the clathrate system. (A) Schematic representation of CO$_2$ clathrate hydrate structures, showing amorphous (left) and crystalline (right) phases. This representation shows the arrangement of water cages, with CO$_2$ guest molecules omitted to reveal the host architecture. (B) Fractional overlap between a reference path and the successive element of the ensemble as a function of the estimated number of force evaluations in application to a CO$_2$ clathrate hydrate system. Shaded areas represent the 10th and 90th percentiles of the distribution. The vertical dotted lines indicate the average amount of force evaluations needed to fully replace an initial path. (C) Transitions between crystalline (upper level) and amorphous (lower level) channels for all replicas. Each point represents an accepted trial within a replica, with the channel state being determined at the trajectory endpoint. Results are shown for the Always Accepting Algorithm (AAA-TPS) (left), and standard one-way shooting algorithm (right). For the one-way shooting algorithm, rejected trajectories are indicated by a red marker. (D) Histogram of the cage count densities of the $5^{12}6^2$ cage (left illustration) and $4^15^{10}6^2$ cage (right illustration) for the whole transition paths for the Always Accepting Algorithm (AAA-TPS). The red line indicates an equal number of both cage types. Two distinct trajectories are highlighted, with progression in time indicated by a color gradient from dark to light.
• Figure S1: This visualization outlines the three-step procedure for determining the MCG order parameter, which quantifies the size of the largest cluster. (a) Guest Condition: Guest molecules (spheres) within pairwise distance $R_g^{cut}$ = 9 (dashed blue circles) form preliminary edges (black lines). (b) Water Condition: Each edge requires at least $N_w = 5$ water molecules simultaneously located within $R_w^{cut}$ = 6 of both guests (transparent spheres) and inside overlapping conical volumes spanning $45^{\circ}$ radially outwards from the inter-guest axis (dark blue shading). (c) Edges are retained only when sufficient bridging waters exist. Guests with at least $N$ connections ($N=1$ for MCG-1) become MCG monomers. The largest connected cluster (red circles) defines the order parameter.
• Figure S2: (a) Geometric representation of the two most relevant cages in the system, $5^{12}6^2$ and $4^15^{10}6^2$, with nodes representing the water molecules and (b) the cups making up these cages.