From 235bb036ce0ca23fac4ce33074f5ddb6a40bba12 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Sat, 23 Jan 2021 18:26:24 -0500 Subject: [PATCH] added electronLabel for rxn family and recipe --- rmgpy/data/kinetics/family.py | 35 ++++++++++++-- rmgpy/data/kinetics/familyTest.py | 47 +++++++++++++++++++ .../groups.py | 2 +- 3 files changed, 79 insertions(+), 5 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 3aa7582941..8f49c4d8db 100644 --- a/rmgpy/data/kinetics/family.py +++ b/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): diff --git a/rmgpy/data/kinetics/familyTest.py b/rmgpy/data/kinetics/familyTest.py index a5b34cecb7..497176538c 100644 --- a/rmgpy/data/kinetics/familyTest.py +++ b/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): """ diff --git a/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/groups.py b/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/groups.py index d7170d3eb7..82047ffe80 100644 --- a/rmgpy/test_data/testing_database/kinetics/families/Surface_Proton_Electron_Reduction_Alpha/groups.py +++ b/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'], ])