diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst
index 1d1a1c0ee..61c5dd18c 100644
--- a/docs/source/pythonapi/deplete.rst
+++ b/docs/source/pythonapi/deplete.rst
@@ -59,6 +59,7 @@ data, such as number densities and reaction rates for each material.
    OperatorResult
    ReactionRates
    Results
+   ResultsList
    TransportOperator
 
 Each of the integrator functions also relies on a number of "helper" functions
@@ -71,4 +72,3 @@ as follows:
 
    integrator.CRAM16
    integrator.CRAM48
-   integrator.save_results
diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py
index 33b6d9af1..e49c4e69c 100644
--- a/openmc/deplete/__init__.py
+++ b/openmc/deplete/__init__.py
@@ -20,5 +20,5 @@ from .operator import *
 from .reaction_rates import *
 from .abc import *
 from .results import *
+from .results_list import *
 from .integrator import *
-from .utilities import *
diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py
index 607650dc6..cf8caffdf 100644
--- a/openmc/deplete/integrator/__init__.py
+++ b/openmc/deplete/integrator/__init__.py
@@ -8,4 +8,3 @@ The integrator subcomponents.
 from .cecm import *
 from .cram import *
 from .predictor import *
-from .save_results import *
diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py
index 9a07cd198..61b58d0b9 100644
--- a/openmc/deplete/integrator/cecm.py
+++ b/openmc/deplete/integrator/cecm.py
@@ -4,7 +4,7 @@ import copy
 from collections.abc import Iterable
 
 from .cram import deplete
-from .save_results import save_results
+from ..results import Results
 
 
 def cecm(operator, timesteps, power, print_out=True):
@@ -51,20 +51,20 @@ def cecm(operator, timesteps, power, print_out=True):
         for i, (dt, p) in enumerate(zip(timesteps, power)):
             # Get beginning-of-timestep reaction rates
             x = [copy.deepcopy(vec)]
-            results = [operator(x[0], p)]
+            op_results = [operator(x[0], p)]
 
             # Deplete for first half of timestep
-            x_middle = deplete(chain, x[0], results[0], dt/2, print_out)
+            x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out)
 
             # Get middle-of-timestep reaction rates
             x.append(x_middle)
-            results.append(operator(x_middle, p))
+            op_results.append(operator(x_middle, p))
 
             # Deplete for full timestep using beginning-of-step materials
-            x_end = deplete(chain, x[0], results[1], dt, print_out)
+            x_end = deplete(chain, x[0], op_results[1], dt, print_out)
 
             # Create results, write to disk
-            save_results(operator, x, results, [t, t + dt], i)
+            Results.save(operator, x, op_results, [t, t + dt], i)
 
             # Advance time, update vector
             t += dt
@@ -72,7 +72,7 @@ def cecm(operator, timesteps, power, print_out=True):
 
         # Perform one last simulation
         x = [copy.deepcopy(vec)]
-        results = [operator(x[0], power[-1])]
+        op_results = [operator(x[0], power[-1])]
 
         # Create results, write to disk
-        save_results(operator, x, results, [t, t], len(timesteps))
+        Results.save(operator, x, op_results, [t, t], len(timesteps))
diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py
index 444ca7baa..9be992c16 100644
--- a/openmc/deplete/integrator/predictor.py
+++ b/openmc/deplete/integrator/predictor.py
@@ -4,7 +4,7 @@ import copy
 from collections.abc import Iterable
 
 from .cram import deplete
-from .save_results import save_results
+from ..results import Results
 
 
 def predictor(operator, timesteps, power, print_out=True):
@@ -46,13 +46,13 @@ def predictor(operator, timesteps, power, print_out=True):
         for i, (dt, p) in enumerate(zip(timesteps, power)):
             # Get beginning-of-timestep reaction rates
             x = [copy.deepcopy(vec)]
-            results = [operator(x[0], p)]
+            op_results = [operator(x[0], p)]
 
             # Create results, write to disk
-            save_results(operator, x, results, [t, t + dt], i)
+            Results.save(operator, x, op_results, [t, t + dt], i)
 
             # Deplete for full timestep
-            x_end = deplete(chain, x[0], results[0], dt, print_out)
+            x_end = deplete(chain, x[0], op_results[0], dt, print_out)
 
             # Advance time, update vector
             t += dt
