Python API 

fastmdxplora.simulation.metad_surface.compute_surface_2d(hills_path, colvar_values=None, *, points=80, minimum_recrossings=4, periodic=(False, False), temperature_K=300.0, names=('cv1', 'cv2'))[source]

A two-dimensional surface, judged one dimension at a time.

The gates are the one-dimensional ones, applied to the free energy along each variable with the other integrated out. That is the whole point of doing it this way: a run can fill a torsion thoroughly while the distance it was also biasing never left one basin, and a verdict on the surface as a whole reports either a pass that hides the second coordinate or a failure that buries the first. The refusal, when there is one, names the dimension.

colvar_values is the trajectory of both variables, shape (frames, 2). Without it no recrossing count can be made in either dimension, and that is said rather than skipped.

The default grid is coarser than the one-dimensional default because the cost is the square of it: 80 by 80 is 6,400 points where 200 by 200 is 40,000, and the second resolves nothing a metadynamics surface supports.

Parameters:

hills_path (str | Path)
colvar_values (ndarray | None)
points (int)
minimum_recrossings (int)
periodic (tuple[bool, bool])
temperature_K (float)
names (tuple[str, str])

Return type:

fastmdxplora.simulation.metad_surface.marginal_profile(surface, axis, *, temperature_K=300.0)[source]

The free energy along one variable, with the others integrated out.

Not the minimum along the other coordinate, which is the projection people usually draw: that reports the bottom of the valley rather than how much room there is in it, so a broad shallow channel and a narrow deep one come out alike. Integrating the populations keeps the width, which is what the free energy along a single coordinate means.

Parameters:

surface (ndarray)
axis (int)
temperature_K (float)

Return type:

fastmdxplora.simulation.metad_surface.read_hills(path)[source]

Read a PLUMED HILLS file, in as many variables as it holds.

PLUMED writes a #! FIELDS header naming every column, and that header is what makes the file readable for more than one variable: the layout is time, then one centre per variable, then one sigma per variable, then height and the bias factor, so a two-variable file has the height in column six where a one-variable file has the bias factor. Reading by position, which is what this did, silently takes a sigma for a height on any run that biases two things.

The sigma columns are found by name, so the count of variables comes from the file rather than from an assumption about it. A file with no header is read by position as one variable, which is what such a file used to be.

Parameters:: path (str | Path)
Return type:: Hills

fastmdxplora.simulation.metad_surface.recrossings(values, *, low, high)[source]

How many times the variable travelled from one side to the other.

Counted with a hysteresis band: a crossing is only counted once the variable has reached the far region, so a run rattling around a threshold does not accumulate crossings it did not make. That is the difference between a barrier observed several times and a coordinate jittering at the top of one.

Parameters:

values (ndarray)
low (float)
high (float)

Return type:

int

fastmdxplora.simulation.metad_surface.surface_from_hills(hills, grid, *, upto=None, periodic=False)[source]

The free energy on grid, from the first upto hills.

The bias is the sum of the deposited Gaussians. For a well-tempered run the bias approaches -(1 - 1/γ) F, so the free energy is the bias scaled by γ/(γ - 1) and negated. Without tempering the bias approaches -F directly, and the scaling would divide by zero, so that case is its own.

Shifted so the lowest point is zero, because only differences mean anything: the absolute value of a free energy from metadynamics is set by where the sum happened to start.

Parameters:

hills (Hills)
grid (ndarray)
upto (int | None)
periodic (bool)

Return type:

fastmdxplora.simulation.metad_surface.surface_from_hills_nd(hills, axes, *, upto=None, periodic=None, chunk=512)[source]

The free energy on a grid of two or more variables.

The same arithmetic as the one-dimensional case, with the Gaussian a product over dimensions. Summed in chunks of hills because the direct form allocates one number per hill per grid point: a 200 by 200 grid and ten thousand hills is four hundred million floats, which is several gigabytes for an intermediate nobody reads.

Periodicity is per dimension. A run biasing a torsion against a distance is periodic in one and not the other, and treating either the way the other wants gives a surface that is wrong at one edge.

Parameters:

hills (Hills)
axes (tuple[ndarray, ...])
upto (int | None)
periodic (tuple[bool, ...] | None)
chunk (int)

Return type:

A free energy along a coordinate, from equilibrium sampling at each point.

Steered molecular dynamics drags a system along a coordinate and reports the work, which depends on how fast it was dragged. Umbrella sampling does the opposite: it holds the system at a series of positions, lets each equilibrate, and recombines the sampling into a potential of mean force. Nothing is hurried, so nothing is dissipated, and the result is a free energy rather than an upper bound on one.

It works only if adjacent windows sample overlapping ranges. The recombination stitches histograms together, and where two neighbours never visit the same value there is nothing to stitch: the free energy on one side cannot be placed relative to the other, and the curve through that gap is interpolation dressed as a measurement. Overlap is therefore checked and a gap is reported rather than bridged.

That check is the reason this module exists rather than a call to an external WHAM program. Every implementation computes the same PMF; few of them say when the windows could not support one.

Where the starting structures come from matters. Windows started from a single structure are strained at the far end of the range, and the strain relaxes into the sampling as drift. The usual source is a steered run: pull once, take a frame near each window’s centre, and each window begins near where it will sit.

class fastmdxplora.simulation.umbrella.UmbrellaPlan(windows, collective_variable, equilibration_fraction=0.2, minimum_overlap=0.03, minimum_samples=200)[source]

Bases: object

Where the windows are, and how hard each holds.

Parameters:

windows (tuple[Window, ...])
collective_variable (str)
equilibration_fraction (float)
minimum_overlap (float)
minimum_samples (int)

as_record()[source]

Return type:: dict[str, Any]

collective_variable: str

equilibration_fraction: float = 0.2: Steps to discard at the start of each window before it counts. A window begins away from where it will settle, and counting that approach as sampling biases the histogram towards where it started. The fraction of each window’s sampling discarded before its histogram is built. A window begins away from where it will settle and the approach is not sampling, so some must go; how much is a judgement about how long a window takes to settle, which depends on the barrier and the force constant, so it belongs to whoever is making the claim rather than to this file. A fifth is a common choice and the default.

minimum_overlap: float = 0.03: How much two neighbours must share before their free energies can be placed relative to one another. Three per cent is enough to stitch and thin: on a real study, pairs at seven per cent passed while a reader might reasonably want fifteen. It is a judgement about how much evidence a joint needs, so it belongs to whoever is making the claim.

minimum_samples: int = 200: How many values a window must have recorded before its histogram means anything. An overlap is the area two histograms share, and a histogram from tens of points is mostly noise – so a run short enough to be a smoke test will produce overlaps, and gaps, that are arithmetic rather than evidence. Like the overlap threshold, it is a judgement, so a study can set its own.

windows: tuple[Window, ...]

class fastmdxplora.simulation.umbrella.Window(index, centre, force_constant)[source]

Bases: object

One position along the coordinate, held there by a spring.

Parameters:

index (int)
centre (float)
force_constant (float)

centre: float

force_constant: float

index: int

plumed_lines(cv_lines)[source]

Parameters:: cv_lines (list[str])
Return type:: str

fastmdxplora.simulation.umbrella.collect_samples(directories, *, equilibration_fraction=0.2)[source]

Read each window’s sampling back from its COLVAR file.

Takes the directories rather than working them out. The first version guessed – <output>/window_00/simulation/COLVAR – and the runs are actually at <output>/runs/window-00/simulation/COLVAR: under a runs directory, with the identifier slugged. Two mistakes in one path, neither visible until a real study finished and found nothing. The caller knows where it put things.

The first part of each window is discarded. A window begins away from where it will settle, and counting the approach as sampling biases the histogram towards where the run started – which is the one place the free energy is guaranteed not to be flat.

Parameters:

directories (dict[int, Any])
equilibration_fraction (float)

Return type:

dict[int, ndarray]

fastmdxplora.simulation.umbrella.compute_pmf(samples, plan, *, temperature_K=300.0, bins=60, minimum_overlap=None, bootstrap_resamples=200, bootstrap_seed=0, _edges=None)[source]

A potential of mean force, or a refusal saying why not.

samples maps a window’s index to the collective-variable values it sampled, after equilibration has been discarded.

The overlap between neighbours is checked first. Where a pair does not share ground, the free energy on one side cannot be placed relative to the other, and a curve drawn through the gap is interpolation presented as a measurement.

Parameters:

samples (dict[int, ndarray])
plan (UmbrellaPlan)
temperature_K (float)
bins (int)
minimum_overlap (float | None)
bootstrap_resamples (int)
bootstrap_seed (int)
_edges (ndarray | None)

Return type:

fastmdxplora.simulation.umbrella.expand_umbrella(config)[source]

Turn a config with an umbrella block into one with a run per window.

An umbrella job is one system held at many positions, which is the shape the batch machinery already runs – so the windows become systems entries differing in the position each holds, and the scheduling, parallelism and per-GPU pinning come from the code that already does those things.

Returns the config unchanged where there is no umbrella block.

Parameters:: config (dict[str, Any])
Return type:: dict[str, Any]

fastmdxplora.simulation.umbrella.overlap_between(a, b, bins=50)[source]

How much two windows’ sampling shares ground, from 0 to 1.

The overlap coefficient: the area shared by two normalised histograms. Zero means they never visited the same value, and no recombination can place one relative to the other.

Parameters:

a (ndarray)
b (ndarray)
bins (int)

Return type:

float

fastmdxplora.simulation.umbrella.plan_from_expanded(config)[source]

Rebuild the window set from a config that has already been expanded.

Returns None where this is not an umbrella study.

Parameters:: config (dict[str, Any])
Return type:: UmbrellaPlan | None

fastmdxplora.simulation.umbrella.plan_windows(spec)[source]

Read an umbrella block into a set of windows.

Positions may be given explicitly, or as a range and a count. The count is the decision that matters: too few and adjacent windows do not overlap, which no amount of sampling repairs.

Parameters:: spec (dict[str, Any])
Return type:: UmbrellaPlan

fastmdxplora.simulation.umbrella.windows_as_sweep(plan)[source]

The windows as a sweep, which is the shape the batch machinery runs.

An umbrella job is one system at many restraint positions, which is the same shape as a parameter sweep – so the runs are expanded and scheduled by the machinery that already does that, rather than by a second one.

Parameters:: plan (UmbrellaPlan)
Return type:: list[dict[str, Any]]

Pulling a system along a coordinate, and what that does and does not give.

Some things do not happen on their own within reach of a simulation. A ligand may take milliseconds to leave a pocket where a run lasts microseconds. Steered molecular dynamics attaches a spring to a collective variable and moves the spring’s anchor, dragging the system along whether or not it wants to go.

What this gives you is a pathway, not a free energy. The work done pulling depends on how fast you pull: drag a ligand out in a nanosecond and most of the work goes into pushing water aside and straining the protein, not into breaking the interactions you meant to measure. That dissipated work does not cancel, and a single fast pull overestimates the barrier, sometimes by a lot.

Jarzynski’s equality recovers a free energy from an ensemble of pulls, and the average is dominated by the rare low-work trajectories – so it needs many repeats and converges badly when the pulling is fast. FastMDXplora does not claim a free energy from a steered run, and reports the work so the claim can be made deliberately by somebody who has done the repeats.

What steered MD is genuinely good for is generating starting structures. Pull once, take frames along the way, and each becomes a window for umbrella sampling – which does give a free energy, from equilibrium sampling at each position rather than from work done in a hurry. That is the standard use, and it is why this exists before umbrella sampling does.

class fastmdxplora.simulation.steered.SteeredPlan(cv, to_value, from_value=None, force_constant=2000.0, steps=500000)[source]

Bases: object

A pull: from where, to where, how hard, and over how long.

Parameters:

cv (MetadynamicsPlan)
to_value (float)
from_value (float | None)
force_constant (float)
steps (int)

as_record()[source]

Return type:: dict[str, Any]

cv: MetadynamicsPlan: The coordinate, resolved to atoms, reusing the metadynamics layer.

force_constant: float = 2000.0

from_value: float | None = None

rate_per_ns(timestep_fs)[source]

How fast the anchor moves, in units of the variable per ns.

The number that decides whether the work means anything. Reported rather than checked against a threshold, because what counts as slow depends on the coordinate and the system – but an order of magnitude is usually obvious.

Parameters:: timestep_fs (float)
Return type:: float | None

steps: int = 500000: Steps over which the anchor travels. Everything about how much this run means depends on this being large.

to_value: float

fastmdxplora.simulation.steered.build_steered_script(plan, reference_pdb=None)[source]

The PLUMED input for a pull.

Written out rather than hidden, because the pulling rate is the thing that decides whether the result means anything and somebody checking a number should be able to read it.

Parameters:

plan (SteeredPlan)
reference_pdb (str | None)

Return type:

str

fastmdxplora.simulation.steered.plan_steered(spec, topology, *, temperature_K=300.0, ligand_resname=None)[source]

Read a steered-MD block, reusing the collective-variable machinery.

Parameters:

spec (dict[str, Any])
topology (Any)
temperature_K (float)
ligand_resname (str | None)

Return type:

SteeredPlan

Holding parts of a system still while the rest settles.

A structure that has just been minimised is not at equilibrium. Heating it lets the solvent find its arrangement, and it also lets the solute move – side chains relax into the vacuum the crystal packing left, a ligand drifts out of the pose that was measured, a membrane thins around a protein that has not yet found its depth. The conventional remedy is to hold the solute in place while the solvent equilibrates around it, and then let go in stages.

Without that, a run reaches production having already lost the arrangement it started from, and the trajectory answers a question about a structure nobody determined. This software had no restraints at all until now: it went from minimisation to unrestrained dynamics in one step.

Four kinds are implemented, which is what the protocols in common use ask for:

position – hold atoms near where they are, the workhorse of equilibration
distance – hold two atoms, or two groups, at a separation
angle – hold three atoms at an angle
torsion – hold four atoms at a dihedral

Each is a harmonic penalty: the force grows with the square of the departure, so a restraint is a spring rather than a wall. That matters for what a restrained run means. A constrained atom cannot move; a restrained one can, and the restraint says how much it cost. Reporting a restrained trajectory as though it were free is a claim the simulation does not support, so the restraints in force at each stage are recorded with the results.

Restraints do not survive into production. A biased production run measures the bias. They are released before it, and a run that keeps them must say so.

class fastmdxplora.simulation.restraints.AngleRestraint(selection: 'str', force_constant: 'float' = 1000.0, target: 'float | None' = None, kind: 'str' = 'angle')[source]

Bases: Restraint

Parameters:

selection (str)
force_constant (float)
target (float | None)
kind (str)

kind: str = 'angle': Set by each subclass. Declared last so a subclass overriding it does not leave a defaulted field ahead of a required one.

class fastmdxplora.simulation.restraints.DistanceRestraint(selection: 'str', force_constant: 'float' = 1000.0, target: 'float | None' = None, kind: 'str' = 'distance')[source]

Bases: Restraint

Parameters:

selection (str)
force_constant (float)
target (float | None)
kind (str)

kind: str = 'distance': Set by each subclass. Declared last so a subclass overriding it does not leave a defaulted field ahead of a required one.

class fastmdxplora.simulation.restraints.PositionRestraint(selection: 'str', force_constant: 'float' = 1000.0, target: 'float | None' = None, kind: 'str' = 'position')[source]

Bases: Restraint

Parameters:

selection (str)
force_constant (float)
target (float | None)
kind (str)

kind: str = 'position': Set by each subclass. Declared last so a subclass overriding it does not leave a defaulted field ahead of a required one.

class fastmdxplora.simulation.restraints.ReleaseSchedule(steps=(1000.0, 500.0, 100.0, 0.0))[source]

Bases: object

How a restraint weakens as equilibration proceeds.

Letting go all at once undoes the point of restraining: the solute is released into a solvent arrangement that formed around a structure held rigid, and the sudden freedom shows up as a jump in energy and a lurch in the structure. Standard protocols step the force down instead – 1000, 500, 100, 0 in kJ/mol/nm² is the shape of it – so each stage starts from something the previous one prepared.

The steps are spread across the equilibration stages. The last is applied to the end of NPT, and production always runs at zero unless somebody asks otherwise and is told what that means.

Parameters:: steps (tuple[float, ...])

as_record()[source]

Return type:: dict[str, Any]

force_at(fraction)[source]

The force constant at a point through equilibration.

fraction runs from 0 at the start to 1 at the end.

Parameters:: fraction (float)
Return type:: float

steps: tuple[float, ...] = (1000.0, 500.0, 100.0, 0.0)

class fastmdxplora.simulation.restraints.Restraint(selection, force_constant=1000.0, target=None, kind='position')[source]

Bases: object

What is held, and how firmly.

force_constant is in kJ/mol/nm² for position and distance restraints and kJ/mol/rad² for angle and torsion. The units differ because the coordinate does: a spring on a length and a spring on an angle are not measured in the same thing, and quietly using one number for both is how an angle restraint ends up a thousand times too weak.

Parameters:

selection (str)
force_constant (float)
target (float | None)
kind (str)

as_record()[source]

Return type:: dict[str, Any]

force_constant: float = 1000.0

kind: str = 'position': Set by each subclass. Declared last so a subclass overriding it does not leave a defaulted field ahead of a required one.

selection: str

target: float | None = None: Where to hold it. None means “where it is now”, which is what a position restraint almost always wants.

property units: str

class fastmdxplora.simulation.restraints.TorsionRestraint(selection: 'str', force_constant: 'float' = 1000.0, target: 'float | None' = None, kind: 'str' = 'torsion')[source]

Bases: Restraint

Parameters:

selection (str)
force_constant (float)
target (float | None)
kind (str)

kind: str = 'torsion': Set by each subclass. Declared last so a subclass overriding it does not leave a defaulted field ahead of a required one.

fastmdxplora.simulation.restraints.build_restraint_forces(omm, topology, positions, restraints)[source]

Turn restraints into OpenMM forces, paired with the parameter that scales each one.

The parameter is what makes staged release possible: a force added to a System cannot be removed once the context exists, but a global parameter can be set to zero, which is the same thing to every atom in the system and much cheaper than rebuilding.

Parameters:

omm (Any)
topology (Any)
positions (Any)
restraints (list[Restraint])

Return type:

list[tuple[Any, str]]

fastmdxplora.simulation.restraints.parse_restraints(spec)[source]

Read restraints from a config block.

Accepts the short form – a selection string, meaning position restraints on it at the default force – and the long form, a list of blocks. The short form is what an equilibration usually wants and should not require a paragraph to say.

Parameters:: spec (Any)
Return type:: list[Restraint]

A binding free energy from a potential of mean force along a distance.

The well depth of a PMF is not a binding free energy. Turning one into the other needs the standard state, because a dissociation constant is defined against a reference concentration and a curve is not: at 1 M the ligand has 1661 cubic Angstroms to itself, and the free energy of binding is the cost of giving that up. cite{gilson1997}

K = 4 pi * integral over the bound range of exp(-beta [A(r) - c]) dr dG = -kT ln(K / V0), V0 = 1661 A^3

with c the constant that A(r) + 2kT ln r settles to in bulk.

Which free energy A(r) is matters, and getting it wrong doubles an entropy. Umbrella sampling recombined by histogram gives the free energy of the radial distribution, in which the 4 pi r^2 volume element is already present: in bulk, where nothing is interacting, A(r) still falls as -2kT ln r simply because a larger shell holds more places to be. The true potential of mean force between the two centres is flat there. Writing the integral above against the distribution form is what makes the r^2 cancel, and the same expression applied to a Jacobian-removed PMF counts the translational entropy twice.

A binding free energy needs the run to have reached bulk. The reference A(r_bulk) is only a reference if the ligand is free there, and that is checkable rather than assumed: in bulk the curve must fall as -2kT ln r. A run whose windows stopped inside the interaction still produces a smooth PMF and a plausible number, and the number is wrong by however much of the well was left outside. Where the outer range does not have that shape, no binding free energy is reported.

The bound state is a choice, and its size is reported. Where the well ends is a definition, not a measurement. The same curve integrated to different cutoffs gives different answers, so the sensitivity across a range of reasonable cutoffs is reported beside the value: a well that is deep and narrow barely moves, and a shallow one moves a great deal, which is the reader’s cue about how much the definition is doing.

fastmdxplora.simulation.binding.STANDARD_VOLUME_NM3 = 1.660539: The volume one molecule has to itself at 1 mol/L, in nm^3. 1661 cubic Angstroms, the standard state a dissociation constant is quoted against.

fastmdxplora.simulation.binding.binding_free_energy(coordinate, free_energy_kjmol, *, temperature_K=300.0, bound_cutoff_nm=None, bulk_fraction=0.25)[source]

A standard-state binding free energy, or a refusal saying why not.

coordinate is the ligand-site distance in nm and free_energy_kjmol the recombined free energy along it, as compute_pmf returns them.

bound_cutoff_nm is where the bound state is taken to end. Left out, it is placed at the first point beyond the minimum where the curve has come within kT of its bulk value, which is a defensible reading of “no longer interacting” and is reported alongside the answer.

Parameters:

coordinate (Any)
free_energy_kjmol (Any)
temperature_K (float)
bound_cutoff_nm (float | None)
bulk_fraction (float)

Return type:

Analyses

Each is a class with the same shape: options in the constructor, compute() for the numbers, run() to write the data, the figure and the record of what it did.

Abstract base class for all FastMDXplora analysis modules.

Every analysis (RMSD, RMSF, Rg, …) subclasses Analysis and implements compute() and plot(). The base class handles everything else: output directory creation, options-manifest serialization, atom selection resolution, status tracking, and the run() convenience method that does compute() -> save_data() -> plot() -> save_figure() in order.

The contract is deliberately small:

compute(traj) returns a Python object (usually a numpy array or pandas DataFrame). It must be deterministic, side-effect-free, and inexpensive to call again with the same input.

plot(result, ax) draws onto a matplotlib Axes. It must not call plt.show() or close the figure — the caller controls that.

save_data(result, path) writes the computed result to disk. The default implementation handles numpy arrays and DataFrames; analyses with non-tabular output can override it.

Outputs land in <output_dir>/<analysis_name>/:

<output_dir>/
└── rmsd/
    ├── rmsd.dat       # the numerical data
    ├── rmsd.png       # the figure
    └── options.json   # parameter manifest for this analysis

class fastmdxplora.analysis.base.Analysis(*, selection=None, output_dir=None, title=None, xlabel=None, ylabel=None, figsize=None, xunit=None, **options)[source]

Bases: ABC

Base class for a single trajectory analysis module.

Subclasses must:

Set the class attribute name to a short identifier (e.g. "rmsd").
Implement compute() which takes an MDTraj trajectory and returns the analysis result.
Implement plot() which renders self.result onto an Axes.

Subclasses may override save_data() and default_selection.

Parameters:

selection (str | None)
output_dir (str | Path | None)
title (str | None)
xlabel (str | None)
ylabel (str | None)
figsize (tuple[float, float] | None)
xunit (str | None)
options (Any)

abstractmethod compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = '': Human-readable description used in figure titles.

figure_title()[source]

Title shown above the figure.

If the user passed title= at construction, that wins. Otherwise the subclass’s description is used; failing that, the analysis name in uppercase.

Return type:: str

findings: dict[str, Any]: What the analysis worked out while running, as opposed to what it was told. Recorded beside the options and kept out of them, because a report lists the options.

frame_axis(traj)[source]

Return (x_values, x_label) for a time-series plot.

Resolution order for the unit:

User-supplied xunit at construction ("ns", "ps", or "frames").
"ns" if the trajectory carries usable timing information (traj.time or traj.timestep).
"frames" as the always-safe fallback.

Parameters:

traj (mdtraj.Trajectory)

Returns:

x (np.ndarray, shape (n_frames,)) – The numerical values for the x axis.
label (str) – A pre-formatted axis label, e.g. "Time (ns)" or "Frame".

Return type:

Notes

MDTraj stores traj.time in picoseconds. When the trajectory lacks timing (e.g., loaded from a PDB without a timestep), the method falls back to frame indices.

honours_selection: bool = True: Whether selection means anything for this analysis. An analysis that decides its own atoms – a protein-ligand measure works out both sides from the ligand’s residue name, dihedrals from the backbone – has nothing to apply it to, and offering a control that does nothing is worse than not offering one: it looks like it worked.

min_atoms_to_align: int = 0

How many atoms this analysis’s selection must match before a rigid body superposition onto it is defined. Three is the minimum for a rotation; below that there is no unique answer.

Alanine dipeptide has one CA, which is the default selection for RMSD and RMSF, and MDTraj responded by printing “UNCONVERGED ROTATION MATRIX. RETURNING IDENTITY” once per frame from its C extension and returning distances measured against no alignment at all. Thousands of lines of it, and a column of numbers that looked like results.

name: str = 'analysis': Short identifier used in output paths and the manifest.

abstractmethod plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (Any)
ax (Axes)

Return type:

None

requires_ligand: bool = False: True for analyses that only apply to protein-ligand complexes (e.g. ligand pose RMSD). The orchestrator runs these automatically when a ligand is present and skips them otherwise. They can still be explicitly requested via include.

reweightable: tuple[str | None, str] | None = None

The per-frame scalar whose ensemble average means something, as (column, label). column names a column when compute() returns a DataFrame and is None when it returns a bare per-frame array. None for the attribute itself – the default – means no weighted average of this analysis is offered.

On a metadynamics run the trajectory is not a Boltzmann ensemble, so every mean reported from it is an average over a distribution the bias flattened on purpose. Declaring this lets that mean be recomputed against the deposited bias and reported beside the raw one. Analyses whose result is not one number per frame leave it unset: reweighting a clustering or a projection is a harder question than a weighted mean, and the report says so rather than guessing.

reweightable_populations: bool = False

Whether compute() returns per-frame categorical labels, either as an array or as a mapping of method name to array. A population is a weighted count of an indicator, so it reweights exactly as a mean does, and how often a state is visited is the thing a biased run distorts most.

What this does not correct is which states exist. The clustering was performed on the biased frames, so the groupings themselves are shaped by where the bias sent the system; reweighting says how often each was really visited, not that the right ones were found.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)
Return type:: AnalysisResult

save_data(result, path)[source]

Write the computed result to path.

Default behaviour:

1-D / 2-D numpy arrays: np.savetxt (whitespace-delimited).
pandas DataFrames: to_csv (comma-delimited).
Anything else: subclass must override.

Returns the path actually written.

Parameters:

result (Any)
path (Path)

Return type:

select_atoms(traj)[source]

Resolve selection to atom indices on a given trajectory.

Returns the full atom index array when selection is None. Raises ValueError if the selection matches zero atoms.

Parameters:: traj (Trajectory)
Return type:: ndarray

time_series: bool = False

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

class fastmdxplora.analysis.base.AnalysisResult(name, status, data=None, output_dir=None, figure_path=None, data_path=None, options_path=None, artifacts=<factory>, message='', started_at='', finished_at='')[source]

Bases: object

Lightweight record of one analysis invocation.

Returned by Analysis.run() and aggregated by the orchestrator. Includes both the computed data (for in-memory consumers) and the on-disk artifact paths (for report generation and provenance).

Parameters:

name (str)
status (str)
data (Any)
output_dir (Path | None)
figure_path (Path | None)
data_path (Path | None)
options_path (Path | None)
artifacts (list[Path])
message (str)
started_at (str)
finished_at (str)

artifacts: list[Path]

data: Any = None

data_path: Path | None = None

figure_path: Path | None = None

finished_at: str = ''

message: str = ''

name: str

options_path: Path | None = None

output_dir: Path | None = None

started_at: str = ''

status: str

to_dict()[source]

Serialize for inclusion in the analysis manifest.

Return type:: dict[str, Any]

fastmdxplora.analysis.base.superposed(traj, *, frame=0, atom_indices=None)[source]

Align a copy, and drop the box that no longer describes it.

Two hazards, both silent, both met on the same run.

mdtraj’s superpose rotates coordinates in place and returns the same object. An analysis that aligned therefore left every later analysis in the same run reading rotated coordinates, so a measure’s result depended on which other measures had run before it. Nothing in the record could show this: every setting was identical either way.

And rotation does not rotate unitcell_vectors. Minimum-image distances computed afterwards map atoms through a box that no longer corresponds to the frame, which does not fail – it answers, wrongly. On a 20 ns trypsin-benzamidine run the same ligand-protein pair measured 1.64 nm before alignment and 1.83 nm after, and pl_interactions reported 252 hydrophobic contacts in company against 10 alone.

So: align a copy, and remove the box. An analysis that wants periodic distances must take them from the unaligned trajectory, where the box is still true. Absent is better than stale, because stale is the one a caller cannot detect.

Analysis-level orchestrator.

AnalysisOrchestrator is the analysis-phase counterpart to the project-level fastmdxplora.FastMDXplora class. Its responsibility is to coordinate the individual Analysis modules: discover what’s available, validate the user’s options, execute the chosen subset in order, capture results and errors, and write a single phase-level manifest.

Architecturally the orchestrator follows a seven-phase pipeline (Aina & Kwan, JCC 2026):

Discovery — what analyses are available?

Validation — does the user’s options dict have the right shape?

Planning — apply include/exclude to produce the execution list.

Defaults — merge per-analysis defaults under user overrides.

Filtering — match kwargs to each analysis’s constructor signature.

Execution — run each analysis sequentially, catching errors.

Consolidation — write the manifest, return the result dict.

In FastMDXplora the orchestrator is constructed by the project-level FastMDXplora.analyze() method, but it also supports direct use as a standalone class (see fastmdxplora.AnalysisOrchestrator).

class fastmdxplora.analysis.orchestrator.AnalysisOrchestrator(trajectory, topology=None, *, output_dir=None, selection=None, scope='solute', ligand_resname=None, stride=None, first=None, last=None)[source]

Bases: object

Coordinate the execution of trajectory analysis modules.

The orchestrator loads the trajectory once at construction (matching the standard pattern) and holds it on self.traj. Subsequent calls to run() operate on that loaded trajectory.

Parameters:

trajectory (path, list of paths, or glob) – Trajectory file(s) to analyze. Passed verbatim to load_trajectory().
topology (path, optional) – Topology file. If omitted, auto-resolution is attempted (see load_trajectory()).
output_dir (path, optional) – Where to write per-analysis subdirectories. Defaults to ./fastmdx_analysis_<timestamp>.
selection (str, optional) – Default MDTraj selection string applied to every analysis that does not override it.
stride (int, optional) – Frame-selection parameters applied at load time.
first (int, optional) – Frame-selection parameters applied at load time.
last (int, optional) – Frame-selection parameters applied at load time.
scope (str)
ligand_resname (str | None)

Examples

Run all registered analyses with defaults:

from fastmdxplora.analysis import AnalysisOrchestrator

ao = AnalysisOrchestrator("traj.dcd", topology="top.pdb")
results = ao.run()

Selectively run RMSD and Rg with custom RMSD options:

results = ao.run(
    include=["rmsd", "rg"],
    options={"rmsd": {"ref": 0, "selection": "name CA"}},
)

Exclude expensive analyses on a quick first pass:

results = ao.run(exclude=["cluster", "dimred"])

analyze(*, include=None, exclude=None, options=None)[source]

Alias for run().

Parameters:

include (list[str] | None)
exclude (list[str] | None)
options (dict[str, dict[str, Any]] | None)

Return type:

dict[str, AnalysisResult]

run(*, include=None, exclude=None, options=None)[source]

Execute the planned analyses against self.traj.

Parameters:

include (list of str, optional) – Subset of analysis names to run. Mutually exclusive with exclude.
exclude (list of str, optional) – Subset of analysis names to skip.
options (dict, optional) – Per-analysis keyword arguments. Keys are analysis names (e.g. "rmsd"); values are dicts forwarded to the analysis constructor. Unrecognized kwargs are silently dropped so the orchestrator can be safely called with a superset of options.

Returns:

Mapping from analysis name to result, in execution order. Also stored on self.results.

Return type:

dict[str, AnalysisResult]

fastmdxplora.analysis.orchestrator.available_analyses()[source]

Return the names of all registered analyses, in registration order.

Return type:: tuple[str, …]

fastmdxplora.analysis.orchestrator.get_analysis_class(name)[source]

Look up a registered analysis class by name.

Parameters:: name (str)
Return type:: type[Analysis]

fastmdxplora.analysis.orchestrator.register_analysis(name, cls)[source]

Register an analysis class under a short name.

Called by each analysis module at import time. Idempotent — re-registering the same name with the same class is a no-op; re-registering a different class raises ValueError.

Parameters:

name (str)
cls (type[Analysis])

Return type:

None

Root-Mean-Square Deviation (RMSD).

Computes the per-frame RMSD of an MD trajectory against a chosen reference frame, with optional rigid-body alignment (Kabsch superposition) applied prior to the distance calculation. The atom subset used for both the alignment and the RMSD calculation is controlled by the selection attribute; by default the alpha-carbon backbone is used, which is the convention for protein conformational analysis.

The output figure is a time-series of RMSD vs. simulation time (or frame index if no timestep is available), with the chosen reference frame marked. Output data is a single-column rmsd.dat of RMSD values in nanometers (MDTraj’s native unit).

References

The implementation delegates to MDTraj’s mdtraj.rmsd(), which uses the QCP method of Theobald (Acta Cryst. A, 2005) — an O(N) algorithm for minimum RMSD without explicit eigendecomposition.

class fastmdxplora.analysis.rmsd.RMSD(*, ref=0, align=True, **kwargs)[source]

Bases: Analysis

Per-frame root-mean-square deviation.

Parameters:

ref (int, default 0) – Reference frame index for the RMSD calculation. Negative indices count from the end (-1 = last frame).
align (bool, default True) – If True (recommended), structurally align each frame to the reference before computing the distance. Without alignment, the result includes rigid-body rotation/translation which is rarely the quantity of interest.
selection (str, optional) – MDTraj atom selection string. Defaults to "name CA" (alpha carbons) for protein analysis. For all atoms, pass selection="all".
**kwargs – Standard base-class options (output_dir, title, xlabel, ylabel, figsize, xunit).

Examples

Default (CA atoms, reference = frame 0, aligned):

rmsd = RMSD()
rmsd.run(trajectory)

RMSD against the last frame, on heavy atoms:

rmsd = RMSD(ref=-1, selection="not element H")

Compute only, plot separately:

rmsd = RMSD()
data = rmsd.compute(trajectory)   # 1-D array in nanometers

The output file is rmsd.dat (single column, nm).

compute(traj)[source]

Compute the per-frame RMSD.

Returns:: RMSD in nanometers (MDTraj convention).
Return type:: np.ndarray of shape (n_frames,)
Parameters:: traj (Trajectory)

default_selection: str | None = 'name CA': Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Root-mean-square deviation': Human-readable description used in figure titles.

frame_axis_for_plot(result, traj)[source]

Return x-axis values + label, using the cached trajectory.

Parameters:

result (ndarray)
traj (Trajectory | None)

Return type:

min_atoms_to_align: int = 3: A superposition needs three atoms to be defined.

name: str = 'rmsd': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

reweightable: tuple[str | None, str] | None = (None, 'RMSD (nm)')

The per-frame scalar whose ensemble average means something, as (column, label). column names a column when compute() returns a DataFrame and is None when it returns a bare per-frame array. None for the attribute itself – the default – means no weighted average of this analysis is offered.

On a metadynamics run the trajectory is not a Boltzmann ensemble, so every mean reported from it is an average over a distribution the bias flattened on purpose. Declaring this lets that mean be recomputed against the deposited bias and reported beside the raw one. Analyses whose result is not one number per frame leave it unset: reweighting a clustering or a projection is a harder question than a weighted mean, and the report says so rather than guessing.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)

time_series: bool = True

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

Root-Mean-Square Fluctuation (RMSF).

Per-residue (default) or per-atom RMSF over the trajectory. The trajectory is first superposed onto a reference (frame 0 or a user-chosen reference) using the selected atom subset to remove rigid-body motion; the per-atom RMSF is then the standard deviation of each atom’s position around its mean. The per-residue RMSF reduces per-atom RMSF to one value per residue by averaging over the atoms that belong to each residue.

This is the standard “flexibility profile” plot used in nearly every MD publication — peaks indicate flexible loops/termini, troughs indicate rigid secondary structure.

class fastmdxplora.analysis.rmsf.RMSF(*, ref=0, per_residue=True, **kwargs)[source]

Bases: Analysis

Per-residue root-mean-square fluctuation.

Parameters:

ref (int, default 0) – Reference frame for the alignment (superposition) step. The choice affects only the bookkeeping; fluctuations are measured relative to each atom’s mean position over the full trajectory, which is invariant under rigid-body alignment.
per_residue (rmsf_nm) when) – If True, collapse the per-atom RMSF down to one value per residue by averaging over the residue’s atoms. If False, return the per-atom array (one value per selected atom).
selection (str, optional) – MDTraj atom selection. Defaults to "name CA" (alpha carbons) for protein analysis.
**kwargs – Standard base-class options.
Output
------
(residue_index (A two-column rmsf.dat)
per_residue
True (is)
(atom_index (or)
False. (rmsf_nm) when)

Examples

Standard per-residue plot for proteins:

rmsf = RMSF()
rmsf.run(trajectory)

Per-atom on backbone heavy atoms:

rmsf = RMSF(per_residue=False, selection="backbone and not element H")

compute(traj)[source]

Compute the RMSF.

Returns:: Two columns: index (residue or atom number) and RMSF in nm.
Return type:: np.ndarray, shape (N, 2)
Parameters:: traj (Trajectory)

default_selection: str | None = 'name CA': Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Root-mean-square fluctuation': Human-readable description used in figure titles.

min_atoms_to_align: int = 3: A superposition needs three atoms to be defined.

name: str = 'rmsf': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

Radius of Gyration (Rg).

Per-frame radius of gyration, a measure of overall molecular size and compactness. For a single-chain protein, Rg typically tracks folding state — unfolded conformations have higher Rg, compact native states have lower Rg.

The formula:

Rg(t) = sqrt( sum_i m_i * |r_i(t) - r_cm(t)|^2  /  sum_i m_i )

The radius of gyration is mass-weighted, which is what GROMACS’s gyrate, cpptraj’s radgyr and a published Rg all report:

Rg = sqrt( sum_i m_i |r_i - R|^2 / sum_i m_i ),   R = sum_i m_i r_i / sum_i m_i

It is computed here rather than by mdtraj.compute_rg(), which weights every atom equally unless told otherwise – and which, when given masses, still measures from the geometric centre rather than the centre of mass. Weighting equally counts each hydrogen for as much as each carbon, which on a protein is a few per cent away from the mass-weighted value and enough to disagree with a number someone is comparing against. Pass mass_weighted=False for the unweighted quantity.

class fastmdxplora.analysis.rg.Rg(*, by_chain=False, mass_weighted=True, **kwargs)[source]

Bases: Analysis

Per-frame radius of gyration.

Parameters:

mass_weighted (bool, default True) – Weight each atom by its mass and measure from the centre of mass, as the conventional definition does. False weights every atom equally.
by_chain (bool, default False) – If True, compute Rg separately for each chain in the topology (in addition to the whole-system Rg). Useful for multi-chain complexes where you want to track the compactness of each subunit.
selection (str, optional) – MDTraj atom selection. Defaults to None which means “use all atoms” — appropriate for Rg of the entire system. For a protein Rg in a solvated system, pass selection="protein".
**kwargs – Standard base-class options.
Output
------
frame) (Single-column rg.dat (Rg in nm per)
when (or multi-column)
by_chain=True (frame, total, chain0, chain1, ....)

compute(traj)[source]

Compute Rg per frame.

Returns:: If by_chain=False: shape (n_frames,), Rg in nm. If by_chain=True: shape (n_frames, 1+n_chains), columns are [Rg_total, Rg_chain0, Rg_chain1, ...].
Return type:: np.ndarray
Parameters:: traj (Trajectory)

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Radius of gyration': Human-readable description used in figure titles.

frame_axis_for_plot(traj, result)[source]

Parameters:

traj (Trajectory | None)
result (ndarray)

Return type:

name: str = 'rg': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

reweightable: tuple[str | None, str] | None = (None, 'Radius of gyration (nm)')

The per-frame scalar whose ensemble average means something, as (column, label). column names a column when compute() returns a DataFrame and is None when it returns a bare per-frame array. None for the attribute itself – the default – means no weighted average of this analysis is offered.

On a metadynamics run the trajectory is not a Boltzmann ensemble, so every mean reported from it is an average over a distribution the bias flattened on purpose. Declaring this lets that mean be recomputed against the deposited bias and reported beside the raw one. Analyses whose result is not one number per frame leave it unset: reweighting a clustering or a projection is a harder question than a weighted mean, and the report says so rather than guessing.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)

time_series: bool = True

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

Solvent-Accessible Surface Area (SASA).

Per-frame SASA computed with the Shrake-Rupley rolling-sphere algorithm (MDTraj’s mdtraj.shrake_rupley()), for the whole molecule, for each residue per frame, or as each residue’s mean over the run – which is the summary that says which residues are buried. Outputs the total SASA time series and, optionally, a per-residue heatmap that shows which residues become exposed/buried over the simulation.

SASA is a sensitive probe of conformational changes that involve burial or exposure of hydrophobic surfaces — it can detect folding/unfolding events, partial unfolding of loops, and binding/unbinding transitions that don’t necessarily show up in RMSD.

References

Shrake, A.; Rupley, J. A. J. Mol. Biol. 1973, 79, 351.

fastmdxplora.analysis.sasa.ATTEMPTS = 5

CI returned answers truncated on both attempts. How often a call is truncated, and whether attempts are independent, is not yet known – see the characterisation test in the suite – so this is a number chosen to be cheap rather than one derived from a measured rate.

Type:: How many times to ask before giving up. Two was not enough

class fastmdxplora.analysis.sasa.SASA(*, mode='total', probe_radius=0.14, n_sphere_points=960, **kwargs)[source]

Bases: Analysis

Solvent-accessible surface area.

Parameters:

mode ({"total", "residue", "average_residue"}, default "total") – "total" returns one value per frame (sum over all atoms). "residue" returns a per-residue SASA matrix (n_frames × n_residues).
probe_radius (float, default 0.14) – Probe (solvent) radius in nm. The default 0.14 nm is the water radius and is the standard choice for biomolecular SASA.
n_sphere_points (int, default 960) – Number of points on the unit sphere for the Shrake-Rupley rolling ball. Higher is more accurate but slower. 960 is MDTraj’s default and provides ~1% precision.
**kwargs – Standard base-class options.
Output
------
``frame (sasa.dat — CSV. Either) –
or (sasa_nm2`` (total))
``frame –
residue (sasa_nm2`` (per)
residue
format). (long)
(residue). (sasa.png — Time series (total) or heatmap)

compute(traj)[source]

Compute SASA per frame.

Returns:: mode="total": columns frame, sasa_nm2. mode="residue": columns frame, residue, sasa_nm2 (long form).
Return type:: pandas.DataFrame
Parameters:: traj (Trajectory)

default_selection: str | None = 'protein': A surface accessible to solvent, computed with the solvent present, is occluded by the very water whose access it measures: the number is not large and slow to reach, it is wrong. A run through the orchestrator never met this because the scope selection resolves to the solute, but a direct call from a notebook did, and paid for it twice – on a small test with 1,500 waters the whole system took 5.9 times as long as the protein alone, and a solvated box carries ten times that many.

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Solvent-accessible surface area': Human-readable description used in figure titles.

frame_axis_for_plot(traj, n_points)[source]

Parameters:

traj (Trajectory | None)
n_points (int)

Return type:

mode: str

n_sphere_points: int

name: str = 'sasa': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (DataFrame)
ax (Axes)

Return type:

None

probe_radius: float

reweightable: tuple[str | None, str] | None = ('sasa_nm2', 'SASA (nm²)')

The per-frame scalar whose ensemble average means something, as (column, label). column names a column when compute() returns a DataFrame and is None when it returns a bare per-frame array. None for the attribute itself – the default – means no weighted average of this analysis is offered.

On a metadynamics run the trajectory is not a Boltzmann ensemble, so every mean reported from it is an average over a distribution the bias flattened on purpose. Declaring this lets that mean be recomputed against the deposited bias and reported beside the raw one. Analyses whose result is not one number per frame leave it unset: reweighting a clustering or a projection is a harder question than a weighted mean, and the report says so rather than guessing.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)

save_data(result, path)[source]

Write the table, and for a per-residue run the average beside it.

That run already contains every number the average needs, so computing the surface a second time to get it would cost minutes for arithmetic. Version 1 wrote all three outputs from one run; this writes two, and the third mode exists for anyone who wants only the summary.

Parameters:: result (DataFrame)
Return type:: Any

time_series: bool = True

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

fastmdxplora.analysis.sasa.UNWRITTEN_MARGIN = 0.4: How far above the run’s usual proportion of zeros a frame must sit before it is taken as unwritten. The fault puts most of a row at zero at once, against a median of none; a real protein sits at a steady tenth or so, because the residues that are buried stay buried.

fastmdxplora.analysis.sasa.VALID_MODES = ('total', 'residue', 'average_residue'): What a SASA run reports. total is the whole molecule per frame, residue every residue per frame, and average_residue each residue’s mean over the run – which is the summary somebody reads to find out what is buried.

Secondary structure assignment.

Per-residue secondary structure across the trajectory using DSSP (Kabsch & Sander algorithm via MDTraj). Produces two outputs:

A time-series heatmap showing the secondary structure of each residue at each frame (residue × frame matrix, colored by DSSP code).

The DSSP codes as a CSV (one row per frame, columns are residues).

DSSP codes used (MDTraj’s “simplified” 3-state output by default):

H : helix (3-10, alpha, pi)
E : strand / extended (beta-sheet)
C : coil (everything else)

The classic “ribbon-plot timeline” emerging from this is one of the most informative single figures in MD trajectory analysis — it shows fold stability, secondary-structure transitions, and termini fraying at a glance.

References

Kabsch, W.; Sander, C. Biopolymers 1983, 22, 2577.

class fastmdxplora.analysis.ss.SS(*, simplified=True, **kwargs)[source]

Bases: Analysis

Per-residue secondary structure via DSSP.

Parameters:

simplified (bool, default True) – If True, use MDTraj’s three-letter simplification (H/E/C). If False, use the full eight-letter DSSP alphabet (H/E/B/G/I/T/S/C) which is then folded down to three classes for the figure but preserved in the saved data.
**kwargs – Standard base-class options.
Output
------
matrix. (ss.dat — CSV with the per-frame DSSP code)
frame) (ss.png — Heatmap (residue ×)
class. (colored by structure)

compute(traj)[source]

Run DSSP per frame.

Returns:: Shape (n_frames, n_residues). Cell values are single-letter DSSP codes. Column names are residue resSeq numbers (PDB numbering when available, else topology indices).
Return type:: pandas.DataFrame
Parameters:: traj (Trajectory)

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Secondary structure (DSSP)': Human-readable description used in figure titles.

name: str = 'ss': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (DataFrame)
ax (Axes)

Return type:

None

save_data(result, path)[source]

Write the computed result to path.

Default behaviour:

1-D / 2-D numpy arrays: np.savetxt (whitespace-delimited).
pandas DataFrames: to_csv (comma-delimited).
Anything else: subclass must override.

Returns the path actually written.

Parameters:: result (DataFrame)
Return type:: Any

Backbone dihedrals and the Ramachandran plot.

Computes the backbone phi (C-N-Cα-C), psi (N-Cα-C-N) and omega (Cα-C-N-Cα) dihedral angles for every protein residue across the trajectory, and produces a Ramachandran plot — the joint distribution of phi/psi pairs coloured by frequency.

Omega is the peptide bond itself, near 180 degrees in almost every residue; the exceptions are the finding, a cis bond most often before a proline. Which angles are measured is a setting.

Output dihedrals.dat is a CSV with one row per (frame, residue) combination plus columns for phi and psi (degrees). Output figure is the 2-D Ramachandran scatter/density plot.

References

The phi/psi assignment follows the IUPAC convention. MDTraj’s mdtraj.compute_phi() / mdtraj.compute_psi() are used; both return angles in radians, which we convert to degrees for the standard Ramachandran display range of (-180°, 180°).

class fastmdxplora.analysis.dihedrals.Dihedrals(*, density=True, angles=('phi', 'psi', 'omega'), bins=72, **kwargs)[source]

Bases: Analysis

Backbone phi/psi dihedrals and the Ramachandran plot.

Parameters:

density (bool, default True) – If True, render the Ramachandran plot as a 2-D histogram (heatmap) showing point density. If False, render as a scatter plot. Density is more readable for long trajectories with many points; scatter is better for short trajectories or when you want to see each sample.
angles (sequence of {"phi", "psi", "omega"}, default all three) – Which backbone torsions to measure. Phi and psi are the Ramachandran pair. Omega is the peptide bond itself, close to 180 degrees in almost every residue – and the exceptions are the finding: a cis peptide bond near zero, most often before a proline, and the departures from planarity that a strained fold produces.
bins (int, default 72) – Number of bins along each axis for the density plot. The Ramachandran range is 360°, so the default 72 bins → 5° resolution.
**kwargs – Standard base-class options. selection is ignored here because the dihedrals are inherently a backbone property; MDTraj’s functions handle the residue iteration internally.
Output
------
``frame (dihedrals.dat — CSV with columns) –
residue
phi_deg
psi_deg``.
default). (dihedrals.png — Ramachandran plot (density heatmap by)

angles: tuple[str, ...]

bins: int

compute(traj)[source]

Compute backbone phi/psi for every (frame, residue) pair.

Returns:: Long format with columns: frame, residue, phi_deg, psi_deg. Residues where either phi or psi cannot be computed (first/last residues, chain breaks) are dropped from the table.
Return type:: pandas.DataFrame
Parameters:: traj (Trajectory)

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

density: bool

description: str = 'Backbone dihedrals (Ramachandran)': Human-readable description used in figure titles.

honours_selection: bool = False: This works out its own atoms, so a general selection has nothing to apply to.

name: str = 'dihedrals': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (DataFrame)
ax (Axes)

Return type:

None

save_data(result, path)[source]

CSV with frame, residue, phi_deg, psi_deg columns.

Parameters:: result (DataFrame)
Return type:: Any

fastmdxplora.analysis.dihedrals.VALID_ANGLES = ('phi', 'psi', 'omega'): Which backbone torsions can be measured. Named here so a form can offer them: the same constant the analysis validates against.

Q-value: fraction of native contacts retained per frame.

The Q-value (or Q-fraction) is the canonical folding-state metric in protein dynamics. For each frame, it reports the fraction of the reference structure’s residue-residue contacts that are still present. Q ≈ 1 means the fold is intact; Q ≈ 0 means it is fully unfolded.

Both the contacts and the measure follow Best, Hummer, & Eaton (PNAS 2013). The set S is over pairs of heavy atoms whose residues lie at least min_seq_separation apart in sequence (default 4, which excludes covalent and local contacts that do not probe the tertiary fold) and which are within a cutoff (default 0.45 nm) in the reference frame.

That S is over atom pairs, not residue pairs, is the whole of the definition and not a detail. Counting one closest-heavy distance per residue pair is a different measure: it gives every contacting pair of residues the same weight whether they touch at one atom or at eight, and it produces a different |S|, a different denominator, and different numbers. Both are defensible measures; only one is the published one, and Q is quoted across papers as though it were a single quantity. Measured against the reference implementation shipped with MDTraj on a peeling hairpin, the residue-pair reading of this formula stood 0.263 away from the published one at its worst frame, on a scale that runs zero to one. The scheme option selects between them and defaults to the paper.

Each native contact then contributes through a switching function rather than a threshold:

Q(X) = 1/|S| * sum over (i,j) in S of
       1 / (1 + exp[beta * (r_ij(X) - lambda * r0_ij)])

with beta = 50 per nm and lambda = 1.8. What that buys is a contact judged against the distance it had natively, not against one number shared by every pair: a contact formed at 0.20 nm and one formed at 0.44 nm are each allowed to stretch by the same factor before they stop counting, and each stops counting gradually rather than at a step. A single threshold applied to every pair reports a fold coming apart while the paper’s measure still reads it as intact.

References

Best, R. B.; Hummer, G.; Eaton, W. A. Native contacts determine protein folding mechanisms in atomistic simulations. PNAS 2013, 110, 17874.

class fastmdxplora.analysis.qvalue.QValue(*, ref=0, cutoff=0.45, beta=50.0, lambda_factor=1.8, min_seq_separation=4, scheme='heavy-atom-pairs', **kwargs)[source]

Bases: Analysis

Fraction of native contacts retained per frame.

Parameters:

ref (int, default 0) – Reference frame defining the “native” state. The contacts present in this frame become the denominator of the Q calculation.
scheme ({"heavy-atom-pairs", "residue-closest-heavy"}, default) – “heavy-atom-pairs” Which set S to build. heavy-atom-pairs is the published definition: every pair of heavy atoms within cutoff in the reference, subject to the sequence separation. residue-closest-heavy uses one closest-heavy distance per residue pair instead, which is a coarser measure that is not comparable with published Q values.
cutoff (float, default 0.45) – Heavy-atom-contact cutoff in nm, used to decide which pairs are natively in contact in the reference frame. It does not decide whether a contact is present later: that is what the switching function is for.
beta (float, default 50.0) – Steepness of the switching function, in 1/nm. The paper’s value. Larger values approach a hard threshold.
lambda_factor (float, default 1.8) – How far a contact may stretch relative to its own native distance before it stops counting. The paper’s value.
min_seq_separation (int, default 4) – Minimum |i - j| in sequence for a pair to be considered. The default 4 excludes local contacts (i±1, i±2, i±3) that don’t probe the global fold.
selection (str, default "protein") –
Which atoms may form contacts. Hydrogens are excluded from S whatever is selected, because the published definition is over heavy atoms, so the useful choices narrow it further:
- "protein" (the default) gives the published measure, every heavy atom of the chain.
- "backbone" reports on the fold’s topology alone and is insensitive to side-chain repacking, which is what makes it the choice when the question is whether the chain has changed course rather than whether a core has loosened.
- "name CA" gives the coarse-grained measure familiar from Go-model work. It is a different quantity from the other two and is not comparable with a Q quoted from an all-atom study.
The choice moves the number as much as the scheme does, so it belongs in the record with it, and it is written to options.json for that reason.
**kwargs – Standard base-class options.
Output
------
[0 (qvalue.dat — single-column file of Q per frame (range)
1]).
frame/time. (qvalue.png — Time series of Q vs.)

SCHEMES = ('heavy-atom-pairs', 'residue-closest-heavy'): The two readings of S, and what each one counts.

compute(traj)[source]

Compute Q per frame.

Returns:: Q values in [0, 1]. NaN if the reference has zero contacts (the calculation is undefined — e.g. an unfolded reference).
Return type:: np.ndarray of shape (n_frames,)
Parameters:: traj (Trajectory)

default_selection: str | None = 'protein': Q is a statement about a fold, and solvent has none. Left at the whole trajectory, S came out of a solvated system as 98% water-water pairs on a small test case, and the number reported would have been how much of the water’s starting arrangement survived. A full run passes a scope selection and never saw this; a direct call from a notebook is the ordinary way to meet it.

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Native contact fraction (Q)': Human-readable description used in figure titles.

frame_axis_for_plot(traj, n_points)[source]

Parameters:

traj (Trajectory | None)
n_points (int)

Return type:

name: str = 'qvalue': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

requires_tertiary_structure = True

The shortest chain this can say anything about.

Q measures how much of a fold is intact, from contacts between residues far enough apart in sequence to probe tertiary structure. A chain shorter than min_seq_separation + 1 residues has no such pair, and no fold for them to describe. Raising there recorded a failed analysis for a peptide that simply has no tertiary structure – the same category error requires_water exists to avoid, where “there is no water here” was a failed phase rather than a question that did not apply.

reweightable: tuple[str | None, str] | None = (None, 'Fraction of native contacts')

The per-frame scalar whose ensemble average means something, as (column, label). column names a column when compute() returns a DataFrame and is None when it returns a bare per-frame array. None for the attribute itself – the default – means no weighted average of this analysis is offered.

On a metadynamics run the trajectory is not a Boltzmann ensemble, so every mean reported from it is an average over a distribution the bias flattened on purpose. Declaring this lets that mean be recomputed against the deposited bias and reported beside the raw one. Analyses whose result is not one number per frame leave it unset: reweighting a clustering or a projection is a harder question than a weighted mean, and the report says so rather than guessing.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)

time_series: bool = True

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

fastmdxplora.analysis.qvalue.native_contact_pairs(traj, *, ref=0, cutoff=0.45, min_seq_separation=4, atom_indices=None)[source]

The published S: heavy-atom pairs in contact in the reference frame.

Returns the pairs in traj’s own atom indexing, together with the distance each pair had in the reference.

This is module-level rather than a method because two callers need the same answer: the analysis that reports Q, and the collective variable that biases it. A CV built from a separately written contact set would bias one quantity and report another, and the disagreement would look like a sampling problem rather than a definition problem.

Built from a neighbour search rather than from every pair of heavy atoms. The reference implementation enumerates all combinations and then discards those beyond the cutoff, which on a real protein means holding millions of pairs to keep some thousands; only pairs already within the cutoff can be native contacts, so the search asks for those directly and the enumeration never happens.

Parameters:

traj (Trajectory)
ref (int)
cutoff (float)
min_seq_separation (int)
atom_indices (Any)

Return type:

tuple[ndarray, ndarray]

Backbone N–H generalised order parameters, S^2.

The quantity NMR relaxation measures for each amide in a protein: how much of a bond vector’s orientation survives the fast internal motion, on a scale where one is rigid and zero is freely reorienting. It is the most direct comparison a trajectory has against a solution measurement, because both describe the same picosecond-to-nanosecond reorientation of the same bond, and unlike crystallographic B-factors it is not damped by a lattice.

S^2 is the plateau of the internal second-rank correlation function,

C(t) = <P2(u(0) . u(t))>, P2(x) = (3x^2 - 1) / 2

with the global tumbling removed. Computed here in the closed form that plateau takes once the internal motion has decorrelated,

S^2 = 3/2 * sum_ab <u_a u_b>^2 - 1/2

over the Cartesian components of the unit bond vector. cite{lipari1982} The two routes agree to numerical precision on a trajectory long enough for the correlation function to have a plateau at all, and the test suite checks them against each other rather than trusting either alone.

Global tumbling has to go first, and how it goes changes the answer. The measurement is of motion relative to the molecule, so the trajectory is superposed before the vectors are taken. The atoms used for that superposition are a choice, recorded with the result: aligning on a flexible terminus drags the frame around with it and depresses S^2 everywhere else.

A short trajectory reports S^2 too high, not too noisy. Motion slower than the run is motion the run never saw, and unsampled motion looks like rigidity. That failure is one-sided, so it cannot be spotted by looking at the scatter. It is checked here by computing the same order parameters from each half of the trajectory: if the halves disagree, the number is reported with a statement that it is an upper bound rather than presented as a measurement.

References

Lipari, G.; Szabo, A. Model-free approach to the interpretation of nuclear magnetic resonance relaxation in macromolecules. 1. J. Am. Chem. Soc. 1982, 104, 4546.

fastmdxplora.analysis.order_parameters.AMIDE_H_NAMES = ('H', 'HN'): Amide hydrogen names in the force fields this software builds systems with. H is the AMBER and PDB v3 name; HN is CHARMM’s; H1 appears on the first residue, where it is one of three and not an amide proton at all, which is why the N-terminus is excluded below.

fastmdxplora.analysis.order_parameters.HALVES_TOLERANCE = 0.02: How far the two halves of a trajectory may disagree before the order parameters are called an upper bound rather than a measurement. Chosen against what the comparison is for: published S^2 sets are quoted to about 0.02, so a disagreement larger than that is larger than the difference anyone would be trying to detect.

class fastmdxplora.analysis.order_parameters.OrderParameters(*, align_selection='name CA', ref=0, **kwargs)[source]

Bases: Analysis

Backbone N–H order parameters against residue number.

Parameters:

align_selection (str, default "name CA") – Atoms used to remove global tumbling. Recorded with the result, because it is a choice that changes the answer.
ref (int, default 0) – Frame the superposition is made onto.
**kwargs – Standard base-class options.
Output
------
columns (order_parameters.dat -- two)
S^2. (residue number and)
number. (order_parameters.png -- S^2 against residue)

compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = 'protein': Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Backbone N--H order parameters (S^2)': Human-readable description used in figure titles.

min_atoms_to_align: int = 3: Superposition needs a frame to define, and three atoms is the fewest that define one. Declared rather than only checked, because the orchestrator reads it at plan time and leaves the analysis out of a molecule too small to align, instead of running it to failure.

name: str = 'order_parameters': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

requires_amide_hydrogens = True: A system without amide hydrogens poses no question here, rather than posing one this fails to answer. See the water gate in the orchestrator for the category this belongs to.

time_series: bool = False

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

fastmdxplora.analysis.order_parameters.amide_pairs(topology, atom_indices=None)[source]

Backbone N and its amide hydrogen, per residue that has one.

Returns (nitrogen index, hydrogen index, residue index).

Proline is absent by chemistry rather than by convention: its backbone nitrogen is in the ring and carries no hydrogen, so there is no vector to measure and no experimental value to compare against. The first residue is absent for a different reason, that its nitrogen carries a charged amino group whose hydrogens are not the amide proton the measurement is about.

Parameters:

topology (Topology)
atom_indices (Any)

Return type:

list[tuple[int, int, int]]

fastmdxplora.analysis.order_parameters.correlation_plateau(vectors, *, lag_fraction=0.5)[source]

S^2 the long way, as the tail of C(t) = <P2(u(0).u(t))>.

Kept because it is the definition, and the closed form above is an identity that holds only once the internal motion has decorrelated. Where the two disagree, the trajectory is too short for the plateau to exist, which is worth knowing and is what the test suite uses this for.

Parameters:

vectors (ndarray)
lag_fraction (float)

Return type:

fastmdxplora.analysis.order_parameters.order_parameters(vectors)[source]

S^2 from unit bond vectors of shape (frames, bonds, 3).

The closed form of the correlation-function plateau. Summing the nine products of components rather than fitting a decay avoids choosing where a plateau begins, which on a noisy correlation function is a choice that moves the answer.

Parameters:: vectors (ndarray)
Return type:: ndarray

Simulated fluctuations against crystallographic B-factors.

The oldest comparison between a trajectory and an experiment, and the one most often quoted without its caveats. A refined B-factor and a simulated RMSF are related through

B = (8 pi^2 / 3) <u^2>,

so a per-residue B converts to the fluctuation amplitude a crystal implies, in the same units the trajectory reports. What comes out is a correlation, and this analysis reports it with the reasons it is not an accuracy.

A crystallographic B is not only motion. It absorbs static disorder across the molecules in the lattice, the refinement’s own restraints, and whatever TLS or occupancy model was used; and the lattice itself damps the loop excursions a solution trajectory is free to make. The result is that B-factors bound loop amplitudes from below, so a simulation that agrees everywhere may be reporting a protein held too tightly, and one that exceeds them in loops may be right. The correlation is worth having and a regression slope is not, which is why only the first is reported.

mdtraj does not carry B-factors, so they are read from the deposited file rather than from the trajectory’s topology. That also keeps the comparison against the structure as deposited rather than against the prepared system, whose B column is whatever the preparation left there.

class fastmdxplora.analysis.bfactor_comparison.BFactorComparison(*, structure=None, align_selection='name CA', **kwargs)[source]

Bases: Analysis

Per-residue RMSF against the fluctuations a crystal implies.

Parameters:

structure (str, optional) – The deposited file to read B-factors from. Discovered from the run directory when not given.
align_selection (str, default "name CA") – Atoms used to remove rigid-body motion before fluctuations are measured, as in the RMSF analysis.
**kwargs – Standard base-class options.
Output
------
number (bfactor_comparison.dat -- residue)
nm (simulated RMSF in)
and
B-factor. (the RMSF implied by the deposited)

compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = 'protein': Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Simulated RMSF against crystallographic B-factors': Human-readable description used in figure titles.

min_atoms_to_align: int = 3

three atoms is the fewest that define a frame, and declaring it lets the orchestrator leave this out of a molecule too small rather than run it to failure.

Type:: As for every analysis that superposes

name: str = 'bfactor_comparison': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

requires_crystallographic_bfactors = True: Only where a deposited structure with real B-factors was found. A run built from a generated or minimised coordinate file poses no question here rather than failing to answer one.

time_series: bool = False

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

fastmdxplora.analysis.bfactor_comparison.B_TO_MSF = 0.037995443865876666: B = (8 pi^2 / 3) <u^2>, so <u^2> = 3B / (8 pi^2), in the file’s units of square Angstroms.

fastmdxplora.analysis.bfactor_comparison.bfactors_from_pdb(path)[source]

Per-residue CA B-factors, keyed by (chain order, residue number).

Read by column rather than by splitting on whitespace, because the PDB format is fixed-width and a five-figure atom serial or a four-character residue name runs its fields together: split on spaces and the B column becomes whatever happened to be next to it.

Chains are keyed by the order their identifiers first appear, which is what survives preparation: PDBFixer keeps chain order while it may not keep the identifiers themselves.

Parameters:: path (str | Path)
Return type:: dict[tuple[int, int], float]

fastmdxplora.analysis.bfactor_comparison.has_crystallographic_bfactors(path)[source]

Whether a file carries B-factors that mean anything.

A minimised or generated structure writes zeros, and a comparison against a column of zeros is not a comparison. One non-zero value is not enough either: a file where a handful of atoms carry a placeholder would pass, so this asks that most of them do.

Parameters:: path (str | Path)
Return type:: bool

Density, energy and temperature, from the record the run already kept.

The simulation phase writes a state table every few hundred steps: step, time, potential and kinetic energy, temperature, volume and density. Nothing read it back. It is the only place the ensemble itself is visible, and the quantities in it are the ones with published values to check against: the density of a water model is a number with one right answer, and a mean temperature that misses the thermostat’s setpoint says the run was not doing what the configuration said.

Each column is treated as the correlated time series it is, with the same settling and effective-sample machinery every other observable here gets, so a density arrives with an error that reflects how many independent observations stand behind it rather than how many lines were written.

