From 41f65babc223c2ff8ce4dc4fbcfc568ec31eb1cc Mon Sep 17 00:00:00 2001 From: ssun30 Date: Tue, 1 Jul 2025 15:55:40 -0400 Subject: [PATCH 1/2] [Solvation] Added excluded species and libraries to solvation --- rmgpy/rmg/input.py | 57 +++++++++++++++++++++++++++------------------- 1 file changed, 34 insertions(+), 23 deletions(-) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index c2bebbe985..6fc2956ec2 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -134,7 +134,7 @@ def catalyst_properties(bindingEnergies=None, metal_db.load(os.path.join(settings['database.directory'], 'surface')) if metal and (bindingEnergies or surfaceSiteDensity): - raise InputError("In catalyst_properties section you should only specify a 'metal' shortcut " + raise InputError("In catalyst_properties section you should only specify a 'metal' shortcut " "or the surfaceSiteDensity and bindingEnergies, but not both.") if metal: @@ -151,7 +151,7 @@ def catalyst_properties(bindingEnergies=None, logging.info("Using default binding energies, Pt(111)") else: rmg.binding_energies = convert_binding_energies(bindingEnergies) - + if surfaceSiteDensity is None: rmg.surface_site_density = metal_db.get_surface_site_density("Pt111") @@ -171,8 +171,8 @@ def catalyst_properties(bindingEnergies=None, def convert_binding_energies(binding_energies): """ Process binding_energies dictionary from the input file - - It converts the values into Energy quantities, and checks that + + It converts the values into Energy quantities, and checks that all elements C,H,O, and N are present. :param binding_energies: a dictionary of element symbol: binding energy pairs @@ -406,7 +406,7 @@ def simple_reactor(temperature, logging.debug(" {0}".format(const_spc)) if const_spc not in species_dict: raise InputError('Species {0} not found in the input file'.format(const_spc)) - + if not isinstance(T, list): sensitivityTemperature = T if not isinstance(P, list): @@ -512,9 +512,9 @@ def constant_V_ideal_gas_reactor(temperature, termination.append(TerminationRateRatio(terminationRateRatio)) if len(termination) == 0: raise InputError('No termination conditions specified for reaction system #{0}.'.format(len(rmg.reaction_systems) + 2)) - + initial_cond = initialMoleFractions - initial_cond["T"] = T + initial_cond["T"] = T initial_cond["P"] = P system = ConstantVIdealGasReactor(rmg.reaction_model.core.phase_system,rmg.reaction_model.edge.phase_system,initial_cond,termination) system.T = Quantity(T) @@ -805,7 +805,7 @@ def constant_T_V_liquid_reactor(temperature, if outlet_volumetric_flow_rate: if outlet_volumetric_flow_rate != outletVolumetricFlowRate: raise InputError(' Inconsistent residence time and inlet volumetric flow rate') - + outlet_volumetric_flow_rate = Quantity(outletVolumetricFlowRate).value_si if inletConcentrations: @@ -907,7 +907,7 @@ def constant_T_V_liquid_reactor(temperature, evap_cond_conditions[key] = item evap_cond_conditions["P"] = vapor_pressure evap_cond_conditions["T"] = initial_conditions["T"] - + system = ConstantTVLiquidReactor(rmg.reaction_model.core.phase_system, rmg.reaction_model.edge.phase_system, initial_conditions, @@ -1162,19 +1162,30 @@ def simulator(atol, rtol, sens_atol=1e-6, sens_rtol=1e-4): rmg.simulator_settings_list.append(SimulatorSettings(atol, rtol, sens_atol, sens_rtol)) -def solvation(solvent,solventData=None): +def solvation(solvent, solventData=None, excludedSpecies=None, excludedLibraries=None): # If solvation module in input file, set the RMG solvent variable - #either a string corresponding to the solvent database or a olvent object - if isinstance(solvent, str): - rmg.solvent = solvent - else: - raise InputError("Solvent not specified properly, solvent must be string") + # either a string corresponding to the solvent database or a solvent object + # The user can define a list of species labels or thermo library names + # to skip solvation correction. This will be useful when explicit solvation + # effects of a library species are already accounted for. + if not isinstance(solvent, str): + raise InputError("Solvent must be a string corresponding to a solvent in the database.") - if isinstance(solventData, SolventData) or solventData is None: - rmg.solvent_data = solventData - else: - raise InputError("Solvent not specified properly, solventData must be None or SolventData object") + if not (isinstance(solventData, SolventData) or solventData is None): + raise InputError("solventData must be a SolventData object or None.") + + if excludedSpecies is not None: + if not (isinstance(excludedSpecies, list) and all(isinstance(s, str) for s in excludedSpecies)): + raise InputError("excludedSpecies must be a list of strings or None.") + + if excludedLibraries is not None: + if not (isinstance(excludedLibraries, list) and all(isinstance(l, str) for l in excludedLibraries)): + raise InputError("excludedLibraries must be a list of strings or None.") + rmg.solvent = solvent + rmg.solvent_data = solventData + rmg.solvation_excluded_species = excludedSpecies if excludedSpecies is not None else [] + rmg.solvation_excluded_libraries = excludedLibraries if excludedLibraries is not None else [] def model(toleranceMoveToCore=None, toleranceRadMoveToCore=np.inf, toleranceMoveEdgeReactionToCore=np.inf, toleranceKeepInEdge=0.0, @@ -1191,7 +1202,7 @@ def model(toleranceMoveToCore=None, toleranceRadMoveToCore=np.inf, toleranceTransitoryDict={}, transitoryStepPeriod=20, toleranceReactionToCoreDeadendRadical=0.0): """ - How to generate the model. `toleranceMoveToCore` must be specified. + How to generate the model. `toleranceMoveToCore` must be specified. toleranceMoveReactionToCore and toleranceReactionInterruptSimulation refers to an additional criterion for forcing an edge reaction to be included in the core by default this criterion is turned off Other parameters are optional and control the pruning. @@ -1521,7 +1532,7 @@ def set_global_rmg(rmg0): def read_input_file(path, rmg0): """ - Read an RMG input file at `path` on disk into the :class:`RMG` object + Read an RMG input file at `path` on disk into the :class:`RMG` object `rmg`. """ global rmg, species_dict, mol_to_frag @@ -1620,7 +1631,7 @@ def read_input_file(path, rmg0): def read_thermo_input_file(path, rmg0): """ - Read an thermo estimation input file at `path` on disk into the :class:`RMG` object + Read an thermo estimation input file at `path` on disk into the :class:`RMG` object `rmg`. """ @@ -1679,7 +1690,7 @@ def read_thermo_input_file(path, rmg0): def save_input_file(path, rmg): """ - Save an RMG input file at `path` on disk from the :class:`RMG` object + Save an RMG input file at `path` on disk from the :class:`RMG` object `rmg`. """ From 8852ecf629f0132c05fe609b47a6b099c55325ae Mon Sep 17 00:00:00 2001 From: ssun30 Date: Mon, 7 Jul 2025 20:10:50 -0400 Subject: [PATCH 2/2] [Solvation] Added logging info when reading solvation inputs --- rmgpy/rmg/input.py | 24 +++++++++++++++++++++++- rmgpy/rmg/model.py | 6 ++++-- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 6fc2956ec2..1f71eb6239 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -1187,6 +1187,19 @@ def solvation(solvent, solventData=None, excludedSpecies=None, excludedLibraries rmg.solvation_excluded_species = excludedSpecies if excludedSpecies is not None else [] rmg.solvation_excluded_libraries = excludedLibraries if excludedLibraries is not None else [] + # Log species and libraries excluded from solvation corrections + if rmg.solvation_excluded_species: + logging.info( + f"Solvation corrections will be skipped for the following species: " + f"{', '.join(rmg.solvation_excluded_species)}" + ) + if rmg.solvation_excluded_libraries: + logging.info( + f"Solvation corrections will be skipped for species from the following thermo libraries: " + f"{', '.join(rmg.solvation_excluded_libraries)}" + ) + + def model(toleranceMoveToCore=None, toleranceRadMoveToCore=np.inf, toleranceMoveEdgeReactionToCore=np.inf, toleranceKeepInEdge=0.0, toleranceInterruptSimulation=1.0, @@ -1769,7 +1782,16 @@ def save_input_file(path, rmg): f.write(')\n\n') if rmg.solvent: - f.write("solvation(\n solvent = '{0!s}'\n)\n\n".format(rmg.solvent)) + f.write("solvation(\n") + f.write(" solvent = '{0!s}',\n".format(rmg.solvent)) + + if rmg.solvation_excluded_species: + f.write(" excludedSpecies = {0},\n".format(repr(rmg.solvation_excluded_species))) + + if rmg.solvation_excluded_libraries: + f.write(" excludedLibraries = {0},\n".format(repr(rmg.solvation_excluded_libraries))) + + f.write(")\n\n") # Simulator tolerances f.write('simulator(\n') diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index d677755b07..19ae9e35d4 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -237,6 +237,8 @@ def __init__(self, core=None, edge=None, surface=None): self.new_surface_spcs_loss = set() self.new_surface_rxns_loss = set() self.solvent_name = "" + self.solvation_excluded_species = [] + self.solvation_excluded_libraries = [] self.surface_site_density = None self.unrealgroups = [ Group().from_adjacency_list( @@ -557,7 +559,7 @@ def make_new_reaction(self, forward, check_existing=True, generate_thermo=True, elif isinstance(forward, LibraryReaction) and forward.is_surface_reaction(): # do fix the library reaction barrier if this is scaled from another metal if any(['Binding energy corrected by LSR' in x.thermo.comment for x in forward.reactants + forward.products]): - forward.fix_barrier_height(solvent=self.solvent_name) + forward.fix_barrier_height(solvent=self.solvent_name) elif forward.kinetics.solute: forward.apply_solvent_correction(solvent=self.solvent_name) if self.pressure_dependence and forward.is_unimolecular(): @@ -1601,7 +1603,7 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False, requires_rms=F self.new_reaction_list = [] self.new_species_list = [] edge_species_to_move = [] - + num_old_core_species = len(self.core.species) num_old_core_reactions = len(self.core.reactions)