topo.utils package¶

Submodules¶

topo.utils.config module¶

Read a simulation control file (md.ini) into a structured configuration.

The control file is an INI file: a flat list of key = value lines. Parsing it used to be a ~100-line block copy-pasted into every run_simulation.py. This module centralizes that logic so a run script can simply do:

import topo

cfg = topo.read_simulation_config("md.ini")
model = topo.models.buildCoarseGrainModel(cfg.pdb_file, **cfg.build_kwargs())
integrator = openmm.LangevinIntegrator(cfg.ref_t, cfg.tau_t, cfg.dt)

All parsed values are returned on a SimulationConfig dataclass, with OpenMM units already applied where appropriate (dt, ref_t, tau_t, ref_p).

class topo.utils.config.SimulationConfig(md_steps: int = 10000, dt: ~typing.Any = <factory>, nstxout: int = 5000, nstlog: int = 5000, nstchk: int | None = 5000, nstcomm: int | None = None, log_precision: int | None = 4, log_width: int | None = 14, model: str = 'topo', constraints: ~typing.Any = 'AllBonds', constraint_tolerance: float = 1e-05, tcoupl: bool = True, ref_t: ~typing.Any = 300.0, tau_t: ~typing.Any | None = None, anneal: bool = False, t_high: ~typing.Any | None = None, anneal_steps: int = 0, anneal_ramp: str = 'jump', anneal_ramp_steps: int = 0, anneal_ramp_increments: int = 20, pcoupl: bool = False, ref_p: ~typing.Any = 1.0, frequency_p: int = 25, pbc: bool = False, box_dimension: ~typing.List[float] | None = None, pdb_file: str | None = None, init_position: str | None = None, domain_def: str | None = None, stride_output_file: str | None = None, output_dir: str = 'traj', outname: str = 'traj', checkpoint: str | None = None, device: str = 'CPU', ppn: int = 1, n_copies: int = 1, copy_shift: float = 5.0, restart: bool = False, minimize: bool = True, config_file: str | None = None, quiet: bool = False)[source]¶

Bases: object

Parsed contents of a simulation control file (md.ini).

anneal: bool = False¶

anneal_ramp: str = 'jump'¶

anneal_ramp_increments: int = 20¶

anneal_ramp_steps: int = 0¶

anneal_steps: int = 0¶

box_dimension: List[float] | None = None¶

build_kwargs() → dict[source]¶

Keyword arguments for topo.models.buildCoarseGrainModel().

Always passes minimize, model and box_dimension; passes domain_def / stride_output_file only when they are set, so the builder’s own defaults apply otherwise. On restart the build-time energy check is skipped (check_forces=False) because the loaded checkpoint state, not the input structure, is what gets simulated.

checkpoint: str | None = None¶

checkpoint_path() → str[source]¶: Resolved checkpoint path: the explicit checkpoint option if given, otherwise <output_dir>/<outname>.chk.

config_file: str | None = None¶

constraint_tolerance: float = 1e-05¶

constraints: Any = 'AllBonds'¶

copy_shift: float = 5.0¶

device: str = 'CPU'¶

domain_def: str | None = None¶

dt: Any¶

frequency_p: int = 25¶

init_position: str | None = None¶

log_precision: int | None = 4¶

log_width: int | None = 14¶

make_platform()[source]¶

Build the OpenMM (platform, properties) pair selected by device.

device = 'GPU' -> CUDA (mixed precision, device 0); anything else -> CPU using ppn threads.

md_steps: int = 10000¶

minimize: bool = True¶

model: str = 'topo'¶

n_copies: int = 1¶

nstchk: int | None = 5000¶

nstcomm: int | None = None¶

nstlog: int = 5000¶

nstxout: int = 5000¶

outname: str = 'traj'¶

output_dir: str = 'traj'¶

output_path(suffix: str = '') → str[source]¶

Path for a generated output file: <output_dir>/<outname><suffix>.

Examples: output_path('.dcd') -> traj/traj.dcd; output_path('_multi.psf') -> traj/traj_multi.psf.

Built with pathlib.Path but returned as str so it can be passed directly to OpenMM/parmed writers (some of which special-case str and would mishandle a raw Path).

pbc: bool = False¶

pcoupl: bool = False¶

pdb_file: str | None = None¶

ppn: int = 1¶

prepare_output_dir() → None[source]¶: Ensure output_dir exists. Path.mkdir(parents=True, exist_ok=True) creates any missing parents and is a no-op if the folder already exists (no manual “does it exist?” check needed).

quench_steps() → int[source]¶

Total steps in the quench phase (0 when anneal is off).

The quench is a separate phase from production: it holds at t_high for anneal_steps and, for a linear ramp, additionally cools over anneal_ramp_steps. Production then runs md_steps at ref_t.

quiet: bool = False¶

ref_p: Any = 1.0¶

ref_t: Any = 300.0¶

restart: bool = False¶

stride_output_file: str | None = None¶

t_high: Any = None¶

tau_t: Any = None¶

tcoupl: bool = True¶

total_steps() → int[source]¶: Grand total integration steps across the quench + production phases.

topo.utils.config._unshift_parsing_error(exc: ParsingError, config_file: str) → ParsingError[source]¶

Rebuild a configparser.ParsingError with honest line numbers.

Parameters:

exc (configparser.ParsingError) – The error raised against the shifted text.
config_file (str) – Path to report as the source.

Returns:

The same errors, renumbered to the user’s file.

Return type:

