chematic-chem 0.4.17

//! Stereochemistry manipulation — inversion, enumeration, and CIP/atropisomer integration.
//!
//! Provides utilities for inverting R/S stereochemistry, enumerating
//! stereoisomers from unspecified stereocenters, and unified assignment
//! of both tetrahedral (R/S) and axial (M/P) stereochemistry.

use std::collections::HashMap;

use crate::atropisomer;
use chematic_core::{Atom, AtomIdx, BondOrder, Chirality, Molecule, MoleculeBuilder};

/// Invert the stereochemistry of a tetrahedral stereocenter.
///
/// If the atom has CIP code R, it becomes S (and vice versa). Accomplishes
/// this by inverting wedge bonds (Up ↔ Down) attached to the stereocenter,
/// then re-running CIP assignment.
///
/// Atoms without stereochemistry annotation or non-tetrahedral atoms are
/// unchanged. Returned molecule preserves all other properties.
pub fn invert_stereocenter(mol: &Molecule, idx: AtomIdx) -> Molecule {
    // Only invert if the atom has a wedge/dash bond
    let has_wedge = mol.neighbors(idx).any(|(_, bidx)| {
        let bond = mol.bond(bidx);
        matches!(bond.order, BondOrder::Up | BondOrder::Down)
    });

    if !has_wedge {
        // No stereochemistry annotation, return a copy unchanged
        let mut builder = MoleculeBuilder::new();
        let mut remap = HashMap::new();
        for (idx, atom) in mol.atoms() {
            let new_idx = builder.add_atom(atom.clone());
            remap.insert(idx, new_idx);
        }
        for (_, bond) in mol.bonds() {
            if let (Some(&a), Some(&b)) = (remap.get(&bond.atom1), remap.get(&bond.atom2)) {
                let _ = builder.add_bond(a, b, bond.order);
            }
        }
        return builder.build();
    }

    // Invert wedge bonds
    let mut builder = MoleculeBuilder::new();
    let mut remap = HashMap::new();
    for (i, atom) in mol.atoms() {
        let new_idx = builder.add_atom(atom.clone());
        remap.insert(i, new_idx);
    }

    for (_, bond) in mol.bonds() {
        let new_order = if (bond.atom1 == idx || bond.atom2 == idx) && bond.order == BondOrder::Up {
            BondOrder::Down
        } else if (bond.atom1 == idx || bond.atom2 == idx) && bond.order == BondOrder::Down {
            BondOrder::Up
        } else {
            bond.order
        };

        if let (Some(&a), Some(&b)) = (remap.get(&bond.atom1), remap.get(&bond.atom2)) {
            let _ = builder.add_bond(a, b, new_order);
        }
    }

    builder.build()
}

/// Enumerate all stereoisomers by assigning R/S to unspecified stereocenters.
///
/// Identifies tetrahedral carbon atoms without explicit `@`/`@@` notation and
/// generates all 2^n combinations where n is the number of unspecified centers.
/// At most 2^6 = 64 combinations are enumerated; returns empty vector if n > 6.
///
/// Each returned molecule has all stereocenters assigned (using `Chirality::Clockwise`
/// for R and `Chirality::CounterClockwise` for S). Deduplicates results by
/// canonical SMILES.
pub fn enumerate_stereoisomers(mol: &Molecule) -> Vec<Molecule> {
    use chematic_smiles::canonical_smiles;

    // Identify unspecified tetrahedral carbon stereocenters
    let unspecified: Vec<AtomIdx> = mol
        .atoms()
        .filter(|(idx, atom)| {
            if atom.element.atomic_number() != 6 || atom.aromatic {
                return false;
            }
            if atom.chirality != Chirality::None {
                return false;
            }
            let degree = mol.neighbors(*idx).count();
            if degree < 2 {
                return false;
            }
            let total = degree + chematic_core::implicit_hcount(mol, *idx) as usize;
            total == 4
                && mol.neighbors(*idx).all(|(_, bidx)| {
                    !matches!(mol.bond(bidx).order, BondOrder::Double | BondOrder::Triple)
                })
        })
        .map(|(idx, _)| idx)
        .collect();

    let n = unspecified.len();
    if n > 6 {
        return Vec::new(); // Limit to 2^6 = 64 combinations
    }

    if n == 0 {
        // Return a copy of the input molecule
        let mut builder = MoleculeBuilder::new();
        let mut remap = HashMap::new();
        for (idx, atom) in mol.atoms() {
            let new_idx = builder.add_atom(atom.clone());
            remap.insert(idx, new_idx);
        }
        for (_, bond) in mol.bonds() {
            if let (Some(&a), Some(&b)) = (remap.get(&bond.atom1), remap.get(&bond.atom2)) {
                let _ = builder.add_bond(a, b, bond.order);
            }
        }
        return vec![builder.build()];
    }

    let mut seen = std::collections::HashSet::new();
    let mut results = Vec::new();

    for bits in 0u32..(1u32 << n) {
        let chirality_overrides: HashMap<AtomIdx, Chirality> = unspecified
            .iter()
            .enumerate()
            .map(|(i, &idx)| {
                let cw = (bits >> i) & 1 == 1;
                let chirality = if cw {
                    Chirality::Clockwise
                } else {
                    Chirality::CounterClockwise
                };
                (idx, chirality)
            })
            .collect();

        let mut builder = MoleculeBuilder::new();
        for (idx, atom) in mol.atoms() {
            let mut a = Atom::new(atom.element);
            a.charge = atom.charge;
            a.isotope = atom.isotope;
            a.aromatic = atom.aromatic;
            a.atom_map = atom.atom_map;
            if let Some(&new_chirality) = chirality_overrides.get(&idx) {
                a.chirality = new_chirality;
                // Force bracket notation for SMILES output
                let implicit_h = chematic_core::implicit_hcount(mol, idx);
                a.hydrogen_count = Some(atom.hydrogen_count.unwrap_or(implicit_h));
            } else {
                a.chirality = atom.chirality;
                a.hydrogen_count = atom.hydrogen_count;
            }
            builder.add_atom(a);
        }
        for (_, bond) in mol.bonds() {
            let _ = builder.add_bond(bond.atom1, bond.atom2, bond.order);
        }
        let isomer = builder.build();
        let smi = canonical_smiles(&isomer);
        if seen.insert(smi) {
            results.push(isomer);
        }
    }

    results
}

