Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix lone pair thermo #808

Closed
wants to merge 3 commits into from
Closed

Fix lone pair thermo #808

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a83ef1b
Unit test for lone pair/biradical thermo.
bslakman Nov 1, 2016
9651845
Move transformLonePairs() method from solvation
bslakman Nov 1, 2016
449315d
Transform lone pairs into biradicals before deciding to use HBI or not.
bslakman Nov 1, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 1 addition & 49 deletions rmgpy/data/solvation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):

Expand Down Expand Up @@ -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
Expand Down
14 changes: 9 additions & 5 deletions rmgpy/data/thermo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down
35 changes: 35 additions & 0 deletions rmgpy/data/thermoTest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
49 changes: 49 additions & 0 deletions rmgpy/molecule/molecule.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down