configparser.ParsingError

topo.utils.config.read_ini(config_file: str, preserve_case: bool = False) → ConfigParser[source]¶

Parse a topo control file into a configparser.ConfigParser.

A topo control file is a flat list of key = value lines. Section headers are optional and their names carry no meaning: every setting in the file is read into one flat set, whether it is written under one header, several, or none at all. So no setting can go unread because of where it sits in the file – a key that is written is a key that is applied.

Parameters:

config_file (str) – Path to the control file (md.ini, optimize.ini, csp.ini, …). Inline comments starting with # or ; are ignored.
preserve_case (bool, optional (default: False)) – If True, keep option names as written instead of lower-casing them (configparser’s default). Needed when the keys are round-tripped back out to another file, as in topo.optimize.

Returns:

Every setting in the file, collected into one section. The section is always present, and is empty only for a file that sets nothing.

Return type:

configparser.ConfigParser

Raises:

FileNotFoundError – If config_file does not exist.
ValueError – If a key is set more than once. There is no principled way to pick a winner, and silently applying one of two values would change a run’s physics without saying so, so a repeated key is an error rather than a last-one-wins.
configparser.ParsingError – If a line is neither a comment nor a key = value pair.

topo.utils.config.read_simulation_config(config_file: str, verbose: bool = True) → SimulationConfig[source]¶

Parse a simulation control file into a SimulationConfig.

Parameters:

config_file (str) – Path to the control file (e.g. md.ini): a flat key = value list (see read_ini()). Inline comments starting with # or ; are ignored. Underscores in md_steps (e.g. 500_000) are allowed.
verbose (bool, optional (default: True)) – If True, echo each parsed setting to stdout (the original behaviour of the inline parser). Set False for a quiet read.

Returns:

All settings with OpenMM units applied where relevant. Options absent from the file fall back to the dataclass defaults.

Return type:

SimulationConfig

Notes

Behavioural details preserved from the original inline parser:

ref_t/tau_t get units only when tcoupl is on.
box_dimension accepts a scalar L (cubic [L, L, L]) or an explicit [x, y, z] list, and is None when pbc is off.
pcoupl requires pbc (asserted).
restart forces minimize off.

topo.utils.config.strtobool(value)[source]¶

Convert a truthy/falsy string to 0/1.

Local replacement for distutils.util.strtobool (distutils was removed in Python 3.12). Accepts the same vocabulary: true/yes/on/1 -> 1, false/no/off/0 -> 0 (case-insensitive); anything else raises ValueError.

topo.utils.multichain module¶

Replicate a single-chain coarse-grained system into many non-interacting copies inside one OpenMM System.

Coarse-grained protein models are tiny (a few hundred beads), so a single chain badly underuses a GPU. Packing N independent copies into one simulation fills the device and yields N trajectories per run for better sampling. The copies must not interact: bonded terms are duplicated per copy, and every CustomNonbondedForce is restricted to intra-copy interactions with addInteractionGroup so copy i never sees copy j.

Typical use:

import topo
cg = topo.models.buildCoarseGrainModel("protein.pdb", domain_def="domain.yaml")
system, topology, positions = topo.make_noninteracting_copies(
    cg.system, cg.topology, cg.positions, n_copies=20)
# ... build a Simulation from system/topology/positions, run, then split the
# multi-chain DCD into per-chain trajectories for analysis.

After the run, the atoms of copy k are the contiguous block [k*n_atoms : (k+1)*n_atoms] (see split_indices()), which makes splitting the trajectory trivial.

topo.utils.multichain._split_cli()[source]¶: CLI: python -m topo.utils.multichain — split a combined multi-copy DCD.

topo.utils.multichain.make_noninteracting_copies(system: ~openmm.openmm.System, topology: ~openmm.app.topology.Topology, positions, n_copies: int, shift=5.0 nm)[source]¶

Convenience wrapper: replicate system, topology and positions into n_copies non-interacting copies in one call.

Return type:: (full_system, full_topology, full_positions)

topo.utils.multichain.replicate_positions(template_positions, n_copies: int, shift=5.0 nm)[source]¶: Replicate positions n_copies times, translating each copy along x by shift so the chains start in separate regions of space (they do not interact, but separation keeps the start clean and aids visualization).

topo.utils.multichain.replicate_system_intra_only(template_system: System, n_copies: int) → System[source]¶

Replicate template_system into n_copies independent copies.

Particles and constraints are duplicated per copy; each supported force is rebuilt with per-copy atom-index offsets. CustomNonbondedForce objects get one interaction group per copy so copies never interact with one another.

Supported forces: HarmonicBondForce, CustomAngleForce, PeriodicTorsionForce, CustomTorsionForce, CustomNonbondedForce. CMMotionRemover is skipped (add a fresh one to the full system if wanted); any other force type raises NotImplementedError.

topo.utils.multichain.replicate_topology(template_topology: Topology, n_copies: int) → Topology[source]¶: Replicate an OpenMM Topology into n_copies chains.

topo.utils.multichain.split_chains(combined_dcd, output_paths, *, chunk_size: int = 1000, center: bool = True)[source]¶

Split a combined multi-copy DCD into one DCD per copy.

The combined trajectory’s atoms are assumed to be n_copies contiguous equal-size blocks (the layout produced by make_noninteracting_copies() / split_indices()), so copy k is atoms [k*m : (k+1)*m].

