Merge pull request ReactionMechanismGenerator#2706 from ReactionMechanismGenerator/vdW_bonds

Vdw bonds

Author: JacksonBurns

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()
-            reason = fails_species_constraints(struct)
-            if reason:
+            if (reason := fails_species_constraints(struct)):
                 raise ForbiddenStructureException(
                     "Species constraints forbids product species {0}. Please "
                     "reformulate constraints, or explicitly "
@@ -1676,11 +1675,17 @@ def is_molecule_forbidden(self, molecule):
             return True
 
         # forbid vdw multi-dentate molecules for surface families
-        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" in self.label.lower() and molecule.is_multidentate():
+            if "vdwbidentate" in self.label.lower():
+                # Within vdWBidentate families, allow vdW in
+                # multi-dentate molecules if at least one bond to the surface
+                # is covalent.
+                if not molecule.has_covalent_surface_bond():
+                    return True
+            else:
+                # for all other families, forbid multi-dentate molecules with any vdW bonds
+                if molecule.has_vdw_surface_bond():
+                    return True
 
         return False
 

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

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:

rmgpy/molecule/molecule.pxd

@@ -206,6 +206,10 @@ cdef class Molecule(Graph):
 
     cpdef int number_of_surface_sites(self) except -1
 
+    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)

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,
@@ -1219,6 +1219,36 @@ 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(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 has_vdw_surface_bond(self):
+        """
+        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):
         """
         Returns ``True`` iff the molecule contains an 'X' surface site.
@@ -1273,9 +1303,14 @@ 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():
+            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)
@@ -3048,9 +3083,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.vertices 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):
@@ -3076,6 +3115,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,
@@ -3114,6 +3154,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 = '*0'
                 else:
                     raise NotImplementedError("Can't remove surface bond of type {}".format(bond.order))
             desorbed_molecule.remove_atom(site)

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)
 
@@ -61,4 +62,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)
+cpdef list find_adsorbate_conjugate_delocalization_paths(Vertex atom1)
+
+cpdef list find_formate_delocalization_paths(Vertex atom1)

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():
@@ -535,3 +553,40 @@ 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:
+
+    - [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.
+
+    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=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():
+            if atom2.is_oxygen() and bond12.is_van_der_waals():
+                for atom3, bond23 in atom2.edges.items():
+                    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():
+                                    if atom5.is_surface_site() and bond45.is_single():
+                                        paths.append([atom1, atom2, atom3, atom4, atom5, bond12, bond23, bond34, bond45])
+    return paths

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_adsorbate_formate_resonance_structures(Graph mol)

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_adsorbate_formate_resonance_structures,
         ]
     else:
         # If the molecule is aromatic, then radical resonance has already been considered
@@ -124,6 +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_adsorbate_formate_resonance_structures)
     return method_list
 
 
@@ -1257,3 +1259,43 @@ def generate_adsorbate_conjugate_resonance_structures(mol):
                                 else:
                                     structures.append(structure)
     return structures
+
+
+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
+    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)
+    """
+    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() and bond23.is_double() and bond34.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

test/database/databaseTest.py

@@ -1577,7 +1577,22 @@ def make_error_message(reactants, message=""):
                 output += "\n" + s.to_adjacency_list(label=s.to_smiles())
             return output
 
-        if len(sample_reactants) == 1 == len(family.forward_template.reactants):
+        expected_reactants = [str(r) for r in family.forward_template.reactants]
+        roots_with_samples = [str(k) for k in sample_reactants.keys()]
+        roots_without_samples = [r for r in expected_reactants if r not in roots_with_samples]
+        if len(sample_reactants) != len(family.forward_template.reactants):
+            # One or more reactant roots produced no usable sample molecule (every candidate was
+            # forbidden by is_molecule_forbidden, or raised UnexpectedChargeError/
+            # ImplicitBenzeneError during make_sample_molecule). Record a descriptive error so it
+            # is logged below alongside any per-sample errors, instead of raising opaquely here.
+            test1.append(
+                f"In family {family_name}, {len(roots_without_samples)} reactant root(s) produced "
+                f"no usable sample molecule: {roots_without_samples}. The family template expects "
+                f"{len(expected_reactants)} reactant(s) {expected_reactants}; only these root(s) "
+                f"yielded samples: {roots_with_samples}. Check the group definitions for the "
+                f"missing reactant root(s)."
+            )
+        elif len(sample_reactants) == 1:
             reactants = list(sample_reactants.values())[0]
             for reactant in reactants:
                 try:
@@ -1670,7 +1685,13 @@ def make_error_message(reactants, message=""):
                     species = rmgpy.species.Species(index=1, molecule=[molecule])
                     species.generate_resonance_structures()
         else:
-            raise ValueError(f"Family had {len(sample_reactants)} reactants?: " f"{', '.join(map(str,sample_reactants.keys())) }")
+            # Reactant count matches the template but is not 1, 2, or 3 (RMG only supports up to
+            # trimolecular). This is not expected for any well-formed family.
+            raise ValueError(
+                f"In family {family_name}, the number of sampled reactant roots "
+                f"({len(sample_reactants)}) matches the template reactant count but is not 1, 2, "
+                f"or 3, which is unexpected: {roots_with_samples}."
+            )
 
         # print out entries skipped from exception we can't currently handle
         if skipped: