Bi nitrite and 2ME CRN

HiPRGen/gen_euvl_testset.py

@@ -0,0 +1,72 @@
+print("Importing required files...")
+from pymatgen.core.structure import Molecule
+from pymatgen.analysis.graphs import MoleculeGraph
+from pymatgen.analysis.local_env import OpenBabelNN
+from monty.serialization import loadfn, dumpfn
+from mol_entry import MoleculeEntry
+from initial_state import find_mol_entry_by_entry_id
+import copy
+print("Done!")
+
+#list of mol_ids for species you want to find with each mol_id occupying its own line. Must have an empty line
+#at the end of the file, but empty lines will not be copied.
+print("Opening and reading species id file...")
+with open('atom_mapping_set.txt') as mol_id_doc: #open the text file and copy each line as the element of a list, removing '\n'
+    mol_id_list = []
+    for line in mol_id_doc:
+        if line.strip(): #excludes empty lines
+            mol_id_list.append(line[0:len(line)-1])
+print("Done!")
+
+print("Opening json file...")
+mol_json = "032923.json"
+database_entries = loadfn(mol_json) #loads the json file to a list of dictionaries
+print("Done!")
+
+print("Reading molecule entries from json file...") #converts dictionary from loaded json file to a list
+mol_entries= [MoleculeEntry.from_mp_doc(e) for e in database_entries]
+print("Done!")
+
+final_mol_ids = []
+
+print('Copying desired molecules to a new list...')
+for mol_id in mol_id_list: #takes the list of ids, finds their corresponding molecules in the json file 
+    for mol_entry in mol_entries:
+        m_id = mol_entry.entry_id
+        if m_id == mol_id:
+            graph = mol_entry.mol_graph
+            for entry in mol_entries: #(as well as differently charged variants) and adds them to the list
+                molecule_graph = entry.mol_graph #tests if other charges found
+                if molecule_graph.isomorphic_to(graph):
+                    final_mol_ids.append(entry.entry_id)
+
+to_remove = ['nbo', 'alpha', 'beta']
+
+test_set = []
+for entry in database_entries:
+    if entry['molecule_id'] in final_mol_ids:
+        new_entry = copy.deepcopy(entry)
+        for name in entry.keys():
+            for remove in to_remove:
+                if remove in name:
+                    val = new_entry.pop(name)
+        test_set.append(new_entry)
+
+dumpfn(test_set, 'atom_mapping_set.json') 
+
+new_mol_json = loadfn('atom_mapping_set.json')
+new_mol_entries= [MoleculeEntry.from_mp_doc(e) for e in new_mol_json]
+
+new_id_list = []
+
+for molecule in new_mol_entries:
+    new_id_list.append(molecule.entry_id)
+
+new_id_set = set(new_id_list)
+mol_id_set = set(mol_id_list)
+differences = mol_id_set - new_id_set
+if differences == set():
+    print('All desired entries successfully copied')
+else:
+    for elem in differences:
+        print('entry ', elem, ' not copied successfully')

HiPRGen/initial_state.py

@@ -1,7 +1,6 @@
 from pymatgen.core.structure import Molecule
 from pymatgen.analysis.graphs import MoleculeGraph
-from pymatgen.analysis.local_env import OpenBabelNN
-from pymatgen.analysis.fragmenter import metal_edge_extender
+from pymatgen.analysis.local_env import OpenBabelNN, metal_edge_extender#, oxygen_edge_extender
 import sqlite3
 
 
@@ -15,7 +14,8 @@ def find_mol_entry_from_xyz_and_charge(mol_entries, xyz_file_path, charge):
     )
 
     # correction to the molecule graph
-    target_mol_graph = metal_edge_extender(target_mol_graph)
+    target_mol_graph = metal_edge_extender(target_mol_graph, cutoff=3.4, metals=None,coordinators=("O"))
+    #target_mol_graph = metal_edge_extender(target_mol_graph)
 
     match = False
     index = -1
