diff --git a/cmeutils/__init__.py b/cmeutils/__init__.py
index 0d2eabd..6de519e 100644
--- a/cmeutils/__init__.py
+++ b/cmeutils/__init__.py
@@ -1,4 +1,10 @@
 from . import gsd_utils, structure
 from .__version__ import __version__
 
-__all__ = ["__version__", "gsd_utils", "structure", "dynamics", "visualize"]
+__all__ = [
+    "__version__",
+    "gsd_utils",
+    "structure",
+    "dynamics",
+    "visualize",
+]
diff --git a/cmeutils/__version__.py b/cmeutils/__version__.py
index 3e8d9f9..5b60188 100644
--- a/cmeutils/__version__.py
+++ b/cmeutils/__version__.py
@@ -1 +1 @@
-__version__ = "1.4.0"
+__version__ = "1.5.0"
diff --git a/cmeutils/geometry.py b/cmeutils/geometry.py
index c28eec0..7c72362 100644
--- a/cmeutils/geometry.py
+++ b/cmeutils/geometry.py
@@ -1,5 +1,49 @@
 import numpy as np
 from numpy.linalg import svd
+from rowan import vector_vector_rotation
+
+
+def get_quaternions(n_views=20):
+    """Generate quaternions over the specified number of views.
+
+    The first (n_view - 3) views will be the views even distributed on a sphere,
+    while the last three views will be the face-on, edge-on, and corner-on
+    views, respectively.
+
+    These quaternions are useful as input to `view_orientation` kwarg in
+    `freud.diffraction.Diffractometer.compute`.
+
+    Parameters
+    ----------
+    n_views : int, default 20
+        The number of views to compute.
+
+    Returns
+    -------
+    list of numpy.ndarray
+        Quaternions as (4,) arrays.
+    """
+    if n_views <= 3 or not isinstance(n_views, int):
+        raise ValueError("Please set n_views to an integer greater than 3.")
+    # Calculate points for even distribution on a sphere
+    ga = np.pi * (3 - 5**0.5)
+    theta = ga * np.arange(n_views - 3)
+    z = np.linspace(1 - 1 / (n_views - 3), 1 / (n_views - 3), n_views - 3)
+    radius = np.sqrt(1 - z * z)
+    points = np.zeros((n_views, 3))
+    points[:-3, 0] = radius * np.cos(theta)
+    points[:-3, 1] = radius * np.sin(theta)
+    points[:-3, 2] = z
+
+    # face on
+    points[-3] = np.array([0, 0, 1])
+    # edge on
+    points[-2] = np.array([0, 1, 1])
+    # corner on
+    points[-1] = np.array([1, 1, 1])
+
+    unit_z = np.array([0, 0, 1])
+    return [vector_vector_rotation(i, unit_z) for i in points]
 
 
 def get_backbone_vector(coordinates):
diff --git a/cmeutils/gsd_utils.py b/cmeutils/gsd_utils.py
index 7bb8bc0..ae88af5 100644
--- a/cmeutils/gsd_utils.py
+++ b/cmeutils/gsd_utils.py
@@ -8,6 +8,24 @@
 from boltons.setutils import IndexedSet
 
 
