From 06f30f6800b8453699d8c53d43f203138dad2468 Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Sun, 6 Jun 2021 12:31:04 -0700 Subject: [PATCH 01/17] categorization script --- src/mrnet/utils/reaction_category.py | 289 +++++++++++++++++++++++++++ src/mrnet/utils/visualization.py | 2 +- 2 files changed, 290 insertions(+), 1 deletion(-) create mode 100644 src/mrnet/utils/reaction_category.py diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py new file mode 100644 index 00000000..204a3242 --- /dev/null +++ b/src/mrnet/utils/reaction_category.py @@ -0,0 +1,289 @@ +# coding: utf-8 +# Copyright (c) Pymatgen Development Team. +# Distributed under the terms of the MIT License. +from ast import literal_eval, operator +from functools import reduce +import mrnet +from mrnet.stochastic.analyze import SimulationAnalyzer, NetworkUpdater + +__author__ = "Hetal D. Patel" +__maintainer__ = "Hetal D. Patel" +__email__ = "hpatel@lbl.gov" +__version__ = "0.1" +__status__ = "Alpha" +__date__ = "May, 2021" + + + + +def reaction_category_RNMC(network_folder, entriesbox, reaction_ids): + + category_dict = {"single_element_reactant":[], "redox":[],"Li_hopping":[], "Li_coord_change":[], + "H_or_F_abstraction":[], "actual_bond_changes": [], "non_local_reaction_center": []} + sa = SimulationAnalyzer(network_folder, entriesbox) + for rxn_id in reaction_ids: + r_info = sa.index_to_reaction(rxn_id) + reactants = [] + products = [] + for r in r_info["reactants"]: + reactants.append(entriesbox.entries_list[r]) + for p in r_info["products"]: + products.append(entriesbox.entries_list[p]) + r = [reactants, products] + rxn_type = reaction_category(r) + category_dict[rxn_type].append(rxn_id) + + +def update_rates_RNMC(network_folder, category_dict, rates=None): + + update = [] + if rates is None: + rates = {"redox": 0.1, "Li_hopping": 0.1, "Li_coord_change": 0.1, "H_or_F_abstraction": 0.1} + + for rxn_type in rates: + rate = rates[rxn_type] + for rxn_id in category_dict[rxn_type]: + update.append((rxn_id, rate)) + + network_updater = NetworkUpdater(network_folder) + network_updater.update_rates(update) + +def reaction_category(r): #r = [[reactnat molecule entries], [product moleucle entries]] + + + all_reactants = [] + react_global_bonds = [] + prod_global_bonds = [] + r_charge = 0 + p_charge = 0 + single_elem_react_or_prod = False + single_elemt_Li = False + LiF = False + Li_hopping = False + + for i in r[0]: + r_charge = r_charge + i.charge + all_reactants = all_reactants + i.species + if i.formula == "F1 Li1": + LiF = True + elif len(i.species) == 1: + single_elem_react_or_prod = True + if i.formula == "Li1": + single_elemt_Li = True + for i in r[1]: + p_charge = p_charge + i.charge + + if i.formula == "F1 Li1": + LiF = True + elif len(i.species) == 1: + print(i.species) + single_elem_react_or_prod = True + if i.formula == "Li1": + single_elemt_Li = True + + + + + mapping_dict = mrnet.core.reactions.get_reaction_atom_mapping(r[0], r[1]) + print("reactant speices", all_reactants) + all_reactant_species = all_reactants + for r_id, react in enumerate(r[0]): + for b in react.bonds: + s = (mapping_dict[0][r_id][b[0]], mapping_dict[0][r_id][b[1]]) + s = sorted(s) + react_global_bonds.append(s) + for p_id, prod in enumerate(r[1]): + for b in prod.bonds: + s = (mapping_dict[1][p_id][b[0]], mapping_dict[1][p_id][b[1]]) + s = sorted(s) + prod_global_bonds.append(s) + print(react_global_bonds, prod_global_bonds) + r_set = set(map(lambda x: repr(x), react_global_bonds)) + p_set = set(map(lambda x: repr(x), prod_global_bonds)) + diff_r_p = list(map(lambda y: literal_eval(y), r_set - p_set)) + diff_p_r = list(map(lambda y: literal_eval(y), p_set - r_set)) + print("diff_p-r", diff_p_r) + print("diff_r-p", diff_r_p) + + if diff_r_p != []: + r_p = reduce(operator.concat, diff_r_p) + else: + r_p = [] + if diff_p_r != []: + p_r = reduce(operator.concat, diff_p_r) + else: + p_r = [] + + print("!!", r_p, p_r) + + # if len(list((set(r_p) & set(p_r)))) != 0: + # print("reaction center behaving") + Li_ind = [i for i, x in enumerate(all_reactant_species) if x == "Li"] + # if Li_ind != []: + + if r_charge != p_charge: + print("redox") + + + + elif r_p == [] or p_r == []: + print("either only bond forming or breaking") + combined_diff = diff_p_r + diff_r_p + print(combined_diff) + if combined_diff == []: + print("UNSURE") + elif list(set(combined_diff[0]).intersection(*combined_diff)) == []: # [(10,1), (10,9), (10,3)] + print("either only forming or breaking bonds during the reaction AND not following rxn center") + print("non_local_reaction_center") + if LiF: + print( + "LiF forming - either only forming or breaking bonds during the reaction AND not following rxn center") + print("LiF forming - non_local_reaction_center") + + elif single_elem_react_or_prod: + print("single element bond forming or breaking") + if single_elemt_Li: + print("Li coordination bond Li+A <> LiA") + else: + print("for AutoTS") + + + elif (set(Li_ind) & set(r_p)) or (set(Li_ind) & set(p_r)): # [10,4,5] + # print((set(Li_ind) & set(r_p)),(set(Li_ind) & set(p_r))) + if len(r[0]) == 1 and len(r[1]) == 1: + print("change in corrdination within a molecule") + elif LiF: + print("LiF coordinating") + else: + print("for AutoTS") + else: + print("coord and covalent bond break/form") + + + else: + print("bonds breaking AND forming", single_elem_react_or_prod, Li_ind) + + if len(list((set(r_p) & set(p_r)))) == 0: # reaction center + # list(set(diff_r_p[0]).intersection(*diff_r_p)) == [] or list(set(diff_p_r[0]).intersection(*diff_p_r)) + # == []: + # print("change in corrdination within a molecule") + print("non_local_reaction_center") + + + elif (set(Li_ind) & set(r_p)) and ( + set(Li_ind) & set(p_r)): # some sort of Li bonding AND Li breaking is involved + Li_changing = (set(Li_ind) & set(r_p) & set(p_r)) + if len(r[0]) == 1 and len(r[1]) == 1: + print("change in corrdination within a molecule") + elif "F" in all_reactant_species: + F_ind = [i for i, x in enumerate(all_reactant_species) if x == "F"] + + edge = [[x, y] for x in F_ind for y in Li_changing] + for e in edge: + e.sort() + if e in react_global_bonds and e in prod_global_bonds: + Li_hopping = True + if Li_hopping: + print("LiF hopping") + else: + print("Li_hopping - inside") + else: + print("Li_hopping") + + + + elif (set(Li_ind) & set(r_p)) or (set(Li_ind) & set(p_r)): + if single_elemt_Li: + # print((set(Li_ind) & set(r_p)),(set(Li_ind) & set(p_r))) + print("coord and covalent bond break/form") + else: + print("for AutoTS") + + elif single_elem_react_or_prod: + print("single element bond forming or breaking") + if single_elemt_Li: + print("Li coordination bond Li+A <> LiA") + else: + print("for AutoTS") + else: + print("for AutoTS") + + + + + + all_reactant_species = [] + react_global_bonds = [] + prod_global_bonds = [] + r_charge = 0 + p_charge = 0 + single_elem_reactant = False + + for i in r[0]: + r_charge = r_charge + i.charge + all_reactant_species = all_reactant_species + i.species + if len(i.species) == 1: + single_elem_reactant = True + for i in r[1]: + p_charge = p_charge + i.charge + + if single_elem_reactant: + return "single_element_reactant" + elif r_charge != p_charge: + return "redox" + else: + mapping_dict = mrnet.core.reactions.get_reaction_atom_mapping(r[0], r[1]) + for r_id, react in enumerate(r[0]): + for b in react.bonds: + s = (mapping_dict[0][r_id][b[0]], mapping_dict[0][r_id][b[1]]) + s = sorted(s) + react_global_bonds.append(s) + for p_id, prod in enumerate(r[1]): + for b in prod.bonds: + s = (mapping_dict[1][p_id][b[0]], mapping_dict[1][p_id][b[1]]) + s = sorted(s) + prod_global_bonds.append(s) + r_set = set(map(lambda x: repr(x), react_global_bonds)) + p_set = set(map(lambda x: repr(x), prod_global_bonds)) + diff_r_p = list(map(lambda y: literal_eval(y), r_set - p_set)) + diff_p_r = list(map(lambda y: literal_eval(y), p_set - r_set)) + print("p-r", diff_p_r) + print("r-p", diff_r_p) + if diff_r_p != []: + r_p = reduce(operator.concat, diff_r_p) + else: + r_p = [] + if diff_p_r != []: + p_r = reduce(operator.concat, diff_p_r) + else: + p_r = [] + + if r_p == [] or p_r == []: + print("either only forming or breaking bonds during the reaction") + if list(set(diff_r_p[0]).intersection(*diff_r_p)) == []: + print("either only forming or breaking bonds during the reaction AND not following rxn center") + return "non_local_reaction_center" + elif len(list((set(r_p) & set(p_r)))) != 0: + print("reaction center behaving") + Li_ind = [i for i, x in enumerate(all_reactant_species) if x == "Li"] + if Li_ind != []: + if (set(Li_ind) & set(r_p)) and (set(Li_ind) & set(p_r)): + print("in both - so hopping") + return "Li_hopping" + elif (set(Li_ind) & set(r_p)) or (set(Li_ind) & set(p_r)): + print("only in one - so coord") + return "Li_coord_change" + elif "H" in all_reactant_species or "F" in all_reactant_species: + print("actual bond change?") + if len(set(r_p) & set(p_r)) == 1 and all_reactant_species[list(set(r_p) & set(p_r))[0]] == "H" or all_reactant_species[ + list(set(r_p) & set(p_r))[0]] == "F": + print("H or F abstraction") + return "H_or_F_abstraction" + else: + return "actual_bond_changes" + # print(r_p, p_r, list(set(r_p) & set(p_r))[0]) + + + + + diff --git a/src/mrnet/utils/visualization.py b/src/mrnet/utils/visualization.py index 0208073d..699a17f3 100644 --- a/src/mrnet/utils/visualization.py +++ b/src/mrnet/utils/visualization.py @@ -26,7 +26,7 @@ def visualize_molecule_entry(molecule_entry, path): output the resulting pdf to path """ - atom_colors = {"O": "red", "H": "gray", "C": "black", "Li": "purple"} + atom_colors = {"O": "red", "H": "gray", "C": "black", "Li": "purple", "F": "green", "P": "orange"} graph = deepcopy(molecule_entry.graph).to_undirected() From 3d9aa38433d2a31904bd5331960cad706b9fe7c8 Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Sun, 6 Jun 2021 12:38:50 -0700 Subject: [PATCH 02/17] updates to category scripts --- src/mrnet/utils/reaction_category.py | 250 ++++++++++----------------- 1 file changed, 90 insertions(+), 160 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 204a3242..7118c87f 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -3,8 +3,11 @@ # Distributed under the terms of the MIT License. from ast import literal_eval, operator from functools import reduce +import math + import mrnet from mrnet.stochastic.analyze import SimulationAnalyzer, NetworkUpdater +from mrnet.utils.constants import KB, ROOM_TEMP, PLANCK __author__ = "Hetal D. Patel" __maintainer__ = "Hetal D. Patel" @@ -14,12 +17,17 @@ __date__ = "May, 2021" - - def reaction_category_RNMC(network_folder, entriesbox, reaction_ids): - - category_dict = {"single_element_reactant":[], "redox":[],"Li_hopping":[], "Li_coord_change":[], - "H_or_F_abstraction":[], "actual_bond_changes": [], "non_local_reaction_center": []} + """ + Method to categorize reactions from RNMC based on reaction ids + :param network_folder: path to network folder + :param entriesbox: EntriesBox instance of the entries + :param reaction_ids: list of RNMC reaction ids + :return: category_dict: dict with key being the type of reactions and item being + a list of reaction ids + """ + + category_dict = {} sa = SimulationAnalyzer(network_folder, entriesbox) for rxn_id in reaction_ids: r_info = sa.index_to_reaction(rxn_id) @@ -31,25 +39,47 @@ def reaction_category_RNMC(network_folder, entriesbox, reaction_ids): products.append(entriesbox.entries_list[p]) r = [reactants, products] rxn_type = reaction_category(r) - category_dict[rxn_type].append(rxn_id) - - -def update_rates_RNMC(network_folder, category_dict, rates=None): - + if rxn_type not in category_dict.keys(): + category_dict[rxn_type] = [] + else: + category_dict[rxn_type].append(rxn_id) + return category_dict + + +def update_rates_RNMC(network_folder, category_dict, barrier=None): + """ + Method to update rates of reactions based on the type of reactions and its reaction id + :param network_folder: path to network folder + :param category_dict: category_dict: dict with key being the type of reactions and item + being a list of reaction ids + :param rates: dict with key being the category of reactions to udpate barriers for + and the value being the barrier + value to update to + """ update = [] - if rates is None: - rates = {"redox": 0.1, "Li_hopping": 0.1, "Li_coord_change": 0.1, "H_or_F_abstraction": 0.1} - - for rxn_type in rates: - rate = rates[rxn_type] + if barrier is None: + barrier = { + "Li_hopping": -0.24, + "Li_coord_change": -0.24, + } + kT = KB * ROOM_TEMP + max_rate = kT / PLANCK + for rxn_type in barrier: + rate = max_rate * math.exp(barrier[rxn_type] / kT) for rxn_id in category_dict[rxn_type]: update.append((rxn_id, rate)) network_updater = NetworkUpdater(network_folder) network_updater.update_rates(update) -def reaction_category(r): #r = [[reactnat molecule entries], [product moleucle entries]] +def reaction_category(r): + """ + Mehtod to categroize a single reaction + :param r: list of reactant and product MoleculeEntry [[reactnat molecule entries], + [product moleucle entries]] + :return: string indicating type of reaction + """ all_reactants = [] react_global_bonds = [] @@ -81,11 +111,7 @@ def reaction_category(r): #r = [[reactnat molecule entries], [product moleucle e if i.formula == "Li1": single_elemt_Li = True - - - mapping_dict = mrnet.core.reactions.get_reaction_atom_mapping(r[0], r[1]) - print("reactant speices", all_reactants) all_reactant_species = all_reactants for r_id, react in enumerate(r[0]): for b in react.bonds: @@ -97,13 +123,10 @@ def reaction_category(r): #r = [[reactnat molecule entries], [product moleucle e s = (mapping_dict[1][p_id][b[0]], mapping_dict[1][p_id][b[1]]) s = sorted(s) prod_global_bonds.append(s) - print(react_global_bonds, prod_global_bonds) r_set = set(map(lambda x: repr(x), react_global_bonds)) p_set = set(map(lambda x: repr(x), prod_global_bonds)) diff_r_p = list(map(lambda y: literal_eval(y), r_set - p_set)) diff_p_r = list(map(lambda y: literal_eval(y), p_set - r_set)) - print("diff_p-r", diff_p_r) - print("diff_r-p", diff_r_p) if diff_r_p != []: r_p = reduce(operator.concat, diff_r_p) @@ -114,176 +137,83 @@ def reaction_category(r): #r = [[reactnat molecule entries], [product moleucle e else: p_r = [] - print("!!", r_p, p_r) - - # if len(list((set(r_p) & set(p_r)))) != 0: - # print("reaction center behaving") Li_ind = [i for i, x in enumerate(all_reactant_species) if x == "Li"] - # if Li_ind != []: - - if r_charge != p_charge: - print("redox") - + if r_charge != p_charge: # redox + return "redox" - elif r_p == [] or p_r == []: - print("either only bond forming or breaking") + elif r_p == [] or p_r == []: # bonds are either only forming or only breaking combined_diff = diff_p_r + diff_r_p - print(combined_diff) - if combined_diff == []: - print("UNSURE") - elif list(set(combined_diff[0]).intersection(*combined_diff)) == []: # [(10,1), (10,9), (10,3)] - print("either only forming or breaking bonds during the reaction AND not following rxn center") - print("non_local_reaction_center") + if ( + combined_diff == [] + ): # there is no difference in reactant and product graphs, so change transfer maybe + # occuring + return "uncategorized" + elif ( + list(set(combined_diff[0]).intersection(*combined_diff)) == [] + ): # reaction center rule not satisfied if LiF: - print( - "LiF forming - either only forming or breaking bonds during the reaction AND not following rxn center") - print("LiF forming - non_local_reaction_center") + return "non_local_reaction_center_LiF_formning" + else: + return "non_local_reaction_center" elif single_elem_react_or_prod: - print("single element bond forming or breaking") if single_elemt_Li: - print("Li coordination bond Li+A <> LiA") + return "Li_coordination_Li+A_to_LiA" else: - print("for AutoTS") - + return "AutoTS" - elif (set(Li_ind) & set(r_p)) or (set(Li_ind) & set(p_r)): # [10,4,5] - # print((set(Li_ind) & set(r_p)),(set(Li_ind) & set(p_r))) + elif (set(Li_ind) & set(r_p)) or ( + set(Li_ind) & set(p_r) + ): # bond changes involve Li if len(r[0]) == 1 and len(r[1]) == 1: - print("change in corrdination within a molecule") + return "coordination_change_within_molecule" # EC monodentate <> EC bidentate elif LiF: - print("LiF coordinating") + return "LiF_coordinating" # A + LiF <> ALiF else: - print("for AutoTS") + print("AutoTS") else: - print("coord and covalent bond break/form") - - - else: - print("bonds breaking AND forming", single_elem_react_or_prod, Li_ind) - - if len(list((set(r_p) & set(p_r)))) == 0: # reaction center - # list(set(diff_r_p[0]).intersection(*diff_r_p)) == [] or list(set(diff_p_r[0]).intersection(*diff_p_r)) - # == []: - # print("change in corrdination within a molecule") - print("non_local_reaction_center") + print( + "coord_and_covalent_bond_changes" + ) # ex. Li coordination causes covalent bond breakage + else: # bonds are being broken AND formed + if len(list((set(r_p) & set(p_r)))) == 0: # check for reaction center + return "non_local_reaction_center" elif (set(Li_ind) & set(r_p)) and ( - set(Li_ind) & set(p_r)): # some sort of Li bonding AND Li breaking is involved - Li_changing = (set(Li_ind) & set(r_p) & set(p_r)) + set(Li_ind) & set(p_r) + ): # some sort of Li bonding AND Li breaking + Li_changing = set(Li_ind) & set(r_p) & set(p_r) if len(r[0]) == 1 and len(r[1]) == 1: - print("change in corrdination within a molecule") + return "coordination_change_within_molecule" elif "F" in all_reactant_species: F_ind = [i for i, x in enumerate(all_reactant_species) if x == "F"] - edge = [[x, y] for x in F_ind for y in Li_changing] for e in edge: e.sort() if e in react_global_bonds and e in prod_global_bonds: Li_hopping = True if Li_hopping: - print("LiF hopping") + return "LiF_hopping" # ALiF + B <> A + BLiF else: - print("Li_hopping - inside") + print("Li_hopping") # LiA + B <> A + LiB else: - print("Li_hopping") - + print("Li_hopping") # LiA + B <> A + LiB - - elif (set(Li_ind) & set(r_p)) or (set(Li_ind) & set(p_r)): + elif (set(Li_ind) & set(r_p)) or ( + set(Li_ind) & set(p_r) + ): # some sort of Li bonding OR breaking if single_elemt_Li: - # print((set(Li_ind) & set(r_p)),(set(Li_ind) & set(p_r))) - print("coord and covalent bond break/form") + return "coord_and_covalent_bond_changes" # ex. Li coordination + # causes covalent bond breakage else: - print("for AutoTS") + return "AutoTS" elif single_elem_react_or_prod: - print("single element bond forming or breaking") if single_elemt_Li: - print("Li coordination bond Li+A <> LiA") + return "Li_coordination_Li+A_to_LiA" else: - print("for AutoTS") - else: - print("for AutoTS") - - - - - - all_reactant_species = [] - react_global_bonds = [] - prod_global_bonds = [] - r_charge = 0 - p_charge = 0 - single_elem_reactant = False - - for i in r[0]: - r_charge = r_charge + i.charge - all_reactant_species = all_reactant_species + i.species - if len(i.species) == 1: - single_elem_reactant = True - for i in r[1]: - p_charge = p_charge + i.charge - - if single_elem_reactant: - return "single_element_reactant" - elif r_charge != p_charge: - return "redox" + return "AutoTS" else: - mapping_dict = mrnet.core.reactions.get_reaction_atom_mapping(r[0], r[1]) - for r_id, react in enumerate(r[0]): - for b in react.bonds: - s = (mapping_dict[0][r_id][b[0]], mapping_dict[0][r_id][b[1]]) - s = sorted(s) - react_global_bonds.append(s) - for p_id, prod in enumerate(r[1]): - for b in prod.bonds: - s = (mapping_dict[1][p_id][b[0]], mapping_dict[1][p_id][b[1]]) - s = sorted(s) - prod_global_bonds.append(s) - r_set = set(map(lambda x: repr(x), react_global_bonds)) - p_set = set(map(lambda x: repr(x), prod_global_bonds)) - diff_r_p = list(map(lambda y: literal_eval(y), r_set - p_set)) - diff_p_r = list(map(lambda y: literal_eval(y), p_set - r_set)) - print("p-r", diff_p_r) - print("r-p", diff_r_p) - if diff_r_p != []: - r_p = reduce(operator.concat, diff_r_p) - else: - r_p = [] - if diff_p_r != []: - p_r = reduce(operator.concat, diff_p_r) - else: - p_r = [] - - if r_p == [] or p_r == []: - print("either only forming or breaking bonds during the reaction") - if list(set(diff_r_p[0]).intersection(*diff_r_p)) == []: - print("either only forming or breaking bonds during the reaction AND not following rxn center") - return "non_local_reaction_center" - elif len(list((set(r_p) & set(p_r)))) != 0: - print("reaction center behaving") - Li_ind = [i for i, x in enumerate(all_reactant_species) if x == "Li"] - if Li_ind != []: - if (set(Li_ind) & set(r_p)) and (set(Li_ind) & set(p_r)): - print("in both - so hopping") - return "Li_hopping" - elif (set(Li_ind) & set(r_p)) or (set(Li_ind) & set(p_r)): - print("only in one - so coord") - return "Li_coord_change" - elif "H" in all_reactant_species or "F" in all_reactant_species: - print("actual bond change?") - if len(set(r_p) & set(p_r)) == 1 and all_reactant_species[list(set(r_p) & set(p_r))[0]] == "H" or all_reactant_species[ - list(set(r_p) & set(p_r))[0]] == "F": - print("H or F abstraction") - return "H_or_F_abstraction" - else: - return "actual_bond_changes" - # print(r_p, p_r, list(set(r_p) & set(p_r))[0]) - - - - - + return "AutoTS" From 5e904f7df878d42ad29170e3399681b98f0c8a11 Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Sun, 6 Jun 2021 12:42:00 -0700 Subject: [PATCH 03/17] black --- src/mrnet/utils/visualization.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/mrnet/utils/visualization.py b/src/mrnet/utils/visualization.py index 699a17f3..2b124029 100644 --- a/src/mrnet/utils/visualization.py +++ b/src/mrnet/utils/visualization.py @@ -26,7 +26,14 @@ def visualize_molecule_entry(molecule_entry, path): output the resulting pdf to path """ - atom_colors = {"O": "red", "H": "gray", "C": "black", "Li": "purple", "F": "green", "P": "orange"} + atom_colors = { + "O": "red", + "H": "gray", + "C": "black", + "Li": "purple", + "F": "green", + "P": "orange", + } graph = deepcopy(molecule_entry.graph).to_undirected() From 0b2eeacb848ed598fac289809dea71b0838d6f4a Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Tue, 22 Jun 2021 11:53:49 -0700 Subject: [PATCH 04/17] update --- src/mrnet/utils/reaction_category.py | 33 +++++++++++++++++++++------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 7118c87f..631d5c46 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -4,7 +4,6 @@ from ast import literal_eval, operator from functools import reduce import math - import mrnet from mrnet.stochastic.analyze import SimulationAnalyzer, NetworkUpdater from mrnet.utils.constants import KB, ROOM_TEMP, PLANCK @@ -46,7 +45,7 @@ def reaction_category_RNMC(network_folder, entriesbox, reaction_ids): return category_dict -def update_rates_RNMC(network_folder, category_dict, barrier=None): +def update_rates_RNMC(network_folder, category_dict, barrier_dict=None): """ Method to update rates of reactions based on the type of reactions and its reaction id :param network_folder: path to network folder @@ -57,22 +56,40 @@ def update_rates_RNMC(network_folder, category_dict, barrier=None): value to update to """ update = [] - if barrier is None: - barrier = { + if barrier_dict is None: + barrier_dict = { "Li_hopping": -0.24, - "Li_coord_change": -0.24, } kT = KB * ROOM_TEMP max_rate = kT / PLANCK - for rxn_type in barrier: - rate = max_rate * math.exp(barrier[rxn_type] / kT) - for rxn_id in category_dict[rxn_type]: + rate_dict = {} + for category, barrier in barrier_dict.items(): + rate_dict[category] = max_rate * math.exp(-barrier_dict[barrier] / kT) + for category, rate in rate_dict.items(): + for rxn_id in category_dict[category]: update.append((rxn_id, rate)) network_updater = NetworkUpdater(network_folder) network_updater.update_rates(update) +def update_rates_specific_rxn(network_folder, reactions_barriers): + """ + Method to update rates for specific reactions + :param network_folder: path to network folder + :param reactions_barriers: list of tuples with (reaction_id, barrier) + """ + + update = [] + for r in reactions_barriers: + r_id = r[0] + r_barrier = r[1] + update.append((r_id, r_barrier)) + + network_updater = NetworkUpdater(network_folder) + network_updater.update_rates(update) + + def reaction_category(r): """ Mehtod to categroize a single reaction From 1059d90faba21ba437fa8770af1331d6d839a11e Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Tue, 22 Jun 2021 13:25:54 -0700 Subject: [PATCH 05/17] udpate2 --- src/mrnet/utils/reaction_category.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 631d5c46..075148fd 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -16,7 +16,7 @@ __date__ = "May, 2021" -def reaction_category_RNMC(network_folder, entriesbox, reaction_ids): +def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids): """ Method to categorize reactions from RNMC based on reaction ids :param network_folder: path to network folder @@ -27,9 +27,8 @@ def reaction_category_RNMC(network_folder, entriesbox, reaction_ids): """ category_dict = {} - sa = SimulationAnalyzer(network_folder, entriesbox) for rxn_id in reaction_ids: - r_info = sa.index_to_reaction(rxn_id) + r_info = simulation_analyzer.index_to_reaction(rxn_id) reactants = [] products = [] for r in r_info["reactants"]: @@ -45,7 +44,7 @@ def reaction_category_RNMC(network_folder, entriesbox, reaction_ids): return category_dict -def update_rates_RNMC(network_folder, category_dict, barrier_dict=None): +def update_rates_RNMC(network_folder_to_udpate, category_dict, barrier_dict=None): """ Method to update rates of reactions based on the type of reactions and its reaction id :param network_folder: path to network folder @@ -69,11 +68,11 @@ def update_rates_RNMC(network_folder, category_dict, barrier_dict=None): for rxn_id in category_dict[category]: update.append((rxn_id, rate)) - network_updater = NetworkUpdater(network_folder) + network_updater = NetworkUpdater(network_folder_to_udpate) network_updater.update_rates(update) -def update_rates_specific_rxn(network_folder, reactions_barriers): +def update_rates_specific_rxn(network_folder_to_update, reactions_barriers): """ Method to update rates for specific reactions :param network_folder: path to network folder @@ -86,7 +85,7 @@ def update_rates_specific_rxn(network_folder, reactions_barriers): r_barrier = r[1] update.append((r_id, r_barrier)) - network_updater = NetworkUpdater(network_folder) + network_updater = NetworkUpdater(network_folder_to_update) network_updater.update_rates(update) From cbf6295cb71f0f1fe7392d0a73f94358ee0dd93b Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Tue, 22 Jun 2021 14:20:45 -0700 Subject: [PATCH 06/17] udpate3 --- src/mrnet/utils/reaction_category.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 075148fd..566810df 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -1,8 +1,9 @@ # coding: utf-8 # Copyright (c) Pymatgen Development Team. # Distributed under the terms of the MIT License. -from ast import literal_eval, operator +from ast import literal_eval from functools import reduce +import operator import math import mrnet from mrnet.stochastic.analyze import SimulationAnalyzer, NetworkUpdater @@ -187,9 +188,9 @@ def reaction_category(r): elif LiF: return "LiF_coordinating" # A + LiF <> ALiF else: - print("AutoTS") + return("AutoTS") else: - print( + return( "coord_and_covalent_bond_changes" ) # ex. Li coordination causes covalent bond breakage @@ -213,9 +214,9 @@ def reaction_category(r): if Li_hopping: return "LiF_hopping" # ALiF + B <> A + BLiF else: - print("Li_hopping") # LiA + B <> A + LiB + return("Li_hopping") # LiA + B <> A + LiB else: - print("Li_hopping") # LiA + B <> A + LiB + return("Li_hopping") # LiA + B <> A + LiB elif (set(Li_ind) & set(r_p)) or ( set(Li_ind) & set(p_r) From 882a89d04bf5113514f6e60f98618ccf88e6d36f Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Tue, 22 Jun 2021 16:25:32 -0700 Subject: [PATCH 07/17] update4 --- src/mrnet/utils/reaction_category.py | 38 +++++++++++++++++++++------- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 566810df..7d813b6c 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -39,12 +39,34 @@ def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids): r = [reactants, products] rxn_type = reaction_category(r) if rxn_type not in category_dict.keys(): - category_dict[rxn_type] = [] + category_dict[rxn_type] = [rxn_id] else: category_dict[rxn_type].append(rxn_id) return category_dict +def reaction_extraction_from_pathway( + simulation_analyzer, target_id, num_paths=100, sort_by="weight" +): + + reaction_dict = simulation_analyzer.extract_reaction_pathways(target_id)[target_id] + paths = list(reaction_dict.keys()) + if sort_by == "weight": + paths_sorted = sorted( + paths, + key=lambda x: (reaction_dict[x]["weight"], reaction_dict[x]["frequency"]), + ) + else: + paths_sorted = sorted( + paths, + key=lambda x: (reaction_dict[x]["frequency"], reaction_dict[x]["weight"]), + ) + paths_sorted = [list(x) for x in paths_sorted] + top_paths_sorted = [list(x) for x in paths_sorted][0:num_paths] + all_reactions = reduce(operator.concat, top_paths_sorted) + return all_reactions + + def update_rates_RNMC(network_folder_to_udpate, category_dict, barrier_dict=None): """ Method to update rates of reactions based on the type of reactions and its reaction id @@ -58,13 +80,13 @@ def update_rates_RNMC(network_folder_to_udpate, category_dict, barrier_dict=None update = [] if barrier_dict is None: barrier_dict = { - "Li_hopping": -0.24, + "Li_hopping": 0.24, } kT = KB * ROOM_TEMP max_rate = kT / PLANCK rate_dict = {} for category, barrier in barrier_dict.items(): - rate_dict[category] = max_rate * math.exp(-barrier_dict[barrier] / kT) + rate_dict[category] = max_rate * math.exp(-barrier / kT) for category, rate in rate_dict.items(): for rxn_id in category_dict[category]: update.append((rxn_id, rate)) @@ -188,11 +210,9 @@ def reaction_category(r): elif LiF: return "LiF_coordinating" # A + LiF <> ALiF else: - return("AutoTS") + return "AutoTS" else: - return( - "coord_and_covalent_bond_changes" - ) # ex. Li coordination causes covalent bond breakage + return "coord_and_covalent_bond_changes" # ex. Li coordination causes covalent bond breakage else: # bonds are being broken AND formed if len(list((set(r_p) & set(p_r)))) == 0: # check for reaction center @@ -214,9 +234,9 @@ def reaction_category(r): if Li_hopping: return "LiF_hopping" # ALiF + B <> A + BLiF else: - return("Li_hopping") # LiA + B <> A + LiB + return "Li_hopping" # LiA + B <> A + LiB else: - return("Li_hopping") # LiA + B <> A + LiB + return "Li_hopping" # LiA + B <> A + LiB elif (set(Li_ind) & set(r_p)) or ( set(Li_ind) & set(p_r) From bd0924df023993b1a6d94ea38dcf7a2f6a2dd419 Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Wed, 23 Jun 2021 11:52:27 -0700 Subject: [PATCH 08/17] udpate5 --- src/mrnet/utils/reaction_category.py | 40 +++++++++++++++++----------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 7d813b6c..6de688ad 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -49,22 +49,32 @@ def reaction_extraction_from_pathway( simulation_analyzer, target_id, num_paths=100, sort_by="weight" ): - reaction_dict = simulation_analyzer.extract_reaction_pathways(target_id)[target_id] - paths = list(reaction_dict.keys()) - if sort_by == "weight": - paths_sorted = sorted( - paths, - key=lambda x: (reaction_dict[x]["weight"], reaction_dict[x]["frequency"]), - ) + simulation_analyzer.extract_reaction_pathways(target_id) + reaction_dict = simulation_analyzer.reaction_pathways_dict[target_id] + if reaction_dict == {}: + return "No paths to " + str(target_id) else: - paths_sorted = sorted( - paths, - key=lambda x: (reaction_dict[x]["frequency"], reaction_dict[x]["weight"]), - ) - paths_sorted = [list(x) for x in paths_sorted] - top_paths_sorted = [list(x) for x in paths_sorted][0:num_paths] - all_reactions = reduce(operator.concat, top_paths_sorted) - return all_reactions + paths = list(reaction_dict.keys()) + if sort_by == "weight": + paths_sorted = sorted( + paths, + key=lambda x: ( + reaction_dict[x]["weight"], + reaction_dict[x]["frequency"], + ), + ) + else: + paths_sorted = sorted( + paths, + key=lambda x: ( + reaction_dict[x]["frequency"], + reaction_dict[x]["weight"], + ), + ) + paths_sorted = [list(x) for x in paths_sorted] + top_paths_sorted = [list(x) for x in paths_sorted][0:num_paths] + all_reactions = reduce(operator.concat, top_paths_sorted) + return all_reactions def update_rates_RNMC(network_folder_to_udpate, category_dict, barrier_dict=None): From 1b34f39f9317b1961283f2c2741575a0dd21922b Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Wed, 23 Jun 2021 12:26:48 -0700 Subject: [PATCH 09/17] dev --- src/mrnet/utils/reaction_category.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 6de688ad..259a9612 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -17,17 +17,17 @@ __date__ = "May, 2021" -def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids): +def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids, category_dict={}): """ Method to categorize reactions from RNMC based on reaction ids :param network_folder: path to network folder :param entriesbox: EntriesBox instance of the entries :param reaction_ids: list of RNMC reaction ids + :param category_dict: if pre-existing dictionary with categorization exist :return: category_dict: dict with key being the type of reactions and item being a list of reaction ids """ - category_dict = {} for rxn_id in reaction_ids: r_info = simulation_analyzer.index_to_reaction(rxn_id) reactants = [] @@ -42,6 +42,10 @@ def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids): category_dict[rxn_type] = [rxn_id] else: category_dict[rxn_type].append(rxn_id) + + for rxn_type in category_dict: + category_dict[rxn_type] = list(set(category_dict[rxn_type])) + return category_dict @@ -73,7 +77,8 @@ def reaction_extraction_from_pathway( ) paths_sorted = [list(x) for x in paths_sorted] top_paths_sorted = [list(x) for x in paths_sorted][0:num_paths] - all_reactions = reduce(operator.concat, top_paths_sorted) + all_reactions = list(set(reduce(operator.concat, top_paths_sorted))) + return all_reactions From 1ee0527021a965cd30030a7e6b0d64cef67e644f Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Wed, 23 Jun 2021 13:00:40 -0700 Subject: [PATCH 10/17] dev --- src/mrnet/utils/reaction_category.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 259a9612..f8192f07 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -160,7 +160,6 @@ def reaction_category(r): if i.formula == "F1 Li1": LiF = True elif len(i.species) == 1: - print(i.species) single_elem_react_or_prod = True if i.formula == "Li1": single_elemt_Li = True From 3e0561cb79244debd17253044bd072c176942eab Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Wed, 23 Jun 2021 14:15:27 -0700 Subject: [PATCH 11/17] dev --- src/mrnet/utils/reaction_category.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index f8192f07..73d0a9c9 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -8,6 +8,7 @@ import mrnet from mrnet.stochastic.analyze import SimulationAnalyzer, NetworkUpdater from mrnet.utils.constants import KB, ROOM_TEMP, PLANCK +from mrnet.utils.visualization import generate_latex_header, generate_latex_footer __author__ = "Hetal D. Patel" __maintainer__ = "Hetal D. Patel" @@ -81,6 +82,23 @@ def reaction_extraction_from_pathway( return all_reactions +def generate_categorization_report(sa, reports_folder, category_dict, categories_to_print=[]): + + + with open( + reports_folder + "/categorization.tex", + "w", + ) as f: + generate_latex_header(f) + for rxn_type, reactions in category_dict.items(): + f.write("\n\n\n") + f.write("REACTION TYPE COUNT: " + rxn_type + " " + str(len(reactions))) + if rxn_type in categories_to_print: + f.write("\n\n\n") + for reaction_index in reactions: + sa.latex_emit_reaction(f, reaction_index) + generate_latex_footer(f) + def update_rates_RNMC(network_folder_to_udpate, category_dict, barrier_dict=None): """ From f0b00a3bc753aea062b9f6d6a7ed01cd9b30d237 Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Wed, 23 Jun 2021 15:27:03 -0700 Subject: [PATCH 12/17] bug in corrd_covalent categorization + report fix --- src/mrnet/utils/reaction_category.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 73d0a9c9..1d7ee7ca 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -92,7 +92,8 @@ def generate_categorization_report(sa, reports_folder, category_dict, categories generate_latex_header(f) for rxn_type, reactions in category_dict.items(): f.write("\n\n\n") - f.write("REACTION TYPE COUNT: " + rxn_type + " " + str(len(reactions))) + str_type = str("Reaction Type: " + rxn_type + " count: " + str(len(reactions))).replace("_", " ") + f.write(str_type) if rxn_type in categories_to_print: f.write("\n\n\n") for reaction_index in reactions: @@ -244,7 +245,7 @@ def reaction_category(r): else: return "AutoTS" else: - return "coord_and_covalent_bond_changes" # ex. Li coordination causes covalent bond breakage + return "AutoTS" else: # bonds are being broken AND formed if len(list((set(r_p) & set(p_r)))) == 0: # check for reaction center From 710fc21c842eef221ed25c6f57f36b26b80c0b03 Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Fri, 25 Jun 2021 13:39:24 -0700 Subject: [PATCH 13/17] Li hopping updates+dG based barriers --- src/mrnet/utils/reaction_category.py | 52 ++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 14 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 1d7ee7ca..9cc5fc44 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -49,6 +49,22 @@ def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids, catego return category_dict +def dG_based_barrier_update_RNMC_(simulation_analyzer, network_folder_to_udpate, reaction_ids,barrier=0.24, + abs_cutoff=0.08): + + kT = KB * ROOM_TEMP + max_rate = kT / PLANCK + rate = max_rate * math.exp(-barrier / kT) + update = [] + for rxn_id in reaction_ids: + r_info = simulation_analyzer.index_to_reaction(rxn_id) + dG = r_info['free_energy'] + if abs(dG) <= abs_cutoff: + update.append((rxn_id, rate)) + + network_updater = NetworkUpdater(network_folder_to_udpate) + network_updater.update_rates(update) + def reaction_extraction_from_pathway( simulation_analyzer, target_id, num_paths=100, sort_by="weight" @@ -84,7 +100,8 @@ def reaction_extraction_from_pathway( def generate_categorization_report(sa, reports_folder, category_dict, categories_to_print=[]): - + if categories_to_print == []: + categories_to_print = list(category_dict.keys()) with open( reports_folder + "/categorization.tex", "w", @@ -129,24 +146,22 @@ def update_rates_RNMC(network_folder_to_udpate, category_dict, barrier_dict=None network_updater.update_rates(update) -def update_rates_specific_rxn(network_folder_to_update, reactions_barriers): +def supress_rxns_RNMC(network_folder_to_update, reactions_ids): """ - Method to update rates for specific reactions - :param network_folder: path to network folder - :param reactions_barriers: list of tuples with (reaction_id, barrier) + Method to add high barriers to reactions + :param network_folder_to_update: path to network folder + :param reactions_ids: list of reaction ids """ update = [] - for r in reactions_barriers: - r_id = r[0] - r_barrier = r[1] - update.append((r_id, r_barrier)) + for r_id in reactions_ids: + update.append((r_id, 0)) network_updater = NetworkUpdater(network_folder_to_update) network_updater.update_rates(update) -def reaction_category(r): +def reaction_category(r, initial_mol_molentries=None): """ Mehtod to categroize a single reaction :param r: list of reactant and product MoleculeEntry [[reactnat molecule entries], @@ -175,7 +190,6 @@ def reaction_category(r): single_elemt_Li = True for i in r[1]: p_charge = p_charge + i.charge - if i.formula == "F1 Li1": LiF = True elif len(i.species) == 1: @@ -225,7 +239,7 @@ def reaction_category(r): list(set(combined_diff[0]).intersection(*combined_diff)) == [] ): # reaction center rule not satisfied if LiF: - return "non_local_reaction_center_LiF_formning" + return "non_local_reaction_center_LiF_forming" else: return "non_local_reaction_center" @@ -267,9 +281,19 @@ def reaction_category(r): if Li_hopping: return "LiF_hopping" # ALiF + B <> A + BLiF else: - return "Li_hopping" # LiA + B <> A + LiB + if initial_mol_molentries is not None: + if set(r[0]).intersection(set(initial_mol_molentries)) != set() and set(r[1]).intersection(set( + initial_mol_molentries)) != set(): + return "Li_hopping_initial_mol_in_reactant_and_product" + else: + return "Li_hopping" # LiA + B <> A + LiB else: - return "Li_hopping" # LiA + B <> A + LiB + if initial_mol_molentries is not None: + if set(r[0]).intersection(set(initial_mol_molentries)) != set() and set(r[1]).intersection(set( + initial_mol_molentries)) != set(): + return "Li_hopping_initial_mol_in_reactant_and_product" + else: + return "Li_hopping" # LiA + B <> A + LiB elif (set(Li_ind) & set(r_p)) or ( set(Li_ind) & set(p_r) From 5cf03e84222af95e447ef844478d774ec08badaf Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Fri, 25 Jun 2021 15:39:51 -0700 Subject: [PATCH 14/17] bug --- src/mrnet/utils/reaction_category.py | 48 ++++++++++++++++++++-------- 1 file changed, 35 insertions(+), 13 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 9cc5fc44..d1f391ba 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -18,7 +18,13 @@ __date__ = "May, 2021" -def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids, category_dict={}): +def reaction_category_RNMC( + simulation_analyzer, + entriesbox, + reaction_ids, + category_dict={}, + initial_mol_molentries=None, +): """ Method to categorize reactions from RNMC based on reaction ids :param network_folder: path to network folder @@ -38,7 +44,7 @@ def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids, catego for p in r_info["products"]: products.append(entriesbox.entries_list[p]) r = [reactants, products] - rxn_type = reaction_category(r) + rxn_type = reaction_category(r, initial_mol_molentries) if rxn_type not in category_dict.keys(): category_dict[rxn_type] = [rxn_id] else: @@ -49,8 +55,14 @@ def reaction_category_RNMC(simulation_analyzer, entriesbox, reaction_ids, catego return category_dict -def dG_based_barrier_update_RNMC_(simulation_analyzer, network_folder_to_udpate, reaction_ids,barrier=0.24, - abs_cutoff=0.08): + +def dG_based_barrier_update_RNMC_( + simulation_analyzer, + network_folder_to_udpate, + reaction_ids, + barrier=0.24, + abs_cutoff=0.08, +): kT = KB * ROOM_TEMP max_rate = kT / PLANCK @@ -58,7 +70,7 @@ def dG_based_barrier_update_RNMC_(simulation_analyzer, network_folder_to_udpate, update = [] for rxn_id in reaction_ids: r_info = simulation_analyzer.index_to_reaction(rxn_id) - dG = r_info['free_energy'] + dG = r_info["free_energy"] if abs(dG) <= abs_cutoff: update.append((rxn_id, rate)) @@ -98,18 +110,23 @@ def reaction_extraction_from_pathway( return all_reactions -def generate_categorization_report(sa, reports_folder, category_dict, categories_to_print=[]): + +def generate_categorization_report( + sa, reports_folder, category_dict, categories_to_print=[] +): if categories_to_print == []: categories_to_print = list(category_dict.keys()) with open( - reports_folder + "/categorization.tex", - "w", + reports_folder + "/categorization.tex", + "w", ) as f: generate_latex_header(f) for rxn_type, reactions in category_dict.items(): f.write("\n\n\n") - str_type = str("Reaction Type: " + rxn_type + " count: " + str(len(reactions))).replace("_", " ") + str_type = str( + "Reaction Type: " + rxn_type + " count: " + str(len(reactions)) + ).replace("_", " ") f.write(str_type) if rxn_type in categories_to_print: f.write("\n\n\n") @@ -282,15 +299,20 @@ def reaction_category(r, initial_mol_molentries=None): return "LiF_hopping" # ALiF + B <> A + BLiF else: if initial_mol_molentries is not None: - if set(r[0]).intersection(set(initial_mol_molentries)) != set() and set(r[1]).intersection(set( - initial_mol_molentries)) != set(): + if ( + set(r[0]).intersection(set(initial_mol_molentries)) != set() + and set(r[1]).intersection(set(initial_mol_molentries)) + != set() + ): return "Li_hopping_initial_mol_in_reactant_and_product" else: return "Li_hopping" # LiA + B <> A + LiB else: if initial_mol_molentries is not None: - if set(r[0]).intersection(set(initial_mol_molentries)) != set() and set(r[1]).intersection(set( - initial_mol_molentries)) != set(): + if ( + set(r[0]).intersection(set(initial_mol_molentries)) != set() + and set(r[1]).intersection(set(initial_mol_molentries)) != set() + ): return "Li_hopping_initial_mol_in_reactant_and_product" else: return "Li_hopping" # LiA + B <> A + LiB From a68f338f43d5b27dd01297348f3f761f90c037d7 Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Fri, 25 Jun 2021 16:34:49 -0700 Subject: [PATCH 15/17] bug --- src/mrnet/utils/reaction_category.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index d1f391ba..7b5f4473 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -56,9 +56,9 @@ def reaction_category_RNMC( return category_dict -def dG_based_barrier_update_RNMC_( +def dG_based_barrier_update_RNMC( simulation_analyzer, - network_folder_to_udpate, + network_folder_to_update, reaction_ids, barrier=0.24, abs_cutoff=0.08, @@ -74,8 +74,8 @@ def dG_based_barrier_update_RNMC_( if abs(dG) <= abs_cutoff: update.append((rxn_id, rate)) - network_updater = NetworkUpdater(network_folder_to_udpate) - network_updater.update_rates(update) + network_updater = NetworkUpdater(network_folder_to_update) + network_updater.update_rates(update) def reaction_extraction_from_pathway( From 33c96df315fbcdf0c35e706aed03fc3add7957ba Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Fri, 25 Jun 2021 17:03:31 -0700 Subject: [PATCH 16/17] li hopping bug --- src/mrnet/utils/reaction_category.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 7b5f4473..9b9bd307 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -305,6 +305,8 @@ def reaction_category(r, initial_mol_molentries=None): != set() ): return "Li_hopping_initial_mol_in_reactant_and_product" + else: + return "Li_hopping" # LiA + B <> A + LiB else: return "Li_hopping" # LiA + B <> A + LiB else: @@ -314,6 +316,8 @@ def reaction_category(r, initial_mol_molentries=None): and set(r[1]).intersection(set(initial_mol_molentries)) != set() ): return "Li_hopping_initial_mol_in_reactant_and_product" + else: + return "Li_hopping" # LiA + B <> A + LiB else: return "Li_hopping" # LiA + B <> A + LiB From 148eee5d40e3a40342d8e5c3919550ed1f5b039e Mon Sep 17 00:00:00 2001 From: Hetal Patel Date: Thu, 1 Jul 2021 16:32:16 -0700 Subject: [PATCH 17/17] cleaned up --- src/mrnet/utils/reaction_category.py | 33 ++++++++++++++++++++++++---- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/src/mrnet/utils/reaction_category.py b/src/mrnet/utils/reaction_category.py index 9b9bd307..e7b8a8b7 100644 --- a/src/mrnet/utils/reaction_category.py +++ b/src/mrnet/utils/reaction_category.py @@ -31,7 +31,7 @@ def reaction_category_RNMC( :param entriesbox: EntriesBox instance of the entries :param reaction_ids: list of RNMC reaction ids :param category_dict: if pre-existing dictionary with categorization exist - :return: category_dict: dict with key being the type of reactions and item being + :return: category_dict: dictionary with key being the type of reactions and item being a list of reaction ids """ @@ -63,6 +63,14 @@ def dG_based_barrier_update_RNMC( barrier=0.24, abs_cutoff=0.08, ): + """ + Method to identify Li_hopping reactions with a dG within the absolate value of the cutoff. + :param simulation_analyzer: SimulationAnalyzer object + :param network_folder_to_update: path to network folder with rn.sqlite to update + :param reaction_ids: list of Li_hopping reaction ids + :param barrier: new barrier to set to the Li_hopping reactions with dG within the absolate value of the cuttoff + :param abs_cutoff: absolate dG value within which the barrier should be modified + """ kT = KB * ROOM_TEMP max_rate = kT / PLANCK @@ -70,7 +78,7 @@ def dG_based_barrier_update_RNMC( update = [] for rxn_id in reaction_ids: r_info = simulation_analyzer.index_to_reaction(rxn_id) - dG = r_info["free_energy"] + dG = r_info["dG"] if abs(dG) <= abs_cutoff: update.append((rxn_id, rate)) @@ -81,6 +89,14 @@ def dG_based_barrier_update_RNMC( def reaction_extraction_from_pathway( simulation_analyzer, target_id, num_paths=100, sort_by="weight" ): + """ + Method to extract reaction ids from the pathways + :param simulation_analyzer: SimulationAnalyzer object + :param target_id: speices id to path find to + :param num_paths: number of paths to extract reactions from + :param sort_by: sort top pathways by frequency or weight + :return: all_reactions: list of reaction ids + """ simulation_analyzer.extract_reaction_pathways(target_id) reaction_dict = simulation_analyzer.reaction_pathways_dict[target_id] @@ -114,6 +130,13 @@ def reaction_extraction_from_pathway( def generate_categorization_report( sa, reports_folder, category_dict, categories_to_print=[] ): + """ + Method to generate categorization report + :param sa: SimulationAnalyzer object + :param reports_folder: path to the save the report to + :param category_dict: dictionary with reaction category as the key and list of reactions as the value + :param categories_to_print: list of reaction types to print, if "[]" all reaction types will be printed + """ if categories_to_print == []: categories_to_print = list(category_dict.keys()) @@ -138,10 +161,10 @@ def generate_categorization_report( def update_rates_RNMC(network_folder_to_udpate, category_dict, barrier_dict=None): """ Method to update rates of reactions based on the type of reactions and its reaction id - :param network_folder: path to network folder + :param network_folder_to_udpate: path to network folder :param category_dict: category_dict: dict with key being the type of reactions and item being a list of reaction ids - :param rates: dict with key being the category of reactions to udpate barriers for + :param barrier_dict: dict with key being the category of reactions to udpate barriers for and the value being the barrier value to update to """ @@ -183,6 +206,8 @@ def reaction_category(r, initial_mol_molentries=None): Mehtod to categroize a single reaction :param r: list of reactant and product MoleculeEntry [[reactnat molecule entries], [product moleucle entries]] + :param initial_mol_molentries: list of inital molecule ids, if None Li_hopping reactions with initial + molecules both in reactant and prodcuts will not be identified as such. :return: string indicating type of reaction """