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/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..d0a1566ab9 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 ################################################################################ @@ -239,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]) @@ -1981,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() @@ -2356,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) @@ -2492,8 +2490,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 @@ -2508,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: @@ -2532,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] @@ -2549,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/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..0ebfbfc54f 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 ################################################################################ @@ -607,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() @@ -624,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) @@ -639,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 @@ -660,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]) @@ -689,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() @@ -706,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 @@ -746,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: @@ -768,8 +757,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/molecule/__init__.py b/rmgpy/molecule/__init__.py new file mode 100644 index 0000000000..f355d0cf1d --- /dev/null +++ b/rmgpy/molecule/__init__.py @@ -0,0 +1,36 @@ +#!/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 .adjlist import * +from .atomtype import * +from .element import * +from .molecule import * +from .group import * +from .molecule_draw import drawMolecule 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/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 62% rename from rmgpy/graph.pxd rename to rmgpy/molecule/graph.pxd index 2dc0db7b7c..b90fe675b5 100644 --- a/rmgpy/graph.pxd +++ b/rmgpy/molecule/graph.pxd @@ -26,14 +26,21 @@ 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 Vertex copy(self) - 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 resetConnectivityValues(self) @@ -45,32 +52,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=?) @@ -84,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) @@ -106,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/graph.py b/rmgpy/molecule/graph.py similarity index 52% rename from rmgpy/graph.py rename to rmgpy/molecule/graph.py index ff5cdb42e5..584c489ab6 100644 --- a/rmgpy/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,18 +66,33 @@ 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, + 'terminal': self.terminal, + 'mapping': self.mapping, } return (Vertex, (), d) def __setstate__(self, d): + self.edges = d['edges'] self.connectivity1 = d['connectivity1'] 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): """ @@ -102,6 +118,8 @@ def resetConnectivityValues(self): self.connectivity2 = -1 self.connectivity3 = -1 self.sortingLabel = -1 + self.terminal = False + self.mapping = None def getVertexConnectivityValue(vertex): """ @@ -128,14 +146,25 @@ 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 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): """ @@ -153,6 +182,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 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.') + ################################################################################ class Graph: @@ -166,44 +208,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 +263,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 +271,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. + Remove the specified `edge` 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 +293,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) - for vertex1 in self.vertices: - for vertex2 in self.edges[vertex1]: - if deep: + if deep: + other.addVertex(vertex.copy()) + else: + edges = vertex.edges + other.addVertex(vertex) + vertex.edges = edges + if deep: + for vertex1 in self.vertices: + for vertex2 in vertex1.edges: 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 +345,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 +363,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 +394,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): @@ -394,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): """ @@ -409,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): """ @@ -422,9 +460,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 +471,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 +518,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] - - #chainLabels=range(len(self.keys())) - #print "Starting at %s in graph: %s"%(self.keys().index(startingVertex),chainLabels) - - cycleList = self.__exploreCyclesRecursively(chain, cycleList) - return cycleList + return self.__exploreCyclesRecursively([startingVertex], []) - 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 +555,22 @@ 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, vertices=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) + vertices = graph.vertices[:] # 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 +586,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 +600,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,17 +620,21 @@ 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) + # 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 cycleList[i]] + return cycleList def isMappingValid(self, other, mapping): @@ -632,17 +643,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,112 +669,101 @@ 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 ################################################################################ -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. - # Initialize callDepth with the size of the largest graph - # Each recursive call to __VF2_match will decrease it by one; + * 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 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! - 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: @@ -767,48 +771,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 @@ -826,172 +827,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 diff --git a/rmgpy/group.pxd b/rmgpy/molecule/group.pxd similarity index 83% rename from rmgpy/group.pxd rename to rmgpy/molecule/group.pxd index fe102a7398..fe38fe0626 100644 --- a/rmgpy/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -24,7 +24,7 @@ # ################################################################################ -from graph cimport Vertex, Edge, Graph +from .graph cimport Vertex, Edge, Graph ################################################################################ @@ -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,16 +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 fromAdjacencyList(str adjlist, bint group=?, bint addH=?) - -cpdef toAdjacencyList(Graph molecule, str label=?, bint group=?, bint removeH=?) + cpdef list findSubgraphIsomorphisms(self, Graph other, dict initialMap=?) diff --git a/rmgpy/group.py b/rmgpy/molecule/group.py similarity index 68% rename from rmgpy/group.py rename to rmgpy/molecule/group.py index a4b4111212..27f258bcc9 100644 --- a/rmgpy/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 ################################################################################ @@ -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=''): @@ -89,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): """ @@ -309,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): """ @@ -435,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): """ @@ -500,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): """ @@ -524,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): @@ -534,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): @@ -545,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 @@ -596,7 +610,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 @@ -604,7 +619,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): """ @@ -670,276 +686,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, 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, 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.pxd b/rmgpy/molecule/molecule.pxd similarity index 66% rename from rmgpy/molecule.pxd rename to rmgpy/molecule/molecule.pxd index 0b13537de7..9a8cc9a2e8 100644 --- a/rmgpy/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 ################################################################################ @@ -36,16 +36,15 @@ 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 - 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) @@ -57,21 +56,23 @@ cdef class Atom(Vertex): ################################################################################ +cpdef object SMILEwriter + 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 ################################################################################ @@ -82,7 +83,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) @@ -94,7 +95,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) @@ -106,39 +107,35 @@ cdef class Molecule(Graph): cpdef deleteHydrogens(self) - cpdef makeHydrogensImplicit(self) - - cpdef makeHydrogensExplicit(self) - 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 other0, dict initialMap=?) + cpdef bint isIsomorphic(self, Graph other, dict initialMap=?) except -2 - cpdef tuple findIsomorphism(self, Graph other0, 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) - 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) @@ -158,20 +155,12 @@ 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) 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.py b/rmgpy/molecule/molecule.py similarity index 60% rename from rmgpy/molecule.py rename to rmgpy/molecule/molecule.py index b420300424..9f10233524 100644 --- a/rmgpy/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -39,9 +39,10 @@ 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 +import openbabel +from .graph import Vertex, Edge, Graph +from .group import GroupAtom, GroupBond, Group, ActionError +from .atomtype import AtomType, atomTypes, getAtomType ################################################################################ @@ -56,7 +57,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 +66,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 +74,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 @@ -83,35 +82,37 @@ 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}, implicitHydrogens={3}, charge={4}, label="{5}")'.format(self.element, self.radicalElectrons, self.spinMultiplicity, self.implicitHydrogens, 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, '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): """ A helper function used when unpickling an object. """ + self.edges = d['edges'] self.connectivity1 = d['connectivity1'] self.connectivity2 = d['connectivity2'] self.connectivity3 = d['connectivity3'] @@ -127,6 +128,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 @@ -141,7 +145,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 +195,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 @@ -292,27 +295,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): """ @@ -342,7 +353,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): """ @@ -429,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 @@ -439,7 +453,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 +460,11 @@ class Molecule(Graph): `InChI` string representing the molecular structure. """ - def __init__(self, atoms=None, bonds=None, implicitH=False, symmetry=1, SMILES='', InChI=''): - Graph.__init__(self, atoms, bonds) - self.implicitHydrogens = implicitH + def __init__(self, atoms=None, symmetry=1, SMILES='', InChI=''): + Graph.__init__(self, atoms) 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,28 +482,24 @@ 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.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): """ @@ -527,13 +535,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): """ @@ -558,7 +566,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): @@ -570,7 +577,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): @@ -579,7 +586,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): @@ -590,16 +597,14 @@ 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 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) @@ -613,92 +618,12 @@ 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: 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 @@ -706,7 +631,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): """ @@ -748,7 +673,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 +681,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 +699,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 +719,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 +738,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): @@ -858,12 +749,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): """ @@ -889,7 +780,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 +788,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 +801,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 +814,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. @@ -937,7 +828,6 @@ def fromOBMol(self, obmol, implicitH=False): cython.declare(atom=Atom, atom1=Atom, atom2=Atom, bond=Bond) self.vertices = [] - self.edges = {} # Add hydrogen atoms to complete molecule if needed obmol.AddHydrogens() @@ -969,9 +859,8 @@ 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] = {} # Add bonds by iterating again through atoms for j in range(0, i): @@ -986,19 +875,13 @@ def fromOBMol(self, obmol, implicitH=False): 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() self.updateAtomTypes() - # Make hydrogens implicit to conserve memory - if implicitH: self.makeHydrogensImplicit() - return self def fromAdjacencyList(self, adjlist): @@ -1007,10 +890,10 @@ 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() - self.makeHydrogensImplicit() return self def toCML(self): @@ -1028,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() @@ -1080,38 +962,29 @@ 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(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() atoms = self.vertices - bonds = self.edges obmol = openbabel.OBMol() for atom in atoms: @@ -1119,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: @@ -1129,19 +1002,14 @@ 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() + from .adjlist import toAdjacencyList + result = toAdjacencyList(self.vertices, label=label, group=False, removeH=removeH) return result def isLinear(self): @@ -1150,7 +1018,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: @@ -1164,18 +1032,15 @@ 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(): + 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 - implicitH = self.implicitHydrogens - self.makeHydrogensExplicit() 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: @@ -1187,11 +1052,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): @@ -1201,397 +1064,26 @@ 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]) + atom1.implicitHydrogens > 1 and len(self.edges[atom2]) + atom2.implicitHydrogens > 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 - 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. """ - - implicitH = self.implicitHydrogens - self.makeHydrogensExplicit() isomers = [self] @@ -1612,10 +1104,6 @@ def generateResonanceIsomers(self): isomers.append(newIsomer) # Move to next resonance isomer index += 1 - - if implicitH: - for isomer in isomers: - isomer.makeHydrogensImplicit() return isomers @@ -1669,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/rmgpy/molecule_draw.py b/rmgpy/molecule/molecule_draw.py similarity index 99% rename from rmgpy/molecule_draw.py rename to rmgpy/molecule/molecule_draw.py index e64bc170cf..225732d3b3 100644 --- a/rmgpy/molecule_draw.py +++ b/rmgpy/molecule/molecule_draw.py @@ -92,7 +92,6 @@ import logging from numpy.linalg import LinAlgError -from rmgpy.molecule import * ################################################################################ @@ -544,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) @@ -589,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 @@ -1133,20 +1132,22 @@ 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. # 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, # 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() @@ -1239,6 +1240,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/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..8a263af027 --- /dev/null +++ b/rmgpy/molecule/symmetry.py @@ -0,0 +1,425 @@ +#!/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 atom.edges.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 + molecule0 = molecule + molecule = molecule0.copy(True) + atom = molecule.vertices[molecule0.vertices.index(atom)] + 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 = atom1.edges[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.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 + 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.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 + 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.edges[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 + bondlist = [] + for atom1, atom2 in bonds: + 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.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) + + # 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 + + for atom in molecule.vertices: + if not molecule.isAtomInCycle(atom): + symmetryNumber *= calculateAtomSymmetryNumber(molecule, atom) + + 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) + + #if molecule.isCyclic(): + # symmetryNumber *= calculateCyclicSymmetryNumber(molecule) + + return symmetryNumber 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/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): 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: diff --git a/setup.py b/setup.py index e53486a2c3..64d0e5ba19 100644 --- a/setup.py +++ b/setup.py @@ -57,12 +57,15 @@ 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=['.']), + Extension('rmgpy.molecule.symmetry', ['rmgpy/molecule/symmetry.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=['.']), 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 94% rename from unittest/atomtypeTest.py rename to unittest/molecule/atomtypeTest.py index 9a30cacd1b..b5b6a7046d 100644 --- a/unittest/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/elementTest.py b/unittest/molecule/elementTest.py similarity index 85% rename from unittest/elementTest.py rename to unittest/molecule/elementTest.py index 3cbc4cb9bb..dcc25cf513 100644 --- a/unittest/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/graphTest.py b/unittest/molecule/graphTest.py similarity index 60% rename from unittest/graphTest.py rename to unittest/molecule/graphTest.py index cf2b3ae7d2..46e9326997 100644 --- a/unittest/graphTest.py +++ b/unittest/molecule/graphTest.py @@ -3,7 +3,7 @@ import unittest -from rmgpy.graph import * +from rmgpy.molecule.graph import * ################################################################################ @@ -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)) - ################################################################################ diff --git a/unittest/groupTest.py b/unittest/molecule/groupTest.py similarity index 89% rename from unittest/groupTest.py rename to unittest/molecule/groupTest.py index 8f8bbd263c..da45a4fad2 100644 --- a/unittest/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 ################################################################################ @@ -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)) diff --git a/unittest/moleculeTest.py b/unittest/molecule/moleculeTest.py similarity index 62% rename from unittest/moleculeTest.py rename to unittest/molecule/moleculeTest.py index 3eb263a272..945678fd79 100644 --- a/unittest/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 ################################################################################ @@ -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,24 +243,9 @@ 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) - 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.implicitHydrogens, atom.implicitHydrogens) - self.assertEqual(self.atom.charge, atom.charge) - self.assertEqual(self.atom.label, atom.label) - ################################################################################ class TestBond(unittest.TestCase): @@ -279,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): @@ -287,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: @@ -298,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: @@ -309,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: @@ -320,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: @@ -331,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': @@ -346,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': @@ -362,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) @@ -376,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) @@ -390,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) @@ -403,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) @@ -416,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) @@ -430,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) @@ -444,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)) @@ -459,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: @@ -482,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): @@ -499,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) @@ -567,32 +537,29 @@ 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)) - 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') 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()) @@ -625,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))) @@ -663,8 +629,6 @@ def testSubgraphIsomorphismAgain(self): 4 H 0 {1,S} """) - molecule.makeHydrogensExplicit() - labeled1 = molecule.getLabeledAtoms().values()[0] labeled2 = group.getLabeledAtoms().values()[0] @@ -672,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))) @@ -703,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))) @@ -725,24 +687,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)) @@ -764,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): """ @@ -778,46 +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)) - - 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) + self.assertFalse(molecule.isBondInCycle(bond)) def testFromSMILESH(self): """ @@ -902,14 +812,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 +909,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..df08c9a253 --- /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.atoms: + for atom2 in atom1.bonds: + 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.atoms: + for atom2 in atom1.bonds: + 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.atoms: + for atom2 in atom1.bonds: + 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.atoms: + for atom2 in atom1.bonds: + 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.atoms: + for atom2 in atom1.bonds: + 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 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 *