EnsembleQL

EnsembleQL is an open-source, declarative temporal query engine for molecular-dynamics trajectories. It lets computational biologists and chemists ask when molecular behavior occurs and how events relate in time, instead of rebuilding each analysis as a bespoke array-processing script.

trajectory -> observables -> predicates -> events -> temporal relationships

The first milestone focuses on transient contacts in intrinsically disordered proteins (IDRs). Its execution engine is C++20; pybind11 exposes a compact Python and CLI interface.

Why events, not only averages?

Frame-wise contact probabilities can erase order. A conventional analysis might say:

R17-D42 = 32%
R17-E53 = 29%

Those values cannot distinguish independent contacts from a directed interaction switch. EnsembleQL can instead report:

R17-D42 -> R17-E53 switching events: 14
median transition gap: 0.8 ns

Events are maximal contiguous intervals over which a frame predicate is true. Temporal operators act on those intervals independently of the molecular observable that produced them.

Quick start

Install a released Linux or macOS wheel from PyPI:

python -m pip install ensembleql

Release wheels include the Chemfiles backend. The commands below build from a source checkout.

Build and test the dependency-free native core:

cmake -S . -B build -DENSEMBLEQL_BUILD_PYTHON=OFF
cmake --build build -j
ctest --test-dir build --output-on-failure

Enable additional trajectory formats—including XTC, TRR, DCD, Amber NetCDF, GRO, LAMMPS trajectories, and TNG—through an installed chemfiles library, or fetch the pinned stable release during configuration:

cmake -S . -B build -DENSEMBLEQL_FETCH_CHEMFILES=ON
cmake --build build -j

Use -DENSEMBLEQL_REQUIRE_CHEMFILES=ON when configuration should fail rather than produce a built-in-only build. Python exposes eql.chemfiles_backend_available() and eql.supported_trajectory_extensions() for capability checks.

To include the backend in an editable Python installation:

CMAKE_ARGS="-DENSEMBLEQL_FETCH_CHEMFILES=ON -DENSEMBLEQL_REQUIRE_CHEMFILES=ON" \
  python -m pip install -e .

For the Python package (Python 3.11+), install into an isolated environment. The build installs its pybind11/scikit-build dependencies:

python -m venv .venv
source .venv/bin/activate
python -m pip install -e .
pytest

Run the included IDR contact-switching example:

from pathlib import Path
import ensembleql as eql

root = Path("examples/idr_contact_switching")
traj = eql.load(root / "switching.xyz", topology=root / "switching.pdb")
events = traj.query("""
    FIND CONTACT(resid 17, resid 42)
    FOLLOWED_BY CONTACT(resid 17, resid 53)
    WITHIN 3ns;
""")
print(events)

Two end-to-end research fixtures are included: IDR contact switching and a peptide–surface adsorption mechanism. Both are intentionally small enough to audit frame by frame; they demonstrate query semantics rather than supply physical cutoff recommendations.

PDB files can be queried directly: a file without MODEL records is one frame, while each MODEL block in an ensemble is a frame. When a trajectory has no timestamps, provide an explicit fallback spacing. Embedded timestamps always take precedence:

traj = eql.load("trajectory.xyz", topology="structure.pdb", default_timestep="2fs")

The synthetic trajectory includes an R17-D42 contact, a one-frame boundary where both contacts exist, then an R17-E53 contact. EnsembleQL returns one directed switch from 0 to 4000 ps with a zero-ps transition gap. This preserves timing and direction that two contact probabilities do not.

DSL

Implemented queries include:

FIND CONTACT(resid 17, resid 42);
FIND CONTACT(resid 17, resid 42, cutoff=0.35nm);
FIND CONTACT(resid 17, resid 42) FOR >= 5ns;
FIND DISTANCE(resid 17, resid 42) < 0.8nm;
FIND RG(protein) < 2.0nm FOR >= 10ns;
FIND CONTACT_COUNT(resid 1:20, resid 40:60) >= 4;
FIND RG(protein, mass_weighted=true) < 2.0nm;
FIND CONTACT_COUNT(resid 1:20, resid 40:60, mode=residue) >= 4;
FIND CONTACT(resid 17, resid 42) OVERLAPS CONTACT(resid 17, resid 53);
FIND HBOND(name N, name O, distance=0.35nm, min_angle=150deg);
FIND DIHEDRAL(name C1, name N2, name CA2, name C2) < -30deg;
FIND HELIX(resid 20:32, minimum_fraction=0.7);
FIND SASA(protein, probe=0.14nm, points=192) < 40nm2;
FIND RMSD(protein) < 0.2nm;
FIND COORDINATION_NUMBER(resname ZN, name O, cutoff=0.3nm) >= 4;
FIND SALT_BRIDGE(resname ARG and name CZ, resname ASP and name CG);
FIND AROMATIC_STACKING(resid 10, resid 25);
FIND SURFACE_DISTANCE(protein, resname SUR) < 0.4nm;
FIND ORIENTATION(resid 1, resid 10, axis=z) < 30deg;
FIND CONTACT(resid 17, resid 42) REPEATS >= 3 WITHIN 20ns;

FOLLOWED_BY, WITHIN, IMMEDIATELY_FOLLOWED_BY, BEFORE, PRECEDES, AFTER, OVERLAPS, DURING, UNTIL, REPEATS, FOR, AND, and OR have AST nodes. Parentheses can group temporal expressions. The parser only creates an AST; the planner resolves selections and deduplicates required observables; the engine then evaluates every required predicate in a single streaming traversal and retains events rather than a full boolean time series.

Selections support resid 17, resid 17:25, name CA, resname ARG, chain A, protein, and hydrogen. They compose with case-insensitive and, or, not, and parentheses; and binds more tightly than or.

CONTACT and CONTACT_COUNT use a 0.45 nm default cutoff. Override it with a dimensionally checked option such as cutoff=4A or cutoff=0.35nm. CONTACT_COUNT accepts mode=atom (the default) or mode=residue; options can appear in either order.

CLI

ensembleql query \
  --topology examples/idr_contact_switching/switching.pdb \
  --trajectory examples/idr_contact_switching/switching.xyz \
  --file examples/idr_contact_switching/query.eql

ensembleql query --topology structure.pdb --trajectory trajectory.xyz \
  --default-timestep 2fs \
  --query "FIND CONTACT(resid 17, resid 42) FOR >= 2ns;" --format json

Output formats are table, csv, and json. EventResults.to_dataframe() returns a pandas DataFrame when pandas is installed and otherwise returns a list of records. EventResults also provides recurrence statistics, explicitly normalized event frequencies, conditional probabilities, transition matrices, motif and recurring-subsequence counts, temporal clusters, event graphs, and aggregated state-transition networks.

Inspect a query before reading trajectory frames:

ensembleql explain \
  --topology examples/idr_contact_switching/switching.pdb \
  --file examples/idr_contact_switching/query.eql

The explanation reports resolved selection expressions, canonical deduplicated observables, frame predicates, temporal operations, and the execution-plan tree. Use --format json for machine-readable output. Python provides the same information through traj.explain(query) or eql.explain(query, topology="structure.pdb").

Each Python Trajectory retains parsed and topology-resolved plans by exact query text. Repeated query() and explain() calls reuse them. Inspect traj.cached_plan_count or call traj.clear_plan_cache() when managing a long-lived interactive session.

Scientific definitions and units

