From d7f7fcbd7ea7c975d53c1a82cf14dc31732b3996 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 13:55:04 -0400 Subject: [PATCH 01/52] Update RMG example 1,3-hexadiene with multiplicity label --- examples/rmg/1,3-hexadiene/input.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/examples/rmg/1,3-hexadiene/input.py b/examples/rmg/1,3-hexadiene/input.py index 28615b51a8..e73b25ba81 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 0 0 {2,S} + 2 H 0 0 {1,S} """), ) species( label='N2', + multiplicity = 1, reactive=False, structure=InChI("InChI=1/N2/c1-2"), ) From 302f9efa4563299b32903dedf48db6ab502e7278 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 13:56:12 -0400 Subject: [PATCH 02/52] Update RMG example c3h4 with multiplicity label --- examples/rmg/c3h4/input.py | 7 +++++++ 1 file changed, 7 insertions(+) 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"), ) From 9f67f0874914ba0a64640031978812282da29a20 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 13:56:43 -0400 Subject: [PATCH 03/52] Update RMG example ch3no2 with multiplicity label --- examples/rmg/ch3no2/input.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/examples/rmg/ch3no2/input.py b/examples/rmg/ch3no2/input.py index 749cd7d7e5..e9219a096d 100644 --- a/examples/rmg/ch3no2/input.py +++ b/examples/rmg/ch3no2/input.py @@ -23,6 +23,7 @@ # List of species species( label='CH3NO2', + multiplicity = 1, reactive=True, structure=adjacencyList( """ @@ -38,6 +39,7 @@ species( label='O2', + multiplicity = 3, reactive=True, structure=adjacencyList( """ @@ -48,11 +50,12 @@ species( label='N2', + multiplicity = 1, reactive=True, structure=adjacencyList( """ - 1 N 1 1 {2,T} - 2 N 1 1 {1,T} + 1 N 0 1 {2,T} + 2 N 0 1 {1,T} """), ) From 7406845ae686531fcd9117d246172a20377cb794 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 13:59:07 -0400 Subject: [PATCH 04/52] Update RMG example diesel with multiplicity label --- examples/rmg/diesel/input.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) 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 From 141803dff118ab57b819c8b12135b732475ebf36 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:03:35 -0400 Subject: [PATCH 05/52] Update RMG example e85 with multiplicity label --- examples/rmg/e85/input.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) 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]"), ) From a536e371f84d5bd7924e4f1b81ff60d65c3d6a94 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:03:52 -0400 Subject: [PATCH 06/52] Update RMG example liquid_phase with multiplicity label --- examples/rmg/liquid_phase/input.py | 7 +++++++ 1 file changed, 7 insertions(+) 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]"), ) From 8a5f90cce3030264e0eb9d020148aac9dba5712c Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:04:06 -0400 Subject: [PATCH 07/52] Update RMG example methylformate with multiplicity label --- examples/rmg/methylformate/input.py | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/examples/rmg/methylformate/input.py b/examples/rmg/methylformate/input.py index ccd3989236..47833dce9d 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 3 0 {2,S} + 2 H 0 0 {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, ) From 44b86f63e3c6c599f9438c7037fc757e1df53b9e Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:04:18 -0400 Subject: [PATCH 08/52] Update RMG example minimal with multiplicity label --- examples/rmg/minimal/input.py | 1 + 1 file changed, 1 insertion(+) 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"), ) From bab0769c461510b97e83ea58ffad4f8dcba500ce Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:04:39 -0400 Subject: [PATCH 09/52] Update RMG example minimal_sensitivity with multiplicity label --- examples/rmg/minimal_sensitivity/input.py | 6 ++++++ 1 file changed, 6 insertions(+) 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"), ) From b7b4725d24b9f19a2600da2a3246c8924143ae8a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:04:49 -0400 Subject: [PATCH 10/52] Update RMG example TEOS with multiplicity label --- examples/rmg/TEOS/input.py | 3 +++ 1 file changed, 3 insertions(+) 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"), ) From 795b41b6b736180dc7b478f130ae2a1300d73547 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:08:11 -0400 Subject: [PATCH 11/52] Add multiplicity label to loadEntry for transport libraries --- rmgpy/data/transport.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/rmgpy/data/transport.py b/rmgpy/data/transport.py index 0b72b2e5fa..4c081dc6c3 100644 --- a/rmgpy/data/transport.py +++ b/rmgpy/data/transport.py @@ -93,6 +93,7 @@ def __init__(self, label='', name='', shortDesc='', longDesc=''): def loadEntry(self, index, label, + multiplicity, molecule, transport, reference=None, @@ -100,10 +101,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, From 22399065b6b6da83c159d3d3451a05a784546f1c Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:09:03 -0400 Subject: [PATCH 12/52] Add multiplicity label to loadEntry for thermo libraries and groups --- rmgpy/data/thermo.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 384ea5dc46..86ba7a953a 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -184,10 +184,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 +218,7 @@ def __init__(self, label='', name='', shortDesc='', longDesc=''): def loadEntry(self, index, label, + multiplicity, molecule, thermo, reference=None, @@ -226,6 +228,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 +236,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 +286,7 @@ def loadEntry(self, label, group, thermo, + multiplicity = [1,2,3,4,5], reference=None, referenceType='', shortDesc='', @@ -290,9 +296,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, From ab77c8ae3231bf4641c722fa795eda2ca170652d Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:11:33 -0400 Subject: [PATCH 13/52] Add multiplicity label to loadEntry for statmech depository --- rmgpy/data/statmech.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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, From f67bab42d5bd632cbea835f782580c13987a90a2 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:13:24 -0400 Subject: [PATCH 14/52] Add multiplicity label to class Entry and loadEntry for forbidden structures --- rmgpy/data/base.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 6092a91eca..e62b60b70b 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(), From 9e6ac85c34ebfdf5fc1244aadccb831aec2c2dc5 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:14:57 -0400 Subject: [PATCH 15/52] Add multiplicity label to loadEntry for kinetics depository --- rmgpy/data/kinetics/depository.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/rmgpy/data/kinetics/depository.py b/rmgpy/data/kinetics/depository.py index 2e03d0418c..a12a9bd7e1 100644 --- a/rmgpy/data/kinetics/depository.py +++ b/rmgpy/data/kinetics/depository.py @@ -137,13 +137,13 @@ def loadEntry(self, rank=None, ): - reactants = [Molecule().fromAdjacencyList(reactant1)] - if reactant2 is not None: reactants.append(Molecule().fromAdjacencyList(reactant2)) - if reactant3 is not None: reactants.append(Molecule().fromAdjacencyList(reactant3)) + reactants = [Molecule().fromAdjacencyList(reactant1, maxMultiplicity=True)] + if reactant2 is not None: reactants.append(Molecule().fromAdjacencyList(reactant2, maxMultiplicity=True)) + if reactant3 is not None: reactants.append(Molecule().fromAdjacencyList(reactant3, maxMultiplicity=True)) - products = [Molecule().fromAdjacencyList(product1)] - if product2 is not None: products.append(Molecule().fromAdjacencyList(product2)) - if product3 is not None: products.append(Molecule().fromAdjacencyList(product3)) + products = [Molecule().fromAdjacencyList(product1, maxMultiplicity=True)] + if product2 is not None: products.append(Molecule().fromAdjacencyList(product2, maxMultiplicity=True)) + if product3 is not None: products.append(Molecule().fromAdjacencyList(product3, maxMultiplicity=True)) reaction = Reaction(reactants=reactants, products=products, degeneracy=degeneracy, duplicate=duplicate, reversible=reversible) From 09fc6521d2ad6a1ce477809dd3363fa0b814ad37 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:24:42 -0400 Subject: [PATCH 16/52] Add multiplicity label to loadEntry for kinetic groups --- rmgpy/data/kinetics/groups.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) 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) From 00db033c792e7a4e3f9288f2ce4814abae9976dc Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:27:10 -0400 Subject: [PATCH 17/52] Add multiplicity label in loadEntry for kinetic libraries --- rmgpy/data/kinetics/library.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index d1c0a5c0ca..e21692f0a8 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].strip()), molecule=[Molecule().fromAdjacencyList(reactant1)])] + if reactant2 is not None: reactants.append(Species(label=reactant2.strip().splitlines()[0].strip(), multiplicity=int(reactant2.strip().splitlines()[1].strip()), molecule=[Molecule().fromAdjacencyList(reactant2)])) + if reactant3 is not None: reactants.append(Species(label=reactant3.strip().splitlines()[0].strip(), multiplicity=int(reactant3.strip().splitlines()[1].strip()), 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].strip()), molecule=[Molecule().fromAdjacencyList(product1)])] + if product2 is not None: products.append(Species(label=product2.strip().splitlines()[0].strip(), multiplicity=int(product2.strip().splitlines()[1].strip()), molecule=[Molecule().fromAdjacencyList(product2)])) + if product3 is not None: products.append(Species(label=product3.strip().splitlines()[0].strip(), multiplicity=int(product3.strip().splitlines()[1].strip()), molecule=[Molecule().fromAdjacencyList(product3)])) comment = "Reaction and kinetics from {0}.".format(self.label) if shortDesc.strip(): From f1572c0a12ee2e508f500c7d1df9fe2fc805a398 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:40:27 -0400 Subject: [PATCH 18/52] Add multiplicity label to group.py - remove multiplicity from GroupAtom - add multiplicity to Group - update __gainRadical() and __loseRadical() - update equivalent() and isSpecificCaseOf() --- rmgpy/molecule/group.pxd | 3 ++- rmgpy/molecule/group.py | 53 ++++++++++++++-------------------------- 2 files changed, 21 insertions(+), 35 deletions(-) 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() From a58ba6997a1869ca71b7d21436a51ea5299a978a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:49:29 -0400 Subject: [PATCH 19/52] Add printMultiplicity=True to toAdjacencyList() in saveSpeciesDictionary() --- rmgpy/chemkin.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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))) From b260eb5a365303898c37c0cdbed40f592ac3cbc1 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:51:01 -0400 Subject: [PATCH 20/52] Add multiplicity label to molecule.py - remove multiplicity from Atom - add multiplicity to Molecule --- rmgpy/molecule/molecule.pxd | 6 +-- rmgpy/molecule/molecule.py | 75 +++++++++++++++++++------------------ 2 files changed, 41 insertions(+), 40 deletions(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 3e474908b8..d5d67d7848 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 @@ -169,7 +169,7 @@ cdef class Molecule(Graph): cpdef fromRDKitMol(self, rdkitmol) - cpdef fromAdjacencyList(self, str adjlist, bint saturateH=?) + cpdef fromAdjacencyList(self, str adjlist, bint saturateH=?, maxMultiplicity=?) cpdef fromXYZ(self, numpy.ndarray atomicNums, numpy.ndarray coordinates) @@ -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..ddbb8b6df4 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=0, 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 @@ -151,7 +149,7 @@ def __reduce__(self): 'sortingLabel': self.sortingLabel, 'atomType': self.atomType.label if self.atomType else None, } - 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): """ @@ -189,7 +187,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 +195,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 +216,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 +226,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 +248,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 +295,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 +306,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 +571,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) @@ -786,6 +769,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 +938,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 +999,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 @@ -1102,6 +1090,16 @@ def fromSMILES(self, smilesstr): # RDKit improperly handles helium and returns it in a triplet state self.fromAdjacencyList('1 He 0 1') return self + elif smilesstr == 'C#O': + # carbon monoxide + self.fromAdjacencyList( + """ + CO + 1 + 1 C 0 1 {2,T} + 2 O 0 1 {1,T} + """) + return self else: rdkitmol = Chem.MolFromSmiles(smilesstr) @@ -1137,7 +1135,7 @@ 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) @@ -1145,18 +1143,20 @@ def fromRDKitMol(self, rdkitmol): number = rdkitatom.GetAtomicNum() element = elements.getElement(number) - # Process spin multiplicity + # Process radicalElectrons 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 + + self.multiplicity = spinMultiplicity # Process charge charge = rdkitatom.GetFormalCharge() - atom = Atom(element, radicalElectrons, spinMultiplicity, charge, '', 0) + atom = Atom(element, radicalElectrons, charge, '', 0) self.vertices.append(atom) # Add bonds by iterating again through atoms @@ -1184,7 +1184,7 @@ def fromRDKitMol(self, rdkitmol): return self newStyleAdjMatcher = re.compile('\s*\d+\s+(\*\d*\s+)?[A-Za-z]+\s+\S+\s+\d').match - def fromAdjacencyList(self, adjlist, saturateH=False): + def fromAdjacencyList(self, adjlist, saturateH=False, maxMultiplicity=False): """ Convert a string adjacency list `adjlist` to a molecular structure. Skips the first line (assuming it's a label) unless `withLabel` is @@ -1194,8 +1194,9 @@ def fromAdjacencyList(self, adjlist, saturateH=False): 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, multiplicity = fromAdjacencyList(adjlist, False, saturateH=saturateH) + if maxMultiplicity: self.multiplicity=multiplicity self.updateConnectivityValues() self.updateAtomTypes() @@ -1445,12 +1446,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=True) return result def isLinear(self): From d72366695c9df976e86e82b161f1270de3f4fa13 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:52:40 -0400 Subject: [PATCH 21/52] Update to- and fromAdjacencyList() for multiplicity in Molecule instead of Atom --- rmgpy/molecule/adjlist.py | 74 ++++++++++++++++++++++----------------- 1 file changed, 41 insertions(+), 33 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 01cc58da1e..08d7f0ee8b 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,12 @@ 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 len(lines[0].split()) == 1: + multiplicity = int(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 @@ -111,7 +118,7 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): index += 1 # Next is the electron state - radicalElectrons = []; spinMultiplicity = [] + radicalElectrons = [] elecState = data[index].upper() if elecState[0] == '{': elecState = elecState[1:-1].split(',') @@ -119,33 +126,19 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): elecState = [elecState] for e in elecState: if e == '0': - radicalElectrons.append(0); spinMultiplicity.append(1) + radicalElectrons.append(0) elif e == '1': - radicalElectrons.append(1); spinMultiplicity.append(2) + radicalElectrons.append(1) 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) + radicalElectrons.append(2) 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) + radicalElectrons.append(3) 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) + radicalElectrons.append(4) elif e == 'X': radicalElectrons.extend([0,1,2,2]) - spinMultiplicity.extend([1,2,1,3]) + else: + raise InvalidAdjacencyListError('Invalid number of radicals.') index += 1 # Next number defines the number of lone electron pairs (if provided) @@ -156,14 +149,16 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): if lpState[0] != '{': if lpState == '0': lonePairElectrons = 0 - if lpState == '1': + elif lpState == '1': lonePairElectrons = 1 - if lpState == '2': + elif lpState == '2': lonePairElectrons = 2 - if lpState == '3': + elif lpState == '3': lonePairElectrons = 3 - if lpState == '4': + elif lpState == '4': lonePairElectrons = 4 + else: + raise InvalidAdjacencyListError('Invalid number of lone electron pairs.') index += 1 else: lonePairElectrons = -1 @@ -172,9 +167,9 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): # 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, radicalElectrons, [0 for e in radicalElectrons], label, [lonePairElectrons]) else: - atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, label, lonePairElectrons) + atom = Atom(atomType[0], radicalElectrons[0], 0, label, lonePairElectrons) # Add the atom to the list atoms.append(atom) @@ -290,7 +285,18 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): print adjlist raise - return atoms + if not group: + if multiplicity == -1: + nRad = 0 + for atom in atoms: + nRad += atom.radicalElectrons + multiplicity = nRad + 1 + + return atoms, multiplicity + + else: + + return atoms ################################################################################ @@ -323,7 +329,7 @@ def getElectronState(radicalElectrons, spinMultiplicity): 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 +343,8 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai pass if label: adjlist += label + '\n' + + if printMultiplicity: adjlist += str(multiplicity) + '\n' # Determine the numbers to use for each atom atomNumbers = {} @@ -360,15 +368,15 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False, removeLonePai atomTypes[atom] = '{{{0}}}'.format(','.join([a.label for a in atom.atomType])) # Electron state(s) if len(atom.radicalElectrons) == 1: - atomElectronStates[atom] = getElectronState(atom.radicalElectrons[0], atom.spinMultiplicity[0]) + 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])) 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)) + atomElectronStates[atom] = '{0}'.format(str(atom.radicalElectrons)) if not removeLonePairs: # Lone Pair(s) atomLonePairs[atom] = atom.lonePairs From 13a539365cc00c3f940302ec9eb8ea28136605ca Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:56:10 -0400 Subject: [PATCH 22/52] Use new molecule.multiplicity attribute in qm.molecule - instead of triplet assumption --- rmgpy/qm/molecule.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From b84567d57e57f108b407e167e21989ab14aae414 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 14:56:58 -0400 Subject: [PATCH 23/52] Add multiplicity to read input file --- rmgpy/rmg/input.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From ed1b8c704db22230388fb697ff9f416827f8e6f1 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 15:06:47 -0400 Subject: [PATCH 24/52] Add copy() to reaction.py to make deep copies of reaction objects --- rmgpy/reaction.pxd | 2 ++ rmgpy/reaction.py | 26 ++++++++++++++++++++++++++ 2 files changed, 28 insertions(+) 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 ################################################################################ From 4b7d957f05e62ff276fd60dd9145a92ef7ad9dca Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 15:12:05 -0400 Subject: [PATCH 25/52] Add multiplicity attribute to Species.py - add condition that given multiplicity must be possible for given number of unpaired electrons (radicals), throws SpeciesError if not - include multiplicity in isIsomorphic() --- rmgpy/species.pxd | 1 + rmgpy/species.py | 24 +++++++++++++++++++----- 2 files changed, 20 insertions(+), 5 deletions(-) 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..49e837f99e 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,15 @@ 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[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 +111,12 @@ def __init__(self, index=-1, label='', thermo=None, conformer=None, self.polarizability = polarizability self.Zrot = Zrot self.energyTransferModel = energyTransferModel + + # Check if multiplicity is possible + 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): """ @@ -189,10 +202,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 From 871c7124c2a1452240f6f89e8158fefb9682b5d3 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 15:16:05 -0400 Subject: [PATCH 26/52] Add species multiplicity model.py - add multiplicity attribute to class Species - add multiplicity to QM thermo database - update makeNewSpecies() --- rmgpy/rmg/model.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) 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]) From 2f4d65174c71a5c7366a612af1bc6d800b0894ad Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 14 Mar 2014 15:20:04 -0400 Subject: [PATCH 27/52] Update family.py for new species multiplicity and search for spin allowed reactions -__generateProductStructures uses angular moment addition theorem to find spin allowed reactions, further improvement and validation might be required --- rmgpy/data/kinetics/family.py | 134 +++++++++++++++++++++------------- 1 file changed, 82 insertions(+), 52 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 300f783e7c..73e204e20d 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -618,13 +618,17 @@ 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() +<<<<<<< HEAD self.forbidden.loadEntry(label=label, group=group, shortDesc=shortDesc, longDesc=longDesc) +======= + self.forbidden.loadEntry(label=label, group=group, multiplicity = multiplicity, shortDesc=shortDesc, longDesc=longDesc, history=history) +>>>>>>> Update family.py for new species multiplicity and search for spin def saveEntry(self, f, entry): """ @@ -834,7 +838,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 +884,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,6 +1138,7 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp # Generate other possible electronic states +<<<<<<< HEAD electronicStructuresList1 = [] electronicStructuresList2 = [] @@ -1142,12 +1147,17 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp struct1a.updateAtomTypes() electronicStructuresList1.append(struct1a) atoms1 = struct1.getRadicalAtoms() +======= + productStructuresList = [] + totalSpin = [] # total spin times 2 +>>>>>>> Update family.py for new species multiplicity and search for spin - 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] +<<<<<<< HEAD if atom1.label != '' and radical1 > 1 and radical1 < 4: if radical1 == 2 and spin1 == 3: @@ -1172,34 +1182,36 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp break else: electronicStructuresList1.append(struct1a) +======= + elif len(reactantStructures) == 2: +>>>>>>> Update family.py for new species multiplicity and search for spin + + spin1 = (reactantStructures[0].multiplicity-1.0)/2.0 + spin2 = (reactantStructures[1].multiplicity-1.0)/2.0 + + count = 0.0 - elif radical1 == 4: + 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) +<<<<<<< HEAD for electronicStructures in electronicStructuresList1: if electronicStructures.isIsomorphic(struct1a): break @@ -1221,10 +1233,19 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp atoms2 = struct2.getRadicalAtoms() for atom2 in atoms2: +======= + if not self.isMoleculeForbidden(struct): + productStructuresList.append([struct]) + + count += 1.0 + + elif len(productStructures) == 2: +>>>>>>> Update family.py for new species multiplicity and search for spin - radical2 = atom2.radicalElectrons - spin2 = atom2.spinMultiplicity + maxSpin1 = productStructures[0].getRadicalCount()/2.0 + maxSpin2 = productStructures[1].getRadicalCount()/2.0 +<<<<<<< HEAD if atom2.label != '' and radical2 > 1 and radical2 < 4: if radical2 == 2 and spin2 == 3: @@ -1249,30 +1270,23 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp break else: electronicStructuresList2.append(struct2a) +======= + count1 = 0.0 +>>>>>>> Update family.py for new species multiplicity and search for spin - 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: +<<<<<<< HEAD atom2.setSpinMultiplicity(3) struct2b = struct2.copy(True) struct2b.updateAtomTypes() @@ -1300,6 +1314,22 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp for structa in electronicStructuresList1: if not (self.isMoleculeForbidden(structa)): productStructuresList.append([structa]) +======= + 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 + + count1 += 1.0 +>>>>>>> Update family.py for new species multiplicity and search for spin return productStructuresList @@ -1407,9 +1437,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] From 712fecbe44c125f1cc778a4071567512c377fb10 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 25 Apr 2014 16:49:07 -0400 Subject: [PATCH 28/52] Change toAdjacencyList() to new adjacency list style --- rmgpy/molecule/adjlist.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 08d7f0ee8b..45451b8dd5 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -359,10 +359,11 @@ def toAdjacencyList(atoms, multiplicity=0, label=None, group=False, removeH=Fals 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])) @@ -376,10 +377,11 @@ def toAdjacencyList(atoms, multiplicity=0, label=None, group=False, removeH=Fals # Atom type atomTypes[atom] = '{0}'.format(atom.element.symbol) # Electron state(s) - atomElectronStates[atom] = '{0}'.format(str(atom.radicalElectrons)) - if not removeLonePairs: - # Lone Pair(s) - atomLonePairs[atom] = atom.lonePairs + atomElectronStates[atom] = '{0}'.format(getElectronState(atom.radicalElectrons, atom.spinMultiplicity)) + # Lone Pair(s) + atomLonePairs[atom] = atom.lonePairs + # Lone Pair(s) + atomCharge[atom] = atom.charge # Determine field widths atomNumberWidth = max([len(s) for s in atomNumbers.values()]) + 1 @@ -387,7 +389,8 @@ def toAdjacencyList(atoms, multiplicity=0, label=None, group=False, removeH=Fals 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: @@ -399,11 +402,12 @@ def toAdjacencyList(atoms, multiplicity=0, label=None, group=False, removeH=Fals 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) + # Radical(s) + adjlist += 'R{0:<{1:d}}'.format(atomElectronStates[atom], atomElectronStateWidth) + # Lone Pair(s) + adjlist += ' L{0:>{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth) + # Charge(s) + adjlist += ' C{0:>{1:d}}'.format(atomCharge[atom], atomChargeWidth) # Bonds list atoms2 = atom.bonds.keys() From d4e0ce4db3b864dafebb0108a88d2192f2df65bd Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 16:52:51 -0400 Subject: [PATCH 29/52] Add multiplicity label to saveEntry() for transport libraries. --- rmgpy/data/transport.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/transport.py b/rmgpy/data/transport.py index 4c081dc6c3..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') From c36551bb665b4da784f2cfae6b30d7403b9d9d5d Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 16:56:08 -0400 Subject: [PATCH 30/52] Add multiplicity label to saveEntry() for thermo libraries. And set printMultiplicity in toAdjacencyList() to False to not print multiplicity as part of the adjacency list --- rmgpy/data/thermo.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index 86ba7a953a..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') From a73f602363ad4deaaf5c976b5ef5e06ca2c116bc Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 16:57:28 -0400 Subject: [PATCH 31/52] Add multiplicity label to saveEntry() for forbiddenstructures.py --- rmgpy/data/base.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index e62b60b70b..2031c6bd61 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -1289,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') From 2faa10c28e9cbc096480c94c397afb566bba509c Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 17:01:29 -0400 Subject: [PATCH 32/52] Replace strip() with split() to read adjacency lists in kinetics library Since the introduction of the multiplicity inside adjacency lists for the kinetics libraries we have to split the multiplicity line (multiplicity n) into its two components where the second component ([1]) is the actual multiplicity number. --- rmgpy/data/kinetics/library.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index e21692f0a8..068fbec2ad 100644 --- a/rmgpy/data/kinetics/library.py +++ b/rmgpy/data/kinetics/library.py @@ -275,13 +275,13 @@ def loadEntry(self, ): - reactants = [Species(label=reactant1.strip().splitlines()[0].strip(), multiplicity=int(reactant1.strip().splitlines()[1].strip()), molecule=[Molecule().fromAdjacencyList(reactant1)])] - if reactant2 is not None: reactants.append(Species(label=reactant2.strip().splitlines()[0].strip(), multiplicity=int(reactant2.strip().splitlines()[1].strip()), molecule=[Molecule().fromAdjacencyList(reactant2)])) - if reactant3 is not None: reactants.append(Species(label=reactant3.strip().splitlines()[0].strip(), multiplicity=int(reactant3.strip().splitlines()[1].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(), multiplicity=int(product1.strip().splitlines()[1].strip()), molecule=[Molecule().fromAdjacencyList(product1)])] - if product2 is not None: products.append(Species(label=product2.strip().splitlines()[0].strip(), multiplicity=int(product2.strip().splitlines()[1].strip()), molecule=[Molecule().fromAdjacencyList(product2)])) - if product3 is not None: products.append(Species(label=product3.strip().splitlines()[0].strip(), multiplicity=int(product3.strip().splitlines()[1].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(): From 8aedab8af7fc50ec8ad7f40243bb0a775e6cdb02 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 17:04:03 -0400 Subject: [PATCH 33/52] Set printMultiplicity label in toAdjacencyList() explicitely to True This makes it easier to understand the difference to other saveEntry() functions --- rmgpy/data/kinetics/common.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/common.py b/rmgpy/data/kinetics/common.py index f6efba6062..09962118c0 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)) From d714d753c488370201b1dd9141a2087ff7ea9119 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 17:05:54 -0400 Subject: [PATCH 34/52] Define type of maxMultiplicity in fromAdjacencyList() --- rmgpy/molecule/molecule.pxd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index d5d67d7848..7c5d934a18 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -169,7 +169,7 @@ cdef class Molecule(Graph): cpdef fromRDKitMol(self, rdkitmol) - cpdef fromAdjacencyList(self, str adjlist, bint saturateH=?, maxMultiplicity=?) + cpdef fromAdjacencyList(self, str adjlist, bint saturateH=?, bint maxMultiplicity=?) cpdef fromXYZ(self, numpy.ndarray atomicNums, numpy.ndarray coordinates) From 154eb1c44a3280fb84e0a98ebc01568e75a376e5 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 17:18:50 -0400 Subject: [PATCH 35/52] Prepare molecule.py for new adjacency list style todjacencyList(): -pass printMultiplicity argument actually to toAdjacencyList() from adjlist.py fromAdjacencyList(): - remove newStyleAdjacencyListMatcher(), all adjacency list are now new style - remove maxMultiplicity argument, multiplicity has now to be explicitely defined --- rmgpy/molecule/molecule.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index ddbb8b6df4..97848c400c 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -108,7 +108,7 @@ class Atom(Vertex): e.g. ``atom.symbol`` instead of ``atom.element.symbol``. """ - def __init__(self, element=None, radicalElectrons=0, charge=0, label='', lonePairs=0, coords=numpy.array([])): + def __init__(self, element=None, radicalElectrons=0, charge=0, label='', lonePairs=None, coords=numpy.array([])): Vertex.__init__(self) if isinstance(element, str): self.element = elements.__dict__[element] @@ -1088,16 +1088,21 @@ 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 - 1 - 1 C 0 1 {2,T} - 2 O 0 1 {1,T} + multiplicity 1 + 1 C U0 L1 {2,T} + 2 O U0 L1 {1,T} """) return self @@ -1183,20 +1188,15 @@ def fromRDKitMol(self, rdkitmol): return self - newStyleAdjMatcher = re.compile('\s*\d+\s+(\*\d*\s+)?[A-Za-z]+\s+\S+\s+\d').match - def fromAdjacencyList(self, adjlist, saturateH=False, maxMultiplicity=False): + def fromAdjacencyList(self, adjlist, saturateH=False): """ Convert a string adjacency list `adjlist` to a molecular structure. Skips the first line (assuming it's a label) unless `withLabel` is ``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, multiplicity = fromAdjacencyList(adjlist, False, saturateH=saturateH) - if maxMultiplicity: self.multiplicity=multiplicity + self.vertices, self.multiplicity = fromAdjacencyList(adjlist, group=False, saturateH=saturateH) self.updateConnectivityValues() self.updateAtomTypes() @@ -1451,7 +1451,7 @@ def toAdjacencyList(self, label='', removeH=False, removeLonePairs=False, printM Convert the molecular structure to a string adjacency list. """ from .adjlist import toAdjacencyList - result = toAdjacencyList(self.vertices, self.multiplicity, label=label, group=False, removeH=removeH, removeLonePairs=removeLonePairs, printMultiplicity=True) + result = toAdjacencyList(self.vertices, self.multiplicity, label=label, group=False, removeH=removeH, removeLonePairs=removeLonePairs, printMultiplicity=printMultiplicity) return result def isLinear(self): From af961ac6aa0ab5933874accf83dad67d06a7cda9 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 17:20:05 -0400 Subject: [PATCH 36/52] Remove maxMultiplicity argument - multiplicity has now to be explicitely defined --- rmgpy/data/kinetics/depository.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/rmgpy/data/kinetics/depository.py b/rmgpy/data/kinetics/depository.py index a12a9bd7e1..2e03d0418c 100644 --- a/rmgpy/data/kinetics/depository.py +++ b/rmgpy/data/kinetics/depository.py @@ -137,13 +137,13 @@ def loadEntry(self, rank=None, ): - reactants = [Molecule().fromAdjacencyList(reactant1, maxMultiplicity=True)] - if reactant2 is not None: reactants.append(Molecule().fromAdjacencyList(reactant2, maxMultiplicity=True)) - if reactant3 is not None: reactants.append(Molecule().fromAdjacencyList(reactant3, maxMultiplicity=True)) + reactants = [Molecule().fromAdjacencyList(reactant1)] + if reactant2 is not None: reactants.append(Molecule().fromAdjacencyList(reactant2)) + if reactant3 is not None: reactants.append(Molecule().fromAdjacencyList(reactant3)) - products = [Molecule().fromAdjacencyList(product1, maxMultiplicity=True)] - if product2 is not None: products.append(Molecule().fromAdjacencyList(product2, maxMultiplicity=True)) - if product3 is not None: products.append(Molecule().fromAdjacencyList(product3, maxMultiplicity=True)) + products = [Molecule().fromAdjacencyList(product1)] + if product2 is not None: products.append(Molecule().fromAdjacencyList(product2)) + if product3 is not None: products.append(Molecule().fromAdjacencyList(product3)) reaction = Reaction(reactants=reactants, products=products, degeneracy=degeneracy, duplicate=duplicate, reversible=reversible) From 49bd3dd225b871262ce21738a1f4abc4dc0c6a5a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 17:22:28 -0400 Subject: [PATCH 37/52] Modify adjlist.py to read and write the new-style adjacency list. --- rmgpy/molecule/adjlist.py | 254 ++++++++++++++++++++++++-------------- 1 file changed, 158 insertions(+), 96 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 45451b8dd5..c91b69c095 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -72,8 +72,9 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): raise InvalidAdjacencyListError('No atoms specified in adjacency list.') # Skip the second line if it contains a multiplicity - if len(lines[0].split()) == 1: - multiplicity = int(lines.pop(0)) + 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.') @@ -117,59 +118,119 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): atomType = [atomType] index += 1 - # Next is the electron state - radicalElectrons = [] - elecState = data[index].upper() - if elecState[0] == '{': - elecState = elecState[1:-1].split(',') - else: - elecState = [elecState] - for e in elecState: - if e == '0': - radicalElectrons.append(0) - elif e == '1': - radicalElectrons.append(1) - elif e == '2': - radicalElectrons.append(2) - elif e == '3': - radicalElectrons.append(3) - elif e == '4': - radicalElectrons.append(4) - elif e == 'X': - radicalElectrons.extend([0,1,2,2]) + # Next the number of unpaired electrons + unpairedElectrons = [] + uState = data[index].upper() + if uState[0] == 'U': + if uState[1] == '{': + uState = uState[2:-1].split(',') else: - raise InvalidAdjacencyListError('Invalid number of radicals.') - index += 1 + 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: + 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 - elif lpState == '1': - lonePairElectrons = 1 - elif lpState == '2': - lonePairElectrons = 2 - elif lpState == '3': - lonePairElectrons = 3 - elif lpState == '4': - lonePairElectrons = 4 + lpState = data[index].upper() + if lpState[0] == 'L': + if lpState[1] == '{': + lpState = lpState[2:-1].split(',') else: - raise InvalidAdjacencyListError('Invalid number of lone electron pairs.') + 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 + lonePairs.append(0) else: - lonePairElectrons = -1 + 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, [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], 0, label, lonePairElectrons) + atom = Atom(atomType[0], unpairedElectrons[0], partialCharges[0], label, lonePairs[0]) # Add the atom to the list atoms.append(atom) @@ -228,21 +289,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 @@ -250,7 +315,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: @@ -286,48 +351,31 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): raise if not group: + nRad = 0 + for atom in atoms: + nRad += atom.radicalElectrons + if multiplicity == -1: - nRad = 0 - for atom in atoms: - nRad += atom.radicalElectrons 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, multiplicity=0, label=None, group=False, removeH=False, removeLonePairs=False, printMultiplicity=False): """ @@ -344,7 +392,7 @@ def toAdjacencyList(atoms, multiplicity=0, label=None, group=False, removeH=Fals if label: adjlist += label + '\n' - if printMultiplicity: adjlist += str(multiplicity) + '\n' + if printMultiplicity: adjlist += 'multiplicity ' + str(multiplicity) + '\n' # Determine the numbers to use for each atom atomNumbers = {} @@ -367,20 +415,27 @@ def toAdjacencyList(atoms, multiplicity=0, label=None, group=False, removeH=Fals 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] = str(atom.radicalElectrons[0]) else: 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)) - # Lone Pair(s) + # Unpaired Electron(s) + atomElectronStates[atom] = '{0}'.format(atom.radicalElectrons) + # Lone Electron Pair(s) atomLonePairs[atom] = atom.lonePairs - # Lone Pair(s) + # Partial Charge(s) atomCharge[atom] = atom.charge # Determine field widths @@ -402,12 +457,19 @@ def toAdjacencyList(atoms, multiplicity=0, label=None, group=False, removeH=Fals adjlist += '{0:<{1:d}}'.format(atomLabels[atom], atomLabelWidth) # Atom type(s) adjlist += '{0:<{1:d}}'.format(atomTypes[atom], atomTypeWidth) - # Radical(s) - adjlist += 'R{0:<{1:d}}'.format(atomElectronStates[atom], atomElectronStateWidth) - # Lone Pair(s) - adjlist += ' L{0:>{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth) - # Charge(s) - adjlist += ' C{0:>{1:d}}'.format(atomCharge[atom], atomChargeWidth) + # 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() From 8db74f9e1a8299dcf509845af378b484e6873da4 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 16 May 2014 17:46:27 -0400 Subject: [PATCH 38/52] Updating examples with new adjacency list style --- examples/cantherm/networks/acetyl+O2/input.py | 9 ++++++++ examples/cantherm/networks/n-butanol/input.py | 4 ++++ examples/generateReactions/input.py | 3 +++ examples/measure/acetylO2/input.py | 9 ++++++++ examples/rmg/1,3-hexadiene/input.py | 4 ++-- examples/rmg/ch3no2/input.py | 22 +++++++++---------- examples/rmg/methylformate/input.py | 4 ++-- examples/thermoEstimator/input.py | 5 +++++ 8 files changed, 45 insertions(+), 15 deletions(-) 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 e73b25ba81..397faee05d 100644 --- a/examples/rmg/1,3-hexadiene/input.py +++ b/examples/rmg/1,3-hexadiene/input.py @@ -32,8 +32,8 @@ reactive=True, structure=adjacencyList( """ - 1 H 0 0 {2,S} - 2 H 0 0 {1,S} + 1 H U0 L0 {2,S} + 2 H U0 L0 {1,S} """), ) species( diff --git a/examples/rmg/ch3no2/input.py b/examples/rmg/ch3no2/input.py index e9219a096d..d7c3f1edfd 100644 --- a/examples/rmg/ch3no2/input.py +++ b/examples/rmg/ch3no2/input.py @@ -27,13 +27,13 @@ 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} """), ) @@ -43,8 +43,8 @@ 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} """), ) @@ -54,8 +54,8 @@ reactive=True, structure=adjacencyList( """ - 1 N 0 1 {2,T} - 2 N 0 1 {1,T} + 1 N U0 L1 {2,T} + 2 N U0 L1 {1,T} """), ) diff --git a/examples/rmg/methylformate/input.py b/examples/rmg/methylformate/input.py index 47833dce9d..77ee25e61e 100644 --- a/examples/rmg/methylformate/input.py +++ b/examples/rmg/methylformate/input.py @@ -33,8 +33,8 @@ reactive=True, structure=adjacencyList( """ - 1 C 3 0 {2,S} - 2 H 0 0 {1,S} + 1 C U3 L0 {2,S} + 2 H U0 L0 {1,S} """), ) species( 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'), ) From 5b2145c7f673607cd4fef5850fddf4cde2fa5344 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 11:03:26 -0400 Subject: [PATCH 39/52] Updating existing unit tests for new adjacency list style --- rmgpy/data/solvationTest.py | 10 +- rmgpy/data/thermoTest.py | 72 +++++++------- rmgpy/molecule/groupTest.py | 59 +++++------ rmgpy/molecule/moleculeTest.py | 176 +++++++++++++++------------------ rmgpy/pdep/networkTest.py | 4 + 5 files changed, 155 insertions(+), 166 deletions(-) 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/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/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/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index 48f5b98cd0..6ced8eb185 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) @@ -469,12 +457,12 @@ 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} +1 *2 C U1 L0 E0 {2,D} {3,S} +2 *1 O U0 L2 E0 {1,D} +3 C U0 L0 E0 {1,S} {4,S} {5,S} {6,S} +4 H U0 L0 E0 {3,S} +5 H U0 L0 E0 {3,S} +6 H U0 L0 E0 {3,S} """ self.molecule = Molecule().fromAdjacencyList(self.adjlist) @@ -557,19 +545,16 @@ def testFromAdjacencyList(self): self.assertTrue(atom1.label == '*2') 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.radicalElectrons == 0) - self.assertTrue(atom2.spinMultiplicity == 1) self.assertTrue(atom2.charge == 0) self.assertTrue(atom3.label == '') self.assertTrue(atom3.element.symbol == 'C') self.assertTrue(atom3.radicalElectrons == 0) - self.assertTrue(atom3.spinMultiplicity == 1) self.assertTrue(atom3.charge == 0) self.assertTrue(bond12.isDouble()) @@ -579,14 +564,18 @@ def testToAdjacencyList(self): """ Test the Molecule.toAdjacencyList() method. """ - adjlist = self.molecule.toAdjacencyList(removeH=False) + adjlist = self.molecule.toAdjacencyList(removeH=False,printMultiplicity=False) self.assertEqual(adjlist.strip(), self.adjlist.strip()) def testFromOldAdjacencyList(self): """ Test we can read things with implicit hydrogens. """ - adjList = "1 O 0" # should be Water + adjList = """ + 1 O U0 L2 {2,S} {3,S} + 2 H U0 L0 {1,S} + 3 H U0 L0 {1,S} + """ # should be Water molecule = Molecule().fromAdjacencyList(adjList, saturateH=True) # only works with saturateH=True self.assertEqual(molecule.getFormula(),'H2O') @@ -605,8 +594,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 +610,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 +654,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 +691,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 +794,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 +803,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 +818,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','CN'] for s in test_strings: molecule = Molecule(SMILES=s) self.assertEqual(s,molecule.toSMILES()) @@ -848,8 +836,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 +847,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, From 16a0cea19c3dbd263c85ab701735f746bcb1d34b Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 11:07:41 -0400 Subject: [PATCH 40/52] Remove maxMultiplicity argument - multiplicity has to be defined for each entry which makes maxMultiplicity unnecessary --- rmgpy/molecule/molecule.pxd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 7c5d934a18..e18e5482cc 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -169,7 +169,7 @@ cdef class Molecule(Graph): cpdef fromRDKitMol(self, rdkitmol) - cpdef fromAdjacencyList(self, str adjlist, bint saturateH=?, bint maxMultiplicity=?) + cpdef fromAdjacencyList(self, str adjlist, bint saturateH=?) cpdef fromXYZ(self, numpy.ndarray atomicNums, numpy.ndarray coordinates) From 31322a51e685dab13a894d43ce24334fdeedaaad Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 11:09:20 -0400 Subject: [PATCH 41/52] Updating existing unit test for new adjacency list style --- rmgpy/molecule/atomtypeTest.py | 108 ++++++++++++++++----------------- 1 file changed, 54 insertions(+), 54 deletions(-) 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): From ca65b4d872c3cc35f66e8285a8c30d91faaf6e72 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 11:32:35 -0400 Subject: [PATCH 42/52] Update fromRDKitMol() for the new adjacency list style - this also updates all the functions based on fromEDKitMol() like fromSMILES() - all these molecule input functions are based on the assumption of maximum multiplicity is the most stable multiplicity (multiplicity=radcial count + 1) --- rmgpy/molecule/molecule.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 97848c400c..98b49a03fe 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -108,7 +108,7 @@ class Atom(Vertex): e.g. ``atom.symbol`` instead of ``atom.element.symbol``. """ - def __init__(self, element=None, radicalElectrons=0, charge=0, label='', lonePairs=None, 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] @@ -1147,19 +1147,10 @@ def fromRDKitMol(self, rdkitmol): # Use atomic number as key for element number = rdkitatom.GetAtomicNum() element = elements.getElement(number) - - # Process radicalElectrons - 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 - - self.multiplicity = spinMultiplicity # Process charge charge = rdkitatom.GetFormalCharge() + radicalElectrons = rdkitatom.GetNumRadicalElectrons() atom = Atom(element, radicalElectrons, charge, '', 0) self.vertices.append(atom) @@ -1186,6 +1177,11 @@ 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 def fromAdjacencyList(self, adjlist, saturateH=False): From d4dd5006e6d136b019e0086ed8c3aaa37c905589 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 11:40:47 -0400 Subject: [PATCH 43/52] Change treatment of missing lone pair card for molecules a missing lone pair card (L) leads to the assumption that there are no lone pairs, therefore L0 for groups a missing L leads to the assumption that it was not defined therefore None --- rmgpy/molecule/adjlist.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index c91b69c095..740b843cea 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -177,9 +177,15 @@ def fromAdjacencyList(adjlist, group=False, saturateH=False): raise InvalidAdjacencyListError('Number of lone electron pairs not recognized.') index += 1 else: - lonePairs.append(0) + if group: + lonePairs.append(None) + else: + lonePairs.append(0) else: - lonePairs.append(0) + if group: + lonePairs.append(None) + else: + lonePairs.append(0) # Next the number of partial charges (if provided) partialCharges = [] From 5be8ba113f43c40adea45df25a84f9db4ec38498 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 11:41:38 -0400 Subject: [PATCH 44/52] Update existing unit test for new adjacency list style --- rmgpy/reactionTest.py | 3 +++ rmgpy/speciesTest.py | 2 ++ rmgpy/transportDataTest.py | 29 +++++++++++++++-------------- 3 files changed, 20 insertions(+), 14 deletions(-) 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/speciesTest.py b/rmgpy/speciesTest.py index f20020d6e6..f7f03aee0c 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)'), @@ -86,6 +87,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) From a37b47d6e6fa0f6e7d6a04cab761b34a6862eaa6 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 11:43:19 -0400 Subject: [PATCH 45/52] Prevent species generation from crashing if molecule is None and add multiplicity to class representation string --- rmgpy/species.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/rmgpy/species.py b/rmgpy/species.py index 49e837f99e..bce224e763 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -97,8 +97,9 @@ def __init__(self, index=-1, label='', multiplicity=102, thermo=None, conformer= dipoleMoment=None, polarizability=None, Zrot=None, energyTransferModel=None, reactive=True): self.index = index - if label != '' and molecule[0].getRadicalCount()>1: - if '_' not in label: label += '_({0})' .format(_multiplicity_labels[multiplicity]) + if label != '' and molecule is not None: + if molecule[0].getRadicalCount()>1: + if '_' not in label: label += '_({0})' .format(_multiplicity_labels[multiplicity]) self.label = label self.multiplicity = multiplicity self.thermo = thermo @@ -113,10 +114,11 @@ def __init__(self, index=-1, label='', multiplicity=102, thermo=None, conformer= self.energyTransferModel = energyTransferModel # Check if multiplicity is possible - 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)) + if molecule is not None: + 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): """ @@ -126,6 +128,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]) From 472b828fc0547e0b10c8f10f2042755aef1e1d65 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 11:43:54 -0400 Subject: [PATCH 46/52] Update existing unit tests for the new adjacency list style --- rmgpy/solver/simpleTest.py | 5 +++++ 1 file changed, 5 insertions(+) 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)")) ) From 6ac5b9c95e87718d1e7a99b9a6a653ff1c050d9d Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 15:26:21 -0400 Subject: [PATCH 47/52] Remove merge conflicts --- rmgpy/data/kinetics/family.py | 125 +--------------------------------- 1 file changed, 1 insertion(+), 124 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 73e204e20d..05d6311d23 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -624,11 +624,7 @@ def loadForbidden(self, label, group, multiplicity = [1,2,3,4,5], shortDesc='', """ if not self.forbidden: self.forbidden = ForbiddenStructures() -<<<<<<< HEAD - self.forbidden.loadEntry(label=label, group=group, shortDesc=shortDesc, longDesc=longDesc) -======= - self.forbidden.loadEntry(label=label, group=group, multiplicity = multiplicity, shortDesc=shortDesc, longDesc=longDesc, history=history) ->>>>>>> Update family.py for new species multiplicity and search for spin + self.forbidden.loadEntry(label=label, group=group, multiplicity = multiplicity, shortDesc=shortDesc, longDesc=longDesc) def saveEntry(self, f, entry): """ @@ -1138,53 +1134,15 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp # Generate other possible electronic states -<<<<<<< HEAD - electronicStructuresList1 = [] - electronicStructuresList2 = [] - - struct1 = productStructures[0] - struct1a = struct1.copy(True) - struct1a.updateAtomTypes() - electronicStructuresList1.append(struct1a) - atoms1 = struct1.getRadicalAtoms() -======= productStructuresList = [] totalSpin = [] # total spin times 2 ->>>>>>> Update family.py for new species multiplicity and search for spin # implement Angular Momentum Addition Theorem if len(reactantStructures) == 1: totalSpin = [(reactantStructures[0].multiplicity-1.0)/2.0] -<<<<<<< HEAD - 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: ->>>>>>> Update family.py for new species multiplicity and search for spin spin1 = (reactantStructures[0].multiplicity-1.0)/2.0 spin2 = (reactantStructures[1].multiplicity-1.0)/2.0 @@ -1211,68 +1169,17 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp struct.multiplicity = int((maxSpin1-count)*2.0+1.0) -<<<<<<< HEAD - for electronicStructures in electronicStructuresList1: - if electronicStructures.isIsomorphic(struct1a): - break - else: - electronicStructuresList1.append(struct1a) - - 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: -======= if not self.isMoleculeForbidden(struct): productStructuresList.append([struct]) count += 1.0 elif len(productStructures) == 2: ->>>>>>> Update family.py for new species multiplicity and search for spin maxSpin1 = productStructures[0].getRadicalCount()/2.0 maxSpin2 = productStructures[1].getRadicalCount()/2.0 -<<<<<<< HEAD - 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 ->>>>>>> Update family.py for new species multiplicity and search for spin while (maxSpin1-count1) >= 0.0: @@ -1286,35 +1193,6 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp if (maxSpin1-count1+maxSpin2-count2-count) in totalSpin: -<<<<<<< HEAD - atom2.setSpinMultiplicity(3) - struct2b = struct2.copy(True) - struct2b.updateAtomTypes() - - for electronicStructures in electronicStructuresList2: - if electronicStructures.isIsomorphic(struct2a): - break - else: - electronicStructuresList2.append(struct2a) - - 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]) -======= struct1 = productStructures[0].copy(deep=True) struct2 = productStructures[1].copy(deep=True) @@ -1329,7 +1207,6 @@ def __generateProductStructures(self, reactantStructures, maps, forward, failsSp count2 += 1.0 count1 += 1.0 ->>>>>>> Update family.py for new species multiplicity and search for spin return productStructuresList From c2199b9154637c486873de44123b0a01c7aa9904 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 19 May 2014 15:27:24 -0400 Subject: [PATCH 48/52] Add write multiplicity label to kinetics groups --- rmgpy/data/kinetics/common.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rmgpy/data/kinetics/common.py b/rmgpy/data/kinetics/common.py index 09962118c0..33f14ca671 100644 --- a/rmgpy/data/kinetics/common.py +++ b/rmgpy/data/kinetics/common.py @@ -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()) From cb0e8680cb3422dc69608fcacba6cfa4f219f42a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 20 May 2014 11:01:05 -0400 Subject: [PATCH 49/52] Add multiplicity and lone pairs to functions required for pickling Atoms and molecule classes --- rmgpy/molecule/molecule.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 98b49a03fe..d4d409a622 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -148,6 +148,7 @@ 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.charge, self.label), d) @@ -161,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 @@ -603,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 From b999f1b473c3abf64f0942a56105b33f8f79b494 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 20 May 2014 11:03:23 -0400 Subject: [PATCH 50/52] Add additional unit test for handling the new adjacency list styles --- rmgpy/molecule/moleculeTest.py | 284 +++++++++++++++++++++++++++------ 1 file changed, 232 insertions(+), 52 deletions(-) diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index 6ced8eb185..d79898e2bb 100644 --- a/rmgpy/molecule/moleculeTest.py +++ b/rmgpy/molecule/moleculeTest.py @@ -456,45 +456,73 @@ class TestMolecule(unittest.TestCase): """ def setUp(self): - self.adjlist = """ -1 *2 C U1 L0 E0 {2,D} {3,S} -2 *1 O U0 L2 E0 {1,D} -3 C U0 L0 E0 {1,S} {4,S} {5,S} {6,S} -4 H U0 L0 E0 {3,S} -5 H U0 L0 E0 {3,S} -6 H U0 L0 E0 {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 @@ -503,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()) @@ -516,68 +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.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.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.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,printMultiplicity=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 U0 L2 {2,S} {3,S} - 2 H U0 L0 {1,S} - 3 H U0 L0 {1,S} - """ # 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): """ @@ -818,7 +998,7 @@ def testSMILES(self): Test that we can generate a few SMILES strings as expected """ import rmgpy.molecule - 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','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()) From ce3c8d7d57f13522594578ee927a497c8bfdfe7a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 20 May 2014 11:05:02 -0400 Subject: [PATCH 51/52] Prevent species creation from crashing of molecule is an empty vector and add multiplicity to pickling functions --- rmgpy/species.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/rmgpy/species.py b/rmgpy/species.py index bce224e763..8ae80ac36c 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -98,8 +98,9 @@ def __init__(self, index=-1, label='', multiplicity=102, thermo=None, conformer= energyTransferModel=None, reactive=True): self.index = index if label != '' and molecule is not None: - if molecule[0].getRadicalCount()>1: - if '_' not in label: label += '_({0})' .format(_multiplicity_labels[multiplicity]) + 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 @@ -115,10 +116,11 @@ def __init__(self, index=-1, label='', multiplicity=102, thermo=None, conformer= # Check if multiplicity is possible if molecule is not None: - 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)) + 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): """ @@ -159,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 From 19f75020ad435e82a557d7a180de0006ae871357 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 20 May 2014 11:17:01 -0400 Subject: [PATCH 52/52] Add multiplicity to species unit test for pickling species --- rmgpy/speciesTest.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rmgpy/speciesTest.py b/rmgpy/speciesTest.py index f7f03aee0c..2d9d9038be 100644 --- a/rmgpy/speciesTest.py +++ b/rmgpy/speciesTest.py @@ -62,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))