Simulation outputs¶

A PyRETIS run leaves behind a directory of plain-text files, and pyretis analyse turns those into a report of tables and figures. This page explains what each file contains and how to read the report – in particular, how to tell a converged simulation from one that has not sampled enough.

What a run writes ¶

A path-sampling run produces one numbered directory per path ensemble, plus a few files at the top level:

retis.toml            the input you supplied
output.toml           the input (defaults resolved, the step size of an
                      OpenMM integrator included) + the restart state
                      and the settings that state was sampled with
output.toml.prev      the run file as it was before its last write
pyretis.log           the run log (see the [output] log keywords)
000/                  the [0^-] ensemble
001/                  the [0^+] ensemble
002/ 003/ ...         the [1^+], [2^+], ... ensembles
report/               written by "pyretis analyse", not by the run

A biretis run numbers its ensembles in another order: 000 and 001 are the [0^-] of A and of B, then come the ensembles of A, those of B and the shared ensembles (see the biretis task). Inside each ensemble directory:

pathensemble.txt      which paths this ensemble HELD, and with what weight
moves.txt             every move this ensemble ATTEMPTED, and how it ended
order.txt             the order parameter along each held path
energy.txt            the energies along each held path
accepted/<n>/         the live paths, the initial paths, the last replaced ones
archive/<n>/          trajectories kept for the long term
rejected/             trials kept by [output] keep_rejected_status
engine_logs/<seq>/    an external engine's logs, kept per [output] engine_log
ensemble.restart      written once, when the initial paths are made

Those first two files answer two different questions, and it is worth being clear about which is which from the start:

File	Answers
`pathensemble.txt`	What was this ensemble holding? One line per path per cycle, with the share of the ensemble that path carried. The crossing probability – and so the rate – is computed from this.
`moves.txt`	What did the sampler try, and did it work? One line per attempted move, accepted or rejected. The acceptance ratios are counted from this, and it is the only record of a rejection.

Neither can be derived from the other. A rejected move never becomes a held path, so it appears only in moves.txt; and under infinite swapping a path can arrive in an ensemble with no move attempted there at all, so it appears only in pathensemble.txt.

The numbering is by slot: 000 is the minus ensemble [0^-], 001 is [0^+], and so on outwards. A single TIS run (a tis task with three interfaces) has no minus ensemble: it writes its one directory under the number [tis] ensemble_number gives it, 002 when the input sets none. A pptis run without zero_left or flux numbers only the windows it builds: 000 is [0^+] with zero_ensemble = true and [1^+] without it, and pyretis analyse labels each directory by that role.

Note

order.txt and energy.txt are written unless the [output] keyword order-file or energy-file is 0 or less; both default to 1. archive/ appears once a replaced path has moved there, which a short run may not reach, rejected/ only with keep_rejected_status, and engine_logs/ only with an external engine. Everything else in the list is always produced.

pathensemble.txt – what each ensemble held ¶

This is the central file: it records the path the ensemble was holding, cycle by cycle, and it is what the crossing-probability analysis reads.

Every line is an occupancy record – which path the ensemble held this cycle, and what share of the ensemble it carried. One path produces one line per cycle it survives, so there are always at least as many lines as cycles. Under infinite swapping there are usually many more: a cycle can spread an ensemble over several live paths at once, each holding a fraction (26839 lines after 20000 cycles in the worked example below).

When a cycle appears more than once, the ensemble was holding more than one path that cycle, and their shares – the FracInc column – sum to one:

#     Step ... Weight            PathNumber  HA-weight  FracInc
        11 ... 2.000000000e+00           14  1.0        0.5
        11 ... 2.000000000e+00            8  1.0        0.5

Attempted moves are not in this file. They are in moves.txt beside it.

#     Step    No.-acc  No.-shoot l m r  Length Acc Mc            Min-O ...
         1          1          0 L M L      68 ACC ld -9.001880000e-01 ...
         3          3          1 L M L     252 ACC sh -9.003640815e-01 ...

The columns are:

