Calculating the rate¶

Everything a path-sampling run does is in service of one number: the rate constant \(k_{AB}\) at which the system leaves state \(A\) for state \(B\). This page explains how that number is put together – the identity the methods rest on, where each factor comes from in a PyRETIS run, and how the routes differ between TIS, RETIS, infinite swapping and the (RE)PPTIS family.

The identity behind every method ¶

A rare transition is rare because the system almost never gets from \(A\) to \(B\). Waiting for it is hopeless; the trick every method here uses is to factorise the event into pieces that are each common enough to measure.

Put an ordered set of interfaces \(\lambda_0 < \lambda_1 < \dots < \lambda_n\) between the two states, with \(\lambda_0\) on the edge of \(A\) and \(\lambda_n\) on the edge of \(B\). Then

\[k_{AB} = f_A \times P_A(\lambda_n | \lambda_0)\]

where

\(f_A\) is the initial flux: how often, per unit time, the system leaves \(A\) through \(\lambda_0\) at all. This is a frequent event, so plain molecular dynamics measures it.
\(P_A(\lambda_n | \lambda_0)\) is the crossing probability: given that a trajectory has just crossed \(\lambda_0\) on its way out, the probability it reaches \(\lambda_n\) before returning to \(A\). This is the rare part.

The crossing probability is itself factorised, which is the whole point:

\[P_A(\lambda_n | \lambda_0) = \prod_{i=0}^{n-1} P_A(\lambda_{i+1} | \lambda_i)\]

Each factor is the probability of getting from one interface to the next – a conditional probability of order 0.1, which ordinary sampling can measure. Their product can be \(10^{-7}\) or smaller without any single simulation ever having to wait that long. Each ensemble \([i^+]\) samples one factor.

That is why a rate of \(10^{-7}\) comes out of runs of a few thousand cycles, and also why the errors multiply: see the worked example in A worked example: the same simulation, short and long, where seven factors with 7–34 % relative error each give a rate with 54 %.

Where each factor comes from in a run ¶

The initial flux ¶

The flux is measured from the two innermost ensembles, [0^-] and [0^+], which between them cover a full round trip in and out of \(A\). PyRETIS computes it as the reciprocal of the mean time for that round trip:

\[f_A = \frac{1}{(\langle L_{0^-}\rangle - 2 + \langle L_{0^+}\rangle - 2) \, \Delta t}\]

where \(\langle L \rangle\) is the mean accepted path length in frames and \(\Delta t\) the time between frames (the engine time step times subcycles). The two subtracted frames are the crossing points shared with the interface, counted once rather than twice.

This is the cheap half of the calculation: the flux is usually well determined long before the crossing probability is, which is why a short run reports a sensible \(f_A\) next to a useless \(P_A\).

A pure task = "md-flux" run measures only this, and by a different route: it counts effective crossings of each interface during plain MD and divides by time. There the time in state \(A\) is cut into windows of skipcross steps and the flux computed in each, which is what gives the estimate an error bar. The window changes the error, not the flux.

The crossing probability ¶

Each ensemble’s own crossing probability curve is a histogram over the highest order parameter its accepted paths reached. Turning the set of curves into one number is called matching: each ensemble’s curve is scaled so that it continues its neighbour’s where they overlap, and the product of the per-ensemble factors is accumulated as the curves are joined. The result is the matched-probability figure, and its value at the outermost interface is \(P_A(\lambda_n|\lambda_0)\).

The rate ¶

\[k_{AB} = f_A \times P_A, \qquad \frac{\sigma_k}{k} = \sqrt{\left(\frac{\sigma_f}{f}\right)^2 + \left(\frac{\sigma_P}{P}\right)^2}\]

The relative errors add in quadrature, so the rate is never more precise than its worse factor – in practice always the crossing probability.

Note

The rate carries the engine’s time unit: a report from an external engine says 1/gromacs or similar. Convert before comparing.

How the methods differ ¶

All of them sample paths and all of them end at \(k = f_A P_A\). They differ in which ensembles they define and how the per-ensemble probabilities are combined.

