From b972628a64ba54f110b271d2646c70cb622ff587 Mon Sep 17 00:00:00 2001
From: donerancl <donera@mit.edu>
Date: Mon, 10 Feb 2025 13:54:36 -0500
Subject: [PATCH 1/6] fixed functions in fragment_utils.py

---
 rmgpy/molecule/fragment_utils.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/rmgpy/molecule/fragment_utils.py b/rmgpy/molecule/fragment_utils.py
index 25fe5e5cb5..201f63ec5f 100644
--- a/rmgpy/molecule/fragment_utils.py
+++ b/rmgpy/molecule/fragment_utils.py
@@ -98,7 +98,7 @@ def match_concentrations_with_same_sums(conc1, conc2, rtol=1e-6):
     seq1 = [tup[1] for tup in conc1]
     seq2 = [tup[1] for tup in conc2]
 
-    matches_seq = FragList.match_sequences(seq1, seq2, rtol)
+    matches_seq = match_sequences(seq1, seq2, rtol)
 
     matches_conc = []
     for match_seq in matches_seq:
@@ -184,7 +184,7 @@ def match_concentrations_with_different_sums(conc1, conc2):
     # let matches_conc match with remaining seq2
     elif pin1 == len(seq1) and pin2 < len(seq2):
         remain_conc2 = [(labels2[pin2], val2)] + conc2[(pin2 + 1):]
-        matches_conc = FragList.match_concentrations_with_different_sums(
+        matches_conc = match_concentrations_with_different_sums(
             matches_conc, remain_conc2
         )
 
@@ -216,7 +216,7 @@ def flatten(combo):
     return_list = []
     for i in combo:
         if isinstance(i, tuple):
-            return_list.extend(FragList.flatten(i))
+            return_list.extend(flatten(i))
         else:
             return_list.append(i)
     return return_list
@@ -278,10 +278,10 @@ def merge_frag_list(to_be_merged):
         frag2 = to_be_merged[-1].smiles  # last fragment in list
 
         if 'R' in frag1 and 'L' in frag2:
-            newfrag = FragList.merge_frag_to_frag(frag1, frag2, 'L')
+            newfrag = merge_frag_to_frag(frag1, frag2, 'L')
 
         elif 'L' in frag1 and 'R' in frag2:
-            newfrag = FragList.merge_frag_to_frag(frag1, frag2, 'R')
+            newfrag = merge_frag_to_frag(frag1, frag2, 'R')
 
         # warn user if last two fragments in list cannot be merged (no R/L
         # combo to be made)
@@ -290,7 +290,7 @@ def merge_frag_list(to_be_merged):
                 frag1, frag2))
 
             if 'L' in frag1 and 'L' in frag2:
