chematic-chem 0.1.32

//! Molecular descriptor functions for drug-likeness and physical property estimation.
//!
//! All functions accept a `&Molecule` reference.  Molecules with aromatic bonds
//! (SMILES lowercase notation) are kekulized internally where hydrogen counts
//! are required; the caller's molecule is never mutated.

use std::collections::HashSet;

use chematic_core::{AtomIdx, BondOrder, BondIdx, Element, Molecule, implicit_hcount};
use chematic_perception::find_sssr;

/// True if `idx` has a double bond to any neighbor whose atomic number equals `target_an`.
fn has_double_bond_to(mol: &Molecule, idx: AtomIdx, target_an: u8) -> bool {
    mol.neighbors(idx).any(|(nb, bidx)| {
        mol.bond(bidx).order == BondOrder::Double
            && mol.atom(nb).element.atomic_number() == target_an
    })
}

/// Count double bonds from `idx` to neighbors whose atomic number equals `target_an`.
fn count_double_bonds_to(mol: &Molecule, idx: AtomIdx, target_an: u8) -> usize {
    mol.neighbors(idx)
        .filter(|&(nb, bidx)| {
            mol.bond(bidx).order == BondOrder::Double
                && mol.atom(nb).element.atomic_number() == target_an
        })
        .count()
}

/// True if any neighbor of `idx` is aromatic.
fn has_aromatic_neighbor(mol: &Molecule, idx: AtomIdx) -> bool {
    mol.neighbors(idx).any(|(nb, _)| mol.atom(nb).aromatic)
}

/// True if any neighbor of `idx` is an aromatic carbon.
fn has_aromatic_carbon_neighbor(mol: &Molecule, idx: AtomIdx) -> bool {
    mol.neighbors(idx)
        .any(|(nb, _)| mol.atom(nb).aromatic && mol.atom(nb).element.atomic_number() == 6)
}

/// Average atomic mass table.
/// Falls back to `atomic_number as f64` for unlisted elements.
fn avg_mass(element: Element) -> f64 {
    match element.atomic_number() {
        1  => 1.008,   // H
        2  => 4.003,   // He
        3  => 6.941,   // Li
        4  => 9.012,   // Be
        5  => 10.811,  // B
        6  => 12.011,  // C
        7  => 14.007,  // N
        8  => 15.999,  // O
        9  => 18.998,  // F
        10 => 20.180,  // Ne
        11 => 22.990,  // Na
        12 => 24.305,  // Mg
        13 => 26.982,  // Al
        14 => 28.086,  // Si
        15 => 30.974,  // P
        16 => 32.065,  // S
        17 => 35.453,  // Cl
        18 => 39.948,  // Ar
        19 => 39.098,  // K
        20 => 40.078,  // Ca
        33 => 74.922,  // As
        34 => 78.971,  // Se
        35 => 79.904,  // Br
        53 => 126.904, // I
        n  => n as f64,
    }
}

/// Monoisotopic (most-abundant isotope) mass table.
/// Falls back to `atomic_number as f64` for unlisted elements.
fn mono_mass(element: Element) -> f64 {
    match element.atomic_number() {
        1  => 1.00783,   // H  (1H)
        6  => 12.0000,   // C  (12C)
        7  => 14.0031,   // N  (14N)
        8  => 15.9949,   // O  (16O)
        9  => 18.9984,   // F  (19F)
        14 => 27.9769,   // Si (28Si)
        15 => 30.9738,   // P  (31P)
        16 => 31.9721,   // S  (32S)
        17 => 34.9689,   // Cl (35Cl)
        35 => 78.9183,   // Br (79Br)
        34 => 79.9165,   // Se (80Se)
        53 => 126.9045,  // I  (127I)
        n  => n as f64,
    }
}

// ---------------------------------------------------------------------------
// 1. Molecular weight
// ---------------------------------------------------------------------------

/// Compute the average molecular weight (Da).
///
/// Sums the average atomic mass of all heavy atoms plus each atom's implicit
/// hydrogen contribution (1.008 Da per H).
pub fn molecular_weight(mol: &Molecule) -> f64 {
    let mut mw = 0.0f64;
    for (idx, atom) in mol.atoms() {
        mw += avg_mass(atom.element);
        let h = implicit_hcount(mol, idx);
        mw += h as f64 * 1.008;
    }
    mw
}

// ---------------------------------------------------------------------------
// 2. Exact mass (monoisotopic)
// ---------------------------------------------------------------------------

/// Compute the monoisotopic (exact) mass (Da).
///
/// Uses the most-abundant isotope for each element, or the atom's explicit
/// isotope label (as an integer approximation) when set.
/// Implicit hydrogens use the ¹H monoisotopic mass (1.00783).
pub fn exact_mass(mol: &Molecule) -> f64 {
    let mut mass = 0.0f64;
    for (idx, atom) in mol.atoms() {
        let m = match atom.isotope {
            Some(iso) => iso as f64,
            None => mono_mass(atom.element),
        };
        mass += m;
        let h = implicit_hcount(mol, idx);
        mass += h as f64 * 1.00783;
    }
    mass
}

// ---------------------------------------------------------------------------
// 3. Heavy atom count
// ---------------------------------------------------------------------------

/// Count non-hydrogen heavy atoms.
///
/// Hydrogen atoms are normally implicit in chematic, but some molecules may
/// carry explicit H atoms in the graph (e.g. from bracket notation `[H]`).
/// Those are excluded from the heavy-atom count.
pub fn heavy_atom_count(mol: &Molecule) -> usize {
    mol.atoms()
        .filter(|(_, atom)| atom.element != Element::H)
        .count()
}

// ---------------------------------------------------------------------------
// 4. Hydrogen bond donor count
// ---------------------------------------------------------------------------

/// Count hydrogen bond donors (N-H or O-H groups).
///
/// Each heavy atom with element N or O that has at least one attached H
/// counts as one donor (not per H — donors are counted per heavy atom).
pub fn hbd_count(mol: &Molecule) -> usize {
    mol.atoms()
        .filter(|(idx, atom)| {
            let an = atom.element.atomic_number();
            (an == 7 || an == 8) && implicit_hcount(mol, *idx) > 0
        })
        .count()
}

// ---------------------------------------------------------------------------
// 5. Hydrogen bond acceptor count (Ertl / RDKit-aligned definition)
// ---------------------------------------------------------------------------

/// Count hydrogen bond acceptors using the Ertl (2000) definition as implemented
/// by RDKit's `rdMolDescriptors.CalcNumHBA`.
///
/// Counts N, O, and divalent S atoms, with the following exclusions:
/// - Aromatic N with H (pyrrole-type `[nH]`): lone pair participates in aromaticity.
/// - Non-aromatic N bonded to C=O (amide N): lone pair delocalized into carbonyl.
/// - O with H bonded to a C=O carbon (carboxylic/ester OH).
/// - O with H bonded to oxidized S with S=O (sulfonic/sulfonamide acid OH).
/// - Oxidized S (degree > 2 or has S=O bonds): lone pair engaged in S=O resonance.
pub fn hba_count(mol: &Molecule) -> usize {
    mol.atoms()
        .filter(|(idx, atom)| {
            let an = atom.element.atomic_number();
            if an == 7 {
                // Nitrogen: charged N (N+ in nitro, quaternary, n+ in thiazolium) is never HBA.
                if atom.charge != 0 {
                    return false;
                }
                let h = implicit_hcount(mol, *idx);
                if atom.aromatic {
                    // [nH] (pyrrole-type aromatic N) is NOT an HBA
                    h == 0
                } else {
                    // Non-aromatic N: exclude amide N (bonded to C=O)
                    !neighbor_has_carbonyl(mol, *idx)
                }
            } else if an == 8 {
                // Oxygen: exclude acid OH bonded to C=O or to oxidized S with S=O
                let h = implicit_hcount(mol, *idx);
                if h > 0 {
                    !neighbor_has_carbonyl(mol, *idx)
                        && !neighbor_is_oxidized_sulfur(mol, *idx)
                } else {
                    true
                }
            } else if an == 16 {
                // Sulfur (Ertl definition includes divalent S with free lone pair)
                if atom.aromatic {
                    // Aromatic S (thiophene-type): count if uncharged
                    atom.charge == 0
                } else {
                    // Non-aromatic S: count only if divalent (X2) and not oxidized (no S=O)
                    let degree = mol.neighbors(*idx).count();
                    let total_valence = degree + implicit_hcount(mol, *idx) as usize;
                    atom.charge == 0
                        && total_valence == 2
                        && !has_double_bond_to(mol, *idx, 8)
                }
            } else {
                false
            }
        })
        .count()
}

