Frag debug rebased 3

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

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:

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

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

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
 

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

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,15 +816,20 @@ 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)
             atomLonePairs[atom] = str(atom.lonePairs)
             # 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

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:

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