Density means nothing at constant volume. Under NVT the box does not move, so the density is a constant, its variance is zero, and a mean with an error bar on it would be a statement about arithmetic rather than about the system. The ensemble is read from the record rather than assumed, and where the volume never changed this says so instead of quoting a spread of zero as a precise measurement.

fastmdxplora.analysis.thermodynamics.COLUMNS: dict[str, tuple[str, str]] = {'density': ('Density', 'g/mL'), 'kinetic_energy': ('Kinetic Energy', 'kJ/mol'), 'potential_energy': ('Potential Energy', 'kJ/mol'), 'temperature': ('Temperature', 'K'), 'total_energy': ('Total Energy', 'kJ/mol'), 'volume': ('Box Volume', 'nm^3')}: The columns worth reporting, and what each is called in the record OpenMM writes. Matched on a substring because the header carries units (“Density (g/mL)”) and the units are part of the name rather than something to parse off it.

class fastmdxplora.analysis.thermodynamics.Thermodynamics(*, state_csv=None, **kwargs)[source]

Bases: Analysis

Ensemble observables from the simulation’s own state record.

Parameters:

state_csv (str, optional) – The state record to read. Discovered beside the run when not given, which is where the simulation phase writes it.
**kwargs – Standard base-class options.
Output
------
observable (thermodynamics.dat -- one row per)
mean (with its settled)

:param : :param the error on it: :param and the number of independent observations behind it.:

compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

description: str = 'Density, energy and temperature from the state record': Human-readable description used in figure titles.

honours_selection: bool = False: The state record is of the whole box, so an atom selection would describe a different quantity than the one written.

name: str = 'thermodynamics': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

requires_state_record = True: Only where the simulation phase left a state record.

save_data(result, path)[source]

Write the computed result to path.

Default behaviour:

1-D / 2-D numpy arrays: np.savetxt (whitespace-delimited).
pandas DataFrames: to_csv (comma-delimited).
Anything else: subclass must override.

Returns the path actually written.

Parameters:

result (ndarray)
path (Path)

Return type:

time_series: bool = False

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

fastmdxplora.analysis.thermodynamics.VOLUME_VARIES_ABOVE = 1e-06: How much the box must vary before the run counts as constant-pressure. A relative standard deviation; anything below this is a fixed box with floating-point noise on it.

fastmdxplora.analysis.thermodynamics.read_state_table(path)[source]

The state record, by column, with the units left in the header.

Read with the csv module rather than by splitting on commas: OpenMM’s header quotes its field names, and several of them contain a comma inside the quotes.

Parameters:: path (str | Path)
Return type:: dict[str, ndarray]

The radial distribution function between two selections.

How the density of one group varies with distance from the other, relative to what it would be if the two ignored each other. A value of one means no structure at that separation; the first peak is the first solvation shell, and where it sits and how tall it stands are quantities with published values for every common water model, which makes this one of the few things a simulation reports that can be checked against a number somebody else measured.

A radial distribution function needs a box. The normalisation divides by the bulk density, and the bulk density is the pair count divided by the volume the system occupies. A trajectory carrying no unit cell has no such volume, and the curve that comes out of assuming one is a histogram wearing the units of a g(r).

And it stops at half the box. Beyond that separation the minimum-image convention no longer supplies a complete shell: part of every sphere lies outside the periodic cell and is counted from the wrong image or not at all, so g(r) sags towards zero for reasons that have nothing to do with the liquid. The curve stays smooth and plausible while it does this, which is why the range is capped here rather than left to the person reading the plot.

fastmdxplora.analysis.rdf.MAX_PAIRS = 400000: Pairs beyond this many are subsampled before the histogram. A protein against every water oxygen is tens of millions of pairs per frame, which is minutes of arithmetic for a curve that a random tenth of them resolves to within the width of its own line.

class fastmdxplora.analysis.rdf.RadialDistribution(*, selection_a='protein', selection_b='water and name O', r_max=None, bin_width=0.005, **kwargs)[source]

Bases: Analysis

g(r) between two atom selections.

Parameters:

selection_a (str, default "protein") – The first group. Distances are measured from these atoms.
selection_b (str, default "water and name O") – The second group. Water oxygens by default, which is the pairing a solvation shell is usually read from.
r_max (float, optional) – Where to stop, in nm. Capped at half the smallest box dimension, and defaults to it.
bin_width (float, default 0.005) – Histogram bin width in nm.
**kwargs – Standard base-class options.
Output
------
columns (rdf.dat -- two)
g(r). (separation in nm and)

compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Radial distribution function between two selections': Human-readable description used in figure titles.

honours_selection: bool = False: The two selections are the analysis’s own, and a scope selection would name a third thing that is neither of them.

name: str = 'rdf': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

requires_periodic_box = True: Without a unit cell there is no bulk density to normalise against.

requires_water = True

The default pair is the solute against water oxygens, and a default run saves the solute alone (save_selection), so this analysis could not succeed on a default configuration: two defaults asking for different things. water_sites already declares this and rdf did not, which is the whole of the defect – the gate existed, this analysis was simply not behind it.

The gate is coarser than the question. Whether this analysis needs water depends on the selections it is given, and the plan is built before options are merged, so a class attribute is all that can be declared. The cost falls on one case: a solvent-free trajectory with custom selections naming neither group as water is now skipped by default rather than run. That case needs include: [rdf], which is honoured as written. The case it fixes is every default run.

time_series: bool = False

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

Conformational clustering.

Clusters frames of the trajectory by structural similarity, producing one cluster-membership labeling per requested method. Three methods are supported:

k-means (default, fast, requires choosing n_clusters)

hierarchical (agglomerative, Ward linkage by default)

dbscan (density-based, no pre-specified cluster count)

All methods operate on the pairwise RMSD distance matrix (computed via MDTraj’s QCP algorithm), which is the standard featurization for conformational clustering. Outputs per method:

cluster_<method>.dat — per-frame integer cluster labels.

cluster_<method>.png — cluster labels as a function of time.

cluster_<method>_counts.png — cluster population bar chart.

cluster_hierarchical_dendrogram.png — hierarchical dendrogram when hierarchical clustering is requested and SciPy is available.

hierarchical_distance_matrix.npy and hierarchical_linkage.npy — reproducibility data for dashboard/report-native dendrogram rendering.

Because this analysis produces multiple files per run, it overrides the base class’s save_data() and _do_plot() methods.

References

Daura, X. et al. J. Mol. Graph. Model. 1999, 18, 122 (RMSD-based MD clustering). Lloyd, S. IEEE Trans. Inf. Theory 1982, 28, 129 (k-means).

class fastmdxplora.analysis.cluster.Cluster(*, methods=('kmeans', 'hierarchical'), features='rmsd', n_clusters=5, eps=0.2, min_samples=5, linkage='average', random_state=42, n_init=10, **kwargs)[source]

Bases: Analysis

Conformational clustering by pairwise RMSD.

Parameters:

methods (sequence of str, default ("kmeans", "hierarchical")) – Which clustering algorithms to run. Each produces its own output files. Valid values: "kmeans", "hierarchical", "dbscan". Two are run by default because they disagree in useful ways: k-means insists every frame joins a cluster, while hierarchical linkage shows how the clusters nest, and a conformational split that both find is worth more than one only either sees.
n_clusters (int, default 5) – Number of clusters (used by k-means and hierarchical).
eps (float, default 0.2) – DBSCAN distance threshold in nm.
min_samples (int, default 5) – DBSCAN minimum samples per cluster.
features ({"rmsd", "coordinates"}, default "rmsd") – What the frames are compared in. "rmsd" measures every pair with its own optimal superposition, so the distance does not depend on how the molecule happens to be placed. "coordinates" superposes every frame onto the first and compares coordinates directly, scaled so a distance is still an RMSD in nm – the cheaper approximation, exact only where one alignment serves every pair. The choice changes the answer more than any parameter here does, so it is worth stating.
random_state (int, default 42) – The seed k-means starts from. K-means finds a local optimum, so a different seed can find a different clustering: one that survives a change of seed is a finding, and one that does not is an artefact of where the algorithm happened to start. Fixed by default so a run repeats, and settable so that can be tested – it was written into the code, which made every run agree and hid the question.
n_init (int, default 10) – How many starts k-means makes, keeping the best. More costs time and buys robustness against a poor start.
linkage ({"ward", "complete", "average", "single"}, default "average") – Hierarchical linkage method. "ward" needs the frames as points rather than as distances: in the coordinate space they are, and from an RMSD matrix a classical MDS embedding stands in for them.
selection (str, optional) – MDTraj atom selection for the RMSD calculation. Defaults to "name CA" (CA-only is fast and capture the global fold well).
**kwargs – Standard base-class options.
Output
------
method (Per) –
- cluster_<method>.dat — CSV with frame, cluster columns.
- cluster_<method>.png — Cluster timeline figure.
<output_dir>/cluster/ (in) –
- cluster_<method>.dat — CSV with frame, cluster columns.
- cluster_<method>.png — Cluster timeline figure.

compute(traj)[source]

Run all requested clustering methods.

Returns:: Maps method name → array of per-frame cluster labels (int). DBSCAN’s -1 label indicates “noise” (unclustered frames).
Return type:: dict
Parameters:: traj (Trajectory)

default_selection: str | None = 'name CA': Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

description: str = 'Conformational clustering': Human-readable description used in figure titles.

eps: float

features: str

linkage: str

methods: list[str]

min_atoms_to_align: int = 3: A superposition needs three atoms to be defined. Without this, MDTraj returns identity rotations and the frames are compared unaligned – a real run clustered a capped alanine and found one distinct cluster where five were asked for, and reported ok.

min_samples: int

n_clusters: int

n_init: int

name: str = 'cluster': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (dict[str, ndarray])
ax (Axes)

Return type:

None

random_state: int

reweightable_populations: bool = True: Populations are weighted counts, so a biased run’s cluster occupancies are recoverable. Which clusters exist is not: see the base class.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)
Return type:: AnalysisResult

fastmdxplora.analysis.cluster.VALID_FEATURES = ('rmsd', 'coordinates')

What the frames are compared in.

rmsd measures every pair with its own optimal superposition, so the distance is invariant to how the molecule happens to be placed. ward needs a Euclidean space and reaches one through a classical MDS embedding.

coordinates superposes every frame onto the first and compares the coordinates directly, which is the cheaper approximation and the space ward and k-means were defined in. The coordinates are scaled by 1/sqrt(n_atoms), so a distance in this space is an RMSD measured against that one common alignment rather than a pairwise one: eps keeps its meaning in nm, and the two feature spaces stay comparable.

Which to use is a real choice and not a detail. A common alignment is exact only where the optimal pairwise superposition is the same one; across a large conformational change it is not, and the two will disagree.

Dimensionality reduction.

Projects the high-dimensional configuration space of an MD trajectory down to two or three dimensions for visualization. Four methods – PCA, MDS, t-SNE and UMAP – which answer different questions: PCA finds the directions of largest variance in the coordinates, while MDS finds an arrangement preserving the distances between frames, and that distance is RMSD.

PCA (default): principal component analysis on the aligned Cartesian coordinates. Linear, fast, decomposes the variance into orthogonal collective modes. Standard in MD analysis.

t-SNE: t-distributed stochastic neighbor embedding. Non-linear, preserves local neighborhood structure, useful for visualizing metastable basins. Stochastic — set random_state for reproducibility.

UMAP (optional): uniform manifold approximation and projection. Non-linear, generally faster than t-SNE, also preserves global structure better. Requires the optional umap-learn package.

Each method produces one 2-D scatter dimred_<method>.png colored by frame index (the “trajectory trace” visualization), plus a data file dimred_<method>.dat with the projected coordinates.

References

Amadei, A.; Linssen, A. B. M.; Berendsen, H. J. C. Proteins 1993, 17, 412 (PCA). van der Maaten, L.; Hinton, G. J. Mach. Learn. Res. 2008, 9, 2579 (t-SNE).

class fastmdxplora.analysis.dimred.DimRed(*, methods=('pca',), n_components=2, perplexity=30.0, n_neighbors=15, min_dist=0.1, random_state=42, **kwargs)[source]

Bases: Analysis

Dimensionality reduction on the trajectory.

Parameters:

methods (sequence of str, default ("pca",)) – Which methods to run. Choices: "pca", "tsne", "umap".
n_components (int, default 2) – Dimensionality of the embedding. For visualization keep at 2 (or 3).
perplexity (float, default 30.0) – t-SNE perplexity parameter. Roughly the effective number of neighbors each point is balanced against; 5-50 is typical.
n_neighbors (int, default 15) – UMAP neighborhood size.
min_dist (float, default 0.1) – UMAP minimum distance between embedded points.
random_state (int, default 42) – Random seed for stochastic methods (t-SNE, UMAP).
selection (str, optional) – MDTraj atom selection used to flatten coordinates. Defaults to "name CA" (CA-only is a standard featurization for protein DimRed).
**kwargs – Standard base-class options.
Output
------
method (Per) –
- dimred_<method>.dat — CSV with frame + component columns.
- dimred_<method>.png — 2-D scatter colored by frame index.
<output_dir>/dimred/ (in) –
- dimred_<method>.dat — CSV with frame + component columns.
- dimred_<method>.png — 2-D scatter colored by frame index.

compute(traj)[source]

Run all requested DimRed methods.

Returns:: Maps method name → (n_frames, n_components) embedding array.
Return type:: dict
Parameters:: traj (Trajectory)

default_selection: str | None = 'name CA': Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

description: str = 'Dimensionality reduction': Human-readable description used in figure titles.

min_atoms_to_align: int = 3: A superposition needs three atoms to be defined. Without this, MDTraj returns identity rotations and the frames are compared unaligned – a real run clustered a capped alanine and found one distinct cluster where five were asked for, and reported ok.

name: str = 'dimred': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (dict[str, ndarray])
ax (Axes)

Return type:

None

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)
Return type:: AnalysisResult

Hydrogen bond analysis.

Identifies hydrogen bonds across the trajectory using either the Baker-Hubbard or Wernet-Nilsson geometric criteria, and produces two outputs: a per-frame H-bond count time series, and a long-form table of which donor/acceptor pairs participated in H-bonds, with occupancy fractions.

\*\*Baker-Hubbard** (default) — D–H···A angle > 120°, H···A distance < 2.5 Å,

applied with a 10% occupancy threshold by default. Standard choice for

protein backbone H-bonds.

\*\*Wernet-Nilsson** — Geometric criterion designed for water; the cutoff

distance is dynamically adjusted by the D–H–A angle. Useful for protein-

water and water-water bonds in solvated systems.

References

Baker, E.; Hubbard, R. Prog. Biophys. Mol. Biol. 1984, 44, 97. Wernet, P. et al. Science 2004, 304, 995.

class fastmdxplora.analysis.hbonds.HBonds(*, method='baker_hubbard', freq=0.1, candidate_freq=0.0, count_multiplier=1, periodic=True, distance_cutoff=0.25, angle_cutoff=120.0, sidechain_only=False, exclude_water=True, **kwargs)[source]

Bases: Analysis

Hydrogen bond identification and counting.

Parameters:

method ({"baker_hubbard", "wernet_nilsson"}, default "baker_hubbard") – Geometric criterion. Baker-Hubbard is the conventional choice for protein backbone; Wernet-Nilsson is better for water hydrogen bonding.
freq (float, default 0.1) – Occupancy above which a bond is called persistent. Every bond is counted in the per-frame series whatever its occupancy; this decides only how many are reported as persistent, alongside the total number found, in n_persistent_bonds. Has no effect on Wernet-Nilsson.
candidate_freq (float, default 0.0) – Occupancy threshold Baker-Hubbard applies when proposing which bonds to evaluate. Zero proposes every bond seen in any frame, which is what a per-frame count needs: a bond present in five per cent of frames is present in those frames, and a threshold applied here would drop it from all of them. Raise it only to restrict the series to bonds that persist. Has no effect on Wernet-Nilsson.
distance_cutoff (float, default 0.25) – Hydrogen-to-acceptor distance in nm, for baker_hubbard. The published Baker-Hubbard value, and settable because other work uses others: 0.35 nm between the heavy atoms is the more common convention and is a different measurement, not a looser one.
angle_cutoff (float, default 120.0) – Donor-hydrogen-acceptor angle in degrees, for baker_hubbard.
sidechain_only (bool, default False) – Count only bonds involving a side chain. A backbone hydrogen bond holds the fold together; a side-chain one is what a substitution can change, and mixing them answers neither question.
exclude_water (bool, default True) – Leave out bonds to water. A solvated trajectory has far more of those than anything else, and counting them buries the protein’s own.
periodic (bool, default True) – Measure across the periodic boundary when the trajectory carries a unit cell. A solvated trajectory is not always imaged, and a molecule split across the boundary looks far from what it is touching. Where there is no unit cell this makes no difference.
count_multiplier (int, default 1) – Multiplies the per-frame count, and is a compatibility device rather than a convention. MDTraj enumerates each hydrogen bond once, as one donor-hydrogen-acceptor triplet, so a multiplier of 2 reports twice the number of bonds present. It exists to reproduce the counts published by version 1 and warns when it is not 1. Must be at least 1. Has no effect on Wernet-Nilsson.
**kwargs – Standard base-class options.
Output
------
counts. (hbonds.dat — CSV with frame-by-frame H-bond)
frame. (hbonds.png — Time-series of H-bond count per)

Notes

The compute() method returns a pandas DataFrame: frame, n_hbonds with one row per frame.

compute(traj)[source]

Compute per-frame H-bond counts.

Returns:: Columns frame, n_hbonds. One row per frame.
Return type:: pandas.DataFrame
Parameters:: traj (Trajectory)

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Hydrogen bonds': Human-readable description used in figure titles.

frame_axis_for_plot(traj, n_points)[source]

Parameters:

traj (Trajectory | None)
n_points (int)

Return type:

name: str = 'hbonds': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (DataFrame)
ax (Axes)

Return type:

None

reweightable: tuple[str | None, str] | None = ('n_hbonds', 'Hydrogen bonds')

The per-frame scalar whose ensemble average means something, as (column, label). column names a column when compute() returns a DataFrame and is None when it returns a bare per-frame array. None for the attribute itself – the default – means no weighted average of this analysis is offered.

On a metadynamics run the trajectory is not a Boltzmann ensemble, so every mean reported from it is an average over a distribution the bias flattened on purpose. Declaring this lets that mean be recomputed against the deposited bias and reported beside the raw one. Analyses whose result is not one number per frame leave it unset: reweighting a clustering or a projection is a harder question than a weighted mean, and the report says so rather than guessing.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)

save_data(result, path)[source]

Write the computed result to path.

Default behaviour:

1-D / 2-D numpy arrays: np.savetxt (whitespace-delimited).
pandas DataFrames: to_csv (comma-delimited).
Anything else: subclass must override.

Returns the path actually written.

Parameters:: result (DataFrame)
Return type:: Any

time_series: bool = True

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

Waters that stay put, and the difference between a site and a molecule.

Most water in a simulation is bulk: it arrives, it leaves, it means nothing. But some positions are occupied throughout — a water wedged between a ligand and a backbone carbonyl, bridging a hydrogen bond that neither could make alone. Those waters are part of the binding site, and a medicinal chemist displacing one pays for it in entropy or gains from it in affinity.

Finding them means clustering water oxygen positions across the trajectory and asking which clusters are occupied often. That much is standard.

What is not standard, and matters, is the distinction this analysis makes. A cluster occupied in ninety per cent of frames can be either of two things:

one water molecule that stayed — bound, with a residence time, and displacing it means displacing that molecule
a position that many waters passed through — a structural site, where the geometry favours a water but no particular water is held

Those are different findings with different consequences, and a cluster occupancy alone cannot tell them apart. Both are reported: how often the site is occupied, and by how many distinct molecules.

Only waters near the solute are considered. Bulk water clusters beautifully and means nothing — with enough frames, every position in the box has been occupied. The cutoff is what makes the result about the protein rather than about the density of water.

class fastmdxplora.analysis.water_sites.WaterSites(*, site_selection='auto', ligand_resname=None, cutoff_nm=0.5, eps_nm=0.12, minimum_occupancy=0.5, minimum_samples=5, maximum_radius_nm=0.25, **kwargs)[source]

Bases: Analysis

Positions where a water sits for much of the trajectory.

Parameters:

site_selection (str, default "auto") – What the waters must be near. "auto" uses the ligand where there is one and the protein otherwise, because a binding-site water is the usual question; anything else is an atom selection.
ligand_resname (str, optional) – Which residue is the ligand. Supplied by the run when there is one; it is a fact about the structure rather than a setting.
cutoff_nm (float, default 0.5) – How near a water must come to be considered at all, in nm. Waters further out are bulk, and bulk clusters beautifully while meaning nothing.
eps_nm (float, default 0.12) – How close two observed positions must be to belong to the same site, in nm. Roughly the width of a water’s thermal motion within one site.
minimum_occupancy (float, default 0.5) – The fraction of frames a site must be occupied to be reported. Below this it is passing traffic rather than a site.
minimum_samples (int, default 5) – How many observed positions a cluster needs before it counts as one. Guards against a handful of coincidental positions becoming a “site”.
maximum_radius_nm (float, default 0.25) – How far a cluster may spread and still be a site, in nm. Clustering links neighbours through neighbours, so without this the whole first hydration shell of a protein chains into one “site” occupied in every frame.
kwargs (Any)

COLUMNS = ('site', 'x', 'y', 'z', 'occupancy', 'n_observations', 'n_distinct_waters', 'longest_stay_frames', 'interpretation'): The shape of the result, whether or not anything was found. An empty frame with no columns raises when sorted, which is how “no site met the threshold” became a crash rather than a finding.

absent_because = 'this trajectory holds no water. `simulation.save_selection` defaults to `not water`, which is nine tenths of the file and nothing this analysis can work without: set it to `all` for a study whose subject is the solvent. An implicit-solvent run has no water to save.': What to do about it, when a run has none. Left out of a planned set rather than failed, but a study that asked for water sites and got no directory deserves the reason: since save_selection defaults to leaving the solvent out of the trajectory, the commonest cause is not that the run had no water but that it was not saved.

compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

description: str = 'Water positions occupied through the trajectory': Human-readable description used in figure titles.

