Allow addLabelsForReaction to properly handle identical references

If a reaction object, which had Molecule or Species objects with
identical references, was input to addLabelsForReaction, the
conflicting references would not allow the same molecule to hold
multiple labels, leading to errors.

This creates a deep copy of the Species object when identical
references are found, since seprate references are necessary for
the labeling.

This commit also adds three test methods for addLabelsForReaction
which previously had no test methods.

Author: goldmanm

rmgpy/data/kinetics/family.py

@@ -2132,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]))
@@ -2152,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

@@ -897,3 +897,128 @@ def test_generate_reactions_from_libraries2(self):
         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)
+