OpenMKM-IO: Thermo YAML Creation bug

Version of pMuTT 1.2.21

Describe the bug List of issues when making a YAML thermo file using the write_thermo_yaml method.

a) Change the default thermo model keyword for interfaces/surfaces with or without lateral interaction to use "surface-lateral-interaction" instead of the "ideal-surface" model which has issues in OpenMKM. https://github.com/VlachosGroup/pMuTT/blob/f57936ed8ad9ec1770f824f469bcea725cdbd4dd/pmutt/omkm/phase.py#L333 Resolution : Use yaml_dict['thermo'] = 'surface-lateral-interaction'

b) Use the latest thermo model for bulk phases class StoichSolid developed for OpenMKM. https://github.com/VlachosGroup/pMuTT/blob/f57936ed8ad9ec1770f824f469bcea725cdbd4dd/pmutt/omkm/phase.py#L101 Resolution : Use 'thermo': 'ref-state-fixed-stoichiometry'

c) write_thermo_yaml doesn't accept the option to activate or de-activate motz-wize corrections for adsorption/desorption reactions based on Sticking Coeeficient. https://github.com/VlachosGroup/pMuTT/blob/f57936ed8ad9ec1770f824f469bcea725cdbd4dd/pmutt/io/omkm.py#L173 Possible Resolution : Add the following lines

--- "a/.\\pMuTT\\pmutt\\io\\omkm.py"
+++ "b/.\\pMuTT-modified\\io\\omkm.py"
@@ -174,6 +174,7 @@ def write_thermo_yaml(phases=None, species=None, reactions=None,
                       lateral_interactions=None, units=None,
                       filename=None, T=300., P=1., newline='\n',
                       ads_act_method='get_H_act',
+                      use_motz_wise='False',
                       yaml_options={'default_flow_style': None, 'indent': 2,
                                     'sort_keys': False, 'width': 79}):
     """Writes the units, phases, species, lateral interactions, reactions and
@@ -246,6 +247,7 @@ def write_thermo_yaml(phases=None, species=None, reactions=None,
             if reaction.id is None:
                 reaction.id = 'r_{:04d}'.format(i)
                 i += 1
+            reaction.use_motz_wise = use_motz_wise
             # Write reaction
             reaction_dict = reaction.to_omkm_yaml(units=units, T=T)
             reactions_out.append(reaction_dict)

d) Adsorbate-Adsorbate lateral interaction units are not written correctly in the YAML file. Currently, pMuTT writes "kcal" instead of "kcal/mol" Related file: https://github.com/VlachosGroup/pMuTT/blob/f57936ed8ad9ec1770f824f469bcea725cdbd4dd/pmutt/mixture/cov.py#L207 Possible Resolution

--- "a/.\\pMuTT\\pmutt\\mixture\\cov.py"
+++ "b/.\\pMuTT-modified\\mixture\\cov.py"
@@ -204,13 +204,15 @@ class PiecewiseCovEffect(_ModelBase):
                                    self.intervals, slopes, self.name))
         return lat_inter_str

-    def to_omkm_yaml(self, energy_unit=None, units=None):
+    def to_omkm_yaml(self, energy_unit='kcal', quantity_unit='mol', units=None):
         """Writes the object in Cantera's YAML format.

         Parameters
         ----------
             energy_unit : str, optional
-                Unit to use for energy. Default is 'cal/mol'
+                Energy unit for slopes. Default is 'kcal'
+            quantity_unit : str, optional
+                Quantity unit for slopes. Default is 'mol'
             units : :class:`~pmutt.omkm.units.Units` object
                 If specified, `energy_unit` is overwritten. Default is None.
         Returns
@@ -218,11 +220,15 @@ class PiecewiseCovEffect(_ModelBase):
             yaml_dict : dict
                 Dictionary compatible with Cantera's YAML format
         """
+        if units is not None:
+            energy_unit = units.energy
+            quantity_unit = units.quantity
+        final = '{}/{}'.format(energy_unit, quantity_unit)
         yaml_dict = {}
         yaml_dict['species'] = [self.name_i, self.name_j]
         yaml_dict['coverage-threshold'] = self.intervals
         '''Assign slope'''
-        strength_param = _Param('strength', self.slopes, '_energy')
+        strength_param = _Param('strength', self.slopes, final)
         _assign_yaml_val(strength_param, yaml_dict, units)
         yaml_dict['id'] = self.name
         return yaml_dict

e) pMuTT doesn't consider explicitly defined E_a value in the spreadsheet, and continues to calculate Ea from either the TS species (if defined) or from the heat of reaction. Resolution: Need to add the following checks for explicitly defined Ea.

--- "a/.\\pMuTT\\pmutt\\omkm\\reaction.py"
+++ "b/.\\pMuTT-modified\\omkm\\reaction.py"
@@ -530,6 +530,13 @@ class SurfaceReaction(Reaction):
             length_unit = units.length
             act_energy_unit = units.act_energy

+        if self.Ea is None:
+            act_val = None
+        else:
+            act_val = c.convert_unit(self.Ea,
+                                     initial='kcal/mol',
+                                     final=act_energy_unit)
+
         yaml_dict = {}
         # Assign reaction name
         yaml_dict['equation'] = self.to_string(stoich_space=True,
@@ -543,8 +550,11 @@ class SurfaceReaction(Reaction):
             A_param = _Param('A', self.sticking_coeff, None)

             # Activation energy
-            act_method = getattr(self, ads_act_method)
-            act_val = act_method(units=act_energy_unit, T=T, P=P)
+            # If activation energy not specified, calculate using requested
+            # method of activation
+            if act_val is None:
+                act_method = getattr(self, ads_act_method)
+                act_val = act_method(units=act_energy_unit, T=T, P=P)

             # Sticking-species
             for species in self.reactants:
@@ -570,7 +580,8 @@ class SurfaceReaction(Reaction):
             # TODO Check units for A. Should have time dependence.
             A_param = _Param('A', A, None)
             # Activation energy
-            act_val = self.get_G_act(units=act_energy_unit, T=T, P=P)
+            if act_val is None:
+                act_val = self.get_G_act(units=act_energy_unit, T=T, P=P)

         # Assign activation energy, beta and pre-exponential factor
         rate_constant_dict = {}

Additional Context for E_a Care must be taken to check all possible scenarios of writing the activation energy when reading from an Excel file. Things like

What takes precedence? Explciit Ea or TS species or Heat or Reaction.
What happens when multiple possible methods of E_a are possible.
Make sure that the correct method is used to calculate the value i.e. act_method and the ads_act_method.

VlachosGroup / pMuTT

OpenMKM-IO: Thermo YAML Creation bug #199