Skip to content

Fix to solvation thermo for lone pairs #364

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

Merged
connie merged 7 commits into from
Apr 1, 2015
Merged

Fix to solvation thermo for lone pairs #364

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
d088397
Added citations and fixed typos in liquid-phase documentation.
bslakman Jul 28, 2014
387f009
Fixed lone pairs bug in Platts group additivity for solvation.
bslakman Mar 30, 2015
388add4
Added a radical correction lookup for solvation thermo.
bslakman Mar 30, 2015
be578cf
Change thermo library commentary for liquid phase.
connie Mar 31, 2015
52e02a8
Added work in progress test for Ammonia.
connie Mar 31, 2015
1d14375
Added an amide to solvation test.
connie Mar 31, 2015
c9625a5
Remove WIP on ammonia solvation test; test now passes.
bslakman Apr 1, 2015
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
42 changes: 27 additions & 15 deletions documentation/source/users/rmg/liquids.rst
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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