adding resonance for adsorbates

rmgpy/data/thermo.py

@@ -1542,77 +1542,96 @@
         Returns a :class:`ThermoData` object, with no Cp0 or CpInf
         """
 
+        # define the comparison function to find the lowest energy
+
+        def lowest_energy(species):
+            if hasattr(species.thermo, 'H298'):
+                return species.thermo.H298.value_si
+            else:
+                return species.thermo.get_enthalpy(298.0)
+
         if species.is_surface_site():
             raise DatabaseError("Can't estimate thermo of vacant site. Should be in library (and should be 0).")
 
         logging.debug("Trying to generate thermo for surface species using first of %d resonance isomer(s):",
                       len(species.molecule))
-        molecule = species.molecule[0]
-        # store any labeled atoms to reapply at the end
-        labeled_atoms = molecule.get_all_labeled_atoms()
-        molecule.clear_labeled_atoms()
-        logging.debug("Before removing from surface:\n" + molecule.to_adjacency_list())
-        # only want/need to do one resonance structure,
-        # because will need to regenerate others in gas phase
-        dummy_molecules = molecule.get_desorbed_molecules()
-        for mol in dummy_molecules:
-            mol.clear_labeled_atoms()
-        if len(dummy_molecules) == 0:
-            raise RuntimeError(f"Cannot get thermo for gas-phase molecule. No valid dummy molecules from original molecule:\n{molecule.to_adjacency_list()}")
+
+        resonance_data = []
+
+        # iterate over all resonance structures of the species
+
+        for i in range(len(species.molecule)):
+            molecule = species.molecule[i]
+            # store any labeled atoms to reapply at the end
+            labeled_atoms = molecule.get_all_labeled_atoms()
+            molecule.clear_labeled_atoms()
+            logging.debug("Before removing from surface:\n" + molecule.to_adjacency_list())
+            # get all desorbed molecules for a given adsorbate
+            dummy_molecules = molecule.get_desorbed_molecules()
+            for mol in dummy_molecules:
+                mol.clear_labeled_atoms()
+            if len(dummy_molecules) == 0:
+                raise RuntimeError(f"Cannot get thermo for gas-phase molecule. No valid dummy molecules from original molecule:\n{molecule.to_adjacency_list()}")
 
-        # if len(molecule) > 1, it will assume all resonance structures have already been generated when it tries to generate them, so evaluate each configuration separately and pick the lowest energy one by H298 value
-        gas_phase_species_from_libraries = []
-        gas_phase_species_estimates = []
-        for dummy_molecule in dummy_molecules:
-            dummy_species = Species()
-            dummy_species.molecule = [dummy_molecule]
-            dummy_species.generate_resonance_structures()
-            dummy_species.thermo = self.get_thermo_data(dummy_species)
-            if dummy_species.thermo.label:
-                gas_phase_species_from_libraries.append(dummy_species)
-            else:
-                gas_phase_species_estimates.append(dummy_species)
+            # if len(molecule) > 1, it will assume all resonance structures have already been generated when it tries to generate them, so evaluate each configuration separately and pick the lowest energy one by H298 value
+            gas_phase_species_from_libraries = []
+            gas_phase_species_estimates = []
+            for dummy_molecule in dummy_molecules:
+                dummy_species = Species()
+                dummy_species.molecule = [dummy_molecule]
+                dummy_species.generate_resonance_structures()
+                dummy_species.thermo = self.get_thermo_data(dummy_species)
+                if dummy_species.thermo.label:
+                    gas_phase_species_from_libraries.append(dummy_species)
+                else:
+                    gas_phase_species_estimates.append(dummy_species)
 
-        # define the comparison function to find the lowest energy
-        def lowest_energy(species):
-            if hasattr(species.thermo, 'H298'):
-                return species.thermo.H298.value_si
+            if gas_phase_species_from_libraries:
+                gas_phase_species = min(gas_phase_species_from_libraries, key=lowest_energy)
             else:
-                return species.thermo.get_enthalpy(298.0)
+                gas_phase_species = min(gas_phase_species_estimates, key=lowest_energy)
 
-        if gas_phase_species_from_libraries:
-            species = min(gas_phase_species_from_libraries, key=lowest_energy)
-        else:
-            species = min(gas_phase_species_estimates, key=lowest_energy)
+            thermo = gas_phase_species.thermo
+            thermo.comment = f"Gas phase thermo for {thermo.label or gas_phase_species.molecule[0].to_smiles()} from {thermo.comment}. Adsorption correction:"
+            logging.debug("Using thermo from gas phase for species {}\n".format(gas_phase_species.label) + repr(thermo))
 
-        thermo = species.thermo
-        thermo.comment = f"Gas phase thermo for {thermo.label or species.molecule[0].to_smiles()} from {thermo.comment}. Adsorption correction:"
-        logging.debug("Using thermo from gas phase for species {}\n".format(species.label) + repr(thermo))
+            if not isinstance(thermo, ThermoData):
+                thermo = thermo.to_thermo_data()
+                find_cp0_and_cpinf(gas_phase_species, thermo)
 
-        if not isinstance(thermo, ThermoData):
-            thermo = thermo.to_thermo_data()
-            find_cp0_and_cpinf(species, thermo)
+            # Get the adsorption energy
+            # Create the ThermoData object
 
-        # Get the adsorption energy
-        # Create the ThermoData object
-        adsorption_thermo = ThermoData(
-            Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
-            Cpdata=([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "J/(mol*K)"),
-            H298=(0.0, "kJ/mol"),
-            S298=(0.0, "J/(mol*K)"),
-        )
+            adsorption_thermo = ThermoData(
+                Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"),
+                Cpdata=([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], "J/(mol*K)"),
+                H298=(0.0, "kJ/mol"),
+                S298=(0.0, "J/(mol*K)"),
+            )
 
-        surface_sites = molecule.get_surface_sites()
-        try:
-            self._add_adsorption_correction(adsorption_thermo, self.groups['adsorptionPt111'], molecule, surface_sites)
-        except (KeyError, DatabaseError):
-            logging.error("Couldn't find in adsorption thermo database:")
-            logging.error(molecule)
-            logging.error(molecule.to_adjacency_list())
-            raise
+            surface_sites = molecule.get_surface_sites()
+            try:
+                self._add_adsorption_correction(adsorption_thermo, self.groups['adsorptionPt111'], molecule, surface_sites)
+            except (KeyError, DatabaseError):
+                logging.error("Couldn't find in adsorption thermo database:")
+                logging.error(molecule)
+                logging.error(molecule.to_adjacency_list())
+                raise
+
+            # (group_additivity=True means it appends the comments)
+            add_thermo_data(thermo, adsorption_thermo, group_additivity=True)
+
+            resonance_data.append(thermo)
 
-        # (group_additivity=True means it appends the comments)
-        add_thermo_data(thermo, adsorption_thermo, group_additivity=True)
+        # Get the enthalpy of formation of every adsorbate at 298 K and determine the resonance structure with the lowest enthalpy of formation
+        # We assume that the lowest enthalpy of formation is the correct thermodynamic property for the adsorbate
+
+        enthalpy298 = []
+
+        for i in range(len(resonance_data)):
+            enthalpy298.append(resonance_data[i].get_enthalpy(298))
+
+        thermo = resonance_data[np.argmin(enthalpy298)]
 
         if thermo.label:
             thermo.label += 'X' * len(surface_sites)
@@ -2849,7 +2868,6 @@
         except ValueError:
             logging.info('Fail to generate inchi/smiles for species below:\n{0}'.format(species.to_adjacency_list()))
 
-
 def find_cp0_and_cpinf(species, heat_capacity):
     """
     Calculate the Cp0 and CpInf values, and add them to the HeatCapacityModel object.