Memory-bounded streaming. The combined DCD is read in chunks of chunk_size frames and each chunk is written straight to the per-copy files, so peak memory is chunk_size * n_atoms * 3 floats regardless of how long the trajectory is — a combined DCD too large to fit in memory is handled fine. A single pass over the input writes all copies.

Coordinates stay in Å throughout (mdtraj’s low-level DCD read and write are both in Å), so no unit conversion is needed or applied.

Parameters:

combined_dcd (str or Path) – The combined multi-copy DCD.
output_paths (sequence of str or Path) – One output DCD path per copy; n_copies = len(output_paths). Parent directories are created as needed.
chunk_size (int, optional) – Frames read/written per chunk (default 1000).
center (bool, optional) – If True (default), translate each copy to its own centre of geometry per frame, so every per-copy trajectory is centred at the origin and the inter-copy copy_shift offset is removed. Set False to preserve the raw coordinates exactly.

Returns:

(n_copies, atoms_per_copy, n_frames)

Return type:

tuple[int, int, int]

topo.utils.multichain.split_indices(n_atoms_single: int, n_copies: int) → List[Tuple[int, int]][source]¶: Return [(start, stop), …] atom-index ranges, one per copy.

topo.utils.nonbonded module¶

Build non-bonded interaction matrices for the TOPO coarse-grained model.

This module constructs distance and energy matrices used by OpenMM for structure-based (native) contacts and repulsive non-native contacts. Inputs: PDB structure, STRIDE hydrogen-bond output, and domain definition YAML.

Main workflow¶

Hydrogen bonds (STRIDE): parsed from stride output; 1 H-bond -> 0.75 kcal/mol, 2+ H-bonds capped at 1.5 kcal/mol per pair.
Backbone-sidechain (BS) and sidechain-sidechain (SS) contacts: distance-based (default cutoff 4.5 Angstrom); BS energy 0.37 kcal/mol; SS energies from BT potential and domain scaling.
Domain scaling: YAML defines intra_domains (residue lists + nscale) and optional inter_domains (pair nscales). nscale is the native-contact scaling factor applied to the sidechain-sidechain well depths (not an absolute energy – the contacts still interact at nscale = 1.0). Single-domain proteins can omit inter_domains.
Non-native contacts: repulsive well with per-residue Rmin/2 from the nearest non-contact CA-CA distance (sum rule); energy NON_NATIVE_KJ.

Typical use¶

>>> from topo.utils.nonbonded import build_nonbonded_interaction
>>> rmin_matrix, energy_matrix = build_nonbonded_interaction(
...     "protein.pdb", "domain.yaml", "stride.dat"
... )
>>> # rmin_matrix (well positions), energy_matrix in nm and kJ/mol for OpenMM

topo.utils.nonbonded._norm_chain(chain) → str[source]¶

Normalize a chain identifier to a canonical string.

Residues are keyed by (chain, resid) so that structures with multiple chains (where the same residue number can appear in more than one chain) are handled correctly. Different sources spell a blank chain differently: MDAnalysis reports '' or a space, while STRIDE prints '-'. All of these are normalized to the empty string so the two sources agree.

topo.utils.nonbonded._residue_chain(res) → str[source]¶

Return the normalized chain ID for an MDAnalysis residue.

Prefers the PDB chainID (matching STRIDE’s chain column); falls back to segid if chainID is unavailable.

topo.utils.nonbonded._scaling_matrix_from_config(domain_to_residues: Dict, intra_nscales: Dict, inter_nscales: Dict) → ndarray[source]¶

Build the SS-energy scaling matrix from already-parsed domain config.

Split out of get_scaling_ss_matrix() so build_nonbonded_interaction() can parse the domain YAML once (for both scaling and the disorder section) rather than reading it twice. See get_scaling_ss_matrix() for the returned matrix.

topo.utils.nonbonded.apply_disorder(rmin_matrix: ndarray, eps_ij: ndarray, rmin_2_nm: ndarray, dis_idx: ndarray, idr_scale: float, index_to_resname: Dict[int, str], ss_interaction_energy: ndarray, eps_gen_kj: float = 0.0)[source]¶

Transform a folded non-bonded build into an IDR-aware one (spec §2.3).

Runs at the very end of build_nonbonded_interaction(), after the folded build is complete and everything is in nm. It overwrites – in place – only the entries that touch a disordered residue, leaving folded-folded pairs and every folded residue’s Karanicolas-Brooks radius byte-identical to a folded-only run.

Three pair classes result (spec §2.1): folded-folded (untouched, Go), IDR-IDR (max(NON_NATIVE_KJ, eps_gen_kj + idr_scale * eps_BT)), and IDR-folded (excluded-volume only). eps_gen_kj is an optional additive, sequence-independent generic-cohesion term (0 by default -> the original idr_scale-only well). Each IDR residue’s radius is switched to the transferable per-AA sigma-radius rvdw = IDR_RMIN_2_NM / 2^(1/6) (folded residues keep their K-B Rmin/2), and every IDR-involving well is placed at the bare sum of those per-bead radii. So IDR-IDR comes out rvdw[i] + rvdw[j] – the O’Brien generic-bt convention, ~11% tighter than the bare Rmin sum – while IDR-folded is Rmin/2[folded] + rvdw[IDR], keeping the folded bead on its native radius rather than shrinking it in cross pairs.

Overriding the per-residue rmin_2_nm array (not merely the pair matrix) is a requirement: the same array feeds the NC<->ribosome excluded volume in CSP (spec §2.6), so both excluded-volume channels see the same rvdw for IDR beads.

