Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Chemical identity comparison for Molecules and Species #1329

Merged
merged 20 commits into from May 22, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
3da3c0e
Refactor some class properties using @property decorator
mliu49 Jan 25, 2018
db1e492
Refactor Molecule SMILES and InChI attributes
mliu49 Jan 25, 2018
aeaf3b5
Add read-only inchi property to Species
mliu49 Jan 25, 2018
3594fc5
Add read-only fingerprint property to Species
mliu49 Jan 25, 2018
010513e
Add read-only multiplicity property to Species
mliu49 Feb 5, 2018
1b9884b
Add unit tests for new Species properties
mliu49 Jan 26, 2018
3d3473f
Add strict argument for isomorphism comparison
mliu49 Aug 27, 2018
da870ca
Add strict argument to Species.isIsomorphic
mliu49 Apr 3, 2019
e36c78e
Remove generate_res argument from Species.isIsomorphic
mliu49 Apr 4, 2019
4422e9c
Rename isomorphic_species_lists to same_species_lists
mliu49 Apr 3, 2019
f2a2296
Add strict argument to Reaction.isIsomorphic
mliu49 Apr 3, 2019
75da1f3
Refactor CERM.checkForExistingSpecies using strict=False isomorphism
mliu49 Apr 4, 2019
de012ad
Refactor product checking in __generateReactions
mliu49 Sep 13, 2018
9e8a491
Revise test for prod_resonance option for generating reactions
mliu49 Sep 13, 2018
7b4604e
Add strict option to Graph.isMappingValid
mliu49 Apr 4, 2019
9052975
Add strict argument to isIdentical methods
mliu49 Apr 4, 2019
78f8a2d
Do not generate resonance structures for degeneracy determination
mliu49 Apr 4, 2019
4084b1e
Make sure selected molecule is reactive
mliu49 Oct 30, 2018
27baa4e
Enable Species instantiation by SMILES or InChI argument
mliu49 Apr 4, 2019
b3ff5c5
Fix reaction degeneracy bug with keep_isomorphic argument
mliu49 Apr 4, 2019
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
9 changes: 4 additions & 5 deletions rmgpy/data/kinetics/common.py
Expand Up @@ -198,7 +198,7 @@ def ensure_independent_atom_ids(input_species, resonance=True):
Modifies the list in place (replacing :class:`Molecule` with :class:`Species`).
Returns None.
"""
ensure_species(input_species, resonance=resonance)
ensure_species(input_species, resonance=resonance, keep_isomorphic=True)
# Method to check that all species' atom ids are different
def independent_ids():
num_atoms = 0
Expand All @@ -214,7 +214,7 @@ def independent_ids():
logging.debug('identical atom ids found between species. regenerating')
for species in input_species:
unreactive_mol_list = [mol for mol in species.molecule if not mol.reactive]
mol = species.molecule[0]
mol = [mol for mol in species.molecule if mol.reactive][0] # Choose first reactive molecule
mol.assignAtomIDs()
species.molecule = [mol]
# Remake resonance structures with new labels
Expand Down Expand Up @@ -275,7 +275,6 @@ def find_degenerate_reactions(rxn_list, same_reactants=None, template=None, kine
# with degenerate transition states
sorted_rxns = []
for rxn0 in selected_rxns:
# find resonance structures for rxn0
rxn0.ensure_species()
if len(sorted_rxns) == 0:
# This is the first reaction, so create a new sublist
Expand All @@ -288,9 +287,9 @@ def find_degenerate_reactions(rxn_list, same_reactants=None, template=None, kine
identical = False
sameTemplate = True
for rxn in sub_list:
isomorphic = rxn0.isIsomorphic(rxn, checkIdentical=False, checkTemplateRxnProducts=True)
isomorphic = rxn0.isIsomorphic(rxn, checkIdentical=False, strict=False, checkTemplateRxnProducts=True)
if isomorphic:
identical = rxn0.isIsomorphic(rxn, checkIdentical=True, checkTemplateRxnProducts=True)
identical = rxn0.isIsomorphic(rxn, checkIdentical=True, strict=False, checkTemplateRxnProducts=True)
if identical:
# An exact copy of rxn0 is already in our list, so we can move on
break
Expand Down
10 changes: 5 additions & 5 deletions rmgpy/data/kinetics/database.py
Expand Up @@ -41,7 +41,7 @@
StickingCoefficientBEP, SurfaceArrhenius, SurfaceArrheniusBEP
from rmgpy.molecule import Molecule, Group
from rmgpy.species import Species
from rmgpy.reaction import Reaction, isomorphic_species_lists
from rmgpy.reaction import Reaction, same_species_lists
from rmgpy.data.base import LogicNode

from .family import KineticsFamily
Expand Down Expand Up @@ -620,11 +620,11 @@ def getForwardReactionForFamilyEntry(self, entry, family, thermoDatabase):
# Remove from that set any reactions that don't produce the desired reactants and products
forward = []; reverse = []
for rxn in generatedReactions:
if (isomorphic_species_lists(reaction.reactants, rxn.reactants)
and isomorphic_species_lists(reaction.products, rxn.products)):
if (same_species_lists(reaction.reactants, rxn.reactants)
and same_species_lists(reaction.products, rxn.products)):
forward.append(rxn)
if (isomorphic_species_lists(reaction.reactants, rxn.products)
and isomorphic_species_lists(reaction.products, rxn.reactants)):
if (same_species_lists(reaction.reactants, rxn.products)
and same_species_lists(reaction.products, rxn.reactants)):
reverse.append(rxn)

# We should now know whether the reaction is given in the forward or
Expand Down
29 changes: 10 additions & 19 deletions rmgpy/data/kinetics/family.py
Expand Up @@ -44,7 +44,7 @@
from rmgpy.constraints import failsSpeciesConstraints
from rmgpy.data.base import Database, Entry, LogicNode, LogicOr, ForbiddenStructures,\
getAllCombinations
from rmgpy.reaction import Reaction, isomorphic_species_lists
from rmgpy.reaction import Reaction, same_species_lists
from rmgpy import settings
from rmgpy.reaction import Reaction
from rmgpy.kinetics import Arrhenius, SurfaceArrhenius,\
Expand Down Expand Up @@ -1570,7 +1570,7 @@ def __createReaction(self, reactants, products, is_forward):
"""

