diff --git a/model.py b/model.py
index 1860c70..d15d094 100644
--- a/model.py
+++ b/model.py
@@ -219,7 +219,7 @@ def satisfiable(logic: Logic, model: Model, interpretation: Dict[Operation, Mode
 
 
 
-def model_closure(initial_set: Set[ModelValue], mfunctions: Set[ModelFunction], forbidden_element: Optional[ModelValue]) -> Set[ModelValue]:
+def model_closure(initial_set: Set[ModelValue], mfunctions: Set[ModelFunction], top: Optional[ModelValue], bottom: Optional[ModelValue]) -> Set[ModelValue]:
     """
     Given an initial set of model values and a set of model functions,
     compute the complete set of model values that are closed
@@ -230,7 +230,7 @@ def model_closure(initial_set: Set[ModelValue], mfunctions: Set[ModelFunction],
     closure_set: Set[ModelValue] = initial_set
     last_new: Set[ModelValue] = initial_set
     changed: bool = True
-    forbidden_found = False
+    topbottom_found = False
 
     while changed:
         changed = False
@@ -256,11 +256,14 @@ def model_closure(initial_set: Set[ModelValue], mfunctions: Set[ModelFunction],
 
                     # Optimization: Break out of computation
                     # early when top or bottom element is foun
-                    if forbidden_element is not None and element == forbidden_element:
-                        forbidden_found = True
+                    if top is not None and element == top:
+                        topbottom_found = True
+                    if bottom is not None and element == bottom:
+                        topbottom_found = True
+                    if topbottom_found:
                         break
 
-                if forbidden_found:
+                if topbottom_found:
                     break
 
                 # We don't need to compute the arguments
@@ -288,14 +291,17 @@ def model_closure(initial_set: Set[ModelValue], mfunctions: Set[ModelFunction],
 
                         # Optimization: Break out of computation
                         # early when top or bottom element is foun
-                        if forbidden_element is not None and element == forbidden_element:
-                            forbidden_found = True
+                        if top is not None and element == top:
+                            topbottom_found = True
+                        if bottom is not None and element == bottom:
+                            topbottom_found = True
+                        if topbottom_found:
                             break
 
-                    if forbidden_found:
+                    if topbottom_found:
                         break
 
-                if forbidden_found:
+                if topbottom_found:
                     break
 
 
@@ -303,7 +309,7 @@ def model_closure(initial_set: Set[ModelValue], mfunctions: Set[ModelFunction],
         changed = len(new_elements) > 0
         last_new = new_elements
 
-        if forbidden_found:
+        if topbottom_found:
             break
 
     return closure_set
diff --git a/vsp.py b/vsp.py
index 7420b68..5c00aff 100644
--- a/vsp.py
+++ b/vsp.py
@@ -145,28 +145,30 @@ def has_vsp(model: Model, interpretation: Dict[Operation, ModelFunction]) -> VSP
             continue
 
         # NOTE: Optimization
-        # If the left subalgebra contains bottom
-        # or the right subalgebra contains top
-        # skip this pair
-        if top is not None and top in ys:
+        # if either subalgebra contains top or bottom, move
+        # onto the next pair
+        if top is not None and (top in xs or top in ys):
             continue
-        if bottom is not None and bottom in xs:
+        if bottom is not None and (bottom in xs or bottom in ys):
             continue
 
         # Compute the closure of all operations
         # with just the xs
-        carrier_set_left: Set[ModelValue] = model_closure(xs, model.logical_operations, bottom)
+        carrier_set_left: Set[ModelValue] = model_closure(xs, model.logical_operations, top, bottom)
 
         # Save to cache
         if cached_xs[0] is not None and not cached_ys[1]:
             closure_cache.append((orig_xs, carrier_set_left))
 
+        if top is not None and top in carrier_set_left:
+            continue
         if bottom is not None and bottom in carrier_set_left:
             continue
 
+
         # Compute the closure of all operations
         # with just the ys
-        carrier_set_right: Set[ModelValue] = model_closure(ys, model.logical_operations, top)
+        carrier_set_right: Set[ModelValue] = model_closure(ys, model.logical_operations, top, bottom)
 
         # Save to cache
         if cached_ys[0] is not None and not cached_ys[1]:
@@ -174,6 +176,8 @@ def has_vsp(model: Model, interpretation: Dict[Operation, ModelFunction]) -> VSP
 
         if top is not None and top in carrier_set_right:
             continue
+        if bottom is not None and bottom in carrier_set_right:
+            continue
 
         # If the carrier set intersects, then move on to the next
         # subalgebra