+def snapshot_to_graph(snap):
+    """Convert a HOOMD snapshot into a networkx bond graph.
+
+    Parameters
+    ----------
+    snap : gsd.hoomd.Frame, required
+        Snapshot (frame) containing particle positions and bond groups.
+    """
+    graph = nx.Graph()
+    positions = snap.particles.position
+    type_ids = snap.particles.typeid
+    types = snap.particles.types
+    for i in range(snap.particles.N):
+        graph.add_node(i, type=types[type_ids[i]], pos=positions[i])
+    graph.add_edges_from(snap.bonds.group)
+    return graph
+
+
 def frame_to_freud_system(frame, ref_length=None):
     """Creates a freud system given a gsd.hoomd.Frame.
 
@@ -335,6 +353,50 @@ def identify_snapshot_connections(snapshot):
     return snapshot
 
 
+def get_centers(gsdfile, new_gsdfile):
+    """Create a gsd file containing the molecule centers from an existing gsd
+     file.
+
+
+    This function calculates the centers of a trajectory given a GSD file
+    and stores them into a new GSD file just for centers. By default it will
+    calculate the centers of an entire trajectory.
+
+    Parameters
+    ----------
+    gsdfile : str
+        Filename of the GSD trajectory.
+    new_gsdfile : str
+        Filename of new GSD for centers.
+    """
+    with (
+        gsd.hoomd.open(new_gsdfile, "w") as new_traj,
+        gsd.hoomd.open(gsdfile, "r") as traj,
+    ):
+        snap = traj[0]
+        cluster_idx = get_molecule_cluster(snap=snap)
+        for snap in traj:
+            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
+            unwrapped_positions = f_box.unwrap(
+                snap.particles.position, snap.particles.image
+            )
+            uw_centers = []
+            for i in range(max(cluster_idx) + 1):
+                cluster_uw_pos = unwrapped_positions[np.where(cluster_idx == i)]
+                uw_centers.append(np.mean(cluster_uw_pos, axis=0))
+            uw_centers = np.stack(uw_centers)
+            new_snap.particles.position = f_box.wrap(uw_centers)
+            new_snap.particles.N = len(uw_centers)
+            new_snap.particles.types = ["A"]
+            new_snap.particles.image = f_box.get_images(uw_centers)
+            new_snap.particles.typeid = np.zeros(len(uw_centers))
+            new_snap.validate()
+            new_traj.append(new_snap)
+
+
 def _fill_connection_info(snapshot, connections, type_):
     p_types = snapshot.particles.types
     p_typeid = snapshot.particles.typeid
diff --git a/cmeutils/structure.py b/cmeutils/structure.py
index 66467d5..889cbc3 100644
--- a/cmeutils/structure.py
+++ b/cmeutils/structure.py
@@ -1,17 +1,18 @@
+import itertools
 import warnings
 
 import freud
 import gsd
 import gsd.hoomd
+import networkx as nx
 import numpy as np
-from rowan import vector_vector_rotation
 
 from cmeutils.geometry import (
     angle_between_vectors,
     dihedral_angle,
     get_plane_normal,
 )
-from cmeutils.gsd_utils import frame_to_freud_system, get_molecule_cluster
+from cmeutils.gsd_utils import frame_to_freud_system, snapshot_to_graph
 from cmeutils.plotting import get_histogram
 
 
@@ -406,69 +407,32 @@ def dihedral_distribution(
         return np.array(dihedrals)
 
 
-def get_quaternions(n_views=20):
-    """Get the quaternions for the specified number of views.
-
-    The first (n_view - 3) views will be the views even distributed on a sphere,
-    while the last three views will be the face-on, edge-on, and corner-on
-    views, respectively.
-
-    These quaternions are useful as input to `view_orientation` kwarg in
-    `freud.diffraction.Diffractometer.compute`.
-
-    Parameters
-    ----------
-    n_views : int, default 20
-        The number of views to compute.
-
-    Returns
-    -------
-    list of numpy.ndarray
-        Quaternions as (4,) arrays.
-    """
-    if n_views <= 3 or not isinstance(n_views, int):
-        raise ValueError("Please set n_views to an integer greater than 3.")
-    # Calculate points for even distribution on a sphere
-    ga = np.pi * (3 - 5**0.5)
-    theta = ga * np.arange(n_views - 3)
-    z = np.linspace(1 - 1 / (n_views - 3), 1 / (n_views - 3), n_views - 3)
-    radius = np.sqrt(1 - z * z)
-    points = np.zeros((n_views, 3))
-    points[:-3, 0] = radius * np.cos(theta)
-    points[:-3, 1] = radius * np.sin(theta)
-    points[:-3, 2] = z
-
-    # face on
-    points[-3] = np.array([0, 0, 1])
-    # edge on
-    points[-2] = np.array([0, 1, 1])
-    # corner on
-    points[-1] = np.array([1, 1, 1])
-
-    unit_z = np.array([0, 0, 1])
-    return [vector_vector_rotation(i, unit_z) for i in points]
-
-
 def gsd_rdf(
     gsdfile,
-    A_name,
-    B_name,
+    A_name=None,
+    B_name=None,
     start=0,
     stop=-1,
     stride=1,
     r_max=None,
     r_min=0,
     bins=100,
-    exclude_bonded=True,
+    exclude_bond_depth=None,
+    exclude_all_bonded=False,
 ):
-    """Compute intermolecular RDF from a GSD file.
+    """Uses freud's RDF module to calculate an RDF averaged over a GSD file.
 
-    This function calculates the radial distribution function given a GSD file
-    and the names of the particle types. By default it will calculate the RDF
-    for the entire trajectory.
+    Notes
+    -----
+    This method lets you set a bond-depth exclusion to prevent neighbors of the
+    same molecule to be used in the RDF calculations. Bond depth is counted in
+    terms of steps away on a connected bond graph.
 
-    It is assumed that the bonding, number of particles, and simulation box do
-    not change during the simulation.
+    For example, ``exclude_bond_depth=1`` excludes pairs directly bonded together,
+    ``exclude_bond_depth=2`` excludes both pairs (i, j) and (i, k) in (i-j-k), and so on.
+    This can be set to any positive integer value, and isn't limited to particles belonging
+    to the same bond, angle or dihedral. To exclude all intramolecular pairs use
+    the ``exclude_all_bonded`` parameter instead.
 
     Parameters
     ----------
@@ -486,71 +450,124 @@ def gsd_rdf(
     stride : int, default = 1
         The stride size when iterating through start:stop
     r_max : float
-        Maximum radius of RDF. If None, half of the maximum box size is used.
-        (default -1)
+        Maximum radius of RDF. If None, half of the minimum box size is used.
     r_min : float
         Minimum radius of RDF. (default 0)
     bins : int
         Number of bins to use when calculating the RDF. (default 100)
-    exclude_bonded : bool
-        Whether to remove particles in same molecule from the neighbor list.
-        (default True)
+    exclude_bond_depth : int, optional (default None)
+        Excludes all pairs within a depth (distance on a bond graph)
+        from the RDF calculation.
+    exclude_all_bonded : bool, optional (default False)
+        Excludes all pairs belonging to the same molecule from the
+        RDF calculation.
 
     Returns
     -------
-    (freud.density.RDF, float)
+    tuple(rdf, rdf_correction) : (freud.density.RDF, float)
+        Access r values with ``rdf.bin_centers`` and g(r) with ``rdf.rdf``
+        rdf_correction is always 1 unless ``exclude_bond_depth`` or ``exclude_all_bonded`` are used.
+        This corrects the g(r) normalization to account for excluded pairs.
+        To include this in the results, g(r) = rdf.rdf * rdf_correction
     """