Column	Meaning
`Step`	The simulation cycle this row describes.
`No.-acc`	Running count of accepted paths in this ensemble: the initial path, then every accepted trial.
`No.-shoot`	Running count of accepted shooting trials in this ensemble. Both counters advance per attempted move, never per row, so at cycle `N` `No.-shoot` is at most `N` and `No.-acc` at most `N + 1`.
`l m r`	Which interfaces the path touches: whether it starts (`l`), crosses the middle (`m`) and ends (`r`) on the left (`L`) or right (`R`) side. `L M R` is a path that crossed and went on to the next state; `L M L` came back.
`Length`	Number of phase points in the path.
`Acc`	The status: always `ACC` in a file PyRETIS 4 writes, because a row is a path the ensemble held (see the warning below). A PyRETIS 3 file has a row per attempted move and holds the rejection codes as well; see Path types and rejection reasons.
`Mc`	The move that generated this path – `sh` shooting (biased shooting included), `tr` time reversal, `mw` multiresolution wire fencing, `ld` an initial path of the run, and so on (see Move types). PyRETIS 4 writes `ld` for an initial path, kicked or loaded; PyRETIS 3 wrote `ki` for a kicked one. The two swaps are different moves and are labelled differently: `s+`/`s-` is the zero swap between `[0^-]` and `[0^+]` (the RETIS one, or the QuanTIS one when `quantis = true`), and `p+`/`p-` is the PPTIS/REPPTIS swap between adjacent ensembles. Both integrate new trajectory, the PPTIS one a different amount each cycle; what each attempt integrated is the `MD-steps` column of `moves.txt`. A `retis` run attempts no other swap: the positive ensembles exchange their paths through infinite swapping, which has no code, so a path the zero swap made keeps its `s-` in every ensemble it moves to (see the note below). A PyRETIS 3 file writes `s+`/`s-` for every swap: the zero swap, the RETIS swaps between `[i^+]` and `[(i+1)^+]`, which exchanged whole paths and integrated nothing, and the REPPTIS swaps. A move that ran no dynamics at all cannot turn an initialisation path into a sampled one, so its product keeps the source path’s code (`ld`, `ki` or `is`) instead of its own: a time reversal of a loaded path is that same path read backwards, a mirror is it seen from the other side, and a swap or target swap that only copied frames is it rearranged. A restart that changes `[tis] high_accept` stamps `hc` on each path an ensemble holds instead: it was accepted under the acceptance rule the continuation no longer runs. Those four codes are what holds a path out of the production statistics. The test for the first three is whether the move integrated any frame for the product; `hc` clears once the ensemble accepts a path under the new rule.
`Min-O`, `Max-O`	The smallest and largest order-parameter value along the path. `Max-O` is the one that decides whether the path reached the next interface, so it is the number the crossing probability is built from.
`Idx-Min`, `Idx-Max`	The indices of the frames where those extrema occur.
`O-shoot`	The order parameter at the shooting point. Wire fencing shoots each of its jumps from its own point and records the last one. Multiresolution wire fencing (`mwf`) shoots from many points and writes `9.0e+03` here, which the analysis does not read as a shooting point.
`Idx-sh`, `Idx-shN`	For `sh`, the index of the shooting point in the old and in the new path. The other moves use them for their own counts: stone skipping (`ss`) its one-step crossing tries, web throwing (`wt`) its accepted segments, wire fencing (`wf`) its usable segments and multiresolution wire fencing (`mwf`) its successful sets, each followed by the path length.
`Weight`	The statistical weight of the row: `HA-weight` divided by `FracInc`.
`PathNumber`	The run-wide identifier of the path, so a path can be followed across cycles and ensembles.
`HA-weight`	The high-acceptance weight, for the moves that use one.
`FracInc`	The share of the ensemble this row’s path carried this cycle. The shares of the rows of one cycle sum to one, and a row whose `HA-weight` is zero carries its share here.

The last three columns are new in PyRETIS 4; a file PyRETIS 3 wrote ends at Weight (see Coming from PyRETIS 3).

Warning

Every line in this file is an accepted path, because the file says what the ensemble held and a rejected move never becomes a held path. Counting these lines therefore tells you nothing about the acceptance rate – it would report every move as accepted. The rejections are in moves.txt.

