diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 63a4ad32..50e9aeb3 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -29,7 +29,7 @@ repos:
         name: isort (python)
         args:
           [ --profile=black, --line-length=80 ]
-        exclude: 'flowermd/tests/assets/.* '
+        exclude: 'flowermd/tests/assets/.* | */__init__.py'
 
   - repo: https://github.com/pycqa/flake8
     rev: 7.0.0
diff --git a/flowermd/__init__.py b/flowermd/__init__.py
index e0995a2c..5af4d2ba 100644
--- a/flowermd/__init__.py
+++ b/flowermd/__init__.py
@@ -1,5 +1,6 @@
 """flowerMD package."""
 
+from .internal.units import Units
 from .base import (
     CoPolymer,
     Lattice,
diff --git a/flowermd/base/simulation.py b/flowermd/base/simulation.py
index 04ba2d44..68320387 100644
--- a/flowermd/base/simulation.py
+++ b/flowermd/base/simulation.py
@@ -11,7 +11,8 @@
 import numpy as np
 import unyt as u
 
-from flowermd.internal import validate_ref_value
+from flowermd import Units
+from flowermd.internal import validate_unit
 from flowermd.utils.actions import StdOutLogger, UpdateWalls
 from flowermd.utils.base_types import HOOMDThermostats
 
@@ -102,6 +103,7 @@ def __init__(
         ]
         self.integrator = None
         self._dt = dt
+        self._kT = None
         self._reference_values = dict()
         self._reference_values = reference_values
         if rigid_constraint and not isinstance(
@@ -198,15 +200,13 @@ def reference_length(self, length):
 
         Parameters
         ----------
-        length : string or unyt.unyt_quantity, required
+        length : reference length * `flowermd.Units`, required
             The reference length of the system.
-            It can be provided in the following forms:
-            1) A string with the format of "value unit", for example "1 nm".
-            2) A unyt.unyt_quantity object with the correct dimension. For
-            example, unyt.unyt_quantity(1, "nm").
+            It can be provided in the following form of:
+            value * `flowermd.Units`, for example 1 * `flowermd.Units.angstrom`.
 
         """
-        validated_length = validate_ref_value(length, u.dimensions.length)
+        validated_length = validate_unit(length, u.dimensions.length)
         self._reference_values["length"] = validated_length
 
     @reference_energy.setter
@@ -215,15 +215,13 @@ def reference_energy(self, energy):
 
         Parameters
         ----------
-        energy : string or unyt.unyt_quantity, required
+        energy : reference energy * `flowermd.Units`, required
             The reference energy of the system.
-            It can be provided in the following forms:
-            1) A string with the format of "value unit", for example "1 kJ/mol".
-            2) A unyt.unyt_quantity object with the correct dimension. For
-            example, unyt.unyt_quantity(1, "kJ/mol").
+            It can be provided in the following form of:
+            value * `flowermd.Units`, for example 1 * `flowermd.Units.kcal/mol`.
 
         """
-        validated_energy = validate_ref_value(energy, u.dimensions.energy)
+        validated_energy = validate_unit(energy, u.dimensions.energy)
         self._reference_values["energy"] = validated_energy
 
     @reference_mass.setter
@@ -232,15 +230,12 @@ def reference_mass(self, mass):
 
         Parameters
         ----------
-        mass : string or unyt.unyt_quantity, required
+        mass : reference mass * `flowermd.Units`, required
             The reference mass of the system.
-            It can be provided in the following forms:
-            1) A string with the format of "value unit", for example "1 amu".
-            2) A unyt.unyt_quantity object with the correct dimension. For
-            example, unyt.unyt_quantity(1, "amu").
-
+            It can be provided in the following form of:
+            value * `flowermd.Units`, for example 1 * `flowermd.Units.amu`.
         """
-        validated_mass = validate_ref_value(mass, u.dimensions.mass)
+        validated_mass = validate_unit(mass, u.dimensions.mass)
         self._reference_values["mass"] = validated_mass
 
     @reference_values.setter
@@ -376,19 +371,79 @@ def real_timestep(self):
         if self._reference_values.get("mass"):
             mass = self._reference_values["mass"].to("kg")
         else:
-            mass = 1 * u.kg
+            mass = 1 * Units.kg
         if self._reference_values.get("length"):
             dist = self.reference_length.to("m")
         else:
-            dist = 1 * u.m
+            dist = 1 * Units.m
         if self._reference_values.get("energy"):
             energy = self.reference_energy.to("J")
         else:
-            energy = 1 * u.J
+            energy = 1 * Units.J
         tau = (mass * (dist**2)) / energy
         timestep = self.dt * (tau**0.5)
         return timestep
 
+    @property
+    def real_temperature(self):
+        """The temperature of the simulation in Kelvin."""
+        if not self._kT:
+            raise ValueError(
+                "Temperature is not set. Please specify the temperature when "
+                "running the simulation, using one of the following run"
+                " methods: `run_nvt`, `run_npt`, `run_update_volume`."
+            )
+        if self._reference_values.get("energy"):
+            energy = self.reference_energy.to("J")
+        else:
+            energy = 1 * Units.J
+        temperature = (self._kT * energy) / u.boltzmann_constant_mks
+        return temperature
+
+    def _temperature_to_kT(self, temperature):
+        """Convert temperature to kT."""
+        if self._reference_values.get("energy"):
+            energy = self.reference_energy.to("J")
+        else:
+            energy = 1 * Units.J
+        if isinstance(temperature, u.unyt_array) or isinstance(
+            temperature, u.unyt_quantity
+        ):
+            temperature = temperature.to("K")
+        else:
+            temperature = temperature * Units.K
+        kT = (temperature * u.boltzmann_constant_mks) / energy
+        return float(kT)
+
+    def _setup_temperature(self, temperature):
+        if isinstance(temperature, (float, int)):
+            # assuming temperature is kT
+            return temperature
+        else:
+            return self._temperature_to_kT(
+                validate_unit(temperature, u.dimensions.temperature)
+            )
+
+    def _time_length_to_n_steps(self, time_length):
+        """Convert time length to number of steps."""
+        if isinstance(time_length, u.unyt_array) or isinstance(
+            time_length, u.unyt_quantity
+        ):
+            time_length = time_length.to("s")
+        else:
+            time_length = time_length * Units.s
+        real_timestep = self.real_timestep.to("s")
+        return int(time_length / real_timestep)
+
+    def _setup_n_steps(self, duration):
+        if isinstance(duration, int):
+            # assuming duration is num steps
+            return duration
+        else:
+            return self._time_length_to_n_steps(
+                validate_unit(duration, u.dimensions.time)
+            )
+
     @property
     def integrate_group(self):
         """The group of particles to apply the integrator to.
@@ -623,10 +678,10 @@ def remove_walls(self, wall_axis):
     def run_update_volume(
         self,
         final_box_lengths,
-        n_steps,
-        period,
-        kT,
+        temperature,
         tau_kt,
+        duration,
+        period,
         thermalize_particles=True,
         write_at_start=True,
     ):
@@ -643,19 +698,23 @@ def run_update_volume(
         ----------
         final_box_lengths : np.ndarray or unyt.array.unyt_array, shape=(3,), required # noqa: E501
             The final box edge lengths in (x, y, z) order.
-        n_steps : int, required
-            Number of steps to run during volume update.
-        period : int, required
-            The number of steps ran between each box update iteration.
-        kT : float or hoomd.variant.Ramp, required
-            The temperature to use during volume update.
+        temperature : flowermd.utils.units or float or int, required
+            The temperature to use during volume update. If no unit is provided,
+            the temperature is assumed to be kT (temperature times Boltzmann
+            constant).
         tau_kt : float, required
             Thermostat coupling period (in simulation time units).
+        duration : int or flowermd.utils.units, required
+            The number of steps or time length to run the simulation. If no unit
+            is provided, the time is assumed to be the number of steps.
+        period : int, required
+            The number of steps ran between each box update iteration.
         write_at_start : bool, default True
             When set to True, triggers writers that evaluate to True
             for the initial step to execute before the next simulation
             time step.
 
+
         Examples
         --------
         In this example, a low density system is initialized with `Pack`
@@ -671,24 +730,29 @@ def run_update_volume(
             pps_mols = PPS(num_mols=20, lengths=15)
             pps_system = Pack(
                 molecules=[pps_mols],
-                force_field=OPLS_AA_PPS(),
-                r_cut=2.5,
                 density=0.5,
+            )
+            pps_system.apply_forcefield(
+                r_cut=2.5,
+                force_field=OPLS_AA_PPS(),
                 auto_scale=True,
                 scale_charges=True
             )
-            sim = Simulation(
-                initial_state=pps_system.hoomd_snapshot,
-                forcefield=pps_system.hoomd_forcefield
-            )
+            sim = Simulation.from_system(pps_system)
             target_box = flowermd.utils.get_target_box_mass_density(
                 density=1.1 * unyt.g/unyt.cm**3, mass=sim.mass.to("g")
             )
             sim.run_update_volume(
-                n_steps=1e4, kT=1.0, tau_kt=1.0, final_box_lengths=target_box
+                final_box_lengths=target_box,
+                temperature=1.0,
+                tau_kt=1.0,
+                duration=1e4,
+                period=100
             )
 
         """
+        self._kT = self._setup_temperature(temperature)
+        _n_steps = self._setup_n_steps(duration)
         if self.reference_length and hasattr(final_box_lengths, "to"):
             ref_unit = self.reference_length.units
             final_box_lengths = final_box_lengths.to(ref_unit)
@@ -701,7 +765,7 @@ def run_update_volume(
         )
         resize_trigger = hoomd.trigger.Periodic(period)
         box_ramp = hoomd.variant.Ramp(
-            A=0, B=1, t_start=self.timestep, t_ramp=int(n_steps)
+            A=0, B=1, t_start=self.timestep, t_ramp=int(_n_steps)
         )
         initial_box = self.state.box
 
@@ -716,13 +780,13 @@ def run_update_volume(
             integrator_method=hoomd.md.methods.ConstantVolume,
             method_kwargs={
                 "thermostat": self._initialize_thermostat(
-                    {"kT": kT, "tau": tau_kt}
+                    {"kT": self._kT, "tau": tau_kt}
                 ),
                 "filter": self.integrate_group,
             },
         )
         if thermalize_particles:
-            self._thermalize_system(kT)
+            self._thermalize_system(self._kT)
 
         if self._wall_forces:
             wall_update = UpdateWalls(sim=self)
@@ -730,20 +794,20 @@ def run_update_volume(
                 trigger=resize_trigger, action=wall_update
             )
             self.operations.updaters.append(wall_updater)
-        std_out_logger = StdOutLogger(n_steps=n_steps, sim=self)
+        std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self)
         std_out_logger_printer = hoomd.update.CustomUpdater(
             trigger=hoomd.trigger.Periodic(self._std_out_freq),
             action=std_out_logger,
         )
         self.operations.updaters.append(std_out_logger_printer)
-        self.run(steps=n_steps + 1, write_at_start=write_at_start)
+        self.run(steps=_n_steps + 1, write_at_start=write_at_start)
         self.operations.updaters.remove(std_out_logger_printer)
         self.operations.updaters.remove(box_resizer)
 
     def run_langevin(
         self,
-        n_steps,
-        kT,
+        duration,
+        temperature,
         tally_reservoir_energy=False,
         default_gamma=1.0,
         default_gamma_r=(1.0, 1.0, 1.0),
@@ -754,10 +818,13 @@ def run_langevin(
 
         Parameters
         ----------
-        n_steps : int, required
-            Number of steps to run the simulation.
-        kT : int or hoomd.variant.Ramp, required
-            The temperature to use during the simulation.
+        duration : int or flowermd.utils.units, required
+            The number of steps or time length to run the simulation. If no unit
+            is provided, the time is assumed to be the number of steps.
+        temperature : flowermd.utils.units or float or int, required
+            The temperature to use during the simulation. If no unit is
+            provided, the temperature is assumed to be kT (temperature times
+            Boltzmann constant).
         tally_reservoir_energy : bool, default False
             When set to True, energy exchange between the thermal reservoir
              and the particles is tracked.
@@ -773,34 +840,36 @@ def run_langevin(
             time step.
 
         """
+        self._kT = self._setup_temperature(temperature)
+        _n_steps = self._setup_n_steps(duration)
         self.set_integrator_method(
             integrator_method=hoomd.md.methods.Langevin,
             method_kwargs={
                 "filter": self.integrate_group,
-                "kT": kT,
+                "kT": self._kT,
                 "tally_reservoir_energy": tally_reservoir_energy,
                 "default_gamma": default_gamma,
                 "default_gamma_r": default_gamma_r,
             },
         )
         if thermalize_particles:
-            self._thermalize_system(kT)
-        std_out_logger = StdOutLogger(n_steps=n_steps, sim=self)
+            self._thermalize_system(self._kT)
+        std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self)
         std_out_logger_printer = hoomd.update.CustomUpdater(
             trigger=hoomd.trigger.Periodic(self._std_out_freq),
             action=std_out_logger,
         )
         self.operations.updaters.append(std_out_logger_printer)
-        self.run(steps=n_steps, write_at_start=write_at_start)
+        self.run(steps=_n_steps, write_at_start=write_at_start)
         self.operations.updaters.remove(std_out_logger_printer)
 
     def run_NPT(
         self,
-        n_steps,
-        kT,
         pressure,
-        tau_kt,
         tau_pressure,
+        temperature,
+        tau_kt,
+        duration,
         couple="xyz",
         box_dof=[True, True, True, False, False, False],
         rescale_all=False,
@@ -812,16 +881,19 @@ def run_NPT(
 
         Parameters
         ----------
-        n_steps: int, required
-            Number of steps to run the simulation.
-        kT: int or hoomd.variant.Ramp, required
-            The temperature to use during the simulation.
         pressure: int or hoomd.variant.Ramp, required
             The pressure to use during the simulation.
-        tau_kt: float, required
-            Thermostat coupling period (in simulation time units).
         tau_pressure: float, required
             Barostat coupling period.
+        temperature: flowermd.utils.units or float or int, required
+            The temperature to use during the simulation. If no unit is
+            provided, the temperature is assumed to be kT (temperature times
+            Boltzmann constant).
+        tau_kt: float, required
+            Thermostat coupling period (in simulation time units).
+        duration: int or flowermd.utils.units, required
+            The number of steps or time length to run the simulation. If no unit
+            is provided, the time is assumed to be the number of steps.
         couple: str, default "xyz"
             Couplings of diagonal elements of the stress tensor/
         box_dof: list of bool;
@@ -839,6 +911,8 @@ def run_NPT(
             time step.
 
         """
+        self._kT = self._setup_temperature(temperature)
+        _n_steps = self._setup_n_steps(duration)
         self.set_integrator_method(
             integrator_method=hoomd.md.methods.ConstantPressure,
             method_kwargs={
@@ -850,26 +924,26 @@ def run_NPT(
                 "gamma": gamma,
                 "filter": self.integrate_group,
                 "thermostat": self._initialize_thermostat(
-                    {"kT": kT, "tau": tau_kt}
+                    {"kT": self._kT, "tau": tau_kt}
                 ),
             },
         )
         if thermalize_particles:
-            self._thermalize_system(kT)
-        std_out_logger = StdOutLogger(n_steps=n_steps, sim=self)
+            self._thermalize_system(self._kT)
+        std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self)
         std_out_logger_printer = hoomd.update.CustomUpdater(
             trigger=hoomd.trigger.Periodic(self._std_out_freq),
             action=std_out_logger,
         )
         self.operations.updaters.append(std_out_logger_printer)
-        self.run(steps=n_steps, write_at_start=write_at_start)
+        self.run(steps=_n_steps, write_at_start=write_at_start)
         self.operations.updaters.remove(std_out_logger_printer)
 
     def run_NVT(
         self,
-        n_steps,
-        kT,
+        temperature,
         tau_kt,
+        duration,
         thermalize_particles=True,
         write_at_start=True,
     ):
@@ -877,12 +951,15 @@ def run_NVT(
 
         Parameters
         ----------
-        n_steps: int, required
-            Number of steps to run the simulation.
-        kT: int or hoomd.variant.Ramp, required
-            The temperature to use during the simulation.
+        temperature: flowermd.utils.units or float or int, required
+            The temperature to use during the simulation. If no unit is
+            provided, the temperature is assumed to be kT (temperature times
+            Boltzmann constant).
         tau_kt: float, required
             Thermostat coupling period (in simulation time units).
+        duration: int or flowermd.utils.units, required
+            The number of steps or time length to run the simulation. If no unit
+            is provided, the time is assumed to be the number of steps.
         thermalize_particles: bool, default True
             When set to True, assigns random velocities to all particles.
         write_at_start : bool, default True
@@ -890,56 +967,89 @@ def run_NVT(
             for the initial step to execute before the next simulation
             time step.
 
+        Examples
+        --------
+        In this example, a system is initialized with `Pack` and a simulation
+        is run in the NVT ensemble at 300 K for 1 ns.
+
+        ::
+
+            import unyt
+            from flowermd.base import Pack, Simulation
+            from flowermd.library import PPS, OPLS_AA_PPS
+
+            pps_mols = PPS(num_mols=20, lengths=15)
+            pps_system = Pack(
+                molecules=[pps_mols],
+                density=0.5,
+            )
+            pps_system.apply_forcefield(
+                r_cut=2.5,
+                force_field=OPLS_AA_PPS(),
+                auto_scale=True,
+                scale_charges=True
+            )
+            sim = Simulation.from_system(pps_system)
+            sim.run_NVT(
+                temperature=300 * flowermd.utils.units.K,
+                tau_kt=1.0,
+                duration= 1 * flowermd.utils.units.ns,
+            )
+
         """
+        self._kT = self._setup_temperature(temperature)
+        _n_steps = self._setup_n_steps(duration)
         self.set_integrator_method(
             integrator_method=hoomd.md.methods.ConstantVolume,
             method_kwargs={
                 "thermostat": self._initialize_thermostat(
-                    {"kT": kT, "tau": tau_kt}
+                    {"kT": self._kT, "tau": tau_kt}
                 ),
                 "filter": self.integrate_group,
             },
         )
         if thermalize_particles:
-            self._thermalize_system(kT)
-        std_out_logger = StdOutLogger(n_steps=n_steps, sim=self)
+            self._thermalize_system(self._kT)
+        std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self)
         std_out_logger_printer = hoomd.update.CustomUpdater(
             trigger=hoomd.trigger.Periodic(self._std_out_freq),
             action=std_out_logger,
         )
         self.operations.updaters.append(std_out_logger_printer)
-        self.run(steps=n_steps, write_at_start=write_at_start)
+        self.run(steps=_n_steps, write_at_start=write_at_start)
         self.operations.updaters.remove(std_out_logger_printer)
 
-    def run_NVE(self, n_steps, write_at_start=True):
+    def run_NVE(self, duration, write_at_start=True):
         """Run the simulation in the NVE ensemble.
 
         Parameters
         ----------
-        n_steps: int, required
-            Number of steps to run the simulation.
+        duration : int or flowermd.utils.units, required
+            The number of steps or time length to run the simulation. If no unit
+            is provided, the time is assumed to be the number of steps.
         write_at_start : bool, default True
             When set to True, triggers writers that evaluate to True
             for the initial step to execute before the next simulation
             time step.
 
         """
+        _n_steps = self._setup_n_steps(duration)
         self.set_integrator_method(
             integrator_method=hoomd.md.methods.ConstantVolume,
             method_kwargs={"filter": self.integrate_group},
         )
-        std_out_logger = StdOutLogger(n_steps=n_steps, sim=self)
+        std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self)
         std_out_logger_printer = hoomd.update.CustomUpdater(
             trigger=hoomd.trigger.Periodic(self._std_out_freq),
             action=std_out_logger,
         )
         self.operations.updaters.append(std_out_logger_printer)
-        self.run(steps=n_steps, write_at_start=write_at_start)
+        self.run(steps=_n_steps, write_at_start=write_at_start)
         self.operations.updaters.remove(std_out_logger_printer)
 
     def run_displacement_cap(
         self,
-        n_steps,
+        duration,
         maximum_displacement=1e-3,
         write_at_start=True,
     ):
@@ -953,17 +1063,18 @@ def run_displacement_cap(
 
         Parameters
         ----------
-        n_steps : int, required
-            Number of steps to run the simulation.
+        duration : int or flowermd.utils.units, required
+            The number of steps or time length to run the simulation. If no unit
+            is provided, the time is assumed to be the number of steps.
         maximum_displacement : float, default 1e-3
             Maximum displacement per step (length)
-
         write_at_start : bool, default True
             When set to True, triggers writers that evaluate to True
             for the initial step to execute before the next simulation
             time step.
 
         """
+        _n_steps = self._setup_n_steps(duration)
         self.set_integrator_method(
             integrator_method=hoomd.md.methods.DisplacementCapped,
             method_kwargs={
@@ -971,30 +1082,44 @@ def run_displacement_cap(
                 "maximum_displacement": maximum_displacement,
             },
         )
-        std_out_logger = StdOutLogger(n_steps=n_steps, sim=self)
+        std_out_logger = StdOutLogger(n_steps=_n_steps, sim=self)
         std_out_logger_printer = hoomd.update.CustomUpdater(
             trigger=hoomd.trigger.Periodic(self._std_out_freq),
             action=std_out_logger,
         )
         self.operations.updaters.append(std_out_logger_printer)
-        self.run(steps=n_steps, write_at_start=write_at_start)
+        self.run(steps=_n_steps, write_at_start=write_at_start)
         self.operations.updaters.remove(std_out_logger_printer)
 
-    def temperature_ramp(self, n_steps, kT_start, kT_final):
+    def temperature_ramp(
+        self,
+        duration,
+        temperature_start,
+        temperature_final,
+    ):
         """Create a temperature ramp.
 
         Parameters
         ----------
-        n_steps : int, required
-            The number of steps to ramp the temperature over.
-        kT_start : float, required
-            The starting temperature.
-        kT_final : float, required
-            The final temperature.
+        duration : int or flowermd.utils.units, required
+            The number of steps or time length to run the simulation. If no unit
+            is provided, the time is assumed to be the number of steps.
+        temperature_start : unyt.unyt_quantity or float, required
+            The initial temperature. If no unit is provided, the temperature is
+            assumed to be kT (temperature times Boltzmann constant).
+        temperature_final : unyt.unyt_quantity or float, required
+            The final temperature. If no unit is provided, the temperature is
+            assumed to be kT (temperature times Boltzmann constant).
 
         """
+        _kT_start = self._setup_temperature(temperature_start)
+        _kT_final = self._setup_temperature(temperature_final)
+        _n_steps = self._setup_n_steps(duration)
         return hoomd.variant.Ramp(
-            A=kT_start, B=kT_final, t_start=self.timestep, t_ramp=int(n_steps)
+            A=_kT_start,
+            B=_kT_final,
+            t_start=self.timestep,
+            t_ramp=int(_n_steps),
         )
 
     def pickle_forcefield(self, file_path="forcefield.pickle"):
diff --git a/flowermd/base/system.py b/flowermd/base/system.py
index 6416be53..d986e477 100644
--- a/flowermd/base/system.py
+++ b/flowermd/base/system.py
@@ -14,7 +14,7 @@
 
 from flowermd.base.forcefield import BaseHOOMDForcefield, BaseXMLForcefield
 from flowermd.base.molecule import Molecule
-from flowermd.internal import check_return_iterable, validate_ref_value
+from flowermd.internal import check_return_iterable, validate_unit
 from flowermd.internal.exceptions import ForceFieldError, MoleculeLoadError
 from flowermd.utils import (
     get_target_box_mass_density,
@@ -208,12 +208,10 @@ def reference_length(self, length):
 
         Parameters
         ----------
-        length : string or unyt.unyt_quantity, required
+        length : reference length * `flowermd.Units`, required
             The reference length of the system.
-            It can be provided in the following forms:
-            1) A string with the format of "value unit", for example "1 nm".
-            2) A unyt.unyt_quantity object with the correct dimension. For
-            example, unyt.unyt_quantity(1, "nm").
+            It can be provided in the following form of:
+            value * `flowermd.Units`, for example 1 * `flowermd.Units.angstrom`.
 
         """
         if self.auto_scale:
@@ -222,7 +220,7 @@ def reference_length(self, length):
                 "Setting reference length manually disables auto "
                 "scaling."
             )
-        validated_length = validate_ref_value(length, u.dimensions.length)
+        validated_length = validate_unit(length, u.dimensions.length)
         self._reference_values["length"] = validated_length
 
     @reference_energy.setter
@@ -231,12 +229,10 @@ def reference_energy(self, energy):
 
         Parameters
         ----------
-        energy : string or unyt.unyt_quantity, required
+        energy : reference energy * `flowermd.Units`, required
             The reference energy of the system.
-            It can be provided in the following forms:
-            1) A string with the format of "value unit", for example "1 kJ/mol".
-            2) A unyt.unyt_quantity object with the correct dimension. For
-            example, unyt.unyt_quantity(1, "kJ/mol").
+            It can be provided in the following form of:
+            value * `flowermd.Units`, for example 1 * `flowermd.Units.kcal/mol`.
 
         """
         if self.auto_scale:
@@ -245,7 +241,7 @@ def reference_energy(self, energy):
                 "Setting reference energy manually disables auto "
                 "scaling."
             )
-        validated_energy = validate_ref_value(energy, u.dimensions.energy)
+        validated_energy = validate_unit(energy, u.dimensions.energy)
         self._reference_values["energy"] = validated_energy
 
     @reference_mass.setter
@@ -254,13 +250,10 @@ def reference_mass(self, mass):
 
         Parameters
         ----------
-        mass : string or unyt.unyt_quantity, required
+        mass : reference mass * `flowermd.Units`, required
             The reference mass of the system.
-            It can be provided in the following forms:
-            1) A string with the format of "value unit", for example "1 amu".
-            2) A unyt.unyt_quantity object with the correct dimension. For
-            example, unyt.unyt_quantity(1, "amu").
-
+            It can be provided in the following form of:
+            value * `flowermd.Units`, for example 1 * `flowermd.Units.amu`.
         """
         if self.auto_scale:
             warnings.warn(
@@ -268,7 +261,7 @@ def reference_mass(self, mass):
                 "Setting reference mass manually disables auto "
                 "scaling."
             )