+    if any([A_name, B_name]) and not all([A_name, B_name]):
+        raise ValueError(
+            "If A_name or B_name is given, the other must be defined as well. "
+            "To calculate an RDF between the same bead type, set A_name and B_name equal to the same value. "
+            "To calculate an RDF between all possible pairs, leave both as ``None``. "
+        )
+
+    if all([exclude_bond_depth, exclude_all_bonded]):
+        raise ValueError(
+            "Only use one of excluded_bond_depth and exclude_all_bonded."
+        )
+
     with gsd.hoomd.open(gsdfile, mode="r") as trajectory:
-        snap = trajectory[0]
+        # Grab the first snapshot for some book-keeping.
+        snap = trajectory[start]
 
+        # Use a value just less than half the minimum box length.
+        # TODO: Iterate through all boxes, use the smallest one? Edge cases might arise with NPT sims
         if r_max is None:
-            # Use a value just less than half the maximum box length.
             r_max = np.nextafter(
-                np.max(snap.configuration.box[:3]) * 0.49, 0, dtype=np.float32
+                np.min(snap.configuration.box[:3]) * 0.49, 0, dtype=np.float32
             )
 
         rdf = freud.density.RDF(bins=bins, r_max=r_max, r_min=r_min)
 
-        type_A = snap.particles.typeid == snap.particles.types.index(A_name)
-        type_B = snap.particles.typeid == snap.particles.types.index(B_name)
+        # Filter particles by type A and type B
+        if A_name is not None and B_name is not None:
+            type_A = snap.particles.typeid == snap.particles.types.index(A_name)
+            type_B = snap.particles.typeid == snap.particles.types.index(B_name)
+            # These 2 *_indices variables store global particle indices (Before filtering by type)
+            # If excluding by bond depth, these need to be passed into filter_nlist
+            type_A_indices = np.where(type_A)[0]
+            type_B_indices = np.where(type_B)[0]
+            exclude_ii = A_name == B_name
+        else:  # Use all particles for this RDF
+            type_A = type_B = np.ones(
+                snap.particles.N, dtype=bool
+            )  # Array of True at all indices
+            type_A_indices = type_B_indices = np.arange(snap.particles.N)
+            exclude_ii = True
+
+        # Build up pair exclusions if exclude_bond_depth or exclude_all_bonded
+        # Reuse these for each frame's RDF
+        # If the bonding topology isn't changing, we only need to get bond graph and excluded pairs once
+        rdf_correction = 1
+        if exclude_bond_depth or exclude_all_bonded:
+            max_idx = snap.particles.N
+            bond_graph = snapshot_to_graph(snap)
+            # This gives a seuquence of tuples [(1, 4), (5, 8)...(i, j)]
+            excluded_pairs = get_excluded_pairs(
+                bond_graph, exclude_bond_depth, exclude_all_bonded
+            )
+            # Map information of sequence of tuples to an array of unique ints.
+            # This is used for faster filtering in filter_nlist() (vectorized instead of for loop)
+            excluded_pairs_encoded = np.array(
+                [i * max_idx + j for i, j in excluded_pairs]
+            )
 