This also explains a path that repeats over several cycles: when a move is rejected the ensemble keeps what it had, and that path is recorded again next cycle. Those repeats are not noise – they carry real statistical weight, since a path held for five cycles counts five times. The rejection that caused them is recorded once, in moves.txt.

Two lines sharing a cycle are not a repeat: that is one ensemble holding two paths at once, each with its share.

Note

Mc here is how the held path was originally generated, which for a path that arrived by a swap is the move that created it in its previous ensemble. For what was attempted here, read moves.txt.

Coming from PyRETIS 3 ¶

The first 17 columns are the PyRETIS 3 columns, in the same byte layout, and PathNumber, HA-weight and FracInc follow them. A reader that stops after Weight parses the file. What a row means is different, and that decides what can be computed from it:

	PyRETIS 3	PyRETIS 4
a row is	the move attempted in that cycle, rejections included	a path the ensemble held in that cycle; always `ACC`
rows per cycle	one	one per path held, each with its share in `FracInc`
`Weight`	the high-acceptance weight	`HA-weight / FracInc`
`Mc`	the move attempted in this ensemble	the move that generated the held path, possibly in another ensemble and in an earlier cycle

Infinite swapping lets an ensemble hold several paths at once, which one row per attempted move cannot express. Weight is defined so that the PyRETIS 3 estimator – each row contributes 1 / Weight – builds the occupancy-weighted histogram: PyRETIS 3’s crossing-probability analysis run on these files gives PyRETIS 4’s crossing probabilities once the rows of initialisation paths (ld, ki, is, hc) are left out, as PyRETIS 4 leaves them out. PyRETIS 3 starts an ensemble’s statistics at its first accepted shooting move, and so reports nothing for an ensemble whose accepted paths all arrived by swaps.

Read from moves.txt instead:

acceptance ratios – every row here is accepted, so counting rows reports every move as accepted;
per-move statistics – Mc is how the held path was made, not what this ensemble attempted;
MD cost – MD-steps records what each attempt integrated, the rejected ones included.

Skip cycles by Step, not by rows: a cycle can have several rows. pyretis analyse -skipb / -skipe count cycles.

moves.txt – what the sampler attempted ¶

One line per attempted Monte Carlo move, accepted or rejected. This is the only record of a rejection, and what the acceptance ratios in both pyretis analyse and pyretis status are counted from.

#     Step   Move   Status    Length              Min-O              Max-O            O-shoot   Idx-sh   With   ->Path  MD-steps
         1     sh      ACC        17   -1.095714675e-01   -7.417634393e-02   -8.921102800e-02        1      -        3        16
         2     s-      ACC         4   -1.029370000e-01   -9.127384900e-02                  -        -    000        5         2
         4     tr      ACC         4   -1.029370000e-01   -9.127384900e-02                  -        -      -        7         0
         5     sh      BTL        11   -9.277306100e-02   -6.659124237e-02   -9.277306100e-02        1      -        -        10
        12     tr      BWI       114   -1.106760990e-01    1.024782825e-01                  -        -      -        -         0

Read line by line it is the story of the run: cycle 1 shot and was accepted, producing path 3; cycle 2 swapped with ensemble 000; cycle 5 shot from order parameter -0.0928 and was rejected for being too long backwards (BTL), so no path was produced – but it still integrated 10 frames of MD, which MD-steps records.