Parameters:

rmin_matrix (np.ndarray, shape (n, n)) – Folded pairwise well positions (nm). Modified in place.
eps_ij (np.ndarray, shape (n, n)) – Folded pairwise well depths (kJ/mol). Modified in place.
rmin_2_nm (np.ndarray, shape (n,)) – Per-residue radius (nm). Modified in place: IDR entries are overridden with the per-AA sigma-radius rvdw = IDR_RMIN_2_NM / 2^(1/6); folded entries keep K-B Rmin/2. (So the array mixes conventions by bead class, both feeding a bare sum.)
dis_idx (np.ndarray of int) – 0-based indices of the disordered residues.
idr_scale (float) – IDR-IDR attraction scale (0 -> self-avoiding chain).
index_to_resname (dict) – 0-based index -> three-letter residue name (from get_residue_mapping()).
ss_interaction_energy (np.ndarray, shape (n, n)) – The RAW per-pair BT energy matrix (kJ/mol) from get_ss_interaction_energy() (abs + shift + kcal->kJ already applied); NOT the domain-scaled version – the IDR well must not inherit nscale.
eps_gen_kj (float, optional) – Additive, sequence-independent generic-cohesion depth (kJ/mol) added to every IDR-IDR well before the excluded-volume floor. Defaults to 0.0 at this level (the neutral, idr_scale-only transform); a domain_def that omits eps_gen_kj: gets the calibrated DEFAULT_EPS_GEN_KJ from read_yaml_config() instead, so user-facing runs default to generic cohesion being ON.

Returns:

rmin_matrix, eps_ij, rmin_2_nm – The same three arrays, now IDR-aware.

Return type:

np.ndarray

topo.utils.nonbonded.build_hb_contact_matrix(hb_pairs: List[Tuple], key_to_index: Dict[Tuple[str, int], int], n_residues: int) → ndarray[source]¶

Build the hydrogen-bond contact matrix from a list of H-bond pairs.

Counts how many H-bonds exist between each residue pair and stores 0, 1, or 2 (values greater than 2 are capped at 2). The model uses 0.75 kcal/mol per single H-bond and 1.5 kcal/mol total for multiple H-bonds; capping at 2 enforces that.

Parameters:

hb_pairs (list of tuple) – List of canonical pairs from parse_hydrogen_bonds(). Each pair is ((chain1, resid1), (chain2, resid2)).
key_to_index (dict) – Maps (chain, resid) -> 0-based residue index, from get_residue_mapping(). Used to place each H-bond at the correct matrix index regardless of chain (multi-chain safe).
n_residues (int) – Number of residues (determines matrix shape).

Returns:

Symmetric matrix. Entry [i, j] is 0, 1, or 2 (number of H-bonds between residue index i and j, capped at 2). Indices are the 0-based residue order from get_residue_mapping().

Return type:

np.ndarray, shape (n_residues, n_residues), dtype int

Notes

H-bond (chain, resid) keys are resolved to matrix indices via key_to_index. A pair whose residue is not found in the mapping (e.g. a hetero residue STRIDE reported but that is not in the protein selection) is skipped with a warning rather than mis-placed.

topo.utils.nonbonded.build_nonbonded_interaction(pdb_file: str, domain_def: str | None = None, stride_output_file: str | None = None, return_rmin_2: bool = False)[source]¶

Build the well-position (Rmin) and energy matrices for TOPO non-bonded contacts.

Combines hydrogen bonds (STRIDE), backbone–sidechain and sidechain–sidechain contacts (distance cutoffs), and domain-based scaling. Native contacts get their well at the measured CA–CA distance with energies from H-bond, BS, and scaled SS terms; non-native pairs get their well at the sum-rule collision distance Rmin/2_i + Rmin/2_j and a small repulsive energy (non_native). Output matrices are (n_residues × n_residues), symmetric, in nm and kJ/mol for use with OpenMM. With return_rmin_2=True also returns the per-residue Rmin/2 array (nm).

Parameters:

pdb_file (str) – Path to PDB structure (all-atom; used for residue list and CA distances).
domain_def (str, optional) – Path to domain YAML (intra_domains, n_residues, optional inter_domains). If None, all residues are treated as a single domain (scaling 1.0 everywhere).
stride_output_file (str, optional) – Path to STRIDE output for hydrogen bond list. If None, the function checks whether the stride program is available; if yes, it runs stride -h pdb_file, caches the output to {prefix}_stride.dat (prefix taken from the PDB filename), and parses that file. STRIDE success is judged by output content (presence of DNR records), not the process exit code, since some STRIDE builds return non-zero even on success. If stride is not found, raises an error (provide a precomputed STRIDE file or install STRIDE).

Returns:

rmin_matrix (np.ndarray, shape (n_residues, n_residues)) – Pairwise well position Rmin (nm). Native: CA–CA distance; non-native: Rmin/2_i + Rmin/2_j (sum rule).
energy_matrix (np.ndarray, shape (n_residues, n_residues)) – Pairwise well depth in kJ/mol. Native: sum of H-bond (0.75/1.5 kcal/mol), backbone–sidechain (0.37 kcal/mol), and scaled SS; non-native: NON_NATIVE_KJ.
rmin_2_nm (np.ndarray, shape (n_residues,)) – Only when return_rmin_2=True (the return is then a 3-tuple): the per-residue Rmin/2 array (nm) the sum rule was built from.

Raises:

