From 2cc43801d0ef344edbdf2255360c76ae4a978f23 Mon Sep 17 00:00:00 2001 From: Corin Wagen Date: Thu, 2 Jul 2026 16:00:50 -0400 Subject: [PATCH] add constraints --- openconf/__init__.py | 13 +- openconf/config.py | 133 ++++++++++++++++++++- openconf/constraints.py | 235 ++++++++++++++++++++++++++++++++++++- openconf/propose/hybrid.py | 45 ++++--- tests/test_constrained.py | 91 ++++++++++++++ 5 files changed, 498 insertions(+), 19 deletions(-) diff --git a/openconf/__init__.py b/openconf/__init__.py index c1c5c63..6c4ad17 100644 --- a/openconf/__init__.py +++ b/openconf/__init__.py @@ -6,7 +6,15 @@ generate_conformers, generate_conformers_from_pose, ) -from .config import ConformerConfig, ConformerPreset, ConstraintSpec, preset_config +from .config import ( + AngleConstraintSpec, + BondConstraintSpec, + ConformerConfig, + ConformerPreset, + ConstraintSpec, + TorsionConstraintSpec, + preset_config, +) from .dedupe import prism_dedupe from .exceptions import OpenConfError, OpenConfRuntimeError, OpenConfValueError from .io import mol_to_smiles, read_sdf, read_xyz, smiles_to_mol, write_sdf, write_xyz @@ -15,6 +23,8 @@ from .torsionlib import TorsionLibrary, TorsionRule __all__ = [ + "AngleConstraintSpec", + "BondConstraintSpec", "ConformerConfig", "ConformerEnsemble", "ConformerPreset", @@ -26,6 +36,7 @@ "RDKitMMFFMinimizer", "Rotor", "RotorModel", + "TorsionConstraintSpec", "TorsionLibrary", "TorsionRule", "build_rotor_model", diff --git a/openconf/config.py b/openconf/config.py index d3b01b3..2548517 100644 --- a/openconf/config.py +++ b/openconf/config.py @@ -31,6 +31,10 @@ ) _SUPPORTED_PARENT_STRATEGIES = frozenset({"softmax", "uniform", "best"}) _SUPPORTED_FINAL_SELECTIONS = frozenset({"energy", "diverse"}) +_DEFAULT_INTERNAL_FORCE_CONSTANT = 100000.0 +_DEFAULT_BOND_TOLERANCE = 0.01 +_DEFAULT_ANGLE_TOLERANCE_DEG = 0.5 +_DEFAULT_TORSION_TOLERANCE_DEG = 0.5 def _require_at_least(name: str, value: float, minimum: float) -> None: @@ -69,13 +73,113 @@ def _validate_move_probs(move_probs: dict[str, float]) -> None: raise OpenConfValueError("move_probs must sum to a positive value.") +def _require_atom_indices(name: str, atoms: tuple[int, ...], count: int) -> None: + if len(atoms) != count: + raise OpenConfValueError(f"{name} must contain exactly {count} atom indices.") + if any(atom < 0 for atom in atoms): + raise OpenConfValueError(f"{name} atom indices must be >= 0.") + if len(set(atoms)) != count: + raise OpenConfValueError(f"{name} atom indices must be distinct.") + + +@dataclass(frozen=True) +class BondConstraintSpec: + """Specification for bond-distance constraint. + + Omit `distance` to freeze current distance from first input conformer. + + Attributes: + atom_i: first atom index + atom_j: second atom index + distance: target distance in Angstrom, or None to use reference conformer + tolerance: allowed half-width around target distance in Angstrom + force_constant: force constant in kcal/mol/A^2 + """ + + atom_i: int + atom_j: int + distance: float | None = None + tolerance: float = _DEFAULT_BOND_TOLERANCE + force_constant: float = _DEFAULT_INTERNAL_FORCE_CONSTANT + + def __post_init__(self) -> None: + _require_atom_indices("BondConstraintSpec", (self.atom_i, self.atom_j), 2) + _require_optional_at_least("distance", self.distance, 0.0) + _require_at_least("tolerance", self.tolerance, 0.0) + _require_greater_than("force_constant", self.force_constant, 0.0) + + +@dataclass(frozen=True) +class AngleConstraintSpec: + """Specification for valence-angle constraint. + + Omit `angle_deg` to freeze current angle from first input conformer. + + Attributes: + atom_i: first atom index + atom_j: central atom index + atom_k: third atom index + angle_deg: target angle in degrees, or None to use reference conformer + tolerance_deg: allowed half-width around target angle in degrees + force_constant: force constant in kcal/mol/rad^2 + """ + + atom_i: int + atom_j: int + atom_k: int + angle_deg: float | None = None + tolerance_deg: float = _DEFAULT_ANGLE_TOLERANCE_DEG + force_constant: float = _DEFAULT_INTERNAL_FORCE_CONSTANT + + def __post_init__(self) -> None: + _require_atom_indices("AngleConstraintSpec", (self.atom_i, self.atom_j, self.atom_k), 3) + _require_optional_at_least("angle_deg", self.angle_deg, 0.0) + if self.angle_deg is not None and self.angle_deg > 180.0: + raise OpenConfValueError(f"angle_deg must be <= 180.0, got {self.angle_deg}.") + _require_at_least("tolerance_deg", self.tolerance_deg, 0.0) + _require_greater_than("force_constant", self.force_constant, 0.0) + + +@dataclass(frozen=True) +class TorsionConstraintSpec: + """Specification for dihedral-angle constraint. + + Omit `dihedral_deg` to freeze current torsion from first input conformer. + + Attributes: + atom_i: first atom index + atom_j: second atom index + atom_k: third atom index + atom_l: fourth atom index + dihedral_deg: target dihedral in degrees, or None to use reference conformer + tolerance_deg: allowed half-width around target dihedral in degrees + force_constant: force constant in kcal/mol/rad^2 + """ + + atom_i: int + atom_j: int + atom_k: int + atom_l: int + dihedral_deg: float | None = None + tolerance_deg: float = _DEFAULT_TORSION_TOLERANCE_DEG + force_constant: float = _DEFAULT_INTERNAL_FORCE_CONSTANT + + def __post_init__(self) -> None: + _require_atom_indices("TorsionConstraintSpec", (self.atom_i, self.atom_j, self.atom_k, self.atom_l), 4) + if self.dihedral_deg is not None and not -360.0 <= self.dihedral_deg <= 360.0: + raise OpenConfValueError(f"dihedral_deg must be between -360.0 and 360.0, got {self.dihedral_deg}.") + _require_at_least("tolerance_deg", self.tolerance_deg, 0.0) + _require_greater_than("force_constant", self.force_constant, 0.0) + + @dataclass class ConstraintSpec: - """Specification for positional constraints during conformer generation. + """Specification for geometry constraints during conformer generation. Used for FEP-style analogue generation where an MCS-aligned pose is provided and a subset of atoms (the core scaffold) must remain fixed while terminal - groups are explored. + groups are explored. Can also freeze or target internal coordinates using + bond, angle, and torsion constraint specifications. Attributes: constrained_atoms: Atom indices that must not move. These are indices @@ -85,14 +189,37 @@ class ConstraintSpec: position_force_constant: MMFF force constant (kcal/mol/Ų) for the harmonic position restraints applied to constrained atoms. Default 1000.0 is very stiff and effectively freezes the core. + bond_constraints: bond-distance constraints + angle_constraints: valence-angle constraints + torsion_constraints: dihedral-angle constraints """ - constrained_atoms: frozenset[int] + constrained_atoms: frozenset[int] = frozenset() position_force_constant: float = 1000.0 + bond_constraints: tuple[BondConstraintSpec, ...] = () + angle_constraints: tuple[AngleConstraintSpec, ...] = () + torsion_constraints: tuple[TorsionConstraintSpec, ...] = () def __post_init__(self) -> None: + if any(atom < 0 for atom in self.constrained_atoms): + raise OpenConfValueError("constrained_atoms values must be >= 0.") _require_greater_than("position_force_constant", self.position_force_constant, 0.0) + @property + def requires_reference_geometry(self) -> bool: + """Return whether constraint targets must be read from input conformer.""" + return bool( + self.constrained_atoms + or any(constraint.distance is None for constraint in self.bond_constraints) + or any(constraint.angle_deg is None for constraint in self.angle_constraints) + or any(constraint.dihedral_deg is None for constraint in self.torsion_constraints) + ) + + @property + def constrained_rotor_atoms(self) -> frozenset[int]: + """Atom indices whose movement should suppress torsion moves.""" + return self.constrained_atoms + @dataclass class ConformerConfig: diff --git a/openconf/constraints.py b/openconf/constraints.py index 4a44fc8..4df2faa 100644 --- a/openconf/constraints.py +++ b/openconf/constraints.py @@ -1,9 +1,13 @@ """Reusable geometry constraints for conformer relaxation.""" from dataclasses import dataclass +from typing import Protocol import numpy as np from rdkit import Chem +from rdkit.Chem import rdMolTransforms + +from .exceptions import OpenConfValueError _DEFAULT_POSITION_FORCE_CONSTANT = 1000.0 _METAL_POSITION_FORCE_CONSTANT = 1e4 @@ -11,6 +15,33 @@ _METAL_LIGAND_DISTANCE_FORCE_CONSTANT = 100000.0 +class _BondConstraintSpecLike(Protocol): + atom_i: int + atom_j: int + distance: float | None + tolerance: float + force_constant: float + + +class _AngleConstraintSpecLike(Protocol): + atom_i: int + atom_j: int + atom_k: int + angle_deg: float | None + tolerance_deg: float + force_constant: float + + +class _TorsionConstraintSpecLike(Protocol): + atom_i: int + atom_j: int + atom_k: int + atom_l: int + dihedral_deg: float | None + tolerance_deg: float + force_constant: float + + @dataclass(frozen=True) class PositionConstraint: """Harmonic position constraint on an atom. @@ -45,6 +76,50 @@ class DistanceConstraint: force_constant: float +@dataclass(frozen=True) +class AngleConstraint: + """Harmonic angle-window constraint between three atoms. + + Attributes: + atom_i: first atom index + atom_j: central atom index + atom_k: third atom index + min_angle_deg: lower angle bound in degrees + max_angle_deg: upper angle bound in degrees + force_constant: force constant in kcal/mol/rad^2 + """ + + atom_i: int + atom_j: int + atom_k: int + min_angle_deg: float + max_angle_deg: float + force_constant: float + + +@dataclass(frozen=True) +class TorsionConstraint: + """Harmonic torsion-window constraint between four atoms. + + Attributes: + atom_i: first atom index + atom_j: second atom index + atom_k: third atom index + atom_l: fourth atom index + min_dihedral_deg: lower dihedral bound in degrees + max_dihedral_deg: upper dihedral bound in degrees + force_constant: force constant in kcal/mol/rad^2 + """ + + atom_i: int + atom_j: int + atom_k: int + atom_l: int + min_dihedral_deg: float + max_dihedral_deg: float + force_constant: float + + @dataclass(frozen=True) class ConstraintModel: """Geometry constraints shared by all minimization paths. @@ -52,10 +127,21 @@ class ConstraintModel: Attributes: position_constraints: atom position constraints distance_constraints: atom-pair distance constraints + angle_constraints: atom-angle constraints + torsion_constraints: atom-torsion constraints """ position_constraints: tuple[PositionConstraint, ...] = () distance_constraints: tuple[DistanceConstraint, ...] = () + angle_constraints: tuple[AngleConstraint, ...] = () + torsion_constraints: tuple[TorsionConstraint, ...] = () + + @property + def has_constraints(self) -> bool: + """Return whether model contains any force-field constraints.""" + return bool( + self.position_constraints or self.distance_constraints or self.angle_constraints or self.torsion_constraints + ) @property def constrained_atoms(self) -> frozenset[int]: @@ -100,6 +186,110 @@ def from_atom_positions( ) return cls(position_constraints=tuple(constraints)) + @classmethod + def from_internal_coordinates( + cls, + mol: Chem.Mol, + bond_constraints: tuple[_BondConstraintSpecLike, ...] = (), + angle_constraints: tuple[_AngleConstraintSpecLike, ...] = (), + torsion_constraints: tuple[_TorsionConstraintSpecLike, ...] = (), + ) -> "ConstraintModel": + """Build internal-coordinate constraints, using first conformer as reference when needed. + + Args: + mol: molecule containing optional reference conformer + bond_constraints: bond constraint specifications + angle_constraints: angle constraint specifications + torsion_constraints: torsion constraint specifications + + Returns: + Constraint model with distance, angle, and torsion constraints + """ + if not bond_constraints and not angle_constraints and not torsion_constraints: + return cls.empty() + + conf = mol.GetConformer(mol.GetConformers()[0].GetId()) if mol.GetNumConformers() else None + distances: list[DistanceConstraint] = [] + angles: list[AngleConstraint] = [] + torsions: list[TorsionConstraint] = [] + + for constraint in bond_constraints: + distance = constraint.distance + if distance is None: + if conf is None: + msg = "Bond constraint without explicit distance requires molecule with reference conformer." + raise OpenConfValueError(msg) + distance = float( + conf.GetAtomPosition(int(constraint.atom_i)).Distance(conf.GetAtomPosition(int(constraint.atom_j))) + ) + distances.append( + DistanceConstraint( + atom_i=int(constraint.atom_i), + atom_j=int(constraint.atom_j), + min_distance=max(0.0, float(distance) - float(constraint.tolerance)), + max_distance=float(distance) + float(constraint.tolerance), + force_constant=float(constraint.force_constant), + ) + ) + + for constraint in angle_constraints: + angle = constraint.angle_deg + if angle is None: + if conf is None: + msg = "Angle constraint without explicit angle requires molecule with reference conformer." + raise OpenConfValueError(msg) + angle = float( + rdMolTransforms.GetAngleDeg( + conf, + int(constraint.atom_i), + int(constraint.atom_j), + int(constraint.atom_k), + ) + ) + angles.append( + AngleConstraint( + atom_i=int(constraint.atom_i), + atom_j=int(constraint.atom_j), + atom_k=int(constraint.atom_k), + min_angle_deg=float(angle) - float(constraint.tolerance_deg), + max_angle_deg=float(angle) + float(constraint.tolerance_deg), + force_constant=float(constraint.force_constant), + ) + ) + + for constraint in torsion_constraints: + dihedral = constraint.dihedral_deg + if dihedral is None: + if conf is None: + msg = "Torsion constraint without explicit dihedral requires molecule with reference conformer." + raise OpenConfValueError(msg) + dihedral = float( + rdMolTransforms.GetDihedralDeg( + conf, + int(constraint.atom_i), + int(constraint.atom_j), + int(constraint.atom_k), + int(constraint.atom_l), + ) + ) + torsions.append( + TorsionConstraint( + atom_i=int(constraint.atom_i), + atom_j=int(constraint.atom_j), + atom_k=int(constraint.atom_k), + atom_l=int(constraint.atom_l), + min_dihedral_deg=float(dihedral) - float(constraint.tolerance_deg), + max_dihedral_deg=float(dihedral) + float(constraint.tolerance_deg), + force_constant=float(constraint.force_constant), + ) + ) + + return cls( + distance_constraints=tuple(distances), + angle_constraints=tuple(angles), + torsion_constraints=tuple(torsions), + ) + @classmethod def from_metal_shell( cls, @@ -167,9 +357,25 @@ def combine(self, other: "ConstraintModel") -> "ConstraintModel": for c in other.distance_constraints } ) + angles = {(c.atom_i, c.atom_j, c.atom_k, c.min_angle_deg, c.max_angle_deg): c for c in self.angle_constraints} + angles.update( + {(c.atom_i, c.atom_j, c.atom_k, c.min_angle_deg, c.max_angle_deg): c for c in other.angle_constraints} + ) + torsions = { + (c.atom_i, c.atom_j, c.atom_k, c.atom_l, c.min_dihedral_deg, c.max_dihedral_deg): c + for c in self.torsion_constraints + } + torsions.update( + { + (c.atom_i, c.atom_j, c.atom_k, c.atom_l, c.min_dihedral_deg, c.max_dihedral_deg): c + for c in other.torsion_constraints + } + ) return ConstraintModel( position_constraints=tuple(positions[idx] for idx in sorted(positions)), distance_constraints=tuple(distances[key] for key in sorted(distances)), + angle_constraints=tuple(angles[key] for key in sorted(angles)), + torsion_constraints=tuple(torsions[key] for key in sorted(torsions)), ) def reset_positions(self, mol: Chem.Mol, conf_id: int) -> None: @@ -233,7 +439,7 @@ def add_constraints_to_force_field(ff: object, constraints: ConstraintModel, fam constraints: constraints to apply family: force-field family, either `"MMFF"` or `"UFF"` """ - if not constraints.position_constraints and not constraints.distance_constraints: + if not constraints.has_constraints: return position_method = getattr(ff, f"{family}AddPositionConstraint", None) @@ -252,3 +458,30 @@ def add_constraints_to_force_field(ff: object, constraints: ConstraintModel, fam float(constraint.max_distance), float(constraint.force_constant), ) + + angle_method = getattr(ff, f"{family}AddAngleConstraint", None) + if angle_method is not None: + for constraint in constraints.angle_constraints: + angle_method( + int(constraint.atom_i), + int(constraint.atom_j), + int(constraint.atom_k), + False, + float(constraint.min_angle_deg), + float(constraint.max_angle_deg), + float(constraint.force_constant), + ) + + torsion_method = getattr(ff, f"{family}AddTorsionConstraint", None) + if torsion_method is not None: + for constraint in constraints.torsion_constraints: + torsion_method( + int(constraint.atom_i), + int(constraint.atom_j), + int(constraint.atom_k), + int(constraint.atom_l), + False, + float(constraint.min_dihedral_deg), + float(constraint.max_dihedral_deg), + float(constraint.force_constant), + ) diff --git a/openconf/propose/hybrid.py b/openconf/propose/hybrid.py index 77644d8..d088a14 100644 --- a/openconf/propose/hybrid.py +++ b/openconf/propose/hybrid.py @@ -144,8 +144,18 @@ def __init__( if constraint_spec is not None else ConstraintModel.empty() ) - self.constraint_model = metal_constraints.combine(pose_constraints) - self._has_position_constraints = bool(self.constraint_model.position_constraints) + internal_constraints = ( + ConstraintModel.from_internal_coordinates( + mol, + bond_constraints=constraint_spec.bond_constraints, + angle_constraints=constraint_spec.angle_constraints, + torsion_constraints=constraint_spec.torsion_constraints, + ) + if constraint_spec is not None + else ConstraintModel.empty() + ) + self.constraint_model = metal_constraints.combine(pose_constraints).combine(internal_constraints) + self._has_constraints = self.constraint_model.has_constraints self.fast_minimizer = get_minimizer( config.minimizer, @@ -294,14 +304,19 @@ def generate_seeds(self, n_seeds: int, prune_rms_thresh: float | None = None) -> max_its = self.config.seed_minimization_iters nthreads = int(self.config.num_threads or 0) - if mmff_props is not None: + if mmff_props is not None and not self.constraint_model.has_constraints: minimize_start = time.perf_counter() energies = minimize_confs_mmff(self.mol, mmff_props, conf_ids, max_its, nthreads) self._add_time_stat("seed_minimization_time_s", time.perf_counter() - minimize_start) seed_results = list(zip(conf_ids, energies, strict=True)) else: minimize_start = time.perf_counter() - seed_results = [(cid, self._minimize_uff_single(self.mol, cid, max_its)) for cid in conf_ids] + seed_results = [ + (cid, self.fast_minimizer.minimize(self.mol, cid)) + if mmff_props is not None + else (cid, self._minimize_uff_single(self.mol, cid, max_its)) + for cid in conf_ids + ] self._add_time_stat("seed_minimization_time_s", time.perf_counter() - minimize_start) # Supplementary seeds covering cis/trans families that ETKDG undersamples @@ -540,14 +555,14 @@ def _select_move_type(self, step: int) -> str: forced = resolve_forced_move( step, self.config.shake_period, - constrained=self._has_position_constraints, + constrained=self._has_constraints, ) if forced is not None: return forced probs = resolve_move_probabilities( self._current_move_probs, - constrained=self._has_position_constraints, + constrained=self._has_constraints, has_ring_flips=bool(self.rotor_model.ring_flips), has_crankshaft=bool(self._moves.crankable_rings), has_kic=bool(self._moves.macro_kic_data), @@ -777,9 +792,9 @@ def propose_batch(self, pool: ConformerPool, step: int) -> list[tuple[int, float staging props (fast_dielectric applied), then transfers accepted (finite-energy) conformers back to self.mol. - When constraint_spec is set, falls back to sequential per-conformer - minimization with MMFF position restraints (MMFFOptimizeMoleculeConfs - does not support custom force field terms). + When constraints are set, falls back to sequential per-conformer + minimization with explicit force-field terms because + MMFFOptimizeMoleculeConfs does not support custom constraints. Args: pool: conformer pool for parent selection @@ -789,7 +804,7 @@ def propose_batch(self, pool: ConformerPool, step: int) -> list[tuple[int, float Accepted conformer IDs, energies, and sources """ # Constraint mode: per-conformer minimization with explicit restraints. - if self.constraint_spec is not None or self.constraint_model.position_constraints: + if self.constraint_model.has_constraints: results: list[tuple[int, float, str]] = [] for i in range(self.config.minimize_batch_size): result = self._propose_constrained(pool, step + i) @@ -891,7 +906,7 @@ def full_refine_final_constrained( Returns: Refined energies in kcal/mol aligned to `final_ids` """ - assert self.constraint_spec is not None + assert self.constraint_model.has_constraints # Keep only finals before refining final_set = set(final_ids) @@ -951,11 +966,13 @@ def run_hybrid_generation( stats = new_generation_stats() if config.collect_stats else {} constraint_spec = config.constraint_spec has_metal_input = any(_is_metal(atom) for atom in mol.GetAtoms()) and mol.GetNumConformers() > 0 - use_input_seed = constraint_spec is not None or has_metal_input + use_input_seed = ( + constraint_spec.requires_reference_geometry if constraint_spec is not None else False + ) or has_metal_input # Filter rotors before building the proposer so _rotor_angles is computed # only for free rotors. - if constraint_spec is not None: + if constraint_spec is not None and constraint_spec.constrained_rotor_atoms: rotor_model = filter_constrained_rotors(rotor_model, constraint_spec.constrained_atoms) effective_config, tuned_defaults_applied = _resolve_runtime_tuned_config(config, rotor_model) @@ -1126,7 +1143,7 @@ def run_hybrid_generation( if effective_config.do_final_refine: final_refine_start = time.perf_counter() - if constraint_spec is not None: + if proposer.constraint_model.has_constraints: final_energies = proposer.full_refine_final_constrained( mol, final_ids, effective_config.max_minimization_iters, dielectric=effective_config.final_dielectric ) diff --git a/tests/test_constrained.py b/tests/test_constrained.py index f9b53ce..6ed9b1f 100644 --- a/tests/test_constrained.py +++ b/tests/test_constrained.py @@ -357,6 +357,97 @@ def test_global_shake_suppressed_in_constrained_mode(): assert "global_shake" not in move_types, f"global_shake appeared in constrained mode: {move_types}" +# --------------------------------------------------------------------------- +# Internal-coordinate constraints +# --------------------------------------------------------------------------- + + +def test_internal_coordinate_constraints_freeze_reference_geometry(): + """Bond, angle, and torsion constraints can freeze input conformer values.""" + from rdkit.Chem import rdMolTransforms + + from openconf import ( + AngleConstraintSpec, + BondConstraintSpec, + ConformerConfig, + ConstraintSpec, + TorsionConstraintSpec, + generate_conformers, + ) + + mol = Chem.AddHs(Chem.MolFromSmiles("CCCC")) + AllChem.EmbedMolecule(mol, randomSeed=0) + conf = mol.GetConformer(0) + ref_distance = conf.GetAtomPosition(1).Distance(conf.GetAtomPosition(2)) + ref_angle = rdMolTransforms.GetAngleDeg(conf, 0, 1, 2) + ref_torsion = rdMolTransforms.GetDihedralDeg(conf, 0, 1, 2, 3) + + config = ConformerConfig( + max_out=3, + pool_max=10, + n_steps=8, + minimize_batch_size=1, + random_seed=0, + constraint_spec=ConstraintSpec( + bond_constraints=(BondConstraintSpec(1, 2),), + angle_constraints=(AngleConstraintSpec(0, 1, 2),), + torsion_constraints=(TorsionConstraintSpec(0, 1, 2, 3),), + ), + ) + ensemble = generate_conformers(mol, config=config, add_hs=False) + + assert ensemble.n_conformers > 0 + for record in ensemble.records: + trial_conf = ensemble.mol.GetConformer(record.conf_id) + distance = trial_conf.GetAtomPosition(1).Distance(trial_conf.GetAtomPosition(2)) + angle = rdMolTransforms.GetAngleDeg(trial_conf, 0, 1, 2) + torsion = rdMolTransforms.GetDihedralDeg(trial_conf, 0, 1, 2, 3) + + assert abs(distance - ref_distance) <= 0.02 + assert abs(angle - ref_angle) <= 1.0 + assert abs(((torsion - ref_torsion + 180.0) % 360.0) - 180.0) <= 1.0 + + +def test_explicit_internal_coordinate_constraints_do_not_require_reference_conformer(): + """Explicit bond, angle, and torsion targets work without input coordinates.""" + from rdkit.Chem import rdMolTransforms + + from openconf import ( + AngleConstraintSpec, + BondConstraintSpec, + ConformerConfig, + ConstraintSpec, + TorsionConstraintSpec, + generate_conformers, + ) + + config = ConformerConfig( + max_out=2, + pool_max=10, + n_seeds=3, + n_steps=0, + minimize_batch_size=1, + random_seed=0, + constraint_spec=ConstraintSpec( + bond_constraints=(BondConstraintSpec(1, 2, distance=1.54, tolerance=0.02),), + angle_constraints=(AngleConstraintSpec(0, 1, 2, angle_deg=112.0, tolerance_deg=2.0),), + torsion_constraints=(TorsionConstraintSpec(0, 1, 2, 3, dihedral_deg=180.0, tolerance_deg=2.0),), + ), + ) + ensemble = generate_conformers("CCCC", config=config) + + assert ensemble.n_conformers > 0 + for record in ensemble.records: + conf = ensemble.mol.GetConformer(record.conf_id) + distance = conf.GetAtomPosition(1).Distance(conf.GetAtomPosition(2)) + angle = rdMolTransforms.GetAngleDeg(conf, 0, 1, 2) + torsion = rdMolTransforms.GetDihedralDeg(conf, 0, 1, 2, 3) + + assert abs(distance - 1.54) <= 0.04 + assert abs(angle - 112.0) <= 3.0 + assert abs(((torsion - 180.0 + 180.0) % 360.0) - 180.0) <= 3.0 + + # --------------------------------------------------------------------------- # SDF output # ---------------------------------------------------------------------------