Skip to content

Commit

Permalink
Merge pull request #364 from bslakman/master
Browse files Browse the repository at this point in the history 
Fix to solvation thermo for lone pairs
  • Loading branch information
@connie
connie committed Apr 1, 2015
2 parents 1df1c3f + c9625a5 commit efb4d1d
Show file tree
Hide file tree
Showing 4 changed files with 241 additions and 90 deletions.
42 changes: 27 additions & 15 deletions documentation/source/users/rmg/liquids.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ dipole interactions related to the polarizability of the solvent) [Vitha2006]_,
`lL` term accounts for the contribution from cavity formation and dispersion (dispersion interactions are
known to scale with solute volume [Vitha2006]_, [Abraham1999]_. The `eE` term, like the `sS` term,
accounts for residual contributions from dipolarity/polarizability related interactions for solutes
whose blend of dipolarity/polarizability differs from that implicitly built into the `S` parameter [Vitha2006]_, [Abraham1990]_, [Jalan2010]_.
whose blend of dipolarity/polarizability differs from that implicitly built into the `S` parameter [Vitha2006]_, [Abraham1999]_, [Jalan2010]_.
The `aA` and `bB` terms account for the contribution of hydrogen bonding between the solute and
the surrounding solvent molecules. H-bonding interactions require two terms as the solute (or solvent)
can act as acceptor (donor) and vice versa. The descriptor `A` is a measure of the solute's ability
Expand Down Expand Up @@ -249,18 +249,30 @@ This is an example of an input file for a liquid-phase system::
saveSimulationProfiles=True,
)

.. [Vitha2006] \ Vitha2006
.. [Abrham1999] \ Abraham1999
.. [Jalan2010] \ Jalan2010
.. [Poole2009] \ Poole2009
.. [Platts1999] \ Platts1999
.. [Mintz2007] \ Mintz2007
.. [Mintz2007a] \ Mintz2007a
.. [Mintz2007b] \ Mintz2007b
.. [Mintz2007c] \ Mintz2007c
.. [Mintz2007d] \ Mintz2007d
.. [Mintz2008] \ Mintz2008
.. [Mintz2008a] \ Mintz2008a
.. [Mintz2009] \ Mintz2009
.. [Rice1985] \ Rice1985
.. [Vitha2006] \ M. Vitha and P.W. Carr. "The chemical interpretation and practice of linear solvation energy relationships in chromatography." *J. Chromatogr. A.* **1126(1-2)**, p. 143-194 (2006).
.. [Abraham1999] \ M.H. Abraham et al. "Correlation and estimation of gas-chloroform and water-chloroformpartition coefficients by a linear free energy relationship method." *J. Pharm. Sci.* **88(7)**, p. 670-679 (1999).
.. [Jalan2010] \ A. Jalan et al. "Predicting solvation energies for kinetic modeling." *Annu. Rep.Prog. Chem., Sect. C* **106**, p. 211-258 (2010).
.. [Poole2009] \ C.F. Poole et al. "Determination of solute descriptors by chromatographic methods." *Anal. Chim. Acta* **652(1-2)** p. 32-53 (2009).
.. [Platts1999] \ J. Platts and D. Butina. "Estimation of molecular linear free energy relation descriptorsusing a group contribution approach." *J. Chem. Inf. Comput. Sci.* **39**, p. 835-845 (1999).
.. [Mintz2007] \ C. Mintz et al. "Enthalpy of solvation correlations for gaseous solutes dissolved inwater and in 1-octanol based on the Abraham model." *J. Chem. Inf. Model.* **47(1)**, p. 115-121 (2007).
.. [Mintz2007a] \ C. Mintz et al. "Enthalpy of solvation corrections for gaseous solutes dissolved in benzene and in alkane solvents based on the Abraham model." *QSAR Comb. Sci.* **26(8)**, p. 881-888 (2007).
.. [Mintz2007b] \ C. Mintz et al. "Enthalpy of solvation correlations for gaseous solutes dissolved in toluene and carbon tetrachloride based on the Abraham model." *J. Sol. Chem.* **36(8)**, p. 947-966 (2007).
.. [Mintz2007c] \ C. Mintz et al. "Enthalpy of solvation correlations for gaseous solutes dissolved indimethyl sulfoxide and propylene carbonate based on the Abraham model." *Thermochim. Acta* **459(1-2)**, p, 17-25 (2007).
.. [Mintz2007d] \ C. Mintz et al. "Enthalpy of solvation correlations for gaseous solutes dissolved inchloroform and 1,2-dichloroethane based on the Abraham model." *Fluid Phase Equilib.* **258(2)**, p. 191-198 (2007).
.. [Mintz2008] \ C. Mintz et al. "Enthalpy of solvation correlations for gaseous solutes dissolved inlinear alkanes (C5-C16) based on the Abraham model." *QSAR Comb. Sci.* **27(2)**, p. 179-186 (2008).
.. [Mintz2008a] \ C. Mintz et al. "Enthalpy of solvation correlations for gaseous solutes dissolved inalcohol solvents based on the Abraham model." *QSAR Comb. Sci.* **27(5)**, p. 627-635 (2008).
.. [Mintz2009] \ C. Mintz et al. "Enthalpy of solvation correlations for organic solutes and gasesdissolved in acetonitrile and acetone." *Thermochim. Acta* **484(1-2)**, p. 65-69 (2009).
.. [Rice1985] \ S.A. Rice. "Diffusion-limited reactions". In *Comprehensive Chemical Kinetics*, EditorsC.H. Bamford, C.F.H. Tipper and R.G. Compton. **25**, (1985).
220 changes: 151 additions & 69 deletions rmgpy/data/solvation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@
from copy import deepcopy
from base import Database, Entry, makeLogicNode, DatabaseError

