From 47bf97b879e8c1b63e4f263d80fcc98405f1b27a Mon Sep 17 00:00:00 2001 From: Bonhyeok Koo Date: Tue, 27 May 2025 18:28:16 -0400 Subject: [PATCH 1/6] WIP: Add two ML options for solvation enthalpy and entropy correction --- rmgpy/data/solvation.py | 35 ++++++++++++++++++++++++++++++++++- rmgpy/thermo/thermoengine.py | 28 +++++++++++++++++++++++++++- 2 files changed, 61 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index b95f708f02a..a723ac34e8d 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -2221,13 +2221,46 @@ 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 + 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_solvation_correction_ml_LSER(self, solute_mol, solvent_mol): + """ + Given a solute_mol and solvent_mol 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 + predicted from the ML model. + """ + correction = SolvationCorrection(0.0, 0.0, 0.0) + + ### ML MODEL HERE ### + + correction.enthalpy = 0 + correction.gibbs = 0 + correction.entropy = 0 + return correction + + def get_solvation_correction_ml(self, solute_mol, solvent_mol): + """ + Given a solute_mol and solvent_mol object, calculates the enthalpy, entropy, + and Gibbs free energy of solvation at 298 K using a machine learning model. + Returns a SolvationCorrection object. + """ + correction = SolvationCorrection(0.0, 0.0, 0.0) + + ### ML MODEL HERE ### + + correction.enthalpy = 0 + correction.gibbs = 0 + correction.entropy = 0 + 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 diff --git a/rmgpy/thermo/thermoengine.py b/rmgpy/thermo/thermoengine.py index e8840770a57..715d5ea5698 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,7 +67,32 @@ 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) + + if False: # Need to add option variable "use_ml_solvation" to use ML solvation correction + if False: # Need to add option variable "use_ml_LSER_parameters" to use ML calculated parameters + + # Case 1. Use LSER with parameters from database to predict solvation correction + solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data) + else: + # Case 2. Use LSER with parameters from ML model to predict solvation correction + 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) + solvation_correction = solvation_database.get_solvation_correction_ml_LSER(solute_mol, solvent_mol) + else: + # Case 3. Use direct ML model to predict solvation correction + 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) + + solvation_correction = solvation_database.get_solvation_correction_ml(solute_mol, solvent_mol) + # 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) From f08c78cecb14c0ff689820eb3f1b1df466100fe9 Mon Sep 17 00:00:00 2001 From: Bonhyeok Koo Date: Wed, 28 May 2025 17:14:10 -0400 Subject: [PATCH 2/6] Add ML solvation correction with fallback to LSER --- rmgpy/data/solvation.py | 61 ++++++++++++++++++------------------ rmgpy/rmg/input.py | 15 +++++++++ rmgpy/thermo/thermoengine.py | 44 +++++++++++++------------- 3 files changed, 67 insertions(+), 53 deletions(-) diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index a723ac34e8d..7b149ccf1e9 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 ################################################################################ @@ -2231,37 +2232,6 @@ def get_solvation_correction(self, solute_data, solvent_data): correction.entropy = self.calc_s(correction.gibbs, correction.enthalpy) return correction - def get_solvation_correction_ml_LSER(self, solute_mol, solvent_mol): - """ - Given a solute_mol and solvent_mol 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 - predicted from the ML model. - """ - correction = SolvationCorrection(0.0, 0.0, 0.0) - - ### ML MODEL HERE ### - - correction.enthalpy = 0 - correction.gibbs = 0 - correction.entropy = 0 - return correction - - def get_solvation_correction_ml(self, solute_mol, solvent_mol): - """ - Given a solute_mol and solvent_mol object, calculates the enthalpy, entropy, - and Gibbs free energy of solvation at 298 K using a machine learning model. - Returns a SolvationCorrection object. - """ - correction = SolvationCorrection(0.0, 0.0, 0.0) - - ### ML MODEL HERE ### - - correction.enthalpy = 0 - correction.gibbs = 0 - correction.entropy = 0 - 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. @@ -2484,3 +2454,32 @@ def check_solvent_in_initial_species(self, rmg, solvent_structure): 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 dummy class for machine learning-based solvation correction. + """ + def __init__(self, model_path: str): + # ex) model_path = "RMG-database/input/thermo/ml/solvation" + self.model_path = model_path + self.model = self.load_model(model_path) + + 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"))" + + def get_solvation_correction(self, solute_mol, solvent_mol): + """ + Given a solute_mol and solvent_mol object, calculates the enthalpy, entropy, + and Gibbs free energy of solvation at 298 K using a machine learning model. + Returns a SolvationCorrection object. + """ + # Need to implement: solute_mol, solvent_mol → feature vector → ML prediction + print(f"[NOTICE] Dummy ML model utilized") + enthalpy = 0.0 #self.model.predict([...])[0] + gibbs = 0.0 + entropy = 0.0 + + from rmgpy.data.solvation import SolvationCorrection + return SolvationCorrection(enthalpy, gibbs, entropy) \ No newline at end of file diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 117c19644cf..c6fe6d8c04f 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(use_ml_solvation=True, name='solvation'): + """ + Enable the use of machine learning model for solvation correction + """ + + from rmgpy.data.solvation import MLSolvation + if use_ml_solvation: + model_path = os.path.join(settings['database.directory'], 'thermo', 'ml', name) # Need to add ML model at RMG-database/thermo/ml/solvation + if not os.path.exists(model_path): + raise InputError('Cannot find ML models folder {}'.format(model_path)) + rmg.ml_solvation = MLSolvation(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/thermo/thermoengine.py b/rmgpy/thermo/thermoengine.py index 715d5ea5698..e686e35cfc2 100644 --- a/rmgpy/thermo/thermoengine.py +++ b/rmgpy/thermo/thermoengine.py @@ -68,22 +68,10 @@ def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''): solvation_database = get_db('solvation') solute_data = solvation_database.get_solute_data(spc) - if False: # Need to add option variable "use_ml_solvation" to use ML solvation correction - if False: # Need to add option variable "use_ml_LSER_parameters" to use ML calculated parameters - - # Case 1. Use LSER with parameters from database to predict solvation correction - solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data) - else: - # Case 2. Use LSER with parameters from ML model to predict solvation correction - 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) - solvation_correction = solvation_database.get_solvation_correction_ml_LSER(solute_mol, solvent_mol) - else: - # Case 3. Use direct ML model to predict solvation correction + try: + from rmgpy.rmg.input import get_input + + ml_solvation = get_input("ml_solvation") solvent_species = solvation_database.get_solvent_structure(solvent_name)[0] solvent_smiles = solvent_species.smiles solvent_mol = Chem.MolFromSmiles(solvent_smiles) @@ -91,12 +79,24 @@ def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''): solute_smiles = spc.smiles solute_mol = Chem.MolFromSmiles(solute_smiles) - solvation_correction = solvation_database.get_solvation_correction_ml(solute_mol, solvent_mol) - - # 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}' + solvation_correction = ml_solvation.get_solvation_correction(solute_mol, solvent_mol) + + wilhoit.S0.value_si += solvation_correction.entropy + wilhoit.H0.value_si += solvation_correction.enthalpy + wilhoit.comment += ( + f" + Solvation correction (H={solvation_correction.enthalpy / 1e3:+.0f} kJ/mol; " + f"S={solvation_correction.entropy:+.0f} J/mol/K) using ML model" + ) + except Exception as e: + logging.warning("ML solvation correction not used: " + str(e)) + # fallback to LSER (default) + solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data) + wilhoit.S0.value_si += solvation_correction.entropy + wilhoit.H0.value_si += solvation_correction.enthalpy + wilhoit.comment += ( + f" + Solvation correction (H={solvation_correction.enthalpy / 1e3:+.0f} kJ/mol; " + f"S={solvation_correction.entropy:+.0f} J/mol/K) using LSER" + ) # Compute E0 by extrapolation to 0 K if spc.conformer is None: From 845c96e9ca7b802c5f58308ae6486bf42fe7a426 Mon Sep 17 00:00:00 2001 From: Bonhyeok Koo Date: Wed, 28 May 2025 17:31:28 -0400 Subject: [PATCH 3/6] Fix typo in the comment in input.py --- rmgpy/rmg/input.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index c6fe6d8c04f..a508d029c80 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -1312,7 +1312,7 @@ def ml_solvation(use_ml_solvation=True, name='solvation'): from rmgpy.data.solvation import MLSolvation if use_ml_solvation: - model_path = os.path.join(settings['database.directory'], 'thermo', 'ml', name) # Need to add ML model at RMG-database/thermo/ml/solvation + 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)) rmg.ml_solvation = MLSolvation(model_path) From 29a7c767a2441349ebe613917db8122100433185 Mon Sep 17 00:00:00 2001 From: Bonhyeok Koo Date: Wed, 4 Jun 2025 17:37:11 -0400 Subject: [PATCH 4/6] WIP: Apply T-dependent solvation correction --- rmgpy/data/solvation.py | 126 ++++++++++++++++++++++++++++++++++++---- rmgpy/reaction.py | 26 +++++++++ rmgpy/species.py | 48 +++++++++++++++ 3 files changed, 188 insertions(+), 12 deletions(-) diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index 7b149ccf1e9..56af39cb9c8 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -47,6 +47,7 @@ from rmgpy.data.thermo import is_aromatic_ring, is_bicyclic, find_aromatic_bonds_from_sub_molecule, \ convert_ring_to_sub_molecule, is_ring_partial_matched, bicyclic_decomposition_for_polyring, \ split_bicyclic_into_single_rings, saturate_ring_bonds +from rmgpy.data import get_db from rdkit import Chem @@ -2222,7 +2223,7 @@ 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. + object Note: This method utilizes the LSER method for solvation correction with parameters from the RMG-database. """ @@ -2433,6 +2434,20 @@ def get_Kfactor_parameters(self, delG298, delH298, delS298, solvent_name, T_tran return kfactor_parameters + def generate_solvation_model(self, species, solvent_name, T_dep=False, method="SoluteGC"): + """ + Returns a TDepModel or StaticModel depending on T_dep flag. + method: 'SoluteGC', 'SoluteML', or "DirectML' + """ + + # NEED TO IMPLEMENT: Extract T_dep and method option from the input file + if T_dep: + # A, B, C, D is calculated (Based on three methods) + return TDepModel(A, B, C, D, solvent_name) + else: + # dH298, dS298 is calculated (Based on three methods) + return StaticModel(dH298, dG298) + 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 @@ -2459,9 +2474,16 @@ class MLSolvation: """ A dummy class for machine learning-based solvation correction. """ - def __init__(self, model_path: str): + def __init__(self, model_path: str, T_dependent: bool = False, method: str = "SoluteGC"): + """ + model path: path to the directory where the machine learning model is stored # ex) model_path = "RMG-database/input/thermo/ml/solvation" + T_dependent: whether T_dependent solvation correction is used + method: the method used for solvation correction, e.g., "SoluteGC", "SoluteML" or "DirectML" + """ self.model_path = model_path + self.T_dependent = T_dependent + self.method = method self.model = self.load_model(model_path) def load_model(self, model_path): @@ -2471,15 +2493,95 @@ def load_model(self, model_path): def get_solvation_correction(self, solute_mol, solvent_mol): """ - Given a solute_mol and solvent_mol object, calculates the enthalpy, entropy, - and Gibbs free energy of solvation at 298 K using a machine learning model. - Returns a SolvationCorrection object. + Based on T_dependent flag and method, returns solvation correction + """ + if self.T_dependent: + solvation_database = get_db('solvation') + if self.method == "SoluteGC": + return solvation_database.get_T_dep_solvation_energy_from_soluteGC(solute_mol, solvent_mol, T) + elif self.method == "SoluteML": + return solvation_database.get_T_dep_solvation_energy_from_soluteML(solute_mol, solvent_mol, T, self.model) + elif self.method == "DirectML": + return solvation_database.get_T_dep_solvation_energy_from_directML(solute_mol, solvent_mol, T, self.model) + else: + return self.get_solvation_correction_at_298K(solute_mol, solvent_mol) + + +class TDepModel: + def __init__(self, A, B, C, D, solvent_name, T_trans_factor=0.75): + """ + solvent_name (str): name of the solvent that is used in CoolProp + A, B, C, D (float): coefficients for the K-factor piecewise function + T_trans_factor (float): fraction of Tc for transition between equations + """ + self.A = A + self.B = B + self.C = C + self.D = D + self.T_trans_factor = T_trans_factor + self.solvent_name = solvent_name + + def get_Kfactor(self, T): """ - # Need to implement: solute_mol, solvent_mol → feature vector → ML prediction - print(f"[NOTICE] Dummy ML model utilized") - enthalpy = 0.0 #self.model.predict([...])[0] - gibbs = 0.0 - entropy = 0.0 + 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 + 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: + 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_free_energy_of_solvation(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. + """ + solvent_name = self.solvent_name + 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 + +class StaticModel: + def __init__(self, dH298, dG298): + self.dH298 = dH298 + self.dG298 = dG298 - from rmgpy.data.solvation import SolvationCorrection - return SolvationCorrection(enthalpy, gibbs, entropy) \ No newline at end of file + def get_free_energy_of_solvation(self, T): + delG = self.dG298 + (T - 298) * self.dS298 + return delG \ No newline at end of file diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index b8ab3951c8c..0326f23155b 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -757,6 +757,11 @@ 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) + + if True: # NEED TO IMPLEMENT: if solvation correction is needed (if solvation module exists?) + ddGsolv = self.get_solvation_free_energy_of_reaction(T) + dGrxn += ddGsolv # add solvation free energy if available + 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 +826,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(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: + print(f"Error in reactant {reactant.name}: {e}") + raise + + for product in self.products: + try: + ddGsolv += product.get_free_energy_of_solvation(T) + except Exception as e: + print(f"Error in product {product.name}: {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/species.py b/rmgpy/species.py index c17bdd182a4..fd865ab9eec 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -605,6 +605,21 @@ 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 the specified + temperature `T` in K. + """ + Gsolv = 0.0 + #if self.has_solvation_thermo(): + # Gsolv = self.get_solvation_thermo_data().get_free_energy_of_solvation(T) + #else: + # raise Exception('Unable to calculate solvation free energy for species {0!r}: ' + # 'no solvationthermo data available.'.format(self.label)) + + Gsolv = self.get_solvation_thermo_data().get_free_energy_of_solvation(T) + return Gsolv def get_sum_of_states(self, e_list): """ @@ -821,6 +836,39 @@ def get_thermo_data(self, solvent_name=''): return self.thermo + def get_solvation_thermo_data(self, solvent_name=''): + """ + Returns a solvationthermo object (TDepModel or StaticModel) of the current Species object. + """ + if self.solvationthermo is not None: + return self.solvationthermo # If already exists, return it + + solvation_database = get_db('solvation') + if solvation_database is None: + raise ValueError("Solvation database is not loaded") + + self.solvationthermo = solvation_database.generate_solvation_model( + species = self, + solvent_name = solvent_name or self.solvent, + ) + + return self.solvationthermo + + #if False: # Need to connect with T-dep option in the input file + # # Need to calculate A~D (Based on three options) + # self.solvationthermo = TDepModel( + # A=1.1, B=2.2, C=3.3, D=4.4, solvent_name=solvent_name or self.solvent + # ) + #else: + # # Need to calculate dH298, dG298 (Based on three options) + # + # + # self.solvationthermo = StaticModel( + # dH298=-20000.0, dG298=-10000.0 / 298.15 + # ) + + #return self.solvationthermo + def generate_transport_data(self): """ Generate the transport_data parameters for the species. From 66af8893d562bcbb91afd763fa8140fea7f77288 Mon Sep 17 00:00:00 2001 From: Bonhyeok Koo Date: Tue, 24 Jun 2025 17:18:49 -0400 Subject: [PATCH 5/6] Apply solvation correction ddGsolv to dGrxn --- rmgpy/data/solvation.py | 233 ++++++++++++++++++++++++++++------- rmgpy/reaction.py | 18 ++- rmgpy/rmg/input.py | 12 +- rmgpy/rmg/main.py | 78 ++++++------ rmgpy/species.pxd | 1 + rmgpy/species.py | 56 +++------ rmgpy/thermo/thermoengine.py | 71 ++++++----- 7 files changed, 305 insertions(+), 164 deletions(-) diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index 56af39cb9c8..da78d6f280d 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -47,7 +47,6 @@ from rmgpy.data.thermo import is_aromatic_ring, is_bicyclic, find_aromatic_bonds_from_sub_molecule, \ convert_ring_to_sub_molecule, is_ring_partial_matched, bicyclic_decomposition_for_polyring, \ split_bicyclic_into_single_rings, saturate_ring_bonds -from rmgpy.data import get_db from rdkit import Chem @@ -464,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): """ @@ -2231,6 +2231,23 @@ def get_solvation_correction(self, solute_data, solvent_data): 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_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) + return correction def get_Kfactor(self, delG298, delH298, delS298, solvent_name, T): @@ -2434,20 +2451,6 @@ def get_Kfactor_parameters(self, delG298, delH298, delS298, solvent_name, T_tran return kfactor_parameters - def generate_solvation_model(self, species, solvent_name, T_dep=False, method="SoluteGC"): - """ - Returns a TDepModel or StaticModel depending on T_dep flag. - method: 'SoluteGC', 'SoluteML', or "DirectML' - """ - - # NEED TO IMPLEMENT: Extract T_dep and method option from the input file - if T_dep: - # A, B, C, D is calculated (Based on three methods) - return TDepModel(A, B, C, D, solvent_name) - else: - # dH298, dS298 is calculated (Based on three methods) - return StaticModel(dH298, dG298) - 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 @@ -2472,54 +2475,184 @@ def check_solvent_in_initial_species(self, rmg, solvent_structure): class MLSolvation: """ - A dummy class for machine learning-based solvation correction. + A class for solvation correction. """ - def __init__(self, model_path: str, T_dependent: bool = False, method: str = "SoluteGC"): + 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_dependent: whether T_dependent solvation correction is used + 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_dependent = T_dependent + self.T_dep = T_dep self.method = method - self.model = self.load_model(model_path) + self.model = self.load_model(model_path) if method != "SoluteGC" else None 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"))" - def get_solvation_correction(self, solute_mol, solvent_mol): + 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) + + 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): """ - Based on T_dependent flag and method, returns solvation correction + Returns solvation correction at 298K based on 3 methods: SoluteGC, SoluteML, or DirectML. """ - if self.T_dependent: - solvation_database = get_db('solvation') - if self.method == "SoluteGC": - return solvation_database.get_T_dep_solvation_energy_from_soluteGC(solute_mol, solvent_mol, T) - elif self.method == "SoluteML": - return solvation_database.get_T_dep_solvation_energy_from_soluteML(solute_mol, solvent_mol, T, self.model) + # 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_T_dep_solvation_energy_from_directML(solute_mol, solvent_mol, T, self.model) + return solvation_database.get_solvation_correction_from_directML(solute_mol, solvent_mol, self.model) else: - return self.get_solvation_correction_at_298K(solute_mol, solvent_mol) + 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. + + Given delG298 (J/mol), delH298 (J/mol), delS298 (J/mol/K), and solvent_name, + it finds the fitted K-factor parameters for the solvent-solute pair. + The parameters (A, B, C, D) are determined by enforcing the smooth continuity of the piecewise + functions at transition temperature and by making K-factor match in value and temperature + gradient at 298 K with those estimated from Abraham and Mintz LSERs. + The four equations are listed here: + 1) Match in value of K-factor with the estimate from the Abraham LSER at 298 K + A + B(1-Tr)^0.355 + Cexp(1-Tr)(Tr)^0.59 = Tr*ln(K-factor) + 2) Match in temperature gradient of K-factor with the estimate from the Mintz LSER at 298 K + -0.355B / Tc (1-Tr)^-0.645 + Cexp(1-Tr) / Tc * (0.59(Tr)^-0.41 - (Tr)^0.59) = d(Tr*ln(K-factor))/dT + 3) Continuity of the piecewise function at T_transition: + A + B(1-Tr)^0.355 + Cexp(1-Tr)(Tr)^0.59 = D(rho_l / rho_c - 1) + 4) Continuous temperature gradient of the piecewise function at T_transition + -0.355B / Tc (1-Tr)^-0.645 + Cexp(1-Tr) / Tc * (0.59(Tr)^-0.41 - (Tr)^0.59) = D / rho_c * d(rho_l)/dT + The conversion between dGsolv estimate from the Abraham and K-factor is shown below: + dGsolv = RTln(K-factor * rho_g / rho_l) + where rho_g is the saturated gas phase density of the solvent. + See the RMG documentation "Liquid Phase Systems" section on a temperature-dependent model for more details. + + 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_trans_factor (float): a temperature [K] for transitioning from the first piecewise function to the + 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 + rho_c = PropsSI('rhomolar_critical', solvent_name) # the critical density of the solvent, in mol/m^3 + + # Generate Amatrix and bvector for Ax = b + Amatrix = np.zeros((4, 4)) + bvec = np.zeros((4, 1)) + # 1. Tr*ln(K-factor) value at T = 298 K + rho_g_298 = get_gas_saturation_density(solvent_name, 298) + rho_l_298 = get_liquid_saturation_density(solvent_name, 298) + K298 = math.exp(delG298 / (298 * constants.R)) / rho_g_298 * rho_l_298 # K-factor + x298 = 298. / Tc * math.log(K298) # Tr*ln(K-factor), in K + Amatrix[0][0] = 1 + Amatrix[0][1] = (1 - 298 / Tc) ** 0.355 + Amatrix[0][2] = math.exp(1 - 298 / Tc) * (298 / Tc) ** 0.59 + Amatrix[0][3] = 0 + bvec[0] = x298 + # 2. d(Tr*ln(K-factor)) / dT at T = 298. Use finite difference method to get the temperature gradient from + # delG, delH, and delS at 298 K + T2 = 299 + delG_T2 = delH298 - delS298 * T2 + rho_g_T2 = get_gas_saturation_density(solvent_name, T2) + rho_l_T2 = get_liquid_saturation_density(solvent_name, T2) + K_T2 = math.exp(delG_T2 / (T2 * constants.R)) / rho_g_T2 * rho_l_T2 + x_T2 = T2 / Tc * math.log(K_T2) # Tln(K-factor) at 299 K, in K + slope298 = (x_T2 - x298) / (T2 - 298) + Amatrix[1][0] = 0 + Amatrix[1][1] = -0.355 / Tc * ((1 - 298 / Tc) ** (-0.645)) + Amatrix[1][2] = 1 / Tc * math.exp(1 - 298 / Tc) * (0.59 * (298 / Tc) ** (-0.41) - (298 / Tc) ** 0.59) + Amatrix[1][3] = 0 + bvec[1] = slope298 + # 3. Tln(K-factor) continuity at T = T_transition + rho_l_Ttran = get_liquid_saturation_density(solvent_name, T_transition) + Amatrix[2][0] = 1 + Amatrix[2][1] = (1 - T_transition / Tc) ** 0.355 + Amatrix[2][2] = math.exp(1 - T_transition / Tc) * (T_transition / Tc) ** 0.59 + Amatrix[2][3] = -(rho_l_Ttran - rho_c) / rho_c + bvec[2] = 0 + # 4. d(Tln(K-factor)) / dT smooth transition at T = T_transition + T3 = T_transition + 1 + rho_l_T3 = get_liquid_saturation_density(solvent_name, T3) + Amatrix[3][0] = 0 + Amatrix[3][1] = -0.355 / Tc * ((1 - T_transition / Tc) ** (-0.645)) + Amatrix[3][2] = 1 / Tc * math.exp(1 - T_transition / Tc) * ( + 0.59 * (T_transition / Tc) ** (-0.41) - (T_transition / Tc) ** 0.59) + Amatrix[3][3] = - ((rho_l_T3 - rho_l_Ttran) / rho_c / (T3 - T_transition)) + bvec[3] = 0 + # solve for the parameters + param, residues, ranks, s = np.linalg.lstsq(Amatrix, bvec, rcond=None) + # store the results in kfactor_parameters class + kfactor_parameters = KfactorParameters() + 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 class TDepModel: - def __init__(self, A, B, C, D, solvent_name, T_trans_factor=0.75): + def __init__(self, kfactor_parameters, solvent_name, T_trans_factor=0.75): """ - solvent_name (str): name of the solvent that is used in CoolProp + 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 """ - self.A = A - self.B = B - self.C = C - self.D = D + 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): """ @@ -2538,7 +2671,7 @@ def get_Kfactor(self, T): DatabaseError: if the given solvent_name is not available in CoolProp. """ - solvent_name = self.solvent_name + solvent_name = self.solvent_name_in_coolprop if solvent_name is not None: Tc = get_critical_temperature(solvent_name) if T < Tc: @@ -2557,8 +2690,7 @@ def get_Kfactor(self, T): f"is not available for the given solvent name: {solvent_name}") return Kfactor - - def get_free_energy_of_solvation(self, T): + 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. @@ -2567,21 +2699,30 @@ def get_free_energy_of_solvation(self, T): 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 + 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 + 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, dG298): + def __init__(self, dH298, dS298): self.dH298 = dH298 - self.dG298 = dG298 + self.dS298 = dS298 def get_free_energy_of_solvation(self, T): - delG = self.dG298 + (T - 298) * self.dS298 + 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 0326f23155b..767f16a293b 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -758,9 +758,15 @@ 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) - if True: # NEED TO IMPLEMENT: if solvation correction is needed (if solvation module exists?) - ddGsolv = self.get_solvation_free_energy_of_reaction(T) - dGrxn += ddGsolv # add solvation free energy if available + # 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 @@ -826,7 +832,7 @@ 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(T): + 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. """ @@ -835,14 +841,14 @@ def get_solvation_free_energy_of_reaction(T): try: ddGsolv -= reactant.get_free_energy_of_solvation(T) except Exception as e: - print(f"Error in reactant {reactant.name}: {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: - print(f"Error in product {product.name}: {e}") + logging.error("Problem with product {!r} in reaction {!s}: {}".format(product.label, self, e)) raise return ddGsolv diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index a508d029c80..50930673a00 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -1305,17 +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(use_ml_solvation=True, name='solvation'): - """ - Enable the use of machine learning model for solvation correction - """ +def ml_solvation(T_dep=False, method="SoluteGC", name='solvation'): from rmgpy.data.solvation import MLSolvation - if use_ml_solvation: + 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)) - rmg.ml_solvation = MLSolvation(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, diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index fb4ae33ef8b..892f6091990 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 0e30ce2b0e9..09911d927cd 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 fd865ab9eec..009ab356874 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): """ @@ -608,19 +616,10 @@ def get_free_energy(self, T): def get_free_energy_of_solvation(self, T): """ - Return the Gibbs free energy of solvation in J/mol for the species at the specified - temperature `T` in K. + Return the Gibbs free energy of solvation in J/mol for the species at temperature T [K]. """ - Gsolv = 0.0 - #if self.has_solvation_thermo(): - # Gsolv = self.get_solvation_thermo_data().get_free_energy_of_solvation(T) - #else: - # raise Exception('Unable to calculate solvation free energy for species {0!r}: ' - # 'no solvationthermo data available.'.format(self.label)) - - Gsolv = self.get_solvation_thermo_data().get_free_energy_of_solvation(T) - return Gsolv - + 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` @@ -836,38 +835,13 @@ def get_thermo_data(self, solvent_name=''): return self.thermo - def get_solvation_thermo_data(self, solvent_name=''): + def get_solvation_thermo_data(self): """ Returns a solvationthermo object (TDepModel or StaticModel) of the current Species object. """ - if self.solvationthermo is not None: - return self.solvationthermo # If already exists, return it - - solvation_database = get_db('solvation') - if solvation_database is None: - raise ValueError("Solvation database is not loaded") - - self.solvationthermo = solvation_database.generate_solvation_model( - species = self, - solvent_name = solvent_name or self.solvent, - ) - - return self.solvationthermo - - #if False: # Need to connect with T-dep option in the input file - # # Need to calculate A~D (Based on three options) - # self.solvationthermo = TDepModel( - # A=1.1, B=2.2, C=3.3, D=4.4, solvent_name=solvent_name or self.solvent - # ) - #else: - # # Need to calculate dH298, dG298 (Based on three options) - # - # - # self.solvationthermo = StaticModel( - # dH298=-20000.0, dG298=-10000.0 / 298.15 - # ) - - #return self.solvationthermo + if self.has_solvation_thermo(): + return self.solvationthermo + raise Exception("No solvationthermo available for species {}".format(self.label)) def generate_transport_data(self): """ diff --git a/rmgpy/thermo/thermoengine.py b/rmgpy/thermo/thermoengine.py index e686e35cfc2..aca51fce89e 100644 --- a/rmgpy/thermo/thermoengine.py +++ b/rmgpy/thermo/thermoengine.py @@ -68,35 +68,35 @@ def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''): solvation_database = get_db('solvation') solute_data = solvation_database.get_solute_data(spc) - try: - from rmgpy.rmg.input import get_input + # try: + # from rmgpy.rmg.input import get_input - ml_solvation = get_input("ml_solvation") - 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) - - solvation_correction = ml_solvation.get_solvation_correction(solute_mol, solvent_mol) - - wilhoit.S0.value_si += solvation_correction.entropy - wilhoit.H0.value_si += solvation_correction.enthalpy - wilhoit.comment += ( - f" + Solvation correction (H={solvation_correction.enthalpy / 1e3:+.0f} kJ/mol; " - f"S={solvation_correction.entropy:+.0f} J/mol/K) using ML model" - ) - except Exception as e: - logging.warning("ML solvation correction not used: " + str(e)) - # fallback to LSER (default) - solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data) - wilhoit.S0.value_si += solvation_correction.entropy - wilhoit.H0.value_si += solvation_correction.enthalpy - wilhoit.comment += ( - f" + Solvation correction (H={solvation_correction.enthalpy / 1e3:+.0f} kJ/mol; " - f"S={solvation_correction.entropy:+.0f} J/mol/K) using LSER" - ) + # ml_solvation = get_input("ml_solvation") + # 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) + + # solvation_correction = ml_solvation.get_solvation_correction(solute_mol, solvent_mol) + + # wilhoit.S0.value_si += solvation_correction.entropy + # wilhoit.H0.value_si += solvation_correction.enthalpy + # wilhoit.comment += ( + # f" + Solvation correction (H={solvation_correction.enthalpy / 1e3:+.0f} kJ/mol; " + # f"S={solvation_correction.entropy:+.0f} J/mol/K) using ML model" + # ) + # except Exception as e: + # logging.warning("ML solvation correction not used: " + str(e)) + # # fallback to LSER (default) + # solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data) + # wilhoit.S0.value_si += solvation_correction.entropy + # wilhoit.H0.value_si += solvation_correction.enthalpy + # wilhoit.comment += ( + # f" + Solvation correction (H={solvation_correction.enthalpy / 1e3:+.0f} kJ/mol; " + # f"S={solvation_correction.entropy:+.0f} J/mol/K) using LSER" + # ) # Compute E0 by extrapolation to 0 K if spc.conformer is None: @@ -208,3 +208,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 From 0ce7b7c20adef34cb3c5e5ebc1f9859e839c9605 Mon Sep 17 00:00:00 2001 From: Bonhyeok Koo Date: Wed, 25 Jun 2025 10:42:33 -0400 Subject: [PATCH 6/6] Remove redundant methods and update solvationTest --- rmgpy/data/solvation.py | 201 ------------------------------- rmgpy/thermo/thermoengine.py | 30 ----- test/rmgpy/data/solvationTest.py | 90 ++++++++------ 3 files changed, 54 insertions(+), 267 deletions(-) diff --git a/rmgpy/data/solvation.py b/rmgpy/data/solvation.py index da78d6f280d..4f6345d97bf 100644 --- a/rmgpy/data/solvation.py +++ b/rmgpy/data/solvation.py @@ -2250,207 +2250,6 @@ def get_solvation_correction_from_directML(self, solute_mol, solvent_mol, model) 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. - - """ - - 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)) - 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. - - """ - - 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. - - 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. - """ - - 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 - - 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. - - Given delG298 (J/mol), delH298 (J/mol), delS298 (J/mol/K), and solvent_name, - it finds the fitted K-factor parameters for the solvent-solute pair. - The parameters (A, B, C, D) are determined by enforcing the smooth continuity of the piecewise - functions at transition temperature and by making K-factor match in value and temperature - gradient at 298 K with those estimated from Abraham and Mintz LSERs. - The four equations are listed here: - 1) Match in value of K-factor with the estimate from the Abraham LSER at 298 K - A + B(1-Tr)^0.355 + Cexp(1-Tr)(Tr)^0.59 = Tr*ln(K-factor) - 2) Match in temperature gradient of K-factor with the estimate from the Mintz LSER at 298 K - -0.355B / Tc (1-Tr)^-0.645 + Cexp(1-Tr) / Tc * (0.59(Tr)^-0.41 - (Tr)^0.59) = d(Tr*ln(K-factor))/dT - 3) Continuity of the piecewise function at T_transition: - A + B(1-Tr)^0.355 + Cexp(1-Tr)(Tr)^0.59 = D(rho_l / rho_c - 1) - 4) Continuous temperature gradient of the piecewise function at T_transition - -0.355B / Tc (1-Tr)^-0.645 + Cexp(1-Tr) / Tc * (0.59(Tr)^-0.41 - (Tr)^0.59) = D / rho_c * d(rho_l)/dT - The conversion between dGsolv estimate from the Abraham and K-factor is shown below: - dGsolv = RTln(K-factor * rho_g / rho_l) - where rho_g is the saturated gas phase density of the solvent. - See the RMG documentation "Liquid Phase Systems" section on a temperature-dependent model for more details. - - 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_trans_factor (float): a temperature [K] for transitioning from the first piecewise function to the - second piecewise function. - - """ - - Tc = get_critical_temperature(solvent_name) - T_transition = Tc * T_trans_factor # T_trans_factor is empirically set to 0.75 by default - rho_c = PropsSI('rhomolar_critical', solvent_name) # the critical density of the solvent, in mol/m^3 - - # Generate Amatrix and bvector for Ax = b - Amatrix = np.zeros((4, 4)) - bvec = np.zeros((4, 1)) - # 1. Tr*ln(K-factor) value at T = 298 K - rho_g_298 = get_gas_saturation_density(solvent_name, 298) - rho_l_298 = get_liquid_saturation_density(solvent_name, 298) - K298 = math.exp(delG298 / (298 * constants.R)) / rho_g_298 * rho_l_298 # K-factor - x298 = 298. / Tc * math.log(K298) # Tr*ln(K-factor), in K - Amatrix[0][0] = 1 - Amatrix[0][1] = (1 - 298 / Tc) ** 0.355 - Amatrix[0][2] = math.exp(1 - 298 / Tc) * (298 / Tc) ** 0.59 - Amatrix[0][3] = 0 - bvec[0] = x298 - # 2. d(Tr*ln(K-factor)) / dT at T = 298. Use finite difference method to get the temperature gradient from - # delG, delH, and delS at 298 K - T2 = 299 - delG_T2 = delH298 - delS298 * T2 - rho_g_T2 = get_gas_saturation_density(solvent_name, T2) - rho_l_T2 = get_liquid_saturation_density(solvent_name, T2) - K_T2 = math.exp(delG_T2 / (T2 * constants.R)) / rho_g_T2 * rho_l_T2 - x_T2 = T2 / Tc * math.log(K_T2) # Tln(K-factor) at 299 K, in K - slope298 = (x_T2 - x298) / (T2 - 298) - Amatrix[1][0] = 0 - Amatrix[1][1] = -0.355 / Tc * ((1 - 298 / Tc) ** (-0.645)) - Amatrix[1][2] = 1 / Tc * math.exp(1 - 298 / Tc) * (0.59 * (298 / Tc) ** (-0.41) - (298 / Tc) ** 0.59) - Amatrix[1][3] = 0 - bvec[1] = slope298 - # 3. Tln(K-factor) continuity at T = T_transition - rho_l_Ttran = get_liquid_saturation_density(solvent_name, T_transition) - Amatrix[2][0] = 1 - Amatrix[2][1] = (1 - T_transition / Tc) ** 0.355 - Amatrix[2][2] = math.exp(1 - T_transition / Tc) * (T_transition / Tc) ** 0.59 - Amatrix[2][3] = -(rho_l_Ttran - rho_c) / rho_c - bvec[2] = 0 - # 4. d(Tln(K-factor)) / dT smooth transition at T = T_transition - T3 = T_transition + 1 - rho_l_T3 = get_liquid_saturation_density(solvent_name, T3) - Amatrix[3][0] = 0 - Amatrix[3][1] = -0.355 / Tc * ((1 - T_transition / Tc) ** (-0.645)) - Amatrix[3][2] = 1 / Tc * math.exp(1 - T_transition / Tc) * ( - 0.59 * (T_transition / Tc) ** (-0.41) - (T_transition / Tc) ** 0.59) - Amatrix[3][3] = - ((rho_l_T3 - rho_l_Ttran) / rho_c / (T3 - T_transition)) - bvec[3] = 0 - # solve for the parameters - param, residues, ranks, s = np.linalg.lstsq(Amatrix, bvec, rcond=None) - # store the results in kfactor_parameters class - kfactor_parameters = KfactorParameters() - 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 - - return kfactor_parameters - 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 diff --git a/rmgpy/thermo/thermoengine.py b/rmgpy/thermo/thermoengine.py index aca51fce89e..54507d8f14a 100644 --- a/rmgpy/thermo/thermoengine.py +++ b/rmgpy/thermo/thermoengine.py @@ -68,36 +68,6 @@ def process_thermo_data(spc, thermo0, thermo_class=NASA, solvent_name=''): solvation_database = get_db('solvation') solute_data = solvation_database.get_solute_data(spc) - # try: - # from rmgpy.rmg.input import get_input - - # ml_solvation = get_input("ml_solvation") - # 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) - - # solvation_correction = ml_solvation.get_solvation_correction(solute_mol, solvent_mol) - - # wilhoit.S0.value_si += solvation_correction.entropy - # wilhoit.H0.value_si += solvation_correction.enthalpy - # wilhoit.comment += ( - # f" + Solvation correction (H={solvation_correction.enthalpy / 1e3:+.0f} kJ/mol; " - # f"S={solvation_correction.entropy:+.0f} J/mol/K) using ML model" - # ) - # except Exception as e: - # logging.warning("ML solvation correction not used: " + str(e)) - # # fallback to LSER (default) - # solvation_correction = solvation_database.get_solvation_correction(solute_data, solvent_data) - # wilhoit.S0.value_si += solvation_correction.entropy - # wilhoit.H0.value_si += solvation_correction.enthalpy - # wilhoit.comment += ( - # f" + Solvation correction (H={solvation_correction.enthalpy / 1e3:+.0f} kJ/mol; " - # f"S={solvation_correction.entropy:+.0f} J/mol/K) using LSER" - # ) - # Compute E0 by extrapolation to 0 K if spc.conformer is None: spc.conformer = Conformer() diff --git a/test/rmgpy/data/solvationTest.py b/test/rmgpy/data/solvationTest.py index 72d0926d552..7055065b66f 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):