Interpreting the results¶

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 resolved settings + the state to restart from
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

and inside each ensemble directory:

pathensemble.txt      the paths this ensemble held, cycle by cycle
trials.txt            one row per attempted move, accepted or not
order.txt             the order parameter along each of those paths
energy.txt            the energies along each of those paths
accepted/<n>/         the trajectory files of the live paths
archive/<n>/          trajectories kept for the long term

The numbering is by slot: 000 is the minus ensemble [0^-], 001 is [0^+], and so on outwards. A run without a minus ensemble (single-ensemble TIS) names its one directory after its interface instead.

Note

order.txt and energy.txt are written only when the corresponding [output] keyword asks for them. 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.

A row is written for each path that carried weight in this ensemble during that cycle. Usually that is one row per cycle, but under infinite swapping a cycle can distribute weight over more than one live path, so an ensemble’s file may hold somewhat more rows than the run has cycles – 26839 rows after 20000 cycles in the worked example below. The Step column is the cycle, so rows sharing a Step belong to the same cycle.

#     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.
`No.-shoot`	Running count of shooting moves in this ensemble.
`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: `ACC`, or a rejection code. See Path types and rejection reasons.
`Mc`	The move that generated this path – `sh` shooting, `tr` time reversal, `ld` loaded, `ki` kicked, and so on.
`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.
`Idx-sh`, `Idx-shN`	The index of the shooting point in the old and in the new path.
`Weight`	The statistical weight of the row (see below).
`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.

The last two columns are written by the scheduler and are absent from files produced by much older versions; everything before them is the long-standing format.

Warning

Every row in this file is an accepted path. The file records what each ensemble holds, not what the sampler tried: a rejected trial never becomes the held path, so it has no row here. Counting the rows of this file therefore tells you nothing about the acceptance rate – it would report every move as accepted. Use trials.txt for that.

The same distinction explains a repeated row: when a trial is rejected, the ensemble keeps its previous path, and that path is recorded again for the new cycle. Runs of identical rows are the normal signature of rejections, and they carry real statistical weight – a path held for five cycles counts five times.

trials.txt – what the sampler attempted ¶

One row per attempted move, accepted or rejected:

#     Step   Move   Status   Old path   New path     Length   Min-op   Max-op
         3     s-      ACC          1          9        252 -0.900364 -0.875088
        12     sh      ACC          9         14        252 -0.900364 -0.875088

Move is the move that was attempted and Status how it ended (ACC or a rejection code). New path is the path number the move produced, or -1 when the trial was rejected and the ensemble kept what it had.

This is the only file that records rejections, and it is what the acceptance ratios in the report are computed from. A cycle in which two ensembles were picked together is a swap and appears as s+ in the lower-index partner and s- in the higher 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. When a path is replaced it moves to archive/<n>/, and only every archive_every-th path is kept there – otherwise a long run would store every trajectory it ever generated. Set [output] archive_every = 1 to keep them all.

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. The example above is taken from a deliberately short run and shows exactly that: see Is it converged?.

Three tables follow.

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, the number of analysed cycles, the shooting and swapping acceptance ratios and the average path length. How to read it is described under Acceptance ratios; the short version is that a ratio shown as n/a was never attempted rather than never accepted, and that the swapping column is n/a for the outer ensembles by design, because their exchange is handled by the infinite-swapping weights rather than by an accept/reject step.

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.
`NNN_lpath`	The distribution of path lengths. A long tail towards the maximum length means trajectories are struggling to terminate; if paths pile up at `maxlength` they are being truncated, which biases the sampling.
`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.

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-prun`	The running average of the overall result. Read it exactly like `NNN_prun`: flat means converged.
`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-prun 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.
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.
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 is measured perfectly well, because it only needs the innermost ensemble. 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 everything else looks healthy: the per-ensemble table of the same run has shooting acceptance ratios between 0.5 and 0.9, and sensible path lengths. Test 1 fails; tests 4 and 5 pass.

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. Note what did not change much: the flux, which was already well determined after 200 cycles (2.8 % error here). All the 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.

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.