@@ -60,7 +60,7 @@ def predictor(operator, timesteps, power, print_out=True):
 
         # Perform one last simulation
         x = [copy.deepcopy(vec)]
-        results = [operator(x[0], power[-1])]
+        op_results = [operator(x[0], power[-1])]
 
         # Create results, write to disk
-        save_results(operator, x, results, [t, t], len(timesteps))
+        Results.save(operator, x, op_results, [t, t], len(timesteps))
diff --git a/openmc/deplete/integrator/save_results.py b/openmc/deplete/integrator/save_results.py
deleted file mode 100644
index 98245fdea..000000000
--- a/openmc/deplete/integrator/save_results.py
+++ /dev/null
@@ -1,42 +0,0 @@
-""" Generic result saving code for integrators.
-
-"""
-from ..results import Results, write_results
-
-
-def save_results(op, x, op_results, t, step_ind):
-    """Creates and writes depletion results to disk
-
-    Parameters
-    ----------
-    op : openmc.deplete.TransportOperator
-        The operator used to generate these results.
-    x : list of list of numpy.array
-        The prior x vectors.  Indexed [i][cell] using the above equation.
-    op_results : list of openmc.deplete.OperatorResult
-        Results of applying transport operator
-    t : list of float
-        Time indices.
-    step_ind : int
-        Step index.
-
-    """
-    # Get indexing terms
-    vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
-
-    # Create results
-    stages = len(x)
-    results = Results()
-    results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
-
-    n_mat = len(burn_list)
-
-    for i in range(stages):
-        for mat_i in range(n_mat):
-            results[i, mat_i, :] = x[i][mat_i][:]
-
-    results.k = [r.k for r in op_results]
-    results.rates = [r.rates for r in op_results]
-    results.time = t
-
-    write_results(results, "depletion_results.h5", step_ind)
diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py
index b7e5623c7..ab740e61e 100644
--- a/openmc/deplete/results.py
+++ b/openmc/deplete/results.py
@@ -11,7 +11,6 @@ import h5py
 
 from . import comm, have_mpi
 from .reaction_rates import ReactionRates
-from openmc.checkvalue import check_filetype_version
 
 _VERSION_RESULTS = (1, 0)
 
@@ -146,6 +145,27 @@ class Results(object):
         # Create storage array
         self.data = np.zeros((stages, self.n_mat, self.n_nuc))
 
+    def export_to_hdf5(self, filename, step):
+        """Export results to an HDF5 file
+
+        Parameters
+        ----------
+        filename : str
+            The filename to write to
+        step : int
+            What step is this?
+
+        """
+        if have_mpi and h5py.get_config().mpi:
+            kwargs = {'driver': 'mpio', 'comm': comm}
+        else:
+            kwargs = {}
+
+        kwargs['mode'] = "w" if step == 0 else "a"
+
+        with h5py.File(filename, **kwargs) as handle:
+            self._to_hdf5(handle, step)
+
     def _write_hdf5_metadata(self, handle):
         """Writes result metadata in HDF5 file
 
@@ -217,7 +237,7 @@ class Results(object):
 
         handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
 
-    def to_hdf5(self, handle, index):
+    def _to_hdf5(self, handle, index):
         """Converts results object into an hdf5 object.
 
         Parameters
@@ -338,49 +358,40 @@ class Results(object):
 
         return results
 
+    @staticmethod
+    def save(op, x, op_results, t, step_ind):
+        """Creates and writes depletion results to disk
 
-def write_results(result, filename, step):
-    """Outputs result to an HDF5 file.
+        Parameters
+        ----------
+        op : openmc.deplete.TransportOperator
+            The operator used to generate these results.
+        x : list of list of numpy.array
+            The prior x vectors.  Indexed [i][cell] using the above equation.
+        op_results : list of openmc.deplete.OperatorResult
+            Results of applying transport operator
+        t : list of float
+            Time indices.
+        step_ind : int
+            Step index.
 
-    Parameters
-    ----------
-    result : Results
-        Object to be stored in a file.
-    filename : String
-        Target filename.
-    step : int
-        What step is this?
+        """
+        # Get indexing terms
+        vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
 
-    """
-    if have_mpi and h5py.get_config().mpi:
-        kwargs = {'driver': 'mpio', 'comm': comm}
-    else:
-        kwargs = {}
+        # Create results
+        stages = len(x)
+        results = Results()
+        results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
 
-    kwargs['mode'] = "w" if step == 0 else "a"
+        n_mat = len(burn_list)
 
