chematic_chem/
workflow.rs

1//! High-level, MCP-ready workflows built from the lower-level chematic APIs.
2//!
3//! These functions intentionally keep I/O out of the core logic. They return
4//! structured Rust values that can be serialized by a future MCP/CLI/WASM layer.
5
6use std::fmt;
7
8use chematic_core::{Atom, AtomIdx, BondOrder, Element, Molecule, MoleculeBuilder};
9use chematic_fp::{atom_pair_fp, ecfp4, maccs};
10use chematic_smarts::{
11    AtomPrimitive, AtomQuery, BondPrimitive, BondQuery, McsConfig, QueryMolecule,
12    find_mcs_with_config,
13};
14
15use crate::{
16    bertz_ct, detect_named_functional_groups, exact_mass, formal_charge_sum, fsp3, ghose_passes,
17    hbd_count, heavy_atom_count, identify_functional_groups, labute_asa, logp_and_mr,
18    molecular_weight, murcko_scaffold, num_heteroatoms, num_stereocenters,
19    pains_passes_and_matches, qed_with_bundle, reos_passes, ring_bundle, sa_score_with_bundle,
20    tpsa, wiener_index,
21};
22
23/// Error type for high-level workflow APIs.
24#[derive(Debug, Clone, PartialEq, Eq)]
25pub enum WorkflowError {
26    /// A SMILES input failed to parse.
27    SmilesParse {
28        /// Index of the input in the caller-provided list.
29        index: usize,
30        /// Original SMILES string.
31        smiles: String,
32        /// Parser error message.
33        message: String,
34    },
35    /// A molecule exceeded the configured atom-count limit.
36    TooManyAtoms {
37        /// Index of the input in the caller-provided list.
38        index: usize,
39        /// Original SMILES string.
40        smiles: String,
41        /// Actual atom count.
42        atom_count: usize,
43        /// Configured limit.
44        max_atoms: usize,
45    },
46    /// The caller provided too many molecules for the requested workflow.
47    TooManyMolecules {
48        /// Actual input count.
49        count: usize,
50        /// Configured limit.
51        max_molecules: usize,
52    },
53    /// The workflow requires at least two valid molecules.
54    NeedAtLeastTwoMolecules,
55}
56
57impl fmt::Display for WorkflowError {
58    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
59        match self {
60            Self::SmilesParse {
61                index,
62                smiles,
63                message,
64            } => {
65                write!(
66                    f,
67                    "SMILES parse failed at index {index} ({smiles:?}): {message}"
68                )
69            }
70            Self::TooManyAtoms {
71                index,
72                smiles,
73                atom_count,
74                max_atoms,
75            } => {
76                write!(
77                    f,
78                    "molecule at index {index} ({smiles:?}) has {atom_count} atoms; max is {max_atoms}"
79                )
80            }
81            Self::TooManyMolecules {
82                count,
83                max_molecules,
84            } => {
85                write!(f, "input has {count} molecules; max is {max_molecules}")
86            }
87            Self::NeedAtLeastTwoMolecules => write!(f, "at least two molecules are required"),
88        }
89    }
90}
91
92impl std::error::Error for WorkflowError {}
93
94/// Safety limits for LLM/tool-facing workflows.
95#[derive(Debug, Clone, Copy, PartialEq, Eq)]
96pub struct WorkflowLimits {
97    /// Maximum molecules accepted by multi-molecule workflows.
98    pub max_molecules: usize,
99    /// Maximum atom count accepted for any single molecule.
100    pub max_atoms: usize,
101}
102
103impl Default for WorkflowLimits {
104    fn default() -> Self {
105        Self {
106            max_molecules: 256,
107            max_atoms: 512,
108        }
109    }
110}
111
112/// Scalar descriptors and structural counts for a molecule.
113#[derive(Debug, Clone, PartialEq)]
114#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
115pub struct DescriptorSummary {
116    pub molecular_weight: f64,
117    pub exact_mass: f64,
118    pub tpsa: f64,
119    pub logp: f64,
120    pub molar_refractivity: f64,
121    pub hbd: usize,
122    pub hba: usize,
123    pub rotatable_bonds: usize,
124    pub heavy_atom_count: usize,
125    pub ring_count: usize,
126    pub num_heteroatoms: usize,
127    pub num_stereocenters: usize,
128    pub num_spiro_atoms: usize,
129    pub num_bridgehead_atoms: usize,
130    pub fsp3: f64,
131    pub qed: f64,
132    pub sa_score: f64,
133    pub formal_charge_sum: i32,
134    pub labute_asa: f64,
135    pub bertz_ct: f64,
136    pub wiener_index: f64,
137}
138
139/// Common oral drug-likeness and structural-alert filter outcomes.
140#[derive(Debug, Clone, PartialEq, Eq)]
141#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
142pub struct FilterSummary {
143    pub lipinski_passes: bool,
144    pub veber_passes: bool,
145    pub egan_passes: bool,
146    pub ghose_passes: bool,
147    pub reos_passes: bool,
148    pub pains_passes: bool,
149    pub pains_alerts: Vec<String>,
150}
151
152/// Functional-group label with atom indices.
153#[derive(Debug, Clone, PartialEq, Eq)]
154#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
155pub struct FunctionalGroupSummary {
156    pub name: String,
157    pub atom_indices: Vec<usize>,
158}
159
160/// Named group label with atom indices.
161#[derive(Debug, Clone, PartialEq, Eq)]
162#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
163pub struct NamedGroupSummary {
164    pub name: String,
165    pub atom_indices: Vec<usize>,
166}
167
168/// Complete single-molecule report for tool-facing integrations.
169#[derive(Debug, Clone, PartialEq)]
170#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
171pub struct MoleculeReport {
172    pub input_smiles: String,
173    pub canonical_smiles: String,
174    pub formula: String,
175    pub murcko_scaffold_smiles: Option<String>,
176    pub descriptors: DescriptorSummary,
177    pub filters: FilterSummary,
178    pub functional_groups: Vec<FunctionalGroupSummary>,
179    pub named_groups: Vec<NamedGroupSummary>,
180}
181
182/// Options for single-molecule reports.
183#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
184pub struct ReportOptions {
185    pub limits: WorkflowLimits,
186}
187
188/// Similarity values for a pair of molecules.
189#[derive(Debug, Clone, PartialEq)]
190#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
191pub struct SimilaritySummary {
192    pub ecfp4_tanimoto: f64,
193    pub maccs_tanimoto: f64,
194    pub atom_pair_tanimoto: f64,
195}
196
197/// Pairwise comparison entry.
198#[derive(Debug, Clone, PartialEq)]
199#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
200pub struct PairwiseComparison {
201    pub left_index: usize,
202    pub right_index: usize,
203    pub similarities: SimilaritySummary,
204}
205
206/// Descriptor delta between two molecules (`right - left`).
207#[derive(Debug, Clone, PartialEq)]
208#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
209pub struct DescriptorDelta {
210    pub left_index: usize,
211    pub right_index: usize,
212    pub molecular_weight: f64,
213    pub exact_mass: f64,
214    pub tpsa: f64,
215    pub logp: f64,
216    pub hbd: isize,
217    pub hba: isize,
218    pub rotatable_bonds: isize,
219    pub qed: f64,
220    pub sa_score: f64,
221}
222
223/// Multi-molecule comparison report.
224#[derive(Debug, Clone, PartialEq)]
225#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
226pub struct MoleculeComparison {
227    pub reports: Vec<MoleculeReport>,
228    pub pairwise: Vec<PairwiseComparison>,
229    pub descriptor_deltas: Vec<DescriptorDelta>,
230    pub mcs_smiles: Option<String>,
231}
232
233/// Options for multi-molecule comparison.
234#[derive(Debug, Clone)]
235pub struct CompareOptions {
236    pub limits: WorkflowLimits,
237    pub mcs_config: McsConfig,
238}
239
240impl Default for CompareOptions {
241    fn default() -> Self {
242        Self {
243            limits: WorkflowLimits::default(),
244            mcs_config: McsConfig {
245                timeout_ms: Some(250),
246                ..McsConfig::default()
247            },
248        }
249    }
250}
251
252/// Per-record result for batch screening.
253#[derive(Debug, Clone, PartialEq)]
254#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
255pub struct ScreeningRecord {
256    pub input_index: usize,
257    pub input_smiles: String,
258    pub report: Option<MoleculeReport>,
259    pub error: Option<String>,
260}
261
262/// Batch screening report.
263#[derive(Debug, Clone, PartialEq)]
264#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
265pub struct ScreeningReport {
266    pub records: Vec<ScreeningRecord>,
267    /// Original input indices selected by MaxMin diversity picking.
268    pub maxmin_picks: Vec<usize>,
269    /// Original input indices grouped by Butina clustering.
270    pub butina_clusters: Vec<Vec<usize>>,
271}
272
273/// Options for batch screening.
274#[derive(Debug, Clone, Copy, PartialEq)]
275pub struct ScreenOptions {
276    pub limits: WorkflowLimits,
277    pub maxmin_pick_count: usize,
278    pub butina_cutoff: f64,
279}
280
281impl Default for ScreenOptions {
282    fn default() -> Self {
283        Self {
284            limits: WorkflowLimits::default(),
285            maxmin_pick_count: 0,
286            butina_cutoff: 0.65,
287        }
288    }
289}
290
291/// Build a complete report from a SMILES string using default options.
292pub fn molecule_report(smiles: &str) -> Result<MoleculeReport, WorkflowError> {
293    molecule_report_with_options(smiles, &ReportOptions::default())
294}
295
296/// Build a complete report from a SMILES string.
297pub fn molecule_report_with_options(
298    smiles: &str,
299    options: &ReportOptions,
300) -> Result<MoleculeReport, WorkflowError> {
301    let mol = parse_checked(0, smiles, options.limits.max_atoms)?;
302    Ok(report_for_molecule(smiles, &mol))
303}
304
305/// Compare two or more SMILES strings using default options.
306pub fn compare_molecules(smiles: &[&str]) -> Result<MoleculeComparison, WorkflowError> {
307    compare_molecules_with_options(smiles, &CompareOptions::default())
308}
309
310/// Compare two or more SMILES strings.
311pub fn compare_molecules_with_options(
312    smiles: &[&str],
313    options: &CompareOptions,
314) -> Result<MoleculeComparison, WorkflowError> {
315    if smiles.len() < 2 {
316        return Err(WorkflowError::NeedAtLeastTwoMolecules);
317    }
318    if smiles.len() > options.limits.max_molecules {
319        return Err(WorkflowError::TooManyMolecules {
320            count: smiles.len(),
321            max_molecules: options.limits.max_molecules,
322        });
323    }
324
325    let mols = parse_many_checked(smiles, options.limits.max_atoms)?;
326    let reports = smiles
327        .iter()
328        .zip(mols.iter())
329        .map(|(s, mol)| report_for_molecule(s, mol))
330        .collect::<Vec<_>>();
331
332    let mut pairwise = Vec::new();
333    let mut descriptor_deltas = Vec::new();
334    for i in 0..mols.len() {
335        for j in (i + 1)..mols.len() {
336            pairwise.push(PairwiseComparison {
337                left_index: i,
338                right_index: j,
339                similarities: similarity_summary(&mols[i], &mols[j]),
340            });
341            descriptor_deltas.push(descriptor_delta(i, j, &reports[i], &reports[j]));
342        }
343    }
344
345    let mol_refs = mols.iter().collect::<Vec<_>>();
346    let qmol = find_mcs_with_config(&mol_refs, &options.mcs_config);
347    let mcs_smiles = query_molecule_to_smiles(&qmol);
348
349    Ok(MoleculeComparison {
350        reports,
351        pairwise,
352        descriptor_deltas,
353        mcs_smiles,
354    })
355}
356
357/// Screen a SMILES list using default options. Invalid records are retained
358/// as per-record errors instead of failing the whole batch.
359pub fn screen_smiles(smiles: &[&str]) -> ScreeningReport {
360    screen_smiles_with_options(smiles, &ScreenOptions::default())
361}
362
363/// Screen a SMILES list. Invalid records are retained as per-record errors
364/// instead of failing the whole batch.
365pub fn screen_smiles_with_options(smiles: &[&str], options: &ScreenOptions) -> ScreeningReport {
366    let mut records = Vec::with_capacity(smiles.len());
367    let mut valid_mols = Vec::new();
368    let mut valid_original_indices = Vec::new();
369
370    for (idx, smi) in smiles.iter().enumerate() {
371        if idx >= options.limits.max_molecules {
372            records.push(ScreeningRecord {
373                input_index: idx,
374                input_smiles: (*smi).to_string(),
375                report: None,
376                error: Some(format!(
377                    "input index exceeds max_molecules={}",
378                    options.limits.max_molecules
379                )),
380            });
381            continue;
382        }
383
384        match parse_checked(idx, smi, options.limits.max_atoms) {
385            Ok(mol) => {
386                let report = report_for_molecule(smi, &mol);
387                valid_original_indices.push(idx);
388                valid_mols.push(mol);
389                records.push(ScreeningRecord {
390                    input_index: idx,
391                    input_smiles: (*smi).to_string(),
392                    report: Some(report),
393                    error: None,
394                });
395            }
396            Err(err) => records.push(ScreeningRecord {
397                input_index: idx,
398                input_smiles: (*smi).to_string(),
399                report: None,
400                error: Some(err.to_string()),
401            }),
402        }
403    }
404
405    let maxmin_picks = if options.maxmin_pick_count == 0 || valid_mols.is_empty() {
406        Vec::new()
407    } else {
408        crate::maxmin_picks(&valid_mols, options.maxmin_pick_count, |a, b| {
409            ecfp4(a).tanimoto(&ecfp4(b))
410        })
411        .into_iter()
412        .map(|valid_idx| valid_original_indices[valid_idx])
413        .collect()
414    };
415
416    let butina_clusters = if valid_mols.is_empty() {
417        Vec::new()
418    } else {
419        crate::butina_cluster(&valid_mols, options.butina_cutoff, |a, b| {
420            ecfp4(a).tanimoto(&ecfp4(b))
421        })
422        .into_iter()
423        .map(|cluster| {
424            cluster
425                .into_iter()
426                .map(|valid_idx| valid_original_indices[valid_idx])
427                .collect()
428        })
429        .collect()
430    };
431
432    ScreeningReport {
433        records,
434        maxmin_picks,
435        butina_clusters,
436    }
437}
438
439fn parse_many_checked(smiles: &[&str], max_atoms: usize) -> Result<Vec<Molecule>, WorkflowError> {
440    smiles
441        .iter()
442        .enumerate()
443        .map(|(idx, smi)| parse_checked(idx, smi, max_atoms))
444        .collect()
445}
446
447fn parse_checked(index: usize, smiles: &str, max_atoms: usize) -> Result<Molecule, WorkflowError> {
448    let mol = chematic_smiles::parse(smiles).map_err(|e| WorkflowError::SmilesParse {
449        index,
450        smiles: smiles.to_string(),
451        message: e.to_string(),
452    })?;
453    if mol.atom_count() > max_atoms {
454        return Err(WorkflowError::TooManyAtoms {
455            index,
456            smiles: smiles.to_string(),
457            atom_count: mol.atom_count(),
458            max_atoms,
459        });
460    }
461    Ok(mol)
462}
463
464fn report_for_molecule(input_smiles: &str, mol: &Molecule) -> MoleculeReport {
465    let scaffold = murcko_scaffold(mol);
466    let murcko_scaffold_smiles = if scaffold.atom_count() == 0 {
467        None
468    } else {
469        Some(chematic_smiles::canonical_smiles(&scaffold))
470    };
471
472    // Compute all ring-derived descriptors with a single find_sssr call.
473    let rb = ring_bundle(mol);
474    let mw = molecular_weight(mol);
475    // Compute LogP and MR together to share the 117-pattern Crippen SMARTS pass.
476    let (logp, mr) = logp_and_mr(mol);
477    let tpsa_val = tpsa(mol);
478    // Single explicit-H + SSSR + 480-pattern pass for both PAINS flag and alert names.
479    let (pains_ok, pains_alert_names) = pains_passes_and_matches(mol);
480
481    MoleculeReport {
482        input_smiles: input_smiles.to_string(),
483        canonical_smiles: chematic_smiles::canonical_smiles(mol),
484        formula: mol.total_formula(),
485        murcko_scaffold_smiles,
486        descriptors: DescriptorSummary {
487            molecular_weight: mw,
488            exact_mass: exact_mass(mol),
489            tpsa: tpsa_val,
490            logp,
491            molar_refractivity: mr,
492            hbd: hbd_count(mol),
493            hba: rb.hba_count,
494            rotatable_bonds: rb.rotatable_bond_count,
495            heavy_atom_count: heavy_atom_count(mol),
496            ring_count: rb.ring_count,
497            num_heteroatoms: num_heteroatoms(mol),
498            num_stereocenters: num_stereocenters(mol),
499            num_spiro_atoms: rb.num_spiro_atoms,
500            num_bridgehead_atoms: rb.num_bridgehead_atoms,
501            fsp3: fsp3(mol),
502            qed: qed_with_bundle(mol, &rb),
503            sa_score: sa_score_with_bundle(mol, &rb),
504            formal_charge_sum: formal_charge_sum(mol),
505            labute_asa: labute_asa(mol),
506            bertz_ct: bertz_ct(mol),
507            wiener_index: wiener_index(mol),
508        },
509        filters: FilterSummary {
510            // Inline only the filters that call rotatable_bond_count or hba_count (→ find_sssr).
511            // Other filters (egan, ghose, reos) don't use ring data — call them directly.
512            lipinski_passes: mw <= 500.0
513                && hbd_count(mol) <= 5
514                && rb.hba_count <= 10
515                && logp <= 5.0,
516            veber_passes: rb.rotatable_bond_count <= 10 && tpsa_val <= 140.0,
517            egan_passes: tpsa_val <= 131.6 && logp <= 5.88,
518            ghose_passes: ghose_passes(mol),
519            reos_passes: reos_passes(mol),
520            pains_passes: pains_ok,
521            pains_alerts: pains_alert_names.into_iter().map(str::to_string).collect(),
522        },
523        functional_groups: identify_functional_groups(mol)
524            .into_iter()
525            .map(|fg| FunctionalGroupSummary {
526                name: fg.atom_types,
527                atom_indices: fg.atom_indices,
528            })
529            .collect(),
530        named_groups: detect_named_functional_groups(mol)
531            .into_iter()
532            .map(|ng| NamedGroupSummary {
533                name: ng.name.to_string(),
534                atom_indices: ng.atoms.into_iter().map(|idx| idx.0 as usize).collect(),
535            })
536            .collect(),
537    }
538}
539
540fn similarity_summary(a: &Molecule, b: &Molecule) -> SimilaritySummary {
541    SimilaritySummary {
542        ecfp4_tanimoto: ecfp4(a).tanimoto(&ecfp4(b)),
543        maccs_tanimoto: maccs(a).tanimoto(&maccs(b)),
544        atom_pair_tanimoto: atom_pair_fp(a).tanimoto(&atom_pair_fp(b)),
545    }
546}
547
548fn descriptor_delta(
549    left_index: usize,
550    right_index: usize,
551    left: &MoleculeReport,
552    right: &MoleculeReport,
553) -> DescriptorDelta {
554    DescriptorDelta {
555        left_index,
556        right_index,
557        molecular_weight: right.descriptors.molecular_weight - left.descriptors.molecular_weight,
558        exact_mass: right.descriptors.exact_mass - left.descriptors.exact_mass,
559        tpsa: right.descriptors.tpsa - left.descriptors.tpsa,
560        logp: right.descriptors.logp - left.descriptors.logp,
561        hbd: right.descriptors.hbd as isize - left.descriptors.hbd as isize,
562        hba: right.descriptors.hba as isize - left.descriptors.hba as isize,
563        rotatable_bonds: right.descriptors.rotatable_bonds as isize
564            - left.descriptors.rotatable_bonds as isize,
565        qed: right.descriptors.qed - left.descriptors.qed,
566        sa_score: right.descriptors.sa_score - left.descriptors.sa_score,
567    }
568}
569
570fn query_molecule_to_smiles(qmol: &QueryMolecule) -> Option<String> {
571    if qmol.atoms.is_empty() {
572        return None;
573    }
574
575    let mut aromatic_atoms = vec![false; qmol.atoms.len()];
576    for (atom_idx, neighbors) in qmol.adj.iter().enumerate() {
577        for (bond_idx, neighbor_idx) in neighbors {
578            if matches!(
579                qmol.bonds[*bond_idx].query,
580                BondQuery::Primitive(BondPrimitive::Aromatic)
581            ) {
582                aromatic_atoms[atom_idx] = true;
583                aromatic_atoms[*neighbor_idx] = true;
584            }
585        }
586    }
587
588    let mut builder = MoleculeBuilder::new();
589    for (idx, qa) in qmol.atoms.iter().enumerate() {
590        let elem = match &qa.query {
591            AtomQuery::Primitive(AtomPrimitive::AtomicNum(n)) => {
592                Element::from_atomic_number(*n).unwrap_or(Element::C)
593            }
594            _ => Element::C,
595        };
596        let mut atom = Atom::new(elem);
597        atom.aromatic = aromatic_atoms[idx];
598        builder.add_atom(atom);
599    }
600
601    for (atom_idx, neighbors) in qmol.adj.iter().enumerate() {
602        for (bond_idx, neighbor_idx) in neighbors {
603            if atom_idx < *neighbor_idx {
604                let order = match &qmol.bonds[*bond_idx].query {
605                    BondQuery::Primitive(BondPrimitive::Double) => BondOrder::Double,
606                    BondQuery::Primitive(BondPrimitive::Triple) => BondOrder::Triple,
607                    BondQuery::Primitive(BondPrimitive::Aromatic) => BondOrder::Aromatic,
608                    _ => BondOrder::Single,
609                };
610                let _ = builder.add_bond(
611                    AtomIdx(atom_idx as u32),
612                    AtomIdx(*neighbor_idx as u32),
613                    order,
614                );
615            }
616        }
617    }
618
619    Some(chematic_smiles::canonical_smiles(&builder.build()))
620}
621
622#[cfg(test)]
623mod tests {
624    use super::*;
625
626    #[test]
627    fn molecule_report_aspirin_has_core_fields() {
628        let report = molecule_report("CC(=O)Oc1ccccc1C(=O)O").unwrap();
629
630        assert_eq!(report.formula, "C9H8O4");
631        assert_eq!(report.descriptors.heavy_atom_count, 13);
632        assert!(report.descriptors.molecular_weight > 180.0);
633        assert!(report.descriptors.tpsa > 60.0);
634        assert!(report.filters.lipinski_passes);
635        assert_eq!(report.murcko_scaffold_smiles.as_deref(), Some("c1ccccc1"));
636    }
637
638    #[test]
639    fn molecule_report_invalid_smiles_returns_structured_error() {
640        let err = molecule_report("C1CC").unwrap_err();
641        assert!(matches!(err, WorkflowError::SmilesParse { index: 0, .. }));
642    }
643
644    #[test]
645    fn compare_molecules_returns_pairwise_and_mcs() {
646        let comparison = compare_molecules(&["c1ccccc1", "Cc1ccccc1"]).unwrap();
647
648        assert_eq!(comparison.reports.len(), 2);
649        assert_eq!(comparison.pairwise.len(), 1);
650        assert_eq!(comparison.descriptor_deltas.len(), 1);
651        assert!(comparison.pairwise[0].similarities.ecfp4_tanimoto > 0.0);
652        assert_eq!(comparison.mcs_smiles.as_deref(), Some("c1ccccc1"));
653    }
654
655    #[test]
656    fn compare_molecules_needs_two_inputs() {
657        let err = compare_molecules(&["CC"]).unwrap_err();
658        assert_eq!(err, WorkflowError::NeedAtLeastTwoMolecules);
659    }
660
661    #[test]
662    fn screen_smiles_keeps_invalid_records_and_original_indices() {
663        let options = ScreenOptions {
664            maxmin_pick_count: 2,
665            ..ScreenOptions::default()
666        };
667        let report = screen_smiles_with_options(&["CC", "C1CC", "c1ccccc1"], &options);
668
669        assert_eq!(report.records.len(), 3);
670        assert!(report.records[0].report.is_some());
671        assert!(report.records[1].report.is_none());
672        assert!(report.records[1].error.is_some());
673        assert!(report.records[2].report.is_some());
674        assert!(report.maxmin_picks.iter().all(|idx| *idx == 0 || *idx == 2));
675        assert!(
676            report
677                .butina_clusters
678                .iter()
679                .flatten()
680                .all(|idx| *idx == 0 || *idx == 2)
681        );
682    }
683
684    #[test]
685    fn limits_reject_large_molecule() {
686        let options = ReportOptions {
687            limits: WorkflowLimits {
688                max_atoms: 1,
689                ..WorkflowLimits::default()
690            },
691        };
692        let err = molecule_report_with_options("CC", &options).unwrap_err();
693        assert!(matches!(
694            err,
695            WorkflowError::TooManyAtoms { atom_count: 2, .. }
696        ));
697    }
698
699    #[test]
700    fn molecule_report_sulfur_compound() {
701        let report = molecule_report("c1ccccc1S(=O)(=O)N").unwrap();
702        assert_eq!(report.descriptors.num_heteroatoms, 4); // S, O, O, N
703        assert!(report.descriptors.molecular_weight > 150.0);
704    }
705
706    #[test]
707    fn molecule_report_halogenated() {
708        let report = molecule_report("ClC(Br)(F)I").unwrap();
709        assert_eq!(report.descriptors.heavy_atom_count, 5); // C, Cl, Br, F, I
710        assert!(report.descriptors.num_heteroatoms == 4); // Cl, Br, F, I (not C)
711    }
712
713    #[test]
714    fn molecule_report_complex_aromatic() {
715        // Quinoline (bicyclic aromatic)
716        let report = molecule_report("c1ccc2ncccc2c1").unwrap();
717        assert_eq!(report.descriptors.heavy_atom_count, 10); // 9 C + 1 N
718        assert!(report.descriptors.ring_count >= 2);
719        assert!(report.filters.lipinski_passes);
720    }
721
722    #[test]
723    fn molecule_report_large_valid_molecule() {
724        // Taxol-like structure (simplified, ~50 atoms)
725        let report = molecule_report("CC(=O)Oc1ccccc1C(=O)N[C@@H]1C[C@H]2CC(C)(C)[C@@H](O)C[C@]2(OC(=O)C(C)C)[C@]1(O)C(=O)c1ccccc1").unwrap();
726        assert!(report.descriptors.heavy_atom_count > 30);
727        assert!(report.descriptors.num_stereocenters > 0);
728    }
729
730    #[test]
731    fn molecule_report_charge_species() {
732        let report = molecule_report("[Na+].[Cl-]").unwrap();
733        assert_eq!(report.formula, "ClNa");
734        assert_eq!(report.descriptors.formal_charge_sum, 0); // +1 + -1 = 0
735    }
736
737    #[test]
738    fn compare_molecules_identical_molecules() {
739        let comparison = compare_molecules(&["c1ccccc1", "c1ccccc1"]).unwrap();
740        let sim = comparison.pairwise[0].similarities.ecfp4_tanimoto;
741        assert!(
742            (sim - 1.0).abs() < 1e-6,
743            "identical molecules should have ~100% similarity"
744        );
745    }
746
747    #[test]
748    fn screen_smiles_empty_batch() {
749        let report = screen_smiles(&[]);
750        assert_eq!(report.records.len(), 0);
751        assert_eq!(report.maxmin_picks.len(), 0);
752    }
753
754    #[test]
755    fn molecule_report_aromatic_nitrogen() {
756        // Pyridine
757        let report = molecule_report("c1ccncc1").unwrap();
758        assert_eq!(report.descriptors.ring_count, 1);
759        assert!(report.descriptors.hba > 0); // N is an acceptor
760    }
761}
chematic_chem/workflow.rs

chematic_chem/
workflow.rs
