Merge pull request #1276 from ReactionMechanismGenerator/fix_calculate_degeneracy_bug

Update calculateDegeneracy to deal with identical reactants

Author: mliu49

rmgpy/data/kinetics/common.py

@@ -156,8 +156,8 @@ def filter_reactions(reactants, products, reactionList):
     warnings.warn("The filter_reactions method is no longer used and may be removed in a future version.", DeprecationWarning)
 
     # Convert from molecules to species and generate resonance isomers.
-    reactants = ensure_species(reactants, resonance=True)
-    products = ensure_species(products, resonance=True)
+    ensure_species(reactants, resonance=True)
+    ensure_species(products, resonance=True)
 
     reactions = reactionList[:]
 
@@ -197,11 +197,10 @@ def filter_reactions(reactants, products, reactionList):
 
 def ensure_species(input_list, resonance=False, keepIsomorphic=False):
     """
-    Given an input list of molecules or species, return a list with only
-    species objects.
+    The input list of :class:`Species` or :class:`Molecule` objects is modified
+    in place to only have :class:`Species` objects. Returns None.
     """
-    output_list = []
-    for item in input_list:
+    for index, item in enumerate(input_list):
         if isinstance(item, Molecule):
             new_item = Species(molecule=[item])
         elif isinstance(item, Species):
@@ -210,9 +209,7 @@ def ensure_species(input_list, resonance=False, keepIsomorphic=False):
             raise TypeError('Only Molecule or Species objects can be handled.')
         if resonance:
             new_item.generate_resonance_structures(keepIsomorphic=keepIsomorphic)
-        output_list.append(new_item)
-
-    return output_list
+        input_list[index] = new_item
 
 
 def generate_molecule_combos(input_species):
@@ -231,13 +228,15 @@ def generate_molecule_combos(input_species):
 
 def ensure_independent_atom_ids(input_species, resonance=True):
     """
-    Given a list or tuple of :class:`Species` objects, ensure that atom ids are
-    independent across all of the species. Optionally, the `resonance` argument
-    can be set to False to not generate resonance structures.
+    Given a list or tuple of :class:`Species` or :class:`Molecule` objects,
+    ensure that atom ids are independent.
+    The `resonance` argument can be set to False to not generate
+    resonance structures.
 
-    Modifies the input species in place, nothing is returned.
+    Modifies the list in place (replacing :class:`Molecule` with :class:`Species`).
+    Returns None.
     """
-
+    ensure_species(input_species, resonance=resonance)
     # Method to check that all species' atom ids are different
     def independent_ids():
         num_atoms = 0

rmgpy/data/kinetics/database.py

@@ -364,7 +364,7 @@ def generate_reactions(self, reactants, products=None, only_families=None, reson
         reactionList = []
         if only_families is None:
             reactionList.extend(self.generate_reactions_from_libraries(reactants, products))
-        reactionList.extend(self.generate_reactions_from_families(reactants, products, only_families=None, resonance=True))
+        reactionList.extend(self.generate_reactions_from_families(reactants, products, only_families=None, resonance=resonance))
         return reactionList
 
     def generate_reactions_from_libraries(self, reactants, products=None):
@@ -386,7 +386,7 @@ def generate_reactions_from_library(self, library, reactants, products=None):
         provided `reactants`, which can be either :class:`Molecule` objects or
         :class:`Species` objects.
         """
-        reactants = ensure_species(reactants)
+        ensure_species(reactants)
 
         reaction_list = []
         for entry in library.entries.values():
@@ -434,10 +434,9 @@ def generate_reactions_from_families(self, reactants, products=None, only_famili
             elif reactants[0].isIsomorphic(reactants[1]):
                 same_reactants = True
 
-        # Convert to Species objects if necessary
-        reactants = ensure_species(reactants)
-
-        # Label reactant atoms for proper degeneracy calculation
+        # Label reactant atoms for proper degeneracy calculation (cannot be in tuple)
+        if isinstance(reactants, tuple):
+            reactants = list(reactants)
         ensure_independent_atom_ids(reactants, resonance=resonance)
 
         combos = generate_molecule_combos(reactants)

rmgpy/data/kinetics/family.py

@@ -46,7 +46,8 @@
 from rmgpy.molecule.resonance import generate_aromatic_resonance_structures
 from rmgpy.species import Species
 
-from .common import saveEntry, ensure_species, find_degenerate_reactions, generate_molecule_combos
+from .common import saveEntry, ensure_species, find_degenerate_reactions, generate_molecule_combos,\
+                    ensure_independent_atom_ids
 from .depository import KineticsDepository
 from .groups import KineticsGroups
 from .rules import KineticsRules
@@ -1545,22 +1546,28 @@ def calculateDegeneracy(self, reaction):
         `ignoreSameReactants= True` to this method.
         """
         reaction.degeneracy = 1
-
-        # find combinations of resonance isomers
-        specReactants = ensure_species(reaction.reactants, resonance=True, keepIsomorphic=True)
-        molecule_combos = generate_molecule_combos(specReactants)
+        # Check if the reactants are the same
+        # If they refer to the same memory address, then make a deep copy so
+        # they can be manipulated independently
+        reactants = reaction.reactants
+        same_reactants = False
+        if len(reactants) == 2:
+            if reactants[0] is reactants[1]:
+                reactants[1] = reactants[1].copy(deep=True)
+                same_reactants = True
+            elif reactants[0].isIsomorphic(reactants[1]):
+                same_reactants = True
+
+        # Label reactant atoms for proper degeneracy calculation
+        ensure_independent_atom_ids(reactants, resonance=True)
+        molecule_combos = generate_molecule_combos(reactants)
 
         reactions = []
         for combo in molecule_combos:
             reactions.extend(self.__generateReactions(combo, products=reaction.products, forward=True))
 
-        # Check if the reactants are the same
-        sameReactants = False
-        if len(specReactants) == 2 and specReactants[0].isIsomorphic(specReactants[1]):
-            sameReactants = True
-
         # remove degenerate reactions
-        reactions = find_degenerate_reactions(reactions, sameReactants, kinetics_family=self)
+        reactions = find_degenerate_reactions(reactions, same_reactants, kinetics_family=self)
 
         # remove reactions with different templates (only for TemplateReaction)
         if isinstance(reaction, TemplateReaction):
@@ -1692,7 +1699,7 @@ def __generateReactions(self, reactants, products=None, forward=True, prod_reson
         # If products is given, remove reactions from the reaction list that
         # don't generate the given products
         if products is not None:
-            products = ensure_species(products, resonance=prod_resonance)
+            ensure_species(products, resonance=prod_resonance)
 
             rxnList0 = rxnList[:]
             rxnList = []
@@ -2125,6 +2132,15 @@ def addAtomLabelsForReaction(self, reaction, output_with_resonance = True):
         # make sure we start with reaction with species objects
         reaction.ensure_species(reactant_resonance=False, product_resonance=False)
 
+        reactants = reaction.reactants
+        products = reaction.products
+        # ensure all species are independent references
+        if len(reactants + products) > len(set([id(s) for s in reactants + products])):
+            logging.debug('Copying reactants and products for reaction {} since they have identical species references'.format(reaction))
+            # not all species are independent
+            reactants = [s.copy(deep=True) for s in reactants]
+            products = [s.copy(deep=True) for s in products]
+
         # get all possible pairs of resonance structures
         reactant_pairs = list(itertools.product(*[s.molecule for s in reaction.reactants]))
         product_pairs = list(itertools.product(*[s.molecule for s in reaction.products]))
@@ -2145,20 +2161,22 @@ def addAtomLabelsForReaction(self, reaction, output_with_resonance = True):
 
         # place the molecules in reaction's species object
         # this prevents overwriting of attributes of species objects by this method
-        for species in reaction.products:
+        for index, species in enumerate(products):
             for labeled_molecule in labeled_products:
                 if species.isIsomorphic(labeled_molecule):
                     species.molecule = [labeled_molecule]
+                    reaction.products[index] = species
                     break
             else:
-                raise ActionError('Could not find isomorphic molecule to fit the original product {}'.format(species))
-        for species in reaction.reactants:
+                raise ActionError('Could not find isomorphic molecule to fit the original product {} from reaction {}'.format(species, reaction))
+        for index, species in enumerate(reactants):
             for labeled_molecule in labeled_reactants:
                 if species.isIsomorphic(labeled_molecule):
                     species.molecule = [labeled_molecule]
+                    reaction.reactants[index] = species
                     break
             else:
-                raise ActionError('Could not find isomorphic molecule to fit the original reactant {}'.format(species))
+                raise ActionError('Could not find isomorphic molecule to fit the original reactant {} from reaction {}'.format(species, reaction))
 
         if output_with_resonance:
             # convert the molecules to species objects with resonance structures

rmgpy/data/kinetics/kineticsTest.py

@@ -688,6 +688,21 @@ def test_ensure_independent_atom_ids(self):
         # checks second resonance structure id
         self.assertNotEqual(s2.molecule[1].atoms[0].id, -1)
 
+    def test_ensure_independent_atom_ids_no_resonance(self):
+        """
+        Ensure ensure_independent_atom_ids does not generate resonance
+        """
+        s1 = Species().fromSMILES('CCC')
+        s2 = Species().fromSMILES('C=C[CH]C')
+        self.assertEqual(s2.molecule[0].atoms[0].id, -1)
+
+        ensure_independent_atom_ids([s1, s2],resonance=False)
+        # checks resonance structures
+        self.assertEqual(len(s2.molecule),1)
+        # checks that atom ids are changed
+        for atom in s2.molecule[0].atoms:
+            self.assertNotEqual(atom.id, -1)
+
     def testSaveEntry(self):
         """
         tests that save entry can run
@@ -839,25 +854,171 @@ def test_generate_reactions_from_families_product_resonance(self):
         self.assertEqual(len(reaction_list), 1)
         self.assertEqual(reaction_list[0].degeneracy, 2)
 
-        reaction_list = self.database.kinetics.generate_reactions_from_families(reactants, products, only_families=['H_Abstraction'], resonance=False)
 
+
+    def test_generate_reactions_from_families_product_resonance2(self):
+        """Test that we can specify the no product resonance structure when generating reactions"""
+        reactants = [
+            Molecule().fromSMILES('CCC=C'),
+            Molecule().fromSMILES('[H]'),
+        ]
+        products = [
+            Molecule().fromSMILES('CC=C[CH2]'),
+            Molecule().fromSMILES('[H][H]'),
+        ]
+
+        reaction_list = self.database.kinetics.generate_reactions_from_families(reactants, products, only_families=['H_Abstraction'], resonance=False)
         self.assertEqual(len(reaction_list), 0)
 
+        self.assertTrue(isinstance(products[0],Species))
+        self.assertEqual(len(products[0].molecule),1)
+
     def test_generate_reactions_from_libraries(self):
         """Test that we can generate reactions from libraries"""
         reactants = [
             Molecule().fromSMILES('CC=O'),
             Molecule().fromSMILES('[H]'),
         ]
-        products = [
-            Molecule().fromSMILES('[CH2]C=O'),
-            Molecule().fromSMILES('[H][H]'),
-        ]
 
         reaction_list = self.database.kinetics.generate_reactions_from_libraries(reactants)
 
         self.assertEqual(len(reaction_list), 3)
 
+    def test_generate_reactions_from_libraries2(self):
+        """Test that we can generate reactions from libraries specifying products"""
+        reactants = [
+            Molecule().fromSMILES('CC=O'),
+            Molecule().fromSMILES('[H]'),
+        ]
+        products = [
+            Molecule().fromSMILES('[CH2]C=O'),
+            Molecule().fromSMILES('[H][H]'),
+        ]
         reaction_list_2 = self.database.kinetics.generate_reactions_from_libraries(reactants, products)
 
         self.assertEqual(len(reaction_list_2), 1)
+
+    def test_add_atom_labels_for_reaction(self):
+        """Test that addAtomLabelsForReaction can identify reactions with resonance
+        The molecule [CH]=C=C has resonance in this reaction"""
+        from rmgpy.data.rmg import getDB
+        reactants = [
+            Molecule().fromSMILES('C=C=C'),
+            Molecule().fromSMILES('[CH]=C=C'),
+        ]
+        products = [
+            Molecule().fromSMILES('C#C[CH2]'),
+            Molecule().fromSMILES('C#CC'),
+        ]
+        reaction = TemplateReaction(reactants =reactants,
+                                    products = products,
+                                    family = 'H_Abstraction')
+        reaction.ensure_species(reactant_resonance=True, product_resonance=True)
+        family = getDB('kinetics').families['H_Abstraction']
+        family.addAtomLabelsForReaction(reaction, output_with_resonance=False)
+
+        # test that the reaction has labels
+        found_labels = []
+        for species in reaction.reactants:
+            for atom in species.molecule[0].atoms:
+                if atom.label != '':
+                    found_labels.append(atom.label)
+        self.assertEqual(len(found_labels), 3)
+        self.assertIn('*1',found_labels)
+        self.assertIn('*2',found_labels)
+        self.assertIn('*3',found_labels)
+
+        # test for the products too
+        found_labels = []
+        for species in reaction.products:
+            for atom in species.molecule[0].atoms:
+                if atom.label != '':
+                    found_labels.append(atom.label)
+        self.assertEqual(len(found_labels), 3)
+        self.assertIn('*1',found_labels)
+        self.assertIn('*2',found_labels)
+        self.assertIn('*3',found_labels)
+
+    def test_add_atom_labels_for_reaction_2(self):
+        """Test that addAtomLabelsForReaction can identify reactions with identical references
+        The molecule [CH]=C=C has resonance in this reaction"""
+        from rmgpy.data.rmg import getDB
+        s1 = Species().fromSMILES('C=C=C')
+        s2 = Species().fromSMILES('C=C=[CH]')
+        s3 = Species().fromSMILES('C#CC')
+        s2.generate_resonance_structures()
+        reactants = [s1,s2]
+        products = [s2,s3]
+        reaction = TemplateReaction(reactants =reactants,
+                                    products = products,
+                                    family = 'H_Abstraction')
+        family = getDB('kinetics').families['H_Abstraction']
+        print reaction.reactants
+        print reaction.products
+        family.addAtomLabelsForReaction(reaction, output_with_resonance=False)
+
+        # test that the reaction has labels
+        found_labels = []
+        for species in reaction.reactants:
+            for atom in species.molecule[0].atoms:
+                if atom.label != '':
+                    found_labels.append(atom.label)
+        self.assertEqual(len(found_labels), 3,'wrong number of labels found {0}'.format(found_labels))
+        self.assertIn('*1',found_labels)
+        self.assertIn('*2',found_labels)
+        self.assertIn('*3',found_labels)
+
+        # test for the products too
+        found_labels = []
+        for species in reaction.products:
+            for atom in species.molecule[0].atoms:
+                if atom.label != '':
+                    found_labels.append(atom.label)
+        self.assertEqual(len(found_labels), 3)
+        self.assertIn('*1',found_labels)
+        self.assertIn('*2',found_labels)
+        self.assertIn('*3',found_labels)
+
+    def test_add_atom_labels_for_reaction_3(self):
+        """Test that addAtomLabelsForReaction can identify reactions with resonance and isotopes"""
+        from rmgpy.data.rmg import getDB
+        mr0 = Molecule().fromAdjacencyList('1    C u0 p0 c0 i13 {3,D} {4,S} {5,S}\n2 *1 C u0 p0 c0 {3,D} {6,S} {7,S}\n3    C u0 p0 c0 {1,D} {2,D}\n4    H u0 p0 c0 {1,S}\n5    H u0 p0 c0 {1,S}\n6    H u0 p0 c0 {2,S}\n7 *4 H u0 p0 c0 {2,S}\n')
+        mr1a = Molecule().fromAdjacencyList('multiplicity 2\n1    C u0 p0 c0 i13 {2,D} {4,S} {5,S}\n2    C u0 p0 c0 {1,D} {3,D}\n3 *1 C u1 p0 c0 {2,D} {6,S}\n4    H u0 p0 c0 {1,S}\n5    H u0 p0 c0 {1,S}\n6    H u0 p0 c0 {3,S}\n')
+        mr1b = Molecule().fromAdjacencyList('multiplicity 2\n1    C u1 p0 c0 i13 {2,S} {4,S} {5,S}\n2    C u0 p0 c0 {1,S} {3,T}\n3 *1 C u0 p0 c0 {2,T} {6,S}\n4    H u0 p0 c0 {1,S}\n5    H u0 p0 c0 {1,S}\n6    H u0 p0 c0 {3,S}\n')
+        mp1a = Molecule().fromAdjacencyList('multiplicity 2\n1    C u0 p0 c0 {2,D} {4,S} {5,S}\n2    C u0 p0 c0 {1,D} {3,D}\n3 *1 C u1 p0 c0 i13 {2,D} {6,S}\n4    H u0 p0 c0 {1,S}\n5    H u0 p0 c0 {1,S}\n6    H u0 p0 c0 {3,S}\n')
+        mp1b = Molecule().fromAdjacencyList('multiplicity 2\n1    C u1 p0 c0 {2,S} {4,S} {5,S}\n2    C u0 p0 c0 {1,S} {3,T}\n3 *1 C u0 p0 c0 i13 {2,T} {6,S}\n4    H u0 p0 c0 {1,S}\n5    H u0 p0 c0 {1,S}\n6    H u0 p0 c0 {3,S}\n')
+        s1 = Species(molecule = [mr0])
+        s2 = Species(molecule = [mr1a,mr1b])
+        s3 = Species(molecule = [mp1a,mp1b])
+        reactants = [s1,s2]
+        products = [s1,s3]
+        reaction = TemplateReaction(reactants =reactants,
+                                    products = products,
+                                    family = 'H_Abstraction')
+        family = getDB('kinetics').families['H_Abstraction']
+        print reaction.reactants
+        print reaction.products
+        family.addAtomLabelsForReaction(reaction, output_with_resonance=False)
+
+        # test that the reaction has labels
+        found_labels = []
+        for species in reaction.reactants:
+            for atom in species.molecule[0].atoms:
+                if atom.label != '':
+                    found_labels.append(atom.label)
+        self.assertEqual(len(found_labels), 3,'wrong number of labels found {0}'.format(found_labels))
+        self.assertIn('*1',found_labels)
+        self.assertIn('*2',found_labels)
+        self.assertIn('*3',found_labels)
+
+        # test for the products too
+        found_labels = []
+        for species in reaction.products:
+            for atom in species.molecule[0].atoms:
+                if atom.label != '':
+                    found_labels.append(atom.label)
+        self.assertEqual(len(found_labels), 3)
+        self.assertIn('*1',found_labels)
+        self.assertIn('*2',found_labels)
+        self.assertIn('*3',found_labels)
+