diff --git a/mbuild/packing.py b/mbuild/packing.py
index e3c3647a5..71dd98e4e 100644
--- a/mbuild/packing.py
+++ b/mbuild/packing.py
@@ -668,6 +668,7 @@ def solvate(
     fix_orientation=False,
     temp_file=None,
     update_port_locations=False,
+    center_solute=True,
 ):
     """Solvate a compound in a box of solvent using PACKMOL.
 
@@ -701,6 +702,8 @@ def solvate(
     update_port_locations : bool, default=False
         After packing, port locations can be updated, but since compounds can be
         rotated, port orientation may be incorrect.
+    center_solute : bool, optional, default=True
+        Move solute center of mass to the center of the `mb.Box` used.
 
     Returns
     -------
@@ -724,8 +727,10 @@ def solvate(
     box_mins = np.asarray(min_tmp) * 10
     box_maxs = np.asarray(max_tmp) * 10
     overlap *= 10
-    center_solute = (box_maxs + box_mins) / 2
-
+    if center_solute:
+        center_solute = (box_maxs + box_mins) / 2
+    else:
+        center_solute = solute.pos * 10
     # Apply edge buffer
     box_maxs = np.subtract(box_maxs, edge * 10)
     box_mins = np.add(box_mins, edge * 10)
diff --git a/mbuild/tests/test_packing.py b/mbuild/tests/test_packing.py
index 5d20067ec..c1f79ec32 100644
--- a/mbuild/tests/test_packing.py
+++ b/mbuild/tests/test_packing.py
@@ -214,9 +214,24 @@ def test_fill_box_multiple(self, ethane, h2o):
 
     def test_solvate(self, ethane, h2o):
         n_solvent = 100
-        solvated = mb.solvate(ethane, h2o, n_solvent=n_solvent, box=[4, 4, 4])
-        assert solvated.n_particles == 8 + n_solvent * 3
-        assert solvated.n_bonds == 7 + n_solvent * 2
+        ethane_pos = ethane.pos
+        solvated1 = mb.solvate(
+            ethane, h2o, n_solvent=n_solvent, box=[4, 4, 4], center_solute=True
+        )
+        solvated2 = mb.solvate(
+            ethane, h2o, n_solvent=n_solvent, box=[4, 4, 4], center_solute=False
+        )
+
+        print(ethane.xyz)
+        print(solvated2.children[0].xyz)
+
+        assert (
+            solvated1.n_particles == solvated2.n_particles == 8 + n_solvent * 3
+        )
+        assert solvated1.n_bonds == solvated2.n_bonds == 7 + n_solvent * 2
+        assert not np.isclose(solvated1.children[0].pos, ethane_pos).all()
+        assert np.isclose(solvated2.children[0].pos, ethane_pos).all()
+
         assert ethane.parent == None
 
     def test_solvate_multiple(self, methane, ethane, h2o):