TMP

arc/processor_test.py

@@ -5,10 +5,13 @@
 This module contains unit tests for the arc.processor module
 """
 
+import os
 import unittest
 
 import arc.processor as processor
-from arc.species.species import ARCSpecies
+from arc.common import ARC_PATH
+from arc.reaction import ARCReaction
+from arc.species import ARCSpecies
 
 
 class TestProcessor(unittest.TestCase):
@@ -30,14 +33,35 @@ def setUpClass(cls):
         cls.ch4.bdes = [(1, 2)]
 
     def test_process_bdes(self):
-        """Test the process_bdes method"""
+        """Test the process_bdes() method"""
         bde_report = processor.process_bdes(label='CH4',
                                             species_dict={'CH4': self.ch4,
                                                           'NH3': self.nh3,
                                                           'H': self.h,
                                                           'CH4_BDE_1_2_A': self.ch4_bde_1_2_a})
         self.assertEqual(bde_report, {(1, 2): 50})
 
+    def test_compare_rates(self):
+        """Test the compare_rates() method"""
+        rxn_1 = ARCReaction(r_species=[self.ch4, self.h],
+                            p_species=[ARCSpecies(label='CH3', smiles='[CH3]'),
+                                       ARCSpecies(label='H2', smiles='[H][H]')],
+                            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': 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'),
+                                                       )
+        print('\n********')
+        print(f'reactions_to_compare: {reactions_to_compare}')
+        raise
+
+
+
 
 if __name__ == '__main__':
     unittest.main(testRunner=unittest.TextTestRunner(verbosity=2))

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']
 
@@ -55,17 +54,20 @@ def get_rate_coefficients(reaction_list: List[dict]) -> List[dict]:
     Returns:
         List[dict]: The list of reactions with rate coefficients.
     """
+    print('Loading RMG database...')
     rmgdb = load_rmg_database()
     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.
@@ -75,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.
@@ -88,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
@@ -158,16 +171,16 @@ def get_kinetics_from_reactions(reactions: List[Reaction]) -> List[dict]:
             'comment': rxn.comment,
             'A': rxn.kinetics.A.value_si if hasattr(rxn.kinetics, 'A') else None,
             'n': rxn.kinetics.n.value_si if hasattr(rxn.kinetics, 'n') else None,
-            'Ea': rxn.kinetics.Ea.value_si if hasattr(rxn.kinetics, 'Ea') else None,
-            'T_min': rxn.kinetics.Tmin.value_si if hasattr(rxn.kinetics, 'Tmin') else None,
-            'T_max': rxn.kinetics.Tmax.value_si if hasattr(rxn.kinetics, 'Tmax') else None,
+            'Ea': rxn.kinetics.Ea.value_si * 0.001 if hasattr(rxn.kinetics, 'Ea') else None,  # kJ/mol
+            'T_min': rxn.kinetics.Tmin.value_si if hasattr(rxn.kinetics, 'Tmin') and rxn.kinetics.Tmin is not None else None,
+            'T_max': rxn.kinetics.Tmax.value_si if hasattr(rxn.kinetics, 'Tmax') and rxn.kinetics.Tmax is not None else None,
         })
     return kinetics_list
 
 
 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``.
@@ -266,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/scripts/rmg_thermo.py

@@ -51,6 +51,7 @@ def get_thermo(species_list: List[dict]) -> List[dict]:
     Returns:
         List[dict]: A list of species dictionaries with thermo properties.
     """
+    print('Loading RMG database...')
     rmgdb = load_rmg_database()
     for i in range(len(species_list)):
         spc = Species(label=species_list[i]['label'])