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/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 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. 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.