from rmgpy.molecule import Molecule, Atom, Bond, Group
from rmgpy.molecule import Molecule, Atom, Bond, Group, atomTypes

################################################################################

Expand Down Expand Up @@ -484,13 +484,14 @@ def loadGroups(self, path):
Load the solute database from the given `path` on disk, where `path`
points to the top-level folder of the solute database.
Two sets of groups for additivity, atom-centered ('abraham') and non atom-centered
('nonacentered').
Three sets of groups for additivity, atom-centered ('abraham'), non atom-centered
('nonacentered'), and radical corrections ('radical')
"""
logging.info('Loading Platts additivity group database from {0}...'.format(path))
self.groups = {}
self.groups['abraham'] = SoluteGroups(label='abraham').load(os.path.join(path, 'abraham.py' ), self.local_context, self.global_context)
self.groups['nonacentered'] = SoluteGroups(label='nonacentered').load(os.path.join(path, 'nonacentered.py' ), self.local_context, self.global_context)
self.groups['radical'] = SoluteGroups(label='radical').load(os.path.join(path, 'radical.py' ), self.local_context, self.global_context)

def save(self, path):
"""
Expand Down Expand Up @@ -519,6 +520,7 @@ def saveGroups(self, path):
if not os.path.exists(path): os.mkdir(path)
self.groups['abraham'].save(os.path.join(path, 'abraham.py'))
self.groups['nonacentered'].save(os.path.join(path, 'nonacentered.py'))
self.groups['radical'].save(os.path.join(path, 'radical.py'))

def loadOld(self, path):
"""
Expand Down Expand Up @@ -672,7 +674,117 @@ def getSoluteDataFromGroups(self, species):
soluteData.comment = "Average of {0}".format(" and ".join(comments))

return soluteData

def saturateRadicals(self, molecule):
saturatedStruct = molecule.copy(deep=True)

# Saturate structure by replacing all radicals with bonds to
# hydrogen atoms
added = {}
for atom in saturatedStruct.atoms:
for i in range(atom.radicalElectrons):
H = Atom('H')
bond = Bond(atom, H, 'S')
saturatedStruct.addAtom(H)
saturatedStruct.addBond(bond)
if atom not in added:
added[atom] = []
added[atom].append([H, bond])
atom.decrementRadical()

# Update the atom types of the saturated structure (not sure why
# this is necessary, because saturating with H shouldn't be
# changing atom types, but it doesn't hurt anything and is not
# very expensive, so will do it anyway)
saturatedStruct.updateConnectivityValues()
saturatedStruct.sortVertices()
saturatedStruct.updateAtomTypes()
saturatedStruct.updateLonePairs()
saturatedStruct.multiplicity = 1

return saturatedStruct, added

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.updateConnectivityValues()
saturatedStruct.sortVertices()
saturatedStruct.updateAtomTypes()
saturatedStruct.updateLonePairs()
saturatedStruct.updateMultiplicity()

return saturatedStruct, addedToPairs

def removeHBonding(self, saturatedStruct, addedToRadicals, addedToPairs, soluteData):

# Remove hydrogen bonds and restore the radical
for atom in addedToRadicals:
for H, bond in addedToRadicals[atom]:
saturatedStruct.removeBond(bond)
saturatedStruct.removeAtom(H)
atom.incrementRadical()

# Change transformed lone pairs back
for atom in addedToPairs:
if addedToPairs[atom] > 0:
for pair in range(1, addedToPairs[atom]):
saturatedStruct.decrementRadical()
saturatedStruct.decrementRadical()
saturatedStruct.incrementLonePairs()

# Update Abraham 'A' H-bonding parameter for unsaturated struct
for atom in saturatedStruct.atoms:
# Iterate over heavy (non-hydrogen) atoms
if atom.isNonHydrogen() and atom.radicalElectrons > 0:
for electron in range(1, atom.radicalElectrons):
# Get solute data for radical group
try:
self.__addGroupSoluteData(soluteData, self.groups['radical'], saturatedStruct, {'*':atom})
except KeyError: pass

return soluteData

def estimateSoluteViaGroupAdditivity(self, molecule):
"""
Return the set of Abraham solute parameters corresponding to a given
Expand All @@ -693,73 +805,43 @@ def estimateSoluteViaGroupAdditivity(self, molecule):
L = 0.13,
A = 0.003
)

