Merge pull request #2136 from cat_fixes for Catalysis fixes

Catalysis fixes - vdw radicals, bidentates, resonance

vdw radicals - We previously implemented a multiplicity check to ensure products in the reverse direction match the multiplicity of the template reactants . However, we are getting vdw products that do not obey multiplicity restriction in the forward direction (#2122) , so we also need to check this direction as well.

radicals on X - if resonance generation puts radicals/lone pairs/or a formal charge on X, RMG crashes. I added a commit on this PR to filter out these structures

bidentate vdw - We were getting bidentate vdw structures. This PR makes them forbidden for surface families

rmgpy/data/kinetics/family.py

@@ -1521,6 +1521,20 @@ def apply_recipe(self, reactant_structures, forward=True, unique=True):
                 atom_labels['*2'].label = '*3'
                 atom_labels['*3'].label = '*2'
 
+            elif label == 'surface_abstraction':
+                atom_labels['*1'].label = '*3'
+                atom_labels['*2'].label = '*5'
+                atom_labels['*3'].label = '*1'
+                atom_labels['*5'].label = '*2'
+
+            elif label == 'surface_abstraction_single_vdw':
+                # *3 migrates from *2-*1 to *4-*5
+                # so swap *1 with *5, swap *2 with *4
+                atom_labels['*1'].label = '*5'
+                atom_labels['*5'].label = '*1'
+                atom_labels['*2'].label = '*4'
+                atom_labels['*4'].label = '*2'
+
         if not forward:
             template = self.reverse_template
             product_num = self.reactant_num or len(template.products)
@@ -1602,12 +1616,16 @@ def apply_recipe(self, reactant_structures, forward=True, unique=True):
                 lowest_labels.append(min(labels))
             product_structures = [s for _, s in sorted(zip(lowest_labels, product_structures))]
 
-        # If applying the family in reverse and the template reactants restrict multiplicity,
-        # we need to make sure that a reverse product matches the multiplicity-constrained
-        # template reactant because the forward template products do not enforce the
-        # multiplicity restriction.
-        if not forward and isinstance(reactant_structure, Molecule):
-            for template in self.forward_template.reactants:
+        # If the template restricts multiplicity, we need to make sure that
+        # a product matches the multiplicity-constrained template
+        # because the template does not ensure that the
+        # multiplicity restriction is obeyed
+        if isinstance(reactant_structure, Molecule):
+            if forward:
+                template_groups = self.forward_template.products
+            else:
+                template_groups = self.forward_template.reactants
+            for template in template_groups:
                 # iterate through the template reactants and check to see if they have a multiplicity constraint
                 if isinstance(template.item, Group):
                     if template.item.multiplicity != []:
@@ -1696,6 +1714,13 @@ def is_molecule_forbidden(self, molecule):
         if self.forbidden is not None and self.forbidden.is_molecule_forbidden(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
+
         return False
 
     def _create_reaction(self, reactants, products, is_forward):

rmgpy/data/kinetics/familyTest.py

@@ -867,6 +867,19 @@ def test_debug_forbidden_reverse_rxn(self, mock_logging):
                       'used as a reactant in the reverse direction.'),
         ])
 
+    def test_molecule_forbidden(self):
+
+        forbidden_mol = Molecule(smiles='*CC.[*]') # vdw bidentate
+
+        mol1 = Molecule(smiles='*CC*') # bidentate
+        mol2 = Molecule(smiles='C.*') # vdw
+        mol3 = Molecule(smiles='CC*') # chemisorbed
+
+        fam = self.database.kinetics.families['Surface_Dissociation_vdW']
+        self.assertTrue(fam.is_molecule_forbidden(forbidden_mol))
+        for allowed_mol in (mol1, mol2, mol3):
+            self.assertFalse(fam.is_molecule_forbidden(allowed_mol))
+
     def test_add_atom_labels_for_reaction(self):
         """Test that we can add atom labels to an existing reaction"""
         reactants = [Species().from_smiles('C=C'), Species().from_smiles('[OH]')]

rmgpy/molecule/resonance.py

@@ -182,11 +182,12 @@ def generate_resonance_structures(mol, clar_structures=True, keep_isomorphic=Fal
                       " definition.".format(mol.to_adjacency_list()))
         raise
     if mol.get_net_charge() != 0:
-        raise ValueError("Got the following structure:\nSMILES: {0}\nAdjacencyList:\n{1}\nNet charge: {2}\n\n"
-                         "Currently RMG cannot process charged species correctly."
-                         "\nIf this structure was entered in SMILES, try using the adjacencyList format for an"
-                         " unambiguous definition.".format(mol.to_smiles(), mol.to_adjacency_list(), mol.get_net_charge()))
-
+        # logging.info("Got the following structure:\nSMILES: {0}\nAdjacencyList:\n{1}\nNet charge: {2}\n\n"
+        #                  "Currently RMG cannot process charged species correctly."
+        #                  "\nIf this structure was entered in SMILES, try using the adjacencyList format for an"
+        #                  " unambiguous definition. "
+        #                  "Returning the input mol".format(mol.to_smiles(), mol.to_adjacency_list(), mol.get_net_charge()))
+        return [mol]
     if not mol.reactive:
         raise ResonanceError('Can only generate resonance structures for reactive molecules! Got the following '
                              'unreactive structure:\n{0}Reactive = {1}'.format(mol.to_adjacency_list(), mol.reactive))
@@ -253,7 +254,8 @@ def _generate_resonance_structures(mol_list, method_list, keep_isomorphic=False,
         copy                if False, append new resonance structures to input list (default)
                             if True, make a new list with all of the resonance structures
     """
-    cython.declare(index=cython.int, molecule=Molecule, new_mol_list=list, new_mol=Molecule, mol=Molecule)
+    cython.declare(index=cython.int, molecule=Molecule, new_mol_list=list, new_mol=Molecule, mol=Molecule, input_charge=cython.int,
+                    x=Atom)
 
     if copy:
         # Make a copy of the list so we don't modify the input list
@@ -299,11 +301,28 @@ def _generate_resonance_structures(mol_list, method_list, keep_isomorphic=False,
         index += 1
 
     # check net charge
-    for mol in mol_list:
-        if mol.get_net_charge() != 0:
-            raise ResonanceError('Resonance generation gave a net charged molecule:\n{0}'
-                                 'Ions are not yet supported in RMG.'.format(
-                mol.to_adjacency_list()))
+    input_charge = mol_list[0].get_net_charge()
+
+    for mol in mol_list[1:]:
+        if mol.get_net_charge() != input_charge:
+            mol_list.remove(mol)
+            logging.debug('Resonance generation created a molecule %s with a net charge of %d '
+                            'which does not match the input mol charge of %d.\n'
+                            'Removing it from resonance structures',mol.smiles,mol.get_net_charge(),input_charge)
+        if mol.contains_surface_site():
+            for x in [atom for atom in mol.atoms if atom.is_surface_site()]:
+                if x.radical_electrons != 0:
+                    mol_list.remove(mol)
+                    logging.debug('Resonance generation created a molecule %s with %d radicals on %s.\n'
+                    'Removing it from resonance structures',mol.smiles,x.radical_electrons,x.symbol)
+                elif x.lone_pairs != 0:
+                    mol_list.remove(mol)
+                    logging.debug('Resonance generation created a molecule %s with %d lone pairs on %s.\n'
+                    'Removing it from resonance structures',mol.smiles,x.lone_pairs,x.symbol)
+                elif x.charge != 0:
+                    mol_list.remove(mol)
+                    logging.debug('Resonance generation created a molecule %s with a charge of %d on %s.\n'
+                    'Removing it from resonance structures',mol.smiles,x.charge,x.symbol)
 
     return mol_list
 

rmgpy/molecule/resonanceTest.py

@@ -1419,6 +1419,13 @@ def test_surface_o(self):
 2 O u0 p2 c0 {1,D}"""))
         self.assertEquals(len(mol_list), 1)
 
+    def test_surface_radical(self):
+        """Test resonance structure generation for radical on surface
+
+        Should not put radical on X atom"""
+        mol_list = generate_resonance_structures(Molecule(smiles='*=C[CH]C'))
+        self.assertEquals(len(mol_list), 1)
+
     def test_sulfur_triple_bond(self):
         """
         Test the prevention of S#S formation through the find_lone_pair_multiplebond_paths and