TMP

arc/plotter.py

@@ -28,7 +28,8 @@
 from rmgpy.quantity import ScalarQuantity
 from rmgpy.species import Species
 
-from arc.common import (extremum_list,
+from arc.common import (calculate_arrhenius_rate_coefficient,
+                        extremum_list,
                         get_angle_in_180_range,
                         get_close_tuple,
                         get_logger,
@@ -529,7 +530,6 @@ def draw_kinetics_plots(rxn_list: list,  # todo
     T_min = ScalarQuantity(value=T_min[0], units=T_min[1])
     T_max = ScalarQuantity(value=T_max[0], units=T_max[1])
     temperatures = np.linspace(T_min.value_si, T_max.value_si, T_count)
-    pressure = 1e7  # Pa  (=100 bar)
 
     pp = None
     if path is not None:
@@ -540,6 +540,7 @@ def draw_kinetics_plots(rxn_list: list,  # todo
 
     for rxn in rxn_list:
         if rxn.kinetics is not None:
+            print(f'Plotting kinetics for {rxn.label} with kinetics {rxn.kinetics}')
             reaction_order = len(rxn.reactants)
             units = ''
             conversion_factor = {1: 1, 2: 1e6, 3: 1e12}
@@ -549,40 +550,46 @@ 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 = [rxn.kinetics.get_rate_coefficient(T, pressure) * conversion_factor[reaction_order] for T in temperatures]
+            arc_k = [calculate_arrhenius_rate_coefficient(A=rxn.kinetics['A'] * conversion_factor[reaction_order],
+                                                          n=rxn.kinetics['n'],
+                                                          Ea=rxn.kinetics['Ea'],
+                                                          T=T,
+                                                          ) for T in temperatures]
             rmg_rxns = list()
             for kinetics in rxn.rmg_kinetics:
-                if kinetics['T_max'] is None or kinetics['T_min'] is None:
+                if kinetics.get('T_max', None) is None or kinetics.get('T_min', None) is None:
                     temps = temperatures
                 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')
                 rmg_rxns.append({'label': kinetics['comment'],
                                  'T': temps,
-                                 'k': [calculate_arrhenius_rate_coefficient(kinetics['A'], kinetics['n'], kinetics['Ea'], T)
+                                 'k': [calculate_arrhenius_rate_coefficient(A=kinetics['A'],
+                                                                            n=kinetics['n'],
+                                                                            Ea=kinetics['Ea'],
+                                                                            T=T,
+                                                                            )
                                        for T in temperatures]})
+            print(f'\ndrawing arc k: {rxn.kinetics}\nand rmg ks: {rxn.rmg_kinetics}\n\n')
             _draw_kinetics_plots(rxn.label, arc_k, temperatures, rmg_rxns, units, pp)
 
     if path is not None:
         pp.close()
 
 
-def calculate_arrhenius_rate_coefficient(A: float, n: float, Ea: float, T: float) -> float:
+def _draw_kinetics_plots(rxn_label, arc_k, temperature, rmg_rxns, units, pp, max_rmg_rxns=5):
     """
-    Calculate the Arrhenius rate coefficient.
+    Draw the kinetics plots.
 
     Args:
-        A (float): Pre-exponential factor.
-        n (float): Temperature exponent.
-        Ea (float): Activation energy in J/mol.
-        T (float): Temperature in Kelvin.
-
-    Returns:
-        float: The rate coefficient at the specified temperature.
+        rxn_label (str): The reaction label.
+        arc_k (list): The ARC rate coefficients.
+        temperature (np.ndarray): The temperatures.
+        rmg_rxns (list): The RMG rate coefficients.
+        units (str): The units of the rate coefficients.
+        pp (PdfPages): Used for storing the image as a multipage PDF file.
+        max_rmg_rxns (int, optional): The maximum number of RMG rate coefficients to plot.
     """
-    return A * (T ** n) * math.exp(-1 * Ea / (R * T))
-
-
-def _draw_kinetics_plots(rxn_label, arc_k, temperature, rmg_rxns, units, pp, max_rmg_rxns=5):
     kinetics_library_priority = ['primaryH2O2', 'Klippenstein_Glarborg2016', 'primaryNitrogenLibrary',
                                  'primarySulfurLibrary', 'N-S_interactions', 'NOx2018',
                                  'Nitrogen_Dean_and_Bozzelli', 'FFCM1(-)', 'JetSurF2.0']
@@ -606,7 +613,7 @@ def _draw_kinetics_plots(rxn_label, arc_k, temperature, rmg_rxns, units, pp, max
                 inverse_temp = [1000 / t for t in rmg_rxn['T']]
                 ax.semilogy(inverse_temp, rmg_rxn['k'], label=rmg_rxn['label'])
                 plotted_rmg_rxns += 1
-    for rmg_rxn in rmg_rxns:
+    for rmg_rxn in rmg_rxns:  # todo: it plots it twice, why?
         if plotted_rmg_rxns <= max_rmg_rxns:
             inverse_temp = [1000 / t for t in rmg_rxn['T']]
             ax.semilogy(inverse_temp, rmg_rxn['k'], label=rmg_rxn['label'])

arc/processor.py

@@ -2,6 +2,7 @@
 Processor module for computing thermodynamic properties and rate coefficients using statistical mechanics.
 """
 
+import math
 import os
 import shutil
 from typing import Optional
@@ -17,6 +18,8 @@
 
 THERMO_SCRIPT_PATH = os.path.join(ARC_PATH, 'arc', 'scripts', 'rmg_thermo.py')
 KINETICS_SCRIPT_PATH = os.path.join(ARC_PATH, 'arc', 'scripts', 'rmg_kinetics.py')
+R = 8.31446261815324  # J/(mol*K)
+EA_UNIT_CONVERSION = {'J/mol': 1, 'kJ/mol': 1e+3, 'cal/mol': 4.184, 'kcal/mol': 4.184e+3}
 
 
 def process_arc_project(thermo_adapter: str,
@@ -237,12 +240,12 @@ def compare_thermo(species_for_thermo_lib: list,
         plotter.draw_thermo_parity_plots(species_list=species_to_compare, path=output_directory)
 
 
-def compare_rates(rxns_for_kinetics_lib: list,  # todo: draw kinetics plot need updaye
+def compare_rates(rxns_for_kinetics_lib: list,
                   output_directory: str,
                   T_min: tuple = None,
                   T_max: tuple = None,
                   T_count: int = 50,
-                  ) -> None:
+                  ) -> list:
     """
     Compare the calculated rates with RMG's estimations or libraries.
 
@@ -252,33 +255,38 @@ def compare_rates(rxns_for_kinetics_lib: list,  # todo: draw kinetics plot need
         T_min (tuple, optional): The minimum temperature for kinetics computations, e.g., (500, 'K').
         T_max (tuple, optional): The maximum temperature for kinetics computations, e.g., (3000, 'K').
         T_count (int, optional): The number of temperature points between ``T_min`` and ``T_max``.
+
+    Returns:
+        list: Reactions for which a rate was both calculated and estimated.
+              Returning this list for testing purposes.
     """
     reactions_to_compare = list()  # reactions for which a rate was both calculated and estimated.
     reactions_kinetics_path = os.path.join(output_directory, 'RMG_kinetics.yml')
     save_yaml_file(path=reactions_kinetics_path,
-                     content=[{'label': rxn.label,
-                              'reactants': [spc.label for spc in rxn.r_species],
-                              'products': [spc.label for spc in rxn.p_species],
-                              'dh_rxn298': rxn.dh_rxn298} for rxn in rxns_for_kinetics_lib])
+                   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])
     command = f'python {KINETICS_SCRIPT_PATH} {reactions_kinetics_path}'
-    execute_command(command=command, no_fail=True)
-    reactions_list = read_yaml_file(path=reactions_kinetics_path)
-    for original_rxn, rmg_rxn in zip(rxns_for_kinetics_lib, reactions_list):
+    out, err = execute_command(command=command, no_fail=True)
+    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()
-        for kinetics_entry in rmg_rxn['kinetics']:
+        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'],
                                                   'n': kinetics_entry['n'],
                                                   'Ea': kinetics_entry['Ea'],
                                                   'comment': kinetics_entry['comment'],
                                                   })
-                reactions_to_compare.append(original_rxn)
+        reactions_to_compare.append(original_rxn)
     if reactions_to_compare:
         plotter.draw_kinetics_plots(reactions_to_compare,
                                     T_min=T_min,
                                     T_max=T_max,
                                     T_count=T_count,
                                     path=output_directory)
+    return reactions_to_compare
 
 
 def compare_transport(species_for_transport_lib: list,

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,34 @@ 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': 1e13, 'n': 0.5, 'Ea': 150},
+                            )
+        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},
+                            )
+        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/reaction/reaction.py

@@ -53,6 +53,8 @@ class ARCReaction(object):
         preserve_param_in_scan (list, optional): Entries are length two iterables of atom indices (1-indexed)
                                                  between which distances and dihedrals of these pivots must be
                                                  preserved. Used for identification of rotors which break a TS.
+        kinetics (Dict[str, float], optional): The high pressure limit rate coefficient calculated by ARC.
+                                               Keys are 'A' in cm-s-mol units, 'n', and 'Ea' in kJ/mol.
 
     Attributes:
         label (str): The reaction's label in the format `r1 + r2 <=> p1 + p2`
@@ -65,9 +67,11 @@ class ARCReaction(object):
         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.
-        kinetics (Arrhenius): The high pressure limit rate coefficient calculated by ARC.
-        rmg_kinetics (List[dict]): The Arrhenius kinetics from RMG's libraries and families. Each dict has 'A', 'n', amd 'Ea'
-                                   parameters as keys, and a 'comment' key with a description of the source of the kinetics.
+        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.
+                                               Each dict has 'A' in cm-s-mol units, 'n', and 'Ea' in kJ/mol as keys,
+                                               and a 'comment' key with a description of the source of the kinetics.
         long_kinetic_description (str): A description for the species entry in the thermo library outputted.
         ts_xyz_guess (list): A list of TS XYZ user guesses, each in a string format.
         multiplicity (int): The reaction surface multiplicity. A trivial guess will be made unless provided.
@@ -96,12 +100,13 @@ def __init__(self,
                  reaction_dict: Optional[dict] = None,
                  species_list: Optional[List[ARCSpecies]] = None,
                  preserve_param_in_scan: Optional[list] = None,
+                 kinetics: Dict[str, float] = None,
                  ):
         self.arrow = ' <=> '
         self.plus = ' + '
         self.r_species = r_species or list()
         self.p_species = p_species or list()
-        self.kinetics = None
+        self.kinetics = kinetics
         self.rmg_kinetics = None
         self.long_kinetic_description = ''
         self._family = None

arc/scripts/rmg_kinetics.py

@@ -55,6 +55,7 @@ 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']],
@@ -158,9 +159,9 @@ 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
 

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'])

arc/testing/process_kinetics/RMG_kinetics.yml

@@ -0,0 +1,86 @@
+- dh_rxn298: null
+  kinetics:
+  - A: 0.4781000000000001
+    Ea: 40.116192000000005
+    T_max: null
+    T_min: null
+    comment: 'Library: FFCM1(-)'
+    kinetics: Arrhenius(A=(478100,'cm^3/(mol*s)'), n=2.5, Ea=(9588,'cal/mol'), T0=(1,'K'))
+    n: 2.5
+  - A: 0.0041
+    Ea: 36.630919999999996
+    T_max: null
+    T_min: null
+    comment: 'Library: NOx2018'
+    kinetics: Arrhenius(A=(4100,'cm^3/(mol*s)'), n=3.156, Ea=(8755,'cal/mol'), T0=(1,'K'))
+    n: 3.156
+  - A: 0.40000000000000013
+    Ea: 41.84
+    T_max: 1500.0
+    T_min: 300.0
+    comment: 'Family: H_Abstraction'
+    kinetics: |-
+      Arrhenius(A=(400000,'cm^3/(mol*s)'), n=0, Ea=(41.84,'kJ/mol'), T0=(1,'K'), Tmin=(300,'K'), Tmax=(1500,'K'), comment="""Estimated using template [X_H;Y_rad_birad_trirad_quadrad] for rate rule [C_methane;H_rad]
+      Euclidian distance = 2.8284271247461903
+      family: H_Abstraction""")
+    n: 0.0
+  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
+  kinetics:
+  - A: 20000000000.0
+    Ea: 104.60000000000002
+    T_max: 1500.0
+    T_min: 300.0
+    comment: 'Family: intra_H_migration'
+    kinetics: |-
+      Arrhenius(A=(2e+10,'s^-1'), n=0, Ea=(104.6,'kJ/mol'), T0=(1,'K'), Tmin=(300,'K'), Tmax=(1500,'K'), comment="""Estimated using template [RnH;Y_rad_out;XH_out] for rate rule [R2H_S;C_rad_out_2H;Cs_H_out_H/NonDeC]
+      Euclidian distance = 4.69041575982343
+      family: intra_H_migration""")
+    n: 0.0
+  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}