/// True if any neighbor of `idx` is a carbon that itself has a double bond to oxygen
/// (i.e., a carbonyl carbon).
fn neighbor_has_carbonyl(mol: &Molecule, idx: AtomIdx) -> bool {
    mol.neighbors(idx).any(|(nb_idx, _)| {
        mol.atom(nb_idx).element.atomic_number() == 6 && has_double_bond_to(mol, nb_idx, 8)
    })
}

/// True if any neighbor of `idx` is a sulfur atom that itself has a S=O double bond
/// (i.e., a sulfoxide, sulfone, or sulfonate S). Used to exclude S–OH from HBA count.
fn neighbor_is_oxidized_sulfur(mol: &Molecule, idx: AtomIdx) -> bool {
    mol.neighbors(idx).any(|(nb_idx, _)| {
        mol.atom(nb_idx).element.atomic_number() == 16 && has_double_bond_to(mol, nb_idx, 8)
    })
}

/// True if any neighbor of `idx` is a carbon that has a C=N double bond
/// (i.e., an imine/amidine/guanidinium carbon).
// ---------------------------------------------------------------------------
// 6. Rotatable bond count
// ---------------------------------------------------------------------------

/// Count rotatable bonds.
///
/// A bond is rotatable when all of the following hold:
/// - It is a single bond (or a stereo bond Up/Down, which is single).
/// - Neither endpoint is terminal (degree > 1 in the heavy-atom graph).
/// - It is not part of any ring (SSSR membership).
/// - It is not an amide bond: if one atom is N and the other is C,
///   and that C has any double bond to an O, the bond is excluded.
pub fn rotatable_bond_count(mol: &Molecule) -> usize {
    let ring_bond_set = ring_bond_indices(mol);

    mol.bonds()
        .filter(|(bidx, bond)| {
            // Stereo bonds Up/Down are also single.
            let is_single = matches!(
                bond.order,
                BondOrder::Single | BondOrder::Up | BondOrder::Down
            );
            is_single
                && !ring_bond_set.contains(bidx)
                && mol.degree(bond.atom1) > 1
                && mol.degree(bond.atom2) > 1
                && !is_amide_bond(mol, bond.atom1, bond.atom2)
        })
        .count()
}

/// Indices of all bonds participating in at least one SSSR ring.
fn ring_bond_indices(mol: &Molecule) -> HashSet<BondIdx> {
    let mut set = HashSet::new();
    for ring in find_sssr(mol).rings() {
        for i in 0..ring.len() {
            let a = ring[i];
            let b = ring[(i + 1) % ring.len()];
            if let Some((bidx, _)) = mol.bond_between(a, b) {
                set.insert(bidx);
            }
        }
    }
    set
}

/// True if the bond between `a` and `b` is an amide-like C-N bond
/// (one atom is N, the other is C with a double bond to O).
fn is_amide_bond(mol: &Molecule, a: AtomIdx, b: AtomIdx) -> bool {
    let an_a = mol.atom(a).element.atomic_number();
    let an_b = mol.atom(b).element.atomic_number();
    let c_idx = match (an_a, an_b) {
        (6, 7) => a,
        (7, 6) => b,
        _ => return false,
    };
    has_double_bond_to(mol, c_idx, 8)
}

// ---------------------------------------------------------------------------
// 7. TPSA
// ---------------------------------------------------------------------------

/// Compute the topological polar surface area (Ų) using the Ertl (2000) table.
///
/// Reference: P. Ertl, B. Rohde, P. Selzer, J. Med. Chem. 2000, 43, 3714-3717.
pub fn tpsa(mol: &Molecule) -> f64 {
    let mut psa = 0.0f64;

    for (idx, atom) in mol.atoms() {
        let an = atom.element.atomic_number();
        let is_aromatic = atom.aromatic;
        let h = implicit_hcount(mol, idx);

        let contribution = match an {
            // Nitrogen
            // Reference values verified against RDKit's _CalcTPSAContribs (2024.x):
            //   [n;H0;X2] (pyridine-type, 2 bonds)           = 12.89 Ų
            //   [n;H1;X2] (pyrrole-type NH, 2 bonds)         = 15.79 Ų  ← RDKit updated value
            //   [n;H0;X≥3] (N-substituted aromatic N, 3+ bonds) = 4.93 Ų
            //   aliphatic NH2                                 = 26.02 Ų
            //   aliphatic NH (secondary)                      = 12.03 Ų
            //   aliphatic N (tertiary, no H)                  =  3.24 Ų
            7 => {
                if is_aromatic {
                    // Count the number of heavy-atom bonds on this aromatic N.
                    let degree = mol.neighbors(idx).count();
                    if h > 0 {
                        15.79 // [nH] pyrrole-type (RDKit uses 15.79, not Ertl 2000's 13.97)
                    } else if degree >= 3 {
                        // N-substituted aromatic N (3+ bonds): neutral → 4.93 Ų
                        // Quaternary aromatic N+ (thiazolium [n+], charge=1) → 3.88 Ų
                        // Confirmed: thiamine [n+]2csc(CCO)c2C, RDKit 3.88 Ų
                        if atom.charge > 0 { 3.88 } else { 4.93 }
                    } else {
                        12.89 // [n;X2]: pyridine-type aromatic N
                    }
                } else {
                    // aliphatic N
                    if atom.charge == 1 {
                        // Nitro group [N+](=O)[O-]: Ertl 2000 value = 41.44 Ų
                        // (the =O and O- oxygens contribute 0 for this group)
                        let (has_oxo, has_o_minus) = mol.neighbors(idx).fold(
                            (false, false),
                            |(oxo, om), (nb, bidx)| {
                                let nb_atom = mol.atom(nb);
                                let is_o = nb_atom.element.atomic_number() == 8;
                                (
                                    oxo || (is_o && mol.bond(bidx).order == BondOrder::Double),
                                    om || (is_o && nb_atom.charge == -1),
                                )
                            },
                        );
                        if has_oxo && has_o_minus { 41.44 } else { 3.24 }
                    } else if h >= 2 {
                        26.02 // NH2
                    } else if h == 1 {
                        // Ertl 2000: sp2 imine N-H (C=N-H in amidine/guanidinium) has
                        // the same TPSA as terminal =NH without H: 23.79 Ų.
                        // Regular secondary amine N-H (sp3, no double bond from N): 12.03 Ų.
                        let is_imine_nh = mol.neighbors(idx).any(|(nb, bidx)| {
                            mol.bond(bidx).order == BondOrder::Double
                                && mol.atom(nb).element.atomic_number() == 6
                        });
                        if is_imine_nh { 23.79 } else { 12.03 }
                    } else {
                        // h=0: tertiary N or bridged/ring imine
                        // Ertl 2000 distinguishes:
                        //   ring/bridged C=N-C (degree≥2, diazepam ring, imidazoline): 12.89 Ų
                        //   tertiary amine (no double bond): 3.24 Ų
                        let is_imine = mol.neighbors(idx).any(|(nb, bidx)| {
                            mol.bond(bidx).order == BondOrder::Double
                                && mol.atom(nb).element.atomic_number() == 6
                        });
                        if is_imine { 12.89 } else { 3.24 }
                    }
                }
            }
            // Oxygen
            8 => {
                if is_aromatic {
                    13.14
                } else if h > 0 {
                    20.23 // OH
                } else {
                    // [O-] in nitro group ([N+](=O)[O-]): contribution absorbed into N+.
                    let is_nitro_o_minus = atom.charge == -1
                        && mol.neighbors(idx).any(|(nb, _)| {
                            mol.atom(nb).element.atomic_number() == 7
                                && mol.atom(nb).charge == 1
                        });
                    if is_nitro_o_minus {
                        0.0
                    } else {
                        // S=O / P=O / N=O contributions are assigned to the heteroatom in
                        // Ertl 2000, so the doubly-bonded O is 0 for those.
                        let dbl_neighbor_an = mol
                            .neighbors(idx)
                            .find(|&(_, bidx)| mol.bond(bidx).order == BondOrder::Double)
                            .map(|(nei, _)| mol.atom(nei).element.atomic_number());
                        match dbl_neighbor_an {
                            Some(6) => 17.07, // carbonyl C=O
                            Some(_) => 0.0,   // S=O, P=O, N=O — handled by the other atom
                            None    => 9.23,  // ether O
                        }
                    }
                }
            }
            // Sulfur — Ertl 2000 atom-type contributions:
            //   aromatic S (thiophene)   = 28.24 Ų
            //   SH (thiol)               = 38.80 Ų
            //   thioether (S, 0 oxo)     = 25.30 Ų
            //   sulfoxide  (S, 1 oxo)    = 36.28 Ų  (S=O O counted as 0)
            //   sulfone/sulfonyl (2+ oxo)= 42.52 Ų  (each S=O O counted as 0)
            16 => {
                if is_aromatic {
                    28.24
                } else if h > 0 {
                    38.80 // S-H (thiol)
                } else {
                    match count_double_bonds_to(mol, idx, 8) {
                        0 => 25.30, // thioether / ring S
                        1 => 36.28, // sulfoxide
                        _ => 42.52, // sulfone, sulfonyl
                    }
                }
            }
            // Phosphorus — Ertl 2000:
            //   P=O present (phosphate, phosphonate): 26.88 Ų
            //   P=O absent (phosphine, phosphite):    34.14 Ų
            15 => {
                if !is_aromatic {
                    let has_oxo = mol.neighbors(idx).any(|(nb, bidx)| {
                        mol.bond(bidx).order == BondOrder::Double
                            && mol.atom(nb).element.atomic_number() == 8
                    });
                    if has_oxo { 26.88 } else { 34.14 }
                } else {
                    0.0
                }
            }
            _ => 0.0,
        };

        psa += contribution;
    }

    psa
}