-                newfrag = FragList.merge_frag_to_frag(
+                newfrag = merge_frag_to_frag(
                     frag1.replace('L', 'R'), frag2, 'L')
         if len(to_be_merged) > 2:
             cut = len(to_be_merged) - 2

From 650e116e8644abfc21baf09cb4ede30169c83297 Mon Sep 17 00:00:00 2001
From: donerancl <donera@mit.edu>
Date: Mon, 2 Jun 2025 12:54:41 -0400
Subject: [PATCH 2/6] moved ammonia project PR changes to new branch after
 update of main to python 3.9

---
 rmgpy/data/kinetics/family.py | 43 +++++++++++++++++++++++++++++++++++
 rmgpy/molecule/molecule.py    | 15 ++++++++++++
 rmgpy/molecule/util.py        | 39 +++++++++++++++++++++++++++++++
 rmgpy/rmg/model.py            |  4 +++-
 rmgpy/rmg/pdep.py             | 12 ++--------
 5 files changed, 102 insertions(+), 11 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 102c5f13ff..d949eaf48d 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -214,6 +214,32 @@ def copy(self):
 
         return other
 
+    def check_if_spin_allowed(self):
+        # get the combined spin for reactants and products
+        reactants_combined_spin, products_combined_spin = self.calculate_combined_spin()
+        # check if there are any matches for combined spin between reactants and products
+        if reactants_combined_spin.intersection(products_combined_spin) != set([]):
+            return True
+        else:
+            logging.debug(f"Reactants combined spin is {reactants_combined_spin}, but the products combined spin is {products_combined_spin}")
+            return False
+        
+    def calculate_combined_spin(self):
+        if len(self.reactants) == 1:
+            reactant_combined_spin = {self.reactants[0].multiplicity}
+        elif len(self.reactants) == 2:
+            reactant_spin_string =  "+".join(sorted([str(reactant.multiplicity) for reactant in self.reactants]))
+            reactant_combined_spin = allowed_spin[reactant_spin_string]
+        else:
+            return None
+        if len(self.products) == 1:
+            product_combined_spin = {self.products[0].multiplicity}
+        elif len(self.products) == 2:
+            product_spin_string = "+".join(sorted([str(product.multiplicity) for product in self.products]))
+            product_combined_spin = allowed_spin[product_spin_string]
+        else:
+            return None
+        return reactant_combined_spin, product_combined_spin
     def apply_solvent_correction(self, solvent):
         """
         apply kinetic solvent correction in this case the parameters are dGTSsite instead of GTS
@@ -4902,3 +4928,20 @@ def get_site_solute_data(rxn):
         return site_data
     else:
         return None
+
+allowed_spin_violation_families =['1,4_Cyclic_birad_scission']
+
+allowed_spin = {
+    "1+1": set([1]),
+    "1+2": set([2]),
+    "1+3": set([3]),
+    "1+4": set([4]),
+    "1+5": set([5]),
+    "2+2": set([1,3]),
+    "2+3": set([2,4]),
+    "2+4": set([3,5]),
+    "2+5": set([4,6]),
+    "3+3": set([1,3,5]),
+    "3+4": set([2,4,6]),
+    "3+5": set([3,5,7]),
+}
\ No newline at end of file
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 1a32f07add..1cd0f3e515 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -60,6 +60,7 @@
 from rmgpy.molecule.kekulize import kekulize
 from rmgpy.molecule.pathfinder import find_shortest_path
 from rmgpy.molecule.fragment import CuttingLabel
+from rmgpy.molecule.util import generate_closed_shell_singlet, generate_singlet_diradicals
 
 ################################################################################
 
@@ -3002,6 +3003,20 @@ def get_desorbed_molecules(self):
 
         return desorbed_molecules
 
+    def update_to_closed_shell_singlet(self):Add commentMore actions
+        for i in range(len(self.atoms)):
+                self.atoms[i].id = i
+        assert self.multiplicity == 1 and self.get_radical_count()>0
+        radical_center_ids = [x.id for x in self.atoms if x.radical_electrons > 0]
+        # remove radicals from radical centers (2)
+        for radical_center_id in radical_center_ids:
+            self.atoms[radical_center_id].decrement_radical()
+        # add removed radicals (2) to one of the radical sites as a lone pair (1)
+        self.atoms[radical_center_ids[0]].increment_lone_pairs()
+        # pick the best resonance structure
+        self.generate_resonance_structures()
+        self.reactive = True
+
 # this variable is used to name atom IDs so that there are as few conflicts by 
 # using the entire space of integer objects
 atom_id_counter = -2 ** 15
diff --git a/rmgpy/molecule/util.py b/rmgpy/molecule/util.py
index 46f8c3504c..0e3dc04171 100644
--- a/rmgpy/molecule/util.py
+++ b/rmgpy/molecule/util.py
@@ -179,3 +179,42 @@ def swap(to_be_swapped, sample):
     new_partner = (to_be_swapped - sample).pop()
 
     return central, original, new_partner
+
+def generate_closed_shell_singlet(m: Molecule):Add commentMore actions
+    for i in range(len(m.atoms)):
+            m.atoms[i].id = i
+    assert m.multiplicity == 1 and m.get_radical_count()>0
+    radical_center_ids = [x.id for x in m.atoms if x.radical_electrons > 0]
+    # remove radicals from radical centers (2)
+    for radical_center_id in radical_center_ids:
+        m.atoms[radical_center_id].decrement_radical()
+    # add removed radicals (2) to one of the radical sites as a lone pair (1)
+    m.atoms[radical_center_ids[0]].increment_lone_pairs()
+    # pick the best resonance structure
+    print([x.smiles for x in m.generate_resonance_structures()])
+    return m.generate_resonance_structures()[0]
+
+def generate_singlet_diradicals(m: Molecule):
+    for i in range(len(m.atoms)):
+        m.atoms[i].id = i
+    
+    singlet_diradicals = [] 
+    for edge in m.get_all_edges():
+        
+        M = m.copy(deep=True) 
+        if edge.get_order_num() > 1: # find a pi bond
+            atom1_id = edge.atom1.id
+            atom2_id = edge.atom2.id
+            M.atoms[atom1_id].increment_radical()  # add a radical to each atom of the pi bond
+            M.atoms[atom2_id].increment_radical()
+            M.get_bond(M.atoms[atom1_id], M.atoms[atom2_id]).decrement_order() # remove 1 pi bond
+            potential_singlet_diradicals = M.generate_resonance_structures()  # generate resonance structures
+            
+            for potential_singlet_diradical in potential_singlet_diradicals: # find all resonance structures with non-neighboring radical sites
+                radical_center_ids = sorted([x.id for x in potential_singlet_diradical.atoms if x.radical_electrons==1]) 
+                potential_singlet_diradical_edges = potential_singlet_diradical.get_all_edges()
+                potential_singlet_diradical_edge_ids = [sorted([x.atom1.id, x.atom2.id]) for x in potential_singlet_diradical_edges]
+                if radical_center_ids not in potential_singlet_diradical_edge_ids:
+                    if potential_singlet_diradical not in singlet_diradicals:
+                        singlet_diradicals.append(potential_singlet_diradical)Add commentMore actions
+    return singlet_diradicals
\ No newline at end of file
diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index d677755b07..1dac846b8a 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -864,7 +864,9 @@ def process_new_reactions(self, new_reactions, new_species, pdep_network=None, g
                 pdep = rxn.generate_high_p_limit_kinetics()
             elif any([any([x.is_subgraph_isomorphic(q) for q in self.unrealgroups]) for y in rxn.reactants + rxn.products for x in y.molecule]):
                 pdep = False
-
+            elif not rxn.reversible:
+                pdep = False
+                
             # If pressure dependence is on, we only add reactions that are not unimolecular;
             # unimolecular reactions will be added after processing the associated networks
             if not pdep:
diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py
index 834f0f7243..508a1892d6 100644
--- a/rmgpy/rmg/pdep.py
+++ b/rmgpy/rmg/pdep.py
@@ -816,14 +816,6 @@ def update(self, reaction_model, pdep_settings, requires_rms=False):
                     spec.conformer = Conformer(E0=spec.get_thermo_data().E0)
                 E0.append(spec.conformer.E0.value_si)
 
-        # Use the average E0 as the reference energy (`energy_correction`) for the network
-        # The `energy_correction` will be added to the free energies and enthalpies for each
-        # configuration in the network.
-        energy_correction = -np.array(E0).mean()
-        for spec in self.reactants + self.products + self.isomers:
-            spec.energy_correction = energy_correction
-        self.energy_correction = energy_correction
-
         # Determine transition state energies on potential energy surface
         # In the absence of any better information, we simply set it to
         # be the reactant ground-state energy + the activation energy
@@ -851,9 +843,9 @@ def update(self, reaction_model, pdep_settings, requires_rms=False):
                                 'type "{2!s}".'.format(rxn, self.index, rxn.kinetics.__class__))
             rxn.fix_barrier_height(force_positive=True)
             if rxn.network_kinetics is None:
-                E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.kinetics.Ea.value_si + energy_correction
+                E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.kinetics.Ea.value_si 
             else:
-                E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.network_kinetics.Ea.value_si + energy_correction
+                E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.network_kinetics.Ea.value_si 
             rxn.transition_state = rmgpy.species.TransitionState(conformer=Conformer(E0=(E0 * 0.001, "kJ/mol")))
 
         # Set collision model

From 15cdb8caf4f06c2a8ddc68a4670b1fb750e85088 Mon Sep 17 00:00:00 2001
From: donerancl <donera@mit.edu>
Date: Mon, 23 Jun 2025 11:23:11 -0400
Subject: [PATCH 3/6] pdep spin conservation updates

---
 rmgpy/data/kinetics/family.py | 32 ++++++++++++++++++++--------
 rmgpy/molecule/molecule.py    | 15 --------------
 rmgpy/molecule/util.py        | 39 -----------------------------------
 rmgpy/rmg/model.py            |  1 -
 4 files changed, 23 insertions(+), 64 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index d949eaf48d..5d63574db1 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -1704,7 +1704,7 @@ def is_molecule_forbidden(self, molecule):
 
         return False
 
-    def _create_reaction(self, reactants, products, is_forward):
+    def _create_reaction(self, reactants, products, is_forward, check_spin = True):
         """
         Create and return a new :class:`Reaction` object containing the
         provided `reactants` and `products` as lists of :class:`Molecule`
@@ -1739,7 +1739,11 @@ def _create_reaction(self, reactants, products, is_forward):
         for key, species_list in zip(['reactants', 'products'], [reaction.reactants, reaction.products]):
             for species in species_list:
                 reaction.labeled_atoms[key] = dict(reaction.labeled_atoms[key], **species.get_all_labeled_atoms())
-
+        if check_spin:
+            if not reaction.check_if_spin_allowed():
+                logging.info("Did not create the following reaction, which violates conservation of spin...")
+                logging.info(str(reaction))
+                return None
         return reaction
 
     def _match_reactant_to_template(self, reactant, template_reactant):
@@ -1802,6 +1806,7 @@ def generate_reactions(self, reactants, products=None, prod_resonance=True, dele
             specified reactants and products within this family.
             Degenerate reactions are returned as separate reactions.
         """
+        check_spin = True
         reaction_list = []
 
         # Forward direction (the direction in which kinetics is defined)
@@ -1989,7 +1994,7 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
                 specified reactants and products within this family.
             Degenerate reactions are returned as separate reactions.
         """
-
+        check_spin = True
         rxn_list = []
 
         # Wrap each reactant in a list if not already done (this is done to
@@ -2045,7 +2050,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
                             pass
                         else:
                             if product_structures is not None:
-                                rxn = self._create_reaction(reactant_structures, product_structures, forward)
+                                if self.label in allowed_spin_violation_families:
+                                    check_spin = False
+                                rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin)
                                 if rxn:
                                     rxn_list.append(rxn)
         # Bimolecular reactants: A + B --> products
@@ -2088,7 +2095,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
                                     pass
                                 else:
                                     if product_structures is not None:
-                                        rxn = self._create_reaction(reactant_structures, product_structures, forward)
+                                        if self.label in allowed_spin_violation_families:
+                                            check_spin = False
+                                        rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin)
                                         if rxn:
                                             rxn_list.append(rxn)
 