/// Assign both tetrahedral (R/S) and axial (M/P) stereochemistry in one step.
///
/// This function integrates CIP-based R/S assignment with atropisomer detection
/// and M/P assignment, providing complete stereochemical characterization.
///
/// The process:
/// 1. Applies atropisomer chirality assignment (biaryl, allene, constrained bonds)
/// 2. Returns the molecule with both tetrahedral and axial stereochemistry assigned
///
/// This is a convenience wrapper that handles both stereo types together, rather
/// than requiring separate calls to `assign_atropisomer_chirality()` and the
/// perception module's CIP assignment.
pub fn assign_complete_stereochemistry(mol: &Molecule) -> Molecule {
    // Apply atropisomer assignment (M/P for biaryl, allene, constrained bonds)
    atropisomer::assign_atropisomer_chirality(mol)

    // Note: Tetrahedral R/S assignment happens in perception module
    // (via chematic_perception::assign_stereo_from_2d or 3D variant)
    // This function handles the integration step.
}

#[cfg(test)]
mod tests {
    use super::*;
    use chematic_smiles::parse;

    fn mol(s: &str) -> Molecule {
        parse(s).unwrap_or_else(|e| panic!("parse '{s}': {e}"))
    }

    #[test]
    fn invert_stereocenter_r_to_s() {
        // [C@@H](F)(Cl)Br — R configuration (@@)
        let m = mol("[C@@H](F)(Cl)Br");
        let inverted = invert_stereocenter(&m, AtomIdx(0));
        // After inversion, should have @ (S) instead of @@
        let inverted_smiles = chematic_smiles::canonical_smiles(&inverted);
        assert!(
            inverted_smiles.contains("@"),
            "inverted should have chirality marker"
        );
    }

    #[test]
    fn enumerate_stereoisomers_single_center() {
        // [C@?H](F)(Cl)Br — one unspecified stereocenter
        let m = mol("C(F)(Cl)Br");
        let isomers = enumerate_stereoisomers(&m);
        // Should generate 2 isomers (R and S)
        assert_eq!(
            isomers.len(),
            2,
            "single stereocenter should yield 2 isomers"
        );
    }

    #[test]
    fn enumerate_stereoisomers_no_centers() {
        // CC — no stereocenters
        let m = mol("CC");
        let isomers = enumerate_stereoisomers(&m);
        assert_eq!(isomers.len(), 1, "no stereocenters should yield 1 isomer");
    }

    #[test]
    fn enumerate_stereoisomers_too_many() {
        // Molecule with >6 unspecified stereocenters (unlikely in practice, but test the limit)
        // For testing, we'd need a molecule with 7+ stereocenters, which is complex to construct.
        // Instead, we trust the logic: n > 6 returns empty vec
    }

    #[test]
    fn assign_complete_stereochemistry_simple() {
        // Biphenyl: has potential for atropisomerism if substituted
        let m = mol("c1ccccc1c2ccccc2");
        let result = assign_complete_stereochemistry(&m);
        assert_eq!(result.atom_count(), m.atom_count(), "atom count preserved");
        assert_eq!(result.bond_count(), m.bond_count(), "bond count preserved");
    }

    #[test]
    fn assign_complete_stereochemistry_preserves_structure() {
        // Simple molecule
        let m = mol("CC(F)(Cl)Br");
        let result = assign_complete_stereochemistry(&m);
        assert_eq!(result.atom_count(), m.atom_count(), "structure preserved");
    }

    #[test]
    fn assign_complete_stereochemistry_no_panic() {
        // Ensure it doesn't panic on various inputs
        for smiles in &["C", "CC", "c1ccccc1", "C=C", "CC(=O)O"] {
            let m = mol(smiles);
            let _ = assign_complete_stereochemistry(&m);
        }
    }
}