Partial path TIS (REPPTIS) in a 1D potential¶
This example shows how to run REPPTIS – Replica Exchange Partial Path Transition Interface Sampling – on the 1D double-well benchmark. REPPTIS cuts the sampled trajectories short: instead of requiring a path to connect the two stable states, each path only has to cross from one interface neighbourhood to the next. The short partial paths are then combined with replica exchange between the partial-path ensembles, which makes it possible to reach permeation timescales that ordinary RETIS cannot. The method and its memory-reduction / replica-exchange scheme are described by Vervust, Zhang, van Erp and Ghysels [1].
The closed-form kinetics (mean first passage time, flux and rate) that can be recovered from the short partial paths via a Markov-state-model analysis are derived by Vervust, Wils, Safaei, Zhang and Ghysels [2]; REPPTIS for permeation builds on the permeability-from-(RE)TIS formalism of Ghysels, Roet, Davoudi and van Erp [3].
Verification status: smoke – see Example test status.
Tutorial quick start¶
- Best starting point:
examples/tutorials/path_sampling/1D-double-well/repptis/. - Edit first:
task = "repptis"and theinterfaceslist in the[simulation]section; the[retis]swapfreq/swapsimuloptions control the partial-path replica exchange. - Run:
pyretisrun -i repptis.toml -pfrom the tutorial folder. - Analyse:
pyretisanalyse -i repptis.tomlafter the run has created its output. - Expected output: the standard per-ensemble folders
(
000…00N) withpathensemble.txt/order.txt, where the partial paths are shorter than the full RETIS paths. - Related checks:
examples/tests/test-internal/partial-path/(internal Langevin engine) andexamples/tests/test-gromacs/test-repptis/(real GROMACS); see Example test status.
Selecting the partial-path task¶
REPPTIS is requested by setting task = "repptis" in the
simulation section. Apart from the task
name, the input file is structured exactly like a RETIS input – the same
interfaces, [engine], [orderparameter] and [retis]
sections apply:
[simulation]
task = "repptis"
steps = 20000
interfaces = [
-0.9,
-0.8,
-0.7,
-0.6,
-0.5,
-0.4,
-0.3,
1.0,
]
The interfaces list defines the partial-path ensembles in the same
way as for RETIS: one \([0^{-}]\) ensemble and one \([i^{+}]\)
ensemble per intermediate interface. The difference is in how far each
trajectory is propagated – a REPPTIS path is terminated as soon as it
reaches a neighbouring interface, rather than being extended all the way
to a stable state.
For a permeation setup, combine task = "repptis" with
permeability = True and zero_left as shown in the
permeability example; the
examples/tests/test-internal/partial-path/ fixture runs exactly
this combination on a flat 1D potential.
Running the simulation¶
Running REPPTIS works the same way as RETIS. Below is the complete input file for the 1D double-well:
Show/hide the full input file »
potential = [
{ class = "DoubleWell", a = 1.0, b = 2.0, c = 0.0 },
]
[simulation]
task = "repptis"
steps = 20000
interfaces = [
-0.9,
-0.8,
-0.7,
-0.6,
-0.5,
-0.4,
-0.3,
1.0,
]
[system]
units = "reduced"
dimensions = 1
temperature = 0.07
[engine]
class = "Langevin"
timestep = 0.002
gamma = 0.3
high_friction = false
seed = 0
[box]
periodic = [
false,
]
[particles]
name = [
"Ar",
]
ptype = [
0,
]
[particles.position]
input_file = "initial.xyz"
[particles.velocity]
generate = "maxwell"
momentum = false
seed = 0
[particles.mass]
Ar = 1.0
[forcefield]
description = "1D double well"
[orderparameter]
class = "Position"
dim = "x"
index = 0
periodic = false
[output]
trajectory-file = 100
energy-file = 100
order-file = 100
[tis]
freq = 0.5
maxlength = 20000
allowmaxlength = false
zero_momentum = false
rescale_energy = false
sigma_v = -1
seed = 0
[initial-path]
method = "kick"
kick-from = "initial"
[retis]
swapfreq = 0.5
nullmoves = true
swapsimul = true
The initial configuration
initial.xyz
is included in the tutorial folder. The simulation is executed with:
pyretisrun -i repptis.toml -p
The -p option displays a progress bar.
Recovering kinetics from partial paths¶
Because REPPTIS paths are cut short, the crossing probability,
mean first passage time, flux and rate cannot be read off directly the
way they are for full RETIS paths. They are instead reconstructed with a
Markov-state-model (MSM) analysis of the partial-path ensembles, using
the closed formulas of Vervust, Wils, Safaei, Zhang and Ghysels [2].
PyRETIS provides that analysis in
pyretis.analysis.repptis_msm: the per-ensemble local
crossing probabilities (the LML / LMR / RML / RMR path-type fractions)
are assembled into a transition matrix with
construct_M(), the overall
crossing probability follows from
global_pcross_msm(), and the
mean first passage time – hence the flux and rate – from
mfpt_to_first_last_state().
The runnable tutorial msm_kinetics.py in this folder walks
through the whole reconstruction. It verifies the exact symmetric
limit pcross = 1 / (2 * (N - 1)) and then computes the crossing
probability, MFPT and rate for a representative asymmetric model:
python msm_kinetics.py
Show/hide the MSM tutorial script »
#!/usr/bin/env python
# Copyright (c) 2026, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""Reconstruct REPPTIS kinetics with the PyRETIS MSM analysis.
A REPPTIS (replica exchange partial path TIS) simulation samples
*short* partial paths, so the overall crossing probability, mean first
passage time (MFPT), flux and rate cannot be read off the path
ensembles directly. They are reconstructed with a Markov-state-model
(MSM) built from the per-ensemble *local* crossing probabilities,
following Vervust, Wils, Safaei, Zhang and Ghysels, J. Chem. Theory
Comput. (2026), https://doi.org/10.1021/acs.jctc.5c01498.
This tutorial demonstrates that analysis with
:py:mod:`pyretis.analysis.repptis_msm`. The starting point is the set
of local crossing probabilities for every partial-path ensemble:
* ``p_mp`` -- LMR: enters from the Left, exits to the Right
(a forward crossing),
* ``p_mm`` -- LML: enters from the Left, returns Left,
* ``p_pm`` -- RML: enters from the Right, exits Left,
* ``p_pp`` -- RMR: enters from the Right, returns Right.
In a real workflow these four numbers per ensemble come from analysing
a REPPTIS run (see the REPPTIS tutorial ``repptis.toml`` in this folder
and ``pyretis.analysis.analyse_repptis_ensemble``). Here we use a
representative set for the 1D double well so the tutorial runs in a
fraction of a second, and we cross-check against the exact symmetric
limit.
Run with::
python msm_kinetics.py
"""
import numpy as np
from pyretis.analysis.repptis_msm import (
construct_M,
construct_tau_vector,
global_pcross_msm,
mfpt_to_first_last_state,
)
def report_crossing_probability(p_mm, p_mp, p_pm, p_pp, n_interfaces):
"""Build the MSM and print the global crossing probability.
Parameters
----------
p_mm, p_mp, p_pm, p_pp : list of float
The per-ensemble local crossing probabilities (length
``n_interfaces - 1``) for the LML, LMR, RML and RMR path types.
n_interfaces : int
The number of interfaces.
Returns
-------
matrix : numpy.ndarray
The MSM transition matrix (reused for the kinetics below).
pcross : float
The global A-to-B crossing probability.
"""
matrix = construct_M(p_mm, p_mp, p_pm, p_pp, n_interfaces)
_, _, y1, _ = global_pcross_msm(matrix)
pcross = float(y1[0, 0])
print(f" global crossing probability P(A->B) = {pcross:.6e}")
return matrix, pcross
def report_kinetics(matrix, pcross, tau_mm, tau_mp, tau_pm, tau_pp,
n_interfaces, flux):
"""Print the rate and the mean first passage time.
The TIS/RETIS rate constant is the product of the initial flux
through the first interface and the overall crossing probability,
``k = flux * pcross``. The mean first passage time (MFPT) is a
separate reconstructed observable, obtained from the MSM together
with the per-ensemble transition times.
Parameters
----------
matrix : numpy.ndarray
The MSM transition matrix from
:py:func:`report_crossing_probability`.
pcross : float
The global A-to-B crossing probability.
tau_mm, tau_mp, tau_pm, tau_pp : list of float
The per-ensemble transition times (length ``n_interfaces``)
for the four path types, in simulation time units.
n_interfaces : int
The number of interfaces.
flux : float
The initial flux through the first interface (crossings per
unit time), as obtained from the ``[0^-]`` flux analysis.
Returns
-------
rate : float
The reconstructed rate constant ``flux * pcross``.
mfpt : float
The reconstructed mean first passage time, in time units.
"""
rate = flux * pcross
print(f" reconstructed rate k = flux * pcross = {rate:.6e}")
n_states = 4 * n_interfaces - 5
tau_vector = construct_tau_vector(
n_interfaces, n_states, tau_mm, tau_mp, tau_pm, tau_pp,
)
# The MFPT machinery splits each visit into a "before", "middle"
# and "after" contribution. For this illustrative example we feed
# the same per-state transition time to each part; a full analysis
# supplies the three measured distributions separately.
_, _, h1, _ = mfpt_to_first_last_state(
matrix, tau_vector, tau_vector, tau_vector,
)
mfpt = float(h1[0, 0])
print(f" mean first passage time (steps) = {mfpt:.6e}")
return rate, mfpt
def main():
"""Run the REPPTIS MSM kinetics demonstration."""
# ----------------------------------------------------------------
# 1) A symmetric, memoryless REPPTIS: every local crossing
# probability is 1/2. The MSM crossing probability then has the
# exact closed form 1 / (2 * (N - 1)), which we verify.
# ----------------------------------------------------------------
n_interfaces = 5
half = [0.5] * (n_interfaces - 1)
print("Symmetric memoryless model (all local probabilities = 0.5):")
_, pcross = report_crossing_probability(
half, half, half, half, n_interfaces,
)
exact = 1.0 / (2.0 * (n_interfaces - 1))
print(f" exact 1/(2*(N-1)) = {exact:.6e}")
assert np.isclose(pcross, exact), "symmetric limit not reproduced"
# ----------------------------------------------------------------
# 2) A representative asymmetric set, as would come from analysing
# a REPPTIS run of the 1D double well (forward crossings get
# harder towards the barrier, so p_mp decreases).
# ----------------------------------------------------------------
p_mp = [0.45, 0.35, 0.30, 0.40] # LMR (forward)
p_mm = [1.0 - value for value in p_mp] # LML (return Left)
p_pp = [0.80, 0.75, 0.70, 0.65] # RMR (return Right)
p_pm = [1.0 - value for value in p_pp] # RML (cross Left)
print("\nRepresentative asymmetric REPPTIS model:")
matrix, pcross = report_crossing_probability(
p_mm, p_mp, p_pm, p_pp, n_interfaces,
)
# ----------------------------------------------------------------
# 3) Combine the crossing probability with the per-ensemble
# transition times and the initial flux to get the rate.
# ----------------------------------------------------------------
tau_mm = [12.0, 14.0, 16.0, 18.0, 20.0]
tau_mp = [40.0, 45.0, 50.0, 55.0, 60.0]
tau_pm = [38.0, 42.0, 46.0, 50.0, 54.0]
tau_pp = [10.0, 11.0, 12.0, 13.0, 14.0]
flux = 1.0e-3 # crossings per step, from the [0^-] flux analysis
print("\nKinetics from the asymmetric model:")
report_kinetics(matrix, pcross, tau_mm, tau_mp, tau_pm, tau_pp,
n_interfaces, flux)
print("\nDone. In a real study the local crossing probabilities "
"and transition\ntimes are measured from a REPPTIS run "
"(see repptis.toml in this folder),\nand the flux from the "
"[0^-] ensemble analysis.")
if __name__ == "__main__":
main()
From a run to the rate (end to end)¶
repptis_to_kinetics.py closes the loop from sampling to
theory: it analyses an actual REPPTIS run for the per-ensemble local
crossing probabilities (LML / LMR / RML / RMR, via
analyse_repptis_ensemble() and
local_crossing_probabilities())
and combines them into the overall crossing probability two
independent ways – the closed-form REPPTIS recursion
(global_crossing_probabilities_from_local())
and the Markov state model
(construct_M() +
global_pcross_msm()) – which
agree. The full workflow is:
pyretisrun -i repptis.toml -p # sample the partial paths
python repptis_to_kinetics.py # reconstruct the kinetics
If run without simulation output present, the script falls back to a small set of illustrative local probabilities so it still demonstrates the analysis.
Show/hide the end-to-end script »
#!/usr/bin/env python
# Copyright (c) 2026, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""From a REPPTIS simulation to its kinetics, end to end.
This tutorial closes the loop from *sampling* to *theory*: it takes the
per-ensemble output of a REPPTIS run, extracts the local crossing
probabilities of each ensemble, and combines them into the overall
crossing probability two independent ways -- the closed-form REPPTIS
recursion and the Markov-state-model -- which must agree. It follows the
analysis of Vervust, Wils, Safaei, Zhang and Ghysels, J. Chem. Theory
Comput. (2026), https://doi.org/10.1021/acs.jctc.5c01498.
Workflow
--------
1. Run the REPPTIS simulation in this folder::
pyretisrun -i repptis.toml -p
This fills the per-ensemble directories ``000/``, ``001/``, ... with
``pathensemble.txt`` files.
2. Reconstruct the kinetics::
python repptis_to_kinetics.py
If no run output is found, the script falls back to a small set of
*illustrative* local crossing probabilities so that it still
demonstrates the analysis (and runs in the tutorial smoke suite); the
banner makes clear which mode is in effect.
"""
import os
import numpy as np
from pyretis.analysis import (
analyse_repptis_ensemble,
construct_M,
global_crossing_probabilities_from_local,
global_pcross_msm,
)
from pyretis.core.pathensemble import generate_ensemble_name
from pyretis.inout.formats.pathensemble import PathEnsembleFile
from pyretis.inout.settings import parse_settings_file
from pyretis.setup.createsimulation import create_ensembles
INPUT_FILE = "repptis.toml"
ANALYSIS_SETTINGS = {
"ngrid": 1000,
"maxblock": 100,
"bins": 100,
"blockskip": 1,
"skip": 0,
}
# Illustrative per-ensemble local crossing probabilities, used only when
# no run output is present (so the tutorial still demonstrates the
# analysis). Order: LMR (forward), LML (return left), RMR (return
# right), RML (cross left), one entry per [i^+] ensemble.
EXAMPLE_LOCAL_PROBS = {
"LMR": [0.45, 0.35, 0.30, 0.40],
"LML": [0.55, 0.65, 0.70, 0.60],
"RMR": [0.80, 0.75, 0.70, 0.65],
"RML": [0.20, 0.25, 0.30, 0.35],
}
def collect_local_probs_from_run():
"""Analyse each ``[i^+]`` ensemble of a run for its local probs.
Returns
-------
probs : dict or None
Lists of per-ensemble local probabilities keyed ``LMR`` /
``LML`` / ``RMR`` / ``RML``, or None if the run output is
missing.
n_interfaces : int
The number of interfaces (needed to size the MSM).
"""
settings = parse_settings_file(INPUT_FILE)
settings.setdefault("analysis", {}).update(ANALYSIS_SETTINGS)
ensembles = create_ensembles(settings)
n_interfaces = len(settings["simulation"]["interfaces"])
probs = {"LMR": [], "LML": [], "RMR": [], "RML": []}
for ensemble in ensembles:
path_ensemble = ensemble["path_ensemble"]
number = path_ensemble.ensemble_number
if number == 0: # the [0^-] ensemble is not an [i^+-] ensemble
continue
filename = os.path.join(
generate_ensemble_name(number), "pathensemble.txt"
)
if not os.path.isfile(filename):
return None, n_interfaces
ens_file = PathEnsembleFile(
filename, "r",
ensemble_settings={"ensemble_number": number,
"interfaces": path_ensemble.interfaces},
)
result = analyse_repptis_ensemble(
ens_file,
{"tis": {"detect": path_ensemble.interfaces[-1],
"ensemble_number": number},
"analysis": ANALYSIS_SETTINGS},
)
local = result["local_probs"]
for key in ("LMR", "LML", "RMR", "RML"):
probs[key].append(local[key])
return probs, n_interfaces
def report_crossing_probability(probs, n_interfaces):
"""Print the global crossing probability from the local probs.
Computes it twice -- via the closed-form REPPTIS recursion and via
the Markov state model -- and checks the two agree.
Parameters
----------
probs : dict
Per-ensemble local probabilities keyed ``LMR`` / ``LML`` /
``RMR`` / ``RML``.
n_interfaces : int
The number of interfaces.
"""
p_mp, p_mm = probs["LMR"], probs["LML"]
p_pp, p_pm = probs["RMR"], probs["RML"]
print("Per-ensemble local crossing probabilities:")
for i, (mp, mm, pp, pm) in enumerate(zip(p_mp, p_mm, p_pp, p_pm)):
print(f" [{i + 1}^+] LMR={mp:.3f} LML={mm:.3f} "
f"RMR={pp:.3f} RML={pm:.3f}")
recursion = global_crossing_probabilities_from_local(
p_mp, p_mm, p_pp, p_pm
)
if np.any(np.isnan(recursion[0])):
print("\nSome ensemble was not sampled in both directions "
"(nan local probability); collect more REPPTIS cycles "
"before reconstructing the kinetics.")
return
p_cross = recursion[2]
print(f"\nGlobal crossing probability (recursion) = {p_cross[-1]:.6e}")
# The MSM uses N-1 local-prob entries per type.
matrix = construct_M(p_mm[:n_interfaces - 1], p_mp[:n_interfaces - 1],
p_pm[:n_interfaces - 1], p_pp[:n_interfaces - 1],
n_interfaces)
_, _, y1, _ = global_pcross_msm(matrix)
print(f"Global crossing probability (MSM) = {float(y1[0, 0]):.6e}")
print("The two independent estimates agree:",
np.isclose(p_cross[-1], float(y1[0, 0])))
def main():
"""Reconstruct the REPPTIS kinetics from a run (or example probs)."""
probs, n_interfaces = collect_local_probs_from_run()
if probs is None:
print("No REPPTIS run output found in this folder; "
"demonstrating with illustrative local probabilities.")
print("Run 'pyretisrun -i repptis.toml -p' first to analyse a "
"real simulation.\n")
probs = EXAMPLE_LOCAL_PROBS
n_interfaces = len(probs["LMR"]) + 1
else:
print("Analysed the REPPTIS run output in this folder.\n")
report_crossing_probability(probs, n_interfaces)
if __name__ == "__main__":
main()
Tested by¶
The REPPTIS task is represented by
examples/tutorials/path_sampling/1D-double-well/repptis/ and
checked by examples/tests/test-internal/partial-path/ (internal
engine) and examples/tests/test-gromacs/test-repptis/ (GROMACS).
The MSM kinetics tutorial msm_kinetics.py runs in the
tutorial smoke suite (examples/run-tutorials.sh, exercised by
test-heavy) and the underlying analysis is unit-tested in
test/analysis/test_repptis_msm.py. Use the fixtures as the
short regression checks and this page as the user-facing guide.
References¶
| [1] | 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 |
| [2] | (1, 2) 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). |
| [3] | 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 |