diff --git a/cmeutils/dynamics.py b/cmeutils/dynamics.py
index a0d164b..32b6a4a 100644
--- a/cmeutils/dynamics.py
+++ b/cmeutils/dynamics.py
@@ -30,7 +30,7 @@ def msd_from_gsd(
         Choose from "window" or "direct". See Freud for the differences
         https://freud.readthedocs.io/en/latest/modules/msd.html#freud.msd.MSD
     """
-    with gsd.hoomd.open(gsdfile, "rb") as trajectory:
+    with gsd.hoomd.open(gsdfile, "r") as trajectory:
         init_box = trajectory[start].configuration.box
         final_box = trajectory[stop].configuration.box
         assert all(
diff --git a/cmeutils/gsd_utils.py b/cmeutils/gsd_utils.py
index 569f094..c1ba685 100644
--- a/cmeutils/gsd_utils.py
+++ b/cmeutils/gsd_utils.py
@@ -21,12 +21,12 @@ def get_type_position(
     ----------
     typename : str or list of str
         Name of particles of which to get the positions (found in
-        gsd.hoomd.Snapshot.particles.types)
+        gsd.hoomd.Frame.particles.types)
         If you want the positions of multiple types, pass in a list
         e.g., ['ca', 'c3']
     gsd_file : str, default None
         Filename of the gsd trajectory
-    snap : gsd.hoomd.Snapshot, default None
+    snap : gsd.hoomd.Frame, default None
         Trajectory snapshot
     gsd_frame : int, default -1
         Frame number to get positions from. Supports negative indexing.
@@ -70,7 +70,7 @@ def get_all_types(gsd_file=None, snap=None, gsd_frame=-1):
     ----------
     gsd_file : str, default None
         Filename of the gsd trajectory
-    snap : gsd.hoomd.Snapshot, default None
+    snap : gsd.hoomd.Frame, default None
         Trajectory snapshot
     gsd_frame : int, default -1
         Frame number to get positions from. Supports negative indexing.
@@ -94,7 +94,7 @@ def get_molecule_cluster(gsd_file=None, snap=None, gsd_frame=-1):
     ----------
     gsd_file : str, default None
         Filename of the gsd trajectory
-    snap : gsd.hoomd.Snapshot, default None
+    snap : gsd.hoomd.Frame, default None
         Trajectory snapshot.
     gsd_frame : int, default -1
         Frame number of gsd_file to use to compute clusters.
@@ -125,13 +125,13 @@ def _validate_inputs(gsd_file, snap, gsd_frame):
     if gsd_file:
         assert isinstance(gsd_frame, int)
         try:
-            with gsd.hoomd.open(name=gsd_file, mode="rb") as f:
+            with gsd.hoomd.open(name=gsd_file, mode="r") as f:
                 snap = f[gsd_frame]
         except Exception as e:
             print("Unable to open the gsd_file")
             raise e
     elif snap:
-        assert isinstance(snap, gsd.hoomd.Snapshot)
+        assert isinstance(snap, gsd.hoomd.Frame)
     return snap
 
 
@@ -148,15 +148,15 @@ def snap_delete_types(snap, delete_types):
 
     Parameters
     ----------
-    snap : gsd.hoomd.Snapshot
+    snap : gsd.hoomd.Frame
         The snapshot to read in
 
     Returns
     -------
-    gsd.hoomd.Snapshot
+    gsd.hoomd.Frame
         The new snapshot with particles deleted.
     """
-    new_snap = gsd.hoomd.Snapshot()
+    new_snap = gsd.hoomd.Frame()
     delete_ids = [snap.particles.types.index(i) for i in delete_types]
     selection = np.where(~np.isin(snap.particles.typeid, delete_ids))[0]
     new_snap.particles.N = len(selection)
@@ -225,7 +225,7 @@ def update_rigid_snapshot(snapshot, mb_compound):
 
     Parameters
     ----------
-    snapshot : gsd.hoomd.Snapshot
+    snapshot : gsd.hoomd.Frame
         The snapshot returned from create_hoomd_forcefield
         or create_hoomd_simulation in mBuild
     mb_compound : mbuild.Compound, required
@@ -307,7 +307,7 @@ def ellipsoid_gsd(gsd_file, new_file, lpar, lperp):
         Value of lperp of the ellipsoids
 
     """
-    with gsd.hoomd.open(new_file, "wb") as new_t:
+    with gsd.hoomd.open(new_file, "w") as new_t:
         with gsd.hoomd.open(gsd_file) as old_t:
             for snap in old_t:
                 snap.particles.type_shapes = [
@@ -352,7 +352,7 @@ def xml_to_gsd(xmlfile, gsdfile):
             overwrite=True,
         )
         hoomd.util.unquiet_status()
-        with gsd.hoomd.open(f.name) as t, gsd.hoomd.open(gsdfile, "wb") as newt:
+        with gsd.hoomd.open(f.name) as t, gsd.hoomd.open(gsdfile, "w") as newt:
             snap = t[0]
             bonds = snap.bonds.group
             bonds = bonds[np.lexsort((bonds[:, 1], bonds[:, 0]))]
diff --git a/cmeutils/structure.py b/cmeutils/structure.py
index 56a73b3..0b38457 100644
--- a/cmeutils/structure.py
+++ b/cmeutils/structure.py
@@ -37,7 +37,7 @@ def angle_distribution(
         Filename of the GSD trajectory.
     A_name, B_name, C_name : str
         Name(s) of particles that form the angle triplet
-        (found in gsd.hoomd.Snapshot.particles.types)
+        (found in gsd.hoomd.Frame.particles.types)
         They must be given in the same order as they form the angle
     start : int
         Starting frame index for accumulating bond lengths.
@@ -78,7 +78,7 @@ def angle_distribution(
     elif degrees and theta_max is None:
         theta_max = 180
 
-    trajectory = gsd.hoomd.open(gsd_file, mode="rb")
+    trajectory = gsd.hoomd.open(gsd_file, mode="r")
     name = "-".join([A_name, B_name, C_name])
     name_rev = "-".join([C_name, B_name, A_name])
 
@@ -149,7 +149,7 @@ def bond_distribution(
         Filename of the GSD trajectory.
     A_name, B_name : str
         Name(s) of particles that form the bond pair
-        (found in gsd.hoomd.Snapshot.particles.types)
+        (found in gsd.hoomd.Frame.particles.types)
     start : int
         Starting frame index for accumulating bond lengths.
         Negative numbers index from the end. (default 0)
@@ -179,7 +179,7 @@ def bond_distribution(
         If histogram is True, returns a 2D array of bin centers and bin heights.
 
     """
-    trajectory = gsd.hoomd.open(gsd_file, mode="rb")
+    trajectory = gsd.hoomd.open(gsd_file, mode="r")
     name = "-".join([A_name, B_name])
     name_rev = "-".join([B_name, A_name])
 
@@ -243,7 +243,7 @@ def dihedral_distribution(
         Filename of the GSD trajectory.
     A_name, B_name, C_name, D_name: str
         Name(s) of particles that form the dihedral quadruplett
-        (found in gsd.hoomd.Snapshot.particles.types)
+        (found in gsd.hoomd.Frame.particles.types)
         They must be given in the same order as they form the dihedral
     start : int
         Starting frame index for accumulating bond lengths.
@@ -273,7 +273,7 @@ def dihedral_distribution(
         If histogram is True, returns a 2D array of bin centers and bin heights.
 
     """
-    trajectory = gsd.hoomd.open(gsd_file, mode="rb")
+    trajectory = gsd.hoomd.open(gsd_file, mode="r")
     name = "-".join([A_name, B_name, C_name, D_name])
     name_rev = "-".join([D_name, C_name, B_name, A_name])
 
@@ -391,7 +391,7 @@ def gsd_rdf(
         Filename of the GSD trajectory.
     A_name, B_name : str
         Name(s) of particles between which to calculate the RDF (found in
-        gsd.hoomd.Snapshot.particles.types)
+        gsd.hoomd.Frame.particles.types)
     start : int
         Starting frame index for accumulating the RDF. Negative numbers index
         from the end. (default 0)
@@ -413,7 +413,7 @@ def gsd_rdf(
     -------
     (freud.density.RDF, float)
     """
-    with gsd.hoomd.open(gsdfile, mode="rb") as trajectory:
+    with gsd.hoomd.open(gsdfile, mode="r") as trajectory:
         snap = trajectory[0]
 
         if r_max is None:
@@ -482,13 +482,13 @@ def get_centers(gsdfile, new_gsdfile):
     new_gsdfile : str
         Filename of new GSD for centers.
     """
-    with gsd.hoomd.open(new_gsdfile, "wb") as new_traj, gsd.hoomd.open(
-        gsdfile, "rb"
+    with gsd.hoomd.open(new_gsdfile, "w") as new_traj, gsd.hoomd.open(
+        gsdfile, "r"
     ) as traj:
         snap = traj[0]
         cluster_idx = gsd_utils.get_molecule_cluster(snap=snap)
         for snap in traj:
-            new_snap = gsd.hoomd.Snapshot()
+            new_snap = gsd.hoomd.Frame()
             new_snap.configuration.box = snap.configuration.box
             f_box = freud.box.Box.from_box(snap.configuration.box)
             # Use the freud box to unwrap the particle positions
@@ -640,7 +640,7 @@ def all_atom_rdf(
     -------
     freud.density.RDF
     """
-    with gsd.hoomd.open(gsdfile, mode="rb") as trajectory:
+    with gsd.hoomd.open(gsdfile, mode="r") as trajectory:
         snap = trajectory[start]
         if r_max is None:
             # Use a value just less than half the maximum box length.
diff --git a/cmeutils/tests/base_test.py b/cmeutils/tests/base_test.py
index 9551568..f9a93ee 100644
--- a/cmeutils/tests/base_test.py
+++ b/cmeutils/tests/base_test.py
@@ -31,13 +31,13 @@ def gsdfile_images(self, tmp_path):
 
     @pytest.fixture
     def snap(self, gsdfile):
-        with gsd.hoomd.open(name=gsdfile, mode="rb") as f:
+        with gsd.hoomd.open(name=gsdfile, mode="r") as f:
             snap = f[-1]
         return snap
 
     @pytest.fixture
     def snap_bond(self, gsdfile_bond):
-        with gsd.hoomd.open(name=gsdfile_bond, mode="rb") as f:
+        with gsd.hoomd.open(name=gsdfile_bond, mode="r") as f:
             snap = f[-1]
         return snap
 
@@ -72,7 +72,7 @@ def p3ht_fresnel(self):
 
 def create_frame(i, add_bonds, images, seed=42):
     np.random.seed(seed)
-    s = gsd.hoomd.Snapshot()
+    s = gsd.hoomd.Frame()
     s.configuration.step = i
     s.particles.N = 5
     s.particles.types = ["A", "B"]
@@ -93,7 +93,7 @@ def create_frame(i, add_bonds, images, seed=42):
 
 
 def create_gsd(filename, add_bonds=False, images=False):
-    with gsd.hoomd.open(name=filename, mode="wb") as f:
+    with gsd.hoomd.open(name=filename, mode="w") as f:
         f.extend(
             [
                 create_frame(i, add_bonds=add_bonds, images=images)
diff --git a/cmeutils/tests/test_gsd.py b/cmeutils/tests/test_gsd.py
index a463c1a..7691a35 100644
--- a/cmeutils/tests/test_gsd.py
+++ b/cmeutils/tests/test_gsd.py
@@ -33,7 +33,7 @@
 class TestGSD(BaseTest):
     def test_ellipsoid_gsd(self, butane_gsd):
         ellipsoid_gsd(butane_gsd, "ellipsoid.gsd", 0.5, 1.0)
-        with gsd.hoomd.open(name="ellipsoid.gsd", mode="rb") as f:
+        with gsd.hoomd.open(name="ellipsoid.gsd", mode="r") as f:
             snap = f[-1]
         assert snap.particles.type_shapes[0]["type"] == "Ellipsoid"
 
diff --git a/cmeutils/visualize.py b/cmeutils/visualize.py
index dce6a8f..41e9def 100644
--- a/cmeutils/visualize.py
+++ b/cmeutils/visualize.py
@@ -31,13 +31,13 @@ def __init__(
         gsd_file : str, required
             Path to a GSD file to load
         frame : int, optional, default 0
-            The frame of the GSD file to load the gsd.hoomd.Snapshot
+            The frame of the GSD file to load the gsd.hoomd.Frame
         view_axis : np.ndarray (3,), optional, default (1, 0, 0)
             Sets the fresnel.camera attributes of camrea position and direction
         color_dict : dict, optional, default None
             Set colors for particle types
         diameter_scale : float, optional default 0.30
-            Scale the diameter values stored in gsd.hoomd.Snapshot
+            Scale the diameter values stored in gsd.hoomd.Frame
         height : float, optional default 10
             Sets the fresnel.camera height attriubute
             Acts like a zoom where larger values zooms out
@@ -62,7 +62,7 @@ def __init__(
         upwrap_positions: bool, optional, default False
             If True, the particle positions are unwrapped in the image
             This requires the GSD file snapshot contain accurate values for
-            gsd.hoomd.Snapshot.particles.image
+            gsd.hoomd.Frame.particles.image
         device, fresnel.Device(), optional
             Set the device to be used by the scene and in rendering.
         show_box: bool, optional, default True
@@ -113,7 +113,7 @@ def frame(self, frame):
 
     @property
     def snapshot(self):
-        """gsd.hoomd.Snapshot loaded from the GSD file.
+        """gsd.hoomd.Frame loaded from the GSD file.
         The snapshot loaded depends on FresnelGSD.frame
         """
         return self._snapshot
@@ -154,7 +154,7 @@ def set_type_color(self, particle_type, color):
         """
         if particle_type not in set(self.particle_types):
             raise ValueError(
-                f"Particle type of {particle_type} is not in the Snapshot"
+                f"Particle type of {particle_type} is not in the Frame"
             )
         self._color_dict[particle_type] = color
 
@@ -163,7 +163,7 @@ def unwrap_positions(self):
         """If set to True, then positions of the particles are
         unwrapped before creating the image.
         This requires that the GSD file snapshots contain accurate
-        image values (gsd.hoomd.Snapshot.particles.image
+        image values (gsd.hoomd.Frame.particles.image
 
         """
         return self._unwrap_positions
@@ -290,7 +290,7 @@ def positions(self):
     @property
     def radius(self):
         """Sets the size of the particles.
-        Determined by the gsd.hoomd.Snapshot.particles.diameter
+        Determined by the gsd.hoomd.Frame.particles.diameter
         values and FresnelGSD.diameter_scale
         """
         return self.snapshot.particles.diameter * self.diameter_scale
diff --git a/environment-dev.yml b/environment-dev.yml
index c8d9243..582cb24 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -3,8 +3,8 @@ channels:
   - conda-forge
 dependencies:
   - fresnel
-  - freud
-  - gsd=2.9
+  - freud>=2.13.1
+  - gsd>=3.0
   - hoomd<4.0=*cpu*
   - mbuild
   - numpy
diff --git a/environment.yml b/environment.yml
index 519a469..ef6cecc 100644
--- a/environment.yml
+++ b/environment.yml
@@ -3,8 +3,8 @@ channels:
   - conda-forge
 dependencies:
   - fresnel
-  - freud
-  - gsd=2.9
+  - freud>=2.13.1
+  - gsd>=3.0
   - hoomd<4.0=*cpu*
   - mbuild
   - numpy