diff --git a/arc/plotter.py b/arc/plotter.py index 493cb9c2a0..79f5273959 100644 --- a/arc/plotter.py +++ b/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() diff --git a/arc/processor.py b/arc/processor.py index 2f70d0e425..6ab371950b 100644 --- a/arc/processor.py +++ b/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'], diff --git a/arc/processor_test.py b/arc/processor_test.py index ce0c229a57..f2c138132a 100644 --- a/arc/processor_test.py +++ b/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'), diff --git a/arc/reaction/reaction.py b/arc/reaction/reaction.py index e454709a75..9a73174d8a 100644 --- a/arc/reaction/reaction.py +++ b/arc/reaction/reaction.py @@ -66,7 +66,7 @@ class ARCReaction(object): r_species (List[ARCSpecies]): A list of reactants :ref:`ARCSpecies ` objects. p_species (List[ARCSpecies]): A list of products :ref:`ARCSpecies ` objects. ts_species (ARCSpecies): The :ref:`ARCSpecies ` 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() diff --git a/arc/scripts/rmg_kinetics.py b/arc/scripts/rmg_kinetics.py index 52f6804497..48461a0fe9 100644 --- a/arc/scripts/rmg_kinetics.py +++ b/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 diff --git a/arc/testing/process_kinetics/RMG_kinetics.yml b/arc/testing/process_kinetics/RMG_kinetics.yml index 91e268eae8..d81dc0b6d3 100644 --- a/arc/testing/process_kinetics/RMG_kinetics.yml +++ b/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 diff --git a/arc/testing/process_kinetics/RMG_kinetics_in.yml b/arc/testing/process_kinetics/RMG_kinetics_in.yml new file mode 100644 index 0000000000..9385d9bb93 --- /dev/null +++ b/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} diff --git a/arc/testing/process_kinetics/rate_plots.pdf b/arc/testing/process_kinetics/rate_plots.pdf index 35c0cdb86e..dc55dd97a2 100644 Binary files a/arc/testing/process_kinetics/rate_plots.pdf and b/arc/testing/process_kinetics/rate_plots.pdf differ