honours_selection: bool = False: It works out its own atoms from the site selection, so a general selection has nothing to apply to.

name: str = 'water_sites': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (DataFrame)
ax (Any)

Return type:

None

requires_water = True: Runs by default only where there is water. An implicit-solvent run, or a trajectory stripped of solvent to save space, has none – and a refusal there is correct but should not fail the phase.

Ligand pose RMSD (RMSD of the ligand after protein alignment).

This is the headline protein-ligand stability metric: it measures whether the ligand stays in its binding pose over the trajectory. Each frame is rigidly aligned onto a reference using the protein atoms (so protein tumbling is removed), and then the RMSD is computed on the ligand atoms of the already-aligned coordinates. A low, flat profile means the ligand holds its pose; a rising profile means it is drifting or unbinding.

This differs from the standard RMSD, which aligns and measures on the same atom set. Here alignment (protein) and measurement (ligand) use different selections, which is the correct way to ask “how much has the ligand moved relative to the protein”.

Output is a single-column ligand_rmsd.dat of RMSD values in nanometers, and a time-series figure.

class fastmdxplora.analysis.ligand_rmsd.LigandRMSD(*, ligand_resname=None, align_selection='protein and name CA', ref=0, **kwargs)[source]

Bases: Analysis

Per-frame RMSD of the ligand after aligning on the protein.

Parameters:

ligand_resname (str) – Residue name of the ligand (e.g. "LIG"). Required — this analysis only makes sense for a protein-ligand complex. The orchestrator supplies it automatically from the setup manifest.
align_selection (str, default "protein and name CA") – Atom selection used for the rigid-body alignment (the receptor frame). Cα atoms are the standard, robust choice.
ref (int, default 0) – Reference frame. Negative indices count from the end.
**kwargs – Standard base-class options.

Notes

The selection attribute is not used for the measurement here (the measured atoms are always the ligand); alignment is controlled by align_selection.

compute(traj)[source]

Compute per-frame ligand RMSD after protein alignment.

Returns:: Ligand RMSD in nanometers.
Return type:: np.ndarray of shape (n_frames,)
Parameters:: traj (Trajectory)

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Ligand pose RMSD (after protein alignment)': Human-readable description used in figure titles.

frame_axis_for_plot(result, traj)[source]

Parameters:

result (ndarray)
traj (Trajectory | None)

Return type:

honours_selection: bool = False: This works out its own atoms, so a general selection has nothing to apply to.

name: str = 'ligand_rmsd': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

requires_ligand: bool = True: True for analyses that only apply to protein-ligand complexes (e.g. ligand pose RMSD). The orchestrator runs these automatically when a ligand is present and skips them otherwise. They can still be explicitly requested via include.

reweightable: tuple[str | None, str] | None = (None, 'Ligand RMSD (nm)')

The per-frame scalar whose ensemble average means something, as (column, label). column names a column when compute() returns a DataFrame and is None when it returns a bare per-frame array. None for the attribute itself – the default – means no weighted average of this analysis is offered.

On a metadynamics run the trajectory is not a Boltzmann ensemble, so every mean reported from it is an average over a distribution the bias flattened on purpose. Declaring this lets that mean be recomputed against the deposited bias and reported beside the raw one. Analyses whose result is not one number per frame leave it unset: reweighting a clustering or a projection is a harder question than a weighted mean, and the report says so rather than guessing.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)

time_series: bool = True

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

Ligand RMSF (per-atom fluctuation of the ligand after protein alignment).

After removing protein rigid-body motion (alignment on the protein), this measures how much each ligand atom fluctuates about its mean position. It reports the ligand’s internal flexibility in the pocket: which parts of the ligand are rigid and which sample multiple positions. Complements the ligand pose RMSD (overall displacement) with a per-atom flexibility profile.

Outputs ligand_rmsf.dat with columns (atom_serial, rmsf_nm) and a bar plot.

class fastmdxplora.analysis.ligand_rmsf.LigandRMSF(*, ligand_resname=None, align_selection='protein and name CA', ref=0, **kwargs)[source]

Bases: Analysis

Per-atom RMSF of the ligand after aligning on the protein.

Parameters:

ligand_resname (str) – Ligand residue name (e.g. "LIG"). Supplied by the orchestrator.
align_selection (str, default "protein and name CA") – Atoms used for the rigid-body alignment (the receptor frame).
ref (int, default 0) – Reference frame for the alignment.
**kwargs – Standard base-class options.

compute(traj)[source]

Compute per-ligand-atom RMSF.

Returns:: Columns: atom serial, RMSF in nm.
Return type:: np.ndarray, shape (n_ligand_atoms, 2)
Parameters:: traj (Trajectory)

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Ligand RMSF (per-atom flexibility, after protein alignment)': Human-readable description used in figure titles.

honours_selection: bool = False: This works out its own atoms, so a general selection has nothing to apply to.

name: str = 'ligand_rmsf': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (ndarray)
ax (Axes)

Return type:

None

requires_ligand: bool = True: True for analyses that only apply to protein-ligand complexes (e.g. ligand pose RMSD). The orchestrator runs these automatically when a ligand is present and skips them otherwise. They can still be explicitly requested via include.

Protein-ligand contacts.

Two complementary, commonly-reported views of how the protein engages the ligand over a trajectory:

Per-frame contact count — the number of protein residues with any heavy atom within cutoff of the ligand, frame by frame. A quick stability signal (a stable pose keeps a roughly constant contact count).
Per-residue contact frequency — for each protein residue, the fraction of frames in which it contacts the ligand. This is the binding-site “interaction fingerprint”: the residues with high frequency line the pocket.

Contacts are defined at the residue level: a residue is “in contact” in a frame if any of its atoms is within cutoff nm of any ligand atom.

Outputs pl_contacts.dat (per-frame count time series) and, alongside it, pl_contacts_per_residue.csv (residue, frequency). The figure shows the per-residue frequency fingerprint (the more informative of the two).

class fastmdxplora.analysis.contacts.Contacts(*, ligand_resname=None, cutoff=0.4, protein_selection='protein', periodic=True, **kwargs)[source]

Bases: Analysis

Protein-ligand contacts: per-frame count and per-residue frequency.

Parameters:

ligand_resname (str) – Ligand residue name (e.g. "LIG"). Supplied automatically by the orchestrator from the setup manifest.
cutoff (float, default 0.4) – Contact distance cutoff in nm (0.4 nm = 4 Angstrom, a standard heavy-atom contact threshold).
periodic (bool, default True) – Measure distances across the periodic boundary when the trajectory carries a unit cell. A solvated trajectory is not always imaged, and a molecule split across the boundary then looks far from everything it is actually touching – a bound ligand can report no contacts at all. Where there is no unit cell this makes no difference. False measures plain distances regardless.
protein_selection (str, default "protein") – The other side of the interaction. Usually the whole protein, and worth changing in three cases: a complex where only one chain matters ("protein and chainid 0"), a domain rather than the whole fold ("protein and resid 40 to 180"), or a receptor that is not a protein at all – "nucleic" for a ligand bound to DNA or RNA, which is a real and common case that the default name does not suggest.
**kwargs – Standard base-class options.

compute(traj)[source]

Compute the per-frame contact count. The per-residue frequency is computed alongside and stashed for save_data/plot.

Returns:: Columns frame, n_contacts (one row per frame).
Return type:: pandas.DataFrame
Parameters:: traj (Trajectory)

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Protein-ligand contacts (count + per-residue frequency)': Human-readable description used in figure titles.

honours_selection: bool = False: This works out its own atoms, so a general selection has nothing to apply to.

name: str = 'pl_contacts': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (DataFrame)
ax (Axes)

Return type:

None

requires_ligand: bool = True: True for analyses that only apply to protein-ligand complexes (e.g. ligand pose RMSD). The orchestrator runs these automatically when a ligand is present and skips them otherwise. They can still be explicitly requested via include.

save_data(result, path)[source]

Write the computed result to path.

Default behaviour:

1-D / 2-D numpy arrays: np.savetxt (whitespace-delimited).
pandas DataFrames: to_csv (comma-delimited).
Anything else: subclass must override.

Returns the path actually written.

Parameters:: result (DataFrame)
Return type:: Any

Protein-ligand hydrogen bonds.

Counts hydrogen bonds formed specifically between the protein and the ligand, frame by frame — the directional polar interactions that anchor the ligand in the pocket. This is distinct from the general hbonds analysis (which counts all hydrogen bonds within the selection); here every reported bond has one partner in the protein and the other in the ligand.

A per-frame H-bond list is computed with Wernet-Nilsson (which returns per-frame donor-H-acceptor triplets), and each triplet is kept only if it bridges protein and ligand. Outputs pl_hbonds.dat (frame, n_hbonds).

class fastmdxplora.analysis.pl_hbonds.ProteinLigandHBonds(*, ligand_resname=None, protein_selection='protein', **kwargs)[source]

Bases: Analysis

Per-frame count of protein-ligand hydrogen bonds.

Parameters:

ligand_resname (str) – Ligand residue name (e.g. "LIG"). Supplied by the orchestrator.
protein_selection (str, default "protein") – The other side of the interaction. Usually the whole protein, and worth changing in three cases: a complex where only one chain matters ("protein and chainid 0"), a domain rather than the whole fold ("protein and resid 40 to 180"), or a receptor that is not a protein at all – "nucleic" for a ligand bound to DNA or RNA, which is a real and common case that the default name does not suggest.
**kwargs – Standard base-class options.

compute(traj)[source]

Compute per-frame protein-ligand H-bond counts.

Returns:: Columns frame, n_hbonds (one row per frame).
Return type:: pandas.DataFrame
Parameters:: traj (Trajectory)

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Protein-ligand hydrogen bonds': Human-readable description used in figure titles.

honours_selection: bool = False: This works out its own atoms, so a general selection has nothing to apply to.

name: str = 'pl_hbonds': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (DataFrame)
ax (Axes)

Return type:

None

requires_ligand: bool = True: True for analyses that only apply to protein-ligand complexes (e.g. ligand pose RMSD). The orchestrator runs these automatically when a ligand is present and skips them otherwise. They can still be explicitly requested via include.

run(traj)[source]

Compute, plot, and save in one call.

This is the orchestrator’s standard entry point. Returns an AnalysisResult regardless of success — check result.status for "ok" vs "error".

Parameters:: traj (Trajectory)

save_data(result, path)[source]

Write the computed result to path.

Default behaviour:

1-D / 2-D numpy arrays: np.savetxt (whitespace-delimited).
pandas DataFrames: to_csv (comma-delimited).
Anything else: subclass must override.

Returns the path actually written.

Parameters:: result (DataFrame)
Return type:: Any

Typed protein-ligand interactions, frame by frame.

Counting contacts says how much of the protein a ligand touches. This says what is holding it: a salt bridge a charge change would destroy, a hydrophobic packing that tolerates one, a hydrogen bond to a backbone that a substitution cannot reach. That difference is what a medicinal chemist is asking.

Eight interaction types, each implemented against a published criterion named in the rule’s own docstring, in interactions. What can take part is chemistry, and that is settled first: see ligand_chemistry, which records how confidently it knew.

The occupancy of each interaction is reported with the observation behind it. A contact present in three frames of five hundred and one present in four hundred and fifty are both “present”; only one means anything, and a fraction alone does not say which is which.

fastmdxplora.analysis.pl_interactions.ALL_KINDS = ('hydrophobic', 'hydrogen_bond', 'salt_bridge', 'pi_stacking', 'pi_cation', 'halogen_bond', 'metal_coordination', 'water_bridge')

what holds the ligand most often first. At module level so it can be the declared default of the kinds option rather than something filled in from None afterwards.

Type:: Every rule, in the order a reader would want them

class fastmdxplora.analysis.pl_interactions.ProteinLigandInteractions(*, ligand_resname=None, protein_selection='protein', ligand_chemistry=None, ligand_net_charge=None, kinds=('hydrophobic', 'hydrogen_bond', 'salt_bridge', 'pi_stacking', 'pi_cation', 'halogen_bond', 'metal_coordination', 'water_bridge'), minimum_occupancy=0.1, periodic=True, **kwargs)[source]

Bases: Analysis

What holds the ligand in place, and how well each contact is observed.

Parameters:

ligand_resname (str) – Ligand residue name. Supplied by the orchestrator.
protein_selection (str, default "protein") – The other side of the interaction. Usually the whole protein, and worth changing in three cases: a complex where only one chain matters ("protein and chainid 0"), a domain rather than the whole fold ("protein and resid 40 to 180"), or a receptor that is not a protein at all – "nucleic" for a ligand bound to DNA or RNA, which is a real and common case that the default name does not suggest.
ligand_chemistry (path, optional) – An SDF stating the ligand’s chemistry. Where a trajectory came from elsewhere and its residue name is not one the Chemical Component Dictionary knows, the bond orders otherwise have to be inferred from the coordinates – which is a guess, and a wrong bond order moves a hydrogen and invents or destroys a hydrogen bond.
ligand_net_charge (int, optional) – The ligand’s net charge, where you know it. It matters for the interactions that are claims about charge: perceiving it from coordinates is ambiguous more often than not, and guanidinium comes out as an anion if the first balancing charge is taken.
kinds (sequence of str, default all eight) – Which interactions to look for. Declared as the tuple rather than filled in from None, so anything reading the signature – a form drawing a control, the help, a config template – can say what it would do.
minimum_occupancy (float, default 0.1) – How often an interaction must appear before it counts towards a binding mode. Below this, a single fleeting contact splits one arrangement into two modes that are really one.
periodic (bool, default True) – Measure across the periodic boundary where the trajectory carries a unit cell.
**kwargs – Standard base-class options.
Outputs
-------
pl_interactions.dat – One row per interaction, with its occupancy, the number of frames it was present, and the number of separate times it formed.
pl_interactions.png – Occupancy per interaction, drawn so that thinly observed contacts are visibly thinly observed.

ALL_KINDS = ('hydrophobic', 'hydrogen_bond', 'salt_bridge', 'pi_stacking', 'pi_cation', 'halogen_bond', 'metal_coordination', 'water_bridge')

compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str | None

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str | None

description: str = 'Protein-ligand interactions by type': Human-readable description used in figure titles.

honours_selection: bool = False: This works out its own atoms, so a general selection has nothing to apply to.

name: str = 'pl_interactions': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:: result (DataFrame)
Return type:: None

requires_ligand: bool = True: True for analyses that only apply to protein-ligand complexes (e.g. ligand pose RMSD). The orchestrator runs these automatically when a ligand is present and skips them otherwise. They can still be explicitly requested via include.

save_data(result, path)[source]

The pair table, and the exact residue table beside it.

A residue’s occupancy is the union of its pairs’ frames, which the pair table cannot express: read from it, a residue touched through eight atoms is somewhere between the largest single pair and the sum of them all. The union is knowable at the point the contacts are counted, so it is written rather than left to be estimated.

Parameters:

result (DataFrame)
path (Path)

Return type:

How often an interaction was there, and how much watching that rests on.

A contact seen in 3 frames of 500 and a contact seen in 450 are both “present”. Only one of them means anything, and reporting both as an occupancy – 0.6 per cent and 90 per cent – hides that the first is three observations and the second is four hundred and fifty.

That distinction is the whole of this module. Occupancy is easy; saying how much observation it rests on is what stops a number being read as more than it is. The same mistake in a different guise ran through this software’s hydrogen bond count for a long time: bonds present in few frames were dropped from every frame, and the plot said one thing while the trajectory said another.

Two further cautions are built in rather than left to the reader.

Frames are not independent. Consecutive frames of a trajectory are correlated – a contact present at 10 ps is very likely present at 10.1 ps – so the naive standard error of a mean over frames is too small, often by a large factor. The number of independent observations is closer to the number of times the contact formed or broke than to the number of frames, and that is what is reported.

A transition rate needs transitions. A matrix built from three observed switches is arithmetic, not kinetics. Where too few were seen, the count is given and the rate is not.

class fastmdxplora.analysis.interaction_summary.Occupancy(kind, ligand_atom, protein_atom, frames_present, frames_total, episodes)[source]

Bases: object

How often one interaction was present, and how well that is known.

Parameters:

kind (str)
ligand_atom (int)
protein_atom (int)
frames_present (int)
frames_total (int)
episodes (int)

as_record()[source]

Return type:: dict[str, Any]

episodes: int: How many times it appeared after being absent. Consecutive frames are correlated, so this is much closer to the number of independent observations than the frame count is.

property fraction: float

frames_present: int: Frames in which it was present, of the frames examined.

frames_total: int

property is_well_sampled: bool

Whether the fraction rests on enough independent observation.

A contact that formed once and stayed was observed once, however many frames it then occupied.

kind: str

ligand_atom: int

protein_atom: int

property uncertainty: float

Standard error of the fraction, counting episodes not frames.

Using frames would give a number several times too small: a contact present in 450 consecutive frames has not been measured 450 times.

fastmdxplora.analysis.interaction_summary.binding_modes(contacts, n_frames, *, minimum_occupancy=0.1)[source]

Which combinations of interactions occurred together, and how often.

A binding mode is the set of interactions present in a frame. Frames sharing a set are in the same mode, and the modes are what a ligand moves between.

Interactions too rare to be part of a mode are left out first, or every frame becomes its own mode: a single fleeting contact would otherwise split one arrangement into two.

Parameters:

contacts (Any)
n_frames (int)
minimum_occupancy (float)

Return type:

fastmdxplora.analysis.interaction_summary.mode_transitions(per_frame, *, minimum_transitions=10)[source]

How the ligand moved between binding modes, where that can be said.

A transition matrix from a handful of observed switches is arithmetic rather than kinetics: three switches give a probability with an uncertainty larger than itself. So the switches are counted first, and the probabilities are reported only where there were enough of them.

This is the check other tools do not make. Computing the matrix is easy; knowing whether the trajectory supports it is the part that decides whether the answer means anything.

Parameters:

per_frame (Any)
minimum_transitions (int)

Return type:

fastmdxplora.analysis.interaction_summary.occupancies(contacts, n_frames)[source]

One entry per interaction, with the observation behind it.

Interactions are grouped by what they are and which atoms they join, so a hydrogen bond that came and went is one entry with an occupancy rather than many entries with none.

Parameters:

contacts (Any)
n_frames (int)

Return type:

list[Occupancy]

fastmdxplora.analysis.interaction_summary.residue_occupancies(contacts, n_frames, label_of)[source]

Exact per-residue occupancy: the union of its pairs’ frames.

A residue usually touches a ligand through several atom pairs, and how often the residue is in contact is the number of frames in which any of them is. That is a union, and a union cannot be recovered from a table of per-pair occupancies: pairs that fire in the same frames give the maximum, pairs that never coincide give the sum, and any real case lies between. Exported as a bracket, comparing this software’s numbers against another tool’s meant comparing an interval against a point, and residues touched through many atoms carried intervals wide enough to swallow the difference under test.

The masks that answer it exist while the contacts are being counted and were being discarded there. This keeps them long enough to take the union, so a residue occupancy leaves this software as a number rather than a range.

label_of maps a protein atom index to the residue label it belongs to.

Parameters:

contacts (Any)
n_frames (int)
label_of (Any)

Return type:

list[dict[str, Any]]

What holds a ligand in place, one interaction type at a time.

Counting contacts says how much of the protein the ligand touches. It does not say what is holding it: a salt bridge that a charge change would destroy, or hydrophobic packing that tolerates one. That difference is what a medicinal chemist is asking, and it is the difference between counting and describing.

Each rule here is implemented against a published criterion, named in its own docstring, and the thresholds are settings rather than constants – because the published values disagree, and a reader should be able to see which was used. Where a widely-used tool departs from the literature, the departure is recorded rather than quietly adopted.

Geometry only. Which atoms can take part is chemistry, and that is settled before this runs: see ligand_chemistry, which also records how confidently.

class fastmdxplora.analysis.interactions.Contact(kind, frame, ligand_atom, protein_atom, distance_nm, angle_deg=None)[source]

Bases: object

One interaction, in one frame.

Parameters:

kind (str)
frame (int)
ligand_atom (int)
protein_atom (int)
distance_nm (float)
angle_deg (float | None)

angle_deg: float | None = None: Set where the rule has one; None where it does not.

distance_nm: float

frame: int

kind: str

ligand_atom: int: Atom indices into the trajectory, ligand side first.

protein_atom: int

fastmdxplora.analysis.interactions.donors_and_acceptors(topology, atom_indices)[source]

Which atoms can donate a hydrogen bond, and which can accept one.

A donor is a nitrogen, oxygen or sulphur with a hydrogen bonded to it; an acceptor is any of those three. This is the criterion Baker and Hubbard used (Prog Biophys Mol Biol 44:97, 1984) and it is what every tool surveyed uses, because with explicit hydrogens present it needs no perception: the hydrogen is either bonded to the nitrogen or it is not.

Returns donors as (heavy, hydrogen) pairs, because the angle is measured at the hydrogen and a nitrogen with two hydrogens can donate twice.

Parameters:

topology (Any)
atom_indices (Any)

Return type:

tuple[list[tuple[int, int]], list[int]]

fastmdxplora.analysis.interactions.halogen_bonds(traj, chemistry, ligand_indices, protein_indices, *, distance_nm=0.35, donor_angle_deg=(130.0, 180.0), include_fluorine=False, periodic=True)[source]

A halogen on the ligand donating to an acceptor on the protein.

The halogen’s sigma-hole points along the C-X axis, so the interaction is directional: the carbon-halogen-acceptor angle has to be near straight. ProLIF requires 130 to 180 degrees within 3.5 A; PLIP uses 165 plus or minus 30 within 4.0 A. The narrower distance is used because it is the one the sigma-hole picture supports, and both are settings.

Fluorine is not counted unless asked for. See _HALOGENS for why, and for which tool disagrees.

Only the ligand donates. Proteins carry no halogens unless somebody has modified them, and a modified residue is something to say so about rather than to guess at.

Parameters:

traj (Any)
chemistry (Any)
ligand_indices (Any)
protein_indices (Any)
distance_nm (float)
donor_angle_deg (tuple[float, float])
include_fluorine (bool)
periodic (bool)

Return type:

