From 34ba00ce3bfa8d16d8713338ffeb67c4daf2c694 Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Fri, 12 Jul 2024 16:13:08 -0400
Subject: [PATCH 01/23] draw vdW bond, the same way we do H bonds

combine drawing of vdW and H bonds

Two commits by  Bjarne Kreitz <bjarne_kreitz@brown.edu>
Merged by Richard West
---
 rmgpy/molecule/converter.py | 4 ++--
 rmgpy/molecule/draw.py      | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/rmgpy/molecule/converter.py b/rmgpy/molecule/converter.py
index ae005cb25b6..5142f916527 100644
--- a/rmgpy/molecule/converter.py
+++ b/rmgpy/molecule/converter.py
@@ -99,10 +99,10 @@ def to_rdkit_mol(mol, remove_h=True, return_mapping=False, sanitize=True,
             label_dict[index] = atom.label
 
     rd_bonds = Chem.rdchem.BondType
-    # no vdW bond in RDKit, so "ZERO" or "OTHER" might be OK
+    # no vdW bond in RDKit, so use UNSPECIFIED
     orders = {'S': rd_bonds.SINGLE, 'D': rd_bonds.DOUBLE,
               'T': rd_bonds.TRIPLE, 'B': rd_bonds.AROMATIC,
-              'Q': rd_bonds.QUADRUPLE, 'vdW': rd_bonds.ZERO,
+              'Q': rd_bonds.QUADRUPLE, 'vdW': rd_bonds.UNSPECIFIED,
               'H': rd_bonds.HYDROGEN, 'R': rd_bonds.UNSPECIFIED,
               None: rd_bonds.UNSPECIFIED}
     # Add the bonds
diff --git a/rmgpy/molecule/draw.py b/rmgpy/molecule/draw.py
index c6931e46ea6..377ca12a39c 100644
--- a/rmgpy/molecule/draw.py
+++ b/rmgpy/molecule/draw.py
@@ -1215,7 +1215,7 @@ def _render_bond(self, atom1, atom2, bond, cr):
                 dv *= 1.6
                 self._draw_line(cr, x1 - du, y1 - dv, x2 - du, y2 - dv)
                 self._draw_line(cr, x1 + du, y1 + dv, x2 + du, y2 + dv, dashed=True)
-            elif bond.is_hydrogen_bond():
+            elif bond.is_hydrogen_bond() or bond.is_van_der_waals():
                 # Draw a dashed line
                 self._draw_line(cr, x1, y1, x2, y2, dashed=True, dash_sizes=[0.5, 3.5])
             else:

From 96c5ea4e6d336b7d6d6b087a7c6b1b4f2c5a8ff6 Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Mon, 15 Jul 2024 11:09:46 -0400
Subject: [PATCH 02/23] Adjust remove vdW bond, using new
 has_covalent_surface_bond function

Combination of two commits by Bjarne Kreitz <bjarne_kreitz@brown.edu>
Merged by Richard West.
---
 rmgpy/molecule/molecule.py | 27 ++++++++++++++++++++++++---
 1 file changed, 24 insertions(+), 3 deletions(-)

diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 1d4f4752e42..f379d6ce606 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -1208,6 +1208,19 @@ def has_bond(self, atom1, atom2):
         """
         return self.has_edge(atom1, atom2)
 
+    def has_covalent_surface_bond(self):
+        """
+        Return True if any bond in this molecule connects a surface site (X) via a covalent bond.
+        """
+        cython.declare(bond=Bond)
+        for bond in self.get_all_edges():
+            if bond.is_van_der_waals():
+                continue
+            atom1, atom2 = bond.atom1, bond.atom2
+            if atom1.is_surface_site() or atom2.is_surface_site():
+                return True
+        return False
+
     def contains_surface_site(self):
         """
         Returns ``True`` iff the molecule contains an 'X' surface site.
@@ -1266,9 +1279,17 @@ def remove_van_der_waals_bonds(self):
         Remove all van der Waals bonds.
         """
         cython.declare(bond=Bond)
-        for bond in self.get_all_edges():
-            if bond.is_van_der_waals():
-                self.remove_bond(bond)
+        if self.is_multidentate():
+            if self.has_covalent_surface_bond():
+                return #preserve the remaining vdW bonds for this structure
+            else:
+                for bond in self.get_all_edges():
+                    if bond.is_van_der_waals():
+                        self.remove_bond(bond)
+        else:
+            for bond in self.get_all_edges():
+                if bond.is_van_der_waals():
+                    self.remove_bond(bond)
 
     def sort_atoms(self):
         """

From 219d9bbe0875bb742e7e7e03fe76cc302405d8f3 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 18 May 2026 23:20:43 -0400
Subject: [PATCH 03/23] Make has_covalent_surface_bond more efficient.

Rather than generate all the bonds, each with two ends (which happens
by iterating over the atoms) and then iterating over the bonds, and
seeing if either end is an X.
We instead iterate over the atoms (there are fewer) and only if one is
an X do we check its bonds. As soon as one is found that's not vdW,
we return True.  Also added a Cython declaration.
---
 rmgpy/molecule/molecule.pxd |  2 ++
 rmgpy/molecule/molecule.py  | 13 ++++++-------
 2 files changed, 8 insertions(+), 7 deletions(-)

diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd
index a7f13e6ce42..969e56e5b20 100644
--- a/rmgpy/molecule/molecule.pxd
+++ b/rmgpy/molecule/molecule.pxd
@@ -183,6 +183,8 @@ cdef class Molecule(Graph):
     cpdef bint is_proton(self)
 
     cpdef bint contains_surface_site(self)
+
+    cpdef bint has_covalent_surface_bond(self)
     
     cpdef bint is_surface_site(self)
 
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index f379d6ce606..82aa51dbf8d 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -1212,13 +1212,12 @@ def has_covalent_surface_bond(self):
         """
         Return True if any bond in this molecule connects a surface site (X) via a covalent bond.
         """
-        cython.declare(bond=Bond)
-        for bond in self.get_all_edges():
-            if bond.is_van_der_waals():
-                continue
-            atom1, atom2 = bond.atom1, bond.atom2
-            if atom1.is_surface_site() or atom2.is_surface_site():
-                return True
+        cython.declare(atom=Atom, bond=Bond)
+        for atom in self.atoms:
+            if atom.is_surface_site():
+                for bond in atom.bonds.values():
+                    if not bond.is_van_der_waals():
+                            return True
         return False
 
     def contains_surface_site(self):

From 7e7e6218a91c1aeac18ec2bb6f3930522515ea62 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 18 May 2026 23:27:01 -0400
Subject: [PATCH 04/23] Add test_has_covalent_surface_bond to test
 has_covalent_surface_bond.

---
 test/rmgpy/molecule/moleculeTest.py | 38 +++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/test/rmgpy/molecule/moleculeTest.py b/test/rmgpy/molecule/moleculeTest.py
index 2c3896148de..33970460366 100644
--- a/test/rmgpy/molecule/moleculeTest.py
+++ b/test/rmgpy/molecule/moleculeTest.py
@@ -3068,6 +3068,44 @@ def test_remove_van_der_waals_bonds(self):
         mol.remove_van_der_waals_bonds()
         assert len(mol.get_all_edges()) == 1
 
+    def test_has_covalent_surface_bond(self):
+        """Test Molecule.has_covalent_surface_bond() distinguishes vdW from covalent X bonds."""
+        # X present but only physisorbed via a vdW bond
+        vdw_only = Molecule().from_adjacency_list(
+            """
+1 X u0 p0 c0 {2,vdW}
+2 H u0 p0 c0 {1,vdW} {3,S}
+3 H u0 p0 c0 {2,S}
+"""
+        )
+        assert not vdw_only.has_covalent_surface_bond()
+
+        # X covalently bonded (chemisorbed)
+        chemisorbed = Molecule().from_adjacency_list(
+            """
+1 H u0 p0 c0 {2,S}
+2 X u0 p0 c0 {1,S}
+"""
+        )
+        assert chemisorbed.has_covalent_surface_bond()
+
+        # Two X atoms: one vdW, one covalent
+        mixed = Molecule().from_adjacency_list(
+            """
+1 X u0 p0 c0 {3,S}
+2 X u0 p0 c0 {4,vdW}
+3 C u0 p0 c0 {1,S} {4,S} {5,S} {6,S}
+4 H u0 p0 c0 {2,vdW} {3,S}
+5 H u0 p0 c0 {3,S}
+6 H u0 p0 c0 {3,S}
+"""
+        )
+        assert mixed.has_covalent_surface_bond()
+
+        # No surface sites at all
+        gas = Molecule().from_smiles("CCO")
+        assert not gas.has_covalent_surface_bond()
+
     def test_get_relevant_cycles(self):
         """
         Test the Molecule.get_relevant_cycles() raises correct error after deprecation.

From c5ea3f75892e201444640b22075d432def33f8b1 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 18 May 2026 23:32:08 -0400
Subject: [PATCH 05/23] Simplify remove_van_der_waals_bonds

We don't need to check if it's multidentate and THEN check
if there's a covalent surface bond.
Knowing that there's a covalent surface bond will suffice,
and still only involves looping through the atoms once each.
(same as counting the X atoms).

And if we haven't returned early, we should remove the bonds.
---
 rmgpy/molecule/molecule.py | 16 +++++-----------
 1 file changed, 5 insertions(+), 11 deletions(-)

diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 82aa51dbf8d..35eb9a05ceb 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -1278,17 +1278,11 @@ def remove_van_der_waals_bonds(self):
         Remove all van der Waals bonds.
         """
         cython.declare(bond=Bond)
-        if self.is_multidentate():
-            if self.has_covalent_surface_bond():
-                return #preserve the remaining vdW bonds for this structure
-            else:
-                for bond in self.get_all_edges():
-                    if bond.is_van_der_waals():
-                        self.remove_bond(bond)
-        else:
-            for bond in self.get_all_edges():
-                if bond.is_van_der_waals():
-                    self.remove_bond(bond)
+        if self.has_covalent_surface_bond():
+            return # preserve any vdW bonds if there's also a covalent X
+        for bond in self.get_all_edges():
+            if bond.is_van_der_waals():
+                self.remove_bond(bond)
 
     def sort_atoms(self):
         """

From 713131e3ac6f99060786eadba7f097f05fb738a9 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Mon, 18 May 2026 23:39:20 -0400
Subject: [PATCH 06/23] Accelerate is_multidentate.

The optimized is_multidentate now short-circuits on the second
surface site, and avoids building an intermediate list.
---
 rmgpy/molecule/molecule.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 35eb9a05ceb..479ca979be3 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -3018,9 +3018,13 @@ def is_multidentate(self):
         Return ``True`` if the adsorbate contains at least two binding sites,
         or ``False`` otherwise.
         """
-        cython.declare(atom=Atom)
-        if len([atom for atom in self.atoms if atom.is_surface_site()])>=2:
-            return True
+        cython.declare(atom=Atom, found_one=cython.bint)
+        found_one = False
+        for atom in self.atoms:
+            if atom.is_surface_site():
+                if found_one:
+                    return True
+                found_one = True
         return False
 
     def get_adatoms(self):

From ebc926145c661b023bcec670b55fe2b1c16da7fb Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Wed, 17 Jul 2024 09:03:32 -0400
Subject: [PATCH 07/23] change is_molecule_forbidden

---
 rmgpy/data/kinetics/family.py | 14 ++++++++++----
 rmgpy/molecule/molecule.py    |  4 +++-
 2 files changed, 13 insertions(+), 5 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 506b4e915f5..87fd92765a4 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -1676,11 +1676,17 @@ def is_molecule_forbidden(self, molecule):
             return True
 
         # forbid vdw multi-dentate molecules for surface families
+        surface_sites = []
         if "surface" in self.label.lower():
-            if molecule.get_num_atoms('X') > 1:
-                for atom in molecule.atoms:
-                    if atom.atomtype.label == 'Xv':
-                        return True
+            if "surface_monodentate_to_vdw_bidentate" in self.label.lower() and molecule.get_num_atoms('X') > 1:
+                surface_sites = [atom.atomtype.label for atom in molecule.atoms if 'X' in atom.atomtype.label]
+                if all(site == 'Xv' for site in surface_sites):
+                    return True
+            else:
+                if molecule.get_num_atoms('X') > 1:
+                    for atom in molecule.atoms:
+                        if atom.atomtype.label == 'Xv':
+                            return True
 
         return False
 
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 479ca979be3..a23546a28c5 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -67,7 +67,7 @@
 def _skip_first(in_tuple):
     return in_tuple[1:]
 
-bond_orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
+bond_orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5, 'vdW': 0}
 
 globals().update({
     'bond_orders': bond_orders,
@@ -3088,6 +3088,8 @@ def get_desorbed_molecules(self):
                     bonded_atom.increment_radical()
                     bonded_atom.increment_lone_pairs()
                     bonded_atom.label = '*4'
+                elif bond.is_van_der_waals():
+                    bonded_atom.label = '*5'
                 else:
                     raise NotImplementedError("Can't remove surface bond of type {}".format(bond.order))
             desorbed_molecule.remove_atom(site)

From 2cf12a2758ce393e0004b49adf8dbc6ecbd4bbcd Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Tue, 12 May 2026 21:05:54 -0400
Subject: [PATCH 08/23] adding find_formate_delocalization_paths to pathfinder

---
 rmgpy/molecule/pathfinder.pxd |  4 +++-
 rmgpy/molecule/pathfinder.py  | 28 ++++++++++++++++++++++++++++
 2 files changed, 31 insertions(+), 1 deletion(-)

diff --git a/rmgpy/molecule/pathfinder.pxd b/rmgpy/molecule/pathfinder.pxd
index 469417983fb..d34ae224e2f 100644
--- a/rmgpy/molecule/pathfinder.pxd
+++ b/rmgpy/molecule/pathfinder.pxd
@@ -61,4 +61,6 @@ cpdef bint is_atom_able_to_lose_lone_pair(Vertex atom)
 
 cpdef list find_adsorbate_delocalization_paths(Vertex atom1)
 
-cpdef list find_adsorbate_conjugate_delocalization_paths(Vertex atom1)
\ No newline at end of file
+cpdef list find_adsorbate_conjugate_delocalization_paths(Vertex atom1)
+
+cpdef list find_formate_delocalization_paths(Vertex atom1)
\ No newline at end of file
diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py
index 8bce591a4dd..9e88cebf52b 100644
--- a/rmgpy/molecule/pathfinder.py
+++ b/rmgpy/molecule/pathfinder.py
@@ -535,3 +535,31 @@ def find_adsorbate_conjugate_delocalization_paths(atom1):
                                     if atom5.is_surface_site():
                                         paths.append([atom1, atom2, atom3, atom4, atom5, bond12, bond23, bond34, bond45])
     return paths
+
+def find_formate_delocalization_paths(atom1):
+    """
+    Find all resonance structures which have a bonding configuration X...O-C-O-X.
+    Examples:
+
+    - XOC(H)XO/XOC(H)XO, where X is the surface site. The adsorption site X
+      is always placed on the left-hand side of the adatom and every adatom
+      is bonded to only one surface site X.
+
+    In this transition atom1 and atom5 are surface sites while atom2
+    and atom4 are oxygen and atom3 is a carbon atom.
+    """
+
+    cython.declare(paths=list, atom2=Vertex, atom3=Vertex, atom4=Vertex, atom5=Vertex, bond12=Edge, bond23=Edge, bond34=Edge, bond45=Edge)
+
+    paths = []
+    if atom1.is_surface_site():
+        for atom2, bond12 in atom1.edges.items():
+            if atom2.is_oxygen() and bond12.is_van_der_waals():
+                for atom3, bond23 in atom2.edges.items():
+                    if atom3.is_carbon():
+                        for atom4, bond34 in atom3.edges.items():
+                            if atom2 is not atom4 and atom4.is_oxygen():
+                                for atom5, bond45 in atom4.edges.items():
+                                    if atom5.is_surface_site():
+                                        paths.append([atom1, atom2, atom3, atom4, atom5, bond12, bond23, bond34, bond45])
+    return paths
\ No newline at end of file

From 7aa5f704d5caf67a4c566bfd5c0026aefe88913e Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 11:31:33 -0400
Subject: [PATCH 09/23] Fix some Cython declarations in pathfinder.py

If we declare these things as Atoms and Bonds then Cython can
directly access the methods like .is_surface_site() and .is_van_der_waals()

I noticed the pxd declares a whole ton of things as Vertex
instead of Atom. This should probably also be fixed, if I'm right.
---
 rmgpy/molecule/pathfinder.pxd |  3 ++-
 rmgpy/molecule/pathfinder.py  | 37 ++++++++++++++++++++++++++++++-----
 2 files changed, 34 insertions(+), 6 deletions(-)

diff --git a/rmgpy/molecule/pathfinder.pxd b/rmgpy/molecule/pathfinder.pxd
index d34ae224e2f..1403d046bcd 100644
--- a/rmgpy/molecule/pathfinder.pxd
+++ b/rmgpy/molecule/pathfinder.pxd
@@ -26,6 +26,7 @@
 ###############################################################################
 
 from .graph cimport Vertex, Edge, Graph
+from .molecule cimport Atom, Bond, Molecule
 
 cpdef list find_butadiene(Vertex start, Vertex end)
 
@@ -63,4 +64,4 @@ cpdef list find_adsorbate_delocalization_paths(Vertex atom1)
 
 cpdef list find_adsorbate_conjugate_delocalization_paths(Vertex atom1)
 
-cpdef list find_formate_delocalization_paths(Vertex atom1)
\ No newline at end of file
+cpdef list find_formate_delocalization_paths(Vertex atom1)
diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py
index 9e88cebf52b..33756be3ebd 100644
--- a/rmgpy/molecule/pathfinder.py
+++ b/rmgpy/molecule/pathfinder.py
@@ -36,7 +36,7 @@
 
 import cython
 
-from rmgpy.molecule.molecule import Atom
+from rmgpy.molecule.molecule import Atom, Bond
 from rmgpy.molecule.graph import Vertex, Edge
 
 def find_butadiene(start, end):
@@ -494,7 +494,15 @@ def find_adsorbate_delocalization_paths(atom1):
     In this transition atom1 and atom4 are surface sites while atom2 and atom3
     are carbon or nitrogen atoms.
     """
-    cython.declare(paths=list, atom2=Vertex, atom3=Vertex, atom4=Vertex, bond12=Edge, bond23=Edge, bond34=Edge)
+    cython.declare(
+        paths=list,
+        atom2=Atom,
+        atom3=Atom,
+        atom4=Atom,
+        bond12=Bond,
+        bond23=Bond,
+        bond34=Bond,
+    )
 
     paths = []
     if atom1.is_surface_site():
@@ -521,7 +529,17 @@ def find_adsorbate_conjugate_delocalization_paths(atom1):
     and atom4 are carbon or nitrogen atoms.
     """
 
-    cython.declare(paths=list, atom2=Vertex, atom3=Vertex, atom4=Vertex, atom5=Vertex, bond12=Edge, bond23=Edge, bond34=Edge, bond45=Edge)
+    cython.declare(
+        paths=list,
+        atom2=Atom,
+        atom3=Atom,
+        atom4=Atom,
+        atom5=Atom,
+        bond12=Bond,
+        bond23=Bond,
+        bond34=Bond,
+        bond45=Bond,
+    )
 
     paths = []
     if atom1.is_surface_site():
@@ -549,8 +567,17 @@ def find_formate_delocalization_paths(atom1):
     and atom4 are oxygen and atom3 is a carbon atom.
     """
 
-    cython.declare(paths=list, atom2=Vertex, atom3=Vertex, atom4=Vertex, atom5=Vertex, bond12=Edge, bond23=Edge, bond34=Edge, bond45=Edge)
-
+    cython.declare(
+        paths=list,
+        atom2=Atom,
+        atom3=Atom,
+        atom4=Atom,
+        atom5=Atom,
+        bond12=Bond,
+        bond23=Bond,
+        bond34=Bond,
+        bond45=Bond,
+    )
     paths = []
     if atom1.is_surface_site():
         for atom2, bond12 in atom1.edges.items():

From 485cffa60639470e3189315f433268d08784f2f1 Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Tue, 12 May 2026 21:13:11 -0400
Subject: [PATCH 10/23] add resonance structure generation for formate to
 resonance.py

---
 rmgpy/molecule/resonance.pxd |  2 ++
 rmgpy/molecule/resonance.py  | 40 ++++++++++++++++++++++++++++++++++++
 2 files changed, 42 insertions(+)

diff --git a/rmgpy/molecule/resonance.pxd b/rmgpy/molecule/resonance.pxd
index b95c0b8ca0b..e0b0381d646 100644
--- a/rmgpy/molecule/resonance.pxd
+++ b/rmgpy/molecule/resonance.pxd
@@ -73,3 +73,5 @@ cpdef list generate_adsorbate_shift_down_resonance_structures(Graph mol)
 cpdef list generate_adsorbate_shift_up_resonance_structures(Graph mol)
 
 cpdef list generate_adsorbate_conjugate_resonance_structures(Graph mol)
+
+cpdef list generate_formate_resonance_structures(Graph mol)
diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py
index 577d7088931..622e419921c 100644
--- a/rmgpy/molecule/resonance.py
+++ b/rmgpy/molecule/resonance.py
@@ -124,6 +124,7 @@ def populate_resonance_algorithms(features=None):
             method_list.append(generate_adsorbate_shift_down_resonance_structures)
             method_list.append(generate_adsorbate_shift_up_resonance_structures)
             method_list.append(generate_adsorbate_conjugate_resonance_structures)
+            method_list.append(generate_formate_resonance_structures)
     return method_list
 
 
@@ -1257,3 +1258,42 @@ def generate_adsorbate_conjugate_resonance_structures(mol):
                                 else:
                                     structures.append(structure)
     return structures
+
+
+def generate_formate_resonance_structures(mol):
+    """
+    Generate all of the resonance structures formed by the shift of two
+    electrons in a conjugated pi bond system of a bidentate adsorbate
+    with a bridging atom in between.
+
+    Example [X]OC(H)O[X]: [X]~OC(H)O[X] <=> [X]OC(H)O~[X]
+    (where '~' denotes a vdW bond)
+    """
+    cython.declare(structures=list, paths=list, index=cython.int, structure=Graph)
+    cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, atom3=Vertex, atom4=Vertex, atom5=Vertex, bond12=Edge, bond23=Edge, bond34=Edge, bond45=Edge)
+    cython.declare(v1=Vertex, v2=Vertex)
+
+    structures = []
+    if mol.is_multidentate():
+        for atom in mol.vertices:
+            paths = pathfinder.find_formate_delocalization_paths(atom)
+            for atom1, atom2, atom3, atom4, atom5, bond12, bond23, bond34, bond45 in paths:
+                if ((atom2.is_oxygen() and bond12.is_van_der_waals()) and
+                        (atom4.is_oxygen() and atom5.is_surface_site() and bond45.is_single())):
+                    bond12.increment_order()
+                    bond23.decrement_order()
+                    bond34.increment_order()
+                    bond45.decrement_order()
+                    structure = mol.copy(deep=True)
+                    bond12.decrement_order()
+                    bond23.increment_order()
+                    bond34.decrement_order()
+                    bond45.increment_order()
+                    try:
+                        structure.update_atomtypes(log_species=False)
+                    except AtomTypeError:
+                        pass
+                    else:
+                        structures.append(structure)
+
+    return structures

From b918ba8d0cb6386b2cf954f4b6b09c4147bd4e80 Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Tue, 12 May 2026 21:48:33 -0400
Subject: [PATCH 11/23] add label for vdW bond in get_desorbed_molecules

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

diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index a23546a28c5..12fc357161e 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -1275,7 +1275,10 @@ def remove_bond(self, bond):
 
     def remove_van_der_waals_bonds(self):
         """
-        Remove all van der Waals bonds.
+        Remove all van der Waals bonds. For multidentate species,
+        vdW bonds are preserved when there are still other
+        covalent bonds with the surface present. If no covalent surface bonds are present,
+        all vdW bonds are removed.
         """
         cython.declare(bond=Bond)
         if self.has_covalent_surface_bond():
@@ -3050,6 +3053,7 @@ def get_desorbed_molecules(self):
         ``*2`` - double bond
         ``*3`` - triple bond
         ``*4`` - quadruple bond
+        ``*0`` - vdW bond
         """
         cython.declare(desorbed_molecules=list, desorbed_molecule=Molecule, sites_to_remove=list, adsorbed_atoms=list,
                        site=Atom, numbonds=cython.int, bonded_atom=Atom, bond=Bond, i=cython.int, j=cython.int, atom0=Atom,
@@ -3089,7 +3093,7 @@ def get_desorbed_molecules(self):
                     bonded_atom.increment_lone_pairs()
                     bonded_atom.label = '*4'
                 elif bond.is_van_der_waals():
-                    bonded_atom.label = '*5'
+                    bonded_atom.label = '*0'
                 else:
                     raise NotImplementedError("Can't remove surface bond of type {}".format(bond.order))
             desorbed_molecule.remove_atom(site)

From a818d863df59cb263ef2b6b959cbdce2b2209ded Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Mon, 18 May 2026 17:39:08 -0400
Subject: [PATCH 12/23] add more constraints to formate path in resonance.py

---
 rmgpy/molecule/resonance.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py
index 622e419921c..756c5e0a855 100644
--- a/rmgpy/molecule/resonance.py
+++ b/rmgpy/molecule/resonance.py
@@ -1262,9 +1262,9 @@ def generate_adsorbate_conjugate_resonance_structures(mol):
 
 def generate_formate_resonance_structures(mol):
     """
-    Generate all of the resonance structures formed by the shift of two
-    electrons in a conjugated pi bond system of a bidentate adsorbate
-    with a bridging atom in between.
+    Generate all resonance structures formed by the shift of two
+    electrons in a conjugated bonding system of a bidentate adsorbate
+    with a bridging atom in between, where one bond to the surface is vdW.
 
     Example [X]OC(H)O[X]: [X]~OC(H)O[X] <=> [X]OC(H)O~[X]
     (where '~' denotes a vdW bond)
@@ -1279,7 +1279,8 @@ def generate_formate_resonance_structures(mol):
             paths = pathfinder.find_formate_delocalization_paths(atom)
             for atom1, atom2, atom3, atom4, atom5, bond12, bond23, bond34, bond45 in paths:
                 if ((atom2.is_oxygen() and bond12.is_van_der_waals()) and
-                        (atom4.is_oxygen() and atom5.is_surface_site() and bond45.is_single())):
+                        (atom4.is_oxygen() and atom5.is_surface_site() and
+                         bond45.is_single() and bond23.is_double() and bond34.is_single())):
                     bond12.increment_order()
                     bond23.decrement_order()
                     bond34.increment_order()

From 148ae76403903024ec590f36aaf9a261aa2cd2b3 Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Mon, 18 May 2026 17:40:55 -0400
Subject: [PATCH 13/23] clean up pathfinder and add constraints

---
 rmgpy/molecule/pathfinder.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py
index 33756be3ebd..11a65a5c2d2 100644
--- a/rmgpy/molecule/pathfinder.py
+++ b/rmgpy/molecule/pathfinder.py
@@ -556,10 +556,10 @@ def find_adsorbate_conjugate_delocalization_paths(atom1):
 
 def find_formate_delocalization_paths(atom1):
     """
-    Find all resonance structures which have a bonding configuration X...O-C-O-X.
+    Find all resonance structures which have a bonding configuration X~O=C-O-X.
     Examples:
 
-    - XOC(H)XO/XOC(H)XO, where X is the surface site. The adsorption site X
+    - [X]~OC(H)O[X]/[X]OC(H)O~[X], where '~' denotes a vdW bond and X is the surface site. The adsorption site X
       is always placed on the left-hand side of the adatom and every adatom
       is bonded to only one surface site X.
 
@@ -583,10 +583,10 @@ def find_formate_delocalization_paths(atom1):
         for atom2, bond12 in atom1.edges.items():
             if atom2.is_oxygen() and bond12.is_van_der_waals():
                 for atom3, bond23 in atom2.edges.items():
-                    if atom3.is_carbon():
+                    if atom3.is_carbon() and bond23.is_double():
                         for atom4, bond34 in atom3.edges.items():
-                            if atom2 is not atom4 and atom4.is_oxygen():
+                            if atom2 is not atom4 and atom4.is_oxygen() and bond34.is_single():
                                 for atom5, bond45 in atom4.edges.items():
-                                    if atom5.is_surface_site():
+                                    if atom5.is_surface_site() and bond45.is_single():
                                         paths.append([atom1, atom2, atom3, atom4, atom5, bond12, bond23, bond34, bond45])
     return paths
\ No newline at end of file

From a5bf55a4457eb29c094df65467df8d26cbb24dbb Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Mon, 18 May 2026 17:56:02 -0400
Subject: [PATCH 14/23] add features for N species

---
 rmgpy/molecule/pathfinder.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py
index 11a65a5c2d2..dc4d9eeeb59 100644
--- a/rmgpy/molecule/pathfinder.py
+++ b/rmgpy/molecule/pathfinder.py
@@ -583,7 +583,7 @@ def find_formate_delocalization_paths(atom1):
         for atom2, bond12 in atom1.edges.items():
             if atom2.is_oxygen() and bond12.is_van_der_waals():
                 for atom3, bond23 in atom2.edges.items():
-                    if atom3.is_carbon() and bond23.is_double():
+                    if (atom3.is_carbon() or atom3.is_nitrogen()) and bond23.is_double():
                         for atom4, bond34 in atom3.edges.items():
                             if atom2 is not atom4 and atom4.is_oxygen() and bond34.is_single():
                                 for atom5, bond45 in atom4.edges.items():

From a4faa1b4294628a41ac213c83b31041f3ec3fa29 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 08:55:30 -0400
Subject: [PATCH 15/23] Consider generate_formate_resonance_structures when no
 features specified.

Suggested by copilot:

Previously, populate_resonance_algorithms() adds generate_formate_resonance_structures only in the features['is_multidentate'] branch, but it is not included in the default (features=None) method list. This means call sites that iterate populate_resonance_algorithms() without passing features (e.g., the isomorphic-resonance enumeration later in this module) will never generate the new formate resonance structures. Consider adding the new generator to the features is None list as well to keep behavior consistent.
---
 rmgpy/molecule/resonance.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py
index 756c5e0a855..d9d17b33341 100644
--- a/rmgpy/molecule/resonance.py
+++ b/rmgpy/molecule/resonance.py
@@ -96,7 +96,8 @@ def populate_resonance_algorithms(features=None):
             generate_clar_structures,
             generate_adsorbate_shift_down_resonance_structures,
             generate_adsorbate_shift_up_resonance_structures,
-            generate_adsorbate_conjugate_resonance_structures
+            generate_adsorbate_conjugate_resonance_structures,
+            generate_formate_resonance_structures,
         ]
     else:
         # If the molecule is aromatic, then radical resonance has already been considered

From da8327355e9c69132b4f1237515b2033c192a82d Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 08:56:53 -0400
Subject: [PATCH 16/23] Rename generate_formate_resonance_structures to
 generate_adsorbate_formate_resonance_structures

Just so all the adsorbate things start with similar names.
---
 rmgpy/molecule/resonance.pxd | 2 +-
 rmgpy/molecule/resonance.py  | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/rmgpy/molecule/resonance.pxd b/rmgpy/molecule/resonance.pxd
index e0b0381d646..47aa3c01c18 100644
--- a/rmgpy/molecule/resonance.pxd
+++ b/rmgpy/molecule/resonance.pxd
@@ -74,4 +74,4 @@ cpdef list generate_adsorbate_shift_up_resonance_structures(Graph mol)
 
 cpdef list generate_adsorbate_conjugate_resonance_structures(Graph mol)
 
-cpdef list generate_formate_resonance_structures(Graph mol)
+cpdef list generate_adsorbate_formate_resonance_structures(Graph mol)
diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py
index d9d17b33341..23fc598684f 100644
--- a/rmgpy/molecule/resonance.py
+++ b/rmgpy/molecule/resonance.py
@@ -97,7 +97,7 @@ def populate_resonance_algorithms(features=None):
             generate_adsorbate_shift_down_resonance_structures,
             generate_adsorbate_shift_up_resonance_structures,
             generate_adsorbate_conjugate_resonance_structures,
-            generate_formate_resonance_structures,
+            generate_adsorbate_formate_resonance_structures,
         ]
     else:
         # If the molecule is aromatic, then radical resonance has already been considered
@@ -125,7 +125,7 @@ def populate_resonance_algorithms(features=None):
             method_list.append(generate_adsorbate_shift_down_resonance_structures)
             method_list.append(generate_adsorbate_shift_up_resonance_structures)
             method_list.append(generate_adsorbate_conjugate_resonance_structures)
-            method_list.append(generate_formate_resonance_structures)
+            method_list.append(generate_adsorbate_formate_resonance_structures)
     return method_list
 
 
@@ -1261,7 +1261,7 @@ def generate_adsorbate_conjugate_resonance_structures(mol):
     return structures
 
 
-def generate_formate_resonance_structures(mol):
+def generate_adsorbate_formate_resonance_structures(mol):
     """
     Generate all resonance structures formed by the shift of two
     electrons in a conjugated bonding system of a bidentate adsorbate

From 1d93792c6bfa9f13e489c80edc99e4cbf2cb384c Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 09:11:55 -0400
Subject: [PATCH 17/23] Test the new remove_van_der_waals_bonds behaviour.

Should NOT remove vdW if there's also a covalent bond to the surface.
---
 test/rmgpy/molecule/moleculeTest.py | 45 ++++++++++++++++++++++++++++-
 1 file changed, 44 insertions(+), 1 deletion(-)

diff --git a/test/rmgpy/molecule/moleculeTest.py b/test/rmgpy/molecule/moleculeTest.py
index 33970460366..5a439016f99 100644
--- a/test/rmgpy/molecule/moleculeTest.py
+++ b/test/rmgpy/molecule/moleculeTest.py
@@ -3056,7 +3056,7 @@ def test_count_aromatic_rings(self):
 
         assert result == [2, 1, 0]
 
-    def test_remove_van_der_waals_bonds(self):
+    def test_remove_van_der_waals_bond(self):
         """Test we can remove a van-der-Waals bond"""
         adjlist = """
 1 X  u0 p0 c0 {2,vdW}
@@ -3068,6 +3068,49 @@ def test_remove_van_der_waals_bonds(self):
         mol.remove_van_der_waals_bonds()
         assert len(mol.get_all_edges()) == 1
 
+    def test_remove_van_der_waals_bonds_bidentate(self):
+        """vdW bonds are preserved on bidentates that also have a covalent X bond, removed otherwise."""
+        # Bidentate ethyl-like fragment: one C–X covalent, one C~X vdW. The vdW should be preserved.
+        mixed = Molecule().from_adjacency_list(
+            """
+1 X u0 p0 c0 {3,S}
+2 X u0 p0 c0 {4,vdW}
+3 C u0 p0 c0 {1,S} {4,S} {5,S} {6,S}
+4 C u0 p0 c0 {2,vdW} {3,S} {7,S} {8,S} {9,S}
+5 H u0 p0 c0 {3,S}
+6 H u0 p0 c0 {3,S}
+7 H u0 p0 c0 {4,S}
+8 H u0 p0 c0 {4,S}
+9 H u0 p0 c0 {4,S}
+"""
+        )
+        assert mixed.has_covalent_surface_bond()
+        n_edges_before = len(mixed.get_all_edges())
+        mixed.remove_van_der_waals_bonds()
+        assert len(mixed.get_all_edges()) == n_edges_before  # nothing removed
+        assert any(bond.is_van_der_waals() for bond in mixed.get_all_edges())
+
+        # Bidentate physisorbed: both C~X bonds are vdW. All vdW bonds should be removed.
+        vdw_only = Molecule().from_adjacency_list(
+            """
+1 X u0 p0 c0 {3,vdW}
+2 X u0 p0 c0 {4,vdW}
+3 C u0 p0 c0 {1,vdW} {4,S} {5,S} {6,S} {7,S}
+4 C u0 p0 c0 {2,vdW} {3,S} {8,S} {9,S} {10,S}
+5 H u0 p0 c0 {3,S}
+6 H u0 p0 c0 {3,S}
+7 H u0 p0 c0 {3,S}
+8 H u0 p0 c0 {4,S}
+9 H u0 p0 c0 {4,S}
+10 H u0 p0 c0 {4,S}
+"""
+        )
+        assert not vdw_only.has_covalent_surface_bond()
+        n_edges_before = len(vdw_only.get_all_edges())
+        vdw_only.remove_van_der_waals_bonds()
+        assert len(vdw_only.get_all_edges()) == n_edges_before - 2
+        assert not any(bond.is_van_der_waals() for bond in vdw_only.get_all_edges())
+
     def test_has_covalent_surface_bond(self):
         """Test Molecule.has_covalent_surface_bond() distinguishes vdW from covalent X bonds."""
         # X present but only physisorbed via a vdW bond

From b91e790e99c92c3d8c07c6ee13a4114bf0ecf13c Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 09:54:05 -0400
Subject: [PATCH 18/23] Add comments describing current implementation of
 is_molecule_forbidden

For checking multi-dentate vdW things.
I'm about to refactor, and want to first document what the current code does,
then change how we do it in the next commit.
---
 rmgpy/data/kinetics/family.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 87fd92765a4..70a0401d588 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -1678,11 +1678,15 @@ def is_molecule_forbidden(self, molecule):
         # forbid vdw multi-dentate molecules for surface families
         surface_sites = []
         if "surface" in self.label.lower():
+            # Within the surface_monodentate_to_vdw_bidentate, allow (don't forbid)
+            # vdW in multi-dentate molecules if at least one bond to the surface
+            # is covalent (not vdW).
             if "surface_monodentate_to_vdw_bidentate" in self.label.lower() and molecule.get_num_atoms('X') > 1:
                 surface_sites = [atom.atomtype.label for atom in molecule.atoms if 'X' in atom.atomtype.label]
                 if all(site == 'Xv' for site in surface_sites):
                     return True
             else:
+                # for all other families, forbid multi-dentate molecules with any vdW bonds
                 if molecule.get_num_atoms('X') > 1:
                     for atom in molecule.atoms:
                         if atom.atomtype.label == 'Xv':

From 75a48965ed19e1ee2c14d08fadce4f62e7f96c73 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 10:33:01 -0400
Subject: [PATCH 19/23] Refactor is_molecule_forbidden checks in reaction
 family.

Lifts is_multidentate() into the outer guard (short-circuits on the second surface atom,
so it's cheaper than get_num_atoms('X') > 1 which counts all atoms).
Replaces the surface_sites list + all(... == 'Xv') check with
`not has_covalent_surface_bond()`.
Drops the now-unused surface_sites = [] initializer.
---
 rmgpy/data/kinetics/family.py | 21 +++++++++------------
 1 file changed, 9 insertions(+), 12 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 70a0401d588..57ceb3de89d 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -1676,21 +1676,18 @@ def is_molecule_forbidden(self, molecule):
             return True
 
         # forbid vdw multi-dentate molecules for surface families
-        surface_sites = []
-        if "surface" in self.label.lower():
-            # Within the surface_monodentate_to_vdw_bidentate, allow (don't forbid)
-            # vdW in multi-dentate molecules if at least one bond to the surface
-            # is covalent (not vdW).
-            if "surface_monodentate_to_vdw_bidentate" in self.label.lower() and molecule.get_num_atoms('X') > 1:
-                surface_sites = [atom.atomtype.label for atom in molecule.atoms if 'X' in atom.atomtype.label]
-                if all(site == 'Xv' for site in surface_sites):
+        if "surface" in self.label.lower() and molecule.is_multidentate():
+            if "surface_monodentate_to_vdw_bidentate" in self.label.lower():
+                # Within the surface_monodentate_to_vdw_bidentate family, allow
+                # (don't forbid) vdW in multi-dentate molecules if at least one
+                # bond to the surface is covalent (not vdW).
+                if not molecule.has_covalent_surface_bond():
                     return True
             else:
                 # for all other families, forbid multi-dentate molecules with any vdW bonds
-                if molecule.get_num_atoms('X') > 1:
-                    for atom in molecule.atoms:
-                        if atom.atomtype.label == 'Xv':
-                            return True
+                for atom in molecule.atoms:
+                    if atom.atomtype.label == 'Xv':
+                        return True
 
         return False
 

From 3865d753de910eb8912555139ca8a9dd771dfe42 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 10:44:47 -0400
Subject: [PATCH 20/23] Add has_vdw_surface_bond method on Molecule.

Also add tests for it, and use it in reaction family is_molecule_forbidden method.
---
 rmgpy/data/kinetics/family.py       |  5 ++--
 rmgpy/molecule/molecule.pxd         |  5 +++-
 rmgpy/molecule/molecule.py          | 12 +++++++++
 test/rmgpy/molecule/moleculeTest.py | 38 +++++++++++++++++++++++++++++
 4 files changed, 56 insertions(+), 4 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 57ceb3de89d..8c291cc9c98 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -1685,9 +1685,8 @@ def is_molecule_forbidden(self, molecule):
                     return True
             else:
                 # for all other families, forbid multi-dentate molecules with any vdW bonds
-                for atom in molecule.atoms:
-                    if atom.atomtype.label == 'Xv':
-                        return True
+                if molecule.has_vdw_surface_bond():
+                    return True
 
         return False
 
diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd
index 969e56e5b20..3316e51b100 100644
--- a/rmgpy/molecule/molecule.pxd
+++ b/rmgpy/molecule/molecule.pxd
@@ -185,7 +185,10 @@ cdef class Molecule(Graph):
     cpdef bint contains_surface_site(self)
 
     cpdef bint has_covalent_surface_bond(self)
-    
+
+    cpdef bint has_vdw_surface_bond(self)
+
+
     cpdef bint is_surface_site(self)
 
     cpdef remove_atom(self, Atom atom)
diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index 12fc357161e..d762d177cd7 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -1220,6 +1220,18 @@ def has_covalent_surface_bond(self):
                             return True
         return False
 
+    def has_vdw_surface_bond(self):
+        """
+        Return True if any bond in this molecule connects a surface site (X) via a van der Waals bond.
+        """
+        cython.declare(atom=Atom, bond=Bond)
+        for atom in self.atoms:
+            if atom.is_surface_site():
+                for bond in atom.bonds.values():
+                    if bond.is_van_der_waals():
+                        return True
+        return False
+
     def contains_surface_site(self):
         """
         Returns ``True`` iff the molecule contains an 'X' surface site.
diff --git a/test/rmgpy/molecule/moleculeTest.py b/test/rmgpy/molecule/moleculeTest.py
index 5a439016f99..78fd97fe69d 100644
--- a/test/rmgpy/molecule/moleculeTest.py
+++ b/test/rmgpy/molecule/moleculeTest.py
@@ -3149,6 +3149,44 @@ def test_has_covalent_surface_bond(self):
         gas = Molecule().from_smiles("CCO")
         assert not gas.has_covalent_surface_bond()
 
+    def test_has_vdw_surface_bond(self):
+        """Test Molecule.has_vdw_surface_bond() detects any vdW X bond."""
+        # X present but only physisorbed via a vdW bond
+        vdw_only = Molecule().from_adjacency_list(
+            """
+1 X u0 p0 c0 {2,vdW}
+2 H u0 p0 c0 {1,vdW} {3,S}
+3 H u0 p0 c0 {2,S}
+"""
+        )
+        assert vdw_only.has_vdw_surface_bond()
+
+        # X covalently bonded (chemisorbed) — no vdW bond
+        chemisorbed = Molecule().from_adjacency_list(
+            """
+1 H u0 p0 c0 {2,S}
+2 X u0 p0 c0 {1,S}
+"""
+        )
+        assert not chemisorbed.has_vdw_surface_bond()
+
+        # Two X atoms: one vdW, one covalent
+        mixed = Molecule().from_adjacency_list(
+            """
+1 X u0 p0 c0 {3,S}
+2 X u0 p0 c0 {4,vdW}
+3 C u0 p0 c0 {1,S} {4,S} {5,S} {6,S}
+4 H u0 p0 c0 {2,vdW} {3,S}
+5 H u0 p0 c0 {3,S}
+6 H u0 p0 c0 {3,S}
+"""
+        )
+        assert mixed.has_vdw_surface_bond()
+
+        # No surface sites at all
+        gas = Molecule().from_smiles("CCO")
+        assert not gas.has_vdw_surface_bond()
+
     def test_get_relevant_cycles(self):
         """
         Test the Molecule.get_relevant_cycles() raises correct error after deprecation.

From 3dc9bef3be571e08b2f13be138b95f925d21e460 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 22:16:11 -0400
Subject: [PATCH 21/23] Fix bug and add tests: unbonded X gives True for
 has_vdw_surface_bond

The has_vdw_surface_bond now ALSO returns True for has LACK OF
a bond (but an X and one other thing).

I guess if you were to give it an adjacency list with just X atoms,
this would return True, which may seem weird.
(But so would the old test which was just for the presence of an Xv atom type)
---
 rmgpy/molecule/molecule.py          |  8 +++++++-
 test/rmgpy/molecule/moleculeTest.py | 14 ++++++++++++++
 2 files changed, 21 insertions(+), 1 deletion(-)

diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py
index d762d177cd7..276f4ef3c54 100644
--- a/rmgpy/molecule/molecule.py
+++ b/rmgpy/molecule/molecule.py
@@ -1222,14 +1222,20 @@ def has_covalent_surface_bond(self):
 
     def has_vdw_surface_bond(self):
         """
-        Return True if any bond in this molecule connects a surface site (X) via a van der Waals bond.
+        Return True if any bond in this molecule connects a surface site (X)
+        via a van der Waals bond, or there's a surface site with no bonds
+        (but at least one other atom in the molecule).
         """
         cython.declare(atom=Atom, bond=Bond)
         for atom in self.atoms:
             if atom.is_surface_site():
+                if not atom.bonds:  # if there are no bonds at all
+                    if len(self.atoms) > 1:  # and there's something besides the surface site
+                        return True  # then treat as vdW bonded
                 for bond in atom.bonds.values():
                     if bond.is_van_der_waals():
                         return True
+                
         return False
 
     def contains_surface_site(self):
diff --git a/test/rmgpy/molecule/moleculeTest.py b/test/rmgpy/molecule/moleculeTest.py
index 78fd97fe69d..d0c31862e70 100644
--- a/test/rmgpy/molecule/moleculeTest.py
+++ b/test/rmgpy/molecule/moleculeTest.py
@@ -3187,6 +3187,20 @@ def test_has_vdw_surface_bond(self):
         gas = Molecule().from_smiles("CCO")
         assert not gas.has_vdw_surface_bond()
 
+        # An unbonded X atom counts as a vdW surface bond
+        unbonded = Molecule().from_adjacency_list(
+            """
+1 X u0 p0 c0
+2 H u0 p0 c0 {3,S}
+3 H u0 p0 c0 {2,S}
+"""
+        )
+        assert unbonded.has_vdw_surface_bond()
+
+        # vacant site alone is not a vdW surface bond
+        vacant = Molecule().from_adjacency_list("1 X u0 p0 c0")
+        assert not vacant.has_vdw_surface_bond()
+
     def test_get_relevant_cycles(self):
         """
         Test the Molecule.get_relevant_cycles() raises correct error after deprecation.

From d13f019ed43db969df06264369d6f4cc531b61b7 Mon Sep 17 00:00:00 2001
From: Richard West <r.west@northeastern.edu>
Date: Tue, 19 May 2026 22:17:13 -0400
Subject: [PATCH 22/23] Avoid calling fails_species_constraints() twice in a
 row.

Since Python 3.8 we can do this "walrus" operator :=
---
 rmgpy/data/kinetics/family.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 8c291cc9c98..9e8fea3e575 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -1655,8 +1655,7 @@ def _generate_product_structures(self, reactant_structures, maps, forward, relab
         for struct in product_structures:
             if self.is_molecule_forbidden(struct):
                 raise ForbiddenStructureException()
-            if fails_species_constraints(struct):
-                reason = fails_species_constraints(struct)
+            if (reason := fails_species_constraints(struct)):
                 raise ForbiddenStructureException(
                     "Species constraints forbids product species {0}. Please "
                     "reformulate constraints, or explicitly "

From ddbe82ac997b37b1ce4ad5061d68c38c2bd98ada Mon Sep 17 00:00:00 2001
From: Bjarne Kreitz <bjarne_kreitz@brown.edu>
Date: Wed, 20 May 2026 21:50:13 -0400
Subject: [PATCH 23/23] restructure is_molecule_forbbiden

---
 rmgpy/data/kinetics/family.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py
index 9e8fea3e575..33d7cd603e3 100644
--- a/rmgpy/data/kinetics/family.py
+++ b/rmgpy/data/kinetics/family.py
@@ -1676,7 +1676,13 @@ def is_molecule_forbidden(self, molecule):
 
         # forbid vdw multi-dentate molecules for surface families
         if "surface" in self.label.lower() and molecule.is_multidentate():
-            if "surface_monodentate_to_vdw_bidentate" in self.label.lower():
+            allowed_vdw_families = [
+                "surface_monodentate_to_vdw_bidentate",
+                "surface_dissociation_vdw_bidentate",
+                "surface_dissociation_vdw_bidentate_beta",
+                "surface_abstraction_vdw_bidentate_beta"
+            ]
+            if any(name in self.label.lower() for name in allowed_vdw_families):
                 # Within the surface_monodentate_to_vdw_bidentate family, allow
                 # (don't forbid) vdW in multi-dentate molecules if at least one
                 # bond to the surface is covalent (not vdW).