-        validated_mass = validate_ref_value(mass, u.dimensions.mass)
+        validated_mass = validate_unit(mass, u.dimensions.mass)
         self._reference_values["mass"] = validated_mass
 
     @reference_values.setter
@@ -404,6 +397,23 @@ def to_gsd(self, file_name):
         with gsd.hoomd.open(file_name, "w") as traj:
             traj.append(self.hoomd_snapshot)
 
+    def save_reference_values(self, file_path="reference_values.pickle"):
+        """Save the reference values of the system to a pickle file.
+
+        Parameters
+        ----------
+        file_path : str, default "reference_values.pickle"
+            The path to save the pickle file to.
+
+        """
+        if not self.reference_values:
+            raise ValueError(
+                "Reference values have not been set. "
+                "See System.reference_values"
+            )
+        f = open(file_path, "wb")
+        pickle.dump(self.reference_values, f)
+
     def _convert_to_gmso(self):
         """Convert the mbuild system to a gmso system."""
         topology = from_mbuild(self.system)
diff --git a/flowermd/internal/__init__.py b/flowermd/internal/__init__.py
index d449c3b3..68f103a2 100644
--- a/flowermd/internal/__init__.py
+++ b/flowermd/internal/__init__.py
@@ -1,2 +1,3 @@
 from .ff_utils import xml_to_gmso_ff
