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 with a single [OPTIONS] section. 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 = 1000, dt: ~typing.Any = <factory>, nstxout: int = 10, nstlog: int = 10, nstchk: int | None = None, 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 = 1000
minimize: bool = True
model: str = 'topo'
n_copies: int = 1
nstchk: int | None = None
nstcomm: int | None = None
nstlog: int = 10
nstxout: int = 10
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.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). Must contain an [OPTIONS] section. 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

  1. 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.

  2. 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.

  3. 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.

  4. Non-native contacts: repulsive well with per-residue Rmin/2 from the nearest non-contact CA-CA distance (sum rule); energy ENERGY_PARAMS[‘non_native’].

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.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: ENERGY_PARAMS[‘non_native’].

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.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).

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][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 keys: intra_domains, n_residues. Optional: inter_domains (omit for single-domain proteins; then inter_nscales will be empty). 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 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.

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

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

topo.utils.read_parms module

topo.utils.read_parms.load_dihedral_parameters()[source]

Load the Karanicolas dihedral parameters bundled with the package.

Reads data/karanicolas_dihe_parm.json from the directory containing this module.

Returns:

The parsed dihedral-parameter table from the JSON file.

Return type:

dict

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 with [run_start], [cuda_metadata] (GPU runs only), and [run_end] sections.

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_run_end(path: str, *, simulation, start_epoch: float, final_structure=None) → None[source]

Append the [run_end] provenance section (wall-clock timing + final state).

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) → None[source]

Write the [run_start] provenance section, overwriting any existing file.

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 count, selected platform), and — for GPU runs — a [cuda_metadata] 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).

topo.utils.runinfo.write_section(path: str, section_name: str, kv_pairs: dict, mode: str = 'a') → None[source]

Write one metadata section ([name] header + key: value lines).

Module contents