diff --git a/src/Auto3D/ASE/geometry.py b/src/Auto3D/ASE/geometry.py index 8ab14fbc..c09e938a 100644 --- a/src/Auto3D/ASE/geometry.py +++ b/src/Auto3D/ASE/geometry.py @@ -20,11 +20,10 @@ from Auto3D.torch_config import TorchConfig, configure_torch from Auto3D.utils.atomic_io import atomic_write_path from Auto3D.utils.energy import E_TOT_HARTREE_PROP, E_TOT_PROP +from Auto3D.utils.output_guard import check_output_not_input, check_output_overwrite from Auto3D.utils.validation import ( check_engine_supports_molecules, check_gpu_requested, - check_output_not_input, - check_output_overwrite, ) __all__ = ["opt_geometry"] diff --git a/src/Auto3D/ASE/thermo.py b/src/Auto3D/ASE/thermo.py index c535dbe6..9b0cc3dc 100644 --- a/src/Auto3D/ASE/thermo.py +++ b/src/Auto3D/ASE/thermo.py @@ -42,11 +42,10 @@ from Auto3D.torch_config import TorchConfig, configure_torch from Auto3D.utils.energy import hartree2ev from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_not_input, check_output_overwrite from Auto3D.utils.validation import ( check_engine_supports_molecules, check_gpu_requested, - check_output_not_input, - check_output_overwrite, ) __all__ = ["calc_thermo"] @@ -686,7 +685,11 @@ def vib_hessian(mol: Chem.Mol, ase_calculator, model, # mass weighting, and the rotational partition function silently disagree # for isotopically labeled input. atoms = mol2atoms(mol, positions=positions) - atoms.set_calculator(ase_calculator) + # atoms.set_calculator() is deprecated since ase 3.22.1 in favor of the + # `.calc` attribute (Minor 6); `pyproject.toml` pins no ase upper bound + # and globally ignores DeprecationWarning, so removal would otherwise + # land as a silent-until-runtime AttributeError with no advance warning. + atoms.calc = ase_calculator charge = rdmolops.GetFormalCharge(mol) # get the Hessian @@ -1194,7 +1197,9 @@ def do_mol_thermo(mol: Chem.Mol, # them via the explicit `positions=` argument, not from mol's conformer), # so nothing here depends on mol's conformer being in sync yet. coord = atoms.get_positions() - vib = vib_hessian(mol, atoms.get_calculator(), model, device, + # atoms.get_calculator() is deprecated since ase 3.22.1; use `.calc` + # (Minor 6, same rationale as the set_calculator() call above). + vib = vib_hessian(mol, atoms.calc, model, device, model_name=model_name, positions=coord) e = atoms.get_potential_energy() geometry = _detect_geometry(atoms) @@ -1691,7 +1696,9 @@ def calc_thermo(path: str, model_name: str, mol_info_func=None, atoms = mol2atoms(mol) calculator.set_charge(charge) - atoms.set_calculator(calculator) + # atoms.set_calculator() is deprecated since ase 3.22.1; use `.calc` + # (Minor 6, same rationale as vib_hessian's call above). + atoms.calc = calculator if mol_info_func is None: idx = mol.GetProp("_Name").strip() diff --git a/src/Auto3D/SPE.py b/src/Auto3D/SPE.py index 09ad68ca..043f5fe6 100644 --- a/src/Auto3D/SPE.py +++ b/src/Auto3D/SPE.py @@ -13,11 +13,10 @@ from Auto3D.models.preflight import resolve_engine_name from Auto3D.torch_config import TorchConfig, configure_torch from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_not_input, check_output_overwrite from Auto3D.utils.validation import ( check_engine_supports_molecules, check_gpu_requested, - check_output_not_input, - check_output_overwrite, ) logger = get_logger(__name__) diff --git a/src/Auto3D/batch_opt/optimization_engine.py b/src/Auto3D/batch_opt/optimization_engine.py index 289a8f2c..ec6f8426 100644 --- a/src/Auto3D/batch_opt/optimization_engine.py +++ b/src/Auto3D/batch_opt/optimization_engine.py @@ -272,7 +272,7 @@ def _step_active_subset( not_oscillating = oscillating_count < patience # Combine the convergence criteria. An `& ~energy_converged` term stood - # here until 4.0.0; `energy_converged` required `fmax < opttol` while + # here until 3.0.0; `energy_converged` required `fmax < opttol` while # `not_converged_post1` is `fmax > opttol`, so the term was the identity # of `&` wherever it was consulted and false-dominated elsewhere -- it # could never change an outcome, including at the `fmax == opttol` @@ -373,7 +373,7 @@ def n_steps( 2. Oscillation detection: drops structures that don't improve for patience steps There is deliberately no energy-stability criterion. One existed until - 4.0.0 but could never fire: it required ``fmax < opttol``, which is + 3.0.0 but could never fire: it required ``fmax < opttol``, which is exactly the condition under which the force criterion has already stopped the structure, so the term was the identity of ``&`` at every element (see ``test_convergence_outcome_never_depends_on_energy_stability``). diff --git a/src/Auto3D/cli/commands/config.py b/src/Auto3D/cli/commands/config.py index 59dd350c..52414a58 100644 --- a/src/Auto3D/cli/commands/config.py +++ b/src/Auto3D/cli/commands/config.py @@ -119,7 +119,7 @@ def execute_config_init( Both refusals below are configuration problems, so both leave through ``handle_error`` at exit **2**. The overwrite refusal in particular had to move: the CLI's own ``-o`` overwrite gate - (``utils.validation.check_output_overwrite``) raises ``ConfigurationError`` + (``utils.output_guard.check_output_overwrite``) raises ``ConfigurationError`` and exits 2, and the CHANGELOG described the two as printing "the same message" -- while this one hard-coded ``SystemExit(1)``, so a script branching on 2 saw one of them and not the other. diff --git a/src/Auto3D/cli/commands/properties.py b/src/Auto3D/cli/commands/properties.py index 24a73a60..ecac1ab9 100644 --- a/src/Auto3D/cli/commands/properties.py +++ b/src/Auto3D/cli/commands/properties.py @@ -18,11 +18,8 @@ from Auto3D.cli.errors import handle_error from Auto3D.exceptions import ConfigurationError, DependencyError from Auto3D.models.preflight import resolve_engine_name -from Auto3D.utils.validation import ( - check_gpu_requested, - check_output_not_input, - check_output_overwrite, -) +from Auto3D.utils.output_guard import check_output_not_input, check_output_overwrite +from Auto3D.utils.validation import check_gpu_requested # Engine names offered for shell completion. Free-form registry names and custom # model paths are also accepted -- each command below validates them with diff --git a/src/Auto3D/config.py b/src/Auto3D/config.py index 189eefef..22e2f51c 100644 --- a/src/Auto3D/config.py +++ b/src/Auto3D/config.py @@ -484,7 +484,7 @@ class OptimizationConfig: memory but may be faster.""" # There is deliberately no energy_tol/energy_patience here. Both existed - # until 4.0.0 and reached an optimizer criterion that could never fire + # until 3.0.0 and reached an optimizer criterion that could never fire # (audit M1), so they were knobs that changed nothing. Removed rather than # kept as inert configuration. diff --git a/src/Auto3D/constants.py b/src/Auto3D/constants.py index 5047c3b0..698a60f3 100644 --- a/src/Auto3D/constants.py +++ b/src/Auto3D/constants.py @@ -71,7 +71,7 @@ DEFAULT_PATIENCE = 250 # Steps before dropping oscillating conformer DEFAULT_BATCHSIZE_ATOMS = 1024 # Atoms per batch for GPU optimization DEFAULT_ENERGY_CLUSTER_WINDOW = 0.1 # eV, for RMSD clustering -# DEFAULT_ENERGY_TOL / DEFAULT_ENERGY_PATIENCE were removed in 4.0.0 along with +# DEFAULT_ENERGY_TOL / DEFAULT_ENERGY_PATIENCE were removed in 3.0.0 along with # the optimizer's energy-stability criterion, which could never fire (audit M1). # Tuning the tolerance -- for fp32 noise or anything else -- changed nothing. DEFAULT_RANDOM_SEED = 42 # Default random seed for reproducibility diff --git a/src/Auto3D/filtering.py b/src/Auto3D/filtering.py index b4af3b98..7bbe02db 100644 --- a/src/Auto3D/filtering.py +++ b/src/Auto3D/filtering.py @@ -1,7 +1,26 @@ #!/usr/bin/env python -"""Optimized conformer filtering with hierarchical RMSD comparison.""" +"""The one conformer filter: duplicate removal with hierarchical RMSD comparison. + +Auto3D carried two implementations of this until 3.0.0 -- this energy-clustered +one and a legacy all-pairs ``filter_unique``, selected by +``ConformerRanker(use_optimized_filtering=...)``. They applied the identical +duplicate criterion and were kept side by side so each could act as the other's +oracle, but a boolean kwarg that swaps one filter for another is two behaviors +to keep in step, and they had already drifted on malformed input: the legacy one +tolerated a record with no usable ``E_tot`` while this one raised ``KeyError`` +from its sort key. The survivor tolerates it (see :func:`_energy_sort_key`), and +the flag and the duplicate implementation are gone. + +The filter reports **why** it dropped things, not just what survived +(:class:`FilterResult`). ``ranking`` used to log "No structure converged" for a +species whose conformers were all dropped for *stereochemistry* -- a message +that sent the reader to the optimizer settings for a problem in the input's +stereo definitions. +""" from __future__ import annotations +from dataclasses import dataclass + from rdkit import Chem from rdkit.Chem import rdMolAlign @@ -12,67 +31,187 @@ ) from Auto3D.utils.connectivity import check_connectivity from Auto3D.utils.convergence import converged_or_unfiltered -from Auto3D.utils.energy import e_tot_ev, try_e_tot_ev +from Auto3D.utils.energy import try_e_tot_ev from Auto3D.utils.stereo_check import species_key, stereo_preserved +__all__ = [ + "DROP_REASONS", + "FilterResult", + "filter_conformers", + "filter_unique_optimized", +] + +#: Every reason a conformer can leave the selected set, in report order. +#: +#: The authoritative vocabulary: :class:`FilterResult` refuses a count keyed by +#: anything else, so a producer that invents a reason fails at construction +#: rather than contributing a silently unlabeled drop to a diagnostic. +#: ``energy_window`` is produced by ``ranking.ConformerRanker.top_window`` +#: rather than by this module -- it is a selection criterion, not a validity +#: one -- but it belongs in the same vocabulary because it reaches the user +#: through the same message. +DROP_REASONS: tuple[str, ...] = ( + "unparsed", + "unconverged", + "stereochemistry", + "connectivity", + "duplicate", + "energy_window", +) -def filter_unique_optimized( +#: Human-readable phrase per reason, used by :meth:`FilterResult.summary`. +_REASON_PHRASES: dict[str, str] = { + "unparsed": "unparseable by RDKit", + "unconverged": "marked Converged=false", + "stereochemistry": "changed stereochemistry during optimization", + "connectivity": "have broken or newly formed bonds", + "duplicate": "duplicates of a kept conformer", + "energy_window": "outside the energy window", +} + + +@dataclass(frozen=True) +class FilterResult: + """What survived filtering, and a count of what did not, by reason. + + Deliberately tiny: a list and a ``{reason: count}`` dict. The alternative + considered -- attaching a reason to each dropped molecule -- would keep + every rejected conformer alive for the duration of a chunk, which is the + memory the filter exists to release. + + Args: + kept: Surviving molecules, sorted by energy (lowest first); records + with no usable energy sort last. + dropped: Count per reason. Keys must come from :data:`DROP_REASONS`; + reasons that did not fire may be omitted or present as 0. + + Raises: + ValueError: ``dropped`` carries a key outside :data:`DROP_REASONS`. + """ + + kept: list[Chem.Mol] + dropped: dict[str, int] + + def __post_init__(self) -> None: + unknown = sorted(set(self.dropped) - set(DROP_REASONS)) + if unknown: + raise ValueError( + f"unknown filter drop reason(s) {unknown}; expected one of " + f"{list(DROP_REASONS)}" + ) + + @property + def reasons(self) -> tuple[str, ...]: + """The reasons that actually fired, in :data:`DROP_REASONS` order.""" + return tuple(r for r in DROP_REASONS if self.dropped.get(r)) + + def summary(self) -> str: + """One phrase per reason that fired, e.g. ``"2 marked Converged=false"``. + + Empty string when nothing was dropped, so a caller can test it for + truth rather than special-casing a placeholder. + """ + return ", ".join( + f"{self.dropped[reason]} {_REASON_PHRASES[reason]}" + for reason in self.reasons + ) + + +def _energy_sort_key(mol: Chem.Mol) -> tuple[bool, float]: + """Sort key that tolerates a record with no usable ``E_tot``. + + Returns ``(False, energy_ev)`` for a record that has one and + ``(True, 0.0)`` for a record that does not, so the energy-less records land + **after** every record that has an energy, whatever their values. + + The 0.0 is a tie-break placeholder, never an energy: reading a missing + ``E_tot`` as 0.0 on its own would sort a garbage record ahead of every + genuine structure (real ``E_tot`` values are large and negative), making it + the reference conformer ``E_rel`` is measured from and the single structure + a ``k=1`` request returns. Among themselves, energy-less records keep their + input order (``list.sort`` is stable), which is the only ordering there is + any evidence for. + """ + energy = try_e_tot_ev(mol) + if energy is None: + return (True, 0.0) + return (False, energy) + + +def filter_conformers( mols: list[Chem.Mol], + *, rmsd_threshold: float = DEFAULT_RMSD_THRESHOLD, energy_cluster_window: float = DEFAULT_ENERGY_CLUSTER_WINDOW, -) -> list[Chem.Mol]: +) -> FilterResult: """Remove duplicate conformers, skipping RMSD comparisons that cannot match. - Sorts by energy and only RMSD-compares molecules close enough in energy to be - duplicates at all, which avoids the O(n^2) comparisons of the legacy - :func:`filter_unique` below **without changing which molecules - survive** -- the partitioning is chosen so that no duplicate pair can be - separated by it. See the comment on the split rule below for why that - holds; it did not hold before 4.0.0. + Sorts by energy and only RMSD-compares molecules close enough in energy to + be duplicates at all, which avoids comparing all pairs **without changing + which molecules survive** -- the partitioning is chosen so that no duplicate + pair can be separated by it. See the comment on the split rule below for why + that holds; it did not hold before 3.0.0. A pair counts as a duplicate only when all three of these agree: the two are the same compound (:func:`Auto3D.utils.stereo_check.species_key`), their heavy-atom RMSD is under ``rmsd_threshold``, and their energies agree within - ``DEFAULT_DUPLICATE_ENERGY_TOL``. + ``DEFAULT_DUPLICATE_ENERGY_TOL`` (or at least one of them has no usable + energy, in which case that term cannot apply and RMSD alone decides). Args: mols: List of RDKit Mol objects with 'E_tot' (Hartree) and, optionally, 'Converged' properties. A record whose 'Converged' property is explicitly false is dropped; a record without the property is kept (it is not filtered on convergence). Records marked - 'Stereo_changed' are excluded. + 'Stereo_changed' are excluded. ``None`` entries -- what + ``SDMolSupplier`` yields for a record RDKit cannot parse -- are + counted and skipped. rmsd_threshold: RMSD threshold for considering structures similar (Angstrom). energy_cluster_window: Energy width (eV) below which molecules are compared to each other. A performance knob only: values below the duplicate energy tolerance cannot shrink the comparison set, because that tolerance is the floor at which a pair can still be a duplicate. The stored Hartree energies are converted to eV on read - (``Auto3D.utils.energy.e_tot_ev``), so the unit is as documented. + (``Auto3D.utils.energy``), so the unit is as documented. Returns: - List of unique molecules, sorted by energy (lowest first). + A :class:`FilterResult` whose ``kept`` list is sorted by energy (lowest + first) and whose ``dropped`` counts say why the rest are missing. """ + dropped: dict[str, int] = {} + + def _drop(reason: str) -> None: + dropped[reason] = dropped.get(reason, 0) + 1 + # Filter converged structures with valid connectivity. A record with no # 'Converged' property is not filtered on convergence (see # Auto3D.utils.convergence): treating its absence as failure deleted every # record of any SDF batchopt did not write. - valid_mols = [] + # + # Checked in this order, one reason attributed per record, so a structure + # that fails several is reported under the first -- the same short-circuit + # order the single `and` chain here used to have, hence the same verdicts. + valid_mols: list[Chem.Mol] = [] for mol in mols: if mol is None: - continue - if ( - converged_or_unfiltered(mol) - and stereo_preserved(mol) - and check_connectivity(mol) - ): + _drop("unparsed") + elif not converged_or_unfiltered(mol): + _drop("unconverged") + elif not stereo_preserved(mol): + _drop("stereochemistry") + elif not check_connectivity(mol): + _drop("connectivity") + else: valid_mols.append(mol) if not valid_mols: - return [] + return FilterResult(kept=[], dropped=dropped) # Sort by energy. E_tot is stored in Hartree; energy_cluster_window and # the duplicate tolerance below are both in eV, so convert on read. - valid_mols.sort(key=e_tot_ev) + # Records with no usable energy sort last -- see _energy_sort_key. + valid_mols.sort(key=_energy_sort_key) + energies = [try_e_tot_ev(mol) for mol in valid_mols] # Partition the energy-sorted list into runs, and only RMSD-compare within a # run. Where a run may end is a correctness question, not a tuning one. @@ -97,22 +236,34 @@ def filter_unique_optimized( # splitting on LARGER gaps is always safe (it only merges runs and compares # more pairs). It therefore stays a performance knob and can no longer become # a correctness hole. Runs are no longer width-bounded, so a dense energy - # ladder degrades to the O(n^2) of the legacy `filter_unique` -- the price of - # the guarantee, on conformer counts that are tens per species. - split_gap = max(DEFAULT_DUPLICATE_ENERGY_TOL, energy_cluster_window) - clusters: list[list[Chem.Mol]] = [] - current_cluster: list[Chem.Mol] = [valid_mols[0]] - previous_e = e_tot_ev(valid_mols[0]) - - for mol in valid_mols[1:]: - e = e_tot_ev(mol) - if e - previous_e <= split_gap: - current_cluster.append(mol) - else: - clusters.append(current_cluster) - current_cluster = [mol] - previous_e = e - clusters.append(current_cluster) + # ladder degrades to comparing all pairs -- the price of the guarantee, on + # conformer counts that are tens per species. + # + # That entire argument rests on every record HAVING an energy. As soon as one + # does not, no gap proves anything about it: `_filter_within_cluster` falls + # back to RMSD-only for a pair where either side's energy is missing, so such + # a record can be a duplicate of any same-species structure at any energy. + # The honest response is to stop partitioning and compare all pairs, which is + # exactly what the legacy all-pairs filter did with this input -- so the + # survivor keeps its verdicts. Malformed input is rare and reaches here only + # through a direct API call (`ConformerRanker` refuses a record with no + # 'E_tot' up front), so the quadratic cost is paid where it is warranted. + if any(energy is None for energy in energies): + clusters: list[list[Chem.Mol]] = [valid_mols] + else: + split_gap = max(DEFAULT_DUPLICATE_ENERGY_TOL, energy_cluster_window) + clusters = [] + current_cluster: list[Chem.Mol] = [valid_mols[0]] + previous_e = energies[0] + + for mol, e in zip(valid_mols[1:], energies[1:], strict=True): + if e - previous_e <= split_gap: + current_cluster.append(mol) + else: + clusters.append(current_cluster) + current_cluster = [mol] + previous_e = e + clusters.append(current_cluster) # Filter unique within each cluster unique_mols: list[Chem.Mol] = [] @@ -120,7 +271,32 @@ def filter_unique_optimized( unique_in_cluster = _filter_within_cluster(cluster, rmsd_threshold) unique_mols.extend(unique_in_cluster) - return unique_mols + n_duplicates = len(valid_mols) - len(unique_mols) + if n_duplicates: + dropped["duplicate"] = n_duplicates + return FilterResult(kept=unique_mols, dropped=dropped) + + +def filter_unique_optimized( + mols: list[Chem.Mol], + rmsd_threshold: float = DEFAULT_RMSD_THRESHOLD, + energy_cluster_window: float = DEFAULT_ENERGY_CLUSTER_WINDOW, +) -> list[Chem.Mol]: + """The surviving molecules from :func:`filter_conformers`, nothing else. + + Kept as the public name it has always been, for callers that want the list + and not the drop counts. Prefer :func:`filter_conformers` when the reason + something is missing has to reach a user. + + Returns: + List of unique molecules, sorted by energy (lowest first); records with + no usable energy sort last. + """ + return filter_conformers( + mols, + rmsd_threshold=rmsd_threshold, + energy_cluster_window=energy_cluster_window, + ).kept def _filter_within_cluster( @@ -206,105 +382,3 @@ def _mol_energy(mol: Chem.Mol) -> float | None: """ return try_e_tot_ev(mol) - -def filter_unique(mols: list[Chem.Mol], crit: float = DEFAULT_RMSD_THRESHOLD) -> list[Chem.Mol]: - """Remove structures that are very similar and remove unconverged structures. - - The legacy all-pairs filter, kept beside :func:`filter_unique_optimized` - (which supersedes it) because ``ConformerRanker(use_optimized_filtering=False)`` - still selects it and because it is the oracle the optimized filter is - compared against. - - This function filters a list of molecules to keep only unique, converged structures. - It first removes unconverged structures and those with invalid connectivity, - then removes similar structures based on RMSD comparison. - - Args: - mols: List of RDKit molecule objects, optionally carrying a 'Converged' - property. A record whose 'Converged' is explicitly false is - dropped; a record without the property is kept (not filtered on - convergence). Records marked 'Stereo_changed' are excluded. - crit: RMSD threshold for considering two structures as identical. - Structures with RMSD below this value are considered duplicates. - Defaults to DEFAULT_RMSD_THRESHOLD (0.3 Angstroms). - - Returns: - List of unique, converged molecules with valid connectivity. - - Example: - >>> from rdkit import Chem - >>> from rdkit.Chem import AllChem - >>> mol = Chem.MolFromSmiles("CCO") - >>> mol = Chem.AddHs(mol) - >>> AllChem.EmbedMolecule(mol, randomSeed=42) - 0 - >>> mol.SetProp("Converged", "true") - >>> filter_unique([mol], crit=0.3) # Returns list with 1 molecule - [...] - """ - # Remove structures that explicitly failed to converge. A record with no - # 'Converged' property is NOT filtered on convergence -- see - # Auto3D.utils.convergence for why absence is not failure. - mols_: list[Chem.Mol] = [] - for mol in mols: - convergence_flag = converged_or_unfiltered(mol) - has_valid_bonds = check_connectivity(mol) - if convergence_flag and has_valid_bonds and stereo_preserved(mol): - mols_.append(mol) - mols = mols_ - - # Remove similar structures. Strip Hs once per molecule (O(n)) instead of on - # both sides of every comparison (O(n^2)); GetBestRMS on no-H forms is - # symmetric so results are unchanged. The ORIGINAL (H-explicit) mols are - # returned; no-H forms are comparison-only. - # - # Heavy-atom RMSD alone collapses conformers that differ only in an O-H / N-H - # rotor orientation. Guard with an energy check: a pair counts as duplicate - # only when the RMSD is below ``crit`` AND the energies agree within - # DEFAULT_DUPLICATE_ENERGY_TOL (eV; 'E_tot' is stored in Hartree and is - # converted on read by Auto3D.utils.energy). Mols without a usable 'E_tot' - # fall back to RMSD-only (energy guard cannot apply). - unique_mols: list[Chem.Mol] = [] - unique_noH: list[Chem.Mol] = [] - unique_energies: list[float | None] = [] - unique_species: list[str] = [] - for mol_i in mols: - mol_i_noH = Chem.RemoveHs(mol_i) - # E_tot is stored in Hartree; DEFAULT_DUPLICATE_ENERGY_TOL is in eV. - e_i: float | None = try_e_tot_ev(mol_i) - species_i = species_key(mol_i) - unique = True - for mol_j_noH, e_j, species_j in zip( - unique_noH, unique_energies, unique_species, strict=True - ): - # Two different compounds are never duplicates of each other, however - # close their geometries. All stereoisomers of one input share a - # ranking group, and two ring diastereomers can sit below the default - # 0.3 A threshold, so without this the RMSD test could delete one of - # them (Auto3D.utils.stereo_check.species_key). Checked before the - # RMSD call it makes unnecessary. - if species_i != species_j: - continue - try: - # temporary bug fix for https://github.com/rdkit/rdkit/issues/6826 - # removing Hs speeds up the calculation - rmsd = rdMolAlign.GetBestRMS(mol_i_noH, mol_j_noH) - except RuntimeError: - # Incomparable pair: treat as distinct (not a duplicate) so the - # conformer is kept. Using 0 would make it look like a perfect - # duplicate and drop a genuinely distinct structure. - rmsd = float("inf") - energy_close = ( - e_i is None - or e_j is None - or abs(e_i - e_j) < DEFAULT_DUPLICATE_ENERGY_TOL - ) - if rmsd < crit and energy_close: - unique = False - break - if unique: - unique_mols.append(mol_i) - unique_noH.append(mol_i_noH) - unique_energies.append(e_i) - unique_species.append(species_i) - return unique_mols diff --git a/src/Auto3D/id_mapping.py b/src/Auto3D/id_mapping.py index 7b2c45e4..028cf96d 100644 --- a/src/Auto3D/id_mapping.py +++ b/src/Auto3D/id_mapping.py @@ -16,15 +16,19 @@ from rdkit import Chem -from Auto3D.exceptions import ConfigurationError, InputValidationError +from Auto3D.exceptions import InputValidationError from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_overwrite from Auto3D.utils.smi_io import iter_smi_records logger = get_logger(__name__) def encode_ids( - path: str, out_dir: str | os.PathLike[str] | None = None + path: str, + out_dir: str | os.PathLike[str] | None = None, + *, + overwrite: bool = False, ) -> tuple[str, dict[str, int]]: """Encode molecule IDs to numeric indices. @@ -46,6 +50,13 @@ def encode_ids( path: Path to the input .smi or .sdf file. out_dir: Directory to write the encoded file into. Defaults to the input file's own directory. + overwrite: Allow replacing an existing file at the encoded path. + Keyword-only, defaults to False -- Auto3D derives this name, so it + belongs to whoever already owns it. This used to refuse + *unconditionally*, with no way for a caller who meant it to say so; + the keyword and its default now match ``decode_ids`` and + ``tautomer.select_tautomers``, and all three raise through the same + ``utils.output_guard.check_output_overwrite``. Returns: Tuple containing: @@ -54,7 +65,8 @@ def encode_ids( Raises: ValueError: If the input file is neither .smi nor .sdf format. - ConfigurationError: If a file already exists at the encoded path. + ConfigurationError: If a file already exists at the encoded path and + `overwrite` is False. InputValidationError: If a molecule has a missing/blank ID or a duplicate ID is encountered. @@ -72,13 +84,7 @@ def encode_ids( directory = Path(out_dir) if out_dir is not None else path_obj.parent new_path = directory / f"{path_obj.stem}_encoded.{extension}" - if new_path.exists(): - raise ConfigurationError( - f"encode_ids would overwrite the existing file {new_path}. " - "Auto3D writes its encoded copy of the input there; move or " - "rename that file, or pass out_dir to write the encoded copy " - "somewhere else." - ) + check_output_overwrite(new_path, overwrite) if extension == "smi": new_data: list[str] = [] @@ -129,7 +135,9 @@ def encode_ids( return str(new_path), mapping -def decode_ids(path: str, mapping: dict[str, int]) -> str: +def decode_ids( + path: str, mapping: dict[str, int], *, overwrite: bool = False +) -> str: """Decode numeric IDs back to original molecule IDs. For an SDF file with numeric IDs, restores the original IDs using @@ -139,10 +147,21 @@ def decode_ids(path: str, mapping: dict[str, int]) -> str: path: Path to the input SDF file with encoded (numeric) IDs. mapping: Dictionary mapping original IDs to their numeric indices (as returned by encode_ids). + overwrite: Allow replacing an existing file at the decoded path. + Keyword-only, defaults to False: ``_out.`` is a name + Auto3D derives from `path`, so an existing file there belongs to + someone else, and ``Chem.SDWriter`` truncates on open. + ``WorkflowOrchestrator`` calls this on the combined output inside a + job directory it created with a bare ``mkdir()``, where the derived + name is always free. Returns: Path to the new SDF file with decoded IDs (adds '_out' suffix). + Raises: + ConfigurationError: A file already exists at the decoded path and + `overwrite` is False. + Example: >>> mapping = {'mol_A': 0, 'mol_B': 1} >>> output_path = decode_ids("encoded_3d.sdf", mapping) @@ -155,6 +174,8 @@ def decode_ids(path: str, mapping: dict[str, int]) -> str: stem_parts = path_obj.stem.split("_")[:-2] new_stem = "_".join(stem_parts) + "_out" new_path = path_obj.parent / f"{new_stem}.{extension}" + # Before the supplier is opened, so a refusal costs nothing. + check_output_overwrite(new_path, overwrite) suppl = Chem.SDMolSupplier(path, removeHs=False) with Chem.SDWriter(str(new_path)) as w: diff --git a/src/Auto3D/ranking.py b/src/Auto3D/ranking.py index 694513f0..b10086d9 100644 --- a/src/Auto3D/ranking.py +++ b/src/Auto3D/ranking.py @@ -6,8 +6,9 @@ from rdkit import Chem from Auto3D.config import SELECTOR_FIELDS, check_selectors_mutually_exclusive +from Auto3D.constants import DEFAULT_ENERGY_CLUSTER_WINDOW from Auto3D.exceptions import ConfigurationError, InputValidationError -from Auto3D.filtering import filter_unique, filter_unique_optimized +from Auto3D.filtering import FilterResult, filter_conformers from Auto3D.utils.connectivity import check_connectivity from Auto3D.utils.convergence import converged_or_unfiltered, has_convergence_flag from Auto3D.utils.energy import ( @@ -59,6 +60,59 @@ def species_id(name: str) -> str: return name.strip().rsplit("_", 2)[0].strip() +#: Selector field (from :data:`Auto3D.config.SELECTOR_FIELDS`) -> the +#: ``ConformerRanker`` method that implements it. +#: +#: ``run`` used to dispatch with a hand-written ``if self.k: ... elif +#: self.window: ...``, which meant a third selector added to ``SELECTOR_FIELDS`` +#: was accepted by ``Auto3DOptions``, accepted by ``CLIConfig``, accepted by +#: ``check_selectors_mutually_exclusive`` -- and then silently ignored here, +#: falling through to the "Parameter k or window needs to be specified" error +#: even though the user had specified one. The registry plus the import-time +#: equality check below turns that into an ImportError at the moment the field +#: list and this table disagree. +_SELECTORS: dict[str, str] = {"k": "top_k", "window": "top_window"} + + +def _verify_selector_registry( + registry: dict[str, str], fields: tuple[str, ...], cls: type +) -> None: + """Raise unless ``registry`` covers ``fields`` with methods that exist. + + Called at import time, and a plain function rather than inline module code + so the failure it guards against can be tested by calling it with a bad + registry -- instead of reloading the module and hoping the reload's own + machinery does not mask the error. + + Args: + registry: The selector-field -> method-name table (:data:`_SELECTORS`). + fields: The authoritative field list (``config.SELECTOR_FIELDS``). + cls: The class expected to implement each mapped method. + + Raises: + ImportError: The registry's field set differs from ``fields``, or a + mapped method is missing from ``cls``. + """ + if set(registry) != set(fields): + raise ImportError( + "Auto3D.ranking._SELECTORS is out of step with " + f"Auto3D.config.SELECTOR_FIELDS: registry has {sorted(registry)}, " + f"config declares {sorted(fields)}. config.py owns the field list; " + "add the missing selector here and give ConformerRanker a method " + "that implements it (a selector nothing dispatches on is accepted " + "by every config class and then silently ignored)." + ) + for field, method in registry.items(): + # Checked as well as the names, because a typo in a method name passes + # the set comparison above and then raises AttributeError from inside + # `run` -- after the whole input has been read and grouped. + if not callable(getattr(cls, method, None)): + raise ImportError( + f"Auto3D.ranking._SELECTORS maps selector {field!r} to " + f"{method!r}, which is not a method of {cls.__name__}." + ) + + class ConformerRanker: """Select conformers that satisfy user-defined energy criteria. @@ -72,12 +126,12 @@ class ConformerRanker: k: Select top-k structures per SMILES. False to disable. window: Select structures within window kcal/mol of minimum. False to disable. - use_optimized_filtering: If True (default), use energy-clustered - RMSD filtering which has O(n*k) complexity instead of O(n^2). - Set to False to use legacy O(n^2) filtering. - energy_cluster_window: Energy window in eV for clustering molecules - before RMSD comparison. Only used when use_optimized_filtering - is True. Default is 0.1 eV. + energy_cluster_window: Energy width in eV below which molecules are + RMSD-compared to each other. A performance knob only -- see + ``Auto3D.filtering.filter_conformers``, which is the single + conformer filter now that ``use_optimized_filtering`` and the + legacy all-pairs implementation it selected are gone. Defaults to + ``Auto3D.constants.DEFAULT_ENERGY_CLUSTER_WINDOW``. overwrite: Allow writing over an existing `out_path`. Defaults to True, which is the historical behavior every caller was written against -- including Auto3D's own pipeline, which always writes @@ -92,8 +146,7 @@ def __init__( threshold: float, k: int | bool = False, window: float | bool = False, - use_optimized_filtering: bool = True, - energy_cluster_window: float = 0.1, + energy_cluster_window: float = DEFAULT_ENERGY_CLUSTER_WINDOW, overwrite: bool = True, ) -> None: # Same C14 guard the three API entry points run. ConformerRanker is a @@ -117,29 +170,62 @@ def __init__( self.threshold = threshold self.k = k self.window = window - self.use_optimized_filtering = use_optimized_filtering self.energy_cluster_window = energy_cluster_window - - def _filter_mols(self, mols: list[Chem.Mol]) -> list[Chem.Mol]: + #: Run-level drop tally, keyed by ``Auto3D.filtering.DROP_REASONS``. + #: ``top_k``/``top_window`` add their per-species counts here (see + #: ``_account``) so ``run``'s "selected 0 structures" warning can state + #: the actual reasons instead of listing the ones it might have been. + #: Reset at the top of ``run``; harmless (and still accurate) for a + #: caller who invokes ``top_k``/``top_window`` directly. + self._drop_totals: dict[str, int] = {} + + def _account(self, result: FilterResult) -> FilterResult: + """Add ``result``'s drop counts to this run's tally, and return it.""" + for reason, count in result.dropped.items(): + if count: + self._drop_totals[reason] = self._drop_totals.get(reason, 0) + count + return result + + def _filter_mols(self, mols: list[Chem.Mol]) -> FilterResult: """Filter molecules to remove duplicates based on RMSD. - Uses either the optimized energy-clustered approach or the legacy - O(n^2) comparison, depending on the use_optimized_filtering setting. - Args: mols: List of RDKit Mol objects to filter. Returns: - List of unique molecules, sorted by energy. + A ``FilterResult``: unique molecules sorted by energy, plus the + per-reason drop counts the "nothing survived" messages below need. """ - if self.use_optimized_filtering: - return filter_unique_optimized( - mols, - rmsd_threshold=self.threshold, - energy_cluster_window=self.energy_cluster_window, - ) + return filter_conformers( + mols, + rmsd_threshold=self.threshold, + energy_cluster_window=self.energy_cluster_window, + ) + + def _log_nothing_selected(self, name: str, result: FilterResult) -> None: + """Say WHY a species contributed no structure to the output. + + This used to be an unconditional "No structure converged for X." -- + which is what a reader saw when every conformer of X was dropped for + *stereochemistry* (an optimization that inverted a center, so the + geometry no longer matches the title) or for *connectivity* (a + structure that fell apart). Both point at the input or the chemistry, + and both were reported as an optimizer convergence problem, sending the + reader to ``--opt-steps`` and ``--convergence-threshold`` for something + neither could fix. + + The literal wording survives for the one case where it was true -- + convergence being the sole reason -- because it is the message users + and their log-scraping scripts have matched on since Auto3D 1.x. When + the sole reason is convergence there is nothing to add; when it is not, + the reasons are named. + """ + if set(result.dropped) <= {"unconverged"}: + logger.info(f"No structure converged for {name}.") else: - return filter_unique(mols, crit=self.threshold) + logger.info( + "No structure selected for %s: %s.", name, result.summary() + ) def top_k(self, df_group: pd.DataFrame, k: int = 1) -> list[Chem.Mol]: """Select top-k lowest-energy structures from a group. @@ -159,28 +245,43 @@ def top_k(self, df_group: pd.DataFrame, k: int = 1) -> list[Chem.Mol]: df2 = df_group.sort_values(by=['energies']) # Optimization: when k=1, skip the expensive RMSD dedup but still - # apply connectivity validation. Return the lowest-energy conformer - # that passes check_connectivity (no broken/formed bonds); if none - # pass, return an empty list. + # apply the validity checks. Return the lowest-energy conformer that + # passes them; if none pass, return an empty list. + # + # The predicates and the drop-reason names are deliberately the same + # ones `filter_conformers` uses, so the diagnostic a user gets for an + # empty selection does not depend on which k they asked for. (Dedup is + # the only thing skipped -- it cannot empty a non-empty set anyway.) + # The loop stops at the first survivor, so the counts are partial when + # something IS selected and complete in exactly the case that produces + # a message. if k == 1: out_mols = [] + dropped: dict[str, int] = {} for mol in df2["mols"]: - if stereo_preserved(mol) and check_connectivity(mol): + if not converged_or_unfiltered(mol): + reason = "unconverged" + elif not stereo_preserved(mol): + reason = "stereochemistry" + elif not check_connectivity(mol): + reason = "connectivity" + else: out_mols = [mol] break + dropped[reason] = dropped.get(reason, 0) + 1 + result = self._account(FilterResult(kept=out_mols, dropped=dropped)) else: - out_mols_ = self._filter_mols(list(df2["mols"])) - if k < len(out_mols_): - out_mols = out_mols_[:k] - else: - out_mols = out_mols_ + result = self._account(self._filter_mols(list(df2["mols"]))) + # Truncation to k is selection, not a filter drop: those conformers + # are valid and unique, they just lost the ranking, so they are not + # counted among `dropped` (nothing is "missing" to explain). + out_mols = result.kept[:k] if k < len(result.kept) else result.kept if len(out_mols) == 0: # names[0] is already the group's species id (see species_id() # above) -- no further splitting here, or a disambiguated id # like "KEY_2" would misreport as "KEY" in this message. - name = names[0].strip() - logger.info(f"No structure converged for {name}.") + self._log_nothing_selected(names[0].strip(), result) else: #Adding relative energies # E_tot is stored in Hartree; E_rel(eV) is, as its name says, eV. @@ -210,21 +311,20 @@ def top_window(self, df_group: pd.DataFrame, window: float = 1.0) -> list[Chem.M raise ValueError(f"All molecules must have the same name, got: {set(names)}") df2 = df_group.sort_values(by=['energies']) - out_mols_ = self._filter_mols(list(df2['mols'])) + result = self._account(self._filter_mols(list(df2['mols']))) out_mols = [] - if len(out_mols_) == 0: + if len(result.kept) == 0: # names[0] is already the group's species id -- see the note in # top_k above. - name = names[0].strip() - logger.info(f"No structure converged for {name}.") + self._log_nothing_selected(names[0].strip(), result) else: # `window` was converted from kcal/mol to eV above, and E_tot is # stored in Hartree, so both sides of the comparison are eV here. # Reading the Hartree number as if it were eV is what made the # window 27.2x too wide for an opt_geometry-produced input. - ref_energy = e_tot_ev(out_mols_[0]) - for mol in out_mols_: + ref_energy = e_tot_ev(result.kept[0]) + for mol in result.kept: my_energy = e_tot_ev(mol) rel_energy = my_energy - ref_energy if rel_energy <= window: @@ -232,6 +332,16 @@ def top_window(self, df_group: pd.DataFrame, window: float = 1.0) -> list[Chem.M out_mols.append(mol) else: break + # The window is the one drop reason this method owns, and it is + # merged into the same tally as the filter's own counts rather than + # reported separately, so `run`'s summary is one accounting of the + # whole selection. `break` above is safe because `kept` ascends in + # energy, so every remaining conformer is outside the window too. + n_outside = len(result.kept) - len(out_mols) + if n_outside: + self._account( + FilterResult(kept=out_mols, dropped={"energy_window": n_outside}) + ) return out_mols def run(self) -> list[Chem.Mol]: @@ -271,9 +381,13 @@ def run(self) -> list[Chem.Mol]: {name: getattr(self, name, None) for name in SELECTOR_FIELDS} ) results = [] + # Fresh tally per run, so a ranker reused for a second file does not + # report the first file's drops (see _account). + self._drop_totals = {} mols, names, energies = [], [], [] n_records = 0 + n_unparsed = 0 n_unconverged = 0 n_unflagged = 0 # Context-managed so the SDF file handle is released promptly rather than @@ -281,6 +395,7 @@ def run(self) -> list[Chem.Mol]: with Chem.SDMolSupplier(self.input_path, removeHs=False) as supplier: for position, mol in enumerate(supplier): if mol is None: + n_unparsed += 1 logger.warning( "Skipping record %d of %s: RDKit could not parse it.", position, self.input_path, @@ -327,10 +442,18 @@ def run(self) -> list[Chem.Mol]: for group_name in groups.indices: group = groups.get_group(group_name) - if self.k: - top_results = self.top_k(group, self.k) - elif self.window: - top_results = self.top_window(group, self.window) + # Dispatch through the registry rather than a hand-written + # if/elif chain on the field names, so the set of selectors this + # method honors is the set `Auto3D.config` declares -- checked for + # equality at import (see _SELECTORS). Iterating SELECTOR_FIELDS + # rather than the dict also means config.py owns the PRECEDENCE, + # not this module's dict-literal order; the mutual-exclusivity + # check above has already refused more than one anyway. + for field in SELECTOR_FIELDS: + value = getattr(self, field, None) + if value: + top_results = getattr(self, _SELECTORS[field])(group, value) + break else: raise ConfigurationError('Parameter k or window needs to be ' 'specified. Append "--k=1" if you' @@ -342,12 +465,28 @@ def run(self) -> list[Chem.Mol]: # returns []. WARNING, not INFO: `logging.lastResort` puts WARNING # and above on stderr even for a caller who never ran # configure_logging, which is every direct API caller. + # + # The reasons are the ones that actually fired. This used to read + # "N record(s) are marked Converged=false and the rest were dropped + # by the connectivity, stereochemistry or energy-window filters" -- + # a hand-maintained list of everything it MIGHT have been, which + # left the reader to work out which, and which had to be edited by + # hand every time a filter was added. `_drop_totals` is the tally + # `top_k`/`top_window` fed as they ran; the two counts below are + # the drops `run` itself made, before grouping, which never reach a + # filter and so are not double-counted. + totals = dict(self._drop_totals) + for reason, count in ( + ("unparsed", n_unparsed), ("unconverged", n_unconverged), + ): + if count: + totals[reason] = totals.get(reason, 0) + count + summary = FilterResult(kept=[], dropped=totals).summary() logger.warning( "Selected 0 structures from %d record(s) in %s, so %s is " - "empty: %d record(s) are marked Converged=false and the rest " - "were dropped by the connectivity, stereochemistry or " - "energy-window filters.", - n_records, self.input_path, self.out_path, n_unconverged, + "empty: %s.", + n_records, self.input_path, self.out_path, + summary or "no filter reported a reason", ) with Chem.SDWriter(self.out_path) as f: @@ -372,5 +511,11 @@ def run(self) -> list[Chem.Mol]: return results +# Runs at import: a selector declared in config.py with nothing wired to it +# here, or wired to a method that does not exist, fails the import rather than +# being discovered by a user whose `--k=1` was silently ignored. +_verify_selector_registry(_SELECTORS, SELECTOR_FIELDS, ConformerRanker) + + # Backward compatibility alias ranking = ConformerRanker diff --git a/src/Auto3D/tautomer.py b/src/Auto3D/tautomer.py index 438941e1..a38add99 100644 --- a/src/Auto3D/tautomer.py +++ b/src/Auto3D/tautomer.py @@ -10,13 +10,20 @@ from Auto3D.exceptions import ConfigurationError from Auto3D.utils.energy import e_tot_hartree, hartree2kcalpermol from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_overwrite logger = get_logger(__name__) __all__ = ["select_tautomers", "get_stable_tautomers"] -def select_tautomers(sdf: str, k: int | None = None, window: float | None = None) -> str: +def select_tautomers( + sdf: str, + k: int | None = None, + window: float | None = None, + *, + overwrite: bool = False, +) -> str: """Select and Write the top-k or E <= window tautomers for each input SMILES Only k or window needs to be specified, NOT both. @@ -24,18 +31,23 @@ def select_tautomers(sdf: str, k: int | None = None, window: float | None = None Output: the path of the low-energy tautomer 3D conformers - Warning: - This function writes ``/_top_tautomers.sdf`` with no - overwrite gate: ``Chem.SDWriter`` truncates on open, so a direct call - such as ``select_tautomers("/data/results.sdf", k=1)`` replaces any - existing ``/data/results_top_tautomers.sdf``. Every route Auto3D - itself takes is safe -- ``get_stable_tautomers`` passes ``main()``'s - output, which lives in a job directory created fresh for that run, and - ``auto3d tautomers`` additionally gates its ``-o`` with - ``check_output_overwrite`` -- so this is a hazard for direct API - callers only, the same residual class as - ``Auto3D.utils.smi_io.smiles2smi`` and ``Auto3D.id_mapping.decode_ids``. See - ``docs/superpowers/follow-ups-after-4.0.0-remediation.md``. + Args: + sdf: Path to the SDF this function reads, normally ``main()``'s output. + k: Keep the top-k tautomers per input molecule. + window: Keep tautomers within this kcal/mol window of the most stable + one. Mutually exclusive with `k`. + overwrite: Allow replacing an existing + ``/_top_tautomers.sdf``. Keyword-only, and + defaults to False because that name is one Auto3D *derives* rather + than one the caller chose: ``Chem.SDWriter`` truncates on open, so + ``select_tautomers("/data/results.sdf", k=1)`` used to replace an + existing ``/data/results_top_tautomers.sdf`` with this call's + selection, silently. Auto3D's own routes are unaffected -- + ``get_stable_tautomers`` passes ``main()``'s output, which lives in + a job directory created fresh for that run, and ``auto3d + tautomers`` additionally gates its ``-o`` with + ``check_output_overwrite`` -- so the default only ever fires for a + direct API caller pointing at an occupied path. Note: Tautomers are ranked by the optimized NNP *electronic* energy (``E_tot``) @@ -48,7 +60,8 @@ def select_tautomers(sdf: str, k: int | None = None, window: float | None = None Raises: ConfigurationError: If both k and window are given, if neither is - given, or if k < 1. ``auto3d tautomers`` already rejects the + given, if k < 1, or if the derived output path exists and + `overwrite` is False. ``auto3d tautomers`` already rejects the both-given case in ``execute_tautomers`` before calling this function, but ``select_tautomers``/``get_stable_tautomers`` are also public Python API entry points that can be called directly, @@ -60,6 +73,16 @@ def select_tautomers(sdf: str, k: int | None = None, window: float | None = None if (k is not None) and (k < 1): raise ConfigurationError(f"tauto_k must be >= 1, got {k}") + # Resolved and gated up front, before the input is read and grouped: + # refusing after all the work is done costs the user the run for nothing, + # and the writer at the bottom truncates the moment it opens. + folder = os.path.dirname(sdf) + # splitext (not split(".")) so an input like 'mol.v2.sdf' keeps 'mol.v2' + # instead of collapsing to 'mol' and risking output collisions. + stem = os.path.splitext(os.path.basename(sdf))[0].strip() + output_path = os.path.join(folder, stem + "_top_tautomers.sdf") + check_output_overwrite(output_path, overwrite) + supplier = Chem.SDMolSupplier(sdf, removeHs=False) mols = [m for m in supplier if m is not None] for mol in mols: @@ -101,14 +124,7 @@ def select_tautomers(sdf: str, k: int | None = None, window: float | None = None else: raise ConfigurationError("Either k OR window needs to be specified") results += out_mols - - folder = os.path.dirname(sdf) - # splitext (not split(".")) so an input like 'mol.v2.sdf' keeps 'mol.v2' - # instead of collapsing to 'mol' and risking output collisions. - stem = os.path.splitext(os.path.basename(sdf))[0].strip() - basename = stem + "_top_tautomers.sdf" - output_path = os.path.join(folder, basename) with Chem.SDWriter(output_path) as w: for mol in results: w.write(mol) diff --git a/src/Auto3D/utils/convergence.py b/src/Auto3D/utils/convergence.py index 68b2505b..b8570922 100644 --- a/src/Auto3D/utils/convergence.py +++ b/src/Auto3D/utils/convergence.py @@ -4,8 +4,7 @@ ``Converged`` is written by :mod:`Auto3D.batch_opt.batchopt` for every record it optimizes -- ``"True"`` or ``"False"`` -- and read by the three filters that decide which conformers survive (:class:`Auto3D.ranking.ConformerRanker`, -:func:`Auto3D.filtering.filter_unique_optimized`, -:func:`Auto3D.filtering.filter_unique`). +:func:`Auto3D.filtering.filter_conformers`). All three used to read it as:: diff --git a/src/Auto3D/utils/energy.py b/src/Auto3D/utils/energy.py index 76fc2787..3ebaecfc 100644 --- a/src/Auto3D/utils/energy.py +++ b/src/Auto3D/utils/energy.py @@ -10,8 +10,7 @@ ``E_tot`` in eV and ``ASE/geometry.opt_geometry`` converted the same tag to Hartree afterwards, so the identical property name carried two units depending on which entry point produced the file -- and the five in-package consumers -(``ranking``, ``filtering.filter_unique_optimized``, -``filtering.filter_unique``) all hard-coded +(``ranking``, ``filtering.filter_conformers``) all hard-coded eV. Feeding an ``opt_geometry`` output straight to ``ConformerRanker(window=2.0)`` therefore opened a window 27.2x too wide, kept 3 conformers where 2 belong, reported ``E_rel`` 0.037 kcal/mol where the truth diff --git a/src/Auto3D/utils/geometry.py b/src/Auto3D/utils/geometry.py index ca8c5018..2f89ae9f 100644 --- a/src/Auto3D/utils/geometry.py +++ b/src/Auto3D/utils/geometry.py @@ -82,7 +82,7 @@ def get_rmsd(mol1: Chem.Mol, mol2: Chem.Mol, remove_hs: bool = True) -> float: The RMSD value in Angstroms. Returns ``float("inf")`` if alignment fails (e.g., due to atom mismatch). An incomparable pair is treated as "distinct" rather than "identical", which is the same convention used - by ``filter_unique``; a downstream ``rmsd < threshold`` check therefore + by the duplicate filter; a downstream ``rmsd < threshold`` check therefore keeps the structure instead of dropping it as a false duplicate. Example: @@ -106,6 +106,6 @@ def get_rmsd(mol1: Chem.Mol, mol2: Chem.Mol, remove_hs: bool = True) -> float: # Temporary bug fix for https://github.com/rdkit/rdkit/issues/6826 rmsd = rdMolAlign.GetBestRMS(mol1_proc, mol2_proc) except RuntimeError: - # Incomparable pair: treat as distinct (inf), matching filter_unique. + # Incomparable pair: treat as distinct (inf), as the filter requires. rmsd = float("inf") return float(rmsd) diff --git a/src/Auto3D/utils/output_guard.py b/src/Auto3D/utils/output_guard.py index 0760452e..4d778ecb 100644 --- a/src/Auto3D/utils/output_guard.py +++ b/src/Auto3D/utils/output_guard.py @@ -8,8 +8,9 @@ ``Auto3D.exceptions``, which is what lets the ``.smi``/``.sdf`` writers under ``utils/`` and the top-level ID/layout helpers gate their own output. -``utils/validation.py`` still re-exports both names, because ``SPE``, -``ASE.thermo`` and ``cli.commands.properties`` import them from there. +This is the only path to these two names. ``utils/validation.py`` re-exported +them for one release while call sites moved over; that re-export is gone, so +importing them from ``validation`` now fails rather than quietly working. """ from __future__ import annotations diff --git a/src/Auto3D/utils/smi_io.py b/src/Auto3D/utils/smi_io.py index 4781705d..38b32523 100644 --- a/src/Auto3D/utils/smi_io.py +++ b/src/Auto3D/utils/smi_io.py @@ -23,6 +23,7 @@ from Auto3D.exceptions import InputValidationError from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.output_guard import check_output_overwrite logger = get_logger(__name__) @@ -86,7 +87,7 @@ def iter_smi_records(path, *, on_malformed="skip"): yield line_no, parts[0], parts[1] -def smiles2smi(smiles: list[str], path: str) -> str: +def smiles2smi(smiles: list[str], path: str, *, overwrite: bool = True) -> str: """Convert a list of SMILES strings to a .smi file with InChIKey IDs. Each SMILES string is converted to a molecule, and its InChIKey is computed @@ -96,10 +97,22 @@ def smiles2smi(smiles: list[str], path: str) -> str: Args: smiles: List of SMILES strings to convert. path: Output file path for the .smi file. + overwrite: Allow replacing an existing file at `path`. Keyword-only. + Defaults to **True**, unlike the writers that derive their own + output name (``id_mapping.encode_ids``/``decode_ids``, + ``tautomer.select_tautomers``): the caller named this file, so + refusing it would gate a decision the caller already made. + ``smiles2mols`` writes into a fresh ``TemporaryDirectory`` on every + call and relies on that. Pass False to make the write refuse an + occupied path. Returns: The output file path. + Raises: + ConfigurationError: `path` exists and `overwrite` is False. + InputValidationError: A SMILES string could not be parsed by RDKit. + Example: >>> smiles2smi(["CCO", "CCC"], "molecules.smi") 'molecules.smi' @@ -107,6 +120,11 @@ def smiles2smi(smiles: list[str], path: str) -> str: # CCO LFQSCWFLJHTTHZ-UHFFFAOYSA-N # CCC ATUOYWHBWRKTHZ-UHFFFAOYSA-N """ + # Checked before any SMILES is parsed: refusing after the work is done is a + # worse experience for no gain, and `open(path, "w+")` below truncates on + # open, so by the time the write starts the file is already gone. + check_output_overwrite(path, overwrite) + lines = [] seen_ids: dict[str, int] = {} for idx, smi in enumerate(smiles): diff --git a/src/Auto3D/utils/stereochemistry.py b/src/Auto3D/utils/stereochemistry.py index f6d7e2af..14c2ba51 100644 --- a/src/Auto3D/utils/stereochemistry.py +++ b/src/Auto3D/utils/stereochemistry.py @@ -18,6 +18,7 @@ from Auto3D.utils.atomic_io import atomic_write_path from Auto3D.utils.logging_config import get_logger +from Auto3D.utils.smi_io import iter_smi_records logger = get_logger(__name__) @@ -252,12 +253,16 @@ def remove_enantiomers(inpath: str, out: str) -> dict[str, list[str]]: If enantiomer removal fails for a molecule, all original SMILES are kept and a message is printed. """ - with open(inpath) as f: - data = f.readlines() - + # iter_smi_records is the single parser for this format (M59). The loop + # this replaced built `vals = line.split()` and indexed vals[0]/vals[1] + # with no guard, so a blank line in `inpath` raised a bare IndexError; + # "raise" keeps that fail-fast behavior (as an InputValidationError naming + # the line, not a crash) for any genuinely malformed row while blank lines + # -- which should never appear in a file this module itself wrote, but + # previously crashed the whole function if one did -- are now skipped like + # every other consumer of this format. smiles: dict[str, list[str]] = defaultdict(lambda: []) - for line in data: - vals = line.split() + for _line_no, smi, mol_id in iter_smi_records(inpath, on_malformed="raise"): # Strip only the trailing isomer-index component write_enumerated_smi # appends (rsplit, maxsplit=1), not everything after the first # underscore: an id like "KEY_2" -- smiles2smi's disambiguation of a @@ -265,7 +270,7 @@ def remove_enantiomers(inpath: str, out: str) -> dict[str, list[str]]: # specifically so it is not dropped -- must survive this grouping # intact, or it silently merges back onto "KEY" here before ranking # ever sees it (M17). - smi, name = vals[0].strip(), vals[1].strip().rsplit("_", 1)[0].strip() + name = mol_id.rsplit("_", 1)[0] smiles[name].append(smi) for key, values in smiles.items(): @@ -407,47 +412,35 @@ def create_enantiomer(smi: str) -> str: 'C[C@@H](O)F' """ stereo_info = get_stereo_info(smi) - new_smi = "" keys = list(stereo_info.keys()) - if len(keys) == 0: + if not keys: # No stereo centers to invert return smi - if len(keys) == 1: - key = keys[0] + + # Single pass with a cursor, one key at a time: emit smi[cursor:key] + # verbatim, then the inverted marker, then advance the cursor past the + # marker just consumed. The tail after the last key is emitted once, + # after the loop, from the final cursor position. + # + # This replaces two copies of the same logic that used to disagree in + # shape: a `len(keys) == 1` branch handled one stereo center, and a + # general loop handled two or more by reading `key2`/`val2` (set inside + # the loop's `else` branch) *after* the `for` -- correct only because + # Python does not scope loop variables to the loop body, so `key2` still + # held the last iteration's value. That accident breaks the moment + # `keys` can be empty going into the tail read, which is exactly what the + # `len(keys) == 1` branch existed to avoid (M60). + inverted = {"@": "@@", "@@": "@"} + new_smi = "" + cursor = 0 + for key in keys: val = stereo_info[key] - if val == "@": - new_smi += smi[:key] - new_smi += "@@" - new_smi += smi[(key + 1) :] - elif val == "@@": - new_smi += smi[:key] - new_smi += "@" - new_smi += smi[(key + 2) :] - else: + if val not in inverted: raise ValueError("Invalid %s" % smi) - return new_smi - - for i in range(len(keys)): - if i == 0: - key = keys[i] - new_smi += smi[:key] - else: - key1 = keys[i - 1] - key2 = keys[i] - val1 = stereo_info[key1] - if val1 == "@": - new_smi += "@@" - new_smi += smi[int(key1 + 1) : key2] - elif val1 == "@@": - new_smi += "@" - new_smi += smi[int(key1 + 2) : key2] - val2 = stereo_info[key2] - if val2 == "@": - new_smi += "@@" - new_smi += smi[int(key2 + 1) :] - elif val2 == "@@": - new_smi += "@" - new_smi += smi[int(key2 + 2) :] + new_smi += smi[cursor:key] + new_smi += inverted[val] + cursor = key + len(val) + new_smi += smi[cursor:] return new_smi @@ -499,16 +492,20 @@ def amend_configuration(smis: str) -> dict[str, list[str]]: For input "N=C1OC(CN2CC(C)OC(C)C2)CN1", if some stereo configurations are missing, this function will attempt to add them. """ - with open(smis) as f: - data = f.readlines() + # iter_smi_records is the single parser for this format (M59). The + # `tuple(line.strip().split())` this replaced required EXACTLY 2 tokens + # (a third raised "too many values to unpack") and had no blank/comment + # handling at all; "raise" preserves the fail-fast behavior for a + # genuinely malformed row (now a named InputValidationError instead of a + # bare ValueError) while tolerating a trailing extra column, matching + # every other consumer of this format. dct: dict[str, list[str]] = defaultdict(lambda: []) - for line in data: - smi, idx = tuple(line.strip().split()) + for _line_no, smi, mol_id in iter_smi_records(smis, on_malformed="raise"): # See the matching note in remove_enantiomers: strip only the # trailing isomer-index component, not everything after the first # underscore, so a disambiguated id like "KEY_2" is not merged back # onto "KEY" here (M17). - idx = idx.rsplit("_", 1)[0].strip() + idx = mol_id.rsplit("_", 1)[0] dct[idx].append(smi) for key in dct.keys(): diff --git a/src/Auto3D/utils/validation.py b/src/Auto3D/utils/validation.py index af2656a9..deed6858 100644 --- a/src/Auto3D/utils/validation.py +++ b/src/Auto3D/utils/validation.py @@ -22,16 +22,7 @@ ModelLoadError, ) from Auto3D.utils.logging_config import get_logger - -# The two output guards moved to the leaf module ``utils/output_guard.py`` so -# that a writer needing only them does not import this module's torch. They are -# re-exported here because ``SPE``, ``ASE.thermo`` and -# ``cli.commands.properties`` import them from this path; new call sites should -# name ``Auto3D.utils.output_guard`` directly. -from Auto3D.utils.output_guard import ( # noqa: F401 - check_output_not_input, - check_output_overwrite, -) +from Auto3D.utils.smi_io import iter_smi_records from Auto3D.utils.stereochemistry import count_unspecified_stereo if TYPE_CHECKING: @@ -291,31 +282,15 @@ def check_smi_format(args: Any) -> tuple[bool, list[str]]: """ ANI = True - smiles_all = [] - with open(args.path) as f: - data = f.readlines() - for line in data: - if line.isspace(): - continue - # Skip '#'-prefixed comment lines, matching cli.commands.validate. - # validate_smiles_file and file_ops.iter_smi_records -- all three must - # agree on what counts as a comment vs. data, or `auto3d validate` - # would approve a file this function then rejects (M25). - if line.lstrip().startswith("#"): - continue - # Tolerate ragged rows the way the rest of the pipeline does: the chunk - # loader reads only the first two whitespace columns (usecols=[0, 1]), so - # trailing tokens (e.g. an inline comment column) must not be rejected - # here. split() never yields empty tokens, so a present SMILES/ID is - # guaranteed non-empty. - parts = line.split() - if len(parts) < 2: - raise InputValidationError( - "Each non-blank line must contain a SMILES and an ID separated by " - f"whitespace, but got: {line.strip()!r}" - ) - smiles = parts[0] # parts[1] is the ID; its presence is enforced above - smiles_all.append(smiles) + # iter_smi_records is the single parser for this format (M59): it already + # skips blank lines and '#'-prefixed comments the same way + # cli.commands.validate.validate_smiles_file does, and its on_malformed + # ("raise") gives the same InputValidationError this loop used to raise by + # hand for a line missing an ID, so `auto3d validate` and this check + # cannot silently disagree about what a well-formed line looks like + # (M25). The parser also tolerates ragged rows (extra whitespace columns + # beyond SMILES+ID), matching the chunk loader's usecols=[0, 1]. + smiles_all = [smiles for _line_no, smiles, _id in iter_smi_records(args.path, on_malformed="raise")] logger.info(f"\tThere are {len(smiles_all)} SMILES in the input file {args.path}.") logger.info("\tAll SMILES and IDs are valid.") diff --git a/tests/test_chunk_manager.py b/tests/test_chunk_manager.py index e5a0c960..fe83549f 100644 --- a/tests/test_chunk_manager.py +++ b/tests/test_chunk_manager.py @@ -395,6 +395,147 @@ def test_chunk_size_clamped_to_at_least_one(self, tmp_path): assert 1 <= len(chunk_info) <= 4 +class TestSmiParserCrossAgreement: + """M59: chunk_manager's pandas-based `.smi` reader stays separate from + ``iter_smi_records`` for performance (avoiding a pure-Python per-line loop + over what can be a very large file), but the two must not silently drift + apart on the input they both exist to read: an already-ID-encoded + intermediate `.smi` file with no blank lines and no comments (chunk_manager + calls this "encode_ids semantics" in prepare_chunks' docstring). If either + parser's tokenizing/whitespace-splitting rule changes, this test should + catch the divergence. + """ + + def test_prepare_chunks_agrees_with_iter_smi_records_on_well_formed_input( + self, tmp_path + ): + """Both parsers must extract the same (smiles, id) pairs, in the same + order, from a well-formed encoded .smi file.""" + from Auto3D.utils.smi_io import iter_smi_records + + rows = [ + ("CCO", "0"), + ("CCCO", "1"), + ("c1ccccc1", "2"), + ("C[C@H](O)F", "3"), + ("CCN", "4"), + ] + input_file = tmp_path / "test_encoded.smi" + input_file.write_text("".join(f"{smi} {mol_id}\n" for smi, mol_id in rows)) + + # memory * capacity is chosen well above len(rows) so chunk_size never + # forces a split: exactly one chunk, holding every row in original + # order, so the chunk file's content is directly comparable to + # iter_smi_records' output order (round-robin distribution across + # multiple chunks would otherwise reorder rows relative to the input). + config = Auto3DOptions( + path=str(tmp_path / "test.smi"), k=1, memory=1, capacity=1000 + ) + manager = ChunkManager( + config=config, + input_path=input_file, + input_format="smi", + job_dir=tmp_path, + workflow_logger=None, + ) + + chunk_info = manager.prepare_chunks() + assert len(chunk_info) == 1, "test assumes no chunk splitting" + + chunk_path, _ = chunk_info[0] + pandas_rows = [ + tuple(line.split()) + for line in Path(chunk_path).read_text().splitlines() + if line.strip() + ] + iter_rows = [ + (smi, mol_id) for _line_no, smi, mol_id in iter_smi_records(str(input_file)) + ] + + assert pandas_rows == rows + assert iter_rows == rows + assert pandas_rows == iter_rows + + def test_ragged_extra_column_agreement(self, tmp_path): + """A trailing whitespace-separated column beyond SMILES+ID must be + dropped identically by both readers (chunk_manager's usecols=[0, 1] + vs. iter_smi_records taking only parts[0]/parts[1]).""" + from Auto3D.utils.smi_io import iter_smi_records + + input_file = tmp_path / "test_encoded.smi" + input_file.write_text("CCO 0 inline_comment_column\nCCCO 1\n") + + config = Auto3DOptions( + path=str(tmp_path / "test.smi"), k=1, memory=1, capacity=1000 + ) + manager = ChunkManager( + config=config, + input_path=input_file, + input_format="smi", + job_dir=tmp_path, + workflow_logger=None, + ) + chunk_info = manager.prepare_chunks() + assert len(chunk_info) == 1 + + chunk_path, _ = chunk_info[0] + pandas_rows = [ + tuple(line.split()) + for line in Path(chunk_path).read_text().splitlines() + if line.strip() + ] + iter_rows = [ + (smi, mol_id) for _line_no, smi, mol_id in iter_smi_records(str(input_file)) + ] + + assert pandas_rows == iter_rows == [("CCO", "0"), ("CCCO", "1")] + + def test_documented_divergence_comment_lines(self, tmp_path): + """Documents a real, deliberate divergence rather than papering over + it: iter_smi_records treats a '#'-prefixed line as a comment and + skips it (matching cli.commands.validate, per M25); chunk_manager's + pd.read_csv has no `comment=` parameter and reads it as a data row. + + This is not a bug this task fixes (chunk_manager stays a separate, + faster reader by design) -- it is recorded here so that a future + change to either parser's comment handling is a deliberate, + visible decision rather than a silent behavior change caught only in + production. If this test starts failing because someone taught + pd.read_csv to skip '#' lines, that is progress: update the + assertion, don't just delete the test. + """ + from Auto3D.utils.smi_io import iter_smi_records + + input_file = tmp_path / "test_encoded.smi" + input_file.write_text("CCO 0\n# 1 2\nCCCO 3\n") + + config = Auto3DOptions( + path=str(tmp_path / "test.smi"), k=1, memory=1, capacity=1000 + ) + manager = ChunkManager( + config=config, + input_path=input_file, + input_format="smi", + job_dir=tmp_path, + workflow_logger=None, + ) + chunk_info = manager.prepare_chunks() + chunk_path, _ = chunk_info[0] + pandas_rows = [ + tuple(line.split()) + for line in Path(chunk_path).read_text().splitlines() + if line.strip() + ] + iter_rows = [ + (smi, mol_id) for _line_no, smi, mol_id in iter_smi_records(str(input_file)) + ] + + # pandas reads the '#' line as data; iter_smi_records skips it. + assert pandas_rows == [("CCO", "0"), ("#", "1"), ("CCCO", "3")] + assert iter_rows == [("CCO", "0"), ("CCCO", "3")] + assert pandas_rows != iter_rows + + class TestLogging: """Tests for ChunkManager logging.""" diff --git a/tests/test_cli_errors.py b/tests/test_cli_errors.py index 84b55324..d67a6a0c 100644 --- a/tests/test_cli_errors.py +++ b/tests/test_cli_errors.py @@ -288,7 +288,7 @@ def _flat(text: str) -> str: def test_overwrite_refusal_does_not_suggest_config_init(capsys, tmp_path): """Raised by the real guard, not hand-built: this pins the raise site's choice of hint together with handle_error's presentation of it.""" - from Auto3D.utils.validation import check_output_overwrite + from Auto3D.utils.output_guard import check_output_overwrite existing = tmp_path / "precious.sdf" existing.write_text("x") diff --git a/tests/test_cli_property_commands.py b/tests/test_cli_property_commands.py index 66a58d22..0ed86409 100644 --- a/tests/test_cli_property_commands.py +++ b/tests/test_cli_property_commands.py @@ -481,7 +481,7 @@ def test_tautomers_refuses_output_equal_to_input(smi): # `auto3d energy junk.sdf --no-gpu -o precious.sdf` used to exit 0, print # "Wrote precious.sdf", and leave precious.sdf at 0 bytes. `config init` has # had -f/--force since it shipped; these four commands did not. The guard -# itself (`Auto3D.utils.validation.check_output_overwrite`) is exercised per +# itself (`Auto3D.utils.output_guard.check_output_overwrite`) is exercised per # API function in tests/test_durability.py; what is pinned here is the CLI # half -- that each command actually *passes* its flag down, which is the part # a refactor drops silently. diff --git a/tests/test_durability.py b/tests/test_durability.py index a8a3c818..ed8f3c39 100644 --- a/tests/test_durability.py +++ b/tests/test_durability.py @@ -510,7 +510,7 @@ def test_guard_compares_real_paths_not_strings(self, job_dir, monkeypatch): a guard that always raised would satisfy every test above. """ from Auto3D.exceptions import ConfigurationError - from Auto3D.utils.validation import check_output_not_input + from Auto3D.utils.output_guard import check_output_not_input sdf = job_dir / "mols.sdf" _write_sdf(sdf, ["m1"]) @@ -555,7 +555,7 @@ def test_a_hardlink_to_the_input_is_refused(self, job_dir): directly on this ext4 box, but it takes the identical samefile path. """ from Auto3D.exceptions import ConfigurationError - from Auto3D.utils.validation import check_output_not_input + from Auto3D.utils.output_guard import check_output_not_input sdf = job_dir / "mols.sdf" _write_sdf(sdf, ["a"]) @@ -572,7 +572,7 @@ def test_a_hardlink_to_the_input_is_refused(self, job_dir): def test_a_genuinely_different_output_is_still_allowed(self, job_dir): """Negative control: without this, a guard that always raised would satisfy every other test in this class.""" - from Auto3D.utils.validation import check_output_not_input + from Auto3D.utils.output_guard import check_output_not_input sdf = job_dir / "mols.sdf" _write_sdf(sdf, ["a"]) @@ -976,7 +976,7 @@ def test_conformer_ranker_refuses_an_existing_output(self, job_dir): def test_the_guard_permits_everything_it_should(self, job_dir): """Negative controls for the shared function itself.""" - from Auto3D.utils.validation import check_output_overwrite + from Auto3D.utils.output_guard import check_output_overwrite existing = job_dir / "there.sdf" existing.write_text("x") @@ -992,7 +992,7 @@ def test_the_guard_accepts_str_and_path_alike(self, job_dir): an `os.path.join` string. A guard that only handled one would silently no-op for the other.""" from Auto3D.exceptions import ConfigurationError - from Auto3D.utils.validation import check_output_overwrite + from Auto3D.utils.output_guard import check_output_overwrite existing = job_dir / "there.sdf" existing.write_text("x") diff --git a/tests/test_filter_unique.py b/tests/test_filter_unique.py deleted file mode 100644 index e3ce2d35..00000000 --- a/tests/test_filter_unique.py +++ /dev/null @@ -1,268 +0,0 @@ -#!/usr/bin/env python -"""Tests for the legacy all-pairs ``Auto3D.filtering.filter_unique``. - -Kept as its own file because ``filter_unique`` is scheduled for removal once -``filter_unique_optimized`` is the single filter; deleting it then means -deleting this file, not surgery on a shared one. -""" -from __future__ import annotations - -import pytest # noqa: F401 -from rdkit import Chem -from rdkit.Chem import AllChem - -from Auto3D.filtering import filter_unique - - -class TestFilterUnique: - """Test the filter_unique function for RMSD-based duplicate filtering.""" - - def test_filter_identical_conformers(self): - """Test that identical conformers are filtered to one.""" - - mol = Chem.MolFromSmiles("CCO") - mol = Chem.AddHs(mol) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") - - # Create identical copies - mol2 = Chem.Mol(mol) - mol2.SetProp("Converged", "true") - - mols = [mol, mol2] - unique_mols = filter_unique(mols, crit=0.3) - - # Should only keep one - assert len(unique_mols) == 1 - - def test_same_geometry_different_energy_kept(self): - """Identical geometry but distinct E_tot must be kept (energy guard). - - Heavy-atom RMSD ~= 0 but the two are distinct minima (the O-H rotamer - case); the energy guard must stop them collapsing into one. - """ - - mol = Chem.MolFromSmiles("CCO") - mol = Chem.AddHs(mol) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") - mol.SetProp("E_tot", "-10.0") - - mol2 = Chem.Mol(mol) # identical geometry - mol2.SetProp("Converged", "true") - mol2.SetProp("E_tot", "-10.5") # |dE| >> tol - - unique_mols = filter_unique([mol, mol2], crit=0.3) - assert len(unique_mols) == 2 - - def test_missing_energy_falls_back_to_rmsd_only(self): - """Without E_tot the energy guard cannot apply -> RMSD-only dedup. - - Preserves the legacy behavior for callers that do not set E_tot. - """ - - mol = Chem.MolFromSmiles("CCO") - mol = Chem.AddHs(mol) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") # no E_tot set - - mol2 = Chem.Mol(mol) - mol2.SetProp("Converged", "true") - - unique_mols = filter_unique([mol, mol2], crit=0.3) - assert len(unique_mols) == 1 - - def test_filter_different_conformers(self): - """Test that different conformers are kept.""" - - mol1 = Chem.MolFromSmiles("CCCCCC") # Hexane - flexible - mol1 = Chem.AddHs(mol1) - AllChem.EmbedMolecule(mol1, randomSeed=42) - mol1.SetProp("Converged", "true") - - mol2 = Chem.MolFromSmiles("CCCCCC") - mol2 = Chem.AddHs(mol2) - AllChem.EmbedMolecule(mol2, randomSeed=123) - mol2.SetProp("Converged", "true") - - # Generate very different conformers by using different seeds - # and moving atoms around - conf = mol2.GetConformer() - pos = conf.GetAtomPosition(0) - conf.SetAtomPosition(0, (pos.x + 0.5, pos.y, pos.z)) - - mols = [mol1, mol2] - unique_mols = filter_unique(mols, crit=0.3) - - # Should keep both (or at least not crash) - assert len(unique_mols) >= 1 - - def test_two_diastereomers_are_never_merged(self): - """A distinct compound must survive dedup, however close its geometry. - - The same guarantee ``tests/test_filtering.py`` asserts for - ``filter_unique_optimized``, asserted here because this is the other - duplicate filter and it applies the identical RMSD-plus-energy criterion. - Fixing one path and not the other would leave the defect reachable - through ``ConformerRanker(use_optimized_filtering=False)``. - - cis/trans-4-tert-butylcyclohexanol: heavy-atom RMSD between the two - diastereomers was measured at 0.300 A, i.e. at the 0.3 A default - threshold. ``crit`` is opened wide here so RMSD and energy both say - "duplicate" and the only thing that can keep the pair apart is the fact - that they are different compounds. - """ - from Auto3D.utils.energy import set_e_tot_from_ev - - def build(smiles: str) -> Chem.Mol: - mol = Chem.AddHs(Chem.MolFromSmiles(smiles)) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") - set_e_tot_from_ev(mol, -10.0) # identical energies - return mol - - cis = build("O[C@H]1CC[C@@H](CC1)C(C)(C)C") - trans = build("O[C@H]1CC[C@H](CC1)C(C)(C)C") - assert Chem.MolToSmiles(cis) != Chem.MolToSmiles(trans), "test premise" - - assert len(filter_unique([cis, trans], crit=10.0)) == 2, ( - "the legacy filter merged two distinct diastereomers, so an input " - "molecule vanished from the output with no record" - ) - - def test_duplicate_conformers_of_one_stereoisomer_still_collapse(self): - """The other half: the stereo guard must narrow dedup, not disable it.""" - from Auto3D.utils.energy import set_e_tot_from_ev - - def build() -> Chem.Mol: - mol = Chem.AddHs(Chem.MolFromSmiles("O[C@H]1CC[C@@H](CC1)C(C)(C)C")) - AllChem.EmbedMolecule(mol, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol) - mol.SetProp("Converged", "true") - set_e_tot_from_ev(mol, -10.0) - return mol - - assert len(filter_unique([build(), build()], crit=10.0)) == 1, ( - "duplicate conformers of one stereoisomer survived, so the stereo " - "guard has switched dedup off rather than narrowing it" - ) - - def test_filter_unconverged_removed(self): - """Test that unconverged structures are removed.""" - - mol1 = Chem.MolFromSmiles("CCO") - mol1 = Chem.AddHs(mol1) - AllChem.EmbedMolecule(mol1, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol1) - mol1.SetProp("Converged", "true") - - mol2 = Chem.MolFromSmiles("CCO") - mol2 = Chem.AddHs(mol2) - AllChem.EmbedMolecule(mol2, randomSeed=123) - mol2.SetProp("Converged", "false") # Not converged - - mols = [mol1, mol2] - unique_mols = filter_unique(mols, crit=0.3) - - # Only converged one should remain - assert len(unique_mols) == 1 - assert unique_mols[0].GetProp("Converged").lower() == "true" - - def test_filter_empty_list(self): - """Test filtering empty list returns empty list.""" - - unique_mols = filter_unique([], crit=0.3) - assert len(unique_mols) == 0 - - def test_filter_custom_threshold(self): - """Test that custom RMSD threshold works.""" - - mol1 = Chem.MolFromSmiles("CCO") - mol1 = Chem.AddHs(mol1) - AllChem.EmbedMolecule(mol1, randomSeed=42) - AllChem.MMFFOptimizeMolecule(mol1) - mol1.SetProp("Converged", "true") - - mol2 = Chem.Mol(mol1) - mol2.SetProp("Converged", "true") - - mols = [mol1, mol2] - - # With very small threshold, might keep both - unique_mols_small = filter_unique(mols, crit=0.0001) - # With large threshold, definitely keep only one - unique_mols_large = filter_unique(mols, crit=10.0) - - # Large threshold should definitely merge identical mols - assert len(unique_mols_large) == 1 - # A tighter threshold can never merge MORE than a looser one -- the - # discarded half of this test's own computation, now actually checked. - assert len(unique_mols_small) >= len(unique_mols_large) - - def test_filter_unique_removehs_is_linear_and_nondestructive(self, monkeypatch): - """Legacy filter_unique strips Hs once per molecule (not per comparison) and - returns the originals with explicit H + exact positions intact.""" - import numpy as np - from rdkit import Chem - from rdkit.Chem import AllChem - - from Auto3D import filtering - - base = Chem.AddHs(Chem.MolFromSmiles("CCCCO")) - cids = AllChem.EmbedMultipleConfs(base, numConfs=5, randomSeed=1) - mols = [] - for cid in cids: - m = Chem.Mol(base, confId=int(cid)) - m.SetProp("Converged", "true") - mols.append(m) - n_atoms = base.GetNumAtoms() - orig_pos = {id(m): m.GetConformer().GetPositions().copy() for m in mols} - - calls = {"n": 0} - real_removehs = filtering.Chem.RemoveHs - - def counting(mol, *a, **k): - calls["n"] += 1 - return real_removehs(mol, *a, **k) - - monkeypatch.setattr(filtering.Chem, "RemoveHs", counting) - - result = filtering.filter_unique(mols, crit=0.01) - assert calls["n"] == len(mols) # once per input, never per pair - assert len(result) == len(mols) - for m in result: - assert m.GetNumAtoms() == n_atoms - assert any(a.GetAtomicNum() == 1 for a in m.GetAtoms()) - assert np.array_equal(m.GetConformer().GetPositions(), orig_pos[id(m)]) - - def test_rmsd_failure_keeps_both(self, monkeypatch): - """An incomparable pair (RMSD raises) must NOT be treated as a duplicate. - - When GetBestRMS raises RuntimeError, filter_unique must treat the pair - as distinct (rmsd = inf) and keep both, mirroring the fix already in - filtering._filter_within_cluster. The previous behavior (rmsd = 0) - made distinct conformers look like perfect duplicates and dropped one. - """ - from Auto3D import filtering - - def make(name): - m = Chem.AddHs(Chem.MolFromSmiles("CCO")) - AllChem.EmbedMolecule(m, randomSeed=abs(hash(name)) % 1000) - AllChem.MMFFOptimizeMolecule(m) - m.SetProp("_Name", name) - m.SetProp("Converged", "true") - return m - - def boom(*args, **kwargs): - raise RuntimeError("GetBestRMS failed") - - # filter_unique calls rdMolAlign.GetBestRMS via the filtering module. - monkeypatch.setattr(filtering.rdMolAlign, "GetBestRMS", boom) - - mols = [make("a"), make("b")] - unique_mols = filtering.filter_unique(mols, crit=0.3) - assert len(unique_mols) == 2 # incomparable pair must NOT be dropped diff --git a/tests/test_filtering.py b/tests/test_filtering.py index 6974ca7c..759858d0 100644 --- a/tests/test_filtering.py +++ b/tests/test_filtering.py @@ -1,12 +1,20 @@ #!/usr/bin/env python -"""Tests for optimized RMSD filtering with energy clustering.""" +"""Tests for the single conformer filter (RMSD dedup with energy clustering).""" from __future__ import annotations +import os + import pytest from rdkit import Chem from rdkit.Chem import AllChem -from Auto3D.filtering import filter_unique_optimized, _filter_within_cluster +from Auto3D.filtering import ( + DROP_REASONS, + FilterResult, + _filter_within_cluster, + filter_conformers, + filter_unique_optimized, +) from Auto3D.utils.energy import set_e_tot_from_ev @@ -331,38 +339,43 @@ def test_small_energy_window_creates_separate_clusters(self): assert len(result) == 2 +def _energyless(smiles: str, seed: int = 42, name: str = "") -> Chem.Mol: + """An embedded, converged conformer carrying NO 'E_tot' property. + + This is what an SDF Auto3D's optimizer did not write can look like: a + hand-built conformer set, or an export that names its energy field + something else. ``ConformerRanker`` refuses such a record up front + (``InputValidationError``), so this shape only reaches the filter through + a direct API call -- which is exactly the caller the two filters used to + disagree for. + """ + mol = Chem.AddHs(Chem.MolFromSmiles(smiles)) + AllChem.EmbedMolecule(mol, randomSeed=seed) + AllChem.MMFFOptimizeMolecule(mol) + mol.SetProp("Converged", "true") + if name: + mol.SetProp("_Name", name) + assert not mol.HasProp("E_tot"), "helper premise: no energy property" + return mol + + class TestMissingEnergyPropertyMustNotCrash: - """filter_unique_optimized must tolerate a record with no 'E_tot', the - way the legacy ``filtering.filter_unique`` already does. - - KNOWN DEFECT (found during cluster E brainstorming, not fixed by this - lane): ``filtering.py:75`` sorts the valid-mols list by - ``Auto3D.utils.energy.e_tot_ev``, which RAISES (KeyError/ValueError) for - a molecule with no usable 'E_tot' property. ``_filter_within_cluster``'s - own energy guard, two dozen lines later in the same module, instead uses - the tolerant ``try_e_tot_ev`` and treats a missing energy as "fall back - to RMSD only". ``filtering.filter_unique`` (the OTHER conformer - filter, sharing the same duplicate criterion since 4.0.1) also uses - ``try_e_tot_ev`` throughout and does not crash on this input. So the two - filters diverge on malformed input: the same list of mols that - ``filter_unique`` happily filters crashes ``filter_unique_optimized``. - - This matters here specifically because cluster B5 is about to delete one - of the two filters, and the survivor is the stricter (crashing) one -- - fixing filtering.py's sort key to use ``try_e_tot_ev``, matching its own - energy guard and the legacy filter, is what should make this pass. + """The one conformer filter must tolerate a record with no 'E_tot'. + + ``filtering.py`` used to sort the valid-mols list by + ``Auto3D.utils.energy.e_tot_ev``, which RAISES (KeyError/ValueError) for a + molecule with no usable 'E_tot'. ``_filter_within_cluster``'s own energy + guard, two dozen lines later in the same module, instead used the tolerant + ``try_e_tot_ev`` and treated a missing energy as "fall back to RMSD only", + as did the legacy all-pairs ``filter_unique`` throughout. So the same list + of mols one filter happily filtered crashed the other -- and the survivor + of the two was the crashing one. + + Tolerating the record must not mean *inventing* an energy for it: a + missing 'E_tot' read as 0.0 would sort a garbage record to the front of a + list of negative energies and hand it to ``top_k`` as the global minimum. """ - @pytest.mark.xfail( - strict=True, - reason=( - "filtering.py:75 sorts by e_tot_ev (raises KeyError for a mol " - "with no 'E_tot' property) instead of the tolerant try_e_tot_ev " - "that _filter_within_cluster's own energy guard and the legacy " - "filtering.filter_unique both use -- the two conformer " - "filters disagree on malformed input (cluster E brainstorm defect)." - ), - ) def test_missing_e_tot_property_does_not_crash(self): mol_no_energy = Chem.AddHs(Chem.MolFromSmiles("CCO")) AllChem.EmbedMolecule(mol_no_energy, randomSeed=42) @@ -380,20 +393,116 @@ def test_missing_e_tot_property_does_not_crash(self): # cannot be deduped by energy -- it must survive alongside the other. assert len(result) == 2 + def test_a_garbled_e_tot_is_tolerated_the_same_way(self): + """``try_e_tot_ev`` swallows ValueError too, so a non-numeric property + must take the same path as an absent one rather than crash.""" + mol_garbled = _energyless("CCO", seed=42) + mol_garbled.SetProp("E_tot", "not-a-number") + + mol_with_energy = _create_mol_with_energy("CC", -10.0) + + result = filter_unique_optimized( + [mol_garbled, mol_with_energy], rmsd_threshold=0.3 + ) + assert len(result) == 2 + + def test_an_energyless_record_sorts_last_whatever_the_real_energies_are(self): + """A tolerated missing energy must not be *compared* as an energy. + + The tolerant path needs some placeholder to sort on; the danger is that + the placeholder takes part in the comparison. A record with no ``E_tot`` + that sorts as if its energy were 0.0 lands ahead of every genuine + structure whose energy is above 0.0 -- and the first element of the + filter's output is what ``top_k``/``top_window`` treat as the global + minimum: the reference ``E_rel`` is measured from, and the single + structure a ``k=1`` request returns. + + ``E_tot`` is a user-supplied SD property, so "every real energy is a + large negative number" is a convention, not a guarantee -- hence the + deliberately positive energy below. The record with no energy must sort + last regardless of where the real energies fall relative to the + placeholder. + """ + low = _create_mol_with_energy("CCO", -100.0) + high = _create_mol_with_energy("CCCO", -50.0) + positive = _create_mol_with_energy("CCCCCCO", +25.0) + unknown = _energyless("CCCCO", seed=7) + + result = filter_unique_optimized( + [unknown, positive, high, low], rmsd_threshold=0.3 + ) + + assert len(result) == 4 + assert [Chem.MolToSmiles(Chem.RemoveHs(m)) for m in result] == [ + Chem.MolToSmiles(Chem.RemoveHs(low)), + Chem.MolToSmiles(Chem.RemoveHs(high)), + Chem.MolToSmiles(Chem.RemoveHs(positive)), + Chem.MolToSmiles(Chem.RemoveHs(unknown)), + ], "the record with no energy must sort after every record that has one" + + def test_two_energyless_duplicates_still_collapse(self): + """The inverse assertion, and the reason the two above are safe. + + A filter that declared every energy-less record distinct would pass + both tests above while silently switching duplicate removal off for + this whole class of input. Two identical geometries with no energy at + all must still collapse to one: the energy guard cannot apply, so RMSD + alone decides -- which is precisely what the legacy all-pairs filter + did. + """ + first = _energyless("CCO", seed=42) + second = _energyless("CCO", seed=42) + + result = filter_unique_optimized([first, second], rmsd_threshold=0.3) -class TestFilterUniqueBehavior: - """Tests verifying behavior matches original filter_unique.""" + assert len(result) == 1, ( + "two bit-identical energy-less conformers both survived, so " + "tolerating a missing energy has disabled dedup for this input " + "rather than falling back to RMSD only" + ) - def test_matches_original_for_simple_case(self): - """Should produce same results as original filter_unique for basic cases. + def test_an_energyless_record_is_compared_against_one_that_has_energy(self): + """Energy-less records are compared against EVERYTHING, not just each + other. - Both filters operate on conformers of the *same* molecule (that is the - real Auto3D contract). Use genuinely distinct conformers of one molecule - so RMSD is well-defined and comparable across both implementations. + No energy gap can prove a pair is not a duplicate when one side has no + energy, so the partitioning that makes the filter sub-quadratic has no + licence to separate such a record from anything. The legacy all-pairs + filter compared it to every survivor; so must this one. """ - from Auto3D.filtering import filter_unique + with_energy = _create_mol_with_energy("CCO", -100.0) + without = _energyless("CCO", seed=42) + # Same compound, same embedding seed and force field -> RMSD ~= 0. + assert Chem.MolToSmiles(Chem.RemoveHs(with_energy)) == Chem.MolToSmiles( + Chem.RemoveHs(without) + ), "test premise: same compound" - def conformer(seed: float, energy_ev: float) -> Chem.Mol: + result = filter_unique_optimized( + [with_energy, without], rmsd_threshold=0.3, energy_cluster_window=0.01 + ) + + assert len(result) == 1, ( + "a record with no energy escaped comparison against a duplicate " + "that has one, because the energy partitioning separated them" + ) + + +class TestTheSurvivingFilterKeepsTheLegacyVerdicts: + """Values recorded from the legacy all-pairs ``filter_unique`` before it was + deleted (cluster B5 phase 4a). + + Auto3D carried two conformer filters with the same duplicate criterion, each + acting as the other's oracle, until 4.1.0. These cases were run against BOTH + and are asserted here as literals so the surviving filter's verdicts stay + pinned now that there is nothing left to compare against. + """ + + def test_distinct_conformers_of_one_molecule_all_survive(self): + """Three genuinely different conformers of a flexible chain: 3 kept. + + Both filters returned 3 for this input. + """ + def conformer(seed: int, energy_ev: float) -> Chem.Mol: m = Chem.AddHs(Chem.MolFromSmiles("CCCCCCO")) # flexible chain AllChem.EmbedMolecule(m, randomSeed=seed) AllChem.MMFFOptimizeMolecule(m) @@ -403,15 +512,121 @@ def conformer(seed: float, energy_ev: float) -> Chem.Mol: mols = [conformer(42, -12.0), conformer(7, -11.0), conformer(123, -10.0)] - original_result = filter_unique(mols, crit=0.3) - optimized_result = filter_unique_optimized( - mols, - rmsd_threshold=0.3, - energy_cluster_window=100.0 # Large window = single cluster = same behavior + assert len(filter_unique_optimized(mols, rmsd_threshold=0.3)) == 3 + # A single cluster (huge window) must give the same answer -- that + # equivalence is the whole justification for the energy partitioning. + assert len( + filter_unique_optimized( + mols, rmsd_threshold=0.3, energy_cluster_window=100.0 + ) + ) == 3 + + def test_a_malformed_mixed_list_keeps_the_recorded_three(self): + """The input the two filters used to DISAGREE about. + + An energy-bearing conformer, an energy-less duplicate of it, an + energy-less record of a different compound, and a distinct third + compound. The legacy filter kept three of the four -- merging the + energy-less ethanol into the one that has an energy -- and so must this + one. + """ + a = _create_mol_with_energy("CCO", -100.0) + a.SetProp("_Name", "ethanol_with_energy") + b = _energyless("CCO", seed=42, name="ethanol_no_energy") + c = _energyless("CCCCO", seed=11, name="butanol_no_energy") + d = _create_mol_with_energy("CCCCCCO", -80.0) + d.SetProp("_Name", "heptanol_with_energy") + + kept = filter_unique_optimized([a, b, c, d], rmsd_threshold=0.3) + + assert {m.GetProp("_Name") for m in kept} == { + "ethanol_with_energy", + "butanol_no_energy", + "heptanol_with_energy", + } + + def test_an_rmsd_threshold_sweep_straddling_a_measured_pair(self): + """A threshold below the pair's actual RMSD keeps both; above, merges. + + Recorded from both filters. Sweeping across the *measured* RMSD -- not + a fixed number -- is what makes this fail for a filter that ignores + ``rmsd_threshold`` altogether, which an equal-length comparison at a + single threshold would not. + """ + from rdkit.Chem import rdMolAlign + + first = _create_mol_with_energy("CCCCCCCC", -100.0) # octane + second = Chem.Mol(first) + AllChem.EmbedMolecule(second, randomSeed=99) + AllChem.MMFFOptimizeMolecule(second) + set_e_tot_from_ev(second, -100.0) # energies agree -> RMSD decides + second.SetProp("Converged", "true") + + rmsd = rdMolAlign.GetBestRMS(Chem.RemoveHs(first), Chem.RemoveHs(second)) + assert rmsd > 0.05, "test premise: the conformers must be distinct" + + for crit, expected in ((rmsd / 2, 2), (rmsd * 2, 1)): + kept = filter_unique_optimized([first, second], rmsd_threshold=crit) + assert len(kept) == expected, f"at rmsd_threshold={crit}" + + def test_two_energyless_conformers_dedup_by_rmsd_alone(self): + """Ported from the legacy filter's own suite. + + Neither record carries ``E_tot``, so the energy half of the duplicate + criterion cannot apply and RMSD alone decides. The legacy filter kept + one; so must this one. + """ + mol = _energyless("CCO", seed=42) + duplicate = Chem.Mol(mol) + duplicate.SetProp("Converged", "true") + + assert len(filter_unique_optimized([mol, duplicate], rmsd_threshold=0.3)) == 1 + + +class TestConvergencePropertyAbsenceFiltersLikeTrue: + """A whole SDF that carries no 'Converged' property must filter exactly as + the same records marked Converged=True do. + + Ported from the legacy filter's suite (it lived beside the validation tests + because ``filter_unique`` used to live in ``utils/validation.py``). Only + ``batchopt`` writes ``Converged``; an ``opt_geometry`` output, an + ORCA/Gaussian export or a hand-built conformer set carries none, and + treating that as "did not converge" deleted every record. + """ + + _SDF = os.path.join( + os.path.dirname(os.path.abspath(__file__)), "files", "example.sdf" + ) + + def _mols(self) -> list[Chem.Mol]: + supp = Chem.SDMolSupplier(self._SDF, removeHs=False) + return [mol for mol in supp if mol is not None] + + def test_an_unflagged_file_keeps_what_a_flagged_one_keeps(self): + flagged = self._mols() + for mol in flagged: + mol.SetProp("Converged", "True") + expected = len(filter_unique_optimized(flagged, rmsd_threshold=0.3)) + assert expected >= 1, "test premise: the flagged file must keep something" + + unflagged = self._mols() + for mol in unflagged: + mol.ClearProp("Converged") + assert not mol.HasProp("Converged") + + result = filter_unique_optimized(unflagged, rmsd_threshold=0.3) + assert len(result) == expected, ( + f"{len(unflagged)} record(s) with no 'Converged' property kept " + f"{len(result)}, but the same records marked Converged=True keep " + f"{expected}" ) - # Same molecule, well-defined RMSD: both implementations must agree. - assert len(original_result) == len(optimized_result) + def test_an_explicit_false_still_empties_the_selection(self): + """The inverse: absence is not failure, but a stated failure is.""" + mols = self._mols() + for mol in mols: + mol.SetProp("Converged", "False") + assert filter_unique_optimized(mols, rmsd_threshold=0.3) == [] def test_filter_within_cluster_removehs_is_linear_and_nondestructive(monkeypatch): @@ -482,3 +697,101 @@ def boom(*a, **k): cluster = [make("a", -1.0), make("b", -0.9)] kept = _filter_within_cluster(cluster, rmsd_threshold=0.3) assert len(kept) == 2 # incomparable pair must NOT be dropped + + +class TestTheFilterSaysWhyItDroppedThings: + """``filter_conformers`` reports a count per reason, not just a survivor list. + + Returning a bare list is what let ``ranking`` tell a user "No structure + converged" for a species whose conformers were every one of them dropped + for *stereochemistry* -- a message that points at ``--opt-steps`` and + ``--convergence-threshold`` for a problem neither can fix. + """ + + def test_each_reason_is_counted_under_its_own_name(self): + good = _create_mol_with_energy("CCO", -100.0) + unconverged = _create_mol_with_energy("CCO", -99.0, converged=False) + stereo_changed = _create_mol_with_energy("CCO", -98.0) + stereo_changed.SetProp("Stereo_changed", "true") + broken = _create_mol_with_energy("CC", -97.0) + conf = broken.GetConformer() + pos = conf.GetAtomPosition(0) + conf.SetAtomPosition(0, (pos.x + 5.0, pos.y, pos.z)) + duplicate = _create_mol_with_energy("CCO", -100.0) # same as `good` + + result = filter_conformers( + [None, good, unconverged, stereo_changed, broken, duplicate], + rmsd_threshold=0.3, + ) + + assert [m is good for m in result.kept] == [True] + assert result.dropped == { + "unparsed": 1, + "unconverged": 1, + "stereochemistry": 1, + "connectivity": 1, + "duplicate": 1, + } + + def test_nothing_dropped_reports_nothing(self): + """The inverse: a clean input must not manufacture a reason. + + A result object that always carried a non-empty ``dropped`` would make + every ranking message name a cause that did not happen. + """ + result = filter_conformers( + [_create_mol_with_energy("CCO", -100.0)], rmsd_threshold=0.3 + ) + assert result.dropped == {} + assert result.reasons == () + assert result.summary() == "" + + def test_summary_names_every_reason_that_fired_in_declared_order(self): + result = FilterResult( + kept=[], + dropped={"duplicate": 3, "unconverged": 2, "stereochemistry": 1}, + ) + assert result.reasons == ("unconverged", "stereochemistry", "duplicate") + assert result.summary() == ( + "2 marked Converged=false, " + "1 changed stereochemistry during optimization, " + "3 duplicates of a kept conformer" + ) + + def test_zero_counts_are_not_reported_as_reasons(self): + result = FilterResult(kept=[], dropped={"unconverged": 0, "duplicate": 2}) + assert result.reasons == ("duplicate",) + assert result.summary() == "2 duplicates of a kept conformer" + + def test_an_unknown_reason_is_refused_at_construction(self): + """DROP_REASONS is the vocabulary, enforced. + + Without this, a producer misspelling a reason contributes a drop that + ``summary()`` silently omits, so a user is told fewer conformers went + missing than actually did. + """ + with pytest.raises(ValueError, match="unknown filter drop reason"): + FilterResult(kept=[], dropped={"unconvrged": 1}) + + def test_every_declared_reason_has_a_phrase(self): + """A reason with no phrase would raise KeyError from inside summary() + -- while reporting a diagnostic, which is the worst place to fail.""" + for reason in DROP_REASONS: + assert FilterResult(kept=[], dropped={reason: 1}).summary() + + def test_truncation_is_not_a_drop_reason(self): + """`k` cutting the list short is selection, not a filter drop. + + Nothing is missing to explain: those conformers are valid and unique, + they lost the ranking. Counting them would make every ``k=1`` run report + drops it should not. + """ + assert "truncated" not in DROP_REASONS + distinct = [ + _create_mol_with_energy("CCO", -100.0), + _create_mol_with_energy("CCCO", -90.0), + _create_mol_with_energy("CCCCO", -80.0), + ] + result = filter_conformers(distinct, rmsd_threshold=0.3) + assert len(result.kept) == 3 + assert result.dropped == {} diff --git a/tests/test_id_mapping.py b/tests/test_id_mapping.py index 313fe297..800086c2 100644 --- a/tests/test_id_mapping.py +++ b/tests/test_id_mapping.py @@ -133,6 +133,14 @@ def test_encode_ids_refuses_to_overwrite_an_existing_file(self, tmp_path): `out_dir` below), but this check keeps the guarantee attached to the function itself, so a caller taking the default location cannot reintroduce the defect. + + The refusal now comes from the shared + `Auto3D.utils.output_guard.check_output_overwrite` (hence the + "already exists" wording) rather than a bespoke message here, so + `encode_ids`, `decode_ids` and `tautomer.select_tautomers` state the + same policy the same way -- and `overwrite=True` can lift it, which the + unconditional refusal this replaced offered no way to do. + `tests/test_output_overwrite_gates.py` covers both directions. """ from Auto3D.exceptions import ConfigurationError @@ -141,7 +149,7 @@ def test_encode_ids_refuses_to_overwrite_an_existing_file(self, tmp_path): users_file = tmp_path / "mols_encoded.smi" users_file.write_bytes(b"IRREPLACEABLE USER DATA\n") - with pytest.raises(ConfigurationError, match="would overwrite"): + with pytest.raises(ConfigurationError, match="already exists"): encode_ids(str(p)) assert users_file.read_bytes() == b"IRREPLACEABLE USER DATA\n" diff --git a/tests/test_output_overwrite_gates.py b/tests/test_output_overwrite_gates.py new file mode 100644 index 00000000..f6ee127c --- /dev/null +++ b/tests/test_output_overwrite_gates.py @@ -0,0 +1,259 @@ +#!/usr/bin/env python +"""Overwrite gates on the public writers that derive their own output names. + +``Chem.SDWriter(path)`` and ``open(path, "w")`` truncate on open, and four +public functions used to do that with no consent gate at all: + +* ``Auto3D.tautomer.select_tautomers`` -- ``select_tautomers("/data/results.sdf", + k=1)`` replaced ``/data/results_top_tautomers.sdf``, a name it invented, with + this call's selection. +* ``Auto3D.id_mapping.decode_ids`` -- same shape, for ``_out.sdf``. +* ``Auto3D.utils.smi_io.smiles2smi`` -- the caller named the file here, so the + gate defaults *open*; what matters is that it can be closed. +* ``Auto3D.id_mapping.encode_ids`` -- refused unconditionally, with no way for a + caller to say yes. + +The policy: permissive where the caller named the file, restrictive where Auto3D +invented the name. Every gate is keyword-only, and every one is asserted in both +directions -- a gate that refused legitimate writes would satisfy the "refuses" +half of each pair and break the pipeline. +""" +from __future__ import annotations + +import pytest +from rdkit import Chem +from rdkit.Chem import AllChem + +from Auto3D.exceptions import ConfigurationError +from Auto3D.id_mapping import decode_ids, encode_ids +from Auto3D.tautomer import select_tautomers +from Auto3D.utils.energy import set_e_tot_from_ev +from Auto3D.utils.smi_io import smiles2smi + +PRECIOUS = b"IRREPLACEABLE USER DATA\n" + + +def _tautomer_sdf(path, names=("mol@taut1", "mol@taut2")) -> str: + """An SDF shaped the way ``select_tautomers`` expects: ``id@tautN`` names.""" + with Chem.SDWriter(str(path)) as writer: + for i, name in enumerate(names): + mol = Chem.AddHs(Chem.MolFromSmiles("CCO")) + AllChem.EmbedMolecule(mol, randomSeed=42 + i) + mol.SetProp("_Name", name) + set_e_tot_from_ev(mol, -10.0 - i) + writer.write(mol) + return str(path) + + +def _encoded_sdf(path) -> str: + """An SDF shaped the way ``decode_ids`` expects: numeric name + ID prop.""" + with Chem.SDWriter(str(path)) as writer: + mol = Chem.AddHs(Chem.MolFromSmiles("CCO")) + AllChem.EmbedMolecule(mol, randomSeed=42) + mol.SetProp("_Name", "0") + mol.SetProp("ID", "0_0_0") + writer.write(mol) + return str(path) + + +class TestSelectTautomers: + """Auto3D invents ``_top_tautomers.sdf``, so the gate defaults shut.""" + + def test_it_refuses_to_replace_an_existing_top_tautomers_file(self, tmp_path): + """The concrete hazard from the follow-ups list.""" + sdf = _tautomer_sdf(tmp_path / "results.sdf") + derived = tmp_path / "results_top_tautomers.sdf" + derived.write_bytes(PRECIOUS) + + with pytest.raises(ConfigurationError, match="already exists"): + select_tautomers(sdf, k=1) + + assert derived.read_bytes() == PRECIOUS + + def test_it_writes_when_the_derived_name_is_free(self, tmp_path): + """Inverse: the ordinary case must still go through untouched. + + Every route Auto3D itself takes lands here -- ``get_stable_tautomers`` + passes ``main()``'s output from a job directory created fresh for that + run -- so a gate that refused unconditionally would break the whole + tautomer pipeline while passing the test above. + """ + sdf = _tautomer_sdf(tmp_path / "results.sdf") + out = select_tautomers(sdf, k=1) + + assert out == str(tmp_path / "results_top_tautomers.sdf") + assert len( + [m for m in Chem.SDMolSupplier(out, removeHs=False) if m is not None] + ) == 1 + + def test_overwrite_true_replaces_it(self, tmp_path): + """Inverse: the gate is a consent gate the caller can lift.""" + sdf = _tautomer_sdf(tmp_path / "results.sdf") + derived = tmp_path / "results_top_tautomers.sdf" + derived.write_bytes(PRECIOUS) + + out = select_tautomers(sdf, k=1, overwrite=True) + + assert out == str(derived) + assert derived.read_bytes() != PRECIOUS + + def test_overwrite_is_keyword_only(self, tmp_path): + """``select_tautomers(sdf, k, window)`` is the documented positional + signature; a fourth positional must not silently become ``overwrite``.""" + sdf = _tautomer_sdf(tmp_path / "results.sdf") + with pytest.raises(TypeError): + select_tautomers(sdf, 1, None, True) + + +class TestDecodeIds: + """Auto3D invents ``_out.sdf``, so the gate defaults shut.""" + + def test_it_refuses_to_replace_an_existing_out_file(self, tmp_path): + sdf = _encoded_sdf(tmp_path / "mols_3d_encoded.sdf") + derived = tmp_path / "mols_out.sdf" + derived.write_bytes(PRECIOUS) + + with pytest.raises(ConfigurationError, match="already exists"): + decode_ids(sdf, {"mol_a": 0}) + + assert derived.read_bytes() == PRECIOUS + + def test_it_writes_when_the_derived_name_is_free(self, tmp_path): + """Inverse: this is the call ``WorkflowOrchestrator`` makes on every run. + + It writes into the job directory it created with a bare ``mkdir()`` + moments earlier, so the derived name is always free there -- but only if + the gate permits a free name. + """ + sdf = _encoded_sdf(tmp_path / "mols_3d_encoded.sdf") + out = decode_ids(sdf, {"mol_a": 0}) + + assert out == str(tmp_path / "mols_out.sdf") + written = [ + m for m in Chem.SDMolSupplier(out, removeHs=False) if m is not None + ] + assert [m.GetProp("_Name") for m in written] == ["mol_a"] + + def test_overwrite_true_replaces_it(self, tmp_path): + sdf = _encoded_sdf(tmp_path / "mols_3d_encoded.sdf") + derived = tmp_path / "mols_out.sdf" + derived.write_bytes(PRECIOUS) + + out = decode_ids(sdf, {"mol_a": 0}, overwrite=True) + + assert out == str(derived) + assert derived.read_bytes() != PRECIOUS + + def test_overwrite_is_keyword_only(self, tmp_path): + sdf = _encoded_sdf(tmp_path / "mols_3d_encoded.sdf") + with pytest.raises(TypeError): + decode_ids(sdf, {"mol_a": 0}, True) + + +class TestEncodeIds: + """Already refused unconditionally; now it refuses *consistently*. + + Same keyword, same default, same exception and same message as the other + three -- and, unlike before, a caller who means it can say so. + """ + + def test_it_still_refuses_by_default(self, tmp_path): + smi = tmp_path / "mols.smi" + smi.write_text("CCO a\n") + derived = tmp_path / "mols_encoded.smi" + derived.write_bytes(PRECIOUS) + + with pytest.raises(ConfigurationError, match="already exists"): + encode_ids(str(smi)) + + assert derived.read_bytes() == PRECIOUS + + def test_overwrite_true_replaces_it(self, tmp_path): + """The half that was impossible before: there was no way to consent.""" + smi = tmp_path / "mols.smi" + smi.write_text("CCO a\n") + derived = tmp_path / "mols_encoded.smi" + derived.write_bytes(PRECIOUS) + + new_path, mapping = encode_ids(str(smi), overwrite=True) + + assert new_path == str(derived) + assert mapping == {"a": 0} + assert derived.read_text() == "CCO 0\n" + + def test_it_writes_when_the_derived_name_is_free(self, tmp_path): + """Inverse: the pipeline's own call, which must keep working.""" + smi = tmp_path / "mols.smi" + smi.write_text("CCO a\nCCC b\n") + new_path, mapping = encode_ids(str(smi)) + assert mapping == {"a": 0, "b": 1} + assert new_path == str(tmp_path / "mols_encoded.smi") + + def test_overwrite_is_keyword_only(self, tmp_path): + smi = tmp_path / "mols.smi" + smi.write_text("CCO a\n") + with pytest.raises(TypeError): + encode_ids(str(smi), None, True) + + +class TestSmiles2Smi: + """The caller named this file, so the gate defaults OPEN. + + ``smiles2mols`` writes it into a ``TemporaryDirectory`` on every call, and a + caller who passes an explicit path has already chosen it. Defaulting shut + here would be a gate on a decision the caller already made. + """ + + def test_the_default_still_overwrites(self, tmp_path): + """Inverse first, because this is the back-compat guarantee.""" + out = tmp_path / "mols.smi" + out.write_bytes(PRECIOUS) + + assert smiles2smi(["CCO"], str(out)) == str(out) + assert out.read_bytes() != PRECIOUS + assert out.read_text().startswith("CCO ") + + def test_overwrite_false_refuses(self, tmp_path): + out = tmp_path / "mols.smi" + out.write_bytes(PRECIOUS) + + with pytest.raises(ConfigurationError, match="already exists"): + smiles2smi(["CCO"], str(out), overwrite=False) + + assert out.read_bytes() == PRECIOUS + + def test_overwrite_false_still_writes_a_free_path(self, tmp_path): + out = tmp_path / "mols.smi" + assert smiles2smi(["CCO"], str(out), overwrite=False) == str(out) + assert out.read_text().startswith("CCO ") + + def test_overwrite_is_keyword_only(self, tmp_path): + with pytest.raises(TypeError): + smiles2smi(["CCO"], str(tmp_path / "mols.smi"), False) + + +class TestTheGateRefusesBeforeDoingTheWork: + """A gate applied just before the write is a gate that wasted the run. + + ``select_tautomers`` reads and groups the whole input before it opens the + writer; refusing only there means the user waits for the work, then loses + it. The check must come first. + """ + + def test_select_tautomers_checks_before_reading_the_input( + self, tmp_path, monkeypatch + ): + import Auto3D.tautomer as tautomer + + sdf = _tautomer_sdf(tmp_path / "results.sdf") + (tmp_path / "results_top_tautomers.sdf").write_bytes(PRECIOUS) + + def _never(*args, **kwargs): + raise AssertionError( + "select_tautomers read its input before checking the output path" + ) + + monkeypatch.setattr(tautomer.Chem, "SDMolSupplier", _never) + + with pytest.raises(ConfigurationError, match="already exists"): + select_tautomers(sdf, k=1) diff --git a/tests/test_padding_invariance.py b/tests/test_padding_invariance.py index a96eb50b..cbe54d4a 100644 --- a/tests/test_padding_invariance.py +++ b/tests/test_padding_invariance.py @@ -2,9 +2,11 @@ Every engine masks padded slots differently: AIMNet2 pads species with 0 and relies on Z=0 being unused, while ANI2x/ANI2xt pad with -1 and rely on that -being torchani's masked-atom sentinel. batch_opt/ANI2xt_no_rep.py:167-172 -documents that second assumption as unverified. These tests verify it (audit -M32, C13). +being torchani's masked-atom sentinel. ``ANI2xt.forward`` in +``batch_opt/ANI2xt_no_rep.py`` documents that second assumption (species +== -1 surviving the periodic-table remap unchanged and relying on +TorchANI's masked-atom convention) as depended-upon but not independently +verified there. These tests verify it (audit M32, C13). """ from __future__ import annotations @@ -30,7 +32,8 @@ class TestPaddingInvariance: # padded-vs-solo AIMNet2 invariant tests/test_model_adapter.py:249-250 # already asserts with atol=1e-2 eV, and ANI2xt's float32 output caps # usable precision at ~float32 ULP (~4e-3 eV) at typical total-energy - # magnitudes per src/Auto3D/batch_opt/ANI2xt_no_rep.py:148-155. 1e-6 would + # magnitudes per the float32-precision Note in ``ANI2xt.forward``'s + # docstring (``src/Auto3D/batch_opt/ANI2xt_no_rep.py``). 1e-6 would # demand sub-ULP reproducibility and flake on a correct model. ANI2x # (torchani, periodic-table indexing) shares ANI2xt's float32 output and # the same -1 species_pad convention, so it gets the same 1e-3 budget diff --git a/tests/test_ranking.py b/tests/test_ranking.py index de7be5f3..795d5caf 100644 --- a/tests/test_ranking.py +++ b/tests/test_ranking.py @@ -47,10 +47,10 @@ def _write_mols_to_sdf(mols: list[Chem.Mol], filepath: str) -> None: class TestConformerRankerWithOptimizedFiltering: - """Tests for ConformerRanker with optimized filtering.""" + """Tests for ConformerRanker's conformer filtering.""" - def test_ranker_with_optimized_filtering_default(self, tmp_path): - """ConformerRanker should use optimized filtering by default.""" + def test_ranker_deduplicates_identical_conformers(self, tmp_path): + """ConformerRanker dedups through the single conformer filter.""" from Auto3D.ranking import ConformerRanker # Create test molecules - all with same SMILES root name. @@ -72,82 +72,66 @@ def test_ranker_with_optimized_filtering_default(self, tmp_path): k=5, ) - # Default should use optimized filtering - assert ranker.use_optimized_filtering is True - results = ranker.run() # With close energies, all in same cluster, should deduplicate to 1 assert len(results) == 1 - def test_ranker_with_legacy_filtering_fallback(self, tmp_path): - """ConformerRanker should support legacy filtering when explicitly requested.""" - from Auto3D.ranking import ConformerRanker + def test_use_optimized_filtering_is_gone_not_silently_ignored(self, tmp_path): + """The flag that selected between two filter implementations is removed. - # Same structure and near-identical energy (within the duplicate energy - # tolerance) -> one unique structure. - mol1 = _create_mol_with_energy("C", -10.0, "mol_1") - mol2 = _create_mol_with_energy("C", -10.005, "mol_2") + There were two conformer filters with the same duplicate criterion, + chosen by ``use_optimized_filtering``, and they had drifted on malformed + input. One filter now, and the flag must raise rather than be swallowed + by ``**kwargs`` -- a caller passing ``False`` deserves to learn the + legacy path is gone instead of silently getting the other one. + """ + from Auto3D.ranking import ConformerRanker input_path = str(tmp_path / "input.sdf") - output_path = str(tmp_path / "output.sdf") - _write_mols_to_sdf([mol1, mol2], input_path) - - ranker = ConformerRanker( - input_path=input_path, - out_path=output_path, - threshold=0.3, - k=5, - use_optimized_filtering=False, - ) - - assert ranker.use_optimized_filtering is False - - results = ranker.run() - # Same behavior, one unique structure - assert len(results) == 1 - - def test_ranker_optimized_vs_legacy_produce_same_results(self, tmp_path): - """Optimized and legacy filtering should produce equivalent results. - - With a large energy_cluster_window, optimized should behave like legacy. + _write_mols_to_sdf([_create_mol_with_energy("C", -10.0, "mol_1")], input_path) + + with pytest.raises(TypeError, match="use_optimized_filtering"): + ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / "output.sdf"), + threshold=0.3, + k=5, + use_optimized_filtering=False, + ) + + def test_a_single_cluster_gives_the_same_answer(self, tmp_path): + """A window wide enough to make one cluster must not change the result. + + That equivalence is the entire justification for partitioning the energy + axis at all, and it used to be asserted by comparing against the legacy + all-pairs filter (deleted in 4.1.0). Asserted directly now. """ from Auto3D.ranking import ConformerRanker # Identical molecules (same structure AND near-identical energy, within # the duplicate energy tolerance) - these should be deduplicated. - mol1 = _create_mol_with_energy("CCCC", -10.0, "a_1") - mol2 = _create_mol_with_energy("CCCC", -10.005, "a_2") # same structure & energy - mol3 = _create_mol_with_energy("CCCC", -10.008, "a_3") # same structure & energy + mols = [ + _create_mol_with_energy("CCCC", -10.0, "a_1"), + _create_mol_with_energy("CCCC", -10.005, "a_2"), + _create_mol_with_energy("CCCC", -10.008, "a_3"), + ] input_path = str(tmp_path / "input.sdf") - output_optimized = str(tmp_path / "output_optimized.sdf") - output_legacy = str(tmp_path / "output_legacy.sdf") - _write_mols_to_sdf([mol1, mol2, mol3], input_path) - - # Test with optimized filtering - use large energy window to match legacy behavior - ranker_optimized = ConformerRanker( - input_path=input_path, - out_path=output_optimized, - threshold=0.3, - k=5, - use_optimized_filtering=True, - energy_cluster_window=100.0, # Large window = single cluster = legacy behavior - ) - results_optimized = ranker_optimized.run() + _write_mols_to_sdf(mols, input_path) - # Test with legacy filtering - ranker_legacy = ConformerRanker( - input_path=input_path, - out_path=output_legacy, - threshold=0.3, - k=5, - use_optimized_filtering=False, - ) - results_legacy = ranker_legacy.run() + counts = [] + for name, window in (("default", None), ("single_cluster", 100.0)): + kwargs = {} if window is None else {"energy_cluster_window": window} + counts.append(len(ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / f"output_{name}.sdf"), + threshold=0.3, + k=5, + **kwargs, + ).run())) - # Should have same number of results - all identical molecules deduplicated to 1 - assert len(results_optimized) == len(results_legacy) - assert len(results_optimized) == 1 # All identical molecules should be deduplicated + assert counts[0] == counts[1] + assert counts[0] == 1 # all identical molecules deduplicate to one def test_energy_cluster_window_parameter(self, tmp_path): """Ranker should accept energy_cluster_window parameter for optimized filtering.""" @@ -198,7 +182,6 @@ def test_top_k_with_optimized_filtering(self, tmp_path): out_path=output_path, threshold=0.3, k=2, - use_optimized_filtering=True, ) df = pd.DataFrame({ @@ -367,7 +350,6 @@ def test_top_window_with_optimized_filtering(self, tmp_path): out_path=output_path, threshold=0.3, window=1.0, # 1 kcal/mol window - use_optimized_filtering=True, ) df = pd.DataFrame({ @@ -789,3 +771,453 @@ def test_enumerate_isomer_false_returns_both_molecules(self, tmp_path): ] assert {m.GetProp("_Name") for m in written} == {key, key_2} assert os.path.getsize(output) > 0 + + +class TestNothingSelectedSaysWhy: + """"No structure converged" used to be the message for every empty group. + + ``ranking`` logged it whether the conformers were dropped for convergence, + for stereochemistry (an optimization that inverted a center, so the geometry + no longer matches the title) or for connectivity (a structure that fell + apart). Two of those three point at the input or the chemistry, and both + were reported as an optimizer convergence problem -- sending the reader to + ``--opt-steps`` and ``--convergence-threshold`` for something neither can + fix. + """ + + @staticmethod + def _group(mols: list[Chem.Mol], name: str = "probe"): + """A one-species ranking group, shaped the way ``run`` builds them.""" + import pandas as pd + + return pd.DataFrame({ + "names": [name] * len(mols), + "energies": [e_tot_ev(m) for m in mols], + "mols": mols, + }) + + @staticmethod + def _ranker(tmp_path, **kwargs): + from Auto3D.ranking import ConformerRanker + + return ConformerRanker( + input_path=str(tmp_path / "in.sdf"), + out_path=str(tmp_path / "out.sdf"), + threshold=0.3, + **kwargs, + ) + + @staticmethod + def _stereo_changed(energy: float, name: str) -> Chem.Mol: + from Auto3D.utils.stereo_check import STEREO_CHANGED_PROP + + mol = _create_mol_with_energy("C/C=C/CCO", energy, name) + mol.SetProp(STEREO_CHANGED_PROP, "true") + return mol + + @staticmethod + def _messages(caplog) -> list[str]: + return [r.getMessage() for r in caplog.records] + + def test_a_stereo_dropped_species_is_not_called_unconverged(self, tmp_path, caplog): + import logging + + mols = [self._stereo_changed(-10.0, "probe_0_0"), + self._stereo_changed(-9.0, "probe_0_1")] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=5) == [] + + messages = self._messages(caplog) + assert not any("No structure converged" in m for m in messages), ( + f"conformers dropped for stereochemistry were reported as an " + f"optimizer convergence failure: {messages}" + ) + assert any( + "probe" in m and "stereochemistry" in m for m in messages + ), f"the real reason was never named: {messages}" + + def test_the_k1_fast_path_reports_the_same_reason(self, tmp_path, caplog): + """k=1 bypasses the RMSD dedup entirely; the diagnostic must not + depend on which k the user asked for.""" + import logging + + mols = [self._stereo_changed(-10.0, "probe_0_0")] + ranker = self._ranker(tmp_path, k=1) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=1) == [] + + messages = self._messages(caplog) + assert not any("No structure converged" in m for m in messages), messages + assert any("stereochemistry" in m for m in messages), messages + + def test_a_connectivity_dropped_species_names_connectivity(self, tmp_path, caplog): + import logging + + from Auto3D.utils.connectivity import check_connectivity + + def broken(energy: float, name: str) -> Chem.Mol: + mol = _create_mol_with_energy("CC", energy, name) + conf = mol.GetConformer() + pos = conf.GetAtomPosition(0) + conf.SetAtomPosition(0, (pos.x + 5.0, pos.y, pos.z)) + assert check_connectivity(mol) is False, "test premise" + return mol + + mols = [broken(-10.0, "probe_0_0"), broken(-9.0, "probe_0_1")] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=5) == [] + + messages = self._messages(caplog) + assert not any("No structure converged" in m for m in messages), messages + assert any("broken or newly formed bonds" in m for m in messages), messages + + def test_the_literal_survives_when_convergence_is_the_sole_reason( + self, tmp_path, caplog + ): + """The inverse, and the reason the assertions above are safe. + + A change that simply stopped emitting "No structure converged" would + satisfy every test above while destroying the message users and their + log-scraping scripts have matched on since Auto3D 1.x. When convergence + IS the sole reason, the exact wording must still appear. + """ + import logging + + mols = [ + _create_mol_with_energy("CCO", -10.0, "probe_0_0", converged=False), + _create_mol_with_energy("CCO", -9.0, "probe_0_1", converged=False), + ] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=5) == [] + + assert any( + m == "No structure converged for probe." for m in self._messages(caplog) + ), self._messages(caplog) + + def test_the_literal_is_not_emitted_alongside_another_reason( + self, tmp_path, caplog + ): + """"Only when that is the sole reason": a mixed group must not claim + convergence.""" + import logging + + mols = [ + _create_mol_with_energy("C/C=C/CCO", -10.0, "probe_0_0", converged=False), + self._stereo_changed(-9.0, "probe_0_1"), + ] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert ranker.top_k(self._group(mols), k=5) == [] + + messages = self._messages(caplog) + assert not any("No structure converged" in m for m in messages), messages + # Both reasons are named, so the reader sees the whole accounting. + assert any( + "Converged=false" in m and "stereochemistry" in m for m in messages + ), messages + + def test_a_successful_selection_logs_no_complaint(self, tmp_path, caplog): + """The other inverse: a group that DOES select must stay silent. + + A message emitted unconditionally would pass the "names the reason" + tests above and spam every ordinary run. + """ + import logging + + mols = [_create_mol_with_energy("CCO", -10.0, "probe_0_0")] + ranker = self._ranker(tmp_path, k=5) + + with caplog.at_level(logging.INFO, logger="Auto3D.ranking"): + assert len(ranker.top_k(self._group(mols), k=5)) == 1 + + messages = self._messages(caplog) + assert not any("No structure" in m for m in messages), messages + + def test_top_window_merges_the_window_into_the_same_accounting(self, tmp_path): + """The energy window is the one drop reason ``top_window`` owns. + + It goes into the same run-level tally as the filter's own counts, so + ``run``'s summary is one accounting of the whole selection rather than + two partial ones. + """ + ranker = self._ranker(tmp_path, window=1.0) + + # Two distinct compounds, 5 eV apart: far outside a 1 kcal/mol window. + mols = [ + _create_mol_with_energy("CCO", -10.0, "probe_0_0"), + _create_mol_with_energy("CCCCCCO", -5.0, "probe_0_1"), + ] + kept = ranker.top_window(self._group(mols), window=1.0) + + assert len(kept) == 1, "the second conformer is outside the window" + assert ranker._drop_totals == {"energy_window": 1} + + def test_a_wide_window_records_no_window_drop(self, tmp_path): + """The inverse: a window nothing falls outside of must not report one.""" + ranker = self._ranker(tmp_path, window=1000.0) + mols = [ + _create_mol_with_energy("CCO", -10.0, "probe_0_0"), + _create_mol_with_energy("CCCCCCO", -5.0, "probe_0_1"), + ] + assert len(ranker.top_window(self._group(mols), window=1000.0)) == 2 + assert ranker._drop_totals == {} + + def test_the_run_level_warning_names_the_reasons_that_fired( + self, tmp_path, caplog + ): + """``run``'s "selected 0 structures" warning used to list every reason + it MIGHT have been. + + The text was "N record(s) are marked Converged=false and the rest were + dropped by the connectivity, stereochemistry or energy-window filters" + -- a hand-maintained disjunction that left the reader to work out which + of the three actually happened, on the one message that reaches a direct + API caller's stderr. Here nothing is unconverged and everything is + stereo-changed, so naming convergence would be wrong. + """ + import logging + + from Auto3D.ranking import ConformerRanker + + mols = [self._stereo_changed(-10.0, "probe_0_0"), + self._stereo_changed(-9.0, "probe_0_1")] + input_path = str(tmp_path / "in.sdf") + _write_mols_to_sdf(mols, input_path) + + ranker = ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / "out.sdf"), + threshold=0.3, + k=5, + ) + with caplog.at_level(logging.WARNING, logger="Auto3D.ranking"): + assert ranker.run() == [] + + warnings = [ + r.getMessage() for r in caplog.records if r.levelno >= logging.WARNING + ] + assert any("Selected 0 structures from 2 record(s)" in m for m in warnings), ( + warnings + ) + assert any( + "changed stereochemistry during optimization" in m for m in warnings + ), warnings + assert not any("Converged=false" in m for m in warnings), ( + f"nothing was unconverged, yet convergence was named: {warnings}" + ) + + def test_the_run_level_warning_still_names_convergence_when_it_applies( + self, tmp_path, caplog + ): + """The inverse: records dropped before grouping (for convergence, or + because RDKit could not parse them) are counted in the same tally.""" + import logging + + from Auto3D.ranking import ConformerRanker + + mols = [ + _create_mol_with_energy("CCO", -10.0, "probe_0_0", converged=False), + _create_mol_with_energy("CCO", -9.0, "probe_0_1", converged=False), + ] + input_path = str(tmp_path / "in.sdf") + _write_mols_to_sdf(mols, input_path) + + ranker = ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / "out.sdf"), + threshold=0.3, + k=5, + ) + with caplog.at_level(logging.WARNING, logger="Auto3D.ranking"): + assert ranker.run() == [] + + warnings = [ + r.getMessage() for r in caplog.records if r.levelno >= logging.WARNING + ] + assert any("2 marked Converged=false" in m for m in warnings), warnings + + def test_the_tally_is_reset_between_runs(self, tmp_path): + """A ranker reused for a second file must not report the first's drops.""" + from Auto3D.ranking import ConformerRanker + + # Stereo-changed records pass run()'s own convergence check and are + # dropped by the FILTER, so they land in the tally -- unlike an + # unconverged record, which run() drops before grouping and counts + # separately. Getting that wrong makes this test vacuous. + dirty = str(tmp_path / "dirty.sdf") + _write_mols_to_sdf( + [self._stereo_changed(-10.0, "a_0_0"), + self._stereo_changed(-9.0, "a_0_1")], + dirty, + ) + clean = str(tmp_path / "clean.sdf") + _write_mols_to_sdf([_create_mol_with_energy("CCO", -10.0, "b_0_0")], clean) + + ranker = ConformerRanker( + input_path=dirty, out_path=str(tmp_path / "out.sdf"), threshold=0.3, k=5, + ) + assert ranker.run() == [] + assert ranker._drop_totals == {"stereochemistry": 2}, "test premise" + + ranker.input_path = clean + assert len(ranker.run()) == 1 + assert ranker._drop_totals == {}, ( + "the second run reported the first run's drops" + ) + + +class TestSelectorDispatchRegistry: + """``run`` dispatches through a registry checked against ``config.py``. + + It used to be a hand-written ``if self.k: ... elif self.window: ...``. A + third selector added to ``Auto3D.config.SELECTOR_FIELDS`` would then be + accepted by ``Auto3DOptions``, accepted by ``CLIConfig``, accepted by + ``check_selectors_mutually_exclusive`` -- and silently ignored here, falling + through to "Parameter k or window needs to be specified" even though the + user had specified one. + """ + + def test_the_registry_matches_the_authoritative_field_list(self): + from Auto3D.config import SELECTOR_FIELDS + from Auto3D.ranking import _SELECTORS + + assert set(_SELECTORS) == set(SELECTOR_FIELDS) + + def test_every_mapped_method_exists(self): + from Auto3D.ranking import _SELECTORS, ConformerRanker + + for field, method in _SELECTORS.items(): + assert callable(getattr(ConformerRanker, method, None)), field + + def test_a_selector_in_config_with_nothing_wired_to_it_is_refused(self): + """The point of the check: adding a field to config without wiring a + method here must be impossible to miss. + + This is the exact call ``Auto3D.ranking`` makes at import, with the + field list a developer would have just extended. + """ + from Auto3D.ranking import ( + _SELECTORS, + ConformerRanker, + _verify_selector_registry, + ) + + with pytest.raises(ImportError, match="out of step with"): + _verify_selector_registry( + _SELECTORS, ("k", "window", "percentile"), ConformerRanker + ) + + def test_a_registry_entry_this_module_does_not_implement_is_refused(self): + """A typo in a method name passes the set comparison, then raises + AttributeError from inside ``run`` -- after the whole input has been + read and grouped.""" + from Auto3D.ranking import ConformerRanker, _verify_selector_registry + + with pytest.raises(ImportError, match="not a method of ConformerRanker"): + _verify_selector_registry( + {"k": "top_k", "window": "top_windwo"}, + ("k", "window"), + ConformerRanker, + ) + + def test_the_real_registry_passes_its_own_check(self): + """The inverse: a check that refused everything would satisfy both + tests above and make ``import Auto3D.ranking`` impossible -- so assert + the shipped configuration is accepted.""" + from Auto3D.config import SELECTOR_FIELDS + from Auto3D.ranking import ( + _SELECTORS, + ConformerRanker, + _verify_selector_registry, + ) + + _verify_selector_registry(_SELECTORS, SELECTOR_FIELDS, ConformerRanker) + + def test_a_third_selector_is_dispatched_not_ignored(self, tmp_path, monkeypatch): + """The defect the registry exists to prevent, exercised end to end. + + With the old hand-written ``if self.k: ... elif self.window: ...``, a + selector added to ``SELECTOR_FIELDS`` and wired here was still ignored: + ``run`` consulted only the two names baked into the chain, so a user who + specified the new selector got "Parameter k or window needs to be + specified" for a parameter they had specified. Reverting ``run`` to that + chain must fail this test. + """ + import Auto3D.ranking as ranking + from Auto3D.ranking import ConformerRanker + + calls = [] + + def top_percentile(self, df_group, percentile): + calls.append(percentile) + selected = list(df_group["mols"])[:1] + for mol in selected: + # Every real selector sets this; run() converts it on the way out. + mol.SetProp("E_rel(eV)", "0.0") + return selected + + monkeypatch.setattr( + ranking, "SELECTOR_FIELDS", ("k", "window", "percentile"), raising=True + ) + monkeypatch.setattr( + ranking, + "_SELECTORS", + {"k": "top_k", "window": "top_window", "percentile": "top_percentile"}, + ) + monkeypatch.setattr( + ConformerRanker, "top_percentile", top_percentile, raising=False + ) + + mols = [_create_mol_with_energy("CCO", -10.0, "probe_0_0")] + input_path = str(tmp_path / "in.sdf") + _write_mols_to_sdf(mols, input_path) + + ranker = ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / "out.sdf"), + threshold=0.3, + ) + ranker.percentile = 90.0 + + results = ranker.run() + + assert calls == [90.0], ( + "the third selector was never dispatched to; run() consulted a " + "hard-coded list of selector names instead of the registry" + ) + assert len(results) == 1 + + def test_k_routes_to_top_k_and_window_to_top_window(self, tmp_path): + """The registry must actually be what dispatch consults.""" + from Auto3D.ranking import ConformerRanker + + mols = [_create_mol_with_energy("CCO", -10.0, "probe_0_0")] + input_path = str(tmp_path / "in.sdf") + _write_mols_to_sdf(mols, input_path) + + for field, expected, value in (("k", "top_k", 1), ("window", "top_window", 5.0)): + called = [] + ranker = ConformerRanker( + input_path=input_path, + out_path=str(tmp_path / f"out_{field}.sdf"), + threshold=0.3, + **{field: value}, + ) + for method in ("top_k", "top_window"): + original = getattr(ranker, method) + + def spy(*args, _m=method, _o=original, _log=called, **kwargs): + _log.append(_m) + return _o(*args, **kwargs) + + setattr(ranker, method, spy) + ranker.run() + assert called == [expected], f"{field} dispatched to {called}" diff --git a/tests/test_stereo_postopt.py b/tests/test_stereo_postopt.py index 1bee6be3..8ed2d8df 100644 --- a/tests/test_stereo_postopt.py +++ b/tests/test_stereo_postopt.py @@ -14,7 +14,6 @@ from Auto3D.filtering import filter_unique_optimized from Auto3D.ranking import ConformerRanker -from Auto3D.filtering import filter_unique from Auto3D.utils.stereo_check import ( STEREO_CHANGED_PROP, apply_optimized_coords, @@ -176,12 +175,22 @@ def test_filter_unique_optimized_drops_the_changed_record(self): assert len(result) == 1, f"expected only the preserved record: {len(result)}" assert result[0].GetProp("E_tot") == "-1.0" - def test_filter_unique_drops_the_changed_record(self): - kept = _optimized(-1.0, changed=False) - dropped = _optimized(-2.0, changed=True) - result = filter_unique([dropped, kept], crit=0.3) - assert len(result) == 1, f"expected only the preserved record: {len(result)}" - assert result[0].GetProp("E_tot") == "-1.0" + def test_the_filter_reports_stereochemistry_as_the_drop_reason(self): + """Not just that it dropped, but that it says why. + + Until 4.1.0 the two conformer filters returned a bare list, so + ``ranking`` reported a stereo-changed species as "No structure + converged" -- pointing the reader at the optimizer settings for a + problem in the input's stereo definitions. + """ + from Auto3D.filtering import filter_conformers + + result = filter_conformers( + [_optimized(-2.0, changed=True), _optimized(-1.0, changed=False)], + rmsd_threshold=0.3, + ) + assert result.dropped == {"stereochemistry": 1} + assert result.reasons == ("stereochemistry",) def test_top_k_one_skips_the_changed_lowest_energy_record(self): """k=1 takes a fast path that bypasses the RMSD filters entirely.""" @@ -207,4 +216,3 @@ def test_unmarked_records_still_survive_every_filter(self): """No regression for molecules that never went through the check.""" mols = [_optimized(-1.0, changed=None), _optimized(-2.0, changed=None)] assert len(filter_unique_optimized(mols, rmsd_threshold=0.3)) == 2 - assert len(filter_unique(mols, crit=0.3)) == 2 diff --git a/tests/test_utils_stereochemistry.py b/tests/test_utils_stereochemistry.py index 0a66bfc5..01c9dea9 100644 --- a/tests/test_utils_stereochemistry.py +++ b/tests/test_utils_stereochemistry.py @@ -9,6 +9,7 @@ import pytest from rdkit import Chem +from Auto3D.exceptions import InputValidationError from Auto3D.utils.stereochemistry import ( amend_configuration, amend_configuration_w, @@ -202,6 +203,30 @@ def test_multiple_centers_all_inverted(self): for orig, res in zip(orig_info, result_info): assert orig != res + def test_three_centers_all_inverted(self): + """M60 regression: 3+ centers used to be handled by a loop reading a + variable (key2) set inside the loop's else-branch and read again + after the loop -- correct only by Python's lack of block scoping. + A single-pass rewrite must still invert every center in order.""" + smi = "C[C@H](O)[C@@H](F)[C@H](Cl)Br" + result = create_enantiomer(smi) + assert result == "C[C@@H](O)[C@H](F)[C@@H](Cl)Br" + orig_info = list(get_stereo_info(smi).values()) + result_info = list(get_stereo_info(result).values()) + assert len(result_info) == len(orig_info) == 3 + for orig, res in zip(orig_info, result_info): + assert orig != res + + def test_four_centers_all_inverted(self): + """Same as above, one more center, to rule out an off-by-one at the + boundary between the removed len(keys)==1 special case and the + general multi-key path.""" + smi = "C[C@H](O)[C@@H](F)[C@H](Cl)[C@@H](Br)I" + result = create_enantiomer(smi) + assert result == "C[C@@H](O)[C@H](F)[C@@H](Cl)[C@H](Br)I" + result_info = list(get_stereo_info(result).values()) + assert len(result_info) == 4 + class TestCheckValue: """Tests for the check_value() function.""" @@ -281,6 +306,32 @@ def test_remove_enantiomers_no_stereo(self): os.unlink(inpath) os.unlink(outpath) + def test_remove_enantiomers_tolerates_blank_lines(self, tmp_path): + """M59: switched to iter_smi_records, which skips blank/comment lines. + + The old hand-rolled parser did `vals = line.split()` then indexed + vals[0]/vals[1] unconditionally, so a blank line raised a bare + IndexError and aborted the whole function. + """ + inpath = tmp_path / "in.smi" + inpath.write_text("CCO mol1_1\n\n# a comment\nCCCO mol2_1\n") + outpath = tmp_path / "out.smi" + + result = remove_enantiomers(str(inpath), str(outpath)) + + assert "mol1" in result + assert "mol2" in result + + def test_remove_enantiomers_rejects_missing_id(self, tmp_path): + """A non-blank, non-comment line with no ID must still fail loudly, + just as an InputValidationError rather than a bare IndexError.""" + inpath = tmp_path / "in.smi" + inpath.write_text("CCO\n") + outpath = tmp_path / "out.smi" + + with pytest.raises(InputValidationError): + remove_enantiomers(str(inpath), str(outpath)) + class TestAmendConfiguration: """Tests for the amend_configuration() function.""" @@ -317,6 +368,39 @@ def test_amend_configuration_w_writes_file(self): finally: os.unlink(path) + def test_amend_configuration_tolerates_blank_and_comment_lines(self, tmp_path): + """M59: switched to iter_smi_records. The old parser did + `tuple(line.strip().split())`, which raises "not enough values to + unpack" on a blank line -- aborting the whole function -- and "too + many values to unpack" on a line with a third whitespace column. + """ + path = tmp_path / "in.smi" + path.write_text( + "C[C@H](O)F mol_1\n\n# comment\nC[C@@H](O)F mol_2\n" + ) + + result = amend_configuration(str(path)) + + assert "mol" in result + + def test_amend_configuration_tolerates_extra_column(self, tmp_path): + """A trailing whitespace-separated column beyond SMILES+ID must not + raise, matching every other consumer of this format.""" + path = tmp_path / "in.smi" + path.write_text("C[C@H](O)F mol_1 extra_column\n") + + result = amend_configuration(str(path)) + + assert "mol" in result + + def test_amend_configuration_rejects_missing_id(self, tmp_path): + """A non-blank, non-comment line with no ID must still fail loudly.""" + path = tmp_path / "in.smi" + path.write_text("C[C@H](O)F\n") + + with pytest.raises(InputValidationError): + amend_configuration(str(path)) + class TestIntegration: """Integration tests using real molecular examples.""" diff --git a/tests/test_utils_validation.py b/tests/test_utils_validation.py index a1643907..67d86e15 100644 --- a/tests/test_utils_validation.py +++ b/tests/test_utils_validation.py @@ -6,7 +6,6 @@ from rdkit import Chem from Auto3D.config import Auto3DOptions -from Auto3D.filtering import filter_unique from Auto3D.utils.connectivity import check_connectivity from Auto3D.utils.validation import ( check_input, @@ -150,113 +149,6 @@ def test_check_connectivity_example_sdf(self): assert check_connectivity(mol) is True -class TestFilterUnique: - """Tests for filter_unique function.""" - - def test_filter_unique_removes_unconverged(self): - """Test that filter_unique removes unconverged structures.""" - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - mols = [mol for mol in supp if mol is not None] - - # Set all mols as unconverged - for mol in mols: - mol.SetProp("Converged", "False") - - result = filter_unique(mols) - assert len(result) == 0 - - def test_filter_unique_keeps_converged(self): - """Test that filter_unique keeps converged structures.""" - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - mols = [mol for mol in supp if mol is not None] - - # Set all mols as converged - for mol in mols: - mol.SetProp("Converged", "True") - - result = filter_unique(mols) - # Should keep at least one unique structure - assert len(result) >= 1 - - def test_filter_unique_removes_duplicates(self): - """Test that filter_unique removes similar structures.""" - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - mols = [mol for mol in supp if mol is not None] - - if len(mols) > 0: - # Duplicate the first molecule - mol = mols[0] - mol.SetProp("Converged", "True") - duplicate = Chem.RWMol(mol) - duplicate.SetProp("Converged", "True") - - result = filter_unique([mol, duplicate], crit=0.3) - # Only one should remain - assert len(result) == 1 - - def test_filter_unique_keeps_records_with_no_converged_property(self): - """Absence of the property is not a failed optimization. - - Only ``batchopt`` writes ``Converged``; an ``opt_geometry`` output, an - ORCA/Gaussian export or a hand-built conformer set carries none. - Treating that as "did not converge" deleted every record. The - consequence asserted here is that such a file filters exactly the same - as one whose records all say Converged=True. - """ - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - flagged = [mol for mol in supp if mol is not None] - for mol in flagged: - mol.SetProp("Converged", "True") - expected = len(filter_unique(flagged)) - assert expected >= 1, "test premise: the flagged file must keep something" - - supp = Chem.SDMolSupplier(path_example_sdf, removeHs=False) - unflagged = [mol for mol in supp if mol is not None] - for mol in unflagged: - mol.ClearProp("Converged") - assert not mol.HasProp("Converged") - - result = filter_unique(unflagged) - assert len(result) == expected, ( - f"{len(unflagged)} record(s) with no 'Converged' property kept " - f"{len(result)}, but the same records marked Converged=True keep " - f"{expected}" - ) - - def test_filter_unique_custom_threshold(self): - """A tighter RMSD threshold must keep MORE structures than a looser - one -- not merely "at least as many", which passes even when the two - thresholds produce identical results (as the fixture in - ``path_example_sdf`` does: two molecules of different sizes, so - ``species_key`` alone already keeps both regardless of ``crit``). - - Constructs two conformers of ONE molecule whose RMSD sits strictly - between the two thresholds, so equality cannot pass silently. - """ - mol1 = Chem.AddHs(Chem.MolFromSmiles("CCCCCCCC")) # octane: flexible - mol2 = Chem.Mol(mol1) - from rdkit.Chem import AllChem, rdMolAlign - AllChem.EmbedMolecule(mol1, randomSeed=1) - AllChem.EmbedMolecule(mol2, randomSeed=99) - AllChem.MMFFOptimizeMolecule(mol1) - AllChem.MMFFOptimizeMolecule(mol2) - mol1.SetProp("Converged", "True") - mol2.SetProp("Converged", "True") - - rmsd = rdMolAlign.GetBestRMS(Chem.RemoveHs(mol1), Chem.RemoveHs(mol2)) - assert rmsd > 0.05, "test premise: conformers must be geometrically distinct" - - crit_strict = rmsd / 2 # below the actual RMSD -> kept separate - crit_lenient = rmsd * 2 # above the actual RMSD -> merged - - result_strict = filter_unique([mol1, mol2], crit=crit_strict) - result_lenient = filter_unique([mol1, mol2], crit=crit_lenient) - - assert len(result_strict) == 2, "strict threshold must not merge distinct conformers" - assert len(result_lenient) == 1, "lenient threshold must merge near-identical conformers" - assert len(result_strict) > len(result_lenient) - - def _options(**overrides): """Build the ``Auto3DOptions`` ``check_valid_configuration`` now takes.