#![forbid(unsafe_code)]
use std::collections::HashSet;
use chematic_core::{AtomIdx, BondIdx, BondOrder, Molecule, MoleculeBuilder, implicit_hcount};
#[derive(Clone, Copy, PartialEq, Eq)]
enum BondOrderMatch {
Single,
Double,
#[allow(dead_code)]
Any,
}
impl BondOrderMatch {
fn matches(self, order: BondOrder) -> bool {
match self {
BondOrderMatch::Single => {
matches!(order, BondOrder::Single | BondOrder::Up | BondOrder::Down)
}
BondOrderMatch::Double => matches!(order, BondOrder::Double),
BondOrderMatch::Any => true,
}
}
}
struct TautomerRule {
#[allow(dead_code)]
name: &'static str,
donor_elem: u8,
bridge_elem: Option<u8>,
acceptor_elem: u8,
donor_bridge_order: BondOrderMatch,
bridge_acceptor_order: BondOrderMatch,
prefer_forward: bool,
}
static RULES: &[TautomerRule] = &[
TautomerRule {
name: "keto-enol",
donor_elem: 8, bridge_elem: Some(6), acceptor_elem: 6,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: true,
},
TautomerRule {
name: "amide-iminol",
donor_elem: 7, bridge_elem: Some(6), acceptor_elem: 8,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "iminol-amide",
donor_elem: 8, bridge_elem: Some(6), acceptor_elem: 7,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: true,
},
TautomerRule {
name: "imine-enamine",
donor_elem: 7, bridge_elem: Some(6), acceptor_elem: 6,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: true,
},
TautomerRule {
name: "1,3-N-to-O",
donor_elem: 7, bridge_elem: None, acceptor_elem: 8,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-N-to-N",
donor_elem: 7, bridge_elem: None, acceptor_elem: 7,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "thioamide",
donor_elem: 7, bridge_elem: Some(6), acceptor_elem: 16,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "thio-iminol-amide",
donor_elem: 16, bridge_elem: Some(6), acceptor_elem: 7,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: true,
},
TautomerRule {
name: "thio-keto-enol",
donor_elem: 16, bridge_elem: Some(6), acceptor_elem: 6,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: true,
},
TautomerRule {
name: "thio-enol-ketone",
donor_elem: 8, bridge_elem: Some(6), acceptor_elem: 16,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-N-to-S",
donor_elem: 7, bridge_elem: None, acceptor_elem: 16,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-S-to-O",
donor_elem: 16, bridge_elem: None, acceptor_elem: 8,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-S-to-N",
donor_elem: 16, bridge_elem: None, acceptor_elem: 7,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-O-to-S",
donor_elem: 8, bridge_elem: None, acceptor_elem: 16,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-S-to-S",
donor_elem: 16, bridge_elem: None, acceptor_elem: 16,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-O-to-N-any-bridge",
donor_elem: 8, bridge_elem: None, acceptor_elem: 7,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: true,
},
TautomerRule {
name: "1,3-O-to-O-any-bridge",
donor_elem: 8, bridge_elem: None, acceptor_elem: 8,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-N-to-C-any-bridge",
donor_elem: 7, bridge_elem: None, acceptor_elem: 6,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: true,
},
TautomerRule {
name: "1,3-C-to-O-any-bridge",
donor_elem: 6, bridge_elem: None, acceptor_elem: 8,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
TautomerRule {
name: "1,3-C-to-N-any-bridge",
donor_elem: 6, bridge_elem: None, acceptor_elem: 7,
donor_bridge_order: BondOrderMatch::Single,
bridge_acceptor_order: BondOrderMatch::Double,
prefer_forward: false,
},
];
fn h_assignment(mol: &Molecule) -> Vec<Option<u32>> {
(0..mol.atom_count())
.map(|i| mol.atom(AtomIdx(i as u32)).hydrogen_count.map(|h| h as u32))
.collect()
}
fn find_direct_aromatic_matches(mol: &Molecule) -> Vec<(AtomIdx, AtomIdx)> {
let mut pairs = Vec::new();
for (d, _) in mol.atoms() {
let da = mol.atom(d);
if !da.aromatic || da.element.atomic_number() != 7 {
continue;
}
if da.hydrogen_count.map_or(true, |h| h == 0) {
continue;
}
for (a, _) in mol.neighbors(d) {
let aa = mol.atom(a);
if aa.aromatic && (aa.element.atomic_number() == 7 || aa.element.atomic_number() == 8) {
pairs.push((d, a));
}
}
}
pairs
}
fn transfer_hydrogen_aromatic(mol: &Molecule, donor: AtomIdx, acceptor: AtomIdx) -> Option<Molecule> {
let donor_h = mol.atom(donor).hydrogen_count?;
if donor_h == 0 {
return None;
}
let acceptor_h = mol.atom(acceptor).hydrogen_count.unwrap_or(0);
let mut builder = MoleculeBuilder::new();
for i in 0..mol.atom_count() {
let idx = AtomIdx(i as u32);
let mut atom = mol.atom(idx).clone();
if idx == donor {
atom.hydrogen_count = Some(donor_h - 1);
} else if idx == acceptor {
atom.hydrogen_count = Some(acceptor_h + 1);
}
builder.add_atom(atom);
}
for i in 0..mol.bond_count() {
let bidx = BondIdx(i as u32);
let b = mol.bond(bidx);
builder.add_bond(b.atom1, b.atom2, b.order)
.expect("transfer_hydrogen_aromatic: bond from a valid molecule must be re-addable");
}
Some(builder.build())
}
fn clone_mol(mol: &Molecule) -> Molecule {
let mut builder = MoleculeBuilder::new();
for i in 0..mol.atom_count() {
builder.add_atom(mol.atom(AtomIdx(i as u32)).clone());
}
for i in 0..mol.bond_count() {
let b = mol.bond(BondIdx(i as u32));
builder.add_bond(b.atom1, b.atom2, b.order)
.expect("clone_mol: bond from a valid molecule must be re-addable");
}
builder.build()
}
const FNV1A_OFFSET: u64 = 0xcbf29ce484222325;
const FNV1A_PRIME: u64 = 0x100000001b3;
fn mol_fingerprint(mol: &Molecule) -> u64 {
let mut atoms: Vec<(u8, i8, u32)> = (0..mol.atom_count())
.map(|i| {
let idx = AtomIdx(i as u32);
let a = mol.atom(idx);
let bos: u32 = mol.neighbors(idx)
.map(|(_, bidx)| mol.bond(bidx).order.order_int() as u32)
.sum();
(a.element.atomic_number(), a.charge, bos)
})
.collect();
atoms.sort();
let mut hash = FNV1A_OFFSET;
for (an, ch, bos) in atoms {
hash ^= an as u64;
hash = hash.wrapping_mul(FNV1A_PRIME);
hash ^= (ch as u8 as u64).wrapping_add(128);
hash = hash.wrapping_mul(FNV1A_PRIME);
hash ^= bos as u64;
hash = hash.wrapping_mul(FNV1A_PRIME);
}
hash
}
fn find_matches(mol: &Molecule, rule: &TautomerRule) -> Vec<(AtomIdx, AtomIdx, AtomIdx)> {
let mut matches = Vec::new();
for i in 0..mol.atom_count() {
let d = AtomIdx(i as u32);
let donor_atom = mol.atom(d);
if donor_atom.element.atomic_number() != rule.donor_elem {
continue;
}
if implicit_hcount(mol, d) == 0 {
continue;
}
for (b, db_bidx) in mol.neighbors(d) {
if !rule.donor_bridge_order.matches(mol.bond(db_bidx).order) {
continue;
}
if let Some(br_elem) = rule.bridge_elem {
if mol.atom(b).element.atomic_number() != br_elem {
continue;
}
}
for (a, ba_bidx) in mol.neighbors(b) {
if a == d {
continue;
}
if !rule.bridge_acceptor_order.matches(mol.bond(ba_bidx).order) {
continue;
}
if mol.atom(a).element.atomic_number() != rule.acceptor_elem {
continue;
}
matches.push((d, b, a));
}
}
}
matches
}
fn transfer_hydrogen(
mol: &Molecule,
donor: AtomIdx,
bridge: AtomIdx,
acceptor: AtomIdx,
) -> Option<Molecule> {
let (db_bidx, _) = mol.bond_between(donor, bridge)?;
let (ba_bidx, _) = mol.bond_between(bridge, acceptor)?;
let mut builder = MoleculeBuilder::new();
for i in 0..mol.atom_count() {
let idx = AtomIdx(i as u32);
let mut atom = mol.atom(idx).clone();
if let Some(h) = atom.hydrogen_count {
if idx == donor {
if h == 0 { return None; }
atom.hydrogen_count = Some(h - 1);
} else if idx == acceptor {
atom.hydrogen_count = Some(h.saturating_add(1));
}
}
builder.add_atom(atom);
}
for i in 0..mol.bond_count() {
let bidx = BondIdx(i as u32);
let b = mol.bond(bidx);
let order = match bidx {
x if x == db_bidx => BondOrder::Double,
x if x == ba_bidx => BondOrder::Single,
_ => b.order,
};
builder.add_bond(b.atom1, b.atom2, order).ok()?;
}
Some(builder.build())
}
fn apply_first_match(mol: &Molecule, rule: &TautomerRule) -> Option<Molecule> {
find_matches(mol, rule)
.into_iter()
.find_map(|(d, b, a)| transfer_hydrogen(mol, d, b, a))
}
fn apply_all_matches(mol: &Molecule, rule: &TautomerRule) -> Vec<Molecule> {
find_matches(mol, rule)
.into_iter()
.filter_map(|(d, b, a)| transfer_hydrogen(mol, d, b, a))
.collect()
}
pub fn canonical_tautomer(mol: &Molecule) -> Molecule {
const MAX_ITER: usize = 16;
let mut current = clone_mol(mol);
let mut seen = HashSet::new();
seen.insert(mol_fingerprint(¤t));
for _ in 0..MAX_ITER {
let mut changed = false;
for rule in RULES.iter().filter(|r| r.prefer_forward) {
if let Some(next) = apply_first_match(¤t, rule) {
let fp = mol_fingerprint(&next);
if !seen.contains(&fp) {
seen.insert(fp);
current = next;
changed = true;
break;
}
}
}
if !changed {
break;
}
}
let mut candidates: Vec<Molecule> = vec![clone_mol(¤t)];
for (d, a) in find_direct_aromatic_matches(¤t) {
if let Some(t) = transfer_hydrogen_aromatic(¤t, d, a) {
candidates.push(t);
}
}
if candidates.len() > 1 {
current = candidates
.into_iter()
.min_by_key(|t| h_assignment(t))
.unwrap();
}
current
}
pub fn enumerate_tautomers(mol: &Molecule) -> Vec<Molecule> {
const MAX_TAUTOMERS: usize = 32;
let mut result = vec![clone_mol(mol)];
let mut seen = HashSet::new();
seen.insert(mol_fingerprint(mol));
let mut h_seen: HashSet<Vec<Option<u32>>> = HashSet::new();
h_seen.insert(h_assignment(mol));
let mut frontier = vec![clone_mol(mol)];
while !frontier.is_empty() && result.len() < MAX_TAUTOMERS {
let current = frontier.remove(0);
for rule in RULES.iter() {
for next in apply_all_matches(¤t, rule) {
let fp = mol_fingerprint(&next);
if !seen.contains(&fp) {
seen.insert(fp);
h_seen.insert(h_assignment(&next));
frontier.push(clone_mol(&next));
result.push(next);
if result.len() >= MAX_TAUTOMERS {
break;
}
}
}
if result.len() >= MAX_TAUTOMERS {
break;
}
}
for (d, a) in find_direct_aromatic_matches(¤t) {
if result.len() >= MAX_TAUTOMERS {
break;
}
if let Some(next) = transfer_hydrogen_aromatic(¤t, d, a) {
let ha = h_assignment(&next);
if !h_seen.contains(&ha) {
h_seen.insert(ha);
seen.insert(mol_fingerprint(&next));
frontier.push(clone_mol(&next));
result.push(next);
}
}
}
}
result
}
#[cfg(test)]
mod tests {
use super::*;
use chematic_smiles::parse;
#[test]
fn test_canonical_no_tautomers() {
let mol = parse("CCO").unwrap();
let t = canonical_tautomer(&mol);
assert_eq!(t.atom_count(), mol.atom_count());
}
#[test]
fn test_canonical_idempotent() {
let mol = parse("CC=O").unwrap(); let t1 = canonical_tautomer(&mol);
let t2 = canonical_tautomer(&t1);
assert_eq!(mol_fingerprint(&t1), mol_fingerprint(&t2));
}
#[test]
fn test_enumerate_single_no_match() {
let mol = parse("C").unwrap(); let tautomers = enumerate_tautomers(&mol);
assert_eq!(tautomers.len(), 1); }
#[test]
fn test_enumerate_cap() {
let mol = parse("CC(=O)CC(=O)C").unwrap(); let tautomers = enumerate_tautomers(&mol);
assert!(tautomers.len() <= 32);
assert!(tautomers.len() >= 1);
}
#[test]
fn test_enumerate_vinyl_alcohol() {
let mol = parse("OC=C").unwrap();
let tautomers = enumerate_tautomers(&mol);
assert!(tautomers.len() >= 2, "Expected >= 2 tautomers for vinyl alcohol, got {}", tautomers.len());
}
#[test]
fn test_canonical_amide_unchanged() {
let mol = parse("CC(=O)N").unwrap();
let t = canonical_tautomer(&mol);
assert_eq!(mol_fingerprint(&t), mol_fingerprint(&mol));
}
#[test]
fn test_canonical_acetylacetone_stable() {
let mol = parse("CC(=O)CC(=O)C").unwrap();
let t = canonical_tautomer(&mol);
assert!(t.atom_count() > 0);
}
#[test]
fn test_enumerate_includes_original() {
let mol = parse("CC=O").unwrap();
let tautomers = enumerate_tautomers(&mol);
assert_eq!(mol_fingerprint(&tautomers[0]), mol_fingerprint(&mol));
}
#[test]
fn test_canonical_keto_unchanged() {
let mol = parse("CC=O").unwrap();
let t = canonical_tautomer(&mol);
assert_eq!(mol_fingerprint(&t), mol_fingerprint(&mol));
}
#[test]
fn test_canonical_acetylacetone_enol() {
let enol = parse("CC(O)=CC(=O)C").unwrap();
let keto = parse("CC(=O)CC(=O)C").unwrap();
let t_enol = canonical_tautomer(&enol);
let t_keto = canonical_tautomer(&keto);
assert!(t_enol.atom_count() > 0);
assert!(t_keto.atom_count() > 0);
}
#[test]
fn test_enumerate_pyrazole_12_shift() {
let mol = parse("c1cc[nH]n1").unwrap();
let tautomers = enumerate_tautomers(&mol);
assert!(
tautomers.len() >= 2,
"Expected >= 2 tautomers for pyrazole, got {}",
tautomers.len()
);
}
#[test]
fn test_canonical_pyrazole_normalization() {
use chematic_smiles::canonical_smiles;
let n1h = parse("c1cc[nH]n1").unwrap();
let tautomers = enumerate_tautomers(&n1h);
let n2h = tautomers
.iter()
.find(|t| h_assignment(t) != h_assignment(&n1h))
.expect("enumerate_tautomers should produce N2H tautomer of pyrazole");
assert_eq!(
canonical_smiles(&canonical_tautomer(&n1h)),
canonical_smiles(&canonical_tautomer(n2h)),
"canonical_tautomer should normalize N1H and N2H to the same form"
);
}
}