from __future__ import annotations
import argparse
import copy
import csv
import json
import sys
from collections import deque
from concurrent.futures import ProcessPoolExecutor
from multiprocessing import get_context
from pathlib import Path
from typing import Dict, List, Optional, Tuple
import numpy as np
from tqdm import tqdm
from pycurves_lib.core.curves_dataclasses import MolecularStructure
from pycurves_lib.io.curves_mol_loader import MolecularLoader
from pycurves_lib.io.curves_output import CurvesOutputFormatter, _json_dumps, _to_jsonable
from pycurves_lib.curves_wrapper import CurvesWrapper
from pycurves_lib.cli.pycurves_cli_options import (
add_pycurves_analysis_options,
pycurves_runner_kwargs,
)
from pycurves_lib.md.trajectory_loader import TrajectoryLoader
from pycurves_lib.topology.frame_annotations import infer_frame_pair_observations
from pycurves_lib.md.trajectory_statistics import (
backbone_conformer_population_tables,
circular_degree_summary,
is_circular_degree_column,
linear_summary,
)
class MDTrajectoryAnalyzer:
"""Run pyCurves over selected trajectory frames and summarize named records."""
def __init__(
self,
topology_file: str,
trajectory_file: Optional[str] = None,
inpfile: Optional[str] = None,
output_dir: str = ".",
frame_convention: str = "standard",
axis_convention: str = "global",
axis_weighting: Optional[bool] = None,
continuous_strands: bool = False,
altloc: Optional[str] = None,
fit_override: Optional[bool] = None,
grv_override: Optional[bool] = None,
mini_override: Optional[bool] = None,
comb_override: Optional[bool] = None,
ends_override: Optional[bool] = None,
topology_mode: str = "reference",
duplex_only: bool = False,
):
self.topology_file = topology_file
self.trajectory_file = trajectory_file
self.output_dir = output_dir
self.frame_convention, self.axis_convention = CurvesWrapper.normalize_conventions(
frame_convention,
axis_convention,
)
self.continuous_strands = continuous_strands
self.duplex_only = duplex_only
self.altloc = MolecularLoader.normalize_altloc(altloc)
self.fit_override = fit_override
self.grv_override = grv_override
self.mini_override = mini_override
self.comb_override = comb_override
self.ends_override = ends_override
self.axis_weighting = axis_weighting
topology_mode = str(topology_mode).strip().lower()
if topology_mode not in {"reference", "annotate"}:
raise ValueError("topology_mode must be one of: reference, annotate")
self.topology_mode = topology_mode
self.reference_topology_file = self._reference_topology(topology_file, trajectory_file, output_dir)
self.template_molecule = self._load_template_molecule(self.reference_topology_file, self.altloc)
self.runner_kwargs = {
"continuous_strands": continuous_strands,
"duplex_only": duplex_only,
"altloc": self.altloc,
"frame_convention": self.frame_convention,
"axis_convention": self.axis_convention,
"axis_weighting": axis_weighting,
"fit_override": fit_override,
"grv_override": grv_override,
"mini_override": mini_override,
"comb_override": comb_override,
"ends_override": ends_override,
}
self.reference_runner = CurvesWrapper(
pdbfile=self.reference_topology_file,
inpfile=inpfile,
output_dir=output_dir,
**self.runner_kwargs,
)
self.inpfile = self.reference_runner.inpfile
self.generated_inpfiles = self.reference_runner.generated_inpfiles
def run(
self,
frame_selector,
selection: Dict,
mode: str = "summary",
mini: Optional[bool] = None,
verbose: bool = False,
warm_start: bool = True,
axis_continuity: bool = True,
workers: int = 1,
) -> Dict:
if workers <= 0:
raise ValueError("workers must be positive")
if self.axis_convention == "local":
mini = False
warm_start = False
if workers > 1 and warm_start:
raise ValueError("workers > 1 requires warm_start=False")
if workers > 1 and axis_continuity:
raise ValueError("workers > 1 requires axis_continuity=False")
frame_payloads = []
table_records: Dict[str, List[Dict]] = {}
processed = 0
effective_mini = False
effective_axis_weighting = False
prev_helical = None
axis_sign_reference = None
iterator = TrajectoryLoader.iter_frames(
self.topology_file,
self.trajectory_file,
frame_selector,
)
def collect_frame(frame_index, frame_time, dataframes) -> None:
nonlocal processed
if mode in {"summary", "both"}:
for table_name, rows in dataframes.items():
if isinstance(rows, list):
for row in rows:
row = dict(row)
row["frame"] = frame_index
row["time"] = frame_time
table_records.setdefault(table_name, []).append(row)
if mode in {"per-frame", "both"}:
frame_payloads.append({
"frame": frame_index,
"time": frame_time,
"dataframes": dataframes,
})
processed += 1
if workers == 1:
for frame in tqdm(iterator, desc="Processing frames"):
(
dataframes,
effective_mini,
current_helical,
axis_direction_signs,
effective_axis_weighting,
) = self._analyze_frame(
frame.coordinates,
mini,
verbose,
prev_helical=prev_helical if warm_start else None,
axis_sign_reference=axis_sign_reference if axis_continuity else None,
capture_helical=warm_start,
)
if current_helical is not None:
prev_helical = current_helical
if axis_continuity and axis_sign_reference is None and axis_direction_signs:
axis_sign_reference = np.asarray(axis_direction_signs, dtype=int)
collect_frame(frame.index, frame.time, dataframes)
else:
pending = deque()
with ProcessPoolExecutor(
max_workers=workers,
mp_context=get_context("spawn"),
initializer=_initialize_md_worker,
initargs=(self,),
) as executor, tqdm(desc="Processing frames") as progress:
for frame in iterator:
pending.append(executor.submit(
_process_md_frame_worker,
frame.coordinates,
frame.index,
frame.time,
mini,
verbose,
))
if len(pending) >= workers:
(
frame_index,
frame_time,
dataframes,
effective_mini,
effective_axis_weighting,
) = pending.popleft().result()
collect_frame(frame_index, frame_time, dataframes)
progress.update()
while pending:
(
frame_index,
frame_time,
dataframes,
effective_mini,
effective_axis_weighting,
) = pending.popleft().result()
collect_frame(frame_index, frame_time, dataframes)
progress.update()
if processed == 0:
raise ValueError("No trajectory frames matched the requested frame selection.")
payload = {
"program": "pyCurves",
"format": "pycurves-trajectory-slim-v2",
"inputs": {
"topology_file": self.topology_file,
"reference_topology_file": self.reference_topology_file,
"trajectory_file": self.trajectory_file,
"inpfile": self.inpfile,
"generated_inpfiles": self.generated_inpfiles,
},
"analysis_options": {
"continuous_strands": self.continuous_strands,
"duplex_only": getattr(self, "duplex_only", False),
"altloc": self.altloc or "first",
"fit": self.fit_override,
"grooves": self.grv_override,
"mini": effective_mini,
"comb": self.comb_override,
"ends": self.ends_override,
"axis_convention": self.axis_convention,
"axis_weighting": effective_axis_weighting,
"topology_mode": self.topology_mode,
"workers": workers,
},
"selection": {
**selection,
"processed_frames": processed,
"warm_start": bool(warm_start and effective_mini),
"axis_continuity": bool(axis_continuity),
},
"frame_convention": {
"name": self.frame_convention,
"compatible_with": ["Curves+", "3DNA", "x3dna"] if self.frame_convention == "standard" else ["Curves 5.3"],
"axis_convention": self.axis_convention,
},
}
if mode in {"per-frame", "both"}:
payload["frames"] = frame_payloads
if mode in {"summary", "both"}:
payload["summary"] = self._summarize_tables(
table_records,
total_frames=processed,
)
return payload
def _analyze_frame(
self,
coordinates: np.ndarray,
mini: Optional[bool],
verbose: bool,
prev_helical: Optional[np.ndarray] = None,
axis_sign_reference: Optional[np.ndarray] = None,
capture_helical: bool = False,
) -> Tuple[Dict, bool, Optional[np.ndarray], List[int], bool]:
molecule = self._molecule_for_frame(coordinates)
runner = self.reference_runner
runner.analyze_molecule(
molecule,
mini=mini,
verbose=verbose,
prev_opt_helical=prev_helical,
axis_sign_reference=axis_sign_reference,
)
if self.topology_mode == "annotate":
reference_rows = runner.ctx.annotations.get(
"base_pair_annotations", []
)
runner.ctx.annotations["frame_base_pair_observations"] = (
infer_frame_pair_observations(molecule, reference_rows)
)
effective_mini = bool(runner.ctx.cfg.mini)
current_helical = (
runner.ctx.params.helical.copy()
if capture_helical and effective_mini
else None
)
axis_direction_signs: List[int] = []
if hasattr(runner.calc, "axis_direction_sign"):
axis_direction_signs = [
int(value)
for value in np.asarray(
runner.calc.axis_direction_sign[:runner.ctx.nst]
).tolist()
]
formatter = CurvesOutputFormatter(runner, fast_groove=True)
dataframes = self._normalize_frame_dataframes(formatter._build_dataframes())
return (
dataframes,
effective_mini,
current_helical,
axis_direction_signs,
bool(runner.ctx.cfg.axis_weighting),
)
[docs]
def write_csv(self, payload: Dict, prefix: str) -> None:
prefix_path = Path(prefix)
prefix_path.parent.mkdir(parents=True, exist_ok=True)
if "summary" in payload:
for name, rows in payload["summary"].items():
self._write_rows_csv(f"{prefix}_{name}_summary.csv", rows)
if "frames" in payload:
per_table: Dict[str, List[Dict]] = {}
for frame in payload["frames"]:
for name, rows in frame["dataframes"].items():
if isinstance(rows, list):
for row in rows:
row = dict(row)
row["frame"] = frame["frame"]
row["time"] = frame["time"]
per_table.setdefault(name, []).append(row)
for name, rows in per_table.items():
self._write_rows_csv(f"{prefix}_{name}_frames.csv", rows)
@staticmethod
def _load_template_molecule(topology_file: str, altloc=None) -> MolecularStructure:
holder = type("MoleculeHolder", (), {"molecule": MolecularStructure()})()
MolecularLoader.load(topology_file, holder, altloc=altloc)
return holder.molecule
@staticmethod
def _reference_topology(topology_file: str, trajectory_file: Optional[str], output_dir: str) -> str:
if trajectory_file is not None:
return topology_file
path = Path(topology_file)
if path.suffix.lower() not in {".pdb", ".brk"}:
return topology_file
output_root = Path(output_dir)
output_root.mkdir(parents=True, exist_ok=True)
ref_path = output_root / f"{path.stem}_first_model.pdb"
if ref_path.exists():
return str(ref_path)
lines = []
in_model = False
saw_model = False
with path.open("r", encoding="utf-8", errors="ignore") as handle:
for line in handle:
record = line[:6]
if record.startswith("MODEL"):
saw_model = True
in_model = True
continue
if record.startswith("ENDMDL") and in_model:
break
if saw_model and not in_model:
continue
if record in {"ATOM ", "HETATM", "TITLE ", "CRYST1", "TER "}:
lines.append(line)
if saw_model and lines:
ref_path.write_text("".join(lines) + "END\n", encoding="utf-8")
return str(ref_path)
return topology_file
def _molecule_for_frame(self, coordinates: np.ndarray) -> MolecularStructure:
coordinates = np.asarray(coordinates, dtype=float)
if coordinates.shape != self.template_molecule.coordinates.shape:
source_indices = self.template_molecule.source_atom_indices
if (source_indices is not None and coordinates.ndim == 2 and coordinates.shape[1] == 3
and len(source_indices) == self.template_molecule.coordinates.shape[0]
and len(source_indices) > 0 and int(np.max(source_indices)) < coordinates.shape[0]):
coordinates = coordinates[source_indices]
if coordinates.shape != self.template_molecule.coordinates.shape:
raise ValueError(
"Trajectory frame atom count does not match the topology molecule: "
f"{coordinates.shape[0]} vs {self.template_molecule.coordinates.shape[0]}"
)
molecule = copy.copy(self.template_molecule)
molecule.coordinates = coordinates.copy()
return molecule
@staticmethod
def _normalize_frame_dataframes(dataframes: Dict) -> Dict:
normalized = {}
for name, rows in dataframes.items():
if name == "groove" and isinstance(rows, dict) and "data" in rows:
normalized[name] = MDTrajectoryAnalyzer._flatten_groove_records(rows)
else:
normalized[name] = rows
return normalized
@staticmethod
def _flatten_groove_records(groove: Dict) -> List[Dict]:
flat_rows: List[Dict] = []
for level, level_data in groove.get("data", {}).items():
for sub_level, sub_data in level_data.get("sub_levels", {}).items():
row = {
"atom_defining_backbone": groove.get("atom_defining_backbone", ""),
"total_levels": groove.get("levels", 0),
"total_sub_levels": groove.get("sub_levels", 0),
"level": int(level),
"base_pair": level_data.get("base_pair", ""),
"sub_level": int(sub_level),
}
for key, value in sub_data.items():
if key != "geometry":
row[key] = value
flat_rows.append(row)
return flat_rows
@staticmethod
def _summarize_tables(
tables: Dict[str, List[Dict]],
total_frames: Optional[int] = None,
) -> Dict[str, List[Dict]]:
summaries = {}
for name, rows in tables.items():
if name == "base_pair_observations":
summaries[name] = MDTrajectoryAnalyzer._pairing_profile_rows(
rows,
total_frames=total_frames,
)
continue
if name == "annotations":
# These are categorical event/detail tables. Treating residue,
# strand, status, and diagnostic fields as generic grouping keys
# produces misleading columns such as strand_1_mean. Pair-state
# populations are reported by base_pair_observations instead.
continue
if not rows:
summaries[name] = []
continue
all_cols = list(dict.fromkeys(col for row in rows for col in row))
identifier_cols = {
"strand", "level", "residue_id", "subunit", "sub_level",
"total_levels", "total_sub_levels", "partner_strand",
"sequence_index", "next_level", "fit_atom_count",
}
numeric_cols = []
for col in all_cols:
if col in {"frame", "time"} or col in identifier_cols:
continue
values = [row.get(col) for row in rows if row.get(col) is not None]
if values and all(isinstance(v, (int, float, np.integer, np.floating)) for v in values):
numeric_cols.append(col)
key_cols = [col for col in all_cols if col not in numeric_cols and col not in {"frame", "time"}]
if not numeric_cols:
summaries[name] = []
continue
groups = {}
if key_cols:
for row in rows:
key = tuple(MDTrajectoryAnalyzer._summary_key_value(row.get(col)) for col in key_cols)
groups.setdefault(key, []).append(row)
else:
groups[()] = rows
table_summary = []
for key, group_rows in groups.items():
out = {col: key[idx] for idx, col in enumerate(key_cols)}
out["count"] = len(group_rows)
for col in numeric_cols:
vals = np.asarray(
[row.get(col) for row in group_rows if row.get(col) is not None],
dtype=float,
)
if MDTrajectoryAnalyzer._is_circular_degree_column(col):
stats = circular_degree_summary(vals)
else:
stats = linear_summary(vals)
out[f"{col}_count"] = int(np.count_nonzero(np.isfinite(vals)))
out[f"{col}_mean"] = stats.mean
out[f"{col}_stddev"] = stats.stddev
table_summary.append(out)
summaries[name] = table_summary
if "backbone" in tables:
summaries.update(backbone_conformer_population_tables(tables["backbone"]))
return summaries
@staticmethod
def _pairing_profile_rows(
rows: List[Dict],
total_frames: Optional[int] = None,
) -> List[Dict]:
"""Summarize categorical pair states without averaging identifiers."""
if not rows:
return []
if total_frames is None:
frame_ids = {
row.get("frame")
for row in rows
if row.get("frame") is not None
}
total_frames = len(frame_ids) if frame_ids else 0
total_frames = max(int(total_frames or 0), 0)
by_pair: Dict[str, List[Dict]] = {}
for row in rows:
pair_id = str(row.get("pair_id") or "").strip()
if not pair_id:
pair_id = (
f"{row.get('residue_1', '')}:"
f"{row.get('residue_2', '')}"
)
by_pair.setdefault(pair_id, []).append(row)
profile: List[Dict] = []
status_order = {"present": 0, "absent": 1, "uncertain": 2}
for pair_id in sorted(by_pair):
pair_rows = by_pair[pair_id]
first = pair_rows[0]
reference_pair = bool(first.get("reference_pair", True))
state_counts: Dict[tuple, int] = {}
for row in pair_rows:
state = (
str(row.get("pair_status") or "uncertain"),
str(row.get("pairing_mode") or ""),
str(row.get("candidate_mode") or ""),
str(row.get("classification_status") or "unassigned"),
str(row.get("observed_lw_family") or ""),
str(row.get("helical_context") or ""),
)
state_counts[state] = state_counts.get(state, 0) + 1
missing_count = max(total_frames - len(pair_rows), 0)
if missing_count:
# Newly discovered pairs have no row in frames where the
# coordinate-only inferrer did not detect them.
missing_status = "uncertain" if reference_pair else "absent"
missing_state = (
missing_status,
"",
"",
"unassigned",
"",
"",
)
state_counts[missing_state] = (
state_counts.get(missing_state, 0) + missing_count
)
for state, frame_count in sorted(
state_counts.items(),
key=lambda item: (
status_order.get(item[0][0], 99),
item[0][1],
item[0][2],
item[0][4],
item[0][5],
),
):
(
pair_status,
pairing_mode,
candidate_mode,
classification_status,
observed_lw_family,
helical_context,
) = state
profile.append({
"pair_id": pair_id,
"reference_pair": reference_pair,
"level": first.get("level"),
"residue_1": first.get("residue_1"),
"residue_2": first.get("residue_2"),
"pair_status": pair_status,
"pairing_mode": pairing_mode,
"candidate_mode": candidate_mode,
"classification_status": classification_status,
"observed_lw_family": observed_lw_family,
"helical_context": helical_context,
"frame_count": frame_count,
"frame_fraction": (
frame_count / total_frames if total_frames else None
),
})
return profile
@staticmethod
def _is_circular_degree_column(column_name: str) -> bool:
return is_circular_degree_column(column_name)
@staticmethod
def _summary_key_value(value):
"""Return a hashable grouping key for JSON-like table identifier values."""
if isinstance(value, (list, dict)):
return json.dumps(_to_jsonable(value), sort_keys=True)
if isinstance(value, np.ndarray):
return json.dumps(_to_jsonable(value), sort_keys=True)
if isinstance(value, np.generic):
return value.item()
return value
@staticmethod
def _write_rows_csv(path: str, rows: List[Dict]) -> None:
fieldnames = list(dict.fromkeys(key for row in rows for key in row.keys()))
with Path(path).open("w", encoding="utf-8", newline="") as handle:
if not fieldnames:
handle.write("")
return
writer = csv.DictWriter(handle, fieldnames=fieldnames, extrasaction="ignore")
writer.writeheader()
writer.writerows(rows)
_WORKER_MD_ANALYZER: Optional[MDTrajectoryAnalyzer] = None
def _initialize_md_worker(analyzer: MDTrajectoryAnalyzer) -> None:
global _WORKER_MD_ANALYZER
_WORKER_MD_ANALYZER = analyzer
def _process_md_frame_worker(
coordinates: np.ndarray,
frame_index: int,
frame_time: Optional[float],
mini: Optional[bool],
verbose: bool,
):
if _WORKER_MD_ANALYZER is None: # pragma: no cover - process-pool guard
raise RuntimeError("MD worker was not initialized.")
dataframes, effective_mini, _, _, axis_weighting = (
_WORKER_MD_ANALYZER._analyze_frame(
coordinates,
mini,
verbose,
)
)
return frame_index, frame_time, dataframes, effective_mini, axis_weighting
def make_frame_selector(spec: Optional[str], start: Optional[int], stop: Optional[int], step: int):
if spec:
exact = set()
ranges = []
explicit_indices = set()
all_ranges_finite = True
range_specs = []
for part in spec.split(","):
part = part.strip()
if not part:
continue
if ":" in part:
values = [int(v) if v else None for v in part.split(":")]
while len(values) < 3:
values.append(None)
s, e, st = values
st = st or 1
range_start = 0 if s is None else s
ranges.append((range_start, e, st))
range_specs.append((range_start, e, st))
if e is None:
all_ranges_finite = False
else:
explicit_indices.update(range(range_start, e, st))
else:
value = int(part)
exact.add(value)
explicit_indices.add(value)
def selected(index: int) -> bool:
if index in exact:
return True
for s, e, st in ranges:
if index >= s and (e is None or index < e) and (index - s) % st == 0:
return True
return False
if all_ranges_finite:
selected.explicit_indices = sorted(explicit_indices)
if not exact and len(range_specs) == 1:
selected.mdtraj_range = range_specs[0]
if hasattr(selected, "explicit_indices"):
delattr(selected, "explicit_indices")
return selected, {"frames": spec}
s = 0 if start is None else start
e = stop
def selected(index: int) -> bool:
return index >= s and (e is None or index < e) and (index - s) % step == 0
if e is not None:
selected.mdtraj_range = (s, e, step)
return selected, {"start": s, "stop": e, "step": step}
def collect_available_frame_indices(topology_file: str, trajectory_file: Optional[str]) -> List[int]:
return [frame.index for frame in TrajectoryLoader.iter_frames(topology_file, trajectory_file)]
def analyze_trajectory(
topology_file: str,
trajectory_file: Optional[str] = None,
inpfile: Optional[str] = None,
output_dir: str = ".",
frames: Optional[str] = None,
start: Optional[int] = None,
stop: Optional[int] = None,
step: int = 1,
mode: str = "per-frame",
frame_convention: str = "standard",
axis_convention: str = "global",
axis_weighting: Optional[bool] = None,
continuous_strands: bool = False,
altloc: Optional[str] = None,
fit: Optional[bool] = None,
grooves: Optional[bool] = None,
mini: Optional[bool] = None,
comb: Optional[bool] = None,
ends: Optional[bool] = None,
topology_mode: str = "reference",
verbose: bool = False,
warm_start: bool = True,
axis_continuity: bool = True,
duplex_only: bool = False,
workers: int = 1,
) -> Dict:
"""Run pyCurves trajectory analysis from Python and return the JSON-like payload.
This is the notebook-friendly equivalent of ``pycurves-md``. Use
``mode="per-frame"`` or ``mode="both"`` when the payload will be converted
to long-form DataFrames with ``pycurves_md_plot.extract_block``.
"""
if mode not in {"per-frame", "summary", "both"}:
raise ValueError("mode must be one of: per-frame, summary, both")
if step <= 0:
raise ValueError("step must be positive")
frame_selector, selection = make_frame_selector(frames, start, stop, step)
analyzer = MDTrajectoryAnalyzer(
topology_file=topology_file,
trajectory_file=trajectory_file,
inpfile=inpfile,
output_dir=output_dir,
frame_convention=frame_convention,
axis_convention=axis_convention,
axis_weighting=axis_weighting,
continuous_strands=continuous_strands,
duplex_only=duplex_only,
altloc=altloc,
fit_override=fit,
grv_override=grooves,
mini_override=mini,
comb_override=comb,
ends_override=ends,
topology_mode=topology_mode,
)
return analyzer.run(
frame_selector=frame_selector,
selection=selection,
mode=mode,
mini=mini,
verbose=verbose,
warm_start=warm_start,
axis_continuity=axis_continuity,
workers=workers,
)
def main() -> None:
from pycurves_lib import __version__
parser = argparse.ArgumentParser(description="Run pyCurves over an MD trajectory.")
parser.add_argument("--version", action="version", version=f"%(prog)s {__version__}")
parser.add_argument("topology", help="Topology/reference structure file, usually PDB or CIF.")
parser.add_argument("trajectory", nargs="?", help="Trajectory file. Omit for multi-model PDB input.")
parser.add_argument("--inp", help="Existing Curves .inp file. If omitted, inferred from the topology.")
parser.add_argument("--output-dir", default=".", help="Directory for auto-generated .inp files.")
parser.add_argument("--output-file", default="pycurves_trajectory.json", help="JSON file or CSV prefix.")
parser.add_argument("--mode", choices=["per-frame", "summary", "both"], default="summary")
parser.add_argument("--format", choices=["json", "csv"], default="json")
parser.add_argument("--frames", help="Frame list/ranges, e.g. '0,10,20:100:5'. Overrides start/stop/step.")
parser.add_argument("--start", type=int, help="First frame index to include.")
parser.add_argument("--stop", type=int, help="Stop before this frame index.")
parser.add_argument("--step", type=int, default=1, help="Frame stride.")
add_pycurves_analysis_options(parser)
parser.add_argument(
"--topology-mode",
choices=["reference", "annotate"],
default="reference",
help=(
"Use the reference pair map only, or infer coordinate-supported "
"pair presence and pairing mode independently for every frame."
),
)
parser.add_argument(
"--no-warm-start",
action="store_true",
help="Do not seed each frame's optimizer from the previous frame. Slower, but useful for diagnosing axis branch flips.",
)
parser.add_argument(
"--no-axis-continuity",
action="store_true",
help="Do not keep the Curves global-axis direction signs aligned to the first processed frame.",
)
parser.add_argument(
"--workers",
type=int,
default=1,
help=(
"Independent frame worker processes (default: 1). Values above 1 "
"require --no-warm-start and --no-axis-continuity."
),
)
parser.add_argument("--verbose", action="store_true", help="Print per-frame pyCurves logs.")
args = parser.parse_args()
frame_selector, selection = make_frame_selector(args.frames, args.start, args.stop, args.step)
analyzer = MDTrajectoryAnalyzer(
topology_file=args.topology,
trajectory_file=args.trajectory,
inpfile=args.inp,
output_dir=args.output_dir,
topology_mode=args.topology_mode,
**pycurves_runner_kwargs(args),
)
try:
payload = analyzer.run(
frame_selector=frame_selector,
selection=selection,
mode=args.mode,
mini=args.mini,
verbose=args.verbose,
warm_start=not args.no_warm_start,
axis_continuity=not args.no_axis_continuity,
workers=args.workers,
)
except (ValueError, ImportError, NotImplementedError) as exc:
raise SystemExit(str(exc)) from exc
if args.format == "json":
Path(args.output_file).write_text(_json_dumps(payload), encoding="utf-8")
else:
analyzer.write_csv(payload, args.output_file.removesuffix(".csv"))
generated = analyzer.generated_inpfiles
if generated:
print("Generated input file(s):", file=sys.stderr)
for path in generated:
print(f" {path}", file=sys.stderr)
print(f"Processed {payload['selection']['processed_frames']} frame(s).", file=sys.stderr)
if __name__ == "__main__":
main()