Merge pull request #2511 to add resonance for surface adsorbates

Adsorbate molecules can exhibit resonance similar to gas-phase molecules. 
For example, the bidentate adsorbate XCXC (where X is a surface site) has the following three configurations:
[X]C#C[X] <-> [X]#CC#[X] <-> [X]=C=C=[X]

In this pull request, @bjkreitz  added functions for the generation of resonance structures for adsorbates.

Resonance occurs only for adsorbates that have at least two binding sites, also known as bidentate adsorbates.
Therefore, a function was added to molecule.py to check if an adsorbate has two binding sites.
The pathfinder.py was updated with functions to find the structures of adsorbates which can eventually have resonance structures. The structures that are searched for in the adsorbates are X-C-C-X and X-C-C-C-X.
Five types of resonance for adsorbates were added to resonance.py. The identified and implemented resonance types are as follows:

Shift 2 electrons from a C=/#C bond to the X-C bond (generate_adsorbate_shift_down_resonance_structures)
Shift 2 electrons from a X=/#C bond to a C-C bond (generate_adsorbate_shift_up_resonance_structures)
Shift 4 electrons from a C#C bond to a X-C bond (generate_adsorbate_double_shift_down_resonance_structures)
Shift 4 electrons from a X#C bond to a C-C bond (generate_adsorbate_double_shift_up_resonance_structures)
Shift 2 electrons in a conjugate pi system for bridged X-C-C-C-X adsorbates (generate_adsorbate_conjugate_resonance_structures)

(this merge commit message is based on the pull request description, which may be out of date. 
See the pull request at #2511 for details)

rmgpy/data/thermo.py

@@ -1560,91 +1560,112 @@ def get_thermo_data_for_surface_species(self, species):
         Returns a :class:`ThermoData` object, with no Cp0 or CpInf
         """
 
+        # define the comparison function to find the lowest energy
+        def species_enthalpy(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 molecule in species.molecule:
+            # 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=species_enthalpy)
             else:
-                return species.thermo.get_enthalpy(298.0)
+                gas_phase_species = min(gas_phase_species_estimates, key=species_enthalpy)
+
+            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 %s: %r", gas_phase_species.label, thermo)
+
+            if not isinstance(thermo, ThermoData):
+                thermo = thermo.to_thermo_data()
+                find_cp0_and_cpinf(gas_phase_species, thermo)
+
+            # 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)"),
+            )
 
-        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)
+            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
 
-        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))
+            # (group_additivity=True means it appends the comments)
+            add_thermo_data(thermo, adsorption_thermo, group_additivity=True)
 
-        if not isinstance(thermo, ThermoData):
-            thermo = thermo.to_thermo_data()
-            find_cp0_and_cpinf(species, thermo)
+            # if the molecule had labels, reapply them
+            for label, atom in labeled_atoms.items():
+                if isinstance(atom,list):
+                    for a in atom:
+                        a.label = label
+                else:
+                    atom.label = label
 
-        # 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)"),
-        )
+            # a tuple of molecule and its thermo
+            resonance_data.append((molecule, thermo))
 
-        surface_sites = molecule.get_surface_sites()
-        try:
-            self._add_adsorption_correction(adsorption_thermo, self.groups[self.adsorption_groups], 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
+        # 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, and the preferred representation.
+
+        resonance_data = sorted(resonance_data, key=lambda x: x[1].H298.value_si)
+
+        # reorder the resonance structures (molecules) by their H298
+        species.molecule = [mol for mol, thermo in resonance_data]
 
-        # (group_additivity=True means it appends the comments)
-        add_thermo_data(thermo, adsorption_thermo, group_additivity=True)
+        thermo = resonance_data[0][1]
 
         if thermo.label:
             thermo.label += 'X' * len(surface_sites)
 
         find_cp0_and_cpinf(species, thermo)
 
-        # if the molecule had labels, reapply them 
-        for label,atom in labeled_atoms.items():
-            if isinstance(atom,list):
-                for a in atom:
-                    a.label = label
-            else:
-                atom.label = label
-
         return thermo
 
     def _add_adsorption_correction(self, adsorption_thermo, adsorption_groups, molecule, surface_sites):
@@ -2574,7 +2595,7 @@ def _add_group_thermo_data(self, thermo_data, database, molecule, atom):
         node = node0
         while node is not None and node.data is None:
             node = node.parent
-        if node is None: 
+        if node is None:
             raise DatabaseError(f'Unable to determine thermo parameters for atom {atom} in molecule {molecule}: '
                                 f'no data for node {node0} or any of its ancestors in database {database.label}.\n' +
                                 molecule.to_adjacency_list())
@@ -2872,7 +2893,6 @@ def register_in_central_thermo_db(self, species):
         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.
@@ -2883,3 +2903,4 @@ def find_cp0_and_cpinf(species, heat_capacity):
     if heat_capacity.CpInf is None:
         cp_inf = species.calculate_cpinf()
         heat_capacity.CpInf = (cp_inf, "J/(mol*K)")
+

rmgpy/molecule/fragment.py

@@ -1240,6 +1240,19 @@ def check_in_ring(self, rd_mol, mapped_atom_idx):
                 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
+
     def from_smiles_like_string(self, smiles_like_string):
         smiles = smiles_like_string
 

rmgpy/molecule/molecule.pxd

@@ -301,6 +301,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

@@ -883,6 +883,13 @@ def is_quadruple(self):
         """
         return self.is_order(4)
 
+    def is_double_or_triple(self):
+        """
+        Return ``True`` if the bond represents a double or triple bond or ``False``
+        if not.
+        """
+        return self.is_order(2) or self.is_order(3)
+
     def is_benzene(self):
         """
         Return ``True`` if the bond represents a benzene bond or ``False`` if
@@ -2864,6 +2871,16 @@ def get_surface_sites(self):
         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,58 @@ 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 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 or nitrogen atoms.
+    """
+    cython.declare(paths=list, atom2=Vertex, atom3=Vertex, atom4=Vertex, bond12=Edge, bond23=Edge, bond34=Edge)
+
+    paths = []
+    if atom1.is_surface_site():
+        for atom2, bond12 in atom1.edges.items():
+            if atom2.is_carbon() or atom2.is_nitrogen():
+                for atom3, bond23 in atom2.edges.items():
+                    if atom3.is_carbon() or atom3.is_nitrogen():
+                        for atom4, bond34 in atom3.edges.items():
+                            if atom4.is_surface_site():
+                                paths.append([atom1, atom2, atom3, atom4, bond12, bond23, bond34])
+    return paths
+
+
+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 or nitrogen atoms.
+    """
+
+    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_carbon() or atom2.is_nitrogen():
+                for atom3, bond23 in atom2.edges.items():
+                    if atom3.is_carbon() or atom3.is_nitrogen():
+                        for atom4, bond34 in atom3.edges.items():
+                            if atom2 is not atom4 and (atom4.is_carbon() or atom4.is_nitrogen()):
+                                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

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)