if sum([atom.radicalElectrons for atom in molecule.atoms]) > 0: # radical species

# Make a copy of the structure so we don't change the original
saturatedStruct = molecule.copy(deep=True)

# Saturate structure by replacing all radicals with bonds to
# hydrogen atoms
added = {}
for atom in saturatedStruct.atoms:
for i in range(atom.radicalElectrons):
H = Atom('H')
bond = Bond(atom, H, 'S')
saturatedStruct.addAtom(H)
saturatedStruct.addBond(bond)
if atom not in added:
added[atom] = []
added[atom].append([H, bond])
atom.decrementRadical()

# Update the atom types of the saturated structure (not sure why
# this is necessary, because saturating with H shouldn't be
# changing atom types, but it doesn't hurt anything and is not
# very expensive, so will do it anyway)
saturatedStruct.updateConnectivityValues()
saturatedStruct.sortVertices()
saturatedStruct.updateAtomTypes()
saturatedStruct.updateLonePairs()
saturatedStruct.multiplicity = 1

# Get solute descriptor estimates for saturated form of structure
soluteData = self.estimateSoluteViaGroupAdditivity(saturatedStruct)
assert soluteData is not None, "Solute data of saturated {0} of molecule {1} is None!".format(saturatedStruct, molecule)

# For each radical site, get radical correction
# Only one radical site should be considered at a time; all others
# should be saturated with hydrogen atoms
for atom in added:

# Remove the added hydrogen atoms and bond and restore the radical
for H, bond in added[atom]:
saturatedStruct.removeBond(bond)
saturatedStruct.removeAtom(H)
atom.incrementRadical()

saturatedStruct.updateConnectivityValues()

else: # non-radical species
# Generate estimate of solute data
for atom in molecule.atoms:
# Iterate over heavy (non-hydrogen) atoms
if atom.isNonHydrogen():
# Get initial solute data from main group database. Every atom must
# be found in the main abraham database
try:
self.__addGroupSoluteData(soluteData, self.groups['abraham'], molecule, {'*':atom})
except KeyError:
logging.error("Couldn't find in main abraham database:")
logging.error(molecule)
logging.error(molecule.toAdjacencyList())
raise
# Get solute data for non-atom centered groups (being found in this group
# database is optional)
try:
self.__addGroupSoluteData(soluteData, self.groups['nonacentered'], molecule, {'*':atom})
except KeyError: pass

addedToRadicals = {} # Dictionary of key = atom, value = dictionary of {H atom: bond}
addedToPairs = {} # Dictionary of key = atom, value = # lone pairs changed
saturatedStruct = molecule.copy(deep=True)

# Convert lone pairs to radicals, then saturate with H.

# 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)

# Now saturate radicals with H
if sum([atom.radicalElectrons for atom in saturatedStruct.atoms]) > 0: # radical species
saturatedStruct, addedToRadicals = self.saturateRadicals(saturatedStruct)

