Add handling for dissociative desorption reactions to flux pair generation

rmgpy/reaction.py

@@ -1500,14 +1500,15 @@ def generate_pairs(self):
         performing flux analysis. The exact procedure for doing so depends on
         the reaction type:
 
-        =================== =============== ========================================
-        Reaction type       Template        Resulting pairs
-        =================== =============== ========================================
-        Isomerization       A     -> C      (A,C)
-        Dissociation        A     -> C + D  (A,C), (A,D)
-        Association         A + B -> C      (A,C), (B,C)
-        Bimolecular         A + B -> C + D  (A,C), (B,D) *or* (A,D), (B,C)
-        =================== =============== ========================================
+        ======================= ==================== ========================================
+        Reaction type           Template             Resulting pairs
+        ======================= ==================== ========================================
+        Isomerization             A      -> C          (A,C)
+        Dissociation              A      -> C + D      (A,C), (A,D)
+        Association               A + B  -> C          (A,C), (B,C)
+        Bimolecular               A + B  -> C + D      (A,C), (B,D) *or* (A,D), (B,C)
+        Dissociative Adsorption   A + 2X -> CX + DX    (A,CX), (A,DX), (X,CX), (X,DX)
+        ======================= ==================== ========================================
 
         There are a number of ways of determining the correct pairing for
         bimolecular reactions. Here we try a simple similarity analysis by comparing
@@ -1521,6 +1522,29 @@ def generate_pairs(self):
             for reactant in self.reactants:
                 for product in self.products:
                     self.pairs.append((reactant, product))
+        elif (len(self.reactants) == 3 and len(self.products) == 2 and \
+              len([r for r in self.reactants if r.is_surface_site()]) == 2 and \
+              len([r for r in self.reactants if not r.contains_surface_site()]) == 1) or \
+             (len(self.products) == 3 and len(self.reactants) == 2 and \
+              len([p for p in self.products if p.is_surface_site()]) == 2 and \
+              len([p for p in self.products if not p.contains_surface_site()]) == 1):
+            # Dissociative adsorption case
+            if len(self.reactants) == 3:
+                gas_reactant = [r for r in self.reactants if not r.is_surface_site()][0]
+                for product in self.products:   
+                    self.pairs.append((gas_reactant, product))
+                site1 = [r for r in self.reactants if r.is_surface_site()][0]
+                site2 = [r for r in self.reactants if r.is_surface_site()][1]
+                self.pairs.append((site1, self.products[0]))
+                self.pairs.append((site2, self.products[1]))
+            else:
+                gas_product = [p for p in self.products if not p.is_surface_site()][0]
+                for reactant in self.reactants:   
+                    self.pairs.append((reactant, gas_product))
+                site1 = [p for p in self.products if p.is_surface_site()][0]
+                site2 = [p for p in self.products if p.is_surface_site()][1]
+                self.pairs.append((self.reactants[0], site1))
+                self.pairs.append((self.reactants[1], site2))
 
         else:  # this is the bimolecular case
             reactants = self.reactants[:]

test/rmgpy/reactionTest.py

@@ -785,6 +785,68 @@ def setup_class(self):
         )
         self.reaction4_pairs = [(PO3, HOPO2), (H2, H_atom)]
 