Column	Meaning
`Step`	The cycle the move was attempted in.
`Move`	The move that actually ran: `sh` shooting, `bias` biased shooting, `wf` wire fencing, `mwf` multiresolution wire fencing, `ss` stone skipping, `wt` web throwing, `tr` time reversal, `mr` mirror, `ts` target swap, `s+`/`s-` the zero swap between `[0^-]` and `[0^+]`, `p+`/`p-` the PPTIS/REPPTIS swap with the neighbour above/below. `bias` and `mwf` are written in full here; `pathensemble.txt` writes their paths as `sh` and `mw`. The initial paths of the run are not an attempted move and have no line.
`Status`	`ACC`, or the code saying why it was rejected. The codes are listed in Path types and rejection reasons.
`Length`	Length of the trial path, in frames.
`Min-O` / `Max-O`	The lowest and highest order parameter the trial reached. For a rejection these say how close it got.
`O-shoot`	The order parameter of the shooting point. Written for the shooting-like moves only (`sh`, `ss`, `wt`, `wf`, `bias`); `-` for the others. Wire fencing records the point of its last jump.
`Idx-sh`	For `sh` and `bias`, the frame the shooting point was taken from; for `ss` the one-step crossing tries, for `wt` the accepted segments, and for `wf` the usable segments.
`With`	The ensemble this move swapped with; `-` when it was not a swap.
`->Path`	The path the move produced, or `-` when it was rejected – a rejected move produces no path.
`MD-steps`	The frames this attempt integrated in this ensemble: its trial, accepted or rejected, and every sub-trajectory the move integrated and discarded, such as the segments of wire fencing and stone skipping. The starting frame of each propagation is handed over rather than integrated, so it is not counted. A move that runs no dynamics – time reversal, mirror – records `0`. For a swap, each partner’s line records what its own leg integrated. This is the MD the ensemble spent: `pyretis analyse` writes its distribution to `report/NNN_lpath.txt.gz` (third and fourth columns), beside that of the lengths of the paths the ensemble held. The report’s `Avg. path length` and the flux come from those held paths, the `Length` of `pathensemble.txt`.

Important

A column that does not apply to a move is written -, never 0. A zero in O-shoot would read as “the shooting point was at zero”, which for a swap or a time reversal is not imprecise but false. Read - as “this move has no such quantity”.

Because a swap touches two ensembles, it is recorded in both: under the + code (s+ or p+) in the lower-index partner and the - code in the higher one, each naming the other under With. The same event, seen from each side.

Navigating it ¶

The file is plain text with fixed columns, so the ordinary tools work. The questions worth asking, and how:

Question	Command
Acceptance rate per move	`awk '!/^#/{n[$2]++; if($3=="ACC") a[$2]++} END{for(m in n) print m, a[m]"/"n[m]}' moves.txt`
Why are moves rejected?	`awk '!/^#/ && $3!="ACC" {print $3}' moves.txt \| sort \| uniq -c \| sort -rn`
Is the sampler stuck?	`awk '!/^#/{print $3}' moves.txt \| uniq -c` – long runs of one code stand out
Which shooting points failed?	`awk '!/^#/ && $2=="sh" && $3!="ACC" {print $7}' moves.txt`
Are swaps happening?	`awk '!/^#/ && $2 ~ /^[sp][+-]$/ {n++} END{print n+0}' moves.txt`
What made path 144?	`awk '!/^#/ && $10==144' moves.txt`
Everything across all ensembles	`awk '!/^#/ && $3!="ACC" {print FILENAME, $3}' 0*/moves.txt`

The fourth of those is worth dwelling on. The distribution of shooting points that failed is the sharpest diagnostic for a badly placed interface: if the rejections cluster near one end of the path, the shooting is being asked to do something the interface spacing makes unlikely. Simply running pyretis status gives the summary; this gives the detail behind it.

Note

Older runs are read too. Before this split, the attempted moves were written into pathensemble.txt itself with a trailing Trial column (0 occupancy, 1 attempt), and before that into a separate trials.txt. All three layouts still analyse, so a run made with any version can be re-analysed with this one.

order.txt and energy.txt – along the path ¶

Both files are a sequence of blocks, one per pathensemble.txt row, in the same order and labelled with the same cycle:

# Cycle: 1, status: ACC, move: ('ld', nan, 0, 0)
#     Time       Orderp
         0    -0.900188
         1    -0.899927

Time is the frame index along the path, not a physical time. energy.txt has the same block structure with potential, kinetic, total energy and temperature columns; a term the engine does not report is written as nan rather than guessed.

Because the blocks are paired one-to-one with the rows of pathensemble.txt, the two files can be read side by side – that is how the analysis reconstructs the order parameter of a given cycle.

