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.

In one command ¶

Point the analysis at the same input file the simulation used, from the directory the simulation ran in:

pyretis analyse -i retis.toml

That reads the per-ensemble output in place, prints the headline numbers to the screen and to pyretisanalyse.log, and writes a full report:

report/
    <name>_report_cycles-NNNNNNNNN.rst    the report, reStructuredText
    <name>_report_cycles-NNNNNNNNN.html   the same, for a browser
    <name>_report_cycles-NNNNNNNNN.tex    the same, LaTeX
    <name>_report_cycles-NNNNNNNNN.pdf    the same, if pdflatex is present
    *.png                                 every figure
    *.txt.gz                              the numbers behind each figure

The cycle count is in the file name, so analysing the same run again after more sampling does not overwrite the earlier report. Which formats are written is set by report. The report always goes to report/ in the directory the analysis is run from: report-dir is not read.

Before it lists the report files, the screen output gives the three numbers this page is about:

RETIS Crossing probability: 3.908476360e-07
(Relative error: 53.824318315 %)
Initial flux (units 1/reduced): 0.271443162
(Relative error: 2.832354761 %)
Rate constant (units 1/reduced): 1.060929183e-07
(Relative error: 53.898789185 %)

The report opens with the same three numbers, as \(P_{\text{cross}}\), \(f_A\) and \(k_{AB}\) with their relative errors.

The rest of this page explains where each of those comes from, and what to do when one of them looks wrong.

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.

A run has one \(\Delta t\) for all its engines. The ensembles of an engine pool ([simulation] ensemble_engines) and the [0^-] of a quantis run, with its [engine0], take their MD from different engine sections, and under infinite swapping a path held by one ensemble can have been generated by the engine of another ensemble; a zero swap builds a path from segments of two engines. Path lengths in frames convert to time with one time per step only, so a run whose engine sections have different times per step is refused before it starts. The analysis takes \(\Delta t\) from [engine]: timestep times subcycles, with the time step of the lammps.in of a LAMMPS engine that leaves it out (see The time step of LAMMPS). An OpenMM run records the step size of its integrator as the timestep of [engine] in output.toml, which the analysis of that file reads (see the OpenMM time step). This is the \(\Delta t\) of every engine section of the run, for the flux of the matched analysis and of the WHAM analysis alike. When the time step of [engine] cannot be read, the WHAM analysis reports its rate per step and the matched analysis asks for [engine] timestep. An OpenMM input that states no timestep is analysed through the output.toml of its run, which records it. The output.toml of a run of an earlier PyRETIS records none: the analysis names it, and timestep = <the step size of the integrator, in ps> in [engine] of that file gives it.

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`	With two, or four or more interfaces: one path ensemble per interface, the `[0^-]` included, each sampled independently. With three interfaces: the one ensemble of those interfaces.	With the `[0^-]`: the flux, the crossing probability and the rate, formed as for `retis`. A single TIS run (three interfaces) has no `[0^-]` ensemble and so no flux: it reports a crossing probability and needs the flux from a separate `md-flux` run.
`retis`	The same ensembles, `[0^-]` included, sampled together by infinite swapping: the zero swap between `[0^-]` and `[0^+]` is accepted or rejected, and the other ensembles share the live paths, each path’s occupancy spread over the ensembles it is valid in.	Flux and crossing probability come from the one run, so it reports the rate directly. This is the usual choice. A second estimator (WHAM) is also available – see below.
`biretis`	Two ladders of `retis` ensembles, one for the paths that start in \(A\) and one for the paths that start in \(B\), each with its own `[0^-]`, which meet in shared ensembles that hold the paths of both states; all sampled together by infinite swapping.	Both rates, \(k_{AB} = f_A P_A\) and \(k_{BA} = f_B P_B\), and \(K = k_{AB} / k_{BA}\), from the one run, with the matched estimator only, in `report/biretis_report.txt` (see the biretis task).
`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, rather than by a plain product. From Python, the same MSM also yields mean first-passage times (`mfpt_to_first_last_state()`); the report does not show them.

Note

pyretis run runs every path-sampling task on the same infinite-swapping scheduler (The Runner section (infinite swapping)); what the task selects is the ensemble layout, the move set and the estimator. task = "infinite_swapping" belongs to the deprecated legacy-runner schema: an input of that schema that carries it is converted when it is read (Migration guide: legacy commands and inputs to the unified CLI), and an input of this version that sets it is refused, with a message that names task = "retis" to set instead.

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" without a skip is a cheap cross-check, since both read the same files; with a skip set, the two do not cover the same window (see Choosing the rate-estimation method).

Running the analysis ¶

pyretis analyse -i input.toml

The input is the same file the simulation used. The analysis reads the per-ensemble output in place and writes report/.

Note

The output.toml the run writes records the input of the run, under the names of the input and with its defaults resolved, and the [current] state of the run. pyretis analyse -i output.toml reads the input sections from it and gives the analysis of the input, the step size of an OpenMM integrator among them; the output.toml of a run of a legacy-runner input records its canonical form. The output.toml of a run made with an earlier PyRETIS, a legacy-runner run with an engine pool among them, holds the settings in the form the sampler ran them instead. pyretis analyse -i output.toml reads them under the names of the input, and a setting such a file does not hold, the [analysis] section among them, takes its default, which a warning says. Such a file of a tis run of one ensemble is refused: it does not say at which interface that ensemble counts a crossing ([tis] detect). So is one whose settings cannot all be read back under the names of the input. Analyse those runs with their input file, in the directory of the run. An output.toml of a tis run of several ensembles is read like any other; each of its ensembles counts its crossings at the interface after its middle one.

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, by the matched and the WHAM analysis alike; `-skipb` overrides it.
`skip_final_cycles`	Cycles discarded at the end, by the matched and the WHAM analysis alike; `-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 how many cycles each option dropped, 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.