+        # Reaction definitions for testing flux pair generation
+        # Isomerization
+        PC4H9 = Species().from_smiles("[CH2]CCC")
+        SC4H9 = Species().from_smiles("C[CH]CCC")
+        self.isomerization_reaction = Reaction(
+            reactants=[PC4H9],
+            products=[SC4H9],
+        )
+        self.isomerization_reaction_pairs = [(PC4H9, SC4H9)]
+
+        # Association
+        propyl = Species().from_smiles("[CH2]CC")
+        methylene = Species().from_smiles("[CH2]")
+        self.association_reaction = Reaction(
+            reactants=[propyl, methylene],
+            products=[PC4H9],
+        )
+        self.association_reaction_pairs = [(propyl, PC4H9), (methylene, PC4H9)]
+
+        # Dissociation
+        self.dissociation_reaction = Reaction(
+            reactants=[PC4H9],
+            products=[propyl, methylene],
+        )
+        self.dissociation_reaction_pairs = [(PC4H9, propyl), (PC4H9, methylene)]
+
+        # Bimolecular
+        butane = Species().from_smiles("CCCC")
+        oh = Species().from_smiles("[OH]")
+        h2o = Species().from_smiles("O")
+        self.bimolecular_reaction = Reaction(
+            reactants=[butane, oh],
+            products=[SC4H9, h2o],
+        )
+        self.bimolecular_reaction_pairs = [(butane, SC4H9), (oh, h2o)]
+
+        # Surface dissociative adsorption
+        NH3 = Species().from_smiles("N")
+        X = Species().from_adjacency_list("1 X u0 p0")
+        HX = Species().from_adjacency_list(
+"""
+1 H u0 p0 c0 {2,S}
+2 X u0 p0 c0 {1,S}
+""")
+        H2NX = Species().from_adjacency_list(
+"""
+1 N u0 p1 c0 {2,S} {3,S} {4,S}
+2 H u0 p0 c0 {1,S}
+3 H u0 p0 c0 {1,S}
+4 X u0 p0 c0 {1,S}
+""")
+        self.surf_dissociation_reaction = Reaction(
+            reactants=[NH3, X, X],
+            products=[H2NX, HX],
+        )
+        self.surf_dissociation_reaction_pairs = [
+            (NH3, H2NX),
+            (NH3, HX),
+            (X, H2NX),
+            (X, HX),
+        ]
+
     def test_is_isomerization(self):
         """
         Test the Reaction.is_isomerization() method.
@@ -1824,6 +1886,57 @@ def test_phosphorus_reaction_pairs(self):
         assert len(self.reaction4.pairs[1]) == 2
         assert self.reaction4.pairs == self.reaction4_pairs
 
+    def test_generate_pairs(self):
+        """
+        This method tests that the correct reaction pairs are being generated for a reaction
+        """
+
+        def pairs_equal(pair1, pair2):
+            if len(pair1) != len(pair2):
+                return False
+            if len(pair1) != 2:
+                return False
+            return pair1[0].is_isomorphic(pair2[0]) and pair1[1].is_isomorphic(pair2[1]) or \
+                     pair1[0].is_isomorphic(pair2[1]) and pair1[1].is_isomorphic(pair2[0])
+        def list_of_pairs_equal(list1, list2):
+            if len(list1) != len(list2):
+                return False
+            for pair1 in list1:
+                found = False
+                for pair2 in list2:
+                    if pairs_equal(pair1, pair2):
+                        found = True
+                        list2.remove(pair2)
+                        break
+                if not found:
+                    return False
+            return True
+
+        # Test Isomerization: A -> C
+        self.isomerization_reaction.generate_pairs()
+        assert len(self.isomerization_reaction.pairs) == 1
+        assert list_of_pairs_equal(self.isomerization_reaction.pairs, self.isomerization_reaction_pairs)
+
+        # Test Dissociation A -> C + D
+        self.dissociation_reaction.generate_pairs()
+        assert len(self.dissociation_reaction.pairs) == 2
+        assert list_of_pairs_equal(self.dissociation_reaction.pairs, self.dissociation_reaction_pairs)
+
+        # Test Association A + B -> C
+        self.association_reaction.generate_pairs()
+        assert len(self.association_reaction.pairs) == 2
+        assert list_of_pairs_equal(self.association_reaction.pairs, self.association_reaction_pairs)
+
+        # Test Bimolecular A + B -> C + D
+        self.bimolecular_reaction.generate_pairs()
+        assert len(self.bimolecular_reaction.pairs) == 2
+        assert list_of_pairs_equal(self.bimolecular_reaction.pairs, self.bimolecular_reaction_pairs)
+
+        # Test Dissociative Adsorption A + 2X -> CX + DX
+        self.surf_dissociation_reaction.generate_pairs()
+        assert len(self.surf_dissociation_reaction.pairs) == 4
+        assert list_of_pairs_equal(self.surf_dissociation_reaction.pairs, self.surf_dissociation_reaction_pairs)
+
 
 class TestReactionToCantera:
     """