Skip to content

Commit

Permalink
adding potential and ElectrodeReactor to main.py and model.py
Browse files Browse the repository at this point in the history
  • Loading branch information
@davidfarinajr @ssun30
davidfarinajr authored and ssun30 committed Mar 29, 2024
1 parent 08f0f33 commit 6e04e13
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 15 deletions.
57 changes: 44 additions & 13 deletions rmgpy/rmg/main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@
from rmgpy.rmg.settings import ModelSettings
from rmgpy.solver.base import TerminationTime, TerminationConversion
from rmgpy.solver.simple import SimpleReactor
from rmgpy.solver.electrode import ElectrodeReactor
from rmgpy.stats import ExecutionStatsWriter
from rmgpy.thermo.thermoengine import submit
from rmgpy.tools.plot import plot_sensitivity
Expand Down Expand Up @@ -170,6 +171,8 @@ def __init__(self, input_file=None, output_directory=None, profiler=None, stats_
self.Tmax = 0.0
self.Pmin = 0.0
self.Pmax = 0.0
self.potential_min = 0.0
self.potential_max = 0.0
self.database = None

def clear(self):
Expand Down Expand Up @@ -762,6 +765,8 @@ def execute(self, initialize=True, **kwargs):
try:
self.Pmin = min([x.Prange[0].value_si if hasattr(x, "Prange") and x.Prange else x.P.value_si for x in self.reaction_systems])
self.Pmax = max([x.Prange[1].value_si if hasattr(x, "Prange") and x.Prange else x.P.value_si for x in self.reaction_systems])
self.potential_min = min([x.potential_range[0].value_si if x.potential_range else x.potential.value_si for x in self.reaction_systems])
self.potential_max = max([x.potential_range[1].value_si if x.potential_range else x.potential.value_si for x in self.reaction_systems])
except AttributeError:
pass

Expand All @@ -775,6 +780,16 @@ def execute(self, initialize=True, **kwargs):
self.rmg_memories.append(RMG_Memory(reaction_system, self.balance_species))
self.rmg_memories[index].generate_cond()
log_conditions(self.rmg_memories, index)
conditions = self.rmg_memories[index].get_cond()
temperature = potential = None
if isinstance(reaction_system, ElectrodeReactor):
if conditions:
potential = conditions.get('potential')
temperature = conditions.get('T')
if not potential:
potential = reaction_system.potential.value_si
if not temperature:
temperature = reaction_system.T.value_si

# Update react flags
if self.filter_reactions:
Expand Down Expand Up @@ -817,12 +832,12 @@ def execute(self, initialize=True, **kwargs):
logging.info("Generating initial reactions...")

# React core species to enlarge edge
self.reaction_model.enlarge(
react_edge=True,
unimolecular_react=self.unimolecular_react,
bimolecular_react=self.bimolecular_react,
trimolecular_react=self.trimolecular_react,
)
self.reaction_model.enlarge(react_edge=True,
unimolecular_react=self.unimolecular_react,
bimolecular_react=self.bimolecular_react,
trimolecular_react=self.trimolecular_react,
temperature=temperature,
potential=potential)

if not np.isinf(self.model_settings_list[0].thermo_tol_keep_spc_in_edge):
self.reaction_model.set_thermodynamic_filtering_parameters(
Expand Down Expand Up @@ -882,6 +897,15 @@ def execute(self, initialize=True, **kwargs):
objects_to_enlarge = []

conditions = self.rmg_memories[index].get_cond()
if isinstance(reaction_system, ElectrodeReactor):
if conditions:
potential = conditions.get('potential')
temperature = conditions.get('T')
if not potential:
potential = reaction_system.potential.value_si
if not temperature:
temperature = reaction_system.T.value_si

if conditions and self.solvent:
T = conditions["T"]
# Set solvent viscosity
Expand Down Expand Up @@ -970,9 +994,9 @@ def execute(self, initialize=True, **kwargs):
self.make_seed_mech() # Just in case the user wants to restart from this
raise

log_conditions(self.rmg_memories, index)
self.rmg_memories[index].add_t_conv_N(t, x, len(obj))
self.rmg_memories[index].generate_cond()
log_conditions(self.rmg_memories, index)

reactor_done = self.reaction_model.add_new_surface_objects(obj, new_surface_species, new_surface_reactions, reaction_system)

Expand Down Expand Up @@ -1089,12 +1113,12 @@ def execute(self, initialize=True, **kwargs):

old_edge_size = len(self.reaction_model.edge.reactions)
old_core_size = len(self.reaction_model.core.reactions)
self.reaction_model.enlarge(
react_edge=True,
unimolecular_react=self.unimolecular_react,
bimolecular_react=self.bimolecular_react,
trimolecular_react=self.trimolecular_react,
)
self.reaction_model.enlarge(react_edge=True,
unimolecular_react=self.unimolecular_react,
bimolecular_react=self.bimolecular_react,
trimolecular_react=self.trimolecular_react,
potential=potential,
temperature=temperature)

if old_edge_size != len(self.reaction_model.edge.reactions) or old_core_size != len(self.reaction_model.core.reactions):
reactor_done = False
Expand Down Expand Up @@ -2204,6 +2228,9 @@ def __init__(self, reaction_system, bspc):
if hasattr(reaction_system, "Prange") and isinstance(reaction_system.Prange, list):
Prange = reaction_system.Prange
self.Ranges["P"] = [np.log(P.value_si) for P in Prange]
if hasattr(reaction_system, 'potential_range') and isinstance(reaction_system.potential_range, list):
potential_range = reaction_system.potential_range
self.Ranges['potential'] = [V.value_si for V in potential_range]
if hasattr(reaction_system, "initial_mole_fractions"):
if bspc:
self.initial_mole_fractions = deepcopy(reaction_system.initial_mole_fractions)
Expand Down Expand Up @@ -2320,6 +2347,8 @@ def obj(y):
new_cond = {yk: yf[i] * (self.Ranges[yk][1] - self.Ranges[yk][0]) + self.Ranges[yk][0] for i, yk in enumerate(ykey)}
if "P" in list(new_cond.keys()):
new_cond["P"] = np.exp(new_cond["P"])
if 'potential' in list(new_cond.keys()):
new_cond['potential'] = new_cond['potential']

if hasattr(self, "initial_mole_fractions"):
for key in self.initial_mole_fractions.keys():
Expand Down Expand Up @@ -2349,6 +2378,8 @@ def log_conditions(rmg_memories, index):
s += "T = {0} K, ".format(item)
elif key == "P":
s += "P = {0} bar, ".format(item / 1.0e5)
elif key == 'potential':
s += 'potential = {0} V, '.format(item)
else:
s += key.label + " = {0}, ".format(item)

Expand Down
35 changes: 33 additions & 2 deletions rmgpy/rmg/model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@
from rmgpy.data.rmg import get_db
from rmgpy.display import display
from rmgpy.exceptions import ForbiddenStructureException
from rmgpy.kinetics import KineticsData, Arrhenius
from rmgpy.kinetics import KineticsData, Arrhenius, SurfaceChargeTransfer
from rmgpy.quantity import Quantity
from rmgpy.reaction import Reaction
from rmgpy.rmg.pdep import PDepReaction, PDepNetwork
Expand Down Expand Up @@ -589,7 +589,9 @@ def make_new_pdep_reaction(self, forward):

return forward

def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bimolecular_react=None, trimolecular_react=None):
def enlarge(self, new_object=None, react_edge=False,
unimolecular_react=None, bimolecular_react=None, trimolecular_react=None,
temperature=None, potential=None):
"""
Enlarge a reaction model by processing the objects in the list `new_object`.
If `new_object` is a
Expand Down Expand Up @@ -704,6 +706,35 @@ def enlarge(self, new_object=None, react_edge=False, unimolecular_react=None, bi
if not np.isinf(self.thermo_tol_keep_spc_in_edge) and self.new_species_list != []:
self.thermo_filter_species(self.new_species_list)

# Generate kinetics of new reactions
if self.new_reaction_list:
logging.info('Generating kinetics for new reactions...')
for reaction in self.new_reaction_list:
# If the reaction already has kinetics (e.g. from a library),
# assume the kinetics are satisfactory
if reaction.kinetics is None:
self.apply_kinetics_to_reaction(reaction)

# For new reactions, convert ArrheniusEP to Arrhenius, and fix barrier heights.
# self.new_reaction_list only contains *actually* new reactions, all in the forward direction.
for reaction in self.new_reaction_list:
# convert KineticsData to Arrhenius forms
if isinstance(reaction.kinetics, KineticsData) and not isinstance(reaction.kinetics, SurfaceChargeTransfer):
reaction.kinetics = reaction.kinetics.to_arrhenius()

if isinstance(reaction.kinetics, SurfaceChargeTransfer):
reaction.set_reference_potential(temperature)

# correct barrier heights of estimated kinetics
if isinstance(reaction, TemplateReaction) or isinstance(reaction,
DepositoryReaction): # i.e. not LibraryReaction
reaction.fix_barrier_height() # also converts ArrheniusEP to Arrhenius.

if self.pressure_dependence and reaction.is_unimolecular():
# If this is going to be run through pressure dependence code,
# we need to make sure the barrier is positive.
reaction.fix_barrier_height(force_positive=True)

# Update unimolecular (pressure dependent) reaction networks
if self.pressure_dependence:
# Recalculate k(T,P) values for modified networks
Expand Down

0 comments on commit 6e04e13

Please sign in to comment.