fastmdxplora.analysis.interactions.hydrogen_bonds(traj, ligand_indices, protein_indices, *, distance_nm=0.35, angle_deg=120.0, periodic=True)[source]

Hydrogen bonds between ligand and protein, in every frame.

The criterion is the literature standard: donor to acceptor within 3.5 A, and the donor-hydrogen-acceptor angle above 120 degrees (Baker & Hubbard 1984; McDonald & Thornton, J Mol Biol 238:777, 1994).

PLIP uses 4.1 A and 100 degrees, which it says is deliberate – refined against low-resolution crystal structures where hydrogen positions are not known. That reasoning does not apply here: these frames have hydrogens placed by the force field, so the angle is measured rather than inferred, and the stricter criterion is the one the measurement supports. PLIP’s values remain reachable by setting them.

Both directions are found: the ligand donating to the protein and the protein donating to the ligand are different interactions, and a ligand that can only accept is a fact about the ligand worth seeing.

Parameters:

traj (Any)
ligand_indices (Any)
protein_indices (Any)
distance_nm (float)
angle_deg (float)
periodic (bool)

Return type:

fastmdxplora.analysis.interactions.hydrophobic_atoms(topology, atom_indices)[source]

Carbons whose only neighbours are carbon or hydrogen.

The criterion PLIP states and ProLIF’s SMARTS encodes. A carbon bonded to nitrogen, oxygen or fluorine is polarised enough that treating it as hydrophobic would count a polar contact twice – once here and once as a hydrogen bond.

Parameters:

topology (Any)
atom_indices (Any)

Return type:

list[int]

fastmdxplora.analysis.interactions.hydrophobic_contacts(traj, ligand_indices, protein_indices, *, distance_nm=0.4, periodic=True)[source]

Hydrophobic contacts between ligand and protein, in every frame.

A pair of hydrophobic carbons within 4.0 A, which is PLIP’s threshold; ProLIF uses 4.5 A. There is no angle: hydrophobic association is entropic rather than directional, so there is no geometry to require, which is why every tool uses a distance alone and a generous one.

The count is large by nature – PLIP notes it can exceed every other type combined – so it is reported per atom pair and left to the caller to reduce. Reducing here would bake in one view of which contact represents a residue, and the reduction PLIP applies is one choice among several.

Parameters:

traj (Any)
ligand_indices (Any)
protein_indices (Any)
distance_nm (float)
periodic (bool)

Return type:

fastmdxplora.analysis.interactions.ligand_aromatic_rings(chemistry, atom_indices)[source]

Aromatic rings in the ligand, from its resolved chemistry.

Aromaticity is a chemical fact rather than a geometric one – a flat ring of carbons is not necessarily aromatic – so it comes from the chemistry, which is why that has to be resolved first.

Parameters:

chemistry (Any)
atom_indices (Any)

Return type:

list[list[int]]

fastmdxplora.analysis.interactions.ligand_charged_groups(chemistry, atom_indices)[source]

Charged groups on the ligand, from its resolved chemistry.

Formal charges come from the chemistry, not from the coordinates, which is why it has to be resolved first. A carboxylate carries -1 across two oxygens whichever one the file happens to mark, so the charge is spread over the group it is delocalised across.

Parameters:

chemistry (Any)
atom_indices (Any)

Return type:

tuple[list[list[int]], list[list[int]]]

fastmdxplora.analysis.interactions.metal_coordination(traj, ligand_indices, protein_indices, *, distance_nm=0.3, periodic=True)[source]

A metal ion coordinated by the ligand, or coordinating it.

A donor atom within 3.0 A of the metal, which is PLIP’s threshold. No angle: the geometry of a metal centre is a property of the whole coordination shell rather than of any one contact, and PLIP fits the shell to known geometries afterwards. That fitting is not done here, because which targets are superfluous to a coordination number is a judgement, and reporting each contact leaves it visible.

The metal may be on either side. An ion in the structure is often neither protein nor ligand in the way a selection divides them, so both directions are searched.

Parameters:

traj (Any)
ligand_indices (Any)
protein_indices (Any)
distance_nm (float)
periodic (bool)

Return type:

fastmdxplora.analysis.interactions.pi_cation(traj, chemistry, ligand_indices, protein_indices, *, distance_nm=0.6, offset_nm=0.2, allow_ambiguous_charge=False, periodic=True)[source]

A positive charge sitting over an aromatic ring, in every frame.

PLIP’s criterion: the charge within 6.0 A of the ring centre, and the offset from the ring’s axis below 2.0 A. The offset matters more here than the distance: a cation beside a ring at 5 A is not interacting with its face, and distance alone cannot tell the two apart.

Refuses on an undetermined ligand charge for the same reason a salt bridge does – it is a claim about charge.

Parameters:

traj (Any)
chemistry (Any)
ligand_indices (Any)
protein_indices (Any)
distance_nm (float)
offset_nm (float)
allow_ambiguous_charge (bool)
periodic (bool)

Return type:

fastmdxplora.analysis.interactions.pi_stacking(traj, chemistry, ligand_indices, protein_indices, *, distance_nm=0.55, angle_tolerance_deg=30.0, offset_nm=0.2, periodic=True)[source]

Aromatic rings stacked on each other, in every frame.

PLIP’s criterion: ring centres within 5.5 A, the angle between the planes within 30 degrees of parallel or of perpendicular, and the offset between the centres below 2.0 A – about the radius of benzene plus a little.

Both arrangements are reported and named. Parallel stacking and the T-shaped, edge-to-face arrangement are different interactions with different geometries, and a ligand that stacks one way and not the other is telling you something about the pocket.

Parameters:

traj (Any)
chemistry (Any)
ligand_indices (Any)
protein_indices (Any)
distance_nm (float)
angle_tolerance_deg (float)
offset_nm (float)
periodic (bool)

Return type:

fastmdxplora.analysis.interactions.protein_aromatic_rings(topology, atom_indices)[source]

Aromatic rings in the protein, from the residues that have them.

Parameters:

topology (Any)
atom_indices (Any)

Return type:

list[list[int]]

fastmdxplora.analysis.interactions.protein_charged_groups(topology, atom_indices)[source]

Charged side chains, as groups of atoms rather than single atoms.

A carboxylate’s charge is shared between two oxygens and a guanidinium’s across three nitrogens, so the distance that matters is to the group’s centre, not to whichever atom happens to be nearest. Measuring to the nearest atom would make the same salt bridge look shorter from one side than the other.

Parameters:

topology (Any)
atom_indices (Any)

Return type:

tuple[list[list[int]], list[list[int]]]

fastmdxplora.analysis.interactions.residues_not_covered(topology, atom_indices)[source]

Residues in a selection that the charge and ring tables do not know.

Those tables are the twenty standard amino acids, which is why they need no perception. The cost is that anything else falls through them silently: point this at DNA and the hydrogen bonds come out right, because they are found from elements and bonds, while the salt bridges and the stacking come out as zero – though a phosphate is charged and a nucleobase is aromatic.

Zero is an answer. “These residues were not examined for charge or aromaticity” is a different one, and the true one.

Returned as a count per residue name so the caller can report it rather than deciding on its own whether it matters: a single modified residue in a large protein is a footnote, and a selection made entirely of nucleotides is not.

Parameters:

topology (Any)
atom_indices (Any)

Return type:

dict[str, int]

fastmdxplora.analysis.interactions.salt_bridges(traj, chemistry, ligand_indices, protein_indices, *, distance_nm=0.45, allow_ambiguous_charge=False, periodic=True)[source]

Salt bridges between ligand and protein, in every frame.

Opposite charges within 4.5 A, which is ProLIF’s threshold; PLIP uses 5.5 A. There is no angle: the interaction is electrostatic and has no preferred direction, which is why every tool surveyed uses a distance alone.

Refuses where the ligand’s charge was not determined. A salt bridge is a claim about charge, and the charge is exactly what perception from coordinates is worst at: guanidinium is +1, and -1 also balances, so a guessed charge can make a cation look like an anion and invent the bridge it was supposed to detect. Where the chemistry was resolved rather than perceived this does not arise; where it was not, state the charge or pass allow_ambiguous_charge knowing what it means.

Parameters:

traj (Any)
chemistry (Any)
ligand_indices (Any)
protein_indices (Any)
distance_nm (float)
allow_ambiguous_charge (bool)
periodic (bool)

Return type:

fastmdxplora.analysis.interactions.water_bridges(traj, ligand_indices, protein_indices, water_indices, *, min_distance_nm=0.25, max_distance_nm=0.41, omega_deg=(71.0, 140.0), periodic=True)[source]

A water molecule hydrogen-bonded to both the ligand and the protein.

PLIP’s criterion: the water oxygen between 2.5 and 4.1 A of a polar atom on each side, and the angle at the water – between the two partners, measured at its oxygen – between 71 and 140 degrees. The lower bound matters as much as the upper: a water in line with both is not bridging them, it is simply between them.

Only single-water bridges are found. Two waters can bridge a gap, and three, and at some chain length the claim stops meaning anything about binding; PLIP draws the line at one and the same line is drawn here. Where a longer chain matters, it is a different question that deserves asking directly rather than falling out of this.

Waters have to be given. Which oxygens count as solvent is a selection, and a run that stripped its waters has none to offer – an empty result there means the trajectory holds no water, not that no bridges formed.

Parameters:

traj (Any)
ligand_indices (Any)
protein_indices (Any)
water_indices (Any)
min_distance_nm (float)
max_distance_nm (float)
omega_deg (tuple[float, float])
periodic (bool)

Return type:

The potential of mean force from an umbrella study.

The free energy along the pulled coordinate, assembled from the umbrella windows. Minima are the states the system prefers; the height between two minima is the cost of crossing, in kJ/mol, and is the number an umbrella study is run to obtain. Windows must overlap for the profile to mean anything: where they do not, the gap is reported as a gap rather than bridged, and bins nobody sampled stay blank in the drawing rather than becoming zeroes that would read as a spurious minimum. Only an umbrella study has windows to read; on any other run this has nothing to say, and says so.

This reads what the simulation phase already computed rather than recomputing it – the stitching is delicate, and doing it twice would invite two answers to the same question.

(Why this analysis exists at all: for an umbrella run the free energy along the coordinate is the study, and it was once written to pmf.json and stopped there – no figure, no manifest entry, no mention in the report, while sixteen lesser curves each got all three.)

class fastmdxplora.analysis.pmf.PMF(*, selection=None, output_dir=None, title=None, xlabel=None, ylabel=None, figsize=None, xunit=None, **options)[source]

Bases: Analysis

The free energy along an umbrella study’s coordinate.

Parameters:

selection (str | None)
output_dir (Path)
title (str | None)
xlabel (str | None)
ylabel (str | None)
figsize (tuple[float, float] | None)
xunit (str | None)
options (dict[str, Any])

compute(traj)[source]

Read the curve, ignoring the trajectory.

The trajectory of any one window is not the study: it is a system held at one position by a spring, and its distribution says nothing about the free energy on its own.

Parameters:: traj (Any)
Return type:: dict[str, Any]

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str

description: str = 'Potential of mean force': Human-readable description used in figure titles.

honours_selection: bool = False

the coordinate was chosen when the windows were planned, and a selection here would describe a different question.

Type:: There is nothing to select

name: str = 'pmf': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (dict[str, Any])
ax (Axes)

Return type:

None

requires_umbrella = True: Only where an umbrella study was run, which is what pmf.json records.

save_data(result, path)[source]

Two columns, with unsampled bins left as nan rather than zeroed.

nan rather than blank, because a blank field makes the row two-column in some places and one-column in others, and every reader of a whitespace table then disagrees about what it holds.

Parameters:

result (dict[str, Any])
path (Path)

Return type:

time_series: bool = False: Not per frame. The x axis is the collective variable, not time.

The free energy surface from a metadynamics run.

The free energy over the chosen collective variable, reconstructed from the bias the run deposited. Read it like a landscape: basins are the states the system visits, and the walls between them are the barriers, in kJ/mol. The surface arrives with its own convergence evidence – basin transitions counted, drift measured – and a run whose bias has not settled still gets its picture, drawn and clearly labelled provisional, with a note beside it saying exactly what is missing. Only a metadynamics run deposits a bias to read; elsewhere this has nothing to say, and says so.

The frames themselves are not averaged: a metadynamics trajectory is deliberately not a Boltzmann ensemble – that is the point of the method – so the surface comes from the hills the simulation phase recorded, and this reads that record rather than trying to recompute it, which would invite two answers to the same question.

(Why this analysis exists at all: the surface was once written to metadynamics_surface.json and stopped there – no figure, no manifest entry, no mention in the report – the same gap the umbrella study had.)

class fastmdxplora.analysis.metad_surface.MetadynamicsSurface(*, selection=None, output_dir=None, title=None, xlabel=None, ylabel=None, figsize=None, xunit=None, **options)[source]

Bases: Analysis

The free energy along a metadynamics run’s collective variable.

Parameters:

selection (str | None)
output_dir (Path)
title (str | None)
xlabel (str | None)
ylabel (str | None)
figsize (tuple[float, float] | None)
xunit (str | None)
options (dict[str, Any])

compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str

description: str = 'Metadynamics free energy surface': Human-readable description used in figure titles.

honours_selection: bool = False: The coordinate was chosen when the bias was planned; a selection here would describe a different question.

name: str = 'metad_surface': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (dict[str, Any])
ax (Axes)

Return type:

None

requires_metadynamics = True: Only where a metadynamics run produced a record.

save_data(result, path)[source]

Write the computed result to path.

Default behaviour:

1-D / 2-D numpy arrays: np.savetxt (whitespace-delimited).
pandas DataFrames: to_csv (comma-delimited).
Anything else: subclass must override.

Returns the path actually written.

Parameters:

result (dict[str, Any])
path (Path)

Return type:

time_series: bool = False

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

The work done by a steered pull.

The cumulative work done on the system as the anchor moves, drawn against the pulled distance. The shape of the curve matters more than the total: a pull that accumulates work smoothly has met resistance all the way, while one that accumulates it in a step has snapped past something – an unbinding event, a broken contact – and the total is the same in both cases, so the total alone hides which happened. Only a steered run records a pull to read; elsewhere this has nothing to say, and says so.

What this is not is a free energy. The work depends on how fast the anchor moved – a fast pull works against solvent and strain as well as the interactions being measured. Jarzynski’s equality recovers a free energy from many pulls, and its average is dominated by rare low-work trajectories, so it needs many more than feels reasonable. One pull gives a pathway, and the figure says so.

(Why this analysis exists at all: the work was once written to steered_work.json and stopped there – no figure, no manifest entry, no mention in the report – the same gap the umbrella study and the metadynamics surface had.)

class fastmdxplora.analysis.steered_work.SteeredWork(*, selection=None, output_dir=None, title=None, xlabel=None, ylabel=None, figsize=None, xunit=None, **options)[source]

Bases: Analysis

Work against the coordinate, for a steered run.

Parameters:

selection (str | None)
output_dir (Path)
title (str | None)
xlabel (str | None)
ylabel (str | None)
figsize (tuple[float, float] | None)
xunit (str | None)
options (dict[str, Any])

compute(traj)[source]

Compute the analysis. Must be deterministic and side-effect-free.

Parameters:: traj (mdtraj.Trajectory) – The trajectory to analyze, already sliced/strided as the user requested. The analysis should respect self.selection if relevant.
Returns:: The analysis result. Most commonly a 1-D or 2-D NumPy array or a pandas DataFrame. The same object is later passed to plot() and save_data().
Return type:: Any

default_selection: str | None = None: Default atom selection (MDTraj selection language). None means “use the whole trajectory”. Subclasses override when an analysis only makes sense on a subset of atoms (e.g. RMSF on CA atoms).

default_xlabel()[source]

X-axis label when the user has not overridden it.

Override in subclasses to set a domain-specific default. Returning None means “leave whatever the plot() method set”. User-supplied xlabel= at construction always wins regardless.

Return type:: str

default_ylabel()[source]

Y-axis label when the user has not overridden it. See default_xlabel().

Return type:: str

description: str = 'Work done by the pull': Human-readable description used in figure titles.

honours_selection: bool = False: Whether selection means anything for this analysis. An analysis that decides its own atoms – a protein-ligand measure works out both sides from the ligand’s residue name, dihedrals from the backbone – has nothing to apply it to, and offering a control that does nothing is worse than not offering one: it looks like it worked.

name: str = 'steered_work': Short identifier used in output paths and the manifest.

plot(result, ax)[source]

Render result onto ax. Do not call plt.show().

Parameters:

result (dict[str, Any])
ax (Axes)

Return type:

None

requires_steered = True: Only where a steered run produced a record.

save_data(result, path)[source]

Write the computed result to path.

Default behaviour:

1-D / 2-D numpy arrays: np.savetxt (whitespace-delimited).
pandas DataFrames: to_csv (comma-delimited).
Anything else: subclass must override.

Returns the path actually written.

Parameters:

result (dict[str, Any])
path (Path)

Return type:

time_series: bool = False

Whether compute returns one value per frame.

A quantity measured every frame has a mean, and a mean is not a measurement until two things are known: whether the system had settled by the time the averaging started, and how many independent observations the average rests on. Both are recorded automatically for an analysis that says yes here.

Declared rather than inferred from the array’s length. A per-atom result on a trajectory that happens to have as many frames as atoms would otherwise be summarised as though it were a time series, and the numbers would look right.

Recovering unbiased averages from a metadynamics trajectory.

A metadynamics run deliberately distorts the ensemble: that is how it escapes minima. A state the bias filled early is visited more often than equilibrium would give, and one filled late less. An average over the frames is therefore an average over a distribution nobody wanted, and reporting it as a property of the system is wrong in a way that looks entirely plausible.

The distortion is known, which is what makes this recoverable. Each frame was sampled under a bias V(s) at its own value of the collective variable, so weighting it by exp(V(s)/kT) undoes the tilt and the weighted average is the unbiased one.

Two things make that statement less simple than it looks.

The bias grows. V is not one function but a sequence of them: the hills deposited by frame ten are not the hills deposited by frame ten thousand. Weighting every frame by the final bias is the common shortcut and it is wrong early in a run, where the bias the frame actually experienced was much smaller. This uses the bias as it stood when each frame was written, summing only the hills deposited before it.

Well-tempered runs converge to a scaled free energy, not to -F: the bias approaches -(1 - 1/gamma) F. Tiwary and Parrinello’s estimator handles this with a time-dependent offset c(t); what is implemented here is the simpler form that holds once the bias has settled, and the caller is told when it has not. A surface still filling gives weights that are only approximately right, which is worth having and worth saying.

Tiwary, P.; Parrinello, M. A time-independent free energy estimator for metadynamics. J Phys Chem B 2015, 119, 736. Branduardi, D.; Bussi, G.; Parrinello, M. Metadynamics with adaptive Gaussians. J Chem Theory Comput 2012, 8, 2247.

fastmdxplora.analysis.reweight.KB_KJ_PER_MOL_K = 0.008314462618: Boltzmann’s constant in kJ/mol/K, as the rest of the package uses it.

class fastmdxplora.analysis.reweight.Weights(values, effective_sample_size, settled, note='')[source]

Bases: object

Per-frame weights, and what they can and cannot be trusted for.

Parameters:

values (ndarray)
effective_sample_size (float)
settled (bool)
note (str)

effective_sample_size: float

How many independent frames the weighted average really rests on.

The Kish estimate: (sum w)^2 / sum w^2. A weighted mean over a thousand frames whose weight is concentrated in five of them is a mean over five, and quoting it as a thousand overstates it by a factor of fourteen. This is the number that says whether a reweighted average means anything.

note: str = ''

settled: bool

Whether the bias had stopped growing when the run ended.

Weights from a surface still filling are approximately right rather than right, because the simple estimator assumes a converged bias.

property usable_fraction: float: The effective sample size as a fraction of the frames.

values: ndarray: One weight per frame, normalised to sum to the frame count.

fastmdxplora.analysis.reweight.bias_at_each_frame(hills_times_ps, hills_centres, hills_sigmas, hills_heights, frame_times_ps, frame_values)[source]

The bias each frame actually felt, from the hills laid down before it.

Not the final bias applied to every frame. A frame written in the first picosecond was sampled under almost no bias, and weighting it as though it had felt the whole of a run’s deposition inflates it by orders of magnitude – which shows up as an effective sample size of one or two.

Parameters:

hills_times_ps (ndarray)
hills_centres (ndarray)
hills_sigmas (ndarray)
hills_heights (ndarray)
frame_times_ps (ndarray)
frame_values (ndarray)

Return type:

fastmdxplora.analysis.reweight.read_colvar(path)[source]

Times and collective-variable values, from a PLUMED COLVAR.

Parameters:: path (str | Path)
Return type:: tuple[ndarray, ndarray]

fastmdxplora.analysis.reweight.weighted_mean(values, weights)[source]

A weighted average, and nothing more clever than that.

Parameters:

values (ndarray)
weights (Weights)

Return type:

float

fastmdxplora.analysis.reweight.weighted_standard_deviation(values, weights)[source]

The spread about the weighted mean, on the effective sample size.

Dividing by the frame count would understate it: a thousand frames whose weight sits in fifty of them carry the uncertainty of fifty.

Parameters:

values (ndarray)
weights (Weights)

Return type:

float

fastmdxplora.analysis.reweight.weights_for_run(simulation_dir, frame_times_ps, *, temperature_K=300.0)[source]

Weights for a run’s trajectory frames, or None where there are none.

frame_times_ps are the trajectory’s own times, which is what makes this correct rather than approximately correct: the collective variable is recorded on PLUMED’s stride and the trajectory on its own, and the two need not coincide. The variable is interpolated onto the frame times rather than assumed to line up with them.

Parameters:

simulation_dir (str | Path)
frame_times_ps (ndarray)
temperature_K (float)

Return type:

Weights | None

fastmdxplora.analysis.reweight.weights_from_bias(bias_kjmol, *, temperature_K=300.0, settled=True)[source]

Turn a per-frame bias into weights that undo it.

Normalised by the largest bias before exponentiating, because exp of a hundred kilojoules over RT overflows a float and returns infinity for every frame – from which the weighted mean is a nan and the effective sample size a nan, and nothing says why.

Parameters:

bias_kjmol (ndarray)
temperature_K (float)
settled (bool)