-        if exclude_bonded:
-            molecules = get_molecule_cluster(snap=snap)
-            molecules_A = molecules[type_A]
-            molecules_B = molecules[type_B]
+            n_excluded = len(excluded_pairs)
+            if A_name == B_name or not any(
+                [A_name, B_name]
+            ):  # Using same type or all particles
+                n_total_pairs = (
+                    len(type_A_indices) * (len(type_A_indices) - 1) / 2
+                )
+            else:  # RDF is not between same types, or using all particles
+                n_total_pairs = len(type_A_indices) * len(type_B_indices)
+            # Overwrite default value only if using exclude_bond_depth
+            if n_total_pairs == n_excluded:
+                warnings.warn(
+                    "Exclusions resulted in no pairs being used to calculate this RDF."
+                )
+            else:
+                rdf_correction = n_total_pairs / (n_total_pairs - n_excluded)
 
         for snap in trajectory[start:stop:stride]:
-            A_pos = snap.particles.position[type_A]
-            if A_name == B_name:
-                B_pos = A_pos
-                exclude_ii = True
-                ab_ratio = 1
-            else:
-                B_pos = snap.particles.position[type_B]
-                exclude_ii = False
-                ab_ratio = len(A_pos) / len(B_pos)
+            A_xyz = snap.particles.position[type_A_indices]
+            B_xyz = snap.particles.position[type_B_indices]
 
+            # Build up the complete neighborlist
             box = snap.configuration.box
-            system = (box, A_pos)
+            system = (box, A_xyz)
             aq = freud.locality.AABBQuery.from_system(system)
