#![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.is_none_or(|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
&& 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()
}
#[derive(Debug, Clone)]
pub struct TautomerConfig {
pub max_iter: usize,
pub max_tautomers: usize,
pub enabled_rules: Vec<usize>,
}
impl Default for TautomerConfig {
fn default() -> Self {
Self {
max_iter: 16,
max_tautomers: 32,
enabled_rules: Vec::new(),
}
}
}
impl TautomerConfig {
pub fn rule_count() -> usize {
RULES.len()
}
pub fn rule_names() -> Vec<&'static str> {
RULES.iter().map(|r| r.name).collect()
}
pub fn keto_enol_only() -> Self {
Self {
enabled_rules: vec![0],
..Self::default()
}
}
}
fn active_rules(config: &TautomerConfig) -> Vec<&'static TautomerRule> {
RULES
.iter()
.enumerate()
.filter_map(|(i, r)| {
if config.enabled_rules.is_empty() || config.enabled_rules.contains(&i) {
Some(r)
} else {
None
}
})
.collect()
}
pub fn canonical_tautomer(mol: &Molecule) -> Molecule {
canonical_tautomer_with_config(mol, &TautomerConfig::default())
}
pub fn canonical_tautomer_with_config(mol: &Molecule, config: &TautomerConfig) -> Molecule {
let mut current = clone_mol(mol);
let mut seen = HashSet::new();
seen.insert(mol_fingerprint(¤t));
for _ in 0..config.max_iter {
let mut changed = false;
for rule in active_rules(config)
.into_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(h_assignment).unwrap();
}
current
}
pub fn enumerate_tautomers(mol: &Molecule) -> Vec<Molecule> {
enumerate_tautomers_with_config(mol, &TautomerConfig::default())
}
pub fn enumerate_tautomers_with_config(mol: &Molecule, config: &TautomerConfig) -> Vec<Molecule> {
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() < config.max_tautomers {
let current = frontier.remove(0);
for rule in active_rules(config).into_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() >= config.max_tautomers {
break;
}
}
}
if result.len() >= config.max_tautomers {
break;
}
}
for (d, a) in find_direct_aromatic_matches(¤t) {
if result.len() >= config.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.is_empty());
}
#[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"
);
}
#[test]
fn test_config_default_same_as_no_config() {
use chematic_smiles::canonical_smiles;
let mol = parse("OC=C").unwrap(); let a = canonical_tautomer(&mol);
let b = canonical_tautomer_with_config(&mol, &TautomerConfig::default());
assert_eq!(
canonical_smiles(&a),
canonical_smiles(&b),
"default config should match canonical_tautomer"
);
}
#[test]
fn test_config_max_iter_one_limits_convergence() {
let mol = parse("OC=C").unwrap();
let config = TautomerConfig {
max_iter: 1,
..TautomerConfig::default()
};
let _ = canonical_tautomer_with_config(&mol, &config);
}
#[test]
fn test_config_max_tautomers_caps_enumerate() {
let mol = parse("CC(=O)CC(=O)C").unwrap(); let config = TautomerConfig {
max_tautomers: 2,
..TautomerConfig::default()
};
let tautomers = enumerate_tautomers_with_config(&mol, &config);
assert_eq!(
tautomers.len(),
2,
"max_tautomers=2 should return exactly 2"
);
}
#[test]
fn test_config_enabled_rules_subset() {
let mol = parse("OC=C").unwrap();
let config = TautomerConfig::keto_enol_only();
let result = canonical_tautomer_with_config(&mol, &config);
assert!(result.atom_count() > 0);
}
#[test]
fn test_config_empty_enabled_rules_equals_all() {
use chematic_smiles::canonical_smiles;
let mol = parse("OC=C").unwrap();
let all = canonical_tautomer_with_config(&mol, &TautomerConfig::default());
let explicit_empty = canonical_tautomer_with_config(
&mol,
&TautomerConfig {
enabled_rules: vec![],
..TautomerConfig::default()
},
);
assert_eq!(
canonical_smiles(&all),
canonical_smiles(&explicit_empty),
"empty enabled_rules should equal all rules"
);
}
#[test]
fn test_rule_count_and_names() {
let count = TautomerConfig::rule_count();
let names = TautomerConfig::rule_names();
assert!(count > 0);
assert_eq!(names.len(), count);
assert!(!names[0].is_empty());
}
}