diff --git a/mbuild/path/__init__.py b/mbuild/path/__init__.py
index cbc39663e..2fa2bce29 100644
--- a/mbuild/path/__init__.py
+++ b/mbuild/path/__init__.py
@@ -8,4 +8,5 @@
     Path,
     Spiral2D,
     StraightLine,
+    ZigZag,
 )
diff --git a/mbuild/utils/volumes.py b/mbuild/path/constraints.py
similarity index 94%
rename from mbuild/utils/volumes.py
rename to mbuild/path/constraints.py
index 04b315d0a..00224b743 100644
--- a/mbuild/utils/volumes.py
+++ b/mbuild/path/constraints.py
@@ -27,13 +27,20 @@ def is_inside(self, points, buffer):
         """Return True if point satisfies constraint (inside), else False."""
         raise NotImplementedError("Must be implemented in subclasses")
 
-    def sample_candidates(self, points, n_candiates, buffer):
+    def sample_candidates(self, points, n_candiates, buffer, region_mins, region_maxs):
         """Sample the volume for canidate points sorted by lowest local density."""
         raise NotImplementedError("Must be implemented in subclasses")
 
-    def find_low_density_points(self, points, n_candidates, buffer, k=10):
+    def find_low_density_points(
+        self, points, n_candidates, buffer, k=10, region_mins=None, region_maxs=None
+    ):
         low_density_points = self.sample_candidates(
-            points=points, n_candidates=n_candidates, buffer=buffer, k=k
+            points=points,
+            n_candidates=n_candidates,
+            buffer=buffer,
+            region_mins=region_mins,
+            region_maxs=region_maxs,
+            k=k,
         )
         return low_density_points
 
@@ -89,7 +96,9 @@ def is_inside(self, points, buffer):
             pbc=self.pbc,
         )
 
-    def sample_candidates(self, points, n_candidates, buffer, k=10):
+    def sample_candidates(
+        self, points, n_candidates, buffer, k=10, region_mins=None, region_maxs=None
+    ):
         """Generate candidate points uniformly distributed inside the box,
         optionally ranked by lowest local density around existing points.
 
@@ -116,8 +125,10 @@ def sample_candidates(self, points, n_candidates, buffer, k=10):
             the nearest neighbor (lowest local density) appear first.
         """
         # Create random candidates inside the box to test and sample from
+        sample_mins = np.asarray(region_mins) if region_mins else self.mins
+        sample_maxs = np.asarray(region_maxs) if region_maxs else self.maxs
         candidates = np.random.uniform(
-            self.mins + buffer, self.maxs - buffer, size=(n_candidates, 3)
+            sample_mins + buffer, sample_maxs - buffer, size=(n_candidates, 3)
         )
         if points is None or len(points) == 0:
             return candidates
diff --git a/mbuild/path/path.py b/mbuild/path/path.py
index 30ea86665..b3b2dca53 100644
--- a/mbuild/path/path.py
+++ b/mbuild/path/path.py
@@ -11,8 +11,8 @@
 from scipy.interpolate import interp1d
 
 from mbuild import Compound
+from mbuild.path.constraints import CuboidConstraint, CylinderConstraint
 from mbuild.path.path_utils import check_path, random_coordinate
-from mbuild.utils.volumes import CuboidConstraint, CylinderConstraint
 
 logger = logging.getLogger(__name__)
 
diff --git a/mbuild/polymer.py b/mbuild/polymer.py
index f4e9a3321..0b42a82d9 100644
--- a/mbuild/polymer.py
+++ b/mbuild/polymer.py
@@ -358,7 +358,10 @@ def build(
                             v1 = this_part_head.pos - this_part_tail.pos
                             v1 /= np.linalg.norm(v1)
                             # v2 is the vector between the previous site pos and next site pos
-                            v2 = coordinates[site_count + 1] - coordinates[site_count - 1]
+                            v2 = (
+                                coordinates[site_count + 1]
+                                - coordinates[site_count - 1]
+                            )
                             v2 /= np.linalg.norm(v2)
                             normal = np.cross(v1, v2)
                             angle = np.arccos(
diff --git a/mbuild/simulation.py b/mbuild/simulation.py
index 08d562aed..9a0cda2d8 100644
--- a/mbuild/simulation.py
+++ b/mbuild/simulation.py
@@ -412,7 +412,7 @@ def hoomd_cap_displacement(
     )
     sim.set_integrator(method=displacement_capped, dt=dt)
     sim.run(n_steps)
-    sim.operations.integrator = None 
+    sim.operations.integrator = None
     sim.update_positions()
     sim._update_snapshot()
 
@@ -528,7 +528,7 @@ def hoomd_fire(
         sim.operations.integrator.reset()
 
     # Update particle positions, save latest state point snapshot
-    sim.operations.integrator = None 
+    sim.operations.integrator = None
     sim._update_snapshot()
     sim.update_positions()
 
diff --git a/mbuild/tests/test_path.py b/mbuild/tests/test_path.py
index 2586689be..cf0af8a84 100644
--- a/mbuild/tests/test_path.py
+++ b/mbuild/tests/test_path.py
@@ -12,6 +12,11 @@
     Spiral2D,
     StraightLine,
 )
+from mbuild.path.constraints import (
+    CuboidConstraint,
+    CylinderConstraint,
+    SphereConstraint,
+)
 from mbuild.path.path_utils import (
     local_density,
     target_density,
@@ -24,11 +29,6 @@
 )
 from mbuild.tests.base_test import BaseTest
 from mbuild.utils.geometry import bounding_box
-from mbuild.utils.volumes import (
-    CuboidConstraint,
-    CylinderConstraint,
-    SphereConstraint,
-)
 
 
 class TestPaths(BaseTest):
diff --git a/mbuild/tests/test_termination.py b/mbuild/tests/test_termination.py
index bd25195ff..bea793857 100644
--- a/mbuild/tests/test_termination.py
+++ b/mbuild/tests/test_termination.py
@@ -7,6 +7,7 @@
     HardSphereRandomWalk,
 )
 from mbuild.path.bias import TargetCoordinate
+from mbuild.path.constraints import CuboidConstraint
 from mbuild.path.termination import (
     ContourLength,
     EndToEndDistance,
@@ -19,9 +20,6 @@
     WithinCoordinate,
 )
 from mbuild.tests.base_test import BaseTest, radius_of_gyration
-from mbuild.utils.volumes import (
-    CuboidConstraint,
-)
 
 
 class TestTermination(BaseTest):