# Make sure the products are in fact different than the reactants
if isomorphic_species_lists(reactants, products):
if same_species_lists(reactants, products):
return None

# Create and return template reaction object
Expand Down Expand Up @@ -1679,6 +1679,8 @@ def addReverseAttribute(self, rxn, react_non_reactive=True):
elif rxn.products[1].isIsomorphic(rxn.products[2]):
sameReactants = 2

ensure_independent_atom_ids(rxn.products)

reactionList = self.__generateReactions([spc.molecule for spc in rxn.products],
products=rxn.reactants, forward=True,
react_non_reactive=react_non_reactive)
Expand Down Expand Up @@ -1713,7 +1715,7 @@ def addReverseAttribute(self, rxn, react_non_reactive=True):
else:
logging.error("Still experiencing error: Expecting one matching reverse reaction, not {0} in reaction family {1} for forward reaction {2}.\n".format(len(reactions), self.label, str(rxn)))
raise KineticsError("Did not find reverse reaction in reaction family {0} for reaction {1}.".format(self.label, str(rxn)))
elif len(reactions) > 1 and not all([reactions[0].isIsomorphic(other, checkTemplateRxnProducts=True) for other in reactions]):
elif len(reactions) > 1 and not all([reactions[0].isIsomorphic(other, strict=False, checkTemplateRxnProducts=True) for other in reactions]):
logging.error("Expecting one matching reverse reaction. Recieved {0} reactions with multiple non-isomorphic ones in reaction family {1} for forward reaction {2}.\n".format(len(reactions), self.label, str(rxn)))
logging.info("Found the following reverse reactions")
for rxn0 in reactions:
Expand Down Expand Up @@ -2129,24 +2131,13 @@ def generate_products_and_reactions(order):
# If products is given, remove reactions from the reaction list that
# don't generate the given products
if products is not None:
ensure_species(products, resonance=prod_resonance)

rxnList0 = rxnList[:]
rxnList = []
for reaction in rxnList0:

products0 = reaction.products[:] if forward else reaction.reactants[:]

# For aromatics, generate aromatic resonance structures to accurately identify isomorphic species
if prod_resonance:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why can we get rid of this stuff?

for i, product in enumerate(products0):
if product.isCyclic:
aromaticStructs = generate_optimal_aromatic_resonance_structures(product)
if aromaticStructs:
products0[i] = aromaticStructs[0]

# Skip reactions that don't match the given products
if isomorphic_species_lists(products, products0):
products0 = reaction.products if forward else reaction.reactants
# Only keep reactions which give the requested products
# If prod_resonance=True, then use strict=False to consider all resonance structures
if same_species_lists(products, products0, strict=not prod_resonance):
rxnList.append(reaction)

# Determine the reactant-product pairs to use for flux analysis
Expand Down Expand Up @@ -2590,7 +2581,7 @@ def getLabeledReactantsAndProducts(self, reactants, products):
pass
else:
if product_structures is not None:
if isomorphic_species_lists(list(products), list(product_structures)):
if same_species_lists(list(products), list(product_structures)):
return reactant_structures, product_structures
else:
continue
Expand Down
18 changes: 13 additions & 5 deletions rmgpy/data/kinetics/kineticsTest.py
Expand Up @@ -590,6 +590,19 @@ def test_degeneracy_multiple_resonance_different_template(self):

self.assertFalse(reaction_list[0].duplicate)

def test_degeneracy_resonance_keep_isomorphic(self):
"""Test that we get the correct degeneracy for [CH2]C=C[CH2] + [H].

Incorrect results would be obtained if isomorphic resonance structures are not kept."""
family_label = 'R_Recombination'
reactants = ['[CH2]C=C[CH2]', '[OH]']
products = ['[CH2]C(O)C=C']

correct_rxn_num = 1
correct_degeneracy = {2}

self.assert_correct_reaction_degeneracy(reactants, correct_rxn_num, correct_degeneracy, family_label, products)


class TestKineticsCommentsParsing(unittest.TestCase):

Expand Down Expand Up @@ -915,8 +928,6 @@ def test_generate_reactions_from_families_product_resonance(self):
self.assertEqual(len(reaction_list), 1)
self.assertEqual(reaction_list[0].degeneracy, 2)



