Studying permeability with PyRETIS

This example shows how to set up a permeability simulation with PyRETIS. For further details on the derivation of the formulas and description of the Monte Carlo moves, please read and cite the permeability from (RE)TIS paper [1].

This example uses three non-interacting particles on a flat potential. The runnable inputs are bundled under examples/tutorials/path_sampling/internal/permeability/ and consist of retis_perm.toml, flat_potential.py, and initial.xyz.

Verification status: smoke – see Tutorial map.

Tutorial quick start

• Best starting point: examples/tutorials/path_sampling/internal/permeability/.

• Edit first: the simulation section for permeability, zero_left, and interfaces; the orderparameter section for the tracked particle and mirror position.

• Run: pyretis run -i retis_perm.toml -p.

• Analyse: pyretis analyse -i retis_perm.toml after the run.

• Expected output: standard RETIS output plus permeability-specific analysis values in report/.

• Related checks: no dedicated heavy-test fixture yet; use this page as the setup reference and Example test status for current coverage.

New simulation settings

The simulation section has a few extra options:

[simulation]
task = "retis"
steps = 50
interfaces = [
    -0.1,
    0.0,
    0.1,
]
zero_left = -0.2
permeability = true

Here we have:

• zero_left: Tells PyRETIS that the [0^-] ensemble has a left boundary that is not located at -inf.

• permeability: If True, any path in the [0^-] ensemble that starts and ends at one of the interfaces is accepted. If False, any path that hits the zero_left interface will be rejected, leading to incorrect flux calculations. This option also triggers pyretis analyse to calculate \(\xi\), \(\frac{\tau}{dz}\) and the permeability.

The mirror and target-swap moves change the state of the order function: the molecule it tracks and whether its coordinate is mirrored. That state belongs to the path the move produced. Every path records the state it was measured with, each move measures its new frames in the state of the path it continues, and a swap carries the state along with the path into its new ensemble, so all frames of a stored path are measured on one molecule and on one side of the mirror. A path that records no state, such as a loaded initial path, was measured in the configured state: the orderparameter index, unmirrored. The states of the live paths are saved in output.toml under [current] order_state, and a restart gives each path its own state back. No setting is needed for this; the PyRETIS 3 keyword swap_attributes is accepted, has no effect, and logs a warning.

New TIS settings

The new mirror and target-swap moves add a few options to the tis section:

mirror_freq = 0.1
target_freq = 0.1
target_indices = [
    0,
    1,
    2,
]

Here we have:

• mirror_freq: The probability of attempting the mirror move in the [0^-] ensemble.

• target_freq: The probability of attempting the target-swap move in the [0^-] ensemble.

• target_indices: A list of atom indices. The target swap is only attempted between these atoms. Make sure that the original orderparameter.index is included in this list.

New order parameter class

This simulation can be run using the new order parameter class pyretis_permeability. This class is a subclass of the pyretis_position order parameter, but alters the output depending on mirror_pos, relative and offset.

[orderparameter]
class = "pyretis_permeability"
dim = "x"
index = 0
offset = 0
relative = false
mirror_pos = -0.15

Here we have:

• dim: The same as for the class Position.

• index: The index of the particle that will be tracked at the start of the simulation. This attribute will be changed by the target-swap move.

• offset: This order parameter adds an offset to the value of compute_s() before wrapping it into the periodic box. This alters the resulting OP value, but all values will fall within the box vectors. If you want to alter the box vectors instead and do not have access to them, you can use pyretis_permeabilityminusoffset, which subtracts the offset after wrapping, before returning the value.

• relative: If True, the output is mapped relative to the box vector (between 0 and 1). Both offset and mirror_pos should be defined relative to this box vector as well.

• mirror_pos: The position of the mirror plane, before applying the offset. For the current implementation, this must be set halfway between the 0-R and 0-L interfaces.

The permeability classes call the function compute_s() before applying the offset and mirror. For the base class this calls the compute function of Position. To use this with your own OP, you can make a subclass of Permeability (the Python class) and override the self.compute_s() function to return your own custom OP before applying the offset and mirroring.

Output of the new moves

The new moves also lead to some new possible responses in moves.txt, the record of every attempted move. For the mirror move, which is always accepted subject to the constraint on mirror_pos, this is just a new move type called mr.

For the target_swap move:

• A new generated label ts to indicate target-swap.

• A new rejection reason: TSS, which means there are no valid indices to swap to.

• Another new rejection reason TSA, which is a rejection based on the Monte Carlo acceptance.

Another changed response is BTS (backward too short), which is a more common rejection for the [0^-] <-> [0^+] swap. This indicates that the attempted trajectory in [0^-] ended at the L interface, so we do not attempt to extend it into the [0^+] ensemble.

New analysis options

Adding permeability = true to the simulation settings triggers pyretis analyse to also calculate and plot \(\xi\), \(\frac{\tau}{dz}\) and the permeability. This follows the formulas described in the paper [1]. For the calculation of \(\tau\) to work, a reference region has to be chosen. This is done by adding tau_ref_bin to the Analysis section of the TOML file.

[analysis]
tau_ref_bin = [
    -0.175,
    -0.125,
]

This value can be changed in retis_perm.toml and the analysis rerun with another reference region. The analysis also plots a 10-bin histogram of the [0^-] region to help the user select a flat histogram region in this space.

Tested by

The reference run examples/tests/test-internal/permeability/ runs this set-up and compares its output files against a committed reference. test/integration/test_permeability_order_state.py runs the target-swap and mirror moves at a high rate on the input in examples/tests/test-internal/permeability/targetswap-validation/ and checks that every stored path is measured on one molecule and on one side of the mirror, and that a run restarted part-way matches an uninterrupted one.

Note

For slow permeants where even RETIS struggles to reach the permeation timescale, the partial-path variant REPPTIS combines short paths with replica exchange; see the REPPTIS example [2] and its Markov-state-model kinetics analysis [3].

References

[1] (1,2)

A. Ghysels, S. Roet, S. Davoudi and T. S. van Erp, Exact non-Markovian permeability from rare event simulations, Phys. Rev. Research 3, 033068 (2021), https://doi.org/10.1103/PhysRevResearch.3.033068

[2]

W. Vervust, D. T. Zhang, T. S. van Erp and A. Ghysels, Path sampling with memory reduction and replica exchange to reach long permeation timescales, Biophys. J. 122, 2960-2972 (2023), https://doi.org/10.1016/j.bpj.2023.02.021

[3]

W. Vervust, E. Wils, S. Safaei, D. T. Zhang and A. Ghysels, Estimating full path lengths and kinetics from partial path transition interface sampling simulations, J. Chem. Theory Comput. (2026), https://doi.org/10.1021/acs.jctc.5c01498 (arXiv:2602.12835).