cross.txt – crossings in a flux run ¶

A flux simulation (task = "md-flux") writes three columns: the step, the interface number, and the direction of the crossing (+ or -). The interfaces are numbered from 1, not from 0.

The trajectory stores ¶

Trajectories are kept in two tiers. accepted/<n>/ holds the paths the ensembles are currently working with, the initial paths of the run, and the paths replaced most recently; each later replacement moves the oldest of those to archive/<n>/. The default is archive_every = 1: every replaced path is kept, so a long external-engine run stores every trajectory it ever generated. Raise [output] archive_every to keep only every n-th replaced path when disk space matters.

A path directory holds the frame files of the path and its traj.txt, which names the file, the frame index and the velocity direction of each phase point, next to its order.txt and energy.txt. A frame file keeps the name the engine gave it. Two files of one path can share that name – an engine numbers its files per process, and a PPTIS or target swap keeps frames of an earlier path – and the later one is then stored with _<k> before its extension (21_internal_1.xyz), the name traj.txt gives it.

The frame index of a frame of a LAMMPS dump (.lammpstrj) is its time step, the number below its ITEM: TIMESTEP line: both LAMMPS engines dump frame k of a path at the time step k * subcycles. The frame index of the other formats is the frame number in the file, from 0.

Reading the report ¶

pyretis analyse -i retis.toml writes report/ containing the same report in several formats (reStructuredText, HTML, LaTeX/PDF) plus the figures as PNG and the underlying numbers as .txt.gz. The headline result is at the top:

* The crossing probability:
  P_cross = 0.000000000e+00  ±  nan %
* The initial flux (unit: 1/reduced):
  f_A = 0.417660720  ±  6.431378533 %
* The rate constant (unit: 1/reduced):
  k_AB = 0.000000000e+00  ±  nan %

The rate constant is the product of the flux out of the initial state and the crossing probability, so an unusable value in either one makes the rate unusable (Calculating the rate derives that product and shows where each factor is measured). The example above is taken from a deliberately short run and shows exactly that: see Is it converged?.

The report of a biretis run is the text file report/biretis_report.txt, with the fluxes, the crossing probabilities and the rates of both directions, and $K$ (see the biretis task).

After its top block, the retis report has three tables.

Interfaces simply restates the interface positions the analysis used, so a report is self-contained. Check it first if numbers look unexpected – it is where a mismatched input file shows up.

Crossing probabilities gives, per ensemble, the probability that a path reaching that interface goes on to reach the next one, with its error estimate. These are the numbers whose product is the overall crossing probability. Individual values in a healthy simulation are typically a few tens of percent; a value of exactly 1.000000 with zero error, or exactly 0, means that ensemble sampled nothing informative.

Pathensemble data reports, per ensemble, TIS cycles, four measures of how the sampling is doing, and the average length of the paths the ensemble held, each weighted by its share. TIS cycles is the number of pathensemble.txt rows analysed: one for each path the ensemble held in a cycle, leaving out the rows of the initial paths and those outside -skipb / -skipe. Under infinite swapping it is often more than the number of cycles (26839 for [0^+] after 20000 cycles in the worked example below). Acceptance ratios explains the four measures in full; in short:

Shoot acc. ratio – the fraction of ordinary shooting trials that were accepted.
Swap acc. ratio – the same for swapping moves. In a retis run it reads n/a for the outer ensembles by design, because only the zero swap is an accept/reject move there; the rest exchange through the infinite-swapping weights. A repptis run proposes the PPTIS swap between adjacent ensembles and reports its ratio in every ensemble that attempted one.
HA eff. – the weight a high-acceptance move (wire fencing, stone skipping) collects per attempt. Those moves are built to be accepted, so counting acceptances says nothing; the weight does.
Exchange – how much of the ensemble’s occupancy is handed to another path per cycle. This is what answers, for every ensemble, the question the swap column can only answer for the zero swap.

A ratio shown as n/a means the move was never attempted, not that it was never accepted. A row that reads n/a throughout is an ensemble with no sampled path to analyse; the analysis log says No accepted paths to analyse in ensemble for it.

The figures ¶