Return type:

Weights

Reporting a biased run’s analyses as the unbiased averages they are not.

Sixteen analyses run on the production trajectory and each reports a mean. On a metadynamics run every one of those means is an average over a distribution the bias deliberately flattened, and reported without qualification it reads as a measurement of the system. The report already says so in its methods. This computes the correction it points at.

Only metadynamics gets weights, and the other two biased methods are left alone deliberately rather than overlooked. An umbrella window is a separate simulation held where it was put; what combines the windows is the potential of mean force, not a weighted average across them. A steered pull is not an equilibrium ensemble at all, so there is no set of weights that turns it into one. Both are stated in the methods text and neither is fixable here.

What comes out is both numbers side by side – the raw average over the biased frames and the reweighted one – with the effective sample size that says how much the second rests on. A reweighted mean over a thousand frames whose weight sits in five of them is a mean over five, and quoting it without that number is the failure mode this is built to avoid.

fastmdxplora.analysis.reweighted_averages.C_OF_T_CHECKPOINTS = 200: How many points the offset is evaluated at before being interpolated onto the frames. c(t) is smooth in time – it tracks the filling of the surface, not the motion of the system – so a couple of hundred is ample, and evaluating it per frame would cost a great deal for no change in the answer.

fastmdxplora.analysis.reweighted_averages.C_OF_T_GRID_POINTS = 320: Points across the collective variable for the two integrals below.

fastmdxplora.analysis.reweighted_averages.FALLBACK_TEMPERATURE_K = 300.0: If the run reached equilibrium temperature somewhere other than 300 K and nothing recorded it, the weights would be silently wrong by the ratio of the temperatures. Falling back is still better than refusing, but it is written into the record so it is not mistaken for a reading.

fastmdxplora.analysis.reweighted_averages.METHOD_MARKERS = (('umbrella', 'umbrella.plumed'), ('steered', 'steered.plumed'), ('metadynamics', 'metadynamics.plumed')): What each method leaves beside the run. The script is written under the method’s own name, so it says which bias was applied even where the result file is missing because the run stopped early.

fastmdxplora.analysis.reweighted_averages.NOT_REWEIGHTABLE = {'steered': 'Production was a steered pull, which is not an equilibrium ensemble: the system was dragged along the coordinate rather than sampling it. These averages describe the pulling, and no reweighting recovers an equilibrium average from a single non-equilibrium trajectory.', 'umbrella': 'This run is one window of an umbrella study, held at its own position on the coordinate by a harmonic restraint. Its averages describe a system held there and are not measurements of the unrestrained system, nor comparable between windows, which differ because the restraints differ. What combines the windows is the potential of mean force, not an average of any quantity across them.'}: Why no set of weights recovers an equilibrium average from these two. Both are stated in the report’s methods text; these are the same claims, put where the numbers are.

fastmdxplora.analysis.reweighted_averages.USABLE_EFFECTIVE_FRAMES = 50.0: Below this many effective frames a reweighted average is reported with a warning attached. The number is not a threshold for correctness – the estimate is what it is – but for whether it carries enough frames to be worth reading, and fifty is where a mean stops being dominated by a handful of them.

fastmdxplora.analysis.reweighted_averages.before_deposition(hills, times_ps)[source]

Frame times nudged so a hill laid at time t is not felt at time t.

PLUMED prints the bias for a step before depositing that step’s hill, and HILLS and COLVAR usually share a stride, so counting a hill as felt at its own deposition time is wrong on every row. On a real run it put the reconstruction out by 1.200 kJ/mol against PLUMED’s own record – exactly one hill at the configured height, which is what made it identifiable.

The nudge is a small fraction of the deposition interval, which is the smallest time difference that matters here, so it moves each frame off the boundary without moving it past the previous hill.

Parameters:

hills (Any)
times_ps (Any)

Return type:

fastmdxplora.analysis.reweighted_averages.biasing_method(output_dir)[source]

Which method biased this run, from what it left on disk.

Parameters:: output_dir (Path)
Return type:: str | None

fastmdxplora.analysis.reweighted_averages.c_of_t(hills, times_ps, temperature_K, periodic=False)[source]

The offset that makes a running bias usable for reweighting.

The bias grows as hills accumulate, so exp(V(s,t)/RT) grows with time wherever the system happens to be. Weighting by it alone therefore ranks frames by when they were written: on a four-nanosecond well-tempered run the last fifth of the frames carried the entire weight, and a reweighted average over five hundred frames rested on seven. That is not a bias that has not settled – it happens to a fully converged run – and no warning about convergence covers it.

Tiwary and Parrinello’s c(t) is the term that removes it: the free-energy offset of the biased ensemble at time t, so that V - c(t) measures where the system is on the coordinate rather than how late it is in the run. For a well-tempered run with bias factor y,

c(t) = RT ln [ int ds e^{(y/(y-1)) V(s,t)/RT}
/ int ds e^{(1/(y-1)) V(s,t)/RT} ]

and the y -> infinity limit of that, which is the ratio to a flat integral, is the form for a run that was not tempered.

Parameters:

hills (Any)
times_ps (ndarray)
temperature_K (float)
periodic (bool | tuple[bool, ...])

Return type:

fastmdxplora.analysis.reweighted_averages.deposited_heights(hills)[source]

The heights actually added to the bias, from what HILLS records.

For a well-tempered run PLUMED does not store the height it deposited. It stores that height multiplied by y/(y-1), where y is the bias factor, so that summing the file gives the free energy directly – which is the convention metad_surface relies on and must keep.

Reconstructing the bias needs it undone. Summing the stored heights overstates the bias by y/(y-1): 11.1% at a bias factor of 10, which is what a real run showed against PLUMED’s own record of the same quantity. That error does not cancel between V and c(t) – both scale by the same factor, so their difference scales too, and since the weights go as exp((V - c(t))/RT) a scaled exponent is an effective-temperature error. It sharpens the weights, understates the effective sample size, and biases every average that rests on them.

A run that was not tempered stores what it deposited, and needs nothing.

Parameters:: hills (Any)
Return type:: ndarray

fastmdxplora.analysis.reweighted_averages.felt_bias(hills, heights, *, times, values, periodic)[source]

The bias each frame felt, measured the short way round on a circle.

A frame at -175 degrees is fifteen degrees from a hill at +170, not three hundred and forty-five, and summed straight it feels none of that hill – so its weight ignores the bias that was actually there.

Done by summing the plain calculation over the coordinate’s periodic images rather than by teaching the weighting maths about circles. A periodic Gaussian is the sum over images, so this is not an approximation of the right answer, it is the right answer; and with a hill width of a fifth of a radian the next image sits eighteen widths away and contributes 1e-70, so one image either side is the whole of it. reweight.py holds generic weighting arithmetic and PLUMED’s conventions belong here.

Parameters:

hills (Any)
heights (ndarray)
times (Any)
values (Any)
periodic (bool | tuple[bool, ...])

Return type:

fastmdxplora.analysis.reweighted_averages.frame_labels(result, n_frames)[source]

Per-frame categorical labels an analysis reported, keyed by method.

Clustering returns a mapping of method name to labels because it runs several and they disagree in useful ways; a single array is accepted too so this does not assume that shape.

Parameters:

result (Any)
n_frames (int)

Return type:

dict[str, ndarray]

fastmdxplora.analysis.reweighted_averages.frame_series(result, spec, n_frames)[source]

The per-frame scalar an analysis reported, or None if it has none.

The length check is the point. SASA in per-residue mode returns a row per residue per frame, and averaging that against per-frame weights would line the two up by position and produce a number that is not wrong so much as meaningless.

Parameters:

result (Any)
spec (tuple[str | None, str])
n_frames (int)

Return type:

ndarray | None

fastmdxplora.analysis.reweighted_averages.populations(labels, weights)[source]

How often each state was really visited, against how often it appeared.

A population is the mean of an indicator, so it reweights exactly as any other average does – and it is the quantity a bias distorts most, since escaping a well is what the bias is for.

Parameters:

labels (ndarray)
weights (Weights)

Return type:

list[dict[str, Any]]

fastmdxplora.analysis.reweighted_averages.read_colvar(path)[source]

Read a PLUMED COLVAR file by its declared field names.

The header names the columns, and reading by name rather than position matters here: a metadynamics COLVAR is time, cv, metad.bias but a steered one is time, cv, pull.work, pull.bias, and column two means something different in each.

Parameters:: path (str | Path)
Return type:: dict[str, ndarray]

fastmdxplora.analysis.reweighted_averages.reweight_results(results, registry, *, n_frames, frame_times_ps, output_dir)[source]

Recompute every reweightable average against the metadynamics bias.

Returns the record written to disk, or None where no bias applies and there is nothing to correct.

Parameters:

results (dict[str, Any])
registry (dict[str, type])
n_frames (int)
frame_times_ps (Any)
output_dir (Path)

Return type:

dict[str, Any] | None

fastmdxplora.analysis.reweighted_averages.weights_for_run(output_dir, frame_times_ps)[source]

Per-frame weights undoing the metadynamics bias, where one was applied.

Returns (None, provenance) for any run that was not biased by metadynamics – plain MD, a steered pull, an umbrella window – because none of those has a set of weights that recovers an equilibrium average.

The bias each frame felt is rebuilt from the hills deposited before it, using that frame’s own value of the collective variable interpolated from COLVAR. Interpolating the variable rather than PLUMED’s bias column is deliberate: the variable is a smooth physical coordinate, while the bias is a sum of narrow Gaussians that interpolation would badly misstate between rows.

Parameters:

output_dir (Path)
frame_times_ps (Any)

Return type:

tuple[Weights | None, dict[str, Any]]

Statistics and provenance

Where a run settled, and how many independent samples it actually holds.

At the top level rather than under analysis, because it is not an analysis: nothing registers it, it produces no figure, and both the analyses and the report ask it the same question. It sat in the analysis package briefly and the report imported it from there, which reads as though measuring a correlation were something analyses do and reports borrow.

Ten analyses averaged over the whole production run without asking either question. A mean root-mean-square deviation of 2.3 Å is not a measurement until two things are known about it: whether the system had stopped changing by the time the averaging started, and how many independent observations the average rests on.

The first. A structure that has just been minimised, heated and pressure- equilibrated is still relaxing when production begins. Averaging from the first frame averages the relaxation together with the equilibrium, and the answer depends on how long the run was rather than on the system.

The second. Frames are not independent. A trajectory written every picosecond from a system whose fluctuations decorrelate over a hundred picoseconds has a hundred times fewer independent samples than it has frames, and a standard error computed as if each frame counted is wrong by a factor of ten. That is how a difference between two systems becomes significant on paper without being real.

Both follow from one quantity. The statistical inefficiency g is the number of frames per independent sample, so a series of n frames carries n / g of them. Chodera’s method chooses where to start averaging by maximising that count: discard too little and the relaxation is still in the average, discard too much and there is nothing left to average.

Chodera, J. D. A simple method for automated equilibration detection in molecular simulations. J. Chem. Theory Comput. 2016, 12, 1799-1805.

class fastmdxplora.statistics.Settled(discard, inefficiency, effective_samples, mean, standard_error, standard_deviation)[source]

Bases: object

What a series supports, once the relaxation is out of it.

Parameters:

discard (int)
inefficiency (float)
effective_samples (float)
mean (float)
standard_error (float)
standard_deviation (float)

as_record()[source]

Return type:: dict[str, Any]

discard: int: Frames discarded before averaging.

effective_samples: float: Independent samples in what is left.

inefficiency: float: Frames per independent sample. One means every frame counts.

mean: float

standard_deviation: float: The spread of the series itself, which is a property of the system rather than of how long it was watched.

standard_error: float: The uncertainty on the mean, from the effective samples rather than from the frame count. Using the frame count is what makes a difference look significant when it is not.

fastmdxplora.statistics.correlation_is_resolved(series)[source]

Whether the series is long enough to measure its own correlation time.

Checked by halving it: an inefficiency the series can resolve does not change much when half the frames are taken away, and one it cannot moves a great deal. On an AR(1) series with a true inefficiency of 2000, four thousand frames gave 361 – which is not an error to tolerate but a number with the wrong meaning, since the independent-sample count built from it said eleven when the truth was two.

The obvious check does not work and it is worth saying why. A sample autocorrelation sums to roughly -1/2 whatever the series, so it goes negative on its own: on that series it crossed zero at lag 334 while the real correlation there was still 0.7. Asking where the correlation decayed answers a question about the estimator rather than the run.

Parameters:: series (ndarray)
Return type:: bool

fastmdxplora.statistics.detect_equilibration(series, *, steps=40)[source]

Where to start averaging, and what is left once you do.

Returns the number of frames to discard, the statistical inefficiency of what remains, and the effective sample count. The discard point is the one maximising that count, which is the trade Chodera’s method makes explicit: keeping the relaxation costs independence, and discarding it costs frames.

Candidate points are strided rather than exhaustive, because the count varies smoothly with where the average starts and evaluating every frame is quadratic for an answer no better.

Parameters:

series (ndarray)
steps (int)

Return type:

tuple[int, float, float]

fastmdxplora.statistics.statistical_inefficiency(series)[source]

Frames per independent sample: g = 1 + 2 sum (1 - t/n) C(t).

C is the normalised fluctuation autocorrelation. The sum is truncated at the first non-positive C, which is the standard convention: past that point the estimates are noise, and summing them adds variance rather than information.

A constant series has no fluctuations to correlate, so g is one: every frame agrees, and there is nothing for a correlation time to describe.

Parameters:: series (ndarray)
Return type:: float

fastmdxplora.statistics.summarise(series, *, minimum_effective_samples=10.0)[source]

A mean with an honest error on it, or a reason there is not one.

The refusal is about independence, not length: a long run of highly correlated frames can hold fewer independent samples than a short run of uncorrelated ones, and it is the second number that decides what the mean is worth.

Parameters:

series (ndarray)
minimum_effective_samples (float)

Return type:

tuple[Settled | None, str | None]

Which code a run was actually made from.

A manifest recorded the version string and nothing else. That string is written by setuptools-scm at install time, so an editable install carries whatever the version was when pip install -e . was last run and drifts silently from then on. A real study came back stamped 2.3.0 for a run that used a feature 2.3.0 did not have: the manifest named a version in which the run could not have happened, and it is the number the report’s reproducibility section prints.

Where the package is imported from a source checkout, the commit says what the version string cannot. Where it is not – a conda or PyPI install – there is no checkout to ask, and the version string is the whole answer because the distribution was built from a tag.

Two decisions are worth stating, because the alternatives are defensible:

The checkout is found from the package, not the working directory. A run started inside some other repository would otherwise record that repository’s commit, which is worse than recording nothing: it is a precise claim about the wrong code.

A dirty tree is reported, not refused. Refusing would block the runs developers make all day, and this project’s habit is to say what a result cannot support rather than to withhold it. A commit beside uncommitted changes does not describe the code that ran, so the flag is what makes the commit honest rather than decorative.

fastmdxplora.provenance.described_structure(record)[source]

One line a person can read, or None where there is nothing to say.

Parameters:: record (dict[str, Any] | None)
Return type:: str | None

fastmdxplora.provenance.source_checkout()[source]

The checkout this package was imported from, or None if installed.

Walks up from the package rather than from the working directory: a run started inside another repository would otherwise report that repository’s commit, which is a precise claim about the wrong code.

Return type:: Path | None

fastmdxplora.provenance.source_provenance()[source]

The commit a run was made from, or None where there is no checkout.

Returns {"commit": ..., "dirty": bool}, and branch where the checkout is on one. dirty is the important field: with uncommitted changes the commit does not describe the code that ran, and saying so is what keeps the commit from being decorative.

Return type:: dict[str, Any] | None

fastmdxplora.provenance.structure_provenance(given, form, path)[source]

Which structure a run was made from, past the point the path says.

The manifest recorded the string somebody typed – 4hhb_cleaned.pdb, ../prep/final.pdb – and nothing else. That is an answer for about a week. The file gets moved, the directory tidied, a second copy made under a name that sorts better, and the run no longer says which structure it used. The report is built from the same field, so a methods section ends up stating that coordinates came from a filename, which is not something a reader can check.

Three things make it checkable, none of which costs anything beside a run of any length:

The digest, which names the bytes whatever became of the path. Two runs agree or they do not, and a file edited since stops matching the run that used it.

The entry the file names itself. A local structure is usually a deposited one that has been through a preparation step, and it keeps the header saying which. That answers “which entry” in exactly the case where the path has stopped answering it.

When, because deposited entries are revised. A run made before a revision used different coordinates from one made after it, and neither the identifier nor the filename records which side of it a run falls on.

Returns None where there is no file to describe – a sequence input, or a fetch that failed – because a record of nothing is worse than no record.

Parameters:

given (str)
form (str)
path (Path | str)

Return type:

dict[str, Any] | None

Why each step happens, said while it happens.

Molecular dynamics has a lot of steps that are obvious once you know them and opaque before that. Why is the protein put in a box of water? Why is the box bigger than the protein? Why heat it before running it? Why is there a separate stage where the volume can change?

A pipeline that does all of this silently is faster to use and teaches nothing. Somebody running their first simulation ends up with a trajectory they cannot defend, because they cannot say why any of it was done.

So FastMDXplora says what it is doing and why, and cites something worth reading where there is one. It is on by default and –no-explain turns it off, because the fourth time through it is noise.

The explanations are keyed, not matched. A call site names the explanation it wants, so an explanation cannot drift onto the wrong step and a step cannot quietly lose its explanation – both are checked. Free-text matching would have made the first failure silent.

fastmdxplora.explain.EXPLANATIONS: dict[str, Explanation] = {'convergence': Explanation(why='Frames are not independent observations. Consecutive frames are nearly the same structure, so a thousand frames may hold only a handful of independent samples -- and an error bar computed as though they were all independent is too small by a large factor. What is reported is how much the trajectory actually supports.', reference='Flyvbjerg & Petersen, Error estimates on averages of correlated data, J Chem Phys 1989'), 'ensemble_choice': Explanation(why='Which ensemble production runs in is a choice, and the usual one is NPT: it matches the constant-pressure conditions an experiment is done under, and the box is free to respond if the system changes shape. NVT production is also legitimate -- it is cheaper, and it removes volume fluctuation from anything sensitive to it -- but only at a density you know is right, and the way you learn that is to run NPT first and take the average box size from it. Going straight to NVT is not the same choice: it fixes the box at whatever solvation produced and simulates there for every step, which is the one option nobody intends.', reference="Braun et al., Best Practices for Foundations in Molecular Simulations, LiveCoMS 2019, Fig. 'Suggested equilibration workflow' (doi:10.33011/livecoms.1.1.5957)"), 'heterogens': Explanation(why='A crystal structure contains more than the protein. Some of it is biology -- a bound ligand, a structural metal -- and some is the chemistry that made the crystal grow, like glycerol or buffer molecules. Simulating the second kind wastes computation and can hold the protein in a shape the crystal imposed. Each one is classified rather than kept or dropped wholesale.', reference='Davis et al., The crystallographic heterogen problem, Acta Cryst D 2008'), 'interactions': Explanation(why='Counting contacts says how much of the protein a ligand touches. Typing them says what holds it there -- a salt bridge a charge change would destroy, or a hydrophobic packing that tolerates one. Those suggest different next experiments, which is why each contact is classified rather than counted.', reference='Adasme et al., PLIP 2021, Nucleic Acids Res 2021'), 'ligand_chemistry': Explanation(why="A structure file gives a ligand's atoms and where they are, but not its bond orders, its aromaticity or its charge -- and a force field needs all three. Getting a bond order wrong moves a hydrogen, and a moved hydrogen invents or destroys a hydrogen bond. The chemistry is looked up rather than guessed wherever it can be.", reference='Westbrook et al., The Chemical Component Dictionary, Bioinformatics 2015'), 'ligand_parameters': Explanation(why="The protein force field knows the twenty amino acids and nothing else, so a bound ligand has no parameters until somebody makes them. OpenFF generates them from the ligand's chemistry, which is why the chemistry had to be settled first.", reference='Qiu et al., OpenFF 2.0.0 Sage, J Chem Theory Comput 2021'), 'membrane': Explanation(why='A membrane protein simulated in water is not the protein: the hydrophobic belt that normally sits in the bilayer is exposed to solvent, and the helices splay apart. The lipids are packed around it so the protein sees the environment it evolved in.', reference='Lomize et al., OPM database and PPM server, Nucleic Acids Res 2012'), 'membrane_barostat': Explanation(why='A bilayer must be free to change thickness independently of its area, so pressure is coupled in the membrane plane and along the normal separately. An ordinary barostat scales all three directions together, which squeezes the membrane and gives the wrong area per lipid -- the number membrane simulations are validated against.', reference='Chow & Ferguson, Isothermal-isobaric molecular dynamics, Comput Phys Commun 1995'), 'metadynamics': Explanation(why='Some things are too slow to see by waiting -- a ligand leaving a pocket might take milliseconds, and a simulation runs for microseconds. Metadynamics fills in the energy landscape along the coordinate you named as the run proceeds, pushing the system out of wells it has already visited so it explores instead of sitting still. The free energy is recovered from the bias that was added.', reference='Barducci et al., Well-tempered metadynamics, Phys Rev Lett 2008'), 'minimize': Explanation(why='The starting structure has strain in it -- atoms slightly too close, bonds slightly too long -- from the experiment, from adding hydrogens, and from dropping the protein into water. At the temperature of a simulation that strain becomes violent motion. Minimisation walks the structure downhill to a nearby arrangement with no such forces in it, before anything moves.', reference=None), 'npt': Explanation(why="P replaces V: the box is now free to change size, and settles to the density real water has at this temperature and pressure. This is not a formality. Solvation leaves a gap around the solute, and in a small box that gap is a large share of the volume: measured here, a solute at 1.0 to 1.2 nm of padding packs near 0.90 g/mL against water's 1.0, while the same solvation at 2.0 nm reaches 0.96. Only a barostat closes it. A box that short has voids in it, which is wrong for anything you measure and a route to the run falling over. How far off your own box is gets reported rather than assumed.", reference='Aqvist et al., Molecular dynamics simulations of water and biomolecules with a Monte Carlo constant pressure algorithm, Chem Phys Lett 2004 (doi:10.1016/j.cplett.2003.12.039)'), 'nvt': Explanation(why='N, V and T are what is held fixed: the number of atoms, the volume of the box, and the temperature. The system starts at zero temperature with the atoms sitting still, and this brings it to the temperature you asked for while the box stays the size solvation made it. Temperature first and density second, because two things equilibrating at once is harder to diagnose when it goes wrong.', reference='Braun et al., Best Practices for Foundations in Molecular Simulations, LiveCoMS 2019 (doi:10.33011/livecoms.1.1.5957)'), 'production': Explanation(why='This is the part that is analysed. Everything before it was getting the system into a state worth measuring; from here the trajectory is a sample of how the system behaves at equilibrium, and the frames written now are the ones every later number comes from.', reference=None), 'protonation': Explanation(why='X-rays do not see hydrogens, so a crystal structure has none, and a simulation needs every one. Which histidines are protonated, and whether a glutamate is charged, depends on the local environment and on pH -- and those choices change the hydrogen bonding that holds the structure together. They are decided here rather than left to a default.', reference='Olsson et al., PROPKA3, J Chem Theory Comput 2011'), 'restraints': Explanation(why='A minimised structure is not at equilibrium, and heating it lets the solute move as well as the solvent -- side chains relax into the space crystal packing left, and a ligand drifts out of the pose that was measured. Holding the solute while the water arranges itself around it, then letting go in stages, means production starts from the structure somebody determined.', reference='Roe & Brooks, A protocol for preparing explicitly solvated systems, J Chem Phys 2020'), 'solvation': Explanation(why='Proteins fold the way they do because of water -- the hydrophobic effect is a property of the solvent, not of the protein. Simulated in vacuum, a structure collapses in on itself. The box is padded well beyond the protein so it never interacts with its own periodic image, and ions are added to physiological concentration because charges in real cells are screened.', reference='Jorgensen et al., Comparison of simple potential functions for water, J Chem Phys 1983')}: Keyed by step, so a call site asks for the explanation it means.