-from .utils import check_return_iterable, validate_ref_value
+from .units import Units
+from .utils import check_return_iterable, validate_unit
diff --git a/flowermd/internal/exceptions.py b/flowermd/internal/exceptions.py
index 8ce79f6d..90d39905 100644
--- a/flowermd/internal/exceptions.py
+++ b/flowermd/internal/exceptions.py
@@ -53,7 +53,7 @@ def __init__(self, msg):
         super().__init__(msg)
 
 
-class ReferenceUnitError(Exception):
+class UnitError(Exception):
     def __init__(self, msg):
         super().__init__(msg)
 
diff --git a/flowermd/internal/units.py b/flowermd/internal/units.py
new file mode 100644
index 00000000..c80a43e6
--- /dev/null
+++ b/flowermd/internal/units.py
@@ -0,0 +1,67 @@
+"""FlowerMD Unit class."""
+
+import unyt as u
+
+
+class Units:
+    """FlowerMD Unit class.
+
+    Example usage:
+    --------------
+
+    ::
+        length = 1.0 * flowermd.Units.angstrom
+        energy = 1.0 * flowermd.Units.kcal_mol
+        mass = 1.0 * flowermd.Units.amu
+        time = 1.0 * flowermd.Units.ps
+        temperature = 1.0 * flowermd.Units.K
+
+
+    """
+
+    # length units
+    m = u.m
+    meter = u.m
+    cm = u.cm
+    centimeter = u.cm
+    nm = u.nm
+    nanometer = u.nm
+    pm = u.pm
+    picometer = u.pm
+    ang = u.angstrom
+    angstrom = u.angstrom
+
+    # energy units
+    J = u.J
+    Joule = u.J
+    kJ = u.kJ
+    cal = u.cal
+    kcal = u.kcal
+    kcal_mol = u.kcal / u.mol
+    kJ_mol = u.kJ / u.mol
+
+    # mass units
+    g = u.g
+    gram = u.g
+    amu = u.amu
+    kg = u.kg
+
+    mol = u.mol
+
+    # time units
+    s = u.s
+    second = u.s
+    ns = u.ns
+    nanosecond = u.ns
+    ps = u.ps
+    picosecond = u.ps
+    fs = u.fs
+    femtosecond = u.fs
+
+    # temperature units
+    K = u.K
+    Kelvin = u.K
+    C = u.degC
+    Celsius = u.degC
+    F = u.degF
+    Fahrenheit = u.degF
diff --git a/flowermd/internal/utils.py b/flowermd/internal/utils.py
index a8f6bea1..cd652b13 100644
--- a/flowermd/internal/utils.py
+++ b/flowermd/internal/utils.py
@@ -1,6 +1,6 @@
 import unyt as u
 
-from flowermd.internal.exceptions import ReferenceUnitError
+from flowermd.internal.exceptions import UnitError
 
 """utils.py
    Internal utility methods for flowerMD.
@@ -19,77 +19,41 @@ def check_return_iterable(obj):
         return [obj]
 
 
-def validate_ref_value(ref_value, dimension):
-    """Validate the reference value and checks the unit dimension.
-    This function validates the reference value. The reference value can be
-    provided in three ways:
-        1. An unyt_quantity instance.
-        2. A string with the value and unit , for example "1.0 g".
-        3. A string with the value and unit separated by a "/", for example
-        "1.0 kcal/mol".
+def validate_unit(value, dimension):
+    """Validate the unit and checks the unit dimension.
+
     Parameters
     ----------
-    ref_value : unyt_quantity or str; required
-        The reference value.
+    value : `unit value * flowermd.Units`; required
+        The unit value to be validated.
     dimension : unyt_dimension; required
-        The dimension of the reference value.
-
-    Returns
-    -------
-    The validated reference value as an unyt.unyt_quantity instance.
+        The dimension of the unit.
     """
 