@@ -2860,3 +2878,4 @@
     if heat_capacity.CpInf is None:
         cp_inf = species.calculate_cpinf()
         heat_capacity.CpInf = (cp_inf, "J/(mol*K)")
+

rmgpy/molecule/fragment.py

@@ -129,6 +129,9 @@
     def is_surface_site(self):
         return False
 
+    def is_multidentate(self):
+        return False
+
     def is_silicon(self):
         return False
 
@@ -2078,6 +2081,22 @@
                 return True
         return False
 
+    def is_multidentate(self):
+        """
+        Return ``True`` if the adsorbate contains at least two binding sites,
+        or ``False`` otherwise.
+        """
+
+        surface_sites = 0
+        for atom in self.vertices:
+            if atom.is_surface_site():
+                surface_sites+=1
+
+        if surface_sites>=2:
+            return True
+
+        return False
+
 # this variable is used to name atom IDs so that there are as few conflicts by 
 # using the entire space of integer objects
 atom_id_counter = -2**15

rmgpy/molecule/molecule.pxd

@@ -289,6 +289,8 @@ cdef class Molecule(Graph):
 
     cpdef list get_adatoms(self)
 
+    cpdef bint is_multidentate(self)
+
     cpdef list get_desorbed_molecules(self)
 
 cdef atom_id_counter

rmgpy/molecule/molecule.py

@@ -2759,6 +2759,16 @@
         cython.declare(atom=Atom)
         return [atom for atom in self.atoms if atom.is_surface_site()]
 
+    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
+        return False
+
     def get_adatoms(self):
         """
         Get a list of adatoms in the molecule.

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)

rmgpy/molecule/pathfinder.py

@@ -480,3 +480,57 @@
     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 multidentate adsorbates which have a bonding configuration X-C-C-X.
+
+    Examples:
+
+    - XCXC, XCHXCH, XCXCH, 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 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 atom3.is_carbon():
+                        for atom4, bond34 in atom3.edges.items():
+                            if atom4.is_surface_site():
+                                path.append([atom1, atom2, atom3, atom4, bond12, bond23, bond34])
+
+    return path
+
+def find_adsorbate_conjugate_delocalization_paths(atom1):
+    """
+    Find all multidentate adsorbates which have a bonding configuration X-C-C-C-X.
+
+    Examples:
+
+    - XCHCHXCH/XCHCHXC, 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, 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():
+            if atom2.is_carbon():
+                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_carbon():
+                                for atom5, bond45 in atom4.edges.items():
+                                    if atom5.is_surface_site():
+                                        path.append([atom1, atom2, atom3, atom4, atom5, bond12, bond23, bond34, bond45])
+
+    return path

rmgpy/molecule/resonance.pxd

@@ -63,3 +63,9 @@ 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_conjugate_resonance_structures(Graph mol)