From a97d2d3e528d2d9ec9819c7af680f5670d001d5c Mon Sep 17 00:00:00 2001 From: bjkreitz Date: Mon, 31 Jul 2023 11:39:41 -0400 Subject: [PATCH] adding resonance for adsorbates --- rmgpy/molecule/molecule.pxd | 2 + rmgpy/molecule/molecule.py | 10 ++ rmgpy/molecule/pathfinder.pxd | 4 + rmgpy/molecule/pathfinder.py | 48 ++++++++ rmgpy/molecule/resonance.pxd | 10 ++ rmgpy/molecule/resonance.py | 199 +++++++++++++++++++++++++++++++++- 6 files changed, 272 insertions(+), 1 deletion(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 3d5062f8bb3..1e8f9e1939a 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -289,6 +289,8 @@ cdef class Molecule(Graph): cpdef list get_adatoms(self) + cpdef bint is_bidentate(self) + cpdef list get_desorbed_molecules(self) cdef atom_id_counter diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index dfaeb70c89c..daba584e1e1 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -2759,6 +2759,16 @@ def get_surface_sites(self): cython.declare(atom=Atom) return [atom for atom in self.atoms if atom.is_surface_site()] + def is_bidentate(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 + return False + def get_adatoms(self): """ Get a list of adatoms in the molecule. diff --git a/rmgpy/molecule/pathfinder.pxd b/rmgpy/molecule/pathfinder.pxd index 28238823ed5..5091e442d01 100644 --- a/rmgpy/molecule/pathfinder.pxd +++ b/rmgpy/molecule/pathfinder.pxd @@ -58,3 +58,7 @@ cpdef list find_N5dc_radical_delocalization_paths(Vertex atom1) cpdef bint is_atom_able_to_gain_lone_pair(Vertex atom) 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 diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py index a115b00d84f..d88e9db2057 100644 --- a/rmgpy/molecule/pathfinder.py +++ b/rmgpy/molecule/pathfinder.py @@ -480,3 +480,51 @@ def is_atom_able_to_lose_lone_pair(atom): return (((atom.is_nitrogen() or atom.is_sulfur()) and atom.lone_pairs in [1, 2, 3]) or (atom.is_oxygen() and atom.lone_pairs in [2, 3]) or atom.is_carbon() and atom.lone_pairs == 1) + +def find_adsorbate_delocalization_paths(atom1): + """ + Find all bidentate adsorbates which have a bonding configuration X-C-C-X. + Examples: + + - XCXC, XCHXCH, XCXCH, where X is the surface site + + In this transition atom1 and atom4 are surface sites while atom2 and atom3 are carbon atoms. + """ + cython.declare(paths=list, atom2=Vertex, atom3=Vertex, atom4=Vertex, bond12=Edge, bond23=Edge, bond34=Edge) + + path=[] + if atom1.is_surface_site(): + for atom2, bond12 in atom1.edges.items(): + if atom2.is_carbon(): + for atom3, bond23 in atom2.edges.items(): + if atom1 is not atom3 and atom3.is_carbon(): + for atom4, bond34 in atom3.edges.items(): + if atom2 is not atom4 and atom4.is_surface_site(): + path.append([atom1, atom2, atom3, atom4, bond12, bond23, bond34]) + return path + return path + +def find_adsorbate_conjugate_delocalization_paths(atom1): + """ + Find all bidentate adsorbates which have a bonding configuration X-C-C-C-X. + Examples: + + - XCHCHXCH/XCHCHXC, where X is the surface site + + In this transition atom1 and atom5 are surface sites while atom2, atom3, and atom4 are carbon atoms. + """ + + cython.declare(paths=list, atom2=Vertex, atom3=Vertex, atom4=Vertex, atom5=Vertex, bond12=Edge, bond23=Edge, bond34=Edge, bond45=Edge) + + path=[] + if atom1.is_surface_site(): + for atom2, bond12 in atom1.edges.items(): + for atom3, bond23 in atom2.edges.items(): + if atom1 is not atom3 and atom3.is_carbon(): + for atom4, bond34 in atom3.edges.items(): + if atom2 is not atom4 and atom4.is_carbon(): + for atom5, bond45 in atom4.edges.items(): + if atom3 is not atom5 and atom5.is_surface_site(): + path.append([atom1, atom2, atom3, atom4, atom5, bond12, bond23, bond34, bond45]) + return path + return path \ No newline at end of file diff --git a/rmgpy/molecule/resonance.pxd b/rmgpy/molecule/resonance.pxd index a9110880662..9651f0a695c 100644 --- a/rmgpy/molecule/resonance.pxd +++ b/rmgpy/molecule/resonance.pxd @@ -63,3 +63,13 @@ cpdef list generate_clar_structures(Graph mol, bint save_order=?) cpdef list _clar_optimization(Graph mol, list constraints=?, max_num=?, save_order=?) cpdef list _clar_transformation(Graph mol, list aromatic_ring) + +cpdef list generate_adsorbate_shift_down_resonance_structures(Graph mol) + +cpdef list generate_adsorbate_shift_up_resonance_structures(Graph mol) + +cpdef list generate_adsorbate_double_shift_up_resonance_structures(Graph mol) + +cpdef list generate_adsorbate_double_shift_down_resonance_structures(Graph mol) + +cpdef list generate_adsorbate_conjugate_resonance_structures(Graph mol) \ No newline at end of file diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py index dd9c2726ee2..78308efd3db 100644 --- a/rmgpy/molecule/resonance.py +++ b/rmgpy/molecule/resonance.py @@ -49,6 +49,12 @@ - ``generate_kekule_structure``: generate a single Kekule structure for an aromatic compound (single/double bond form) - ``generate_opposite_kekule_structure``: for monocyclic aromatic species, rotate the double bond assignment - ``generate_clar_structures``: generate all structures with the maximum number of pi-sextet assignments +- Multidentate adsorbates only + - ``generate_adsorbate_shift_down_resonance_structures``: shift 2 electrons from a C=/#C bond to the X-C bond + - ``generate_adsorbate_shift_up_resonance_structures``: shift 2 electrons from a X=/#C bond to a C-C bond + - ``generate_adsorbate_double_shift_down_resonance_structures``: shift 4 elecrons from a C#C bond to a X-C bond + - ``generate_adsorbate_double_shift_up_resonance_structures``: shift 4 electrons from a X#C bond to a C-C bond + - ``generate_adsorbate_conjugate_resonance_structures``: shift 2 electrons in a conjugate pi system for bridged X-C-C-C-X adsorbates """ import logging @@ -86,6 +92,11 @@ def populate_resonance_algorithms(features=None): generate_aryne_resonance_structures, generate_kekule_structure, generate_clar_structures, + generate_adsorbate_shift_down_resonance_structures, + generate_adsorbate_shift_up_resonance_structures, + generate_adsorbate_double_shift_down_resonance_structures, + generate_adsorbate_double_shift_up_resonance_structures, + generate_adsorbate_conjugate_resonance_structures ] else: # If the molecule is aromatic, then radical resonance has already been considered @@ -109,7 +120,12 @@ def populate_resonance_algorithms(features=None): # solution. A more holistic approach would be to identify these cases in generate_resonance_structures, # and pass a list of forbidden atom ID's to find_lone_pair_multiple_bond_paths. method_list.append(generate_lone_pair_multiple_bond_resonance_structures) - + if features['is_bidentate']: + method_list.append(generate_adsorbate_shift_down_resonance_structures) + method_list.append(generate_adsorbate_shift_up_resonance_structures) + method_list.append(generate_adsorbate_double_shift_down_resonance_structures) + method_list.append(generate_adsorbate_double_shift_up_resonance_structures) + method_list.append(generate_adsorbate_conjugate_resonance_structures) return method_list @@ -129,6 +145,7 @@ def analyze_molecule(mol, save_order=False): 'is_aryl_radical': False, 'hasNitrogenVal5': False, 'hasLonePairs': False, + 'is_bidentate':mol.is_bidentate(), } if features['is_cyclic']: @@ -1116,3 +1133,183 @@ def _clar_transformation(mol, aromatic_ring): for bond in bond_list: bond.order = 1.5 + + +def generate_adsorbate_shift_down_resonance_structures(mol): + """ + Generate all of the resonance structures formed by the shift a pi bond between two C-C atoms to both X-C bonds. + Example XCXCH: [X]C=C=[X] <=> [X]=CC#[X] + (where '=' denotes a double bond, '#' denotes a triple bond) + """ + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, atom3=Vertex, atom4=Vertex, bond12=Edge, bond23=Edge, bond34=Edge) + cython.declare(v1=Vertex, v2=Vertex) + + structures = [] + if mol.is_bidentate(): + for atom in mol.vertices: + # print(atom) + paths = pathfinder.find_adsorbate_delocalization_paths(atom) + for atom1, atom2, atom3, atom4, bond12, bond23, bond34 in paths: + if bond23.is_single(): + pass + else: + bond12.increment_order() + bond23.decrement_order() + bond34.increment_order() + structure = mol.copy(deep=True) + bond12.decrement_order() + bond23.increment_order() + bond34.decrement_order() + try: + structure.update_atomtypes(log_species=False) + except AtomTypeError: + pass + else: + structures.append(structure) + return structures + +def generate_adsorbate_shift_up_resonance_structures(mol): + """ + Generate all of the resonance structures formed by the shift of two electrons from X-C bonds to increase the bond + order between two C-C atoms by 1. + Example XCXCH: [X]=CC#[X] <=> [X]C=C=[X] + (where '=' denotes a double bond, '#' denotes a triple bond) + """ + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, atom3=Vertex, atom4=Vertex, bond12=Edge, bond23=Edge, bond34=Edge) + cython.declare(v1=Vertex, v2=Vertex) + + structures = [] + if mol.is_bidentate(): + for atom in mol.vertices: + paths = pathfinder.find_adsorbate_delocalization_paths(atom) + for atom1, atom2, atom3, atom4, bond12, bond23, bond34 in paths: + if (bond12.is_double() and bond23.is_single() and bond34.is_double()) or (bond23.is_double() and bond12.is_double() and bond34.is_double()) or (bond12.is_triple() and bond23.is_single() and bond34.is_triple()) or (bond12.is_triple() and bond23.is_single() and bond34.is_double()): + bond12.decrement_order() + bond23.increment_order() + bond34.decrement_order() + structure = mol.copy(deep=True) + bond12.increment_order() + bond23.decrement_order() + bond34.increment_order() + try: + structure.update_atomtypes(log_species=False) + except AtomTypeError: + pass + else: + structures.append(structure) + return structures + +def generate_adsorbate_double_shift_up_resonance_structures(mol): + """ + Generate all of the resonance structures formed by the shift of four electrons from X-C bonds to increase the bond + order between two C-C atoms by 2. + Example XCXC: [X]#CC#[X] <=> [X]C#C[X] + (where '#' denotes a triple bond) + """ + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, atom3=Vertex, atom4=Vertex, bond12=Edge, bond23=Edge, bond34=Edge) + cython.declare(v1=Vertex, v2=Vertex) + + structures = [] + if mol.is_bidentate(): + for atom in mol.vertices: + paths = pathfinder.find_adsorbate_delocalization_paths(atom) + for atom1, atom2, atom3, atom4, bond12, bond23, bond34 in paths: + if bond23.is_single() and bond12.is_triple() and bond34.is_triple(): + bond12.decrement_order() + bond12.decrement_order() + bond23.increment_order() + bond23.increment_order() + bond34.decrement_order() + bond34.decrement_order() + structure = mol.copy(deep=True) + bond12.increment_order() + bond12.increment_order() + bond23.decrement_order() + bond23.decrement_order() + bond34.increment_order() + bond34.increment_order() + try: + structure.update_atomtypes(log_species=False) + except AtomTypeError: + pass + else: + structures.append(structure) + return structures + +def generate_adsorbate_double_shift_down_resonance_structures(mol): + """ + Generate all of the resonance structures formed by the shift of four electrons from a C#C bond to increase the bond + order between two X-C atoms by 2. + Example XCXC: [X]C#C[X] <=> [X]#CC#[X] + (where '#' denotes a triple bond) + """ + cython.declare(structures=list, paths=list, index=cython.int, structure=Graph) + cython.declare(atom=Vertex, atom1=Vertex, atom2=Vertex, atom3=Vertex, atom4=Vertex, bond12=Edge, bond23=Edge, bond34=Edge) + cython.declare(v1=Vertex, v2=Vertex) + + structures = [] + if mol.is_bidentate(): + for atom in mol.vertices: + paths = pathfinder.find_adsorbate_delocalization_paths(atom) + for atom1, atom2, atom3, atom4, bond12, bond23, bond34 in paths: + if bond23.is_triple() and bond12.is_single() and bond34.is_single(): + bond12.increment_order() + bond12.increment_order() + bond23.decrement_order() + bond23.decrement_order() + bond34.increment_order() + bond34.increment_order() + structure = mol.copy(deep=True) + bond12.decrement_order() + bond12.decrement_order() + bond23.increment_order() + bond23.increment_order() + bond34.decrement_order() + bond34.decrement_order() + try: + structure.update_atomtypes(log_species=False) + except AtomTypeError: + pass + else: + structures.append(structure) + return structures + +def generate_adsorbate_conjugate_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 C atom in between. + Example XCHCHXC: [X]#CC=C[X] <=> [X]=C=CC=[X] + (where '#' denotes a triple bond, '=' denotes a double 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_bidentate(): + for atom in mol.vertices: + paths = pathfinder.find_adsorbate_conjugate_delocalization_paths(atom) + for atom1, atom2, atom3, atom4, atom45, bond12, bond23, bond34, bond45 in paths: + if bond23.is_single(): + if (bond12.is_double() or bond12.is_triple()): + if (bond34.is_double() or bond34.is_triple()): + if (bond45.is_single() or bond45.is_double()): + bond12.decrement_order() + bond23.increment_order() + bond34.decrement_order() + bond45.increment_order() + structure = mol.copy(deep=True) + bond12.increment_order() + bond23.decrement_order() + bond34.increment_order() + bond45.decrement_order() + try: + structure.update_atomtypes(log_species=False) + except AtomTypeError: + pass + else: + structures.append(structure) + return structures \ No newline at end of file