diff --git a/examples/cantherm/networks/acetyl+O2/input.py b/examples/cantherm/networks/acetyl+O2/input.py index 1a06738e86..6dde39e9e3 100644 --- a/examples/cantherm/networks/acetyl+O2/input.py +++ b/examples/cantherm/networks/acetyl+O2/input.py @@ -16,6 +16,7 @@ species( label = 'acetylperoxy', + multiplicity = 2, structure = SMILES('CC(=O)O[O]'), E0 = (-34.6,'kcal/mol'), modes = [ @@ -38,6 +39,7 @@ species( label = 'hydroperoxylvinoxy', + multiplicity = 2, structure = SMILES('[CH2]C(=O)OO'), E0 = (-32.4,'kcal/mol'), modes = [ @@ -61,6 +63,7 @@ species( label = 'acetyl', + multiplicity = 2, structure = SMILES('C[C]=O'), E0 = (0.0,'kcal/mol'), #(-20.5205,"kJ/mol") modes = [ @@ -75,6 +78,7 @@ species( label = 'oxygen', + multiplicity = 3, structure = SMILES('[O][O]'), E0 = (0.0,'kcal/mol'), #(-5.74557,"kJ/mol") modes = [ @@ -88,6 +92,7 @@ species( label = 'ketene', + multiplicity = 1, structure = SMILES('C=C=O'), E0 = (-6.6,'kcal/mol'), #modes = [ @@ -101,6 +106,7 @@ species( label = 'lactone', + multiplicity = 1, structure = SMILES('C1OC1(=O)'), E0 = (-30.8,'kcal/mol'), #modes = [ @@ -116,6 +122,7 @@ species( label = 'hydroxyl', + multiplicity = 2, structure = SMILES('[OH]'), E0 = (0.0,'kcal/mol'), #modes = [ @@ -129,6 +136,7 @@ species( label = 'hydroperoxyl', + multiplicity = 2, structure = SMILES('O[O]'), E0 = (0.0,'kcal/mol'), #modes = [ @@ -142,6 +150,7 @@ species( label = 'nitrogen', + multiplicity = 1, structure = SMILES('N#N'), molecularWeight = (28.04,"g/mol"), collisionModel = TransportData(sigma=(3.70,'angstrom'), epsilon=(94.9,'K')), diff --git a/examples/cantherm/networks/n-butanol/input.py b/examples/cantherm/networks/n-butanol/input.py index ec48996eca..f76c823984 100644 --- a/examples/cantherm/networks/n-butanol/input.py +++ b/examples/cantherm/networks/n-butanol/input.py @@ -15,6 +15,7 @@ species( label = 'n-C4H10', + multiplicity = 1, structure = SMILES('CCCCO'), E0 = (-317.807,'kJ/mol'), modes = [ @@ -35,6 +36,7 @@ species( label = 'C4H8', + multiplicity = 1, structure = SMILES('C=CCC'), E0 = (-17.8832,'kJ/mol'), modes = [ @@ -50,6 +52,7 @@ species( label = 'H2O', + multiplicity = 1, structure = SMILES('O'), E0 = (-269.598,'kJ/mol'), modes = [ @@ -63,6 +66,7 @@ species( label = "bath_gas", + multiplicity = 1, E0 = (0,'kJ/mol'), molecularWeight = (28.04,"g/mol"), collisionModel = TransportData(sigma=(3.41,"angstrom"), epsilon=(124,"K")), diff --git a/examples/generateReactions/input.py b/examples/generateReactions/input.py index f0d2dff35b..c37f253000 100644 --- a/examples/generateReactions/input.py +++ b/examples/generateReactions/input.py @@ -11,18 +11,21 @@ # List of species species( label='ethane', + multiplicity = 1, reactive=True, structure=SMILES("CC"), ) species( label='H', + multiplicity = 2, reactive=True, structure=SMILES("[H]"), ) species( label='butane', + multiplicity = 1, reactive=True, structure=SMILES("CCCC"), ) diff --git a/examples/measure/acetylO2/input.py b/examples/measure/acetylO2/input.py index dc1b25c184..91d0a76424 100644 --- a/examples/measure/acetylO2/input.py +++ b/examples/measure/acetylO2/input.py @@ -16,6 +16,7 @@ species( label='acetylperoxy', + multiplicity = 2, SMILES='CC(=O)O[O]', E0=(-34.6,'kcal/mol'), states=States( @@ -39,6 +40,7 @@ species( label='hydroperoxylvinoxy', + multiplicity = 2, SMILES='[CH2]C(=O)OO', E0=(-32.4,'kcal/mol'), states=States( @@ -63,6 +65,7 @@ species( label='acetyl', + multiplicity = 1, SMILES='C[C]=O', E0=(0.0,'kcal/mol'), states=States( @@ -84,6 +87,7 @@ species( label='oxygen', + multiplicity = 3, SMILES='[O][O]', E0=(0.0,'kcal/mol'), states=States( @@ -102,30 +106,35 @@ species( label='ketene', + multiplicity = 1, SMILES='C=C=O', E0=(-6.6,'kcal/mol'), ) species( label='lactone', + multiplicity = 1, SMILES='C1OC1(=O)', E0=(-30.8,'kcal/mol'), ) species( label='hydroxyl', + multiplicity = 2, SMILES='[OH]', E0=(0.0,'kcal/mol'), ) species( label='hydroperoxyl', + multiplicity = 2, SMILES='O[O]', E0=(0.0,'kcal/mol'), ) species( label='nitrogen', + multiplicity = 1, SMILES='N#N', lennardJones=LennardJones(sigma=(3.70,'angstrom'), epsilon=(94.9,'K')), collisionModel = SingleExponentialDown( diff --git a/examples/rmg/1,3-hexadiene/input.py b/examples/rmg/1,3-hexadiene/input.py index 28615b51a8..397faee05d 100644 --- a/examples/rmg/1,3-hexadiene/input.py +++ b/examples/rmg/1,3-hexadiene/input.py @@ -8,28 +8,37 @@ kineticsEstimator = 'rate rules', ) +# Constraints on generated species +generatedSpeciesConstraints( + maximumRadicalElectrons = 2, +) + # List of species species( label='HXD13', + multiplicity = 1, reactive=True, structure=SMILES("C=CC=CCC"), ) species( label='CH4', + multiplicity = 1, reactive=True, structure=SMILES("C"), ) species( label='H2', + multiplicity = 1, reactive=True, structure=adjacencyList( """ - 1 H 0 {2,S} - 2 H 0 {1,S} + 1 H U0 L0 {2,S} + 2 H U0 L0 {1,S} """), ) species( label='N2', + multiplicity = 1, reactive=False, structure=InChI("InChI=1/N2/c1-2"), ) diff --git a/examples/rmg/TEOS/input.py b/examples/rmg/TEOS/input.py index 44fecc5b47..0f04878f81 100644 --- a/examples/rmg/TEOS/input.py +++ b/examples/rmg/TEOS/input.py @@ -11,16 +11,19 @@ # List of species species( label='TEOS', + multiplicity = 1, reactive=True, structure=InChI("InChI=1/C8H20O4Si/c1-5-9-13(10-6-2,11-7-3)12-8-4/h5-8H2,1-4H3"), ) species( label='EtOH', + multiplicity = 1, reactive=True, structure=InChI("InChI=1/C2H6O/c1-2-3/h3H,2H2,1H3"), ) species( label='Ar', + multiplicity = 1, reactive=False, structure=InChI("InChI=1/Ar"), ) diff --git a/examples/rmg/c3h4/input.py b/examples/rmg/c3h4/input.py index ffc50b4d3b..9a737939c1 100644 --- a/examples/rmg/c3h4/input.py +++ b/examples/rmg/c3h4/input.py @@ -8,19 +8,26 @@ kineticsEstimator = 'rate rules', ) +generatedSpeciesConstraints( + maximumRadicalElectrons = 4, +) + # List of species species( label='CH2', reactive=True, + multiplicity = 3, structure=SMILES("[CH2]"), ) species( label='C2H2', + multiplicity = 1, reactive=True, structure=SMILES("C#C"), ) species( label='N2', + multiplicity = 1, reactive=False, structure=InChI("InChI=1/N2/c1-2"), ) diff --git a/examples/rmg/ch3no2/input.py b/examples/rmg/ch3no2/input.py index 749cd7d7e5..d7c3f1edfd 100644 --- a/examples/rmg/ch3no2/input.py +++ b/examples/rmg/ch3no2/input.py @@ -23,36 +23,39 @@ # List of species species( label='CH3NO2', + multiplicity = 1, reactive=True, structure=adjacencyList( """ - 1 C 0 0 {2,S} {3,S} {4,S} {5,S} - 2 H 0 0 {1,S} - 3 H 0 0 {1,S} - 4 H 0 0 {1,S} - 5 N 0 0 {1,S} {6,D} {7,S} - 6 O 0 2 {5,D} - 7 O 0 3 {5,S} + 1 C U0 L0 {2,S} {3,S} {4,S} {5,S} + 2 H U0 L0 {1,S} + 3 H U0 L0 {1,S} + 4 H U0 L0 {1,S} + 5 N U0 L0 {1,S} {6,D} {7,S} + 6 O U0 L2 {5,D} + 7 O U0 L3 {5,S} """), ) species( label='O2', + multiplicity = 3, reactive=True, structure=adjacencyList( """ - 1 O 1 2 {2,S} - 2 O 1 2 {1,S} + 1 O U1 L2 {2,S} + 2 O U1 L2 {1,S} """), ) species( label='N2', + multiplicity = 1, reactive=True, structure=adjacencyList( """ - 1 N 1 1 {2,T} - 2 N 1 1 {1,T} + 1 N U0 L1 {2,T} + 2 N U0 L1 {1,T} """), ) diff --git a/examples/rmg/diesel/input.py b/examples/rmg/diesel/input.py index 1b3a010db4..f616f9168c 100644 --- a/examples/rmg/diesel/input.py +++ b/examples/rmg/diesel/input.py @@ -11,43 +11,51 @@ # List of species species( label='n-decylbz', + multiplicity = 1, reactive=True, structure=SMILES("CCCCCCCCCCc1ccccc1"), ) species( label='n-C11', + multiplicity = 1, reactive=True, structure=SMILES("CCCCCCCCCCC"), ) species( label='n-C13', + multiplicity = 1, reactive=True, structure=SMILES("CCCCCCCCCCCCC"), ) species( label='n-C16', + multiplicity = 1, reactive=True, structure=SMILES("CCCCCCCCCCCCCCCC"), ) species( label='n-C19', + multiplicity = 1, reactive=True, structure=SMILES("CCCCCCCCCCCCCCCCCCC"), ) species( label='n-C21', + multiplicity = 1, reactive=True, structure=SMILES("CCCCCCCCCCCCCCCCCCCCC"), ) species( label='1M-napthalene', + multiplicity = 1, reactive=True, structure=SMILES("Cc1cccc2ccccc12"), ) species( label='O2', + multiplicity = 3, reactive=True, - structure=SMILES("O=O"), + structure=SMILES("[O][O]"), ) # Reaction systems diff --git a/examples/rmg/e85/input.py b/examples/rmg/e85/input.py index 6ffb04c275..34b29799ee 100644 --- a/examples/rmg/e85/input.py +++ b/examples/rmg/e85/input.py @@ -8,34 +8,45 @@ kineticsEstimator = 'rate rules', ) +# Constraints on generated species +generatedSpeciesConstraints( + maximumRadicalElectrons = 2, +) + # List of species species( label='O2', # oxygen + multiplicity = 3, reactive=True, - structure=SMILES("O=O"), + structure=SMILES("[O][O]"), ) species( label='C8H18i', # isooctane + multiplicity = 1, reactive=True, structure=SMILES("CC(C)CC(C)(C)C"), ) species( label='C2H6On', # ethanol + multiplicity = 1, reactive=True, structure=SMILES("CCO"), ) species( label='C7H8t', # toluene + multiplicity = 1, reactive=True, structure=SMILES("Cc1ccccc1"), ) species( label='C6H12n', # hex-1-ene + multiplicity = 1, reactive=True, structure=SMILES("CCCCC=C"), ) species( label='Ar', # argon + multiplicity = 1, reactive=False, structure=SMILES("[Ar]"), ) diff --git a/examples/rmg/liquid_phase/input.py b/examples/rmg/liquid_phase/input.py index c47247abdf..d9893e9a64 100644 --- a/examples/rmg/liquid_phase/input.py +++ b/examples/rmg/liquid_phase/input.py @@ -8,15 +8,22 @@ kineticsEstimator = 'rate rules', ) +# Constraints on generated species +generatedSpeciesConstraints( + maximumRadicalElectrons = 3, +) + # List of species species( label='octane', + multiplicity = 1, reactive=True, structure=SMILES("C(CCCCC)CC"), ) species( label='oxygen', + multiplicity = 3, reactive=True, structure=SMILES("[O][O]"), ) diff --git a/examples/rmg/methylformate/input.py b/examples/rmg/methylformate/input.py index ccd3989236..77ee25e61e 100644 --- a/examples/rmg/methylformate/input.py +++ b/examples/rmg/methylformate/input.py @@ -11,80 +11,95 @@ # List of species species( label='Mfmt', + multiplicity = 1, reactive=True, structure=SMILES("COC=O"), ) species( label='O2', + multiplicity = 3, reactive=True, structure=SMILES("[O][O]"), ) species( label='C2H', + multiplicity = 2, reactive=True, structure=SMILES("C#[C]"), ) species( label='CH', + multiplicity = 4, reactive=True, structure=adjacencyList( """ - 1 C 3 {2,S} - 2 H 0 {1,S} + 1 C U3 L0 {2,S} + 2 H U0 L0 {1,S} """), ) species( label='H2O', + multiplicity = 1, reactive=True, structure=SMILES("O"), ) species( label='H2', + multiplicity = 1, reactive=True, structure=SMILES("[H][H]"), ) species( label='CO', + multiplicity = 1, reactive=True, - structure=SMILES("[C]=O"), + structure=SMILES("C#O"), ) species( label='CO2', + multiplicity = 1, reactive=True, structure=SMILES("C(=O)=O"), ) species( label='CH4', + multiplicity = 1, reactive=True, structure=SMILES("C"), ) species( label='CH3', + multiplicity = 2, reactive=True, structure=SMILES("[CH3]"), ) species( label='CH3OH', + multiplicity = 1, reactive=True, structure=SMILES("CO"), ) species( label='C2H4', + multiplicity = 1, reactive=True, structure=SMILES("C=C"), ) species( label='C2H2', + multiplicity = 1, reactive=True, structure=SMILES("C#C"), ) species( label='CH2O', + multiplicity = 1, reactive=True, structure=SMILES("C=O"), ) species( label='CH3CHO', + multiplicity = 1, reactive=True, structure=SMILES("CC=O"), ) @@ -93,6 +108,7 @@ # Bath gas species( label='Ar', + multiplicity = 1, reactive=False, structure=InChI("InChI=1S/Ar"), ) @@ -157,4 +173,5 @@ drawMolecules=False, generatePlots=False, saveConcentrationProfiles=False, + saveEdgeSpecies=True, ) diff --git a/examples/rmg/minimal/input.py b/examples/rmg/minimal/input.py index 93dd6d1c26..b4f31d77de 100644 --- a/examples/rmg/minimal/input.py +++ b/examples/rmg/minimal/input.py @@ -11,6 +11,7 @@ # List of species species( label='ethane', + multiplicity = 1, reactive=True, structure=SMILES("CC"), ) diff --git a/examples/rmg/minimal_sensitivity/input.py b/examples/rmg/minimal_sensitivity/input.py index c1c57f0936..f8aefedd3e 100644 --- a/examples/rmg/minimal_sensitivity/input.py +++ b/examples/rmg/minimal_sensitivity/input.py @@ -8,9 +8,15 @@ kineticsEstimator = 'rate rules', ) +# Constraints on generated species +generatedSpeciesConstraints( + maximumRadicalElectrons = 2, +) + # List of species species( label='ethane', + multiplicity = 1, reactive=True, structure=SMILES("CC"), ) diff --git a/examples/thermoEstimator/input.py b/examples/thermoEstimator/input.py index 6f1bb29a48..8b838d8662 100644 --- a/examples/thermoEstimator/input.py +++ b/examples/thermoEstimator/input.py @@ -4,22 +4,27 @@ species( label='DIPK', + multiplicity = 1, structure=SMILES("CC(C)C(=O)C(C)C"), ) species( label='O2', + multiplicity = 3, structure=SMILES("[O][O]"), ) species( label='R_tert', + multiplicity = 2, structure=SMILES("CC(C)C(=O)[C](C)C"), ) species( label='R_pri', + multiplicity = 2, structure=SMILES("CC(C)C(=O)C(C)[CH2]"), ) species( label='Cineole', + multiplicity = 1, structure=SMILES('CC12CCC(CC1)C(C)(C)O2'), ) diff --git a/rmgpy/chemkin.py b/rmgpy/chemkin.py index 019d8920de..ed5e359c43 100644 --- a/rmgpy/chemkin.py +++ b/rmgpy/chemkin.py @@ -1556,7 +1556,7 @@ def saveSpeciesDictionary(path, species): with open(path, 'w') as f: for spec in species: try: - f.write(spec.molecule[0].toAdjacencyList(label=getSpeciesIdentifier(spec), removeH=False)) + f.write(spec.molecule[0].toAdjacencyList(label=getSpeciesIdentifier(spec), removeH=False, printMultiplicity=True)) f.write('\n') except: raise ChemkinError('Ran into error saving dictionary for species {0}. Please check your files.'.format(getSpeciesIdentifier(spec))) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 6092a91eca..2031c6bd61 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -68,6 +68,7 @@ class Entry: =================== ======================================================== `index` A unique nonnegative integer index for the entry `label` A unique string identifier for the entry (or '' if not used) + `multiplicity` The multiplicity of this species, multiplicity = 2*total_spin+1 `item` The item that this entry represents `parent` The parent of the entry in the hierarchy (or ``None`` if not used) `children` A list of the children of the entry in the hierarchy (or ``None`` if not used) @@ -84,6 +85,7 @@ class Entry: def __init__(self, index=-1, label='', + multiplicity = 100, item=None, parent=None, children=None, @@ -96,6 +98,7 @@ def __init__(self, ): self.index = index self.label = label + self.multiplicity = multiplicity self.item = item self.parent = parent self.children = children or [] @@ -1250,7 +1253,7 @@ def saveOld(self, path): """ self.saveOldDictionary(path) - def loadEntry(self, label, molecule=None, group=None, shortDesc='', longDesc=''): + def loadEntry(self, label, multiplicity = 200, molecule=None, group=None, shortDesc='', longDesc=''): """ Load an entry from the forbidden structures database. This method is automatically called during loading of the forbidden structures @@ -1271,6 +1274,7 @@ def loadEntry(self, label, molecule=None, group=None, shortDesc='', longDesc='') item = Group().fromAdjacencyList(group) self.entries[label] = Entry( label = label, + multiplicity = multiplicity, item = item, shortDesc = shortDesc, longDesc = longDesc.strip(), @@ -1285,6 +1289,7 @@ def saveEntry(self, f, entry, name='entry'): f.write('{0}(\n'.format(name)) f.write(' label = "{0}",\n'.format(entry.label)) + f.write(' multiplicity = {0},\n'.format(entry.multiplicity)) if isinstance(entry.item, Molecule): f.write(' molecule = \n') diff --git a/rmgpy/data/kinetics/common.py b/rmgpy/data/kinetics/common.py index f6efba6062..33f14ca671 100644 --- a/rmgpy/data/kinetics/common.py +++ b/rmgpy/data/kinetics/common.py @@ -143,7 +143,7 @@ def sortEfficiencies(efficiencies0): elif isinstance(reactant, Species): f.write(' reactant{0:d} = \n'.format(i+1)) f.write('"""\n') - f.write(reactant.molecule[0].toAdjacencyList(label=reactant.label, removeH=False)) + f.write(reactant.molecule[0].toAdjacencyList(label=reactant.label, removeH=False, printMultiplicity=True)) f.write('""",\n') elif isinstance(reactant, Group): f.write(' group{0:d} = \n'.format(i+1)) @@ -170,6 +170,7 @@ def sortEfficiencies(efficiencies0): if not entry.item.reversible: f.write(' reversible = {0!r},\n'.format(entry.item.reversible)) elif isinstance(entry.item, Group): + f.write(' multiplicity = {0},\n'.format(entry.multiplicity)) f.write(' group = \n') f.write('"""\n') f.write(entry.item.toAdjacencyList()) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 300f783e7c..05d6311d23 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -618,13 +618,13 @@ def loadRecipe(self, actions): assert action[0] in ['CHANGE_BOND','FORM_BOND','BREAK_BOND','GAIN_RADICAL','LOSE_RADICAL','GAIN_PAIR','LOSE_PAIR'] self.forwardRecipe.addAction(action) - def loadForbidden(self, label, group, shortDesc='', longDesc=''): + def loadForbidden(self, label, group, multiplicity = [1,2,3,4,5], shortDesc='', longDesc=''): """ Load information about a forbidden structure. """ if not self.forbidden: self.forbidden = ForbiddenStructures() - self.forbidden.loadEntry(label=label, group=group, shortDesc=shortDesc, longDesc=longDesc) + self.forbidden.loadEntry(label=label, group=group, multiplicity = multiplicity, shortDesc=shortDesc, longDesc=longDesc) def saveEntry(self, f, entry): """ @@ -834,7 +834,7 @@ def addKineticsRulesFromTrainingSet(self, thermoDatabase=None): reverse_entries = [] for entry in entries: try: - template = self.getReactionTemplate(deepcopy(entry.item)) + template = self.getReactionTemplate(entry.item.copy()) except UndeterminableKineticsError: # Some entries might be stored in the reverse direction for # this family; save them so we can try this @@ -880,11 +880,11 @@ def addKineticsRulesFromTrainingSet(self, thermoDatabase=None): # Estimate the thermo for the reactants and products 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] + item.reactants = [Species(multiplicity=m.multiplicity,molecule=[m]) for m in item.reactants] for reactant in item.reactants: reactant.generateResonanceIsomers() reactant.thermo = thermoDatabase.getThermoData(reactant) - item.products = [Species(molecule=[m]) for m in item.products] + item.products = [Species(multiplicity=m.multiplicity,molecule=[m]) for m in item.products] for product in item.products: product.generateResonanceIsomers() product.thermo = thermoDatabase.getThermoData(product) @@ -1134,172 +1134,79 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp # Generate other possible electronic states - electronicStructuresList1 = [] - electronicStructuresList2 = [] - - struct1 = productStructures[0] - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() - electronicStructuresList1.append(struct1a) - atoms1 = struct1.getRadicalAtoms() + productStructuresList = [] + totalSpin = [] # total spin times 2 - for atom1 in atoms1: + # implement Angular Momentum Addition Theorem + if len(reactantStructures) == 1: - radical1 = atom1.radicalElectrons - spin1 = atom1.spinMultiplicity + totalSpin = [(reactantStructures[0].multiplicity-1.0)/2.0] - if atom1.label != '' and radical1 > 1 and radical1 < 4: - - if radical1 == 2 and spin1 == 3: - atom1.setSpinMultiplicity(1) - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() - elif radical1 == 2 and spin1 == 1: - atom1.setSpinMultiplicity(3) - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() - elif radical1 == 3 and spin1 == 4: - atom1.setSpinMultiplicity(2) - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() - elif radical1 == 3 and spin1 == 2: - atom1.setSpinMultiplicity(4) - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() - - for electronicStructures in electronicStructuresList1: - if electronicStructures.isIsomorphic(struct1a): - break - else: - electronicStructuresList1.append(struct1a) + elif len(reactantStructures) == 2: - elif radical1 == 4: + spin1 = (reactantStructures[0].multiplicity-1.0)/2.0 + spin2 = (reactantStructures[1].multiplicity-1.0)/2.0 + + count = 0.0 + + while (spin1+spin2-count) >= abs(spin1-spin2): - if spin1 == 5: - atom1.setSpinMultiplicity(3) - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() + totalSpin.append(spin1+spin2-count) + + count += 1 + + if len(productStructures) == 1: + + maxSpin1 = productStructures[0].getRadicalCount()/2.0 + + count = 0.0 + + while (maxSpin1-count) >= 0.0: - atom1.setSpinMultiplicity(1) - struct1b = struct1.copy(True) - struct1b.updateAtomTypes() - elif spin1 == 3: - atom1.setSpinMultiplicity(5) - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() + if (maxSpin1-count) in totalSpin: - atom1.setSpinMultiplicity(1) - struct1b = struct1.copy(True) - struct1b.updateAtomTypes() - elif spin1 == 1: - atom1.setSpinMultiplicity(5) - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() + struct = productStructures[0].copy(deep=True) - atom1.setSpinMultiplicity(3) - struct1b = struct1.copy(True) - struct1b.updateAtomTypes() + struct.multiplicity = int((maxSpin1-count)*2.0+1.0) - for electronicStructures in electronicStructuresList1: - if electronicStructures.isIsomorphic(struct1a): - break - else: - electronicStructuresList1.append(struct1a) + if not self.isMoleculeForbidden(struct): + productStructuresList.append([struct]) - for electronicStructures in electronicStructuresList1: - if electronicStructures.isIsomorphic(struct1b): - break - else: - electronicStructuresList1.append(struct1b) - - if len(productStructures) == 2: - - struct2 = productStructures[1] - struct2a = struct2.copy(True) - struct2a.updateAtomTypes() - electronicStructuresList2.append(struct2a) - atoms2 = struct2.getRadicalAtoms() - - for atom2 in atoms2: + count += 1.0 + + elif len(productStructures) == 2: - radical2 = atom2.radicalElectrons - spin2 = atom2.spinMultiplicity + maxSpin1 = productStructures[0].getRadicalCount()/2.0 + maxSpin2 = productStructures[1].getRadicalCount()/2.0 - if atom2.label != '' and radical2 > 1 and radical2 < 4: - - if radical2 == 2 and spin2 == 3: - atom2.setSpinMultiplicity(1) - struct2a = struct2.copy(True) - struct2a.updateAtomTypes() - elif radical2 == 2 and spin2 == 1: - atom2.setSpinMultiplicity(3) - struct2a = struct2.copy(True) - struct2a.updateAtomTypes() - elif radical2 == 3 and spin2 == 4: - atom2.setSpinMultiplicity(2) - struct2a = struct2.copy(True) - struct2a.updateAtomTypes() - elif radical2 == 3 and spin2 == 2: - atom2.setSpinMultiplicity(4) - struct2a = struct2.copy(True) - struct2a.updateAtomTypes() - - for electronicStructures in electronicStructuresList2: - if electronicStructures.isIsomorphic(struct2a): - break - else: - electronicStructuresList2.append(struct2a) + count1 = 0.0 - elif radical2 == 4: + while (maxSpin1-count1) >= 0.0: - if spin2 == 5: - atom2.setSpinMultiplicity(3) - struct2a = struct2.copy(True) - struct2a.updateAtomTypes() + count2 = 0.0 - atom2.setSpinMultiplicity(1) - struct2b = struct2.copy(True) - struct2b.updateAtomTypes() - elif spin2 == 3: - atom2.setSpinMultiplicity(5) - struct2a = struct2.copy(True) - struct2a.updateAtomTypes() + while (maxSpin2-count2) >= 0.0: + + count = 0.0 - atom2.setSpinMultiplicity(1) - struct2b = struct2.copy(True) - struct2b.updateAtomTypes() - elif spin2 == 1: - atom2.setSpinMultiplicity(5) - struct2a = struct2.copy(True) - struct2a.updateAtomTypes() + while (maxSpin1-count1+maxSpin2-count2-count) >= abs((maxSpin1-count1)-(maxSpin2-count2)): + + if (maxSpin1-count1+maxSpin2-count2-count) in totalSpin: - atom2.setSpinMultiplicity(3) - struct2b = struct2.copy(True) - struct2b.updateAtomTypes() - - for electronicStructures in electronicStructuresList2: - if electronicStructures.isIsomorphic(struct2a): - break - else: - electronicStructuresList2.append(struct2a) + struct1 = productStructures[0].copy(deep=True) + struct2 = productStructures[1].copy(deep=True) + + struct1.multiplicity = int((maxSpin1-count1)*2.0+1.0) + struct2.multiplicity = int((maxSpin2-count2)*2.0+1.0) + + if not self.isMoleculeForbidden(struct1) and not self.isMoleculeForbidden(struct2): + productStructuresList.append([struct1,struct2]) + + count += 1.0 + + count2 += 1.0 - for electronicStructures in electronicStructuresList2: - if electronicStructures.isIsomorphic(struct2b): - break - else: - electronicStructuresList2.append(struct2b) - - if len(productStructures) == 2: - - for structa in electronicStructuresList1: - for structb in electronicStructuresList2: - if not (self.isMoleculeForbidden(structa) or self.isMoleculeForbidden(structb)): - productStructuresList.append([structa,structb]) - elif len(productStructures) == 1: - - for structa in electronicStructuresList1: - if not (self.isMoleculeForbidden(structa)): - productStructuresList.append([structa]) + count1 += 1.0 return productStructuresList @@ -1407,9 +1314,9 @@ def generateReactions(self, reactants, failsSpeciesConstraints=None): for reaction in reactionList: moleculeDict = {} for molecule in reaction.reactants: - moleculeDict[molecule] = Species(molecule=[molecule]) + moleculeDict[molecule] = Species(multiplicity=molecule.multiplicity,molecule=[molecule]) for molecule in reaction.products: - moleculeDict[molecule] = Species(molecule=[molecule]) + moleculeDict[molecule] = Species(multiplicity=molecule.multiplicity,molecule=[molecule]) reaction.reactants = [moleculeDict[molecule] for molecule in reaction.reactants] reaction.products = [moleculeDict[molecule] for molecule in reaction.products] reaction.pairs = [(moleculeDict[reactant],moleculeDict[product]) for reactant, product in reaction.pairs] diff --git a/rmgpy/data/kinetics/groups.py b/rmgpy/data/kinetics/groups.py index 3662c1f9b9..fcffbb21c5 100644 --- a/rmgpy/data/kinetics/groups.py +++ b/rmgpy/data/kinetics/groups.py @@ -81,14 +81,16 @@ def __init__(self, def __repr__(self): return ''.format(self.label) - def loadEntry(self, index, label, group, kinetics, reference=None, referenceType='', shortDesc='', longDesc=''): + def loadEntry(self, index, label, group, kinetics, multiplicity = [1,2,3,4,5], reference=None, referenceType='', shortDesc='', longDesc=''): if group[0:3].upper() == 'OR{' or group[0:4].upper() == 'AND{' or group[0:7].upper() == 'NOT OR{' or group[0:8].upper() == 'NOT AND{': item = makeLogicNode(group) else: item = Group().fromAdjacencyList(group) + item.multiplicity = multiplicity self.entries[label] = Entry( index = index, label = label, + multiplicity = multiplicity, item = item, data = kinetics, reference = reference, @@ -159,10 +161,14 @@ def getReactionTemplate(self, reaction): # template is a list of the actual matched nodes # forwardTemplate is a list of the top level nodes that should be matched if len(template) != len(forwardTemplate): +# print 'len(template):', len(template) +# print 'len(forwardTemplate):', len(forwardTemplate) msg = 'Unable to find matching template for reaction {0} in reaction family {1}.'.format(str(reaction), str(self)) msg += 'Trying to match {0} but matched {1}'.format(str(forwardTemplate),str(template)) +# print 'reactants' # for reactant in reaction.reactants: # print reactant.toAdjacencyList() + '\n' +# print 'products' # for product in reaction.products: # print product.toAdjacencyList() + '\n' raise UndeterminableKineticsError(reaction, message=msg) diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index d1c0a5c0ca..068fbec2ad 100644 --- a/rmgpy/data/kinetics/library.py +++ b/rmgpy/data/kinetics/library.py @@ -274,13 +274,14 @@ def loadEntry(self, longDesc='', ): - reactants = [Species(label=reactant1.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(reactant1)])] - if reactant2 is not None: reactants.append(Species(label=reactant2.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(reactant2)])) - if reactant3 is not None: reactants.append(Species(label=reactant3.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(reactant3)])) + + reactants = [Species(label=reactant1.strip().splitlines()[0].strip(), multiplicity=int(reactant1.strip().splitlines()[1].split()[1]), molecule=[Molecule().fromAdjacencyList(reactant1)])] + if reactant2 is not None: reactants.append(Species(label=reactant2.strip().splitlines()[0].strip(), multiplicity=int(reactant2.strip().splitlines()[1].split()[1]), molecule=[Molecule().fromAdjacencyList(reactant2)])) + if reactant3 is not None: reactants.append(Species(label=reactant3.strip().splitlines()[0].strip(), multiplicity=int(reactant3.strip().splitlines()[1].split()[1]), molecule=[Molecule().fromAdjacencyList(reactant3)])) - products = [Species(label=product1.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(product1)])] - if product2 is not None: products.append(Species(label=product2.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(product2)])) - if product3 is not None: products.append(Species(label=product3.strip().splitlines()[0].strip(), molecule=[Molecule().fromAdjacencyList(product3)])) + products = [Species(label=product1.strip().splitlines()[0].strip(), multiplicity=int(product1.strip().splitlines()[1].split()[1]), molecule=[Molecule().fromAdjacencyList(product1)])] + if product2 is not None: products.append(Species(label=product2.strip().splitlines()[0].strip(), multiplicity=int(product2.strip().splitlines()[1].split()[1]), molecule=[Molecule().fromAdjacencyList(product2)])) + if product3 is not None: products.append(Species(label=product3.strip().splitlines()[0].strip(), multiplicity=int(product3.strip().splitlines()[1].split()[1]), molecule=[Molecule().fromAdjacencyList(product3)])) comment = "Reaction and kinetics from {0}.".format(self.label) if shortDesc.strip(): diff --git a/rmgpy/data/solvationTest.py b/rmgpy/data/solvationTest.py index ee13cf1503..06341d2b58 100644 --- a/rmgpy/data/solvationTest.py +++ b/rmgpy/data/solvationTest.py @@ -47,8 +47,9 @@ def testSoluteGeneration(self): ] for name, smiles, S, B, E, L, A, V in self.testCases: - species = Species(molecule=[Molecule(SMILES=smiles)]) - soluteData = self.database.getSoluteDataFromGroups(Species(molecule=[species.molecule[0]])) + molecule=Molecule(SMILES=smiles) + species = Species(multiplicity=molecule.multiplicity, molecule=[molecule]) + soluteData = self.database.getSoluteDataFromGroups(species) print name, soluteData print self.assertAlmostEqual(soluteData.S, S) print self.assertAlmostEqual(soluteData.B, B) @@ -96,8 +97,9 @@ def testCorrectionGeneration(self): ] for solventName, soluteName, smiles, H, G, S in self.testCases: - species = Species(molecule=[Molecule(SMILES=smiles)]) - soluteData = self.database.getSoluteDataFromGroups(Species(molecule=[species.molecule[0]])) + molecule=Molecule(SMILES=smiles) + species = Species(multiplicity=molecule.multiplicity, molecule=[molecule]) + soluteData = self.database.getSoluteDataFromGroups(species) solventData = self.database.getSolventData(solventName) solvationCorrection = self.database.getSolvationCorrection(soluteData, solventData) #print("Enthalpy of solvation for {0} in {1} is {2} J/mol".format(soluteName, solventName, solvationCorrection.enthalpy)) diff --git a/rmgpy/data/statmech.py b/rmgpy/data/statmech.py index 770f04139d..3d2e6f6533 100644 --- a/rmgpy/data/statmech.py +++ b/rmgpy/data/statmech.py @@ -131,10 +131,14 @@ def loadEntry(self, shortDesc='', longDesc='', ): + + item = Molecule().fromAdjacencyList(molecule) + item.multiplicity = multiplicity + self.entries[label] = Entry( index = index, label = label, - item = Molecule().fromAdjacencyList(molecule), + item = item, data = statmech, reference = reference, referenceType = referenceType, diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 384ea5dc46..30eab72911 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -54,13 +54,14 @@ def saveEntry(f, entry): """ f.write('entry(\n') - f.write(' index = {0:d},\n'.format(entry.index)) - f.write(' label = "{0}",\n'.format(entry.label)) + f.write(' index = {0:d},\n'.format(entry.index)) + f.write(' label = "{0}",\n'.format(entry.label)) + f.write(' multiplicity = {0},\n'.format(entry.multiplicity)) if isinstance(entry.item, Molecule): f.write(' molecule = \n') f.write('"""\n') - f.write(entry.item.toAdjacencyList(removeH=False)) + f.write(entry.item.toAdjacencyList(removeH=False,printMultiplicity=False)) f.write('""",\n') elif isinstance(entry.item, Group): f.write(' group = \n') @@ -184,10 +185,11 @@ class ThermoDepository(Database): def __init__(self, label='', name='', shortDesc='', longDesc=''): Database.__init__(self, label=label, name=name, shortDesc=shortDesc, longDesc=longDesc) - def loadEntry(self, index, label, molecule, thermo, reference=None, referenceType='', shortDesc='', longDesc=''): + def loadEntry(self, index, label, multiplicity, molecule, thermo, reference=None, referenceType='', shortDesc='', longDesc=''): entry = Entry( index = index, label = label, + multiplicity = multiplicity, item = Molecule().fromAdjacencyList(molecule), data = thermo, reference = reference, @@ -217,6 +219,7 @@ def __init__(self, label='', name='', shortDesc='', longDesc=''): def loadEntry(self, index, label, + multiplicity, molecule, thermo, reference=None, @@ -226,6 +229,7 @@ def loadEntry(self, ): molecule = Molecule().fromAdjacencyList(molecule) + molecule.multiplicity = multiplicity # Internal checks for adding entry to the thermo library if label in self.entries.keys(): @@ -233,11 +237,13 @@ def loadEntry(self, for entry in self.entries.values(): if molecule.isIsomorphic(entry.item): - raise DatabaseError('Adjacency list of {0} matches that of existing molecule {1} in thermo library. Please correct your library.'.format(label, entry.label)) + if multiplicity == entry.multiplicity: + raise DatabaseError('Adjacency list and multiplicity of {0} matches that of existing molecule {1} in thermo library. Please correct your library.'.format(label, entry.label)) self.entries[label] = Entry( index = index, label = label, + multiplicity = multiplicity, item = molecule, data = thermo, reference = reference, @@ -281,6 +287,7 @@ def loadEntry(self, label, group, thermo, + multiplicity = [1,2,3,4,5], reference=None, referenceType='', shortDesc='', @@ -290,9 +297,11 @@ def loadEntry(self, item = makeLogicNode(group) else: item = Group().fromAdjacencyList(group) + item.multiplicity = multiplicity self.entries[label] = Entry( index = index, label = label, + multiplicity = multiplicity, item = item, data = thermo, reference = reference, diff --git a/rmgpy/data/thermoTest.py b/rmgpy/data/thermoTest.py index 374ff2de7c..f9614dd27b 100644 --- a/rmgpy/data/thermoTest.py +++ b/rmgpy/data/thermoTest.py @@ -19,8 +19,8 @@ class TestThermoDatabase(unittest.TestCase): # Only load these once to save time database = ThermoDatabase() database.load(os.path.join(settings['database.directory'], 'thermo')) - oldDatabase = ThermoDatabase() - oldDatabase.loadOld(os.path.join(settings['database.directory'], '../output/RMG_database')) +# oldDatabase = ThermoDatabase() +# oldDatabase.loadOld(os.path.join(settings['database.directory'], '../output/RMG_database')) def setUp(self): @@ -29,7 +29,7 @@ def setUp(self): """ self.database = self.__class__.database - self.oldDatabase = self.__class__.oldDatabase +# self.oldDatabase = self.__class__.oldDatabase self.Tlist = [300, 400, 500, 600, 800, 1000, 1500] @@ -65,14 +65,15 @@ def testNewThermoGeneration(self): for smiles, symm, H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500 in self.testCases: Cplist = [Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500] - species = Species(molecule=[Molecule(SMILES=smiles)]) + molecule=Molecule(SMILES=smiles) + species = Species(multiplicity=molecule.multiplicity, molecule=molecule) species.generateResonanceIsomers() species.molecule[0] - thermoData = self.database.getThermoDataFromGroups(Species(molecule=[species.molecule[0]])) + thermoData = self.database.getThermoDataFromGroups(species) molecule = species.molecule[0] for mol in species.molecule[1:]: - thermoData0 = self.database.getAllThermoData(Species(molecule=[mol]))[0][0] - for data in self.database.getAllThermoData(Species(molecule=[mol]))[1:]: + thermoData0 = self.database.getAllThermoData(Species(multiplicity=mol.multiplicity, molecule=[mol]))[0][0] + for data in self.database.getAllThermoData(Species(multiplicity=mol.multiplicity, molecule=[mol]))[1:]: if data.getEnthalpy(298) < thermoData0.getEnthalpy(298): thermoData0 = data if thermoData0.getEnthalpy(298) < thermoData.getEnthalpy(298): @@ -92,14 +93,15 @@ def testSymmetryNumberGeneration(self): to select the stablest resonance isomer. """ for smiles, symm, H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500 in self.testCases: - species = Species(molecule=[Molecule(SMILES=smiles)]) + molecule=Molecule(SMILES=smiles) + species = Species(multiplicity=molecule.multiplicity, molecule=molecule) species.generateResonanceIsomers() - thermoData = self.database.getThermoDataFromGroups(Species(molecule=[species.molecule[0]])) + thermoData = self.database.getThermoDataFromGroups(Species(multiplicity=species.molecule[0].multiplicity, molecule=[species.molecule[0]])) # pick the molecule with lowest H298 molecule = species.molecule[0] for mol in species.molecule[1:]: - thermoData0 = self.database.getAllThermoData(Species(molecule=[mol]))[0][0] - for data in self.database.getAllThermoData(Species(molecule=[mol]))[1:]: + thermoData0 = self.database.getAllThermoData(Species(multiplicity=mol.multiplicity, molecule=[mol]))[0][0] + for data in self.database.getAllThermoData(Species(multiplicity=mol.multiplicity, molecule=[mol]))[1:]: if data.getEnthalpy(298) < thermoData0.getEnthalpy(298): thermoData0 = data if thermoData0.getEnthalpy(298) < thermoData.getEnthalpy(298): @@ -107,30 +109,30 @@ def testSymmetryNumberGeneration(self): molecule = mol self.assertEqual(molecule.calculateSymmetryNumber(), symm, msg="Symmetry number error for {0}".format(smiles)) - @work_in_progress - def testOldThermoGeneration(self): - """ - Test that the old ThermoDatabase generates relatively accurate thermo data. - """ - for smiles, symm, H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500 in self.testCases: - Cplist = [Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500] - species = Species(molecule=[Molecule(SMILES=smiles)]) - species.generateResonanceIsomers() - thermoData = self.oldDatabase.getThermoData(Species(molecule=[species.molecule[0]])) - molecule = species.molecule[0] - for mol in species.molecule[1:]: - thermoData0 = self.oldDatabase.getAllThermoData(Species(molecule=[mol]))[0][0] - for data in self.oldDatabase.getAllThermoData(Species(molecule=[mol]))[1:]: - if data.getEnthalpy(298) < thermoData0.getEnthalpy(298): - thermoData0 = data - if thermoData0.getEnthalpy(298) < thermoData.getEnthalpy(298): - thermoData = thermoData0 - molecule = mol - - self.assertAlmostEqual(H298, thermoData.getEnthalpy(298) / 4184, places=1, msg="H298 error for {0}".format(smiles)) - self.assertAlmostEqual(S298, thermoData.getEntropy(298) / 4.184, places=1, msg="S298 error for {0}".format(smiles)) - for T, Cp in zip(self.Tlist, Cplist): - self.assertAlmostEqual(Cp, thermoData.getHeatCapacity(T) / 4.184, places=1, msg="Cp{1} error for {0}".format(smiles, T)) +# @work_in_progress +# def testOldThermoGeneration(self): +# """ +# Test that the old ThermoDatabase generates relatively accurate thermo data. +# """ +# for smiles, symm, H298, S298, Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500 in self.testCases: +# Cplist = [Cp300, Cp400, Cp500, Cp600, Cp800, Cp1000, Cp1500] +# species = Species(molecule=[Molecule(SMILES=smiles)]) +# species.generateResonanceIsomers() +# thermoData = self.oldDatabase.getThermoData(Species(molecule=[species.molecule[0]])) +# molecule = species.molecule[0] +# for mol in species.molecule[1:]: +# thermoData0 = self.oldDatabase.getAllThermoData(Species(molecule=[mol]))[0][0] +# for data in self.oldDatabase.getAllThermoData(Species(molecule=[mol]))[1:]: +# if data.getEnthalpy(298) < thermoData0.getEnthalpy(298): +# thermoData0 = data +# if thermoData0.getEnthalpy(298) < thermoData.getEnthalpy(298): +# thermoData = thermoData0 +# molecule = mol +# +# self.assertAlmostEqual(H298, thermoData.getEnthalpy(298) / 4184, places=1, msg="H298 error for {0}".format(smiles)) +# self.assertAlmostEqual(S298, thermoData.getEntropy(298) / 4.184, places=1, msg="S298 error for {0}".format(smiles)) +# for T, Cp in zip(self.Tlist, Cplist): +# self.assertAlmostEqual(Cp, thermoData.getHeatCapacity(T) / 4.184, places=1, msg="Cp{1} error for {0}".format(smiles, T)) class TestThermoDatabaseAromatics(TestThermoDatabase): """ diff --git a/rmgpy/data/transport.py b/rmgpy/data/transport.py index 0b72b2e5fa..332065ab22 100644 --- a/rmgpy/data/transport.py +++ b/rmgpy/data/transport.py @@ -52,8 +52,9 @@ def saveEntry(f, entry): database to the file object `f`. """ f.write('entry(\n') - f.write(' index = {0:d},\n'.format(entry.index)) - f.write(' label = "{0}",\n'.format(entry.label)) + f.write(' index = {0:d},\n'.format(entry.index)) + f.write(' label = "{0}",\n'.format(entry.label)) + f.write(' multiplicity = {0},\n'.format(entry.multiplicity)) if isinstance(entry.item, Molecule): f.write(' molecule = \n') @@ -93,6 +94,7 @@ def __init__(self, label='', name='', shortDesc='', longDesc=''): def loadEntry(self, index, label, + multiplicity, molecule, transport, reference=None, @@ -100,10 +102,15 @@ def loadEntry(self, shortDesc='', longDesc='', ): + + item = Molecule().fromAdjacencyList(molecule) + item.multiplicity = multiplicity + self.entries[label] = Entry( index = index, label = label, - item = Molecule().fromAdjacencyList(molecule), + multiplicity =multiplicity, + item = item, data = transport, reference = reference, referenceType = referenceType, diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 01cc58da1e..740b843cea 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -56,6 +56,7 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): atoms = [] atomdict = {} bonds = {} + multiplicity = -1 try: @@ -69,6 +70,13 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): label = lines.pop(0) if len(lines) == 0: raise InvalidAdjacencyListError('No atoms specified in adjacency list.') + + # Skip the second line if it contains a multiplicity + if lines[0].split()[0] == 'multiplicity': + multiplicity = int(lines[0].split()[1]) + lines.pop(0) + if len(lines) == 0: + raise InvalidAdjacencyListError('No atoms specified in adjacency list.') mistake1 = re.compile('\{[^}]*\s+[^}]*\}') # Iterate over the remaining lines, generating Atom or GroupAtom objects @@ -110,71 +118,125 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): atomType = [atomType] index += 1 - # Next is the electron state - radicalElectrons = []; spinMultiplicity = [] - elecState = data[index].upper() - if elecState[0] == '{': - elecState = elecState[1:-1].split(',') + # Next the number of unpaired electrons + unpairedElectrons = [] + uState = data[index].upper() + if uState[0] == 'U': + if uState[1] == '{': + uState = uState[2:-1].split(',') + else: + uState = [uState[1]] + for u in uState: + if u == '0': + unpairedElectrons.append(0) + elif u == '1': + unpairedElectrons.append(1) + elif u == '2': + unpairedElectrons.append(2) + elif u == '3': + unpairedElectrons.append(3) + elif u == '4': + unpairedElectrons.append(4) + elif u == 'x': + unpairedElectrons.append(1) + unpairedElectrons.append(2) + unpairedElectrons.append(3) + unpairedElectrons.append(4) + else: + raise InvalidAdjacencyListError('Number of unpaired electrons not recognized.') + index += 1 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 == '3D': - radicalElectrons.append(3); spinMultiplicity.append(2) - elif e == '3Q': - radicalElectrons.append(3); spinMultiplicity.append(4) - elif e == '4': - radicalElectrons.append(4); spinMultiplicity.append(5) - elif e == '4S': - radicalElectrons.append(4); spinMultiplicity.append(1) - elif e == '4T': - radicalElectrons.append(4); spinMultiplicity.append(3) - elif e == '4V': - radicalElectrons.append(4); spinMultiplicity.append(5) - elif e == 'X': - radicalElectrons.extend([0,1,2,2]) - spinMultiplicity.extend([1,2,1,3]) - index += 1 + raise InvalidAdjacencyListError('Number of unpaired electrons not defined.') - # Next number defines the number of lone electron pairs (if provided) - lonePairElectrons = -1 - # if not group and len(data) > index: + # Next the number of lone electron pairs (if provided) + lonePairs = [] if len(data) > index: - lpState = data[index] - if lpState[0] != '{': - if lpState == '0': - lonePairElectrons = 0 - if lpState == '1': - lonePairElectrons = 1 - if lpState == '2': - lonePairElectrons = 2 - if lpState == '3': - lonePairElectrons = 3 - if lpState == '4': - lonePairElectrons = 4 + lpState = data[index].upper() + if lpState[0] == 'L': + if lpState[1] == '{': + lpState = lpState[2:-1].split(',') + else: + lpState = [lpState[1]] + for l in lpState: + if l == '0': + lonePairs.append(0) + elif l == '1': + lonePairs.append(1) + elif l == '2': + lonePairs.append(2) + elif l == '3': + lonePairs.append(3) + elif l == '4': + lonePairs.append(4) + elif l == 'x': + lonePairs.append(1) + lonePairs.append(2) + lonePairs.append(3) + lonePairs.append(4) + else: + raise InvalidAdjacencyListError('Number of lone electron pairs not recognized.') index += 1 else: - lonePairElectrons = -1 + if group: + lonePairs.append(None) + else: + lonePairs.append(0) else: - lonePairElectrons = -1 + if group: + lonePairs.append(None) + else: + lonePairs.append(0) + + # Next the number of partial charges (if provided) + partialCharges = [] + if len(data) > index: + eState = data[index].upper() + if eState[0] == 'E': + if eState[1] == '{': + eState = eState[2:-1].split(',') + else: + eState = [eState[1:]] + for e in eState: + if e == '0': + partialCharges.append(0) + elif e == '+1': + partialCharges.append(1) + elif e == '+2': + partialCharges.append(2) + elif e == '+3': + partialCharges.append(3) + elif e == '+4': + partialCharges.append(4) + elif e == '-1': + partialCharges.append(-1) + elif e == '-2': + partialCharges.append(-2) + elif e == '-3': + partialCharges.append(-3) + elif e == '-4': + partialCharges.append(-4) + elif e == 'x': + partialCharges.append(1) + partialCharges.append(2) + partialCharges.append(3) + partialCharges.append(4) + partialCharges.append(-1) + partialCharges.append(-2) + partialCharges.append(-3) + partialCharges.append(-4) + else: + raise InvalidAdjacencyListError('Number of partial charges not recognized.') + index += 1 + else: + partialCharges.append(0) + else: + partialCharges.append(0) # Create a new atom based on the above information if group: - atom = GroupAtom(atomType, radicalElectrons, spinMultiplicity, [0 for e in radicalElectrons], label, [lonePairElectrons]) + atom = GroupAtom(atomType, unpairedElectrons, [0 for e in unpairedElectrons], label, lonePairs) else: - atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, label, lonePairElectrons) + atom = Atom(atomType[0], unpairedElectrons[0], partialCharges[0], label, lonePairs[0]) # Add the atom to the list atoms.append(atom) @@ -233,21 +295,25 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): if saturateH: # 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, 'Cl': 1, 'He': 0, 'Ne': 0, 'Ar': 0} + orbitals = {'H': 1, 'C': 4, 'O': 4, 'N': 4, 'S': 4, 'Si': 4, 'Cl': 1, '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] + orbital = orbitals[atom.symbol] except KeyError: - raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown valence for atom "{0}".'.format(atom.symbol)) - radical = atom.radicalElectrons + raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown orbital for atom "{0}".'.format(atom.symbol)) + order = 0 for atom2, bond in atom.bonds.items(): order += orders[bond.order] - count = valence - radical - int(order) + count = orbital - atom.radicalElectrons - atom.lonePairs - int(order) + + if count < 0: + raise InvalidAdjacencyListError('Incorrect electron configuration on atom.') + for i in range(count): - a = Atom('H', 0, 1, 0, '') + a = Atom(element='H', radicalElectrons=0, charge=0, label='', lonePairs=0) b = Bond(atom, a, 'S') newAtoms.append(a) atom.bonds[a] = b @@ -255,7 +321,7 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): atoms.extend(newAtoms) # Calculate the number of lone pair electrons requiring molecule with all hydrogen atoms present - if not group and lonePairElectrons == -1: + if not group and lonePairs is not None: valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0, 'Cl': 1} orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} for atom in atoms: @@ -290,40 +356,34 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): print adjlist raise - return atoms + if not group: + nRad = 0 + for atom in atoms: + nRad += atom.radicalElectrons + + if multiplicity == -1: + multiplicity = nRad + 1 + + n = 0 + while (nRad + 1 - n*2) > 0: + + if (nRad + 1 - n*2) == multiplicity: + break + n=n+1 + else: + print adjlist + raise InvalidAdjacencyListError('Multiplicity not in agreement with total number of radicals.') + + return atoms, multiplicity + else: + + 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 and spinMultiplicity == 2: - electronState = '3D' - elif radicalElectrons == 3 and spinMultiplicity == 4: - electronState = '3Q' - elif radicalElectrons == 4 and spinMultiplicity == 1: - electronState = '4S' - elif radicalElectrons == 4 and spinMultiplicity == 3: - electronState = '4T' - elif radicalElectrons == 4 and spinMultiplicity == 5: - electronState = '4V' - 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, removeLonePairs=False): + +def toAdjacencyList(atoms, multiplicity=0, label=None, group=False, removeH=False, removeLonePairs=False, printMultiplicity=False): """ Convert a chemical graph defined by a list of `atoms` into a string adjacency list. @@ -337,6 +397,8 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai pass if label: adjlist += label + '\n' + + if printMultiplicity: adjlist += 'multiplicity ' + str(multiplicity) + '\n' # Determine the numbers to use for each atom atomNumbers = {} @@ -351,27 +413,36 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai atomTypes = {} atomElectronStates = {} atomLonePairs = {} + atomCharge = {} if group: for atom in atomNumbers: # Atom type(s) - if len(atom.atomType) == 1: + 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) + # Unpaired Electron(s) if len(atom.radicalElectrons) == 1: - atomElectronStates[atom] = getElectronState(atom.radicalElectrons[0], atom.spinMultiplicity[0]) + atomElectronStates[atom] = str(atom.radicalElectrons[0]) else: - atomElectronStates[atom] = '{{{0}}}'.format(','.join([getElectronState(radical, spin) for radical, spin in zip(atom.radicalElectrons, atom.spinMultiplicity)])) + atomElectronStates[atom] = '{{{0}}}'.format(','.join([str(radical) for radical in atom.radicalElectrons])) + # Lone Electron Pair(s) + if atom.lonePairs is None: + atomLonePairs[atom] = None + elif len(atom.lonePairs) == 1: + atomLonePairs[atom] = str(atom.lonePairs[0]) + else: + atomLonePairs[atom] = '{{{0}}}'.format(','.join([str(pair) for pair in atom.lonePairs])) 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)) - if not removeLonePairs: - # Lone Pair(s) - atomLonePairs[atom] = atom.lonePairs + # Unpaired Electron(s) + atomElectronStates[atom] = '{0}'.format(atom.radicalElectrons) + # Lone Electron Pair(s) + atomLonePairs[atom] = atom.lonePairs + # Partial Charge(s) + atomCharge[atom] = atom.charge # Determine field widths atomNumberWidth = max([len(s) for s in atomNumbers.values()]) + 1 @@ -379,7 +450,8 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai if atomLabelWidth > 0: atomLabelWidth += 1 atomTypeWidth = max([len(s) for s in atomTypes.values()]) + 1 atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()]) - atomLonePairWidth = 2 + atomLonePairWidth = 1 + atomChargeWidth = 1 # Assemble the adjacency list for atom in atoms: @@ -391,11 +463,19 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai 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) - if group == False and not removeLonePairs: - # Lone Pair(s) - adjlist += '{0:>{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth) + # Unpaired Electron(s) + adjlist += 'U{0:<{1:d}}'.format(atomElectronStates[atom], atomElectronStateWidth) + # Lone Electron Pair(s) + if atomLonePairs[atom] != 'None': + adjlist += ' L{0:>{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth) + if not group: + # Partial Charge(s) + if atomCharge[atom] > 0: + adjlist += ' E+{0:>{1:d}}'.format(atomCharge[atom], atomChargeWidth) + elif atomCharge[atom] < 0: + adjlist += ' E{0:>{1:d}}'.format(atomCharge[atom], atomChargeWidth) + else: + adjlist += ' E{0:>{1:d}} '.format(atomCharge[atom], atomChargeWidth) # Bonds list atoms2 = atom.bonds.keys() diff --git a/rmgpy/molecule/atomtypeTest.py b/rmgpy/molecule/atomtypeTest.py index 97de266933..01e5ff58df 100644 --- a/rmgpy/molecule/atomtypeTest.py +++ b/rmgpy/molecule/atomtypeTest.py @@ -112,75 +112,75 @@ def setUp(self): self.mol1 = Molecule().fromSMILES('COC(=O)CC=C=CC#C') # self.mol2 = Molecule().fromSMILES('c1ccccc1') ## the fromSMILES method currently Kekulizes, so to test Benzene we use fromAdjacencyList - self.mol2 = Molecule().fromAdjacencyList('''1 C 0 0 {2,B} {6,B} {7,S} - 2 C 0 0 {1,B} {3,B} {8,S} - 3 C 0 0 {2,B} {4,B} {9,S} - 4 C 0 0 {3,B} {5,B} {10,S} - 5 C 0 0 {4,B} {6,B} {11,S} - 6 C 0 0 {1,B} {5,B} {12,S} - 7 H 0 0 {1,S} - 8 H 0 0 {2,S} - 9 H 0 0 {3,S} - 10 H 0 0 {4,S} - 11 H 0 0 {5,S} - 12 H 0 0 {6,S}''') + self.mol2 = Molecule().fromAdjacencyList('''1 C U0 L0 {2,B} {6,B} {7,S} + 2 C U0 L0 {1,B} {3,B} {8,S} + 3 C U0 L0 {2,B} {4,B} {9,S} + 4 C U0 L0 {3,B} {5,B} {10,S} + 5 C U0 L0 {4,B} {6,B} {11,S} + 6 C U0 L0 {1,B} {5,B} {12,S} + 7 H U0 L0 {1,S} + 8 H U0 L0 {2,S} + 9 H U0 L0 {3,S} + 10 H U0 L0 {4,S} + 11 H U0 L0 {5,S} + 12 H U0 L0 {6,S}''') self.mol3 = Molecule().fromSMILES('[H]') self.mol4 = Molecule().fromSMILES( 'O=[Si][Si][Si]=[Si]=[Si][Si]#[Si]SS=S') - self.mol5 = Molecule().fromAdjacencyList('''1 H 0 0 {3,S} - 2 H 0 0 {3,S} - 3 N 0 0 {1,S} {2,S} {4,D} - 4 N 0 2 {3,D}''') + self.mol5 = Molecule().fromAdjacencyList('''1 H U0 L0 {3,S} + 2 H U0 L0 {3,S} + 3 N U0 L0 {1,S} {2,S} {4,D} + 4 N U0 L2 {3,D}''') self.mol6 = Molecule().fromSMILES('[Ar]') self.mol7 = Molecule().fromSMILES('[He]') self.mol8 = Molecule().fromSMILES('[Ne]') - self.mol9 = Molecule().fromAdjacencyList('''1 N 0 1 {2,S} {3,S} {4,S} - 2 H 0 0 {1,S} - 3 H 0 0 {1,S} - 4 H 0 0 {1,S}''') + self.mol9 = Molecule().fromAdjacencyList('''1 N U0 L1 {2,S} {3,S} {4,S} + 2 H U0 L0 {1,S} + 3 H U0 L0 {1,S} + 4 H U0 L0 {1,S}''') - self.mol10 = Molecule().fromAdjacencyList('''1 N 1 1 {2,S} {3,S} - 2 H 0 0 {1,S} - 3 H 0 0 {1,S}''') + self.mol10 = Molecule().fromAdjacencyList('''1 N U1 L1 {2,S} {3,S} + 2 H U0 L0 {1,S} + 3 H U0 L0 {1,S}''') - self.mol11 = Molecule().fromAdjacencyList('''1 N 2 1 {2,S} - 2 H 0 0 {1,S}''') + self.mol11 = Molecule().fromAdjacencyList('''1 N U2 L1 {2,S} + 2 H U0 L0 {1,S}''') - self.mol12 = Molecule().fromAdjacencyList('''1 N 0 1 {2,T} - 2 C 1 0 {1,T}''') + self.mol12 = Molecule().fromAdjacencyList('''1 N U0 L1 {2,T} + 2 C U1 L0 {1,T}''') - self.mol13 = Molecule().fromAdjacencyList('''1 N 0 0 {2,S} {3,S} {4,S} {5,S} - 2 H 0 0 {1,S} - 3 H 0 0 {1,S} - 4 H 0 0 {1,S} - 5 O 0 3 {1,S}''') + self.mol13 = Molecule().fromAdjacencyList('''1 N U0 L0 {2,S} {3,S} {4,S} {5,S} + 2 H U0 L0 {1,S} + 3 H U0 L0 {1,S} + 4 H U0 L0 {1,S} + 5 O U0 L3 {1,S}''') - self.mol14 = Molecule().fromAdjacencyList('''1 N 0 2 {2,D} - 2 N 0 0 {1,D} {3,D} - 3 O 0 2 {2,D}''') + self.mol14 = Molecule().fromAdjacencyList('''1 N U0 L2 {2,D} + 2 N U0 L0 {1,D} {3,D} + 3 O U0 L2 {2,D}''') - self.mol15 = Molecule().fromAdjacencyList('''1 N 0 1 {2,T} - 2 N 0 0 {1,T} {3,S} - 3 O 0 3 {2,S}''') + self.mol15 = Molecule().fromAdjacencyList('''1 N U0 L1 {2,T} + 2 N U0 L0 {1,T} {3,S} + 3 O U0 L3 {2,S}''') - self.mol16 = Molecule().fromAdjacencyList('''1 N 0 1 {2,D} {3,S} - 2 O 0 2 {1,D} - 3 O 1 2 {1,S}''') + self.mol16 = Molecule().fromAdjacencyList('''1 N U0 L1 {2,D} {3,S} + 2 O U0 L2 {1,D} + 3 O U1 L2 {1,S}''') - self.mol17 = Molecule().fromAdjacencyList('''1 N 1 1 {2,D} - 2 O 0 2 {1,D}''') + self.mol17 = Molecule().fromAdjacencyList('''1 N U1 L1 {2,D} + 2 O U0 L2 {1,D}''') - self.mol18 = Molecule().fromAdjacencyList('''1 N 0 0 {2,B} {6,B} {7,S} - 2 C 0 0 {1,B} {3,B} {8,S} - 3 C 0 0 {2,B} {4,B} {9,S} - 4 C 0 0 {3,B} {5,B} {10,S} - 5 C 0 0 {4,B} {6,B} {11,S} - 6 N 0 1 {1,B} {5,B} - 7 O 0 3 {1,S} - 8 H 0 0 {2,S} - 9 H 0 0 {3,S} - 10 H 0 0 {4,S} - 11 H 0 0 {5,S}''') + self.mol18 = Molecule().fromAdjacencyList('''1 N U0 L0 {2,B} {6,B} {7,S} + 2 C U0 L0 {1,B} {3,B} {8,S} + 3 C U0 L0 {2,B} {4,B} {9,S} + 4 C U0 L0 {3,B} {5,B} {10,S} + 5 C U0 L0 {4,B} {6,B} {11,S} + 6 N U0 L1 {1,B} {5,B} + 7 O U0 L3 {1,S} + 8 H U0 L0 {2,S} + 9 H U0 L0 {3,S} + 10 H U0 L0 {4,S} + 11 H U0 L0 {5,S}''') def atomType(self, mol, atomID): diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd index efdc8bdde5..87378d5eb3 100644 --- a/rmgpy/molecule/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -33,7 +33,6 @@ cdef class GroupAtom(Vertex): cdef public list atomType cdef public list radicalElectrons - cdef public list spinMultiplicity cdef public list charge cdef public str label cdef public list lonePairs @@ -80,6 +79,8 @@ cdef class GroupBond(Edge): cdef class Group(Graph): + cdef public list multiplicity + # These read-only attribues act as a "fingerprint" for accelerating # subgraph isomorphism checks cdef public short carbonCount diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index e9eea060ae..891ec5e4c6 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -61,7 +61,6 @@ class GroupAtom(Vertex): =================== =================== ==================================== `atomType` ``list`` The allowed atom types (as :class:`AtomType` objects) `radicalElectrons` ``list`` The allowed numbers of radical electrons (as short integers) - `spinMultiplicity` ``list`` The allowed spin multiplicities (as short integers) `charge` ``list`` The allowed formal charges (as short integers) `label` ``str`` A string label that can be used to tag individual atoms `lonePairs` ``list`` The number of lone electron pairs @@ -69,19 +68,18 @@ class GroupAtom(Vertex): Each list represents a logical OR construct, i.e. an atom will match the group if it matches *any* item in the list. However, the - `radicalElectrons`, `spinMultiplicity`, and `charge` attributes are linked + `radicalElectrons`, and `charge` attributes are linked such that an atom must match values from the same index in each of these in order to match. """ - def __init__(self, atomType=None, radicalElectrons=None, spinMultiplicity=None, charge=None, label='', lonePairs=None): + def __init__(self, atomType=None, radicalElectrons=None, charge=None, label='', lonePairs=None): Vertex.__init__(self) self.atomType = atomType or [] for index in range(len(self.atomType)): if isinstance(self.atomType[index], str): self.atomType[index] = atomTypes[self.atomType[index]] self.radicalElectrons = radicalElectrons or [] - self.spinMultiplicity = spinMultiplicity or [] self.charge = charge or [] self.label = label self.lonePairs = lonePairs or [] @@ -100,7 +98,7 @@ def __reduce__(self): 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, self.lonePairs), d) + return (GroupAtom, (atomType, self.radicalElectrons, self.charge, self.label, self.lonePairs), d) def __setstate__(self, d): """ @@ -132,7 +130,7 @@ def copy(self): Return a deep copy of the :class:`GroupAtom` object. Modifying the attributes of the copy will not affect the original. """ - return GroupAtom(self.atomType[:], self.radicalElectrons[:], self.spinMultiplicity[:], self.charge[:], self.label) + return GroupAtom(self.atomType[:], self.radicalElectrons[:], self.charge[:], self.label) def __changeBond(self, order): """ @@ -191,15 +189,12 @@ def __gainRadical(self, radical): where `radical` specifies the number of radical electrons to add. """ radicalElectrons = [] - spinMultiplicity = [] if any([len(atomType.incrementRadical) == 0 for atomType in self.atomType]): raise ActionError('Unable to update GroupAtom due to GAIN_RADICAL action: Unknown atom type produced from set "{0}".'.format(self.atomType)) - for electron, spin in zip(self.radicalElectrons, self.spinMultiplicity): + for electron in self.radicalElectrons: radicalElectrons.append(electron + radical) - spinMultiplicity.append(spin + radical) - # Set the new radical electron counts and spin multiplicities + # Set the new radical electron counts self.radicalElectrons = radicalElectrons - self.spinMultiplicity = spinMultiplicity def __loseRadical(self, radical): """ @@ -207,28 +202,17 @@ def __loseRadical(self, radical): where `radical` specifies the number of radical electrons to remove. """ radicalElectrons = [] - spinMultiplicity = [] pairs = set() if any([len(atomType.decrementRadical) == 0 for atomType in self.atomType]): raise ActionError('Unable to update GroupAtom due to LOSE_RADICAL action: Unknown atom type produced from set "{0}".'.format(self.atomType)) - for electron, spin in zip(self.radicalElectrons, self.spinMultiplicity): + for electron in self.radicalElectrons: electron = electron - radical if electron < 0: - raise ActionError('Unable to update GroupAtom due to LOSE_RADICAL action: Invalid radical electron set "{0}".'.format(self.radicalElectrons)) - spin = spin - radical - if spin <= 0: - spin += 2 - - pair = (electron,spin) - if pair in pairs: - continue # with next electron,spin pair, so we don't get redundant answers. Otherwise.... - pairs.add(pair) + raise ActionError('Unable to update GroupAtom due to LOSE_RADICAL action: Invalid radical electron set "{0}".'.format(self.radicalElectrons)) radicalElectrons.append(electron) - spinMultiplicity.append(spin) - # Set the new radical electron counts and spin multiplicities + # Set the new radical electron counts self.radicalElectrons = radicalElectrons - self.spinMultiplicity = spinMultiplicity def __gainPair(self, pair): """ @@ -310,14 +294,14 @@ def equivalent(self, other): else: return False # Each free radical electron state in self must have an equivalent in other (and vice versa) - for radical1, spin1 in zip(self.radicalElectrons, self.spinMultiplicity): - for radical2, spin2 in zip(other.radicalElectrons, other.spinMultiplicity): - if radical1 == radical2 and spin1 == spin2: break + for radical1 in self.radicalElectrons: + for radical2 in other.radicalElectrons: + if radical1 == radical2: break else: return False - for radical1, spin1 in zip(other.radicalElectrons, other.spinMultiplicity): - for radical2, spin2 in zip(self.radicalElectrons, self.spinMultiplicity): - if radical1 == radical2 and spin1 == spin2: break + for radical1 in other.radicalElectrons: + for radical2 in self.radicalElectrons: + if radical1 == radical2: break else: return False #Each charge in self must have an equivalent in other (and vice versa) @@ -355,9 +339,9 @@ def isSpecificCaseOf(self, other): else: return False # Each free radical electron state in self must have an equivalent in other (and vice versa) - for radical1, spin1 in zip(self.radicalElectrons, self.spinMultiplicity): # all these must match - for radical2, spin2 in zip(other.radicalElectrons, other.spinMultiplicity): # can match any of these - if radical1 == radical2 and spin1 == spin2: break + for radical1 in self.radicalElectrons: + for radical2 in other.radicalElectrons: + if radical1 == radical2: break else: return False #Each charge in self must have an equivalent in other @@ -514,6 +498,7 @@ class Group(Graph): def __init__(self, atoms=None): Graph.__init__(self, atoms) + self.multiplicity = [1,2,3,4,5] self.updateConnectivityValues() self.updateFingerprint() diff --git a/rmgpy/molecule/groupTest.py b/rmgpy/molecule/groupTest.py index ab4411630b..5a09ceb133 100644 --- a/rmgpy/molecule/groupTest.py +++ b/rmgpy/molecule/groupTest.py @@ -18,7 +18,7 @@ def setUp(self): """ A method called before each unit test in this class. """ - self.atom = GroupAtom(atomType=[atomTypes['Cd']], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + self.atom = GroupAtom(atomType=[atomTypes['Cd']], radicalElectrons=[1], charge=[0], label='*1') def testApplyActionBreakBond(self): """ @@ -26,7 +26,7 @@ def testApplyActionBreakBond(self): """ action = ['BREAK_BOND', '*1', 'S', '*2'] for label, atomType in atomTypes.iteritems(): - atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1') atom = atom0.copy() try: atom.applyAction(action) @@ -34,7 +34,6 @@ def testApplyActionBreakBond(self): for a in atomType.breakBond: self.assertTrue(a in atom.atomType) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons) - self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) except ActionError: @@ -46,7 +45,7 @@ def testApplyActionFormBond(self): """ action = ['FORM_BOND', '*1', 'S', '*2'] for label, atomType in atomTypes.iteritems(): - atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1') atom = atom0.copy() try: atom.applyAction(action) @@ -54,7 +53,6 @@ def testApplyActionFormBond(self): for a in atomType.formBond: self.assertTrue(a in atom.atomType) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons) - self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) except ActionError: @@ -66,7 +64,7 @@ def testApplyActionIncrementBond(self): """ action = ['CHANGE_BOND', '*1', 1, '*2'] for label, atomType in atomTypes.iteritems(): - atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1') atom = atom0.copy() try: atom.applyAction(action) @@ -74,7 +72,6 @@ def testApplyActionIncrementBond(self): for a in atomType.incrementBond: self.assertTrue(a in atom.atomType) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons) - self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) except ActionError: @@ -86,7 +83,7 @@ def testApplyActionDecrementBond(self): """ action = ['CHANGE_BOND', '*1', -1, '*2'] for label, atomType in atomTypes.iteritems(): - atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1') atom = atom0.copy() try: atom.applyAction(action) @@ -94,7 +91,6 @@ def testApplyActionDecrementBond(self): for a in atomType.decrementBond: self.assertTrue(a in atom.atomType) self.assertEqual(atom0.radicalElectrons, atom.radicalElectrons) - self.assertEqual(atom0.spinMultiplicity, atom.spinMultiplicity) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) except ActionError: @@ -106,7 +102,7 @@ def testApplyActionGainRadical(self): """ action = ['GAIN_RADICAL', '*1', 1] for label, atomType in atomTypes.iteritems(): - atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1') atom = atom0.copy() try: atom.applyAction(action) @@ -114,7 +110,6 @@ def testApplyActionGainRadical(self): for a in atomType.incrementRadical: self.assertTrue(a in atom.atomType, "GAIN_RADICAL on {0} gave {1} not {2}".format(atomType, atom.atomType, atomType.incrementRadical)) self.assertEqual(atom0.radicalElectrons, [r - 1 for r in atom.radicalElectrons]) - self.assertEqual(atom0.spinMultiplicity, [s - 1 for s in atom.spinMultiplicity]) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) except ActionError: @@ -126,7 +121,7 @@ def testApplyActionLoseRadical(self): """ action = ['LOSE_RADICAL', '*1', 1] for label, atomType in atomTypes.iteritems(): - atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + atom0 = GroupAtom(atomType=[atomType], radicalElectrons=[1], charge=[0], label='*1') atom = atom0.copy() try: atom.applyAction(action) @@ -134,7 +129,6 @@ def testApplyActionLoseRadical(self): for a in atomType.incrementRadical: self.assertTrue(a in atom.atomType, "LOSE_RADICAL on {0} gave {1} not {2}".format(atomType, atom.atomType, atomType.decrementRadical)) self.assertEqual(atom0.radicalElectrons, [r + 1 for r in atom.radicalElectrons]) - self.assertEqual(atom0.spinMultiplicity, [s + 1 for s in atom.spinMultiplicity]) self.assertEqual(atom0.charge, atom.charge) self.assertEqual(atom0.label, atom.label) except ActionError: @@ -146,8 +140,8 @@ def testEquivalent(self): """ for label1, atomType1 in atomTypes.iteritems(): for label2, atomType2 in atomTypes.iteritems(): - atom1 = GroupAtom(atomType=[atomType1], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') - atom2 = GroupAtom(atomType=[atomType2], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + atom1 = GroupAtom(atomType=[atomType1], radicalElectrons=[1], charge=[0], label='*1') + atom2 = GroupAtom(atomType=[atomType2], radicalElectrons=[1], charge=[0], label='*1') if label1 == label2 or atomType2 in atomType1.generic or atomType1 in atomType2.generic: self.assertTrue(atom1.equivalent(atom2), '{0!s} is not equivalent to {1!s}'.format(atom1, atom2)) self.assertTrue(atom2.equivalent(atom1), '{0!s} is not equivalent to {1!s}'.format(atom2, atom1)) @@ -157,7 +151,7 @@ def testEquivalent(self): # Now see if charge and radical count are checked properly for charge in range(3): for radicals in range(2): - atom3 = GroupAtom(atomType=[atomType1], radicalElectrons=[radicals], spinMultiplicity=[radicals + 1], charge=[charge], label='*1') + atom3 = GroupAtom(atomType=[atomType1], radicalElectrons=[radicals], charge=[charge], label='*1') if radicals == 1 and charge == 0: self.assertTrue(atom1.equivalent(atom3), '{0!s} is not equivalent to {1!s}'.format(atom1, atom3)) self.assertTrue(atom1.equivalent(atom3), '{0!s} is not equivalent to {1!s}'.format(atom3, atom1)) @@ -171,11 +165,11 @@ def testIsSpecificCaseOf(self): """ for label1, atomType1 in atomTypes.iteritems(): for label2, atomType2 in atomTypes.iteritems(): - atom1 = GroupAtom(atomType=[atomType1], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') - atom2 = GroupAtom(atomType=[atomType2], radicalElectrons=[1], spinMultiplicity=[2], charge=[0], label='*1') + atom1 = GroupAtom(atomType=[atomType1], radicalElectrons=[1], charge=[0], label='*1') + atom2 = GroupAtom(atomType=[atomType2], radicalElectrons=[1], charge=[0], label='*1') # And make more generic types of these two atoms - atom1gen = GroupAtom(atomType=[atomType1], radicalElectrons=[0, 1], spinMultiplicity=[1, 2], charge=[0, 1], label='*1') - atom2gen = GroupAtom(atomType=[atomType2], radicalElectrons=[0, 1], spinMultiplicity=[1, 2], charge=[0, 1], label='*1') + atom1gen = GroupAtom(atomType=[atomType1], radicalElectrons=[0, 1], charge=[0, 1], label='*1') + atom2gen = GroupAtom(atomType=[atomType2], radicalElectrons=[0, 1], charge=[0, 1], label='*1') if label1 == label2 or atomType2 in atomType1.generic: self.assertTrue(atom1.isSpecificCaseOf(atom2), '{0!s} is not a specific case of {1!s}'.format(atom1, atom2)) self.assertTrue(atom1.isSpecificCaseOf(atom2gen), '{0!s} is not a specific case of {1!s}'.format(atom1, atom2gen)) @@ -193,7 +187,6 @@ def testCopy(self): 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) @@ -207,7 +200,6 @@ def testPickle(self): 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) @@ -362,9 +354,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} U0 {2,{S,D}} {3,S} +2 *1 {Os,Od} U0 {1,{S,D}} +3 R!H U0 {1,S} """ self.group = Group().fromAdjacencyList(self.adjlist) @@ -447,6 +439,7 @@ def testToAdjacencyList(self): Test the Group.toAdjacencyList() method. """ adjlist = self.group.toAdjacencyList() + print 'adjlist', adjlist self.assertEqual(adjlist.strip(), self.adjlist.strip()) def testIsIsomorphic(self): @@ -454,9 +447,9 @@ def testIsIsomorphic(self): Test the Group.isIsomorphic() method. """ adjlist = """ -1 *1 {Os,Od} 0 {3,{S,D}} -2 R!H 0 {3,S} -3 *2 {Cs,Cd} 0 {1,{S,D}} {2,S} +1 *1 {Os,Od} U0 {3,{S,D}} +2 R!H U0 {3,S} +3 *2 {Cs,Cd} U0 {1,{S,D}} {2,S} """ group = Group().fromAdjacencyList(adjlist) self.assertTrue(self.group.isIsomorphic(group)) @@ -467,9 +460,9 @@ def testFindIsomorphism(self): Test the Group.findIsomorphism() method. """ adjlist = """ -1 *1 {Os,Od} 0 {3,{S,D}} -2 R!H 0 {3,S} -3 *2 {Cs,Cd} 0 {1,{S,D}} {2,S} +1 *1 {Os,Od} U0 {3,{S,D}} +2 R!H U0 {3,S} +3 *2 {Cs,Cd} U0 {1,{S,D}} {2,S} """ group = Group().fromAdjacencyList(adjlist) result = self.group.findIsomorphism(group) @@ -491,7 +484,7 @@ def testIsSubgraphIsomorphic(self): Test the Group.isSubgraphIsomorphic() method. """ adjlist = """ -1 *1 {Cs,Cd} 0 +1 *1 {Cs,Cd} U0 """ group = Group().fromAdjacencyList(adjlist) self.assertTrue(self.group.isSubgraphIsomorphic(group)) @@ -502,7 +495,7 @@ def testFindSubgraphIsomorphisms(self): Test the Group.findSubgraphIsomorphisms() method. """ adjlist = """ -1 *1 {Cs,Cd} 0 +1 *1 {Cs,Cd} U0 """ group = Group().fromAdjacencyList(adjlist) result = self.group.findSubgraphIsomorphisms(group) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 3e474908b8..e18e5482cc 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -38,7 +38,6 @@ cdef class Atom(Vertex): cdef public Element element cdef public short radicalElectrons - cdef public short spinMultiplicity cdef public short charge cdef public str label cdef public AtomType atomType @@ -103,6 +102,7 @@ cdef class Molecule(Graph): cdef public bint implicitHydrogens cdef public int symmetryNumber + cdef public int multiplicity cdef public object rdMol cdef public int rdMolConfId cdef str _fingerprint @@ -185,7 +185,7 @@ cdef class Molecule(Graph): # cpdef tRDKitMol(self) - cpdef toAdjacencyList(self, str label=?, bint removeH=?, bint removeLonePairs=?) + cpdef toAdjacencyList(self, str label=?, bint removeH=?, bint removeLonePairs=?, printMultiplicity=?) cpdef bint isLinear(self) except -2 diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index cd1cac7104..d4d409a622 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -97,7 +97,6 @@ class Atom(Vertex): `atomType` :class:`AtomType` The :ref:`atom type ` `element` :class:`Element` The chemical element the atom represents `radicalElectrons` ``short`` The number of radical electrons - `spinMultiplicity` ``short`` The spin multiplicity of the atom `charge` ``short`` The formal charge of the atom `label` ``str`` A string label that can be used to tag individual atoms `coords` ``numpy array`` The (x,y,z) coordinates in Angstrom @@ -109,14 +108,13 @@ class Atom(Vertex): e.g. ``atom.symbol`` instead of ``atom.element.symbol``. """ - def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge=0, label='', lonePairs=0, coords=numpy.array([])): + def __init__(self, element=None, radicalElectrons=0, charge=0, label='', lonePairs=-100, coords=numpy.array([])): Vertex.__init__(self) if isinstance(element, str): self.element = elements.__dict__[element] else: self.element = element self.radicalElectrons = radicalElectrons - self.spinMultiplicity = spinMultiplicity self.charge = charge self.label = label self.atomType = None @@ -150,8 +148,9 @@ def __reduce__(self): 'connectivity3': self.connectivity3, 'sortingLabel': self.sortingLabel, 'atomType': self.atomType.label if self.atomType else None, + 'lonePairs': self.lonePairs, } - return (Atom, (self.element.symbol, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label), d) + return (Atom, (self.element.symbol, self.radicalElectrons, self.charge, self.label), d) def __setstate__(self, d): """ @@ -163,6 +162,7 @@ def __setstate__(self, d): self.connectivity3 = d['connectivity3'] self.sortingLabel = d['sortingLabel'] self.atomType = atomTypes[d['atomType']] if d['atomType'] else None + self.lonePairs = d['lonePairs'] @property def mass(self): return self.element.mass @@ -189,7 +189,6 @@ def equivalent(self, other): atom = other return (self.element is atom.element and self.radicalElectrons == atom.radicalElectrons and - self.spinMultiplicity == atom.spinMultiplicity and self.charge == atom.charge) elif isinstance(other, GroupAtom): cython.declare(a=AtomType, radical=cython.short, spin=cython.short, charge=cython.short) @@ -198,8 +197,8 @@ def equivalent(self, other): if self.atomType.equivalent(a): break else: return False - for radical, spin in zip(ap.radicalElectrons, ap.spinMultiplicity): - if self.radicalElectrons == radical and self.spinMultiplicity == spin: break + for radical in ap.radicalElectrons: + if self.radicalElectrons == radical: break else: return False for charge in ap.charge: @@ -219,7 +218,7 @@ def isSpecificCaseOf(self, other): if isinstance(other, Atom): return self.equivalent(other) elif isinstance(other, GroupAtom): - cython.declare(atom=GroupAtom, a=AtomType, radical=cython.short, spin=cython.short, charge=cython.short, index=cython.int) + cython.declare(atom=GroupAtom, a=AtomType, radical=cython.short, charge=cython.short, index=cython.int) atom = other if self.atomType is None: return False @@ -229,8 +228,7 @@ def isSpecificCaseOf(self, other): return False for index in range(len(atom.radicalElectrons)): radical = atom.radicalElectrons[index] - spin = atom.spinMultiplicity[index] - if self.radicalElectrons == radical and self.spinMultiplicity == spin: break + if self.radicalElectrons == radical: break else: return False # until we have charges and lone pairs in the group values we neglect them here @@ -252,7 +250,6 @@ def copy(self): a.resetConnectivityValues() a.element = self.element a.radicalElectrons = self.radicalElectrons - a.spinMultiplicity = self.spinMultiplicity a.charge = self.charge a.label = self.label a.atomType = self.atomType @@ -300,16 +297,10 @@ def incrementRadical(self): Update the atom pattern as a result of applying a GAIN_RADICAL action, where `radical` specifies the number of radical electrons to add. """ - # Set the new radical electron counts and spin multiplicities + # Set the new radical electron count self.radicalElectrons += 1 if self.radicalElectrons <= 0: raise ActionError('Unable to update Atom due to GAIN_RADICAL action: Invalid radical electron set "{0}".'.format(self.radicalElectrons)) - elif self.radicalElectrons == 1: - self.spinMultiplicity = 2 - elif self.radicalElectrons == 2: - self.spinMultiplicity = 3 # Assume this always results in the triplet, as they tend to be more stable than the singlet (though there are exceptions!) - else: - self.spinMultiplicity = self.radicalElectrons + 1 def decrementRadical(self): """ @@ -317,18 +308,10 @@ def decrementRadical(self): where `radical` specifies the number of radical electrons to remove. """ cython.declare(radicalElectrons=cython.short) - # Set the new radical electron counts and spin multiplicities + # Set the new radical electron count radicalElectrons = self.radicalElectrons = self.radicalElectrons - 1 if radicalElectrons < 0: raise ActionError('Unable to update Atom due to LOSE_RADICAL action: Invalid radical electron set "{0}".'.format(self.radicalElectrons)) - elif radicalElectrons == 0: - self.spinMultiplicity = 1 - elif radicalElectrons == 1: - self.spinMultiplicity = 2 - elif radicalElectrons == 2: - self.spinMultiplicity = 3 # Assume this always results in the triplet, as they tend to be more stable than the singlet (though there are exceptions!) - else: - self.spinMultiplicity = self.radicalElectrons + 1 def setLonePairs(self,lonePairs): """ @@ -590,15 +573,17 @@ class Molecule(Graph): Attribute Type Description ======================= =========== ======================================== `symmetryNumber` ``int`` The (estimated) external + internal symmetry number of the molecule + `multiplicity` ``int`` The multiplicity of this species, multiplicity = 2*total_spin+1 ======================= =========== ======================================== A new molecule object can be easily instantiated by passing the `SMILES` or `InChI` string representing the molecular structure. """ - def __init__(self, atoms=None, symmetry=1, SMILES='', InChI='', SMARTS = ''): + def __init__(self, atoms=None, symmetry=1, multiplicity=187, SMILES='', InChI='', SMARTS = ''): Graph.__init__(self, atoms) self.symmetryNumber = symmetry + self.multiplicity = multiplicity self._fingerprint = None if SMILES != '': self.fromSMILES(SMILES) elif InChI != '': self.fromInChI(InChI) @@ -620,7 +605,7 @@ def __reduce__(self): """ A helper function used when pickling an object. """ - return (Molecule, (self.vertices, self.symmetryNumber)) + return (Molecule, (self.vertices, self.symmetryNumber, self.multiplicity)) def __getAtoms(self): return self.vertices def __setAtoms(self, atoms): self.vertices = atoms @@ -786,6 +771,7 @@ def copy(self, deep=False): other = cython.declare(Molecule) g = Graph.copy(self, deep) other = Molecule(g.vertices) + other.multiplicity = self.multiplicity return other def merge(self, other): @@ -954,6 +940,9 @@ def isIsomorphic(self, other, initialMap=None): # sufficient!) condition for the associated molecules to be isomorphic if self.getFingerprint() != other.getFingerprint(): return False + # check multiplicity + if self.multiplicity != other.multiplicity: + return False # Do the full isomorphism comparison result = Graph.isIsomorphic(self, other, initialMap) return result @@ -1012,7 +1001,8 @@ def isSubgraphIsomorphic(self, other, initialMap=None): carbonCount < group.carbonCount or nitrogenCount < group.nitrogenCount or oxygenCount < group.oxygenCount or - sulfurCount < group.sulfurCount): + sulfurCount < group.sulfurCount or + self.multiplicity not in group.multiplicity): return False # Do the isomorphism comparison @@ -1100,7 +1090,22 @@ def fromSMILES(self, smilesstr): # Special handling of helium if smilesstr == '[He]': # RDKit improperly handles helium and returns it in a triplet state - self.fromAdjacencyList('1 He 0 1') + self.fromAdjacencyList( + """ + He + multiplicity 1 + 1 He U0 L1 + """) + return self + elif smilesstr == 'C#O': + # carbon monoxide + self.fromAdjacencyList( + """ + CO + multiplicity 1 + 1 C U0 L1 {2,T} + 2 O U0 L1 {1,T} + """) return self else: @@ -1137,26 +1142,19 @@ def fromRDKitMol(self, rdkitmol): rdkitmol = Chem.AddHs(rdkitmol) Chem.rdmolops.Kekulize(rdkitmol, clearAromaticFlags=True) - # iterate though atoms in rdkitmol + # iterate through atoms in rdkitmol for i in range(rdkitmol.GetNumAtoms()): rdkitatom = rdkitmol.GetAtomWithIdx(i) # Use atomic number as key for element number = rdkitatom.GetAtomicNum() element = elements.getElement(number) - - # Process spin multiplicity - radicalElectrons = rdkitatom.GetNumRadicalElectrons() - - # Assume this is always true - # There are cases where 2 radicalElectrons is a singlet, but - # the triplet is often more stable, - spinMultiplicity = radicalElectrons + 1 # Process charge charge = rdkitatom.GetFormalCharge() + radicalElectrons = rdkitatom.GetNumRadicalElectrons() - atom = Atom(element, radicalElectrons, spinMultiplicity, charge, '', 0) + atom = Atom(element, radicalElectrons, charge, '', 0) self.vertices.append(atom) # Add bonds by iterating again through atoms @@ -1181,9 +1179,13 @@ def fromRDKitMol(self, rdkitmol): self.updateLonePairs() self.updateAtomTypes() + # Assume this is always true + # There are cases where 2 radicalElectrons is a singlet, but + # the triplet is often more stable, + self.multiplicity = self.getRadicalCount() + 1 + return self - newStyleAdjMatcher = re.compile('\s*\d+\s+(\*\d*\s+)?[A-Za-z]+\s+\S+\s+\d').match def fromAdjacencyList(self, adjlist, saturateH=False): """ Convert a string adjacency list `adjlist` to a molecular structure. @@ -1191,11 +1193,8 @@ def fromAdjacencyList(self, adjlist, saturateH=False): ``False``. """ from .adjlist import fromAdjacencyList - if not self.newStyleAdjMatcher(adjlist[adjlist.find('1 '):]): - logging.warning("There could be an old-style adjacency list. Please check library and input before proceeding!") - print adjlist - - self.vertices = fromAdjacencyList(adjlist, False, saturateH=saturateH) + + self.vertices, self.multiplicity = fromAdjacencyList(adjlist, group=False, saturateH=saturateH) self.updateConnectivityValues() self.updateAtomTypes() @@ -1445,12 +1444,12 @@ def toRDKitMol(self, removeHs=True, returnMapping=False): return rdkitmol, rdAtomIndices return rdkitmol - def toAdjacencyList(self, label='', removeH=False, removeLonePairs=False): + def toAdjacencyList(self, label='', removeH=False, removeLonePairs=False, printMultiplicity=True): """ Convert the molecular structure to a string adjacency list. """ from .adjlist import toAdjacencyList - result = toAdjacencyList(self.vertices, label=label, group=False, removeH=removeH, removeLonePairs=removeLonePairs) + result = toAdjacencyList(self.vertices, self.multiplicity, label=label, group=False, removeH=removeH, removeLonePairs=removeLonePairs, printMultiplicity=printMultiplicity) return result def isLinear(self): diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index 48f5b98cd0..d79898e2bb 100644 --- a/rmgpy/molecule/moleculeTest.py +++ b/rmgpy/molecule/moleculeTest.py @@ -18,7 +18,7 @@ def setUp(self): """ A method called before each unit test in this class. """ - self.atom = Atom(element=getElement('C'), radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + self.atom = Atom(element=getElement('C'), radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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, charge=0, label='*1') + atom = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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, charge=0, label='*1') + atom = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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, charge=0, label='*1') + atom = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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, charge=0, label='*1') + atom = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=2) if element.symbol == 'O': self.assertTrue(atom.isOxygen()) else: @@ -87,20 +87,16 @@ def testIncrementRadical(self): Test the Atom.incrementRadical() method. """ radicalElectrons = self.atom.radicalElectrons - spinMultiplicity = self.atom.spinMultiplicity self.atom.incrementRadical() self.assertEqual(self.atom.radicalElectrons, radicalElectrons + 1) - self.assertEqual(self.atom.spinMultiplicity, spinMultiplicity + 1) def testDecrementRadical(self): """ Test the Atom.decrementRadical() method. """ radicalElectrons = self.atom.radicalElectrons - spinMultiplicity = self.atom.spinMultiplicity self.atom.decrementRadical() self.assertEqual(self.atom.radicalElectrons, radicalElectrons - 1) - self.assertEqual(self.atom.spinMultiplicity, spinMultiplicity - 1) def testApplyActionBreakBond(self): """ @@ -108,12 +104,11 @@ def testApplyActionBreakBond(self): """ action = ['BREAK_BOND', '*1', 'S', '*2'] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -123,12 +118,11 @@ def testApplyActionFormBond(self): """ action = ['FORM_BOND', '*1', 'S', '*2'] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -138,12 +132,11 @@ def testApplyActionIncrementBond(self): """ action = ['CHANGE_BOND', '*1', 1, '*2'] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -153,12 +146,11 @@ def testApplyActionDecrementBond(self): """ action = ['CHANGE_BOND', '*1', -1, '*2'] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -168,12 +160,11 @@ def testApplyActionGainRadical(self): """ action = ['GAIN_RADICAL', '*1', 1] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -183,12 +174,11 @@ def testApplyActionLoseRadical(self): """ action = ['LOSE_RADICAL', '*1', 1] for element in elementList: - atom0 = Atom(element=element, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + atom0 = Atom(element=element, radicalElectrons=1, charge=0, label='*1', lonePairs=0) 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.charge, atom.charge) self.assertEqual(atom0.label, atom.label) @@ -198,8 +188,8 @@ def testEquivalent(self): """ for index1, element1 in enumerate(elementList[0:10]): for index2, element2 in enumerate(elementList[0:10]): - atom1 = Atom(element=element1, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') - atom2 = Atom(element=element2, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + atom1 = Atom(element=element1, radicalElectrons=1, charge=0, label='*1', lonePairs=0) + atom2 = Atom(element=element2, radicalElectrons=1, charge=0, label='*1', lonePairs=0) if index1 == index2: self.assertTrue(atom1.equivalent(atom2)) self.assertTrue(atom2.equivalent(atom1)) @@ -213,8 +203,8 @@ def testIsSpecificCaseOf(self): """ for index1, element1 in enumerate(elementList[0:10]): for index2, element2 in enumerate(elementList[0:10]): - atom1 = Atom(element=element1, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') - atom2 = Atom(element=element2, radicalElectrons=1, spinMultiplicity=2, charge=0, label='*1') + atom1 = Atom(element=element1, radicalElectrons=1, charge=0, label='*1', lonePairs=0) + atom2 = Atom(element=element2, radicalElectrons=1, charge=0, label='*1', lonePairs=0) if index1 == index2: self.assertTrue(atom1.isSpecificCaseOf(atom2)) else: @@ -228,7 +218,6 @@ def testCopy(self): self.assertEqual(self.atom.element.symbol, atom.element.symbol) self.assertEqual(self.atom.atomType, atom.atomType) self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons) - self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity) self.assertEqual(self.atom.charge, atom.charge) self.assertEqual(self.atom.label, atom.label) @@ -242,7 +231,6 @@ def testPickle(self): self.assertEqual(self.atom.element.symbol, atom.element.symbol) self.assertEqual(self.atom.atomType, atom.atomType) self.assertEqual(self.atom.radicalElectrons, atom.radicalElectrons) - self.assertEqual(self.atom.spinMultiplicity, atom.spinMultiplicity) self.assertEqual(self.atom.charge, atom.charge) self.assertEqual(self.atom.label, atom.label) @@ -468,45 +456,73 @@ class TestMolecule(unittest.TestCase): """ def setUp(self): - self.adjlist = """ -1 *2 C 1 0 {2,D} {3,S} -2 *1 O 0 2 {1,D} -3 C 0 0 {1,S} {4,S} {5,S} {6,S} -4 H 0 0 {3,S} -5 H 0 0 {3,S} -6 H 0 0 {3,S} + self.adjlist_1 = """ +1 *1 C U1 L0 E0 {2,S} {3,S} {4,S} +2 H U0 L0 E0 {1,S} +3 H U0 L0 E0 {1,S} +4 *2 N U0 L0 E+1 {1,S} {5,S} {6,D} +5 O U0 L3 E-1 {4,S} +6 O U0 L2 E0 {4,D} """ - self.molecule = Molecule().fromAdjacencyList(self.adjlist) + self.molecule = [Molecule().fromAdjacencyList(self.adjlist_1)] + + self.adjlist_2 = """ +1 *1 C U1 L0 {2,S} {3,S} {4,S} +2 H U0 L0 {1,S} +3 H U0 L0 {1,S} +4 *2 N U0 L0 {1,S} {5,S} {6,D} +5 O U0 L3 {4,S} +6 O U0 L2 {4,D} + """ + self.molecule.append(Molecule().fromAdjacencyList(self.adjlist_2)) + + self.adjlist_3 = """ +1 *1 C U1 {2,S} {3,S} {4,S} +2 H U0 {1,S} +3 H U0 {1,S} +4 *2 N U0 {1,S} {5,S} {6,D} +5 O U0 {4,S} +6 O U0 {4,D} + """ + self.molecule.append(Molecule().fromAdjacencyList(self.adjlist_3)) + + self.adjlist_4 = """ +1 *1 C U1 L0 {2,S} +2 *2 N U0 L0 {1,S} {3,S} {4,D} +3 O U0 L3 {2,S} +4 O U0 L2 {2,D} + """ + self.molecule.append(Molecule().fromAdjacencyList(self.adjlist_4,saturateH=True)) def testClearLabeledAtoms(self): """ Test the Molecule.clearLabeledAtoms() method. """ - self.molecule.clearLabeledAtoms() - for atom in self.molecule.atoms: + self.molecule[0].clearLabeledAtoms() + for atom in self.molecule[0].atoms: self.assertEqual(atom.label, '') def testContainsLabeledAtom(self): """ Test the Molecule.containsLabeledAtom() method. """ - for atom in self.molecule.atoms: + for atom in self.molecule[0].atoms: if atom.label != '': - self.assertTrue(self.molecule.containsLabeledAtom(atom.label)) - self.assertFalse(self.molecule.containsLabeledAtom('*3')) - self.assertFalse(self.molecule.containsLabeledAtom('*4')) - self.assertFalse(self.molecule.containsLabeledAtom('*5')) - self.assertFalse(self.molecule.containsLabeledAtom('*6')) + self.assertTrue(self.molecule[0].containsLabeledAtom(atom.label)) + self.assertFalse(self.molecule[0].containsLabeledAtom('*3')) + self.assertFalse(self.molecule[0].containsLabeledAtom('*4')) + self.assertFalse(self.molecule[0].containsLabeledAtom('*5')) + self.assertFalse(self.molecule[0].containsLabeledAtom('*6')) def testGetLabeledAtom(self): """ Test the Molecule.getLabeledAtom() method. """ - for atom in self.molecule.atoms: + for atom in self.molecule[0].atoms: if atom.label != '': - self.assertEqual(atom, self.molecule.getLabeledAtom(atom.label)) + self.assertEqual(atom, self.molecule[0].getLabeledAtom(atom.label)) try: - self.molecule.getLabeledAtom('*3') + self.molecule[0].getLabeledAtom('*3') self.fail('Unexpected successful return from Molecule.getLabeledAtom() with invalid atom label.') except ValueError: pass @@ -515,8 +531,8 @@ def testGetLabeledAtoms(self): """ Test the Molecule.getLabeledAtoms() method. """ - labeled = self.molecule.getLabeledAtoms() - for atom in self.molecule.atoms: + labeled = self.molecule[0].getLabeledAtoms() + for atom in self.molecule[0].atoms: if atom.label != '': self.assertTrue(atom.label in labeled) self.assertTrue(atom in labeled.values()) @@ -528,67 +544,220 @@ def testGetFormula(self): """ Test the Molecule.getLabeledAtoms() method. """ - self.assertEqual(self.molecule.getFormula(), 'C2H3O') + self.assertEqual(self.molecule[0].getFormula(), 'CH2NO2') + self.assertEqual(self.molecule[1].getFormula(), 'CH2NO2') + self.assertEqual(self.molecule[2].getFormula(), 'CH2NO2') + self.assertEqual(self.molecule[3].getFormula(), 'CH2NO2') def testRadicalCount(self): """ Test the Molecule.getRadicalCount() method. """ - self.assertEqual( self.molecule.getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule.atoms]) ) + self.assertEqual( self.molecule[0].getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule[0].atoms]) ) + self.assertEqual( self.molecule[1].getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule[1].atoms]) ) + self.assertEqual( self.molecule[2].getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule[2].atoms]) ) + self.assertEqual( self.molecule[3].getRadicalCount(), sum([atom.radicalElectrons for atom in self.molecule[3].atoms]) ) def testGetMolecularWeight(self): """ Test the Molecule.getMolecularWeight() method. """ - self.assertAlmostEqual(self.molecule.getMolecularWeight() * 1000, 43.04, 2) + self.assertAlmostEqual(self.molecule[0].getMolecularWeight() * 1000, 60.03, 2) + self.assertAlmostEqual(self.molecule[1].getMolecularWeight() * 1000, 60.03, 2) + self.assertAlmostEqual(self.molecule[2].getMolecularWeight() * 1000, 60.03, 2) + self.assertAlmostEqual(self.molecule[3].getMolecularWeight() * 1000, 60.03, 2) def testFromAdjacencyList(self): """ Test the Molecule.fromAdjacencyList() method. """ - 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 = atom1.bonds[atom2] - bond13 = atom1.bonds[atom3] + + # molecule 1 + + self.assertTrue(self.molecule[0].multiplicity == 2) + + atom1 = self.molecule[0].atoms[0] + atom2 = self.molecule[0].atoms[3] + atom3 = self.molecule[0].atoms[4] + atom4 = self.molecule[0].atoms[5] + self.assertTrue(self.molecule[0].hasBond(atom2,atom1)) + self.assertTrue(self.molecule[0].hasBond(atom2,atom3)) + self.assertTrue(self.molecule[0].hasBond(atom2,atom4)) + self.assertFalse(self.molecule[0].hasBond(atom1,atom3)) + self.assertFalse(self.molecule[0].hasBond(atom1,atom4)) + bond21 = atom2.bonds[atom1] + bond23 = atom2.bonds[atom3] + bond24 = atom2.bonds[atom4] - self.assertTrue(atom1.label == '*2') + self.assertTrue(atom1.label == '*1') self.assertTrue(atom1.element.symbol == 'C') self.assertTrue(atom1.radicalElectrons == 1) - self.assertTrue(atom1.spinMultiplicity == 2) self.assertTrue(atom1.charge == 0) - self.assertTrue(atom2.label == '*1') - self.assertTrue(atom2.element.symbol == 'O') + self.assertTrue(atom2.label == '*2') + self.assertTrue(atom2.element.symbol == 'N') self.assertTrue(atom2.radicalElectrons == 0) - self.assertTrue(atom2.spinMultiplicity == 1) - self.assertTrue(atom2.charge == 0) + self.assertTrue(atom2.charge == 1) self.assertTrue(atom3.label == '') - self.assertTrue(atom3.element.symbol == 'C') + self.assertTrue(atom3.element.symbol == 'O') self.assertTrue(atom3.radicalElectrons == 0) - self.assertTrue(atom3.spinMultiplicity == 1) - self.assertTrue(atom3.charge == 0) + self.assertTrue(atom3.charge == -1) + + self.assertTrue(atom4.label == '') + self.assertTrue(atom4.element.symbol == 'O') + self.assertTrue(atom4.radicalElectrons == 0) + self.assertTrue(atom4.charge == 0) - self.assertTrue(bond12.isDouble()) - self.assertTrue(bond13.isSingle()) + self.assertTrue(bond21.isSingle()) + self.assertTrue(bond23.isSingle()) + self.assertTrue(bond24.isDouble()) + + + # molecule 2 + + self.assertTrue(self.molecule[1].multiplicity == 2) + + atom1 = self.molecule[1].atoms[0] + atom2 = self.molecule[1].atoms[3] + atom3 = self.molecule[1].atoms[4] + atom4 = self.molecule[1].atoms[5] + self.assertTrue(self.molecule[1].hasBond(atom2,atom1)) + self.assertTrue(self.molecule[1].hasBond(atom2,atom3)) + self.assertTrue(self.molecule[1].hasBond(atom2,atom4)) + self.assertFalse(self.molecule[1].hasBond(atom1,atom3)) + self.assertFalse(self.molecule[1].hasBond(atom1,atom4)) + bond21 = atom2.bonds[atom1] + bond23 = atom2.bonds[atom3] + bond24 = atom2.bonds[atom4] + + self.assertTrue(atom1.label == '*1') + self.assertTrue(atom1.element.symbol == 'C') + self.assertTrue(atom1.radicalElectrons == 1) + self.assertTrue(atom1.charge == 0) + + self.assertTrue(atom2.label == '*2') + self.assertTrue(atom2.element.symbol == 'N') + self.assertTrue(atom2.radicalElectrons == 0) + self.assertTrue(atom2.charge == 1) + + self.assertTrue(atom3.label == '') + self.assertTrue(atom3.element.symbol == 'O') + self.assertTrue(atom3.radicalElectrons == 0) + self.assertTrue(atom3.charge == -1) + + self.assertTrue(atom4.label == '') + self.assertTrue(atom4.element.symbol == 'O') + self.assertTrue(atom4.radicalElectrons == 0) + self.assertTrue(atom4.charge == 0) + + self.assertTrue(bond21.isSingle()) + self.assertTrue(bond23.isSingle()) + self.assertTrue(bond24.isDouble()) + + + # molecule 3 + + self.assertTrue(self.molecule[2].multiplicity == 2) + + atom1 = self.molecule[2].atoms[0] + atom2 = self.molecule[2].atoms[3] + atom3 = self.molecule[2].atoms[4] + atom4 = self.molecule[2].atoms[5] + self.assertTrue(self.molecule[2].hasBond(atom2,atom1)) + self.assertTrue(self.molecule[2].hasBond(atom2,atom3)) + self.assertTrue(self.molecule[2].hasBond(atom2,atom4)) + self.assertFalse(self.molecule[2].hasBond(atom1,atom3)) + self.assertFalse(self.molecule[2].hasBond(atom1,atom4)) + bond21 = atom2.bonds[atom1] + bond23 = atom2.bonds[atom3] + bond24 = atom2.bonds[atom4] + + self.assertTrue(atom1.label == '*1') + self.assertTrue(atom1.element.symbol == 'C') + self.assertTrue(atom1.radicalElectrons == 1) + self.assertTrue(atom1.charge == 0) + + self.assertTrue(atom2.label == '*2') + self.assertTrue(atom2.element.symbol == 'N') + self.assertTrue(atom2.radicalElectrons == 0) + self.assertTrue(atom2.charge == 1) + + self.assertTrue(atom3.label == '') + self.assertTrue(atom3.element.symbol == 'O') + self.assertTrue(atom3.radicalElectrons == 0) + self.assertTrue(atom3.charge == -1) + + self.assertTrue(atom4.label == '') + self.assertTrue(atom4.element.symbol == 'O') + self.assertTrue(atom4.radicalElectrons == 0) + self.assertTrue(atom4.charge == 0) + + self.assertTrue(bond21.isSingle()) + self.assertTrue(bond23.isSingle()) + self.assertTrue(bond24.isDouble()) + + + # molecule 4 + + self.assertTrue(self.molecule[3].multiplicity == 2) + + atom1 = self.molecule[3].atoms[0] + atom2 = self.molecule[3].atoms[1] + atom3 = self.molecule[3].atoms[2] + atom4 = self.molecule[3].atoms[3] + self.assertTrue(self.molecule[3].hasBond(atom2,atom1)) + self.assertTrue(self.molecule[3].hasBond(atom2,atom3)) + self.assertTrue(self.molecule[3].hasBond(atom2,atom4)) + self.assertFalse(self.molecule[3].hasBond(atom1,atom3)) + self.assertFalse(self.molecule[3].hasBond(atom1,atom4)) + bond21 = atom2.bonds[atom1] + bond23 = atom2.bonds[atom3] + bond24 = atom2.bonds[atom4] + + self.assertTrue(atom1.label == '*1') + self.assertTrue(atom1.element.symbol == 'C') + self.assertTrue(atom1.radicalElectrons == 1) + self.assertTrue(atom1.charge == 0) + + self.assertTrue(atom2.label == '*2') + self.assertTrue(atom2.element.symbol == 'N') + self.assertTrue(atom2.radicalElectrons == 0) + self.assertTrue(atom2.charge == 1) + + self.assertTrue(atom3.label == '') + self.assertTrue(atom3.element.symbol == 'O') + self.assertTrue(atom3.radicalElectrons == 0) + self.assertTrue(atom3.charge == -1) + + self.assertTrue(atom4.label == '') + self.assertTrue(atom4.element.symbol == 'O') + self.assertTrue(atom4.radicalElectrons == 0) + self.assertTrue(atom4.charge == 0) + + self.assertTrue(bond21.isSingle()) + self.assertTrue(bond23.isSingle()) + self.assertTrue(bond24.isDouble()) + + def testToAdjacencyList(self): """ Test the Molecule.toAdjacencyList() method. """ - adjlist = self.molecule.toAdjacencyList(removeH=False) - self.assertEqual(adjlist.strip(), self.adjlist.strip()) + adjlist_1 = self.molecule[0].toAdjacencyList(removeH=False,printMultiplicity=False) + self.assertEqual(adjlist_1.strip(), self.adjlist_1.strip()) - def testFromOldAdjacencyList(self): - """ - Test we can read things with implicit hydrogens. - """ - adjList = "1 O 0" # should be Water - molecule = Molecule().fromAdjacencyList(adjList, saturateH=True) # only works with saturateH=True - self.assertEqual(molecule.getFormula(),'H2O') +# def testFromOldAdjacencyList(self): +# """ +# Test we can read things with implicit hydrogens. +# """ +# adjList = """ +# 1 O 0 +# """ # should be Water +# molecule = Molecule().fromAdjacencyList(adjList, saturateH=True) # only works with saturateH=True +# self.assertEqual(molecule.getFormula(),'H2O') def testIsomorphism(self): """ @@ -605,8 +774,8 @@ def testSubgraphIsomorphism(self): """ molecule = Molecule().fromSMILES('C=CC=C[CH]C') group = Group().fromAdjacencyList(""" - 1 Cd 0 {2,D} - 2 Cd 0 {1,D} + 1 Cd U0 {2,D} + 2 Cd U0 {1,D} """) self.assertTrue(molecule.isSubgraphIsomorphic(group)) @@ -621,30 +790,30 @@ def testSubgraphIsomorphism(self): def testSubgraphIsomorphismAgain(self): molecule = Molecule() molecule.fromAdjacencyList(""" - 1 * C 0 {2,D} {7,S} {8,S} - 2 C 0 {1,D} {3,S} {9,S} - 3 C 0 {2,S} {4,D} {10,S} - 4 C 0 {3,D} {5,S} {11,S} - 5 C 0 {4,S} {6,S} {12,S} {13,S} - 6 C 0 {5,S} {14,S} {15,S} {16,S} - 7 H 0 {1,S} - 8 H 0 {1,S} - 9 H 0 {2,S} - 10 H 0 {3,S} - 11 H 0 {4,S} - 12 H 0 {5,S} - 13 H 0 {5,S} - 14 H 0 {6,S} - 15 H 0 {6,S} - 16 H 0 {6,S} + 1 * C U0 {2,D} {7,S} {8,S} + 2 C U0 {1,D} {3,S} {9,S} + 3 C U0 {2,S} {4,D} {10,S} + 4 C U0 {3,D} {5,S} {11,S} + 5 C U0 {4,S} {6,S} {12,S} {13,S} + 6 C U0 {5,S} {14,S} {15,S} {16,S} + 7 H U0 {1,S} + 8 H U0 {1,S} + 9 H U0 {2,S} + 10 H U0 {3,S} + 11 H U0 {4,S} + 12 H U0 {5,S} + 13 H U0 {5,S} + 14 H U0 {6,S} + 15 H U0 {6,S} + 16 H U0 {6,S} """) group = Group() group.fromAdjacencyList(""" - 1 * C 0 {2,D} {3,S} {4,S} - 2 C 0 {1,D} - 3 H 0 {1,S} - 4 H 0 {1,S} + 1 * C U0 {2,D} {3,S} {4,S} + 2 C U0 {1,D} + 3 H U0 {1,S} + 4 H U0 {1,S} """) labeled1 = molecule.getLabeledAtoms().values()[0] @@ -665,22 +834,21 @@ def testSubgraphIsomorphismAgain(self): def testSubgraphIsomorphismManyLabels(self): molecule = Molecule() # specific case (species) molecule.fromAdjacencyList(""" -1 *1 C 1 {2,S} {3,S} {4,S} -2 C 0 {1,S} {3,S} {5,S} {6,S} -3 C 0 {1,S} {2,S} {7,S} {8,S} -4 H 0 {1,S} -5 H 0 {2,S} -6 H 0 {2,S} -7 H 0 {3,S} -8 H 0 {3,S} +1 *1 C U1 {2,S} {3,S} {4,S} +2 C U0 {1,S} {3,S} {5,S} {6,S} +3 C U0 {1,S} {2,S} {7,S} {8,S} +4 H U0 {1,S} +5 H U0 {2,S} +6 H U0 {2,S} +7 H U0 {3,S} +8 H U0 {3,S} """) - print molecule.toAdjacencyList() group = Group() # general case (functional group) group.fromAdjacencyList(""" -1 *1 C 1 {2,S}, {3,S} -2 R!H 0 {1,S} -3 R!H 0 {1,S} +1 *1 C U1 {2,S}, {3,S} +2 R!H U0 {1,S} +3 R!H U0 {1,S} """) labeled1 = molecule.getLabeledAtoms() @@ -703,21 +871,21 @@ def testAdjacencyList(self): Check the adjacency list read/write functions for a full molecule. """ molecule1 = Molecule().fromAdjacencyList(""" - 1 C 0 {2,D} {7,S} {8,S} - 2 C 0 {1,D} {3,S} {9,S} - 3 C 0 {2,S} {4,D} {10,S} - 4 C 0 {3,D} {5,S} {11,S} - 5 C 1 {4,S} {6,S} {12,S} - 6 C 0 {5,S} {13,S} {14,S} {15,S} - 7 H 0 {1,S} - 8 H 0 {1,S} - 9 H 0 {2,S} - 10 H 0 {3,S} - 11 H 0 {4,S} - 12 H 0 {5,S} - 13 H 0 {6,S} - 14 H 0 {6,S} - 15 H 0 {6,S} + 1 C U0 {2,D} {7,S} {8,S} + 2 C U0 {1,D} {3,S} {9,S} + 3 C U0 {2,S} {4,D} {10,S} + 4 C U0 {3,D} {5,S} {11,S} + 5 C U1 {4,S} {6,S} {12,S} + 6 C U0 {5,S} {13,S} {14,S} {15,S} + 7 H U0 {1,S} + 8 H U0 {1,S} + 9 H U0 {2,S} + 10 H U0 {3,S} + 11 H U0 {4,S} + 12 H U0 {5,S} + 13 H U0 {6,S} + 14 H U0 {6,S} + 15 H U0 {6,S} """) molecule2 = Molecule().fromSMILES('C=CC=C[CH]C') self.assertTrue(molecule1.isIsomorphic(molecule2)) @@ -806,7 +974,7 @@ def testRadicalCH(self): """ molecule = Molecule().fromSMILES('[CH]') self.assertEqual(molecule.atoms[0].radicalElectrons, 3) - self.assertEqual(molecule.atoms[0].spinMultiplicity, 4) + self.assertEqual(molecule.multiplicity, 4) self.assertEqual(molecule.getRadicalCount(), 3) def testRadicalCH2(self): @@ -815,7 +983,7 @@ def testRadicalCH2(self): """ molecule = Molecule().fromSMILES('[CH2]') self.assertEqual(molecule.atoms[0].radicalElectrons, 2) - self.assertEqual(molecule.atoms[0].spinMultiplicity, 3) + self.assertEqual(molecule.multiplicity, 3) self.assertEqual(molecule.getRadicalCount(), 2) def testRadicalCH2CH2CH2(self): @@ -830,7 +998,7 @@ def testSMILES(self): Test that we can generate a few SMILES strings as expected """ import rmgpy.molecule - test_strings =['CO', '[C]', '[CH]', 'OO', '[H][H]', '[H]', '[He]', '[O]', 'O', '[CH3]', 'C', '[OH]', 'CCC', 'CC', 'N#N', '[O]O', 'C[CH2]', '[Ar]', 'CCCC','O=C=O','CN'] + test_strings =['C#O', '[C]', '[CH]', 'OO', '[H][H]', '[H]', '[He]', '[O]', 'O', '[CH3]', 'C', '[OH]', 'CCC', 'CC', 'N#N', '[O]O', 'C[CH2]', '[Ar]', 'CCCC','O=C=O','[C]#N'] for s in test_strings: molecule = Molecule(SMILES=s) self.assertEqual(s,molecule.toSMILES()) @@ -848,8 +1016,8 @@ def testAugmentedInChI(self): Test that Augmented InChI generation is printing the /mult layer """ mol = Molecule().fromAdjacencyList(""" - 1 C 1 {2,S} - 2 C 1 {1,S} + 1 C U1 {2,S} + 2 C U1 {1,S} """) self.assertEqual(mol.toAugmentedInChI(), 'InChI=1S/C2H4/c1-2/h1-2H2/mult3') @@ -859,8 +1027,8 @@ def testAugmentedInChIKey(self): Test that Augmented InChI Key generation is printing the mult layer """ mol = Molecule().fromAdjacencyList(""" - 1 C 1 {2,S} - 2 C 1 {1,S} + 1 C U1 {2,S} + 2 C U1 {1,S} """) self.assertEqual(mol.toAugmentedInChIKey(), 'VGGSQFUCUMXWEO-UHFFFAOYSAmult3') diff --git a/rmgpy/pdep/networkTest.py b/rmgpy/pdep/networkTest.py index 51f017d905..475d7432ed 100644 --- a/rmgpy/pdep/networkTest.py +++ b/rmgpy/pdep/networkTest.py @@ -57,6 +57,7 @@ def setUp(self): """ self.nC4H10O = Species( label = 'n-C4H10O', + multiplicity = 1, conformer = Conformer( E0 = (-317.807,'kJ/mol'), modes = [ @@ -78,6 +79,7 @@ def setUp(self): self.nC4H8 = Species( label = 'n-C4H8', + multiplicity = 1, conformer = Conformer( E0 = (-17.8832,'kJ/mol'), modes = [ @@ -94,6 +96,7 @@ def setUp(self): self.H2O = Species( label = 'H2O', + multiplicity = 1, conformer = Conformer( E0 = (-269.598,'kJ/mol'), modes = [ @@ -108,6 +111,7 @@ def setUp(self): self.N2 = Species( label = 'N2', + multiplicity = 1, molecularWeight = (28.04,"g/mol"), transportData=TransportData(sigma=(3.41, "angstrom"), epsilon=(124, "K")), energyTransferModel = None, diff --git a/rmgpy/qm/molecule.py b/rmgpy/qm/molecule.py index c318a4407f..f63390e103 100644 --- a/rmgpy/qm/molecule.py +++ b/rmgpy/qm/molecule.py @@ -42,7 +42,7 @@ def __init__(self, settings, uniqueID, molecule, multiplicity, uniqueIDlong=None #: A short unique ID such as an augmented InChI Key. self.uniqueID = uniqueID self.molecule = molecule - #: The multiplicity, eg. the number of free radicals plus one. + #: The multiplicity self.multiplicity = multiplicity if uniqueIDlong is None: self.uniqueIDlong = uniqueID @@ -236,7 +236,7 @@ def createGeometry(self): """ Creates self.geometry with RDKit geometries """ - multiplicity = sum([i.radicalElectrons for i in self.molecule.atoms]) + 1 + multiplicity = self.molecule.multiplicity self.geometry = Geometry(self.settings, self.uniqueID, self.molecule, multiplicity, uniqueIDlong=self.uniqueIDlong) self.geometry.generateRDKitGeometries() return self.geometry diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 8eeb90797f..230852e40a 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -100,6 +100,8 @@ cdef class Reaction: cpdef generatePairs(self) + cpdef copy(self) + ################################################################################ cdef class ReactionModel: diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 5edb29fb98..1fae86560a 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -45,6 +45,7 @@ import logging import re import os.path +from copy import copy, deepcopy import rmgpy.constants as constants from rmgpy.molecule.molecule import Molecule, Atom @@ -980,6 +981,31 @@ def generate3dTS(self, reactants, products): zCoord[k] = dirVec[k].z*lenVec[k] reactionAxis = [sum(xCoord), sum(yCoord), sum(zCoord)] products[i].reactionAxis = reactionAxis + + def copy(self): + """ + Create a deep copy of the current reaction. + """ + + cython.declare(other=Reaction) + + other = Reaction.__new__(Reaction) + other.index = self.index + other.label = self.label + other.reactants = [] + for reactant in self.reactants: + other.reactants.append(reactant.copy(deep=True)) + other.products = [] + for product in self.products: + other.products.append(product.copy(deep=True)) + other.kinetics = deepcopy(self.kinetics) + other.reversible = self.reversible + other.transitionState = deepcopy(self.transitionState) + other.duplicate = self.duplicate + other.degeneracy = self.degeneracy + other.pairs = deepcopy(self.pairs) + + return other ################################################################################ diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index 0f276e3b52..7a295e2af0 100644 --- a/rmgpy/reactionTest.py +++ b/rmgpy/reactionTest.py @@ -90,6 +90,7 @@ def setUp(self): """ ethylene = Species( label = 'C2H4', + multiplicity = 1, conformer = Conformer( E0 = (44.7127, 'kJ/mol'), modes = [ @@ -117,6 +118,7 @@ def setUp(self): hydrogen = Species( label = 'H', + multiplicity = 2, conformer = Conformer( E0 = (211.794, 'kJ/mol'), modes = [ @@ -131,6 +133,7 @@ def setUp(self): ethyl = Species( label = 'C2H5', + multiplicity = 2, conformer = Conformer( E0 = (111.603, 'kJ/mol'), modes = [ diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 87e24e7b21..474435227f 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -88,9 +88,9 @@ def database( assert isinstance(kineticsFamilies,list), "kineticsFamilies should be either 'default', 'all', 'none', or a list of names eg. ['H_Abstraction','R_Recombination'] or ['!Intra_Disproportionation']." rmg.kineticsFamilies = kineticsFamilies -def species(label, structure, reactive=True): +def species(label, multiplicity, structure, reactive=True): logging.debug('Found {0} species "{1}" ({2})'.format('reactive' if reactive else 'nonreactive', label, structure.toSMILES())) - spec, isNew = rmg.reactionModel.makeNewSpecies(structure, label=label, reactive=reactive) + spec, isNew = rmg.reactionModel.makeNewSpecies(structure, multiplicity, label=label, reactive=reactive) assert isNew, "Species {0} is a duplicate of {1}. Species in input file must be unique".format(label,spec.label) rmg.initialSpecies.append(spec) speciesDict[label] = spec diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 1e252a604b..60ec90ec58 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -61,6 +61,9 @@ from pdep import PDepReaction, PDepNetwork, PressureDependenceError # generateThermoDataFromQM under the Species class imports the qm package +#: This dictionary is used to add multiplicity to species label +_multiplicity_labels = {1:'S',2:'D',3:'T',4:'Q',5:'V',} + ################################################################################ @@ -70,11 +73,11 @@ class Species(rmgpy.species.Species): solventViscosity = None diffusionTemp = None - def __init__(self, index=-1, label='', thermo=None, conformer=None, + def __init__(self, index=-1, label='', multiplicity=103, thermo=None, conformer=None, molecule=None, transportData=None, molecularWeight=None, dipoleMoment=None, polarizability=None, Zrot=None, energyTransferModel=None, reactive=True, coreSizeAtCreation=0): - rmgpy.species.Species.__init__(self, index, label, thermo, conformer, molecule, transportData, molecularWeight, dipoleMoment, polarizability, Zrot, energyTransferModel, reactive) + rmgpy.species.Species.__init__(self, index, label, multiplicity, thermo, conformer, molecule, transportData, molecularWeight, dipoleMoment, polarizability, Zrot, energyTransferModel, reactive) self.coreSizeAtCreation = coreSizeAtCreation def __reduce__(self): @@ -151,7 +154,8 @@ def generateThermoData(self, database, thermoClass=NASA, quantumMechanics=None): if thermo0 is not None: # Write the QM molecule thermo to a library so that can be used in future RMG jobs. quantumMechanics.database.loadEntry(index = len(quantumMechanics.database.entries) + 1, - label = molecule.toSMILES(), + label = molecule.toSMILES() + '_({0})'.format(_multiplicity_labels[molecule.multiplicity]), + multiplicity = molecule.multiplicity, molecule = molecule.toAdjacencyList(), thermo = thermo0, shortDesc = thermo0.comment @@ -372,7 +376,7 @@ def checkForExistingSpecies(self, molecule): # At this point we can conclude that the structure does not exist return False, None - def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True): + def makeNewSpecies(self, object, multiplicity, label='', reactive=True, checkForExisting=True): """ Formally create a new species from the specified `object`, which can be either a :class:`Molecule` object or an :class:`rmgpy.species.Species` @@ -386,6 +390,7 @@ def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True) molecule = object molecule.clearLabeledAtoms() + molecule.multiplicity = multiplicity # If desired, check to ensure that the species is new; return the # existing species if not new @@ -404,7 +409,7 @@ def makeNewSpecies(self, object, label='', reactive=True, checkForExisting=True) # so that we can use the label in file paths label = molecule.toSMILES().replace('/','').replace('\\','') logging.debug('Creating new species {0}'.format(label)) - spec = Species(index=self.speciesCounter+1, label=label, molecule=[molecule], reactive=reactive) + spec = Species(index=self.speciesCounter+1, label=label, multiplicity=multiplicity, molecule=[molecule], reactive=reactive) spec.coreSizeAtCreation = len(self.core.species) spec.generateResonanceIsomers() spec.molecularWeight = Quantity(spec.molecule[0].getMolecularWeight()*1000.,"amu") @@ -513,8 +518,8 @@ def makeNewReaction(self, forward, checkExisting=True): """ # Determine the proper species objects for all reactants and products - reactants = [self.makeNewSpecies(reactant)[0] for reactant in forward.reactants] - products = [self.makeNewSpecies(product)[0] for product in forward.products ] + reactants = [self.makeNewSpecies(reactant,reactant.multiplicity,)[0] for reactant in forward.reactants] + products = [self.makeNewSpecies(product,product.multiplicity,)[0] for product in forward.products ] if forward.pairs is not None: for pairIndex in range(len(forward.pairs)): reactantIndex = forward.reactants.index(forward.pairs[pairIndex][0]) diff --git a/rmgpy/solver/simpleTest.py b/rmgpy/solver/simpleTest.py index 6864c805df..8e1ac6bb47 100644 --- a/rmgpy/solver/simpleTest.py +++ b/rmgpy/solver/simpleTest.py @@ -26,18 +26,22 @@ def testSolve(self): CH4 + C2H5 <=> CH3 + C2H6. """ CH4 = Species( + multiplicity=1, molecule=[Molecule().fromSMILES("C")], thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],"K"), Cpdata=([ 8.615, 9.687,10.963,12.301,14.841,16.976,20.528],"cal/(mol*K)"), H298=(-17.714,"kcal/mol"), S298=(44.472,"cal/(mol*K)")) ) CH3 = Species( + multiplicity=2, molecule=[Molecule().fromSMILES("[CH3]")], thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],"K"), Cpdata=([ 9.397,10.123,10.856,11.571,12.899,14.055,16.195],"cal/(mol*K)"), H298=( 9.357,"kcal/mol"), S298=(45.174,"cal/(mol*K)")) ) C2H6 = Species( + multiplicity=1, molecule=[Molecule().fromSMILES("CC")], thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],"K"), Cpdata=([12.684,15.506,18.326,20.971,25.500,29.016,34.595],"cal/(mol*K)"), H298=(-19.521,"kcal/mol"), S298=(54.799,"cal/(mol*K)")) ) C2H5 = Species( + multiplicity=2, molecule=[Molecule().fromSMILES("C[CH2]")], thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],"K"), Cpdata=([11.635,13.744,16.085,18.246,21.885,24.676,29.107],"cal/(mol*K)"), H298=( 29.496,"kcal/mol"), S298=(56.687,"cal/(mol*K)")) ) @@ -90,6 +94,7 @@ def testSolve(self): # Solve a reaction system and check if the analytical jacobian matches the finite difference jacobian H2 = Species( + multiplicity = 1, molecule=[Molecule().fromSMILES("[H][H]")], thermo=ThermoData(Tdata=([300,400,500,600,800,1000,1500],"K"), Cpdata=([6.89,6.97,6.99,7.01,7.08,7.22,7.72],"cal/(mol*K)"), H298=( 0,"kcal/mol"), S298=(31.23,"cal/(mol*K)")) ) diff --git a/rmgpy/species.pxd b/rmgpy/species.pxd index 6080745d76..8e2fa9f200 100644 --- a/rmgpy/species.pxd +++ b/rmgpy/species.pxd @@ -38,6 +38,7 @@ cdef class Species: cdef public int index cdef public str label + cdef public int multiplicity cdef public HeatCapacityModel thermo cdef public Conformer conformer cdef public object transportData diff --git a/rmgpy/species.py b/rmgpy/species.py index 2ec3caf696..8ae80ac36c 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -49,6 +49,9 @@ import rmgpy.quantity as quantity from rmgpy.molecule import Molecule +#: This dictionary is used to add multiplicity to species label +_multiplicity_labels = {1:'S',2:'D',3:'T',4:'Q',5:'V',} + ################################################################################ class SpeciesError(Exception): @@ -71,6 +74,7 @@ class Species(object): ======================= ==================================================== `index` A unique nonnegative integer index `label` A descriptive string label + `multiplicity` The multiplicity of this species, integer multiplicity = 2*total_spin+1 `thermo` The heat capacity model for the species `conformer` The molecular conformer for the species `molecule` A list of the :class:`Molecule` objects describing the molecular structure @@ -88,12 +92,17 @@ class Species(object): :class:`rmg.model.Species` inherits from this class, and adds some extra methods. """ - def __init__(self, index=-1, label='', thermo=None, conformer=None, + def __init__(self, index=-1, label='', multiplicity=102, thermo=None, conformer=None, molecule=None, transportData=None, molecularWeight=None, dipoleMoment=None, polarizability=None, Zrot=None, energyTransferModel=None, reactive=True): self.index = index + if label != '' and molecule is not None: + if molecule != []: + if molecule[0].getRadicalCount()>1: + if '_' not in label: label += '_({0})' .format(_multiplicity_labels[multiplicity]) self.label = label + self.multiplicity = multiplicity self.thermo = thermo self.conformer = conformer self.molecule = molecule or [] @@ -104,6 +113,14 @@ def __init__(self, index=-1, label='', thermo=None, conformer=None, self.polarizability = polarizability self.Zrot = Zrot self.energyTransferModel = energyTransferModel + + # Check if multiplicity is possible + if molecule is not None: + if molecule != []: + n_rad = molecule[0].getRadicalCount() + if not (n_rad + 1 == multiplicity or n_rad - 1 == multiplicity or n_rad - 3 == multiplicity or n_rad - 5 == multiplicity): + print molecule[0].toAdjacencyList() + raise SpeciesError('Impossible multiplicity for {0}: multiplicity = {1} and number of unpaired electrons = {2}'.format(label,multiplicity,n_rad)) def __repr__(self): """ @@ -113,6 +130,7 @@ def __repr__(self): string = 'Species(' if self.index != -1: string += 'index={0:d}, '.format(self.index) if self.label != -1: string += 'label="{0}", '.format(self.label) + if self.multiplicity != -1: string += 'multiplicity={0}, '.format(self.multiplicity) if self.thermo is not None: string += 'thermo={0!r}, '.format(self.thermo) if self.conformer is not None: string += 'conformer={0!r}, '.format(self.conformer) if len(self.molecule) > 0: string += 'molecule=[{0!r}], '.format(self.molecule[0]) @@ -143,7 +161,7 @@ def __reduce__(self): """ A helper function used when pickling an object. """ - return (Species, (self.index, self.label, self.thermo, self.conformer, self.molecule, self.transportData, self.molecularWeight, self.dipoleMoment, self.polarizability, self.Zrot, self.energyTransferModel, self.reactive)) + return (Species, (self.index, self.label, self.multiplicity, self.thermo, self.conformer, self.molecule, self.transportData, self.molecularWeight, self.dipoleMoment, self.polarizability, self.Zrot, self.energyTransferModel, self.reactive)) def getMolecularWeight(self): return self._molecularWeight @@ -189,10 +207,11 @@ def isIsomorphic(self, other): if molecule.isIsomorphic(other): return True elif isinstance(other, Species): - for molecule1 in self.molecule: - for molecule2 in other.molecule: - if molecule1.isIsomorphic(molecule2): - return True + if self.multiplicity == other.multiplicity: + for molecule1 in self.molecule: + for molecule2 in other.molecule: + if molecule1.isIsomorphic(molecule2): + return True else: raise ValueError('Unexpected value "{0!r}" for other parameter; should be a Molecule or Species object.'.format(other)) return False diff --git a/rmgpy/speciesTest.py b/rmgpy/speciesTest.py index f20020d6e6..2d9d9038be 100644 --- a/rmgpy/speciesTest.py +++ b/rmgpy/speciesTest.py @@ -27,6 +27,7 @@ def setUp(self): self.species = Species( index=1, label='C2H4', + multiplicity = 1, thermo=ThermoData( Tdata=([300.0,400.0,500.0,600.0,800.0,1000.0,1500.0],'K'), Cpdata=([3.0,4.0,5.0,6.0,8.0,10.0,15.0],'cal/(mol*K)'), @@ -61,6 +62,7 @@ def testPickle(self): species = cPickle.loads(cPickle.dumps(self.species,-1)) self.assertEqual(self.species.index, species.index) self.assertEqual(self.species.label, species.label) + self.assertEqual(self.species.multiplicity, species.multiplicity) self.assertEqual(self.species.thermo.H298.value_si, species.thermo.H298.value_si) self.assertEqual(self.species.thermo.H298.units, species.thermo.H298.units) self.assertEqual(len(self.species.conformer.modes), len(species.conformer.modes)) @@ -86,6 +88,7 @@ def testOutput(self): exec('species = {0!r}'.format(self.species)) self.assertEqual(self.species.index, species.index) self.assertEqual(self.species.label, species.label) + self.assertEqual(self.species.multiplicity, species.multiplicity) self.assertEqual(self.species.thermo.H298.value_si, species.thermo.H298.value_si) self.assertEqual(self.species.thermo.H298.units, species.thermo.H298.units) self.assertEqual(len(self.species.conformer.modes), len(species.conformer.modes)) diff --git a/rmgpy/transportDataTest.py b/rmgpy/transportDataTest.py index 12981e4c2e..d8191f55c3 100644 --- a/rmgpy/transportDataTest.py +++ b/rmgpy/transportDataTest.py @@ -252,7 +252,8 @@ def testJoback(self): ['benzene', 'c1ccccc1', None, None, None], ] for name, smiles, sigma, epsilon, comment in self.testCases: - species = Species(molecule=[Molecule(SMILES=smiles)]) + molecule=Molecule(SMILES=smiles) + species = Species(multiplicity=molecule.multiplicity,molecule=[molecule]) transportData, blank, blank2 = self.transportdb.getTransportPropertiesViaGroupEstimates(species) # check Joback worked. # If we don't know what to expect, don't check (just make sure we didn't crash) @@ -266,21 +267,21 @@ def testJoback(self): def testJobackOnBenzeneBonds(self): "Test Joback doesn't crash on Cb desription of beneze" adjlist = """ - 1 C 0 0 {2,D} {6,S} {7,S} - 2 C 0 0 {1,D} {3,S} {8,S} - 3 C 0 0 {2,S} {4,D} {9,S} - 4 C 0 0 {3,D} {5,S} {10,S} - 5 C 0 0 {4,S} {6,D} {11,S} - 6 C 0 0 {1,S} {5,D} {12,S} - 7 H 0 0 {1,S} - 8 H 0 0 {2,S} - 9 H 0 0 {3,S} - 10 H 0 0 {4,S} - 11 H 0 0 {5,S} - 12 H 0 0 {6,S} + 1 C U0 L0 {2,D} {6,S} {7,S} + 2 C U0 L0 {1,D} {3,S} {8,S} + 3 C U0 L0 {2,S} {4,D} {9,S} + 4 C U0 L0 {3,D} {5,S} {10,S} + 5 C U0 L0 {4,S} {6,D} {11,S} + 6 C U0 L0 {1,S} {5,D} {12,S} + 7 H U0 L0 {1,S} + 8 H U0 L0 {2,S} + 9 H U0 L0 {3,S} + 10 H U0 L0 {4,S} + 11 H U0 L0 {5,S} + 12 H U0 L0 {6,S} """ m = Molecule().fromAdjacencyList(adjlist) - species = Species(molecule=[m]) + species = Species(multiplicity=m.multiplicity, molecule=[m]) transportData, blank, blank2 = self.transportdb.getTransportPropertiesViaGroupEstimates(species) self.assertIsNotNone(transportData)