# Saturated structure should now have no unpaired electrons, and only "expected" lone pairs
# based on the valency
for atom in saturatedStruct.atoms:
# Iterate over heavy (non-hydrogen) atoms
if atom.isNonHydrogen():
# Get initial solute data from main group database. Every atom must
# be found in the main abraham database
try:
self.__addGroupSoluteData(soluteData, self.groups['abraham'], saturatedStruct, {'*':atom})
except KeyError:
logging.error("Couldn't find in main abraham database:")
logging.error(saturatedStruct)
logging.error(saturatedStruct.toAdjacencyList())
raise
# Get solute data for non-atom centered groups (being found in this group
# database is optional)
try:
self.__addGroupSoluteData(soluteData, self.groups['nonacentered'], saturatedStruct, {'*':atom})
except KeyError: pass

soluteData = self.removeHBonding(saturatedStruct, addedToRadicals, addedToPairs, soluteData)

return soluteData

def __addGroupSoluteData(self, soluteData, database, molecule, atom):
Expand Down
67 changes: 62 additions & 5 deletions rmgpy/data/solvationTest.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,12 +77,7 @@ def testSoluteGeneration(self):
"Test we can estimate Abraham solute parameters correctly using group contributions"

self.testCases = [
# from RMG-Java test runs by RWest (mostly in agreement with Jalan et. al. supplementary data)
['1,2-ethanediol', 'C(CO)O', 0.823, 0.685, 0.327, 2.572, 0.693, None],
# a nitrogen case
#['acetonitrile', 'CC#N', 0.9, 0.33, 0.237, 1.739, 0.04, None],
# a sulfur case
#['ethanethiol', 'CCS', 0.35, 0.24, 0.392, 2.173, 0, None]
]

for name, smiles, S, B, E, L, A, V in self.testCases:
Expand All @@ -94,6 +89,68 @@ def testSoluteGeneration(self):
self.assertAlmostEqual(soluteData.L, L, places=2)
self.assertAlmostEqual(soluteData.A, A, places=2)

def testLonePairSoluteGeneration(self):
"Test we can obtain solute parameters via group additivity for a molecule with lone pairs"
molecule=Molecule().fromAdjacencyList(
"""
CH2_singlet
multiplicity 1
1 C u0 p1 c0 {2,S} {3,S}
2 H u0 p0 c0 {1,S}
3 H u0 p0 c0 {1,S}
""")
species = Species(molecule=[molecule])
soluteData = self.database.getSoluteDataFromGroups(species)
self.assertTrue(soluteData is not None)

def testSoluteDataGenerationAmmonia(self):
"Test we can obtain solute parameters via group additivity for ammonia"
molecule=Molecule().fromAdjacencyList(
"""
1 N u0 p1 c0 {2,S} {3,S} {4,S}
2 H u0 p0 c0 {1,S}
3 H u0 p0 c0 {1,S}
4 H u0 p0 c0 {1,S}
""")
species = Species(molecule=[molecule])
soluteData = self.database.getSoluteDataFromGroups(species)
self.assertTrue(soluteData is not None)

def testSoluteDataGenerationAmide(self):
"Test that we can obtain solute parameters via group additivity for an amide"
molecule=Molecule().fromAdjacencyList(
"""
1 N u0 p1 {2,S} {3,S} {4,S}
2 H u0 {1,S}
3 C u0 {1,S} {6,S} {7,S} {8,S}
4 C u0 {1,S} {5,D} {9,S}
5 O u0 p2 {4,D}
6 H u0 {3,S}
7 H u0 {3,S}
8 H u0 {3,S}
9 H u0 {4,S}
""")
species = Species(molecule=[molecule])
soluteData = self.database.getSoluteDataFromGroups(species)
self.assertTrue(soluteData is not None)

def testRadicalandLonePairGeneration(self):
"""
Test we can obtain solute parameters via group additivity for a molecule with both lone
pairs and a radical
"""
molecule=Molecule().fromAdjacencyList(
"""
[C]OH
multiplicity 2
1 C u1 p1 c0 {2,S}
2 O u0 p2 c0 {1,S} {3,S}
3 H u0 p0 c0 {2,S}
""")
species = Species(molecule=[molecule])
soluteData = self.database.getSoluteDataFromGroups(species)
self.assertTrue(soluteData is not None)

def testCorrectionGeneration(self):
"Test we can estimate solvation thermochemistry."
self.testCases = [
Expand Down
Loading

0 comments on commit efb4d1d

Please sign in to comment.