// ---------------------------------------------------------------------------
// 8. LogP (Wildman-Crippen, calibrated)
// ---------------------------------------------------------------------------

/// Compute a Wildman-Crippen LogP using a calibrated atom-type table.
///
/// Atom type contributions are derived analytically from the RDKit reference
/// dataset (175 molecules) and confirmed against the Wildman-Crippen 1999 paper.
/// Key improvements over the simplified table:
/// - H atom contributions (H on C=+0.1230, H on N=+0.2142,
///   H on aliphatic-OH=−0.2677, H on carboxylic-OH=+0.2980)
/// - Aromatic C: [cH]=0.1581 vs [c]=0.1441
/// - Aromatic N (both [nH] and [n;H0]): −0.3239 (was +0.2626)
/// - S: thioether=+0.6482, aromatic=+0.6237 (was 0.2432/0.0)
/// - O: alcohol=−0.2893, ether=−0.0684, aromatic=+0.1552, carbonyl=−0.0509
/// - Cl: aromatic=+0.7904, aliphatic=+0.6895
/// Per-atom Crippen LogP contributions (heavy atoms only; H contributions are
/// folded into the heavy atom they are attached to). Index matches mol.atoms().
pub fn logp_crippen_per_atom(mol: &Molecule) -> Vec<f64> {
    mol.atoms().map(|(idx, atom)| {
        let an = atom.element.atomic_number();
        let ar = atom.aromatic;
        let h  = implicit_hcount(mol, idx);
        let heavy = match an {
            6  => crippen_carbon(mol, idx, ar, h),
            7  => crippen_nitrogen(mol, idx, ar),
            8  => crippen_oxygen(mol, idx, ar, h),
            16 => crippen_sulfur(mol, idx, ar),
            9  => crippen_halogen(mol, idx, ar, 0.2761, 0.4202),
            17 => crippen_halogen(mol, idx, ar, 0.7904, 0.6895),
            35 => crippen_halogen(mol, idx, ar, 0.8995, 0.8456),
            53 => crippen_halogen(mol, idx, ar, 0.7416, 0.8857),
            15 => {
                let has_oxo = mol.neighbors(idx).any(|(nb, bidx)| {
                    mol.bond(bidx).order == BondOrder::Double
                        && mol.atom(nb).element.atomic_number() == 8
                });
                if has_oxo { 0.7933 } else { -0.3451 }
            }
            _ => 0.0,
        };
        let h_contrib = if h == 0 { 0.0 } else {
            crippen_hydrogen(mol, idx, an, ar) * h as f64
        };
        heavy + h_contrib
    }).collect()
}

pub fn logp_crippen(mol: &Molecule) -> f64 {
    logp_crippen_per_atom(mol).iter().sum()
}

/// Crippen contribution for Carbon atoms.
fn crippen_carbon(mol: &Molecule, idx: AtomIdx, ar: bool, h: u8) -> f64 {
    if ar {
        // Aromatic C with exocyclic C=heteroatom (e.g., caffeine C2/C6, uracil C2/C4)
        // gets C10 = −0.3800, not the standard aromatic C value.
        let has_exocyclic_heteroatom_double = mol.neighbors(idx).any(|(nb, bidx)| {
            mol.bond(bidx).order == BondOrder::Double
                && !mol.atom(nb).aromatic
                && mol.atom(nb).element.atomic_number() != 6
        });
        if has_exocyclic_heteroatom_double {
            -0.3800
        } else if h > 0 {
            0.1581  // C11 [cH]
        } else {
            // Aromatic C bonded to non-aromatic N (aniline/amide/etc.) gets 0.4619.
            // Confirmed: aniline, diphenylamine, triphenylamine, paracetamol.
            // Does NOT apply to endocyclic aromatic N (pyridine, indole).
            if mol.neighbors(idx).any(|(nb, _)| {
                let a = mol.atom(nb);
                a.element.atomic_number() == 7 && !a.aromatic
            }) {
                return 0.4619;
            }
            // Aryl ether C (Ar-O-R): aromatic C bonded to O (single bond, no H on O).
            // RDKit per-atom analysis confirmed 0.5437 for aryl ether C vs 0.1441 for standard [c].
            let bonded_to_ether_o = mol.neighbors(idx).any(|(nb, bidx)| {
                mol.atom(nb).element.atomic_number() == 8
                    && !mol.atom(nb).aromatic
                    && mol.bond(bidx).order == BondOrder::Single
                    && implicit_hcount(mol, nb) == 0  // O has no H (ether, not phenol)
            });
            if bonded_to_ether_o {
                return 0.5437;  // Aryl ether C
            }
            // Ring-junction C (all neighbors aromatic, ≥2 of them aromatic C)
            // covers C4a/C8a in naphthalene, C3a/C7a in indole, etc.
            // Excludes caffeine C2 (only 1 aromatic C neighbor; 2 are aromatic N).
            let all_aromatic_nbrs = mol.neighbors(idx)
                .filter(|(nb, _)| mol.atom(*nb).aromatic)
                .count();
            let aromatic_c_nbrs = mol.neighbors(idx)
                .filter(|(nb, _)| mol.atom(*nb).aromatic
                    && mol.atom(*nb).element.atomic_number() == 6)
                .count();
            if all_aromatic_nbrs >= 3 && aromatic_c_nbrs >= 2 {
                0.2956  // junction [c]
            } else {
                0.1441  // C12 [c] (substituted or N-rich junction)
            }
        }
    } else {
        let has_double_to_n = has_double_bond_to(mol, idx, 7);
        let has_double_to_heteroatom = has_double_to_n || mol.neighbors(idx).any(|(nb, bidx)| {
            let bo = mol.bond(bidx).order;
            (bo == BondOrder::Double || bo == BondOrder::Triple)
                && mol.atom(nb).element.atomic_number() != 6
                && mol.atom(nb).element.atomic_number() != 7
        });
        let has_double_to_c = mol.neighbors(idx).any(|(nb, bidx)| {
            mol.bond(bidx).order == BondOrder::Double
                && !mol.atom(nb).aromatic
                && mol.atom(nb).element.atomic_number() == 6
        });

        if has_double_to_n {
            // C=N (imine, oxime, guanidine, amidine, nitrile C≡N handled here too).
            // RDKit assigns −0.2783 regardless of aryl context.
            -0.2783
        } else if has_double_to_heteroatom {
            // C=X (X = O, S, etc.) adjacent to aromatic C (Ar-CHO, Ar-COOH, Ar-CO-R)
            // takes a different value than purely aliphatic C=X. Confirmed against
            // benzoic_acid and methyl_benzoate (+0.2574 exact).
            if has_aromatic_carbon_neighbor(mol, idx) {
                -0.1226
            } else {
                -0.3800  // C10: aliphatic C=X (ketone, aldehyde, ester, carboxyl)
            }
        } else if has_double_to_c {
            // Alkene C (C=C) — Wildman-Crippen C5 type: ~+0.2274
            0.2274
        } else {
            // sp3 C: distinguish heteroatom-bonded, benzylic, pure alkyl
            let bonded_to_n = mol.neighbors(idx).any(|(nb, _)| {
                mol.atom(nb).element.atomic_number() == 7
            });
            let bonded_to_heteroatom = bonded_to_n || mol.neighbors(idx).any(|(nb, _)| {
                matches!(mol.atom(nb).element.atomic_number(), 8|9|15|16|17|35|53)
            });
            if bonded_to_heteroatom {
                // sp3 C bonded to N that is also benzylic gets 0.1193.
                // Confirmed: N-benzyl compounds, chlorpromazine side chain.
                if bonded_to_n && has_aromatic_carbon_neighbor(mol, idx) {
                    0.1193
                } else {
                    -0.2035  // C6/C7/C8: sp3 C bonded to N/O/S/halogen
                }
            } else if has_aromatic_carbon_neighbor(mol, idx) {
                // Benzylic C (Wildman-Crippen C25–C28).
                // Per-atom RDKit comparison confirmed the following corrected values:
                match h {
                    3 =>  0.0845,  // CH3-Ar (C25)  — was 0.0764
                    2 => -0.0516,  // CH2-Ar (C26)  — was -0.0597
                    1 =>  0.1193,  // CH-Ar  (C27)  — was -0.1415 (sign reversed!)
                    _ => -0.0967,  // C<-Ar  (C28)  — was -0.2037
                }
            } else {
                // Pure alkyl C: distinguish branching (C3, ≥3 C neighbors) from straight chain (C1/C2)
                // Branching C (isobutane CH, quaternary C): 0.0000 (was incorrectly 0.1441)
                // Straight-chain C (ethane CH3, propane CH2): 0.1441
                // Safety: sp2 alkene C (=C<) intercepted by has_double_to_c branch above.
                let c_nbr_count = mol.neighbors(idx)
                    .filter(|(nb, _)| mol.atom(*nb).element.atomic_number() == 6)
                    .count();
                if c_nbr_count >= 3 { 0.0000 } else { 0.1441 }
            }
        }
    }
}

