.. _examples-permeability:

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
:file:`examples/tutorials/path_sampling/permeability/` and consist of
:file:`retis_perm.rst`, :file:`flat_potential.py`, and
:file:`initial.xyz`.

**Verification status:** smoke -- see :ref:`example-tutorial-map`.

Tutorial quick start
--------------------

* **Best starting point:** :file:`examples/tutorials/path_sampling/permeability/`.
* **Edit first:** the ``Simulation`` section for ``permeability``,
  ``zero_left``, and interfaces; the ``Orderparameter`` section for the
  tracked particle and mirror position.
* **Run:** ``pyretisrun -i retis_perm.rst -p``.
* **Analyse:** ``pyretisanalyse -i out.rst`` after the run has created
  ``out.rst``.
* **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 :ref:`example-test-status` for current
  coverage.

New simulation settings
-----------------------
The ``Simulation`` section has a few extra options:

.. literalinclude:: /_static/examples/permeability/retis_perm.rst
   :language: rst
   :lines: 9-12

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 ``pyretisanalyse``
    to calculate
    :math:`\xi`, :math:`\frac{\tau}{dz}` and the permeability.

  * ``swap_attributes``: A list of ``path_ensemble`` attributes that need to
    be swapped with every REPEX swap. The newly implemented Monte Carlo moves
    alter the order parameter in ``[0^-]``. These altered order parameters
    need to be swapped with the path whenever the path is exchanged between
    ensembles.


New TIS settings
----------------
The new ``mirror`` and ``target-swap`` moves add a few options to the
``TIS`` section:

.. literalinclude:: /_static/examples/permeability/retis_perm.rst
   :language: rst
   :lines: 43-45

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
``Permeability``. This class is a subclass of ``Position``, but alters the
output depending on
``mirror_pos``, ``relative`` and ``offset``.

.. literalinclude:: /_static/examples/permeability/retis_perm.rst
   :language: rst
   :lines: 73-80

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 ``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`` 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 pathensemble.txt.
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
``pyretisanalyse`` to also calculate and plot :math:`\xi`,
:math:`\frac{\tau}{dz}` and the permeability. This follows the formulas
described in the paper [1]_.
For the calculation of :math:`\tau` to work, a reference region has to be
chosen. This is done by adding ``tau_ref_bin`` to the ``Analysis`` section of
the ``.rst`` file.

.. literalinclude:: /_static/examples/permeability/retis_perm.rst
   :language: rst
   :lines: 148-149


This value can be changed in ``out.rst`` and the analysis can then be 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
---------

This tutorial currently documents a specialized setup and is not yet
covered by a dedicated fixture under :file:`examples/tests/`. Before a
release depends on this workflow, add a short permeability check that
runs the input file and verifies that the permeability analysis entries
are produced.


.. 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
   :ref:`REPPTIS example <examples-repptis>` [2]_ and its
   Markov-state-model kinetics analysis [3]_.

References
----------

.. [1] 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).