A. Select a configuration in the bin i . At t > 0 use eqn. 2 (with a shift of index) to weight the configurations according to their previous bin. At t = 0 set all configurations in i as equiprobable. B. Run dynamics of duration τ . C. Calculate λ for the evolved configura- tion and find bin j given λ . D. Associate the evolved configuration to bin j for the next timestep, also recording its origin as i . E. Set M t i,j = M t i,j + 1 /n t $_{i}$. ii. Set sampling rate n$_{i}$ t +1 using eqn. 1. (c) Print the matrix M t and the vector ⃗ d $^{t}$. (d) Print any further observables derived using eqn. 4. (f) Set t = t + τ . (e) Set ⃗ d t + τ = M t ⃗ d $^{t}$. (b) Loop for all i s.t. d t i = 0: i. Repeat n t i times: (a) Set all transition probability estimates M t i,j = 0. 5. Loop for t = 0 to t = ∞ : 4. Prepare a vector ⃗ n 0 giving the number of shots to make from each bin at t = 0. 3. Prepare a vector ⃗ d 0 giving the estimated occupa- tion probability of each bin on λ at t = 0. 2. Find λ for each configuration, and associate it to the appropriate discrete bin w.r.t. λ . 1. Prepare a set of (one or more) configurations of the system at time t = 0. The region B is treated as absorbing in this way: at the end of each timestep, the probability vector entry corresponding to B is set equal to zero and the entry corresponding to the region A is incremented by the flux which has been deleted. No new configurations are actu- ally transferred to A , only some of the ‘probability mass’ tracked by ⃗ d $^{t}$. This ‘short circuit’ of the matrix is equiv- alent to a system with an infinite reservoir of states in A and absorbing boundary conditions at B ; which is the premise normally adopted for FFS. In order to achieve a definition of the forward flux which is strictly independent of τ it is necessary that runs which enter the ‘sink’ region, B , are halted immediately. This may be computationally costly if the coordinate λ is costly to calculate, but cannot be avoided if an accurate flux is desired. If paths were allowed to enter and leave B , this would not correspond to the accepted definition of forward reaction flux as the probability per unit time of a first passage from A to B . The set of starting configurations can be prepared in any way that is of interest. E.g. one could prepare an equilibrium set for a given control parameter and then use S-PRES to perform a quench, i.e. change the param- eter and observe the system dynamics. Or one could start out from a single configuration and observe how trajec- tories diverge from this point. If the system is intended to be set up in a stationary state then a conventional rare event sampling method or an initialising round of S- PRES can be run for whatever length of time is needed to prepare a set of configurations with correct associated weights over a good range of λ . S-PRES can be used in two ways, either to calcu- late the time dependent probability distribution of some static observable or to calculate the time dependent re- action flux φ ( t ) between two specifically chosen ‘source’ and ‘sink’ bins on λ (which in the following we call A and B ). The latter quantity is the non-stationary ana- logue of that which is usually calculated using TPS and FFS methods, and the former of that which is usually calculated via IS techniques. Flux calculations typically require special treatment of boundary conditions, in or- der to remove the effects of granularity in time and in order to create a system which can remain far from equi- librium indefinitely. To aid implementation a step-by-step guide to an ex- ample program is provided. I. Boundary Conditions for Flux Calculations H. S-PRES Algorithm Pseudocode The extent to which a given M t can be transplanted to a different timepoint t or to describe a system which was initialised with different starting conditions, or with a different history of control parameters, must be judged by the practitioner. If the series of M t is stable for suc- cessive timepoints then this indicates that the dynamics, at least with respect to λ , may have converged to some stationary limit. Under the assumption of stationarity in the dynamics, diagonalisation of M t to find the infinite- time distribution with respect to λ can usefully be carried out, as well as other tactics commonly used to condense a description of the kinetics from a Markov-like matrix of estimated transition probabilities [17]. over the ‘parent’ bins i can be carried out without loss of information as in the case of Markov dynamics, meaning that our M t − τ i,j represents an average over the configu- rations in i , rather than having the same value for each configuration in i . In particular, the practitioner should be aware that non-trivial correlations of the process on a time-scale longer than τ are not captured by the descrip- tion in terms of the matrices M $^{t}$. 5