-    with h5py.File(filename, **kwargs) as handle:
-        result.to_hdf5(handle, step)
+        for i in range(stages):
+            for mat_i in range(n_mat):
+                results[i, mat_i, :] = x[i][mat_i][:]
 
+        results.k = [r.k for r in op_results]
+        results.rates = [r.rates for r in op_results]
+        results.time = t
 
-def read_results(filename):
-    """Return a list of Results objects from an HDF5 file.
-
-    Parameters
-    ----------
-    filename : str
-        The filename to read from.
-
-    Returns
-    -------
-    results : list of openmc.deplete.Results
-        The result objects.
-
-    """
-    with h5py.File(str(filename), "r") as fh:
-        check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
-
-        # Get number of results stored
-        n = fh["number"].value.shape[0]
-
-        return [Results.from_hdf5(fh, i) for i in range(n)]
+        results.export_to_hdf5("depletion_results.h5", step_ind)
diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py
new file mode 100644
index 000000000..d3d45e955
--- /dev/null
+++ b/openmc/deplete/results_list.py
@@ -0,0 +1,105 @@
+import h5py
+import numpy as np
+
+from .results import Results, _VERSION_RESULTS
+from openmc.checkvalue import check_filetype_version
+
+
+class ResultsList(list):
+    """A list of openmc.deplete.Results objects
+
+    Parameters
+    ----------
+    filename : str
+        The filename to read from.
+
+    """
+    def __init__(self, filename):
+        super().__init__()
+        with h5py.File(str(filename), "r") as fh:
+            check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
+
+            # Get number of results stored
+            n = fh["number"].value.shape[0]
+
+            for i in range(n):
+                self.append(Results.from_hdf5(fh, i))
+
+    def get_atoms(self, mat, nuc):
+        """Get nuclide concentration over time from a single material
+
+        Parameters
+        ----------
+        mat : str
+            Material name to evaluate
+        nuc : str
+            Nuclide name to evaluate
+
+        Returns
+        -------
+        time : numpy.ndarray
+            Array of times in [s]
+        concentration : numpy.ndarray
+            Total number of atoms for specified nuclide
+
+        """
+        time = np.empty_like(self)
+        concentration = np.empty_like(self)
+
+        # Evaluate value in each region
+        for i, result in enumerate(self):
+            time[i] = result.time[0]
+            concentration[i] = result[0, mat, nuc]
+
+        return time, concentration
+
+    def get_reaction_rate(self, mat, nuc, rx):
+        """Get reaction rate in a single material/nuclide over time
+
+        Parameters
+        ----------
+        mat : str
+            Material name to evaluate
+        nuc : str
+            Nuclide name to evaluate
+        rx : str
+            Reaction rate to evaluate
+
+        Returns
+        -------
+        time : numpy.ndarray
+            Array of times in [s]
+        rate : numpy.ndarray
+            Array of reaction rates
+
+        """
+        time = np.empty_like(self)
+        rate = np.empty_like(self)
+
+        # Evaluate value in each region
+        for i, result in enumerate(self):
+            time[i] = result.time[0]
+            rate[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
+
+        return time, rate
+
+    def get_eigenvalue(self):
+        """Evaluates the eigenvalue from a results list.
+
+        Returns
+        -------
+        time : numpy.ndarray
+            Array of times in [s]
+        eigenvalue : numpy.ndarray
+            k-eigenvalue at each time
+
+        """
+        time = np.empty_like(self)
+        eigenvalue = np.empty_like(self)
+
+        # Get time/eigenvalue at each point
+        for i, result in enumerate(self):
+            time[i] = result.time[0]
+            eigenvalue[i] = result.k[0]
+
+        return time, eigenvalue
diff --git a/openmc/deplete/utilities.py b/openmc/deplete/utilities.py
deleted file mode 100644
index 59d530546..000000000
--- a/openmc/deplete/utilities.py
+++ /dev/null
@@ -1,101 +0,0 @@
-"""The utilities module.
-
-Contains functions that can be used to post-process objects that come out of
-the results module.
-"""
-
-import numpy as np
-
-
-def evaluate_single_nuclide(results, mat, nuc):
-    """Evaluates a single nuclide in a single material from a results list.
-
-    Parameters
-    ----------
-    results : list of results
-        The results to extract data from.  Must be sorted and continuous.
-    mat : str
-        Material name to evaluate
-    nuc : str
-        Nuclide name to evaluate
-
-    Returns
-    -------
-    time : numpy.ndarray
-        Time vector
-    concentration : numpy.ndarray
-        Total number of atoms in the material
-
-    """
-    n_points = len(results)
-    time = np.zeros(n_points)
-    concentration = np.zeros(n_points)
-
-    # Evaluate value in each region
-    for i, result in enumerate(results):
-        time[i] = result.time[0]
-        concentration[i] = result[0, mat, nuc]
-
-    return time, concentration
-
-
-def evaluate_reaction_rate(results, mat, nuc, rx):
-    """Return reaction rate in a single material/nuclide from a results list.
-
-    Parameters
-    ----------
-    results : list of openmc.deplete.Results
-        The results to extract data from.  Must be sorted and continuous.
-    mat : str
-        Material name to evaluate
-    nuc : str
-        Nuclide name to evaluate
-    rx : str
-        Reaction rate to evaluate
-
-    Returns
-    -------
-    time : numpy.ndarray
-        Time vector.
-    rate : numpy.ndarray
-        Reaction rate.
-
-    """
-    n_points = len(results)
-    time = np.zeros(n_points)
-    rate = np.zeros(n_points)
-    # Evaluate value in each region
-    for i, result in enumerate(results):
-        time[i] = result.time[0]
-        rate[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
-
-    return time, rate
-
-
-def evaluate_eigenvalue(results):
-    """Evaluates the eigenvalue from a results list.
-
-    Parameters
-    ----------
-    results : list of openmc.deplete.Results
-        The results to extract data from.  Must be sorted and continuous.
-
-    Returns
-    -------
-    time : numpy.ndarray
-        Time vector.
-    eigenvalue : numpy.ndarray
-        Eigenvalue.
-
-    """
-    n_points = len(results)
-    time = np.zeros(n_points)
-    eigenvalue = np.zeros(n_points)
-
-    # Evaluate value in each region
-    for i, result in enumerate(results):
-
-        time[i] = result.time[0]
-        eigenvalue[i] = result.k[0]
-
-    return time, eigenvalue
diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py
index 706793d93..05fe97a95 100644
--- a/tests/dummy_operator.py
+++ b/tests/dummy_operator.py
@@ -34,19 +34,15 @@ class DummyOperator(TransportOperator):
 
         Returns
         -------
-        k : float
-            Zero.
-        rates : ReactionRates
-            Reaction rates from this simulation.
-        seed : int
-            Zero.
+        openmc.deplete.OperatorResult
+            Result of transport operator
+
         """
+        mats = ["1"]
+        nuclides = ["1", "2"]
+        reactions = ["1"]
 
-        cell_to_ind = {"1" : 0}
-        nuc_to_ind = {"1" : 0, "2" : 1}
-        react_to_ind = {"1" : 0}
-
-        reaction_rates = ReactionRates(cell_to_ind, nuc_to_ind, react_to_ind)
+        reaction_rates = ReactionRates(mats, nuclides, reactions)
 
         reaction_rates[0, 0, 0] = vec[0][0]
         reaction_rates[0, 1, 0] = vec[0][1]
@@ -105,18 +101,18 @@ class DummyOperator(TransportOperator):
         return ["1", "2"]
 
     @property
-    def burn_list(self):
+    def local_mats(self):
         """
-        burn_list : list of str
-            A list of all cell IDs to be burned.  Used for sorting the simulation.
+        local_mats : list of str
+            A list of all material IDs to be burned.  Used for sorting the simulation.
         """
 
         return ["1"]
 
     @property
-    def mat_tally_ind(self):
+    def burnable_mats(self):
         """Maps cell name to index in global geometry."""
-        return {"1": 0}
+        return self.local_mats
 
 
     @property
@@ -125,11 +121,11 @@ class DummyOperator(TransportOperator):
         reaction_rates : ReactionRates
             Reaction rates from the last operator step.
         """
-        cell_to_ind = {"1" : 0}
-        nuc_to_ind = {"1" : 0, "2" : 1}
-        react_to_ind = {"1" : 0}
+        mats = ["1"]
+        nuclides = ["1", "2"]
+        reactions = ["1"]
 
-        return ReactionRates(cell_to_ind, nuc_to_ind, react_to_ind)
+        return ReactionRates(mats, nuclides, reactions)
 
     def initial_condition(self):
         """Returns initial vector.
@@ -153,8 +149,8 @@ class DummyOperator(TransportOperator):
             A list of all nuclide names. Used for sorting the simulation.
         burn_list : list of int
             A list of all cell IDs to be burned.  Used for sorting the simulation.
-        full_burn_dict : OrderedDict of str to int
+        full_burn_list : OrderedDict of str to int
             Maps cell name to index in global geometry.
-        """
 
-        return self.volume, self.nuc_list, self.burn_list, self.mat_tally_ind
+        """
+        return self.volume, self.nuc_list, self.local_mats, self.burnable_mats
diff --git a/tests/regression_tests/test_deplete_full.py b/tests/regression_tests/test_deplete_full.py
index fb31ac0c0..19e2c7664 100644
--- a/tests/regression_tests/test_deplete_full.py
+++ b/tests/regression_tests/test_deplete_full.py
@@ -8,8 +8,6 @@ import numpy as np
 import openmc
 from openmc.data import JOULE_PER_EV
 import openmc.deplete
-from openmc.deplete import results
-from openmc.deplete import utilities
 
 from tests.regression_tests import config
 from .example_geometry import generate_problem
@@ -67,8 +65,8 @@ def test_full(run_in_tmpdir):
         return
 
     # Load the reference/test results
-    res_test = results.read_results(path_test)
-    res_ref = results.read_results(path_reference)
+    res_test = openmc.deplete.ResultsList(path_test)
+    res_ref = openmc.deplete.ResultsList(path_reference)
 
     # Assert same mats
     for mat in res_ref[0].mat_to_ind:
@@ -88,8 +86,8 @@ def test_full(run_in_tmpdir):
     tol = 1.0e-6
     for mat in res_test[0].mat_to_ind:
         for nuc in res_test[0].nuc_to_ind:
-            _, y_test = utilities.evaluate_single_nuclide(res_test, mat, nuc)
-            _, y_old = utilities.evaluate_single_nuclide(res_ref, mat, nuc)
+            _, y_test = res_test.get_atoms(mat, nuc)
+            _, y_old = res_ref.get_atoms(mat, nuc)
 
             # Test each point
             correct = True
diff --git a/tests/unit_tests/test_deplete_cecm.py b/tests/unit_tests/test_deplete_cecm.py
index 6deb6cd3d..466a2eec7 100644
--- a/tests/unit_tests/test_deplete_cecm.py
+++ b/tests/unit_tests/test_deplete_cecm.py
@@ -5,8 +5,6 @@ These tests integrate a simple test problem described in dummy_geometry.py.
 
 from pytest import approx
 import openmc.deplete
-from openmc.deplete import results
-from openmc.deplete import utilities
 
 from tests import dummy_operator
 
@@ -23,10 +21,10 @@ def test_cecm(run_in_tmpdir):
     openmc.deplete.cecm(op, dt, power, print_out=False)
 
     # Load the files
-    res = results.read_results(op.output_dir / "depletion_results.h5")
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
 
-    _, y1 = utilities.evaluate_single_nuclide(res, "1", "1")
-    _, y2 = utilities.evaluate_single_nuclide(res, "1", "2")
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
 
     # Mathematica solution
     s1 = [1.86872629872102, 1.395525772416039]
diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py
index b20990088..a1768625b 100644
--- a/tests/unit_tests/test_deplete_integrator.py
+++ b/tests/unit_tests/test_deplete_integrator.py
@@ -1,4 +1,4 @@
-"""Tests for integrator.py
+"""Tests for saving results
 
 It is worth noting that openmc.deplete.integrate is extremely complex, to the
 point I am unsure if it can be reasonably unit-tested.  For the time being, it
@@ -11,11 +11,11 @@ import os
 from unittest.mock import MagicMock
 
 import numpy as np
-from openmc.deplete import (integrator, ReactionRates, results, comm,
+from openmc.deplete import (ReactionRates, Results, ResultsList, comm,
                             OperatorResult)
 
 
-def test_save_results(run_in_tmpdir):
+def test_results_save(run_in_tmpdir):
     """Test data save module"""
 
     stages = 3
@@ -72,11 +72,11 @@ def test_save_results(run_in_tmpdir):
 
     op_result1 = [OperatorResult(k, rates) for k, rates in zip(eigvl1, rate1)]
     op_result2 = [OperatorResult(k, rates) for k, rates in zip(eigvl2, rate2)]
-    integrator.save_results(op, x1, op_result1, t1, 0)
-    integrator.save_results(op, x2, op_result2, t2, 1)
+    Results.save(op, x1, op_result1, t1, 0)
+    Results.save(op, x2, op_result2, t2, 1)
 
     # Load the files
-    res = results.read_results("depletion_results.h5")
+    res = ResultsList("depletion_results.h5")
 
     for i in range(stages):
         for mat_i, mat in enumerate(burn_list):
diff --git a/tests/unit_tests/test_deplete_predictor.py b/tests/unit_tests/test_deplete_predictor.py
index 0d283855c..50803e508 100644
--- a/tests/unit_tests/test_deplete_predictor.py
+++ b/tests/unit_tests/test_deplete_predictor.py
@@ -5,8 +5,6 @@ These tests integrate a simple test problem described in dummy_geometry.py.
 
 from pytest import approx
 import openmc.deplete
-from openmc.deplete import results
-from openmc.deplete import utilities
 
 from tests import dummy_operator
 
@@ -23,10 +21,10 @@ def test_predictor(run_in_tmpdir):
     openmc.deplete.predictor(op, dt, power, print_out=False)
 
     # Load the files
-    res = results.read_results(op.output_dir / "depletion_results.h5")
+    res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
 
-    _, y1 = utilities.evaluate_single_nuclide(res, "1", "1")
-    _, y2 = utilities.evaluate_single_nuclide(res, "1", "2")
+    _, y1 = res.get_atoms("1", "1")
+    _, y2 = res.get_atoms("1", "2")
 
     # Mathematica solution
     s1 = [2.46847546272295, 0.986431226850467]
diff --git a/tests/regression_tests/test_deplete_utilities.py b/tests/unit_tests/test_deplete_resultslist.py
similarity index 62%
rename from tests/regression_tests/test_deplete_utilities.py
rename to tests/unit_tests/test_deplete_resultslist.py
index 82d4d56a4..aad8cd9f6 100644
--- a/tests/regression_tests/test_deplete_utilities.py
+++ b/tests/unit_tests/test_deplete_resultslist.py
@@ -1,26 +1,22 @@
-""" Tests the utilities classes.
-
-This also tests the results read/write code.
-"""
+"""Tests the ResultsList class"""
 
 from pathlib import Path
 
 import numpy as np
 import pytest
-from openmc.deplete import results
-from openmc.deplete import utilities
+import openmc.deplete
 
 
 @pytest.fixture
 def res():
     """Load the reference results"""
-    filename = Path(__file__).with_name('test_reference.h5')
-    return results.read_results(filename)
+    filename = Path(__file__).parents[1] / 'regression_tests' / 'test_reference.h5'
+    return openmc.deplete.ResultsList(filename)
 
 
-def test_evaluate_single_nuclide(res):
-    """Tests evaluating single nuclide utility code."""
-    t, n = utilities.evaluate_single_nuclide(res, "1", "Xe135")
+def test_get_atoms(res):
+    """Tests evaluating single nuclide concentration."""
+    t, n = res.get_atoms("1", "Xe135")
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
     n_ref = [6.6747328233649218e+08, 3.5421791038348462e+14,
@@ -29,9 +25,9 @@ def test_evaluate_single_nuclide(res):
     np.testing.assert_array_equal(t, t_ref)
     np.testing.assert_array_equal(n, n_ref)
 
-def test_evaluate_reaction_rate(res):
-    """Tests evaluating reaction rate utility code."""
-    t, r = utilities.evaluate_reaction_rate(res, "1", "Xe135", "(n,gamma)")
+def test_get_reaction_rate(res):
+    """Tests evaluating reaction rate."""
+    t, r = res.get_reaction_rate("1", "Xe135", "(n,gamma)")
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
     n_ref = np.array([6.6747328233649218e+08, 3.5421791038348462e+14,
@@ -43,9 +39,9 @@ def test_evaluate_reaction_rate(res):
     np.testing.assert_array_equal(r, n_ref * xs_ref)
 
 
-def test_evaluate_eigenvalue(res):
+def test_get_eigenvalue(res):
     """Tests evaluating eigenvalue."""
-    t, k = utilities.evaluate_eigenvalue(res)
+    t, k = res.get_eigenvalue()
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
     k_ref = [1.181281798790367, 1.1798750921988739, 1.1965943696058159,