adding potential and ElectrodeReactor to main.py and model.py

Author: davidfarinajr

rmgpy/rmg/main.py

@@ -73,6 +73,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
@@ -166,6 +167,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):
@@ -683,6 +686,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
 
@@ -696,6 +701,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:
@@ -709,7 +724,7 @@ def execute(self, initialize=True, **kwargs):
                     atol=self.simulator_settings_list[0].atol,
                     rtol=self.simulator_settings_list[0].rtol,
                     filter_reactions=True,
-                    conditions=self.rmg_memories[index].get_cond(),
+                    conditions=conditions,
                 )
 
                 self.update_reaction_threshold_and_react_flags(
@@ -730,7 +745,9 @@ def execute(self, initialize=True, **kwargs):
             self.reaction_model.enlarge(react_edge=True,
                                         unimolecular_react=self.unimolecular_react,
                                         bimolecular_react=self.bimolecular_react,
-                                        trimolecular_react=self.trimolecular_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(
@@ -791,6 +808,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
@@ -820,7 +846,7 @@ def execute(self, initialize=True, **kwargs):
                                 prune=prune,
                                 model_settings=model_settings,
                                 simulator_settings=simulator_settings,
-                                conditions=self.rmg_memories[index].get_cond()
+                                conditions=conditions
                             )
                         except:
                             logging.error("Model core reactions:")
@@ -833,9 +859,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)
@@ -920,7 +946,9 @@ def execute(self, initialize=True, **kwargs):
                         self.reaction_model.enlarge(react_edge=True,
                                                     unimolecular_react=self.unimolecular_react,
                                                     bimolecular_react=self.bimolecular_react,
-                                                    trimolecular_react=self.trimolecular_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):
@@ -2210,6 +2238,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)
@@ -2324,6 +2355,8 @@ def obj(y):
                        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():
@@ -2353,6 +2386,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)
 

rmgpy/rmg/model.py

@@ -47,7 +47,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
@@ -501,7 +501,8 @@ 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):
+                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
@@ -630,9 +631,13 @@ def enlarge(self, new_object=None, react_edge=False,
         # 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):
+            if isinstance(reaction.kinetics, KineticsData) and not isinstance(reaction.kinetics, SurfaceChargeTransfer):
                 reaction.kinetics = reaction.kinetics.to_arrhenius()
-            #  correct barrier heights of estimated kinetics
+
+            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.