Merge pull request #316 from lpratalimaffei/dev

Dev

mechanalyzer/builder/creckclass.py

@@ -0,0 +1,74 @@
+""" functions for managing creck classes and dictionaries
+"""
+
+import pandas as pd
+
+R = ['R', 'H', 'OH', 'O2', 'O', 'CH3', 'HO2', 'HCO', 'C2H3']
+RSR = ['C5H5', 'C7H7', 'C6H5CH2', 'C6H5O', 'A1O']
+
+
+def build_creckclass_fromdct(rxn_creckclass_dct, classtype_dct = {}):
+    """ build a dataframe with reactiontype, classtype, speciestype, bimoltype for each reaction
+
+    Args:
+        rxn_creckclass_dct (dict{rxn: {'speciestype':speciestype, 
+        'reactiontype':reactiontype}}): assigned species and reaction type
+        classtype_dct (dict{reactiontype(str): classtype(str)}, optional): _description_. Defaults to None.
+    
+    return creckclass_df: dataframe[['classtype','speciestype','reactiontype','bimoltype'][rxn]]
+    """
+
+    creckclass_df = pd.DataFrame.from_dict(rxn_creckclass_dct, orient='index')
+    # add info on bimolecular type
+    for speciestype, df_sptype in creckclass_df.groupby('speciestype'):
+        for idx, reactiontype in df_sptype['reactiontype'].items():
+            bimoltype = set_bimol_type(idx, speciestype, reactiontype)
+            creckclass_df.loc[idx, ['bimoltype']] = bimoltype
+    for reactiontype, df_rxntype in creckclass_df.groupby('reactiontype'):
+        rxns = df_rxntype.index
+        if reactiontype in classtype_dct.keys():
+            cltype = classtype_dct[reactiontype]
+        else:
+            cltype = 'UNSORTED'
+        creckclass_df.loc[rxns, 'classtype'] = cltype
+
+    return creckclass_df
+
+def set_bimol_type(rxn, speciestype, reactiontype):
+    """
+    check if a bimolecular reaction is of type M+M, RSR+M, R+M, R+R, RSR+RSR
+    from speciestype A-M/R/RSR and reactiontype
+    NB for now, relies on rudimental classification. can be upgraded
+    by analysing radical character in autochem.
+    suggestion: check if molecule or radical, then for rsr check presence of multiple rad structures
+    risk: for aromatics, resonance will always be detected
+    """
+    # return unimol for unimolecular reactions
+    if len(rxn[0]) == 1:
+        return 'UNIMOL'
+    # CHECK SPECIES TYPE 1 - CLASSIFIED ACCORDING TO THE SPECIES TYPE FROM THE USER
+    try:
+        species_type_1 = speciestype.split('-')[-1]
+    except AttributeError:
+        return 'UNSORTED'
+
+    # CHECK SPECIES TYPE 2
+    try:
+        species_type_2 = reactiontype.split('_')[1]
+    except AttributeError:
+        return 'UNSORTED'
+
+    if species_type_2 in R:
+        species_type_2 = 'R'
+    elif species_type_2 in RSR:
+        species_type_2 = 'RSR'
+    elif '-' in species_type_2:
+        species_type_2 = species_type_2.split('-')[1]
+    else:
+        species_type_2 = 'NA'
+
+    type_list = [species_type_1, species_type_2]
+    type_list.sort()
+
+    return '+'.join(type_list)
+

mechanalyzer/builder/sort_fct.py

@@ -322,7 +322,12 @@ def sort(self, hierarchy, species_list):
                       'submech', 'submech_prompt', 'submech_keepsubfuel',
                       'submech_deletelarge', 'maxval', 'maxratio',
                       ]
-
+        # check that all criteria in hierarchy are allowed
+        notavail = [hr for hr in hierarchy[:-1] if hr not in criteria_all]
+        assert len(notavail) == 0, (
+            '*Error: classification according to {} not available; criteria are: \n {}'.format(notavail,
+                                                                                               '\n'.join(criteria_all))
+        )
         # check on pes/subpes criteria:
         # if subpes alone, also add pes
         if (('subpes' in hierarchy and 'pes' not in hierarchy) or

mechanalyzer/builder/sorter.py

@@ -15,7 +15,7 @@
     'H5H3ratio': 0,
     'kratio': 1e50, #depends too much on temperature to make it a default that makes sense
     'kabs': 1e50,
-    'lookforpromptchains': 1 
+    'lookforpromptchains': 0 
 }
 
 # Functions to take the mechanism strings (may want to further simplify)
@@ -37,8 +37,8 @@ def sorted_mech(spc_str, mech_str, isolate_spc, sort_lst, spc_therm_dct=None, dc
     srt_mch, rxn_param_dct = _sort_objs(
         spc_str, mech_str, sort_lst, isolate_spc, stereo_optns=stereo_optns)
 
-    pes_groups = None
-    rxns_filter = None
+    pes_groups = ''
+    rxns_filter = ''
     # if prompt groups detected: retrieve grps info
     if 'submech_prompt' in sort_lst and spc_therm_dct and dct_flt_grps:
         DFG.update(dct_flt_grps)

mechanalyzer/calculator/nonboltz.py

@@ -293,7 +293,6 @@ def get_max_reactivity(hot_sp, hot_sp_df, therm_df, T0, Tref):
                      for prds in hot_sp_df['prd_names_lst'].values])
 
     for rxn in hot_sp_df.index:
