Introduction to PyRETIS and rare event methods¶
PyRETIS is a computational library for performing molecular simulations of rare events with a focus on transition interface sampling (TIS)  and replica exchange transition interface sampling (RETIS) .
Rare events are called rare because they happen at time or length scales much longer than we are able to simulate through brute-force simulations. Raindrops form every day, but compared to the motion of water molecules, raindrop formation is rare, and simulating it through brute-force simulations would take forever and a day.
Rare event methods aim to sample these rare events inaccessible through brute-force simulation.
Transition Interface Sampling¶
To explain the basic concepts, we consider an illustrative example: The transition between two stable states, from the reactant (labeled A) to the product (labeled B) as illustrated in Fig. 1.
These two states are defined by a progress coordinate (or order parameter), , where state A is for and state B for . The central quantity calculated in TIS/RETIS simulations is the rate constant for the transition which can be expressed as:
and we see that it contains two parts:
- The initial flux which measures how often trajectories start off at the foot of the reaction barrier from the reaction side . In TIS this is obtained from a molecular dynamics (MD) simulation which in PyRETIS is requested by running a md-flux simulation. Note that this is not needed for the RETIS method as explained.
- The crossing probability of reaching before given that has just been crossed. In PyRETIS this can be accomplished by running a tis or retis simulation.
As this is a rare event, the crossing probability is extremely small and nearly impossible to compute in brute-force simulations. Transition interface sampling divides the region between A and B into sub-regions using interfaces, denoted (see Fig. 2). The first interface, is positioned at the interface defining state A (), and the next interface, is put at such that the probability of reaching from is no longer extremely small. We can continue in this fashion and place more interfaces until we reach . What we effectively have done with these interfaces is to split up the computation of the extremely small crossing probability into the computation of many, not-so small, crossing probabilities:
Here, the intermediate crossing probabilities, , are formally defined as the probability of a path crossing given that it originated from , ended in or , and had at least one crossing with in the past.
The interfaces we have placed defines the so-called path ensembles. A path ensemble comprises all possible trajectories that start at the foot of reaction barrier from the reactant side (), end in it or at the product region () and having reached a certain threshold value () between the start point and the final point. This path ensemble is labelled as . The probabilities, can then be obtained as the fraction of paths in the ensemble that also cross .
What is left now, is to have an efficient way of generating trajectories for the various path ensembles. This is in fact done by making use of a selection of Monte Carlo (MC) moves. For TIS  we choose between two moves:
The shooting move which is adapted from the transition path sampling (TPS) shooting algorithm   to allow variable trajectory length. In this move, we generate a new trajectory from an existing trajectory one by:
- Picking randomly one of the discrete MD steps in the present trajectory.
- Modifying the velocities of this phase point (e.g. by randomly drawing new velocities from a Maxwellian distribution).
- Generating a new trajectory from this new phase point by integrating (i.e. running MD simulations) forward and backward in time until A or B is reached. The new trajectory is then obtained by merging the backward and forward trajectories.
The new trajectory is accepted as part of only if all the following criteria are satisfied:
- A detailed balance condition for the energy and path length. 
- It starts at
- It has at least one crossing with before ending in A or B.
The shooting move gives a much higher chance to generate a valid trajectory at each trial compared to simply starting from a random phase point within the reactant well.
The time reversal move which generates trajectories by simply changing the time direction of a path. 
These two moves are illustrated in Fig. 4.
Replica Exchange Transition Interface Sampling¶
The RETIS method is similar to TIS and employs both the shooting and time reversal moves. In addition, RETIS makes use of the swapping move and defines a new ensemble which consist of trajectories that explore the reactant state (see the illustration in Fig. 3).
The swapping move acts between different path simulations. If two simulations generate simultaneously two paths that are valid for each other’s path ensemble, these two paths can be swapped. The swapping moves increase with negligible extra computational cost the number of accepted paths in the ensembles and decrease significantly the correlations between the consecutive paths within the same ensemble. All the moves used for generating trajectories are illustrated in Fig. 4.
There is one notable exception where the swapping move is more computationally demanding: the swap between the and ensembles requires MD simulations. 
In a TIS simulation, as explained above, we have to perform an extra simulation in order to calculate the initial flux. In RETIS, this initial flux can directly be obtained by:
where is the average path length in the ensemble and the average path length in the ensemble.
|||(1, 2, 3, 4) T. S. van Erp, D. Moroni and P. G. Bolhuis, J. Chem. Phys. 118, 7762 (2003) https://dx.doi.org/10.1063%2F1.1562614|
|||(1, 2) T. S. van Erp, Phys. Rev. Lett. 98, 26830 (2007) https://dx.doi.org/10.1103/PhysRevLett.98.268301|
|||C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, J. Chem. Phys. 108, 1964 (1998) https://dx.doi.org/10.1063/1.475562|
|||P. G. Bolhuis, D. Chandler, C. Dellago and P. L. Geissler, Annu. Rev. Phys. Chem. 53, 291 (2002) https://dx.doi.org/10.1146/annurev.physchem.53.082301.113146|