-            nlist = aq.query(
-                B_pos, {"r_max": r_max, "exclude_ii": exclude_ii}
-            ).toNeighborList()
-
-            if exclude_bonded:
-                pre_filter = len(nlist)
-                nlist.filter(
-                    molecules_A[nlist.point_indices]
-                    != molecules_B[nlist.query_point_indices]
+            query_args = {"r_max": r_max, "exclude_ii": exclude_ii}
+            nlist = aq.query(B_xyz, query_args).toNeighborList()
+
+            # Filter excluded pairs from all pairs in nlist
+            if exclude_bond_depth or exclude_all_bonded:
+                nlist = filter_nlist(
+                    nlist=nlist,
+                    excluded_pairs_encoded=excluded_pairs_encoded,
+                    query_indices=type_A_indices,
+                    point_indices=type_B_indices,
+                    max_idx=snap.particles.N,
                 )
-                post_filter = len(nlist)
 
             rdf.compute(aq, neighbors=nlist, reset=False)
-
-        normalization = post_filter / pre_filter if exclude_bonded else 1
-        normalization *= ab_ratio
-
-        return rdf, normalization
+    return rdf, rdf_correction
 
 
 def structure_factor(
@@ -680,50 +697,6 @@ def diffraction_pattern(
     return dp
 
 
-def get_centers(gsdfile, new_gsdfile):
-    """Create a gsd file containing the molecule centers from an existing gsd
-     file.
-
-
-    This function calculates the centers of a trajectory given a GSD file
-    and stores them into a new GSD file just for centers. By default it will
-    calculate the centers of an entire trajectory.
-
-    Parameters
-    ----------
-    gsdfile : str
-        Filename of the GSD trajectory.
-    new_gsdfile : str
-        Filename of new GSD for centers.
-    """
-    with (
-        gsd.hoomd.open(new_gsdfile, "w") as new_traj,
-        gsd.hoomd.open(gsdfile, "r") as traj,
-    ):
-        snap = traj[0]
-        cluster_idx = get_molecule_cluster(snap=snap)
-        for snap in traj:
-            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
-            unwrapped_positions = f_box.unwrap(
-                snap.particles.position, snap.particles.image
-            )
-            uw_centers = []
-            for i in range(max(cluster_idx) + 1):
-                cluster_uw_pos = unwrapped_positions[np.where(cluster_idx == i)]
-                uw_centers.append(np.mean(cluster_uw_pos, axis=0))
-            uw_centers = np.stack(uw_centers)
-            new_snap.particles.position = f_box.wrap(uw_centers)
-            new_snap.particles.N = len(uw_centers)
-            new_snap.particles.types = ["A"]
-            new_snap.particles.image = f_box.get_images(uw_centers)
-            new_snap.particles.typeid = np.zeros(len(uw_centers))
-            new_snap.validate()
-            new_traj.append(new_snap)
-
-
 def order_parameter(aa_gsd, cg_gsd, mapping, r_max, a_max, large=6, start=-10):
     """Calculate the order parameter of a system.
 
@@ -887,54 +860,83 @@ def concentration_profile(snap, A_indices, B_indices, n_bins=70, box_axis=0):
     return d_profile[:-1], A_count, B_count, total_count
 
 
-def all_atom_rdf(
-    gsdfile,
-    start=0,
-    stop=-1,
-    stride=1,
-    r_max=None,
-    r_min=0,
-    bins=100,
+def get_excluded_pairs(
+    bond_graph, excluded_bond_depth=None, exclude_all_bonded=False
 ):
-    """Compute intermolecular RDF from a GSD file.
+    """Returns a set of (i, j) pairs to exclude based on step distance of a bond graph."""
+    excluded_pairs = set()
+    if excluded_bond_depth:
+        for i in bond_graph.nodes:
+            lengths = nx.single_source_shortest_path_length(
+                bond_graph, i, cutoff=excluded_bond_depth
+            )
+            for j, dist in lengths.items():
+                # use j > i for 2 reasons: skip adding the same pair twice and its used in filter_nlist
+                if j > i and dist <= excluded_bond_depth:
+                    excluded_pairs.add((i, j))
+    elif exclude_all_bonded:
+        for component in nx.connected_components(bond_graph):
+            component = sorted(component)
+            for i, j in itertools.combinations(component, 2):
+                excluded_pairs.add((i, j))
+    return excluded_pairs
+
+
+def filter_nlist(
+    nlist, excluded_pairs_encoded, query_indices, point_indices, max_idx
+):
+    """Filter a freud NeighborList by removing excluded pairs.
 
-    This function calculates the radial distribution function given a GSD file
-    for all atoms. By default it will calculate the RDF
-    for the entire trajectory.
-    It is assumed that the bonding, number of particles, and simulation box do
-    not change during the simulation.
+    Handles the index space mismatch between the NeighborList (which uses
+    local indices into type filtered position arrays) and excluded_pairs
+    (which uses global particle indices from the bond graph).
 
     Parameters
     ----------
-    gsdfile : str
-        Filename of the GSD trajectory.
-    start : int, default 0
-        Starting frame index for accumulating the RDF. Negative numbers index
-        from the end.
-    stop : int, default -1
-        Final frame index for accumulating the RDF. If None, the last frame
-        will be used.
-    stride : int, default 1
-        The stride size when iterating through start:stop
-    r_max : float, default None
-        Maximum radius of RDF. If None, half of the maximum box size is used.
-    r_min : float, default 0
-        Minimum radius of RDF.
-    bins : int, default 100
-        Number of bins to use when calculating the RDF.
+    nlist : freud.locality.NeighborList
+        Neighbor list to filter, as returned by AABBQuery.toNeighborList().
+        Query point and point indices are local (0 indexed within their
+        type subsets).
+    excluded_pairs_encoded : np.ndarray of int, shape (N_excluded,)
+        Encoded global particle index pairs to exclude, precomputed as
+        ``i * max_idx + j`` where ``i < j`` and ``max_idx = N_particles``.
+        Encodes the output of get_excluded_pairs().
+    query_indices : np.ndarray of int, shape (N_query,)
+        Global particle indices of the query points (type A particles).
+        Maps local query index to global particle index.
+    point_indices : np.ndarray of int, shape (N_points,)
+        Global particle indices of the reference points (type B particles).
+        Maps local point index to global particle index.
 
     Returns
     -------
-    freud.density.RDF
+    freud.locality.NeighborList
+        Filtered neighbor list with excluded pairs removed. Local indices
+        are preserved, so it can be passed directly to freud.density.RDF.compute().
     """
-    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.
-            r_max = np.nextafter(
-                np.max(snap.configuration.box[:3]) * 0.5, 0, dtype=np.float32
-            )
-        rdf = freud.density.RDF(bins=bins, r_max=r_max, r_min=r_min)
-        for snap in trajectory[start:stop:stride]:
-            rdf.compute(snap, reset=False)
-    return rdf
+    # Local indices of the neighborlist, since it was created after filtering particle type
+    i_local = nlist.query_point_indices
+    j_local = nlist.point_indices
+
+    # excluded_pair indices are global, they came from the the bond graph of all particles
+    # Convert i_local and j_local to global indices so we can look up against excluded_pairs
+    i_global = query_indices[i_local]
+    j_global = point_indices[j_local]
+
+    # Sort all pairs so that the first index is always < the second.
+    lo = np.minimum(i_global, j_global)
+    hi = np.maximum(i_global, j_global)
+
+    # Encode each pair as a single integer for vectorized lookup.
+    # (lo, hi) = lo * max_idx + hi which is unique as long as hi < max_idx
+    neighbor_pairs_encoded = lo * max_idx + hi
+
+    # vectorized set lookup
+    keep = ~np.isin(neighbor_pairs_encoded, excluded_pairs_encoded)
+    return freud.locality.NeighborList.from_arrays(
+        num_query_points=nlist.num_query_points,
+        num_points=nlist.num_points,
+        query_point_indices=i_local[keep],
+        point_indices=j_local[keep],
+        vectors=nlist.vectors[keep],
+    )
diff --git a/cmeutils/tests/assets/AB-traj.gsd b/cmeutils/tests/assets/AB-traj.gsd
new file mode 100644
index 0000000..58755f5
Binary files /dev/null and b/cmeutils/tests/assets/AB-traj.gsd differ
diff --git a/cmeutils/tests/assets/lj-fluid.gsd b/cmeutils/tests/assets/lj-fluid.gsd
new file mode 100644
index 0000000..8cde5c6
Binary files /dev/null and b/cmeutils/tests/assets/lj-fluid.gsd differ
diff --git a/cmeutils/tests/base_test.py b/cmeutils/tests/base_test.py
index 3d2ea86..0840e18 100644
--- a/cmeutils/tests/base_test.py
+++ b/cmeutils/tests/base_test.py
@@ -50,6 +50,14 @@ def snap_bond(self, gsdfile_bond):
             snap = f[-1]
         return snap
 
+    @pytest.fixture
+    def AB_chain_gsd(self):
+        return path.join(asset_dir, "AB-traj.gsd")
+
+    @pytest.fixture
+    def LJ_gsd(self):
+        return path.join(asset_dir, "lj-fluid.gsd")
+
     @pytest.fixture
     def butane_gsd(self):
         return path.join(asset_dir, "butanes.gsd")
diff --git a/cmeutils/tests/test_geometry.py b/cmeutils/tests/test_geometry.py
index 8335bb5..ad4f706 100644
--- a/cmeutils/tests/test_geometry.py
+++ b/cmeutils/tests/test_geometry.py
@@ -9,6 +9,7 @@
     angle_between_vectors,
     get_backbone_vector,
     get_plane_normal,
+    get_quaternions,
     moit,
     radial_grid_positions,
     spherical_grid_positions,
@@ -38,6 +39,17 @@ def test_moit(self):
         _moit = moit(points=[(-1, 0, 0), (1, 0, 0)], masses=[1, 1])
         assert np.array_equal(_moit, np.array([0, 2.0, 2.0]))
 
+    def test_get_quaternions(self):
+        with pytest.raises(ValueError):
+            get_quaternions(0)
+
+        with pytest.raises(ValueError):
+            get_quaternions(5.3)
+
+        qs = get_quaternions()
+        assert len(qs) == 20
+        assert len(qs[0]) == 4
+
     def test_get_plane_normal(self):
         points = np.array([[1, 0, 0], [0, 1, 0], [-1, 0, 0], [0, -1, 0]])
         ctr, norm = get_plane_normal(points)
diff --git a/cmeutils/tests/test_gsd.py b/cmeutils/tests/test_gsd.py
index 6c926cb..0272ec4 100644
--- a/cmeutils/tests/test_gsd.py
+++ b/cmeutils/tests/test_gsd.py
@@ -12,6 +12,7 @@
     ellipsoid_gsd,
     frame_to_freud_system,
     get_all_types,
+    get_centers,
     get_molecule_cluster,
     get_type_position,
     identify_snapshot_connections,
@@ -32,6 +33,11 @@
 
 
 class TestGSD(BaseTest):
+    def test_get_centers(self, gsdfile):
+        new_gsdfile = "centers.gsd"
+        centers = get_centers(gsdfile, new_gsdfile)
+        assert isinstance(centers, type(None))
+
     def test_frame_to_freud_system(self, butane_gsd):
         with gsd.hoomd.open(butane_gsd) as traj:
             frame = traj[0]
diff --git a/cmeutils/tests/test_structure.py b/cmeutils/tests/test_structure.py
index ff6779d..9c11abb 100644
--- a/cmeutils/tests/test_structure.py
+++ b/cmeutils/tests/test_structure.py
@@ -3,16 +3,15 @@
 import freud
 import numpy as np
 import pytest
+from scipy.integrate import trapezoid
 
+from cmeutils.geometry import get_quaternions
 from cmeutils.structure import (
-    all_atom_rdf,
     angle_distribution,
     bond_distribution,
     concentration_profile,
     diffraction_pattern,
     dihedral_distribution,
-    get_centers,
-    get_quaternions,
     gsd_rdf,
     order_parameter,
     structure_factor,
@@ -296,37 +295,104 @@ def test_structure_factor_bad_method(self, gsdfile_bond):
         with pytest.raises(ValueError):
             structure_factor(gsdfile_bond, k_min=0.2, k_max=5, method="a")
 
-    def test_gsd_rdf(self, gsdfile_bond):
-        rdf_ex, norm = gsd_rdf(gsdfile_bond, "A", "B")
-        rdf_noex, norm2 = gsd_rdf(gsdfile_bond, "A", "B", exclude_bonded=False)
-        assert np.isclose(norm2, 2 / 3, 1e-4)
-        assert not np.array_equal(rdf_noex, rdf_ex)
-
-    def test_gsd_rdf_samename(self, gsdfile_bond):
-        rdf_ex, norm = gsd_rdf(gsdfile_bond, "A", "A")
-        rdf_noex, norm2 = gsd_rdf(gsdfile_bond, "A", "A", exclude_bonded=False)
-        assert norm2 == 1
-        assert not np.array_equal(rdf_noex, rdf_ex)
+    def test_rdf_bad_args(self, AB_chain_gsd):
+        with pytest.raises(ValueError):
+            gsd_rdf(
+                gsdfile=AB_chain_gsd,
+                start=0,
+                stop=10,
+                exclude_bond_depth=2,
+                exclude_all_bonded=True,
+            )
 
-    def test_gsd_rdf_pair_order(self, gsdfile_bond):
-        rdf, norm = gsd_rdf(gsdfile_bond, "A", "B")
-        rdf_y = rdf.rdf * norm
-        rdf2, norm2 = gsd_rdf(gsdfile_bond, "B", "A")
-        rdf_y2 = rdf2.rdf * norm2
+        with pytest.raises(ValueError):
+            gsd_rdf(
+                gsdfile=AB_chain_gsd,
+                A_name="A",
+                start=0,
+                stop=10,
+            )
 
-        for i, j in zip(rdf_y, rdf_y2):
-            assert np.allclose(i, j, atol=1e-4)
+    def test_gsd_rdf(self, AB_chain_gsd):
+        rdf, scale_factor = gsd_rdf(
+            gsdfile=AB_chain_gsd,
+            start=0,
+            stop=10,
+            exclude_bond_depth=0,
+            exclude_all_bonded=False,
+        )
+        assert isinstance(rdf, freud.density.RDF)
+        rdf.rdf
+        assert scale_factor == 1
 
-    def test_get_quaternions(self):
-        with pytest.raises(ValueError):
-            get_quaternions(0)
+    def test_gsd_rdf_exclude_all_bonded(self, AB_chain_gsd):
+        rdf, scale_factor = gsd_rdf(
+            gsdfile=AB_chain_gsd,
+            start=0,
+            stop=10,
+            exclude_all_bonded=True,
+        )
+        assert scale_factor == 1
+        assert np.array_equal(rdf.rdf, np.zeros_like(rdf.rdf))
 
-        with pytest.raises(ValueError):
-            get_quaternions(5.3)
+    def test_gsd_rdf_exclusions(self, AB_chain_gsd):
+        """Exclude bonded neighbor."""
+        rdf, scale_factor = gsd_rdf(
+            gsdfile=AB_chain_gsd,
+            start=0,
+            stop=10,
+            exclude_bond_depth=1,
+            exclude_all_bonded=False,
+        )
+        assert scale_factor != 1
+        # In this GSD file, the bond length is ~1, so the first peak shows up around r=1
+        zero_indices = np.where(rdf.bin_centers < 1.5)[0]
+        zero_array = np.zeros_like(zero_indices)
+        rdf_values = rdf.rdf[zero_indices]
+        assert np.allclose(rdf_values, zero_array)
+
+        rdf, scale_factor = gsd_rdf(
+            gsdfile=AB_chain_gsd,
+            start=0,
+            stop=10,
+            exclude_bond_depth=2,
+            exclude_all_bonded=False,
+        )
+        assert scale_factor != 1
+        # The second peak shows up around r=2, should be gone
+        zero_indices = np.where(rdf.bin_centers < 2.5)[0]
+        zero_array = np.zeros_like(zero_indices)
+        rdf_values = rdf.rdf[zero_indices]
+        assert np.allclose(rdf_values, zero_array)
+
+    def test_gsd_rdf_r_max(self, LJ_gsd):
+        """Test 2 RDFs with different r_cuts. The values of the shared r_cut region should be very close"""
+        rdf, scale_factor = gsd_rdf(
+            gsdfile=LJ_gsd,
+            start=0,
+            stop=10,
+            r_max=3,
+            exclude_bond_depth=0,
+            exclude_all_bonded=False,
+        )
+        rdf2, scale_factor2 = gsd_rdf(
+            gsdfile=LJ_gsd,
+            start=0,
+            stop=10,
+            r_max=4,
+            exclude_bond_depth=0,
+            exclude_all_bonded=False,
+        )
 
-        qs = get_quaternions()
-        assert len(qs) == 20
-        assert len(qs[0]) == 4
+        check_indices = np.where(rdf.bin_centers < 3)[0]
+        check_indices2 = np.where(rdf2.bin_centers < 3)[0]
+        rdf_integral = trapezoid(
+            rdf.rdf[check_indices], rdf.bin_centers[check_indices]
+        )
+        rdf2_integral = trapezoid(
+            rdf2.rdf[check_indices2], rdf2.bin_centers[check_indices2]
+        )
+        assert np.isclose(rdf_integral, rdf2_integral, rtol=0.05)
 
     def test_order_parameter(self, p3ht_gsd, p3ht_cg_gsd, mapping):
         r_max = 2
@@ -339,15 +405,6 @@ def test_order_parameter(self, p3ht_gsd, p3ht_cg_gsd, mapping):
         assert np.isclose(order[0], 0.33125)
         assert len(cl_idx[0]) == 160
 
-    def test_all_atom_rdf(self, gsdfile):
-        rdf = all_atom_rdf(gsdfile)
-        assert isinstance(rdf, freud.density.RDF)
-
-    def test_get_centers(self, gsdfile):
-        new_gsdfile = "centers.gsd"
-        centers = get_centers(gsdfile, new_gsdfile)
-        assert isinstance(centers, type(None))
-
     def test_conc_profiel(self, slab_snapshot):
         A_indices = np.arange(20)
         B_indices = np.arange(20, 40)
diff --git a/environment-dev.yml b/environment-dev.yml
index 73f9889..fea946f 100644
--- a/environment-dev.yml
+++ b/environment-dev.yml
@@ -8,7 +8,7 @@ dependencies:
   - gsd >=3.0
   - hoomd >=5.0
   - mbuild >=1.0
-  - numpy >=2.0,<2.3
+  - numpy >=2.0
   - pip
   - python >= 3.10,<=3.13
   - matplotlib
diff --git a/environment.yml b/environment.yml
index ca08476..28f7bee 100644
--- a/environment.yml
+++ b/environment.yml
@@ -8,7 +8,7 @@ dependencies:
   - gsd >=3.0
   - hoomd >=5.0
   - mbuild >=1.0
-  - numpy >=2.0,<2.3
+  - numpy >=2.0
   - pip
   - python >= 3.10,<=3.13
   - matplotlib