-        print(rxn)
         rcts = hot_sp_df['rct_names_lst'][rxn]
         prds = hot_sp_df['prd_names_lst'][rxn]
         # used to filter out too fast isomerization channels, but causes failure when only unimol channels are present

mechanalyzer/parser/mech.py

@@ -4,7 +4,7 @@
 
 import sys
 import autoparse.pattern as app
-import ioformat.ptt
+from ioformat import ptt
 from ioformat import remove_comment_lines
 from mechanalyzer.parser import ckin_ as ckin
 
@@ -40,41 +40,41 @@ def parse_sort(sort_str):
             sort_str, delim_pattern=app.escape('!'))
 
     # Read and format information from the isolate_submech block
-    spc_block = ioformat.ptt.end_block(sort_str, 'isolate_submech')
+    spc_block = ptt.end_block(sort_str, 'isolate_submech')
 
     if not spc_block: # this should be checked because spc_block might be None - section not mandatory
         spc_lst = []
     elif not str.isspace(spc_block):
-        spc_lst = list(ioformat.ptt.values_from_block(
+        spc_lst = list(ptt.values_from_block(
             spc_block, val_ptt=app.one_or_more(app.CKINSAFE_CHAR)))
     else:
         spc_lst = []
 
     # Read and format information from the sort_mech block
-    isol_block = ioformat.ptt.end_block(sort_str, 'sort_mech')
+    isol_block = ptt.end_block(sort_str, 'sort_mech')
     # main criteria
-    crit_block = ioformat.ptt.paren_blocks(
+    crit_block = ptt.paren_blocks(
         isol_block, key='criteria')
 
     if crit_block:
-        crit_tup = ioformat.ptt.values_from_block(
+        crit_tup = ptt.values_from_block(
             crit_block[0][1], val_ptt=app.one_or_more(app.URLSAFE_CHAR))
     else:
         crit_tup = ()
 
-    head_block = ioformat.ptt.keyword_value_blocks(
+    head_block = ptt.keyword_value_blocks(
         isol_block, key='n_criteria_headers')
     nhead = int(head_block[0][1]) if head_block is not None else 0
 
-    keepbelow = ioformat.ptt.keyword_value_blocks(
+    keepbelow = ptt.keyword_value_blocks(
         isol_block, key='stoich_keepbelow')
     if keepbelow is not None:
         spc_lst += ['keepbelow ' + keepbelow[0][1].strip(),]
-    deleteabove = ioformat.ptt.keyword_value_blocks(
+    deleteabove = ptt.keyword_value_blocks(
         isol_block, key='stoich_deleteabove')
     if deleteabove is not None:
         spc_lst += ['deleteabove ' + deleteabove[0][1].strip(),]
-    singlespecies = ioformat.ptt.keyword_value_blocks(
+    singlespecies = ptt.keyword_value_blocks(
         isol_block, key='singlespecies')
     if singlespecies is not None:
         if singlespecies[0][1].strip() == 'True':
@@ -87,10 +87,10 @@ def parse_sort(sort_str):
     sort_lst = list(sort_tup)
 
     # prompt criteria
-    prompt_block = ioformat.ptt.end_block(sort_str, 'prompt_filter')
+    prompt_block = ptt.end_block(sort_str, 'prompt_filter')
     prompt_filter_dct = {}
     if prompt_block is not None:
-        prompt_block = ioformat.ptt.keyword_value_blocks(
+        prompt_block = ptt.keyword_value_blocks(
             prompt_block)        
         dct_0 = dict(prompt_block)
         for key in dct_0.keys():
@@ -102,3 +102,23 @@ def parse_sort(sort_str):
         print('*ERROR: sort_mech section is not defined')
 
     return spc_lst, sort_lst, prompt_filter_dct
+
+def parse_classtype(classtype_str):
+    """ parse the class type string. useful to process larger class groups within a (CRECK) mech
+
+    Args:
+        classtype_str (str): string containing reaction types classified in larger groups
+                             inline comments should be removed in advance.
+        reaction_typeclass_dct (dct{reactiontype(str): classtype(str)}):
+        dictionary with correspondence between reaction type and class type
+    """
+
+    grp_blocks = ptt.named_end_blocks(classtype_str, 'classtype', footer='classtype')
+    reaction_typeclass_dct = {}
+    for classtype, block in grp_blocks.items():
+        reactiontypes = block.split()
+        if reactiontypes:
+            for reactiontype in reactiontypes:
+                reaction_typeclass_dct[reactiontype] = classtype
+
+    return reaction_typeclass_dct