Fix resonance save order

rmgpy/molecule/filtration.py

@@ -43,6 +43,7 @@
 which is quite like http://www.chem.ucla.edu/~harding/IGOC/R/resonance_contributor_preference_rules.html)
 """
 
+import itertools
 import logging
 
 from rmgpy.exceptions import ResonanceError
@@ -238,16 +239,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 = []  # sortingLabels for radical sites
-        mul_bond_sorting_list = []  # sortingLabels for multiple bind sites in the form of (atom1,atom2) tuples
+        rad_sorting_list = []  # atom indices for radical sites
+        mul_bond_sorting_list = []  # atom indices for multiple bind 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.edges.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:
@@ -272,16 +275,31 @@ 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``.
+
+    Replaces ``Vertex.sorting_label``, which the isomorphism machinery leaves either unset or
+    holding a stale permutation. Position identifies an atom across the structures of one species
+    only while they share an atom order, which ``save_order`` guarantees and ``sort_atoms`` does
+    not; that limitation is inherited, since the sorting labels equalled the position in exactly
+    the case where the structures had been sorted.
+    """
+    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.edges.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]
@@ -340,17 +358,19 @@ 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:
+                # charges in disconnected components, e.g. an ionic pair, have no through-bond distance
+                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),
@@ -363,7 +383,7 @@ def stabilize_charges_by_proximity(mol_list):
                                                   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.add(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:

rmgpy/molecule/resonance.py

@@ -56,6 +56,7 @@
 """
 
 import logging
+from functools import partial
 from operator import attrgetter
 
 import cython
@@ -251,11 +252,13 @@ def generate_resonance_structures(mol, clar_structures=True, keep_isomorphic=Fal
                                            keep_isomorphic=keep_isomorphic,
                                            save_order=save_order)
         if features['isPolycyclicAromatic'] and clar_structures:
-            _generate_resonance_structures(mol_list, [generate_clar_structures],
+            _generate_resonance_structures(mol_list,
+                                           [partial(generate_clar_structures, save_order=save_order)],
                                            keep_isomorphic=keep_isomorphic,
                                            save_order=save_order)
         else:
-            _generate_resonance_structures(mol_list, [generate_aromatic_resonance_structure],
+            _generate_resonance_structures(mol_list,
+                                           [partial(generate_aromatic_resonance_structure, save_order=save_order)],
                                            keep_isomorphic=keep_isomorphic,
                                            save_order=save_order)
 
@@ -282,6 +285,7 @@ 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 input list (default)
                             if True, make a new list with all of the resonance structures
+        save_order          if True, preserve the atom order of the input molecules
     """
     cython.declare(index=cython.int, molecule=Graph, new_mol_list=list, new_mol=Graph, mol=Graph,
                    input_charge=cython.int, x=Vertex)

test/rmgpy/molecule/filtrationTest.py

@@ -158,7 +158,7 @@ def test_radical_site(self):
         ]
 
         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))
         assert len(filtered_list) == 2
@@ -247,7 +247,7 @@ def test_electronegativity(self):
         ]
 
         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))
         assert len(filtered_list) == 4
@@ -334,3 +334,49 @@ def test_aromaticity(self):
 
         filtered_list = aromaticity_filtration(mol_list, analyze_molecule(mol_list[0]))
         assert len(filtered_list) == 3
+
+    def test_charge_filtration_independent_of_atom_order(self):
+        """Test that the charge filtration heuristics do not depend on the atom order
+
+        The heuristics used to key on Vertex.sorting_label, which the isomorphism machinery
+        leaves either unset or holding a stale permutation, so requesting save_order
+        discarded structures -- for the aromatics below, every aromatic structure.
+        """
+        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
+            )
+            assert len(sorted_list) == len(saved_list) == expected
+            for mol_list in (sorted_list, saved_list):
+                assert sum(1 for mol in mol_list if any(bond.is_benzene() for bond in mol.get_all_edges())) == aromatic
+
+    def test_charge_filtration_of_disconnected_ions(self):
+        """Test that salts, whose charges straddle disconnected fragments, survive filtration
+
+        find_shortest_path has no path between fragments, so no opposite-charge distance
+        accumulates. The proximity rule must then compare the like-charge distance against its
+        own maximum: comparing the opposite-charge distance against it pops every structure,
+        because the two are unrelated quantities, and filtration is left with nothing.
+        """
+        for smiles, expected in (
+            ("[Li+].[OH-]", 1),
+            ("[Li+].[O-]C=O", 2),
+            ("[Li+].[Li+].[O-][O-]", 1),
+            ("[Li+].[Li+].[O-]C(=O)[O-]", 3),
+        ):
+            for save_order in (False, True):
+                mol_list = generate_resonance_structures(
+                    Molecule().from_smiles(smiles), keep_isomorphic=True, save_order=save_order
+                )
+                assert len(mol_list) == expected

test/rmgpy/molecule/resonanceTest.py

@@ -1362,6 +1362,128 @@ def test_resonance_without_changing_atom_order2(self):
                 atom2_nb = {nb.id for nb in list(atom2.bonds.keys())}
                 assert atom1_nb == atom2_nb
 