FileNotFoundError – If pdb_file (or stride_output_file when given) cannot be opened.
RuntimeError – If stride_output_file is None and the stride program is not found or not executable.

Example

>>> dist, energy = build_nonbonded_interaction("2ww4.pdb", "domain.yaml", "stride.dat")
>>> dist, energy = build_nonbonded_interaction("2ww4.pdb")  # single domain, run stride if available
>>> dist.shape
(283, 283)

topo.utils.nonbonded.calculate_rmin_2_values(binary_contact_matrix: ndarray, ca_distances: ndarray, n_residues: int) → List[float][source]¶

Compute per-residue collision radius Rmin/2 for non-native contacts.

For residue i, Rmin/2[i] = 0.5 * RMIN_SCALE_FACTOR * (minimum CA-CA distance to residues that are (1) not in contact with i and (2) not within LOCAL_SEPARATION in sequence). The minimum CA-CA distance is taken as the LJ sigma (closest approach); RMIN_SCALE_FACTOR = 2^(1/6) scales it to the full collision diameter Rmin, and the 0.5 gives Rmin/2 – the per-residue collision radius (structure-derived Karanicolas-Brooks). If no such residue j exists, Rmin/2[i] = 0.0. Non-native pairs combine by the SUM rule Rmin_ij = Rmin/2[i] + Rmin/2[j].

Parameters:

binary_contact_matrix (np.ndarray, shape (n_residues, n_residues)) – Binary contact matrix (1 = in contact, 0 = not).
ca_distances (np.ndarray, shape (n_residues, n_residues)) – CA–CA distances (same units as used later; typically Angstrom).
n_residues (int) – Number of residues.

Returns:

Rmin/2[i] per residue i; non-native well position via the sum rule rmin_matrix[i, j] = Rmin/2[i] + Rmin/2[j].

Return type:

list of float

topo.utils.nonbonded.format_residue_ranges(nums: List[int]) → str[source]¶

Compress a list of residue numbers into a compact [a-b, c, d-e] string.

The inverse of parse_residue_list(): consecutive runs collapse to start-end so a long contiguous block prints as one range instead of thousands of numbers (e.g. an IDR covering residues 18..164 prints [18-164]). Duplicates are removed and the output is sorted.

Parameters:: nums (list of int) – Residue numbers (any order, duplicates allowed).
Returns:: e.g. "[18-164, 177-219, 279-363]"; "[]" for an empty list.
Return type:: str

Example

>>> format_residue_ranges([3, 1, 2, 5, 6, 7, 10])
'[1-3, 5-7, 10]'

topo.utils.nonbonded.get_bs_contact_matrix(u: Universe, cutoff: float = 4.5) → ndarray[source]¶

Build backbone–sidechain contact matrix from structure (distance-based).

Two residues are in contact if any backbone atom of one is within cutoff of any sidechain atom of the other (and vice versa), and they are not within LOCAL_SEPARATION along the sequence. The matrix is symmetric and stores the total number of directional contacts (i → j and j → i).

Parameters:

u (mda.Universe) – MDAnalysis universe with protein structure (must have backbone/sidechain).
cutoff (float, optional) – Distance cutoff in Angstrom (default DEFAULT_CUTOFF = 4.5).

Returns:

Symmetric matrix. Entry [i, j] is the number of directional backbone–sidechain contacts between residue index i and j (0, 1, or 2). Indices are 0-based residue order from get_residue_mapping().

Return type:

np.ndarray, shape (n_residues, n_residues), dtype int

Notes

Backbone: protein backbone atoms excluding H. Sidechain: protein non-backbone excluding H. Pairs within abs(resid_i - resid_j) <= LOCAL_SEPARATION are excluded, but only when both residues are in the same chain; residues in different chains are never excluded by sequence separation (multi-chain safe).

Why the count is 0, 1, or 2 (and not just a 0/1 flag). A BS contact is directional: the backbone belongs to one residue and the sidechain to the other, so the roles are not interchangeable. A pair (i, j) therefore admits two independent contacts, and either, both, or neither may be present:

0 — neither backbone reaches the other’s sidechain. (The pair may still be a native contact via an H-bond or an SS contact; see build_nonbonded_interaction().)
1 — exactly one direction: backbone of i is within cutoff of a sidechain atom of j, or backbone of j is within cutoff of a sidechain atom of i.
2 — both directions at once: backbone of i touches sidechain of j and backbone of j touches sidechain of i (mutually interdigitated).

That is exactly what the two-step construction below computes: contacts_bs holds directed (backbone-residue -> sidechain-residue) pairs, which are stored in an asymmetric matrix where M[i, j] = 1 means “backbone of i within cutoff of sidechain of j”. The returned value is the symmetrized count M + M.T, whose entries are the number of directions realized. The energy is linear in that count (0, 1, or 2 x BACKBONE_SIDECHAIN_KJ), so a count of 2 marks a tighter, mutually buried pair and earns twice the well depth.

Two easy confusions worth flagging:

This is not the same rule as the H-bond count, which also tops out at 2 (see get_hb_contact_matrix()) but for an unrelated reason: there the 2 counts two H-bonds reported by STRIDE, not two directions. The BS 2 is a structural maximum, not a cap.
Glycine has no sidechain atoms, so it can only ever be the backbone partner — any BS contact involving a glycine caps at 1 by construction.

topo.utils.nonbonded.get_hb_contact_matrix(stride_output_file: str, key_to_index: Dict[Tuple[str, int], int], n_residues: int) → ndarray[source]¶

