diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index b95f708f02..4f6345d97b 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -48,6 +48,7 @@ convert_ring_to_sub_molecule, is_ring_partial_matched, bicyclic_decomposition_for_polyring, \ split_bicyclic_into_single_rings, saturate_ring_bonds +from rdkit import Chem ################################################################################ @@ -462,10 +463,11 @@ class KfactorParameters(object): y_solute = mole fraction of the solute in a gas phase at equilibrium in a binary dilute mixture x_solute = mole fraction of the solute in a liquid phase at equilibrium in a binary dilute mixture """ - def __init__(self, A=None, B=None, C=None, D=None, T_transition=None): + def __init__(self, A=None, B=None, C=None, D=None, T_transition=None, solvent_name=None): self.lower_T = [A, B, C] self.higher_T = D self.T_transition = T_transition # in K + self.solvent_name = solvent_name # in coolprop format class SoluteData(object): """ @@ -2222,125 +2224,121 @@ def get_solvation_correction(self, solute_data, solvent_data): Given a solute_data and solvent_data object, calculates the enthalpy, entropy, and Gibbs free energy of solvation at 298 K. Returns a SolvationCorrection object + Note: This method utilizes the LSER method for solvation correction with parameters + from the RMG-database. """ correction = SolvationCorrection(0.0, 0.0, 0.0) correction.enthalpy = self.calc_h(solute_data, solvent_data) correction.gibbs = self.calc_g(solute_data, solvent_data) correction.entropy = self.calc_s(correction.gibbs, correction.enthalpy) - return correction - - def get_Kfactor(self, delG298, delH298, delS298, solvent_name, T): - """Returns a K-factor for the input temperature given the solvation properties of a solute in a solvent - at 298 K. - Args: - delG298 (float): solvation free energy at 298 K in J/mol. - delH298 (float): solvation enthalpy at 298 K in J/mol. - delS298 (float): solvation entropy at 298 K in J/mol/K. - solvent_name (str): name of the solvent that is used in CoolProp. - T (float): input temperature in K. - - Returns: - Kfactor (float): K-factor, which is a ratio of the mole fraction of a solute in a gas-phase to - the mole fraction of a solute in a liquid-phase at equilibrium. - - Raises: - InputError: if the input temperature is above the critical temperature of the solvent. - DatabaseError: if the given solvent_name is not available in CoolProp. + return correction + + def get_solvation_correction_from_soluteML(self, solute_mol, solvent_mol, model): + """ + Dummy function to be replaced with soluteML model implementation. + """ + correction = SolvationCorrection(0.0, 0.0, 0.0) + return correction + + def get_solvation_correction_from_directML(self, solute_mol, solvent_mol, model): + """ + Dummy function to be replaced with directML model implementation. """ + correction = SolvationCorrection(0.0, 0.0, 0.0) - if solvent_name is not None: - Tc = get_critical_temperature(solvent_name) - if T < Tc: - kfactor_parameters = self.get_Kfactor_parameters(delG298, delH298, delS298, solvent_name) - A = kfactor_parameters.lower_T[0] - B = kfactor_parameters.lower_T[1] - C = kfactor_parameters.lower_T[2] - D = kfactor_parameters.higher_T - T_transition = kfactor_parameters.T_transition - rho_c = PropsSI('rhomolar_critical', solvent_name) # critical density of the solvent in mol/m^3 - rho_l = get_liquid_saturation_density(solvent_name, T) # saturated liquid phase density of the solvent, in mol/m^3 - if T < T_transition: - Kfactor = math.exp((A + B * (1 - T / Tc) ** 0.355 + C * math.exp(1 - T / Tc) * (T / Tc) ** 0.59) / (T / Tc)) - else: - Kfactor = math.exp(D * (rho_l / rho_c -1) / (T / Tc)) + return correction + + def check_solvent_in_initial_species(self, rmg, solvent_structure): + """ + Given the instance of RMG class and the solvent_structure, it checks whether the solvent is listed as one + of the initial species. + If the SMILES / adjacency list for all the solvents exist in the solvent library, it uses the solvent's + molecular structure to determine whether the species is the solvent or not. + If the solvent library does not have SMILES / adjacency list, then it uses the solvent's string name + to determine whether the species is the solvent or not + """ + for spec in rmg.initial_species: + if solvent_structure is not None: + spec.is_solvent = spec.is_isomorphic(solvent_structure) else: - raise InputError("The input temperature {0} K cannot be greater than " - "or equal to the critical temperature, {1} K".format(T, Tc)) - else: - raise DatabaseError("K-factor calculation or temperature-dependent solvation free energy calculation " - f"is not available for the given solvent name: {solvent_name}") - return Kfactor - - def get_T_dep_solvation_energy_from_LSER_298(self, solute_data, solvent_data, T): - """Returns solvation free energy and K-factor for the input temperature based on the 298 K solvation properties - calculated from the LSER method. - - Args: - solute_data: :class:`SoluteData` object `solute_data`. - solvent_data: :class:`SolventData` object `solvent_data`. - T (float): input temperature in K. - - Returns: - delG (float): solvation free energy at the input temperature in J/mol. - Kfactor (float): K-factor at the input temperature. K-factor is defined as a ratio of the mole fraction - of a solute in a gas-phase to the mole fraction of a solute in a liquid-phase at equilibrium. - - Raises: - DatabaseError: if `solute_data.name_in_coolprop` is None or `solute_data` has any missing Abarham or - Mintz solvent parameters. + spec.is_solvent = rmg.solvent == spec.label + if not any([spec.is_solvent for spec in rmg.initial_species]): + if solvent_structure is not None: + logging.info('One of the initial species must be the solvent') + logging.warning("Solvent is not an initial species") + else: + logging.info('One of the initial species must be the solvent with the same string name') + logging.warning("Solvent is not an initial species with the same string name") +class MLSolvation: + """ + A class for solvation correction. + """ + def __init__(self, T_dep: bool = False, method: str = "SoluteGC", model_path: str = ""): """ + model path: path to the directory where the machine learning model is stored + # ex) model_path = "RMG-database/input/thermo/ml/solvation" + T_dep: whether T_dep solvation correction is used + method: the method used for solvation correction, e.g., "SoluteGC", "SoluteML" or "DirectML" + """ + self.model_path = model_path + self.T_dep = T_dep + self.method = method + self.model = self.load_model(model_path) if method != "SoluteGC" else None - solvent_name = solvent_data.name_in_coolprop - # check whether all solvent parameters exist - solvent_param_list = [solvent_data.s_g, solvent_data.b_g, solvent_data.e_g, - solvent_data.l_g, solvent_data.a_g, solvent_data.c_g, - solvent_data.s_h, solvent_data.b_h, solvent_data.e_h, - solvent_data.l_h, solvent_data.a_h, solvent_data.c_h, - ] - param_found_list = [param is not None for param in solvent_param_list] - - if solvent_name is not None and all(param_found_list): - correction = self.get_solvation_correction(solute_data, solvent_data) - delG298 = correction.gibbs # in J/mol - delH298 = correction.enthalpy # in J/mol - delS298 = correction.entropy # in J/mol/K - else: - raise DatabaseError("K-factor parameter calculation is not available for the solvent " - "whose `name_in_coolprop` is None or that is missing the Abraham and Mintz " - "solvent parameters.") - - return self.get_T_dep_solvation_energy_from_input_298(delG298, delH298, delS298, solvent_name, T) - - def get_T_dep_solvation_energy_from_input_298(self, delG298, delH298, delS298, solvent_name, T): - """Returns solvation free energy and K-factor for the input temperature based on the given solvation properties - values at 298 K. - - Args: - delG298 (float): solvation free energy at 298 K in J/mol. - delH298 (float): solvation enthalpy at 298 K in J/mol. - delS298 (float): solvation entropy at 298 K in J/mol/K. - solvent_name (str): name of the solvent that is used in CoolProp. - T (float): input temperature in K. + def load_model(self, model_path): + # Need to implement model loading code (e.g., joblib.load, torch.load, etc.) + print(f"[NOTICE] Dummy ML model loaded from: {model_path}") + return None # Need to return something like "return joblib.load(os.path.join(model_path, "model.pkl"))" - Returns: - delG (float): solvation free energy at the input temperature in J/mol. - Kfactor (float): K-factor at the input temperature. K-factor is defined as a ratio of the mole fraction - of a solute in a gas-phase to the mole fraction of a solute in a liquid-phase at equilibrium - kH (float): the Henry's law constant at the input temperature. kH is defined as the ratio of the pressure - of a solute in the gas-phase to the concentration of a solute in the liquid-phase at equilibrium. + def generate_solvation_model(self, spc, solvent_name): """ + Returns a TDepModel or StaticModel depending on T_dep flag. + """ + correction = self.get_solvation_correction(spc, solvent_name) - Kfactor = self.get_Kfactor(delG298, delH298, delS298, solvent_name, T) - rho_g = get_gas_saturation_density(solvent_name, T) # saturated gas phase density of the solvent, in mol/m^3 - rho_l = get_liquid_saturation_density(solvent_name, T) # saturated liquid phase density of the solvent, in mol/m^3 - Psat_g = get_gas_saturation_pressure(solvent_name, T) - kH = Kfactor * Psat_g / rho_l # Henry's law constant as a fraction of the gas-phase partial pressure of solute to the liquid-phase concentration of solute - delG = constants.R * T * math.log(Kfactor * rho_g / (rho_l)) # in J/mol - return delG, Kfactor, kH + dH298 = correction.enthalpy # J/mol + dS298 = correction.entropy # J/mol/K + dG298 = correction.gibbs # J/mol + if self.T_dep: + kfactor_parameters = self.get_Kfactor_parameters(dG298, dH298, dS298, solvent_name) + return TDepModel(kfactor_parameters, solvent_name) + else: + return StaticModel(dH298, dS298) + + def get_solvation_correction(self, spc, solvent_name): + """ + Returns solvation correction at 298K based on 3 methods: SoluteGC, SoluteML, or DirectML. + """ + # Get solute_data and solvent_data + from rmgpy.data.rmg import get_db # if we import get_db at the top, it will cause circular import error + solvation_database = get_db('solvation') + solute_data = solvation_database.get_solute_data(spc) + solvent_data = solvation_database.get_solvent_data(solvent_name) + + if self.method == "SoluteGC": + return solvation_database.get_solvation_correction(solute_data, solvent_data) + + elif self.method in ["SoluteML", "DirectML"]: + + # Get solute_mol and solvent_mol + solvent_species = solvation_database.get_solvent_structure(solvent_name)[0] + solvent_smiles = solvent_species.smiles + solvent_mol = Chem.MolFromSmiles(solvent_smiles) + + solute_smiles = spc.smiles + solute_mol = Chem.MolFromSmiles(solute_smiles) + + if self.method == "SoluteML": + return solvation_database.get_solvation_correction_from_soluteML(solute_mol, solvent_mol, self.model) + elif self.method == "DirectML": + return solvation_database.get_solvation_correction_from_directML(solute_mol, solvent_mol, self.model) + else: + raise ValueError(f"Unknown solvation method: {self.method}") + def get_Kfactor_parameters(self, delG298, delH298, delS298, solvent_name, T_trans_factor=0.75): """Returns fitted parameters for the K-factor piecewise functions given the solvation properties at 298 K. @@ -2372,6 +2370,13 @@ def get_Kfactor_parameters(self, delG298, delH298, delS298, solvent_name, T_tran second piecewise function. """ + # Find solvent name in CoolProp database + from rmgpy.data.rmg import get_db # if we import get_db at the top, it will cause circular import error + solvation_database = get_db('solvation') + solvent_data = solvation_database.get_solvent_data(solvent_name) + solvent_name = solvent_data.name_in_coolprop + if solvent_name is None: + raise DatabaseError("Cannot find the solvent name in CoolProp database.") Tc = get_critical_temperature(solvent_name) T_transition = Tc * T_trans_factor # T_trans_factor is empirically set to 0.75 by default @@ -2427,27 +2432,96 @@ def get_Kfactor_parameters(self, delG298, delH298, delS298, solvent_name, T_tran kfactor_parameters.lower_T = [float(param[0]), float(param[1]), float(param[2])] kfactor_parameters.higher_T = float(param[3]) kfactor_parameters.T_transition = T_transition + kfactor_parameters.solvent_name = solvent_name # name in coolprop return kfactor_parameters - def check_solvent_in_initial_species(self, rmg, solvent_structure): + +class TDepModel: + def __init__(self, kfactor_parameters, solvent_name, T_trans_factor=0.75): """ - Given the instance of RMG class and the solvent_structure, it checks whether the solvent is listed as one - of the initial species. - If the SMILES / adjacency list for all the solvents exist in the solvent library, it uses the solvent's - molecular structure to determine whether the species is the solvent or not. - If the solvent library does not have SMILES / adjacency list, then it uses the solvent's string name - to determine whether the species is the solvent or not + solvent_name (str): name of the solvent + solvent_name_in_coolprop (str): name of the solvent in CoolProp database + A, B, C, D (float): coefficients for the K-factor piecewise function + T_trans_factor (float): fraction of Tc for transition between equations """ - for spec in rmg.initial_species: - if solvent_structure is not None: - spec.is_solvent = spec.is_isomorphic(solvent_structure) - else: - spec.is_solvent = rmg.solvent == spec.label - if not any([spec.is_solvent for spec in rmg.initial_species]): - if solvent_structure is not None: - logging.info('One of the initial species must be the solvent') - logging.warning("Solvent is not an initial species") + self.A = kfactor_parameters.lower_T[0] + self.B = kfactor_parameters.lower_T[1] + self.C = kfactor_parameters.lower_T[2] + self.D = kfactor_parameters.higher_T + self.T_trans_factor = T_trans_factor + self.solvent_name = solvent_name + self.solvent_name_in_coolprop = kfactor_parameters.solvent_name # name in coolprop + + def get_Kfactor(self, T): + """ + Returns a K-factor for the input temperature given the solvation properties of a solute in a solvent + at 298 K. + + Args: + T (float): input temperature in K. + + Returns: + Kfactor (float): K-factor, which is a ratio of the mole fraction of a solute in a gas-phase to + the mole fraction of a solute in a liquid-phase at equilibrium. + + Raises: + InputError: if the input temperature is above the critical temperature of the solvent. + DatabaseError: if the given solvent_name is not available in CoolProp. + + """ + solvent_name = self.solvent_name_in_coolprop + if solvent_name is not None: + Tc = get_critical_temperature(solvent_name) + if T < Tc: + T_transition = self.T_trans_factor * Tc + rho_c = PropsSI('rhomolar_critical', solvent_name) # critical density of the solvent in mol/m^3 + rho_l = get_liquid_saturation_density(solvent_name, T) # saturated liquid phase density of the solvent, in mol/m^3 + if T < T_transition: + Kfactor = math.exp((self.A + self.B * (1 - T / Tc) ** 0.355 + self.C * math.exp(1 - T / Tc) * (T / Tc) ** 0.59) / (T / Tc)) + else: + Kfactor = math.exp(self.D * (rho_l / rho_c -1) / (T / Tc)) else: - logging.info('One of the initial species must be the solvent with the same string name') - logging.warning("Solvent is not an initial species with the same string name") + raise InputError("The input temperature {0} K cannot be greater than " + "or equal to the critical temperature, {1} K".format(T, Tc)) + else: + raise DatabaseError("K-factor calculation or temperature-dependent solvation free energy calculation " + f"is not available for the given solvent name: {solvent_name}") + return Kfactor + + def get_T_dep_solvation_energy(self, T): + """Returns solvation free energy for the input temperature based on the given solvation properties + values at 298 K. + + Args: + T (float): input temperature in K. + + Returns: + delG (float): solvation free energy at the input temperature in J/mol. + Kfactor (float): K-factor at the input temperature. K-factor is defined as a ratio of the mole fraction + of a solute in a gas-phase to the mole fraction of a solute in a liquid-phase at equilibrium + kH (float): the Henry's law constant at the input temperature. kH is defined as the ratio of the pressure + of a solute in the gas-phase to the concentration of a solute in the liquid-phase at equilibrium. + """ + solvent_name = self.solvent_name_in_coolprop + Kfactor = self.get_Kfactor(T) + rho_g = get_gas_saturation_density(solvent_name, T) # saturated gas phase density of the solvent, in mol/m^3 + rho_l = get_liquid_saturation_density(solvent_name, T) # saturated liquid phase density of the solvent, in mol/m^3 + Psat_g = get_gas_saturation_pressure(solvent_name, T) + kH = Kfactor * Psat_g / rho_l # Henry's law constant as a fraction of the gas-phase partial pressure of solute to the liquid-phase concentration of solute + delG = constants.R * T * math.log(Kfactor * rho_g / (rho_l)) # in J/mol + return delG, Kfactor, kH + + def get_free_energy_of_solvation(self, T): + delG, Kfactor, kH = self.get_T_dep_solvation_energy(T) + return delG + + +class StaticModel: + def __init__(self, dH298, dS298): + self.dH298 = dH298 + self.dS298 = dS298 + + def get_free_energy_of_solvation(self, T): + delG = self.dH298 - T * self.dS298 + return delG \ No newline at end of file diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index b8ab3951c8..767f16a293 100644 --- a/rmgpy/reaction.py +++ b/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 diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 117c19644c..50930673a0 100644 --- a/rmgpy/rmg/input.py +++ b/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: diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index fb4ae33ef8..892f609199 100644 --- a/rmgpy/rmg/main.py +++ b/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): """ diff --git a/rmgpy/species.pxd b/rmgpy/species.pxd index 0e30ce2b0e..09911d927c 100644 --- a/rmgpy/species.pxd +++ b/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 diff --git a/rmgpy/species.py b/rmgpy/species.py index c17bdd182a..009ab35687 100644 --- a/rmgpy/species.py +++ b/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. diff --git a/rmgpy/thermo/thermoengine.py b/rmgpy/thermo/thermoengine.py index e8840770a5..54507d8f14 100644 --- a/rmgpy/thermo/thermoengine.py +++ b/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)) \ No newline at end of file diff --git a/test/rmgpy/data/solvationTest.py b/test/rmgpy/data/solvationTest.py index 72d0926d55..7055065b66 100644 --- a/test/rmgpy/data/solvationTest.py +++ b/test/rmgpy/data/solvationTest.py @@ -35,6 +35,8 @@ SoluteData, SolvationDatabase, SolventLibrary, + MLSolvation, + TDepModel, get_critical_temperature, get_liquid_saturation_density, get_gas_saturation_density, @@ -45,12 +47,20 @@ from rmgpy.exceptions import InputError import pytest +class DummyRMGDatabase: + """ + Dummy class used in testing to satisfy the get_db('solvation') call within MLSolvation. + """ + def __init__(self, solvation_database): + self.solvation = solvation_database class TestSoluteDatabase: @classmethod def setup_class(cls): cls.database = SolvationDatabase() cls.database.load(os.path.join(settings["database.directory"], "solvation")) + import rmgpy.data.rmg + rmgpy.data.rmg.database = DummyRMGDatabase(cls.database) @classmethod def teardown_class(cls): @@ -377,71 +387,79 @@ def test_Kfactor_parameters(self): species = Species().from_smiles("CCC(C)=O") # 2-Butanone for a solute solute_data = self.database.get_solute_data(species) solvent_data = self.database.get_solvent_data("water") - correction = self.database.get_solvation_correction(solute_data, solvent_data) - delG298 = correction.gibbs # in J/mol - delH298 = correction.enthalpy # in J/mol - delS298 = correction.entropy # in J/mol/K + solvent_name = solvent_data.name_in_coolprop - kfactor_parameters = self.database.get_Kfactor_parameters(delG298, delH298, delS298, solvent_name) - assert round(abs(kfactor_parameters.lower_T[0] - -9.780), 3) == 0 # check up to 3 decimal places - assert round(abs(kfactor_parameters.lower_T[1] - 0.492), 3) == 0 - assert round(abs(kfactor_parameters.lower_T[2] - 10.485), 3) == 0 - assert round(abs(kfactor_parameters.higher_T - 1.147), 3) == 0 - assert round(abs(kfactor_parameters.T_transition - 485.3), 1) == 0 + solvation = MLSolvation(T_dep=True, method="SoluteGC") + solvation_model = solvation.generate_solvation_model(species, solvent_name) + + # Confirm that we get the expected model type + assert isinstance(solvation_model, TDepModel) + + # Extract fitted parameters and check against expected values + assert round(abs(solvation_model.A - -9.780), 3) == 0 # check up to 3 decimal places + assert round(abs(solvation_model.B - 0.492), 3) == 0 + assert round(abs(solvation_model.C - 10.485), 3) == 0 + assert round(abs(solvation_model.D - 1.147), 3) == 0 + + Tc = get_critical_temperature(solvation_model.solvent_name_in_coolprop) + calculated_Tc = solvation_model.T_trans_factor * Tc + assert round(abs(calculated_Tc - 485.3), 1) == 0 + # check that DatabaseError is raised when the solvent's name_in_coolprop is None - solvent_data = self.database.get_solvent_data("chloroform") - solvent_name = solvent_data.name_in_coolprop with pytest.raises(DatabaseError): - self.database.get_Kfactor_parameters( - delG298, - delH298, - delS298, - solvent_name, - ) + solvation.generate_solvation_model(species, "chloroform") def test_Tdep_solvation_calculation(self): """ Test we can calculate the temperature dependent solvation free energy and K-factor using both `get_T_dep_solvation_energy_from_LSER_298` and `get_T_dep_solvation_energy_from_input_298` methods. """ - # First, test `get_T_dep_solvation_energy_from_LSER_298` method. species = Species().from_smiles("CCC1=CC=CC=C1") # ethylbenzene species.generate_resonance_structures() - solute_data = self.database.get_solute_data(species) - solvent_data = self.database.get_solvent_data("benzene") + solvent_name = "benzene" + + solvation = MLSolvation(T_dep=True, method="SoluteGC") + solvation_model = solvation.generate_solvation_model(species, solvent_name) + + # Confirm that we get the expected model type + assert isinstance(solvation_model, TDepModel) + + # Evaluate at T = 500 K T = 500 # in K - delG, Kfactor, kH = self.database.get_T_dep_solvation_energy_from_LSER_298(solute_data, solvent_data, T) + delG, Kfactor, kH = solvation_model.get_T_dep_solvation_energy(T) assert round(abs(Kfactor - 0.403), 3) == 0 assert round(abs(delG / 1000 - -13.59), 2) == 0 # delG is in J/mol + # For temperature greater than or equal to the critical temperature of the solvent, # it should raise InputError T = 1000 with pytest.raises(InputError): - self.database.get_T_dep_solvation_energy_from_LSER_298( - solute_data, - solvent_data, - T, - ) + solvation_model.get_free_energy_of_solvation(T) - # Now test `get_T_dep_solvation_energy_from_input_298` method. + def test_TDepModel(self): delG298 = -23570 # in J/mol delH298 = -40612 # in J/mol delS298 = (delH298 - delG298) / 298 # in J/mol/K solvent_name = "benzene" + + solvation = MLSolvation(T_dep=True, method="SoluteGC") + kfactor_parameters = solvation.get_Kfactor_parameters(delG298, delH298, delS298, solvent_name) + solvation_model = TDepModel( + kfactor_parameters=kfactor_parameters, + solvent_name=solvent_name + ) + + # Test at T = 500 K T = 500 # in K - delG, Kfactor, kH = self.database.get_T_dep_solvation_energy_from_input_298(delG298, delH298, delS298, solvent_name, T) + delG, Kfactor, kH = solvation_model.get_T_dep_solvation_energy(T) assert round(abs(Kfactor - 0.567), 3) == 0 assert round(abs(delG / 1000 - -12.18), 2) == 0 # delG is in J/mol + # test that it raises InputError for T above the critical temperature T = 1000 with pytest.raises(InputError): - self.database.get_T_dep_solvation_energy_from_input_298( - delG298, - delH298, - delS298, - solvent_name, - T, - ) + solvation_model.get_free_energy_of_solvation(T) + @pytest.mark.skip(reason="Skip for Electrochem PR.") def test_initial_species(self):