Stop the charge filtration heuristics keying on sorting_label

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,36 +317,39 @@ 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
     for i, mol in enumerate(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:

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))

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

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):
     """