Build the hydrogen-bond contact matrix directly from a STRIDE output file.

Convenience function that parses the file and builds the matrix. Values are 0, 1, or 2 (multiple H-bonds capped at 2).

Parameters:

stride_output_file (str) – Path to the STRIDE output file.
key_to_index (dict) – Maps (chain, resid) -> 0-based residue index, from get_residue_mapping() (multi-chain safe).
n_residues (int) – Number of residues (must match the system used for STRIDE).

Returns:

Symmetric H-bond contact matrix. See build_hb_contact_matrix().

Return type:

np.ndarray, shape (n_residues, n_residues), dtype int

Example

>>> key_to_index, _, n_residues = get_residue_mapping(u)
>>> hb_matrix = get_hb_contact_matrix("stride.dat", key_to_index, n_residues)
>>> hb_matrix.shape
(283, 283)
>>> (hb_matrix > 1).sum() // 2   # number of pairs with 2 H-bonds (symmetric)
20

topo.utils.nonbonded.get_residue_mapping(universe: Universe) → Tuple[Dict[Tuple[str, int], int], Dict[int, str], int][source]¶

Build residue index and name mappings from an MDAnalysis protein universe.

Uses the order of residues in universe.select_atoms("protein").residues to define 0-based indices. Useful for filling contact/energy matrices by residue index instead of PDB resid.

Parameters:

universe (mda.Universe) – MDAnalysis universe with a loaded protein structure.

Returns:

key_to_index (dict) – Maps (chain, resid) -> 0-based index in residue list, where chain is the normalized chain ID (see _norm_chain()) and resid is the PDB residue ID. Keying on (chain, resid) instead of resid alone is required for multi-chain structures, where the same residue number can appear in more than one chain.
index_to_resname (dict) – Maps 0-based index -> three-letter residue name (e.g. ‘ALA’, ‘GLY’).
n_residues (int) – Total number of protein residues.

Example

>>> u = mda.Universe("protein.pdb")
>>> key_to_idx, idx_to_name, n = get_residue_mapping(u)
>>> key_to_idx[('A', 1)]   # first residue's index (chain A, resid 1)
0
>>> idx_to_name[0]
'MET'

topo.utils.nonbonded.get_scaling_ss_matrix(domain_def: str) → ndarray[source]¶

Build scaling matrix for sidechain–sidechain energies by domain.

Reads domain definitions from YAML and builds an (n × n) matrix where entry [i, j] is the intra-domain nscale if residues i and j are in the same domain, or the inter-domain nscale if they are in different domains (defaulting to 1.0 — no scaling — when no inter nscale is defined for that domain pair).

Parameters:: domain_def (str) – Path to domain definition YAML (see read_yaml_config()).
Returns:: Symmetric matrix. Rows/columns are ordered by sorted residue list (all residues that appear in domain_to_residues). Values are floats (typically 0.0 to 1.0) used to scale SS contact energies.
Return type:: np.ndarray, shape (n, n)

topo.utils.nonbonded.get_ss_contact_matrix(u: Universe, cutoff: float = 4.5) → ndarray[source]¶

Build sidechain–sidechain contact matrix from structure (distance-based).

Two residues are in contact if any sidechain atom of one is within cutoff of any sidechain atom of the other, and they are not within LOCAL_SEPARATION along the sequence. The matrix is binary (0 or 1) and symmetric.

Parameters:

u (mda.Universe) – MDAnalysis universe with protein structure.
cutoff (float, optional) – Distance cutoff in Angstrom (default DEFAULT_CUTOFF = 4.5).

Returns:

Symmetric binary matrix. Entry [i, j] is 1 if residue pair (i, j) has at least one sidechain–sidechain contact within cutoff, else 0. Indices are 0-based residue order from get_residue_mapping().

Return type:

np.ndarray, shape (n_residues, n_residues), dtype int

topo.utils.nonbonded.get_ss_interaction_energy(u: Universe, bt_file: str = 'bt_potential.csv') → ndarray[source]¶

Build residue–residue sidechain interaction energy matrix from BT potential.

Uses the BT potential CSV (residue name vs residue name) and the residue sequence in u to fill an (n_residues × n_residues) matrix of pairwise energies in kJ/mol.

Parameters:

u (mda.Universe) – MDAnalysis universe; residue order determines matrix indices.
bt_file (str, optional) – Passed to load_bt_potential() (default 'bt_potential.csv').

Returns:

Matrix of pairwise sidechain interaction energies in kJ/mol. Entry [i, j] is the BT energy for the residue types at indices i and j.

Return type:

np.ndarray, shape (n_residues, n_residues)

topo.utils.nonbonded.load_bt_potential(bt_file: str = 'bt_potential.csv') → DataFrame[source]¶

Load BT (Betancourt–Thirumalai) potential matrix from CSV.

If bt_file is a filename only (no path), the file is loaded from topo/parameters/data/ so the package works regardless of current working directory. If bt_file is an absolute path, that path is used.

Parameters:: bt_file (str, optional) – Filename (e.g. 'bt_potential.csv') or absolute path to the CSV. Default 'bt_potential.csv' is looked up in topo/parameters/data/.
Returns:: Matrix indexed by residue name (rows and columns). Values are in kJ/mol, computed as KCAL_TO_KJ * |raw_value - bt_shift| from the CSV.
Return type:: pd.DataFrame
Raises:: FileNotFoundError – If the CSV file cannot be found.

Example

