use std::collections::{BTreeMap, BTreeSet, VecDeque};
use crate::{
AtomId, BondDirection, BondId, BondOrder, Molecule, RingFindingError, RingInfo,
ValenceAssignment, ValenceError, ValenceModel, canon_rank::FragmentRankScope,
canon_rank::rank_fragment_atoms_for_kekulize, molecule::TopologyBlock,
read_parts::MoleculeReadParts,
};
#[derive(Debug, Clone, PartialEq, Eq, thiserror::Error)]
pub enum KekulizeError {
#[error("can't kekulize mol. unkekulized atoms: {atoms:?}")]
UnkekulizedAtoms { atoms: Vec<AtomId> },
#[error("non-ring atom {atom} marked aromatic")]
NonRingAromaticAtom { atom: AtomId },
#[error("incomplete RDKit kekulize port for {branch}: {reason}")]
ProtocolDebt {
branch: &'static str,
reason: &'static str,
},
#[error("kekulization changed valence on atom {atom}: {before}!={after}")]
ValenceChanged {
atom: AtomId,
before: i32,
after: i32,
},
#[error("kekulize fragment bitset size mismatch: atoms={atoms}, bonds={bonds}")]
FragmentBitsetSizeMismatch { atoms: usize, bonds: usize },
#[error("canonical rank {kind} symbol size mismatch: expected={expected}, actual={actual}")]
CanonicalRankSymbolSizeMismatch {
kind: &'static str,
expected: usize,
actual: usize,
},
#[error(transparent)]
RingFinding(#[from] RingFindingError),
#[error(transparent)]
Valence(#[from] ValenceError),
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub(crate) struct KekulizeAssignment {
bond_orders: Vec<Option<BondOrder>>,
bond_aromatic_flags: Vec<Option<bool>>,
bond_directions: Vec<Option<BondDirection>>,
atom_aromatic_flags: Vec<Option<bool>>,
atom_explicit_hydrogens: Vec<Option<u8>>,
atom_no_implicit: Vec<Option<bool>>,
discovered_rings: Option<RingInfo>,
}
struct KekulizeCandidateState {
all_atoms: Vec<AtomId>,
candidate_atoms: BTreeSet<AtomId>,
aromatic_edges: Vec<(BondId, AtomId, AtomId)>,
questions: Vec<AtomId>,
done: Vec<AtomId>,
}
#[derive(Debug, Clone, PartialEq, Eq)]
struct KekulizeRing {
atoms: Vec<AtomId>,
bonds: Vec<BondId>,
source_ring: usize,
}
impl KekulizeAssignment {
fn empty(num_atoms: usize, num_bonds: usize) -> Self {
Self::with_sizes(num_atoms, num_bonds)
}
fn with_sizes(num_atoms: usize, num_bonds: usize) -> Self {
Self {
bond_orders: vec![None; num_bonds],
bond_aromatic_flags: vec![None; num_bonds],
bond_directions: vec![None; num_bonds],
atom_aromatic_flags: vec![None; num_atoms],
atom_explicit_hydrogens: vec![None; num_atoms],
atom_no_implicit: vec![None; num_atoms],
discovered_rings: None,
}
}
pub(crate) fn bond_order(&self, bond: BondId) -> Option<BondOrder> {
self.bond_orders[bond.index()]
}
pub(crate) fn bond_aromatic_flag(&self, bond: BondId) -> Option<bool> {
self.bond_aromatic_flags[bond.index()]
}
pub(crate) fn bond_direction(&self, bond: BondId) -> Option<BondDirection> {
self.bond_directions[bond.index()]
}
pub(crate) fn atom_aromatic_flag(&self, atom: AtomId) -> Option<bool> {
self.atom_aromatic_flags[atom.index()]
}
pub(crate) fn atom_explicit_hydrogens(&self, atom: AtomId) -> Option<u8> {
self.atom_explicit_hydrogens[atom.index()]
}
pub(crate) fn atom_no_implicit(&self, atom: AtomId) -> Option<bool> {
self.atom_no_implicit[atom.index()]
}
pub(crate) fn discovered_rings(&self) -> Option<&RingInfo> {
self.discovered_rings.as_ref()
}
#[cfg(test)]
fn is_empty(&self) -> bool {
self.bond_orders.iter().all(Option::is_none)
&& self.bond_aromatic_flags.iter().all(Option::is_none)
&& self.bond_directions.iter().all(Option::is_none)
&& self.atom_aromatic_flags.iter().all(Option::is_none)
&& self.atom_explicit_hydrogens.iter().all(Option::is_none)
&& self.atom_no_implicit.iter().all(Option::is_none)
}
}
pub(crate) fn apply_kekulize_assignment(
topology: &mut TopologyBlock,
assignment: &KekulizeAssignment,
) -> bool {
let mut changed = false;
for bond in &mut topology.bonds {
if let Some(next) = assignment.bond_order(bond.id()) {
if bond.order() != next {
bond.set_order(next);
changed = true;
}
}
if let Some(next) = assignment.bond_aromatic_flag(bond.id()) {
if bond.is_aromatic() != next {
bond.set_aromatic(next);
changed = true;
}
}
if let Some(next) = assignment.bond_direction(bond.id()) {
if bond.direction() != next {
bond.set_direction(next);
changed = true;
}
}
}
for atom in &mut topology.atoms {
if let Some(next) = assignment.atom_aromatic_flag(atom.id()) {
if atom.is_aromatic() != next {
atom.set_aromatic(next);
changed = true;
}
}
if let Some(next) = assignment.atom_explicit_hydrogens(atom.id()) {
if atom.explicit_hydrogens() != next {
atom.set_explicit_hydrogens(next);
changed = true;
}
}
if let Some(next) = assignment.atom_no_implicit(atom.id()) {
if atom.no_implicit() != next {
atom.set_no_implicit(next);
changed = true;
}
}
}
changed
}
fn topology_is_aromatic_atom(topology: &TopologyBlock, atom: AtomId) -> bool {
let atom_data = &topology.atoms[atom.index()];
if atom_data.is_aromatic() {
return true;
}
topology.bonds.iter().any(|bond| {
(bond.begin() == atom || bond.end() == atom)
&& (bond.is_aromatic() || bond.order() == BondOrder::Aromatic)
})
}
fn kekulize_is_aromatic_bond(bond: &crate::Bond) -> bool {
bond.is_aromatic() || bond.order() == BondOrder::Aromatic
}
pub(crate) fn kekulize_assignment(
molecule: &Molecule,
rings: Option<&RingInfo>,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<KekulizeAssignment, KekulizeError> {
kekulize_assignment_from_parts_with_valence(
molecule.topology_block(),
molecule.stereo_groups(),
rings,
molecule.derived_cache().valence.as_ref(),
clear_aromatic_flags,
canonical,
max_backtracks,
)
}
pub(crate) fn kekulize_assignment_from_parts(
topology: &TopologyBlock,
stereo_groups: &[crate::stereo::StereoGroup],
rings: Option<&RingInfo>,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<KekulizeAssignment, KekulizeError> {
kekulize_assignment_from_parts_with_valence(
topology,
stereo_groups,
rings,
None,
clear_aromatic_flags,
canonical,
max_backtracks,
)
}
fn kekulize_assignment_from_parts_with_valence(
topology: &TopologyBlock,
stereo_groups: &[crate::stereo::StereoGroup],
rings: Option<&RingInfo>,
cached_valence: Option<&ValenceAssignment>,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<KekulizeAssignment, KekulizeError> {
let atoms_to_use = vec![true; topology.atoms.len()];
let bonds_to_use = vec![true; topology.bonds.len()];
kekulize_fragment_assignment(
topology,
stereo_groups,
rings,
cached_valence,
&atoms_to_use,
bonds_to_use,
clear_aromatic_flags,
canonical,
max_backtracks,
)
}
pub(crate) fn kekulize_assignment_from_read_parts(
read: MoleculeReadParts<'_>,
rings: Option<&RingInfo>,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<KekulizeAssignment, KekulizeError> {
kekulize_assignment_from_read_parts_with_valence(
read,
rings,
None,
clear_aromatic_flags,
canonical,
max_backtracks,
)
}
pub(crate) fn kekulize_assignment_from_read_parts_with_valence(
read: MoleculeReadParts<'_>,
rings: Option<&RingInfo>,
cached_valence: Option<&ValenceAssignment>,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<KekulizeAssignment, KekulizeError> {
kekulize_assignment_from_parts_with_valence(
read.topology(),
read.stereo_groups(),
rings,
cached_valence,
clear_aromatic_flags,
canonical,
max_backtracks,
)
}
pub(crate) fn kekulize_fragment_assignment(
topology: &TopologyBlock,
stereo_groups: &[crate::stereo::StereoGroup],
rings: Option<&RingInfo>,
cached_valence: Option<&ValenceAssignment>,
atoms_to_use: &[bool],
mut bonds_to_use: Vec<bool>,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<KekulizeAssignment, KekulizeError> {
let trace_kekulize = std::env::var_os("COSMOLKIT_TRACE_KEKULIZE").is_some();
if atoms_to_use.len() != topology.atoms.len() || bonds_to_use.len() != topology.bonds.len() {
return Err(KekulizeError::FragmentBitsetSizeMismatch {
atoms: atoms_to_use.len(),
bonds: bonds_to_use.len(),
});
}
let mut assignment = KekulizeAssignment::with_sizes(topology.atoms.len(), topology.bonds.len());
if !atoms_to_use.iter().any(|selected| *selected) {
return Ok(assignment);
}
let mut found_aromatic = false;
for bond in &topology.bonds {
if bonds_to_use[bond.id().index()] {
if bond
.query()
.is_some_and(crate::valence::has_bond_type_query)
{
bonds_to_use[bond.id().index()] = false;
} else if kekulize_is_aromatic_bond(bond) {
found_aromatic = true;
}
}
}
let pre_kek_valence = if let Some(valence) = cached_valence {
valence.clone()
} else {
crate::valence::assign_valence_with_options_from_parts(
&topology.atoms,
&topology.bonds,
&topology.adjacency,
ValenceModel::RdkitLike,
false,
)?
};
if topology
.atoms
.iter()
.any(|atom| topology_is_aromatic_atom(topology, atom.id()))
{
found_aromatic = true;
}
if !found_aromatic {
return Ok(assignment);
}
let discovered_rings = if rings.is_none() {
Some(crate::rings::find_sssr_from_parts(
topology.atoms.len(),
&topology.bonds,
&topology.adjacency,
)?)
} else {
None
};
let rings = rings
.or(discovered_rings.as_ref())
.expect("Kekulize has either cached or newly discovered ring information");
let atom_ranks = if canonical {
rank_fragment_atoms_for_kekulize(
topology,
stereo_groups,
FragmentRankScope::new(atoms_to_use, &bonds_to_use),
)?
} else {
(0..topology.atoms.len()).collect::<Vec<_>>()
};
if trace_kekulize {
eprintln!("kekulize-fragment atomRanks={atom_ranks:?}");
}
let mut wedged_atoms = BTreeSet::new();
for bond in &topology.bonds {
if bonds_to_use[bond.id().index()]
&& matches!(
bond.direction(),
BondDirection::BeginWedge | BondDirection::BeginDash
)
{
wedged_atoms.insert(bond.begin());
}
}
let kekulize_rings =
ordered_kekulize_rings(topology, rings, atoms_to_use, &bonds_to_use, &wedged_atoms);
let ring_systems = fused_ring_systems(&kekulize_rings);
if trace_kekulize {
eprintln!(
"kekulize-rings arings={:?}",
kekulize_rings
.iter()
.map(|ring| ring
.atoms
.iter()
.map(|atom| atom.index())
.collect::<Vec<_>>())
.collect::<Vec<_>>()
);
eprintln!(
"kekulize-rings brings={:?}",
kekulize_rings
.iter()
.map(|ring| ring
.bonds
.iter()
.map(|bond| bond.index())
.collect::<Vec<_>>())
.collect::<Vec<_>>()
);
eprintln!("kekulize-rings fused={ring_systems:?}");
}
for system in ring_systems {
let ring_bond_set = system
.iter()
.flat_map(|&ring_idx| kekulize_rings[ring_idx].bonds.iter().copied())
.collect::<BTreeSet<_>>();
let ring_atom_set = system
.iter()
.flat_map(|&ring_idx| kekulize_rings[ring_idx].atoms.iter().copied())
.collect::<BTreeSet<_>>();
if ring_bond_set.iter().any(|bond| {
let bond = &topology.bonds[bond.index()];
kekulize_is_aromatic_bond(bond)
}) || ring_atom_set
.iter()
.any(|atom| topology_is_aromatic_atom(topology, *atom))
{
kekulize_fused_assignment(
topology,
rings,
&kekulize_rings,
&system,
&pre_kek_valence,
&atom_ranks,
max_backtracks,
&mut assignment,
)?;
}
}
if clear_aromatic_flags {
let full_molecule_scope = atoms_to_use.iter().all(|selected| *selected)
&& bonds_to_use.iter().all(|selected| *selected);
for bond in &topology.bonds {
if bonds_to_use[bond.id().index()] && bond.is_aromatic() {
assignment.bond_aromatic_flags[bond.id().index()] = Some(false);
}
}
for atom in &topology.atoms {
if atoms_to_use[atom.id().index()] && atom.is_aromatic() {
if full_molecule_scope && rings.num_atom_rings(atom.id()) == 0 {
return Err(KekulizeError::NonRingAromaticAtom { atom: atom.id() });
}
assignment.atom_aromatic_flags[atom.id().index()] = Some(false);
if matches!(atom.atomic_number(), 7 | 15)
&& atom.formal_charge() == 0
&& atom.explicit_hydrogens() == 1
{
assignment.atom_no_implicit[atom.id().index()] = Some(false);
assignment.atom_explicit_hydrogens[atom.id().index()] = Some(0);
}
}
}
}
let mut topology_after = topology.clone();
apply_kekulize_assignment(&mut topology_after, &assignment);
let mut post_kek_valence = pre_kek_valence.clone();
for atom in &topology_after.atoms {
let idx = atom.id().index();
if assignment.atom_explicit_hydrogens[idx].is_some()
|| assignment.atom_no_implicit[idx].is_some()
{
let (explicit_valence, implicit_hydrogens) =
crate::valence::assign_valence_state_for_atom_from_parts(
&topology_after.atoms,
&topology_after.bonds,
&topology_after.adjacency,
atom.id(),
false,
)?;
post_kek_valence.explicit_valence[idx] = explicit_valence;
post_kek_valence.implicit_hydrogens[idx] = implicit_hydrogens;
}
}
for atom in &topology.atoms {
let idx = atom.id().index();
let before =
pre_kek_valence.explicit_valence[idx] + pre_kek_valence.implicit_hydrogens[idx];
let after =
post_kek_valence.explicit_valence[idx] + post_kek_valence.implicit_hydrogens[idx];
if before != after {
return Err(KekulizeError::ValenceChanged {
atom: atom.id(),
before,
after,
});
}
}
assignment.discovered_rings = discovered_rings;
Ok(assignment)
}
fn kekulize_fused_assignment(
topology: &TopologyBlock,
rings: &RingInfo,
kekulize_rings: &[KekulizeRing],
fused_ring_indices: &[usize],
valence: &ValenceAssignment,
atom_ranks: &[usize],
max_backtracks: usize,
assignment: &mut KekulizeAssignment,
) -> Result<(), KekulizeError> {
let state = mark_double_bond_candidates(
topology,
rings,
kekulize_rings,
fused_ring_indices,
valence,
assignment,
)?;
let matched_bonds = if let Some(matched_bonds) = kekulize_worker_matching(
topology,
&state,
&state.done,
atom_ranks,
max_backtracks,
&BTreeSet::new(),
) {
matched_bonds
} else if let Some(matched_bonds) = permute_dummies_and_match(
topology,
&state,
&state.questions,
atom_ranks,
max_backtracks,
) {
matched_bonds
} else {
return Err(KekulizeError::UnkekulizedAtoms {
atoms: state.candidate_atoms.iter().copied().collect(),
});
};
for (bond, _, _) in state.aromatic_edges {
assignment.bond_orders[bond.index()] = Some(if matched_bonds.contains(&bond) {
BondOrder::Double
} else {
BondOrder::Single
});
if matched_bonds.contains(&bond)
&& topology.bonds[bond.index()].direction() != BondDirection::None
{
assignment.bond_directions[bond.index()] = Some(BondDirection::None);
}
}
Ok(())
}
fn mark_double_bond_candidates(
topology: &TopologyBlock,
rings: &RingInfo,
kekulize_rings: &[KekulizeRing],
fused_ring_indices: &[usize],
valence: &ValenceAssignment,
assignment: &mut KekulizeAssignment,
) -> Result<KekulizeCandidateState, KekulizeError> {
let mut all_atoms = Vec::new();
let mut seen_atoms = BTreeSet::new();
for &ring_idx in fused_ring_indices {
for &atom_id in &kekulize_rings[ring_idx].atoms {
if seen_atoms.insert(atom_id) {
all_atoms.push(atom_id);
}
}
}
let all_atom_set = all_atoms.iter().copied().collect::<BTreeSet<_>>();
if !all_atoms.iter().any(|atom| {
topology.atoms[atom.index()].atomic_number() == 0
|| topology_is_aromatic_atom(topology, *atom)
}) {
return Ok(KekulizeCandidateState {
all_atoms,
candidate_atoms: BTreeSet::new(),
aromatic_edges: Vec::new(),
questions: Vec::new(),
done: Vec::new(),
});
}
let fused_source_ring_set = fused_ring_indices
.iter()
.map(|&ring_idx| kekulize_rings[ring_idx].source_ring)
.collect::<BTreeSet<_>>();
let mut is_ring_not_cand = BTreeSet::new();
for &ring_idx in fused_ring_indices {
let mut ring_is_candidate = false;
let source_ring = kekulize_rings[ring_idx].source_ring;
for atom in &rings.atom_rings()[source_ring] {
if topology_is_aromatic_atom(topology, *atom) && rings.num_atom_rings(*atom) == 1 {
ring_is_candidate = true;
break;
}
}
if !ring_is_candidate {
is_ring_not_cand.insert(source_ring);
}
}
let mut candidate_atoms = BTreeSet::<AtomId>::new();
let mut aromatic_edges = Vec::<(BondId, AtomId, AtomId)>::new();
let mut seen_make_single = BTreeSet::<BondId>::new();
let mut questions = Vec::<AtomId>::new();
let mut done = Vec::<AtomId>::new();
for &atom_id in &all_atoms {
let atom = &topology.atoms[atom_id.index()];
if atom.atomic_number() != 0 && !topology_is_aromatic_atom(topology, atom_id) {
done.push(atom_id);
continue;
}
let mut sbo = 0i32;
let mut n_to_ignore = 0usize;
let mut non_ar_non_dummy_nbr = 0usize;
for bond in topology
.bonds
.iter()
.filter(|bond| bond.begin() == atom_id || bond.end() == atom_id)
{
let other = if bond.begin() == atom_id {
bond.end()
} else {
bond.begin()
};
let other_atom = &topology.atoms[other.index()];
if all_atom_set.contains(&other)
&& other_atom.atomic_number() != 0
&& !topology_is_aromatic_atom(topology, other)
{
non_ar_non_dummy_nbr += 1;
}
if kekulize_is_aromatic_bond(bond)
&& matches!(
bond.order(),
BondOrder::Single | BondOrder::Double | BondOrder::Aromatic
)
{
sbo += 1;
if seen_make_single.insert(bond.id()) {
assignment.bond_orders[bond.id().index()] = Some(BondOrder::Single);
aromatic_edges.push((bond.id(), bond.begin(), bond.end()));
}
} else {
let bond_contrib =
crate::valence::bond_valence_contrib(bond, atom_id)?.round() as i32;
sbo += bond_contrib;
if bond_contrib == 0 {
n_to_ignore += 1;
}
}
}
let num_atom_rings = rings.num_atom_rings(atom_id);
let num_non_cand_rings = rings
.atom_members(atom_id)
.iter()
.filter(|ring_idx| {
fused_source_ring_set.contains(ring_idx) && is_ring_not_cand.contains(ring_idx)
})
.count();
if atom.atomic_number() == 0
&& non_ar_non_dummy_nbr < num_atom_rings
&& num_non_cand_rings < num_atom_rings
{
candidate_atoms.insert(atom_id);
questions.push(atom_id);
continue;
}
if atom.atomic_number() == 0 {
continue;
}
sbo += i32::from(atom.explicit_hydrogens()) + valence.implicit_hydrogens[atom_id.index()];
let mut dv = rdkit_kekulize_default_valence(atom.atomic_number())?;
let mut chrg = atom.formal_charge();
if rdkit_is_early_atom(atom.atomic_number()) {
chrg = -chrg;
}
if atom.atomic_number() == 6 && chrg > 0 {
chrg = -chrg;
}
dv += i32::from(chrg);
let tbo =
valence.explicit_valence[atom_id.index()] + valence.implicit_hydrogens[atom_id.index()];
let n_radicals = i32::from(atom.radical_electrons());
let total_degree = i32::try_from(
topology
.bonds
.iter()
.filter(|bond| bond.begin() == atom_id || bond.end() == atom_id)
.count(),
)
.unwrap_or(i32::MAX)
+ valence.implicit_hydrogens[atom_id.index()]
- i32::try_from(n_to_ignore).unwrap_or(i32::MAX);
let valence_list = crate::rdkit_valence_list(atom.atomic_number())?.ok_or(
ValenceError::UnsupportedBranch {
reason: "kekulize valence list atomic number out of range",
},
)?;
let mut vi = 1usize;
while tbo > dv && vi < valence_list.len() && valence_list[vi] > 0 {
dv = valence_list[vi] + i32::from(chrg);
vi += 1;
}
if tbo == 5
&& sbo == 4
&& dv == 3
&& total_degree == 3
&& n_radicals == 0
&& chrg == 0
&& atom.explicit_hydrogens() == 0
&& valence.implicit_hydrogens[atom_id.index()] == 0
&& matches!(atom.atomic_number(), 7 | 15 | 33)
{
dv = 5;
}
if total_degree + n_radicals >= dv {
continue;
}
if dv == (sbo + 1 + n_radicals)
|| (n_radicals == 0 && atom.no_implicit() && dv == (sbo + 2))
{
candidate_atoms.insert(atom_id);
}
}
Ok(KekulizeCandidateState {
all_atoms,
candidate_atoms,
aromatic_edges,
questions,
done,
})
}
fn kekulize_worker_matching(
topology: &TopologyBlock,
state: &KekulizeCandidateState,
initial_done: &[AtomId],
atom_ranks: &[usize],
max_backtracks: usize,
switched_off: &BTreeSet<AtomId>,
) -> Option<BTreeSet<BondId>> {
let trace_kekulize = std::env::var_os("COSMOLKIT_TRACE_KEKULIZE").is_some();
let all_atom_set = state.all_atoms.iter().copied().collect::<BTreeSet<_>>();
let mut aromatic_edge_ids = vec![false; topology.bonds.len()];
for (bond, _, _) in &state.aromatic_edges {
aromatic_edge_ids[bond.index()] = true;
}
let mut adjacency = vec![Vec::<(AtomId, BondId)>::new(); topology.atoms.len()];
for bond in &topology.bonds {
let begin = bond.begin();
let end = bond.end();
if all_atom_set.contains(&begin) && all_atom_set.contains(&end) {
adjacency[begin.index()].push((end, bond.id()));
adjacency[end.index()].push((begin, bond.id()));
}
}
for neighbors in &mut adjacency {
neighbors.sort_by_key(|(atom, _)| (atom_ranks[atom.index()], atom.index()));
}
let sorted_atoms = worker_sorted_atoms(topology, &state.all_atoms, &all_atom_set, atom_ranks);
if trace_kekulize {
eprintln!(
"kekulize-worker sortedAtms={:?}",
sorted_atoms
.iter()
.map(|atom| atom.index())
.collect::<Vec<_>>()
);
eprintln!(
"kekulize-worker allAtms={:?} done={:?} candidates={:?} questions={:?}",
state
.all_atoms
.iter()
.map(|atom| atom.index())
.collect::<Vec<_>>(),
initial_done
.iter()
.map(|atom| atom.index())
.collect::<Vec<_>>(),
state
.candidate_atoms
.iter()
.map(|atom| atom.index())
.collect::<Vec<_>>(),
state
.questions
.iter()
.map(|atom| atom.index())
.collect::<Vec<_>>()
);
}
let mut d_bnd_cands = vec![false; topology.atoms.len()];
for atom in &state.candidate_atoms {
if !switched_off.contains(atom) {
d_bnd_cands[atom.index()] = true;
}
}
let mut d_bnd_adds = vec![false; topology.bonds.len()];
let mut local_bonds_added = vec![false; topology.bonds.len()];
let mut done = initial_done.to_vec();
let mut done_flags = vec![false; topology.atoms.len()];
for atom in &done {
done_flags[atom.index()] = true;
}
let mut astack = VecDeque::<AtomId>::new();
let mut astack_flags = vec![false; topology.atoms.len()];
let mut options = BTreeMap::<AtomId, VecDeque<(AtomId, BondId)>>::new();
let mut last_opt = None::<AtomId>;
let mut btmoves = Vec::<AtomId>::new();
let mut matched = BTreeSet::<BondId>::new();
let mut num_bt = 0usize;
while done.len() < sorted_atoms.len() || !astack.is_empty() {
let curr = if let Some(curr) = astack.pop_front() {
astack_flags[curr.index()] = false;
curr
} else {
sorted_atoms
.iter()
.copied()
.find(|atom| !done_flags[atom.index()])
.expect("starting point not found")
};
done.push(curr);
done_flags[curr.index()] = true;
let c_cand = d_bnd_cands[curr.index()];
let mut opts = if let Some(saved) = options.get(&curr) {
saved.clone()
} else {
let mut lstack = VecDeque::new();
let mut opts_v = Vec::<(AtomId, BondId)>::new();
let mut wedged_opts_v = Vec::<(AtomId, BondId)>::new();
for &(nbr_idx, nbr_bond) in &adjacency[curr.index()] {
if done_flags[nbr_idx.index()] {
continue;
}
if !astack_flags[nbr_idx.index()] {
lstack.push_back(nbr_idx);
astack_flags[nbr_idx.index()] = true;
}
let bond = &topology.bonds[nbr_bond.index()];
if c_cand
&& d_bnd_cands[nbr_idx.index()]
&& (aromatic_edge_ids[nbr_bond.index()]
|| kekulize_is_aromatic_bond(bond)
|| topology.atoms[curr.index()].atomic_number() == 0
|| topology.atoms[nbr_idx.index()].atomic_number() == 0)
{
if matches!(
bond.direction(),
BondDirection::BeginWedge | BondDirection::BeginDash
) {
wedged_opts_v.push((nbr_idx, nbr_bond));
} else {
opts_v.push((nbr_idx, nbr_bond));
}
}
}
let mut computed = VecDeque::new();
computed.extend(opts_v);
computed.extend(wedged_opts_v);
astack.extend(lstack);
computed
};
if trace_kekulize {
eprintln!(
"kekulize-worker step curr={} cCand={} opts={:?} astack={:?} done={:?}",
curr.index(),
c_cand,
opts.iter()
.map(|(atom, bond)| (atom.index(), bond.index()))
.collect::<Vec<_>>(),
astack.iter().map(|atom| atom.index()).collect::<Vec<_>>(),
done.iter().map(|atom| atom.index()).collect::<Vec<_>>()
);
}
if c_cand {
if let Some((ncnd, bond_id)) = opts.pop_front() {
d_bnd_cands[curr.index()] = false;
d_bnd_cands[ncnd.index()] = false;
d_bnd_adds[bond_id.index()] = true;
local_bonds_added[bond_id.index()] = true;
matched.insert(bond_id);
if trace_kekulize {
eprintln!(
"kekulize-worker add curr={} ncnd={} bond={}",
curr.index(),
ncnd.index(),
bond_id.index()
);
}
if options.contains_key(&curr) {
if opts.is_empty() {
options.remove(&curr);
let _ = btmoves.pop();
last_opt = btmoves.last().copied();
} else {
options.insert(curr, opts);
}
} else if !opts.is_empty() {
last_opt = Some(curr);
btmoves.push(curr);
options.insert(curr, opts);
}
} else if topology.atoms[curr.index()].atomic_number() != 0 {
if let Some(last_opt_atom) = last_opt {
if num_bt < max_backtracks {
if trace_kekulize {
eprintln!(
"kekulize-worker backtrack curr={} lastOpt={} numBT={}",
curr.index(),
last_opt_atom.index(),
num_bt
);
}
back_track(
topology,
&mut options,
last_opt_atom,
&mut done,
&mut done_flags,
&mut astack,
&mut astack_flags,
&mut d_bnd_cands,
&mut d_bnd_adds,
&mut matched,
);
num_bt += 1;
} else {
for (bond_idx, was_added) in local_bonds_added.iter().enumerate() {
if *was_added {
matched.remove(&BondId::new(bond_idx));
}
}
return None;
}
} else {
for (bond_idx, was_added) in local_bonds_added.iter().enumerate() {
if *was_added {
matched.remove(&BondId::new(bond_idx));
}
}
return None;
}
}
}
}
Some(matched)
}
fn permute_dummies_and_match(
topology: &TopologyBlock,
state: &KekulizeCandidateState,
questions: &[AtomId],
atom_ranks: &[usize],
max_backtracks: usize,
) -> Option<BTreeSet<BondId>> {
if questions.is_empty() {
return None;
}
for switched_off in question_switch_masks(questions) {
if let Some(matched) = kekulize_worker_matching(
topology,
state,
&[],
atom_ranks,
max_backtracks,
&switched_off,
) {
return Some(matched);
}
}
None
}
fn back_track(
topology: &TopologyBlock,
_options: &mut BTreeMap<AtomId, VecDeque<(AtomId, BondId)>>,
last_opt: AtomId,
done: &mut Vec<AtomId>,
done_flags: &mut [bool],
aqueue: &mut VecDeque<AtomId>,
astack_flags: &mut [bool],
d_bnd_cands: &mut [bool],
d_bnd_adds: &mut [bool],
matched: &mut BTreeSet<BondId>,
) {
let split = done
.iter()
.position(|atom| *atom == last_opt)
.expect("lastOpt must exist in done");
let tdone = done[..split].to_vec();
for atom in done[split..].iter().rev() {
aqueue.push_front(*atom);
}
for atom in done.iter().skip(split) {
done_flags[atom.index()] = false;
}
for atom in done[split..].iter().rev() {
astack_flags[atom.index()] = true;
}
for bond in &topology.bonds {
let bi = bond.id().index();
if d_bnd_adds[bi] {
let aid1 = bond.begin();
let aid2 = bond.end();
if !tdone.contains(&aid1) && !tdone.contains(&aid2) {
d_bnd_adds[bi] = false;
matched.remove(&bond.id());
d_bnd_cands[aid1.index()] = true;
d_bnd_cands[aid2.index()] = true;
}
}
}
*done = tdone;
}
fn worker_sorted_atoms(
topology: &TopologyBlock,
all_atoms: &[AtomId],
all_atom_set: &BTreeSet<AtomId>,
atom_ranks: &[usize],
) -> Vec<AtomId> {
let mut wedge_end_atoms = BTreeSet::new();
for bond in &topology.bonds {
if matches!(
bond.direction(),
BondDirection::BeginWedge | BondDirection::BeginDash
) && all_atom_set.contains(&bond.end())
{
wedge_end_atoms.insert(bond.end());
}
}
let mut sorted_atoms = all_atoms.to_vec();
sorted_atoms.sort_by(|left, right| {
match (
wedge_end_atoms.contains(left),
wedge_end_atoms.contains(right),
) {
(true, false) => std::cmp::Ordering::Less,
(false, true) => std::cmp::Ordering::Greater,
_ => (atom_ranks[left.index()], left.index())
.cmp(&(atom_ranks[right.index()], right.index())),
}
});
sorted_atoms
}
fn question_switch_masks(questions: &[AtomId]) -> Vec<BTreeSet<AtomId>> {
let question_count = questions.len();
(1usize..(1usize << question_count))
.map(|mask| {
questions
.iter()
.enumerate()
.filter_map(|(idx, atom)| ((mask & (1usize << idx)) != 0).then_some(*atom))
.collect()
})
.collect()
}
fn rdkit_kekulize_default_valence(atomic_number: u8) -> Result<i32, ValenceError> {
crate::valence::rdkit_default_valence(atomic_number)
}
#[allow(dead_code)]
fn kekulize_if_possible_assignment(
molecule: &Molecule,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<Option<KekulizeAssignment>, KekulizeError> {
match kekulize_assignment(
molecule,
None,
clear_aromatic_flags,
canonical,
max_backtracks,
) {
Ok(assignment) => Ok(Some(assignment)),
Err(error) if kekulize_if_possible_is_recoverable_failure(&error) => Ok(None),
Err(error) => Err(error),
}
}
#[allow(dead_code)]
fn kekulize_if_possible_value(
molecule: &Molecule,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<(Molecule, bool), KekulizeError> {
match kekulize_if_possible_assignment(
molecule,
clear_aromatic_flags,
canonical,
max_backtracks,
)? {
Some(assignment) => {
let mut kekulized = molecule.clone();
apply_kekulize_assignment(kekulized.topology_block_mut(), &assignment);
Ok((kekulized, true))
}
None => Ok((molecule.clone(), false)),
}
}
fn kekulize_if_possible_is_recoverable_failure(error: &KekulizeError) -> bool {
matches!(
error,
KekulizeError::UnkekulizedAtoms { .. }
| KekulizeError::NonRingAromaticAtom { .. }
| KekulizeError::ValenceChanged { .. }
| KekulizeError::Valence(_)
| KekulizeError::RingFinding(_)
)
}
fn rdkit_is_early_atom(atomic_number: u8) -> bool {
matches!(
atomic_number,
3 | 4
| 5
| 11
| 12
| 13
| 19
| 20
| 21
| 22
| 30
| 31
| 32
| 37
| 38
| 39
| 40
| 41
| 48
| 49
| 50
| 51
| 55
| 56
| 57
| 58
| 59
| 60
| 61
| 72
| 73
| 80
| 81
| 82
| 83
| 87
| 88
| 89
| 90
| 91
| 92
| 93
| 104
| 105
| 106
| 107
| 108
| 109
| 110
| 111
| 112
| 113
| 114
| 115
| 116
| 117
| 118
)
}
fn ordered_kekulize_rings(
topology: &TopologyBlock,
rings: &RingInfo,
atoms_to_use: &[bool],
bonds_to_use: &[bool],
wedged_atoms: &BTreeSet<AtomId>,
) -> Vec<KekulizeRing> {
let mut front = Vec::new();
let mut back = Vec::new();
for ring_idx in 0..rings.num_rings() {
let atoms = &rings.atom_rings()[ring_idx];
let bonds = &rings.bond_rings()[ring_idx];
let contains_non_dummy = atoms
.iter()
.all(|atom| atom.index() < topology.atoms.len() && atoms_to_use[atom.index()])
&& atoms
.iter()
.any(|atom| topology.atoms[atom.index()].atomic_number() != 0);
if !contains_non_dummy {
continue;
}
if !bonds.iter().all(|bond| bonds_to_use[bond.index()]) {
continue;
}
let start_pos = atoms
.iter()
.position(|atom| wedged_atoms.contains(atom))
.unwrap_or(0);
let rotated_atoms = (0..atoms.len())
.map(|offset| atoms[(offset + start_pos) % atoms.len()])
.collect::<Vec<_>>();
let entry = KekulizeRing {
atoms: rotated_atoms,
bonds: bonds.clone(),
source_ring: ring_idx,
};
if start_pos == 0 && !atoms.iter().any(|atom| wedged_atoms.contains(atom)) {
back.push(entry);
} else {
front.push(entry);
}
}
front.extend(back);
front
}
fn fused_ring_systems(rings: &[KekulizeRing]) -> Vec<Vec<usize>> {
let neighbor_map = crate::source_port_helpers::rdkit_make_ring_neighbor_map(
rings.len(),
|index| rings[index].bonds.as_slice(),
0,
0,
);
let mut systems = Vec::<Vec<usize>>::new();
let mut done = vec![false; rings.len()];
let mut current = 0usize;
while current < rings.len() {
let mut fused = Vec::new();
crate::source_port_helpers::rdkit_pick_fused_rings(
current,
&neighbor_map,
&mut fused,
&mut done,
0,
)
.expect("RingUtils map contains every source ring index");
systems.push(fused);
let Some(next) = done.iter().position(|is_done| !*is_done) else {
break;
};
current = next;
}
systems
}
fn fused_ring_systems_stack_order_for_tests(rings: &[KekulizeRing]) -> Vec<Vec<usize>> {
let mut systems = Vec::<Vec<usize>>::new();
let mut seen = vec![false; rings.len()];
for start in 0..rings.len() {
if seen[start] {
continue;
}
let mut stack = vec![start];
let mut system = Vec::new();
seen[start] = true;
while let Some(ring) = stack.pop() {
system.push(ring);
for next in 0..rings.len() {
if !seen[next] && rings_share_bond(rings, ring, next) {
seen[next] = true;
stack.push(next);
}
}
}
systems.push(system);
}
systems
}
fn rings_share_bond(rings: &[KekulizeRing], left: usize, right: usize) -> bool {
rings[left]
.bonds
.iter()
.any(|bond| rings[right].bonds.contains(bond))
}
#[cfg(test)]
mod tests {
use super::*;
use crate::{
AtomSpec, BondOrder, BondQueryPredicate, BondSpec, Element, MoleculeBuilder, QueryNode,
symmetrize_sssr,
};
fn kekulize_fragment_assignment(
molecule: &Molecule,
rings: &RingInfo,
atoms_to_use: &[bool],
bonds_to_use: Vec<bool>,
clear_aromatic_flags: bool,
canonical: bool,
max_backtracks: usize,
) -> Result<KekulizeAssignment, KekulizeError> {
super::kekulize_fragment_assignment(
molecule.topology_block(),
molecule.stereo_groups(),
Some(rings),
None,
atoms_to_use,
bonds_to_use,
clear_aromatic_flags,
canonical,
max_backtracks,
)
}
fn ordered_kekulize_rings(
molecule: &Molecule,
rings: &RingInfo,
atoms_to_use: &[bool],
bonds_to_use: &[bool],
wedged_atoms: &BTreeSet<AtomId>,
) -> Vec<KekulizeRing> {
super::ordered_kekulize_rings(
molecule.topology_block(),
rings,
atoms_to_use,
bonds_to_use,
wedged_atoms,
)
}
#[test]
fn nonaromatic_fused_product_defers_ring_discovery_to_aromaticity_like_rdkit() {
let mut molecule =
Molecule::from_smiles_with_sanitize("Cc1nc2c(nc1C)C(=O)C1=C(C2=O)C2C=CC1CC2", false)
.unwrap();
let explicit_hydrogens = [3, 0, 1, 1, 0, 0, 0, 3, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 1, 1];
let no_implicit = [
true, false, true, true, false, true, false, true, false, true, true, false, false,
true, true, true, false, true, true, true,
];
let bond_orders = [
BondOrder::Single,
BondOrder::Single,
BondOrder::Single,
BondOrder::Single,
BondOrder::Double,
BondOrder::Single,
BondOrder::Single,
BondOrder::Single,
BondOrder::Single,
BondOrder::Double,
BondOrder::Single,
BondOrder::Single,
BondOrder::Double,
BondOrder::Double,
BondOrder::Single,
BondOrder::Single,
BondOrder::Double,
BondOrder::Single,
BondOrder::Double,
BondOrder::Double,
BondOrder::Single,
BondOrder::Single,
BondOrder::Single,
];
{
let topology = molecule.topology_block_mut();
for (atom, (&hydrogens, &no_implicit)) in topology
.atoms
.iter_mut()
.zip(explicit_hydrogens.iter().zip(&no_implicit))
{
atom.set_aromatic(false);
atom.set_explicit_hydrogens(hydrogens);
atom.set_no_implicit(no_implicit);
}
for (bond, &order) in topology.bonds.iter_mut().zip(&bond_orders) {
bond.set_order(order);
bond.set_aromatic(false);
}
}
let valence =
crate::assign_valence_with_options(&molecule, crate::ValenceModel::RdkitLike, false)
.unwrap();
molecule.derived_cache_mut().rings = None;
molecule.derived_cache_mut().valence = Some(valence);
let result = molecule
.sanitize_with_ops(crate::SanitizeOps::KEKULIZE | crate::SanitizeOps::SET_AROMATICITY)
.unwrap();
let rings = result
.derived_cache()
.rings
.as_ref()
.expect("SetAromaticity should initialize symmetric ring information");
assert_eq!(rings.find_type(), crate::RingFindType::SymmSssr);
assert_eq!(rings.num_rings(), 5);
assert_eq!(
result
.atoms()
.iter()
.filter(|atom| atom.is_aromatic())
.map(|atom| atom.id().index())
.collect::<Vec<_>>(),
vec![10, 11, 14, 17, 18, 19]
);
assert_eq!(
result
.bonds()
.iter()
.filter(|bond| bond.is_aromatic())
.map(|bond| bond.id().index())
.collect::<Vec<_>>(),
vec![10, 13, 17, 18, 20, 22]
);
}
fn mark_double_bond_candidates(
molecule: &Molecule,
rings: &RingInfo,
kekulize_rings: &[KekulizeRing],
fused_ring_indices: &[usize],
assignment: &mut KekulizeAssignment,
) -> Result<KekulizeCandidateState, KekulizeError> {
let valence = crate::valence::assign_valence(&molecule, crate::ValenceModel::RdkitLike)?;
super::mark_double_bond_candidates(
molecule.topology_block(),
rings,
kekulize_rings,
fused_ring_indices,
&valence,
assignment,
)
}
fn kekulize_fused_assignment(
molecule: &Molecule,
rings: &RingInfo,
kekulize_rings: &[KekulizeRing],
fused_ring_indices: &[usize],
atom_ranks: &[usize],
max_backtracks: usize,
assignment: &mut KekulizeAssignment,
) -> Result<(), KekulizeError> {
let valence = crate::valence::assign_valence(&molecule, crate::ValenceModel::RdkitLike)?;
super::kekulize_fused_assignment(
molecule.topology_block(),
rings,
kekulize_rings,
fused_ring_indices,
&valence,
atom_ranks,
max_backtracks,
assignment,
)
}
fn kekulize_worker_matching(
molecule: &Molecule,
state: &KekulizeCandidateState,
initial_done: &[AtomId],
atom_ranks: &[usize],
max_backtracks: usize,
switched_off: &BTreeSet<AtomId>,
) -> Option<BTreeSet<BondId>> {
super::kekulize_worker_matching(
molecule.topology_block(),
state,
initial_done,
atom_ranks,
max_backtracks,
switched_off,
)
}
fn permute_dummies_and_match(
molecule: &Molecule,
state: &KekulizeCandidateState,
questions: &[AtomId],
atom_ranks: &[usize],
max_backtracks: usize,
) -> Option<BTreeSet<BondId>> {
super::permute_dummies_and_match(
molecule.topology_block(),
state,
questions,
atom_ranks,
max_backtracks,
)
}
#[allow(clippy::too_many_arguments)]
fn back_track(
molecule: &Molecule,
options: &mut BTreeMap<AtomId, VecDeque<(AtomId, BondId)>>,
last_opt: AtomId,
done: &mut Vec<AtomId>,
done_flags: &mut [bool],
aqueue: &mut VecDeque<AtomId>,
astack_flags: &mut [bool],
d_bnd_cands: &mut [bool],
d_bnd_adds: &mut [bool],
matched: &mut BTreeSet<BondId>,
) {
super::back_track(
molecule.topology_block(),
options,
last_opt,
done,
done_flags,
aqueue,
astack_flags,
d_bnd_cands,
d_bnd_adds,
matched,
);
}
fn worker_sorted_atoms(
molecule: &Molecule,
all_atoms: &[AtomId],
all_atom_set: &BTreeSet<AtomId>,
atom_ranks: &[usize],
) -> Vec<AtomId> {
super::worker_sorted_atoms(
molecule.topology_block(),
all_atoms,
all_atom_set,
atom_ranks,
)
}
fn build_fused_aromatic_naphthalene(atom_is_aromatic: bool) -> Molecule {
let mut builder = MoleculeBuilder::new();
let atoms = (0..10)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(atom_is_aromatic)))
.collect::<Vec<_>>();
for &(begin, end) in &[
(0, 1),
(1, 2),
(2, 3),
(3, 4),
(4, 5),
(5, 0),
(4, 6),
(6, 7),
(7, 8),
(8, 9),
(9, 5),
] {
builder
.add_bond(
BondSpec::new(atoms[begin], atoms[end], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
builder.build().unwrap()
}
fn build_fused_dummy_aromatic_naphthalene() -> Molecule {
let mut builder = MoleculeBuilder::new();
let mut atoms = (0..10)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
atoms[0] = builder.add_atom(AtomSpec::new(Element::DUMMY));
for &(begin, end) in &[
(0, 1),
(1, 2),
(2, 3),
(3, 4),
(4, 5),
(5, 0),
(4, 6),
(6, 7),
(7, 8),
(8, 9),
(9, 5),
] {
builder
.add_bond(
BondSpec::new(atoms[begin], atoms[end], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
builder.build().unwrap()
}
fn fused_ring_inputs(molecule: &Molecule) -> (RingInfo, Vec<KekulizeRing>, Vec<Vec<usize>>) {
let rings = symmetrize_sssr(molecule).unwrap();
let kekulize_rings = ordered_kekulize_rings(
molecule,
&rings,
&vec![true; molecule.num_atoms()],
&vec![true; molecule.num_bonds()],
&std::collections::BTreeSet::new(),
);
let ring_systems = fused_ring_systems(&kekulize_rings);
(rings, kekulize_rings, ring_systems)
}
#[test]
fn shared_ring_utils_kekulize_preserves_recursive_neighbor_order() {
let rings = vec![
KekulizeRing {
atoms: Vec::new(),
bonds: vec![BondId::new(0), BondId::new(1)],
source_ring: 0,
},
KekulizeRing {
atoms: Vec::new(),
bonds: vec![BondId::new(1), BondId::new(2)],
source_ring: 1,
},
KekulizeRing {
atoms: Vec::new(),
bonds: vec![BondId::new(0), BondId::new(3)],
source_ring: 2,
},
KekulizeRing {
atoms: Vec::new(),
bonds: vec![BondId::new(2), BondId::new(4)],
source_ring: 3,
},
];
assert_eq!(fused_ring_systems(&rings), vec![vec![0, 1, 3, 2]]);
}
#[test]
fn kekulize_assignment_returns_empty_for_empty_molecule() {
let molecule = Molecule::new();
let assignment = kekulize_assignment(&molecule, None, true, false, 100).unwrap();
assert!(assignment.is_empty());
}
#[test]
fn kekulize_assignment_kekulizes_benzene_like_aromatic_cycle() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let assignment = kekulize_assignment(&molecule, None, true, false, 100).unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
let single_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Single))
.count();
assert_eq!(double_bonds, 3);
assert_eq!(single_bonds, 3);
assert!(
molecule
.bonds()
.iter()
.all(|bond| assignment.bond_aromatic_flag(bond.id()) == Some(false))
);
assert!(
molecule
.atoms()
.iter()
.all(|atom| assignment.atom_aromatic_flag(atom.id()) == Some(false))
);
}
#[test]
fn kekulize_assignment_preserves_aromatic_flags_when_clear_is_disabled_like_rdkit() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let assignment = kekulize_assignment(&molecule, None, false, false, 100).unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
let single_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Single))
.count();
assert_eq!(double_bonds, 3);
assert_eq!(single_bonds, 3);
assert!(
molecule
.bonds()
.iter()
.all(|bond| assignment.bond_aromatic_flag(bond.id()).is_none())
);
assert!(
molecule
.atoms()
.iter()
.all(|atom| assignment.atom_aromatic_flag(atom.id()).is_none())
);
}
#[test]
fn kekulize_assignment_noncanonical_keeps_atom_index_seed_for_benzene_like_rdkit_core_api() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let assignment = kekulize_assignment(&molecule, None, false, false, 100).unwrap();
let bond_orders = molecule
.bonds()
.iter()
.map(|bond| assignment.bond_order(bond.id()).unwrap())
.collect::<Vec<_>>();
assert_eq!(
bond_orders,
vec![
BondOrder::Double,
BondOrder::Single,
BondOrder::Double,
BondOrder::Single,
BondOrder::Double,
BondOrder::Single
]
);
}
#[test]
fn kekulize_assignment_fragment_returns_empty_when_atoms_to_use_is_empty_like_rdkit() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let assignment = kekulize_fragment_assignment(
&molecule,
&rings,
&[false; 6],
vec![true; molecule.num_bonds()],
true,
false,
100,
)
.unwrap();
assert!(assignment.is_empty());
}
#[test]
fn kekulize_assignment_fragment_allows_partial_scope_non_ring_aromatic_atom_like_rdkit() {
let mut builder = MoleculeBuilder::new();
let a0 = builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true));
let _a1 = builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true));
let molecule = builder.build().unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let assignment = kekulize_fragment_assignment(
&molecule,
&rings,
&[true, false],
vec![],
true,
false,
100,
)
.unwrap();
assert_eq!(assignment.atom_aromatic_flag(a0), Some(false));
}
#[test]
fn kekulize_fragment_rejects_bitset_size_mismatch() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let error = kekulize_fragment_assignment(
&molecule,
&rings,
&[true; 5],
vec![true; molecule.num_bonds()],
true,
false,
100,
)
.unwrap_err();
assert!(matches!(
error,
KekulizeError::FragmentBitsetSizeMismatch { atoms: 5, bonds: 6 }
));
}
#[test]
fn kekulize_fragment_rejects_bond_mask_size_mismatch_like_rdkit() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let error = kekulize_fragment_assignment(
&molecule,
&rings,
&[true; 6],
vec![true; 5],
true,
false,
100,
)
.unwrap_err();
assert!(matches!(
error,
KekulizeError::FragmentBitsetSizeMismatch { atoms: 6, bonds: 5 }
));
}
#[test]
fn kekulize_fragment_filters_out_rings_outside_selected_atom_mask_like_rdkit() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let mut atoms_to_use = vec![true; molecule.num_atoms()];
atoms_to_use[0] = false;
let assignment = kekulize_fragment_assignment(
&molecule,
&rings,
&atoms_to_use,
vec![true; molecule.num_bonds()],
false,
false,
100,
)
.unwrap();
assert!(assignment.bond_orders.iter().all(Option::is_none));
}
#[test]
fn kekulize_fragment_clears_pyrrolic_explicit_hydrogen_state_like_rdkit() {
let molecule = Molecule::from_smiles_with_sanitize("c1cc[nH]c1", false).unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let pyrrolic = molecule
.atoms()
.iter()
.find(|atom| atom.atomic_number() == 7)
.map(|atom| atom.id())
.unwrap();
let assignment = kekulize_fragment_assignment(
&molecule,
&rings,
&vec![true; molecule.num_atoms()],
vec![true; molecule.num_bonds()],
true,
false,
100,
)
.unwrap();
assert_eq!(assignment.atom_aromatic_flag(pyrrolic), Some(false));
assert_eq!(assignment.atom_no_implicit(pyrrolic), Some(false));
assert_eq!(assignment.atom_explicit_hydrogens(pyrrolic), Some(0));
}
#[test]
fn kekulize_assignment_runs_canonical_fragment_ranking_for_plain_aromatic_cycle() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let assignment = kekulize_assignment(&molecule, None, true, true, 100).unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
let single_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Single))
.count();
assert_eq!(double_bonds, 3);
assert_eq!(single_bonds, 3);
}
#[test]
fn kekulize_if_possible_returns_assignment_for_kekulizable_molecule() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let assignment = kekulize_if_possible_assignment(&molecule, true, false, 100)
.unwrap()
.unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
assert_eq!(double_bonds, 3);
}
#[test]
fn kekulize_if_possible_preserves_aromatic_flags_when_clear_is_disabled() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let assignment = kekulize_if_possible_assignment(&molecule, false, false, 100)
.unwrap()
.unwrap();
assert!(
molecule
.bonds()
.iter()
.all(|bond| assignment.bond_aromatic_flag(bond.id()).is_none())
);
assert!(
molecule
.atoms()
.iter()
.all(|atom| assignment.atom_aromatic_flag(atom.id()).is_none())
);
}
#[test]
fn kekulize_if_possible_runs_canonical_fragment_ranking_for_plain_aromatic_cycle() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let assignment = kekulize_if_possible_assignment(&molecule, true, true, 100)
.unwrap()
.unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
assert_eq!(double_bonds, 3);
}
#[test]
fn kekulize_if_possible_assignment_returns_none_for_recoverable_failure() {
let mut builder = MoleculeBuilder::new();
let atoms: Vec<_> = (0..5)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect();
for i in 0..5 {
builder
.add_bond(
BondSpec::new(atoms[i], atoms[(i + 1) % 5], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
let molecule = builder.build().unwrap();
let assignment = kekulize_if_possible_assignment(&molecule, true, false, 100).unwrap();
assert!(assignment.is_none());
}
#[test]
fn kekulize_if_possible_value_returns_kekulized_clone_on_success() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let (kekulized, changed) = kekulize_if_possible_value(&molecule, true, false, 100).unwrap();
assert!(changed);
assert!(
kekulized
.bonds()
.iter()
.any(|bond| bond.order() == BondOrder::Double)
);
assert!(kekulized.bonds().iter().all(|bond| !bond.is_aromatic()));
assert!(
molecule
.bonds()
.iter()
.all(|bond| bond.order() == BondOrder::Aromatic && bond.is_aromatic())
);
}
#[test]
fn kekulize_if_possible_value_returns_original_clone_when_kekulization_fails() {
let mut builder = MoleculeBuilder::new();
let atoms = (0..5)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
for idx in 0..5 {
builder
.add_bond(
BondSpec::new(atoms[idx], atoms[(idx + 1) % 5], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
let molecule = builder.build().unwrap();
let (restored, changed) = kekulize_if_possible_value(&molecule, true, false, 100).unwrap();
assert!(!changed);
assert_eq!(restored.bonds().len(), molecule.bonds().len());
assert!(restored.bonds().iter().zip(molecule.bonds().iter()).all(
|(restored_bond, original_bond)| restored_bond.order() == original_bond.order()
&& restored_bond.is_aromatic() == original_bond.is_aromatic()
));
assert!(
restored
.atoms()
.iter()
.zip(molecule.atoms().iter())
.all(|(restored_atom, original_atom)| restored_atom.is_aromatic()
== original_atom.is_aromatic())
);
}
#[test]
fn ordered_kekulize_rings_prioritizes_wedged_rings_like_rdkit() {
let mut builder = MoleculeBuilder::new();
let ring_one = (0..3)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
let ring_two = (0..3)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
for idx in 0..3 {
builder
.add_bond(
BondSpec::new(ring_one[idx], ring_one[(idx + 1) % 3], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
for idx in 0..3 {
let mut bond =
BondSpec::new(ring_two[idx], ring_two[(idx + 1) % 3], BondOrder::Aromatic)
.with_aromatic(true);
if idx == 0 {
bond = bond.with_direction(crate::BondDirection::BeginWedge);
}
builder.add_bond(bond).unwrap();
}
let molecule = builder.build().unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let wedged_atoms = std::iter::once(ring_two[0]).collect::<std::collections::BTreeSet<_>>();
let ordered = ordered_kekulize_rings(
&molecule,
&rings,
&vec![true; molecule.num_atoms()],
&vec![true; molecule.num_bonds()],
&wedged_atoms,
);
assert_eq!(ordered.len(), 2);
assert_eq!(ordered[0].source_ring, 1);
}
#[test]
fn ordered_kekulize_rings_rotates_wedged_ring_start_atom_like_rdkit() {
let mut builder = MoleculeBuilder::new();
let ring = (0..4)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
for idx in 0..4 {
let mut bond = BondSpec::new(ring[idx], ring[(idx + 1) % 4], BondOrder::Aromatic)
.with_aromatic(true);
if idx == 2 {
bond = bond.with_direction(crate::BondDirection::BeginDash);
}
builder.add_bond(bond).unwrap();
}
let molecule = builder.build().unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let wedged_atoms = std::iter::once(ring[2]).collect::<std::collections::BTreeSet<_>>();
let ordered = ordered_kekulize_rings(
&molecule,
&rings,
&vec![true; molecule.num_atoms()],
&vec![true; molecule.num_bonds()],
&wedged_atoms,
);
assert_eq!(ordered.len(), 1);
assert_eq!(ordered[0].atoms[0], ring[2]);
}
#[test]
fn canonical_kekulize_handles_stereo_ranking_branch() {
let mut builder = MoleculeBuilder::new();
let atoms = (0..6)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
for idx in 0..6 {
builder
.add_bond(
BondSpec::new(atoms[idx], atoms[(idx + 1) % 6], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
let mut molecule = builder.build().unwrap();
molecule.topology_block_mut().atoms[0].set_chiral_tag(crate::ChiralTag::TetrahedralCw);
let assignment = kekulize_assignment(&molecule, None, true, true, 100).unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
assert_eq!(double_bonds, 3);
}
#[test]
fn kekulize_assignment_handles_dummy_atom_permutation_branch() {
let molecule = Molecule::from_smiles_with_sanitize("*1cccc1", false).unwrap();
let assignment = kekulize_assignment(&molecule, None, true, false, 100).unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
let single_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Single))
.count();
assert_eq!(double_bonds, 2, "{:?}", assignment.bond_orders);
assert_eq!(single_bonds, 1, "{:?}", assignment.bond_orders);
}
#[test]
fn kekulize_assignment_detects_aromaticity_from_bond_state_without_atom_flags() {
let mut builder = MoleculeBuilder::new();
let atoms = (0..6)
.map(|_| builder.add_atom(AtomSpec::new(Element::C)))
.collect::<Vec<_>>();
for idx in 0..6 {
builder
.add_bond(
BondSpec::new(atoms[idx], atoms[(idx + 1) % 6], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
let molecule = builder.build().unwrap();
let assignment = kekulize_assignment(&molecule, None, true, false, 100).unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
assert_eq!(double_bonds, 3);
}
#[test]
fn kekulize_assignment_detects_fused_aromaticity_from_bond_state_without_atom_flags() {
let molecule = build_fused_aromatic_naphthalene(false);
let assignment = kekulize_assignment(&molecule, None, true, false, 100).unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
assert_eq!(double_bonds, 5);
}
#[test]
fn kekulize_if_possible_mark_double_bond_candidates_tracks_dummy_atoms_in_fused_systems_like_rdkit()
{
let molecule = build_fused_dummy_aromatic_naphthalene();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
assert!(state.questions.contains(&AtomId::new(10)));
}
#[test]
fn kekulize_if_possible_mark_double_bond_candidates_marks_non_aromatic_atoms_done_inside_fused_systems()
{
let mut molecule = build_fused_aromatic_naphthalene(true);
molecule.topology_block_mut().atoms[0].set_aromatic(false);
for bond_idx in [0usize, 5usize] {
molecule.topology_block_mut().bonds[bond_idx].set_aromatic(false);
molecule.topology_block_mut().bonds[bond_idx].set_order(BondOrder::Single);
}
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
assert!(state.done.contains(&AtomId::new(0)));
assert!(!state.candidate_atoms.contains(&AtomId::new(0)));
}
#[test]
fn kekulize_fused_returns_empty_candidate_state_without_aromatic_or_dummy_atoms() {
let mut molecule = build_fused_aromatic_naphthalene(true);
for atom in molecule.topology_block_mut().atoms.iter_mut() {
atom.set_aromatic(false);
}
for bond in molecule.topology_block_mut().bonds.iter_mut() {
bond.set_aromatic(false);
bond.set_order(BondOrder::Single);
}
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
assert!(state.candidate_atoms.is_empty());
assert!(state.aromatic_edges.is_empty());
assert!(state.questions.is_empty());
assert!(state.done.is_empty());
}
#[test]
fn kekulize_assignment_reports_problem_atoms_for_unmatchable_fused_system() {
let mut builder = MoleculeBuilder::new();
let atoms = (0..5)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
for idx in 0..5 {
builder
.add_bond(
BondSpec::new(atoms[idx], atoms[(idx + 1) % 5], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
let molecule = builder.build().unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let kekulize_rings = ordered_kekulize_rings(
&molecule,
&rings,
&vec![true; molecule.num_atoms()],
&vec![true; molecule.num_bonds()],
&std::collections::BTreeSet::new(),
);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let error = kekulize_fused_assignment(
&molecule,
&rings,
&kekulize_rings,
&[0],
&(0..molecule.num_atoms()).collect::<Vec<_>>(),
100,
&mut assignment,
)
.unwrap_err();
assert!(matches!(
error,
KekulizeError::UnkekulizedAtoms { atoms }
if atoms == vec![
AtomId::new(0),
AtomId::new(1),
AtomId::new(2),
AtomId::new(3),
AtomId::new(4),
]
));
}
#[test]
fn kekulize_fused_assignment_sets_bond_orders_and_clears_directions_like_rdkit() {
let mut molecule = build_fused_aromatic_naphthalene(true);
for bond in molecule.topology_block_mut().bonds.iter_mut() {
bond.set_direction(BondDirection::BeginWedge);
}
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
kekulize_fused_assignment(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&(0..molecule.num_atoms()).collect::<Vec<_>>(),
100,
&mut assignment,
)
.unwrap();
let assigned_orders = molecule
.bonds()
.iter()
.filter_map(|bond| assignment.bond_order(bond.id()))
.collect::<Vec<_>>();
let num_doubles = assigned_orders
.iter()
.filter(|order| **order == BondOrder::Double)
.count();
let num_singles = assigned_orders
.iter()
.filter(|order| **order == BondOrder::Single)
.count();
assert_eq!(num_doubles, 5);
assert_eq!(num_singles, 6);
assert!(molecule.bonds().iter().any(|bond| {
assignment.bond_order(bond.id()) == Some(BondOrder::Double)
&& assignment.bond_direction(bond.id()) == Some(BondDirection::None)
}));
}
#[test]
fn back_track_rolls_back_local_bonds_and_keeps_saved_options_like_rdkit() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let mut options = std::collections::BTreeMap::from([
(
AtomId::new(1),
std::collections::VecDeque::from([(AtomId::new(3), BondId::new(2))]),
),
(
AtomId::new(2),
std::collections::VecDeque::from([(AtomId::new(4), BondId::new(3))]),
),
]);
let mut done = vec![AtomId::new(0), AtomId::new(1), AtomId::new(2)];
let mut done_flags = vec![false; molecule.num_atoms()];
for atom in &done {
done_flags[atom.index()] = true;
}
let mut aqueue = std::collections::VecDeque::new();
let mut astack_flags = vec![false; molecule.num_atoms()];
let mut d_bnd_cands = vec![false; molecule.num_atoms()];
let mut d_bnd_adds = vec![false; molecule.num_bonds()];
d_bnd_adds[0] = true;
d_bnd_adds[1] = true;
let mut matched = [BondId::new(0), BondId::new(1)]
.into_iter()
.collect::<std::collections::BTreeSet<_>>();
back_track(
&molecule,
&mut options,
AtomId::new(1),
&mut done,
&mut done_flags,
&mut aqueue,
&mut astack_flags,
&mut d_bnd_cands,
&mut d_bnd_adds,
&mut matched,
);
assert_eq!(done, vec![AtomId::new(0)]);
assert!(done_flags[0]);
assert!(!done_flags[1]);
assert!(!done_flags[2]);
assert_eq!(
aqueue.into_iter().collect::<Vec<_>>(),
vec![AtomId::new(1), AtomId::new(2)]
);
assert!(astack_flags[1]);
assert!(astack_flags[2]);
assert!(d_bnd_adds[0]);
assert!(!d_bnd_adds[1]);
assert!(matched.contains(&BondId::new(0)));
assert!(!matched.contains(&BondId::new(1)));
assert!(d_bnd_cands[1]);
assert!(d_bnd_cands[2]);
assert_eq!(
options.get(&AtomId::new(1)).cloned(),
Some(std::collections::VecDeque::from([(
AtomId::new(3),
BondId::new(2),
)]))
);
assert_eq!(
options.get(&AtomId::new(2)).cloned(),
Some(std::collections::VecDeque::from([(
AtomId::new(4),
BondId::new(3),
)]))
);
}
#[test]
fn kekulize_if_possible_worker_back_track_keeps_saved_options_for_replay() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let expected = std::collections::VecDeque::from([(AtomId::new(3), BondId::new(2))]);
let mut options = std::collections::BTreeMap::from([
(AtomId::new(1), expected.clone()),
(
AtomId::new(2),
std::collections::VecDeque::from([(AtomId::new(4), BondId::new(3))]),
),
]);
let mut done = vec![AtomId::new(0), AtomId::new(1), AtomId::new(2)];
let mut done_flags = vec![false; molecule.num_atoms()];
for atom in &done {
done_flags[atom.index()] = true;
}
let mut aqueue = std::collections::VecDeque::new();
let mut astack_flags = vec![false; molecule.num_atoms()];
let mut d_bnd_cands = vec![false; molecule.num_atoms()];
let mut d_bnd_adds = vec![false; molecule.num_bonds()];
d_bnd_adds[0] = true;
let mut matched = std::collections::BTreeSet::from([BondId::new(0)]);
back_track(
&molecule,
&mut options,
AtomId::new(1),
&mut done,
&mut done_flags,
&mut aqueue,
&mut astack_flags,
&mut d_bnd_cands,
&mut d_bnd_adds,
&mut matched,
);
assert_eq!(options.get(&AtomId::new(1)), Some(&expected));
assert_eq!(matched, std::collections::BTreeSet::from([BondId::new(0)]));
}
#[test]
fn kekulize_if_possible_permutes_dummy_candidates_after_initial_worker_failure() {
let molecule = Molecule::from_smiles_with_sanitize("*1cccc1", false).unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
let atom_ranks = (0..molecule.num_atoms()).collect::<Vec<_>>();
let initial = kekulize_worker_matching(
&molecule,
&state,
&state.done,
&atom_ranks,
100,
&std::collections::BTreeSet::new(),
);
let permuted =
permute_dummies_and_match(&molecule, &state, &state.questions, &atom_ranks, 100);
assert!(initial.is_none());
assert_eq!(
permuted.as_ref().map(|matched| matched.len()),
Some(2),
"{permuted:?}"
);
}
#[test]
fn permute_dummies_and_match_can_recover_without_backtracking_via_mask_switch() {
let molecule = Molecule::from_smiles_with_sanitize("*1cccc1", false).unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
let atom_ranks = (0..molecule.num_atoms()).collect::<Vec<_>>();
let permuted =
permute_dummies_and_match(&molecule, &state, &state.questions, &atom_ranks, 0);
assert_eq!(permuted.as_ref().map(|matched| matched.len()), Some(2));
}
#[test]
fn permute_dummies_and_match_returns_none_after_exhausting_all_switch_masks() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
let atom_ranks = (0..molecule.num_atoms()).collect::<Vec<_>>();
let permuted =
permute_dummies_and_match(&molecule, &state, &[AtomId::new(0)], &atom_ranks, 100);
assert!(permuted.is_none());
}
#[test]
fn kekulize_worker_matching_succeeds_for_plain_aromatic_cycle() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
let atom_ranks = (0..molecule.num_atoms()).collect::<Vec<_>>();
let matched = kekulize_worker_matching(
&molecule,
&state,
&state.done,
&atom_ranks,
100,
&std::collections::BTreeSet::new(),
);
assert_eq!(matched.as_ref().map(|set| set.len()), Some(3));
}
#[test]
fn kekulize_worker_matching_returns_none_when_backtracks_are_exhausted() {
let molecule = Molecule::from_smiles_with_sanitize("*1cccc1", false).unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
let atom_ranks = (0..molecule.num_atoms()).collect::<Vec<_>>();
let matched = kekulize_worker_matching(
&molecule,
&state,
&state.done,
&atom_ranks,
0,
&std::collections::BTreeSet::new(),
);
assert!(matched.is_none());
}
#[test]
fn kekulize_worker_matching_returns_empty_match_when_candidates_are_switched_off() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
let atom_ranks = (0..molecule.num_atoms()).collect::<Vec<_>>();
let switched_off = state.candidate_atoms.clone();
let matched = kekulize_worker_matching(
&molecule,
&state,
&state.done,
&atom_ranks,
100,
&switched_off,
);
assert_eq!(matched, Some(std::collections::BTreeSet::new()));
}
#[test]
fn permute_dummies_and_match_returns_none_without_questions() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
let atom_ranks = (0..molecule.num_atoms()).collect::<Vec<_>>();
let matched = permute_dummies_and_match(&molecule, &state, &[], &atom_ranks, 100);
assert!(matched.is_none());
}
#[test]
fn kekulize_if_possible_does_not_treat_port_shape_errors_as_recoverable() {
assert!(!kekulize_if_possible_is_recoverable_failure(
&KekulizeError::FragmentBitsetSizeMismatch { atoms: 1, bonds: 2 }
));
assert!(!kekulize_if_possible_is_recoverable_failure(
&KekulizeError::CanonicalRankSymbolSizeMismatch {
kind: "atom",
expected: 1,
actual: 0,
}
));
}
#[test]
fn question_switch_masks_follow_rdkit_bit_order() {
let masks = question_switch_masks(&[AtomId::new(1), AtomId::new(3), AtomId::new(5)]);
assert_eq!(
masks,
vec![
[AtomId::new(1)].into_iter().collect(),
[AtomId::new(3)].into_iter().collect(),
[AtomId::new(1), AtomId::new(3)].into_iter().collect(),
[AtomId::new(5)].into_iter().collect(),
[AtomId::new(1), AtomId::new(5)].into_iter().collect(),
[AtomId::new(3), AtomId::new(5)].into_iter().collect(),
[AtomId::new(1), AtomId::new(3), AtomId::new(5)]
.into_iter()
.collect(),
]
);
}
#[test]
fn question_switch_masks_returns_empty_for_no_questions() {
let masks = question_switch_masks(&[]);
assert!(masks.is_empty());
}
#[test]
fn worker_sorted_atoms_prioritizes_wedge_end_atoms_before_rank() {
let mut builder = MoleculeBuilder::new();
let atoms = (0..3)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
builder
.add_bond(
BondSpec::new(atoms[0], atoms[1], BondOrder::Aromatic)
.with_aromatic(true)
.with_direction(crate::BondDirection::BeginWedge),
)
.unwrap();
builder
.add_bond(BondSpec::new(atoms[1], atoms[2], BondOrder::Aromatic).with_aromatic(true))
.unwrap();
let molecule = builder.build().unwrap();
let all_atoms = vec![atoms[0], atoms[1], atoms[2]];
let all_atom_set = all_atoms
.iter()
.copied()
.collect::<std::collections::BTreeSet<_>>();
let sorted = worker_sorted_atoms(&molecule, &all_atoms, &all_atom_set, &[2, 10, 0]);
assert_eq!(sorted, vec![atoms[1], atoms[2], atoms[0]]);
}
#[test]
fn worker_sorted_atoms_uses_rank_order_when_no_wedge_bias_exists() {
let mut builder = MoleculeBuilder::new();
let atoms = (0..3)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
builder
.add_bond(BondSpec::new(atoms[0], atoms[1], BondOrder::Aromatic).with_aromatic(true))
.unwrap();
builder
.add_bond(BondSpec::new(atoms[1], atoms[2], BondOrder::Aromatic).with_aromatic(true))
.unwrap();
let molecule = builder.build().unwrap();
let all_atoms = vec![atoms[0], atoms[1], atoms[2]];
let all_atom_set = all_atoms
.iter()
.copied()
.collect::<std::collections::BTreeSet<_>>();
let sorted = worker_sorted_atoms(&molecule, &all_atoms, &all_atom_set, &[2, 0, 1]);
assert_eq!(sorted, vec![atoms[1], atoms[2], atoms[0]]);
}
#[test]
fn mark_double_bond_candidates_detects_dummy_and_aromatic_fused_state() {
let molecule = build_fused_dummy_aromatic_naphthalene();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
assert!(!state.aromatic_edges.is_empty());
assert!(!state.candidate_atoms.is_empty());
assert!(!state.questions.is_empty());
assert!(
state
.candidate_atoms
.contains(&AtomId::new(molecule.atoms().len() - 1))
);
}
#[test]
fn mark_double_bond_candidates_excludes_non_aromatic_ring_bond_from_aromatic_edges() {
let mut molecule = build_fused_dummy_aromatic_naphthalene();
molecule.topology_block_mut().bonds[0].set_aromatic(false);
molecule.topology_block_mut().bonds[0].set_order(BondOrder::Single);
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
&ring_systems[0],
&mut assignment,
)
.unwrap();
assert!(
!state
.aromatic_edges
.iter()
.any(|(bond, _, _)| *bond == BondId::new(0))
);
assert_eq!(assignment.bond_order(BondId::new(0)), None);
}
#[test]
fn kekulize_assignment_excludes_bond_type_query_from_bonds_to_use() {
let mut builder = MoleculeBuilder::new();
let a0 = builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true));
let a1 = builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true));
let a2 = builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true));
builder
.add_bond(
BondSpec::new(a0, a1, BondOrder::Aromatic)
.with_aromatic(true)
.with_query(QueryNode::predicate(BondQueryPredicate::OrderIn(vec![
BondOrder::Single,
BondOrder::Aromatic,
]))),
)
.unwrap();
builder
.add_bond(BondSpec::new(a1, a2, BondOrder::Aromatic).with_aromatic(true))
.unwrap();
builder
.add_bond(BondSpec::new(a2, a0, BondOrder::Aromatic).with_aromatic(true))
.unwrap();
let molecule = builder.build().unwrap();
let assignment = kekulize_assignment(&molecule, None, true, false, 100).unwrap();
assert_eq!(assignment.bond_order(BondId::new(0)), None);
assert_eq!(assignment.bond_aromatic_flag(BondId::new(0)), None);
assert_eq!(assignment.bond_aromatic_flag(BondId::new(1)), Some(false));
assert_eq!(assignment.bond_aromatic_flag(BondId::new(2)), Some(false));
}
#[test]
fn mark_double_bond_candidates_keeps_exocyclic_aryl_carbon_kekulizable_like_rdkit() {
let molecule =
Molecule::from_smiles("O=C1N(/N=C(/C)C1=NN/C2=C/C(OC)=CC=C2)C=3C=CC=CC=3").unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let phenyl_system = ring_systems
.iter()
.find(|system| {
system.iter().any(|&ring_idx| {
kekulize_rings[ring_idx]
.atoms
.iter()
.any(|atom| atom.index() == 17)
})
})
.unwrap();
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
phenyl_system,
&mut assignment,
)
.unwrap();
assert!(state.candidate_atoms.contains(&AtomId::new(17)));
assert!(kekulize_assignment(&molecule, None, true, true, 100).is_ok());
}
#[test]
fn mark_double_bond_candidates_excludes_aromatic_nh_in_row_134_like_rdkit() {
let molecule = Molecule::from_smiles(
"CC1=C(/C=C2\\C(C)=C3/C(C)=C(/C=C4\\C(C)=C5/C(CCC(=O)O)=C(C(C)=C5N4)C=C6\\C(CCC(=O)O)=C(C)C(=C6N2)C=C1)N3)C=C(\\C=C)C",
)
.unwrap();
let (rings, kekulize_rings, ring_systems) = fused_ring_inputs(&molecule);
let porphyrin_system = ring_systems
.iter()
.find(|system| {
system.iter().any(|&ring_idx| {
kekulize_rings[ring_idx]
.atoms
.iter()
.any(|atom| matches!(atom.index(), 26 | 39 | 42))
})
})
.unwrap();
let mut assignment = KekulizeAssignment::empty(molecule.num_atoms(), molecule.num_bonds());
let state = mark_double_bond_candidates(
&molecule,
&rings,
&kekulize_rings,
porphyrin_system,
&mut assignment,
)
.unwrap();
assert!(!state.candidate_atoms.contains(&AtomId::new(26)));
assert!(!state.candidate_atoms.contains(&AtomId::new(39)));
assert!(!state.candidate_atoms.contains(&AtomId::new(42)));
}
#[test]
fn kekulize_assignment_handles_row_134_porhyrin_like_rdkit() {
let molecule = Molecule::from_smiles(
"CC1=C(/C=C2\\C(C)=C3/C(C)=C(/C=C4\\C(C)=C5/C(CCC(=O)O)=C(C(C)=C5N4)C=C6\\C(CCC(=O)O)=C(C)C(=C6N2)C=C1)N3)C=C(\\C=C)C",
)
.unwrap();
assert!(kekulize_assignment(&molecule, None, true, true, 100).is_ok());
assert!(kekulize_assignment(&molecule, None, true, false, 100).is_ok());
}
fn kekulize_assignment_handles_pyrrolic_five_membered_ring_like_rdkit() {
let molecule = Molecule::from_smiles("[nH]1cccc1").unwrap();
assert_eq!(molecule.atoms()[0].explicit_hydrogens(), 1);
assert!(kekulize_assignment(&molecule, None, true, true, 100).is_ok());
assert!(kekulize_assignment(&molecule, None, true, false, 100).is_ok());
}
#[test]
fn kekulize_assignment_keeps_non_bond_type_query_in_bonds_to_use() {
let mut builder = MoleculeBuilder::new();
let atoms = (0..6)
.map(|_| builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true)))
.collect::<Vec<_>>();
builder
.add_bond(
BondSpec::new(atoms[0], atoms[1], BondOrder::Aromatic)
.with_aromatic(true)
.with_query(QueryNode::predicate(BondQueryPredicate::IsInRing(true))),
)
.unwrap();
for idx in 1..6 {
builder
.add_bond(
BondSpec::new(atoms[idx], atoms[(idx + 1) % 6], BondOrder::Aromatic)
.with_aromatic(true),
)
.unwrap();
}
let molecule = builder.build().unwrap();
let assignment = kekulize_assignment(&molecule, None, true, false, 100).unwrap();
assert_eq!(assignment.bond_aromatic_flag(BondId::new(0)), Some(false));
}
#[test]
fn kekulize_assignment_rejects_non_ring_aromatic_atom_like_rdkit() {
let mut builder = MoleculeBuilder::new();
builder.add_atom(AtomSpec::new(Element::C).with_aromatic(true));
let molecule = builder.build().unwrap();
let error = kekulize_assignment(&molecule, None, true, false, 100).unwrap_err();
assert!(
matches!(error, KekulizeError::NonRingAromaticAtom { atom } if atom == AtomId::new(0))
);
}
#[test]
fn kekulize_assignment_preserves_total_valence_like_rdkit() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let before =
crate::assign_valence_with_options(&molecule, crate::ValenceModel::RdkitLike, false)
.unwrap();
let assignment = kekulize_assignment(&molecule, None, true, false, 100).unwrap();
let mut kekulized = molecule.clone();
apply_kekulize_assignment(kekulized.topology_block_mut(), &assignment);
let after =
crate::assign_valence_with_options(&kekulized, crate::ValenceModel::RdkitLike, false)
.unwrap();
for atom in molecule.atoms() {
let idx = atom.id().index();
assert_eq!(
before.explicit_valence[idx] + before.implicit_hydrogens[idx],
after.explicit_valence[idx] + after.implicit_hydrogens[idx]
);
}
}
#[test]
fn kekulize_assignment_reuses_provided_ring_info_like_rdkit_entrypoint() {
let molecule = Molecule::from_smiles_with_sanitize("c1ccccc1", false).unwrap();
let rings = symmetrize_sssr(&molecule).unwrap();
let assignment = kekulize_assignment(&molecule, Some(&rings), true, false, 100).unwrap();
let double_bonds = molecule
.bonds()
.iter()
.filter(|bond| assignment.bond_order(bond.id()) == Some(BondOrder::Double))
.count();
assert_eq!(double_bonds, 3);
}
#[test]
fn kekulize_source_markers_do_not_hide_elided_rdkit_lines() {
let source = include_str!("kekulize.rs");
for line in source.lines().filter(|line| line.contains("RDKit")) {
assert!(
!line.contains("..."),
"RDKit source marker must not contain elided source: {line}"
);
}
}
}