TMP

arc/plotter.py

@@ -499,7 +499,7 @@ def draw_parity_plot(var_arc, var_rmg, labels, var_label, var_units, pp=None):
     plt.close(fig=fig)
 
 
-def draw_kinetics_plots(rxn_list: list,  # todo
+def draw_kinetics_plots(rxn_list: list,
                         T_min: Optional[Tuple[float, str]] = None,
                         T_max: Optional[Tuple[float, str]] = None,
                         T_count: int = 50,
@@ -550,7 +550,7 @@ def draw_kinetics_plots(rxn_list: list,  # todo
                 units = r' (cm$^3$/(mol s))'
             elif reaction_order == 3:
                 units = r' (cm$^6$/(mol$^2$ s))'
-            arc_k = [calculate_arrhenius_rate_coefficient(A=rxn.kinetics['A'] * conversion_factor[reaction_order],
+            arc_k = [calculate_arrhenius_rate_coefficient(A=rxn.kinetics['A'],
                                                           n=rxn.kinetics['n'],
                                                           Ea=rxn.kinetics['Ea'],
                                                           T=T,
@@ -562,9 +562,10 @@ def draw_kinetics_plots(rxn_list: list,  # todo
                 else:
                     temps = np.linspace(kinetics['T_min'].value_si, kinetics['T_max'].value_si, T_count)
                 print(f'\n\n------------appending kinetics {kinetics["comment"]} to rmg_rxns')
+                print(f'order and factor: {reaction_order} {conversion_factor[reaction_order]}')
                 rmg_rxns.append({'label': kinetics['comment'],
                                  'T': temps,
-                                 'k': [calculate_arrhenius_rate_coefficient(A=kinetics['A'],
+                                 'k': [calculate_arrhenius_rate_coefficient(A=kinetics['A'] * conversion_factor[reaction_order],
                                                                             n=kinetics['n'],
                                                                             Ea=kinetics['Ea'],
                                                                             T=T,
@@ -598,26 +599,29 @@ def _draw_kinetics_plots(rxn_label, arc_k, temperature, rmg_rxns, units, pp, max
     plt.title(rxn_label)
     inverse_temperature = [1000 / t for t in temperature]
     ax.semilogy(inverse_temperature, arc_k, 'k--', linewidth=2.5, label='ARC')
-    plotted_rmg_rxns = 0
+    plotted_rmg_rxns = list()
     remaining_rmg_rxns = list()
-    for rmg_rxn in rmg_rxns:
+    for i, rmg_rxn in enumerate(rmg_rxns):
         if 'family' in rmg_rxn['label'].lower():
             inverse_temp = [1000 / t for t in rmg_rxn['T']]
+            print(f'plotting 606 rmg_rxn {rmg_rxn["label"]}')
             ax.semilogy(inverse_temp, rmg_rxn['k'], label=rmg_rxn['label'])
-            plotted_rmg_rxns += 1
+            plotted_rmg_rxns.append(i)
         else:
             remaining_rmg_rxns.append(rmg_rxn)
     for priority_lib in kinetics_library_priority:
-        for rmg_rxn in remaining_rmg_rxns:
-            if priority_lib.lower() in rmg_rxn['label'].lower() and plotted_rmg_rxns <= max_rmg_rxns:
+        for i, rmg_rxn in enumerate(remaining_rmg_rxns):
+            if i not in plotted_rmg_rxns and priority_lib.lower() in rmg_rxn['label'].lower() and len(plotted_rmg_rxns) <= max_rmg_rxns:
                 inverse_temp = [1000 / t for t in rmg_rxn['T']]
+                print(f'plotting 615 {priority_lib} rmg_rxn {rmg_rxn["label"]}')
                 ax.semilogy(inverse_temp, rmg_rxn['k'], label=rmg_rxn['label'])
-                plotted_rmg_rxns += 1
-    for rmg_rxn in rmg_rxns:  # todo: it plots it twice, why?
-        if plotted_rmg_rxns <= max_rmg_rxns:
+                plotted_rmg_rxns.append(i)
+    for i, rmg_rxn in enumerate(rmg_rxns):
+        if i not in plotted_rmg_rxns and len(plotted_rmg_rxns) <= max_rmg_rxns:
             inverse_temp = [1000 / t for t in rmg_rxn['T']]
+            print(f'plotting 622 rmg_rxn {rmg_rxn["label"]}')
             ax.semilogy(inverse_temp, rmg_rxn['k'], label=rmg_rxn['label'])
-            plotted_rmg_rxns += 1
+            plotted_rmg_rxns.append(i)
     plt.xlabel(r'1000 / T (K$^-$$^1$)')
     plt.ylabel(f'Rate coefficient{units}')
     plt.legend()

arc/processor.py

@@ -266,12 +266,25 @@ def compare_rates(rxns_for_kinetics_lib: list,
                    content=[{'label': rxn.label,
                              'reactants': [spc.mol.to_adjacency_list() for spc in rxn.r_species],
                              'products': [spc.mol.to_adjacency_list() for spc in rxn.p_species],
-                             'dh_rxn298': rxn.dh_rxn298} for rxn in rxns_for_kinetics_lib])
+                             'dh_rxn298': rxn.dh_rxn298,
+                             'family': rxn.family,
+                             } for rxn in rxns_for_kinetics_lib],
+                   )
     command = f'python {KINETICS_SCRIPT_PATH} {reactions_kinetics_path}'
     out, err = execute_command(command=command, no_fail=True)
+    print('\n\n')
+    print('out:')
+    for o in out:
+        print(o)
+    print(f'err:')
+    for e in err:
+        print(e)
+    print('\n\n\n\n\n')
     reactions_list_w_rmg_kinetics = read_yaml_file(path=reactions_kinetics_path)
     for original_rxn, rxn_w_rmg_kinetics in zip(rxns_for_kinetics_lib, reactions_list_w_rmg_kinetics):
         original_rxn.rmg_kinetics = original_rxn.rmg_kinetics or list()
+        if 'kinetics' not in rxn_w_rmg_kinetics:
+            continue
         for kinetics_entry in rxn_w_rmg_kinetics['kinetics']:
             if 'A' in kinetics_entry and 'n' in kinetics_entry and 'Ea' in kinetics_entry:
                 original_rxn.rmg_kinetics.append({'A': kinetics_entry['A'],

arc/processor_test.py

@@ -46,11 +46,12 @@ def test_compare_rates(self):
         rxn_1 = ARCReaction(r_species=[self.ch4, self.h],
                             p_species=[ARCSpecies(label='CH3', smiles='[CH3]'),
                                        ARCSpecies(label='H2', smiles='[H][H]')],
-                            kinetics={'A': 1e13, 'n': 0.5, 'Ea': 150},
+                            kinetics={'A': 4.79e+05, 'n': 2.5, 'Ea': 40.12},
                             )
         rxn_2 = ARCReaction(r_species=[ARCSpecies(label='nC3H7', smiles='[CH2]CC')],
                             p_species=[ARCSpecies(label='iC3H7', smiles='C[CH]C')],
-                            kinetics={'A': 1e12, 'n': 0.25, 'Ea': 50},
+                            # kinetics={'A': 2.14e13, 'n': -0.40, 'Ea': 163.8},
+                            kinetics={'A': 7.18e5, 'n': 2.05, 'Ea': 151.88},
                             )
         reactions_to_compare = processor.compare_rates(rxns_for_kinetics_lib=[rxn_1, rxn_2],
                                                        output_directory=os.path.join(ARC_PATH, 'arc', 'testing', 'process_kinetics'),

arc/reaction/reaction.py

@@ -66,7 +66,7 @@ class ARCReaction(object):
         r_species (List[ARCSpecies]): A list of reactants :ref:`ARCSpecies <species>` objects.
         p_species (List[ARCSpecies]): A list of products :ref:`ARCSpecies <species>` objects.
         ts_species (ARCSpecies): The :ref:`ARCSpecies <species>` corresponding to the reaction's TS.
-        dh_rxn298 (float): The heat of reaction at 298K.
+        dh_rxn298 (float): The heat of reaction at 298K in J/mol.
         kinetics (Dict[str, float]): The high pressure limit rate coefficient calculated by ARC.
                                      Keys are 'A' in cm-s-mol units, 'n', and 'Ea' in kJ/mol.
         rmg_kinetics (List[Dict[str, float]]): The Arrhenius kinetics from RMG's libraries and families.
@@ -142,10 +142,10 @@ def __init__(self,
                                 f'reactants and {len(self.products)} products for reaction {self.label}.')
         if not isinstance(self.ts_xyz_guess, list):
             self.ts_xyz_guess = [self.ts_xyz_guess]
-        for spc in r_species + p_species:
-            if not isinstance(spc, ARCSpecies):
-                raise InputError(f'All reactants and products must be ARCSpecies objects. Got {spc} which is a '
-                                 f'{type(spc)} object.')
+        # for spc in r_species + p_species:
+        #     if not isinstance(spc, ARCSpecies):
+        #         raise InputError(f'All reactants and products must be ARCSpecies objects. Got {spc} which is a '
+        #                          f'{type(spc)} object.')
         self.remove_dup_species()
         self.check_atom_balance()
 

arc/scripts/rmg_kinetics.py

@@ -7,19 +7,18 @@
 """
 
 import os
-from typing import TYPE_CHECKING, List, Optional, Tuple
+from typing import List, Optional, Tuple
 
 from common import parse_command_line_arguments, read_yaml_file, save_yaml_file
 
 from rmgpy.data.kinetics.common import find_degenerate_reactions
+from rmgpy.data.kinetics.family import KineticsFamily, TemplateReaction
 from rmgpy.data.rmg import RMGDatabase
 from rmgpy import settings as rmg_settings
 from rmgpy.reaction import same_species_lists, Reaction
+from rmgpy.rmg.model import get_family_library_object
 from rmgpy.species import Species
 
-if TYPE_CHECKING:
-    from rmgpy.data.kinetics.family import KineticsFamily
-
 
 DB_PATH = rmg_settings['database.directory']
 
@@ -60,13 +59,15 @@ def get_rate_coefficients(reaction_list: List[dict]) -> List[dict]:
     for i in range(len(reaction_list)):
         rxn = Reaction(reactants=[Species().from_adjacency_list(adjlist) for adjlist in reaction_list[i]['reactants']],
                        products=[Species().from_adjacency_list(adjlist) for adjlist in reaction_list[i]['products']])
-        reaction_list[i]['kinetics'] = determine_rmg_kinetics(rmgdb=rmgdb, reaction=rxn, dh_rxn298=reaction_list[i]['dh_rxn298'])
+        reaction_list[i]['kinetics'] = determine_rmg_kinetics(rmgdb=rmgdb, reaction=rxn, dh_rxn298=reaction_list[i]['dh_rxn298'],
+                                                              family=reaction_list[i]['family'] if 'family' in reaction_list[i] else None)
     return reaction_list
 
 
 def determine_rmg_kinetics(rmgdb: RMGDatabase,
                            reaction: Reaction,
                            dh_rxn298: Optional[float] = None,
+                           family: Optional[str] = None,
                            ) -> List[dict]:
     """
     Determine kinetics for `reaction` (an RMG Reaction object) from RMG's database, if possible.
@@ -76,6 +77,7 @@ def determine_rmg_kinetics(rmgdb: RMGDatabase,
         rmgdb (RMGDatabase): The RMG database instance.
         reaction (Reaction): The RMG Reaction object.
         dh_rxn298 (float, optional): The heat of reaction at 298 K in J/mol.
+        family (str, optional): The RMG family label.
 
     Returns: List[dict]
         All matching RMG reactions kinetics (both libraries and families) as a dict of parameters.
@@ -89,14 +91,24 @@ def determine_rmg_kinetics(rmgdb: RMGDatabase,
                 library_reaction.comment = f'Library: {library.label}'
                 rmg_reactions.append(library_reaction)
                 break
-    # Families:
-    fam_list = loop_families(rmgdb, reaction)
+    # # Families:
+    # family = get_family_library_object(family)
+    # template = family.get_reaction_template(TemplateReaction(reactants=reaction.reactants, products=reaction.products))
+    # # Get the kinetics for the reaction
+    # kinetics, source, entry, is_forward = family.get_kinetics(
+    #     reaction=reaction, template_labels=template, degeneracy=reaction.degeneracy, estimator='rate rules',
+    #     return_all_kinetics=False,
+    # )
+    # print(f'direct kinetics: {kinetics}')
+
+    fam_list = loop_families(rmgdb, reaction)  # todo: n is 0 for H abstraction from CH4, why?
+    dh_rxn298 = dh_rxn298 or get_dh_rxn298(rmgdb=rmgdb, reaction=reaction)  # J/mol
     for family, degenerate_reactions in fam_list:
         for deg_rxn in degenerate_reactions:
             template = family.retrieve_template(deg_rxn.template)
             kinetics = family.estimate_kinetics_using_rate_rules(template)[0]
             kinetics.change_rate(deg_rxn.degeneracy)
-            kinetics = kinetics.to_arrhenius(dh_rxn298 or get_dh_rxn298(rmgdb=rmgdb, reaction=reaction))  # Convert ArrheniusEP to Arrhenius using the dHrxn at 298K
+            kinetics = kinetics.to_arrhenius(dh_rxn298)  # Convert ArrheniusEP to Arrhenius
             deg_rxn.kinetics = kinetics
             deg_rxn.comment = f'Family: {deg_rxn.family}'
             deg_rxn.reactants = reaction.reactants
@@ -168,7 +180,7 @@ def get_kinetics_from_reactions(reactions: List[Reaction]) -> List[dict]:
 
 def loop_families(rmgdb: RMGDatabase,
                   reaction: Reaction,
-                  ) -> List[Tuple['KineticsFamily', list]]:
+                  ) -> List[Tuple[KineticsFamily, list]]:
     """
     Loop through kinetic families and return a list of tuples of (family, degenerate_reactions)
     corresponding to ``reaction``.
@@ -267,7 +279,10 @@ def load_rmg_database() -> RMGDatabase:
     kinetics_libraries = read_yaml_file(path=os.path.join(os.path.dirname(__file__), 'libraries.yaml'))['kinetics']
     thermo_libraries = read_yaml_file(path=os.path.join(os.path.dirname(__file__), 'libraries.yaml'))['thermo']
     rmgdb.load_thermo(path=os.path.join(DB_PATH, 'thermo'), thermo_libraries=thermo_libraries, depository=True, surface=False)
-    rmgdb.load_kinetics(path=os.path.join(DB_PATH, 'kinetics'), reaction_libraries=kinetics_libraries, kinetics_families='default')
+    rmgdb.load_kinetics(path=os.path.join(DB_PATH, 'kinetics'),
+                        reaction_libraries=kinetics_libraries,
+                        kinetics_families='default',
+                        kinetics_depositories=['training'])
     return rmgdb
 
 

arc/testing/process_kinetics/RMG_kinetics.yml

@@ -1,4 +1,5 @@
 - dh_rxn298: null
+  family: H_Abstraction
   kinetics:
   - A: 0.4781000000000001
     Ea: 40.116192000000005
@@ -46,6 +47,7 @@
     multiplicity 2
     1 H u1 p0 c0
 - dh_rxn298: null
+  family: intra_H_migration
   kinetics:
   - A: 20000000000.0
     Ea: 104.60000000000002

arc/testing/process_kinetics/RMG_kinetics_in.yml

@@ -0,0 +1,50 @@
+- dh_rxn298: null
+  label: CH4 + H <=> CH3 + H2
+  products:
+  - |
+    multiplicity 2
+    1 C u1 p0 c0 {2,S} {3,S} {4,S}
+    2 H u0 p0 c0 {1,S}
+    3 H u0 p0 c0 {1,S}
+    4 H u0 p0 c0 {1,S}
+  - |
+    1 H u0 p0 c0 {2,S}
+    2 H u0 p0 c0 {1,S}
+  reactants:
+  - |
+    1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S}
+    2 H u0 p0 c0 {1,S}
+    3 H u0 p0 c0 {1,S}
+    4 H u0 p0 c0 {1,S}
+    5 H u0 p0 c0 {1,S}
+  - |
+    multiplicity 2
+    1 H u1 p0 c0
+- dh_rxn298: null
+  label: nC3H7 <=> iC3H7
+  products:
+  - |
+    multiplicity 2
+    1  C u0 p0 c0 {2,S} {4,S} {5,S} {6,S}
+    2  C u1 p0 c0 {1,S} {3,S} {7,S}
+    3  C u0 p0 c0 {2,S} {8,S} {9,S} {10,S}
+    4  H u0 p0 c0 {1,S}
+    5  H u0 p0 c0 {1,S}
+    6  H u0 p0 c0 {1,S}
+    7  H u0 p0 c0 {2,S}
+    8  H u0 p0 c0 {3,S}
+    9  H u0 p0 c0 {3,S}
+    10 H u0 p0 c0 {3,S}
+  reactants:
+  - |
+    multiplicity 2
+    1  C u1 p0 c0 {2,S} {4,S} {5,S}
+    2  C u0 p0 c0 {1,S} {3,S} {6,S} {7,S}
+    3  C u0 p0 c0 {2,S} {8,S} {9,S} {10,S}
+    4  H u0 p0 c0 {1,S}
+    5  H u0 p0 c0 {1,S}
+    6  H u0 p0 c0 {2,S}
+    7  H u0 p0 c0 {2,S}
+    8  H u0 p0 c0 {3,S}
+    9  H u0 p0 c0 {3,S}
+    10 H u0 p0 c0 {3,S}