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.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, and where they go by report-dir.

The screen output ends with the three numbers this page is about:

* The crossing probability:
  P_cross = 3.908476360e-07  ±  53.824318315 %
* The initial flux (unit: 1/reduced):
  f_A = 0.271443162  ±  2.832354761 %
* The rate constant (unit: 1/reduced):
  k_AB = 1.060929183e-07  ±  53.898789185 %

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 timestep of the lammps.in of a LAMMPS engine that leaves it out. 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.

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, 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. 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 every default resolved, and the [current] state. pyretis analyse -i output.toml reads the input sections from it and gives the analysis of the input; 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 scheduler configuration instead. pyretis analyse -i output.toml reads its settings under the names of the input, through the key table, and a setting it does not hold, the [analysis] section among them, takes its default, which a warning says. The state file of a tis run of one ensemble is refused: the analysis of that ensemble counts a crossing at [tis] detect when the input gives [tis] ensemble_number, and the scheduler configuration holds neither [tis] detect nor whether the input gave the number. Analyse such a run with its input file. The analysis of a tis run of several ensembles counts the crossings of each ensemble at the interface after its middle one, and reads its state file as the state file of any other run.

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.