>>> df = load_bt_potential()  # uses topo/parameters/data/bt_potential.csv
>>> df.loc['ALA', 'GLY']  # energy for ALA–GLY pair in kJ/mol
1.234

topo.utils.nonbonded.parse_hydrogen_bonds(stride_output_file: str) → List[Tuple][source]¶

Parse hydrogen bonds from a STRIDE output file (donor-acceptor pairs).

Only DNR (donor) lines are read so that each physical H-bond is counted once. STRIDE writes each bond twice (DNR and ACC); using DNR only avoids double-counting while still allowing pairs with multiple H-bonds to appear multiple times in the list.

Parameters:: stride_output_file (str) – Path to the STRIDE output file (e.g. stride.dat).
Returns:: Each element is a pair ((chain1, resid1), (chain2, resid2)) with the tuple sorted so the pair is canonical. chain is the normalized chain ID and resid the PDB residue ID. Repeated pairs indicate multiple H-bonds between the same residue pair. The (chain, resid) keys match those produced by get_residue_mapping(), so multi-chain structures (where a residue number repeats across chains) are resolved correctly.
Return type:: list of tuple
Raises:: FileNotFoundError – If the STRIDE file cannot be opened.

Notes

STRIDE line format: DNR RES chain resid seq -> RES chain resid seq, where chain is a chain letter (e.g. A) or - for a blank chain. The regex captures resname, chain, and the first resid for donor and acceptor.

Example

>>> pairs = parse_hydrogen_bonds("stride.dat")
>>> pairs[0]  # e.g. (('A', 3), ('A', 43))
(('A', 3), ('A', 43))

topo.utils.nonbonded.parse_residue_list(residue_items: List) → List[int][source]¶

Parse a list of residue specifiers into a flat list of residue numbers.

Accepts integers and strings. Strings may be a single number (e.g. "5") or a range "start-end" (inclusive), which is expanded to all integers in [start, end].

Parameters:: residue_items (list) – Elements are int or str. Examples: [1, 2, "5-10", 15].
Returns:: Sorted residue numbers are not guaranteed; order follows the input and range expansion.
Return type:: list of int

Example

>>> parse_residue_list([1, "3-5", 7])
[1, 3, 4, 5, 7]

topo.utils.nonbonded.print_pairs_with_multiple_hb(hb_pairs: List[Tuple]) → Dict[Tuple, int][source]¶

Count H-bonds per residue pair, print pairs with more than one, and return counts.

Parameters:: hb_pairs (list of tuple) – List of canonical H-bond pairs from parse_hydrogen_bonds().
Returns:: Maps canonical pair ((chain1, resid1), (chain2, resid2)) to number of H-bonds (before any cap). Useful for logging or downstream analysis.
Return type:: dict

Example

>>> pairs = parse_hydrogen_bonds("stride.dat")
>>> counts = print_pairs_with_multiple_hb(pairs)
Residue pairs with more than 1 hydrogen bond:
  A3 -- A43:  2 H-bonds
  ...
  (total: 20 pairs)
>>> counts[(('A', 3), ('A', 43))]
2

topo.utils.nonbonded.read_yaml_config(filepath: str) → Tuple[Dict, Dict, Dict, Dict | None][source]¶

Read and parse domain definition YAML (intra/inter nscales, residue lists).

nscale is the native-contact scaling factor applied to the sidechain-sidechain contact well depths – a multiplier, not an absolute energy (the contacts still interact at nscale = 1.0). Per-domain values scale intra-domain contacts; per-interface values scale contacts between two domains.

Required key: n_residues. All three sections are optional (spec §2.2): intra_domains (omit for a single unscaled domain), inter_domains (omit for single-domain proteins; then inter_nscales will be empty), and disordered (omit for a fully-folded protein). Residues not listed in any domain are assigned to domain 'X' with intra nscale 1.0 and inter 1.0 to all other domains.

The optional disordered: section marks intrinsically-disordered (IDR) residues; it carries residues: (required, same syntax as a domain’s residue list), idr_scale: (optional, default DEFAULT_IDR_SCALE = 1.0 – the sequence-specific IDR-IDR attraction knob) and eps_gen_kj: (optional, default DEFAULT_EPS_GEN_KJ = 2.25 kJ/mol – an additive, sequence-independent generic-cohesion depth). The two defaults are the jointly SAXS-calibrated pair. It is returned as the fourth value (None when absent) and consumed by apply_disorder() inside build_nonbonded_interaction().

The per-domain scaling key is nscale. The legacy key strength is still accepted as a deprecated alias (a one-time deprecation notice is printed).

Parameters:

filepath (str) – Path to the YAML file (e.g. domain.yaml).

Returns:

domain_to_residues (dict) – Domain name -> list of residue numbers (1-based).
intra_nscales (dict) – Domain name -> float (intra-domain contact nscale).
inter_nscales (dict) – (domain1, domain2) -> float (inter-domain nscale); symmetric keys (d1, d2) and (d2, d1) are both set.
disorder (dict or None) – None when there is no disordered: section; otherwise {'residues': [int, ...], 'idr_scale': float, 'eps_gen_kj': float} (1-based residue numbers).

Raises:

FileNotFoundError – If the YAML file cannot be opened.
KeyError – If a domain entry has neither an nscale nor a legacy strength value.

Example

YAML format:

n_residues: 110
intra_domains:
  A: { residues: [1-50], nscale: 1.0 }
  B: { residues: [51-110], nscale: 1.0 }
inter_domains:
  A-B: 0.5