@@ -2112,8 +2121,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
                                         pass
                                     else:
                                         if product_structures is not None:
-                                            rxn = self._create_reaction(reactant_structures, product_structures,
-                                                                        forward)
+                                            if self.label in allowed_spin_violation_families:
+                                                check_spin = False
+                                            rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin)
                                             if rxn:
                                                 rxn_list.append(rxn)
 
@@ -2166,7 +2176,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
                             pass
                         else:
                             if product_structures is not None:
-                                rxn = self._create_reaction(reactant_structures, product_structures, forward)
+                                if self.label in allowed_spin_violation_families:
+                                    check_spin = False
+                                rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin)
                                 if rxn:
                                     rxn_list.append(rxn)
             else:
@@ -2231,7 +2243,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
                             pass
                         else:
                             if product_structures is not None:
-                                rxn = self._create_reaction(reactant_structures, product_structures, forward)
+                                if self.label in allowed_spin_violation_families:
+                                    check_spin = False
+                                rxn = self._create_reaction(reactant_structures, product_structures, forward, check_spin = check_spin)
                                 if rxn:
                                     rxn_list.append(rxn)
 
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 1cd0f3e515..1a32f07add 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -60,7 +60,6 @@
 from rmgpy.molecule.kekulize import kekulize
 from rmgpy.molecule.pathfinder import find_shortest_path
 from rmgpy.molecule.fragment import CuttingLabel
-from rmgpy.molecule.util import generate_closed_shell_singlet, generate_singlet_diradicals
 
 ################################################################################
 
@@ -3003,20 +3002,6 @@ def get_desorbed_molecules(self):
 
         return desorbed_molecules
 
-    def update_to_closed_shell_singlet(self):Add commentMore actions
-        for i in range(len(self.atoms)):
-                self.atoms[i].id = i
-        assert self.multiplicity == 1 and self.get_radical_count()>0
-        radical_center_ids = [x.id for x in self.atoms if x.radical_electrons > 0]
-        # remove radicals from radical centers (2)
-        for radical_center_id in radical_center_ids:
-            self.atoms[radical_center_id].decrement_radical()
-        # add removed radicals (2) to one of the radical sites as a lone pair (1)
-        self.atoms[radical_center_ids[0]].increment_lone_pairs()
-        # pick the best resonance structure
-        self.generate_resonance_structures()
-        self.reactive = True
-
 # this variable is used to name atom IDs so that there are as few conflicts by 
 # using the entire space of integer objects
 atom_id_counter = -2 ** 15
diff --git a/rmgpy/molecule/util.py b/rmgpy/molecule/util.py
index 0e3dc04171..46f8c3504c 100644
--- a/rmgpy/molecule/util.py
+++ b/rmgpy/molecule/util.py
@@ -179,42 +179,3 @@ def swap(to_be_swapped, sample):
     new_partner = (to_be_swapped - sample).pop()
 
     return central, original, new_partner
