Merge pull request #364 from bslakman/master

Fix to solvation thermo for lone pairs

Author: connie

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

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

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 = [