diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 506b4e915f..33d7cd603e 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 " @@ -1676,11 +1675,23 @@ 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(): + 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). + 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 diff --git a/rmgpy/molecule/converter.py b/rmgpy/molecule/converter.py index ae005cb25b..5142f91652 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 c6931e46ea..377ca12a39 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: diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index a7f13e6ce4..3316e51b10 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -183,7 +183,12 @@ cdef class Molecule(Graph): cpdef bint is_proton(self) 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 1d4f4752e4..276f4ef3c5 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, @@ -1208,6 +1208,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. @@ -1263,9 +1293,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) @@ -3004,9 +3039,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): @@ -3032,6 +3071,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, @@ -3070,6 +3110,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) diff --git a/rmgpy/molecule/pathfinder.pxd b/rmgpy/molecule/pathfinder.pxd index 469417983f..1403d046bc 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) @@ -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) \ No newline at end of file +cpdef list find_adsorbate_conjugate_delocalization_paths(Vertex atom1) + +cpdef list find_formate_delocalization_paths(Vertex atom1) diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py index 8bce591a4d..dc4d9eeeb5 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(): @@ -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 \ No newline at end of file diff --git a/rmgpy/molecule/resonance.pxd b/rmgpy/molecule/resonance.pxd index b95c0b8ca0..47aa3c01c1 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_adsorbate_formate_resonance_structures(Graph mol) diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py index 577d708893..23fc598684 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_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 diff --git a/test/rmgpy/molecule/moleculeTest.py b/test/rmgpy/molecule/moleculeTest.py index 2c3896148d..d0c31862e7 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,139 @@ 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 + 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_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() + + # 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.