def test_generate_reactions_from_families_product_resonance2(self):
"""Test that we can specify the no product resonance structure when generating reactions"""
reactants = [
Expand All @@ -931,9 +942,6 @@ def test_generate_reactions_from_families_product_resonance2(self):
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 = [
Expand Down
8 changes: 4 additions & 4 deletions rmgpy/molecule/graph.pxd
Expand Up @@ -40,7 +40,7 @@ cdef class Vertex(object):

cpdef Vertex copy(self)

cpdef bint equivalent(self, Vertex other) except -2
cpdef bint equivalent(self, Vertex other, bint strict=?) except -2

cpdef bint isSpecificCaseOf(self, Vertex other) except -2

Expand Down Expand Up @@ -110,9 +110,9 @@ cdef class Graph(object):

cpdef restore_vertex_order(self)

cpdef bint isIsomorphic(self, Graph other, dict initialMap=?, bint saveOrder=?) except -2
cpdef bint isIsomorphic(self, Graph other, dict initialMap=?, bint saveOrder=?, bint strict=?) except -2

cpdef list findIsomorphism(self, Graph other, dict initialMap=?, bint saveOrder=?)
cpdef list findIsomorphism(self, Graph other, dict initialMap=?, bint saveOrder=?, bint strict=?)

cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?, bint saveOrder=?) except -2

Expand Down Expand Up @@ -156,6 +156,6 @@ cdef class Graph(object):

cpdef list getLargestRing(self, Vertex vertex)

cpdef bint isMappingValid(self, Graph other, dict mapping, bint equivalent=?) except -2
cpdef bint isMappingValid(self, Graph other, dict mapping, bint equivalent=?, bint strict=?) except -2

cpdef list get_edges_in_cycle(self, list vertices, bint sort=?)
54 changes: 36 additions & 18 deletions rmgpy/molecule/graph.pyx
Expand Up @@ -95,7 +95,7 @@ cdef class Vertex(object):
new = Vertex()
return new

cpdef bint equivalent(self, Vertex other) except -2:
cpdef bint equivalent(self, Vertex other, bint strict=True) except -2:
"""
Return :data:`True` if two vertices `self` and `other` are semantically
equivalent, or :data:`False` if not. You should reimplement this
Expand Down Expand Up @@ -495,20 +495,30 @@ cdef class Graph(object):
else:
self.vertices = self.ordered_vertices

cpdef bint isIsomorphic(self, Graph other, dict initialMap=None, bint saveOrder=False) except -2:
cpdef bint isIsomorphic(self, Graph other, dict initialMap=None, bint saveOrder=False, bint strict=True) except -2:
"""
Returns :data:`True` if two graphs are isomorphic and :data:`False`
otherwise. Uses the VF2 algorithm of Vento and Foggia.

Args:
initialMap (dict, optional): initial atom mapping to use
saveOrder (bool, optional): if ``True``, reset atom order after performing atom isomorphism
strict (bool, optional): if ``False``, perform isomorphism ignoring electrons
"""
return vf2.isIsomorphic(self, other, initialMap, saveOrder=saveOrder)
return vf2.isIsomorphic(self, other, initialMap, saveOrder=saveOrder, strict=strict)

cpdef list findIsomorphism(self, Graph other, dict initialMap=None, bint saveOrder=False):
cpdef list findIsomorphism(self, Graph other, dict initialMap=None, bint saveOrder=False, bint strict=True):
"""
Returns :data:`True` if `other` is subgraph isomorphic and :data:`False`
otherwise, and the matching mapping.
Uses the VF2 algorithm of Vento and Foggia.

Args:
initialMap (dict, optional): initial atom mapping to use
saveOrder (bool, optional): if ``True``, reset atom order after performing atom isomorphism
strict (bool, optional): if ``False``, perform isomorphism ignoring electrons
"""
return vf2.findIsomorphism(self, other, initialMap, saveOrder=saveOrder)
return vf2.findIsomorphism(self, other, initialMap, saveOrder=saveOrder, strict=strict)

cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=None, bint saveOrder=False) except -2:
"""
Expand Down Expand Up @@ -1065,25 +1075,28 @@ cdef class Graph(object):
return longest_cycle


cpdef bint isMappingValid(self, Graph other, dict mapping, bint equivalent=True) except -2:
cpdef bint isMappingValid(self, Graph other, dict mapping, bint equivalent=True, bint strict=True) except -2:
"""
Check that a proposed `mapping` of vertices from `self` to `other`
is valid by checking that the vertices and edges involved in the
mapping are mutually equivalent. If equivalent is true it checks
if atoms and edges are equivalent, if false it checks if they
are specific cases of each other.
mapping are mutually equivalent. If equivalent is ``True`` it checks
if atoms and edges are equivalent, if ``False`` it checks if they
are specific cases of each other. If strict is ``True``, electrons
and bond orders are considered, and ignored if ``False``.
"""
cdef Vertex vertex1, vertex2
cdef list vertices1, vertices2
cdef bint selfHasEdge, otherHasEdge
cdef int i, j

method = 'equivalent' if equivalent else 'isSpecificCaseOf'


# Check that the mapped pairs of vertices compare True
for vertex1, vertex2 in mapping.items():
if not getattr(vertex1,method)(vertex2):
return False
if equivalent:
if not vertex1.equivalent(vertex2, strict=strict):
return False
else:
if not vertex1.isSpecificCaseOf(vertex2):
return False

# Check that any edges connected mapped vertices are equivalent
vertices1 = mapping.keys()
Expand All @@ -1094,10 +1107,15 @@ cdef class Graph(object):
otherHasEdge = other.hasEdge(vertices2[i], vertices2[j])
if selfHasEdge and otherHasEdge:
# Both graphs have the edge, so we must check it for equivalence
edge1 = self.getEdge(vertices1[i], vertices1[j])
edge2 = other.getEdge(vertices2[i], vertices2[j])
if not getattr(edge1,method)(edge2):
return False
if strict:
edge1 = self.getEdge(vertices1[i], vertices1[j])
edge2 = other.getEdge(vertices2[i], vertices2[j])
if equivalent:
if not edge1.equivalent(edge2):
return False
else:
if not edge1.isSpecificCaseOf(edge2):
return False
elif not equivalent and selfHasEdge and not otherHasEdge:
#in the subgraph case self can have edges other doesn't have
continue
Expand Down
6 changes: 3 additions & 3 deletions rmgpy/molecule/group.pxd
Expand Up @@ -63,7 +63,7 @@ cdef class GroupAtom(Vertex):

cpdef applyAction(self, list action)

cpdef bint equivalent(self, Vertex other) except -2
cpdef bint equivalent(self, Vertex other, bint strict=?) except -2

cpdef bint isSpecificCaseOf(self, Vertex other) except -2

Expand Down Expand Up @@ -168,9 +168,9 @@ cdef class Group(Graph):

cpdef update_charge(self)

cpdef bint isIsomorphic(self, Graph other, dict initialMap=?, bint saveOrder=?) except -2
cpdef bint isIsomorphic(self, Graph other, dict initialMap=?, bint saveOrder=?, bint strict=?) except -2

cpdef list findIsomorphism(self, Graph other, dict initialMap=?, bint saveOrder=?)
cpdef list findIsomorphism(self, Graph other, dict initialMap=?, bint saveOrder=?, bint strict=?)

cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?, bint generateInitialMap=?, bint saveOrder=?) except -2

Expand Down