diff --git a/examples/rmg/minimal/input.py b/examples/rmg/minimal/input.py index 4ac1156a39..db17316407 100644 --- a/examples/rmg/minimal/input.py +++ b/examples/rmg/minimal/input.py @@ -12,9 +12,23 @@ species( label='ethane', reactive=True, - structure=SMILES("CC"), + structure=fragment_adj("""1 C u0 p0 c0 {2,S} {4,S} {5,S} {6,S} +2 C u0 p0 c0 {1,S} {3,S} {7,S} {8,S} +3 R u0 p0 c0 {2,S} +4 H u0 p0 c0 {1,S} +5 H u0 p0 c0 {1,S} +6 H u0 p0 c0 {1,S} +7 H u0 p0 c0 {2,S} +8 H u0 p0 c0 {2,S} +"""), ) +# species( +# label='ethane', +# reactive=True, +# structure=SMILES("CC"), +# ) + # Reaction systems simpleReactor( temperature=(1350,'K'), @@ -35,7 +49,7 @@ model( toleranceKeepInEdge=0.0, - toleranceMoveToCore=0.1, + toleranceMoveToCore=0.01, toleranceInterruptSimulation=0.1, maximumEdgeSpecies=100000, ) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 5848fa5a8a..1b5ea9d810 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -755,6 +755,8 @@ def loadSpeciesDictionary(path): mapping the species identifiers to the loaded species. Resonance isomers for each species are automatically generated. """ + from afm.fragment import Fragment + import re speciesDict = {} inerts = [Species().fromSMILES(inert) for inert in ('[He]', '[Ne]', 'N#N', '[Ar]')] @@ -763,7 +765,18 @@ def loadSpeciesDictionary(path): for line in f: if line.strip() == '' and adjlist.strip() != '': # Finish this adjacency list - species = Species().fromAdjacencyList(adjlist) + if len(re.findall(r'([LR]\d?)', adjlist)) != 0: + frag = Fragment().fromAdjacencyList(adjlist) + species = Species(molecule = [frag]) + for label in adjlist.splitlines(): + if label.strip(): + break + else: + label = '' + if len(label.split()) > 0 and not label.split()[0].isdigit(): + species.label = label.strip() + else: + species = Species().fromAdjacencyList(adjlist) species.generate_resonance_structures() label = species.label for inert in inerts: @@ -781,7 +794,18 @@ def loadSpeciesDictionary(path): adjlist += line else: #reach end of file if adjlist.strip() != '': - species = Species().fromAdjacencyList(adjlist) + if len(re.findall(r'([LR]\d?)', adjlist)) != 0: + frag = Fragment().fromAdjacencyList(adjlist) + species = Species(molecule = [frag]) + for label in adjlist.splitlines(): + if label.strip(): + break + else: + label = '' + if len(label.split()) > 0 and not label.split()[0].isdigit(): + species.label = label.strip() + else: + species = Species().fromAdjacencyList(adjlist) species.generate_resonance_structures() label = species.label for inert in inerts: diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 43fa934675..350b8c5bd2 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -275,13 +275,26 @@ def getSpecies(self, path, resonance=True): Load the dictionary containing all of the species in a kinetics library or depository. """ from rmgpy.species import Species + from afm.fragment import Fragment + import re speciesDict = OrderedDict() with open(path, 'r') as f: adjlist = '' for line in f: if line.strip() == '' and adjlist.strip() != '': # Finish this adjacency list - species = Species().fromAdjacencyList(adjlist) + if len(re.findall(r'([LR]\d?)', adjlist)) != 0: + frag = Fragment().fromAdjacencyList(adjlist) + species = Species(molecule = [frag]) + for label in adjlist.splitlines(): + if label.strip(): + break + else: + label = '' + if len(label.split()) > 0 and not label.split()[0].isdigit(): + species.label = label.strip() + else: + species = Species().fromAdjacencyList(adjlist) if resonance: species.generate_resonance_structures() label = species.label diff --git a/rmgpy/data/kinetics/common.py b/rmgpy/data/kinetics/common.py index c4bebaa21e..c0b7f912e1 100644 --- a/rmgpy/data/kinetics/common.py +++ b/rmgpy/data/kinetics/common.py @@ -156,8 +156,9 @@ def ensure_species(input_list, resonance=False, keep_isomorphic=False): The input list of :class:`Species` or :class:`Molecule` objects is modified in place to only have :class:`Species` objects. Returns None. """ + from afm.fragment import Fragment for index, item in enumerate(input_list): - if isinstance(item, Molecule): + if isinstance(item, Molecule) or isinstance(item, Fragment): new_item = Species(molecule=[item]) elif isinstance(item, Species): new_item = item diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 46515e2aa4..c25f058c57 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -63,6 +63,8 @@ from rmgpy.exceptions import InvalidActionError, ReactionPairsError, KineticsError,\ UndeterminableKineticsError, ForbiddenStructureException,\ KekulizationError, ActionError, DatabaseError + +from afm.fragment import Fragment, CuttingLabel import itertools ################################################################################ @@ -1245,10 +1247,13 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True): # Also copy structures so we don't modify the originals # Since the tagging has already occurred, both the reactants and the # products will have tags - if isinstance(reactantStructures[0], Group): + if any(isinstance(reactant, Fragment) for reactant in reactantStructures): + reactantStructure = Fragment() + elif isinstance(reactantStructures[0], Group): reactantStructure = Group() - else: + elif isinstance(reactantStructures[0], Molecule): reactantStructure = Molecule() + for s in reactantStructures: reactantStructure = reactantStructure.merge(s.copy(deep=True)) @@ -1334,7 +1339,7 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True): productStructure = reactantStructure if not productStructure.props['validAromatic']: - if isinstance(productStructure, Molecule): + if isinstance(productStructure, Molecule) or isinstance(productStructure, Fragment): # For molecules, kekulize the product to redistribute bonds appropriately productStructure.kekulize() else: @@ -1443,6 +1448,12 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True): # Split product structure into multiple species if necessary productStructures = productStructure.split() + # check if a Fragment is a Molecule, then change it to Molecule + if isinstance(productStructures[0], Fragment): + for index, product in enumerate(productStructures): + if not any(isinstance(vertex, CuttingLabel) for vertex in product.vertices): + mol = Molecule(atoms=product.vertices) + productStructures[index] = mol # Make sure we've made the expected number of products if len(template.products) != len(productStructures): @@ -1459,7 +1470,10 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True): # Remove vdW bonds for struct in productStructures: - struct.removeVanDerWaalsBonds() + if isinstance(struct, Fragment): + continue + else: + struct.removeVanDerWaalsBonds() # Make sure we don't create a different net charge between reactants and products reactant_net_charge = product_net_charge = 0 @@ -1471,7 +1485,7 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True): # If product structures are Molecule objects, update their atom types # If product structures are Group objects and the reaction is in certain families # (families with charged substances), the charge of structures will be updated - if isinstance(struct, Molecule): + if isinstance(struct, Molecule) or isinstance(struct, Fragment): struct.update() elif isinstance(struct, Group): struct.resetRingMembership() @@ -1511,7 +1525,7 @@ def applyRecipe(self, reactantStructures, forward=True, unique=True): labels = [int(label) for label in labels if label] lowest_labels.append(min(labels)) productStructures = [s for _, s in sorted(zip(lowest_labels, productStructures))] - + # Return the product structures return productStructures diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index f2984a99f1..af3c90c6b2 100644 --- a/rmgpy/data/kinetics/library.py +++ b/rmgpy/data/kinetics/library.py @@ -447,7 +447,8 @@ def load(self, path, local_context=None, global_context=None): raise DatabaseError('RMG does not accept reactions with more than 3 products in its solver. Reaction {0} in kinetics library {1} has {2} reactants.'.format(rxn, self.label, len(rxn.products))) if self.autoGenerated: - self.checkForDuplicates(markDuplicates=True) + # self.checkForDuplicates(markDuplicates=True) + pass else: self.checkForDuplicates() self.convertDuplicatesToMulti() diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 791daadfb0..3a08f3932c 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -56,6 +56,7 @@ def saturate(atoms): ''' newAtoms = [] for atom in atoms: + if not isinstance(atom, Atom): continue try: max_number_of_valence_electrons = PeriodicSystem.valence_electrons[atom.symbol] except KeyError: @@ -639,9 +640,13 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): if group: atom = GroupAtom(atomType, unpairedElectrons, partialCharges, label, lonePairs, props) else: - atom = Atom(atomType[0], unpairedElectrons[0], partialCharges[0], label, lonePairs[0]) - if isotope != -1: - atom.element = getElement(atom.number, isotope) + try: + atom = Atom(atomType[0], unpairedElectrons[0], partialCharges[0], label, lonePairs[0]) + if isotope != -1: + atom.element = getElement(atom.number, isotope) + except KeyError: + from afm.fragment import CuttingLabel + atom = CuttingLabel(name=atomType[0], label=label) # Add the atom to the list atoms.append(atom) @@ -707,11 +712,12 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): if not group: # Molecule consistency check # Electron and valency consistency check for each atom - for atom in atoms: ConsistencyChecker.check_partial_charge(atom) + for atom in atoms: + if isinstance(atom, Atom): + ConsistencyChecker.check_partial_charge(atom) nRad = sum([atom.radicalElectrons for atom in atoms]) - absolute_spin_per_electron = 1/2. - if multiplicity == None: multiplicity = 2* (nRad * absolute_spin_per_electron) + 1 + if multiplicity == None: multiplicity = nRad + 1 ConsistencyChecker.check_multiplicity(nRad, multiplicity) for atom in atoms: ConsistencyChecker.check_hund_rule(atom, multiplicity) @@ -756,7 +762,7 @@ def toAdjacencyList(atoms, multiplicity, label=None, group=False, removeH=False, atomNumbers = {} index = 0 for atom in atoms: - if removeH and atom.element.symbol == 'H' and atom.label == '': continue + if removeH and atom.symbol == 'H' and atom.label == '': continue atomNumbers[atom] = '{0:d}'.format(index + 1) index += 1 @@ -810,7 +816,7 @@ def toAdjacencyList(atoms, multiplicity, label=None, group=False, removeH=False, else: for atom in atomNumbers: # Atom type - atomTypes[atom] = '{0}'.format(atom.element.symbol) + atomTypes[atom] = '{0}'.format(atom.symbol) # Unpaired Electron(s) atomUnpairedElectrons[atom] = '{0}'.format(atom.radicalElectrons) # Lone Electron Pair(s) @@ -818,7 +824,12 @@ def toAdjacencyList(atoms, multiplicity, label=None, group=False, removeH=False, # Partial Charge(s) atomCharge[atom] = '+'+str(atom.charge) if atom.charge > 0 else '' + str(atom.charge) # Isotopes - atomIsotope[atom] = atom.element.isotope + if isinstance(atom, Atom): + atomIsotope[atom] = atom.element.isotope + else: + # cutting labels in + # fragment cases + atomIsotope[atom] = atom.isotope # Determine field widths diff --git a/rmgpy/molecule/converter.py b/rmgpy/molecule/converter.py index e21e418c5c..0df42dee03 100644 --- a/rmgpy/molecule/converter.py +++ b/rmgpy/molecule/converter.py @@ -37,6 +37,7 @@ import cython # Assume that rdkit is installed from rdkit import Chem +from rdkit.Chem import rdchem # Test if openbabel is installed try: import openbabel @@ -66,6 +67,7 @@ def toRDKitMol(mol, removeHs=True, returnMapping=False, sanitize=True): mol.sortAtoms() atoms = mol.vertices rdAtomIndices = {} # dictionary of RDKit atom indices + label_dict = {} # store label of atom for Framgent rdkitmol = Chem.rdchem.EditableMol(Chem.rdchem.Mol()) for index, atom in enumerate(mol.vertices): if atom.element.symbol == 'X': @@ -82,6 +84,13 @@ def toRDKitMol(mol, removeHs=True, returnMapping=False, sanitize=True): pass else: rdAtomIndices[atom] = index + if atom.label: + saved_index = index + label = atom.label + if label in label_dict: + label_dict[label].append(saved_index) + else: + label_dict[label] = [saved_index] rdBonds = Chem.rdchem.BondType orders = {'S': rdBonds.SINGLE, 'D': rdBonds.DOUBLE, 'T': rdBonds.TRIPLE, 'B': rdBonds.AROMATIC, 'Q': rdBonds.QUADRUPLE} @@ -99,6 +108,10 @@ def toRDKitMol(mol, removeHs=True, returnMapping=False, sanitize=True): # Make editable mol into a mol and rectify the molecule rdkitmol = rdkitmol.GetMol() + if label_dict: + for label, ind_list in label_dict.iteritems(): + for ind in ind_list: + Chem.SetSupplementalSmilesLabel(rdkitmol.GetAtomWithIdx(ind), label) if sanitize: Chem.SanitizeMol(rdkitmol) if removeHs: diff --git a/rmgpy/molecule/draw.py b/rmgpy/molecule/draw.py index 3dfea871b6..c6dcff24ce 100644 --- a/rmgpy/molecule/draw.py +++ b/rmgpy/molecule/draw.py @@ -54,7 +54,7 @@ from rmgpy.qm.molecule import Geometry from rdkit.Chem import AllChem -from rmgpy.molecule.molecule import Molecule +from rmgpy.molecule.molecule import Atom, Molecule from numpy.linalg import LinAlgError @@ -163,7 +163,7 @@ def draw(self, molecule, format, target=None): self.implicitHydrogens = {} surfaceSites = [] for atom in self.molecule.atoms: - if atom.isHydrogen() and atom.label == '': atomsToRemove.append(atom) + if isinstance(atom, Atom) and atom.isHydrogen() and atom.label == '': atomsToRemove.append(atom) elif atom.isSurfaceSite(): surfaceSites.append(atom) if len(atomsToRemove) < len(self.molecule.atoms) - len(surfaceSites): for atom in atomsToRemove: @@ -949,7 +949,7 @@ def render(self, cr, offset=None): drawLonePairs = False for atom in atoms: - if atom.isNitrogen(): + if isinstance(atom, Atom) and atom.isNitrogen(): drawLonePairs = True left = 0.0 @@ -1219,7 +1219,8 @@ def __renderAtom(self, symbol, atom, x0, y0, cr, heavyFirst=True, drawLonePairs= boundingRect = [x1, y1, x2, y2] # Set color for text - if atom.element.isotope != -1: cr.set_source_rgba(0.0, 0.5, 0.0, 1.0) + if not isinstance(atom, Atom): cr.set_source_rgba(0.0, 0.5, 0.0, 1.0) + elif atom.element.isotope != -1: cr.set_source_rgba(0.0, 0.5, 0.0, 1.0) elif heavyAtom == 'C': cr.set_source_rgba(0.0, 0.0, 0.0, 1.0) elif heavyAtom == 'N': cr.set_source_rgba(0.0, 0.0, 1.0, 1.0) elif heavyAtom == 'O': cr.set_source_rgba(1.0, 0.0, 0.0, 1.0) @@ -1537,13 +1538,13 @@ def draw(self, reaction, format, path=None): for reactant in reaction.reactants: if isinstance(reactant, Species): molecule = reactant.molecule[0] - elif isinstance(reactant, Molecule): + else: molecule = reactant reactants.append(MoleculeDrawer().draw(molecule, format)) for product in reaction.products: if isinstance(product, Species): molecule = product.molecule[0] - elif isinstance(product, Molecule): + else: molecule = product products.append(MoleculeDrawer().draw(molecule, format)) diff --git a/rmgpy/molecule/filtration.py b/rmgpy/molecule/filtration.py index 4c70ab63dd..27ce12c686 100644 --- a/rmgpy/molecule/filtration.py +++ b/rmgpy/molecule/filtration.py @@ -55,6 +55,12 @@ def filter_structures(mol_list, mark_unreactive=True, allow_expanded_octet=True, We often get too many resonance structures from the combination of all rules, particularly for species containing lone pairs. This function filters them out by minimizing the number of C/N/O/S atoms without a full octet. """ + + from afm.fragment import Fragment + if isinstance(mol_list[0], Fragment): + for mol in mol_list: + mol.update() + if not all([(mol.multiplicity == mol_list[0].multiplicity) for mol in mol_list]): raise ValueError("Cannot filter structures with different multiplicities!") @@ -102,12 +108,13 @@ def get_octet_deviation(mol, allow_expanded_octet=True): if `allow_expanded_octet` is ``True`` (by default), then the function also considers dectet for third row elements (currently sulfur is the only hypervalance third row element in RMG) """ - if not isinstance(mol, Molecule): + from afm.fragment import Fragment, CuttingLabel + if not isinstance(mol, (Molecule, Fragment)): raise TypeError("Octet deviation could only be determined for Molecule objects.") octet_deviation = 0 # This is the overall "score" for the molecule, summed across all non-H atoms for atom in mol.vertices: - if atom.isHydrogen(): + if isinstance(atom, CuttingLabel) or atom.isHydrogen(): continue val_electrons = 2 * (int(atom.getBondOrdersForAtom()) + atom.lonePairs) + atom.radicalElectrons if atom.isCarbon() or atom.isNitrogen() or atom.isOxygen(): diff --git a/rmgpy/molecule/kekulize.pyx b/rmgpy/molecule/kekulize.pyx index 21a257a5bb..ca6d68409d 100644 --- a/rmgpy/molecule/kekulize.pyx +++ b/rmgpy/molecule/kekulize.pyx @@ -52,11 +52,12 @@ and the process continues until the entire molecule can be solved. import logging +from .graph cimport Graph, Vertex, Edge from .molecule cimport Atom, Bond, Molecule from .element import PeriodicSystem from rmgpy.exceptions import KekulizationError, AtomTypeError -cpdef kekulize(Molecule mol): +cpdef kekulize(Graph mol): """ Kekulize an aromatic molecule in place. If the molecule cannot be kekulized, a KekulizationError will be raised. However, the molecule will be left in @@ -67,8 +68,8 @@ cpdef kekulize(Molecule mol): """ cdef list ring, rings, aromatic_rings, resolved_rings cdef set endo_bonds, exo_bonds - cdef Atom atom1, atom2, atom - cdef Bond bond + cdef Vertex atom1, atom2, atom + cdef Edge bond cdef bint aromatic, successful, bridged cdef int itercount, maxiter cdef AromaticRing aromatic_ring diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 70e0d363ab..7c0b14b2dd 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -258,7 +258,7 @@ cdef class Molecule(Graph): cpdef bint atomIDValid(self) - cpdef bint isIdentical(self, Molecule other, bint strict=?) except -2 + cpdef bint isIdentical(self, Graph other, bint strict=?) except -2 cpdef dict enumerate_bonds(self) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index a421ff50d9..d84e7e5650 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1341,7 +1341,7 @@ def isIsomorphic(self, other, initialMap=None, generateInitialMap=False, saveOrd """ # 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(other, Molecule): + if not isinstance(other, Graph): raise TypeError('Got a {0} object for parameter "other", when a Molecule object is required.'.format(other.__class__)) # Do the quick isomorphism comparison using the fingerprint # Two fingerprint strings matching is a necessary (but not @@ -2323,8 +2323,9 @@ def isIdentical(self, other, strict=True): If ``strict=False``, performs the check ignoring electrons and resonance structures. """ cython.declare(atomIndicies=set, otherIndices=set, atomList=list, otherList=list, mapping = dict) + from afm.fragment import Fragment - if not isinstance(other, Molecule): + if not isinstance(other, (Molecule, Fragment)): raise TypeError('Got a {0} object for parameter "other", when a Molecule object is required.'.format(other.__class__)) # Get a set of atom indices for each molecule diff --git a/rmgpy/molecule/pathfinder.pxd b/rmgpy/molecule/pathfinder.pxd index 26ea27c1f7..c33e89cc62 100644 --- a/rmgpy/molecule/pathfinder.pxd +++ b/rmgpy/molecule/pathfinder.pxd @@ -26,13 +26,13 @@ ############################################################################### from .molecule cimport Atom, Bond, Molecule -from .graph cimport Vertex, Edge +from .graph cimport Vertex, Edge, Graph -cpdef list find_butadiene(Atom start, Atom end) +cpdef list find_butadiene(Vertex start, Vertex end) -cpdef list find_butadiene_end_with_charge(Atom start) +cpdef list find_butadiene_end_with_charge(Vertex start) -cpdef list find_allyl_end_with_charge(Atom start) +cpdef list find_allyl_end_with_charge(Vertex start) cpdef list find_shortest_path(Vertex start, Vertex end, list path=*) @@ -42,20 +42,20 @@ cpdef list add_allyls(list path) cpdef list add_inverse_allyls(list path) -cpdef dict compute_atom_distance(list atom_indices, Molecule mol) +cpdef dict compute_atom_distance(list atom_indices, Graph mol) -cpdef list find_allyl_delocalization_paths(Atom atom1) +cpdef list find_allyl_delocalization_paths(Vertex atom1) -cpdef list find_lone_pair_multiple_bond_paths(Atom atom1) +cpdef list find_lone_pair_multiple_bond_paths(Vertex atom1) -cpdef list find_adj_lone_pair_radical_delocalization_paths(Atom atom1) +cpdef list find_adj_lone_pair_radical_delocalization_paths(Vertex atom1) -cpdef list find_adj_lone_pair_multiple_bond_delocalization_paths(Atom atom1) +cpdef list find_adj_lone_pair_multiple_bond_delocalization_paths(Vertex atom1) -cpdef list find_adj_lone_pair_radical_multiple_bond_delocalization_paths(Atom atom1) +cpdef list find_adj_lone_pair_radical_multiple_bond_delocalization_paths(Vertex atom1) -cpdef list find_N5dc_radical_delocalization_paths(Atom atom1) +cpdef list find_N5dc_radical_delocalization_paths(Vertex atom1) -cpdef bint is_atom_able_to_gain_lone_pair(Atom atom) +cpdef bint is_atom_able_to_gain_lone_pair(Vertex atom) -cpdef bint is_atom_able_to_lose_lone_pair(Atom atom) +cpdef bint is_atom_able_to_lose_lone_pair(Vertex atom) diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py index 7d47da6f1f..e74cdfc8e4 100644 --- a/rmgpy/molecule/pathfinder.py +++ b/rmgpy/molecule/pathfinder.py @@ -38,7 +38,7 @@ from Queue import Queue from rmgpy.molecule.molecule import Atom, Bond - +from rmgpy.molecule.graph import Vertex, Edge def find_butadiene(start, end): """ @@ -253,7 +253,7 @@ def find_allyl_delocalization_paths(atom1): """ Find all the delocalization paths allyl to the radical center indicated by `atom1`. """ - cython.declare(paths=list, atom2=Atom, atom3=Atom, bond12=Bond, bond23=Bond) + cython.declare(paths=list, atom2=Vertex, atom3=Vertex, bond12=Edge, bond23=Edge) # No paths if atom1 is not a radical if atom1.radicalElectrons <= 0: @@ -283,7 +283,7 @@ def find_lone_pair_multiple_bond_paths(atom1): - N[N+]([O-])=O <=> N[N+](=O)[O-], these structures are isomorphic but not identical, this transition is important for correct degeneracy calculations """ - cython.declare(paths=list, atom2=Atom, atom3=Atom, bond12=Bond, bond23=Bond) + cython.declare(paths=list, atom2=Vertex, atom3=Vertex, bond12=Edge, bond23=Edge) # No paths if atom1 has no lone pairs, or cannot lose them, or is a carbon atom if atom1.lonePairs <= 0 or not is_atom_able_to_lose_lone_pair(atom1) or atom1.isCarbon(): @@ -332,7 +332,7 @@ def find_adj_lone_pair_radical_delocalization_paths(atom1): (where ':' denotes a lone pair, '.' denotes a radical, '-' not in [] denotes a single bond, '-'/'+' denote charge) The bond between the sites does not have to be single, e.g.: [:O.+]=[::N-] <=> [::O]=[:N.] """ - cython.declare(paths=list, atom2=Atom, bond12=Bond) + cython.declare(paths=list, atom2=Vertex, bond12=Edge) paths = [] if atom1.radicalElectrons >= 1 and\ @@ -367,7 +367,7 @@ def find_adj_lone_pair_multiple_bond_delocalization_paths(atom1): (where ':' denotes a lone pair, '.' denotes a radical, '-' not in [] denotes a single bond, '-'/'+' denote charge) (In direction 1 atom1 a lone pair, in direction 2 atom1 a lone pair) """ - cython.declare(paths=list, atom2=Atom, atom3=Atom, bond12=Bond, bond23=Bond) + cython.declare(paths=list, atom2=Vertex, atom3=Vertex, bond12=Bond, bond23=Edge) paths = [] @@ -409,7 +409,7 @@ def find_adj_lone_pair_radical_multiple_bond_delocalization_paths(atom1): (In direction 1 atom1 a lone pair, gains a radical, and atom2 looses a radical. In direction 2 atom1 a lone pair, looses a radical, and atom2 gains a radical) """ - cython.declare(paths=list, atom2=Atom, atom3=Atom, bond12=Bond, bond23=Bond) + cython.declare(paths=list, atom2=Vertex, atom3=Vertex, bond12=Edge, bond23=Edge) paths = [] @@ -445,7 +445,7 @@ def find_N5dc_radical_delocalization_paths(atom1): In this transition atom1 is the middle N+ (N5dc), atom2 is the radical site, and atom3 is negatively charged A "if atom1.atomType.label == 'N5dc'" check should be done before calling this function """ - cython.declare(paths=list, atom2=Atom, atom3=Atom, bond12=Bond, bond23=Bond) + cython.declare(paths=list, atom2=Vertex, atom3=Vertex, bond12=Edge, bond23=Edge) path = [] diff --git a/rmgpy/molecule/resonance.pxd b/rmgpy/molecule/resonance.pxd index a1670aa925..34d213f477 100644 --- a/rmgpy/molecule/resonance.pxd +++ b/rmgpy/molecule/resonance.pxd @@ -30,36 +30,36 @@ from .molecule cimport Atom, Bond, Molecule cpdef list populate_resonance_algorithms(dict features=?) -cpdef dict analyze_molecule(Molecule mol) +cpdef dict analyze_molecule(Graph mol) -cpdef list generate_resonance_structures(Molecule mol, bint clar_structures=?, bint keep_isomorphic=?, bint filter_structures=?) +cpdef list generate_resonance_structures(Graph mol, bint clar_structures=?, bint keep_isomorphic=?, bint filter_structures=?) cpdef list _generate_resonance_structures(list mol_list, list method_list, bint keep_isomorphic=?, bint copy=?, bint filter_structures=?) -cpdef list generate_allyl_delocalization_resonance_structures(Molecule mol) +cpdef list generate_allyl_delocalization_resonance_structures(Graph mol) -cpdef list generate_lone_pair_multiple_bond_resonance_structures(Molecule mol) +cpdef list generate_lone_pair_multiple_bond_resonance_structures(Graph mol) -cpdef list generate_adj_lone_pair_radical_resonance_structures(Molecule mol) +cpdef list generate_adj_lone_pair_radical_resonance_structures(Graph mol) -cpdef list generate_adj_lone_pair_multiple_bond_resonance_structures(Molecule mol) +cpdef list generate_adj_lone_pair_multiple_bond_resonance_structures(Graph mol) -cpdef list generate_adj_lone_pair_radical_multiple_bond_resonance_structures(Molecule mol) +cpdef list generate_adj_lone_pair_radical_multiple_bond_resonance_structures(Graph mol) -cpdef list generate_N5dc_radical_resonance_structures(Molecule mol) +cpdef list generate_N5dc_radical_resonance_structures(Graph mol) -cpdef list generate_isomorphic_resonance_structures(Molecule mol, bint saturate_h=?) +cpdef list generate_isomorphic_resonance_structures(Graph mol, bint saturate_h=?) -cpdef list generate_optimal_aromatic_resonance_structures(Molecule mol, dict features=?) +cpdef list generate_optimal_aromatic_resonance_structures(Graph mol, dict features=?) -cpdef list generate_aromatic_resonance_structure(Molecule mol, list aromatic_bonds=?, bint copy=?) +cpdef list generate_aromatic_resonance_structure(Graph mol, list aromatic_bonds=?, bint copy=?) -cpdef list generate_aryne_resonance_structures(Molecule mol) +cpdef list generate_aryne_resonance_structures(Graph mol) -cpdef list generate_kekule_structure(Molecule mol) +cpdef list generate_kekule_structure(Graph mol) -cpdef list generate_clar_structures(Molecule mol) +cpdef list generate_clar_structures(Graph mol) -cpdef list _clar_optimization(Molecule mol, list constraints=?, max_num=?) +cpdef list _clar_optimization(Graph mol, list constraints=?, max_num=?) -cpdef list _clar_transformation(Molecule mol, list aromatic_ring) +cpdef list _clar_transformation(Graph mol, list aromatic_ring) diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py index dc342357db..f51ce588a0 100644 --- a/rmgpy/molecule/resonance.py +++ b/rmgpy/molecule/resonance.py @@ -254,7 +254,7 @@ def _generate_resonance_structures(mol_list, method_list, keep_isomorphic=False, copy if False, append new resonance structures to input list (default) if True, make a new list with all of the resonance structures """ - cython.declare(index=cython.int, molecule=Molecule, new_mol_list=list, new_mol=Molecule, mol=Molecule) + cython.declare(index=cython.int, molecule=Graph, new_mol_list=list, new_mol=Graph, mol=Graph) if copy: # Make a copy of the list so we don't modify the input list @@ -315,8 +315,8 @@ def generate_allyl_delocalization_resonance_structures(mol): Biradicals on a single atom are not supported. """ - cython.declare(structures=list, paths=list, index=cython.int, structure=Molecule) - cython.declare(atom=Atom, atom1=Atom, atom2=Atom, atom3=Atom, bond12=Bond, bond23=Bond) + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, atom3=Vertex, bond12=Edge, bond23=Edge) cython.declare(v1=Vertex, v2=Vertex) structures = [] @@ -351,8 +351,8 @@ def generate_lone_pair_multiple_bond_resonance_structures(mol): Examples: aniline (Nc1ccccc1), azide, [:NH2]C=[::O] <=> [NH2+]=C[:::O-] (where ':' denotes a lone pair, '.' denotes a radical, '-' not in [] denotes a single bond, '-'/'+' denote charge) """ - cython.declare(structures=list, paths=list, index=cython.int, structure=Molecule) - cython.declare(atom=Atom, atom1=Atom, atom2=Atom, atom3=Atom, bond12=Bond, bond23=Bond) + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, atom3=Vertex, bond12=Edge, bond23=Edge) cython.declare(v1=Vertex, v2=Vertex) structures = [] @@ -392,8 +392,8 @@ def generate_adj_lone_pair_radical_resonance_structures(mol): NO2 example: O=[:N]-[::O.] <=> O=[N.+]-[:::O-] (where ':' denotes a lone pair, '.' denotes a radical, '-' not in [] denotes a single bond, '-'/'+' denote charge) """ - cython.declare(structures=list, paths=list, index=cython.int, structure=Molecule) - cython.declare(atom=Atom, atom1=Atom, atom2=Atom) + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex) cython.declare(v1=Vertex, v2=Vertex) structures = [] @@ -434,8 +434,8 @@ def generate_adj_lone_pair_multiple_bond_resonance_structures(mol): Here atom1 refers to the N/S/O atom, atom 2 refers to the any R!H (atom2's lonePairs aren't affected) (In direction 1 atom1 a lone pair, in direction 2 atom1 a lone pair) """ - cython.declare(structures=list, paths=list, index=cython.int, structure=Molecule, direction=cython.int) - cython.declare(atom=Atom, atom1=Atom, atom2=Atom, bond12=Bond) + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph, direction=cython.int) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, bond12=Edge) cython.declare(v1=Vertex, v2=Vertex) structures = [] @@ -482,8 +482,8 @@ def generate_adj_lone_pair_radical_multiple_bond_resonance_structures(mol): (In direction 1 atom1 a lone pair, gains a radical, and atom2 looses a radical. In direction 2 atom1 a lone pair, looses a radical, and atom2 gains a radical) """ - cython.declare(structures=list, paths=list, index=cython.int, structure=Molecule, direction=cython.int) - cython.declare(atom=Atom, atom1=Atom, atom2=Atom, bond12=Bond) + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph, direction=cython.int) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, bond12=Edge) cython.declare(v1=Vertex, v2=Vertex) structures = [] @@ -573,7 +573,7 @@ def generate_optimal_aromatic_resonance_structures(mol, features=None): In certain cases where multiple forms have the same number of aromatic rings, multiple structures will be returned. If there's an error (eg. in RDKit) it just returns an empty list. """ - cython.declare(molecule=Molecule, rings=list, aromaticBonds=list, kekuleList=list, maxNum=cython.int, mol_list=list, + cython.declare(molecule=Graph, rings=list, aromaticBonds=list, kekuleList=list, maxNum=cython.int, mol_list=list, new_mol_list=list, ring=list, bond=Bond, order=float, originalBonds=list, originalOrder=list, i=cython.int, counter=cython.int) @@ -722,7 +722,7 @@ def generate_aryne_resonance_structures(mol): """ cython.declare(rings=list, ring=list, new_mol_list=list, bond_list=list, i=cython.int, j=cython.int, bond_orders=str, new_orders=str, - ind=cython.int, bond=Bond, new_mol=Molecule) + ind=cython.int, bond=Edge, new_mol=Graph) rings = mol.getRelevantCycles() rings = [ring for ring in rings if len(ring) == 6] @@ -786,9 +786,10 @@ def generate_kekule_structure(mol): Returns a single Kekule structure as an element of a list of length 1. If there's an error (eg. in RDKit) then it just returns an empty list. """ - cython.declare(atom=Atom, molecule=Molecule) + cython.declare(atom=Vertex, molecule=Graph) for atom in mol.atoms: + if not isinstance(atom, Atom): continue if atom.atomType.label == 'Cb' or atom.atomType.label == 'Cbf': break else: @@ -864,8 +865,8 @@ def generate_clar_structures(mol): Returns a list of :class:`Molecule` objects corresponding to the Clar structures. """ - cython.declare(output=list, mol_list=list, new_mol=Molecule, aromaticRings=list, bonds=list, solution=list, - y=list, x=list, index=cython.int, bond=Bond, ring=list) + cython.declare(output=list, mol_list=list, new_mol=Graph, aromaticRings=list, bonds=list, solution=list, + y=list, x=list, index=cython.int, bond=Edge, ring=list) if not mol.isCyclic(): return [] @@ -932,7 +933,7 @@ def _clar_optimization(mol, constraints=None, max_num=None): Hansen, P.; Zheng, M. The Clar Number of a Benzenoid Hydrocarbon and Linear Programming. J. Math. Chem. 1994, 15 (1), 93–107. """ - cython.declare(molecule=Molecule, aromaticRings=list, exo=list, l=cython.int, m=cython.int, n=cython.int, + cython.declare(molecule=Graph, aromaticRings=list, exo=list, l=cython.int, m=cython.int, n=cython.int, a=list, objective=list, status=cython.int, solution=list, innerSolutions=list) from lpsolve55 import lpsolve @@ -1075,7 +1076,7 @@ def _clar_transformation(mol, aromatic_ring): This function directly modifies the input molecule and does not return anything. """ - cython.declare(bondList=list, i=cython.int, atom1=Atom, atom2=Atom, bond=Bond) + cython.declare(bondList=list, i=cython.int, atom1=Vertex, atom2=Vertex, bond=Edge) bond_list = [] diff --git a/rmgpy/molecule/util.py b/rmgpy/molecule/util.py index 63fca67323..cfba9b222a 100644 --- a/rmgpy/molecule/util.py +++ b/rmgpy/molecule/util.py @@ -35,6 +35,7 @@ def retrieveElementCount(obj): """Converts an (augmented) inchi or Molecule into a dictionary element -> count""" + from afm.fragment import Fragment element_count = {} if isinstance(obj, str): @@ -60,7 +61,7 @@ def retrieveElementCount(obj): return element_count - elif isinstance(obj, Molecule): + elif isinstance(obj, Molecule) or isinstance(obj, Fragment): return obj.get_element_count() else: diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index c4aeac7291..e6d5effd2c 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -28,6 +28,7 @@ cimport rmgpy.constants as constants from rmgpy.species cimport Species, TransitionState from rmgpy.molecule.molecule cimport Atom, Molecule +from rmgpy.molecule.graph cimport Vertex, Graph from rmgpy.molecule.element cimport Element from rmgpy.kinetics.model cimport KineticsModel from rmgpy.kinetics.arrhenius cimport Arrhenius diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 491f4f7824..e8ce736e59 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -968,8 +968,9 @@ def isBalanced(self): each side of the reaction equation, or ``False`` if not. """ from rmgpy.molecule.element import elementList + from afm.fragment import CuttingLabel, Fragment - cython.declare(reactantElements=dict, productElements=dict, molecule=Molecule, atom=Atom, element=Element) + cython.declare(reactantElements=dict, productElements=dict, molecule=Graph, atom=Vertex, element=Element) reactantElements = {}; productElements = {} for element in elementList: @@ -979,18 +980,34 @@ def isBalanced(self): for reactant in self.reactants: if isinstance(reactant, Species): molecule = reactant.molecule[0] + for atom in molecule.atoms: + if not isinstance(atom, CuttingLabel): + reactantElements[atom.element] += 1 elif isinstance(reactant, Molecule): molecule = reactant - for atom in molecule.atoms: - reactantElements[atom.element] += 1 + for atom in molecule.atoms: + if not isinstance(atom, CuttingLabel): + reactantElements[atom.element] += 1 + elif isinstance(reactant, Fragment): + for atom in reactant.atoms: + if not isinstance(atom, CuttingLabel): + reactantElements[atom.element] += 1 for product in self.products: if isinstance(product, Species): molecule = product.molecule[0] + for atom in molecule.atoms: + if not isinstance(atom, CuttingLabel): + productElements[atom.element] += 1 elif isinstance(product, Molecule): molecule = product - for atom in molecule.atoms: - productElements[atom.element] += 1 + for atom in molecule.atoms: + if not isinstance(atom, CuttingLabel): + productElements[atom.element] += 1 + elif isinstance(product, Fragment): + for atom in product.atoms: + if not isinstance(atom, CuttingLabel): + productElements[atom.element] += 1 for element in elementList: if reactantElements[element] != productElements[element]: diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index cf6fb128f0..1a0eba6d87 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -177,6 +177,14 @@ def InChI(string): def adjacencyList(string): return Molecule().fromAdjacencyList(string) +def fragment_adj(string): + from afm.fragment import Fragment + return Fragment().fromAdjacencyList(string) + +def fragment_SMILES(string): + from afm.fragment import Fragment + return Fragment().from_SMILES_like_string(string) + # Reaction systems def simpleReactor(temperature, pressure, @@ -738,6 +746,8 @@ def readInputFile(path, rmg0): 'catalystProperties': catalystProperties, 'species': species, 'SMARTS': SMARTS, + 'fragment_adj': fragment_adj, + 'fragment_SMILES': fragment_SMILES, 'SMILES': SMILES, 'InChI': InChI, 'adjacencyList': adjacencyList, diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 8187362050..58c636c98c 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -311,7 +311,7 @@ def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True, spec.creationIteration = self.iterationNum spec.generate_resonance_structures() spec.molecularWeight = Quantity(spec.molecule[0].getMolecularWeight()*1000.,"amu") - + if generateThermo: self.generateThermo(spec) @@ -823,9 +823,12 @@ def generateThermo(self, spc, rename=False): if not spc.thermo: submit(spc, self.solventName) + # if rename and spc.thermo and spc.thermo.label != '': # check if thermo libraries have a name for it + # logging.info('Species {0} renamed {1} based on thermo library name'.format(spc.label, spc.thermo.label)) + # spc.label = spc.thermo.label if rename and spc.thermo and spc.thermo.label != '': # check if thermo libraries have a name for it - logging.info('Species {0} renamed {1} based on thermo library name'.format(spc.label, spc.thermo.label)) - spc.label = spc.thermo.label + logging.info('Species {0} NOT renamed {1} but get thermo based on thermo library'.format(spc.label, spc.thermo.label)) + spc.label = spc.SMILES spc.generateEnergyTransferModel() diff --git a/rmgpy/species.pxd b/rmgpy/species.pxd index 5611ee4d72..5fed38243b 100644 --- a/rmgpy/species.pxd +++ b/rmgpy/species.pxd @@ -33,6 +33,7 @@ from rmgpy.thermo.model cimport HeatCapacityModel from rmgpy.statmech.conformer cimport Conformer from rmgpy.kinetics.model cimport TunnelingModel from rmgpy.molecule.molecule cimport Atom, Bond, Molecule +from rmgpy.molecule.graph cimport Graph ################################################################################ diff --git a/rmgpy/species.py b/rmgpy/species.py index 15b9d51675..7d55e78abd 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -50,6 +50,7 @@ import rmgpy.quantity as quantity +from rmgpy.molecule.graph import Vertex, Edge, Graph from rmgpy.molecule.molecule import Atom, Bond, Molecule from rmgpy.pdep import SingleExponentialDown from rmgpy.statmech.conformer import Conformer @@ -116,7 +117,13 @@ def __init__(self, index=-1, label='', thermo=None, conformer=None, molecule=Non self.molecule = [Molecule(InChI=InChI)] self._inchi = InChI elif SMILES: - self.molecule = [Molecule(SMILES=SMILES)] + # check it is fragment or molecule + import re + from afm.fragment import Fragment + if re.findall(r'([LR]\d?)', SMILES) != []: # Fragment + self.molecule = [Fragment(SMILES=SMILES)] + else: # Molecule + self.molecule = [Molecule(SMILES=SMILES)] self._smiles = SMILES # Check multiplicity of each molecule is the same @@ -234,7 +241,8 @@ def isIsomorphic(self, other, generateInitialMap=False, strict=True): generateInitialMap (bool, optional): If ``True``, make initial map by matching labeled atoms strict (bool, optional): If ``False``, perform isomorphism ignoring electrons. """ - if isinstance(other, Molecule): + from afm.fragment import Fragment + if isinstance(other, Molecule) or isinstance(other, Fragment): for molecule in self.molecule: if molecule.isIsomorphic(other, generateInitialMap=generateInitialMap, strict=strict): return True @@ -254,10 +262,15 @@ def isIdentical(self, other, strict=True): If ``strict=False``, performs the check ignoring electrons and resonance structures. """ + from afm.fragment import Fragment if isinstance(other, Molecule): for molecule in self.molecule: if molecule.isIdentical(other, strict=strict): return True + elif isinstance(other, Fragment): + for molecule in self.molecule: + if molecule.isIdentical(other, strict=strict): + return True elif isinstance(other, Species): for molecule1 in self.molecule: for molecule2 in other.molecule: @@ -334,12 +347,16 @@ def toCantera(self, useChemkinIdentifier = False): instead. Be sure that species' labels are unique when setting it False. """ import cantera as ct + from afm.fragment import CuttingLabel # Determine the number of each type of element in the molecule elementDict = {} # elementCounts = [0,0,0,0] - for atom in self.molecule[0].atoms: + for vertex in self.molecule[0].vertices: # The atom itself - symbol = atom.element.symbol + if not isinstance(vertex, CuttingLabel): + symbol = vertex.element.symbol + else: # that means this vertex is CuttingLabel + symbol = 'Cl' if symbol not in elementDict: elementDict[symbol] = 1 else: @@ -589,7 +606,7 @@ def has_reactive_molecule(self): """ `True` if the species has at least one reactive molecule, `False` otherwise """ - cython.declare(molecule=Molecule) + cython.declare(molecule=Graph) return any([molecule.reactive for molecule in self.molecule]) def copy(self, deep=False): @@ -687,7 +704,12 @@ def generateTransportData(self): raise #count = sum([1 for atom in self.molecule[0].vertices if atom.isNonHydrogen()]) - self.transportData = transportDB.getTransportProperties(self)[0] + if isinstance(self.molecule[0], Molecule): + self.transportData = transportDB.getTransportProperties(self)[0] + else: + # assume it's a species for Fragment + self.molecule[0].assign_representative_species() + self.transportData = transportDB.getTransportProperties(self.molecule[0].species_repr)[0] def getTransportData(self): diff --git a/rmgpy/thermo/thermoengine.py b/rmgpy/thermo/thermoengine.py index 1176c20437..2e37a01d1f 100644 --- a/rmgpy/thermo/thermoengine.py +++ b/rmgpy/thermo/thermoengine.py @@ -132,8 +132,8 @@ def generateThermoData(spc, thermoClass=NASA, solventName=''): except Exception: logging.debug('Could not obtain the thermo database. Not generating thermo...') return None - - thermo0 = thermodb.getThermoData(spc) + + thermo0 = thermodb.getThermoData(spc) # 1. maybe only submit cyclic core # 2. to help radical prediction, HBI should also @@ -168,8 +168,17 @@ def evaluator(spc, solventName = ''): """ logging.debug("Evaluating spc %s ", spc) - spc.generate_resonance_structures() - thermo = generateThermoData(spc,solventName=solventName) + from rmgpy.molecule import Molecule + + if isinstance(spc.molecule[0], Molecule): + spc.generate_resonance_structures() + thermo = generateThermoData(spc, solventName=solventName) + else: + # assume it's a species for Fragment + spc.molecule[0].assign_representative_species() + spc_repr = spc.molecule[0].species_repr + spc_repr.generate_resonance_structures() + thermo = generateThermoData(spc_repr, solventName=solventName) return thermo