disordered:            # optional; None when absent
  residues: [1-24]
  idr_scale: 1.0        # sequence-specific BT scale (default)
  eps_gen_kj: 2.25      # additive generic cohesion, kJ/mol (default)

>>> dom, intra, inter, disorder = read_yaml_config("domain.yaml")
>>> dom["A"][:3]
[1, 2, 3]
>>> inter[("A", "B")]
0.5

topo.utils.runinfo module¶

Run provenance / metadata logging.

Writes a human-readable *_runinfo.log capturing software versions, hardware, GPU/CUDA details, and timing for a simulation run. This is a pure side channel: nothing here changes the simulation, it only records how and where a run was produced, which is useful for reproducibility and for debugging performance differences between machines.

Typical use from a run script:

import time, topo
start = time.time()
...                                   # build model, set up `simulation`
topo.runinfo.write_run_start(
    path, control_file=cfg.config_file, checkpoint_file=checkpoint,
    restart=restart_active, steps_planned=nsteps_remain,
    simulation=simulation, use_gpu=(cfg.device == 'GPU'), ppn=cfg.ppn)
simulation.step(nsteps_remain)
topo.runinfo.write_run_end(path, simulation=simulation, start_epoch=start)

The log is a simple INI-like text file: a one-line # TOPO run info banner followed by aligned [run], [software], [hardware], [gpu] (GPU runs only) and [result] sections. The first four are written before the run so a crashed run still leaves a “what was attempted” file; [result] is appended at the end.

topo.utils.runinfo._context_from_path(path: str) → str[source]¶: Short banner context from the output dir, e.g. .../run/prod/x.log -> run/prod.

topo.utils.runinfo._human_time(seconds: float) → str[source]¶: Wall-clock as a single human-friendly figure (s / min / h).

topo.utils.runinfo._module_version(module) → str[source]¶: Best-effort version string for a Python module.

topo.utils.runinfo._safe_platform_property(simulation, key: str) → str[source]¶: Best-effort OpenMM platform property lookup.

topo.utils.runinfo.collect_cuda_metadata(simulation) → dict[source]¶

Collect CUDA/device metadata when available.

Combines OpenMM platform properties (device name, driver/runtime version, precision, device index) with a best-effort nvidia-smi query for friendlier diagnostics. Every lookup degrades gracefully to a descriptive “not available” string rather than raising.

topo.utils.runinfo.write_banner(path: str, title: str) → None[source]¶: (Over)write the file with the # TOPO run info banner header.

topo.utils.runinfo.write_run_end(path: str, *, simulation, start_epoch: float, final_structure=None, section_label=None) → None[source]¶

Append the [result] section (status + wall-clock timing + final state).

section_label renames the section for a folded multi-phase file (e.g. "result: stage 2 translocation"); defaults to "result".

Parameters:

path (str) – Output run-info log path (same as passed to write_run_start()).
simulation (openmm.app.Simulation) – The simulation, queried for final step count and simulated time.
start_epoch (float) – time.time() value captured at the start of the run, used to compute elapsed wall-clock time.
final_structure (str, optional) – Path to the PDB of the last conformation written for this run (usable as init_position for a follow-up run).

topo.utils.runinfo.write_run_start(path: str, *, control_file, checkpoint_file, restart, steps_planned, simulation, use_gpu, ppn, coord_source=None, vel_source=None, title=None, append=False, section_label=None) → None[source]¶

Write the banner + [run]/[software]/[hardware] (and [gpu]) provenance sections, overwriting any existing file.

append / section_label support a folded, multi-phase run-info file (one file per residue, one section group per stage – used by the CSP runner). When append is True the banner is not re-written and the run-invariant [software]/[hardware]/[gpu] sections are skipped (they were written by the first phase); only the [<section_label>] run section is appended. section_label names that section (default "run"), so successive stages read as e.g. [run: stage 2 translocation].

Records timing, input/checkpoint paths, planned step count, the source of the initial coordinates/velocities, package versions (Python/NumPy/ParmEd/OpenMM), hardware (host, OS, CPU, selected platform), and — for GPU runs — a friendly accelerator line plus a detailed [gpu] section.

Parameters:

path (str) – Output run-info log path (e.g. "traj/traj_runinfo.log").
control_file (str) – Path to the simulation control file used for this run.
checkpoint_file (str) – Resolved checkpoint path for the run.
restart (bool) – Whether the run actually restarted from a checkpoint.
steps_planned (int) – Number of integration steps this invocation will run.
simulation (openmm.app.Simulation) – The constructed simulation (used to read the selected platform / GPU info).
use_gpu (bool) – Whether the run is using the CUDA platform.
ppn (int) – Requested CPU threads (informational).
coord_source (str, optional) – Human-readable description of where the initial coordinates came from (checkpoint, init_position, or pdb_file).
vel_source (str, optional) – Human-readable description of where the initial velocities came from (checkpoint or a Boltzmann distribution).
title (str, optional) – Short banner label (e.g. "L=8 · stage 3"). When omitted, a context is derived from the output directory (e.g. "L_008/stage_3").

topo.utils.runinfo.write_section(path: str, section_name: str, kv_pairs: dict, mode: str = 'a') → None[source]¶: Append one metadata section ([name] header + aligned key : value lines).

topo.utils package¶

Submodules¶

topo.utils.config module¶

topo.utils.multichain module¶

topo.utils.nonbonded module¶

Main workflow¶

Typical use¶

topo.utils.runinfo module¶

Module contents¶