Per ensemble, named by its directory number:

Figure	What it shows, and what to look for
`NNN_pcross`	The crossing probability as a function of the order parameter. It should fall smoothly and monotonically from 1 at the ensemble’s own interface. Steps and flat plateaus mean few distinct paths contributed – the histogram is being carried by a handful of trajectories.
`NNN_prun`	The running average of the crossing probability against cycle number. This is the primary convergence diagnostic. It should start noisy and flatten out. If it is still drifting at the right edge, the run has not converged, whatever the quoted error says.
`NNN_perror`	Block-error analysis: the estimated error against block length. Correlated data makes the estimate grow with block length until it reaches a plateau; the plateau value is the honest error. If the curve is still rising at the largest block, the run is too short for its own error estimate to be trusted.
`000_lerror`, `001_lerror`	The block-error analysis of the average path length of `[0^-]` and of `[0^+]`. The flux is computed from these two averages, and its error from theirs: read them like `NNN_perror`.
`NNN_lpath`	The distribution of the lengths of the paths the ensemble held, with their average and standard deviation. A long tail towards the maximum length means trajectories are struggling to terminate, and if paths pile up at `maxlength` they are being truncated, which biases the sampling. `NNN_lpath.txt.gz` holds this distribution and, in its third and fourth columns, that of the MD each attempted move integrated in this ensemble – the `MD-steps` of `moves.txt`, rejected trials included – which is what the sampling of this ensemble cost.
`NNN_shoots`	Where the shooting points were taken and which of them were accepted. Accepted shots concentrated in a narrow region mean the ensemble is exploring poorly.

A pptis or repptis report draws NNN_prun_sl and NNN_perror_sl, and NNN_prun_sr and NNN_perror_sr, in place of NNN_prun and NNN_perror: the same running average and block error, for the paths that start at the left and at the right interface of the ensemble.

And for the simulation as a whole:

Figure	What it shows
`total-probability`	Every ensemble’s crossing probability curve on one axis, before matching.
`matched-probability`	The single curve obtained by matching those together. This is the overall crossing probability, and it should decay smoothly over orders of magnitude. A kink at an interface means neighbouring ensembles disagree where they overlap – usually too few paths, or interfaces spaced too widely.
`overall-rrun`	The running estimate of the rate constant against cycle number: the flux times the running overall crossing probability. Read it exactly like `NNN_prun`: flat means converged. The analysis of a `pptis` or `repptis` run draws the running overall crossing probability instead, as `overall-prun`.
`overall-err`	The block-error analysis of the overall result.

Is it converged?¶

No single number answers this. Work through the following, in order – each is a way for a run to be wrong that the others do not catch:

Did the outermost ensemble ever reach the final state? If the overall crossing probability is exactly 0 and the rate is 0 with a nan error, the answer is no, and nothing else in the report matters yet. This is the normal appearance of a run that is simply too short.
Has the running average flattened? Look at overall-rrun and the per-ensemble NNN_prun. A curve still trending at the right edge means the estimate is still moving; the run needs more cycles.

The two below are the same ensemble of the same system, at 200 and at 20000 cycles. On the left the estimate sits at zero until cycle 150 – nothing had crossed yet – and is then climbing steeply when the run ends: whatever number that run reports is the value the curve happened to have reached. On the right it swings wildly for the first few thousand cycles and then settles onto 0.242 and stays there. That settling is what “converged” looks like.
Has the block-error curve reached a plateau? In overall-err and NNN_perror, a curve still rising at the largest block length means the quoted error is an underestimate – the samples are more correlated than the run is long.
Is the matched probability smooth? Kinks at interface positions in matched-probability mean neighbouring ensembles do not agree in their overlap region. A healthy one falls smoothly through orders of magnitude, with the interfaces (dashed) spread evenly along the decay – here from 1 to about $4 \times 10^{-7}$, the crossing probability the run reports. It flattens past the last interface because nothing beyond it was sampled.
Are the per-ensemble numbers sane? Crossing probabilities of exactly 0 or exactly 1, acceptance ratios near zero, or path-length distributions piled up against maxlength all point at a badly configured ensemble rather than at a merely short run.

