Add T-Dependent and ML Options for Solvation Corrections

rmgpy/reaction.py

@@ -757,6 +757,17 @@ def get_equilibrium_constant(self, T, potential=0., type='Kc', surface_site_dens
         )
         # Use free energy of reaction to calculate Ka
         dGrxn = self.get_free_energy_of_reaction(T, potential)
+
+        # Apply solvation correction
+        first_reactant = self.reactants[0] if self.reactants else None
+        if first_reactant and first_reactant.has_solvation_thermo():
+            try:
+                ddGsolv = self.get_solvation_free_energy_of_reaction(T)
+                dGrxn += ddGsolv
+                logging.debug("Applied solvation correction: ΔΔG_solv = {:.2f} J/mol for reaction {!s}".format(ddGsolv, self))
+            except Exception as e:
+                logging.info("Failed to compute solvation correction for reaction {!s}: {}".format(self, e))
+
         K = np.exp(-dGrxn / constants.R / T)
         # Convert Ka to Kc or Kp if specified
         # Assume a pressure of 1e5 Pa for gas phase species
@@ -821,6 +832,27 @@ def get_equilibrium_constant(self, T, potential=0., type='Kc', surface_site_dens
             raise ReactionError('Got equilibrium constant of 0')
         return K
 
+    def get_solvation_free_energy_of_reaction(self, T):
+        """
+        Retrun the solvation free energy of reaction in J/mol evaluated at temperature T in K.
+        """
+        ddGsolv = 0.0
+        for reactant in self.reactants:
+            try:
+                ddGsolv -= reactant.get_free_energy_of_solvation(T)
+            except Exception as e:
+                logging.error("Problem with reactant {!r} in reaction {!s}: {}".format(reactant.label, self, e))
+                raise
+
+        for product in self.products:
+            try:
+                ddGsolv += product.get_free_energy_of_solvation(T)
+            except Exception as e:
+                logging.error("Problem with product {!r} in reaction {!s}: {}".format(product.label, self, e))
+                raise
+
+        return ddGsolv
+
     def get_enthalpies_of_reaction(self, Tlist):
         """
         Return the enthalpies of reaction in J/mol evaluated at temperatures

rmgpy/rmg/input.py

@@ -1305,6 +1305,17 @@ def ml_estimator(thermo=True,
         logging.warning('"onlyCyclics" should be True when "onlyHeterocyclics" is True. '
                         'Machine learning estimator is restricted to only heterocyclic species thermo')
 
+def ml_solvation(T_dep=False, method="SoluteGC", name='solvation'):
+
+    from rmgpy.data.solvation import MLSolvation
+    if method != "SoluteGC":
+        model_path = os.path.join(settings['database.directory'], 'thermo', 'ml', name) # Need to add ML model at RMG-database/input/thermo/ml/solvation
+        if not os.path.exists(model_path):
+            raise InputError('Cannot find ML models folder {}'.format(model_path))
+    else:
+        model_path = "" # SoluteGC does not need model files
+
+    rmg.ml_solvation = MLSolvation(T_dep=T_dep, method=method, model_path=model_path)
 
 def pressure_dependence(
         method,
@@ -1565,6 +1576,7 @@ def read_input_file(path, rmg0):
         'model': model,
         'quantumMechanics': quantum_mechanics,
         'mlEstimator': ml_estimator,
+        'mlSolvation': ml_solvation,
         'pressureDependence': pressure_dependence,
         'options': options,
         'generatedSpeciesConstraints': generated_species_constraints,
@@ -1645,6 +1657,7 @@ def read_thermo_input_file(path, rmg0):
         'adjacencyList': adjacency_list,
         'quantumMechanics': quantum_mechanics,
         'mlEstimator': ml_estimator,
+        'mlSolvation': ml_solvation,
     }
 
     try:
@@ -1851,6 +1864,8 @@ def get_input(name):
             return rmg.quantum_mechanics
         elif name == 'ml_estimator':
             return rmg.ml_estimator, rmg.ml_settings
+        elif name == 'ml_solvation':
+            return rmg.ml_solvation
         elif name == 'thermo_central_database':
             return rmg.thermo_central_database
         else:

rmgpy/rmg/main.py

@@ -1456,47 +1456,51 @@ def check_model(self):
 
         # Check all core reactions (in both directions) for collision limit violation
         violators = []
-        for rxn in self.reaction_model.core.reactions:
-            if rxn.is_surface_reaction():
-                # Don't check collision limits for surface reactions.
-                continue
-            violator_list = rxn.check_collision_limit_violation(t_min=self.Tmin, t_max=self.Tmax, p_min=self.Pmin, p_max=self.Pmax)
-            if violator_list:
-                violators.extend(violator_list)
-        # Whether or not violators were found, rename 'collision_rate_violators.log' if it exists
-        new_file = os.path.join(self.output_directory, "collision_rate_violators.log")
-        old_file = os.path.join(self.output_directory, "collision_rate_violators_OLD.log")
-        if os.path.isfile(new_file):
-            # If there are no violators, yet the violators log exists (probably from a previous run
-            # in the same folder), rename it.
-            if os.path.isfile(old_file):
-                os.remove(old_file)
-            os.rename(new_file, old_file)
-        if violators:
-            logging.info("\n")
-            logging.warning(
-                "{0} CORE reactions violate the collision rate limit!"
-                "\nSee the 'collision_rate_violators.log' for details.\n\n".format(len(violators))
-            )
-            with open(new_file, "w") as violators_f:
-                violators_f.write(
-                    "*** Collision rate limit violators report ***\n"
-                    '"Violation factor" is the ratio of the rate coefficient to the collision limit'
-                    " rate at the relevant conditions\n\n"
+
+        if not self.solvent:
+            for rxn in self.reaction_model.core.reactions:
+                if rxn.is_surface_reaction():
+                    # Don't check collision limits for surface reactions.
+                    continue
+                violator_list = rxn.check_collision_limit_violation(t_min=self.Tmin, t_max=self.Tmax, p_min=self.Pmin, p_max=self.Pmax)
+                if violator_list:
+                    violators.extend(violator_list)
+            # Whether or not violators were found, rename 'collision_rate_violators.log' if it exists
+            new_file = os.path.join(self.output_directory, "collision_rate_violators.log")
+            old_file = os.path.join(self.output_directory, "collision_rate_violators_OLD.log")
+            if os.path.isfile(new_file):
+                # If there are no violators, yet the violators log exists (probably from a previous run
+                # in the same folder), rename it.
+                if os.path.isfile(old_file):
+                    os.remove(old_file)
+                os.rename(new_file, old_file)
+            if violators:
+                logging.info("\n")
+                logging.warning(
+                    "{0} CORE reactions violate the collision rate limit!"
+                    "\nSee the 'collision_rate_violators.log' for details.\n\n".format(len(violators))
                 )
-                for violator in violators:
-                    if isinstance(violator[0].kinetics, (ThirdBody,Troe)):
-                        rxn_string = violator[0].to_chemkin(self.reaction_model.core.species)
-                    else:
-                        rxn_string = violator[0].to_chemkin()
-                    direction = violator[1]
-                    ratio = violator[2]
-                    condition = violator[3]
+                with open(new_file, "w") as violators_f:
                     violators_f.write(
-                        f"{rxn_string}\n" f"Direction: {direction}\n" f"Violation factor: {ratio:.2f}\n" f"Violation condition: {condition}\n\n\n"
+                        "*** Collision rate limit violators report ***\n"
+                        '"Violation factor" is the ratio of the rate coefficient to the collision limit'
+                        " rate at the relevant conditions\n\n"
                     )
+                    for violator in violators:
+                        if isinstance(violator[0].kinetics, (ThirdBody,Troe)):
+                            rxn_string = violator[0].to_chemkin(self.reaction_model.core.species)
+                        else:
+                            rxn_string = violator[0].to_chemkin()
+                        direction = violator[1]
+                        ratio = violator[2]
+                        condition = violator[3]
+                        violators_f.write(
+                            f"{rxn_string}\n" f"Direction: {direction}\n" f"Violation factor: {ratio:.2f}\n" f"Violation condition: {condition}\n\n\n"
+                        )
+            else:
+                logging.info("No collision rate violators found in the model's core.")
         else:
-            logging.info("No collision rate violators found in the model's core.")
+            logging.info("Skipping collision limit checks: liquid phase reactions.")
 
     def initialize_seed_mech(self):
         """

rmgpy/species.pxd

@@ -42,6 +42,7 @@ cdef class Species:
     cdef public int index
     cdef public str label
     cdef public object thermo
+    cdef public object solvationthermo
     cdef public Conformer conformer
     cdef public object transport_data
     cdef public list molecule

rmgpy/species.py

@@ -100,6 +100,7 @@ def __init__(self, index=-1, label='', thermo=None, conformer=None, molecule=Non
         self.index = index
         self.label = label
         self.thermo = thermo
+        self.solvationthermo = None
         self.conformer = conformer
         self.molecule = molecule or []
         self.transport_data = transport_data
@@ -495,6 +496,13 @@ def has_thermo(self):
         ``False`` otherwise.
         """
         return self.thermo is not None
+
+    def has_solvation_thermo(self):
+        """
+        Return ``True`` if the species has solvationthermo, or 
+        ``False`` otherwise.
+        """
+        return self.solvationthermo is not None
 
     def contains_surface_site(self):
         """
@@ -605,7 +613,13 @@ def get_free_energy(self, T):
             raise Exception('Unable to calculate free energy for species {0!r}: '
                             'no thermo or statmech data available.'.format(self.label))
         return G
-
+
+    def get_free_energy_of_solvation(self, T):
+        """
+        Return the Gibbs free energy of solvation in J/mol for the species at temperature T [K].
+        """
+        return self.get_solvation_thermo_data().get_free_energy_of_solvation(T)
+
     def get_sum_of_states(self, e_list):
         """
         Return the sum of states :math:`N(E)` at the specified energies `e_list`
@@ -821,6 +835,14 @@ def get_thermo_data(self, solvent_name=''):
 
         return self.thermo
 
+    def get_solvation_thermo_data(self):
+        """
+        Returns a solvationthermo object (TDepModel or StaticModel) of the current Species object.
+        """
+        if self.has_solvation_thermo():
+            return self.solvationthermo
+        raise Exception("No solvationthermo available for species {}".format(self.label))
+
     def generate_transport_data(self):
         """
         Generate the transport_data parameters for the species.

rmgpy/thermo/thermoengine.py

@@ -39,6 +39,7 @@
 from rmgpy.molecule import Molecule
 from rmgpy.molecule.fragment import Fragment
 
+from rdkit import Chem
 
 def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''):
     """
@@ -66,11 +67,6 @@ def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''):
     if solvent_data and not "Liquid thermo library" in thermo0.comment:
         solvation_database = get_db('solvation')
         solute_data = solvation_database.get_solute_data(spc)
-        solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data)
-        # correction is added to the entropy and enthalpy
-        wilhoit.S0.value_si = (wilhoit.S0.value_si + solvation_correction.entropy)
-        wilhoit.H0.value_si = (wilhoit.H0.value_si + solvation_correction.enthalpy)
-        wilhoit.comment += f' + Solvation correction (H={solvation_correction.enthalpy/1e3:+.0f}kJ/mol;S={solvation_correction.entropy:+.0f}J/mol/K) with {solvent_name} as solvent and solute estimated using {solute_data.comment}'
 
     # Compute E0 by extrapolation to 0 K
     if spc.conformer is None:
@@ -182,3 +178,18 @@ def submit(spc, solvent_name=''):
 
     """
     spc.thermo = evaluator(spc, solvent_name=solvent_name)
+
+    # generate solvationthermo if needed
+    if solvent_name and spc.thermo and "Liquid thermo library" not in spc.thermo.comment:
+        solvation_database = get_db('solvation')
+        if solvation_database:
+            try:
+                from rmgpy.rmg.input import get_input
+                ml_solvation = get_input("ml_solvation")
+                solvationthermo = ml_solvation.generate_solvation_model(
+                    spc = spc,
+                    solvent_name = solvent_name
+                )
+                spc.solvationthermo = solvationthermo
+            except Exception as e:
+                    logging.warning("Failed to generate solvation thermo for {}: {}".format(spc.label, e))