From a83ef1be5d3fc6525076934a23e7f4347cc3bb21 Mon Sep 17 00:00:00 2001 From: Belinda Slakman Date: Tue, 1 Nov 2016 17:33:03 -0400 Subject: [PATCH 1/3] Unit test for lone pair/biradical thermo. If we're estimating via group additivity, a species with lone pair or biradical on the same atom should give the same thermo result. Testing on a cyclopropane birad C1C[C]1. Otherwise, we fall up to more general, erroneous groups since the lone pair form is not a radical, and it's not bonded to hydrogens. --- rmgpy/data/thermoTest.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/rmgpy/data/thermoTest.py b/rmgpy/data/thermoTest.py index da3469478b..db2a96d4b0 100644 --- a/rmgpy/data/thermoTest.py +++ b/rmgpy/data/thermoTest.py @@ -237,6 +237,41 @@ def testParseThermoComments(self): self.assertEqual(source['GAV']['ring'][0][1],-1) # the weight of benzene contribution should be -1 self.assertEqual(source['GAV']['group'][0][1],2) # weight of the group(Cs-CsCsHH) conbtribution should be 2 + def testLonePairThermoGeneration(self): + """ + Test that for the biradical and lone pair form of the same molecule, we + are getting the same thermo data by transforming the lone pair into a biradical. + """ + spc1 = Species(molecule=[Molecule().fromAdjacencyList(""" + multiplicity 3 + 1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} + 2 C u0 p0 c0 {1,S} {3,S} {6,S} {7,S} + 3 C u2 p0 c0 {1,S} {2,S} + 4 H u0 p0 c0 {1,S} + 5 H u0 p0 c0 {1,S} + 6 H u0 p0 c0 {2,S} + 7 H u0 p0 c0 {2,S} + """ + )]) + + spc2 = Species(molecule=[Molecule().fromAdjacencyList(""" + 1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} + 2 C u0 p0 c0 {1,S} {3,S} {6,S} {7,S} + 3 C u0 p1 c0 {1,S} {2,S} + 4 H u0 p0 c0 {1,S} + 5 H u0 p0 c0 {1,S} + 6 H u0 p0 c0 {2,S} + 7 H u0 p0 c0 {2,S} + """ + )]) + + spc1.generateResonanceIsomers() + spc2.generateResonanceIsomers() + thermo1 = self.database.getThermoDataFromGroups(spc1) + thermo2 = self.database.getThermoDataFromGroups(spc2) + + self.assertEqual(thermo1.getEnthalpy(298), thermo2.getEnthalpy(298)) + class TestThermoDatabaseAromatics(TestThermoDatabase): """ Test only Aromatic species. From 96518457620485a135dba4bc17e6fa6f20744b04 Mon Sep 17 00:00:00 2001 From: Belinda Slakman Date: Tue, 1 Nov 2016 16:57:21 -0400 Subject: [PATCH 2/3] Move transformLonePairs() method from solvation to molecule. Since we're now doing it for both gas phase and solvation thermo estimation, and it was more of a molecule method anyway. --- rmgpy/data/solvation.py | 50 +------------------------------------- rmgpy/molecule/molecule.py | 49 +++++++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+), 49 deletions(-) diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index 690bef04e6..5fe5408780 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -701,55 +701,7 @@ def getSoluteDataFromGroups(self, species): soluteData.comment = "Average of {0}".format(" and ".join(comments)) return soluteData - - def transformLonePairs(self, molecule): - """ - Changes lone pairs in a molecule to two radicals for purposes of finding - solute data via group additivity. Transformed for each atom based on valency. - """ - saturatedStruct = molecule.copy(deep=True) - addedToPairs = {} - for atom in saturatedStruct.atoms: - addedToPairs[atom] = 0 - if atom.lonePairs > 0: - charge = atom.charge # Record this so we can conserve it when checking - bonds = saturatedStruct.getBonds(atom) - sumBondOrders = 0 - for key, bond in bonds.iteritems(): - if bond.order == 'S': sumBondOrders += 1 - if bond.order == 'D': sumBondOrders += 2 - if bond.order == 'T': sumBondOrders += 3 - if bond.order == 'B': sumBondOrders += 1.5 # We should always have 2 'B' bonds (but what about Cbf?) - if atomTypes['Val4'] in atom.atomType.generic: # Carbon, Silicon - while(atom.radicalElectrons + charge + sumBondOrders < 4): - atom.decrementLonePairs() - atom.incrementRadical() - atom.incrementRadical() - addedToPairs[atom] += 1 - if atomTypes['Val5'] in atom.atomType.generic: # Nitrogen - while(atom.radicalElectrons + charge + sumBondOrders < 3): - atom.decrementLonePairs() - atom.incrementRadical() - atom.incrementRadical() - addedToPairs[atom] += 1 - if atomTypes['Val6'] in atom.atomType.generic: # Oxygen, sulfur - while(atom.radicalElectrons + charge + sumBondOrders < 2): - atom.decrementLonePairs() - atom.incrementRadical() - atom.incrementRadical() - addedToPairs[atom] += 1 - if atomTypes['Val7'] in atom.atomType.generic: # Chlorine - while(atom.radicalElectrons + charge + sumBondOrders < 1): - atom.decrementLonePairs() - atom.incrementRadical() - atom.incrementRadical() - addedToPairs[atom] += 1 - - saturatedStruct.update() - saturatedStruct.updateLonePairs() - - return saturatedStruct, addedToPairs def removeHBonding(self, saturatedStruct, addedToRadicals, addedToPairs, soluteData): @@ -809,7 +761,7 @@ def estimateSoluteViaGroupAdditivity(self, molecule): # Change lone pairs to radicals based on valency if sum([atom.lonePairs for atom in saturatedStruct.atoms]) > 0: # molecule contains lone pairs - saturatedStruct, addedToPairs = self.transformLonePairs(saturatedStruct) + saturatedStruct, addedToPairs = saturatedStruct.transformLonePairs() # Now saturate radicals with H if sum([atom.radicalElectrons for atom in saturatedStruct.atoms]) > 0: # radical species diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 24de6d9a82..410e8dbaea 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1563,6 +1563,55 @@ def getNetCharge(self): charge += atom.charge return charge + def transformLonePairs(self): + """ + Changes lone pairs in a molecule to two radicals for purposes of finding + solute data or thermo data via group additivity. Transformed for each atom + based on valency. + """ + saturatedStruct = self.copy(deep=True) + addedToPairs = {} + + for atom in saturatedStruct.atoms: + addedToPairs[atom] = 0 + if atom.lonePairs > 0: + charge = atom.charge # Record this so we can conserve it when checking + bonds = saturatedStruct.getBonds(atom) + sumBondOrders = 0 + for key, bond in bonds.iteritems(): + if bond.order == 'S': sumBondOrders += 1 + if bond.order == 'D': sumBondOrders += 2 + if bond.order == 'T': sumBondOrders += 3 + if bond.order == 'B': sumBondOrders += 1.5 # We should always have 2 'B' bonds (but what about Cbf?) + if atomTypes['Val4'] in atom.atomType.generic: # Carbon, Silicon + while(atom.radicalElectrons + charge + sumBondOrders < 4): + atom.decrementLonePairs() + atom.incrementRadical() + atom.incrementRadical() + addedToPairs[atom] += 1 + if atomTypes['Val5'] in atom.atomType.generic: # Nitrogen + while(atom.radicalElectrons + charge + sumBondOrders < 3): + atom.decrementLonePairs() + atom.incrementRadical() + atom.incrementRadical() + addedToPairs[atom] += 1 + if atomTypes['Val6'] in atom.atomType.generic: # Oxygen, sulfur + while(atom.radicalElectrons + charge + sumBondOrders < 2): + atom.decrementLonePairs() + atom.incrementRadical() + atom.incrementRadical() + addedToPairs[atom] += 1 + if atomTypes['Val7'] in atom.atomType.generic: # Chlorine + while(atom.radicalElectrons + charge + sumBondOrders < 1): + atom.decrementLonePairs() + atom.incrementRadical() + atom.incrementRadical() + addedToPairs[atom] += 1 + + saturatedStruct.update() + saturatedStruct.updateLonePairs() + + return saturatedStruct, addedToPairs def saturate(self): """ Saturate the molecule by replacing all radicals with bonds to hydrogen atoms. Changes self molecule object. From 449315d0f36f237afef82855262b2341af7624d0 Mon Sep 17 00:00:00 2001 From: Belinda Slakman Date: Tue, 1 Nov 2016 17:00:10 -0400 Subject: [PATCH 3/3] Transform lone pairs into biradicals before deciding to use HBI or not. Before, we were doing the normal group additivity thermo on these species, but we don't have Benson groups for them. So we were falling up to general/ erroneous nodes. --- rmgpy/data/thermo.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index ae1a5795b7..2e9023f1c1 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -1363,15 +1363,19 @@ def estimateThermoViaGroupAdditivity(self, molecule): # will probably not visit the right atoms, and so will get the thermo wrong molecule.sortAtoms() - if molecule.isRadical(): # radical species - thermoData = self.estimateRadicalThermoViaHBI(molecule, self.computeGroupAdditivityThermo) + # Convert lone pairs to radicals before determining if molecule is a radical + transformed_mol = molecule.copy(deep=True) + transformed_mol, addedToPairs = transformed_mol.transformLonePairs() + + if transformed_mol.isRadical(): # radical species + thermoData = self.estimateRadicalThermoViaHBI(transformed_mol, self.computeGroupAdditivityThermo) return thermoData else: # non-radical species - thermoData = self.computeGroupAdditivityThermo(molecule) + thermoData = self.computeGroupAdditivityThermo(transformed_mol) # Correct entropy for symmetry number - if not 'saturated' in molecule.props: - thermoData.S298.value_si -= constants.R * math.log(molecule.getSymmetryNumber()) + if not 'saturated' in molecule.props: + thermoData.S298.value_si -= constants.R * math.log(transformed_mol.getSymmetryNumber()) return thermoData