diff --git a/documentation/source/users/rmg/input.rst b/documentation/source/users/rmg/input.rst
index 5edf1c572e..67c5f335a4 100644
--- a/documentation/source/users/rmg/input.rst
+++ b/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,
diff --git a/rmgpy/constraints.py b/rmgpy/constraints.py
index d3781004de..548de5ecd5 100644
--- a/rmgpy/constraints.py
+++ b/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:
diff --git a/rmgpy/constraintsTest.py b/rmgpy/constraintsTest.py
index 20479befb9..e4a20560af 100644
--- a/rmgpy/constraintsTest.py
+++ b/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.
diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py
index 8a906a0bcf..50f181ef39 100644
--- a/rmgpy/rmg/input.py
+++ b/rmgpy/rmg/input.py
@@ -813,6 +813,7 @@ def generated_species_constraints(**kwargs):
         'maximumNitrogenAtoms',
         'maximumSiliconAtoms',
         'maximumSulfurAtoms',
+        'maximumSurfaceSites',
         'maximumHeavyAtoms',
         'maximumRadicalElectrons',
         'maximumSingletCarbenes',