-
-def generate_closed_shell_singlet(m: Molecule):Add commentMore actions
-    for i in range(len(m.atoms)):
-            m.atoms[i].id = i
-    assert m.multiplicity == 1 and m.get_radical_count()>0
-    radical_center_ids = [x.id for x in m.atoms if x.radical_electrons > 0]
-    # remove radicals from radical centers (2)
-    for radical_center_id in radical_center_ids:
-        m.atoms[radical_center_id].decrement_radical()
-    # add removed radicals (2) to one of the radical sites as a lone pair (1)
-    m.atoms[radical_center_ids[0]].increment_lone_pairs()
-    # pick the best resonance structure
-    print([x.smiles for x in m.generate_resonance_structures()])
-    return m.generate_resonance_structures()[0]
-
-def generate_singlet_diradicals(m: Molecule):
-    for i in range(len(m.atoms)):
-        m.atoms[i].id = i
-    
-    singlet_diradicals = [] 
-    for edge in m.get_all_edges():
-        
-        M = m.copy(deep=True) 
-        if edge.get_order_num() > 1: # find a pi bond
-            atom1_id = edge.atom1.id
-            atom2_id = edge.atom2.id
-            M.atoms[atom1_id].increment_radical()  # add a radical to each atom of the pi bond
-            M.atoms[atom2_id].increment_radical()
-            M.get_bond(M.atoms[atom1_id], M.atoms[atom2_id]).decrement_order() # remove 1 pi bond
-            potential_singlet_diradicals = M.generate_resonance_structures()  # generate resonance structures
-            
-            for potential_singlet_diradical in potential_singlet_diradicals: # find all resonance structures with non-neighboring radical sites
-                radical_center_ids = sorted([x.id for x in potential_singlet_diradical.atoms if x.radical_electrons==1]) 
-                potential_singlet_diradical_edges = potential_singlet_diradical.get_all_edges()
-                potential_singlet_diradical_edge_ids = [sorted([x.atom1.id, x.atom2.id]) for x in potential_singlet_diradical_edges]
-                if radical_center_ids not in potential_singlet_diradical_edge_ids:
-                    if potential_singlet_diradical not in singlet_diradicals:
-                        singlet_diradicals.append(potential_singlet_diradical)Add commentMore actions
-    return singlet_diradicals
\ No newline at end of file
diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index 1dac846b8a..6cbcf77a4a 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -311,7 +311,6 @@ def make_new_species(self, object, label="", reactive=True, check_existing=True,
             try:
                 mols = molecule.cut_molecule(cut_through=False)
             except AttributeError:
-                # it's Molecule object, change it to Fragment and then cut
                 molecule = Fragment().from_adjacency_list(molecule.to_adjacency_list())
                 mols = molecule.cut_molecule(cut_through=False)
             if len(mols) == 1:

From edc03c4ab9b4a806705faeaab008f4061edc1c96 Mon Sep 17 00:00:00 2001
From: donerancl <donera@mit.edu>
Date: Mon, 23 Jun 2025 17:45:40 -0400
Subject: [PATCH 4/6] debugging pdep

---
 rmgpy/rmg/model.py | 141 ++++++++++++++++++++++++++++++++++++---------
 rmgpy/rmg/pdep.py  | 128 +++++++++++++++++++++++++---------------
 2 files changed, 195 insertions(+), 74 deletions(-)

diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index 6cbcf77a4a..856db801de 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -1945,36 +1945,123 @@ def update_unimolecular_reaction_networks(self):
         # Two partial networks having the same source and containing one or
         # more explored isomers in common must be merged together to avoid
         # double-counting of rates
-        for networks in self.network_dict.values():
-            network_count = len(networks)
-            for index0, network0 in enumerate(networks):
-                index = index0 + 1
-                while index < network_count:
-                    found = False
-                    network = networks[index]
-                    if network0.source == network.source:
-                        # The networks contain the same source, but do they contain any common included isomers (other than the source)?
-                        for isomer in network0.explored:
-                            if isomer != network.source and isomer in network.explored:
-                                # The networks contain an included isomer in common, so we need to merge them
-                                found = True
-                                break
-                    if found:
-                        # The networks contain the same source and one or more common included isomers
-                        # Therefore they need to be merged together
-                        logging.info("Merging PDepNetwork #{0:d} and PDepNetwork #{1:d}".format(network0.index, network.index))
-                        network0.merge(network)
-                        networks.remove(network)
-                        self.network_list.remove(network)
-                        network_count -= 1
-                    else:
-                        index += 1
+        # --- REVISED MERGING LOGIC FOR ROBUSTNESS ---
+
+        # 1. Collect all networks into a single, mutable list to process
+        working_networks = list(self.network_list) 
+
+        # Keep track of the indices of networks that have been merged INTO another network.
+        networks_to_exclude = set()
+
+        # Iterate through the working list.
+        for i in range(len(working_networks)):
+            network0 = working_networks[i]
+            
+            # If this network0 has already been merged into another, skip it.
+            if network0.index in networks_to_exclude:
+                continue
+
+            # Inner loop: 'network' is the potential network to be merged INTO network0
+            for j in range(i + 1, len(working_networks)):
+                network = working_networks[j]
+                
+                # If 'network' has already been merged, skip it.
+                if network.index in networks_to_exclude:
+                    continue
+
+                # Check for shared isomers (excluding source if desired by current logic)
+                shared_isomers = []
+                for isomer in network0.explored:
+                    # If you want to merge if *any* isomer is shared, remove `isomer != network.source`
+                    if isomer != network.source and isomer in network.explored:
+                        shared_isomers.append(isomer)
+                        break 
+
+                if shared_isomers: # 'found' is true if shared_isomers is not empty
+                    # --- LOGGING BEFORE MERGE ---
+                    logging.info("-" * 80)
+                    # Safely get source label(s) for logging
+                    source0_label = [x.label for x in network0.source] if isinstance(network0.source, list) else network0.source.label
+                    source_label = [x.label for x in network.source] if isinstance(network.source, list) else network.source.label
+                    logging.info(f"Attempting to merge PDepNetwork #{network0.index} (source: {source0_label}) and PDepNetwork #{network.index} (source: {source_label})")
+                    logging.info(f"Shared isomers triggering merge: {[s.label for s in shared_isomers]}")
+
+                    def log_species_energies(prefix, network_obj):
+                        """Helper function to log species energies, avoiding duplicates in output."""
+                        # Use a set to collect unique species for logging
+                        species_for_logging = set(network_obj.explored)
+                        if hasattr(network_obj, 'source'):
+                            # Ensure source handling is robust (single species vs. list of species)
+                            if isinstance(network_obj.source, list):
+                                species_for_logging.update(network_obj.source)
+                            else:
+                                species_for_logging.add(network_obj.source)
+                        
+                        # Add products. Assuming 'products' contains ProductSet objects
+                        # and ProductSet has a 'species' attribute (a list/tuple of Species objects).
+                        # If 'products' can contain single Species objects directly, adjust as needed.
+                        for prod_item in network_obj.products:
+                            if hasattr(prod_item, 'species') and prod_item.species is not None:
+                                species_for_logging.update(prod_item.species)
+                            else: # Assume it's a single Species object itself
+                                species_for_logging.add(prod_item)
+
+                        logging.info(f"{prefix} PDepNetwork #{network_obj.index} species energies:")
+                        
+                        # Sort for consistent output
+                        for sp in sorted(list(species_for_logging), key=lambda x: x.label if hasattr(x, 'label') else str(x)):
+                            try:
+                                if isinstance(sp, Species):
+                                    if hasattr(sp, 'conformer') and sp.conformer is not None and \
+                                       hasattr(sp.conformer, 'E0') and sp.conformer.E0 is not None:
+                                        logging.info(f"  {sp.label}: {sp.conformer.E0.value:.2f} {sp.conformer.E0.units}")
+                                    else:
+                                        logging.info(f"  {sp.label}: Conformer E0 not available")
+                                # This `else` block catches objects that are not `Species`
+                                # but might have an `E0` attribute (e.g., your ProductSet if it has one)
+                                else:
+                                    if hasattr(sp, 'label') and hasattr(sp, 'E0') and sp.E0 is not None:
+                                        logging.info(f"  {sp.label}: {sp.E0:.2f}")
+                                    elif hasattr(sp, 'E0') and sp.E0 is not None:
+                                        logging.info(f"  {str(sp)}: {sp.E0:.2f}") # Fallback for non-labeled E0 objects
+                                    else:
+                                        logging.info(f"  {sp.label if hasattr(sp, 'label') else str(sp)}: E0 not available")
+                            except Exception as e:
+                                logging.warning(f"  Error logging energy for {sp.label if hasattr(sp, 'label') else str(sp)}: {e}")
+
+                    log_species_energies("Before merge,", network0)
+                    log_species_energies("Before merge,", network)
+                    # --- END LOGGING BEFORE MERGE ---
+
+                    # Perform the merge. network is absorbed into network0.
+                    network0.merge(network) 
+                    
+                    # Mark 'network' for exclusion from the final list.
+                    networks_to_exclude.add(network.index) 
+
+                    # --- LOGGING AFTER MERGE ---
+                    logging.info(f"After merge, PDepNetwork #{network0.index} (merged network) species energies:")
+                    log_species_energies("After merge,", network0) # Log network0 *after* it has absorbed 'network'
+                    logging.info("-" * 80)
+                    # --- END LOGGING AFTER MERGE ---
+
+        # 2. Build the new self.network_list by excluding merged networks
+        self.network_list = [net for net in working_networks if net.index not in networks_to_exclude]
+
+        # 3. Rebuild self.network_dict from the updated self.network_list
+        self.network_dict.clear()
+        for network in self.network_list:
+            source_key = tuple(network.source) if isinstance(network.source, list) else network.source
+            if source_key not in self.network_dict:
+                self.network_dict[source_key] = []
+            self.network_dict[source_key].append(network)
+
+        # The rest of the method for updating invalid networks
+        # (This part of the count check needs adjustment if network.source can be a list)
+        count = sum([1 for network in self.network_list if not network.valid and not (len(network.explored) == 0 and (isinstance(network.source, list) and len(network.source) > 1 or not isinstance(network.source, list) and network.source is not None and network.source.label is not None))])
 
-        count = sum([1 for network in self.network_list if not network.valid and not (len(network.explored) == 0 and len(network.source) > 1)])
         logging.info("Updating {0:d} modified unimolecular reaction networks (out of {1:d})...".format(count, len(self.network_list)))
 
-        # Iterate over all the networks, updating the invalid ones as necessary
-        # self = reaction_model object
         updated_networks = []
         for network in self.network_list:
             if not network.valid:
diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py
index 508a1892d6..f473d74aa1 100644
--- a/rmgpy/rmg/pdep.py
+++ b/rmgpy/rmg/pdep.py
@@ -596,62 +596,96 @@ def merge(self, other):
         Merge the partial network `other` into this network.
         """
         # Make sure the two partial networks have the same source configuration
-        assert self.source == other.source
-
-        # Merge isomers
-        for isomer in other.isomers:
-            if isomer not in self.isomers:
-                self.isomers.append(isomer)
-        # Merge explored
-        for isomer in other.explored:
-            if isomer not in self.explored:
-                self.explored.append(isomer)
-        # Merge reactants
-        for reactants in other.reactants:
-            if reactants not in self.reactants:
-                self.reactants.append(reactants)
-        # Merge products
-        for products in other.products:
-            if products not in self.products:
-                self.products.append(products)
+        # assert self.source == other.source
+        # Use sets for merging to ensure uniqueness
+        self.isomers = list(set(self.isomers) | set(other.isomers))
+        self.explored = list(set(self.explored) | set(other.explored))
+        self.reactants = list(set(self.reactants) | set(other.reactants))
+        self.products = list(set(self.products) | set(other.products)) # <--- FIX FOR PRODUCTS
 
         # However, products that have been explored are actually isomers
-        # These should be removed from the list of products!
         products_to_remove = []
-        for products in self.products:
-            if len(products.species) == 1 and products.species[0] in self.isomers:
-                products_to_remove.append(products)
-        for products in products_to_remove:
-            self.products.remove(products)
+        for products_set in self.products: # Iterate through the ProductSet objects
+            # Check if all species in the products_set are in the merged isomers list
+            all_species_in_isomer_list = True
+            for sp in products_set.species: # Assuming ProductSet has a 'species' attribute (list of Species)
+                if sp not in self.isomers:
+                    all_species_in_isomer_list = False
+                    break
+            if all_species_in_isomer_list and len(products_set.species) == 1: # Only remove if it's a single isomer
+                products_to_remove.append(products_set)
 
-        # Merge path reactions
+        for products_set in products_to_remove:
+            self.products.remove(products_set)
+
+        # Merge path reactions (needs careful __eq__ for Reaction objects)
+        # Use a temporary set for efficient uniqueness check
+        current_path_reactions_set = set(self.path_reactions)
         for reaction in other.path_reactions:
-            found = False
-            for rxn in self.path_reactions:
-                if reaction.reactants == rxn.reactants and reaction.products == rxn.products:
-                    # NB the isEquivalent() method that used to be on the previous line also checked reverse direction.
-                    # I am not sure which is appropriate 
-                    found = True
-                    break
-            if not found:
-                self.path_reactions.append(reaction)
+            current_path_reactions_set.add(reaction) # Set automatically handles duplicates if Reaction has __eq__/__hash__
+        self.path_reactions = list(current_path_reactions_set)
 
-        # Also merge net reactions (so that when we update the network in the
-        # future, we update the existing net reactions rather than making new ones)
-        # Q: What to do when a net reaction exists in both networks being merged?
+        # Merge net reactions
+        current_net_reactions_set = set(self.net_reactions)
         for reaction in other.net_reactions:
-            found = False
-            for rxn in self.net_reactions:
-                if reaction.reactants == rxn.reactants and reaction.products == rxn.products:
-                    # NB the isEquivalent() method that used to be on the previous line also checked reverse direction.
-                    # I am not sure which is appropriate 
-                    found = True
-                    break
-            if not found:
-                self.net_reactions.append(reaction)
+            current_net_reactions_set.add(reaction)
+        self.net_reactions = list(current_net_reactions_set)
 
-        # Mark this network as invalid
         self.valid = False
+        # # Merge isomers
+        # for isomer in other.isomers:
+        #     if isomer not in self.isomers:
+        #         self.isomers.append(isomer)
+        # # Merge explored
+        # for isomer in other.explored:
+        #     if isomer not in self.explored:
+        #         self.explored.append(isomer)
+        # # Merge reactants
+        # for reactants in other.reactants:
+        #     if reactants not in self.reactants:
+        #         self.reactants.append(reactants)
+        # # Merge products
+        # for products in other.products:
+        #     if products not in self.products:
+        #         self.products.append(products)
+
+        # # However, products that have been explored are actually isomers
+        # # These should be removed from the list of products!
+        # products_to_remove = []
+        # for products in self.products:
+        #     if len(products.species) == 1 and products.species[0] in self.isomers:
+        #         products_to_remove.append(products)
+        # for products in products_to_remove:
+        #     self.products.remove(products)
+
+        # # Merge path reactions
+        # for reaction in other.path_reactions:
+        #     found = False
+        #     for rxn in self.path_reactions:
+        #         if reaction.reactants == rxn.reactants and reaction.products == rxn.products:
+        #             # NB the isEquivalent() method that used to be on the previous line also checked reverse direction.
+        #             # I am not sure which is appropriate 
+        #             found = True
+        #             break
+        #     if not found:
+        #         self.path_reactions.append(reaction)
+
+        # # Also merge net reactions (so that when we update the network in the
+        # # future, we update the existing net reactions rather than making new ones)
+        # # Q: What to do when a net reaction exists in both networks being merged?
+        # for reaction in other.net_reactions:
+        #     found = False
+        #     for rxn in self.net_reactions:
+        #         if reaction.reactants == rxn.reactants and reaction.products == rxn.products:
+        #             # NB the isEquivalent() method that used to be on the previous line also checked reverse direction.
+        #             # I am not sure which is appropriate 
+        #             found = True
+        #             break
+        #     if not found:
+        #         self.net_reactions.append(reaction)
+
+        # # Mark this network as invalid
+        # self.valid = False
 
     def update_configurations(self, reaction_model):
         """

From aac9a5651364c2886993c294a25ca6bcb4ef471c Mon Sep 17 00:00:00 2001
From: donerancl <donera@mit.edu>
Date: Wed, 25 Jun 2025 11:17:33 -0400
Subject: [PATCH 5/6]  revert update_unimolecular_reaction_network

---
 rmgpy/rmg/model.py | 142 +++++++++------------------------------------
 1 file changed, 28 insertions(+), 114 deletions(-)

diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py
index 856db801de..bd4bca92ff 100644
--- a/rmgpy/rmg/model.py
+++ b/rmgpy/rmg/model.py
@@ -1945,123 +1945,36 @@ def update_unimolecular_reaction_networks(self):
         # Two partial networks having the same source and containing one or
         # more explored isomers in common must be merged together to avoid
         # double-counting of rates
-        # --- REVISED MERGING LOGIC FOR ROBUSTNESS ---
-
-        # 1. Collect all networks into a single, mutable list to process
-        working_networks = list(self.network_list) 
-
-        # Keep track of the indices of networks that have been merged INTO another network.
-        networks_to_exclude = set()
-
-        # Iterate through the working list.
-        for i in range(len(working_networks)):
-            network0 = working_networks[i]
-            
-            # If this network0 has already been merged into another, skip it.
-            if network0.index in networks_to_exclude:
-                continue
-
-            # Inner loop: 'network' is the potential network to be merged INTO network0
-            for j in range(i + 1, len(working_networks)):
-                network = working_networks[j]
-                
-                # If 'network' has already been merged, skip it.
-                if network.index in networks_to_exclude:
-                    continue
-
-                # Check for shared isomers (excluding source if desired by current logic)
-                shared_isomers = []
-                for isomer in network0.explored:
-                    # If you want to merge if *any* isomer is shared, remove `isomer != network.source`
-                    if isomer != network.source and isomer in network.explored:
-                        shared_isomers.append(isomer)
-                        break 
-
-                if shared_isomers: # 'found' is true if shared_isomers is not empty
-                    # --- LOGGING BEFORE MERGE ---
-                    logging.info("-" * 80)
-                    # Safely get source label(s) for logging
-                    source0_label = [x.label for x in network0.source] if isinstance(network0.source, list) else network0.source.label
-                    source_label = [x.label for x in network.source] if isinstance(network.source, list) else network.source.label
-                    logging.info(f"Attempting to merge PDepNetwork #{network0.index} (source: {source0_label}) and PDepNetwork #{network.index} (source: {source_label})")
-                    logging.info(f"Shared isomers triggering merge: {[s.label for s in shared_isomers]}")
-
-                    def log_species_energies(prefix, network_obj):
-                        """Helper function to log species energies, avoiding duplicates in output."""
-                        # Use a set to collect unique species for logging
-                        species_for_logging = set(network_obj.explored)
-                        if hasattr(network_obj, 'source'):
-                            # Ensure source handling is robust (single species vs. list of species)
-                            if isinstance(network_obj.source, list):
-                                species_for_logging.update(network_obj.source)
-                            else:
-                                species_for_logging.add(network_obj.source)
-                        
-                        # Add products. Assuming 'products' contains ProductSet objects
-                        # and ProductSet has a 'species' attribute (a list/tuple of Species objects).
-                        # If 'products' can contain single Species objects directly, adjust as needed.
-                        for prod_item in network_obj.products:
-                            if hasattr(prod_item, 'species') and prod_item.species is not None:
-                                species_for_logging.update(prod_item.species)
-                            else: # Assume it's a single Species object itself
-                                species_for_logging.add(prod_item)
-
-                        logging.info(f"{prefix} PDepNetwork #{network_obj.index} species energies:")
-                        
-                        # Sort for consistent output
-                        for sp in sorted(list(species_for_logging), key=lambda x: x.label if hasattr(x, 'label') else str(x)):
-                            try:
-                                if isinstance(sp, Species):
-                                    if hasattr(sp, 'conformer') and sp.conformer is not None and \
-                                       hasattr(sp.conformer, 'E0') and sp.conformer.E0 is not None:
-                                        logging.info(f"  {sp.label}: {sp.conformer.E0.value:.2f} {sp.conformer.E0.units}")
-                                    else:
-                                        logging.info(f"  {sp.label}: Conformer E0 not available")
-                                # This `else` block catches objects that are not `Species`
-                                # but might have an `E0` attribute (e.g., your ProductSet if it has one)
-                                else:
-                                    if hasattr(sp, 'label') and hasattr(sp, 'E0') and sp.E0 is not None:
-                                        logging.info(f"  {sp.label}: {sp.E0:.2f}")
-                                    elif hasattr(sp, 'E0') and sp.E0 is not None:
-                                        logging.info(f"  {str(sp)}: {sp.E0:.2f}") # Fallback for non-labeled E0 objects
-                                    else:
-                                        logging.info(f"  {sp.label if hasattr(sp, 'label') else str(sp)}: E0 not available")
-                            except Exception as e:
-                                logging.warning(f"  Error logging energy for {sp.label if hasattr(sp, 'label') else str(sp)}: {e}")
-
-                    log_species_energies("Before merge,", network0)
-                    log_species_energies("Before merge,", network)
-                    # --- END LOGGING BEFORE MERGE ---
-
-                    # Perform the merge. network is absorbed into network0.
-                    network0.merge(network) 
-                    
-                    # Mark 'network' for exclusion from the final list.
-                    networks_to_exclude.add(network.index) 
-
-                    # --- LOGGING AFTER MERGE ---
-                    logging.info(f"After merge, PDepNetwork #{network0.index} (merged network) species energies:")
-                    log_species_energies("After merge,", network0) # Log network0 *after* it has absorbed 'network'
-                    logging.info("-" * 80)
-                    # --- END LOGGING AFTER MERGE ---
-
-        # 2. Build the new self.network_list by excluding merged networks
-        self.network_list = [net for net in working_networks if net.index not in networks_to_exclude]
-
-        # 3. Rebuild self.network_dict from the updated self.network_list
-        self.network_dict.clear()
-        for network in self.network_list:
-            source_key = tuple(network.source) if isinstance(network.source, list) else network.source
-            if source_key not in self.network_dict:
-                self.network_dict[source_key] = []
-            self.network_dict[source_key].append(network)
-
-        # The rest of the method for updating invalid networks
-        # (This part of the count check needs adjustment if network.source can be a list)
-        count = sum([1 for network in self.network_list if not network.valid and not (len(network.explored) == 0 and (isinstance(network.source, list) and len(network.source) > 1 or not isinstance(network.source, list) and network.source is not None and network.source.label is not None))])
+        for networks in self.network_dict.values():
+            network_count = len(networks)
+            for index0, network0 in enumerate(networks):
+                index = index0 + 1
+                while index < network_count:
+                    found = False
+                    network = networks[index]
+                    if network0.source == network.source:
+                        # The networks contain the same source, but do they contain any common included isomers (other than the source)?
+                        for isomer in network0.explored:
+                            if isomer != network.source and isomer in network.explored:
+                                # The networks contain an included isomer in common, so we need to merge them
+                                found = True
+                                break
+                    if found:
+                        # The networks contain the same source and one or more common included isomers
+                        # Therefore they need to be merged together
+                        logging.info("Merging PDepNetwork #{0:d} and PDepNetwork #{1:d}".format(network0.index, network.index))
+                        network0.merge(network)
+                        networks.remove(network)
+                        self.network_list.remove(network)
+                        network_count -= 1
+                    else:
+                        index += 1
 
+        count = sum([1 for network in self.network_list if not network.valid and not (len(network.explored) == 0 and len(network.source) > 1)])
         logging.info("Updating {0:d} modified unimolecular reaction networks (out of {1:d})...".format(count, len(self.network_list)))
 
+        # Iterate over all the networks, updating the invalid ones as necessary
+        # self = reaction_model object
         updated_networks = []
         for network in self.network_list:
             if not network.valid:
@@ -2118,6 +2031,7 @@ def log_species_energies(prefix, network_obj):
             # Move to the next core reaction
             index += 1
 
+
     def mark_chemkin_duplicates(self):
         """
         Check that all reactions that will appear the chemkin output have been checked as duplicates.

From 960dfb4bf0d396fdfb788556ee27c210747ca7cd Mon Sep 17 00:00:00 2001
From: donerancl <donera@mit.edu>
Date: Wed, 25 Jun 2025 11:22:41 -0400
Subject: [PATCH 6/6]  revrt def merge

---
 rmgpy/rmg/pdep.py | 128 +++++++++++++++++-----------------------------
 1 file changed, 47 insertions(+), 81 deletions(-)

diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py
index f473d74aa1..508a1892d6 100644
--- a/rmgpy/rmg/pdep.py
+++ b/rmgpy/rmg/pdep.py
@@ -596,96 +596,62 @@ def merge(self, other):
         Merge the partial network `other` into this network.
         """
         # Make sure the two partial networks have the same source configuration
-        # assert self.source == other.source
-        # Use sets for merging to ensure uniqueness
-        self.isomers = list(set(self.isomers) | set(other.isomers))
-        self.explored = list(set(self.explored) | set(other.explored))
-        self.reactants = list(set(self.reactants) | set(other.reactants))
-        self.products = list(set(self.products) | set(other.products)) # <--- FIX FOR PRODUCTS
+        assert self.source == other.source
+
+        # Merge isomers
+        for isomer in other.isomers:
+            if isomer not in self.isomers:
+                self.isomers.append(isomer)
+        # Merge explored
+        for isomer in other.explored:
+            if isomer not in self.explored:
+                self.explored.append(isomer)
+        # Merge reactants
+        for reactants in other.reactants:
+            if reactants not in self.reactants:
+                self.reactants.append(reactants)
+        # Merge products
+        for products in other.products:
+            if products not in self.products:
+                self.products.append(products)
 
         # However, products that have been explored are actually isomers
+        # These should be removed from the list of products!
         products_to_remove = []
-        for products_set in self.products: # Iterate through the ProductSet objects
-            # Check if all species in the products_set are in the merged isomers list
-            all_species_in_isomer_list = True
-            for sp in products_set.species: # Assuming ProductSet has a 'species' attribute (list of Species)
-                if sp not in self.isomers:
-                    all_species_in_isomer_list = False
-                    break
-            if all_species_in_isomer_list and len(products_set.species) == 1: # Only remove if it's a single isomer
-                products_to_remove.append(products_set)
-
-        for products_set in products_to_remove:
-            self.products.remove(products_set)
+        for products in self.products:
+            if len(products.species) == 1 and products.species[0] in self.isomers:
+                products_to_remove.append(products)
+        for products in products_to_remove:
+            self.products.remove(products)
 
-        # Merge path reactions (needs careful __eq__ for Reaction objects)
-        # Use a temporary set for efficient uniqueness check
-        current_path_reactions_set = set(self.path_reactions)
+        # Merge path reactions
         for reaction in other.path_reactions:
-            current_path_reactions_set.add(reaction) # Set automatically handles duplicates if Reaction has __eq__/__hash__
-        self.path_reactions = list(current_path_reactions_set)
+            found = False
+            for rxn in self.path_reactions:
+                if reaction.reactants == rxn.reactants and reaction.products == rxn.products:
+                    # NB the isEquivalent() method that used to be on the previous line also checked reverse direction.
+                    # I am not sure which is appropriate 
+                    found = True
+                    break
+            if not found:
+                self.path_reactions.append(reaction)
 
-        # Merge net reactions
-        current_net_reactions_set = set(self.net_reactions)
+        # Also merge net reactions (so that when we update the network in the
+        # future, we update the existing net reactions rather than making new ones)
+        # Q: What to do when a net reaction exists in both networks being merged?
         for reaction in other.net_reactions:
-            current_net_reactions_set.add(reaction)
-        self.net_reactions = list(current_net_reactions_set)
+            found = False
+            for rxn in self.net_reactions:
+                if reaction.reactants == rxn.reactants and reaction.products == rxn.products:
+                    # NB the isEquivalent() method that used to be on the previous line also checked reverse direction.
+                    # I am not sure which is appropriate 
+                    found = True
+                    break
+            if not found:
+                self.net_reactions.append(reaction)
 
+        # Mark this network as invalid
         self.valid = False
-        # # Merge isomers
-        # for isomer in other.isomers:
-        #     if isomer not in self.isomers:
-        #         self.isomers.append(isomer)
-        # # Merge explored
-        # for isomer in other.explored:
-        #     if isomer not in self.explored:
-        #         self.explored.append(isomer)
-        # # Merge reactants
-        # for reactants in other.reactants:
-        #     if reactants not in self.reactants:
-        #         self.reactants.append(reactants)
-        # # Merge products
-        # for products in other.products:
-        #     if products not in self.products:
-        #         self.products.append(products)
-
-        # # However, products that have been explored are actually isomers
-        # # These should be removed from the list of products!
-        # products_to_remove = []
-        # for products in self.products:
-        #     if len(products.species) == 1 and products.species[0] in self.isomers:
-        #         products_to_remove.append(products)
-        # for products in products_to_remove:
-        #     self.products.remove(products)
-
-        # # Merge path reactions
-        # for reaction in other.path_reactions:
-        #     found = False
-        #     for rxn in self.path_reactions:
-        #         if reaction.reactants == rxn.reactants and reaction.products == rxn.products:
-        #             # NB the isEquivalent() method that used to be on the previous line also checked reverse direction.
-        #             # I am not sure which is appropriate 
-        #             found = True
-        #             break
-        #     if not found:
-        #         self.path_reactions.append(reaction)
-
-        # # Also merge net reactions (so that when we update the network in the
-        # # future, we update the existing net reactions rather than making new ones)
-        # # Q: What to do when a net reaction exists in both networks being merged?
-        # for reaction in other.net_reactions:
-        #     found = False
-        #     for rxn in self.net_reactions:
-        #         if reaction.reactants == rxn.reactants and reaction.products == rxn.products:
-        #             # NB the isEquivalent() method that used to be on the previous line also checked reverse direction.
-        #             # I am not sure which is appropriate 
-        #             found = True
-        #             break
-        #     if not found:
-        #         self.net_reactions.append(reaction)
-
-        # # Mark this network as invalid
-        # self.valid = False
 
     def update_configurations(self, reaction_model):
         """