Merge pull request #2069 from ReactionMechanismGenerator/surface_site…

…_constaint

The bidentate families can make molecules with a bunch of surface sites, like 3 or 4 or more. This PR adds a maximumSurfaceSites constraint if you want to limit the number of surface sites that can be in a species in your model. Recommended value is probably 2 (but default if unspecified is infinite.)

documentation/source/users/rmg/input.rst

@@ -909,6 +909,7 @@ all of RMG's reaction families. ::
         maximumSiliconAtoms=2,
         maximumSulfurAtoms=2,
         maximumHeavyAtoms=10,
+        maximumSurfaceSites=2,
         maximumRadicalElectrons=2,
         maximumSingletCarbenes=1,
         maximumCarbeneRadicals=0,

rmgpy/constraints.py

@@ -82,6 +82,11 @@ def fails_species_constraints(species):
         if struct.get_num_atoms('S') > max_sulfur_atoms:
             return True
 
+    max_surface_sites = species_constraints.get('maximumSurfaceSites', -1)
+    if max_surface_sites != -1:
+        if struct.get_num_atoms('X') > max_surface_sites:
+            return True
+
     max_heavy_atoms = species_constraints.get('maximumHeavyAtoms', -1)
     if max_heavy_atoms != -1:
         if struct.get_num_atoms() - struct.get_num_atoms('H') > max_heavy_atoms:

rmgpy/constraintsTest.py

@@ -61,6 +61,7 @@ def setUpClass(cls):
             maximumNitrogenAtoms=1,
             maximumSiliconAtoms=1,
             maximumSulfurAtoms=1,
+            maximumSurfaceSites=2,
             maximumHeavyAtoms=3,
             maximumRadicalElectrons=2,
             maximumSingletCarbenes=1,
@@ -159,6 +160,57 @@ def test_sulfur_constraint(self):
         mol2 = Molecule(smiles='SCS')
         self.assertTrue(fails_species_constraints(mol2))
 
+    def test_surface_site_constraint(self):
+        """
+        Test that we can constrain the max number of surface sites.
+        """
+
+        mol_1site = Molecule().from_adjacency_list("""
+1 O u0 p2 c0 {2,D}
+2 C u0 p0 c0 {1,D} {3,D}
+3 X u0 p0 c0 {2,D}
+""")
+        mol_2site = Molecule().from_adjacency_list("""
+1 C u0 p0 c0 {2,D} {3,D}
+2 C u0 p0 c0 {1,D} {4,D}
+3 X u0 p0 c0 {1,D}
+4 X u0 p0 c0 {2,D}
+""")
+
+        mol_3site_vdW = Molecule().from_adjacency_list("""
+1 C u0 p0 c0 {2,D} {3,D}
+2 C u0 p0 c0 {1,D} {4,D}
+3 X u0 p0 c0 {1,D}
+4 X u0 p0 c0 {2,D}
+6 X u0 p0 c0
+""")
+
+        mol_3site = Molecule().from_adjacency_list("""
+1 C u0 p0 c0 {4,S} {2,D} {7,S}
+2 C u0 p0 c0 {1,D} {3,S} {8,S}
+3 C u0 p0 c0 {2,S} {5,S} {6,S} {9,S}
+4 H u0 p0 c0 {1,S}
+5 H u0 p0 c0 {3,S}
+6 H u0 p0 c0 {3,S}
+7 X u0 p0 c0 {1,S}
+8 X u0 p0 c0 {2,S}
+9 X u0 p0 c0 {3,S}
+""")
+        max_carbon = self.rmg.species_constraints['maximumCarbonAtoms']
+        max_heavy_atoms = self.rmg.species_constraints['maximumHeavyAtoms']
+
+        self.rmg.species_constraints['maximumCarbonAtoms'] = 3
+        self.rmg.species_constraints['maximumHeavyAtoms'] = 6
+
+        self.assertFalse(fails_species_constraints(mol_1site))
+        self.assertFalse(fails_species_constraints(mol_2site))
+
+        self.assertTrue(fails_species_constraints(mol_3site_vdW))
+        self.assertTrue(fails_species_constraints(mol_3site))
+
+        self.rmg.species_constraints['maximumCarbonAtoms'] = max_carbon
+        self.rmg.species_constraints['maximumHeavyAtoms'] = max_heavy_atoms
+
     def test_heavy_constraint(self):
         """
         Test that we can constrain the max number of heavy atoms.

rmgpy/rmg/input.py

@@ -813,6 +813,7 @@ def generated_species_constraints(**kwargs):
         'maximumNitrogenAtoms',
         'maximumSiliconAtoms',
         'maximumSulfurAtoms',
+        'maximumSurfaceSites',
         'maximumHeavyAtoms',
         'maximumRadicalElectrons',
         'maximumSingletCarbenes',