diff --git a/arc/molecule/filtration.py b/arc/molecule/filtration.py index aa7eaf8d00..2a4e863196 100644 --- a/arc/molecule/filtration.py +++ b/arc/molecule/filtration.py @@ -14,6 +14,8 @@ which is quite like http://www.chem.ucla.edu/~harding/IGOC/R/resonance_contributor_preference_rules.html) """ +import itertools + from arc.common import get_logger from arc.exceptions import ResonanceError from arc.molecule.element import PeriodicSystem @@ -206,16 +208,18 @@ def charge_filtration(filtered_list, charge_span_list): filtered_list = [filtered_mol for index, filtered_mol in enumerate(filtered_list) if charge_span_list[index] == min_charge_span] # the minimal charge span layer # Find the radical and multiple bond sites in all filtered_list structures: - rad_sorting_list = [] # sorting_label for radical sites - mul_bond_sorting_list = [] # sorting_label for multiple bond sites in the form of (atom1, atom2) tuples + rad_sorting_list = [] # atom indices for radical sites + mul_bond_sorting_list = [] # atom indices for multiple bond sites in the form of (atom1, atom2) tuples for mol in filtered_list: - for atom in mol.vertices: - if atom.radical_electrons and int(atom.sorting_label) not in rad_sorting_list: - rad_sorting_list.append(int(atom.sorting_label)) + indices = get_atom_indices(mol) + for index1, atom in enumerate(mol.vertices): + if atom.radical_electrons and index1 not in rad_sorting_list: + rad_sorting_list.append(index1) for atom2, bond in atom.bonds.items(): + index2 = indices[id(atom2)] # check if bond is multiple, store only from one side (atom1 < atom2) for consistency - if atom2.sorting_label > atom.sorting_label and bond.is_double() or bond.is_triple(): - mul_bond_sorting_list.append((int(atom.sorting_label), int(atom2.sorting_label))) + if index2 > index1 and (bond.is_double() or bond.is_triple()): + mul_bond_sorting_list.append((index1, index2)) # Find unique radical and multiple bond sites in charged_list and append to unique_charged_list: unique_charged_list = [] for mol in charged_list: @@ -240,16 +244,28 @@ def charge_filtration(filtered_list, charge_span_list): return filtered_list +def get_atom_indices(mol): + """ + Return a mapping of ``id(atom)`` to the atom's position in ``mol.vertices``. + + The mapping identifies an atom consistently across the structures of a single species only + while those structures share a common atom order. + """ + return {id(atom): index for index, atom in enumerate(mol.vertices)} + + def find_unique_sites_in_charged_list(mol, rad_sorting_list, mul_bond_sorting_list): """ A helper function for reactive site discovery in charged species """ - for atom in mol.vertices: - if atom.radical_electrons and int(atom.sorting_label) not in rad_sorting_list: + indices = get_atom_indices(mol) + for index1, atom in enumerate(mol.vertices): + if atom.radical_electrons and index1 not in rad_sorting_list: return [mol] for atom2, bond in atom.bonds.items(): - if (atom2.sorting_label > atom.sorting_label and (bond.is_double() or bond.is_triple()) - and (int(atom.sorting_label), int(atom2.sorting_label)) not in mul_bond_sorting_list + index2 = indices[id(atom2)] + if (index2 > index1 and (bond.is_double() or bond.is_triple()) + and (index1, index2) not in mul_bond_sorting_list and not (atom.is_sulfur() and atom2.is_sulfur())): # We check that both atoms aren't S, otherwise we get [S.-]=[S.+] as a structure of S2 triplet return [mol] @@ -301,6 +317,7 @@ def stabilize_charges_by_proximity(mol_list): """ Only keep structures that obey the charge proximity rule. Opposite charges will be as close as possible to one another, and vice versa. + Charged atoms in different molecular fragments have no path between them and are not counted. """ indices_to_pop = [] charge_distance_list = [] # indices match mol_list @@ -308,29 +325,31 @@ def stabilize_charges_by_proximity(mol_list): # Try finding well-defined pairs of formally-charged atoms to apply the proximity principle # (opposite charges will be as close as possible to one another, and vice versa) cumulative_opposite_charge_distance = cumulative_similar_charge_distance = 0 - for atom1 in mol.vertices: - if atom1.charge: - for atom2 in mol.vertices: - if atom2.charge and atom2.sorting_label > atom1.sorting_label: - # found two charged atoms - if (atom1.charge > 0) ^ (atom2.charge > 0): # xor - # they have opposing signs when ONLY one is positive - cumulative_opposite_charge_distance += len(find_shortest_path(atom1, atom2)) - else: - # they have similar signs - cumulative_similar_charge_distance += len(find_shortest_path(atom1, atom2)) + charged_atoms = [atom for atom in mol.vertices if atom.charge] + for atom1, atom2 in itertools.combinations(charged_atoms, 2): + # found two charged atoms + path = find_shortest_path(atom1, atom2) + if path is None: + continue + if (atom1.charge > 0) ^ (atom2.charge > 0): # xor + # they have opposing signs when ONLY one is positive + cumulative_opposite_charge_distance += len(path) + else: + # they have similar signs + cumulative_similar_charge_distance += len(path) charge_distance_list.append([cumulative_opposite_charge_distance, cumulative_similar_charge_distance]) - min_cumulative_opposite_charge_distance = min([distances[0] for distances in charge_distance_list] - or [0]) # in Python 3 use `min(list, default=0)` + min_cumulative_opposite_charge_distance = min((distances[0] for distances in charge_distance_list), + default=0) for i, distances in enumerate(charge_distance_list): # after generating the charge_distance_list, iterate through it and mark structures to pop if distances[0] > min_cumulative_opposite_charge_distance: indices_to_pop.append(i) - max_cumulative_similar_charge_distance = max([distances[1] for i, distances in - enumerate(charge_distance_list) if i not in indices_to_pop] or [0]) + max_cumulative_similar_charge_distance = max((distances[1] for i, distances in + enumerate(charge_distance_list) if i not in indices_to_pop), + default=0) for i, distances in enumerate(charge_distance_list): - if distances[0] < max_cumulative_similar_charge_distance: + if distances[1] < max_cumulative_similar_charge_distance: indices_to_pop.append(i) for i in reversed(range(len(mol_list))): # pop starting from the end, so indices won't change if i in indices_to_pop: diff --git a/arc/molecule/filtration_test.py b/arc/molecule/filtration_test.py index 6c583fde90..9409b4f2fa 100644 --- a/arc/molecule/filtration_test.py +++ b/arc/molecule/filtration_test.py @@ -4,7 +4,8 @@ import unittest from arc.molecule.filtration import (get_octet_deviation_list, get_octet_deviation, filter_structures, - charge_filtration, get_charge_span_list, aromaticity_filtration) + charge_filtration, get_atom_indices, get_charge_span_list, + aromaticity_filtration) from arc.molecule.molecule import Molecule from arc.molecule.resonance import generate_resonance_structures, analyze_molecule @@ -126,7 +127,7 @@ def test_radical_site(self): Molecule().from_adjacency_list(adj3)] for mol in mol_list: - mol.update() # the charge_filtration uses the atom.sorting_label attribute + mol.update() filtered_list = charge_filtration(mol_list, get_charge_span_list(mol_list)) self.assertEqual(len(filtered_list), 2) @@ -213,7 +214,7 @@ def test_electronegativity(self): Molecule().from_adjacency_list(adj7)] for mol in mol_list: - mol.update() # the charge_filtration uses the atom.sorting_label attribute + mol.update() filtered_list = charge_filtration(mol_list, get_charge_span_list(mol_list)) self.assertEqual(len(filtered_list), 4) @@ -301,6 +302,53 @@ def test_aromaticity(self): filtered_list = aromaticity_filtration(mol_list, analyze_molecule(mol_list[0])) self.assertEqual(len(filtered_list), 3) + def test_get_atom_indices(self): + """Test that atoms are mapped to their position in mol.vertices""" + mol = Molecule().from_smiles('[O]N=O') + indices = get_atom_indices(mol) + self.assertEqual(len(indices), len(mol.vertices)) + for index, atom in enumerate(mol.vertices): + self.assertEqual(indices[id(atom)], index) + + def test_charge_filtration_independent_of_atom_order(self): + """Test that the charge filtration heuristics return the same structures for either atom order""" + for smiles, expected, aromatic in (('[O]c1ccc([N+](=O)[O-])cc1', 11, 2), + ('[CH2]c1ccc([N+](=O)[O-])cc1', 11, 2), + ('[O]c1ccccc1[N+](=O)[O-]', 11, 2), + ('[O]c1cccc([N+](=O)[O-])c1', 9, 2), + ('[O]N=O', 4, 0), + ('C=N[O]', 3, 0), + ('NC=O', 2, 0), + ): + sorted_list = generate_resonance_structures(Molecule().from_smiles(smiles), + keep_isomorphic=True, save_order=False) + saved_list = generate_resonance_structures(Molecule().from_smiles(smiles), + keep_isomorphic=True, save_order=True) + self.assertEqual(len(sorted_list), expected, msg=smiles) + self.assertEqual(len(saved_list), expected, msg=smiles) + for mol_list in (sorted_list, saved_list): + self.assertEqual(sum(1 for mol in mol_list + if any(bond.is_benzene() for bond in mol.get_all_edges())), + aromatic, msg=smiles) + + def test_charge_filtration_of_disconnected_ions(self): + """Test that a species whose charges sit in disconnected components is filtered + + Charged atoms in different molecular fragments have no path between them, so the charge + proximity heuristic does not count them, and a species whose only opposite-charge pairs + are of that kind is not filtered out entirely. + """ + for save_order in (False, True): + for smiles, expected in (('[Li+].[OH-]', 1), + ('[NH4+].[O-]N=O', 1), + ('O=C([O-])[O-].[NH4+].[NH4+]', 1), + ('[Li+].[Li+].[O-][O-]', 1), + ): + mol_list = generate_resonance_structures(Molecule().from_smiles(smiles), + save_order=save_order) + self.assertEqual(len(mol_list), expected, msg=smiles) + self.assertTrue(mol_list[0].get_charge_span()) + if __name__ == '__main__': unittest.main(testRunner=unittest.TextTestRunner(verbosity=2)) diff --git a/arc/molecule/resonance.py b/arc/molecule/resonance.py index 87ee756ddb..1c71ad061e 100644 --- a/arc/molecule/resonance.py +++ b/arc/molecule/resonance.py @@ -271,6 +271,8 @@ def _generate_resonance_structures(mol_list, method_list, keep_isomorphic=False, if True, only remove structures that give is_identical=True copy: if False, append new resonance structures to the input list (default) if True, make a new list with all the resonance structures + save_order: if True, the atom order is maintained; it is forwarded to every method in + ``method_list`` that is listed in ``_SAVE_ORDER_METHODS`` """ cython.declare(index=cython.int, molecule=Graph, new_mol_list=list, new_mol=Graph, mol=Graph, input_charge=cython.int, x=Vertex) @@ -296,7 +298,10 @@ def _generate_resonance_structures(mol_list, method_list, keep_isomorphic=False, charge_span = molecule.get_charge_span() if octet_deviation <= min_octet_deviation + 2 and charge_span <= min_charge_span + 1: for method in method_list: - new_mol_list.extend(method(molecule)) + if method in _SAVE_ORDER_METHODS: + new_mol_list.extend(method(molecule, save_order=save_order)) + else: + new_mol_list.extend(method(molecule)) if octet_deviation < min_octet_deviation: # update min_octet_deviation to make this criterion tighter min_octet_deviation = octet_deviation @@ -960,6 +965,13 @@ def generate_clar_structures(mol, save_order=False): return mol_list +_SAVE_ORDER_METHODS = ( + generate_optimal_aromatic_resonance_structures, + generate_aromatic_resonance_structure, + generate_clar_structures, +) + + # Define this helper function at the module level (outside the cpdef method): def _aromatic_ring_sort_key(ring): sum_ids = 0 diff --git a/arc/molecule/resonance_test.py b/arc/molecule/resonance_test.py index 5668e08378..1200a9cd29 100644 --- a/arc/molecule/resonance_test.py +++ b/arc/molecule/resonance_test.py @@ -1,11 +1,29 @@ #!/usr/bin/env python3 # encoding: utf-8 +import inspect import unittest from arc.molecule.molecule import Molecule -from arc.molecule.resonance import (_clar_optimization, _clar_transformation, generate_clar_structures, - generate_kekule_structure, generate_optimal_aromatic_resonance_structures, generate_resonance_structures) +from arc.molecule.resonance import (_SAVE_ORDER_METHODS, _clar_optimization, _clar_transformation, + generate_aromatic_resonance_structure, generate_clar_structures, generate_kekule_structure, + generate_optimal_aromatic_resonance_structures, generate_resonance_structures, + populate_resonance_algorithms) + + +def resonance_fingerprint(mol): + """ + Return an atom order independent fingerprint of ``mol``, keyed on the atom IDs. + + The fingerprint is a tuple of a sorted tuple of + (id, element symbol, charge, radical electrons, lone pairs) per atom, + and a sorted tuple of (lower id, higher id, bond order) per bond. + """ + atoms = tuple(sorted((atom.id, atom.element.symbol, atom.charge, atom.radical_electrons, atom.lone_pairs) + for atom in mol.atoms)) + bonds = tuple(sorted((min(bond.atom1.id, bond.atom2.id), max(bond.atom1.id, bond.atom2.id), bond.order) + for bond in mol.get_all_edges())) + return atoms, bonds class ResonanceTest(unittest.TestCase): @@ -1160,6 +1178,66 @@ def test_resonance_without_changing_atom_order2(self): atom2_nb = {nb.id for nb in list(atom2.bonds.keys())} self.assertEqual(atom1_nb, atom2_nb) + def test_resonance_of_polycyclic_aromatics_without_changing_atom_order(self): + """Test that save_order is honored by the Clar structures of polycyclic aromatic molecules""" + for smiles, num_structures in [('c1ccc2ccccc2c1', 3), ('C#Cc1ccc2ccccc2c1', 4), ('c1ccc2cc3ccccc3cc2c1', 4)]: + mol = Molecule(smiles=smiles) + mol.assign_atom_ids() + original_order = [atom.id for atom in mol.atoms] + neighbors = [{neighbor.id for neighbor in atom.bonds.keys()} for atom in mol.atoms] + + ordered = generate_resonance_structures(mol.copy(deep=True), save_order=True) + self.assertEqual(len(ordered), num_structures) + for res_mol in ordered: + self.assertEqual([atom.id for atom in res_mol.atoms], original_order) + for index, atom in enumerate(res_mol.atoms): + self.assertEqual(atom.element.symbol, mol.atoms[index].element.symbol) + self.assertEqual({neighbor.id for neighbor in atom.bonds.keys()}, neighbors[index]) + + sorted_ = generate_resonance_structures(mol.copy(deep=True), save_order=False) + self.assertEqual(len(sorted_), num_structures) + self.assertTrue(any([atom.id for atom in res_mol.atoms] != original_order + for res_mol in sorted_)) + for res_mol in ordered: + self.assertTrue(any(res_mol.copy(deep=True).is_isomorphic(other.copy(deep=True)) + for other in sorted_)) + + kept_ordered = generate_resonance_structures(mol.copy(deep=True), keep_isomorphic=True, save_order=True) + kept_sorted = generate_resonance_structures(mol.copy(deep=True), keep_isomorphic=True, save_order=False) + self.assertEqual({resonance_fingerprint(res_mol) for res_mol in kept_ordered}, + {resonance_fingerprint(res_mol) for res_mol in kept_sorted}) + + def test_clar_structures_of_charge_separated_polycyclic_aromatics(self): + """Test that the Clar structures of a charge separated polycyclic aromatic survive filtration""" + mol = Molecule(smiles='[O-][N+](=O)c1cccc2ccccc12') + mol.assign_atom_ids() + + non_isomorphic = generate_resonance_structures(mol.copy(deep=True), keep_isomorphic=False) + identical = generate_resonance_structures(mol.copy(deep=True), keep_isomorphic=True) + + self.assertEqual(len(non_isomorphic), 4) + self.assertEqual(len(identical), 7) + self.assertGreaterEqual(len(identical), len(non_isomorphic)) + + ordered = generate_resonance_structures(mol.copy(deep=True), keep_isomorphic=True, save_order=True) + sorted_ = generate_resonance_structures(mol.copy(deep=True), keep_isomorphic=True, save_order=False) + self.assertEqual({resonance_fingerprint(res_mol) for res_mol in ordered}, + {resonance_fingerprint(res_mol) for res_mol in sorted_}) + + +class SaveOrderMethodsTest(unittest.TestCase): + """ + Contains unit tests for _SAVE_ORDER_METHODS. + """ + + def test_save_order_aware_algorithms(self): + """Test that _SAVE_ORDER_METHODS lists exactly the dispatchable algorithms taking save_order""" + dispatchable = set(populate_resonance_algorithms()) + dispatchable.add(generate_aromatic_resonance_structure) + accepts_save_order = {method for method in dispatchable + if 'save_order' in inspect.signature(method).parameters} + self.assertEqual(set(_SAVE_ORDER_METHODS), accepts_save_order) + class ClarTest(unittest.TestCase): """