@@ -71,8 +71,8 @@ def insert_initial_state(
         mol_entries,
         initial_state_db,
         factor_zero = 1.0,
-        factor_two = 1.0,
-        factor_duplicate = 0.5
+        factor_two = 0.5,
+        factor_duplicate = 1.0,
 ):
     """
     initial state is a dict mapping species ids to counts.

HiPRGen/mc_analysis.py

@@ -527,7 +527,7 @@ def reaction_tally_report(
     for (reaction_index, number) in sorted(
             reaction_tally.items(), key=lambda pair: -pair[1]):
         if number > cutoff:
-            report_generator.emit_text(str(number) + " occourances of:")
+            report_generator.emit_text(str(number) + " occurrences of:")
             report_generator.emit_reaction(
                 network_loader.index_to_reaction(reaction_index),
                 label=str(reaction_index)
@@ -551,7 +551,7 @@ def species_report(network_loader, species_report_path):
         report_generator.emit_text(str(mol.entry_id))
         report_generator.emit_text(
             "formula: " + mol.formula)
-
+        report_generator.emit_text("free energy = " + str(mol.free_energy) + " eV")
         report_generator.emit_molecule(i)
         report_generator.emit_newline()
 
@@ -750,7 +750,7 @@ def export_pathways_to_json(pathfinding, species_id, path):
 
 
 
-def generate_pathway_report(
+def pathway_report(
         pathfinding,
         species_id,
         report_file_path,
@@ -839,7 +839,7 @@ def compute_expected_final_state(self):
         )
 
         for seed in self.network_loader.trajectories:
-            state = np.copy(self.network_loader.initial_state_array)
+            state = np.copy(self.network_loader.initial_state_array[seed])
             for step in self.network_loader.trajectories[seed]:
                 reaction_index = self.network_loader.trajectories[seed][step][0]
                 time = self.network_loader.trajectories[seed][step][1]
@@ -858,6 +858,24 @@ def compute_expected_final_state(self):
         self.expected_final_state = (
             self.expected_final_state / len(self.network_loader.trajectories))
 
+    def compute_trajectory_final_states(self):
+        self.final_states = {}
+        for seed in self.network_loader.trajectories:
+            state = np.copy(self.network_loader.initial_state_array[seed])
+            for step in self.network_loader.trajectories[seed]:
+                reaction_index = self.network_loader.trajectories[seed][step][0]
+                time = self.network_loader.trajectories[seed][step][1]
+                reaction = self.network_loader.index_to_reaction(reaction_index)
+
+                for i in range(reaction['number_of_reactants']):
+                    reactant_index = reaction['reactants'][i]
+                    state[reactant_index] -= 1
+
+                for j in range(reaction['number_of_products']):
+                    product_index = reaction['products'][j]
+                    state[product_index] += 1
+            self.final_states[seed] = state
+
     def compute_production_consumption_info(self):
         self.consuming_reactions = {}
         self.producing_reactions = {}
@@ -888,13 +906,13 @@ def compute_production_consumption_info(self):
                         self.producing_reactions[product_index][reaction_index] += 1
 
     def compute_state_time_series(self, seed):
-        state_dimension_size = len(self.network_loader.initial_state_array)
+        state_dimension_size = len(self.network_loader.initial_state_array[seed])
         step_dimension_size = len(self.network_loader.trajectories[seed])
         time_series = np.zeros(
             (step_dimension_size, state_dimension_size),
             dtype=int)
 
-        state = np.copy(self.network_loader.initial_state_array)
+        state = np.copy(self.network_loader.initial_state_array[seed])
         for step in self.network_loader.trajectories[seed]:
             reaction_index = self.network_loader.trajectories[seed][step][0]
             time = self.network_loader.trajectories[seed][step][1]
@@ -918,9 +936,11 @@ def time_series_graph(
             seeds,
             species_of_interest,
             path,
+            custom_y_max=None,
+            custom_colorstyle_list=None,
             colors = list(mcolors.TABLEAU_COLORS.values()),
-            styles = ['solid', 'dotted', 'dashed', 'dashdot'],
-            internal_index_labels=True
+            styles = ['solid', 'dotted', 'dashed', 'dashdot','solid', 'dotted', 'dashed', 'dashdot','solid', 'dotted', 'dashed', 'dashdot','solid', 'dotted', 'dashed', 'dashdot'],
+            internal_index_labels=True,
     ):
 
 
@@ -952,11 +972,16 @@ def time_series_graph(
 
         line_dict = {}
         i = 0
-        for species_index in species_of_interest:
-            r = i % len(colors)
-            q = i // len(colors)
-            line_dict[species_index] = (colors[r], styles[q])
-            i += 1
+        if custom_colorstyle_list is None:
+            for species_index in species_of_interest:
+                r = i % len(colors)
+                q = i // len(colors)
+                line_dict[species_index] = (colors[r], styles[q])
+                i += 1
+        else:
+            for species_index in species_of_interest:
+                line_dict[species_index] = (custom_colorstyle_list[i][0], custom_colorstyle_list[i][1])
+                i += 1
 
 
         fig, (ax0, ax1, ax2) = plt.subplots(
@@ -965,15 +990,19 @@ def time_series_graph(
             gridspec_kw={'height_ratios':[2,2,1]})
 
         y_max = 0
-        for step in range(total_time_series.shape[0]):
-            for species_index in species_of_interest:
-                y_max = max(y_max, total_time_series[step,species_index])
+        if custom_y_max is None:
+            for step in range(total_time_series.shape[0]):
+                for species_index in species_of_interest:
+                    y_max = max(y_max, total_time_series[step,species_index])
+            y_max += 1
+        else:
+            y_max = custom_y_max
 
         ax0.set_xlim([0,total_time_series.shape[0]])
-        ax0.set_ylim([0,y_max+1])
+        ax0.set_ylim([0,y_max])
 
         ax1.set_xlim([0,total_time_series.shape[0]])
-        ax1.set_ylim([0,(y_max+1)/10])
+        ax1.set_ylim([0,(y_max)/10])
 
 
         ticks = np.arange(0, total_time_series.shape[0])
@@ -1097,8 +1126,8 @@ def sink_filter(self, species_index):
         mol = self.network_loader.mol_entries[species_index]
         if (number_of_consuming_reactions + number_of_producing_reactions > 0  and
             ratio > 1.5 and
-            expected_value > 0.1 and
-            mol.spin_multiplicity == 1):
+            expected_value > 0.1):# and
+            #mol.spin_multiplicity == 1):
             return True
         else:
             return False
@@ -1121,11 +1150,11 @@ def export_consumption_to_json(simulation_replayer, species_index, path):
             }
         reactions[reaction_id] = json_reaction
 
-        dumpfn({
-            'reactions' : reactions,
-            'producing_reactions' : producing_reactions,
-            'consuming_reactions' : consuming_reactions},
-               path)
+    dumpfn({
+        'reactions' : reactions,
+        'producing_reactions' : producing_reactions,
+        'consuming_reactions' : consuming_reactions},
+           path)
 
 def pad_time_series(time_series, max_number_of_steps):
     num_steps = time_series.shape[0]
@@ -1144,15 +1173,22 @@ def pad_time_series(time_series, max_number_of_steps):
 
 
 
-def export_sinks_to_json(simulation_replayer, path):
-    sink_data = simulation_replayer.sink_data
+def export_full_sink_data_to_json(simulation_replayer, path):
     sink_data_json = {}
-    for i in sink_data:
+    for i in simulation_replayer.sink_data: # All molecules formed at least once are in sink_data
         mol = simulation_replayer.network_loader.mol_entries[i]
-        sink_data_json[mol.entry_id] = sink_data[i]
+        sink_data_json[mol.entry_id] = simulation_replayer.sink_data[i]
 
     dumpfn(sink_data_json, path)
 
+def export_sinks_to_json(simulation_replayer, path):
+    sinks_json = {}
+    for i in simulation_replayer.sinks: # Only molecules that pass sink_filter are in sinks
+        mol = simulation_replayer.network_loader.mol_entries[i]
+        sinks_json[mol.entry_id] = simulation_replayer.sink_data[i]
+
+    dumpfn(sinks_json, path)
+
 def export_species_report_to_json(network_loader, path):
     data = {}
     for i in range(network_loader.number_of_species):
@@ -1260,3 +1296,26 @@ def sink_report(
         report_generator.emit_newline()
 
     report_generator.finished()
+
+
+def final_state_report(
+        simulation_replayer,
+        final_state_report_path
+):
+    final_state = []
+    for species_index, value in enumerate(simulation_replayer.expected_final_state):
+        if value > 0.0:
+            final_state.append({"index": species_index, "value": value})
+    sorted_final_state = sorted(final_state, key=lambda x: -x["value"])
+
+    report_generator = ReportGenerator(
+        simulation_replayer.network_loader.mol_entries,
+        final_state_report_path,
+        rebuild_mol_pictures=False)
+
+    for entry in sorted_final_state:
+        report_generator.emit_text("amount: " + str(entry["value"]))
+        report_generator.emit_molecule(entry["index"])
+        report_generator.emit_newline()
+    report_generator.finished()
+