/// Crippen contribution for Nitrogen atoms.
fn crippen_nitrogen(mol: &Molecule, idx: AtomIdx, ar: bool) -> f64 {
    if ar {
        // N11: all aromatic N ([nH] and [n]) use −0.3239.
        // Confirmed from pyridine, pyrrole, imidazole, pyrimidine.
        return -0.3239;
    }
    let h = implicit_hcount(mol, idx);
    let atom = mol.atom(idx);

    // Oxidized N: N=O (nitroso) or N+(=O) (nitro). Priority: check before imine.
    if has_double_bond_to(mol, idx, 8) {
        // N+ in nitro group vs neutral nitroso N
        return if atom.charge > 0 { -0.3396 } else { 0.1836 };
    }
    // Imine/guanidine/amidine N with direct C=N double bond.
    // RDKit: =N (h=0) → 0.1836, =NH (h=1) → 0.0839.
    // Does NOT apply to N merely adjacent to C=N (those use aliphatic values).
    if has_double_bond_to(mol, idx, 6) {
        return match h { 0 => 0.1836, _ => 0.0839 };
    }
    // Aryl N: bonded to aromatic C (aniline, diphenylamine, aryl amide, etc.).
    // RDKit uses same values regardless of whether N is also amide — no amide branch.
    if has_aromatic_carbon_neighbor(mol, idx) {
        return match h {
            0 => -0.4458,   // tertiary aryl N (N-methylaniline type)
            1 => -0.5188,   // secondary aryl NH (aniline NH type)
            _ => -1.0270,   // primary aryl NH2 (aniline NH2 type)
        };
    }
    // Non-aryl N adjacent to carbonyl C (amide, carbamate, urea).
    // Must come AFTER the aryl check: paracetamol's aryl amide N uses aryl values above.
    if neighbor_has_carbonyl(mol, idx) {
        return match h {
            0 => {
                // Urea-type: carbonyl C has another N neighbor (N-CO-N) → 0.0000.
                // Regular tertiary amide (DMF, DMA, BOC): carbonyl C has no other N → -0.3187.
                // Confirmed from dimethyl_urea vs DMF/DMA RDKit per-atom contributions.
                let is_urea_type = mol.neighbors(idx).any(|(cn, _)| {
                    mol.atom(cn).element.atomic_number() == 6
                        && has_double_bond_to(mol, cn, 8)
                        && mol.neighbors(cn).any(|(n2, _)| {
                            mol.atom(n2).element.atomic_number() == 7 && n2 != idx
                        })
                });
                if is_urea_type { 0.0000 } else { -0.3187 }
            }
            _ => -0.7011,  // primary/secondary amide N; confirmed from urea
        };
    }
    // Secondary aliphatic N (h=1) singly adjacent to one guanidine/amidine C=N:
    // N14 type in Wildman-Crippen (-0.335). This covers the chain NH in arginine's
    // guanidine side chain. Does NOT apply to bridge NH between two C=N groups
    // (that is doubly adjacent and treated as regular secondary amine).
    if h == 1 {
        let imine_c_nbrs = mol.neighbors(idx)
            .filter(|(nb, _)| {
                mol.atom(*nb).element.atomic_number() == 6
                    && has_double_bond_to(mol, *nb, 7)
            })
            .count();
        if imine_c_nbrs == 1 {
            return -0.335;
        }
    }
    // Aliphatic N: amide and amine use the same values in Wildman-Crippen.
    // Confirmed: dimethylformamide/acetamide tertiary amide N → -0.3187 (not 0.0).
    match h {
        0 => -0.3187,   // tertiary N (amine or amide)
        1 => -0.7096,   // secondary NH
        _ => -1.0190,   // primary NH2
    }
}

/// Crippen contribution for Oxygen atoms.
fn crippen_oxygen(mol: &Molecule, idx: AtomIdx, ar: bool, h: u8) -> f64 {
    if ar {
        0.1552           // O9: aromatic O (furan); confirmed from furan LogP=1.2796
    } else if h > 0 {
        // OH (alcohol, phenol, carboxylic acid) — H contribution handled separately.
        -0.2893
    } else {
        // Nitro group O (bonded to N+): both =O and -O- of [N+](=O)[O-] get 0.0335.
        if mol.neighbors(idx).any(|(nb, _)| {
            mol.atom(nb).element.atomic_number() == 7 && mol.atom(nb).charge > 0
        }) {
            return 0.0335;
        }
        let is_double_bonded = mol
            .neighbors(idx)
            .any(|(_, bidx)| mol.bond(bidx).order == BondOrder::Double);
        if is_double_bonded {
            -0.0509      // O8: carbonyl =O; confirmed from acetone
        } else {
            // Aryl ether O (Ar-O-R) requires special handling:
            // When ether O is bonded to aromatic C, RDKit uses distinct atomic type.
            // Confirmed from anisole and diphenyl_ether per-atom RDKit analysis.
            let bonded_to_aromatic_c = mol.neighbors(idx).any(|(nb, _)| {
                mol.atom(nb).aromatic && mol.atom(nb).element.atomic_number() == 6
            });
            if bonded_to_aromatic_c {
                -0.4195  // O: aryl ether (Ar-O-R)
            } else {
                // Carbamate/urethane ether O (N-CO-O): the adjacent C has both C=O and N.
                // This is distinct from regular ester O (C-CO-O, which has no N on the C=O carbon).
                // Confirmed from n_boc_piperazine RDKit per-atom contributions.
                let is_carbamate_o = mol.neighbors(idx).any(|(cn, _)| {
                    mol.atom(cn).element.atomic_number() == 6
                        && has_double_bond_to(mol, cn, 8)
                        && mol.neighbors(cn).any(|(n2, _)| {
                            mol.atom(n2).element.atomic_number() == 7
                        })
                });
                if is_carbamate_o { 0.4833 } else { -0.0684 }   // O4/O5: ether O
            }
        }
    }
}

/// Crippen contribution for Sulfur atoms.
fn crippen_sulfur(mol: &Molecule, idx: AtomIdx, ar: bool) -> f64 {
    if ar {
        return 0.6237;   // S3: aromatic S (thiophene); from thiophene LogP=1.7481
    }
    let h = implicit_hcount(mol, idx);
    let oxo_count = count_double_bonds_to(mol, idx, 8);
    if h > 0 && oxo_count == 0 {
        0.3132           // S4: thiol; confirmed from thiophenol, cysteine
    } else {
        match oxo_count {
            0 => 0.6482, // S1: thioether; confirmed from dimethylsulfide, THT
            1 => -0.2854,// S2: sulfoxide; derived from DMSO
            _ => -0.5684,// S3: sulfone; derived from DMSO2
        }
    }
}

/// Crippen contribution for halogens; `ar_val` when on aromatic ring, `al_val` on aliphatic.
fn crippen_halogen(mol: &Molecule, idx: AtomIdx, ar: bool, ar_val: f64, al_val: f64) -> f64 {
    if ar || has_aromatic_neighbor(mol, idx) { ar_val } else { al_val }
}

/// Crippen H-atom contribution per hydrogen on atom `idx`.
fn crippen_hydrogen(mol: &Molecule, idx: AtomIdx, an: u8, ar: bool) -> f64 {
    match an {
        6 => 0.1230,   // H1: H on C (any hybridization); from alkane series
        7 => 0.2142,   // H2: H on N; from pyrrole, imidazole
        8 => {
            if ar {
                0.1125
            } else if has_aromatic_neighbor(mol, idx) {
                // Phenolic OH (no carbonyl) — confirmed: phenol, catechol,
                // resorcinol, hydroquinone, salicylic_acid (via fix1+fix2), dopamine.
                0.1319
            } else if neighbor_has_carbonyl(mol, idx) {
                0.2980 // H3: carboxylic/ester OH; from acetic acid
            } else {
                -0.2677 // H4: aliphatic alcohol OH; from methanol, ethanol
            }
        }
        _ => 0.1125,
    }
}

// ---------------------------------------------------------------------------
// 9. Lipinski Rule of Five
// ---------------------------------------------------------------------------

/// Apply Lipinski's Rule of Five.
///
/// Returns `true` when all four criteria are satisfied:
/// - Molecular weight ≤ 500 Da
/// - H-bond donors ≤ 5
/// - H-bond acceptors ≤ 10
/// - Crippen LogP ≤ 5.0
pub fn lipinski_passes(mol: &Molecule) -> bool {
    molecular_weight(mol) <= 500.0
        && hbd_count(mol) <= 5
        && hba_count(mol) <= 10
        && logp_crippen(mol) <= 5.0
}

// ---------------------------------------------------------------------------
// Fsp3 — fraction of sp3 carbons
// ---------------------------------------------------------------------------

/// Fraction of sp3 carbons: sp3_C / total_C.
///
/// sp3 carbon is defined as a non-aromatic carbon that has no double or triple
/// bond to any neighbour (i.e. hybridisation is effectively sp3).
/// Returns 0.0 if the molecule contains no carbon atoms.
pub fn fsp3(mol: &Molecule) -> f64 {
    let c_total = mol
        .atoms()
        .filter(|(_, a)| a.element.atomic_number() == 6)
        .count();
    if c_total == 0 {
        return 0.0;
    }
    let sp3 = mol
        .atoms()
        .filter(|(idx, a)| {
            a.element.atomic_number() == 6
                && !a.aromatic
                && mol.neighbors(*idx).all(|(_, bidx)| {
                    !matches!(mol.bond(bidx).order, BondOrder::Double | BondOrder::Triple)
                })
        })
        .count();
    sp3 as f64 / c_total as f64
}

// ---------------------------------------------------------------------------
// Aromatic ring count
// ---------------------------------------------------------------------------

/// Number of aromatic rings in the molecule (from SSSR).
///
/// A ring is considered aromatic when every atom in it carries the
/// `aromatic` flag.
pub fn aromatic_ring_count(mol: &Molecule) -> usize {
    find_sssr(mol)
        .rings()
        .iter()
        .filter(|ring| ring.iter().all(|&idx| mol.atom(idx).aromatic))
        .count()
}

// ---------------------------------------------------------------------------
// 11. Formal charge sum
// ---------------------------------------------------------------------------

/// Sum of all formal charges in the molecule.
pub fn formal_charge_sum(mol: &Molecule) -> i32 {
    mol.atoms().map(|(_, a)| a.charge as i32).sum()
}

// ---------------------------------------------------------------------------
// 12. Molar Refractivity (Wildman-Crippen additive model)
//
// Atom-type contributions taken from RDKit's Crippen.txt (Wildman & Crippen 1999).
// ---------------------------------------------------------------------------

fn mr_carbon(mol: &Molecule, idx: AtomIdx, ar: bool, h: u8) -> f64 {
    if ar {
        if h > 0 { 3.35 } else { 3.50 }          // C18 / avg of C19-C25
    } else {
        let has_double_to_heteroatom = mol.neighbors(idx).any(|(nb, bidx)| {
            let bo = mol.bond(bidx).order;
            (bo == BondOrder::Double || bo == BondOrder::Triple)
                && mol.atom(nb).element.atomic_number() != 6
        });
        let has_double_to_c = mol.neighbors(idx).any(|(nb, bidx)| {
            mol.bond(bidx).order == BondOrder::Double
                && !mol.atom(nb).aromatic
                && mol.atom(nb).element.atomic_number() == 6
        });
        if has_double_to_heteroatom {
            5.007   // C5: sp2 C=X (carbonyl, imine, thiocarbonyl, …)
        } else if has_double_to_c {
            3.513   // C6: alkene C
        } else {
            let bonded_to_heteroatom = mol.neighbors(idx).any(|(nb, _)| {
                matches!(mol.atom(nb).element.atomic_number(), 7|8|9|15|16|17|35|53)
            });
            if bonded_to_heteroatom { 2.753 } else { 2.503 }   // C3 / C1
        }
    }
}

fn mr_nitrogen(mol: &Molecule, idx: AtomIdx, ar: bool) -> f64 {
    if ar { return 2.202; }   // N11
    let h = implicit_hcount(mol, idx);
    match h {
        0 => 1.839,   // N7: tertiary amine
        1 => 2.173,   // N2: secondary amine
        _ => 2.262,   // N1: primary amine
    }
}

fn mr_oxygen(mol: &Molecule, idx: AtomIdx, ar: bool, h: u8) -> f64 {
    if ar { return 1.08; }       // O1: aromatic o (furan)
    if h > 0 { return 0.8238; }  // O2: OH
    let is_double = mol
        .neighbors(idx)
        .any(|(_, bidx)| mol.bond(bidx).order == BondOrder::Double);
    if is_double { 0.0 } else { 1.085 }  // O9 carbonyl O / O3 ether O
}

fn mr_sulfur(_mol: &Molecule, _idx: AtomIdx, ar: bool) -> f64 {
    if ar { 6.691 } else { 7.591 }    // S3 aromatic / S1 thioether
}

fn mr_hydrogen(mol: &Molecule, idx: AtomIdx, an: u8, ar: bool) -> f64 {
    match an {
        6 => 1.057,    // H1
        7 => 0.9627,   // H3
        8 => {
            if ar { 1.112 }
            else if neighbor_has_carbonyl(mol, idx) { 1.805 } // H4: COOH / ester OH
            else { 1.395 }  // H2: alcohol / phenol OH
        }
        _ => 1.112,    // HS fallback
    }
}

/// Per-atom Molar Refractivity contributions (Wildman & Crippen 1999).
/// H contributions folded into the attached heavy atom. Index matches mol.atoms().
pub fn mr_per_atom(mol: &Molecule) -> Vec<f64> {
    mol.atoms().map(|(idx, atom)| {
        let an = atom.element.atomic_number();
        let ar = atom.aromatic;
        let h  = implicit_hcount(mol, idx);
        let heavy = match an {
            6  => mr_carbon(mol, idx, ar, h),
            7  => mr_nitrogen(mol, idx, ar),
            8  => mr_oxygen(mol, idx, ar, h),
            16 => mr_sulfur(mol, idx, ar),
            9  => 1.108,
            17 => 5.853,
            35 => 8.927,
            53 => 14.02,
            15 => 6.920,
            _  => 3.243,
        };
        let h_contrib = if h == 0 { 0.0 } else {
            mr_hydrogen(mol, idx, an, ar) * h as f64
        };
        heavy + h_contrib
    }).collect()
}

/// Compute Molar Refractivity using the Wildman-Crippen additive model.
///
/// Uses the same atom-type framework as `logp_crippen` but with MR contributions
/// from Wildman & Crippen 1999 (J. Chem. Inf. Comput. Sci. 39, 868-873).
pub fn molar_refractivity(mol: &Molecule) -> f64 {
    mr_per_atom(mol).iter().sum()
}

// ---------------------------------------------------------------------------
// 13. Basic count descriptors
// ---------------------------------------------------------------------------

/// Number of heteroatoms (non-C, non-H heavy atoms).
pub fn num_heteroatoms(mol: &Molecule) -> usize {
    mol.atoms()
        .filter(|(_, a)| {
            let an = a.element.atomic_number();
            an != 1 && an != 6
        })
        .count()
}

/// Total number of rings (SSSR count).
pub fn ring_count(mol: &Molecule) -> usize {
    find_sssr(mol).rings().len()
}

/// Number of non-aromatic (aliphatic) rings.
pub fn num_aliphatic_rings(mol: &Molecule) -> usize {
    find_sssr(mol)
        .rings()
        .iter()
        .filter(|ring| ring.iter().any(|&idx| !mol.atom(idx).aromatic))
        .count()
}

/// Number of saturated (all sp3) rings.
///
/// A ring is saturated when every atom has no double or triple bonds to
/// neighbors (regardless of aromaticity flag).
pub fn num_saturated_rings(mol: &Molecule) -> usize {
    find_sssr(mol)
        .rings()
        .iter()
        .filter(|ring| {
            ring.iter().all(|&idx| {
                mol.neighbors(idx).all(|(_, bidx)| {
                    !matches!(mol.bond(bidx).order, BondOrder::Double | BondOrder::Triple | BondOrder::Aromatic)
                })
            })
        })
        .count()
}

/// Number of aromatic rings containing at least one heteroatom (N, O, S, P, …).
///
/// Examples: pyridine (1), furan (1), imidazole (1), benzene (0).
pub fn num_aromatic_heterocycles(mol: &Molecule) -> usize {
    find_sssr(mol).rings().iter().filter(|ring| {
        ring.iter().all(|&idx| mol.atom(idx).aromatic)
            && ring.iter().any(|&idx| {
                let an = mol.atom(idx).element.atomic_number();
                an != 6 && an != 1
            })
    }).count()
}

/// Number of non-aromatic rings containing at least one heteroatom.
///
/// A ring is aliphatic when at least one of its atoms is not aromatic.
/// Examples: piperidine (1), morpholine (1), tetrahydrofuran (1).
pub fn num_aliphatic_heterocycles(mol: &Molecule) -> usize {
    find_sssr(mol).rings().iter().filter(|ring| {
        ring.iter().any(|&idx| !mol.atom(idx).aromatic)
            && ring.iter().any(|&idx| {
                let an = mol.atom(idx).element.atomic_number();
                an != 6 && an != 1
            })
    }).count()
}

/// Number of fully saturated rings (no unsaturated bonds) containing at least one heteroatom.
///
/// Examples: piperidine (1), oxetane (1), azetidine (1).
/// A ring with any double, triple, or aromatic bond is not saturated.
pub fn num_saturated_heterocycles(mol: &Molecule) -> usize {
    find_sssr(mol).rings().iter().filter(|ring| {
        ring.iter().all(|&idx| {
            mol.neighbors(idx).all(|(_, bidx)| {
                !matches!(mol.bond(bidx).order,
                    BondOrder::Double | BondOrder::Triple | BondOrder::Aromatic)
            })
        }) && ring.iter().any(|&idx| {
            let an = mol.atom(idx).element.atomic_number();
            an != 6 && an != 1
        })
    }).count()
}

/// Number of spiro atoms.
///
/// A spiro atom belongs to exactly 2 rings and is the sole shared atom between them.
/// Example: spiro[4.5]decane (`C1CCCCC11CCCC1`) has 1 spiro atom.
pub fn num_spiro_atoms(mol: &Molecule) -> usize {
    let sssr = find_sssr(mol);
    let rings = sssr.rings();
    mol.atoms().filter(|(idx, _)| {
        let member: Vec<_> = rings.iter().filter(|r| r.contains(idx)).collect();
        if member.len() != 2 { return false; }
        // Spiro: the two rings share exactly this one atom (no shared bond = not fused).
        member[0].iter().filter(|a| member[1].contains(a)).count() == 1
    }).count()
}

/// Number of bridgehead atoms.
///
/// A bridgehead atom belongs to 2 or more rings and has 3 or more bonds to other
/// ring atoms, where the rings it belongs to share at least one pair of atoms that
/// are NOT directly bonded (distinguishing bridged from fused or spiro systems).
///
/// Example: norbornane (`C1CC2CCC1C2`) has 2 bridgehead atoms.
/// Naphthalene has 0 (the junction atoms are fused — directly bonded to each other).
/// Spiro[4.5]decane has 0 (the spiro center is not bridged).
pub fn num_bridgehead_atoms(mol: &Molecule) -> usize {
    let sssr = find_sssr(mol);
    let rings = sssr.rings();
    mol.atoms().filter(|(idx, _)| {
        if sssr.atoms_in_ring_count(*idx) < 2 { return false; }
        let ring_bonds = mol.neighbors(*idx)
            .filter(|(nb, _)| sssr.contains_atom(*nb))
            .count();
        if ring_bonds < 3 { return false; }
        let member_rings: Vec<_> = rings.iter().filter(|r| r.contains(idx)).collect();
        for i in 0..member_rings.len() {
            for j in (i + 1)..member_rings.len() {
                let shared: Vec<AtomIdx> = member_rings[i].iter()
                    .filter(|a| member_rings[j].contains(a))
                    .copied()
                    .collect();
                if shared.len() < 2 { continue; }
                // Skip when one ring is entirely contained in the other — this is an artifact
                // of the SSSR returning a symmetric-difference ring instead of the minimal ring.
                // Real ring pairs always have |shared| < min(|R_i|, |R_j|).
                if shared.len() == member_rings[i].len() || shared.len() == member_rings[j].len() {
                    continue;
                }
                // If any pair of shared atoms is not directly bonded, this is a bridge junction.
                for a in 0..shared.len() {
                    for b in (a + 1)..shared.len() {
                        if mol.bond_between(shared[a], shared[b]).is_none() {
                            return true;
                        }
                    }
                }
            }
        }
        false
    }).count()
}

/// Number of assigned stereocenters (tetrahedral R/S from CIP assignment).
///
/// Runs `assign_cip` internally. Returns 0 for molecules with no stereo.
pub fn num_stereocenters(mol: &Molecule) -> usize {
    use chematic_core::CipCode;
    crate::cip::assign_cip(mol)
        .assignments
        .iter()
        .filter(|(_, c)| matches!(c, CipCode::R | CipCode::S))
        .count()
}

/// Number of unspecified (undefined) stereocenters.
///
/// Counts sp3 carbons with exactly 4 substituents whose chirality is not
/// specified (no @/@@ in the input SMILES).
pub fn num_unspecified_stereocenters(mol: &Molecule) -> usize {
    use chematic_core::Chirality;
    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();
            let total = degree + implicit_hcount(mol, *idx) as usize;
            total == 4
                && mol.neighbors(*idx).all(|(_, bidx)| {
                    !matches!(mol.bond(bidx).order, BondOrder::Double | BondOrder::Triple)
                })
        })
        .count()
}

// ---------------------------------------------------------------------------
// 14. Drug-likeness filters
// ---------------------------------------------------------------------------

/// Veber (2002) oral bioavailability filter.
///
/// Passes when TPSA ≤ 140 Ų **and** rotatable bonds ≤ 10.
pub fn veber_passes(mol: &Molecule) -> bool {
    tpsa(mol) <= 140.0 && rotatable_bond_count(mol) <= 10
}

/// Egan (2000) absorption/permeability filter ("Egg model").
///
/// Passes when TPSA ≤ 131.6 Ų **and** Crippen LogP ≤ 5.88.
pub fn egan_passes(mol: &Molecule) -> bool {
    tpsa(mol) <= 131.6 && logp_crippen(mol) <= 5.88
}

/// REOS (Rapid Elimination Of Swill) filter.
///
/// Six criteria for hit-identification library quality:
/// MW 200–500, LogP −5 to 5, HBD 0–5, HBA 0–10, formal charge −2 to 2,
/// rotatable bonds 0–8, heavy atoms 15–50.
pub fn reos_passes(mol: &Molecule) -> bool {
    let mw   = molecular_weight(mol);
    let lp   = logp_crippen(mol);
    let hbd  = hbd_count(mol) as i32;
    let hba  = hba_count(mol) as i32;
    let fc   = formal_charge_sum(mol);
    let rotb = rotatable_bond_count(mol) as i32;
    let hac  = heavy_atom_count(mol) as i32;

    (200.0..=500.0).contains(&mw)
        && (-5.0..=5.0).contains(&lp)
        && (0..=5).contains(&hbd)
        && (0..=10).contains(&hba)
        && (-2..=2).contains(&fc)
        && (0..=8).contains(&rotb)
        && (15..=50).contains(&hac)
}

/// Ghose (1999) drug-likeness filter.
///
/// Four criteria: MW 160–480, LogP −0.4 to 5.6, heavy atoms 20–70,
/// Molar Refractivity 40–130.
pub fn ghose_passes(mol: &Molecule) -> bool {
    let mw  = molecular_weight(mol);
    let lp  = logp_crippen(mol);
    let hac = heavy_atom_count(mol) as f64;
    let mr  = molar_refractivity(mol);

    (160.0..=480.0).contains(&mw)
        && (-0.4..=5.6).contains(&lp)
        && (20.0..=70.0).contains(&hac)
        && (40.0..=130.0).contains(&mr)
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

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

    /// Parse a SMILES string, panicking on failure.
    fn mol(smiles: &str) -> Molecule {
        parse(smiles).unwrap_or_else(|e| panic!("failed to parse {smiles:?}: {e}"))
    }

    // Tolerance helpers.
    fn approx(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() <= tol
    }

    fn pct2(a: f64, b: f64) -> bool {
        // within 2% relative, or within 0.05 Da absolute (for very small values)
        approx(a, b, b.abs() * 0.02 + 0.05)
    }

    // -- Test 1: methane molecular weight ------------------------------------
    #[test]
    fn test_mw_methane() {
        let m = mol("C");
        // CH4: 12.011 + 4*1.008 = 16.043
        assert!(pct2(molecular_weight(&m), 16.043), "methane MW = {}", molecular_weight(&m));
    }

    // -- Test 2: water molecular weight -------------------------------------
    #[test]
    fn test_mw_water() {
        let m = mol("O");
        // H2O: 15.999 + 2*1.008 = 18.015
        assert!(pct2(molecular_weight(&m), 18.015), "water MW = {}", molecular_weight(&m));
    }

    // -- Test 3: ethanol molecular weight -----------------------------------
    #[test]
    fn test_mw_ethanol() {
        let m = mol("CCO");
        // C2H6O: 2*12.011 + 6*1.008 + 15.999 = 46.068
        assert!(pct2(molecular_weight(&m), 46.068), "ethanol MW = {}", molecular_weight(&m));
    }

    // -- Test 4: benzene molecular weight -----------------------------------
    #[test]
    fn test_mw_benzene() {
        let m = mol("c1ccccc1");
        // C6H6: 6*12.011 + 6*1.008 = 78.114
        assert!(pct2(molecular_weight(&m), 78.114), "benzene MW = {}", molecular_weight(&m));
    }

    // -- Test 5: aspirin molecular weight -----------------------------------
    #[test]
    fn test_mw_aspirin() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        // C9H8O4: MW ~180.16
        let mw = molecular_weight(&m);
        assert!(approx(mw, 180.16, 1.0), "aspirin MW = {mw}");
    }

    // -- Test 6: methane exact mass -----------------------------------------
    #[test]
    fn test_exact_mass_methane() {
        let m = mol("C");
        // 12C + 4*(1H): 12.0000 + 4*1.00783 = 16.0313
        let em = exact_mass(&m);
        assert!(approx(em, 16.031, 0.01), "methane exact mass = {em}");
    }

    // -- Test 7: benzene heavy atom count -----------------------------------
    #[test]
    fn test_hac_benzene() {
        let m = mol("c1ccccc1");
        assert_eq!(heavy_atom_count(&m), 6);
    }

    // -- Test 8: aspirin heavy atom count -----------------------------------
    #[test]
    fn test_hac_aspirin() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        // C9H8O4: 9 C + 4 O = 13 heavy atoms
        assert_eq!(heavy_atom_count(&m), 13);
    }

    // -- Test 9: ethanol HBD ------------------------------------------------
    #[test]
    fn test_hbd_ethanol() {
        let m = mol("CCO");
        assert_eq!(hbd_count(&m), 1); // one OH
    }

    // -- Test 10: aniline HBD -----------------------------------------------
    #[test]
    fn test_hbd_aniline() {
        let m = mol("Nc1ccccc1");
        assert_eq!(hbd_count(&m), 1); // one NH2
    }

    // -- Test 11: benzene HBD -----------------------------------------------
    #[test]
    fn test_hbd_benzene() {
        let m = mol("c1ccccc1");
        assert_eq!(hbd_count(&m), 0);
    }

    // -- Test 12: ethanol HBA -----------------------------------------------
    #[test]
    fn test_hba_ethanol() {
        let m = mol("CCO");
        assert_eq!(hba_count(&m), 1); // one O
    }

    // -- Test 13: aspirin HBA -----------------------------------------------
    #[test]
    fn test_hba_aspirin() {
        // aspirin: 4 O total, but carboxyl OH is excluded → 3 (RDKit-aligned)
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        assert_eq!(hba_count(&m), 3);
    }

    // -- Test 14: benzene rotatable bonds ------------------------------------
    #[test]
    fn test_rot_benzene() {
        let m = mol("c1ccccc1");
        assert_eq!(rotatable_bond_count(&m), 0);
    }

    // -- Test 15: aspirin rotatable bonds ------------------------------------
    #[test]
    fn test_rot_aspirin() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        // Rotatable: CH3-C(=O), C(=O)-O (ester oxygen), O-aryl-C
        // Non-rotatable: ring bonds, C=O double bonds, terminal CH3 (degree 1)
        // Expected: 3
        let r = rotatable_bond_count(&m);
        assert_eq!(r, 3, "aspirin rotatable bonds = {r}");
    }

    // -- Test 16: water TPSA -------------------------------------------------
    #[test]
    fn test_tpsa_water() {
        let m = mol("O");
        // single O with 2H → 20.23
        let t = tpsa(&m);
        assert!(approx(t, 20.23, 1.0), "water TPSA = {t}");
    }

    // -- Test 17: aniline TPSA -----------------------------------------------
    #[test]
    fn test_tpsa_aniline() {
        let m = mol("Nc1ccccc1");
        // NH2 (aliphatic) → 26.02
        let t = tpsa(&m);
        assert!(approx(t, 26.02, 5.0), "aniline TPSA = {t}");
    }

    // -- Test 18: aspirin Lipinski -------------------------------------------
    #[test]
    fn test_lipinski_aspirin() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        assert!(lipinski_passes(&m));
    }

    // -- Test 19: benzene Lipinski ------------------------------------------
    #[test]
    fn test_lipinski_benzene() {
        let m = mol("c1ccccc1");
        assert!(lipinski_passes(&m));
    }

    // -- Additional tests ---------------------------------------------------

    // Benzene exact mass
    #[test]
    fn test_exact_mass_benzene() {
        let m = mol("c1ccccc1");
        // C6H6: 6*12 + 6*1.00783 = 78.04698
        let em = exact_mass(&m);
        assert!(approx(em, 78.047, 0.05), "benzene exact mass = {em}");
    }

    // Ethanol exact mass
    #[test]
    fn test_exact_mass_ethanol() {
        let m = mol("CCO");
        // C2H6O: 2*12 + 6*1.00783 + 15.9949 = 46.0419
        let em = exact_mass(&m);
        assert!(approx(em, 46.042, 0.05), "ethanol exact mass = {em}");
    }

    // Aspirin logp and Lipinski components
    #[test]
    fn test_logp_aspirin_is_reasonable() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        let lp = logp_crippen(&m);
        // The simplified model gives a low but not absurd value; just check range.
        assert!(lp > -5.0 && lp < 5.0, "aspirin logp = {lp}");
    }

    // Heavy atom count for ethanol
    #[test]
    fn test_hac_ethanol() {
        let m = mol("CCO");
        assert_eq!(heavy_atom_count(&m), 3); // 2 C + 1 O
    }

    // HBA for aniline (one N)
    #[test]
    fn test_hba_aniline() {
        let m = mol("Nc1ccccc1");
        assert_eq!(hba_count(&m), 1); // one N
    }

    // Rotatable bonds for n-butane (single chain)
    #[test]
    fn test_rot_butane() {
        let m = mol("CCCC");
        // C1-C2, C2-C3, C3-C4 → three bonds; C1 has degree 1 (only C2 neighbor), C4 too
        // degree(C1)=1 → excluded; degree(C4)=1 → excluded
        // Only C2-C3 is non-terminal on both sides? Actually C2 has degree 2 (C1,C3),
        // C3 has degree 2 (C2,C4). Bond C2-C3: both non-terminal → rotatable (1).
        // Bond C1-C2: C1 degree 1 → skip. Bond C3-C4: C4 degree 1 → skip.
        assert_eq!(rotatable_bond_count(&m), 1, "n-butane has 1 rotatable bond");
    }

    // TPSA for aspirin (should be > 0)
    #[test]
    fn test_tpsa_aspirin_positive() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        let t = tpsa(&m);
        assert!(t > 0.0, "aspirin TPSA = {t}");
    }

    // -- Fsp3 tests --------------------------------------------------------

    #[test]
    fn test_fsp3_benzene() {
        let m = mol("c1ccccc1");
        // all aromatic C, no sp3
        assert!((fsp3(&m) - 0.0).abs() < 1e-9, "benzene Fsp3 should be 0");
    }

    #[test]
    fn test_fsp3_cyclohexane() {
        let m = mol("C1CCCCC1");
        // all sp3 C
        assert!((fsp3(&m) - 1.0).abs() < 1e-9, "cyclohexane Fsp3 should be 1");
    }

    #[test]
    fn test_fsp3_aspirin() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        // 9 C total: 1 sp3 (methyl), 2 sp2 C=O, 6 aromatic
        // sp3 = 1, total C = 9 → Fsp3 = 1/9 ≈ 0.111
        let f = fsp3(&m);
        assert!(f > 0.05 && f < 0.25, "aspirin Fsp3={f} expected ~0.111");
    }

    #[test]
    fn test_fsp3_no_carbon() {
        let m = mol("[NH4+]");
        assert!((fsp3(&m) - 0.0).abs() < 1e-9, "no-carbon mol Fsp3 should be 0");
    }

    // -- aromatic_ring_count tests -----------------------------------------

    #[test]
    fn test_aromatic_ring_count_benzene() {
        let m = mol("c1ccccc1");
        assert_eq!(aromatic_ring_count(&m), 1);
    }

    #[test]
    fn test_aromatic_ring_count_naphthalene() {
        let m = mol("c1ccc2ccccc2c1");
        assert_eq!(aromatic_ring_count(&m), 2);
    }

    #[test]
    fn test_aromatic_ring_count_cyclohexane() {
        let m = mol("C1CCCCC1");
        assert_eq!(aromatic_ring_count(&m), 0);
    }

    #[test]
    fn test_aromatic_ring_count_aspirin() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        assert_eq!(aromatic_ring_count(&m), 1);
    }

    // -- formal_charge_sum tests -------------------------------------------

    #[test]
    fn test_formal_charge_neutral_aspirin() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        assert_eq!(formal_charge_sum(&m), 0);
    }

    #[test]
    fn test_formal_charge_quaternary_n() {
        // Trimethylammonium: [N+]
        let m = mol("CC[N+](C)(C)C");
        assert_eq!(formal_charge_sum(&m), 1);
    }

    #[test]
    fn test_formal_charge_zwitterion() {
        // Glycine zwitterion: [NH3+]CC(=O)[O-]
        let m = mol("[NH3+]CC(=O)[O-]");
        assert_eq!(formal_charge_sum(&m), 0);
    }

    // -- molar_refractivity tests -------------------------------------------

    #[test]
    fn test_mr_benzene_range() {
        // Benzene MR ≈ 26.0 (RDKit reference)
        let m = mol("c1ccccc1");
        let mr = molar_refractivity(&m);
        assert!(mr > 20.0 && mr < 35.0, "benzene MR={mr:.2}");
    }

    #[test]
    fn test_mr_aspirin_range() {
        // Aspirin MR ≈ 45.5 (Ghose filter: 40-130 → should pass)
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        let mr = molar_refractivity(&m);
        assert!(mr > 35.0 && mr < 65.0, "aspirin MR={mr:.2}");
    }

    #[test]
    fn test_mr_chlorobenzene_higher_than_benzene() {
        // Cl adds ~5.85 to MR
        let m_bz  = mol("c1ccccc1");
        let m_clb = mol("c1ccc(Cl)cc1");
        assert!(
            molar_refractivity(&m_clb) > molar_refractivity(&m_bz),
            "chlorobenzene should have higher MR than benzene"
        );
    }

    // -- Veber filter tests -------------------------------------------------

    #[test]
    fn test_veber_aspirin_passes() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        assert!(veber_passes(&m), "aspirin should pass Veber filter");
    }

    #[test]
    fn test_veber_large_flexible_fails() {
        // A molecule with many rotatable bonds should fail
        let m = mol("CCCCCCCCCCCCC(=O)O");  // myristic acid - 12 rotatable bonds
        let rotb = rotatable_bond_count(&m);
        if rotb > 10 {
            assert!(!veber_passes(&m), "myristic acid (rotb={rotb}) should fail Veber");
        }
    }

    // -- Egan filter tests --------------------------------------------------

    #[test]
    fn test_egan_aspirin_passes() {
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        assert!(egan_passes(&m), "aspirin should pass Egan filter");
    }

    // -- REOS filter tests --------------------------------------------------

    #[test]
    fn test_reos_aspirin_passes() {
        // Aspirin: MW=180, logP~1.2, HBD=1, HBA=3, charge=0, rotb=3, HAC=13
        // HAC=13 < 15 → REOS fails! (aspirin is small)
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        // Just test that we can call it without panicking
        let _ = reos_passes(&m);
    }

    #[test]
    fn test_reos_ibuprofen_passes() {
        // Ibuprofen: MW=206, logP~3.5, HBD=1, HBA=1, charge=0, rotb=4, HAC=13
        // HAC=13 < 15 → likely fails
        let m = mol("CC(C)Cc1ccc(cc1)C(C)C(=O)O");
        let _ = reos_passes(&m);
    }

    #[test]
    fn test_reos_diazepam_passes() {
        // Diazepam: MW~285, logP~2.9, HBD=0, HBA=2, charge=0, rotb=1, HAC~22 — all in range
        let m = mol("CN1C(=O)CN=C(c2ccccc2)c2cc(Cl)ccc21");
        assert!(reos_passes(&m), "diazepam should pass REOS");
    }

    // -- Ghose filter tests -------------------------------------------------

    #[test]
    fn test_ghose_aspirin_range() {
        // Aspirin MW=180 is below Ghose MW lower bound (160 ≤ MW ≤ 480). Should pass MW.
        // HAC=13 < 20 lower bound → Ghose fails for aspirin
        let m = mol("CC(=O)Oc1ccccc1C(=O)O");
        // Just verify it doesn't panic
        let _ = ghose_passes(&m);
    }

    #[test]
    fn test_ghose_ibuprofen_passes() {
        // Ibuprofen: MW=206, logP~3.5, HAC=13 (borderline)
        let m = mol("CC(C)Cc1ccc(cc1)C(C)C(=O)O");
        let _ = ghose_passes(&m);
    }

    // -- Basic count descriptor tests ----------------------------------------

    #[test]
    fn test_num_heteroatoms_aspirin() {
        // Aspirin: 4 O atoms (ester C=O, ester O, carboxylic C=O, carboxylic OH)
        assert_eq!(num_heteroatoms(&mol("CC(=O)Oc1ccccc1C(=O)O")), 4);
    }

    #[test]
    fn test_num_heteroatoms_benzene_zero() {
        assert_eq!(num_heteroatoms(&mol("c1ccccc1")), 0);
    }

    #[test]
    fn test_ring_count_benzene() {
        assert_eq!(ring_count(&mol("c1ccccc1")), 1);
    }

    #[test]
    fn test_ring_count_naphthalene() {
        assert_eq!(ring_count(&mol("c1ccc2ccccc2c1")), 2);
    }

    #[test]
    fn test_ring_count_acyclic_zero() {
        assert_eq!(ring_count(&mol("CCO")), 0);
    }

    #[test]
    fn test_num_saturated_rings_cyclohexane() {
        assert_eq!(num_saturated_rings(&mol("C1CCCCC1")), 1);
    }

    #[test]
    fn test_num_saturated_rings_benzene_zero() {
        assert_eq!(num_saturated_rings(&mol("c1ccccc1")), 0);
    }

    #[test]
    fn test_num_aliphatic_rings_cyclohexane() {
        assert_eq!(num_aliphatic_rings(&mol("C1CCCCC1")), 1);
    }

    #[test]
    fn test_num_aliphatic_rings_benzene_zero() {
        assert_eq!(num_aliphatic_rings(&mol("c1ccccc1")), 0);
    }

    #[test]
    fn test_num_stereocenters_alanine() {
        // L-alanine: 1 R/S center
        assert_eq!(num_stereocenters(&mol("[C@@H](N)(C)C(=O)O")), 1);
    }

    #[test]
    fn test_num_stereocenters_achiral_zero() {
        assert_eq!(num_stereocenters(&mol("CC(=O)O")), 0);
    }


}