class fastmdxplora.explain.Explanation(why, reference=None)[source]

Bases: object

Why a step happens, and where to read more.

Parameters:

why (str)
reference (str | None)

as_text(width=74)[source]

Parameters:: width (int)
Return type:: str

reference: str | None = None: A citation worth following, when one exists.

why: str

fastmdxplora.explain.explain(key)[source]

The explanation for a step, or None where there is not one.

Returning None rather than raising, because a missing explanation should not stop a run – but the guard in the test suite fails on one, so it does not go unnoticed either.

Parameters:: key (str)
Return type:: Explanation | None

What is worth knowing before a run starts, rather than after.

A run says a great deal about what it is doing: which heterogens it discarded, that a force field wants hard truncation, that a metal in a site will not be held there by charge alone. All of it is true and all of it arrives once the run has begun – by which point the person who would have changed something has stopped watching, and on a cluster has gone home.

Most of it is decidable earlier. A structure and a set of settings are enough to say that this protein has a zinc in a site and these parameters will not hold it, or that this cutoff wants a bigger box than this padding will make. Said while somebody is still choosing, the same sentence changes what they do instead of explaining what already happened.

This is deliberately not a validator. A validator answers “will this run”, and these all describe things that run perfectly well and produce a result worth doubting. Nothing here refuses anything.

class fastmdxplora.advisories.Advisory(setting, summary, detail, remedy)[source]

Bases: object

Something worth knowing, and what to do about it.

setting names the field this is about, so an interface can show it beside the control rather than in a list somewhere else – which is the whole point of saying it early.

Parameters:

setting (str)
summary (str)
detail (str)
remedy (str)

as_text()[source]

Return type:: str

detail: str

remedy: str

setting: str

summary: str

fastmdxplora.advisories.DODECAHEDRON_WIDTH_FRACTION = 0.52: How much of solute-plus-padding survives as a dodecahedron’s smallest width. Measured rather than derived: 1.2 nm of padding on a 0.9 nm solute gives 1.70 nm across where the sum says 3.3, and the ratio held near a half at every padding tried.

fastmdxplora.advisories.advise(structure, settings=None)[source]

Everything worth saying about this structure and these settings.

structure is the report from gui.structure_info.count_structure() or anything shaped like it; settings is a config’s setup and simulation blocks flattened together. Both may be absent, and what can be said falls back to what is known.

Parameters:

structure (dict[str, Any] | None)
settings (dict[str, Any] | None)

Return type:

list[Advisory]

Campaigns

Batch orchestration — the single execution path for all runs.

Every FastMDXplora run goes through BatchExplorer. There is no separate “single run” path: a one-system config is simply a batch of one. This keeps one code path, one config shape (systems: is always a list), and one mental model.

Output layout adapts to the run count:

One run → flat, familiar layout written directly to the output directory: output/setup/, output/simulation/, etc., with the usual manifest.json and resolved_config.yml. No runs/ wrapper, no batch_manifest.json.

Many runs → each run in output/runs/<id>/ (a complete study), plus a top-level batch_manifest.json indexing them all.

Execution modes (execution: block):

sequential (default) — one run at a time, in process.

parallel — a process pool of workers runs at once. On GPU, set devices: [0, 1, ...] and each worker is pinned to a distinct device round-robin (one run per GPU), which is the only safe way to parallelize GPU MD — oversubscribing a single GPU is slower than sequential.

Process-based (not thread-based) parallelism is mandatory: OpenMM contexts and the GIL don’t share across threads. Each run is therefore dispatched to a subprocess via a module-level worker function.

class fastmdxplora.batch.explorer.BatchExplorer(config=None, *, config_data=None, output_dir=None, verbose=False, continue_on_error=None, force=False)[source]

Bases: object

Run one or more FastMDXplora studies (systems × sweep).

Parameters:

config (str | os.PathLike) – Path to a YAML config with a systems: list (and optionally sweep: / execution:).
output_dir (str | os.PathLike | None) – Root output directory. One run → written here directly; many runs → each in runs/<id>/. Defaults to a timestamped directory.
verbose (bool) – Forwarded to each run.
continue_on_error (bool | None) – Override the config’s execution.continue_on_error. If None, the config value (default True) is used.
config_data (dict[str, Any] | None)
force (bool)

Examples

>>> BatchExplorer(config="study.yml").run()

dry_run()[source]

Report the plan without executing anything.

Prints each run, its system, swept values, target output directory, and the phases that would execute, then returns a list of RunResult with status "planned" (no phases populated).

Return type:: list[‘RunResult’]

results: list[dict[str, Any]]

run()[source]

Execute every run. Returns per-run result records.

Return type:: list[dict[str, Any]]

fastmdxplora.batch.explorer.HEARTBEAT_SECONDS = 60.0: How long the terminal may say nothing before it says something. A window of a real study runs for hours, and every worker’s output goes to its own log so three of them do not interleave – which left the terminal silent for the whole run, with no way to tell a study that is working from one that has hung.

Batch sweep expansion.

Turns a batch configuration (a list of systems and a sweep of parameter axes) into a flat list of concrete single-run configurations — the Cartesian product of systems × sweep points.

This module is deliberately pure: it does no I/O and constructs no orchestrators. It just computes what runs should happen and with what options, so the expansion logic is trivially testable in isolation.

Sweep axes use dotted keys naming a phase option, e.g. simulation.temperature_K. Multiple axes form a full Cartesian product. Each entry in systems may carry its own per-phase option overrides, which are applied beneath the sweep values (sweep wins, since the sweep is the thing being varied).

Example

Input:

systems:
  - {id: a, system: a.pdb}
  - {id: b, system: b.pdb}
sweep:
  simulation.temperature_K: [300, 310]

produces four runs: a@300, a@310, b@300, b@310.

class fastmdxplora.batch.sweep.RunSpec(run_id, system, options=<factory>, sweep_values=<factory>, system_id='')[source]

Bases: object

One concrete run in a batch.

Parameters:

run_id (str)
system (str)
options (dict[str, dict[str, Any]])
sweep_values (dict[str, Any])
system_id (str)

run_id

Unique, human-readable identifier (also the run subdirectory name).

Type:: str

system

The system input for this run.

Type:: str

options

Per-phase options for this run (merged: base < system < sweep).

Type:: dict[str, dict[str, Any]]

sweep_values

The swept axis values for this run (dotted key -> value), recorded for the batch manifest.

Type:: dict[str, Any]

system_id

The originating system’s id.

Type:: str

options: dict[str, dict[str, Any]]

run_id: str

sweep_values: dict[str, Any]

system: str

system_id: str = ''

to_dict()[source]

Return type:: dict[str, Any]

exception fastmdxplora.batch.sweep.SweepError[source]

Bases: ValueError

Raised for malformed systems/sweep specifications.

fastmdxplora.batch.sweep.expand_runs(*, systems, sweep, base_options=None, base_system=None)[source]

Expand systems × sweep into a flat list of RunSpec.

Parameters:

systems (list of normalized system dicts, or None) – From normalize_systems(). If None, a single implicit system is used from base_system.
sweep (dict of axis -> values, or None) – From normalize_sweep(). If None/empty, no sweep is applied (one run per system).
base_options (dict, optional) – Project-level per-phase options that every run inherits (the top-level phase blocks of the config). Lowest priority.
base_system (str, optional) – Used when systems is None (single implicit system).

Returns:

One entry per (system × sweep-point), in deterministic order: systems outer, sweep axes inner (in declared order).

Return type:

list[RunSpec]

fastmdxplora.batch.sweep.is_batch_config(data)[source]

Return True if a parsed config requests batch mode.

Batch mode is active when the config contains a non-empty systems list or a non-empty sweep mapping.

Parameters:: data (dict[str, Any])
Return type:: bool

fastmdxplora.batch.sweep.normalize_summary_for_validation(data)[source]

Validate the batch keys of a parsed config (raises SweepError).

Checks that systems and sweep are well-formed, and that every sweep axis names a real option in its phase’s schema (so a typo like simulation.temperature_K is caught with a did-you-mean style message rather than silently producing a junk run).

Parameters:: data (dict[str, Any])
Return type:: None

fastmdxplora.batch.sweep.normalize_sweep(raw)[source]

Validate and normalize the sweep mapping.

Each key is a dotted phase.option and each value is a non-empty list of values to try. Returns the same shape with values coerced to lists (a scalar is treated as a one-element axis).

Parameters:: raw (Any)
Return type:: dict[str, list[Any]]

fastmdxplora.batch.sweep.normalize_systems(raw)[source]

Validate and normalize the systems list.

Each entry must be a mapping with at least system (the input). id is optional (defaults to a positional sN / derived name). Any phase-named keys (setup, simulation, …) are treated as per-system option overrides.

Parameters:: raw (Any)
Return type:: list[dict[str, Any]]

One table from a campaign’s members, and what their spread means.

A campaign leaves each member in its own directory with its own analyses. Reading them is arithmetic; knowing what the spread across them measures is not, and it is the whole of why this exists.

Replicas and variants are opposites, and they look identical on disk. Members that differ only by random seed are repeats of one measurement, so the spread of their means is the error on it. Members that differ by system, mutation or parameter are different measurements, so the spread between them is the result. Averaging the second kind reports a mutant series’ biology as noise; quoting the first kind as a finding reports noise as biology. This tells them apart from the campaign’s own sweep record rather than from the shape of the numbers, and says which it decided.

Replicas are what makes a reported uncertainty checkable. Every settled mean carries a standard error estimated from one trajectory’s autocorrelation. Ten trajectories give the same quantity a second way, as the spread of ten means, and the two should agree. Where the single-run estimate is much the smaller, the run was too short to see its own correlation time and the error bars it printed were too tight. That comparison is only available to a seed sweep, which is why it is reported only there.

fastmdxplora.batch.aggregate.SEED_AXES = frozenset({'random_seed', 'simulation.random_seed'}): Sweep axes that make members repeats rather than different studies. Written out rather than matched on the word “seed”, because a study sweeping something merely named seed-like is not replicating anything.

fastmdxplora.batch.aggregate.aggregate_members(batch_dir)[source]

Collect a campaign’s members into one comparison.

Returns the per-member values for every analysis that reported a settled mean, and, where the members are replicas, the comparison between the error each run estimated for itself and the spread the replicas actually show.

Parameters:: batch_dir (str | Path)
Return type:: dict[str, Any]

fastmdxplora.batch.aggregate.read_member_findings(run_dir)[source]

Every analysis’s findings in one member, keyed by analysis name.

Parameters:: run_dir (str | Path)
Return type:: dict[str, dict[str, Any]]

Cross-run comparison report for a sweep / multi-system study.

After a batch of runs completes, this module reads each run’s analysis outputs and produces a single comparison report at the batch root:

<batch_output>/comparison/
overlay_<analysis>.png # all runs’ curves on one axes trend_<analysis>.png # summary scalar vs swept parameter comparison_summary.csv # one row per run, summary scalars comparison_report.md # the written report tying it together

Two complementary views are produced:

Overlays — for per-frame analyses (RMSD, Rg, Q-value, total SASA), every run’s time series is drawn on one set of axes, labelled by its swept value, so divergence across the sweep is visible at a glance.

Trends — each run is reduced to a summary scalar (e.g. mean RMSD over the production trajectory) and plotted against the swept parameter, turning a directory of runs into a structure-property relationship.

The report degrades gracefully: analyses that didn’t run, runs that errored, and sweeps over non-numeric axes are handled without failing — the report simply includes what it can and notes the rest.

fastmdxplora.batch.compare.build_comparison_report(batch_output_dir)[source]

Build the cross-run comparison report for a completed batch.

Parameters:: batch_output_dir (path) – The batch root directory containing batch_manifest.json and a runs/ directory.
Returns:: The path to the comparison report directory, or None if there was nothing to compare (fewer than two successful runs, or no analysis outputs were found).
Return type:: Path or None

fastmdxplora.batch.compare.flabel_of(analysis)[source]

Summary-function label for an analysis (used in table headers).

Parameters:: analysis (str)
Return type:: str

Validation

Measurements of the software’s own behaviour, run against each release. Distinct from the test suite, which asserts: these report rates and comparisons against implementations that share no code with it.

Guardrails measured, rather than asserted.

Every guardrail in this software has a test showing it fires on a case chosen because it should fire. That is sensitivity, and sensitivity alone is not evidence: a checker that refused every study would score perfectly on it. The claim worth making is that the guardrails fire on defects and stay quiet on ordinary work, and the second half needs a corpus of ordinary work with the expected answer written down first.

So this runs two corpora and reports both rates. A case names what it does, what should happen, and why that is the right answer; the harness records what did happen and compares. Nothing here runs dynamics: the guardrails being measured decide before or after a trajectory, not during one, so the whole corpus completes in seconds and can run on every release rather than once before a paper.

Three outcomes are distinguished, because collapsing them would hide the distinction the software is built on. A case may proceed with a number, be refused with a reason, or be qualified – answered, but with a statement attached saying what the answer does not support. A qualification is not a refusal, and counting it as one would make the software look more obstructive than it is; counting it as a clean pass would make it look less careful.

fastmdxplora.validation.corpus.CLEAN: list[Case] = [Case(name='binding free energy from a run that reached bulk', run=<function _complete_pmf>, expect='proceeded', because="the tail follows a free ligand's shape, so the reference is a reference", mentioning=''), Case(name='order parameters on a settled peptide', run=<function _order_parameters_on_a_settled_run>, expect='proceeded', because='the halves agree, so the values are not qualified', mentioning=''), Case(name='radial distribution within the box', run=<function _rdf_within_the_box>, expect='proceeded', because='the whole requested range lies inside half the box', mentioning=''), Case(name='mutation named against the residue that is there', run=<function _mutation_that_matches>, expect='proceeded', because='the structure holds what the mutation says it holds', mentioning=''), Case(name='metadynamics that crossed and settled', run=<function _metadynamics_that_crossed_and_settled>, expect='proceeded', because='the system visited both basins repeatedly and the hills flattened', mentioning=''), Case(name='a run long against its correlation time', run=<function _a_run_long_against_its_correlation>, expect='proceeded', because='the halves agree, so the effective sample count means what it says', mentioning=''), Case(name='a box that already fits its cutoff', run=<function _a_box_that_already_fits>, expect='proceeded', because='no growing was needed, so nothing is reported about it', mentioning=''), Case(name='fraction of native contacts on a chain long enough to have a fold', run=<function _q_on_a_chain_long_enough>, expect='proceeded', because='a hairpin puts residues far apart in sequence within contact distance, which is what Q measures', mentioning=''), Case(name='a two-dimensional surface where both coordinates moved', run=<function _a_surface_where_both_coordinates_moved>, expect='proceeded', because='both variables visited both basins, so neither axis is the shape of the bias', mentioning=''), Case(name='density from a constant-pressure run', run=<function _density_at_constant_pressure>, expect='proceeded', because='the box breathes, so the density is measured rather than set', mentioning=''), Case(name='a save selection that keeps what was asked for', run=<function _a_selection_that_keeps_something>, expect='proceeded', because='the selection matches atoms, so there is nothing to cap or correct', mentioning=''), Case(name='a torsion with two states on the circle', run=<function _a_torsion_with_two_states_on_the_circle>, expect='proceeded', because='two rotamers half a turn apart are two states however the coordinate is read, so there is a barrier to report', mentioning=''), Case(name='crystallographic water retained on request', run=<function _crystallographic_water_retained_on_request>, expect='proceeded', because='the water model is part of the force field, so there is no chemistry to retrieve and nothing to refuse over', mentioning=''), Case(name='a lone ion kept without chemistry', run=<function _a_lone_ion_kept_without_chemistry>, expect='proceeded', because='a monatomic component has no bonds for an SDF to describe and the protein force field carries its parameters', mentioning=''), Case(name='umbrella windows that overlap', run=<function _umbrella_windows_that_overlap>, expect='proceeded', because='neighbours share ground, so the recombination rests on sampling rather than on interpolation across a gap', mentioning='')]: Ordinary studies, where nothing should fire. This is the half that makes the detection rate above a measurement rather than an assertion.

class fastmdxplora.validation.corpus.Case(name, run, expect, because, mentioning='')[source]

Bases: object

One study, and what the software should say about it.

Parameters:

name (str)
run (Callable[[], Any])
expect (str)
because (str)
mentioning (str)

because: str: Why that is the right answer. Recorded so a disagreement can be read without reconstructing the intent from the code.

expect: str: “refused”, “qualified” or “proceeded”.

mentioning: str = '': A fragment the refusal or qualification must contain, so a case cannot pass by refusing for an unrelated reason.

name: str

run: Callable[[], Any]: What the case exercises. Returns whatever the guardrail returns; raising is one of the ways a guardrail can refuse.

fastmdxplora.validation.corpus.DEFECTS: list[Case] = [Case(name='binding free energy from a run that never reached bulk', run=<function _truncated_pmf>, expect='refused', because='the reference the well depth is measured against is not a reference unless the ligand is free there', mentioning='-2kT ln r'), Case(name='order parameters on a structure without hydrogens', run=<function _order_parameters_without_hydrogens>, expect='refused', because='the measurement is of a bond vector and the bond is absent', mentioning='without hydrogens'), Case(name='radial distribution with no periodic box', run=<function _rdf_without_a_box>, expect='refused', because='there is no volume to take a bulk density from', mentioning='no unit cell'), Case(name='mutation named against the wrong residue', run=<function _mutation_against_the_wrong_residue>, expect='refused', because='applying it would simulate a protein nobody chose', mentioning='that position holds'), Case(name='fraction of native contacts on a chain too short to have a fold', run=<function _q_on_a_chain_too_short>, expect='refused', because='no residue pair is far enough apart in sequence to make a tertiary contact', mentioning='min_seq_separation'), Case(name='two-dimensional surface with one coordinate stuck', run=<function _surface_with_a_stuck_dimension>, expect='refused', because='the surface across a coordinate that did not move is the shape of the bias', mentioning='dist'), Case(name='radial distribution asked for beyond half the box', run=<function _rdf_past_half_the_box>, expect='qualified', because='the range is cut to where the shells are complete, and the cut is stated rather than silent', mentioning='half the smallest box'), Case(name='density from a constant-volume run', run=<function _density_at_constant_volume>, expect='qualified', because='the density is a constant the setup fixed, so a mean with an error on it would describe arithmetic', mentioning='constant the setup fixed'), Case(name='metadynamics stopped before any recrossing', run=<function _metadynamics_without_a_recrossing>, expect='refused', because='a barrier the system never crossed is the shape of the bias, not of the landscape', mentioning='cross'), Case(name='a run too short for its own correlation time', run=<function _a_run_too_short_for_its_own_correlation>, expect='qualified', because='the effective sample count is an upper bound, and the error bar printed from it is too tight', mentioning='not long'), Case(name='a box smaller than twice the cutoff', run=<function _a_box_too_small_for_its_cutoff>, expect='qualified', because='growing the padding without limit to reach an impossible cutoff would solvate forever, so it stops and says what it tried', mentioning='stopped'), Case(name='a save selection that would keep nothing', run=<function _a_selection_that_would_save_nothing>, expect='qualified', because='a box of pure water is a legitimate study, so `not water` matching none of it saves everything and names itself', mentioning='matched none'), Case(name='a save selection that will not parse', run=<function _a_malformed_save_selection>, expect='refused', because='a mistake in the study file, fixable in one edit, that would otherwise affect every frame', mentioning='not a selection')]: Studies with one named thing wrong. The expected answer was written before any of them was run.

class fastmdxplora.validation.corpus.Outcome(case: 'str', expected: 'str', observed: 'str', agreed: 'bool', detail: 'str' = '', because: 'str' = '')[source]

Bases: object

Parameters:

case (str)
expected (str)
observed (str)
agreed (bool)
detail (str)
because (str)

agreed: bool

because: str = ''

case: str

detail: str = ''

expected: str

observed: str

fastmdxplora.validation.corpus.run_case(case)[source]

Run one case and say whether the software agreed with the record.

Parameters:: case (Case)
Return type:: Outcome

fastmdxplora.validation.corpus.run_corpus(defects, clean)[source]

Both corpora, with the two rates that matter reported together.

The detection rate alone is what every guardrail test in this repository already shows. The false-refusal rate is the half that makes it a measurement.

Parameters:

defects (list[Case])
clean (list[Case])

Return type: