added electronLabel for rxn family and recipe

rmgpy/data/kinetics/family.py

@@ -57,7 +57,7 @@
 from rmgpy.kinetics import Arrhenius, SurfaceArrhenius, SurfaceArrheniusBEP, StickingCoefficient, \
                            StickingCoefficientBEP, ArrheniusBM, SurfaceChargeTransfer
 from rmgpy.kinetics.uncertainties import RateUncertainty, rank_accuracy_map
-from rmgpy.molecule import Bond, GroupBond, Group, Molecule
+from rmgpy.molecule import Bond, GroupBond, Group, GroupAtom, Molecule, Atom
 from rmgpy.molecule.atomtype import ATOMTYPES
 from rmgpy.reaction import Reaction, same_species_lists
 from rmgpy.species import Species
@@ -227,6 +227,7 @@ class ReactionRecipe(object):
 
     def __init__(self, actions=None):
         self.actions = actions or []
+        self.electron_label = None
 
     def add_action(self, action):
         """
@@ -350,8 +351,17 @@ def _apply(self, struct, forward, unique):
                             atom.apply_action(['LOSE_RADICAL', label, 1])
                         elif (action[0] == 'LOSE_CHARGE' and forward) or (action[0] == 'GAIN_CHARGE' and not forward):
                             atom.apply_action(['LOSE_CHARGE', label, 1])
-                        elif (action[0] == 'LOSE_CHARGE' and forward) or (action[0] == 'GAIN_CHARGE' and not forward):
-                            atom.apply_action(['LOSE_CHARGE', label, 1])
+                        elif (action[0] == 'GAIN_CHARGE' and forward) or (action[0] == 'LOSE_CHARGE' and not forward):
+                            atom.apply_action(['GAIN_CHARGE', label, 1])
+                            if isinstance(struct, Molecule):
+                                electron = Atom('e')
+                                electron.update_charge()
+                                electron.lone_pairs = 0
+                            elif isinstance(struct, Group):
+                                electron = GroupAtom(atomtype=['e'])
+                                #struct = struct.merge(Group(atoms=[GroupAtom(atomtype=['e'])]))
+                            electron.label = self.electron_label
+                            struct.add_atom(electron)
 
             elif action[0] in ['LOSE_PAIR', 'GAIN_PAIR']:
 
@@ -671,6 +681,7 @@ def load(self, path, local_context=None, global_context=None, depository_labels=
         local_context['reverseMap'] = None
         local_context['reactantNum'] = None
         local_context['productNum'] = None
+        local_context['electronLabel'] = None
         local_context['autoGenerated'] = False
         self.groups = KineticsGroups(label='{0}/groups'.format(self.label))
         logging.debug("Loading kinetics family groups from {0}".format(os.path.join(path, 'groups.py')))
@@ -682,6 +693,7 @@ def load(self, path, local_context=None, global_context=None, depository_labels=
 
         self.reactant_num = local_context.get('reactantNum', None)
         self.product_num = local_context.get('productNum', None)
+        self.electron_label = local_context.get('electronLabel', None)
 
         self.auto_generated = local_context.get('autoGenerated', False)
 
@@ -696,6 +708,7 @@ def load(self, path, local_context=None, global_context=None, depository_labels=
             self.forward_template.products = self.forward_template.reactants[:]
             self.reverse_template = None
             self.reverse_recipe = self.forward_recipe.get_reverse()
+            self.reverse_recipe.electron_label = self.electron_label
         else:
             self.reverse = local_context.get('reverse', None)
             self.reversible = True if local_context.get('reversible', None) is None else local_context.get('reversible', None)
@@ -704,9 +717,11 @@ def load(self, path, local_context=None, global_context=None, depository_labels=
                 self.reverse_template = Reaction(reactants=self.forward_template.products,
                                                  products=self.forward_template.reactants)
                 self.reverse_recipe = self.forward_recipe.get_reverse()
+                self.reverse_recipe.electron_label = self.electron_label
                 if self.reverse is None:
                     self.reverse = '{0}_reverse'.format(self.label)
 
+        self.forward_recipe.electron_label = self.electron_label
         self.rules = KineticsRules(label='{0}/rules'.format(self.label))
         logging.debug("Loading kinetics family rules from {0}".format(os.path.join(path, 'rules.py')))
         self.rules.load(os.path.join(path, 'rules.py'), local_context, global_context)
@@ -971,6 +986,8 @@ def save_groups(self, path):
             f.write('reactantNum = {0}\n\n'.format(self.reactant_num))
         if self.product_num is not None:
             f.write('productNum = {0}\n\n'.format(self.product_num))
+        if self.electron_label is not None:
+            f.write('electronLabel = "{0}"\n\n'.format(self.electron_label))
 
         if self.auto_generated is not None:
             f.write('autoGenerated = {0}\n\n'.format(self.auto_generated))
@@ -1309,6 +1326,7 @@ def apply_recipe(self, reactant_structures, forward=True, unique=True):
         # There is some hardcoding of reaction families in this function, so
         # we need the label of the reaction family for this
         label = self.label.lower()
+        remove_electrons = True
 
         # Merge reactant structures into single structure
         # Also copy structures so we don't modify the originals
@@ -1381,8 +1399,14 @@ def apply_recipe(self, reactant_structures, forward=True, unique=True):
 
             # Generate the product structure by applying the recipe
             self.forward_recipe.apply_forward(reactant_structure, unique)
+            actions = [a[0] for a in self.forward_recipe.actions]
+            if "GAIN_CHARGE" in actions:
+                remove_electrons = False
         else:
             self.reverse_recipe.apply_forward(reactant_structure, unique)
+            actions = [a[0] for a in self.reverse_recipe.actions]
+            if "GAIN_CHARGE" in actions:
+                remove_electrons = False
 
         # Now that we have applied the recipe, let's start calling
         # this thing the product_structure (although it's the same object in memory)
@@ -1490,7 +1514,10 @@ def apply_recipe(self, reactant_structures, forward=True, unique=True):
             product_num = self.product_num or len(template.products)
 
         # Split product structure into multiple species if necessary
-        product_structures = [p for p in product_structure.split() if not p.is_electron()]
+        if remove_electrons:
+            product_structures = [p for p in product_structure.split() if not p.is_electron()]
+        else:
+            product_structures = [p for p in product_structure.split()]
 
         # Make sure we've made the expected number of products
         if product_num != len(product_structures):

rmgpy/data/kinetics/familyTest.py

@@ -43,6 +43,7 @@
 from rmgpy.data.thermo import ThermoDatabase
 from rmgpy.molecule import Molecule
 from rmgpy.species import Species
+from rmgpy.reaction import Reaction
 
 
 ###################################################
@@ -624,6 +625,52 @@ def test_r_addition_com(self):
         self.assertEqual(len(products), 1)
         self.assertTrue(expected_products[0].is_isomorphic(products[0]))
 
+    def test_surface_proton_electron_reduction_alpha(self):
+        """
+        Test that the Surface_Proton_Electron_Reduction_Alpha family can successfully match the reaction and returns properly product structures.
+        """
+        family = self.database.families['Surface_Proton_Electron_Reduction_Alpha']
+        m_electron = Molecule().from_smiles("e")
+        m_proton = Molecule().from_smiles("[H+]")
+        m_x = Molecule().from_adjacency_list("1 X u0 p0")
+        m_ch2x = Molecule().from_adjacency_list(
+            """
+            1 C u0 p0 c0 {2,S} {3,S} {4,D}
+            2 H u0 p0 c0 {1,S}
+            3 H u0 p0 c0 {1,S}
+            4 X u0 p0 c0 {1,D}
+            """
+            )
+        m_ch3x = Molecule().from_adjacency_list(
+            """
+            1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+            2 H u0 p0 c0 {1,S}
+            3 H u0 p0 c0 {1,S}
+            4 H u0 p0 c0 {1,S}
+            5 X u0 p0 c0 {1,S}
+            """
+            )
+
+        reactants = [m_electron,m_proton,m_ch2x]
+        expected_products = [m_ch3x]
+
+        labeled_rxn = Reaction(reactants=reactants, products=expected_products)
+        family.add_atom_labels_for_reaction(labeled_rxn)
+        prods = family.apply_recipe([m.molecule[0] for m in labeled_rxn.reactants])
+        self.assertTrue(expected_products[0].is_isomorphic(prods[0]))
+
+        self.assertEqual(len(prods), 1)
+        self.assertTrue(expected_products[0].is_isomorphic(prods[0]))
+        reacts = family.apply_recipe(prods, forward=False)
+        self.assertEqual(len(reacts), 3)
+
+        prods = [Species(molecule=[p]) for p in prods]
+        reacts = [Species(molecule=[r]) for r in reacts]
+
+        fam_rxn = Reaction(reactants=reacts,products=prods)
+
+        self.assertTrue(fam_rxn.is_isomorphic(labeled_rxn))
+
     def test_save_family(self):
         """
 

rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/groups.py

@@ -20,10 +20,10 @@
 
 reactantNum = 3
 productNum = 1
+electronLabel = '*4'
 
 recipe(actions=[
     ['LOSE_CHARGE', '*3', 1],
-    # ['LOSE_RADICAL','*3', 1],
     ['CHANGE_BOND', '*1', -1, '*2'],
     ['FORM_BOND', '*1', 1, '*3'],
 ])