`task`	What it samples	How the rate is formed
`md-flux`	Plain MD in \(A\), recording crossings of each interface (`cross.txt`).	Gives \(f_A\) only – no crossing probability, so no rate. Use it to get the flux cheaply, or to place interfaces.
`tis`	One path ensemble per interface, each sampled independently.	The per-ensemble probabilities are matched and multiplied exactly as above. With no `[0^-]` ensemble there is no flux, so a pure TIS run reports a crossing probability and needs the flux from a separate `md-flux` run.
`retis`	The same ensembles plus `[0^-]`, sampled together, with swap moves exchanging paths between neighbours.	Flux and crossing probability come from the one run, so it reports the rate directly. This is the usual choice.
`infinite_swapping`	The same physics as RETIS, but instead of accepting or rejecting a swap the coordinator spreads each path’s occupancy over the ensembles it is valid in.	Identical sampling target, so the same matched estimator applies. A second estimator (WHAM) is also available – see below.
`pptis` / `repptis`	Partial paths: each ensemble samples only between neighbouring interfaces rather than all the way back to \(A\), which keeps path lengths short in long or diffusive systems.	The local probabilities are assembled into a Markov state model over the interfaces and the global crossing probability is obtained from it (`global_pcross_msm()`), rather than by a plain product. The same MSM also yields mean first-passage times.

Note

Whichever task is written in the input, the sampling goes through one code path – the infinite-swapping coordinator. What the task selects is the ensemble layout, the move set and the estimator.

Two estimators over the same data ¶

For a run that produced per-ensemble output, two independent estimators can be applied to it, chosen with the method keyword:

matched (the default; crossing is accepted as its old name): The point-matching product described above – the classic RETIS estimator, with block-error analysis.
wham: Weighted-histogram analysis over the same paths, which reweights all ensembles together instead of matching them pairwise.
both: Run each and report both.

They are different estimators of the same physical rate, so on a converged run they agree within their statistical error. Disagreement beyond that is a signal – usually that the run is not as converged as the error bars suggest. Running method = "both" is a cheap cross-check, since both read the same files.

Running the analysis ¶

pyretis analyse -i input.toml

The input is the same file the simulation used (or the output.toml the run wrote, which carries the resolved settings). The analysis reads the per-ensemble output in place and writes report/.

The arguments are listed in Description of input arguments for pyretis analyse; the two worth knowing are -skipb and -skipe, which choose the cycle window without editing the input file:

pyretis analyse -i input.toml -skipb 500            # drop equilibration
pyretis analyse -i input.toml -skipb 500 -skipe 100  # and a cut-short tail

Everything else is set in the [analysis] section of the input file. The keywords that change the result rather than the presentation are:

Keyword	Effect on the rate
`method`	Which estimator(s) to use: `matched`, `wham` or `both`.
`skip_initial_cycles`	Cycles discarded at the start; `-skipb` overrides it.
`skip_final_cycles`	Cycles discarded at the end; `-skipe` overrides it.
`maxblock`, `blockskip`	The block lengths used for the error estimate. They do not change the rate, but they do change its quoted error – and the plateau of the block-error curve is what makes that error trustworthy.
`ngrid`, `bins`	Resolution of the crossing-probability grid and of the histograms. Too coarse a grid visibly steps the matched curve.

A run analysed twice with different windows gives two legitimate answers over two different sets of cycles; the log states which cycles were used, and the report is stamped with the cycle count.

Is the number trustworthy?¶

The rate is the most-processed quantity the code produces, so it is the last thing to converge and the first to mislead. Before quoting one:

Work through Is it converged? – especially whether the running average has flattened and the block-error curve has reached its plateau.
Check the flux and the crossing probability separately. They fail differently: the flux is usually fine early, so a suspicious rate is nearly always a crossing-probability problem.
Look at the matched-probability curve for kinks at the interfaces. A kink means neighbouring ensembles disagree where they overlap, and the product across that seam is not to be trusted.
If in doubt, run method = "both" and compare the two estimators.

PyRETIS ships a validation suite (examples/validation/) that runs this whole chain on a one-dimensional double well whose rate is known analytically from Kramers’ theory, and checks that the sampled rate lands on it. It is the reference for “does the machinery still give the right answer” – see that directory’s README.rst.