diff --git a/documentation/source/users/rmg/liquids.rst b/documentation/source/users/rmg/liquids.rst index d13f818613..188121651c 100644 --- a/documentation/source/users/rmg/liquids.rst +++ b/documentation/source/users/rmg/liquids.rst @@ -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 @@ -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). diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index a55076d1a0..35a305f765 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -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 ################################################################################ @@ -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): """ @@ -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): """ @@ -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 @@ -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): diff --git a/rmgpy/data/solvationTest.py b/rmgpy/data/solvationTest.py index cb44ceba33..2f65a6c923 100644 --- a/rmgpy/data/solvationTest.py +++ b/rmgpy/data/solvationTest.py @@ -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: @@ -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 = [ diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index d9aa4142d4..9c0bad5e60 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -688,7 +688,7 @@ def getThermoDataFromLibraries(self, species, trainingSet=None): if thermoData is not None: assert len(thermoData) == 3, "thermoData should be a tuple at this point" if rmgpy.rmg.main.solvent is not None and trainingSet is None: - thermoData[0].comment += 'Thermo library "corrected": ' + label + thermoData[0].comment += 'Thermo library corrected for liquid phase: ' + label else: thermoData[0].comment += 'Thermo library: ' + label return thermoData