+    def test_resonance_without_changing_atom_order3(self):
+        """Test generating resonance structures for polycyclic aromatics without changing the atom order"""
+        mol = Molecule().from_adjacency_list(
+            """
+1  O u0 p2 c0 {2,S} {12,S}
+2  C u0 p0 c0 {1,S} {3,S} {11,D}
+3  C u0 p0 c0 {2,S} {4,D} {13,S}
+4  C u0 p0 c0 {3,D} {5,S} {14,S}
+5  C u0 p0 c0 {4,S} {6,D} {15,S}
+6  C u0 p0 c0 {5,D} {7,S} {11,S}
+7  C u0 p0 c0 {6,S} {8,D} {16,S}
+8  C u0 p0 c0 {7,D} {9,S} {17,S}
+9  C u0 p0 c0 {8,S} {10,D} {18,S}
+10 C u0 p0 c0 {9,D} {11,S} {19,S}
+11 C u0 p0 c0 {2,D} {6,S} {10,S}
+12 H u0 p0 c0 {1,S}
+13 H u0 p0 c0 {3,S}
+14 H u0 p0 c0 {4,S}
+15 H u0 p0 c0 {5,S}
+16 H u0 p0 c0 {7,S}
+17 H u0 p0 c0 {8,S}
+18 H u0 p0 c0 {9,S}
+19 H u0 p0 c0 {10,S}
+"""
+        )
+
+        res_mols = mol.copy(deep=True).generate_resonance_structures(save_order=True)
+
+        # Assign atom ids
+        for molecule in [mol] + res_mols:
+            for idx, atom in enumerate(molecule.atoms):
+                atom.id = idx
+
+        # Compare the atom symbols and their nearest neighbors
+        for res_mol in res_mols:
+            for atom1, atom2 in zip(mol.atoms, res_mol.atoms):
+                assert atom1.element.symbol == atom2.element.symbol
+                atom1_nb = {nb.id for nb in list(atom1.bonds.keys())}
+                atom2_nb = {nb.id for nb in list(atom2.bonds.keys())}
+                assert atom1_nb == atom2_nb
+
+    def test_resonance_with_save_order_keeps_clar_structures(self):
+        """Test that save_order does not discard the Clar structures of a charged polycyclic aromatic
+
+        The charge filtration heuristics compare ``atom.sorting_label``, which is left unset until
+        something sorts the molecule. When the Clar structures were the only ones whose atoms had
+        been sorted, ``stabilize_charges_by_proximity`` popped exactly them. The default path is
+        pinned here too.
+        """
+        adjlist = """
+1  O u0 p3 c-1 {2,S}
+2  N u0 p0 c+1 {1,S} {3,D} {4,S}
+3  O u0 p2 c0 {2,D}
+4  C u0 p0 c0 {2,S} {5,S} {13,D}
+5  C u0 p0 c0 {4,S} {6,D} {14,S}
+6  C u0 p0 c0 {5,D} {7,S} {15,S}
+7  C u0 p0 c0 {6,S} {8,D} {16,S}
+8  C u0 p0 c0 {7,D} {9,S} {13,S}
+9  C u0 p0 c0 {8,S} {10,D} {17,S}
+10 C u0 p0 c0 {9,D} {11,S} {18,S}
+11 C u0 p0 c0 {10,S} {12,D} {19,S}
+12 C u0 p0 c0 {11,D} {13,S} {20,S}
+13 C u0 p0 c0 {4,D} {8,S} {12,S}
+14 H u0 p0 c0 {5,S}
+15 H u0 p0 c0 {6,S}
+16 H u0 p0 c0 {7,S}
+17 H u0 p0 c0 {9,S}
+18 H u0 p0 c0 {10,S}
+19 H u0 p0 c0 {11,S}
+20 H u0 p0 c0 {12,S}
+"""
+
+        def count_aromatic(mol_list):
+            return sum(1 for mol in mol_list if any(bond.is_benzene() for bond in mol.get_all_edges()))
+
+        def distinct(mol_list):
+            # is_isomorphic re-sorts its operands, so compare copies and leave mol_list untouched
+            reps = []
+            for mol in mol_list:
+                candidate = mol.copy(deep=True)
+                if not any(candidate.is_isomorphic(rep) for rep in reps):
+                    reps.append(candidate)
+            return reps
+
+        # 1-nitronaphthalene
+        sorted_list = generate_resonance_structures(
+            Molecule().from_adjacency_list(adjlist), keep_isomorphic=True, save_order=False
+        )
+        saved_list = generate_resonance_structures(
+            Molecule().from_adjacency_list(adjlist), keep_isomorphic=True, save_order=True
+        )
+        unique_list = generate_resonance_structures(
+            Molecule().from_adjacency_list(adjlist), keep_isomorphic=False, save_order=True
+        )
+        unique_sorted_list = generate_resonance_structures(
+            Molecule().from_adjacency_list(adjlist), keep_isomorphic=False, save_order=False
+        )
+
+        assert len(distinct(saved_list)) == len(distinct(sorted_list)) == 4
+        assert count_aromatic(distinct(saved_list)) == count_aromatic(distinct(sorted_list)) == 3
+
+        # the default path: the delocalized structure, the two Clar sextets, and the unreactive input
+        assert len(unique_list) == len(unique_sorted_list) == 4
+        assert sum(1 for mol in unique_list if mol.reactive) == 3
+        assert count_aromatic(unique_list) == 3
+
+    def test_resonance_with_save_order_non_polycyclic_aromatic(self):
+        """Test that save_order does not change the structures of a non-polycyclic aromatic
+
+        RMG counts only six-membered all-carbon rings as aromatic, so these reach
+        generate_aromatic_resonance_structure rather than generate_clar_structures, covering the
+        other save_order-aware dispatch. Indole is bicyclic but only one of its rings counts.
+        """
+        for smiles, expected in (("[CH2]c1ccccc1", 5), ("[O]c1ccc(O)cc1", 5), ("c1ccc2[nH]ccc2c1", 2)):
+            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
+            )
+            assert len(saved_list) == len(sorted_list) == expected
+
     def test_adsorbate_resonance_cc1(self):
         """Test if all three resonance structures for X#CC#X are generated"""
         adjlist = """