Refactoring of Molecule module and change to how graphs are stored.

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):

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 *
 

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):

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 *
 

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):

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 ([email protected]) and the
+#	RMG Team ([email protected])
+#
+#	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