If a run fails only test 2 or 3, it needs more cycles. If it fails 4 or 5, more cycles will not help on their own: the interfaces are probably placed too far apart, or the move set does not suit the system.

A worked example: the same simulation, short and long ¶

The clearest way to learn the difference is to read one system’s report twice. Both runs below are the shipped 1D double-well RETIS example, unchanged except for [simulation] steps.

Too short: 200 cycles ¶

* The crossing probability:
  P_cross = 0.000000000e+00  ±  nan %
* The initial flux (unit: 1/reduced):
  f_A = 0.417660720  ±  6.431378533 %
* The rate constant (unit: 1/reduced):
  k_AB = 0.000000000e+00  ±  nan %

This looks alarming and is merely short. Nothing is broken: the flux has a value, because it only needs the [0^-] and [0^+] ensembles. But in 200 cycles no path ever reached the outermost interface, so the crossing probability is exactly zero, and the rate – their product – with it. The nan is the honest report of an error that cannot be formed from zero successes, not a numerical fault.

The tell that this is a short run and not a broken one is that its numbers rest on a handful of moves. The per-ensemble table of the same run has sensible path lengths, growing outward from 314 frames in [0^+] to 1322 in [6^+], but each ensemble attempted only 6 to 16 shooting moves: four ensembles report a crossing probability of exactly 0, and [5^+] accepted one shot of thirteen (0.077). Test 1 fails. Test 5 fails as well, but on so few moves that it cannot yet tell a badly configured ensemble from a short run: the same input run for 20000 cycles, below, passes it.

Long enough: 20000 cycles ¶

* 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 %

Now there is a rate. The flux moved too, from 0.418 ± 6.4 % at 200 cycles to 0.271 ± 2.8 %: a change of more than five times the error the short run quoted, so a short run’s error estimate can be far too small for the flux as well. All the remaining difficulty is in the crossing probability, and it is worth seeing why:

| Ensemble |   Pcross   |   Error    | Rel. error (%) |
|  [0^+]   |  0.244074  |  0.017055  |    6.987693    |
|  [1^+]   |  0.155206  |  0.018709  |   12.054258    |
|  [2^+]   |  0.124862  |  0.020380  |   16.322180    |
|  [3^+]   |  0.124981  |  0.022884  |   18.309615    |
|  [4^+]   |  0.106953  |  0.027190  |   25.422839    |
|  [5^+]   |  0.188040  |  0.031766  |   16.893095    |
|  [6^+]   |  0.032874  |  0.011243  |   34.198995    |

Each ensemble measures a conditional probability of order 0.1, none of them alarming on its own. The overall crossing probability is their product – which is how a rate of $10^{-7}$ is obtained from simulations that never wait $10^{7}$ steps for anything – and the relative errors compound along the chain, from 7 % at the innermost interface to 34 % at the outermost, giving the 54 % on the final result.

So this run has converged in the sense that matters least and most. It has converged enough to produce a rate, which the 200-cycle run could not; it has not converged to a precise one. If you need better than a factor of two, this system needs more cycles, more interfaces in the outer region, or both. Reading the error is not optional.

A second system: water-dimer dissociation with GROMACS ¶

The double well above is a one-particle toy. The same report on a real MD engine – the shipped water-dimer RETIS example, GROMACS at 250 K, 2000 cycles, order parameter the O–O distance – reads:

* The crossing probability:
  P_cross = 0.116483662  ±  33.757375544 %
* The initial flux (unit: 1/gromacs):
  f_A = 0.155524126  ±  24.403500782 %
* The rate constant (unit: 1/gromacs):
  k_AB = 1.811601973e-02  ±  41.654426584 %

The striking difference is the size of the answer: a crossing probability of 0.12 against the double well’s $4 \times 10^{-7}$. Breaking a water dimer is not a rare event on this order parameter – it is barely activated – so the per-ensemble probabilities run the other way:

| Ensemble |   Pcross   |   Error    | Rel. error (%) |
|  [0^+]   |  0.428929  |  0.121031  |   28.216993    |
|  [1^+]   |  0.558257  |  0.073001  |   13.076596    |
|  [2^+]   |  0.722222  |  0.081805  |   11.326792    |
|  [3^+]   |  0.882557  |  0.043049  |    4.877726    |
|  [4^+]   |  0.878774  |  0.030353  |    3.453983    |
|  [5^+]   |  0.929374  |  0.019230  |    2.069189    |
|  [6^+]   |  0.934468  |  0.018839  |    2.016036    |

The water dimer's matched crossing probability, falling by less than two orders of magnitude.

Put that beside the double well’s curve above and the difference is the whole story: this one spans well under two orders of magnitude, and most of the drop happens between the first two interfaces. They rise outward, 0.43 to 0.93, where the double well’s fell. Once the dimer is stretched past the inner interfaces it comes apart almost every time, so the outer ensembles are easy and the whole difficulty sits at the first interface – which is also where the largest relative error is. Read that as a design signal: this interface set is finer than this system needs at the outside and too coarse at the inside.

The acceptance table tells the same story from the sampling side:

| Ensemble | TIS cycles | Shoot acc. ratio | Swap acc. ratio |  Exchange  |
|  [0^-]   |    1998    |     0.617647     |    0.978632     |  0.000000  |
|  [0^+]   |    7074    |     0.689922     |    0.978632     |  0.242415  |
|  [1^+]   |    9798    |     0.600000     |      n/a        |  0.445224  |
|  [2^+]   |   10054    |     0.460630     |      n/a        |  0.597915  |
|  [3^+]   |   10394    |     0.351852     |      n/a        |  0.711665  |
|  [4^+]   |    9870    |     0.274510     |      n/a        |  0.738265  |
|  [5^+]   |    9178    |     0.178862     |      n/a        |  0.758949  |
|  [6^+]   |    8636    |     0.093496     |      n/a        |  0.757857  |

Shooting acceptance falls steeply outward, 0.69 down to 0.09 – an external engine on a floppy molecule is far harder to shoot in than a Langevin particle in a smooth well, where the same column only fell to 0.64. An acceptance below 0.1 in the outermost ensemble is the number to act on: it says most of the work spent there is thrown away, and that the interfaces are too widely spaced at the outside.

Note also that Exchange rises outward here (0.24 to 0.76) while acceptance falls. That is the infinite-swapping machinery doing its job: where shooting struggles, the ensembles lean more on exchanging occupancy with their neighbours.

Note

The engine unit shows in the report (1/gromacs), because the rate carries whatever time unit the engine works in. Compare rates only after converting to a physical unit.

The health check ¶

The per-ensemble table of the long run is what a well-behaved RETIS simulation looks like:

| Ensemble | TIS cycles | Shoot acc. ratio | Swap acc. ratio | Avg. path length |
|  [0^-]   |   19998    |     0.808360     |    1.000000     |    1.5051e+03    |
|  [0^+]   |   26839    |     0.918100     |    1.000000     |    340.941287    |
|  [1^+]   |   29581    |     0.810756     |      n/a        |    507.382453    |
|  [2^+]   |   27511    |     0.718249     |      n/a        |    608.313264    |
|  [3^+]   |   26428    |     0.686131     |      n/a        |    713.262891    |
|  [4^+]   |   26173    |     0.682670     |      n/a        |    819.305071    |
|  [5^+]   |   26499    |     0.671490     |      n/a        |    954.391600    |
|  [6^+]   |   24253    |     0.636364     |      n/a        |    1.1479e+03    |

Two trends are the signature of a healthy run, and both are visible:

Acceptance falls outwards, 0.92 down to 0.64. Outer ensembles ask more of a trajectory, so more trials are rejected. A flat column – and especially a column of exactly 1.000000 – would mean the rejections are not being counted, not that the sampling is perfect.
Paths grow longer outwards, 341 up to 1148 frames, because a path has further to travel before it can return. Lengths that stop growing, or that pile up at maxlength, mean trajectories are being truncated.

The swap column reads n/a from [1^+] outwards by design, as explained above: only the zero swap is an accept/reject move.