-    def _is_valid_dimension(ref_unit):
-        if ref_unit.dimensions != dimension:
-            raise ReferenceUnitError(
-                f"Invalid unit dimension. The reference "
-                f"value must be in {dimension} "
-                f"dimension."
+    def _sample_unit_str(dimension):
+        if dimension == u.dimensions.temperature:
+            return "temperature = 300 * flowermd.Units.K"
+        elif dimension == u.dimensions.mass:
+            return "mass = 1.0 * flowermd.Units.g"
+        elif dimension == u.dimensions.length:
+            return "length = 1.0 * flowermd.Units.angstrom"
+        elif dimension == u.dimensions.time:
+            return "time = 1.0 * flowermd.Units.ps"
+        elif dimension == u.dimensions.energy:
+            return "energy = 1.0 * flowermd.Units.kcal_mol"
+
+    if isinstance(value, u.unyt_quantity):
+        if value.units.dimensions != dimension:
+            raise UnitError(
+                f"Invalid unit dimension. The unit must be in "
+                f"{dimension} dimension. Check `flowermd.Units` for "
+                f"valid units."
             )
-        return True
-
-    def _parse_and_validate_unit(value, unit_str):
-        if hasattr(u, unit_str):
-            if unit_str == "amu":
-                u_unit = u.Unit("amu")
-            else:
-                u_unit = getattr(u, unit_str)
-            if _is_valid_dimension(u_unit):
-                return float(value) * u_unit
-        # if the unit contains "/" character, for example "g/mol", check if
-        # the unit is a valid unit and has the correct dimension.
-        if len(unit_str.split("/")) == 2:
-            unit1, unit2 = unit_str.split("/")
-            if hasattr(u, unit1) and hasattr(u, unit2):
-                comb_unit = getattr(u, unit1) / getattr(u, unit2)
-                if _is_valid_dimension(comb_unit):
-                    return float(value) * comb_unit
-        raise ReferenceUnitError(
-            f"Invalid reference value. Please provide "
-            f"a reference value with unit of "
-            f"{dimension} dimension."
-        )
-
-    def _is_float(num):
-        try:
-            return float(num)
-        except ValueError:
-            raise ValueError("The reference value is not a number.")
-
-    # if ref_value is an instance of unyt_quantity, check the dimension.
-    if isinstance(ref_value, u.unyt_quantity) and _is_valid_dimension(
-        ref_value.units
-    ):
-        return ref_value
-    # if ref_value is a string, check if it is a number and if it is, check if
-    # the unit exists in unyt and has the correct dimension.
-    elif isinstance(ref_value, str) and len(ref_value.split()) == 2:
-        value, unit_str = ref_value.split()
-        value = _is_float(value)
-        return _parse_and_validate_unit(value, unit_str)
+        return value
     else:
-        raise ReferenceUnitError(
-            f"Invalid reference value. Please provide "
-            f"a reference value with unit of "
-            f"{dimension} dimension."
+        raise UnitError(
+            "The unit value must be provided from the "
+            "`flowermd.Units` class. For example, "
+            f"{_sample_unit_str(dimension)}. Check "
+            "`flowermd.Units` for valid units."
         )
diff --git a/flowermd/tests/base/test_simulation.py b/flowermd/tests/base/test_simulation.py
index a756f808..56aadbf1 100644
--- a/flowermd/tests/base/test_simulation.py
+++ b/flowermd/tests/base/test_simulation.py
@@ -6,9 +6,8 @@
 import hoomd
 import numpy as np
 import pytest
-import unyt as u
 
-from flowermd import Simulation
+from flowermd import Simulation, Units
 from flowermd.base import Pack
 from flowermd.library import OPLS_AA_PPS
 from flowermd.library.forcefields import EllipsoidForcefield
@@ -20,7 +19,7 @@
 class TestSimulate(BaseTest):
     def test_initialize_from_system(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
-        sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500)
+        sim.run_NVT(temperature=1.0, tau_kt=0.01, duration=500)
         assert len(sim.forces) == len(benzene_system.hoomd_forcefield)
         assert sim.reference_values == benzene_system.reference_values
 
@@ -30,7 +29,7 @@ def test_initialize_from_system_separate_ff(
         sim = Simulation.from_system(
             benzene_cg_system, forcefield=cg_single_bead_ff
         )
-        sim.run_NVT(kT=0.1, tau_kt=10, n_steps=500)
+        sim.run_NVT(temperature=0.1, tau_kt=10, duration=500)
 
     def test_initialize_from_system_missing_ff(self, benzene_cg_system):
         with pytest.raises(ValueError):
@@ -70,9 +69,9 @@ def test_set_ref_values(self, benzene_system):
             forcefield=benzene_system.hoomd_forcefield,
         )
         ref_value_dict = {
-            "length": 1 * u.angstrom,
-            "energy": 3.0 * u.kcal / u.mol,
-            "mass": 1.25 * u.Unit("amu"),
+            "length": 1 * Units.angstrom,
+            "energy": 3.0 * Units.kcal_mol,
+            "mass": 1.25 * Units.amu,
         }
         sim.reference_values = ref_value_dict
         assert sim.reference_length == ref_value_dict["length"]
@@ -84,35 +83,55 @@ def test_set_ref_length(self, benzene_system):
             initial_state=benzene_system.hoomd_snapshot,
             forcefield=benzene_system.hoomd_forcefield,
         )
-        sim.reference_length = 1 * u.angstrom
-        assert sim.reference_length == 1 * u.angstrom
+        sim.reference_length = 1 * Units.angstrom
+        assert sim.reference_length == 1 * Units.angstrom
 
     def test_set_ref_energy(self, benzene_system):
         sim = Simulation(
             initial_state=benzene_system.hoomd_snapshot,
             forcefield=benzene_system.hoomd_forcefield,
         )
-        sim.reference_energy = 3.0 * u.kcal / u.mol
-        assert sim.reference_energy == 3.0 * u.kcal / u.mol
+        sim.reference_energy = 3.0 * Units.kcal_mol
+        assert sim.reference_energy == 3.0 * Units.kcal_mol
 
     def test_set_ref_mass(self, benzene_system):
         sim = Simulation(
             initial_state=benzene_system.hoomd_snapshot,
             forcefield=benzene_system.hoomd_forcefield,
         )
-        sim.reference_mass = 1.25 * u.amu
-        assert sim.reference_mass == 1.25 * u.amu
+        sim.reference_mass = 1.25 * Units.amu
+        assert sim.reference_mass == 1.25 * Units.amu
 
     def test_NVT(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
-        sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500)
+        sim.run_NVT(temperature=1.0, tau_kt=0.01, duration=500)
         assert isinstance(sim.method, hoomd.md.methods.ConstantVolume)
 
+    def test_NVT_real_units(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        sim.run_NVT(
+            temperature=35.225 * Units.K, tau_kt=0.01, duration=1 * Units.ps
+        )
+        assert isinstance(sim.method, hoomd.md.methods.ConstantVolume)
+        assert np.isclose(sim._kT, 1.0, atol=1e-1)
+        assert sim.timestep == int(1 * Units.ps / sim.real_timestep)
+
     def test_NPT(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         sim.run_NPT(
-            kT=1.0,
-            n_steps=500,
+            temperature=1.0,
+            duration=500,
+            pressure=0.0001,
+            tau_kt=0.001,
+            tau_pressure=0.01,
+        )
+        assert isinstance(sim.method, hoomd.md.methods.ConstantPressure)
+
+    def test_NPT_real_units(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        sim.run_NPT(
+            temperature=200.0 * Units.K,
+            duration=0.1 * Units.ps,
             pressure=0.0001,
             tau_kt=0.001,
             tau_pressure=0.01,
@@ -121,25 +140,54 @@ def test_NPT(self, benzene_system):
 
     def test_langevin(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
-        sim.run_langevin(n_steps=500, kT=1.0)
+        sim.run_langevin(duration=500, temperature=1.0)
+        assert isinstance(sim.method, hoomd.md.methods.Langevin)
+
+    def test_langevin_real_units(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        sim.run_langevin(
+            duration=0.1 * Units.ps, temperature=10.0 * Units.Celsius
+        )
         assert isinstance(sim.method, hoomd.md.methods.Langevin)
 
     def test_NVE(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
-        sim.run_NVE(n_steps=500)
+        sim.run_NVE(duration=500)
+        assert isinstance(sim.method, hoomd.md.methods.ConstantVolume)
+
+    def test_NVE_real_units(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        sim.run_NVE(duration=0.1 * Units.ps)
         assert isinstance(sim.method, hoomd.md.methods.ConstantVolume)
 
     def test_displacement_cap(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
-        sim.run_displacement_cap(n_steps=500, maximum_displacement=1e-4)
+        sim.run_displacement_cap(duration=500, maximum_displacement=1e-4)
+        assert isinstance(sim.method, hoomd.md.methods.DisplacementCapped)
+
+    def test_displacement_cap_real_units(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        sim.run_displacement_cap(
+            duration=01.0 * Units.ps, maximum_displacement=1e-4
+        )
         assert isinstance(sim.method, hoomd.md.methods.DisplacementCapped)
 
     def test_update_volume_target_box(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         sim.run_update_volume(
-            kT=1.0,
+            temperature=1.0,
+            tau_kt=0.01,
+            duration=500,
+            period=1,
+            final_box_lengths=sim.box_lengths_reduced * 0.5,
+        )
+
+    def test_update_volume_real_units(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        sim.run_update_volume(
+            temperature=100.0 * Units.K,
             tau_kt=0.01,
-            n_steps=500,
+            duration=0.2 * Units.ps,
             period=1,
             final_box_lengths=sim.box_lengths_reduced * 0.5,
         )
@@ -148,9 +196,9 @@ def test_update_volume_walls(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         sim.add_walls(wall_axis=(1, 0, 0), sigma=1.0, epsilon=1.0, r_cut=1.12)
         sim.run_update_volume(
-            kT=1.0,
+            temperature=1.0,
             tau_kt=0.01,
-            n_steps=500,
+            duration=500,
             period=5,
             final_box_lengths=sim.box_lengths_reduced * 0.5,
         )
@@ -164,9 +212,9 @@ def test_update_volume_no_units(self, benzene_system):
         init_box = sim.box_lengths_reduced
         sim.run_update_volume(
             final_box_lengths=init_box / 2,
-            kT=1.0,
+            temperature=1.0,
             tau_kt=0.01,
-            n_steps=500,
+            duration=500,
             period=1,
         )
         assert np.allclose(sim.box_lengths_reduced * 2, init_box)
@@ -174,18 +222,18 @@ def test_update_volume_no_units(self, benzene_system):
     def test_update_volume_density(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         target_box = get_target_box_mass_density(
-            density=0.05 * u.Unit("g") / u.Unit("cm**3"), mass=sim.mass.to(u.g)
+            density=0.05 * Units.g / Units.cm**3, mass=sim.mass.to(Units.g)
         )
         sim.run_update_volume(
-            kT=1.0,
+            temperature=1.0,
             tau_kt=0.01,
-            n_steps=500,
+            duration=500,
             period=1,
             final_box_lengths=target_box,
         )
         assert np.isclose(
-            sim.density.to(u.g / u.cm**3).value,
-            (0.05 * (u.g / u.cm**3)).value,
+            sim.density.to(Units.g / Units.cm**3).value,
+            (0.05 * (Units.g / Units.cm**3)).value,
             atol=1e-4,
         )
 
@@ -193,12 +241,12 @@ def test_update_volume_by_density_factor(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
         init_density = copy.deepcopy(sim.density)
         target_box = get_target_box_mass_density(
-            density=init_density * 5, mass=sim.mass.to(u.g)
+            density=init_density * 5, mass=sim.mass.to(Units.g)
         )
         sim.run_update_volume(
-            kT=1.0,
+            temperature=1.0,
             tau_kt=0.01,
-            n_steps=500,
+            duration=500,
             period=1,
             final_box_lengths=target_box,
         )
@@ -208,18 +256,22 @@ def test_update_volume_by_density_factor(self, benzene_system):
 
     def test_change_methods(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
-        sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)
+        sim.run_NVT(temperature=1.0, tau_kt=0.01, duration=0)
         assert isinstance(sim.method, hoomd.md.methods.ConstantVolume)
         sim.run_NPT(
-            kT=1.0, tau_kt=0.01, tau_pressure=0.1, pressure=0.001, n_steps=0
+            temperature=1.0,
+            tau_kt=0.01,
+            tau_pressure=0.1,
+            pressure=0.001,
+            duration=0,
         )
         assert isinstance(sim.method, hoomd.md.methods.ConstantPressure)
 
     def test_change_dt(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
-        sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)
+        sim.run_NVT(temperature=1.0, tau_kt=0.01, duration=0)
         sim.dt = 0.003
-        sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=0)
+        sim.run_NVT(temperature=1.0, tau_kt=0.01, duration=0)
         assert sim.dt == 0.003
 
     def test_scale_epsilon(self, benzene_system):
@@ -282,7 +334,7 @@ def test_set_integrate_group(self, benzene_system):
         tag_filter = hoomd.filter.Tags([0, 1, 2, 3])
         sim.integrate_group = tag_filter
         assert not isinstance(sim.integrate_group, hoomd.filter.All)
-        sim.run_NVT(n_steps=200, kT=1.0, tau_kt=0.01)
+        sim.run_NVT(duration=200, temperature=1.0, tau_kt=0.01)
 
     def test_pickle_ff(self, benzene_system):
         sim = Simulation.from_system(benzene_system)
@@ -333,7 +385,7 @@ def test_rigid_sim(self):
             forcefield=ellipsoid_ff.hoomd_forces,
             rigid_constraint=rigid,
         )
-        sim.run_NVT(n_steps=0, kT=1.0, tau_kt=sim.dt * 100)
+        sim.run_NVT(duration=0, temperature=1.0, tau_kt=sim.dt * 100)
         assert sim.integrator.integrate_rotational_dof is True
         assert sim.mass_reduced == 800.0
 
@@ -352,7 +404,7 @@ def test_save_restart_gsd(self, benzene_system):
 
     def test_gsd_logger(self, benzene_system):
         sim = Simulation.from_system(benzene_system, gsd_write_freq=1)
-        sim.run_NVT(n_steps=5, kT=1.0, tau_kt=0.001)
+        sim.run_NVT(duration=5, temperature=1.0, tau_kt=0.001)
         sim.operations.writers[-2].flush()
         expected_gsd_quantities = [
             "flowermd/base/simulation/Simulation/timestep",
@@ -420,8 +472,27 @@ def test_bad_ff(self, benzene_system):
 
     def test_flush(self, benzene_system):
         sim = Simulation.from_system(benzene_system, gsd_write_freq=100)
-        sim.run_NVT(kT=1.0, tau_kt=0.01, n_steps=500, write_at_start=False)
+        sim.run_NVT(
+            temperature=1.0, tau_kt=0.01, duration=500, write_at_start=False
+        )
         sim.flush_writers()
         with gsd.hoomd.open("trajectory.gsd") as traj:
             assert len(traj) > 0
         os.remove("trajectory.gsd")
+
+    def test_real_temperature(self, benzene_system):
+        sim = Simulation.from_system(benzene_system)
+        with pytest.raises(ValueError):
+            sim.real_temperature
+        sim.run_NVT(temperature=1.0, tau_kt=0.01, duration=100)
+        assert sim.real_temperature.units == Units.K
+        assert np.isclose(sim.real_temperature, 35.225, atol=1e-4)
+
+    def test_real_temperature_no_energy_units(self, benzene_system):
+        sim = Simulation(
+            initial_state=benzene_system.hoomd_snapshot,
+            forcefield=benzene_system.hoomd_forcefield,
+            reference_values=dict(),
+        )
+        sim.run_NVT(temperature=1e-10, tau_kt=0.01, duration=100)
+        assert np.isclose(sim.real_temperature, 7.2429e12)
diff --git a/flowermd/tests/base/test_system.py b/flowermd/tests/base/test_system.py
index cade5f96..5ec020d3 100644
--- a/flowermd/tests/base/test_system.py
+++ b/flowermd/tests/base/test_system.py
@@ -1,15 +1,15 @@
 import os
+import pickle
 
 import gmso
 import hoomd
 import numpy as np
 import pytest
-import unyt as u
 from cmeutils.geometry import get_backbone_vector
 from unyt import Unit
 
-from flowermd import Lattice, Pack
-from flowermd.internal.exceptions import ForceFieldError, ReferenceUnitError
+from flowermd import Lattice, Pack, Units
+from flowermd.internal.exceptions import ForceFieldError, UnitError
 from flowermd.library import (
     OPLS_AA,
     OPLS_AA_DIMETHYLETHER,
@@ -263,9 +263,9 @@ def test_rebuild_snapshot(self, polyethylene):
         init_snap = system.hoomd_snapshot
         new_snap = system.hoomd_snapshot
         assert init_snap == new_snap
-        system.reference_length = 1 * u.angstrom
-        system.reference_energy = 1 * u.kcal / u.mol
-        system.reference_mass = 1 * u.amu
+        system.reference_length = 1 * Units.angstrom
+        system.reference_energy = 1 * Units.kcal_mol
+        system.reference_mass = 1 * Units.amu
         assert system._snap_refs != system.reference_values
         assert system._ff_refs != system.reference_values
         new_snap = system.hoomd_snapshot
@@ -280,9 +280,9 @@ def test_ref_values_no_autoscale(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        system.reference_length = 1 * u.angstrom
-        system.reference_energy = 1 * u.kcal / u.mol
-        system.reference_mass = 1 * u.amu
+        system.reference_length = 1 * Units.angstrom
+        system.reference_energy = 1 * Units.kcal_mol
+        system.reference_mass = 1 * Units.amu
         assert np.allclose(
             system.hoomd_snapshot.configuration.box[:3],
             system.gmso_system.box.lengths.to("angstrom").value,
@@ -300,34 +300,15 @@ def test_set_ref_values(self, polyethylene):
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
         ref_value_dict = {
-            "length": 1 * u.angstrom,
-            "energy": 3.0 * u.kcal / u.mol,
-            "mass": 1.25 * u.Unit("amu"),
+            "length": 1 * Units.angstrom,
+            "energy": 3.0 * Units.kcal_mol,
+            "mass": 1.25 * Units.amu,
         }
         system.reference_values = ref_value_dict
         assert system.reference_length == ref_value_dict["length"]
         assert system.reference_energy == ref_value_dict["energy"]
         assert system.reference_mass == ref_value_dict["mass"]
 
-    def test_set_ref_values_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        ref_value_dict = {
-            "length": "1 angstrom",
-            "energy": "3 kcal/mol",
-            "mass": "1.25 amu",
-        }
-        system.reference_values = ref_value_dict
-        assert system.reference_length == 1 * u.angstrom
-        assert system.reference_energy == 3 * u.kcal / u.mol
-        assert system.reference_mass == 1.25 * u.amu
-
     def test_set_ref_values_missing_key(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
         system = Pack(
@@ -338,8 +319,8 @@ def test_set_ref_values_missing_key(self, polyethylene):
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
         ref_value_dict = {
-            "length": 1 * u.angstrom,
-            "energy": 3.0 * u.kcal / u.mol,
+            "length": 1 * Units.angstrom,
+            "energy": 3.0 * Units.kcal_mol,
         }
         with pytest.raises(ValueError):
             system.reference_values = ref_value_dict
@@ -354,11 +335,11 @@ def test_set_ref_values_invalid_type(self, polyethylene):
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
         ref_value_dict = {
-            "length": 1 * u.angstrom,
-            "energy": 3.0 * u.kcal / u.mol,
+            "length": 1 * Units.angstrom,
+            "energy": 3.0 * Units.kcal_mol,
             "mass": 1.25,
         }
-        with pytest.raises(ReferenceUnitError):
+        with pytest.raises(UnitError):
             system.reference_values = ref_value_dict
 
     def test_set_ref_values_auto_scale_true(self, polyethylene):
@@ -371,9 +352,9 @@ def test_set_ref_values_auto_scale_true(self, polyethylene):
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True
         )
         ref_value_dict = {
-            "length": 1 * u.angstrom,
-            "energy": 3.0 * u.kcal / u.mol,
-            "mass": 1.25 * u.Unit("amu"),
+            "length": 1 * Units.angstrom,
+            "energy": 3.0 * Units.kcal_mol,
+            "mass": 1.25 * Units.amu,
         }
         with pytest.warns():
             system.reference_values = ref_value_dict
@@ -387,8 +368,8 @@ def test_set_ref_length(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        system.reference_length = 1 * u.angstrom
-        assert system.reference_length == 1 * u.angstrom
+        system.reference_length = 1 * Units.angstrom
+        assert system.reference_length == 1 * Units.angstrom
 
     def test_set_ref_length_invalid_type(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=5)
@@ -399,45 +380,9 @@ def test_set_ref_length_invalid_type(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        with pytest.raises(ReferenceUnitError):
+        with pytest.raises(UnitError):
             system.reference_length = 1.0
 
-    def test_ref_length_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        system.reference_length = "1 angstrom"
-        assert system.reference_length == 1 * u.angstrom
-
-    def test_ref_length_invalid_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_length = "1.0"
-
-    def test_ref_length_invalid_unit_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_length = "1.0 invalid_unit"
-
     def test_ref_length_invalid_dimension(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
         system = Pack(
@@ -447,20 +392,8 @@ def test_ref_length_invalid_dimension(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_length = 1.0 * u.g
-
-    def test_ref_length_invalid_dimension_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_length = "1.0 g"
+        with pytest.raises(UnitError):
+            system.reference_length = 1.0 * Units.g
 
     def test_ref_length_auto_scale_true(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
@@ -471,8 +404,8 @@ def test_ref_length_auto_scale_true(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True
         )
-        system.reference_length = 1 * u.angstrom
-        assert system.reference_length == 1 * u.angstrom
+        system.reference_length = 1 * Units.angstrom
+        assert system.reference_length == 1 * Units.angstrom
 
     def test_set_ref_energy(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
@@ -483,8 +416,8 @@ def test_set_ref_energy(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        system.reference_energy = 1 * u.kcal / u.mol
-        assert system.reference_energy == 1 * u.kcal / u.mol
+        system.reference_energy = 1 * Units.kcal_mol
+        assert system.reference_energy == 1 * Units.kcal_mol
 
     def test_set_ref_energy_invalid_type(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
@@ -495,57 +428,9 @@ def test_set_ref_energy_invalid_type(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        with pytest.raises(ReferenceUnitError):
+        with pytest.raises(UnitError):
             system.reference_energy = 1.0
 
-    def test_ref_energy_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        system.reference_energy = "1 kJ"
-        assert system.reference_energy == 1 * u.kJ
-
-    def test_ref_energy_string_comb_units(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        system.reference_energy = "1 kcal/mol"
-        assert system.reference_energy == 1 * u.kcal / u.mol
-
-    def test_ref_energy_invalid_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_energy = "1.0"
-
-    def test_ref_energy_invalid_unit_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_energy = "1.0 invalid_unit"
-
     def test_ref_energy_invalid_dimension(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
         system = Pack(
@@ -555,20 +440,8 @@ def test_ref_energy_invalid_dimension(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_energy = 1.0 * u.g
-
-    def test_ref_energy_invalid_dimension_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_energy = "1.0 m"
+        with pytest.raises(UnitError):
+            system.reference_energy = 1.0 * Units.g
 
     def test_set_ref_energy_auto_scale_true(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
@@ -579,8 +452,8 @@ def test_set_ref_energy_auto_scale_true(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True
         )
-        system.reference_energy = 1 * u.kcal / u.mol
-        assert system.reference_energy == 1 * u.kcal / u.mol
+        system.reference_energy = 1 * Units.kcal_mol
+        assert system.reference_energy == 1 * Units.kcal_mol
 
     def test_set_ref_mass(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
@@ -592,8 +465,8 @@ def test_set_ref_mass(self, polyethylene):
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
 
-        system.reference_mass = 1.0 * u.amu
-        assert system.reference_mass == 1.0 * u.amu
+        system.reference_mass = 1.0 * Units.amu
+        assert system.reference_mass == 1.0 * Units.amu
 
     def test_set_ref_mass_invalid_type(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
@@ -604,45 +477,9 @@ def test_set_ref_mass_invalid_type(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        with pytest.raises(ReferenceUnitError):
+        with pytest.raises(UnitError):
             system.reference_mass = 1.0
 
-    def test_ref_mass_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        system.reference_mass = "1 g"
-        assert system.reference_mass == 1.0 * u.g
-
-    def test_ref_mass_invalid_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_mass = "1.0"
-
-    def test_ref_mass_invalid_unit_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_mass = "1.0 invalid_unit"
-
     def test_ref_mass_invalid_dimension(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
         system = Pack(
@@ -652,20 +489,8 @@ def test_ref_mass_invalid_dimension(self, polyethylene):
         system.apply_forcefield(
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
         )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_energy = 1.0 * u.m
-
-    def test_ref_mass_invalid_dimension_string(self, polyethylene):
-        polyethylene = polyethylene(lengths=5, num_mols=1)
-        system = Pack(
-            molecules=[polyethylene],
-            density=1.0,
-        )
-        system.apply_forcefield(
-            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=False
-        )
-        with pytest.raises(ReferenceUnitError):
-            system.reference_mass = "1.0 m"
+        with pytest.raises(UnitError):
+            system.reference_energy = 1.0 * Units.m
 
     def test_set_ref_mass_auto_scale_true(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
@@ -677,8 +502,8 @@ def test_set_ref_mass_auto_scale_true(self, polyethylene):
             r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True
         )
 
-        system.reference_mass = 1.0 * u.amu
-        assert system.reference_mass == 1.0 * u.amu
+        system.reference_mass = 1.0 * Units.amu
+        assert system.reference_mass == 1.0 * Units.amu
 
     def test_apply_forcefield_no_forcefield(self, polyethylene):
         polyethylene = polyethylene(lengths=5, num_mols=1)
@@ -891,28 +716,61 @@ def test_pickle_ff_not_ff(self, polyethylene):
         with pytest.raises(ValueError):
             system.pickle_forcefield("forcefield.pickle")
 
+    def test_save_reference_values(self, polyethylene):
+        polyethylene = polyethylene(lengths=5, num_mols=1)
+        system = Pack(
+            molecules=[polyethylene],
+            density=1.0,
+        )
+        system.apply_forcefield(
+            r_cut=2.5, force_field=[OPLS_AA()], auto_scale=True
+        )
+        system.save_reference_values()
+        assert os.path.isfile(
+            os.path.join(os.getcwd(), "reference_values.pickle")
+        )
+        ref_values = pickle.load(
+            open(os.path.join(os.getcwd(), "reference_values.pickle"), "rb")
+        )
+        assert ref_values == system.reference_values
+        os.remove(os.path.join(os.getcwd(), "reference_values.pickle"))
+
+    def test_save_empty_reference_values(self, polyethylene):
+        polyethylene = polyethylene(lengths=5, num_mols=1)
+        system = Pack(
+            molecules=[polyethylene],
+            density=1.0,
+            base_units=dict(),
+        )
+        with pytest.raises(ValueError):
+            system.save_reference_values()
+
     def test_mass_density(self, benzene_molecule):
         benzene_mol = benzene_molecule(n_mols=100)
-        system = Pack(molecules=[benzene_mol], density=0.8 * u.g / u.cm**3)
-        assert system.density.units == u.g / u.cm**3
+        system = Pack(
+            molecules=[benzene_mol], density=0.8 * Units.g / Units.cm**3
+        )
+        assert system.density.units == Units.g / Units.cm**3
 
     def test_number_density(self, benzene_molecule):
         benzene_mol = benzene_molecule(n_mols=100)
-        system = Pack(molecules=[benzene_mol], density=0.8 / u.cm**3)
-        assert system.density.units == u.cm**-3
+        system = Pack(molecules=[benzene_mol], density=0.8 / Units.cm**3)
+        assert system.density.units == Units.cm**-3
 
     def test_bad_density(self, benzene_molecule):
         benzene_mol = benzene_molecule(n_mols=100)
         with pytest.raises(ValueError):
-            Pack(molecules=[benzene_mol], density=0.8 * u.J / u.kg)
+            Pack(molecules=[benzene_mol], density=0.8 * Units.J / Units.kg)
 
     def test_density_warning(self, benzene_molecule):
         benzene_mol = benzene_molecule(n_mols=100)
         with pytest.warns():
             system = Pack(molecules=[benzene_mol], density=0.8)
-            assert system.density.units == u.g / u.cm**3
+            assert system.density.units == Units.g / Units.cm**3
 
     def test_n_particles_no_ff(self, benzene_molecule):
         benzene_mol = benzene_molecule(n_mols=100)
-        system = Pack(molecules=[benzene_mol], density=0.8 * u.g / u.cm**3)
+        system = Pack(
+            molecules=[benzene_mol], density=0.8 * Units.g / Units.cm**3
+        )
         assert system.n_particles == 100 * 12
diff --git a/flowermd/tests/utils/test_utils.py b/flowermd/tests/utils/test_utils.py
index 505eb0d3..26c4d54b 100644
--- a/flowermd/tests/utils/test_utils.py
+++ b/flowermd/tests/utils/test_utils.py
@@ -2,8 +2,9 @@
 import pytest
 import unyt as u
 
-from flowermd.internal import check_return_iterable, validate_ref_value
-from flowermd.internal.exceptions import ReferenceUnitError
+from flowermd import Units
+from flowermd.internal import check_return_iterable, validate_unit
+from flowermd.internal.exceptions import UnitError
 from flowermd.utils import (
     _calculate_box_length,
     get_target_box_mass_density,
@@ -24,37 +25,39 @@ def test_check_return_iterable(self):
             {"test": 1, "test2": 2}
         ]
 
-    def test_validate_ref_value(self):
-        assert validate_ref_value(1.0 * u.g, u.dimensions.mass) == 1.0 * u.g
-        assert validate_ref_value("1.0 g", u.dimensions.mass) == 1.0 * u.g
-        assert validate_ref_value("1.0 kcal/mol", u.dimensions.energy) == (
-            1.0 * u.kcal / u.mol
-        )
-        assert validate_ref_value("1.0 amu", u.dimensions.mass) == 1.0 * u.Unit(
-            "amu"
-        )
-        with pytest.raises(ReferenceUnitError):
-            validate_ref_value("1.0 g", u.dimensions.energy)
-
-        with pytest.raises(ReferenceUnitError):
-            validate_ref_value("1.0 kcal/invalid", u.dimensions.energy)
-
-        with pytest.raises(ReferenceUnitError):
-            validate_ref_value("1.0 invalid", u.dimensions.energy)
-
-        with pytest.raises(ValueError):
-            validate_ref_value("test g", u.dimensions.mass)
+    def test_validate_unit(self):
+        value = 1.0 * Units.g
+        dimension = u.dimensions.mass
+        assert validate_unit(value, dimension) == value
+        value = 1.0 * Units.kcal_mol
+        dimension = u.dimensions.energy
+        assert validate_unit(value, dimension) == value
+        value = 1.0 * Units.angstrom
+        dimension = u.dimensions.length
+        assert validate_unit(value, dimension) == value
+        value = 1.0 * Units.ps
+        dimension = u.dimensions.time
+        assert validate_unit(value, dimension) == value
+        value = 300 * Units.K
+        dimension = u.dimensions.temperature
+        assert validate_unit(value, dimension) == value
+        with pytest.raises(UnitError):
+            validate_unit(1.0 * Units.nm, u.dimensions.mass)
+        with pytest.raises(UnitError):
+            validate_unit(1.0 * Units.g, u.dimensions.length)
+        with pytest.raises(UnitError):
+            validate_unit(1.0, u.dimensions.energy)
 
     def test_target_box_mass_density(self):
-        mass = u.unyt_quantity(4.0, u.g)
-        density = u.unyt_quantity(0.5, u.g / u.cm**3)
+        mass = 4 * Units.g
+        density = 0.5 * (Units.g / Units.cm**3)
         target_box = get_target_box_mass_density(density=density, mass=mass)
         assert target_box[0] == target_box[1] == target_box[2]
-        assert np.array_equal(target_box, np.array([2] * 3) * u.cm)
+        assert np.array_equal(target_box, np.array([2] * 3) * Units.cm)
 
     def test_target_box_one_constraint_mass(self):
-        mass = u.unyt_quantity(4.0, u.g)
-        density = u.unyt_quantity(0.5, u.g / u.cm**3)
+        mass = 4 * Units.g
+        density = 0.5 * Units.g / Units.cm**3
         cubic_box = get_target_box_mass_density(density=density, mass=mass)
         tetragonal_box = get_target_box_mass_density(
             density=density, mass=mass, x_constraint=cubic_box[0] / 2
@@ -64,8 +67,8 @@ def test_target_box_one_constraint_mass(self):
         assert tetragonal_box[0] == cubic_box[0] / 2
 
     def test_target_box_two_constraint_mass(self):
-        mass = u.unyt_quantity(4.0, u.g)
-        density = u.unyt_quantity(0.5, u.g / u.cm**3)
+        mass = 4 * Units.g
+        density = 0.5 * (Units.g / Units.cm**3)
         cubic_box = get_target_box_mass_density(density=density, mass=mass)
         ortho_box = get_target_box_mass_density(
             density=density,
@@ -78,7 +81,7 @@ def test_target_box_two_constraint_mass(self):
         assert ortho_box[0] == cubic_box[0] / 2
 
     def test_target_box_number_density(self):
-        sigma = 1 * u.nm
+        sigma = 1 * Units.nm
         n_beads = 100
         density = 1 / sigma**3
         target_box = get_target_box_number_density(
@@ -88,7 +91,7 @@ def test_target_box_number_density(self):
         assert np.allclose(L**3, 100, atol=1e-8)
 
     def test_target_box_one_constraint_number_density(self):
-        sigma = 1 * u.nm
+        sigma = 1 * Units.nm
         n_beads = 100
         density = 1 / sigma**3
         cubic_box = get_target_box_number_density(
@@ -103,7 +106,7 @@ def test_target_box_one_constraint_number_density(self):
         assert np.allclose(tetragonal_box[1].value, 6.56419787945, atol=1e-5)
 
     def test_target_box_two_constraint_number_density(self):
-        sigma = 1 * u.nm
+        sigma = 1 * Units.nm
         n_beads = 100
         density = 1 / sigma**3
         cubic_box = get_target_box_number_density(
@@ -121,27 +124,27 @@ def test_target_box_two_constraint_number_density(self):
         )
 
     def test_calculate_box_length_bad_args(self):
-        mass_density = 1 * u.g / (u.cm**3)
-        number_density = 1 / (1 * u.nm**3)
+        mass_density = 1 * Units.g / (Units.cm**3)
+        number_density = 1 / (1 * Units.nm**3)
         with pytest.raises(ValueError):
             get_target_box_mass_density(density=number_density, mass=100)
         with pytest.raises(ValueError):
             get_target_box_number_density(density=mass_density, n_beads=100)
 
     def test_calculate_box_length_fixed_l_1d(self):
-        mass = u.unyt_quantity(6.0, u.g)
-        density = u.unyt_quantity(0.5, u.g / u.cm**3)
-        fixed_L = u.unyt_quantity(3.0, u.cm)
+        mass = 6.0 * Units.g
+        density = 0.5 * (Units.g / Units.cm**3)
+        fixed_L = 3.0 * Units.cm
         box_length = _calculate_box_length(
             mass=mass, density=density, fixed_L=fixed_L
         )
-        assert box_length == 2.0 * u.cm
+        assert box_length == 2.0 * Units.cm
 
     def test_calculate_box_length_fixed_l_2d(self):
-        mass = u.unyt_quantity(12.0, u.g)
-        density = u.unyt_quantity(0.5, u.g / u.cm**3)
-        fixed_L = u.unyt_array([3.0, 2.0], u.cm)
+        mass = 12.0 * Units.g
+        density = 0.5 * Units.g / Units.cm**3
+        fixed_L = [3.0, 2.0] * Units.cm
         box_length = _calculate_box_length(
             mass=mass, density=density, fixed_L=fixed_L
         )
-        assert box_length == 4.0 * u.cm
+        assert box_length == 4.0 * Units.cm