• Coordinates and distances are normalized to nm. XYZ coordinates are interpreted as angstroms; PDB is used for topology metadata only.
• Times are normalized to ps. Supported distance units are nm, angstrom, and A; supported time units are fs, ps, ns, and us. Query thresholds and configurable fallback timesteps require explicit units except integer-like counts.
• DISTANCE(A,B) is the minimum distance between distinct atoms in A and B. It uses Euclidean distance without a cell and a nearest-image search for orthorhombic or triclinic periodic cells.
• CONTACT(A,B) is true when that minimum distance is less than or equal to 0.45 nm. CONTACT_COUNT counts unique unordered atom pairs, or unique unordered (chain, resid) pairs with mode=residue, at or below the same inclusive cutoff. Both honor orthorhombic and triclinic periodic boundaries.
• RG(A) is the unweighted root-mean-square distance of selected atom coordinates from their geometric centroid. mass_weighted=true uses standard atomic weights derived from PDB element symbols. On periodic frames, both variants reconstruct the selected atoms by traversing PDB CONECT bonds before calculating the centroid. See observable definitions.
• Event endpoints are the timestamps of the first and last true sampled frames; duration is end - start. Single-sample events therefore have zero observed duration. Irregular timestamps are supported, while duplicate or decreasing timestamps are rejected. FOR >=, contact cutoffs, WITHIN, and interval boundary comparisons are inclusive. See sampled-time event semantics.
• Periodic cells use a nearest-image search over neighboring lattice translations. XYZ comments accept box=20,20,20A; extended XYZ accepts a full Lattice="..." matrix in angstroms. Periodic RG requires the selected atoms to belong to one connected PDB bond component. See periodic-boundary conventions.

Architecture

Public headers separate trajectory/topology I/O, selections, geometry, observables, event extraction, interval algebra, AST parsing, planning, and execution. FrameReader is the backend-neutral streaming interface used by the built-in XYZ/PDB readers and the optional multi-format Chemfiles adapter. Observable classes are likewise independent of the parser.

Contact enumeration uses spatial hashing for sufficiently large non-periodic, orthorhombic, and triclinic selections, while small selections use the reference pairwise kernel. Stateful contact observables reuse Verlet-style candidate lists across frames and rebuild after a half-skin displacement or cell change. With OpenMP available, Cartesian candidate-distance filtering uses an explicit SIMD loop and SASA distributes selected atoms across threads; configure with -DENSEMBLEQL_ENABLE_OPENMP=OFF to force serial kernels. Boolean CONTACT evaluation retains scalar early exit. Other extension points include memory-mapped trajectory readers and additional observables. See the roadmap.

Current limitations

PDB supplies topology metadata, element-derived standard atomic weights, and explicit CONECT bonds. XYZ and single- or multi-model PDB trajectories are always supported; the optional chemfiles backend adds the formats documented in trajectory backends. PDB CRYST1 records preserve orthorhombic or triclinic cells. Chemfiles coordinates and cell vectors are converted from angstroms to nm; numeric time properties are interpreted as ps, while formats without one use the configured fallback. Periodic RG requires connected bond metadata, and EnsembleQL does not yet infer standard-residue bonds. Event intervals use sampled timestamps and therefore do not infer behavior between frames. AND/OR combine frame predicates; temporal relations operate on extracted intervals.

Development

The native target compiles with -Wall -Wextra -Wpedantic; CI additionally enables ENSEMBLEQL_WARNINGS_AS_ERRORS. Enable microbenchmarks with -DENSEMBLEQL_BUILD_BENCHMARKS=ON. The synthetic benchmarks compare optimized and pairwise contact detection at increasing atom counts, exercise neighbor-list reuse and triclinic hashing, report boolean-contact early-exit time, and cover streaming event extraction and temporal joins.

The scientific-validation suite compares every observable family with pinned MDAnalysis, MDTraj, or NumPy calculations on published AdK and membrane/peptide trajectories. Its machine-readable manifests record provenance, checksums, thresholds, tolerances, and expected event intervals; the large generated fixtures remain outside Git.

Maintainer release builds use tested Linux and macOS wheels, Trusted Publishing, and tag/version consistency checks. See the contribution guide, release procedure, and changelog. Research users can cite the software using CITATION.cff.

Contributions should preserve scientific definitions, add boundary-condition tests, and keep file-format backends independent from the engine.