From a31603f1ac9c07d88c7691f1fa58ac78c79899a4 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 09:33:26 -0400 Subject: [PATCH 01/23] Created rmgpy.molecule subpackage. There is enough functionality involving molecules to warrant a full subpackage. Within the subpackage, feel free to split the content into several smaller modules to enable faster partial rebuilds. --- rmgpy/molecule/__init__.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) create mode 100644 rmgpy/molecule/__init__.py diff --git a/rmgpy/molecule/__init__.py b/rmgpy/molecule/__init__.py new file mode 100644 index 0000000000..9a0b0fbe19 --- /dev/null +++ b/rmgpy/molecule/__init__.py @@ -0,0 +1,29 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +################################################################################ +# +# RMG - Reaction Mechanism Generator +# +# Copyright (c) 2002-2009 Prof. William H. Green (whgreen@mit.edu) and the +# RMG Team (rmg_dev@mit.edu) +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the 'Software'), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +# +################################################################################ From 4c8db82244fccadba0c8ffca2d1962a06d391426 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 09:39:38 -0400 Subject: [PATCH 02/23] Moved several modules into the rmgpy.molecule subpackage. As you can see, we already have several modules devoted to working with molecules and molecular graphs. This will make them easier to find and generally clean up the source tree. --- rmgpy/{ => molecule}/atomtype.pxd | 0 rmgpy/{ => molecule}/atomtype.py | 0 rmgpy/{ => molecule}/element.pxd | 0 rmgpy/{ => molecule}/element.py | 0 rmgpy/{ => molecule}/graph.pxd | 0 rmgpy/{ => molecule}/graph.py | 0 rmgpy/{ => molecule}/group.pxd | 0 rmgpy/{ => molecule}/group.py | 0 rmgpy/{ => molecule}/molecule.pxd | 0 rmgpy/{ => molecule}/molecule.py | 0 rmgpy/{ => molecule}/molecule_draw.py | 0 setup.py | 12 +++++++----- 12 files changed, 7 insertions(+), 5 deletions(-) rename rmgpy/{ => molecule}/atomtype.pxd (100%) rename rmgpy/{ => molecule}/atomtype.py (100%) rename rmgpy/{ => molecule}/element.pxd (100%) rename rmgpy/{ => molecule}/element.py (100%) rename rmgpy/{ => molecule}/graph.pxd (100%) rename rmgpy/{ => molecule}/graph.py (100%) rename rmgpy/{ => molecule}/group.pxd (100%) rename rmgpy/{ => molecule}/group.py (100%) rename rmgpy/{ => molecule}/molecule.pxd (100%) rename rmgpy/{ => molecule}/molecule.py (100%) rename rmgpy/{ => molecule}/molecule_draw.py (100%) diff --git a/rmgpy/atomtype.pxd b/rmgpy/molecule/atomtype.pxd similarity index 100% rename from rmgpy/atomtype.pxd rename to rmgpy/molecule/atomtype.pxd diff --git a/rmgpy/atomtype.py b/rmgpy/molecule/atomtype.py similarity index 100% rename from rmgpy/atomtype.py rename to rmgpy/molecule/atomtype.py diff --git a/rmgpy/element.pxd b/rmgpy/molecule/element.pxd similarity index 100% rename from rmgpy/element.pxd rename to rmgpy/molecule/element.pxd diff --git a/rmgpy/element.py b/rmgpy/molecule/element.py similarity index 100% rename from rmgpy/element.py rename to rmgpy/molecule/element.py diff --git a/rmgpy/graph.pxd b/rmgpy/molecule/graph.pxd similarity index 100% rename from rmgpy/graph.pxd rename to rmgpy/molecule/graph.pxd diff --git a/rmgpy/graph.py b/rmgpy/molecule/graph.py similarity index 100% rename from rmgpy/graph.py rename to rmgpy/molecule/graph.py diff --git a/rmgpy/group.pxd b/rmgpy/molecule/group.pxd similarity index 100% rename from rmgpy/group.pxd rename to rmgpy/molecule/group.pxd diff --git a/rmgpy/group.py b/rmgpy/molecule/group.py similarity index 100% rename from rmgpy/group.py rename to rmgpy/molecule/group.py diff --git a/rmgpy/molecule.pxd b/rmgpy/molecule/molecule.pxd similarity index 100% rename from rmgpy/molecule.pxd rename to rmgpy/molecule/molecule.pxd diff --git a/rmgpy/molecule.py b/rmgpy/molecule/molecule.py similarity index 100% rename from rmgpy/molecule.py rename to rmgpy/molecule/molecule.py diff --git a/rmgpy/molecule_draw.py b/rmgpy/molecule/molecule_draw.py similarity index 100% rename from rmgpy/molecule_draw.py rename to rmgpy/molecule/molecule_draw.py diff --git a/setup.py b/setup.py index e53486a2c3..0a0f1c9122 100644 --- a/setup.py +++ b/setup.py @@ -57,12 +57,14 @@ def getMainExtensionModules(): return [ - Extension('rmgpy.atomtype', ['rmgpy/atomtype.py'], include_dirs=['.']), - Extension('rmgpy.element', ['rmgpy/element.py'], include_dirs=['.']), - Extension('rmgpy.graph', ['rmgpy/graph.py'], include_dirs=['.']), - Extension('rmgpy.group', ['rmgpy/group.py'], include_dirs=['.']), + # Molecules and molecular representations + Extension('rmgpy.molecule.atomtype', ['rmgpy/molecule/atomtype.py'], include_dirs=['.']), + Extension('rmgpy.molecule.element', ['rmgpy/molecule/element.py'], include_dirs=['.']), + Extension('rmgpy.molecule.graph', ['rmgpy/molecule/graph.py'], include_dirs=['.']), + Extension('rmgpy.molecule.group', ['rmgpy/molecule/group.py'], include_dirs=['.']), + Extension('rmgpy.molecule.molecule', ['rmgpy/molecule/molecule.py'], include_dirs=['.']), + # Miscellaneous Extension('rmgpy.kinetics', ['rmgpy/kinetics.py'], include_dirs=['.']), - Extension('rmgpy.molecule', ['rmgpy/molecule.py'], include_dirs=['.']), Extension('rmgpy.quantity', ['rmgpy/quantity.py'], include_dirs=['.']), Extension('rmgpy.reaction', ['rmgpy/reaction.py'], include_dirs=['.']), Extension('rmgpy.species', ['rmgpy/species.py'], include_dirs=['.']), From ac12f010f6aa67cc1709da3da1bb2109c20a3274 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 10:00:32 -0400 Subject: [PATCH 03/23] Updated references to moved rmgpy.molecule submodules. Most of the time I think you can just import things from the subpackage directly, e.g. :: from rmgpy.molecule import Molecule, Group, fromAdjacencyList rather than trying to remember which submolecule these are defined in. The only exception would be if you are cimporting from a Cython module; in this case you will need to refer to the full subpackage. --- rmgpy/data/base.py | 3 +-- rmgpy/data/kinetics.py | 3 +-- rmgpy/data/states.py | 2 +- rmgpy/data/thermo.py | 3 +-- rmgpy/molecule/__init__.py | 6 ++++++ rmgpy/molecule/group.pxd | 2 +- rmgpy/molecule/group.py | 4 ++-- rmgpy/molecule/molecule.pxd | 8 ++++---- rmgpy/molecule/molecule.py | 6 +++--- rmgpy/molecule/molecule_draw.py | 3 ++- rmgpy/reaction.pxd | 4 ++-- rmgpy/reaction.py | 2 +- rmgpy/rmg/output.py | 2 +- 13 files changed, 26 insertions(+), 22 deletions(-) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 17d254cb28..9dbe853dc0 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -41,8 +41,7 @@ import codecs from rmgpy.quantity import * -from rmgpy.molecule import Molecule -from rmgpy.group import Group, InvalidAdjacencyListError +from rmgpy.molecule import Molecule, Group, InvalidAdjacencyListError from reference import * diff --git a/rmgpy/data/kinetics.py b/rmgpy/data/kinetics.py index 55ba68c5a6..ca4734d51a 100644 --- a/rmgpy/data/kinetics.py +++ b/rmgpy/data/kinetics.py @@ -40,8 +40,7 @@ from rmgpy.quantity import Quantity from rmgpy.reaction import Reaction, ReactionError from rmgpy.kinetics import * -from rmgpy.group import GroupBond, Group -from rmgpy.molecule import Bond +from rmgpy.molecule import Bond, GroupBond, Group from rmgpy.species import Species ################################################################################ diff --git a/rmgpy/data/states.py b/rmgpy/data/states.py index 250aaf991b..e4331dc81c 100644 --- a/rmgpy/data/states.py +++ b/rmgpy/data/states.py @@ -34,7 +34,7 @@ from rmgpy.quantity import constants from rmgpy.statmech import * -from rmgpy.group import Group +from rmgpy.molecule import Group from base import * diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 39d1af5ad3..affcddfc8d 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -43,8 +43,7 @@ from rmgpy.quantity import constants from rmgpy.thermo import * -from rmgpy.molecule import Molecule, Atom, Bond -from rmgpy.group import Group +from rmgpy.molecule import Molecule, Atom, Bond, Group ################################################################################ diff --git a/rmgpy/molecule/__init__.py b/rmgpy/molecule/__init__.py index 9a0b0fbe19..10474c6ff5 100644 --- a/rmgpy/molecule/__init__.py +++ b/rmgpy/molecule/__init__.py @@ -27,3 +27,9 @@ # DEALINGS IN THE SOFTWARE. # ################################################################################ + +from .atomtype import * +from .element import * +from .molecule import * +from .group import * +from .molecule_draw import drawMolecule diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd index fe102a7398..34387aa48c 100644 --- a/rmgpy/molecule/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -24,7 +24,7 @@ # ################################################################################ -from graph cimport Vertex, Edge, Graph +from .graph cimport Vertex, Edge, Graph ################################################################################ diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index a4b4111212..7c388d00d1 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -35,8 +35,8 @@ import cython -from graph import Vertex, Edge, Graph -from atomtype import atomTypes +from .graph import Vertex, Edge, Graph +from .atomtype import atomTypes ################################################################################ diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 0b13537de7..ced7d6e799 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -24,10 +24,10 @@ # ################################################################################ -from graph cimport Vertex, Edge, Graph -from atomtype cimport AtomType -from group cimport GroupAtom, GroupBond, Group -from element cimport Element +from .graph cimport Vertex, Edge, Graph +from .atomtype cimport AtomType +from .group cimport GroupAtom, GroupBond, Group +from .element cimport Element ################################################################################ diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index b420300424..362349f2fb 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -39,9 +39,9 @@ import os import re import element as elements -from graph import Vertex, Edge, Graph -from group import GroupAtom, GroupBond, Group, ActionError, fromAdjacencyList, toAdjacencyList -from atomtype import AtomType, atomTypes, getAtomType +from .graph import Vertex, Edge, Graph +from .group import GroupAtom, GroupBond, Group, ActionError, fromAdjacencyList, toAdjacencyList +from .atomtype import AtomType, atomTypes, getAtomType ################################################################################ diff --git a/rmgpy/molecule/molecule_draw.py b/rmgpy/molecule/molecule_draw.py index e64bc170cf..7bb762b76a 100644 --- a/rmgpy/molecule/molecule_draw.py +++ b/rmgpy/molecule/molecule_draw.py @@ -92,7 +92,6 @@ import logging from numpy.linalg import LinAlgError -from rmgpy.molecule import * ################################################################################ @@ -1239,6 +1238,8 @@ def drawMolecule(molecule, path=None, surface=''): if __name__ == '__main__': + from rmgpy.molecule.molecule import Molecule + molecule = Molecule() # Test #1: Straight chain backbone, no functional groups diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 7bd21a8dd9..e8d746285b 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -26,8 +26,8 @@ from quantity cimport constants from species cimport Species, TransitionState -from molecule cimport Atom, Molecule -from element cimport Element +from rmgpy.molecule.molecule cimport Atom, Molecule +from rmgpy.molecule.element cimport Element from kinetics cimport KineticsModel, Arrhenius cimport numpy diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index aabbd1adba..3ccb140c4b 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -47,7 +47,7 @@ import os.path from quantity import constants -from molecule import Molecule +from rmgpy.molecule.molecule import Molecule from species import Species from kinetics import Arrhenius, KineticsData, ArrheniusEP, ThirdBody diff --git a/rmgpy/rmg/output.py b/rmgpy/rmg/output.py index 896137dd03..4483bdc5a0 100644 --- a/rmgpy/rmg/output.py +++ b/rmgpy/rmg/output.py @@ -61,7 +61,7 @@ def saveOutputHTML(path, reactionModel): from model import PDepReaction - from rmgpy.molecule_draw import drawMolecule + from rmgpy.molecule import drawMolecule try: import jinja2 except ImportError: From 5021e15d2a9d2f8d217ce5ccb3f4515012dff143 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 10:05:49 -0400 Subject: [PATCH 04/23] Moved rmgpy.molecule unit test modules into molecule subfolder. The idea is that the unit test modules should generally be organized parallel to the main rmgpy modules. --- unittest/molecule/__init__.py | 0 unittest/{ => molecule}/atomtypeTest.py | 0 unittest/{ => molecule}/elementTest.py | 0 unittest/{ => molecule}/graphTest.py | 0 unittest/{ => molecule}/groupTest.py | 0 unittest/{ => molecule}/moleculeTest.py | 0 6 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 unittest/molecule/__init__.py rename unittest/{ => molecule}/atomtypeTest.py (100%) rename unittest/{ => molecule}/elementTest.py (100%) rename unittest/{ => molecule}/graphTest.py (100%) rename unittest/{ => molecule}/groupTest.py (100%) rename unittest/{ => molecule}/moleculeTest.py (100%) diff --git a/unittest/molecule/__init__.py b/unittest/molecule/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/unittest/atomtypeTest.py b/unittest/molecule/atomtypeTest.py similarity index 100% rename from unittest/atomtypeTest.py rename to unittest/molecule/atomtypeTest.py diff --git a/unittest/elementTest.py b/unittest/molecule/elementTest.py similarity index 100% rename from unittest/elementTest.py rename to unittest/molecule/elementTest.py diff --git a/unittest/graphTest.py b/unittest/molecule/graphTest.py similarity index 100% rename from unittest/graphTest.py rename to unittest/molecule/graphTest.py diff --git a/unittest/groupTest.py b/unittest/molecule/groupTest.py similarity index 100% rename from unittest/groupTest.py rename to unittest/molecule/groupTest.py diff --git a/unittest/moleculeTest.py b/unittest/molecule/moleculeTest.py similarity index 100% rename from unittest/moleculeTest.py rename to unittest/molecule/moleculeTest.py From 3f8ecee9ce05c55329bcabddf4bec758ed0cbfed Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 10:11:03 -0400 Subject: [PATCH 05/23] Updated references to moved rmgpy.molecule submodules in unit tests. All of the unit tests in the unittest/molecule subfolder now run as they did at the start of this branch. --- unittest/molecule/atomtypeTest.py | 6 +++--- unittest/molecule/elementTest.py | 10 +++++----- unittest/molecule/graphTest.py | 2 +- unittest/molecule/groupTest.py | 4 ++-- unittest/molecule/moleculeTest.py | 6 +++--- unittest/test.py | 11 ++++++----- 6 files changed, 20 insertions(+), 19 deletions(-) diff --git a/unittest/molecule/atomtypeTest.py b/unittest/molecule/atomtypeTest.py index 9a30cacd1b..b5b6a7046d 100644 --- a/unittest/molecule/atomtypeTest.py +++ b/unittest/molecule/atomtypeTest.py @@ -7,8 +7,8 @@ import unittest -import rmgpy.atomtype -from rmgpy.atomtype import AtomType, getAtomType +import rmgpy.molecule.atomtype +from rmgpy.molecule.atomtype import AtomType, getAtomType ################################################################################ @@ -21,7 +21,7 @@ def setUp(self): """ A function run before each unit test in this class. """ - self.atomType = rmgpy.atomtype.atomTypes['Cd'] + self.atomType = rmgpy.molecule.atomtype.atomTypes['Cd'] def testPickle(self): """ diff --git a/unittest/molecule/elementTest.py b/unittest/molecule/elementTest.py index 3cbc4cb9bb..dcc25cf513 100644 --- a/unittest/molecule/elementTest.py +++ b/unittest/molecule/elementTest.py @@ -7,8 +7,8 @@ import unittest -from rmgpy.element import Element -import rmgpy.element +from rmgpy.molecule.element import Element +import rmgpy.molecule.element ################################################################################ @@ -21,7 +21,7 @@ def setUp(self): """ A function run before each unit test in this class. """ - self.element = rmgpy.element.C + self.element = rmgpy.molecule.element.C def testPickle(self): """ @@ -50,8 +50,8 @@ def testGetElement(self): """ Test the rmgpy.elements.getElement() method. """ - self.assertTrue(rmgpy.element.getElement(6) is self.element) - self.assertTrue(rmgpy.element.getElement('C') is self.element) + self.assertTrue(rmgpy.molecule.element.getElement(6) is self.element) + self.assertTrue(rmgpy.molecule.element.getElement('C') is self.element) ################################################################################ diff --git a/unittest/molecule/graphTest.py b/unittest/molecule/graphTest.py index cf2b3ae7d2..40817f7b6d 100644 --- a/unittest/molecule/graphTest.py +++ b/unittest/molecule/graphTest.py @@ -3,7 +3,7 @@ import unittest -from rmgpy.graph import * +from rmgpy.molecule.graph import * ################################################################################ diff --git a/unittest/molecule/groupTest.py b/unittest/molecule/groupTest.py index 8f8bbd263c..d5140c86cc 100644 --- a/unittest/molecule/groupTest.py +++ b/unittest/molecule/groupTest.py @@ -3,8 +3,8 @@ import unittest -from rmgpy.group import * -from rmgpy.atomtype import atomTypes +from rmgpy.molecule.group import * +from rmgpy.molecule.atomtype import atomTypes ################################################################################ diff --git a/unittest/molecule/moleculeTest.py b/unittest/molecule/moleculeTest.py index 3eb263a272..a6cf957643 100644 --- a/unittest/molecule/moleculeTest.py +++ b/unittest/molecule/moleculeTest.py @@ -3,9 +3,9 @@ import unittest -from rmgpy.molecule import * -from rmgpy.group import Group -from rmgpy.element import getElement, elementList +from rmgpy.molecule.molecule import * +from rmgpy.molecule.group import Group +from rmgpy.molecule.element import getElement, elementList ################################################################################ diff --git a/unittest/test.py b/unittest/test.py index 57ef7ee591..e47c83b011 100644 --- a/unittest/test.py +++ b/unittest/test.py @@ -10,12 +10,13 @@ from cantherm.gaussianTest import * from cantherm.geometryTest import * -from atomtypeTest import * -from elementTest import * -from graphTest import * -from groupTest import * +from molecule.atomtypeTest import * +from molecule.elementTest import * +from molecule.graphTest import * +from molecule.groupTest import * +from molecule.moleculeTest import * + from kineticsTest import * -from moleculeTest import * from quantityTest import * from reactionTest import * from speciesTest import * From 7f8e1d951eacf4a6846bce681b2be3cf34c19ae2 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 10:27:05 -0400 Subject: [PATCH 06/23] Created new rmgpy.molecule.symmetry module for symmetry number calculation. This functionality is complex enough to warrant its own module, especially when we get around to adding QMTP functionality to RMG-Py. --- rmgpy/molecule/molecule.pxd | 10 +- rmgpy/molecule/molecule.py | 378 +------------------------------- rmgpy/molecule/symmetry.pxd | 39 ++++ rmgpy/molecule/symmetry.py | 416 ++++++++++++++++++++++++++++++++++++ setup.py | 1 + 5 files changed, 462 insertions(+), 382 deletions(-) create mode 100644 rmgpy/molecule/symmetry.pxd create mode 100644 rmgpy/molecule/symmetry.py diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index ced7d6e799..4ecb0988b2 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -166,12 +166,4 @@ cdef class Molecule(Graph): cpdef findAllDelocalizationPaths(self, Atom atom1) - cpdef int calculateAtomSymmetryNumber(self, Atom atom) - - cpdef int calculateBondSymmetryNumber(self, Atom atom1, Atom atom2) - - cpdef int calculateAxisSymmetryNumber(self) - - cpdef int calculateCyclicSymmetryNumber(self) - - cpdef int calculateSymmetryNumber(self) + cpdef int calculateSymmetryNumber(self) except -1 diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 362349f2fb..ddf2593ee8 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -42,7 +42,7 @@ from .graph import Vertex, Edge, Graph from .group import GroupAtom, GroupBond, Group, ActionError, fromAdjacencyList, toAdjacencyList from .atomtype import AtomType, atomTypes, getAtomType - + ################################################################################ class Atom(Vertex): @@ -1208,383 +1208,15 @@ def countInternalRotors(self): count += 1 return count - def calculateAtomSymmetryNumber(self, atom): - """ - Return the symmetry number centered at `atom` in the structure. The - `atom` of interest must not be in a cycle. - """ - symmetryNumber = 1 - - single = 0; double = 0; triple = 0; benzene = 0 - numNeighbors = 0 - for bond in self.edges[atom].values(): - if bond.isSingle(): single += 1 - elif bond.isDouble(): double += 1 - elif bond.isTriple(): triple += 1 - elif bond.isBenzene(): benzene += 1 - numNeighbors += 1 - - # If atom has zero or one neighbors, the symmetry number is 1 - if numNeighbors < 2: return symmetryNumber - - # Create temporary structures for each functional group attached to atom - molecule = self.copy() - for atom2 in molecule.bonds[atom].keys(): molecule.removeBond(atom, atom2) - molecule.removeAtom(atom) - groups = molecule.split() - - # Determine equivalence of functional groups around atom - groupIsomorphism = dict([(group, dict()) for group in groups]) - for group1 in groups: - for group2 in groups: - if group1 is not group2 and group2 not in groupIsomorphism[group1]: - groupIsomorphism[group1][group2] = group1.isIsomorphic(group2) - groupIsomorphism[group2][group1] = groupIsomorphism[group1][group2] - elif group1 is group2: - groupIsomorphism[group1][group1] = True - count = [sum([int(groupIsomorphism[group1][group2]) for group2 in groups]) for group1 in groups] - for i in range(count.count(2) / 2): - count.remove(2) - for i in range(count.count(3) / 3): - count.remove(3); count.remove(3) - for i in range(count.count(4) / 4): - count.remove(4); count.remove(4); count.remove(4) - count.sort(); count.reverse() - - if atom.radicalElectrons == 0: - if single == 4: - # Four single bonds - if count == [4]: symmetryNumber *= 12 - elif count == [3, 1]: symmetryNumber *= 3 - elif count == [2, 2]: symmetryNumber *= 2 - elif count == [2, 1, 1]: symmetryNumber *= 1 - elif count == [1, 1, 1, 1]: symmetryNumber *= 1 - elif single == 2: - # Two single bonds - if count == [2]: symmetryNumber *= 2 - elif double == 2: - # Two double bonds - if count == [2]: symmetryNumber *= 2 - elif atom.radicalElectrons == 1: - if single == 3: - # Three single bonds - if count == [3]: symmetryNumber *= 6 - elif count == [2, 1]: symmetryNumber *= 2 - elif count == [1, 1, 1]: symmetryNumber *= 1 - elif atom.radicalElectrons == 2: - if single == 2: - # Two single bonds - if count == [2]: symmetryNumber *= 2 - - return symmetryNumber - - def calculateBondSymmetryNumber(self, atom1, atom2): - """ - Return the symmetry number centered at `bond` in the structure. - """ - bond = self.edges[atom1][atom2] - symmetryNumber = 1 - if bond.isSingle() or bond.isDouble() or bond.isTriple(): - if atom1.equivalent(atom2): - # An O-O bond is considered to be an "optical isomer" and so no - # symmetry correction will be applied - if atom1.atomType == atom2.atomType == 'Os' and \ - atom1.radicalElectrons == atom2.radicalElectrons == 0: - pass - # If the molecule is diatomic, then we don't have to check the - # ligands on the two atoms in this bond (since we know there - # aren't any) - elif len(self.vertices) == 2: - symmetryNumber = 2 - else: - molecule = self.copy() - molecule.removeBond(atom1, atom2) - fragments = molecule.split() - if len(fragments) != 2: return symmetryNumber - - fragment1, fragment2 = fragments - if atom1 in fragment1.atoms: fragment1.removeAtom(atom1) - if atom2 in fragment1.atoms: fragment1.removeAtom(atom2) - if atom1 in fragment2.atoms: fragment2.removeAtom(atom1) - if atom2 in fragment2.atoms: fragment2.removeAtom(atom2) - groups1 = fragment1.split() - groups2 = fragment2.split() - - # Test functional groups for symmetry - if len(groups1) == len(groups2) == 1: - if groups1[0].isIsomorphic(groups2[0]): symmetryNumber *= 2 - elif len(groups1) == len(groups2) == 2: - if groups1[0].isIsomorphic(groups2[0]) and groups1[1].isIsomorphic(groups2[1]): symmetryNumber *= 2 - elif groups1[1].isIsomorphic(groups2[0]) and groups1[0].isIsomorphic(groups2[1]): symmetryNumber *= 2 - elif len(groups1) == len(groups2) == 3: - if groups1[0].isIsomorphic(groups2[0]) and groups1[1].isIsomorphic(groups2[1]) and groups1[2].isIsomorphic(groups2[2]): symmetryNumber *= 2 - elif groups1[0].isIsomorphic(groups2[0]) and groups1[1].isIsomorphic(groups2[2]) and groups1[2].isIsomorphic(groups2[1]): symmetryNumber *= 2 - elif groups1[0].isIsomorphic(groups2[1]) and groups1[1].isIsomorphic(groups2[2]) and groups1[2].isIsomorphic(groups2[0]): symmetryNumber *= 2 - elif groups1[0].isIsomorphic(groups2[1]) and groups1[1].isIsomorphic(groups2[0]) and groups1[2].isIsomorphic(groups2[2]): symmetryNumber *= 2 - elif groups1[0].isIsomorphic(groups2[2]) and groups1[1].isIsomorphic(groups2[0]) and groups1[2].isIsomorphic(groups2[1]): symmetryNumber *= 2 - elif groups1[0].isIsomorphic(groups2[2]) and groups1[1].isIsomorphic(groups2[1]) and groups1[2].isIsomorphic(groups2[0]): symmetryNumber *= 2 - - return symmetryNumber - - def calculateAxisSymmetryNumber(self): - """ - Get the axis symmetry number correction. The "axis" refers to a series - of two or more cumulated double bonds (e.g. C=C=C, etc.). Corrections - for single C=C bonds are handled in getBondSymmetryNumber(). - - Each axis (C=C=C) has the potential to double the symmetry number. - If an end has 0 or 1 groups (eg. =C=CJJ or =C=C-R) then it cannot - alter the axis symmetry and is disregarded:: - - A=C=C=C.. A-C=C=C=C-A - - s=1 s=1 - - If an end has 2 groups that are different then it breaks the symmetry - and the symmetry for that axis is 1, no matter what's at the other end:: - - A\ A\ /A - T=C=C=C=C-A T=C=C=C=T - B/ A/ \B - s=1 s=1 - - If you have one or more ends with 2 groups, and neither end breaks the - symmetry, then you have an axis symmetry number of 2:: - - A\ /B A\ - C=C=C=C=C C=C=C=C-B - A/ \B A/ - s=2 s=2 - """ - - symmetryNumber = 1 - - # List all double bonds in the structure - doubleBonds = [] - for atom1 in self.edges: - for atom2 in self.edges[atom1]: - if self.edges[atom1][atom2].isDouble() and self.vertices.index(atom1) < self.vertices.index(atom2): - doubleBonds.append((atom1, atom2)) - - # Search for adjacent double bonds - cumulatedBonds = [] - for i, bond1 in enumerate(doubleBonds): - atom11, atom12 = bond1 - for bond2 in doubleBonds[i+1:]: - atom21, atom22 = bond2 - if atom11 is atom21 or atom11 is atom22 or atom12 is atom21 or atom12 is atom22: - listToAddTo = None - for cumBonds in cumulatedBonds: - if (atom11, atom12) in cumBonds or (atom21, atom22) in cumBonds: - listToAddTo = cumBonds - if listToAddTo is not None: - if (atom11, atom12) not in listToAddTo: listToAddTo.append((atom11, atom12)) - if (atom21, atom22) not in listToAddTo: listToAddTo.append((atom21, atom22)) - else: - cumulatedBonds.append([(atom11, atom12), (atom21, atom22)]) - - # Also keep isolated double bonds - for bond1 in doubleBonds: - for bonds in cumulatedBonds: - if bond1 in bonds: - break - else: - cumulatedBonds.append([bond1]) - - # For each set of adjacent double bonds, check for axis symmetry - for bonds in cumulatedBonds: - - # Do nothing if less than two cumulated bonds - if len(bonds) < 1: continue - - # Do nothing if axis is in cycle - found = False - for atom1, atom2 in bonds: - if self.isBondInCycle(atom1, atom2): found = True - if found: continue - - # Find terminal atoms in axis - # Terminal atoms labelled T: T=C=C=C=T - axis = [] - for bond in bonds: axis.extend(bond) - terminalAtoms = [] - for atom in axis: - if axis.count(atom) == 1: terminalAtoms.append(atom) - if len(terminalAtoms) != 2: continue - - # Remove axis from (copy of) structure - structure = self.copy() - for atom1, atom2 in bonds: - structure.removeBond(atom1, atom2) - atomsToRemove = [] - for atom in structure.atoms: - if len(structure.bonds[atom]) == 0 and atom not in terminalAtoms: # it's not bonded to anything - atomsToRemove.append(atom) - for atom in atomsToRemove: structure.removeAtom(atom) - - # Split remaining fragments of structure - end_fragments = structure.split() - - # - # there can be two groups at each end A\ /B - # T=C=C=C=T - # A/ \B - - # to start with nothing has broken symmetry about the axis - symmetry_broken=False - end_fragments_to_remove = [] - for fragment in end_fragments: # a fragment is one end of the axis - - # remove the atom that was at the end of the axis and split what's left into groups - terminalAtom = None - for atom in terminalAtoms: - if atom in fragment.atoms: - terminalAtom = atom - fragment.removeAtom(atom) - break - else: - continue - - groups = [] - if len(fragment.atoms) > 0: - groups = fragment.split() - - # If end has only one group then it can't contribute to (nor break) axial symmetry - # Eg. this has no axis symmetry: A-T=C=C=C=T-A - # so we remove this end from the list of interesting end fragments - if len(groups) == 0: - end_fragments_to_remove.append(fragment) - continue # next end fragment - elif len(groups)==1 and terminalAtom.radicalElectrons == 0: - end_fragments_to_remove.append(fragment) - continue # next end fragment - elif len(groups)==1 and terminalAtom.radicalElectrons != 0: - symmetry_broken = True - elif len(groups)==2: - if not groups[0].isIsomorphic(groups[1]): - # this end has broken the symmetry of the axis - symmetry_broken = True - - for fragment in end_fragments_to_remove: - end_fragments.remove(fragment) - - # If there are end fragments left that can contribute to symmetry, - # and none of them broke it, then double the symmetry number - # NB>> This assumes coordination number of 4 (eg. Carbon). - # And would be wrong if we had /B - # =C=C=C=C=T-B - # \B - # (for some T with coordination number 5). - if end_fragments and not symmetry_broken: - symmetryNumber *= 2 - - return symmetryNumber - - def calculateCyclicSymmetryNumber(self): - """ - Get the symmetry number correction for cyclic regions of a molecule. - For complicated fused rings the smallest set of smallest rings is used. - """ - - symmetryNumber = 1 - - # Get symmetry number for each ring in structure - rings = self.getSmallestSetOfSmallestRings() - for ring in rings: - - # Make copy of structure - structure = self.copy() - - # Remove bonds of ring from structure - for i, atom1 in enumerate(ring): - for atom2 in ring[i+1:]: - if structure.hasBond(atom1, atom2): - structure.removeBond(atom1, atom2) - - structures = structure.split() - groups = [] - for struct in structures: - for atom in ring: - if atom in struct.atoms(): struct.removeAtom(atom) - groups.append(struct.split()) - - # Find equivalent functional groups on ring - equivalentGroups = [] - for group in groups: - found = False - for eqGroup in equivalentGroups: - if not found: - if group.isIsomorphic(eqGroup[0]): - eqGroup.append(group) - found = True - if not found: - equivalentGroups.append([group]) - - # Find equivalent bonds on ring - equivalentBonds = [] - for i, atom1 in enumerate(ring): - for atom2 in ring[i+1:]: - if self.hasBond(atom1, atom2): - bond = self.getBond(atom1, atom2) - found = False - for eqBond in equivalentBonds: - if not found: - if bond.equivalent(eqBond[0]): - eqBond.append(group) - found = True - if not found: - equivalentBonds.append([bond]) - - # Find maximum number of equivalent groups and bonds - maxEquivalentGroups = 0 - for groups in equivalentGroups: - if len(groups) > maxEquivalentGroups: - maxEquivalentGroups = len(groups) - maxEquivalentBonds = 0 - for bonds in equivalentBonds: - if len(bonds) > maxEquivalentBonds: - maxEquivalentBonds = len(bonds) - - if maxEquivalentGroups == maxEquivalentBonds == len(ring): - symmetryNumber *= len(ring) - else: - symmetryNumber *= max(maxEquivalentGroups, maxEquivalentBonds) - - print len(ring), maxEquivalentGroups, maxEquivalentBonds, symmetryNumber - - - return symmetryNumber - def calculateSymmetryNumber(self): """ Return the symmetry number for the structure. The symmetry number includes both external and internal modes. """ - symmetryNumber = 1 - - implicitH = self.implicitHydrogens - self.makeHydrogensExplicit() - - for atom in self.vertices: - if not self.isAtomInCycle(atom): - symmetryNumber *= self.calculateAtomSymmetryNumber(atom) - - for atom1 in self.edges: - for atom2 in self.edges[atom1]: - if self.vertices.index(atom1) < self.vertices.index(atom2) and not self.isBondInCycle(atom1, atom2): - symmetryNumber *= self.calculateBondSymmetryNumber(atom1, atom2) - - symmetryNumber *= self.calculateAxisSymmetryNumber() - - #if self.isCyclic(): - # symmetryNumber *= self.calculateCyclicSymmetryNumber() - - self.symmetryNumber = symmetryNumber - - if implicitH: self.makeHydrogensImplicit() - - return symmetryNumber - + from rmgpy.molecule.symmetry import calculateSymmetryNumber + self.symmetryNumber = calculateSymmetryNumber(self) + return self.symmetryNumber + def generateResonanceIsomers(self): """ Generate and return all of the resonance isomers of this molecule. diff --git a/rmgpy/molecule/symmetry.pxd b/rmgpy/molecule/symmetry.pxd new file mode 100644 index 0000000000..642d97491e --- /dev/null +++ b/rmgpy/molecule/symmetry.pxd @@ -0,0 +1,39 @@ +################################################################################ +# +# RMG - Reaction Mechanism Generator +# +# Copyright (c) 2009-2011 by the RMG Team (rmg_dev@mit.edu) +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the 'Software'), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +# +################################################################################ + +from .molecule cimport Atom, Bond, Molecule + +################################################################################ + +cpdef int calculateAtomSymmetryNumber(Molecule molecule, Atom atom) except -1 + +cpdef int calculateBondSymmetryNumber(Molecule molecule, Atom atom1, Atom atom2) except -1 + +cpdef int calculateAxisSymmetryNumber(Molecule molecule) except -1 + +cpdef int calculateCyclicSymmetryNumber(Molecule molecule) except -1 + +cpdef int calculateSymmetryNumber(Molecule molecule) except -1 diff --git a/rmgpy/molecule/symmetry.py b/rmgpy/molecule/symmetry.py new file mode 100644 index 0000000000..aafcc8eea5 --- /dev/null +++ b/rmgpy/molecule/symmetry.py @@ -0,0 +1,416 @@ +#!/usr/bin/env python +# encoding: utf-8 + +################################################################################ +# +# RMG - Reaction Mechanism Generator +# +# Copyright (c) 2009-2011 by the RMG Team (rmg_dev@mit.edu) +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the 'Software'), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +# +################################################################################ + +""" +This module provides functionality for estimating the symmetry number of a +molecule from its chemical graph representation. +""" + +def calculateAtomSymmetryNumber(molecule, atom): + """ + Return the symmetry number centered at `atom` in the structure. The + `atom` of interest must not be in a cycle. + """ + symmetryNumber = 1 + + single = 0; double = 0; triple = 0; benzene = 0 + numNeighbors = 0 + for bond in molecule.edges[atom].values(): + if bond.isSingle(): single += 1 + elif bond.isDouble(): double += 1 + elif bond.isTriple(): triple += 1 + elif bond.isBenzene(): benzene += 1 + numNeighbors += 1 + + # If atom has zero or one neighbors, the symmetry number is 1 + if numNeighbors < 2: return symmetryNumber + + # Create temporary structures for each functional group attached to atom + molecule = molecule.copy() + for atom2 in molecule.bonds[atom].keys(): molecule.removeBond(atom, atom2) + molecule.removeAtom(atom) + groups = molecule.split() + + # Determine equivalence of functional groups around atom + groupIsomorphism = dict([(group, dict()) for group in groups]) + for group1 in groups: + for group2 in groups: + if group1 is not group2 and group2 not in groupIsomorphism[group1]: + groupIsomorphism[group1][group2] = group1.isIsomorphic(group2) + groupIsomorphism[group2][group1] = groupIsomorphism[group1][group2] + elif group1 is group2: + groupIsomorphism[group1][group1] = True + count = [sum([int(groupIsomorphism[group1][group2]) for group2 in groups]) for group1 in groups] + for i in range(count.count(2) / 2): + count.remove(2) + for i in range(count.count(3) / 3): + count.remove(3); count.remove(3) + for i in range(count.count(4) / 4): + count.remove(4); count.remove(4); count.remove(4) + count.sort(); count.reverse() + + if atom.radicalElectrons == 0: + if single == 4: + # Four single bonds + if count == [4]: symmetryNumber *= 12 + elif count == [3, 1]: symmetryNumber *= 3 + elif count == [2, 2]: symmetryNumber *= 2 + elif count == [2, 1, 1]: symmetryNumber *= 1 + elif count == [1, 1, 1, 1]: symmetryNumber *= 1 + elif single == 2: + # Two single bonds + if count == [2]: symmetryNumber *= 2 + elif double == 2: + # Two double bonds + if count == [2]: symmetryNumber *= 2 + elif atom.radicalElectrons == 1: + if single == 3: + # Three single bonds + if count == [3]: symmetryNumber *= 6 + elif count == [2, 1]: symmetryNumber *= 2 + elif count == [1, 1, 1]: symmetryNumber *= 1 + elif atom.radicalElectrons == 2: + if single == 2: + # Two single bonds + if count == [2]: symmetryNumber *= 2 + + return symmetryNumber + +################################################################################ + +def calculateBondSymmetryNumber(molecule, atom1, atom2): + """ + Return the symmetry number centered at `bond` in the structure. + """ + bond = molecule.edges[atom1][atom2] + symmetryNumber = 1 + if bond.isSingle() or bond.isDouble() or bond.isTriple(): + if atom1.equivalent(atom2): + # An O-O bond is considered to be an "optical isomer" and so no + # symmetry correction will be applied + if atom1.atomType == atom2.atomType == 'Os' and \ + atom1.radicalElectrons == atom2.radicalElectrons == 0: + pass + # If the molecule is diatomic, then we don't have to check the + # ligands on the two atoms in this bond (since we know there + # aren't any) + elif len(molecule.vertices) == 2: + symmetryNumber = 2 + else: + molecule = molecule.copy() + molecule.removeBond(atom1, atom2) + fragments = molecule.split() + if len(fragments) != 2: return symmetryNumber + + fragment1, fragment2 = fragments + if atom1 in fragment1.atoms: fragment1.removeAtom(atom1) + if atom2 in fragment1.atoms: fragment1.removeAtom(atom2) + if atom1 in fragment2.atoms: fragment2.removeAtom(atom1) + if atom2 in fragment2.atoms: fragment2.removeAtom(atom2) + groups1 = fragment1.split() + groups2 = fragment2.split() + + # Test functional groups for symmetry + if len(groups1) == len(groups2) == 1: + if groups1[0].isIsomorphic(groups2[0]): symmetryNumber *= 2 + elif len(groups1) == len(groups2) == 2: + if groups1[0].isIsomorphic(groups2[0]) and groups1[1].isIsomorphic(groups2[1]): symmetryNumber *= 2 + elif groups1[1].isIsomorphic(groups2[0]) and groups1[0].isIsomorphic(groups2[1]): symmetryNumber *= 2 + elif len(groups1) == len(groups2) == 3: + if groups1[0].isIsomorphic(groups2[0]) and groups1[1].isIsomorphic(groups2[1]) and groups1[2].isIsomorphic(groups2[2]): symmetryNumber *= 2 + elif groups1[0].isIsomorphic(groups2[0]) and groups1[1].isIsomorphic(groups2[2]) and groups1[2].isIsomorphic(groups2[1]): symmetryNumber *= 2 + elif groups1[0].isIsomorphic(groups2[1]) and groups1[1].isIsomorphic(groups2[2]) and groups1[2].isIsomorphic(groups2[0]): symmetryNumber *= 2 + elif groups1[0].isIsomorphic(groups2[1]) and groups1[1].isIsomorphic(groups2[0]) and groups1[2].isIsomorphic(groups2[2]): symmetryNumber *= 2 + elif groups1[0].isIsomorphic(groups2[2]) and groups1[1].isIsomorphic(groups2[0]) and groups1[2].isIsomorphic(groups2[1]): symmetryNumber *= 2 + elif groups1[0].isIsomorphic(groups2[2]) and groups1[1].isIsomorphic(groups2[1]) and groups1[2].isIsomorphic(groups2[0]): symmetryNumber *= 2 + + return symmetryNumber + +################################################################################ + +def calculateAxisSymmetryNumber(molecule): + """ + Get the axis symmetry number correction. The "axis" refers to a series + of two or more cumulated double bonds (e.g. C=C=C, etc.). Corrections + for single C=C bonds are handled in getBondSymmetryNumber(). + + Each axis (C=C=C) has the potential to double the symmetry number. + If an end has 0 or 1 groups (eg. =C=CJJ or =C=C-R) then it cannot + alter the axis symmetry and is disregarded:: + + A=C=C=C.. A-C=C=C=C-A + + s=1 s=1 + + If an end has 2 groups that are different then it breaks the symmetry + and the symmetry for that axis is 1, no matter what's at the other end:: + + A\ A\ /A + T=C=C=C=C-A T=C=C=C=T + B/ A/ \B + s=1 s=1 + + If you have one or more ends with 2 groups, and neither end breaks the + symmetry, then you have an axis symmetry number of 2:: + + A\ /B A\ + C=C=C=C=C C=C=C=C-B + A/ \B A/ + s=2 s=2 + """ + + symmetryNumber = 1 + + # List all double bonds in the structure + doubleBonds = [] + for atom1 in molecule.edges: + for atom2 in molecule.edges[atom1]: + if molecule.edges[atom1][atom2].isDouble() and molecule.vertices.index(atom1) < molecule.vertices.index(atom2): + doubleBonds.append((atom1, atom2)) + + # Search for adjacent double bonds + cumulatedBonds = [] + for i, bond1 in enumerate(doubleBonds): + atom11, atom12 = bond1 + for bond2 in doubleBonds[i+1:]: + atom21, atom22 = bond2 + if atom11 is atom21 or atom11 is atom22 or atom12 is atom21 or atom12 is atom22: + listToAddTo = None + for cumBonds in cumulatedBonds: + if (atom11, atom12) in cumBonds or (atom21, atom22) in cumBonds: + listToAddTo = cumBonds + if listToAddTo is not None: + if (atom11, atom12) not in listToAddTo: listToAddTo.append((atom11, atom12)) + if (atom21, atom22) not in listToAddTo: listToAddTo.append((atom21, atom22)) + else: + cumulatedBonds.append([(atom11, atom12), (atom21, atom22)]) + + # Also keep isolated double bonds + for bond1 in doubleBonds: + for bonds in cumulatedBonds: + if bond1 in bonds: + break + else: + cumulatedBonds.append([bond1]) + + # For each set of adjacent double bonds, check for axis symmetry + for bonds in cumulatedBonds: + + # Do nothing if less than two cumulated bonds + if len(bonds) < 1: continue + + # Do nothing if axis is in cycle + found = False + for atom1, atom2 in bonds: + if molecule.isBondInCycle(atom1, atom2): found = True + if found: continue + + # Find terminal atoms in axis + # Terminal atoms labelled T: T=C=C=C=T + axis = [] + for bond in bonds: axis.extend(bond) + terminalAtoms = [] + for atom in axis: + if axis.count(atom) == 1: terminalAtoms.append(atom) + if len(terminalAtoms) != 2: continue + + # Remove axis from (copy of) structure + structure = molecule.copy() + for atom1, atom2 in bonds: + structure.removeBond(atom1, atom2) + atomsToRemove = [] + for atom in structure.atoms: + if len(structure.bonds[atom]) == 0 and atom not in terminalAtoms: # it's not bonded to anything + atomsToRemove.append(atom) + for atom in atomsToRemove: structure.removeAtom(atom) + + # Split remaining fragments of structure + end_fragments = structure.split() + + # + # there can be two groups at each end A\ /B + # T=C=C=C=T + # A/ \B + + # to start with nothing has broken symmetry about the axis + symmetry_broken=False + end_fragments_to_remove = [] + for fragment in end_fragments: # a fragment is one end of the axis + + # remove the atom that was at the end of the axis and split what's left into groups + terminalAtom = None + for atom in terminalAtoms: + if atom in fragment.atoms: + terminalAtom = atom + fragment.removeAtom(atom) + break + else: + continue + + groups = [] + if len(fragment.atoms) > 0: + groups = fragment.split() + + # If end has only one group then it can't contribute to (nor break) axial symmetry + # Eg. this has no axis symmetry: A-T=C=C=C=T-A + # so we remove this end from the list of interesting end fragments + if len(groups) == 0: + end_fragments_to_remove.append(fragment) + continue # next end fragment + elif len(groups)==1 and terminalAtom.radicalElectrons == 0: + end_fragments_to_remove.append(fragment) + continue # next end fragment + elif len(groups)==1 and terminalAtom.radicalElectrons != 0: + symmetry_broken = True + elif len(groups)==2: + if not groups[0].isIsomorphic(groups[1]): + # this end has broken the symmetry of the axis + symmetry_broken = True + + for fragment in end_fragments_to_remove: + end_fragments.remove(fragment) + + # If there are end fragments left that can contribute to symmetry, + # and none of them broke it, then double the symmetry number + # NB>> This assumes coordination number of 4 (eg. Carbon). + # And would be wrong if we had /B + # =C=C=C=C=T-B + # \B + # (for some T with coordination number 5). + if end_fragments and not symmetry_broken: + symmetryNumber *= 2 + + return symmetryNumber + +################################################################################ + +def calculateCyclicSymmetryNumber(molecule): + """ + Get the symmetry number correction for cyclic regions of a molecule. + For complicated fused rings the smallest set of smallest rings is used. + """ + + symmetryNumber = 1 + + # Get symmetry number for each ring in structure + rings = molecule.getSmallestSetOfSmallestRings() + for ring in rings: + + # Make copy of structure + structure = molecule.copy() + + # Remove bonds of ring from structure + for i, atom1 in enumerate(ring): + for atom2 in ring[i+1:]: + if structure.hasBond(atom1, atom2): + structure.removeBond(atom1, atom2) + + structures = structure.split() + groups = [] + for struct in structures: + for atom in ring: + if atom in struct.atoms(): struct.removeAtom(atom) + groups.append(struct.split()) + + # Find equivalent functional groups on ring + equivalentGroups = [] + for group in groups: + found = False + for eqGroup in equivalentGroups: + if not found: + if group.isIsomorphic(eqGroup[0]): + eqGroup.append(group) + found = True + if not found: + equivalentGroups.append([group]) + + # Find equivalent bonds on ring + equivalentBonds = [] + for i, atom1 in enumerate(ring): + for atom2 in ring[i+1:]: + if molecule.hasBond(atom1, atom2): + bond = molecule.getBond(atom1, atom2) + found = False + for eqBond in equivalentBonds: + if not found: + if bond.equivalent(eqBond[0]): + eqBond.append(group) + found = True + if not found: + equivalentBonds.append([bond]) + + # Find maximum number of equivalent groups and bonds + maxEquivalentGroups = 0 + for groups in equivalentGroups: + if len(groups) > maxEquivalentGroups: + maxEquivalentGroups = len(groups) + maxEquivalentBonds = 0 + for bonds in equivalentBonds: + if len(bonds) > maxEquivalentBonds: + maxEquivalentBonds = len(bonds) + + if maxEquivalentGroups == maxEquivalentBonds == len(ring): + symmetryNumber *= len(ring) + else: + symmetryNumber *= max(maxEquivalentGroups, maxEquivalentBonds) + + print len(ring), maxEquivalentGroups, maxEquivalentBonds, symmetryNumber + + + return symmetryNumber + +################################################################################ + +def calculateSymmetryNumber(molecule): + """ + Return the symmetry number for the structure. The symmetry number + includes both external and internal modes. + """ + symmetryNumber = 1 + + implicitH = molecule.implicitHydrogens + molecule.makeHydrogensExplicit() + + for atom in molecule.vertices: + if not molecule.isAtomInCycle(atom): + symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) + + for atom1 in molecule.edges: + for atom2 in molecule.edges[atom1]: + if molecule.vertices.index(atom1) < molecule.vertices.index(atom2) and not molecule.isBondInCycle(atom1, atom2): + symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) + + symmetryNumber *= calculateAxisSymmetryNumber(molecule) + + #if molecule.isCyclic(): + # symmetryNumber *= calculateCyclicSymmetryNumber(molecule) + + if implicitH: molecule.makeHydrogensImplicit() + + return symmetryNumber diff --git a/setup.py b/setup.py index 0a0f1c9122..64d0e5ba19 100644 --- a/setup.py +++ b/setup.py @@ -63,6 +63,7 @@ def getMainExtensionModules(): Extension('rmgpy.molecule.graph', ['rmgpy/molecule/graph.py'], include_dirs=['.']), Extension('rmgpy.molecule.group', ['rmgpy/molecule/group.py'], include_dirs=['.']), Extension('rmgpy.molecule.molecule', ['rmgpy/molecule/molecule.py'], include_dirs=['.']), + Extension('rmgpy.molecule.symmetry', ['rmgpy/molecule/symmetry.py'], include_dirs=['.']), # Miscellaneous Extension('rmgpy.kinetics', ['rmgpy/kinetics.py'], include_dirs=['.']), Extension('rmgpy.quantity', ['rmgpy/quantity.py'], include_dirs=['.']), From e9656a3f359b2ace128a4ce795bbd20d6cce5dc5 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 10:42:07 -0400 Subject: [PATCH 07/23] Split symmetry number unit tests to separate module. Again, the idea is to parallel the structure of the modules in the rmgpy.* source tree. Note that, as before, a few of the unit tests fail because they represent tricky edge cases that the current symmetry number algorithms do not yet capture. --- unittest/molecule/moleculeTest.py | 311 --------------------------- unittest/molecule/symmetryTest.py | 344 ++++++++++++++++++++++++++++++ 2 files changed, 344 insertions(+), 311 deletions(-) create mode 100644 unittest/molecule/symmetryTest.py diff --git a/unittest/molecule/moleculeTest.py b/unittest/molecule/moleculeTest.py index a6cf957643..2f40425815 100644 --- a/unittest/molecule/moleculeTest.py +++ b/unittest/molecule/moleculeTest.py @@ -902,14 +902,7 @@ def testAugmentedInChIKey(self): """) self.assertEqual(mol.toAugmentedInChIKey(), 'VGGSQFUCUMXWEO-UHFFFAOYSAmult3') -################################################################################ -class TestMoleculeSymmetry(unittest.TestCase): - """ - Contains unit tests of the methods for computing symmetry numbers for a - given Molecule object. - """ - def testLinearMethane(self): """ Test the Molecule.isLinear() method. @@ -1006,311 +999,7 @@ def testCountInternalRotorsDimethylAcetylene(self): This is a "hard" test that currently fails. """ self.assertEqual(Molecule().fromSMILES('CC#CC').countInternalRotors(), 1) - - def testAtomSymmetryNumberMethane(self): - """ - Test the Molecule.calculateAtomSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('C') - symmetryNumber = 1 - for atom in molecule.atoms: - if not molecule.isAtomInCycle(atom): - symmetryNumber *= molecule.calculateAtomSymmetryNumber(atom) - self.assertEqual(symmetryNumber, 12) - - def testAtomSymmetryNumberMethyl(self): - """ - Test the Molecule.calculateAtomSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('[CH3]') - symmetryNumber = 1 - for atom in molecule.atoms: - if not molecule.isAtomInCycle(atom): - symmetryNumber *= molecule.calculateAtomSymmetryNumber(atom) - self.assertEqual(symmetryNumber, 6) - - def testAtomSymmetryNumberEthane(self): - """ - Test the Molecule.calculateAtomSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('CC') - symmetryNumber = 1 - for atom in molecule.atoms: - if not molecule.isAtomInCycle(atom): - symmetryNumber *= molecule.calculateAtomSymmetryNumber(atom) - self.assertEqual(symmetryNumber, 9) - - def testAtomSymmetryNumberPropane(self): - """ - Test the Molecule.calculateAtomSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('CCC') - symmetryNumber = 1 - for atom in molecule.atoms: - if not molecule.isAtomInCycle(atom): - symmetryNumber *= molecule.calculateAtomSymmetryNumber(atom) - self.assertEqual(symmetryNumber, 18) - - def testAtomSymmetryNumberIsobutane(self): - """ - Test the Molecule.calculateAtomSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('CC(C)C') - symmetryNumber = 1 - for atom in molecule.atoms: - if not molecule.isAtomInCycle(atom): - symmetryNumber *= molecule.calculateAtomSymmetryNumber(atom) - self.assertEqual(symmetryNumber, 81) - def testBondSymmetryNumberEthane(self): - """ - Test the Molecule.calculateBondSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('CC') - symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: - if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): - symmetryNumber *= molecule.calculateBondSymmetryNumber(atom1, atom2) - self.assertEqual(symmetryNumber, 2) - - def testBondSymmetryNumberPropane(self): - """ - Test the Molecule.calculateBondSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('CCC') - symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: - if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): - symmetryNumber *= molecule.calculateBondSymmetryNumber(atom1, atom2) - self.assertEqual(symmetryNumber, 1) - - def testBondSymmetryNumberButane(self): - """ - Test the Molecule.calculateBondSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('CCCC') - symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: - if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): - symmetryNumber *= molecule.calculateBondSymmetryNumber(atom1, atom2) - self.assertEqual(symmetryNumber, 2) - - def testBondSymmetryNumberEthylene(self): - """ - Test the Molecule.calculateBondSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('C=C') - symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: - if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): - symmetryNumber *= molecule.calculateBondSymmetryNumber(atom1, atom2) - self.assertEqual(symmetryNumber, 2) - - def testBondSymmetryNumberAcetylene(self): - """ - Test the Molecule.calculateBondSymmtryNumber() method. - """ - molecule = Molecule().fromSMILES('C#C') - symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: - if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): - symmetryNumber *= molecule.calculateBondSymmetryNumber(atom1, atom2) - self.assertEqual(symmetryNumber, 2) - - def testAxisSymmetryNumberEthylene(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C').calculateAxisSymmetryNumber(), 2) - - def testAxisSymmetryNumberPropadiene(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=C').calculateAxisSymmetryNumber(), 2) - - def testAxisSymmetryNumberButatriene(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=C=C').calculateAxisSymmetryNumber(), 2) - - def testAxisSymmetryNumberButatrienyl(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=C=[CH]').calculateAxisSymmetryNumber(), 1) - - def testAxisSymmetryNumberPropadienyl(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=[C]').calculateAxisSymmetryNumber(), 2) - - def testAxisSymmetryNumber12Butadienyl(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('CC=C=[C]').calculateAxisSymmetryNumber(), 1) - - def testAxisSymmetryNumber12Hexadienyl(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=CCCC').calculateAxisSymmetryNumber(), 1) - - def testAxisSymmetryNumber1(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('CC(C)=C=C(CC)CC').calculateAxisSymmetryNumber(), 2) - - def testAxisSymmetryNumber2(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=C(C(C(C(C=C=C)=C=C)=C=C)=C=C)').calculateAxisSymmetryNumber(), 2) - - def testAxisSymmetryNumber3(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=[C]C(C)(C)[C]=C=C').calculateAxisSymmetryNumber(), 4) - - def testAxisSymmetryNumber4(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=C=O').calculateAxisSymmetryNumber(), 2) - - def testAxisSymmetryNumber5(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('CC=C=C=O').calculateAxisSymmetryNumber(), 1) - - def testAxisSymmetryNumber5(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=C=N').calculateAxisSymmetryNumber(), 1) - - def testAxisSymmetryNumber5(self): - """ - Test the Molecule.calculateAxisSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=C=[N]').calculateAxisSymmetryNumber(), 2) - -# def testCyclicSymmetryNumber(self): -# -# # cyclohexane -# molecule = Molecule().fromInChI('InChI=1/C6H12/c1-2-4-6-5-3-1/h1-6H2') -# molecule.makeHydrogensExplicit() -# symmetryNumber = molecule.calculateCyclicSymmetryNumber() -# self.assertEqual(symmetryNumber, 12) - - def testTotalSymmetryNumberEthane(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('CC').calculateSymmetryNumber(), 18) - - def testTotalSymmetryNumber1(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C=[C]C(C)(C)[C]=C=C').calculateSymmetryNumber(), '???') - - def testTotalSymmetryNumber2(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C(=CC(c1ccccc1)C([CH]CCCCCC)C=Cc1ccccc1)[CH]CCCCCC').calculateSymmetryNumber(), 1) - - def testSymmetryNumberHydroxyl(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('[OH]').calculateSymmetryNumber(), 1) - - def testSymmetryNumberOxygen(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('O=O').calculateAxisSymmetryNumber(), 1) - self.assertEqual(Molecule().fromSMILES('O=O').calculateSymmetryNumber(), 2) - - def testSymmetryNumberDicarbon(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('[C]#[C]').calculateSymmetryNumber(), 2) - - def testSymmetryNumberHydrogen(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('[H][H]').calculateSymmetryNumber(), 2) - - def testSymmetryNumberAcetylene(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C#C').calculateSymmetryNumber(), 2) - - def testSymmetryNumberButadiyne(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C#CC#C').calculateSymmetryNumber(), 2) - - def testSymmetryNumberMethane(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C').calculateSymmetryNumber(), 12) - - def testSymmetryNumberFormaldehyde(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=O').calculateSymmetryNumber(), 2) - - def testSymmetryNumberMethyl(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('[CH3]').calculateSymmetryNumber(), 6) - - def testSymmetryNumberWater(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('O').calculateSymmetryNumber(), 2) - - def testSymmetryNumberEthylene(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=C').calculateSymmetryNumber(), 4) - - def testSymmetryNumberEthenyl(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C=[CH]').calculateSymmetryNumber(), 1) - - def testSymmetryNumberCyclic(self): - """ - Test the Molecule.calculateSymmtryNumber() method. - """ - self.assertEqual(Molecule().fromSMILES('C1=C=C=1').calculateSymmetryNumber(), '6?') - ################################################################################ if __name__ == '__main__': diff --git a/unittest/molecule/symmetryTest.py b/unittest/molecule/symmetryTest.py new file mode 100644 index 0000000000..a3969fad1b --- /dev/null +++ b/unittest/molecule/symmetryTest.py @@ -0,0 +1,344 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- + +import unittest + +from rmgpy.molecule.molecule import * +from rmgpy.molecule.symmetry import * + +################################################################################ + +class TestMoleculeSymmetry(unittest.TestCase): + """ + Contains unit tests of the methods for computing symmetry numbers for a + given Molecule object. + """ + + def testAtomSymmetryNumberMethane(self): + """ + Test the Molecule.calculateAtomSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C') + symmetryNumber = 1 + for atom in molecule.atoms: + if not molecule.isAtomInCycle(atom): + symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) + self.assertEqual(symmetryNumber, 12) + + def testAtomSymmetryNumberMethyl(self): + """ + Test the Molecule.calculateAtomSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('[CH3]') + symmetryNumber = 1 + for atom in molecule.atoms: + if not molecule.isAtomInCycle(atom): + symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) + self.assertEqual(symmetryNumber, 6) + + def testAtomSymmetryNumberEthane(self): + """ + Test the Molecule.calculateAtomSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CC') + symmetryNumber = 1 + for atom in molecule.atoms: + if not molecule.isAtomInCycle(atom): + symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) + self.assertEqual(symmetryNumber, 9) + + def testAtomSymmetryNumberPropane(self): + """ + Test the Molecule.calculateAtomSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CCC') + symmetryNumber = 1 + for atom in molecule.atoms: + if not molecule.isAtomInCycle(atom): + symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) + self.assertEqual(symmetryNumber, 18) + + def testAtomSymmetryNumberIsobutane(self): + """ + Test the Molecule.calculateAtomSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CC(C)C') + symmetryNumber = 1 + for atom in molecule.atoms: + if not molecule.isAtomInCycle(atom): + symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) + self.assertEqual(symmetryNumber, 81) + + def testBondSymmetryNumberEthane(self): + """ + Test the Molecule.calculateBondSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CC') + symmetryNumber = 1 + for atom1 in molecule.bonds: + for atom2 in molecule.bonds[atom1]: + if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): + symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) + self.assertEqual(symmetryNumber, 2) + + def testBondSymmetryNumberPropane(self): + """ + Test the Molecule.calculateBondSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CCC') + symmetryNumber = 1 + for atom1 in molecule.bonds: + for atom2 in molecule.bonds[atom1]: + if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): + symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) + self.assertEqual(symmetryNumber, 1) + + def testBondSymmetryNumberButane(self): + """ + Test the Molecule.calculateBondSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CCCC') + symmetryNumber = 1 + for atom1 in molecule.bonds: + for atom2 in molecule.bonds[atom1]: + if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): + symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) + self.assertEqual(symmetryNumber, 2) + + def testBondSymmetryNumberEthylene(self): + """ + Test the Molecule.calculateBondSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C') + symmetryNumber = 1 + for atom1 in molecule.bonds: + for atom2 in molecule.bonds[atom1]: + if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): + symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) + self.assertEqual(symmetryNumber, 2) + + def testBondSymmetryNumberAcetylene(self): + """ + Test the Molecule.calculateBondSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C#C') + symmetryNumber = 1 + for atom1 in molecule.bonds: + for atom2 in molecule.bonds[atom1]: + if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): + symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) + self.assertEqual(symmetryNumber, 2) + + def testAxisSymmetryNumberEthylene(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 2) + + def testAxisSymmetryNumberPropadiene(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=C') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 2) + + def testAxisSymmetryNumberButatriene(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=C=C') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 2) + + def testAxisSymmetryNumberButatrienyl(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=C=[CH]') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 1) + + def testAxisSymmetryNumberPropadienyl(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=[C]') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 2) + + def testAxisSymmetryNumber12Butadienyl(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CC=C=[C]') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 1) + + def testAxisSymmetryNumber12Hexadienyl(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=CCCC') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 1) + + def testAxisSymmetryNumber1(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CC(C)=C=C(CC)CC') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 2) + + def testAxisSymmetryNumber2(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=C(C(C(C(C=C=C)=C=C)=C=C)=C=C)') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 2) + + def testAxisSymmetryNumber3(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=[C]C(C)(C)[C]=C=C') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 4) + + def testAxisSymmetryNumber4(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=C=O') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 2) + + def testAxisSymmetryNumber5(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('CC=C=C=O') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 1) + + def testAxisSymmetryNumber6(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=C=N') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 1) + + def testAxisSymmetryNumber7(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('C=C=C=[N]') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 2) + + def testAxisSymmetryOxygenSinglet(self): + """ + Test the Molecule.calculateAxisSymmtryNumber() method. + """ + molecule = Molecule().fromSMILES('O=O') + self.assertEqual(calculateAxisSymmetryNumber(molecule), 1) + +# def testCyclicSymmetryNumber(self): +# +# # cyclohexane +# molecule = Molecule().fromInChI('InChI=1/C6H12/c1-2-4-6-5-3-1/h1-6H2') +# molecule.makeHydrogensExplicit() +# symmetryNumber = molecule.calculateCyclicSymmetryNumber() +# self.assertEqual(symmetryNumber, 12) + + def testTotalSymmetryNumberEthane(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('CC').calculateSymmetryNumber(), 18) + + def testTotalSymmetryNumber1(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C=C=[C]C(C)(C)[C]=C=C').calculateSymmetryNumber(), '???') + + def testTotalSymmetryNumber2(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C(=CC(c1ccccc1)C([CH]CCCCCC)C=Cc1ccccc1)[CH]CCCCCC').calculateSymmetryNumber(), 1) + + def testSymmetryNumberHydroxyl(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('[OH]').calculateSymmetryNumber(), 1) + + def testSymmetryNumberOxygen(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('O=O').calculateSymmetryNumber(), 2) + + def testSymmetryNumberDicarbon(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('[C]#[C]').calculateSymmetryNumber(), 2) + + def testSymmetryNumberHydrogen(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('[H][H]').calculateSymmetryNumber(), 2) + + def testSymmetryNumberAcetylene(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C#C').calculateSymmetryNumber(), 2) + + def testSymmetryNumberButadiyne(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C#CC#C').calculateSymmetryNumber(), 2) + + def testSymmetryNumberMethane(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C').calculateSymmetryNumber(), 12) + + def testSymmetryNumberFormaldehyde(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C=O').calculateSymmetryNumber(), 2) + + def testSymmetryNumberMethyl(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('[CH3]').calculateSymmetryNumber(), 6) + + def testSymmetryNumberWater(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('O').calculateSymmetryNumber(), 2) + + def testSymmetryNumberEthylene(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C=C').calculateSymmetryNumber(), 4) + + def testSymmetryNumberEthenyl(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C=[CH]').calculateSymmetryNumber(), 1) + + def testSymmetryNumberCyclic(self): + """ + Test the Molecule.calculateSymmtryNumber() method. + """ + self.assertEqual(Molecule().fromSMILES('C1=C=C=1').calculateSymmetryNumber(), '6?') + +################################################################################ + +if __name__ == '__main__': + unittest.main( testRunner = unittest.TextTestRunner(verbosity=2) ) \ No newline at end of file From 287f161bd72c38780ea7d9fb15081d3feaac0c8b Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 11:17:09 -0400 Subject: [PATCH 08/23] Removed implicit hydrogens functionality from Molecule objects. The intent was to conserve memory in some situations by storing the number of hydrogen atoms adjacent to each heavy atom as an integer, instead of storing the hydrogen atoms explicitly as Atom objects. (The idea came from OpenBabel, which has this functionality.) However, this turned out to be far from nontrivial to implement when dealing with the many graph manipulation requirements of an RMG job, and I don't think the memory savings is worth the trouble any more. There was probably also a time savings in the graph isomorphism evaluation when comparing implicit to implicit, but I think we can get at least some of that savings back in other ways. --- rmgpy/molecule/group.py | 7 +- rmgpy/molecule/molecule.pxd | 17 +-- rmgpy/molecule/molecule.py | 201 ++++-------------------------- rmgpy/molecule/molecule_draw.py | 2 +- rmgpy/molecule/symmetry.py | 5 - unittest/molecule/moleculeTest.py | 90 ++----------- 6 files changed, 48 insertions(+), 274 deletions(-) diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index 7c388d00d1..d035e49a6a 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -70,8 +70,7 @@ class GroupAtom(Vertex): group if it matches *any* item in the list. However, the `radicalElectrons`, `spinMultiplicity`, and `charge` attributes are linked such that an atom must match values from the same index in each of these in - order to match. Unlike an :class:`Atom` object, an :class:`GroupAtom` - cannot store implicit hydrogen atoms. + order to match. """ def __init__(self, atomType=None, radicalElectrons=None, spinMultiplicity=None, charge=None, label=''): @@ -763,7 +762,7 @@ def fromAdjacencyList(adjlist, group=False, addH=False): if group: atom = GroupAtom(atomType, radicalElectrons, spinMultiplicity, [0 for e in radicalElectrons], label) else: - atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, 0, label) + atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, label) # Add the atom to the list atoms.append(atom) @@ -835,7 +834,7 @@ def fromAdjacencyList(adjlist, group=False, addH=False): order += orders[bond.order] count = valence - radical - int(order) for i in range(count): - a = Atom('H', 0, 1, 0, 0, '') + a = Atom('H', 0, 1, 0, '') b = Bond('S') newAtoms.append(a) bonds[atom][a] = b diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 4ecb0988b2..ac5b6cf53d 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -36,7 +36,6 @@ cdef class Atom(Vertex): cdef public Element element cdef public short radicalElectrons cdef public short spinMultiplicity - cdef public short implicitHydrogens cdef public short charge cdef public str label cdef public AtomType atomType @@ -106,10 +105,6 @@ cdef class Molecule(Graph): cpdef deleteHydrogens(self) - cpdef makeHydrogensImplicit(self) - - cpdef makeHydrogensExplicit(self) - cpdef clearLabeledAtoms(self) cpdef bint containsLabeledAtom(self, str label) @@ -118,9 +113,9 @@ cdef class Molecule(Graph): cpdef dict getLabeledAtoms(self) - cpdef bint isIsomorphic(self, Graph other0, dict initialMap=?) + cpdef bint isIsomorphic(self, Graph other, dict initialMap=?) - cpdef tuple findIsomorphism(self, Graph other0, dict initialMap=?) + cpdef tuple findIsomorphism(self, Graph other, dict initialMap=?) cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?) @@ -132,13 +127,13 @@ cdef class Molecule(Graph): cpdef draw(self, str path) - cpdef fromCML(self, str cmlstr, bint implicitH=?) + cpdef fromCML(self, str cmlstr) - cpdef fromInChI(self, str inchistr, bint implicitH=?) + cpdef fromInChI(self, str inchistr) - cpdef fromSMILES(self, str smilesstr, bint implicitH=?) + cpdef fromSMILES(self, str smilesstr) - cpdef fromOBMol(self, obmol, bint implicitH=?) + cpdef fromOBMol(self, obmol) cpdef fromAdjacencyList(self, str adjlist) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index ddf2593ee8..321784decc 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -56,7 +56,6 @@ class Atom(Vertex): `element` :class:`Element` The chemical element the atom represents `radicalElectrons` ``short`` The number of radical electrons `spinMultiplicity` ``short`` The spin multiplicity of the atom - `implicitHydrogens` ``short`` The number of implicit hydrogen atoms bonded to this atom `charge` ``short`` The formal charge of the atom `label` ``str`` A string label that can be used to tag individual atoms =================== =================== ==================================== @@ -66,7 +65,7 @@ class Atom(Vertex): e.g. ``atom.symbol`` instead of ``atom.element.symbol``. """ - def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, implicitHydrogens=0, charge=0, label=''): + def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge=0, label=''): Vertex.__init__(self) if isinstance(element, str): self.element = elements.__dict__[element] @@ -74,7 +73,6 @@ def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, implici self.element = element self.radicalElectrons = radicalElectrons self.spinMultiplicity = spinMultiplicity - self.implicitHydrogens = implicitHydrogens self.charge = charge self.label = label self.atomType = None @@ -93,7 +91,7 @@ def __repr__(self): """ Return a representation that can be used to reconstruct the object. """ - return 'Atom(element="{0}", radicalElectrons={1}, spinMultiplicity={2}, implicitHydrogens={3}, charge={4}, label="{5}")'.format(self.element, self.radicalElectrons, self.spinMultiplicity, self.implicitHydrogens, self.charge, self.label) + return 'Atom(element="{0}", radicalElectrons={1}, spinMultiplicity={2}, charge={3}, label="{4}")'.format(self.element, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label) def __reduce__(self): """ @@ -106,7 +104,7 @@ def __reduce__(self): 'sortingLabel': self.sortingLabel, 'atomType': self.atomType.label if self.atomType else None, } - return (Atom, (self.element.symbol, self.radicalElectrons, self.spinMultiplicity, self.implicitHydrogens, self.charge, self.label), d) + return (Atom, (self.element.symbol, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label), d) def __setstate__(self, d): """ @@ -141,7 +139,6 @@ def equivalent(self, other): return (self.element is atom.element and self.radicalElectrons == atom.radicalElectrons and self.spinMultiplicity == atom.spinMultiplicity and - self.implicitHydrogens == atom.implicitHydrogens and self.charge == atom.charge) elif isinstance(other, GroupAtom): cython.declare(a=AtomType, radical=cython.short, spin=cython.short, charge=cython.short) @@ -192,7 +189,7 @@ def copy(self): Generate a deep copy of the current atom. Modifying the attributes of the copy will not affect the original. """ - a = Atom(self.element, self.radicalElectrons, self.spinMultiplicity, self.implicitHydrogens, self.charge, self.label) + a = Atom(self.element, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label) a.atomType = self.atomType return a @@ -439,7 +436,6 @@ class Molecule(Graph): ======================= =========== ======================================== Attribute Type Description ======================= =========== ======================================== - `implicitHydrogens` ``bool`` ``True`` if the hydrogen atoms are stored implicitly, ``False`` if stored explicity `symmetryNumber` ``int`` The (estimated) external + internal symmetry number of the molecule ======================= =========== ======================================== @@ -447,12 +443,11 @@ class Molecule(Graph): `InChI` string representing the molecular structure. """ - def __init__(self, atoms=None, bonds=None, implicitH=False, symmetry=1, SMILES='', InChI=''): + def __init__(self, atoms=None, bonds=None, symmetry=1, SMILES='', InChI=''): Graph.__init__(self, atoms, bonds) - self.implicitHydrogens = implicitH self.symmetryNumber = symmetry - if SMILES != '': self.fromSMILES(SMILES, implicitH) - elif InChI != '': self.fromInChI(InChI, implicitH) + if SMILES != '': self.fromSMILES(SMILES) + elif InChI != '': self.fromInChI(InChI) def __str__(self): """ @@ -470,7 +465,7 @@ def __reduce__(self): """ A helper function used when pickling an object. """ - return (Molecule, (self.vertices, self.edges, self.implicitHydrogens, self.symmetryNumber)) + return (Molecule, (self.vertices, self.edges, self.symmetryNumber)) def __getAtoms(self): return self.vertices def __setAtoms(self, atoms): self.vertices = atoms @@ -558,7 +553,6 @@ def getMolecularWeight(self): H = elements.getElement('H') for atom in self.vertices: mass += atom.element.mass - mass += atom.implicitHydrogens * H.mass return mass def copy(self, deep=False): @@ -596,10 +590,8 @@ def split(self): def deleteHydrogens(self): """ - Irreversibly delete all non-labeled hydrogens, without incrementing - the "implicitHydrogens" count of the neighbouring atom, or updating - Connectivity Values or the implicitHydrogens flag. If there's nothing - but Hydrogens, it does nothing. + Irreversibly delete all non-labeled hydrogens without updating + connectivity values. If there's nothing but hydrogens, it does nothing. It destroys information; be careful with it. """ cython.declare(atom=Atom, neighbor=Atom, hydrogens=list) @@ -619,86 +611,6 @@ def deleteHydrogens(self): for atom in hydrogens: self.removeAtom(atom) - - def makeHydrogensImplicit(self): - """ - Convert all explicitly stored hydrogen atoms to be stored implicitly, - unless they are labeled atoms (then they remain explicit). - An implicit hydrogen atom is stored on the heavy atom it is connected - to as a single integer counter. This is done to save memory. - """ - - cython.declare(atom=Atom, neighbor=Atom, hydrogens=list) - - # Check that the structure contains at least one heavy atom - for atom in self.vertices: - if not atom.isHydrogen(): - break - else: - # No heavy atoms, so leave explicit - return - - # Count the hydrogen atoms on each non-hydrogen atom and set the - # `implicitHydrogens` attribute accordingly - hydrogens = [] - for atom in self.vertices: - if atom.isHydrogen() and atom.label == '': - neighbor = self.edges[atom].keys()[0] - neighbor.implicitHydrogens += 1 - hydrogens.append(atom) - - # Remove the hydrogen atoms from the structure - for atom in hydrogens: - self.removeAtom(atom) - - # The connectivity values are different in implicit and explicit mode, - # so reset them so they are recomputed when needed - if hydrogens: - self.resetConnectivityValues() - - # Set implicitHydrogens flag to True - self.implicitHydrogens = True - - def makeHydrogensExplicit(self): - """ - Convert all implicitly stored hydrogen atoms to be stored explicitly. - An explicit hydrogen atom is stored as its own atom in the graph, with - a single bond to the heavy atom it is attached to. This consumes more - memory, but may be required for certain tasks (e.g. subgraph matching). - """ - - cython.declare(atom=Atom, H=Atom, bond=Bond, hydrogens=list, numAtoms=cython.short) - - # Create new hydrogen atoms for each implicit hydrogen - hydrogens = [] - for atom in self.vertices: - while atom.implicitHydrogens > 0: - H = Atom(element='H') - bond = Bond(order='S') - hydrogens.append((H, atom, bond)) - atom.implicitHydrogens -= 1 - - # Add the hydrogens to the graph - numAtoms = len(self.vertices) - for H, atom, bond in hydrogens: - self.addAtom(H) - self.addBond(H, atom, bond) - H.atomType = getAtomType(H, {atom:bond}) - # If known, set the connectivity information - H.connectivity1 = 1 - H.connectivity2 = atom.connectivity1 - H.connectivity3 = atom.connectivity2 - H.sortingLabel = numAtoms - numAtoms += 1 - - # The connectivity values are different in implicit and explicit mode, - # so reset them so they are recomputed when needed - if hydrogens: - self.resetConnectivityValues() - - # Set implicitHydrogens flag to False - self.implicitHydrogens = False - def updateAtomTypes(self): """ Iterate through the atoms in the structure, checking their atom types @@ -748,7 +660,7 @@ def getLabeledAtoms(self): labeled[atom.label] = atom return labeled - def isIsomorphic(self, other0, initialMap=None): + def isIsomorphic(self, other, initialMap=None): """ Returns :data:`True` if two graphs are isomorphic and :data:`False` otherwise. The `initialMap` attribute can be used to specify a required @@ -756,27 +668,15 @@ def isIsomorphic(self, other0, initialMap=None): while the atoms of `other` are the values). The `other` parameter must be a :class:`Molecule` object, or a :class:`TypeError` is raised. """ - cython.declare(other=Molecule, selfImplicitH=cython.bint, otherImplicitH=cython.bint) # It only makes sense to compare a Molecule to a Molecule for full # isomorphism, so raise an exception if this is not what was requested - if not isinstance(other0, Molecule): - raise TypeError('Got a {0} object for parameter "other0", when a Molecule object is required.'.format(other0.__class__)) - other = other0 - # Ensure that both self and other have the same implicit hydrogen status - # If not, make them both explicit just to be safe - selfImplicitH = self.implicitHydrogens - otherImplicitH = other.implicitHydrogens - if not selfImplicitH or not otherImplicitH: - self.makeHydrogensExplicit() - other.makeHydrogensExplicit() + if not isinstance(other, Molecule): + raise TypeError('Got a {0} object for parameter "other", when a Molecule object is required.'.format(other.__class__)) # Do the isomorphism comparison result = Graph.isIsomorphic(self, other, initialMap) - # Restore implicit status if needed - if selfImplicitH and not self.implicitHydrogens: self.makeHydrogensImplicit() - if otherImplicitH and not other.implicitHydrogens: other.makeHydrogensImplicit() return result - def findIsomorphism(self, other0, initialMap=None): + def findIsomorphism(self, other, initialMap=None): """ Returns :data:`True` if `other` is isomorphic and :data:`False` otherwise, and the matching mapping. The `initialMap` attribute can be @@ -786,24 +686,12 @@ def findIsomorphism(self, other0, initialMap=None): and the atoms of `other` for the values. The `other` parameter must be a :class:`Molecule` object, or a :class:`TypeError` is raised. """ - cython.declare(other=Molecule, selfImplicitH=cython.bint, otherImplicitH=cython.bint) # It only makes sense to compare a Molecule to a Molecule for full # isomorphism, so raise an exception if this is not what was requested - if not isinstance(other0, Molecule): - raise TypeError('Got a {0} object for parameter "other0", when a Molecule object is required.'.format(other0.__class__)) - other = other0 - # Ensure that both self and other have the same implicit hydrogen status - # If not, make them both explicit just to be safe - selfImplicitH = self.implicitHydrogens - otherImplicitH = other.implicitHydrogens - if not selfImplicitH or not otherImplicitH: - self.makeHydrogensExplicit() - other.makeHydrogensExplicit() + if not isinstance(other, Molecule): + raise TypeError('Got a {0} object for parameter "other", when a Molecule object is required.'.format(other.__class__)) # Do the isomorphism comparison result = Graph.findIsomorphism(self, other, initialMap) - # Restore implicit status if needed - if selfImplicitH and not self.implicitHydrogens: self.makeHydrogensImplicit() - if otherImplicitH and not other.implicitHydrogens: other.makeHydrogensImplicit() return result def isSubgraphIsomorphic(self, other, initialMap=None): @@ -818,13 +706,8 @@ def isSubgraphIsomorphic(self, other, initialMap=None): # isomorphism, so raise an exception if this is not what was requested if not isinstance(other, Group): raise TypeError('Got a {0} object for parameter "other", when a Molecule object is required.'.format(other.__class__)) - # Ensure that self is explicit (assume other is explicit) - implicitH = self.implicitHydrogens - if implicitH: self.makeHydrogensExplicit() # Do the isomorphism comparison result = Graph.isSubgraphIsomorphic(self, other, initialMap) - # Restore implicit status if needed - if implicitH and not self.implicitHydrogens: self.makeHydrogensImplicit() return result def findSubgraphIsomorphisms(self, other, initialMap=None): @@ -842,13 +725,8 @@ def findSubgraphIsomorphisms(self, other, initialMap=None): # isomorphism, so raise an exception if this is not what was requested if not isinstance(other, Group): raise TypeError('Got a {0} object for parameter "other", when a Molecule object is required.'.format(other.__class__)) - # Ensure that self is explicit (assume other is explicit) - implicitH = self.implicitHydrogens - if implicitH: self.makeHydrogensExplicit() # Do the isomorphism comparison result = Graph.findSubgraphIsomorphisms(self, other, initialMap) - # Restore implicit status if needed - if implicitH and not self.implicitHydrogens: self.makeHydrogensImplicit() return result def isAtomInCycle(self, atom): @@ -889,7 +767,7 @@ def _repr_png_(self): return png - def fromCML(self, cmlstr, implicitH=False): + def fromCML(self, cmlstr): """ Convert a string of CML `cmlstr` to a molecular structure. Uses `OpenBabel `_ to perform the conversion. @@ -897,10 +775,10 @@ def fromCML(self, cmlstr, implicitH=False): import pybel cmlstr = cmlstr.replace('\t', '') mol = pybel.readstring('cml', cmlstr) - self.fromOBMol(mol.OBMol, implicitH) + self.fromOBMol(mol.OBMol) return self - def fromInChI(self, inchistr, implicitH=False): + def fromInChI(self, inchistr): """ Convert an InChI string `inchistr` to a molecular structure. Uses `OpenBabel `_ to perform the conversion. @@ -910,10 +788,10 @@ def fromInChI(self, inchistr, implicitH=False): return self.fromAdjacencyList('1 H 1') import pybel mol = pybel.readstring('inchi', inchistr) - self.fromOBMol(mol.OBMol, implicitH) + self.fromOBMol(mol.OBMol) return self - def fromSMILES(self, smilesstr, implicitH=False): + def fromSMILES(self, smilesstr): """ Convert a SMILES string `smilesstr` to a molecular structure. Uses `OpenBabel `_ to perform the conversion. @@ -923,10 +801,10 @@ def fromSMILES(self, smilesstr, implicitH=False): return self.fromAdjacencyList('1 H 1') import pybel mol = pybel.readstring('smiles', smilesstr) - self.fromOBMol(mol.OBMol, implicitH) + self.fromOBMol(mol.OBMol) return self - def fromOBMol(self, obmol, implicitH=False): + def fromOBMol(self, obmol): """ Convert an OpenBabel OBMol object `obmol` to a molecular structure. Uses `OpenBabel `_ to perform the conversion. @@ -969,7 +847,7 @@ def fromOBMol(self, obmol, implicitH=False): # Process charge charge = obatom.GetFormalCharge() - atom = Atom(element, radicalElectrons, spinMultiplicity, 0, charge) + atom = Atom(element, radicalElectrons, spinMultiplicity, charge) self.vertices.append(atom) self.edges[atom] = {} @@ -996,9 +874,6 @@ def fromOBMol(self, obmol, implicitH=False): self.updateConnectivityValues() self.updateAtomTypes() - # Make hydrogens implicit to conserve memory - if implicitH: self.makeHydrogensImplicit() - return self def fromAdjacencyList(self, adjlist): @@ -1010,7 +885,6 @@ def fromAdjacencyList(self, adjlist): self.vertices, self.edges = fromAdjacencyList(adjlist, False, True) self.updateConnectivityValues() self.updateAtomTypes() - self.makeHydrogensImplicit() return self def toCML(self): @@ -1098,14 +972,9 @@ def toOBMol(self): import openbabel - cython.declare(implicitH=cython.bint) cython.declare(atom=Atom, atom1=Atom, bonds=dict, atom2=Atom, bond=Bond) cython.declare(index1=cython.int, index2=cython.int, order=cython.int) - # Make hydrogens explicit while we perform the conversion - implicitH = self.implicitHydrogens - self.makeHydrogensExplicit() - # Sort the atoms before converting to ensure output is consistent # between different runs self.sortAtoms() @@ -1129,19 +998,13 @@ def toOBMol(self): obmol.AssignSpinMultiplicity(True) - # Restore implicit hydrogens if necessary - if implicitH: self.makeHydrogensImplicit() - return obmol def toAdjacencyList(self, label='', removeH=False): """ Convert the molecular structure to a string adjacency list. """ - implicitH = self.implicitHydrogens - self.makeHydrogensExplicit() result = toAdjacencyList(self, label=label, group=False, removeH=removeH) - if implicitH: self.makeHydrogensImplicit() return result def isLinear(self): @@ -1150,7 +1013,7 @@ def isLinear(self): otherwise. """ - atomCount = len(self.vertices) + sum([atom.implicitHydrogens for atom in self.vertices]) + atomCount = len(self.vertices) # Monatomic molecules are definitely nonlinear if atomCount == 1: @@ -1165,15 +1028,12 @@ def isLinear(self): # True if all bonds are double bonds (e.g. O=C=O) allDoubleBonds = True for atom1 in self.edges: - if atom1.implicitHydrogens > 0: allDoubleBonds = False for bond in self.edges[atom1].values(): if not bond.isDouble(): allDoubleBonds = False if allDoubleBonds: return True # True if alternating single-triple bonds (e.g. H-C#C-H) # This test requires explicit hydrogen atoms - implicitH = self.implicitHydrogens - self.makeHydrogensExplicit() for atom in self.vertices: bonds = self.edges[atom].values() if len(bonds)==1: @@ -1187,11 +1047,9 @@ def isLinear(self): break # fail if we haven't continued else: # didn't fail - if implicitH: self.makeHydrogensImplicit() return True # not returned yet? must be nonlinear - if implicitH: self.makeHydrogensImplicit() return False def countInternalRotors(self): @@ -1204,7 +1062,7 @@ def countInternalRotors(self): for atom1 in self.edges: for atom2, bond in self.edges[atom1].iteritems(): if self.vertices.index(atom1) < self.vertices.index(atom2) and bond.isSingle() and not self.isBondInCycle(atom1, atom2): - if len(self.edges[atom1]) + atom1.implicitHydrogens > 1 and len(self.edges[atom2]) + atom2.implicitHydrogens > 1: + if len(self.edges[atom1]) > 1 and len(self.edges[atom2]) > 1: count += 1 return count @@ -1221,9 +1079,6 @@ def generateResonanceIsomers(self): """ Generate and return all of the resonance isomers of this molecule. """ - - implicitH = self.implicitHydrogens - self.makeHydrogensExplicit() isomers = [self] @@ -1244,10 +1099,6 @@ def generateResonanceIsomers(self): isomers.append(newIsomer) # Move to next resonance isomer index += 1 - - if implicitH: - for isomer in isomers: - isomer.makeHydrogensImplicit() return isomers diff --git a/rmgpy/molecule/molecule_draw.py b/rmgpy/molecule/molecule_draw.py index 7bb762b76a..0364473570 100644 --- a/rmgpy/molecule/molecule_draw.py +++ b/rmgpy/molecule/molecule_draw.py @@ -1138,7 +1138,7 @@ def drawMolecule(molecule, path=None, surface=''): # We will delete them from the *copied* list of atoms, and store them here: implicitHydrogensToDraw = {} for atom in molecule.atoms: - implicitHydrogensToDraw[atom] = atom.implicitHydrogens + implicitHydrogensToDraw[atom] = 0 atoms = molecule.atoms[:] # bonds = molecule.bonds.copy() is too shallow for a dict-of-dicts, diff --git a/rmgpy/molecule/symmetry.py b/rmgpy/molecule/symmetry.py index aafcc8eea5..3dd8810576 100644 --- a/rmgpy/molecule/symmetry.py +++ b/rmgpy/molecule/symmetry.py @@ -394,9 +394,6 @@ def calculateSymmetryNumber(molecule): """ symmetryNumber = 1 - implicitH = molecule.implicitHydrogens - molecule.makeHydrogensExplicit() - for atom in molecule.vertices: if not molecule.isAtomInCycle(atom): symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) @@ -411,6 +408,4 @@ def calculateSymmetryNumber(molecule): #if molecule.isCyclic(): # symmetryNumber *= calculateCyclicSymmetryNumber(molecule) - if implicitH: molecule.makeHydrogensImplicit() - return symmetryNumber diff --git a/unittest/molecule/moleculeTest.py b/unittest/molecule/moleculeTest.py index 2f40425815..e73d3a74d1 100644 --- a/unittest/molecule/moleculeTest.py +++ b/unittest/molecule/moleculeTest.py @@ -18,7 +18,7 @@ def setUp(self): """ A method called before each unit test in this class. """ - self.atom = Atom(element=getElement('C'), radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=3, charge=0, label='*1') + self.atom = Atom(element=getElement('C'), radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') def testMass(self): """ @@ -43,7 +43,7 @@ def testIsHydrogen(self): Test the Atom.isHydrogen() method. """ for element in elementList: - atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=0, charge=0, label='*1') + atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') if element.symbol == 'H': self.assertTrue(atom.isHydrogen()) else: @@ -54,7 +54,7 @@ def testIsNonHydrogen(self): Test the Atom.isNonHydrogen() method. """ for element in elementList: - atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=0, charge=0, label='*1') + atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') if element.symbol == 'H': self.assertFalse(atom.isNonHydrogen()) else: @@ -65,7 +65,7 @@ def testIsCarbon(self): Test the Atom.isCarbon() method. """ for element in elementList: - atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=0, charge=0, label='*1') + atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') if element.symbol == 'C': self.assertTrue(atom.isCarbon()) else: @@ -76,7 +76,7 @@ def testIsOxygen(self): Test the Atom.isOxygen() method. """ for element in elementList: - atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=0, charge=0, label='*1') + atom = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') if element.symbol == 'O': self.assertTrue(atom.isOxygen()) else: @@ -108,13 +108,12 @@ def testApplyActionBreakBond(self): """ action = ['BREAK_BOND', '*1', 'S', '*2'] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=3, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') atom = atom0.copy() atom.applyAction(action) self.assertEqual(atom0.element, atom.element) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons) self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(atom0.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -124,13 +123,12 @@ def testApplyActionFormBond(self): """ action = ['FORM_BOND', '*1', 'S', '*2'] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=3, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') atom = atom0.copy() atom.applyAction(action) self.assertEqual(atom0.element, atom.element) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons) self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(atom0.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -140,13 +138,12 @@ def testApplyActionIncrementBond(self): """ action = ['CHANGE_BOND', '*1', 1, '*2'] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=3, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') atom = atom0.copy() atom.applyAction(action) self.assertEqual(atom0.element, atom.element) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons) self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(atom0.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -156,13 +153,12 @@ def testApplyActionDecrementBond(self): """ action = ['CHANGE_BOND', '*1', -1, '*2'] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=3, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') atom = atom0.copy() atom.applyAction(action) self.assertEqual(atom0.element, atom.element) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons) self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(atom0.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -172,13 +168,12 @@ def testApplyActionGainRadical(self): """ action = ['GAIN_RADICAL', '*1', 1] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=3, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') atom = atom0.copy() atom.applyAction(action) self.assertEqual(atom0.element, atom.element) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons - 1) self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity - 1) - self.assertEqual(atom0.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -188,13 +183,12 @@ def testApplyActionLoseRadical(self): """ action = ['LOSE_RADICAL', '*1', 1] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, implicitHydrogens=3, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') atom = atom0.copy() atom.applyAction(action) self.assertEqual(atom0.element, atom.element) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons + 1) self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity + 1) - self.assertEqual(atom0.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -235,7 +229,6 @@ def testCopy(self): self.assertEqual(self.atom.atomType, atom.atomType) self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons) self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(self.atom.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(self.atom.charge, atom.charge) self.assertEqual(self.atom.label, atom.label) @@ -250,7 +243,6 @@ def testPickle(self): self.assertEqual(self.atom.atomType, atom.atomType) self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons) self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(self.atom.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(self.atom.charge, atom.charge) self.assertEqual(self.atom.label, atom.label) @@ -264,7 +256,6 @@ def testOutput(self): self.assertEqual(self.atom.atomType, atom.atomType) self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons) self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(self.atom.implicitHydrogens, atom.implicitHydrogens) self.assertEqual(self.atom.charge, atom.charge) self.assertEqual(self.atom.label, atom.label) @@ -567,7 +558,7 @@ def testFromAdjacencyList(self): """ Test the Molecule.fromAdjacencyList() method. """ - atom1, atom2, atom3 = self.molecule.atoms + atom1, atom2, atom3 = self.molecule.atoms[0:3] self.assertTrue(self.molecule.hasBond(atom1,atom2)) self.assertTrue(self.molecule.hasBond(atom1,atom3)) self.assertFalse(self.molecule.hasBond(atom2,atom3)) @@ -578,21 +569,18 @@ def testFromAdjacencyList(self): self.assertTrue(atom1.element.symbol == 'C') self.assertTrue(atom1.radicalElectrons == 1) self.assertTrue(atom1.spinMultiplicity == 2) - self.assertTrue(atom1.implicitHydrogens == 0) self.assertTrue(atom1.charge == 0) self.assertTrue(atom2.label == '*1') self.assertTrue(atom2.element.symbol == 'O') self.assertTrue(atom2.radicalElectrons == 0) self.assertTrue(atom2.spinMultiplicity == 1) - self.assertTrue(atom2.implicitHydrogens == 0) self.assertTrue(atom2.charge == 0) self.assertTrue(atom3.label == '') self.assertTrue(atom3.element.symbol == 'C') self.assertTrue(atom3.radicalElectrons == 0) self.assertTrue(atom3.spinMultiplicity == 1) - self.assertTrue(atom3.implicitHydrogens == 3) self.assertTrue(atom3.charge == 0) self.assertTrue(bond12.isDouble()) @@ -663,8 +651,6 @@ def testSubgraphIsomorphismAgain(self): 4 H 0 {1,S} """) - molecule.makeHydrogensExplicit() - labeled1 = molecule.getLabeledAtoms().values()[0] labeled2 = group.getLabeledAtoms().values()[0] @@ -725,24 +711,6 @@ def testAdjacencyList(self): 6 C 0 {5,S} """) molecule2 = Molecule().fromSMILES('C=CC=C[CH]C') - - molecule1.makeHydrogensExplicit() - molecule2.makeHydrogensExplicit() - self.assertTrue(molecule1.isIsomorphic(molecule2)) - self.assertTrue(molecule2.isIsomorphic(molecule1)) - - molecule1.makeHydrogensImplicit() - molecule2.makeHydrogensImplicit() - self.assertTrue(molecule1.isIsomorphic(molecule2)) - self.assertTrue(molecule2.isIsomorphic(molecule1)) - - molecule1.makeHydrogensExplicit() - molecule2.makeHydrogensImplicit() - self.assertTrue(molecule1.isIsomorphic(molecule2)) - self.assertTrue(molecule2.isIsomorphic(molecule1)) - - molecule1.makeHydrogensImplicit() - molecule2.makeHydrogensExplicit() self.assertTrue(molecule1.isIsomorphic(molecule2)) self.assertTrue(molecule2.isIsomorphic(molecule1)) @@ -784,40 +752,6 @@ def testIsInCycleCyclohexane(self): self.assertTrue(molecule.isBondInCycle(atom1, atom2)) else: self.assertFalse(molecule.isBondInCycle(atom1, atom2)) - - def testImplicitHydrogens(self): - """ - Test that a molecule can be converted to and from implicit hydrogen - mode with no loss of data. - """ - self.assertTrue(self.molecule.implicitHydrogens) - self.assertEqual(len(self.molecule.atoms), 3) - self.assertEqual(self.molecule.atoms[0].implicitHydrogens, 0) - self.assertEqual(self.molecule.atoms[1].implicitHydrogens, 0) - self.assertEqual(self.molecule.atoms[2].implicitHydrogens, 3) - self.assertEqual(sum([1 for atom in self.molecule.atoms if atom.isHydrogen()]), 0) - self.assertEqual(self.molecule.getFormula(), 'C2H3O') - self.assertAlmostEqual(self.molecule.getMolecularWeight() * 1000, 43.04, 2) - - self.molecule.makeHydrogensExplicit() - self.assertFalse(self.molecule.implicitHydrogens) - self.assertEqual(len(self.molecule.atoms), 6) - self.assertEqual(self.molecule.atoms[0].implicitHydrogens, 0) - self.assertEqual(self.molecule.atoms[1].implicitHydrogens, 0) - self.assertEqual(self.molecule.atoms[2].implicitHydrogens, 0) - self.assertEqual(sum([1 for atom in self.molecule.atoms if atom.isHydrogen()]), 3) - self.assertEqual(self.molecule.getFormula(), 'C2H3O') - self.assertAlmostEqual(self.molecule.getMolecularWeight() * 1000, 43.04, 2) - - self.molecule.makeHydrogensImplicit() - self.assertTrue(self.molecule.implicitHydrogens) - self.assertEqual(len(self.molecule.atoms), 3) - self.assertEqual(self.molecule.atoms[0].implicitHydrogens, 3) - self.assertEqual(self.molecule.atoms[1].implicitHydrogens, 0) - self.assertEqual(self.molecule.atoms[2].implicitHydrogens, 0) - self.assertEqual(sum([1 for atom in self.molecule.atoms if atom.isHydrogen()]), 0) - self.assertEqual(self.molecule.getFormula(), 'C2H3O') - self.assertAlmostEqual(self.molecule.getMolecularWeight() * 1000, 43.04, 2) def testFromSMILESH(self): """ From 6573ee9964d69c0ff1aaee4d023d38e7ab31197c Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 11:39:37 -0400 Subject: [PATCH 09/23] Removed assorted references to implicit/explicit hydrogens. Since this functionality no longer exists, we need to remove these references for RMG-Py to run. --- rmgpy/chemkin.py | 7 ------- rmgpy/data/kinetics.py | 2 -- rmgpy/data/thermo.py | 17 +---------------- rmgpy/rmg/model.py | 22 ---------------------- 4 files changed, 1 insertion(+), 47 deletions(-) diff --git a/rmgpy/chemkin.py b/rmgpy/chemkin.py index bb1d588664..23e37db719 100755 --- a/rmgpy/chemkin.py +++ b/rmgpy/chemkin.py @@ -843,13 +843,6 @@ def writeThermoEntry(species): elementCounts.append(1) else: elementCounts[elements.index(symbol)] += 1 - # Also handle implicit hydrogen atoms - symbol = 'H' - if symbol not in elements: - elements.append(symbol) - elementCounts.append(atom.implicitHydrogens) - else: - elementCounts[elements.index(symbol)] += atom.implicitHydrogens # Remove elements with zero count index = 0 while index < len(elementCounts): diff --git a/rmgpy/data/kinetics.py b/rmgpy/data/kinetics.py index ca4734d51a..07eabf3216 100644 --- a/rmgpy/data/kinetics.py +++ b/rmgpy/data/kinetics.py @@ -2491,8 +2491,6 @@ def __generateReactions(self, reactants, forward=True): for i in range(len(reactants)): for j in range(len(reactants[i])): reactants[i][j] = reactants[i][j].copy(deep=True) - # Each molecule must have explicit hydrogen atoms - reactants[i][j].makeHydrogensExplicit() if forward: template = self.forwardTemplate diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index affcddfc8d..394e6d90ae 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -606,12 +606,7 @@ def getThermoDataFromGroups(self, species): :class:`Species` object `species` by estimation using the group additivity values. If no group additivity values are loaded, a :class:`DatabaseError` is raised. - """ - # Ensure molecules are using explicit hydrogens - implicitH = [mol.implicitHydrogens for mol in species.molecule] - for molecule in species.molecule: - molecule.makeHydrogensExplicit() - + """ thermo = [] for molecule in species.molecule: molecule.clearLabeledAtoms() @@ -623,11 +618,6 @@ def getThermoDataFromGroups(self, species): indices = H298.argsort() species.molecule = [species.molecule[ind] for ind in indices] - implicitH = [implicitH[ind] for ind in indices] - - # Restore implicit hydrogens if necessary - for implicit, molecule in zip(implicitH, species.molecule): - if implicit: molecule.makeHydrogensImplicit() return (thermo[indices[0]], None, None) @@ -638,9 +628,6 @@ def estimateThermoViaGroupAdditivity(self, molecule): additivity values. If no group additivity values are loaded, a :class:`DatabaseError` is raised. """ - implicitH = molecule.implicitHydrogens - molecule.makeHydrogensExplicit() - # For thermo estimation we need the atoms to already be sorted because we # iterate over them; if the order changes during the iteration then we # will probably not visit the right atoms, and so will get the thermo wrong @@ -767,8 +754,6 @@ def estimateThermoViaGroupAdditivity(self, molecule): molecule.calculateSymmetryNumber() thermoData.S298.value -= constants.R * math.log(molecule.symmetryNumber) - if implicitH: molecule.makeHydrogensImplicit() - return thermoData def __addThermoData(self, thermoData1, thermoData2): diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 755c18154a..265e75b2c4 100755 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -72,12 +72,6 @@ def generateThermoData(self, database, thermoClass=MultiNASA): Generates the thermo data for each structure (resonance isomer), picks that with lowest H298 value, and saves it to `self.thermoData`. """ - - # Ensure molecules are using explicit hydrogens - implicitH = [mol.implicitHydrogens for mol in self.molecule] - for molecule in self.molecule: - molecule.makeHydrogensExplicit() - # Get the thermo data for the species from the database thermo0 = database.thermo.getThermoData(self) @@ -112,10 +106,6 @@ def generateThermoData(self, database, thermoClass=MultiNASA): err = math.sqrt(numpy.sum((self.thermo.getHeatCapacities(Tlist) - thermo0.getHeatCapacities(Tlist))**2)/len(Tlist))/constants.R logging.log(logging.WARNING if err > 0.1 else 0, 'Average RMS error in heat capacity fit to {0} = {1:g}*R'.format(self, err)) - # Restore implicit hydrogens if necessary - for implicit, molecule in zip(implicitH, self.molecule): - if implicit: molecule.makeHydrogensImplicit() - return self.thermo def generateStatesData(self, database): @@ -127,10 +117,7 @@ def generateStatesData(self, database): if not self.thermo: raise Exception("Unable to determine states model for species {0}: No thermodynamics model found.".format(self)) molecule = self.molecule[0] - implicitH = molecule.implicitHydrogens - molecule.makeHydrogensExplicit() self.states = database.states.getStatesData(molecule, self.thermo) - if implicitH: molecule.makeHydrogensImplicit() def generateLennardJonesParameters(self): """ @@ -264,7 +251,6 @@ def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True) molecule = object molecule.clearLabeledAtoms() - molecule.makeHydrogensImplicit() # If desired, check to ensure that the species is new; return the # existing species if not new @@ -294,11 +280,6 @@ def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True) else: self.speciesDict[formula] = [spec] - # Store hydrogens implicitly to conserve memory and speed up isomorphism - for mol in spec.molecule: - #mol.updateConnectivityValues() - mol.makeHydrogensImplicit() - self.speciesCounter += 1 # Since the species is new, add it to the list of new species @@ -490,15 +471,12 @@ def react(self, database, speciesA, speciesB=None): for moleculeA in speciesA.molecule: reactionList.extend(database.kinetics.generateReactions([moleculeA])) moleculeA.clearLabeledAtoms() - moleculeA.makeHydrogensImplicit() else: for moleculeA in speciesA.molecule: for moleculeB in speciesB.molecule: reactionList.extend(database.kinetics.generateReactions([moleculeA, moleculeB])) moleculeA.clearLabeledAtoms() - moleculeA.makeHydrogensImplicit() moleculeB.clearLabeledAtoms() - moleculeB.makeHydrogensImplicit() return reactionList def enlarge(self, newObject): From 1033d46bc705d39e65cb9cc37c39156bf3156cc2 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 16:04:04 -0400 Subject: [PATCH 10/23] Now storing dict of edges on Vertex and vertices on Edge. Previously we stored the edges in a dict of dicts on the Graph, and did not store the vertices on the edge at all. However, the convenience of having these attributes on the Vertex and Edge objects outweights the slight increase in memory use. This is a pretty significant change in approach, and many of the methods of the Graph class needed modification as a result. --- rmgpy/molecule/graph.pxd | 21 +- rmgpy/molecule/graph.py | 314 +++++++++++------------ unittest/molecule/graphTest.py | 453 +++++++++++++++++---------------- 3 files changed, 397 insertions(+), 391 deletions(-) diff --git a/rmgpy/molecule/graph.pxd b/rmgpy/molecule/graph.pxd index 2dc0db7b7c..e85721dec1 100644 --- a/rmgpy/molecule/graph.pxd +++ b/rmgpy/molecule/graph.pxd @@ -26,6 +26,8 @@ cdef class Vertex(object): + cdef readonly dict edges + cdef public short connectivity1 cdef public short connectivity2 cdef public short connectivity3 @@ -45,32 +47,37 @@ cpdef short getVertexSortingLabel(Vertex vertex) except -1 # all values should b cdef class Edge(object): - cpdef bint equivalent(Edge self, Edge other) + cdef readonly Vertex vertex1, vertex2 + + cpdef Edge copy(self) + + cpdef bint equivalent(Edge self, Edge other) except -2 + + cpdef bint isSpecificCaseOf(self, Edge other) except -2 - cpdef bint isSpecificCaseOf(self, Edge other) + cpdef Vertex getOtherVertex(self, Vertex vertex) ################################################################################ cdef class Graph: cdef public list vertices - cdef public dict edges cpdef Vertex addVertex(self, Vertex vertex) - cpdef Edge addEdge(self, Vertex vertex1, Vertex vertex2, Edge edge) + cpdef Edge addEdge(self, Edge edge) cpdef dict getEdges(self, Vertex vertex) cpdef Edge getEdge(self, Vertex vertex1, Vertex vertex2) - cpdef bint hasVertex(self, Vertex vertex) + cpdef bint hasVertex(self, Vertex vertex) except -2 - cpdef bint hasEdge(self, Vertex vertex1, Vertex vertex2) + cpdef bint hasEdge(self, Vertex vertex1, Vertex vertex2) except -2 cpdef removeVertex(self, Vertex vertex) - cpdef removeEdge(self, Vertex vertex1, Vertex vertex2) + cpdef removeEdge(self, Edge edge) cpdef Graph copy(self, bint deep=?) diff --git a/rmgpy/molecule/graph.py b/rmgpy/molecule/graph.py index ff5cdb42e5..ac02052bb1 100644 --- a/rmgpy/molecule/graph.py +++ b/rmgpy/molecule/graph.py @@ -58,6 +58,7 @@ class Vertex(object): """ def __init__(self): + self.edges = {} self.resetConnectivityValues() def __reduce__(self): @@ -65,6 +66,7 @@ def __reduce__(self): A helper function used when pickling an object. """ d = { + 'edges': self.edges, 'connectivity1': self.connectivity1, 'connectivity2': self.connectivity2, 'connectivity3': self.connectivity3, @@ -73,6 +75,7 @@ def __reduce__(self): return (Vertex, (), d) def __setstate__(self, d): + self.edges = d['edges'] self.connectivity1 = d['connectivity1'] self.connectivity2 = d['connectivity2'] self.connectivity3 = d['connectivity3'] @@ -128,14 +131,15 @@ class Edge(object): in the derived class. """ - def __init__(self): - pass + def __init__(self, vertex1, vertex2): + self.vertex1 = vertex1 + self.vertex2 = vertex2 def __reduce__(self): """ A helper function used when pickling an object. """ - return (Edge, ()) + return (Edge, (self.vertex1, self.vertex2)) def equivalent(self, other): """ @@ -153,6 +157,19 @@ class if your edges have semantic information. """ return True + def getOtherVertex(self, vertex): + """ + Given a vertex that makes up part of the edge, return the other vertex. + Raise a :class:`ValueError` if the given vertex is not part of the + edge. + """ + if self.vertex1 is not vertex and self.vertex2 is not vertex: + raise ValueError('The given vertex is not one of the vertices of this edge.') + elif self.vertex1 is vertex: + return self.vertex2 + elif self.vertex2 is vertex: + return self.vertex1 + ################################################################################ class Graph: @@ -166,44 +183,48 @@ class Graph: or the :meth:`getEdges` method. """ - def __init__(self, vertices=None, edges=None): + def __init__(self, vertices=None): self.vertices = vertices or [] - self.edges = edges or {} def __reduce__(self): """ A helper function used when pickling an object. """ - return (Graph, (self.vertices, self.edges)) + return (Graph, (self.vertices,)) def addVertex(self, vertex): """ Add a `vertex` to the graph. The vertex is initialized with no edges. """ self.vertices.append(vertex) - self.edges[vertex] = dict() + vertex.edges = dict() return vertex - def addEdge(self, vertex1, vertex2, edge): + def addEdge(self, edge): """ - Add an `edge` to the graph as an edge connecting the two vertices - `vertex1` and `vertex2`. + Add an `edge` to the graph. The two vertices in the edge must already + exist in the graph, or a :class:`ValueError` is raised. """ - self.edges[vertex1][vertex2] = edge - self.edges[vertex2][vertex1] = edge + if edge.vertex1 not in self.vertices or edge.vertex2 not in self.vertices: + raise ValueError('Attempted to add edge between vertices not in the graph.') + edge.vertex1.edges[edge.vertex2] = edge + edge.vertex2.edges[edge.vertex1] = edge return edge def getEdges(self, vertex): """ Return a list of the edges involving the specified `vertex`. """ - return self.edges[vertex] + return vertex.edges def getEdge(self, vertex1, vertex2): """ Returns the edge connecting vertices `vertex1` and `vertex2`. """ - return self.edges[vertex1][vertex2] + try: + return vertex1.edges[vertex2] + except KeyError: + raise ValueError('The specified vertices are not connected by an edge in this graph.') def hasVertex(self, vertex): """ @@ -217,7 +238,7 @@ def hasEdge(self, vertex1, vertex2): Returns ``True`` if vertices `vertex1` and `vertex2` are connected by an edge, or ``False`` if not. """ - return vertex2 in self.edges[vertex1] if vertex1 in self.edges else False + return vertex1 in self.vertices and vertex2 in vertex1.edges def removeVertex(self, vertex): """ @@ -225,21 +246,20 @@ def removeVertex(self, vertex): not remove vertices that no longer have any edges as a result of this removal. """ - for vertex2 in self.vertices: - if vertex2 is not vertex: - if vertex in self.edges[vertex2]: - del self.edges[vertex2][vertex] - del self.edges[vertex] + cython.declare(vertex2=Vertex) + for vertex2 in vertex.edges: + del vertex2.edges[vertex] + vertex.edges = dict() self.vertices.remove(vertex) - def removeEdge(self, vertex1, vertex2): + def removeEdge(self, edge): """ Remove the edge having vertices `vertex1` and `vertex2` from the graph. Does not remove vertices that no longer have any edges as a result of this removal. """ - del self.edges[vertex1][vertex2] - del self.edges[vertex2][vertex1] + del edge.vertex1.edges[edge.vertex2] + del edge.vertex2.edges[edge.vertex1] def copy(self, deep=False): """ @@ -248,43 +268,50 @@ def copy(self, deep=False): If `deep` is ``False`` or not specified, a shallow copy is made: the original vertices and edges are used in the new graph. """ - other = cython.declare(Graph) + cython.declare(other=Graph) + cython.declare(vertex=Vertex, vertex1=Vertex, vertex2=Vertex) + cython.declare(edge=Edge) + cython.declare(edges=dict) + cython.declare(index1=cython.int, index2=cython.int) + other = Graph() for vertex in self.vertices: - other.addVertex(vertex.copy() if deep else vertex) + if deep: + other.addVertex(vertex.copy()) + else: + edges = vertex.edges + other.addVertex(vertex) + vertex.edges = edges for vertex1 in self.vertices: - for vertex2 in self.edges[vertex1]: + for vertex2 in vertex1.edges: if deep: index1 = self.vertices.index(vertex1) index2 = self.vertices.index(vertex2) - other.addEdge(other.vertices[index1], other.vertices[index2], - self.edges[vertex1][vertex2].copy()) - else: - other.addEdge(vertex1, vertex2, self.edges[vertex1][vertex2]) + edge = vertex1.edges[vertex2].copy() + edge.vertex1 = other.vertices[index1] + edge.vertex2 = other.vertices[index2] + other.addEdge(edge) return other def merge(self, other): """ Merge two graphs so as to store them in a single Graph object. """ - + cython.declare(new=Graph) + cython.declare(vertex=Vertex, vertex1=Vertex, vertex2=Vertex) + # Create output graph - new = cython.declare(Graph) new = Graph() # Add vertices to output graph for vertex in self.vertices: + edges = vertex.edges new.addVertex(vertex) + vertex.edges = edges for vertex in other.vertices: + edges = vertex.edges new.addVertex(vertex) - - # Add edges to output graph - for v1 in self.vertices: - for v2 in self.edges[v1]: - new.edges[v1][v2] = self.edges[v1][v2] - for v1 in other.vertices: - for v2 in other.edges[v1]: - new.edges[v1][v2] = other.edges[v1][v2] + vertex.edges = edges return new @@ -293,13 +320,12 @@ def split(self): Convert a single Graph object containing two or more unconnected graphs into separate graphs. """ - + cython.declare(new1=Graph, new2=Graph) + cython.declare(vertex=Vertex, vertex1=Vertex, vertex2=Vertex) + cython.declare(verticesToMove=list) + cython.declare(index=cython.int) + # Create potential output graphs - new1 = cython.declare(Graph) - new2 = cython.declare(Graph) - verticesToMove = cython.declare(list) - index = cython.declare(cython.int) - new1 = self.copy() new2 = Graph() @@ -312,30 +338,20 @@ def split(self): # Iterate until there are no more atoms to move index = 0 while index < len(verticesToMove): - for v2 in self.edges[verticesToMove[index]]: + for v2 in verticesToMove[index].edges: if v2 not in verticesToMove: verticesToMove.append(v2) index += 1 - + # If all atoms are to be moved, simply return new1 if len(new1.vertices) == len(verticesToMove): return [new1] - # Copy to new graph - for vertex in verticesToMove: - new2.addVertex(vertex) - for v1 in verticesToMove: - for v2, edge in new1.edges[v1].iteritems(): - new2.edges[v1][v2] = edge - - # Remove from old graph - for v1 in new2.vertices: - for v2 in new2.edges[v1]: - if v1 in verticesToMove and v2 in verticesToMove: - del new1.edges[v1][v2] + # Copy to new graph and remove from old graph for vertex in verticesToMove: - new1.removeVertex(vertex) - + new2.vertices.append(vertex) + new1.vertices.remove(vertex) + new = [new2] new.extend(new1.split()) return new @@ -353,22 +369,19 @@ def updateConnectivityValues(self): Update the connectivity values for each vertex in the graph. These are used to accelerate the isomorphism checking. """ - - cython.declare(count=cython.short, edges=dict) cython.declare(vertex1=Vertex, vertex2=Vertex) - + cython.declare(count=cython.short) + for vertex1 in self.vertices: - count = len(self.edges[vertex1]) + count = len(vertex1.edges) vertex1.connectivity1 = count for vertex1 in self.vertices: count = 0 - edges = self.edges[vertex1] - for vertex2 in edges: count += vertex2.connectivity1 + for vertex2 in vertex1.edges: count += vertex2.connectivity1 vertex1.connectivity2 = count for vertex1 in self.vertices: count = 0 - edges = self.edges[vertex1] - for vertex2 in edges: count += vertex2.connectivity2 + for vertex2 in vertex1.edges: count += vertex2.connectivity2 vertex1.connectivity3 = count def sortVertices(self): @@ -422,9 +435,10 @@ def findSubgraphIsomorphisms(self, other, initialMap=None): def isCyclic(self): """ - Return :data:`True` if one or more cycles are present in the structure - and :data:`False` otherwise. + Return ``True`` if one or more cycles are present in the graph or + ``False`` otherwise. """ + cython.declare(vertex=Vertex) for vertex in self.vertices: if self.isVertexInCycle(vertex): return True @@ -432,46 +446,46 @@ def isCyclic(self): def isVertexInCycle(self, vertex): """ - Return :data:`True` if `vertex` is in one or more cycles in the graph, - or :data:`False` if not. + Return ``True`` if the given `vertex` is contained in one or more + cycles in the graph, or ``False`` if not. """ - chain = cython.declare(list) - chain = [vertex] - return self.__isChainInCycle(chain) + return self.__isChainInCycle([vertex]) - def isEdgeInCycle(self, vertex1, vertex2): + def isEdgeInCycle(self, edge): """ Return :data:`True` if the edge between vertices `vertex1` and `vertex2` is in one or more cycles in the graph, or :data:`False` if not. """ - cycle_list = self.getAllCycles(vertex1) - for cycle in cycle_list: - if vertex2 in cycle: + cython.declare(cycles=list) + cycles = self.getAllCycles(edge.vertex1) + for cycle in cycles: + if edge.vertex2 in cycle: return True return False def __isChainInCycle(self, chain): """ - Is the `chain` in a cycle? - Returns True/False. - Recursively calls itself + Return ``True`` if the given `chain` of vertices is contained in one + or more cycles or ``False`` otherwise. This function recursively calls + itself. """ - # Note that this function no longer returns the cycle; just True/False - vertex2 = cython.declare(Vertex) - edge = cython.declare(Edge) - found = cython.declare(cython.bint) + cython.declare(vertex1=Vertex, vertex2=Vertex) + cython.declare(edge=Edge) - for vertex2, edge in self.edges[chain[-1]].iteritems(): + vertex1 = chain[-1] + for vertex2 in vertex1.edges: if vertex2 is chain[0] and len(chain) > 2: return True elif vertex2 not in chain: - # make the chain a little longer and explore again + # Make the chain a little longer and explore again chain.append(vertex2) - found = self.__isChainInCycle(chain) - if found: return True - # didn't find a cycle down this path (-vertex2), - # so remove the vertex from the chain - chain.remove(vertex2) + if self.__isChainInCycle(chain): + # We found a cycle, so the return value must be True + return True + else: + # We did not find a cycle down this path, so remove the vertex from the chain + chain.remove(vertex2) + # If we reach this point then we did not find any cycles involving this chain return False def getAllCycles(self, startingVertex): @@ -479,49 +493,31 @@ def getAllCycles(self, startingVertex): Given a starting vertex, returns a list of all the cycles containing that vertex. """ - chain = cython.declare(list) - cycleList = cython.declare(list) - - cycleList=list() - chain = [startingVertex] + return self.__exploreCyclesRecursively([startingVertex], []) - #chainLabels=range(len(self.keys())) - #print "Starting at %s in graph: %s"%(self.keys().index(startingVertex),chainLabels) - - cycleList = self.__exploreCyclesRecursively(chain, cycleList) - return cycleList - - def __exploreCyclesRecursively(self, chain, cycleList): + def __exploreCyclesRecursively(self, chain, cycles): """ - Finds cycles by spidering through a graph. - Give it a chain of atoms that are connected, `chain`, - and a list of cycles found so far `cycleList`. - If `chain` is a cycle, it is appended to `cycleList`. - Then chain is expanded by one atom (in each available direction) - and the function is called again. This recursively spiders outwards - from the starting chain, finding all the cycles. + Search the graph for cycles by recursive spidering. Given a `chain` + (list) of connected atoms and a list of `cycles` found so far, find any + cycles involving the chain of atoms and append them to the list of + cycles. This function recursively calls itself. """ - vertex2 = cython.declare(Vertex) - edge = cython.declare(Edge) - - # chainLabels = cython.declare(list) - # chainLabels=[self.keys().index(v) for v in chain] - # print "found %d so far. Chain=%s"%(len(cycleList),chainLabels) + cython.declare(vertex1=Vertex, vertex2=Vertex) - for vertex2, edge in self.edges[chain[-1]].iteritems(): - # vertex2 will loop through each of the atoms - # that are bonded to the last atom in the chain. + vertex1 = chain[-1] + # Loop over each of the atoms neighboring the last atom in the chain + for vertex2 in vertex1.edges: if vertex2 is chain[0] and len(chain) > 2: - # it is the first atom in the chain - so the chain IS a cycle! - cycleList.append(chain[:]) + # It is the first atom in the chain, so the chain is a cycle! + cycles.append(chain[:]) elif vertex2 not in chain: - # make the chain a little longer and explore again + # Make the chain a little longer and explore again chain.append(vertex2) - cycleList = self.__exploreCyclesRecursively(chain, cycleList) - # any cycles down this path (-vertex2) have now been found, - # so remove the vertex from the chain + cycles = self.__exploreCyclesRecursively(chain, cycles) + # Any cycles down this path have now been found, so remove vertex2 from the chain chain.pop(-1) - return cycleList + # At this point we should have discovered all of the cycles involving the current chain + return cycles def getSmallestSetOfSmallestRings(self): """ @@ -534,27 +530,21 @@ def getSmallestSetOfSmallestRings(self): from a Connection Table." *J. Chem. Inf. Comput. Sci.* **33**, p. 657-662 (1993). """ - - graph = cython.declare(Graph) - done = cython.declare(cython.bint) - verticesToRemove = cython.declare(list) - cycleList = cython.declare(list) - cycles = cython.declare(list) - vertex = cython.declare(Vertex) - rootVertex = cython.declare(Vertex) - found = cython.declare(cython.bint) - cycle = cython.declare(list) - graphs = cython.declare(list) + cython.declare(graph=Graph) + cython.declare(done=cython.bint, found=cython.bint) + cython.declare(cycleList=list, cycles=list, cycle=list, graphs=list, neighbors=list) + cython.declare(verticesToRemove=list) + cython.declare(vertex=Vertex, rootVertex=Vertex) # Make a copy of the graph so we don't modify the original - graph = self.copy() + graph = self.copy(deep=True) # Step 1: Remove all terminal vertices done = False while not done: verticesToRemove = [] - for vertex1, value in graph.edges.iteritems(): - if len(value) == 1: verticesToRemove.append(vertex1) + for vertex in graph.vertices: + if len(vertex.edges) == 1: verticesToRemove.append(vertex) done = len(verticesToRemove) == 0 # Remove identified vertices from graph for vertex in verticesToRemove: @@ -570,8 +560,6 @@ def getSmallestSetOfSmallestRings(self): for vertex in verticesToRemove: graph.removeVertex(vertex) - ### also need to remove EDGES that are not in ring - # Step 3: Split graph into remaining subgraphs graphs = graph.split() @@ -586,22 +574,15 @@ def getSmallestSetOfSmallestRings(self): for vertex in graph.vertices: if rootVertex is None: rootVertex = vertex - elif len(graph.edges[vertex]) < len(graph.edges[rootVertex]): + elif len(vertex.edges) < len(rootVertex.edges): rootVertex = vertex # Get all cycles involving the root vertex cycles = graph.getAllCycles(rootVertex) if len(cycles) == 0: - # this vertex is no longer in a ring. - # remove all its edges - neighbours = graph.edges[rootVertex].keys()[:] - for vertex2 in neighbours: - graph.removeEdge(rootVertex, vertex2) - # then remove it + # This vertex is no longer in a ring, so remove it graph.removeVertex(rootVertex) - #print("Removed vertex that's no longer in ring") - continue # (pick a new root Vertex) -# raise Exception('Did not find expected cycle!') + continue # Keep the smallest of the cycles found above cycle = cycles[0] @@ -613,13 +594,13 @@ def getSmallestSetOfSmallestRings(self): # Remove from the graph all vertices in the cycle that have only two edges verticesToRemove = [] for vertex in cycle: - if len(graph.edges[vertex]) <= 2: + if len(vertex.edges) <= 2: verticesToRemove.append(vertex) if len(verticesToRemove) == 0: # there are no vertices in this cycle that with only two edges - # Remove edge between root vertex and any one vertex it is connected to - graph.removeEdge(rootVertex, graph.edges[rootVertex].keys()[0]) + vertex = rootVertex.edges.keys()[0] + graph.removeEdge(rootVertex.edges[vertex]) else: for vertex in verticesToRemove: graph.removeVertex(vertex) @@ -632,17 +613,21 @@ def isMappingValid(self, other, mapping): is valid by checking that the vertices and edges involved in the mapping are mutually equivalent. """ - #cython.declare(vertex1=Vertex, vertex2=Vertex, vertices1=cython.list, vertices2=cython.list) - #cython.declare(i=cython.int, j=cython.int) + cython.declare(vertex1=Vertex, vertex2=Vertex) + cython.declare(vertices1=list, vertices2=list) + cython.declare(selfHasEdge=cython.bint, otherHasEdge=cython.bint) + cython.declare(i=cython.int, j=cython.int) + # Check that the mapped pairs of vertices are equivalent - for vertex1, vertex2 in mapping.iteritems(): + for vertex1, vertex2 in mapping.items(): if not vertex1.equivalent(vertex2): return False + # Check that any edges connected mapped vertices are equivalent vertices1 = mapping.keys() vertices2 = mapping.values() - for i in range(len(mapping)): - for j in range(i+1, len(mapping)): + for i in range(len(vertices1)): + for j in range(i+1, len(vertices1)): selfHasEdge = self.hasEdge(vertices1[i], vertices1[j]) otherHasEdge = other.hasEdge(vertices2[i], vertices2[j]) if selfHasEdge and otherHasEdge: @@ -654,6 +639,7 @@ def isMappingValid(self, other, mapping): elif selfHasEdge or otherHasEdge: # Only one of the graphs has the edge, so the mapping must be invalid return False + # If we're here then the vertices and edges are equivalent, so the # mapping is valid return True diff --git a/unittest/molecule/graphTest.py b/unittest/molecule/graphTest.py index 40817f7b6d..46e9326997 100644 --- a/unittest/molecule/graphTest.py +++ b/unittest/molecule/graphTest.py @@ -20,50 +20,48 @@ def setUp(self): A function run before each unit test in this class. """ vertices = [Vertex() for i in range(6)] - edges = [Edge() for i in range(5)] + edges = [ + Edge(vertices[0], vertices[1]), + Edge(vertices[1], vertices[2]), + Edge(vertices[2], vertices[3]), + Edge(vertices[3], vertices[4]), + Edge(vertices[4], vertices[5]), + ] self.graph = Graph() - for i in range(6): - self.graph.addVertex(vertices[i]) - self.graph.addEdge(vertices[0], vertices[1], edges[0]) - self.graph.addEdge(vertices[1], vertices[2], edges[1]) - self.graph.addEdge(vertices[2], vertices[3], edges[2]) - self.graph.addEdge(vertices[3], vertices[4], edges[3]) - self.graph.addEdge(vertices[4], vertices[5], edges[4]) + for vertex in vertices: self.graph.addVertex(vertex) + for edge in edges: self.graph.addEdge(edge) - def testAddVertex(self): + def test_addVertex(self): """ Test the Graph.addVertex() method. """ vertex = Vertex() self.graph.addVertex(vertex) self.assertTrue(vertex in self.graph.vertices) - self.assertTrue(vertex in self.graph.edges) - self.assertTrue(self.graph.edges[vertex] == {}) + self.assertTrue(vertex.edges == {}) - def testAddEdge(self): + def test_addEdge(self): """ Test the Graph.addEdge() method. """ - vertex1 = Vertex(); vertex2 = Vertex(); edge = Edge() + vertex1 = Vertex(); vertex2 = Vertex(); edge = Edge(vertex1, vertex2) try: - self.graph.addEdge(vertex1, vertex2, edge) + self.graph.addEdge(edge) self.fail('Added edge between vertices not in graph to graph.') - except KeyError: + except ValueError: pass self.graph.addVertex(vertex1) self.graph.addVertex(vertex2) - self.graph.addEdge(vertex1, vertex2, edge) + self.graph.addEdge(edge) self.assertTrue(vertex1 in self.graph.vertices) - self.assertTrue(vertex1 in self.graph.edges) + self.assertTrue(vertex1 in vertex2.edges) self.assertTrue(vertex2 in self.graph.vertices) - self.assertTrue(vertex2 in self.graph.edges) - self.assertTrue(vertex2 in self.graph.edges[vertex1]) - self.assertTrue(vertex1 in self.graph.edges[vertex2]) - self.assertTrue(self.graph.edges[vertex1][vertex2] is edge) - self.assertTrue(self.graph.edges[vertex2][vertex1] is edge) + self.assertTrue(vertex2 in vertex1.edges) + self.assertTrue(vertex1.edges[vertex2] is edge) + self.assertTrue(vertex2.edges[vertex1] is edge) - def testGetEdge(self): + def test_getEdge(self): """ Test the Graph.getEdge() method. """ @@ -72,17 +70,17 @@ def testGetEdge(self): try: edge = self.graph.getEdge(vertex1, vertex2) self.fail('Returned an edge between vertices that should not be connected in graph.') - except KeyError: + except ValueError: pass vertex1 = self.graph.vertices[2] vertex2 = self.graph.vertices[3] edge = self.graph.getEdge(vertex1, vertex2) self.assertNotEqual(edge, None) self.assertTrue(isinstance(edge, Edge)) - self.assertTrue(self.graph.edges[vertex1][vertex2] is edge) - self.assertTrue(self.graph.edges[vertex2][vertex1] is edge) + self.assertTrue(vertex1.edges[vertex2] is edge) + self.assertTrue(vertex2.edges[vertex1] is edge) - def testGetEdges(self): + def test_getEdges(self): """ Test the Graph.getEdges() method. """ @@ -93,7 +91,7 @@ def testGetEdges(self): self.assertTrue(self.graph.vertices[1] in edges) self.assertTrue(self.graph.vertices[3] in edges) - def testHasVertex(self): + def test_hasVertex(self): """ Test the Graph.hasVertex() method. """ @@ -102,7 +100,7 @@ def testHasVertex(self): for v in self.graph.vertices: self.assertTrue(self.graph.hasVertex(v)) - def testHasEdge(self): + def test_hasEdge(self): """ Test the Graph.hasEdge() method. """ @@ -113,7 +111,7 @@ def testHasEdge(self): vertex2 = self.graph.vertices[3] self.assertTrue(self.graph.hasEdge(vertex1, vertex2)) - def testRemoveVertex(self): + def test_removeVertex(self): """ Test the Graph.removeVertex() method. """ @@ -121,21 +119,114 @@ def testRemoveVertex(self): self.assertTrue(self.graph.hasVertex(vertex)) self.graph.removeVertex(vertex) self.assertFalse(self.graph.hasVertex(vertex)) - self.assertFalse(vertex in self.graph.edges) - for v in self.graph.edges: - self.assertFalse(vertex in self.graph.edges[v]) + for v in self.graph.vertices: + self.assertFalse(vertex in v.edges) - def testRemoveEdge(self): + def test_removeEdge(self): """ Test the Graph.removeEdge() method. """ vertex1 = self.graph.vertices[2] vertex2 = self.graph.vertices[3] self.assertTrue(self.graph.hasEdge(vertex1, vertex2)) - self.graph.removeEdge(vertex1, vertex2) - self.assertFalse(self.graph.hasEdge(vertex1, vertex2)) + edge = self.graph.getEdge(vertex1, vertex2) + self.graph.removeEdge(edge) + self.assertFalse(vertex1 in vertex2.edges) + self.assertFalse(vertex2 in vertex1.edges) + + def test_copy(self): + """ + Test the graph copy function to ensure a complete copy of the graph is + made while preserving vertices and edges. + """ + + vertices = [Vertex() for i in range(6)] + edges = [ + Edge(vertices[0], vertices[1]), + Edge(vertices[1], vertices[2]), + Edge(vertices[2], vertices[3]), + Edge(vertices[3], vertices[4]), + Edge(vertices[4], vertices[5]), + ] + + graph = Graph() + for vertex in vertices: graph.addVertex(vertex) + for edge in edges: graph.addEdge(edge) - def testResetConnectivityValues(self): + graph2 = graph.copy() + for vertex in graph.vertices: + self.assertTrue(graph2.hasVertex(vertex)) + for v1 in graph.vertices: + for v2 in v1.edges: + self.assertTrue(graph2.hasEdge(v1, v2)) + self.assertTrue(graph2.hasEdge(v2, v1)) + self.assertTrue(graph2.isIsomorphic(graph)) + self.assertTrue(graph.isIsomorphic(graph2)) + + def test_split(self): + """ + Test the graph split function to ensure a proper splitting of the graph + is being done. + """ + + vertices = [Vertex() for i in range(6)] + edges = [ + Edge(vertices[0], vertices[1]), + Edge(vertices[1], vertices[2]), + Edge(vertices[2], vertices[3]), + Edge(vertices[4], vertices[5]), + ] + + graph = Graph() + for vertex in vertices: graph.addVertex(vertex) + for edge in edges: graph.addEdge(edge) + + graphs = graph.split() + + self.assertTrue(len(graphs) == 2) + self.assertTrue(len(graphs[0].vertices) == 4 or len(graphs[0].vertices) == 2) + self.assertTrue(len(graphs[0].vertices) + len(graphs[1].vertices) == len(graph.vertices)) + + def test_merge(self): + """ + Test the graph merge function to ensure a proper merging of the graph + is being done. + """ + + vertices1 = [Vertex() for i in range(4)] + edges1 = [ + Edge(vertices1[0], vertices1[1]), + Edge(vertices1[1], vertices1[2]), + Edge(vertices1[2], vertices1[3]), + ] + + vertices2 = [Vertex() for i in range(3)] + edges2 = [ + Edge(vertices2[0], vertices2[1]), + Edge(vertices2[1], vertices2[2]), + ] + + graph1 = Graph() + for vertex in vertices1: graph1.addVertex(vertex) + for edge in edges1: graph1.addEdge(edge) + + graph2 = Graph() + for vertex in vertices2: graph2.addVertex(vertex) + for edge in edges2: graph2.addEdge(edge) + + graph = graph1.merge(graph2) + + self.assertTrue(len(graph1.vertices) + len(graph2.vertices) == len(graph.vertices)) + for vertex1 in vertices1: + self.assertTrue(vertex1 in graph.vertices) + for vertex2 in vertex1.edges: + self.assertTrue(vertex2 in graph.vertices) + for vertex2 in vertices2: + self.assertTrue(vertex2 in graph.vertices) + for vertex1 in vertex2.edges: + self.assertTrue(vertex1 in vertex2.edges) + + def test_resetConnectivityValues(self): """ Test the Graph.resetConnectivityValues() method. """ @@ -145,8 +236,8 @@ def testResetConnectivityValues(self): self.assertEqual(vertex.connectivity2, -1) self.assertEqual(vertex.connectivity3, -1) self.assertEqual(vertex.sortingLabel, -1) - - def testUpdateConnectivityValues(self): + + def test_updateConnectivityValues(self): """ Test the Graph.updateConnectivityValues() method. """ @@ -175,8 +266,8 @@ def testUpdateConnectivityValues(self): self.assertEqual(self.graph.vertices[5].connectivity2, 2) self.assertEqual(self.graph.vertices[5].connectivity3, 3) self.assertEqual(self.graph.vertices[5].sortingLabel, -1) - - def testSortVertices(self): + + def test_sortVertices(self): """ Test the Graph.sortVertices() method. """ @@ -187,36 +278,8 @@ def testSortVertices(self): self.assertTrue(vertex1.connectivity3 >= vertex2.connectivity3) self.assertTrue(vertex1.connectivity2 >= vertex2.connectivity2) self.assertTrue(vertex1.connectivity1 >= vertex2.connectivity1) - - def testCopy(self): - """ - Test the graph copy function to ensure a complete copy of the graph is - made while preserving vertices and edges. - """ - - vertices = [Vertex() for i in range(6)] - edges = [Edge() for i in range(5)] - - graph = Graph() - for vertex in vertices: graph.addVertex(vertex) - graph.addEdge(vertices[0], vertices[1], edges[0]) - graph.addEdge(vertices[1], vertices[2], edges[1]) - graph.addEdge(vertices[2], vertices[3], edges[2]) - graph.addEdge(vertices[3], vertices[4], edges[3]) - graph.addEdge(vertices[4], vertices[5], edges[4]) - - graph2 = graph.copy() - for vertex in graph.vertices: - self.assertTrue(vertex in graph2.edges) - self.assertTrue(graph2.hasVertex(vertex)) - for v1 in graph.vertices: - for v2 in graph.edges[v1]: - self.assertTrue(graph2.hasEdge(v1, v2)) - self.assertTrue(graph2.hasEdge(v2, v1)) - self.assertTrue(graph2.isIsomorphic(graph)) - self.assertTrue(graph.isIsomorphic(graph2)) - - def testVertexConnectivityValues(self): + + def test_vertex_connectivity_values(self): """ Tests the vertex connectivity values as introduced by Morgan (1965). @@ -232,248 +295,198 @@ def testVertexConnectivityValues(self): """ vertices = [Vertex() for i in range(6)] - edges = [Edge() for i in range(5)] - + edges = [ + Edge(vertices[0], vertices[1]), + Edge(vertices[1], vertices[2]), + Edge(vertices[2], vertices[3]), + Edge(vertices[3], vertices[4]), + Edge(vertices[1], vertices[5]), + ] + graph = Graph() for vertex in vertices: graph.addVertex(vertex) - graph.addEdge(vertices[0], vertices[1], edges[0]) - graph.addEdge(vertices[1], vertices[2], edges[1]) - graph.addEdge(vertices[2], vertices[3], edges[2]) - graph.addEdge(vertices[3], vertices[4], edges[3]) - graph.addEdge(vertices[1], vertices[5], edges[4]) + for edge in edges: graph.addEdge(edge) graph.updateConnectivityValues() for i,cv_ in enumerate([1,3,2,2,1,1]): cv = vertices[i].connectivity1 - self.assertEqual(cv, cv_, "On vertex %d got connectivity[0]=%d but expected %d"%(i,cv,cv_)) + self.assertEqual(cv, cv_, "On vertex {0:d} got connectivity[0]={1:d} but expected {2:d}".format(i,cv,cv_)) for i,cv_ in enumerate([3,4,5,3,2,3]): cv = vertices[i].connectivity2 - self.assertEqual(cv, cv_, "On vertex %d got connectivity[1]=%d but expected %d"%(i,cv,cv_)) + self.assertEqual(cv, cv_, "On vertex {0:d} got connectivity[0]={1:d} but expected {2:d}".format(i,cv,cv_)) for i,cv_ in enumerate([4,11,7,7,3,4]): cv = vertices[i].connectivity3 - self.assertEqual(cv, cv_, "On vertex %d got connectivity[2]=%d but expected %d"%(i,cv,cv_)) + self.assertEqual(cv, cv_, "On vertex {0:d} got connectivity[0]={1:d} but expected {2:d}".format(i,cv,cv_)) - def testSplit(self): - """ - Test the graph split function to ensure a proper splitting of the graph - is being done. - """ - - vertices = [Vertex() for i in range(6)] - edges = [Edge() for i in range(4)] - - graph = Graph() - for vertex in vertices: graph.addVertex(vertex) - graph.addEdge(vertices[0], vertices[1], edges[0]) - graph.addEdge(vertices[1], vertices[2], edges[1]) - graph.addEdge(vertices[2], vertices[3], edges[2]) - graph.addEdge(vertices[4], vertices[5], edges[3]) - - graphs = graph.split() - - self.assertTrue(len(graphs) == 2) - self.assertTrue(len(graphs[0].vertices) == 4 or len(graphs[0].vertices) == 2) - self.assertTrue(len(graphs[0].vertices) + len(graphs[1].vertices) == len(graph.vertices)) - - def testMerge(self): - """ - Test the graph merge function to ensure a proper merging of the graph - is being done. - """ - - vertices1 = [Vertex() for i in range(4)] - edges1 = [Edge() for i in range(3)] - - vertices2 = [Vertex() for i in range(3)] - edges2 = [Edge() for i in range(2)] - - graph1 = Graph() - for vertex in vertices1: graph1.addVertex(vertex) - graph1.addEdge(vertices1[0], vertices1[1], edges1[0]) - graph1.addEdge(vertices1[1], vertices1[2], edges1[1]) - graph1.addEdge(vertices1[2], vertices1[3], edges1[2]) - - graph2 = Graph() - for vertex in vertices2: graph2.addVertex(vertex) - graph2.addEdge(vertices2[0], vertices2[1], edges2[0]) - graph2.addEdge(vertices2[1], vertices2[2], edges2[1]) - - graph = graph1.merge(graph2) - - self.assertTrue(len(graph1.vertices) + len(graph2.vertices) == len(graph.vertices)) - for vertex1 in vertices1: - self.assertTrue(vertex1 in graph.edges) - self.assertTrue(len(graph1.edges[vertex1]) == len(graph.edges[vertex1])) - for vertex2 in graph1.edges[vertex1]: - self.assertTrue(vertex2 in graph.edges[vertex1]) - self.assertTrue(graph1.edges[vertex1][vertex2] is graph.edges[vertex1][vertex2]) - for vertex2 in vertices2: - self.assertTrue(vertex2 in graph.edges) - self.assertTrue(len(graph2.edges[vertex2]) == len(graph.edges[vertex2])) - for vertex1 in graph2.edges[vertex2]: - self.assertTrue(vertex1 in graph.edges[vertex2]) - self.assertTrue(graph2.edges[vertex2][vertex1] is graph.edges[vertex2][vertex1]) - - def testIsomorphism(self): + def test_isomorphism(self): """ Check the graph isomorphism functions. """ vertices1 = [Vertex() for i in range(6)] - edges1 = [Edge() for i in range(5)] + edges1 = [ + Edge(vertices1[0], vertices1[1]), + Edge(vertices1[1], vertices1[2]), + Edge(vertices1[2], vertices1[3]), + Edge(vertices1[3], vertices1[4]), + Edge(vertices1[4], vertices1[5]), + ] + vertices2 = [Vertex() for i in range(6)] - edges2 = [Edge() for i in range(5)] + edges2 = [ + Edge(vertices2[0], vertices2[1]), + Edge(vertices2[1], vertices2[2]), + Edge(vertices2[2], vertices2[3]), + Edge(vertices2[3], vertices2[4]), + Edge(vertices2[4], vertices2[5]), + ] graph1 = Graph() for vertex in vertices1: graph1.addVertex(vertex) - graph1.edges[vertices1[0]] = { vertices1[1]: edges1[0] } - graph1.edges[vertices1[1]] = { vertices1[0]: edges1[0], vertices1[2]: edges1[1] } - graph1.edges[vertices1[2]] = { vertices1[1]: edges1[1], vertices1[3]: edges1[2] } - graph1.edges[vertices1[3]] = { vertices1[2]: edges1[2], vertices1[4]: edges1[3] } - graph1.edges[vertices1[4]] = { vertices1[3]: edges1[3], vertices1[5]: edges1[4] } - graph1.edges[vertices1[5]] = { vertices1[4]: edges1[4] } + for edge in edges1: graph1.addEdge(edge) graph2 = Graph() for vertex in vertices2: graph2.addVertex(vertex) - graph2.edges[vertices2[0]] = { vertices2[1]: edges2[4] } - graph2.edges[vertices2[1]] = { vertices2[0]: edges2[4], vertices2[2]: edges2[3] } - graph2.edges[vertices2[2]] = { vertices2[1]: edges2[3], vertices2[3]: edges2[2] } - graph2.edges[vertices2[3]] = { vertices2[2]: edges2[2], vertices2[4]: edges2[1] } - graph2.edges[vertices2[4]] = { vertices2[3]: edges2[1], vertices2[5]: edges2[0] } - graph2.edges[vertices2[5]] = { vertices2[4]: edges2[0] } + for edge in edges2: graph2.addEdge(edge) self.assertTrue(graph1.isIsomorphic(graph2)) self.assertTrue(graph1.isSubgraphIsomorphic(graph2)) self.assertTrue(graph2.isIsomorphic(graph1)) self.assertTrue(graph2.isSubgraphIsomorphic(graph1)) - def testSubgraphIsomorphism(self): + def test_subgraphIsomorphism(self): """ Check the subgraph isomorphism functions. """ vertices1 = [Vertex() for i in range(6)] - edges1 = [Edge() for i in range(5)] + edges1 = [ + Edge(vertices1[0], vertices1[1]), + Edge(vertices1[1], vertices1[2]), + Edge(vertices1[2], vertices1[3]), + Edge(vertices1[3], vertices1[4]), + Edge(vertices1[4], vertices1[5]), + ] vertices2 = [Vertex() for i in range(2)] - edges2 = [Edge() for i in range(1)] + edges2 = [ + Edge(vertices2[0], vertices2[1]), + ] graph1 = Graph() for vertex in vertices1: graph1.addVertex(vertex) - graph1.edges[vertices1[0]] = { vertices1[1]: edges1[0] } - graph1.edges[vertices1[1]] = { vertices1[0]: edges1[0], vertices1[2]: edges1[1] } - graph1.edges[vertices1[2]] = { vertices1[1]: edges1[1], vertices1[3]: edges1[2] } - graph1.edges[vertices1[3]] = { vertices1[2]: edges1[2], vertices1[4]: edges1[3] } - graph1.edges[vertices1[4]] = { vertices1[3]: edges1[3], vertices1[5]: edges1[4] } - graph1.edges[vertices1[5]] = { vertices1[4]: edges1[4] } + for edge in edges1: graph1.addEdge(edge) graph2 = Graph() for vertex in vertices2: graph2.addVertex(vertex) - graph2.edges[vertices2[0]] = { vertices2[1]: edges2[0] } - graph2.edges[vertices2[1]] = { vertices2[0]: edges2[0] } + for edge in edges2: graph2.addEdge(edge) self.assertFalse(graph1.isIsomorphic(graph2)) self.assertFalse(graph2.isIsomorphic(graph1)) self.assertTrue(graph1.isSubgraphIsomorphic(graph2)) - ismatch, mapList = graph1.findSubgraphIsomorphisms(graph2) - self.assertTrue(ismatch) + mapList = graph1.findSubgraphIsomorphisms(graph2) self.assertTrue(len(mapList) == 10) for mapping in mapList: self.assertTrue( graph1.isMappingValid(graph2,mapping) ) self.assertTrue( graph1.isMappingValid(graph2,mapping) ) + + def test_pickle(self): + """ + Test that a Graph object can be successfully pickled and unpickled + with no loss of information. + """ - def testIsCyclic(self): + vertices = [Vertex() for i in range(6)] + edges = [ + Edge(vertices[0], vertices[1]), + Edge(vertices[1], vertices[2]), + Edge(vertices[2], vertices[3]), + Edge(vertices[3], vertices[4]), + Edge(vertices[4], vertices[5]), + ] + + graph0 = Graph() + for vertex in vertices: graph0.addVertex(vertex) + for edge in edges: graph0.addEdge(edge) + graph0.updateConnectivityValues() + + import cPickle + graph = cPickle.loads(cPickle.dumps(graph0)) + + self.assertEqual(len(graph0.vertices), len(graph.vertices)) + for v1, v2 in zip(graph0.vertices, graph.vertices): + self.assertEqual(v1.connectivity1, v2.connectivity1) + self.assertEqual(v1.connectivity2, v2.connectivity2) + self.assertEqual(v1.connectivity3, v2.connectivity3) + self.assertEqual(v1.sortingLabel, v2.sortingLabel) + self.assertEqual(len(v1.edges), len(v2.edges)) + self.assertTrue(graph0.isIsomorphic(graph)) + self.assertTrue(graph.isIsomorphic(graph0)) + + def test_isCyclic(self): """ Test the Graph.isCyclic() method. """ self.assertFalse(self.graph.isCyclic()) - self.graph.addEdge(self.graph.vertices[0], self.graph.vertices[3], Edge()) # To create a cycle + edge = Edge(self.graph.vertices[0], self.graph.vertices[3]) + self.graph.addEdge(edge) # To create a cycle self.assertTrue(self.graph.isCyclic()) - def testIsVertexInCycle(self): + def test_isVertexInCycle(self): """ Test the Graph.isVertexInCycle() method. """ for vertex in self.graph.vertices: self.assertFalse(self.graph.isVertexInCycle(vertex)) - self.graph.addEdge(self.graph.vertices[0], self.graph.vertices[3], Edge()) # To create a cycle + edge = Edge(self.graph.vertices[0], self.graph.vertices[3]) + self.graph.addEdge(edge) # To create a cycle for vertex in self.graph.vertices[0:4]: self.assertTrue(self.graph.isVertexInCycle(vertex)) for vertex in self.graph.vertices[4:]: self.assertFalse(self.graph.isVertexInCycle(vertex)) - def testIsEdgeInCycle(self): + def test_isEdgeInCycle(self): """ Test the Graph.isEdgeInCycle() method. """ - for vertex1 in self.graph.edges: - for vertex2 in self.graph.edges[vertex1]: - self.assertFalse(self.graph.isEdgeInCycle(vertex1, vertex2)) - self.graph.addEdge(self.graph.vertices[0], self.graph.vertices[3], Edge()) # To create a cycle - for vertex1 in self.graph.edges: - for vertex2 in self.graph.edges[vertex1]: + for vertex1 in self.graph.vertices: + for vertex2, edge in vertex1.edges.items(): + self.assertFalse(self.graph.isEdgeInCycle(edge)) + edge = Edge(self.graph.vertices[0], self.graph.vertices[3]) + self.graph.addEdge(edge) # To create a cycle + for vertex1 in self.graph.vertices: + for vertex2, edge in vertex1.edges.items(): if self.graph.vertices.index(vertex1) < 4 and self.graph.vertices.index(vertex2) < 4: - self.assertTrue(self.graph.isEdgeInCycle(vertex1, vertex2)) + self.assertTrue(self.graph.isEdgeInCycle(edge)) else: - self.assertFalse(self.graph.isEdgeInCycle(vertex1, vertex2)) + self.assertFalse(self.graph.isEdgeInCycle(edge)) - def testGetAllCycles(self): + def test_getAllCycles(self): """ Test the Graph.getAllCycles() method. """ cycleList = self.graph.getAllCycles(self.graph.vertices[0]) self.assertEqual(len(cycleList), 0) - self.graph.addEdge(self.graph.vertices[0], self.graph.vertices[3], Edge()) # To create a cycle + edge = Edge(self.graph.vertices[0], self.graph.vertices[3]) + self.graph.addEdge(edge) # To create a cycle cycleList = self.graph.getAllCycles(self.graph.vertices[0]) self.assertEqual(len(cycleList), 2) self.assertEqual(len(cycleList[0]), 4) self.assertEqual(len(cycleList[1]), 4) - def testGetSmallestSetOfSmallestRings(self): + def test_getSmallestSetOfSmallestRings(self): """ Test the Graph.getSmallestSetOfSmallestRings() method. """ cycleList = self.graph.getSmallestSetOfSmallestRings() self.assertEqual(len(cycleList), 0) - self.graph.addEdge(self.graph.vertices[0], self.graph.vertices[3], Edge()) # To create a cycle + edge = Edge(self.graph.vertices[0], self.graph.vertices[3]) + self.graph.addEdge(edge) # To create a cycle cycleList = self.graph.getSmallestSetOfSmallestRings() self.assertEqual(len(cycleList), 1) self.assertEqual(len(cycleList[0]), 4) - - def testPickle(self): - """ - Test that a Graph object can be successfully pickled and unpickled - with no loss of information. - """ - - vertices = [Vertex() for i in range(6)] - edges = [Edge() for i in range(5)] - - graph0 = Graph() - for vertex in vertices: graph0.addVertex(vertex) - graph0.addEdge(vertices[0], vertices[1], edges[0]) - graph0.addEdge(vertices[1], vertices[2], edges[1]) - graph0.addEdge(vertices[2], vertices[3], edges[2]) - graph0.addEdge(vertices[3], vertices[4], edges[3]) - graph0.addEdge(vertices[4], vertices[5], edges[4]) - graph0.updateConnectivityValues() - - import cPickle - graph = cPickle.loads(cPickle.dumps(graph0)) - - self.assertEqual(len(graph0.vertices), len(graph.vertices)) - self.assertEqual(len(graph0.edges), len(graph.edges)) - for v1, v2 in zip(graph0.vertices, graph.vertices): - self.assertEqual(v1.connectivity1, v2.connectivity1) - self.assertEqual(v1.connectivity2, v2.connectivity2) - self.assertEqual(v1.connectivity3, v2.connectivity3) - self.assertEqual(v1.sortingLabel, v2.sortingLabel) - self.assertEqual(len(graph0.edges[v1]), len(graph.edges[v2])) - self.assertTrue(graph0.isIsomorphic(graph)) - self.assertTrue(graph.isIsomorphic(graph0)) - ################################################################################ From 4f103d9f050376861093de51ea2db555e26b526b Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 16:14:49 -0400 Subject: [PATCH 11/23] Added new implementation of VF2 algorithm for graph isomorphism. The new implementation exploits the fact that we are now storing the edges on the Vertex objects to dramatically speed up the isomorphism check by significantly decreasing the number of calls into the slow Python/C API. In earlier tests I was seeing a ~8x speedup in isomorphism evaluation for Graph objects, and about ~2x for Molecule objects. --- rmgpy/molecule/graph.pxd | 43 ++-- rmgpy/molecule/graph.py | 484 +++++++++++++++++++-------------------- 2 files changed, 252 insertions(+), 275 deletions(-) diff --git a/rmgpy/molecule/graph.pxd b/rmgpy/molecule/graph.pxd index e85721dec1..b90fe675b5 100644 --- a/rmgpy/molecule/graph.pxd +++ b/rmgpy/molecule/graph.pxd @@ -28,14 +28,19 @@ cdef class Vertex(object): cdef readonly dict edges + # These attributes are used in the VF2 graph isomorphism algorithm cdef public short connectivity1 cdef public short connectivity2 cdef public short connectivity3 cdef public short sortingLabel + cdef public bint terminal + cdef public Vertex mapping - cpdef bint equivalent(self, Vertex other) + cpdef Vertex copy(self) - cpdef bint isSpecificCaseOf(self, Vertex other) + cpdef bint equivalent(self, Vertex other) except -2 + + cpdef bint isSpecificCaseOf(self, Vertex other) except -2 cpdef resetConnectivityValues(self) @@ -91,21 +96,21 @@ cdef class Graph: cpdef sortVertices(self) - cpdef bint isIsomorphic(self, Graph other, dict initialMap=?) + cpdef bint isIsomorphic(self, Graph other, dict initialMap=?) except -2 - cpdef tuple findIsomorphism(self, Graph other, dict initialMap=?) + cpdef list findIsomorphism(self, Graph other, dict initialMap=?) - cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?) + cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?) except -2 - cpdef tuple findSubgraphIsomorphisms(self, Graph other, dict initialMap=?) + cpdef list findSubgraphIsomorphisms(self, Graph other, dict initialMap=?) - cpdef bint isCyclic(self) + cpdef bint isCyclic(self) except -2 - cpdef bint isVertexInCycle(self, Vertex vertex) + cpdef bint isVertexInCycle(self, Vertex vertex) except -2 - cpdef bint isEdgeInCycle(self, Vertex vertex1, Vertex vertex2) + cpdef bint isEdgeInCycle(self, Edge edge) except -2 - cpdef bint __isChainInCycle(self, list chain) + cpdef bint __isChainInCycle(self, list chain) except -2 cpdef getAllCycles(self, Vertex startingVertex) @@ -113,22 +118,16 @@ cdef class Graph: cpdef getSmallestSetOfSmallestRings(self) - cpdef bint isMappingValid(self, Graph other, dict mapping) + cpdef bint isMappingValid(self, Graph other, dict mapping) except -2 ################################################################################ -cpdef VF2_isomorphism(Graph graph1, Graph graph2, bint subgraph=?, - bint findAll=?, dict initialMap=?) +cdef VF2_isomorphism(Graph graph1, Graph graph2, bint subgraph=?, bint findAll=?, dict initialMapping=?) -cpdef bint __VF2_feasible(Graph graph1, Graph graph2, Vertex vertex1, - Vertex vertex2, dict map21, dict map12, list terminals1, list terminals2, - bint subgraph) except -2 # bint should be 0 or 1 +cpdef bint VF2_feasible(Graph graph1, Graph graph2, Vertex vertex1, Vertex vertex2, bint subgraph) except -2 -cpdef bint __VF2_match(Graph graph1, Graph graph2, dict map21, dict map12, - list terminals1, list terminals2, bint subgraph, bint findAll, - list map21List, list map12List, int call_depth) except -2 # bint should be 0 or 1 +cpdef bint VF2_match(Graph graph1, Graph graph2, bint subgraph, bint findAll, list mappingList, int callDepth) except -2 -cpdef list __VF2_terminals(Graph graph, dict mapping) +cpdef VF2_addToMapping(Vertex vertex1, Vertex vertex2) -cpdef list __VF2_updateTerminals(Graph graph, dict mapping, list old_terminals, - Vertex new_vertex) +cpdef VF2_removeFromMapping(Vertex vertex1, Vertex vertex2) diff --git a/rmgpy/molecule/graph.py b/rmgpy/molecule/graph.py index ac02052bb1..b46d56a3da 100644 --- a/rmgpy/molecule/graph.py +++ b/rmgpy/molecule/graph.py @@ -71,6 +71,8 @@ def __reduce__(self): 'connectivity2': self.connectivity2, 'connectivity3': self.connectivity3, 'sortingLabel': self.sortingLabel, + 'terminal': self.terminal, + 'mapping': self.mapping, } return (Vertex, (), d) @@ -80,6 +82,17 @@ def __setstate__(self, d): self.connectivity2 = d['connectivity2'] self.connectivity3 = d['connectivity3'] self.sortingLabel = d['sortingLabel'] + self.terminal = d['terminal'] + self.mapping = d['mapping'] + + def copy(self): + """ + Return a copy of the vertex. The default implementation assumes that no + semantic information is associated with each vertex, and therefore + simply returns a new :class:`Vertex` object. + """ + new = Vertex() + return new def equivalent(self, other): """ @@ -105,6 +118,8 @@ def resetConnectivityValues(self): self.connectivity2 = -1 self.connectivity3 = -1 self.sortingLabel = -1 + self.terminal = False + self.mapping = None def getVertexConnectivityValue(vertex): """ @@ -141,6 +156,16 @@ def __reduce__(self): """ return (Edge, (self.vertex1, self.vertex2)) + def copy(self): + """ + Return a copy of the edge. The default implementation assumes that no + semantic information is associated with each edge, and therefore + simply returns a new :class:`Edge` object. Note that the vertices are + not copied in this implementation. + """ + new = Edge(self.vertex1, self.vertex2) + return new + def equivalent(self, other): """ Return ``True`` if two edges `self` and `other` are semantically @@ -407,7 +432,7 @@ def isIsomorphic(self, other, initialMap=None): Returns :data:`True` if two graphs are isomorphic and :data:`False` otherwise. Uses the VF2 algorithm of Vento and Foggia. """ - return VF2_isomorphism(self, other, False, False, initialMap)[0] + return VF2_isomorphism(self, other, False, False, initialMap) def findIsomorphism(self, other, initialMap=None): """ @@ -422,7 +447,7 @@ def isSubgraphIsomorphic(self, other, initialMap=None): Returns :data:`True` if `other` is subgraph isomorphic and :data:`False` otherwise. Uses the VF2 algorithm of Vento and Foggia. """ - return VF2_isomorphism(self, other, True, False, initialMap)[0] + return VF2_isomorphism(self, other, True, False, initialMap) def findSubgraphIsomorphisms(self, other, initialMap=None): """ @@ -646,106 +671,95 @@ def isMappingValid(self, other, mapping): ################################################################################ -def VF2_isomorphism(graph1, graph2, subgraph=False, findAll=False, initialMap=None): +class VF2Error(Exception): """ - Determines if two :class:`Graph` objects `graph1` and `graph2` are - isomorphic. A number of options affect how the isomorphism check is - performed: - - * If `subgraph` is ``True``, the isomorphism function will treat `graph2` - as a subgraph of `graph1`. In this instance a subgraph can either mean a - smaller graph (i.e. fewer vertices and/or edges) or a less specific graph. - - * If `findAll` is ``True``, all valid isomorphisms will be found and - returned; otherwise only the first valid isomorphism will be returned. - - * The `initialMap` parameter can be used to pass a previously-established - mapping. This mapping will be preserved in all returned valid - isomorphisms. - - The isomorphism algorithm used is the VF2 algorithm of Vento and Foggia. - The function returns a boolean `isMatch` indicating whether or not one or - more valid isomorphisms have been found, and a list `mapList` of the valid - isomorphisms, each consisting of a dictionary mapping from vertices of - `graph1` to corresponding vertices of `graph2`. + An exception raised if an error occurs within the VF2 graph isomorphism + algorithm. Pass a string describing the error. """ + pass - cython.declare(isMatch=cython.bint, map12List=list, map21List=list) - cython.declare(terminals1=list, terminals2=list, callDepth=cython.int) - cython.declare(vert=Vertex) - - map21List = list() - - # Some quick initial checks to avoid using the full algorithm if the - # graphs are obviously not isomorphic (based on graph size) - if not subgraph: - if len(graph2.vertices) != len(graph1.vertices): - # The two graphs don't have the same number of vertices, so they - # cannot be isomorphic - return False, map21List - elif len(graph1.vertices) == len(graph2.vertices) == 0: - logging.warning("Tried matching empty graphs (returning True)") - # The two graphs don't have any vertices; this means they are - # trivially isomorphic - return True, map21List - else: - if len(graph2.vertices) > len(graph1.vertices): - # The second graph has more vertices than the first, so it cannot be - # a subgraph of the first - return False, map21List - - if initialMap is None: initialMap = {} - map12List = list() +def VF2_isomorphism(graph1, graph2, subgraph=False, findAll=False, initialMapping=None): + """ + Use the VF2 algorithm of Vento and Foggia to evaluate the isomorphism of + the graphs `graph1` and `graph2`. A number of options affect how the + isomorphism check is performed: + + * If `subgraph` is ``True``, the function will determine if `graph2` is a + subgraph of `graph1` instead of a full graph. + * If `findAll` is ``True``, this function returns a list of valid mappings + from `graph1` to `graph2`; each mapping is a ``dict`` with vertices from + `graph1` as the keys and vertices from `graph2` as the values. If + `findAll` is ``False``, this function simply returns ``True`` if a + valid mapping was found, or ``False`` if not. + + * The `initialMapping` parameter is used to specify a mapping of vertices + from `graph1` to those in `subgraph` that are fixed in the isomorphism + check; this mapping will appear in every returned mapping. Note that no + validation of this initial mapping is performed in this function. + + """ + cython.declare(vertex1=Vertex, vertex2=Vertex) + cython.declare(mappingList=list) + cython.declare(callDepth=cython.int) + + # Some quick isomorphism checks based on graph sizes + if not subgraph and len(graph2.vertices) != len(graph1.vertices): + # The two graphs don't have the same number of vertices, so they + # cannot be isomorphic + return list() if findAll else False + elif not subgraph and len(graph2.vertices) == len(graph1.vertices) == 0: + # The two graphs don't have any vertices; this means they are + # trivially isomorphic + return list() if findAll else True + elif subgraph and len(graph2.vertices) > len(graph1.vertices): + # The second graph has more vertices than the first, so it cannot be + # a subgraph of the first + return list() if findAll else False + # Initialize callDepth with the size of the largest graph - # Each recursive call to __VF2_match will decrease it by one; + # Each recursive call to VF2_match will decrease it by one; # when the whole graph has been explored, it should reach 0 # It should never go below zero! - callDepth = min(len(graph1.vertices), len(graph2.vertices)) - len(initialMap) + callDepth = len(graph2.vertices) # Sort the vertices in each graph to make the isomorphism more efficient graph1.sortVertices() graph2.sortVertices() - # Generate initial mapping pairs - # map21 = map to 2 from 1 - # map12 = map to 1 from 2 - map21 = initialMap - map12 = dict([(v,k) for k,v in initialMap.iteritems()]) - - # Generate an initial set of terminals - terminals1 = __VF2_terminals(graph1, map21) - terminals2 = __VF2_terminals(graph2, map12) - - isMatch = __VF2_match(graph1, graph2, map21, map12, \ - terminals1, terminals2, subgraph, findAll, map21List, map12List, callDepth) - - if findAll: - return len(map21List) > 0, map21List - else: - return isMatch, map21 - -def __VF2_feasible(graph1, graph2, vertex1, vertex2, map21, map12, terminals1, - terminals2, subgraph): + # Initialize mapping by clearing any previous mapping information + for vertex1 in graph1.vertices: + vertex1.mapping = None + vertex1.terminal = False + for vertex2 in graph2.vertices: + vertex2.mapping = None + vertex2.terminal = False + # Set the initial mapping if provided + if initialMapping is not None: + for vertex1, vertex2 in initialMapping.items(): + VF2_addToMapping(vertex1, vertex2) + callDepth -= len(initialMapping) + + mappingList = [] + isMatch = VF2_match(graph1, graph2, subgraph, findAll, mappingList, callDepth) + + return mappingList if findAll else isMatch + +def VF2_feasible(graph1, graph2, vertex1, vertex2, subgraph): """ - Returns :data:`True` if two vertices `vertex1` and `vertex2` from graphs - `graph1` and `graph2`, respectively, are feasible matches. `mapping21` and - `mapping12` are the current state of the mapping from `graph1` to `graph2` - and vice versa, respectively. `terminals1` and `terminals2` are lists of - the vertices that are directly connected to the already-mapped vertices. - `subgraph` is :data:`True` if graph2 is to be treated as a potential - subgraph of graph1. i.e. graph1 is a specific case of graph2. - - Uses the VF2 algorithm of Vento and Foggia. The feasibility is assessed - through a series of semantic and structural checks. Only the combination - of the semantic checks and the level 0 structural check are both - necessary and sufficient to ensure feasibility. (This does *not* mean that - vertex1 and vertex2 are always a match, although the level 1 and level 2 - checks preemptively eliminate a number of false positives.) + Return ``True`` if two vertices `vertex1` and `vertex2` from graphs + `graph1` and `graph2`, respectively, are feasible matches. The `subgraph` + flag indicates whether or not to treat `graph2` as a subgraph of `graph1`. + + The feasibility is assessed through a series of semantic and structural + checks. Only the combination of the semantic checks and the level 0 + structural check are both necessary and sufficient to ensure feasibility. + (This does *not* mean that `vertex1` and `vertex2` are always a match, + although the level 1 and level 2 checks preemptively eliminate a number of + false positives.) """ - - cython.declare(vert1=Vertex, vert2=Vertex, edge1=Edge, edge2=Edge, edges1=dict, edges2=dict) - cython.declare(i=cython.int) + cython.declare(vert1=Vertex, vert2=Vertex) + cython.declare(edge1=Edge, edge2=Edge) cython.declare(term1Count=cython.int, term2Count=cython.int, neither1Count=cython.int, neither2Count=cython.int) if not subgraph: @@ -753,48 +767,45 @@ def __VF2_feasible(graph1, graph2, vertex1, vertex2, map21, map12, terminals1, if vertex1.connectivity1 != vertex2.connectivity1: return False if vertex1.connectivity2 != vertex2.connectivity2: return False if vertex1.connectivity3 != vertex2.connectivity3: return False - + # Semantic check #1: vertex1 and vertex2 must be equivalent if subgraph: if not vertex1.isSpecificCaseOf(vertex2): return False else: if not vertex1.equivalent(vertex2): return False - - # Get edges adjacent to each vertex - edges1 = graph1.edges[vertex1] - edges2 = graph2.edges[vertex2] - + # Semantic check #2: adjacent vertices to vertex1 and vertex2 that are # already mapped should be connected by equivalent edges - for vert2 in edges2: - if vert2 in map12: - vert1 = map12[vert2] - if not vert1 in edges1: # atoms not joined in graph1 + for vert2 in vertex2.edges: + if vert2.mapping is not None: + vert1 = vert2.mapping + if vert1 not in vertex1.edges: + # The vertices are joined in graph2, but not in graph1 return False - edge1 = edges1[vert1] - edge2 = edges2[vert2] + edge1 = vertex1.edges[vert1] + edge2 = vertex2.edges[vert2] if subgraph: if not edge1.isSpecificCaseOf(edge2): return False - else: # exact match required + else: if not edge1.equivalent(edge2): return False - # there could still be edges in graph1 that aren't in graph2. - # this is ok for subgraph matching, but not for exact matching + # There could still be edges in graph1 that aren't in graph2; this is okay + # for subgraph matching, but not for exact matching if not subgraph: - for vert1 in edges1: - if vert1 in map21: - vert2 = map21[vert1] - if not vert2 in edges2: return False + for vert1 in vertex1.edges: + if vert1.mapping is not None: + if vert1.mapping not in vertex2.edges: + # The vertices are joined in graph1, but not in graph2 + return False # Count number of terminals adjacent to vertex1 and vertex2 term1Count = 0; term2Count = 0; neither1Count = 0; neither2Count = 0 - - for vert1 in edges1: - if vert1 in terminals1: term1Count += 1 - elif vert1 not in map21: neither1Count += 1 - for vert2 in edges2: - if vert2 in terminals2: term2Count += 1 - elif vert2 not in map12: neither2Count += 1 + for vert1 in vertex1.edges: + if vert1.terminal: term1Count += 1 + elif vert1.mapping is not None: neither1Count += 1 + for vert2 in vertex2.edges: + if vert2.terminal: term2Count += 1 + elif vert2.mapping is not None: neither2Count += 1 # Level 2 look-ahead: the number of adjacent vertices of vertex1 and # vertex2 that are non-terminals must be equal @@ -812,172 +823,139 @@ def __VF2_feasible(graph1, graph2, vertex1, vertex2, map21, map12, terminals1, # Level 0 look-ahead: all adjacent vertices of vertex2 already in the # mapping must map to adjacent vertices of vertex1 - for vert2 in edges2: - if vert2 in map12: - vert1 = map12[vert2] - if vert1 not in edges1: return False + for vert2 in vertex2.edges: + if vert2.mapping is not None: + if vert2.mapping not in vertex1.edges: return False # Also, all adjacent vertices of vertex1 already in the mapping must map to # adjacent vertices of vertex2, unless we are subgraph matching if not subgraph: - for vert1 in edges1: - if vert1 in map21: - vert2 = map21[vert1] - if vert2 not in edges2: return False + for vert1 in vertex1.edges: + if vert1.mapping is not None: + if vert1.mapping not in vertex2.edges: return False - # All of our tests have been passed, so the two vertices are a feasible - # pair + # All of our tests have been passed, so the two vertices are a feasible pair return True -def __VF2_match(graph1, graph2, map21, map12, terminals1, terminals2, subgraph, - findAll, map21List, map12List, callDepth): +def VF2_match(graph1, graph2, subgraph, findAll, mappingList, callDepth): """ - A recursive function used to explore two graphs `graph1` and `graph2` for - isomorphism by attempting to map them to one another. `mapping21` and - `mapping12` are the current state of the mapping from `graph1` to `graph2` - and vice versa, respectively. `terminals1` and `terminals2` are lists of - the vertices that are directly connected to the already-mapped vertices. - `subgraph` is :data:`True` if graph2 is to be treated as a potential - subgraph of graph1. i.e. graph1 is a specific case of graph2. - - If findAll=True then it adds valid mappings to map21List and - map12List, but returns False when done (or True if the initial mapping is complete) - - Uses the VF2 algorithm of Vento and Foggia, which is O(N) in spatial complexity - and O(N**2) (best-case) to O(N! * N) (worst-case) in temporal complexity. + Recursively explore two graphs `graph1` and `graph2` in search of one or + more isomorphism relationships by attempting to map vertices to one + another. The `subgraph` flag indicates whether or not to treat `graph2` as + a subgraph of `graph1`. The `findAll` flag indicates whether to find all + isomorphisms or only the first. The `mappingList` parameter stores the + current list of found mappings. The `callDepth` parameter keeps track of + how many matching pairs of vertices have been identified, and is used to + know when an isomorphism is found. Returns ``True`` if at least one + isomorphism was found or ``False`` if none were found. """ - - cython.declare(vertices1=list, new_terminals1=list, new_terminals2=list) cython.declare(vertex1=Vertex, vertex2=Vertex) - cython.declare(ismatch=cython.bint) - - # Make sure we don't get cause in an infinite recursive loop + cython.declare(mapping=dict) + + # The call depth should never be negative! if callDepth < 0: - logging.error("Recursing too deep. Now {0:d}".format(callDepth)) - if callDepth < -100: - raise Exception("Recursing infinitely deep!") + raise VF2Error('Negative call depth encountered in VF2_match().') # Done if we have mapped to all vertices in graph if callDepth == 0: - if not subgraph: - assert len(map21) == len(graph1.vertices), \ - "Calldepth mismatch: callDepth = {0:d}, len(map21) = {1:d}, len(map12) = {2:d}, len(graph1.vertices) = {3:d}, len(graph2.vertices) = {4:d}".format(callDepth, len(map21), len(map12), len(graph1.vertices), len(graph2.vertices)) - if findAll: - map21List.append(map21.copy()) - map12List.append(map12.copy()) - return True - else: - assert len(map12) == len(graph2.vertices), \ - "Calldepth mismatch: callDepth = {0:d}, len(map21) = {1:d}, len(map12) = {2:d}, len(graph1.vertices) = {3:d}, len(graph2.vertices) = {4:d}".format(callDepth, len(map21), len(map12), len(graph1.vertices), len(graph2.vertices)) - if findAll: - map21List.append(map21.copy()) - map12List.append(map12.copy()) - return True + if findAll: + mapping = {} + for vertex2 in graph2.vertices: + assert vertex2.mapping is not None + assert vertex2.mapping.mapping is vertex2 + mapping[vertex2.mapping] = vertex2 + mappingList.append(mapping) + return True # Create list of pairs of candidates for inclusion in mapping - # Note that the extra Python overhead is not worth making this a standalone - # method, so we simply put it inline here - # If we have terminals for both graphs, then use those as a basis for the - # pairs - if len(terminals1) > 0 and len(terminals2) > 0: - vertices1 = terminals1 - vertex2 = terminals2[0] - # Otherwise construct list from all *remaining* vertices (not matched) + vertices1 = [] + for vertex2 in graph2.vertices: + if vertex2.terminal: + # graph2 has terminals, so graph1 also must have terminals + for vertex1 in graph1.vertices: + if vertex1.terminal: + vertices1.append(vertex1) + break else: - # vertex2 is the lowest-labelled un-mapped vertex from graph2 - # Note that this assumes that graph2.vertices is properly sorted - vertices1 = [] - for vertex1 in graph1.vertices: - if vertex1 not in map21: - vertices1.append(vertex1) - for vertex2 in graph2.vertices: - if vertex2 not in map12: - break - else: - raise Exception("Could not find a pair to propose!") + # graph2 does not have terminals, so graph1 also must not have terminals + vertex2 = graph2.vertices[0] + vertices1 = graph1.vertices[:] for vertex1 in vertices1: - # propose a pairing - if __VF2_feasible(graph1, graph2, vertex1, vertex2, map21, map12, \ - terminals1, terminals2, subgraph): - # Update mapping accordingly - map21[vertex1] = vertex2 - map12[vertex2] = vertex1 - - # update terminals - new_terminals1 = __VF2_updateTerminals(graph1, map21, terminals1, vertex1) - new_terminals2 = __VF2_updateTerminals(graph2, map12, terminals2, vertex2) - + # Propose a pairing + if VF2_feasible(graph1, graph2, vertex1, vertex2, subgraph): + # Add proposed match to mapping + VF2_addToMapping(vertex1, vertex2) # Recurse - ismatch = __VF2_match(graph1, graph2, \ - map21, map12, new_terminals1, new_terminals2, subgraph, findAll, \ - map21List, map12List, callDepth-1) - if ismatch: + isMatch = VF2_match(graph1, graph2, subgraph, findAll, mappingList, callDepth-1) + if isMatch: if not findAll: return True # Undo proposed match - del map21[vertex1] - del map12[vertex2] - # changes to 'new_terminals' will be discarded and 'terminals' is unchanged + VF2_removeFromMapping(vertex1, vertex2) + # None of the proposed matches led to a complete isomorphism, so return False return False -def __VF2_terminals(graph, mapping): +def VF2_addToMapping(vertex1, vertex2): """ - For a given graph `graph` and associated partial mapping `mapping`, - generate a list of terminals, vertices that are directly connected to - vertices that have already been mapped. - - List is sorted (using key=__getSortLabel) before returning. + Add a pair of vertices `vertex1` and `vertex2` to the current mapping, + and update the terminals information for the neighbors of each vertex. """ + cython.declare(v=Vertex) + + # Map the vertices to one another + vertex1.mapping = vertex2 + vertex2.mapping = vertex1 + + # Remove these vertices from the set of terminals + vertex1.terminal = False + vertex2.terminal = False + + # Add any neighboring vertices not already in mapping to terminals + for v in vertex1.edges: + v.terminal = v.mapping is None + for v in vertex2.edges: + v.terminal = v.mapping is None - cython.declare(terminals=list) - terminals = list() - for vertex2 in graph.vertices: - if vertex2 not in mapping: - for vertex1 in mapping: - if vertex2 in graph.edges[vertex1]: - terminals.append(vertex2) - break - return terminals - -def __VF2_updateTerminals(graph, mapping, old_terminals, new_vertex): +def VF2_removeFromMapping(vertex1, vertex2): """ - For a given graph `graph` and associated partial mapping `mapping`, - *updates* a list of terminals, vertices that are directly connected to - vertices that have already been mapped. You have to pass it the previous - list of terminals `old_terminals` and the vertex `vertex` that has been - added to the mapping. Returns a new *copy* of the terminals. + Remove a pair of vertices `vertex1` and `vertex2` from the current mapping, + and update the terminals information for the neighbors of each vertex. """ - - cython.declare(terminals=list, vertex1=Vertex, vertex2=Vertex, edges=dict) - cython.declare(i=cython.int, sorting_label=cython.short, sorting_label2=cython.short) - - # Copy the old terminals, leaving out the new_vertex - terminals = old_terminals[:] - if new_vertex in terminals: terminals.remove(new_vertex) - - # Add the terminals of new_vertex - edges = graph.edges[new_vertex] - for vertex1 in edges: - if vertex1 not in mapping: # only add if not already mapped - # find spot in the sorted terminals list where we should put this vertex - sorting_label = vertex1.sortingLabel - i=0; sorting_label2=-1 # in case terminals list empty - for i in range(len(terminals)): - vertex2 = terminals[i] - sorting_label2 = vertex2.sortingLabel - if sorting_label2 >= sorting_label: - break - # else continue going through the list of terminals - else: # got to end of list without breaking, - # so add one to index to make sure vertex goes at end - i+=1 - if sorting_label2 == sorting_label: # this vertex already in terminals. - continue # try next vertex in graph[new_vertex] - - # insert vertex in right spot in terminals - terminals.insert(i,vertex1) - - return terminals - -################################################################################ + cython.declare(v=Vertex, v2=Vertex) + + # Unmap the vertices from one another + vertex1.mapping = None + vertex2.mapping = None + + # Restore these vertices to the set of terminals + for v in vertex1.edges: + if v.mapping is not None: + vertex1.terminal = True + break + else: + vertex1.terminal = False + for v in vertex2.edges: + if v.mapping is not None: + vertex2.terminal = True + break + else: + vertex2.terminal = False + + # Recompute the terminal status of any neighboring atoms + for v in vertex1.edges: + if v.mapping is not None: continue + for v2 in v.edges: + if v2.mapping is not None: + v.terminal = True + break + else: + v.terminal = False + for v in vertex2.edges: + if v.mapping is not None: continue + for v2 in v.edges: + if v2.mapping is not None: + v.terminal = True + break + else: + v.terminal = False From 5efe1141512d8ccf9291411bfdb722cba414d063 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 16:37:57 -0400 Subject: [PATCH 12/23] Updated symmetry number functions to work with new Graph format. These functions needed quite a bit of work to adjust to storing the bonds on the atoms itself. In particular, we now need to make deep copies of the Molecule object in more places, since adding and removing bonds now modified the Atom objects. This may need more attention in the future, but seems okay for now. --- rmgpy/molecule/symmetry.py | 52 ++++++++++++++++++++----------- unittest/molecule/symmetryTest.py | 20 ++++++------ 2 files changed, 43 insertions(+), 29 deletions(-) diff --git a/rmgpy/molecule/symmetry.py b/rmgpy/molecule/symmetry.py index 3dd8810576..8a263af027 100644 --- a/rmgpy/molecule/symmetry.py +++ b/rmgpy/molecule/symmetry.py @@ -41,7 +41,7 @@ def calculateAtomSymmetryNumber(molecule, atom): single = 0; double = 0; triple = 0; benzene = 0 numNeighbors = 0 - for bond in molecule.edges[atom].values(): + for bond in atom.edges.values(): if bond.isSingle(): single += 1 elif bond.isDouble(): double += 1 elif bond.isTriple(): triple += 1 @@ -52,8 +52,9 @@ def calculateAtomSymmetryNumber(molecule, atom): if numNeighbors < 2: return symmetryNumber # Create temporary structures for each functional group attached to atom - molecule = molecule.copy() - for atom2 in molecule.bonds[atom].keys(): molecule.removeBond(atom, atom2) + molecule0 = molecule + molecule = molecule0.copy(True) + atom = molecule.vertices[molecule0.vertices.index(atom)] molecule.removeAtom(atom) groups = molecule.split() @@ -108,7 +109,7 @@ def calculateBondSymmetryNumber(molecule, atom1, atom2): """ Return the symmetry number centered at `bond` in the structure. """ - bond = molecule.edges[atom1][atom2] + bond = atom1.edges[atom2] symmetryNumber = 1 if bond.isSingle() or bond.isDouble() or bond.isTriple(): if atom1.equivalent(atom2): @@ -123,9 +124,14 @@ def calculateBondSymmetryNumber(molecule, atom1, atom2): elif len(molecule.vertices) == 2: symmetryNumber = 2 else: - molecule = molecule.copy() - molecule.removeBond(atom1, atom2) - fragments = molecule.split() + molecule.removeBond(bond) + structure = molecule.copy(True) + molecule.addBond(bond) + + atom1 = structure.atoms[molecule.atoms.index(atom1)] + atom2 = structure.atoms[molecule.atoms.index(atom2)] + fragments = structure.split() + if len(fragments) != 2: return symmetryNumber fragment1, fragment2 = fragments @@ -149,7 +155,8 @@ def calculateBondSymmetryNumber(molecule, atom1, atom2): elif groups1[0].isIsomorphic(groups2[1]) and groups1[1].isIsomorphic(groups2[0]) and groups1[2].isIsomorphic(groups2[2]): symmetryNumber *= 2 elif groups1[0].isIsomorphic(groups2[2]) and groups1[1].isIsomorphic(groups2[0]) and groups1[2].isIsomorphic(groups2[1]): symmetryNumber *= 2 elif groups1[0].isIsomorphic(groups2[2]) and groups1[1].isIsomorphic(groups2[1]) and groups1[2].isIsomorphic(groups2[0]): symmetryNumber *= 2 - + + return symmetryNumber ################################################################################ @@ -189,9 +196,9 @@ def calculateAxisSymmetryNumber(molecule): # List all double bonds in the structure doubleBonds = [] - for atom1 in molecule.edges: - for atom2 in molecule.edges[atom1]: - if molecule.edges[atom1][atom2].isDouble() and molecule.vertices.index(atom1) < molecule.vertices.index(atom2): + for atom1 in molecule.vertices: + for atom2 in atom1.edges: + if atom1.edges[atom2].isDouble() and molecule.vertices.index(atom1) < molecule.vertices.index(atom2): doubleBonds.append((atom1, atom2)) # Search for adjacent double bonds @@ -228,7 +235,7 @@ def calculateAxisSymmetryNumber(molecule): # Do nothing if axis is in cycle found = False for atom1, atom2 in bonds: - if molecule.isBondInCycle(atom1, atom2): found = True + if molecule.isBondInCycle(atom1.edges[atom2]): found = True if found: continue # Find terminal atoms in axis @@ -241,12 +248,19 @@ def calculateAxisSymmetryNumber(molecule): if len(terminalAtoms) != 2: continue # Remove axis from (copy of) structure - structure = molecule.copy() + bondlist = [] for atom1, atom2 in bonds: - structure.removeBond(atom1, atom2) + bond = atom1.edges[atom2] + bondlist.append(bond) + molecule.removeBond(bond) + structure = molecule.copy(True) + terminalAtoms = [structure.vertices[molecule.vertices.index(atom)] for atom in terminalAtoms] + for bond in bondlist: + molecule.addBond(bond) + atomsToRemove = [] - for atom in structure.atoms: - if len(structure.bonds[atom]) == 0 and atom not in terminalAtoms: # it's not bonded to anything + for atom in structure.vertices: + if len(atom.edges) == 0 and atom not in terminalAtoms: # it's not bonded to anything atomsToRemove.append(atom) for atom in atomsToRemove: structure.removeAtom(atom) @@ -398,9 +412,9 @@ def calculateSymmetryNumber(molecule): if not molecule.isAtomInCycle(atom): symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) - for atom1 in molecule.edges: - for atom2 in molecule.edges[atom1]: - if molecule.vertices.index(atom1) < molecule.vertices.index(atom2) and not molecule.isBondInCycle(atom1, atom2): + for atom1 in molecule.vertices: + for atom2 in atom1.edges: + if molecule.vertices.index(atom1) < molecule.vertices.index(atom2) and not molecule.isBondInCycle(atom1.edges[atom2]): symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) symmetryNumber *= calculateAxisSymmetryNumber(molecule) diff --git a/unittest/molecule/symmetryTest.py b/unittest/molecule/symmetryTest.py index a3969fad1b..df08c9a253 100644 --- a/unittest/molecule/symmetryTest.py +++ b/unittest/molecule/symmetryTest.py @@ -75,8 +75,8 @@ def testBondSymmetryNumberEthane(self): """ molecule = Molecule().fromSMILES('CC') symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: + for atom1 in molecule.atoms: + for atom2 in atom1.bonds: if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) self.assertEqual(symmetryNumber, 2) @@ -87,8 +87,8 @@ def testBondSymmetryNumberPropane(self): """ molecule = Molecule().fromSMILES('CCC') symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: + for atom1 in molecule.atoms: + for atom2 in atom1.bonds: if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) self.assertEqual(symmetryNumber, 1) @@ -99,8 +99,8 @@ def testBondSymmetryNumberButane(self): """ molecule = Molecule().fromSMILES('CCCC') symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: + for atom1 in molecule.atoms: + for atom2 in atom1.bonds: if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) self.assertEqual(symmetryNumber, 2) @@ -111,8 +111,8 @@ def testBondSymmetryNumberEthylene(self): """ molecule = Molecule().fromSMILES('C=C') symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: + for atom1 in molecule.atoms: + for atom2 in atom1.bonds: if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) self.assertEqual(symmetryNumber, 2) @@ -123,8 +123,8 @@ def testBondSymmetryNumberAcetylene(self): """ molecule = Molecule().fromSMILES('C#C') symmetryNumber = 1 - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: + for atom1 in molecule.atoms: + for atom2 in atom1.bonds: if molecule.atoms.index(atom1) < molecule.atoms.index(atom2): symmetryNumber *= calculateBondSymmetryNumber(molecule, atom1, atom2) self.assertEqual(symmetryNumber, 2) From 21f57576841be70359ad861485a8e4a0a5b4c809 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 17:04:48 -0400 Subject: [PATCH 13/23] Added new rmgpy.molecule.adjlist module for working with adjacency lists. This functionality is used by both Molecule and Group objects. Before it was stored in the rmgpy.molecule.group module, a bit of an unfortunate compromise. Now that rmgpy.molecule is a subpackage, we have space to place the adjacency list functionality in a module of its own, as done here. --- rmgpy/molecule/__init__.py | 1 + rmgpy/molecule/adjlist.py | 338 +++++++++++++++++++++++++++++++++++++ rmgpy/molecule/group.pxd | 6 - rmgpy/molecule/group.py | 279 +----------------------------- rmgpy/molecule/molecule.py | 8 +- 5 files changed, 348 insertions(+), 284 deletions(-) create mode 100644 rmgpy/molecule/adjlist.py diff --git a/rmgpy/molecule/__init__.py b/rmgpy/molecule/__init__.py index 10474c6ff5..f355d0cf1d 100644 --- a/rmgpy/molecule/__init__.py +++ b/rmgpy/molecule/__init__.py @@ -28,6 +28,7 @@ # ################################################################################ +from .adjlist import * from .atomtype import * from .element import * from .molecule import * diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py new file mode 100644 index 0000000000..dee9bb47da --- /dev/null +++ b/rmgpy/molecule/adjlist.py @@ -0,0 +1,338 @@ +################################################################################ +# +# ChemPy - A chemistry toolkit for Python +# +# Copyright (c) 2012 by Joshua W. Allen (jwallen@mit.edu) +# +# Permission is hereby granted, free of charge, to any person obtaining a +# copy of this software and associated documentation files (the "Software"), +# to deal in the Software without restriction, including without limitation +# the rights to use, copy, modify, merge, publish, distribute, sublicense, +# and/or sell copies of the Software, and to permit persons to whom the +# Software is furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +# DEALINGS IN THE SOFTWARE. +# +################################################################################ + +""" +This module contains functionality for reading from and writing to the +adjacency list format used by Reaction Mechanism Generator (RMG). +""" + +from .molecule import Atom, Bond +from .group import GroupAtom, GroupBond +#import chempy.molecule.atomtype as atomtypes + +################################################################################ + +class InvalidAdjacencyListError(Exception): + """ + An exception used to indicate that an RMG-style adjacency list is invalid. + Pass a string describing the reason the adjacency list is invalid + """ + pass + +################################################################################ + +def fromAdjacencyList(adjlist, group=False): + """ + Convert a string adjacency list `adjlist` into a set of :class:`Atom` and + :class:`Bond` objects. + """ + atoms = [] + atomdict = {} + bonds = {} + + try: + + adjlist = adjlist.strip() + lines = adjlist.splitlines() + if adjlist == '' or len(lines) == 0: + raise InvalidAdjacencyListError('Empty adjacency list.') + + # Skip the first line if it contains a label + if len(lines[0].split()) == 1: + label = lines.pop(0) + if len(lines) == 0: + raise InvalidAdjacencyListError('No atoms specified in adjacency list.') + + # Iterate over the remaining lines, generating Atom or GroupAtom objects + for line in lines: + + # Sometimes commas are used to delimit bonds in the bond list, + # so replace them just in case + line = line.replace('},{', '} {') + + data = line.split() + + # Skip if blank line + if len(data) == 0: continue + + # First item is index for atom + # Sometimes these have a trailing period (as if in a numbered list), + # so remove it just in case + aid = int(data[0].strip('.')) + + # If second item starts with '*', then atom is labeled + label = ''; index = 1 + if data[1][0] == '*': + label = data[1] + index += 1 + + # Next is the element or functional group element + # A list can be specified with the {,} syntax + atomType = data[index] + if atomType[0] == '{': + atomType = atomType[1:-1].split(',') + else: + atomType = [atomType] + + # Next is the element or atom type + # A list can be specified with the {,} syntax + atomType = data[index] + if atomType[0] == '{': + atomType = atomType[1:-1].split(',') + else: + atomType = [atomType] + index += 1 + + # Next is the electron state + radicalElectrons = []; spinMultiplicity = [] + elecState = data[index].upper() + if elecState[0] == '{': + elecState = elecState[1:-1].split(',') + else: + elecState = [elecState] + for e in elecState: + if e == '0': + radicalElectrons.append(0); spinMultiplicity.append(1) + elif e == '1': + radicalElectrons.append(1); spinMultiplicity.append(2) + elif e == '2': + radicalElectrons.append(2); spinMultiplicity.append(1) + radicalElectrons.append(2); spinMultiplicity.append(3) + elif e == '2S': + radicalElectrons.append(2); spinMultiplicity.append(1) + elif e == '2T': + radicalElectrons.append(2); spinMultiplicity.append(3) + elif e == '3': + radicalElectrons.append(3); spinMultiplicity.append(4) + elif e == '4': + radicalElectrons.append(4); spinMultiplicity.append(5) + index += 1 + + # Create a new atom based on the above information + if group: + atom = GroupAtom(atomType, radicalElectrons, spinMultiplicity, [0 for e in radicalElectrons], label) + else: + atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, label) + + # Add the atom to the list + atoms.append(atom) + atomdict[aid] = atom + + # Process list of bonds + bonds[aid] = {} + for datum in data[index:]: + + # Sometimes commas are used to delimit bonds in the bond list, + # so strip them just in case + datum = datum.strip(',') + + aid2, comma, order = datum[1:-1].partition(',') + aid2 = int(aid2) + if aid == aid2: + raise InvalidAdjacencyListError('Attempted to create a bond between atom {0:d} and itself.'.format(aid)) + + if order[0] == '{': + order = order[1:-1].split(',') + else: + order = [order] + + bonds[aid][aid2] = order + + # Check consistency using bonddict + for atom1 in bonds: + for atom2 in bonds[atom1]: + if atom2 not in bonds: + raise InvalidAdjacencyListError('Atom {0:d} not in bond dictionary.'.format(atom2)) + elif atom1 not in bonds[atom2]: + raise InvalidAdjacencyListError('Found bond between {0:d} and {1:d}, but not the reverse.'.format(atom1, atom2)) + elif bonds[atom1][atom2] != bonds[atom2][atom1]: + raise InvalidAdjacencyListError('Found bonds between {0:d} and {1:d}, but of different orders "{2}" and "{3}".'.format(atom1, atom2, bonds[atom1][atom2], bonds[atom2][atom1])) + + # Convert bonddict to use Atom[group] and Bond[group] objects + atomkeys = atomdict.keys() + atomkeys.sort() + for aid1 in atomkeys: + atomkeys2 = bonds[aid1].keys() + atomkeys2.sort() + for aid2 in atomkeys2: + if aid1 < aid2: + atom1 = atomdict[aid1] + atom2 = atomdict[aid2] + order = bonds[aid1][aid2] + if group: + bond = GroupBond(atom1, atom2, order) + elif len(order) == 1: + bond = Bond(atom1, atom2, order[0]) + else: + raise InvalidAdjacencyListError('Multiple bond orders specified for an atom in a Molecule.') + atom1.edges[atom2] = bond + atom2.edges[atom1] = bond + + # Add explicit hydrogen atoms to complete structure if desired + if not group: + valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0} + orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} + newAtoms = [] + for atom in atoms: + try: + valence = valences[atom.symbol] + except KeyError: + raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown valence for atom "{0}".'.format(atom.symbol)) + radical = atom.radicalElectrons + order = 0 + for atom2, bond in atom.bonds.items(): + order += orders[bond.order] + count = valence - radical - int(order) + for i in range(count): + a = Atom('H', 0, 1, 0, '') + b = Bond(atom, a, 'S') + newAtoms.append(a) + atom.bonds[a] = b + a.bonds[atom] = b + atoms.extend(newAtoms) + + except InvalidAdjacencyListError: + print adjlist + raise + + return atoms + +################################################################################ + +def getElectronState(radicalElectrons, spinMultiplicity): + """ + Return the electron state corresponding to the given number of radical + electrons `radicalElectrons` and spin multiplicity `spinMultiplicity`. + Raise a :class:`ValueError` if the electron state cannot be determined. + """ + electronState = '' + if radicalElectrons == 0: + electronState = '0' + elif radicalElectrons == 1: + electronState = '1' + elif radicalElectrons == 2 and spinMultiplicity == 1: + electronState = '2S' + elif radicalElectrons == 2 and spinMultiplicity == 3: + electronState = '2T' + elif radicalElectrons == 3: + electronState = '3' + elif radicalElectrons == 4: + electronState = '4' + else: + raise ValueError('Unable to determine electron state for {0:d} radical electrons with spin multiplicity of {1:d}.'.format(radicalElectrons, spinMultiplicity)) + return electronState + +def toAdjacencyList(atoms, label=None, group=False, removeH=False): + """ + Convert a chemical graph defined by a list of `atoms` into a string + adjacency list. + """ + adjlist = '' + + # Don't remove hydrogen atoms if the molecule consists only of hydrogen atoms + try: + if removeH and all([atom.element.symbol == 'H' for atom in atoms]): removeH = False + except AttributeError: + pass + + if label: adjlist += label + '\n' + + # Determine the numbers to use for each atom + atomNumbers = {} + index = 0 + for atom in atoms: + if removeH and atom.element.symbol == 'H' and atom.label == '': continue + atomNumbers[atom] = '{0:d}'.format(index + 1) + index += 1 + + atomLabels = {atom: '{0}'.format(atom.label) for atom in atomNumbers} + + atomTypes = {} + atomElectronStates = {} + if group: + for atom in atomNumbers: + # Atom type(s) + if len(atom.atomType) == 1: + atomTypes[atom] = atom.atomType[0].label + else: + atomTypes[atom] = '{{{0}}}'.format(','.join([a.label for a in atom.atomType])) + # Electron state(s) + if len(atom.radicalElectrons) == 1: + atomElectronStates[atom] = getElectronState(atom.radicalElectrons[0], atom.spinMultiplicity[0]) + else: + atomElectronStates[atom] = '{{{0}}}'.format(','.join([getElectronState(radical, spin) for radical, spin in zip(atom.radicalElectrons, atom.spinMultiplicity)])) + else: + for atom in atomNumbers: + # Atom type + atomTypes[atom] = '{0}'.format(atom.element.symbol) + # Electron state(s) + atomElectronStates[atom] = '{0}'.format(getElectronState(atom.radicalElectrons, atom.spinMultiplicity)) + + # Determine field widths + atomNumberWidth = max([len(s) for s in atomNumbers.values()]) + 1 + atomLabelWidth = max([len(s) for s in atomLabels.values()]) + if atomLabelWidth > 0: atomLabelWidth += 1 + atomTypeWidth = max([len(s) for s in atomTypes.values()]) + 1 + atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()]) + + # Assemble the adjacency list + for atom in atoms: + if atom not in atomNumbers: continue + + # Atom number + adjlist += '{0:<{1:d}}'.format(atomNumbers[atom], atomNumberWidth) + # Atom label + adjlist += '{0:<{1:d}}'.format(atomLabels[atom], atomLabelWidth) + # Atom type(s) + adjlist += '{0:<{1:d}}'.format(atomTypes[atom], atomTypeWidth) + # Electron state(s) + adjlist += '{0:<{1:d}}'.format(atomElectronStates[atom], atomElectronStateWidth) + + # Bonds list + atoms2 = atom.bonds.keys() + # sort them the same way as the atoms + atoms2.sort(key=atoms.index) + + for atom2 in atoms2: + if atom2 not in atomNumbers: continue + + bond = atom.bonds[atom2] + adjlist += ' {{{0},'.format(atomNumbers[atom2]) + + # Bond type(s) + if group: + if len(bond.order) == 1: + adjlist += bond.order[0] + else: + adjlist += '{{{0}}}'.format(','.join(bond.order)) + else: + adjlist += bond.order + adjlist += '}' + + # Each atom begins on a new line + adjlist += '\n' + + return adjlist diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd index 34387aa48c..c5d8623952 100644 --- a/rmgpy/molecule/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -113,9 +113,3 @@ cdef class Group(Graph): cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?) cpdef tuple findSubgraphIsomorphisms(self, Graph other, dict initialMap=?) - -################################################################################ - -cpdef fromAdjacencyList(str adjlist, bint group=?, bint addH=?) - -cpdef toAdjacencyList(Graph molecule, str label=?, bint group=?, bint removeH=?) diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index d035e49a6a..7aaf218c3d 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -595,7 +595,8 @@ def fromAdjacencyList(self, adjlist): Skips the first line (assuming it's a label) unless `withLabel` is ``False``. """ - self.vertices, self.edges = fromAdjacencyList(adjlist, group=True, addH=False) + from .adjlist import fromAdjacencyList + self.vertices = fromAdjacencyList(adjlist, group=True) self.updateConnectivityValues() return self @@ -603,7 +604,8 @@ def toAdjacencyList(self, label=''): """ Convert the molecular structure to a string adjacency list. """ - return toAdjacencyList(self, label='', group=True) + from .adjlist import toAdjacencyList + return toAdjacencyList(self.vertices, label='', group=True) def isIsomorphic(self, other, initialMap=None): """ @@ -669,276 +671,3 @@ def findSubgraphIsomorphisms(self, other, initialMap=None): raise TypeError('Got a {0} object for parameter "other", when a Group object is required.'.format(other.__class__)) # Do the isomorphism comparison return Graph.findSubgraphIsomorphisms(self, other, initialMap) - -################################################################################ - -class InvalidAdjacencyListError(Exception): - """ - An exception used to indicate that an RMG-style adjacency list is invalid. - Pass a string giving specifics about the particular exceptional behavior. - """ - pass - -def fromAdjacencyList(adjlist, group=False, addH=False): - """ - Convert a string adjacency list `adjlist` into a set of :class:`Atom` and - :class:`Bond` objects (if `group` is ``False``) or a set of - :class:`GroupAtom` and :class:`GroupBond` objects (if `group` is - ``True``). Only adds hydrogen atoms if `addH` is ``True``. Skips the first - line (assuming it's a label) unless `withLabel` is ``False``. - """ - - from molecule import Atom, Bond - - atoms = []; atomdict = {}; bonds = {} - - try: - - adjlist = adjlist.strip() - if adjlist == '': - raise InvalidAdjacencyListError('Empty adjacency list.') - - lines = adjlist.splitlines() - # Skip the first line if it contains a label - if len(lines) > 0 and len(lines[0].split()) == 1: - label = lines.pop(0) - # Iterate over the remaining lines, generating Atom or GroupAtom objects - if len(lines) == 0: - raise InvalidAdjacencyListError('No atoms specified in adjacency list.') - for line in lines: - - # Sometimes commas are used to delimit bonds in the bond list, - # so replace them just in case - line = line.replace('},{', '} {') - - data = line.split() - - # Skip if blank line - if len(data) == 0: continue - - # First item is index for atom - # Sometimes these have a trailing period (as if in a numbered list), - # so remove it just in case - aid = int(data[0].strip('.')) - - # If second item starts with '*', then atom is labeled - label = ''; index = 1 - if data[1][0] == '*': - label = data[1]; index = 2 - - # Next is the element or functional group element - # A list can be specified with the {,} syntax - atomType = data[index] - if atomType[0] == '{': - atomType = atomType[1:-1].split(',') - else: - atomType = [atomType] - - # Next is the electron state - radicalElectrons = []; spinMultiplicity = [] - elecState = data[index+1].upper() - if elecState[0] == '{': - elecState = elecState[1:-1].split(',') - else: - elecState = [elecState] - for e in elecState: - if e == '0': - radicalElectrons.append(0); spinMultiplicity.append(1) - elif e == '1': - radicalElectrons.append(1); spinMultiplicity.append(2) - elif e == '2': - radicalElectrons.append(2); spinMultiplicity.append(1) - radicalElectrons.append(2); spinMultiplicity.append(3) - elif e == '2S': - radicalElectrons.append(2); spinMultiplicity.append(1) - elif e == '2T': - radicalElectrons.append(2); spinMultiplicity.append(3) - elif e == '3': - radicalElectrons.append(3); spinMultiplicity.append(4) - elif e == '4': - radicalElectrons.append(4); spinMultiplicity.append(5) - - # Create a new atom based on the above information - if group: - atom = GroupAtom(atomType, radicalElectrons, spinMultiplicity, [0 for e in radicalElectrons], label) - else: - atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, label) - - # Add the atom to the list - atoms.append(atom) - atomdict[aid] = atom - - # Process list of bonds - bonds[aid] = {} - for datum in data[index+2:]: - - # Sometimes commas are used to delimit bonds in the bond list, - # so strip them just in case - datum = datum.strip(',') - - aid2, comma, order = datum[1:-1].partition(',') - aid2 = int(aid2) - if aid == aid2: - raise InvalidAdjacencyListError('Attempted to create a bond between atom {0:d} and itself.'.format(aid)) - - if order[0] == '{': - order = order[1:-1].split(',') - else: - order = [order] - - bonds[aid][aid2] = order - - # Check consistency using bonddict - for atom1 in bonds: - for atom2 in bonds[atom1]: - if atom2 not in bonds: - raise InvalidAdjacencyListError('Atom {0:d} not in bond dictionary.'.format(atom2)) - elif atom1 not in bonds[atom2]: - raise InvalidAdjacencyListError('Found bond between {0:d} and {1:d}, but not the reverse.'.format(atom1, atom2)) - elif bonds[atom1][atom2] != bonds[atom2][atom1]: - raise InvalidAdjacencyListError('Found bonds between {0:d} and {1:d}, but of different orders "{2}" and "{3}".'.format(atom1, atom2, bonds[atom1][atom2], bonds[atom2][atom1])) - - # Convert bonddict to use Atom[group] and Bond[group] objects - atomkeys = atomdict.keys() - atomkeys.sort() - for aid1 in atomkeys: - bonds[atomdict[aid1]] = {} - atomkeys2 = bonds[aid1].keys() - atomkeys2.sort() - for aid2 in atomkeys2: - if aid1 < aid2: - order = bonds[aid1][aid2] - if group: - bonds[atomdict[aid1]][atomdict[aid2]] = GroupBond(order) - elif len(order) == 1: - bonds[atomdict[aid1]][atomdict[aid2]] = Bond(order[0]) - else: - raise InvalidAdjacencyListError('Multiple bond orders specified for an atom in a Molecule.') - else: - bonds[atomdict[aid1]][atomdict[aid2]] = bonds[atomdict[aid2]][atomdict[aid1]] - del bonds[aid1] - - # Add explicit hydrogen atoms to complete structure if desired - if addH and not group: - valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0} - orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} - newAtoms = [] - for atom in atoms: - try: - valence = valences[atom.symbol] - except KeyError: - raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown valence for atom "{0}".'.format(atom.symbol)) - radical = atom.radicalElectrons - order = 0 - for atom2, bond in bonds[atom].iteritems(): - order += orders[bond.order] - count = valence - radical - int(order) - for i in range(count): - a = Atom('H', 0, 1, 0, '') - b = Bond('S') - newAtoms.append(a) - bonds[atom][a] = b - bonds[a] = {atom: b} - atoms.extend(newAtoms) - - except InvalidAdjacencyListError: - print adjlist - raise - - return atoms, bonds - -def toAdjacencyList(molecule, label='', group=False, removeH=False): - """ - Convert the `molecule` object to an adjacency list. `group` specifies - whether the graph object is a complete molecule (if ``False``) or a - substructure group (if ``True``). The `label` parameter is an optional - string to put as the first line of the adjacency list; if set to the empty - string, this line will be omitted. If `removeH` is ``True``, hydrogen atoms - (that do not have labels) will not be printed; this is a valid shorthand, - as they can usually be inferred as long as the free electron numbers are - accurate. - """ - - adjlist = '' - - # Don't remove hydrogen atoms if the molecule consists only of hydrogen atoms - try: - if removeH and all([atom.isHydrogen() for atom in molecule.atoms]): removeH = False - except AttributeError: - pass - - if label != '': adjlist += label + '\n' - - molecule.updateConnectivityValues() # so we can sort by them - molecule.sortAtoms() - atoms = molecule.atoms - bonds = molecule.bonds - - # Determine the numbers to use for each atom - atomNumbers = {}; index = 0 - for atom in atoms: - if removeH and atom.isHydrogen() and atom.label=='': continue - atomNumbers[atom] = index + 1 - index += 1 - - for atom in atoms: - if removeH and atom.isHydrogen() and atom.label=='': continue - - # Atom number - adjlist += '{0:<2} '.format(atomNumbers[atom]) - - # Atom label - adjlist += '{0:<2} '.format(atom.label) - - if group: - # Atom type(s) - if len(atom.atomType) == 1: - adjlist += atom.atomType[0].label + ' ' - else: - adjlist += '{{{0}}} '.format(','.join([a.label for a in atom.atomType])) - # Electron state(s) - if len(atom.radicalElectrons) > 1: adjlist += '{' - for radical, spin in zip(atom.radicalElectrons, atom.spinMultiplicity): - if radical == 0: adjlist += '0' - elif radical == 1: adjlist += '1' - elif radical == 2 and spin == 1: adjlist += '2S' - elif radical == 2 and spin == 3: adjlist += '2T' - elif radical == 3: adjlist += '3' - elif radical == 4: adjlist += '4' - if len(atom.radicalElectrons) > 1: adjlist += ',' - if len(atom.radicalElectrons) > 1: adjlist = adjlist[0:-1] + '}' - else: - # Atom type - adjlist += "{0:<5} ".format(atom.symbol) - # Electron state(s) - if atom.radicalElectrons == 0: adjlist += '0' - elif atom.radicalElectrons == 1: adjlist += '1' - elif atom.radicalElectrons == 2 and atom.spinMultiplicity == 1: adjlist += '2S' - elif atom.radicalElectrons == 2 and atom.spinMultiplicity == 3: adjlist += '2T' - elif atom.radicalElectrons == 3: adjlist += '3' - elif atom.radicalElectrons == 4: adjlist += '4' - - # Bonds list - atoms2 = bonds[atom].keys() - # sort them the same way as the atoms - atoms2.sort(key=atoms.index) - - for atom2 in atoms2: - if removeH and atom2.isHydrogen() and atom2.label=='': continue - bond = bonds[atom][atom2] - adjlist += ' {{{0:d},'.format(atomNumbers[atom2]) - - # Bond type(s) - if group: - if len(bond.order) == 1: - adjlist += bond.order[0] - else: - adjlist += '{{{0}}}'.format(','.join(bond.order)) - else: - adjlist += bond.order - adjlist += '}' - - # Each atom begins on a new line - adjlist += '\n' - - return adjlist diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 321784decc..d25e523884 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -40,7 +40,7 @@ import re import element as elements from .graph import Vertex, Edge, Graph -from .group import GroupAtom, GroupBond, Group, ActionError, fromAdjacencyList, toAdjacencyList +from .group import GroupAtom, GroupBond, Group, ActionError from .atomtype import AtomType, atomTypes, getAtomType ################################################################################ @@ -882,7 +882,8 @@ def fromAdjacencyList(self, adjlist): Skips the first line (assuming it's a label) unless `withLabel` is ``False``. """ - self.vertices, self.edges = fromAdjacencyList(adjlist, False, True) + from .adjlist import fromAdjacencyList + self.vertices = fromAdjacencyList(adjlist, False) self.updateConnectivityValues() self.updateAtomTypes() return self @@ -1004,7 +1005,8 @@ def toAdjacencyList(self, label='', removeH=False): """ Convert the molecular structure to a string adjacency list. """ - result = toAdjacencyList(self, label=label, group=False, removeH=removeH) + from .adjlist import toAdjacencyList + result = toAdjacencyList(self.vertices, label=label, group=False, removeH=removeH) return result def isLinear(self): From d07a1df725573d68cb87d51bb9fa960cc34abb5b Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 17:15:53 -0400 Subject: [PATCH 14/23] Updated Atom, Bond, and Molecule to work with new Graph format. As before, this is mostly updating to reflect that the edges are now stored on the Atom objects instead of the Molecule objects. Some of the corresponding unit tests also needed a bit of work. --- rmgpy/molecule/molecule.pxd | 40 ++++++------ rmgpy/molecule/molecule.py | 101 +++++++++++++++--------------- unittest/molecule/moleculeTest.py | 88 ++++++++++---------------- 3 files changed, 104 insertions(+), 125 deletions(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index ac5b6cf53d..1e640c465a 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -40,11 +40,11 @@ cdef class Atom(Vertex): cdef public str label cdef public AtomType atomType - cpdef bint equivalent(self, Vertex other) + cpdef bint equivalent(self, Vertex other) except -2 - cpdef bint isSpecificCaseOf(self, Vertex other) + cpdef bint isSpecificCaseOf(self, Vertex other) except -2 - cpdef Atom copy(self) + cpdef Vertex copy(self) cpdef bint isHydrogen(self) @@ -60,17 +60,17 @@ cdef class Bond(Edge): cdef public str order - cpdef bint equivalent(self, Edge other) + cpdef bint equivalent(self, Edge other) except -2 - cpdef bint isSpecificCaseOf(self, Edge other) + cpdef bint isSpecificCaseOf(self, Edge other) except -2 - cpdef Bond copy(self) + cpdef Edge copy(self) - cpdef bint isSingle(self) + cpdef bint isSingle(self) except -2 - cpdef bint isDouble(self) + cpdef bint isDouble(self) except -2 - cpdef bint isTriple(self) + cpdef bint isTriple(self) except -2 ################################################################################ @@ -81,7 +81,7 @@ cdef class Molecule(Graph): cpdef addAtom(self, Atom atom) - cpdef addBond(self, Atom atom1, Atom atom2, Bond bond) + cpdef addBond(self, Bond bond) cpdef dict getBonds(self, Atom atom) @@ -93,7 +93,7 @@ cdef class Molecule(Graph): cpdef removeAtom(self, Atom atom) - cpdef removeBond(self, Atom atom1, Atom atom2) + cpdef removeBond(self, Bond bond) cpdef sortAtoms(self) @@ -107,23 +107,23 @@ cdef class Molecule(Graph): cpdef clearLabeledAtoms(self) - cpdef bint containsLabeledAtom(self, str label) + cpdef bint containsLabeledAtom(self, str label) except -2 cpdef Atom getLabeledAtom(self, str label) cpdef dict getLabeledAtoms(self) - cpdef bint isIsomorphic(self, Graph other, dict initialMap=?) + cpdef bint isIsomorphic(self, Graph other, dict initialMap=?) except -2 - cpdef tuple findIsomorphism(self, Graph other, dict initialMap=?) + cpdef list findIsomorphism(self, Graph other, dict initialMap=?) - cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?) + cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?) except -2 - cpdef tuple findSubgraphIsomorphisms(self, Graph other, dict initialMap=?) + cpdef list findSubgraphIsomorphisms(self, Graph other, dict initialMap=?) - cpdef bint isAtomInCycle(self, Atom atom) + cpdef bint isAtomInCycle(self, Atom atom) except -2 - cpdef bint isBondInCycle(self, Atom atom1, Atom atom2) + cpdef bint isBondInCycle(self, Bond bond) except -2 cpdef draw(self, str path) @@ -153,9 +153,9 @@ cdef class Molecule(Graph): cpdef toAdjacencyList(self, str label=?, bint removeH=?) - cpdef bint isLinear(self) + cpdef bint isLinear(self) except -2 - cpdef int countInternalRotors(self) + cpdef int countInternalRotors(self) except -2 cpdef getAdjacentResonanceIsomers(self) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index d25e523884..01e28baa75 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -81,23 +81,24 @@ def __str__(self): """ Return a human-readable string representation of the object. """ - return "".format( - str(self.element) + - '.' * self.radicalElectrons + - '+' * self.charge if self.charge > 0 else '-' * -self.charge + return '{0}{1}{2}'.format( + str(self.element), + '.' * self.radicalElectrons, + '+' * self.charge if self.charge > 0 else '-' * -self.charge, ) def __repr__(self): """ Return a representation that can be used to reconstruct the object. """ - return 'Atom(element="{0}", radicalElectrons={1}, spinMultiplicity={2}, charge={3}, label="{4}")'.format(self.element, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label) + return "".format(str(self)) def __reduce__(self): """ A helper function used when pickling an object. """ d = { + 'edges': self.edges, 'connectivity1': self.connectivity1, 'connectivity2': self.connectivity2, 'connectivity3': self.connectivity3, @@ -110,6 +111,7 @@ def __setstate__(self, d): """ A helper function used when unpickling an object. """ + self.edges = d['edges'] self.connectivity1 = d['connectivity1'] self.connectivity2 = d['connectivity2'] self.connectivity3 = d['connectivity3'] @@ -125,6 +127,9 @@ def number(self): return self.element.number @property def symbol(self): return self.element.symbol + @property + def bonds(self): return self.edges + def equivalent(self, other): """ Return ``True`` if `other` is indistinguishable from this atom, or @@ -289,27 +294,35 @@ class Bond(Edge): """ - def __init__(self, order=1): - Edge.__init__(self) + def __init__(self, atom1, atom2, order=1): + Edge.__init__(self, atom1, atom2) self.order = order def __str__(self): """ Return a human-readable string representation of the object. """ - return ''.format(self.order) + return self.order def __repr__(self): """ Return a representation that can be used to reconstruct the object. """ - return 'Bond(order="{0}")'.format(self.order) + return ''.format(self.order) def __reduce__(self): """ A helper function used when pickling an object. """ - return (Bond, (self.order,)) + return (Bond, (self.vertex1, self.vertex2, self.order)) + + @property + def atom1(self): + return self.vertex1 + + @property + def atom2(self): + return self.vertex2 def equivalent(self, other): """ @@ -339,7 +352,7 @@ def copy(self): Generate a deep copy of the current bond. Modifying the attributes of the copy will not affect the original. """ - return Bond(self.order) + return Bond(self.vertex1, self.vertex2, self.order) def isSingle(self): """ @@ -443,8 +456,8 @@ class Molecule(Graph): `InChI` string representing the molecular structure. """ - def __init__(self, atoms=None, bonds=None, symmetry=1, SMILES='', InChI=''): - Graph.__init__(self, atoms, bonds) + def __init__(self, atoms=None, symmetry=1, SMILES='', InChI=''): + Graph.__init__(self, atoms) self.symmetryNumber = symmetry if SMILES != '': self.fromSMILES(SMILES) elif InChI != '': self.fromInChI(InChI) @@ -465,28 +478,24 @@ def __reduce__(self): """ A helper function used when pickling an object. """ - return (Molecule, (self.vertices, self.edges, self.symmetryNumber)) + return (Molecule, (self.vertices, self.symmetryNumber)) def __getAtoms(self): return self.vertices def __setAtoms(self, atoms): self.vertices = atoms atoms = property(__getAtoms, __setAtoms) - def __getBonds(self): return self.edges - def __setBonds(self, bonds): self.edges = bonds - bonds = property(__getBonds, __setBonds) - def addAtom(self, atom): """ Add an `atom` to the graph. The atom is initialized with no bonds. """ return self.addVertex(atom) - def addBond(self, atom1, atom2, bond): + def addBond(self, bond): """ Add a `bond` to the graph as an edge connecting the two atoms `atom1` and `atom2`. """ - return self.addEdge(atom1, atom2, bond) + return self.addEdge(bond) def getBonds(self, atom): """ @@ -522,13 +531,13 @@ def removeAtom(self, atom): """ return self.removeVertex(atom) - def removeBond(self, atom1, atom2): + def removeBond(self, bond): """ Remove the bond between atoms `atom1` and `atom2` from the graph. Does not remove atoms that no longer have any bonds as a result of this removal. """ - return self.removeEdge(atom1, atom2) + return self.removeEdge(bond) def sortAtoms(self): """ @@ -564,7 +573,7 @@ def copy(self, deep=False): """ other = cython.declare(Molecule) g = Graph.copy(self, deep) - other = Molecule(g.vertices, g.edges) + other = Molecule(g.vertices) return other def merge(self, other): @@ -573,7 +582,7 @@ def merge(self, other): object. The merged :class:`Molecule` object is returned. """ g = Graph.merge(self, other) - molecule = Molecule(atoms=g.vertices, bonds=g.edges) + molecule = Molecule(atoms=g.vertices) return molecule def split(self): @@ -584,7 +593,7 @@ def split(self): graphs = Graph.split(self) molecules = [] for g in graphs: - molecule = Molecule(atoms=g.vertices, bonds=g.edges) + molecule = Molecule(atoms=g.vertices) molecules.append(molecule) return molecules @@ -605,7 +614,7 @@ def deleteHydrogens(self): hydrogens = [] for atom in self.vertices: if atom.isHydrogen() and atom.label == '': - neighbor = self.edges[atom].keys()[0] + neighbor = atom.edges.keys()[0] hydrogens.append(atom) # Remove the hydrogen atoms from the structure for atom in hydrogens: @@ -618,7 +627,7 @@ def updateAtomTypes(self): environment) and complete (i.e. are as detailed as possible). """ for atom in self.vertices: - atom.atomType = getAtomType(atom, self.edges[atom]) + atom.atomType = getAtomType(atom, atom.edges) def clearLabeledAtoms(self): """ @@ -736,12 +745,12 @@ def isAtomInCycle(self, atom): """ return self.isVertexInCycle(atom) - def isBondInCycle(self, atom1, atom2): + def isBondInCycle(self, bond): """ Return :data:`True` if the bond between atoms `atom1` and `atom2` is in one or more cycles in the graph, or :data:`False` if not. """ - return self.isEdgeInCycle(atom1, atom2) + return self.isEdgeInCycle(bond) def draw(self, path): """ @@ -815,7 +824,6 @@ def fromOBMol(self, obmol): cython.declare(atom=Atom, atom1=Atom, atom2=Atom, bond=Bond) self.vertices = [] - self.edges = {} # Add hydrogen atoms to complete molecule if needed obmol.AddHydrogens() @@ -849,7 +857,6 @@ def fromOBMol(self, obmol): atom = Atom(element, radicalElectrons, spinMultiplicity, charge) self.vertices.append(atom) - self.edges[atom] = {} # Add bonds by iterating again through atoms for j in range(0, i): @@ -864,11 +871,8 @@ def fromOBMol(self, obmol): elif obbond.IsTriple(): order = 'T' elif obbond.IsAromatic(): order = 'B' - bond = Bond(order) - atom1 = self.vertices[i] - atom2 = self.vertices[j] - self.edges[atom1][atom2] = bond - self.edges[atom2][atom1] = bond + bond = Bond(self.vertices[i], self.vertices[j], order) + self.addBond(bond) # Set atom types and connectivity values self.updateConnectivityValues() @@ -981,7 +985,6 @@ def toOBMol(self): self.sortAtoms() atoms = self.vertices - bonds = self.edges obmol = openbabel.OBMol() for atom in atoms: @@ -989,8 +992,8 @@ def toOBMol(self): a.SetAtomicNum(atom.number) a.SetFormalCharge(atom.charge) orders = {'S': 1, 'D': 2, 'T': 3, 'B': 5} - for atom1, bonds in bonds.iteritems(): - for atom2, bond in bonds.iteritems(): + for atom1 in self.vertices: + for atom2, bond in atom1.edges.iteritems(): index1 = atoms.index(atom1) index2 = atoms.index(atom2) if index1 < index2: @@ -1029,15 +1032,15 @@ def isLinear(self): # True if all bonds are double bonds (e.g. O=C=O) allDoubleBonds = True - for atom1 in self.edges: - for bond in self.edges[atom1].values(): + for atom1 in self.vertices: + for bond in atom1.edges.values(): if not bond.isDouble(): allDoubleBonds = False if allDoubleBonds: return True # True if alternating single-triple bonds (e.g. H-C#C-H) # This test requires explicit hydrogen atoms for atom in self.vertices: - bonds = self.edges[atom].values() + bonds = atom.edges.values() if len(bonds)==1: continue # ok, next atom if len(bonds)>2: @@ -1061,10 +1064,10 @@ def countInternalRotors(self): are considered to be internal rotors. """ count = 0 - for atom1 in self.edges: - for atom2, bond in self.edges[atom1].iteritems(): - if self.vertices.index(atom1) < self.vertices.index(atom2) and bond.isSingle() and not self.isBondInCycle(atom1, atom2): - if len(self.edges[atom1]) > 1 and len(self.edges[atom2]) > 1: + for atom1 in self.vertices: + for atom2, bond in atom1.edges.items(): + if self.vertices.index(atom1) < self.vertices.index(atom2) and bond.isSingle() and not self.isBondInCycle(bond): + if len(atom1.edges) > 1 and len(atom2.edges) > 1: count += 1 return count @@ -1154,12 +1157,12 @@ def findAllDelocalizationPaths(self, atom1): # Find all delocalization paths paths = [] - for atom2, bond12 in self.edges[atom1].iteritems(): + for atom2, bond12 in atom1.edges.items(): # Vinyl bond must be capable of gaining an order if bond12.order in ['S', 'D']: - for atom3, bond23 in self.getBonds(atom2).iteritems(): + for atom3, bond23 in atom2.edges.items(): # Allyl bond must be capable of losing an order without breaking - if atom1 is not atom3 and bond23.order in ['D', 'T']: + if atom1 is not atom3 and (bond23.order == 'D' or bond23.order == 'T'): paths.append([atom1, atom2, atom3, bond12, bond23]) return paths diff --git a/unittest/molecule/moleculeTest.py b/unittest/molecule/moleculeTest.py index e73d3a74d1..945678fd79 100644 --- a/unittest/molecule/moleculeTest.py +++ b/unittest/molecule/moleculeTest.py @@ -246,19 +246,6 @@ def testPickle(self): self.assertEqual(self.atom.charge, atom.charge) self.assertEqual(self.atom.label, atom.label) - def testOutput(self): - """ - Test that we can reconstruct a Atom object from its repr() - output with no loss of information. - """ - exec('atom = {0!r}'.format(self.atom)) - self.assertEqual(self.atom.element.symbol, atom.element.symbol) - self.assertEqual(self.atom.atomType, atom.atomType) - self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons) - self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(self.atom.charge, atom.charge) - self.assertEqual(self.atom.label, atom.label) - ################################################################################ class TestBond(unittest.TestCase): @@ -270,7 +257,7 @@ def setUp(self): """ A method called before each unit test in this class. """ - self.bond = Bond(order='D') + self.bond = Bond(atom1=None, atom2=None, order='D') self.orderList = ['S','D','T','B'] def testIsSingle(self): @@ -278,7 +265,7 @@ def testIsSingle(self): Test the Bond.isSingle() method. """ for order in self.orderList: - bond = Bond(order=order) + bond = Bond(None, None, order=order) if order == 'S': self.assertTrue(bond.isSingle()) else: @@ -289,7 +276,7 @@ def testIsDouble(self): Test the Bond.isDouble() method. """ for order in self.orderList: - bond = Bond(order=order) + bond = Bond(None, None, order=order) if order == 'D': self.assertTrue(bond.isDouble()) else: @@ -300,7 +287,7 @@ def testIsTriple(self): Test the Bond.isTriple() method. """ for order in self.orderList: - bond = Bond(order=order) + bond = Bond(None, None, order=order) if order == 'T': self.assertTrue(bond.isTriple()) else: @@ -311,7 +298,7 @@ def testIsBenzene(self): Test the Bond.isBenzene() method. """ for order in self.orderList: - bond = Bond(order=order) + bond = Bond(None, None, order=order) if order == 'B': self.assertTrue(bond.isBenzene()) else: @@ -322,7 +309,7 @@ def testIncrementOrder(self): Test the Bond.incrementOrder() method. """ for order in self.orderList: - bond = Bond(order=order) + bond = Bond(None, None, order=order) try: bond.incrementOrder() if order == 'S': @@ -337,7 +324,7 @@ def testDecrementOrder(self): Test the Bond.decrementOrder() method. """ for order in self.orderList: - bond = Bond(order=order) + bond = Bond(None, None, order=order) try: bond.decrementOrder() if order == 'D': @@ -353,7 +340,7 @@ def testApplyActionBreakBond(self): """ action = ['BREAK_BOND', '*1', 'S', '*2'] for order0 in self.orderList: - bond0 = Bond(order=order0) + bond0 = Bond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -367,7 +354,7 @@ def testApplyActionFormBond(self): """ action = ['FORM_BOND', '*1', 'S', '*2'] for order0 in self.orderList: - bond0 = Bond(order=order0) + bond0 = Bond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -381,7 +368,7 @@ def testApplyActionIncrementBond(self): """ action = ['CHANGE_BOND', '*1', 1, '*2'] for order0 in self.orderList: - bond0 = Bond(order=order0) + bond0 = Bond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -394,7 +381,7 @@ def testApplyActionDecrementBond(self): """ action = ['CHANGE_BOND', '*1', -1, '*2'] for order0 in self.orderList: - bond0 = Bond(order=order0) + bond0 = Bond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -407,7 +394,7 @@ def testApplyActionGainRadical(self): """ action = ['GAIN_RADICAL', '*1', 1] for order0 in self.orderList: - bond0 = Bond(order=order0) + bond0 = Bond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -421,7 +408,7 @@ def testApplyActionLoseRadical(self): """ action = ['LOSE_RADICAL', '*1', 1] for order0 in self.orderList: - bond0 = Bond(order=order0) + bond0 = Bond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -435,8 +422,8 @@ def testEquivalent(self): """ for order1 in self.orderList: for order2 in self.orderList: - bond1 = Bond(order=order1) - bond2 = Bond(order=order2) + bond1 = Bond(None, None, order=order1) + bond2 = Bond(None, None, order=order2) if order1 == order2: self.assertTrue(bond1.equivalent(bond2)) self.assertTrue(bond2.equivalent(bond1)) @@ -450,8 +437,8 @@ def testIsSpecificCaseOf(self): """ for order1 in self.orderList: for order2 in self.orderList: - bond1 = Bond(order=order1) - bond2 = Bond(order=order2) + bond1 = Bond(None, None, order=order1) + bond2 = Bond(None, None, order=order2) if order1 == order2: self.assertTrue(bond1.isSpecificCaseOf(bond2)) else: @@ -473,14 +460,6 @@ def testPickle(self): bond = cPickle.loads(cPickle.dumps(self.bond)) self.assertEqual(self.bond.order, bond.order) - def testOutput(self): - """ - Test that we can reconstruct a Bond object from its repr() - output with no loss of information. - """ - exec('bond = {0!r}'.format(self.bond)) - self.assertEqual(self.bond.order, bond.order) - ################################################################################ class TestMolecule(unittest.TestCase): @@ -490,9 +469,9 @@ class TestMolecule(unittest.TestCase): def setUp(self): self.adjlist = """ -1 *2 C 1 {2,D} {3,S} -2 *1 O 0 {1,D} -3 C 0 {1,S} +1 *2 C 1 {2,D} {3,S} +2 *1 O 0 {1,D} +3 C 0 {1,S} """ self.molecule = Molecule().fromAdjacencyList(self.adjlist) @@ -562,8 +541,8 @@ def testFromAdjacencyList(self): self.assertTrue(self.molecule.hasBond(atom1,atom2)) self.assertTrue(self.molecule.hasBond(atom1,atom3)) self.assertFalse(self.molecule.hasBond(atom2,atom3)) - bond12 = self.molecule.bonds[atom1][atom2] - bond13 = self.molecule.bonds[atom1][atom3] + bond12 = atom1.bonds[atom2] + bond13 = atom1.bonds[atom3] self.assertTrue(atom1.label == '*2') self.assertTrue(atom1.element.symbol == 'C') @@ -613,8 +592,7 @@ def testSubgraphIsomorphism(self): """) self.assertTrue(molecule.isSubgraphIsomorphic(group)) - match, mapping = molecule.findSubgraphIsomorphisms(group) - self.assertTrue(match) + mapping = molecule.findSubgraphIsomorphisms(group) self.assertTrue(len(mapping) == 4, "len(mapping) = %d, should be = 4" % (len(mapping))) for map in mapping: self.assertTrue(len(map) == min(len(molecule.atoms), len(group.atoms))) @@ -658,8 +636,7 @@ def testSubgraphIsomorphismAgain(self): self.assertTrue(molecule.isSubgraphIsomorphic(group, initialMap)) initialMap = {labeled1: labeled2} - match, mapping = molecule.findSubgraphIsomorphisms(group, initialMap) - self.assertTrue(match) + mapping = molecule.findSubgraphIsomorphisms(group, initialMap) self.assertTrue(len(mapping) == 2, "len(mapping) = %d, should be = 2" % (len(mapping))) for map in mapping: self.assertTrue(len(map) == min(len(molecule.atoms), len(group.atoms))) @@ -689,8 +666,7 @@ def testSubgraphIsomorphismManyLabels(self): initialMap[atom1] = labeled2[label] self.assertTrue(molecule.isSubgraphIsomorphic(group, initialMap)) - match, mapping = molecule.findSubgraphIsomorphisms(group, initialMap) - self.assertTrue(match) + mapping = molecule.findSubgraphIsomorphisms(group, initialMap) self.assertTrue(len(mapping) == 1) for map in mapping: self.assertTrue(len(map) == min(len(molecule.atoms), len(group.atoms))) @@ -732,9 +708,9 @@ def testIsInCycleEthane(self): molecule = Molecule().fromSMILES('CC') for atom in molecule.atoms: self.assertFalse(molecule.isAtomInCycle(atom)) - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: - self.assertFalse(molecule.isBondInCycle(atom1, atom2)) + for atom1 in molecule.atoms: + for atom2, bond in atom1.bonds.items(): + self.assertFalse(molecule.isBondInCycle(bond)) def testIsInCycleCyclohexane(self): """ @@ -746,12 +722,12 @@ def testIsInCycleCyclohexane(self): self.assertFalse(molecule.isAtomInCycle(atom)) elif atom.isCarbon(): self.assertTrue(molecule.isAtomInCycle(atom)) - for atom1 in molecule.bonds: - for atom2 in molecule.bonds[atom1]: + for atom1 in molecule.atoms: + for atom2, bond in atom1.bonds.items(): if atom1.isCarbon() and atom2.isCarbon(): - self.assertTrue(molecule.isBondInCycle(atom1, atom2)) + self.assertTrue(molecule.isBondInCycle(bond)) else: - self.assertFalse(molecule.isBondInCycle(atom1, atom2)) + self.assertFalse(molecule.isBondInCycle(bond)) def testFromSMILESH(self): """ From 744b10ca1e9fa62f072dfd62d22fba53e808be64 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Tue, 5 Jun 2012 17:24:29 -0400 Subject: [PATCH 15/23] Updated GroupAtom, GroupBond, and Group to work with new Graph format. As before, this is mostly updating to reflect that the edges are now stored on the GroupAtom objects instead of the Group objects. Some of the corresponding unit tests also needed a bit of work. --- rmgpy/molecule/group.pxd | 16 ++++---- rmgpy/molecule/group.py | 63 ++++++++++++++++++----------- unittest/molecule/groupTest.py | 74 ++++++++++++---------------------- 3 files changed, 72 insertions(+), 81 deletions(-) diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd index c5d8623952..fe38fe0626 100644 --- a/rmgpy/molecule/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -36,7 +36,7 @@ cdef class GroupAtom(Vertex): cdef public list charge cdef public str label - cpdef copy(self) + cpdef Vertex copy(self) cpdef __changeBond(self, short order) @@ -60,7 +60,7 @@ cdef class GroupBond(Edge): cdef public list order - cpdef copy(self) + cpdef Edge copy(self) cpdef __changeBond(self, short order) @@ -76,7 +76,7 @@ cdef class Group(Graph): cpdef addAtom(self, GroupAtom atom) - cpdef addBond(self, GroupAtom atom1, GroupAtom atom2, GroupBond bond) + cpdef addBond(self, GroupBond bond) cpdef dict getBonds(self, GroupAtom atom) @@ -88,7 +88,7 @@ cdef class Group(Graph): cpdef removeAtom(self, GroupAtom atom) - cpdef removeBond(self, GroupAtom atom1, GroupAtom GroupAtom2) + cpdef removeBond(self, GroupBond bond) cpdef sortAtoms(self) @@ -106,10 +106,10 @@ cdef class Group(Graph): cpdef toAdjacencyList(self, str label=?) - cpdef bint isIsomorphic(self, Graph other, dict initialMap=?) + cpdef bint isIsomorphic(self, Graph other, dict initialMap=?) except -2 - cpdef tuple findIsomorphism(self, Graph other, dict initialMap=?) + cpdef list findIsomorphism(self, Graph other, dict initialMap=?) - cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?) + cpdef bint isSubgraphIsomorphic(self, Graph other, dict initialMap=?) except -2 - cpdef tuple findSubgraphIsomorphisms(self, Graph other, dict initialMap=?) + cpdef list findSubgraphIsomorphisms(self, Graph other, dict initialMap=?) diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index 7aaf218c3d..27f258bcc9 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -88,23 +88,42 @@ def __reduce__(self): """ A helper function used when pickling an object. """ + d = { + 'edges': self.edges, + 'connectivity1': self.connectivity1, + 'connectivity2': self.connectivity2, + 'connectivity3': self.connectivity3, + 'sortingLabel': self.sortingLabel, + } atomType = self.atomType if atomType is not None: atomType = [a.label for a in atomType] - return (GroupAtom, (atomType, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label)) + return (GroupAtom, (atomType, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label), d) + + def __setstate__(self, d): + """ + A helper function used when unpickling an object. + """ + self.edges = d['edges'] + self.connectivity1 = d['connectivity1'] + self.connectivity2 = d['connectivity2'] + self.connectivity3 = d['connectivity3'] + self.sortingLabel = d['sortingLabel'] def __str__(self): """ Return a human-readable string representation of the object. """ - return "".format(self.atomType) + return '[{0}]'.format(','.join([repr(a.label) for a in self.atomType])) def __repr__(self): """ Return a representation that can be used to reconstruct the object. """ - atomType = ','.join(['"{0}"'.format(a.label) for a in self.atomType]) - return "GroupAtom(atomType=[{0}], radicalElectrons={1}, spinMultiplicity={2}, charge={3}, label='{4}')".format(atomType, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label) + return "".format(self) + + @property + def bonds(self): return self.edges def copy(self): """ @@ -308,34 +327,34 @@ class GroupBond(Edge): group if it matches *any* item in the list. """ - def __init__(self, order=None): - Edge.__init__(self) + def __init__(self, atom1, atom2, order=None): + Edge.__init__(self, atom1, atom2) self.order = order or [] def __str__(self): """ Return a human-readable string representation of the object. """ - return "".format(self.order) + return str(self.order) def __repr__(self): """ Return a representation that can be used to reconstruct the object. """ - return "GroupBond(order={0})".format(self.order) + return "".format(self.order) def __reduce__(self): """ A helper function used when pickling an object. """ - return (GroupBond, (self.order,)) + return (GroupBond, (self.vertex1, self.vertex2, self.order)) def copy(self): """ Return a deep copy of the :class:`GroupBond` object. Modifying the attributes of the copy will not affect the original. """ - return GroupBond(self.order[:]) + return GroupBond(self.vertex1, self.vertex2, self.order[:]) def __changeBond(self, order): """ @@ -434,36 +453,32 @@ class Group(Graph): Corresponding alias methods have also been provided. """ - def __init__(self, atoms=None, bonds=None): - Graph.__init__(self, atoms, bonds) + def __init__(self, atoms=None): + Graph.__init__(self, atoms) self.updateConnectivityValues() def __reduce__(self): """ A helper function used when pickling an object. """ - return (Group, (self.vertices, self.edges)) + return (Group, (self.vertices,)) def __getAtoms(self): return self.vertices def __setAtoms(self, atoms): self.vertices = atoms atoms = property(__getAtoms, __setAtoms) - def __getBonds(self): return self.edges - def __setBonds(self, bonds): self.edges = bonds - bonds = property(__getBonds, __setBonds) - def addAtom(self, atom): """ Add an `atom` to the graph. The atom is initialized with no bonds. """ return self.addVertex(atom) - def addBond(self, atom1, atom2, bond): + def addBond(self, bond): """ Add a `bond` to the graph as an edge connecting the two atoms `atom1` and `atom2`. """ - return self.addEdge(atom1, atom2, bond) + return self.addEdge(bond) def getBonds(self, atom): """ @@ -499,13 +514,13 @@ def removeAtom(self, atom): """ return self.removeVertex(atom) - def removeBond(self, atom1, atom2): + def removeBond(self, bond): """ Remove the bond between atoms `atom1` and `atom2` from the graph. Does not remove atoms that no longer have any bonds as a result of this removal. """ - return self.removeEdge(atom1, atom2) + return self.removeEdge(bond) def sortAtoms(self): """ @@ -523,7 +538,7 @@ def copy(self, deep=False): """ other = cython.declare(Group) g = Graph.copy(self, deep) - other = Group(g.vertices, g.edges) + other = Group(g.vertices) return other def merge(self, other): @@ -533,7 +548,7 @@ def merge(self, other): object is returned. """ g = Graph.merge(self, other) - molecule = Group(atoms=g.vertices, bonds=g.edges) + molecule = Group(atoms=g.vertices) return molecule def split(self): @@ -544,7 +559,7 @@ def split(self): graphs = Graph.split(self) molecules = [] for g in graphs: - molecule = Group(atoms=g.vertices, bonds=g.edges) + molecule = Group(atoms=g.vertices) molecules.append(molecule) return molecules diff --git a/unittest/molecule/groupTest.py b/unittest/molecule/groupTest.py index d5140c86cc..da45a4fad2 100644 --- a/unittest/molecule/groupTest.py +++ b/unittest/molecule/groupTest.py @@ -192,19 +192,6 @@ def testPickle(self): self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity) self.assertEqual(self.atom.charge, atom.charge) self.assertEqual(self.atom.label, atom.label) - - def testOutput(self): - """ - Test that we can reconstruct a GroupAtom object from its repr() - output with no loss of information. - """ - exec('atom = {0!r}'.format(self.atom)) - self.assertEqual(len(self.atom.atomType), len(atom.atomType)) - self.assertEqual(self.atom.atomType[0].label, atom.atomType[0].label) - self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons) - self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity) - self.assertEqual(self.atom.charge, atom.charge) - self.assertEqual(self.atom.label, atom.label) ################################################################################ @@ -217,7 +204,7 @@ def setUp(self): """ A method called before each unit test in this class. """ - self.bond = GroupBond(order=['D']) + self.bond = GroupBond(None, None, order=['D']) self.orderList = [['S'], ['D'], ['T'], ['B'], ['S','D'], ['D','S'], ['D','T'], ['S','D','T']] def testApplyActionBreakBond(self): @@ -226,7 +213,7 @@ def testApplyActionBreakBond(self): """ action = ['BREAK_BOND', '*1', 'S', '*2'] for order0 in self.orderList: - bond0 = GroupBond(order=order0) + bond0 = GroupBond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -240,7 +227,7 @@ def testApplyActionFormBond(self): """ action = ['FORM_BOND', '*1', 'S', '*2'] for order0 in self.orderList: - bond0 = GroupBond(order=order0) + bond0 = GroupBond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -254,7 +241,7 @@ def testApplyActionIncrementBond(self): """ action = ['CHANGE_BOND', '*1', 1, '*2'] for order0 in self.orderList: - bond0 = GroupBond(order=order0) + bond0 = GroupBond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -267,7 +254,7 @@ def testApplyActionDecrementBond(self): """ action = ['CHANGE_BOND', '*1', -1, '*2'] for order0 in self.orderList: - bond0 = GroupBond(order=order0) + bond0 = GroupBond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -280,7 +267,7 @@ def testApplyActionGainRadical(self): """ action = ['GAIN_RADICAL', '*1', 1] for order0 in self.orderList: - bond0 = GroupBond(order=order0) + bond0 = GroupBond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -294,7 +281,7 @@ def testApplyActionLoseRadical(self): """ action = ['LOSE_RADICAL', '*1', 1] for order0 in self.orderList: - bond0 = GroupBond(order=order0) + bond0 = GroupBond(None, None, order=order0) bond = bond0.copy() try: bond.applyAction(action) @@ -308,8 +295,8 @@ def testEquivalent(self): """ for order1 in self.orderList: for order2 in self.orderList: - bond1 = GroupBond(order=order1) - bond2 = GroupBond(order=order2) + bond1 = GroupBond(None, None, order=order1) + bond2 = GroupBond(None, None, order=order2) if order1 == order2 or (all([o in order2 for o in order1]) and all([o in order1 for o in order2])): self.assertTrue(bond1.equivalent(bond2)) self.assertTrue(bond2.equivalent(bond1)) @@ -323,8 +310,8 @@ def testIsSpecificCaseOf(self): """ for order1 in self.orderList: for order2 in self.orderList: - bond1 = GroupBond(order=order1) - bond2 = GroupBond(order=order2) + bond1 = GroupBond(None, None, order=order1) + bond2 = GroupBond(None, None, order=order2) if order1 == order2 or all([o in order2 for o in order1]): self.assertTrue(bond1.isSpecificCaseOf(bond2)) else: @@ -347,15 +334,6 @@ def testPickle(self): bond = cPickle.loads(cPickle.dumps(self.bond)) self.assertEqual(len(self.bond.order), len(bond.order)) self.assertEqual(self.bond.order, bond.order) - - def testOutput(self): - """ - Test that we can reconstruct a GroupBond object from its repr() - output with no loss of information. - """ - exec('bond = {0!r}'.format(self.bond)) - self.assertEqual(len(self.bond.order), len(bond.order)) - self.assertEqual(self.bond.order, bond.order) ################################################################################ @@ -366,9 +344,9 @@ class TestGroup(unittest.TestCase): def setUp(self): self.adjlist = """ -1 *2 {Cs,Cd} 0 {2,{S,D}} {3,S} -2 *1 {Os,Od} 0 {1,{S,D}} -3 R!H 0 {1,S} +1 *2 {Cs,Cd} 0 {2,{S,D}} {3,S} +2 *1 {Os,Od} 0 {1,{S,D}} +3 R!H 0 {1,S} """ self.group = Group().fromAdjacencyList(self.adjlist) @@ -426,8 +404,8 @@ def testFromAdjacencyList(self): self.assertTrue(self.group.hasBond(atom1,atom2)) self.assertTrue(self.group.hasBond(atom1,atom3)) self.assertFalse(self.group.hasBond(atom2,atom3)) - bond12 = self.group.bonds[atom1][atom2] - bond13 = self.group.bonds[atom1][atom3] + bond12 = atom1.bonds[atom2] + bond13 = atom1.bonds[atom3] self.assertTrue(atom1.label == '*2') self.assertTrue(atom1.atomType[0].label in ['Cs','Cd']) @@ -477,18 +455,17 @@ def testFindIsomorphism(self): """ group = Group().fromAdjacencyList(adjlist) result = self.group.findIsomorphism(group) - self.assertTrue(result[0]) - self.assertEqual(len(result[1]), 1) - for atom1, atom2 in result[1][0].iteritems(): + self.assertEqual(len(result), 1) + for atom1, atom2 in result[0].items(): self.assertTrue(atom1 in self.group.atoms) self.assertTrue(atom2 in group.atoms) self.assertTrue(atom1.equivalent(atom2)) - for atom3 in self.group.bonds[atom1]: - atom4 = result[1][0][atom3] - self.assertTrue(atom4 in group.bonds[atom2]) + for atom3 in atom1.bonds: + atom4 = result[0][atom3] + self.assertTrue(atom4 in atom2.bonds) self.assertTrue(atom3.equivalent(atom4)) - bond1 = self.group.bonds[atom1][atom3] - bond2 = group.bonds[atom2][atom4] + bond1 = atom1.bonds[atom3] + bond2 = atom2.bonds[atom4] self.assertTrue(bond1.equivalent(bond2)) def testIsSubgraphIsomorphic(self): @@ -511,9 +488,8 @@ def testFindSubgraphIsomorphisms(self): """ group = Group().fromAdjacencyList(adjlist) result = self.group.findSubgraphIsomorphisms(group) - self.assertTrue(result[0]) - self.assertEqual(len(result[1]), 1) - for atom1, atom2 in result[1][0].iteritems(): + self.assertEqual(len(result), 1) + for atom1, atom2 in result[0].iteritems(): self.assertTrue(atom1 in self.group.atoms) self.assertTrue(atom2 in group.atoms) self.assertTrue(atom1.equivalent(atom2)) From 0581d153701651f1da322f16256dff1afe59cc5d Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Wed, 6 Jun 2012 13:48:11 -0400 Subject: [PATCH 16/23] Updated molecule drawing functionality to work with new Graph format. As before, this is mostly updating to reflect that the edges are now stored on the Atom objects instead of the Molecule objects. --- rmgpy/molecule/molecule_draw.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/rmgpy/molecule/molecule_draw.py b/rmgpy/molecule/molecule_draw.py index 0364473570..225732d3b3 100644 --- a/rmgpy/molecule/molecule_draw.py +++ b/rmgpy/molecule/molecule_draw.py @@ -543,7 +543,7 @@ def findLongestPath(chemGraph, atoms0): """ atom1 = atoms0[-1] paths = [atoms0] - for atom2 in chemGraph.bonds[atom1]: + for atom2 in atom1.bonds: if atom2 not in atoms0: atoms = atoms0[:] atoms.append(atom2) @@ -588,7 +588,7 @@ def findBackbone(chemGraph, ringSystems): # Find the terminal atoms - those that only have one explicit bond terminalAtoms = [] for atom in chemGraph.atoms: - if len(chemGraph.bonds[atom]) == 1: + if len(atom.bonds) == 1: terminalAtoms.append(atom) # Starting from each terminal atom, find the longest straight path to @@ -1132,6 +1132,8 @@ def drawMolecule(molecule, path=None, surface=''): except ImportError: print 'Cairo not found; molecule will not be drawn.' return + + molecule = molecule.copy(deep=True) # This algorithm now works with explicit hydrogen atoms on the molecule. # Please ensure all the subroutines do also. @@ -1144,8 +1146,8 @@ def drawMolecule(molecule, path=None, surface=''): # bonds = molecule.bonds.copy() is too shallow for a dict-of-dicts, # so we loop one level deep and copy the inner dicts. bonds = dict() - for atom1,atom2dict in molecule.bonds.iteritems(): - bonds[atom1] = atom2dict.copy() + for atom1 in molecule.atoms: + bonds[atom1] = atom1.edges.copy() cycles = molecule.getSmallestSetOfSmallestRings() From 8a9cb44aa6132f064b65f3f9fdb7e1fea2a8f577 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Wed, 6 Jun 2012 13:53:30 -0400 Subject: [PATCH 17/23] Updated thermo and kinetics database code to work with new Graph format. Mostly this is creating Bond and GroupBond objects with the two atoms specified in the __init__() method. The findIsomorphism() and findSubgraphIsomorphisms() methods also now only return the list of mappings; this caused a few additional changes in the database code. --- rmgpy/data/kinetics.py | 82 ++++++++++++++++++++---------------------- rmgpy/data/thermo.py | 8 ++--- 2 files changed, 43 insertions(+), 47 deletions(-) diff --git a/rmgpy/data/kinetics.py b/rmgpy/data/kinetics.py index 07eabf3216..d0a1566ab9 100644 --- a/rmgpy/data/kinetics.py +++ b/rmgpy/data/kinetics.py @@ -238,15 +238,15 @@ def __apply(self, struct, doForward, unique): atom2.applyAction(['CHANGE_BOND', label1, -info, label2]) bond.applyAction(['CHANGE_BOND', label1, -info, label2]) elif (action[0] == 'FORM_BOND' and doForward) or (action[0] == 'BREAK_BOND' and not doForward): - bond = GroupBond(order=['S']) if pattern else Bond(order='S') - struct.addBond(atom1, atom2, bond) + bond = GroupBond(atom1, atom2, order=['S']) if pattern else Bond(atom1, atom2, order='S') + struct.addBond(bond) atom1.applyAction(['FORM_BOND', label1, info, label2]) atom2.applyAction(['FORM_BOND', label1, info, label2]) elif (action[0] == 'BREAK_BOND' and doForward) or (action[0] == 'FORM_BOND' and not doForward): if not struct.hasBond(atom1, atom2): raise InvalidActionError('Attempted to remove a nonexistent bond.') bond = struct.getBond(atom1, atom2) - struct.removeBond(atom1, atom2) + struct.removeBond(bond) atom1.applyAction(['BREAK_BOND', label1, info, label2]) atom2.applyAction(['BREAK_BOND', label1, info, label2]) @@ -1980,7 +1980,7 @@ def addKineticsRulesFromTrainingSet(self, thermoDatabase=None): data.changeT0(1) # Estimate the thermo for the reactants and products - item = deepcopy(entry.item) + item = Reaction(reactants=[m.copy(deep=True) for m in entry.item.reactants], products=[m.copy(deep=True) for m in entry.item.products]) item.reactants = [Species(molecule=[m]) for m in item.reactants] for reactant in item.reactants: reactant.generateResonanceIsomers() @@ -2355,9 +2355,8 @@ def __matchReactantToTemplate(self, reactant, templateReactant): if isinstance(struct, LogicNode): mappings = [] for child_structure in struct.getPossibleStructures(self.groups.entries): - ismatch, map = reactant.findSubgraphIsomorphisms(child_structure) - if ismatch: mappings.extend(map) - return len(mappings) > 0, mappings + mappings.extend(reactant.findSubgraphIsomorphisms(child_structure)) + return mappings elif isinstance(struct, Group): return reactant.findSubgraphIsomorphisms(struct) @@ -2505,18 +2504,17 @@ def __generateReactions(self, reactants, forward=True): # Iterate over all resonance isomers of the reactant for molecule in reactants[0]: - ismatch, mappings = self.__matchReactantToTemplate(molecule, template.reactants[0]) - if ismatch: - for map in mappings: - reactantStructures = [molecule] - try: - productStructures = self.__generateProductStructures(reactantStructures, [map], forward) - except ForbiddenStructureException: - pass - else: - if productStructures is not None: - rxn = self.__createReaction(reactantStructures, productStructures, forward) - if rxn: rxnList.append(rxn) + mappings = self.__matchReactantToTemplate(molecule, template.reactants[0]) + for map in mappings: + reactantStructures = [molecule] + try: + productStructures = self.__generateProductStructures(reactantStructures, [map], forward) + except ForbiddenStructureException: + pass + else: + if productStructures is not None: + rxn = self.__createReaction(reactantStructures, productStructures, forward) + if rxn: rxnList.append(rxn) # Bimolecular reactants: A + B --> products elif len(reactants) == 2 and len(template.reactants) == 2: @@ -2529,11 +2527,30 @@ def __generateReactions(self, reactants, forward=True): for moleculeB in moleculesB: # Reactants stored as A + B - ismatchA, mappingsA = self.__matchReactantToTemplate(moleculeA, template.reactants[0]) - ismatchB, mappingsB = self.__matchReactantToTemplate(moleculeB, template.reactants[1]) + mappingsA = self.__matchReactantToTemplate(moleculeA, template.reactants[0]) + mappingsB = self.__matchReactantToTemplate(moleculeB, template.reactants[1]) # Iterate over each pair of matches (A, B) - if ismatchA and ismatchB: + for mapA in mappingsA: + for mapB in mappingsB: + reactantStructures = [moleculeA, moleculeB] + try: + productStructures = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward) + except ForbiddenStructureException: + pass + else: + if productStructures is not None: + rxn = self.__createReaction(reactantStructures, productStructures, forward) + if rxn: rxnList.append(rxn) + + # Only check for swapped reactants if they are different + if reactants[0] is not reactants[1]: + + # Reactants stored as B + A + mappingsA = self.__matchReactantToTemplate(moleculeA, template.reactants[1]) + mappingsB = self.__matchReactantToTemplate(moleculeB, template.reactants[0]) + + # Iterate over each pair of matches (A, B) for mapA in mappingsA: for mapB in mappingsB: reactantStructures = [moleculeA, moleculeB] @@ -2546,27 +2563,6 @@ def __generateReactions(self, reactants, forward=True): rxn = self.__createReaction(reactantStructures, productStructures, forward) if rxn: rxnList.append(rxn) - # Only check for swapped reactants if they are different - if reactants[0] is not reactants[1]: - - # Reactants stored as B + A - ismatchA, mappingsA = self.__matchReactantToTemplate(moleculeA, template.reactants[1]) - ismatchB, mappingsB = self.__matchReactantToTemplate(moleculeB, template.reactants[0]) - - # Iterate over each pair of matches (A, B) - if ismatchA and ismatchB: - for mapA in mappingsA: - for mapB in mappingsB: - reactantStructures = [moleculeA, moleculeB] - try: - productStructures = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward) - except ForbiddenStructureException: - pass - else: - if productStructures is not None: - rxn = self.__createReaction(reactantStructures, productStructures, forward) - if rxn: rxnList.append(rxn) - # Remove duplicates from the reaction list index0 = 0 while index0 < len(rxnList): diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 394e6d90ae..0feacbb50a 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -646,9 +646,9 @@ def estimateThermoViaGroupAdditivity(self, molecule): for atom in saturatedStruct.atoms: for i in range(atom.radicalElectrons): H = Atom('H') - bond = Bond('S') + bond = Bond(atom, H, 'S') saturatedStruct.addAtom(H) - saturatedStruct.addBond(atom, H, bond) + saturatedStruct.addBond(bond) if atom not in added: added[atom] = [] added[atom].append([H, bond]) @@ -675,7 +675,7 @@ def estimateThermoViaGroupAdditivity(self, molecule): # Remove the added hydrogen atoms and bond and restore the radical for H, bond in added[atom]: - saturatedStruct.removeBond(atom, H) + saturatedStruct.removeBond(bond) saturatedStruct.removeAtom(H) atom.incrementRadical() @@ -692,7 +692,7 @@ def estimateThermoViaGroupAdditivity(self, molecule): # Re-saturate for H, bond in added[atom]: saturatedStruct.addAtom(H) - saturatedStruct.addBond(atom, H, bond) + saturatedStruct.addBond(bond) atom.decrementRadical() # Subtract the enthalpy of the added hydrogens From 6bbb693662b4e085178cd5ac6c582dbb25dbad1a Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 14 Jun 2012 13:46:02 -0400 Subject: [PATCH 18/23] Reordered a couple of 'if' blocks to be more efficient. Small changes, but I think improvements. One reduces the number of if checks, the other reduces the number of loops. --- rmgpy/molecule/graph.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/rmgpy/molecule/graph.py b/rmgpy/molecule/graph.py index b46d56a3da..b79a46f912 100644 --- a/rmgpy/molecule/graph.py +++ b/rmgpy/molecule/graph.py @@ -188,12 +188,12 @@ def getOtherVertex(self, vertex): Raise a :class:`ValueError` if the given vertex is not part of the edge. """ - if self.vertex1 is not vertex and self.vertex2 is not vertex: - raise ValueError('The given vertex is not one of the vertices of this edge.') - elif self.vertex1 is vertex: + if self.vertex1 is vertex: return self.vertex2 elif self.vertex2 is vertex: return self.vertex1 + else: + raise ValueError('The given vertex is not one of the vertices of this edge.') ################################################################################ @@ -307,9 +307,9 @@ def copy(self, deep=False): edges = vertex.edges other.addVertex(vertex) vertex.edges = edges - for vertex1 in self.vertices: - for vertex2 in vertex1.edges: - if deep: + if deep: + for vertex1 in self.vertices: + for vertex2 in vertex1.edges: index1 = self.vertices.index(vertex1) index2 = self.vertices.index(vertex2) edge = vertex1.edges[vertex2].copy() From 82b6a2b44d6c86affeb5546a76c490c3273270e2 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 14 Jun 2012 13:46:36 -0400 Subject: [PATCH 19/23] Small changes to docstring and comment. I *think* these are more accurate... --- rmgpy/molecule/graph.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/rmgpy/molecule/graph.py b/rmgpy/molecule/graph.py index b79a46f912..bed5ff5ee7 100644 --- a/rmgpy/molecule/graph.py +++ b/rmgpy/molecule/graph.py @@ -279,7 +279,7 @@ def removeVertex(self, vertex): def removeEdge(self, edge): """ - Remove the edge having vertices `vertex1` and `vertex2` from the graph. + Remove the specified `edge` from the graph. Does not remove vertices that no longer have any edges as a result of this removal. """ @@ -697,7 +697,6 @@ def VF2_isomorphism(graph1, graph2, subgraph=False, findAll=False, initialMappin from `graph1` to those in `subgraph` that are fixed in the isomorphism check; this mapping will appear in every returned mapping. Note that no validation of this initial mapping is performed in this function. - """ cython.declare(vertex1=Vertex, vertex2=Vertex) cython.declare(mappingList=list) @@ -717,7 +716,7 @@ def VF2_isomorphism(graph1, graph2, subgraph=False, findAll=False, initialMappin # a subgraph of the first return list() if findAll else False - # Initialize callDepth with the size of the largest graph + # Initialize callDepth with the size of the smallest graph # Each recursive call to VF2_match will decrease it by one; # when the whole graph has been explored, it should reach 0 # It should never go below zero! From d177753ac8fe55423f6063595a327e709cbe1ad0 Mon Sep 17 00:00:00 2001 From: Josh Allen Date: Thu, 14 Jun 2012 15:01:13 -0400 Subject: [PATCH 20/23] Fixed bug in cycle detection algorithm caused by new graph format. The getSmallestSetOfSmallestRings() method generates a copy of the graph before applying the SSSR algorithm so as to not modify the original. Before, this was safe because we could make a shallow copy and reuse the same Vertex and Edge objects, since we didn't store any information about the graph connectivity on these objects. Now that we are, we must make a deep copy of the graph to use for the SSSR algorithm. As a result, we also need to map the vertices of the copy back to those of the original graph before returning. This was causing the molecule drawing to fail for cyclic species; the problem should now be fixed. --- rmgpy/molecule/graph.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/rmgpy/molecule/graph.py b/rmgpy/molecule/graph.py index bed5ff5ee7..1fd64b1a52 100644 --- a/rmgpy/molecule/graph.py +++ b/rmgpy/molecule/graph.py @@ -558,11 +558,12 @@ def getSmallestSetOfSmallestRings(self): cython.declare(graph=Graph) cython.declare(done=cython.bint, found=cython.bint) cython.declare(cycleList=list, cycles=list, cycle=list, graphs=list, neighbors=list) - cython.declare(verticesToRemove=list) + cython.declare(verticesToRemove=list, vertices=list) cython.declare(vertex=Vertex, rootVertex=Vertex) # Make a copy of the graph so we don't modify the original graph = self.copy(deep=True) + vertices = graph.vertices[:] # Step 1: Remove all terminal vertices done = False @@ -630,6 +631,10 @@ def getSmallestSetOfSmallestRings(self): for vertex in verticesToRemove: graph.removeVertex(vertex) + # Map atoms in cycles back to atoms in original graph + for i in range(len(cycleList)): + cycleList[i] = [self.vertices[vertices.index(v)] for v in cycles[i]] + return cycleList def isMappingValid(self, other, mapping): From acccd69fc2800410b1e44e64f01a60038d6be5bb Mon Sep 17 00:00:00 2001 From: Richard West Date: Tue, 19 Jun 2012 12:31:17 -0400 Subject: [PATCH 21/23] Fix ring correction searching in group additive thermo estimator. The new method of storing bonds on atoms requires this change. Closes #81. NB. Ring corrections are searched using a molecule fragment with only the ring and no adjoining atoms. If a ring correction definition depends on ligands, it will not be found. (This is not a change in this commit, just an observation) --- rmgpy/data/thermo.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 0feacbb50a..0ebfbfc54f 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -732,14 +732,17 @@ def estimateThermoViaGroupAdditivity(self, molecule): # each ring one time; this doesn't work yet rings = molecule.getSmallestSetOfSmallestRings() for ring in rings: - # Make a temporary structure containing only the atoms in the ring + # NB. if any of the ring corrections depend on ligands not in the ring, they will not be found! ringStructure = Molecule() - for atom in ring: ringStructure.addAtom(atom) + newAtoms = dict() + for atom in ring: + newAtoms[atom] = atom.copy() + ringStructure.addAtom(newAtoms[atom]) # (addAtom deletes the atom's bonds) for atom1 in ring: for atom2 in ring: if molecule.hasBond(atom1, atom2): - ringStructure.addBond(atom1, atom2, molecule.getBond(atom1, atom2)) + ringStructure.addBond(Bond(newAtoms[atom1], newAtoms[atom2], atom1.bonds[atom2].order )) # Get thermo correction for this ring try: From 108886f4d93a1af578c39db35e91836e3649ef79 Mon Sep 17 00:00:00 2001 From: Richard West Date: Tue, 19 Jun 2012 13:37:53 -0400 Subject: [PATCH 22/23] Fix issue #82 in rmgpy.molecule.graph.Graph.getSmallestSetOfSmallestRings Closes #82. --- rmgpy/molecule/graph.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/molecule/graph.py b/rmgpy/molecule/graph.py index 1fd64b1a52..584c489ab6 100644 --- a/rmgpy/molecule/graph.py +++ b/rmgpy/molecule/graph.py @@ -633,7 +633,7 @@ def getSmallestSetOfSmallestRings(self): # Map atoms in cycles back to atoms in original graph for i in range(len(cycleList)): - cycleList[i] = [self.vertices[vertices.index(v)] for v in cycles[i]] + cycleList[i] = [self.vertices[vertices.index(v)] for v in cycleList[i]] return cycleList From 0bd9b13baa960406dc583b573380461bfe0cb8bb Mon Sep 17 00:00:00 2001 From: Richard West Date: Tue, 19 Jun 2012 19:53:13 -0400 Subject: [PATCH 23/23] Change SMILES writer to remove stereochemistry. We (currently) don't have any stereochemistry, it's just guessed by OpenBabel, and makes the SMILES strings longer and more ugly than they need to be. This removes the @@H signs, and also might make it faster as we make fewer calls through pybel. --- rmgpy/molecule/molecule.pxd | 2 ++ rmgpy/molecule/molecule.py | 20 ++++++++++---------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 1e640c465a..9a8cc9a2e8 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -56,6 +56,8 @@ cdef class Atom(Vertex): ################################################################################ +cpdef object SMILEwriter + cdef class Bond(Edge): cdef public str order diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 01e28baa75..9f10233524 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -39,10 +39,11 @@ import os import re import element as elements +import openbabel from .graph import Vertex, Edge, Graph from .group import GroupAtom, GroupBond, Group, ActionError from .atomtype import AtomType, atomTypes, getAtomType - + ################################################################################ class Atom(Vertex): @@ -439,7 +440,10 @@ def applyAction(self, action): raise ActionError('Unable to update GroupBond: Invalid action {0}.'.format(action)) ################################################################################ - +SMILEwriter = openbabel.OBConversion() +SMILEwriter.SetOutFormat('smi') +SMILEwriter.SetOptions("i",SMILEwriter.OUTOPTIONS) # turn off isomer and stereochemistry information (the @ signs!) + class Molecule(Graph): """ A representation of a molecular structure using a graph data type, extending @@ -907,7 +911,6 @@ def toInChI(self): Convert a molecular structure to an InChI string. Uses `OpenBabel `_ to perform the conversion. """ - import openbabel # This version does not write a warning to stderr if stereochemistry is undefined obmol = self.toOBMol() obConversion = openbabel.OBConversion() @@ -959,24 +962,21 @@ def toAugmentedInChIKey(self): return key+'mult'+str(radicalNumber+1) else: return key - + + def toSMILES(self): """ Convert a molecular structure to an SMILES string. Uses `OpenBabel `_ to perform the conversion. """ - import pybel - mol = pybel.Molecule(self.toOBMol()) - return mol.write('smiles').strip() + mol = self.toOBMol() + return SMILEwriter.WriteString(mol).strip() def toOBMol(self): """ Convert a molecular structure to an OpenBabel OBMol object. Uses `OpenBabel `_ to perform the conversion. """ - - import openbabel - cython.declare(atom=Atom, atom1=Atom, bonds=dict, atom2=Atom, bond=Bond) cython.declare(index1=cython.int, index2=cython.int, order=cython.int)