diff --git a/src/relentless/model/potential/__init__.py b/src/relentless/model/potential/__init__.py
index 7231ccf9..31ddb2a4 100644
--- a/src/relentless/model/potential/__init__.py
+++ b/src/relentless/model/potential/__init__.py
@@ -27,6 +27,18 @@
     BondSpline
 
 
+Angle potentials
+================
+
+.. autosummary::
+    :toctree: generated/
+
+    AngleSpline
+    CosineAngle
+    HarmonicAngle
+    HarmonicCosineAngle
+
+
 Developer classes
 =================
 
@@ -36,6 +48,9 @@
     Potential
     Parameters
     BondedPotential
+    BondedSpline
+    AnglePotential
+    AngleParameters
     BondPotential
     BondParameters
     PairPotential
@@ -43,6 +58,14 @@
 
 """
 
+from .angle import (
+    AngleParameters,
+    AnglePotential,
+    AngleSpline,
+    CosineAngle,
+    HarmonicAngle,
+    HarmonicCosineAngle,
+)
 from .bond import FENEWCA, BondParameters, BondPotential, BondSpline, HarmonicBond
 from .pair import (
     Depletion,
@@ -52,4 +75,4 @@
     PairSpline,
     Yukawa,
 )
-from .potential import BondedPotential, Parameters, Potential
+from .potential import BondedPotential, BondedSpline, Parameters, Potential
diff --git a/src/relentless/model/potential/angle.py b/src/relentless/model/potential/angle.py
new file mode 100644
index 00000000..cd975694
--- /dev/null
+++ b/src/relentless/model/potential/angle.py
@@ -0,0 +1,471 @@
+import numpy
+
+from . import potential
+
+
+class AngleParameters(potential.Parameters):
+    """Parameters for a angle potential."""
+
+    pass
+
+
+class AnglePotential(potential.BondedPotential):
+    r"""Abstract base class for an angle potential.
+
+    The angle potential is defined by the angle between three bonded particles.
+    The angle between particles :math:`i`, :math:`j`, and :math:`k` is defined as:
+
+    .. math::
+
+        \theta_{ijk} = \arccos\left(\frac{\mathbf{r}_{ij} \cdot
+        \mathbf{r}_{jk}}{|\mathbf{r}_{ij}||\mathbf{r}_{jk}|}\right)
+
+    where :math:`\mathbf{r}_{ij}` and :math:`\mathbf{r}_{jk}` are the vectors
+    from particle *i* to particle *j* and from particle *j* to particle *k*,
+    respectively.
+
+    """
+
+    def derivative(self, type_, var, theta):
+        r"""Evaluate derivative with respect to a variable.
+
+        The derivative is evaluated using the :meth:`_derivative` function for all
+        :math:`u_{0,\lambda}(theta)`.
+
+        The derivative will be carried out with respect to ``var`` for all
+        :class:`~relentless.variable.Variable` parameters. The appropriate chain
+        rules are handled automatically. If the potential does not depend on
+        ``var``, the derivative will be zero by definition.
+
+        Parameters
+        ----------
+        _type : tuple[str]
+            The type for which to calculate the derivative.
+        var : :class:`~relentless.variable.Variable`
+            The variable with respect to which the derivative is calculated.
+        theta : float or list
+            The bond distance(s) at which to evaluate the derivative.
+
+        Returns
+        -------
+        float or numpy.ndarray
+            The bond derivative evaluated at ``theta``. The return type is consistent
+            with ``theta``.
+
+        Raises
+        ------
+        ValueError
+            If any value in ``theta`` is negative.
+        TypeError
+            If the parameter with respect to which to take the derivative
+            is not a :class:`~relentless.variable.Variable`.
+
+        """
+
+        return super().derivative(type_=type_, var=var, x=theta)
+
+
+class HarmonicAngle(AnglePotential):
+    r"""Harmonic angle potential.
+
+    .. math::
+
+        u(\theta) = \frac{k}{2} (\theta - \theta_0)^2
+
+    where :math:`\theta` is the angle between three bonded particles. The parameters
+    for each type are:
+
+    +-------------+--------------------------------------------------+
+    | Parameter   | Description                                      |
+    +=============+==================================================+
+    | ``k``       | Spring constant :math:`k`.                       |
+    +-------------+--------------------------------------------------+
+    | ``theta0``  | Minimum-energy angle :math:`\theta_0`.           |
+    +-------------+--------------------------------------------------+
+
+    Parameters
+    ----------
+    types : tuple[str]
+        Types.
+    name : str
+        Unique name of the potential. Defaults to ``__u[id]``, where ``id`` is the
+        unique integer ID of the potential.
+
+    Attributes
+    ----------
+    coeff : :class:`AngleParameters`
+        Parameters of the potential for each type.
+
+    Examples
+    --------
+    Harmonic Angle::
+
+        >>> u = relentless.potential.angle.HarmonicAngle(("A",))
+        >>> u.coeff["A"].update({'k': 1000, 'theta0': 1.9})
+
+    """
+
+    def __init__(self, types, name=None):
+        super().__init__(keys=types, params=("k", "theta0"), name=name)
+
+    def energy(self, type_, theta):
+        """Evaluate angle energy."""
+        params = self.coeff.evaluate(type_)
+        k = params["k"]
+        theta0 = params["theta0"]
+
+        theta, u, s = self._zeros(theta)
+
+        u = 0.5 * k * (theta - theta0) ** 2
+
+        if s:
+            u = u.item()
+        return u
+
+    def force(self, type_, theta):
+        """Evaluate angle force."""
+        params = self.coeff.evaluate(type_)
+        k = params["k"]
+        theta0 = params["theta0"]
+
+        theta, f, s = self._zeros(theta)
+
+        f = -k * (theta - theta0)
+
+        if s:
+            f = f.item()
+        return f
+
+    def _derivative(self, param, theta, k, theta0):
+        r"""Evaluate angle derivative with respect to a variable."""
+        theta, d, s = self._zeros(theta)
+
+        if param == "k":
+            d = (theta - theta0) ** 2 / 2
+        elif param == "theta0":
+            d = -k * (theta - theta0)
+        else:
+            raise ValueError("Unknown parameter")
+
+        if s:
+            d = d.item()
+        return d
+
+
+class HarmonicCosineAngle(AnglePotential):
+    r"""Harmonic cosine angle potential.
+
+    .. math::
+
+        u(\theta) = \frac{k}{2} (\cos \theta - \cos \theta_0)^2
+
+    where :math:`\theta` is the angle between three bonded particles. The parameters
+    for each type are:
+
+    +-------------+--------------------------------------------------+
+    | Parameter   | Description                                      |
+    +=============+==================================================+
+    | ``k``       | Spring constant :math:`k`.                       |
+    +-------------+--------------------------------------------------+
+    | ``theta0``  | Minimum-energy angle :math:`\theta_0`.           |
+    +-------------+--------------------------------------------------+
+
+    Parameters
+    ----------
+    types : tuple[str]
+        Types.
+    name : str
+        Unique name of the potential. Defaults to ``__u[id]``, where ``id`` is the
+        unique integer ID of the potential.
+
+    Attributes
+    ----------
+    coeff : :class:`AngleParameters`
+        Parameters of the potential for each type.
+
+    Examples
+    --------
+    Cosine Squared Angle::
+
+        >>> u = relentless.potential.angle.HarmonicCosineAngle(("A",))
+        >>> u.coeff["A"].update({'k': 1000, 'theta0': 1})
+
+    """
+
+    def __init__(self, types, name=None):
+        super().__init__(keys=types, params=("k", "theta0"), name=name)
+
+    def energy(self, type_, theta):
+        """Evaluate angle energy."""
+        params = self.coeff.evaluate(type_)
+        k = params["k"]
+        theta0 = params["theta0"]
+
+        theta, u, s = self._zeros(theta)
+
+        u = 0.5 * k * (numpy.cos(theta) - numpy.cos(theta0)) ** 2
+
+        if s:
+            u = u.item()
+        return u
+
+    def force(self, type_, theta):
+        """Evaluate angle force."""
+        params = self.coeff.evaluate(type_)
+        k = params["k"]
+        theta0 = params["theta0"]
+
+        theta, f, s = self._zeros(theta)
+
+        f = k * (numpy.cos(theta) - numpy.cos(theta0)) * numpy.sin(theta)
+
+        if s:
+            f = f.item()
+        return f
+
+    def _derivative(self, param, theta, k, theta0):
+        r"""Evaluate angle derivative with respect to a variable."""
+        theta, d, s = self._zeros(theta)
+
+        if param == "k":
+            d = 0.5 * (numpy.cos(theta) - numpy.cos(theta0)) ** 2
+        elif param == "theta0":
+            d = k * (numpy.cos(theta) - numpy.cos(theta0)) * numpy.sin(theta0)
+        else:
+            raise ValueError("Unknown parameter")
+
+        if s:
+            d = d.item()
+        return d
+
+
+class CosineAngle(AnglePotential):
+    r"""Cosine angle potential.
+
+    .. math::
+
+        u(\theta) = k (1 + \cos \theta)
+
+    where :math:`\theta` is the angle between three bonded particles. The parameters
+    for each type are:
+
+    +-------------+--------------------------------------------------+
+    | Parameter   | Description                                      |
+    +=============+==================================================+
+    | ``k``       | Spring constant :math:`k`.                       |
+    +-------------+--------------------------------------------------+
+
+    Parameters
+    ----------
+    types : tuple[str]
+        Types.
+    name : str
+        Unique name of the potential. Defaults to ``__u[id]``, where ``id`` is the
+        unique integer ID of the potential.
+
+    Attributes
+    ----------
+    coeff : :class:`AngleParameters`
+        Parameters of the potential for each type.
+
+    Examples
+    --------
+    Cosine Angle::
+
+        >>> u = relentless.potential.angle.CosineAngle(("A",))
+        >>> u.coeff["A"].update({'k': 1000})
+
+    """
+
+    def __init__(self, types, name=None):
+        super().__init__(keys=types, params=("k",), name=name)
+
+    def energy(self, type_, theta):
+        """Evaluate angle energy."""
+        params = self.coeff.evaluate(type_)
+        k = params["k"]
+
+        theta, u, s = self._zeros(theta)
+
+        u = k * (1 + numpy.cos(theta))
+
+        if s:
+            u = u.item()
+        return u
+
+    def force(self, type_, theta):
+        """Evaluate angle force."""
+        params = self.coeff.evaluate(type_)
+        k = params["k"]
+
+        theta, f, s = self._zeros(theta)
+
+        f = k * numpy.sin(theta)
+
+        if s:
+            f = f.item()
+        return f
+
+    def _derivative(self, param, theta, k):
+        r"""Evaluate angle derivative with respect to a variable."""
+        theta, d, s = self._zeros(theta)
+
+        if param == "k":
+            d = 1 + numpy.cos(theta)
+        else:
+            raise ValueError("Unknown parameter")
+
+        if s:
+            d = d.item()
+        return d
+
+
+class AngleSpline(potential.BondedSpline, AnglePotential):
+    """Spline angle potentials.
+
+    The angle spline potential is defined by interpolation through a set of
+    knot points. The interpolation scheme uses Akima splines.
+
+    Parameters
+    ----------
+    types : tuple[str]
+        Types.
+    num_knots : int
+        Number of knots.
+    mode : str
+        Mode for storing the values of the knots in
+        :class:`~relentless.variable.Variable` that can be optimized. If
+        ``mode='value'``, the knot amplitudes are manipulated directly.
+        If ``mode='diff'``, the amplitude of the *last* knot is fixed, and
+        differences between neighboring knots are manipulated for all other
+        knots. Defaults to ``'diff'``.
+    name : str
+        Unique name of the potential. Defaults to ``__u[id]``, where ``id`` is the
+        unique integer ID of the potential.
+
+    Examples
+    --------
+    The spline potential is setup from a tabulated potential instead of
+    specifying knot parameters directly::
+
+        spline = relentless.potential.angle.AngleSpline(types=("angleA",), num_knots=3)
+        spline.from_array(("angleA"),[0,1,2],[4,2,0])
+
+    However, the knot variables can be iterated over and manipulated directly::
+
+        for r,k in spline.knots("angleA"):
+            k.value = 1.0
+
+    """
+
+    _space_coord_name = "theta"
+
+    def __init__(self, types, num_knots, mode="diff", name=None):
+        super().__init__(types=types, num_knots=num_knots, mode=mode, name=name)
+
+    def from_array(self, types, theta, u):
+        r"""Set up the potential from knot points.
+
+        Parameters
+        ----------
+        types : tuple[str]
+            The type for which to set up the potential.
+        theta : list
+            Position of each knot.
+        u : list
+            Potential energy of each knot.
+
+        Raises
+        ------
+        ValueError
+            If the number of ``theta`` values is not the same as the number of knots.
+        ValueError
+            If the number of ``u`` values is not the same as the number of knots.
+
+        """
+
+        if theta[0] != 0.0 and theta[-1] != numpy.pi:
+            raise ValueError("The first and last knot must be at 0 and pi.")
+        return super().from_array(types=types, x=theta, u=u)
+
+    def energy(self, type_, theta):
+        """Evaluate potential energy.
+
+        Parameters
+        ----------
+        type_
+            Type parametrizing the potential in :attr:`coeff<container>`.
+        theta : float or list
+            Potential energy coordinate.
+
+        Returns
+        -------
+        float or numpy.ndarray
+            The pair energy evaluated at ``theta``. The return type is consistent
+            with ``theta``.
+
+        """
+        return super().energy(type_=type_, x=theta)
+
+    def force(self, type_, theta):
+        """Evaluate force magnitude.
+
+        The force is the (negative) magnitude of the ``theta`` gradient.
+
+        Parameters
+        ----------
+        type_
+            Type parametrizing the potential in :attr:`coeff<container>`.
+        theta : float or list
+            Potential energy coordinate.
+
+        Returns
+        -------
+        float or numpy.ndarray
+            The force evaluated at ``theta``. The return type is consistent
+            with ``theta``.
+
+        """
+        return super().force(type_=type_, x=theta)
+
+    def derivative(self, type_, var, theta):
+        r"""Evaluate derivative with respect to a variable.
+
+        The derivative is evaluated using the :meth:`_derivative` function for all
+        :math:`u_{0,\lambda}(theta)`.
+
+        The derivative will be carried out with respect to ``var`` for all
+        :class:`~relentless.variable.Variable` parameters. The appropriate chain
+        rules are handled automatically. If the potential does not depend on
+        ``var``, the derivative will be zero by definition.
+
+        Parameters
+        ----------
+        _type : tuple[str]
+            The type for which to calculate the derivative.
+        var : :class:`~relentless.variable.Variable`
+            The variable with respect to which the derivative is calculated.
+        theta : float or list
+            The bond distance(s) at which to evaluate the derivative.
+
+        Returns
+        -------
+        float or numpy.ndarray
+            The bond derivative evaluated at ``theta``. The return type is consistent
+            with ``theta``.
+
+        Raises
+        ------
+        ValueError
+            If any value in ``theta`` is negative.
+        TypeError
+            If the parameter with respect to which to take the derivative
+            is not a :class:`~relentless.variable.Variable`.
+
+        """
+
+        return super().derivative(type_=type_, var=var, x=theta)
+
+    def _derivative(self, param, theta, **params):
+        """Evaluate angle derivative with respect to a variable."""
+        return super()._derivative(param=param, x=theta, **params)
diff --git a/src/relentless/model/potential/bond.py b/src/relentless/model/potential/bond.py
index 763b4b9d..00529137 100644
--- a/src/relentless/model/potential/bond.py
+++ b/src/relentless/model/potential/bond.py
@@ -1,8 +1,5 @@
 import numpy
 
-from relentless import math
-from relentless.model import variable
-
 from . import potential
 
 
@@ -15,7 +12,43 @@ class BondParameters(potential.Parameters):
 class BondPotential(potential.BondedPotential):
     r"""Abstract base class for a bond potential."""
 
-    pass
+    def derivative(self, type_, var, r):
+        r"""Evaluate bond derivative with respect to a variable.
+
+        The derivative is evaluated using the :meth:`_derivative` function for all
+        :math:`u_{0,\lambda}(r)`.
+
+        The derivative will be carried out with respect to ``var`` for all
+        :class:`~relentless.variable.Variable` parameters. The appropriate chain
+        rules are handled automatically. If the potential does not depend on
+        ``var``, the derivative will be zero by definition.
+
+        Parameters
+        ----------
+        _type : tuple[str]
+            The type for which to calculate the derivative.
+        var : :class:`~relentless.variable.Variable`
+            The variable with respect to which the derivative is calculated.
+        r : float or list
+            The bond distance(s) at which to evaluate the derivative.
+
+        Returns
+        -------
+        float or numpy.ndarray
+            The bond derivative evaluated at ``r``. The return type is consistent
+            with ``r``.
+
+        Raises
+        ------
+        ValueError
+            If any value in ``r`` is negative.
+        TypeError
+            If the parameter with respect to which to take the derivative
+            is not a :class:`~relentless.variable.Variable`.
+
+        """
+
+        return super().derivative(type_=type_, var=var, x=r)
 
 
 class HarmonicBond(BondPotential):
@@ -270,11 +303,11 @@ def _derivative(self, param, r, k, r0, epsilon, sigma):
         return d
 
 
-class BondSpline(BondPotential):
-    """Spline bond potential.
+class BondSpline(potential.BondedSpline, BondPotential):
+    """Spline bond potentials.
 
-    The bond potential is defined by interpolation through a set of knot points.
-    The interpolation scheme uses Akima splines.
+    The bonded spline potential is defined by interpolation through a set of
+    knot points. The interpolation scheme uses Akima splines.
 
     Parameters
     ----------
@@ -308,35 +341,10 @@ class BondSpline(BondPotential):
 
     """
 
-    valid_modes = ("value", "diff")
+    _space_coord_name = "r"
 
     def __init__(self, types, num_knots, mode="diff", name=None):
-        if isinstance(num_knots, int) and num_knots >= 2:
-            self._num_knots = num_knots
-        else:
-            raise ValueError("Number of spline knots must be an integer >= 2.")
-
-        if mode in self.valid_modes:
-            self._mode = mode
-        else:
-            raise ValueError("Invalid mode, choose from: " + ",".join(self.valid_modes))
-
-        params = []
-        for i in range(self.num_knots):
-            ri, ki = self.knot_params(i)
-            params.append(ri)
-            params.append(ki)
-        super().__init__(keys=types, params=params, name=name)
-
-    @classmethod
-    def from_json(cls, data, name=None):
-        u = super().from_json(data, name)
-        # reset the knot values as variables since they were set as floats
-        for type in u.coeff:
-            for i, (r, k) in enumerate(u.knots(type)):
-                u._set_knot(type, i, r, k)
-
-        return u
+        super().__init__(types=types, num_knots=num_knots, mode=mode, name=name)
 
     def from_array(self, types, r, u):
         r"""Set up the potential from knot points.
@@ -358,203 +366,45 @@ def from_array(self, types, r, u):
             If the number of ``u`` values is not the same as the number of knots.
 
         """
-        # check that r and u have the right shape
-        if len(r) != self.num_knots:
-            raise ValueError("r must have the same length as the number of knots")
-        if len(u) != self.num_knots:
-            raise ValueError("u must have the same length as the number of knots")
-
-        # convert to r,knot form given the mode
-        rs = numpy.asarray(r, dtype=float)
-        ks = numpy.asarray(u, dtype=float)
-        if self.mode == "diff":
-            # difference is next knot minus m y knot,
-            # with last knot fixed at its current value
-            ks[:-1] -= ks[1:]
-
-        # convert knot positions to differences
-        drs = numpy.zeros_like(rs)
-        drs[0] = rs[0]
-        drs[1:] = rs[1:] - rs[:-1]
-
-        for i in range(self.num_knots):
-            self._set_knot(types, i, drs[i], ks[i])
-
-    def to_json(self):
-        data = super().to_json()
-        data["num_knots"] = self.num_knots
-        data["mode"] = self.mode
-        return data
-
-    def knot_params(self, i):
-        r"""Get the parameter names for a given knot.
-
-        Parameters
-        ----------
-        i : int
-            Key for the knot variable.
 
-        Returns
-        -------
-        str
-            The parameter name of the :math:`r` value.
-        str
-            The parameter name of the knot value.
-
-        Raises
-        ------
-        TypeError
-            If the knot key is not an integer.
-
-        """
-        if not isinstance(i, int):
-            raise TypeError("Knots are keyed by integers")
-        return f"dr-{i}", f"{self.mode}-{i}"
-
-    def _set_knot(self, types, i, dr, k):
-        """Set the value of knot variables.
-
-        The meaning of the value of the knot variable is defined by the ``mode``.
-        This method is mostly meant to coerce the knot variable types.
-
-        Parameters
-        ----------
-        types : tuple[str]
-            The type for which to set up the potential.
-        i : int
-            Index of the knot.
-        r : float
-            Relative position of each knot from previous one.
-        u : float
-            Value of the knot variable.
-
-        """
-        if i > 0 and dr <= 0:
-            raise ValueError("Knot spacings must be positive")
-
-        dri, ki = self.knot_params(i)
-        if isinstance(self.coeff[types][dri], variable.IndependentVariable):
-            self.coeff[types][dri].value = dr
-        else:
-            self.coeff[types][dri] = variable.IndependentVariable(dr)
-
-        if isinstance(self.coeff[types][ki], variable.IndependentVariable):
-            self.coeff[types][ki].value = k
-        else:
-            self.coeff[types][ki] = variable.IndependentVariable(k)
+        return super().from_array(types=types, x=r, u=u)
 
     def energy(self, type_, r):
-        params = self.coeff.evaluate(type_)
-        r, u, s = self._zeros(r)
-        u = self._interpolate(params)(r)
-        if s:
-            u = u.item()
-        return u
-
-    def force(self, type_, r):
-        params = self.coeff.evaluate(type_)
-        r, f, s = self._zeros(r)
-        f = -self._interpolate(params).derivative(r, 1)
-        if s:
-            f = f.item()
-        return f
-
-    def _derivative(self, param, r, **params):
-        r, d, s = self._zeros(r)
-        h = 0.001
-
-        if "dr-" in param:
-            f_low = self._interpolate(params)(r)
-            knot_p = params[param]
-            params[param] = knot_p + h
-            f_high = self._interpolate(params)(r)
-            params[param] = knot_p
-            d = (f_high - f_low) / h
-        elif self.mode + "-" in param:
-            # perturb knot param value
-            knot_p = params[param]
-            params[param] = knot_p + h
-            f_high = self._interpolate(params)(r)
-            params[param] = knot_p - h
-            f_low = self._interpolate(params)(r)
-            params[param] = knot_p
-            d = (f_high - f_low) / (2 * h)
-        else:
-            raise ValueError("Parameter cannot be differentiated")
-
-        if s:
-            d = d.item()
-        return d
-
-    def _interpolate(self, params):
-        """Interpolate the knot points into a spline potential.
+        """Evaluate potential energy.
 
         Parameters
         ----------
-        params : dict
-            The knot parameters
+        type_
+            Type parametrizing the potential in :attr:`coeff<container>`.
+        r : float or list
+            Potential energy coordinate.
 
         Returns
         -------
-        :class:`~relentless.math.Interpolator`
-            The interpolated spline potential.
+        float or numpy.ndarray
+            The pair energy evaluated at ``r``. The return type is consistent
+            with ``r``.
 
         """
-        dr = numpy.zeros(self.num_knots)
-        u = numpy.zeros(self.num_knots)
-        for i in range(self.num_knots):
-            dri, ki = self.knot_params(i)
-            dr[i] = params[dri]
-            u[i] = params[ki]
-
-        if numpy.any(dr[1:] <= 0):
-            raise ValueError("Knot differences must be positive")
+        return super().energy(type_=type_, x=r)
 
-        # reconstruct r from differences
-        r = numpy.cumsum(dr)
-
-        # reconstruct the energies from differences, starting from the end
-        if self.mode == "diff":
-            u = numpy.flip(numpy.cumsum(numpy.flip(u)))
-
-        return math.AkimaSpline(x=r, y=u)
-
-    @property
-    def num_knots(self):
-        """int: Number of knots."""
-        return self._num_knots
-
-    @property
-    def mode(self):
-        """str: Spline construction mode."""
-        return self._mode
+    def force(self, type_, r):
+        """Evaluate force magnitude.
 
-    def knots(self, types):
-        r"""Generator for knot points.
+        The force is the (negative) magnitude of the ``x`` gradient.
 
         Parameters
         ----------
-        types : tuple[str]
-            The types for which to retrieve the knot points.
+        type_
+            Type parametrizing the potential in :attr:`coeff<container>`.
+        x : float or list
+            Potential energy coordinate.
 
-        Yields
-        ------
-        :class:`~relentless.variable.Variable`
-            The next :math:`dr` variable in the parameters.
-        :class:`~relentless.variable.Variable`
-            The next knot variable in the parameters.
+        Returns
+        -------
+        float or numpy.ndarray
+            The force evaluated at ``r``. The return type is consistent
+            with ``r``.
 
         """
-        for i in range(self.num_knots):
-            dri, ki = self.knot_params(i)
-            yield self.coeff[types][dri], self.coeff[types][ki]
-
-    @property
-    def design_variables(self):
-        """tuple: Designable variables of the spline."""
-        dvars = []
-        for types in self.coeff:
-            for i, (dr, k) in enumerate(self.knots(types)):
-                if i != self.num_knots - 1:
-                    dvars.append(k)
-        return tuple(dvars)
+        return super().force(type_=type_, x=r)
diff --git a/src/relentless/model/potential/potential.py b/src/relentless/model/potential/potential.py
index 1cc7b9e9..0b7691c6 100644
--- a/src/relentless/model/potential/potential.py
+++ b/src/relentless/model/potential/potential.py
@@ -5,7 +5,7 @@
 
 import numpy
 
-from relentless import collections, mpi
+from relentless import collections, math, mpi
 from relentless.model import variable
 
 
@@ -389,11 +389,11 @@ class BondedPotential(Potential):
     parameter :math:`\lambda`.
     """
 
-    def derivative(self, type_, var, r):
+    def derivative(self, type_, var, x):
         r"""Evaluate bond derivative with respect to a variable.
 
         The derivative is evaluated using the :meth:`_derivative` function for all
-        :math:`u_{0,\lambda}(r)`.
+        :math:`u_{0,\lambda}(x)`.
 
         The derivative will be carried out with respect to ``var`` for all
         :class:`~relentless.variable.Variable` parameters. The appropriate chain
@@ -406,37 +406,35 @@ def derivative(self, type_, var, r):
             The type for which to calculate the derivative.
         var : :class:`~relentless.variable.Variable`
             The variable with respect to which the derivative is calculated.
-        r : float or list
+        x : float or list
             The bond distance(s) at which to evaluate the derivative.
 
         Returns
         -------
         float or numpy.ndarray
-            The bond derivative evaluated at ``r``. The return type is consistent
-            with ``r``.
+            The bond derivative evaluated at ``x``. The return type is consistent
+            with ``x``.
 
         Raises
         ------
         ValueError
-            If any value in ``r`` is negative.
+            If any value in ``x`` is negative.
         TypeError
             If the parameter with respect to which to take the derivative
             is not a :class:`~relentless.variable.Variable`.
-        ValueError
-            If the potential is shifted without setting ``rmax``.
 
         """
         params = self.coeff.evaluate(type_)
-        r, deriv, scalar_r = self._zeros(r)
-        if any(r < 0):
-            raise ValueError("r cannot be negative")
+        x, deriv, scalar_x = self._zeros(x)
+        if any(x < 0):
+            raise ValueError("x cannot be negative")
         if not isinstance(var, variable.Variable):
             raise TypeError(
                 "Parameter with respect to which to take the derivative"
                 " must be a Variable."
             )
 
-        flags = numpy.ones(r.shape[0], dtype=bool)
+        flags = numpy.ones(x.shape[0], dtype=bool)
 
         for p in self.coeff.params:
             # try to take chain rule w.r.t. variable first
@@ -453,19 +451,19 @@ def derivative(self, type_, var, r):
                 continue
 
             # now take the parameter derivative
-            below = numpy.zeros(r.shape[0], dtype=bool)
-            above = numpy.zeros(r.shape[0], dtype=bool)
+            below = numpy.zeros(x.shape[0], dtype=bool)
+            above = numpy.zeros(x.shape[0], dtype=bool)
 
             flags = numpy.logical_and(~below, ~above)
-            deriv[flags] += self._derivative(p, r[flags], **params) * dp_dvar
+            deriv[flags] += self._derivative(p, x[flags], **params) * dp_dvar
 
         # coerce derivative back into shape of the input
-        if scalar_r:
+        if scalar_x:
             deriv = deriv.item()
         return deriv
 
     @abc.abstractmethod
-    def _derivative(self, param, r, **params):
+    def _derivative(self, param, x, **params):
         """Implementation of the parameter derivative function.
 
         This abstract method defines the interface for computing the parameter
@@ -477,7 +475,7 @@ def _derivative(self, param, r, **params):
         ----------
         param : str
             Name of the parameter.
-        r : float or list
+        x : float or list
             The bond distance(s) at which to evaluate the derivative.
         **params : kwargs
             Named parameters of the potential.
@@ -485,8 +483,324 @@ def _derivative(self, param, r, **params):
         Returns
         -------
         float or numpy.ndarray
-            The bond derivative evaluated at ``r``. The return type is consistent
-            with ``r``.
+            The bond derivative evaluated at ``x``. The return type is consistent
+            with ``x``.
 
         """
         pass
+
+
+class BondedSpline(BondedPotential):
+    """Base class for bonded spline potentials.
+
+    The bonded spline potential is defined by interpolation through a set of
+    knot points. The interpolation scheme uses Akima splines.
+
+    This class should not be instantiated directly. Instead, use the appropriate
+    spline type, i.e., :class:`~relentless.model.potential.BondSpline` or
+    :class:`~relentless.model.potential.AngleSpline`.
+
+    Parameters
+    ----------
+    types : tuple[str]
+        Types.
+    num_knots : int
+        Number of knots.
+    mode : str
+        Mode for storing the values of the knots in
+        :class:`~relentless.variable.Variable` that can be optimized. If
+        ``mode='value'``, the knot amplitudes are manipulated directly.
+        If ``mode='diff'``, the amplitude of the *last* knot is fixed, and
+        differences between neighboring knots are manipulated for all other
+        knots. Defaults to ``'diff'``.
+    name : str
+        Unique name of the potential. Defaults to ``__u[id]``, where ``id`` is the
+        unique integer ID of the potential.
+
+    """
+
+    valid_modes = ("value", "diff")
+    _space_coord_name = "x"
+
+    def __init__(self, types, num_knots, mode="diff", name=None):
+        if isinstance(num_knots, int) and num_knots >= 2:
+            self._num_knots = num_knots
+        else:
+            raise ValueError("Number of spline knots must be an integer >= 2.")
+
+        if mode in self.valid_modes:
+            self._mode = mode
+        else:
+            raise ValueError("Invalid mode, choose from: " + ",".join(self.valid_modes))
+
+        params = []
+        for i in range(self.num_knots):
+            xi, ki = self.knot_params(i)
+            params.append(xi)
+            params.append(ki)
+        super().__init__(keys=types, params=params, name=name)
+
+    @classmethod
+    def from_json(cls, data, name=None):
+        u = super().from_json(data, name)
+        # reset the knot values as variables since they were set as floats
+        for type in u.coeff:
+            for i, (x, k) in enumerate(u.knots(type)):
+                u._set_knot(type, i, x, k)
+
+        return u
+
+    def from_array(self, types, x, u):
+        r"""Set up the potential from knot points.
+
+        Parameters
+        ----------
+        types : tuple[str]
+            The type for which to set up the potential.
+        x : list
+            Position of each knot.
+        u : list
+            Potential energy of each knot.
+
+        Raises
+        ------
+        ValueError
+            If the number of ``x`` values is not the same as the number of knots.
+        ValueError
+            If the number of ``u`` values is not the same as the number of knots.
+
+        """
+        # check that r and u have the right shape
+        if len(x) != self.num_knots:
+            raise ValueError("x must have the same length as the number of knots")
+        if len(u) != self.num_knots:
+            raise ValueError("u must have the same length as the number of knots")
+
+        # convert to r,knot form given the mode
+        xs = numpy.asarray(x, dtype=float)
+        ks = numpy.asarray(u, dtype=float)
+        if self.mode == "diff":
+            # difference is next knot minus m y knot,
+            # with last knot fixed at its current value
+            ks[:-1] -= ks[1:]
+
+        # convert knot positions to differences
+        dxs = numpy.zeros_like(xs)
+        dxs[0] = xs[0]
+        dxs[1:] = xs[1:] - xs[:-1]
+
+        for i in range(self.num_knots):
+            self._set_knot(types, i, dxs[i], ks[i])
+
+    def to_json(self):
+        data = super().to_json()
+        data["num_knots"] = self.num_knots
+        data["mode"] = self.mode
+        return data
+
+    def knot_params(self, i):
+        r"""Get the parameter names for a given knot.
+
+        Parameters
+        ----------
+        i : int
+            Key for the knot variable.
+
+        Returns
+        -------
+        str
+            The parameter name of the :math:`x` value.
+        str
+            The parameter name of the knot value.
+
+        Raises
+        ------
+        TypeError
+            If the knot key is not an integer.
+
+        """
+        if not isinstance(i, int):
+            raise TypeError("Knots are keyed by integers")
+        return f"d{self._space_coord_name}-{i}", f"{self.mode}-{i}"
+
+    def _set_knot(self, types, i, dx, k):
+        """Set the value of knot variables.
+
+        The meaning of the value of the knot variable is defined by the ``mode``.
+        This method is mostly meant to coerce the knot variable types.
+
+        Parameters
+        ----------
+        types : tuple[str]
+            The type for which to set up the potential.
+        i : int
+            Index of the knot.
+        x : float
+            Relative position of each knot from previous one.
+        u : float
+            Value of the knot variable.
+
+        """
+        if i > 0 and dx <= 0:
+            raise ValueError("Knot spacings must be positive")
+
+        dxi, ki = self.knot_params(i)
+        if isinstance(self.coeff[types][dxi], variable.IndependentVariable):
+            self.coeff[types][dxi].value = dx
+        else:
+            self.coeff[types][dxi] = variable.IndependentVariable(dx)
+
+        if isinstance(self.coeff[types][ki], variable.IndependentVariable):
+            self.coeff[types][ki].value = k
+        else:
+            self.coeff[types][ki] = variable.IndependentVariable(k)
+
+    def energy(self, type_, x):
+        """Evaluate potential energy.
+
+        Parameters
+        ----------
+        type_
+            Type parametrizing the potential in :attr:`coeff<container>`.
+        x : float or list
+            Potential energy coordinate.
+
+        Returns
+        -------
+        float or numpy.ndarray
+            The pair energy evaluated at ``x``. The return type is consistent
+            with ``x``.
+
+        """
+        params = self.coeff.evaluate(type_)
+        x, u, s = self._zeros(x)
+        u = self._interpolate(params)(x)
+        if s:
+            u = u.item()
+        return u
+
+    def force(self, type_, x):
+        """Evaluate force magnitude.
+
+        The force is the (negative) magnitude of the ``x`` gradient.
+
+        Parameters
+        ----------
+        type_
+            Type parametrizing the potential in :attr:`coeff<container>`.
+        x : float or list
+            Potential energy coordinate.
+
+        Returns
+        -------
+        float or numpy.ndarray
+            The force evaluated at ``x``. The return type is consistent
+            with ``x``.
+
+        """
+        params = self.coeff.evaluate(type_)
+        x, f, s = self._zeros(x)
+        f = -self._interpolate(params).derivative(x, 1)
+        if s:
+            f = f.item()
+        return f
+
+    def _derivative(self, param, x, **params):
+        x, d, s = self._zeros(x)
+        h = 0.001
+
+        if f"d{self._space_coord_name}-" in param:
+            f_low = self._interpolate(params)(x)
+            knot_p = params[param]
+            params[param] = knot_p + h
+            f_high = self._interpolate(params)(x)
+            params[param] = knot_p
+            d = (f_high - f_low) / h
+        elif self.mode + "-" in param:
+            # perturb knot param value
+            knot_p = params[param]
+            params[param] = knot_p + h
+            f_high = self._interpolate(params)(x)
+            params[param] = knot_p - h
+            f_low = self._interpolate(params)(x)
+            params[param] = knot_p
+            d = (f_high - f_low) / (2 * h)
+        else:
+            raise ValueError("Parameter cannot be differentiated")
+
+        if s:
+            d = d.item()
+        return d
+
+    def _interpolate(self, params):
+        """Interpolate the knot points into a spline potential.
+
+        Parameters
+        ----------
+        params : dict
+            The knot parameters
+
+        Returns
+        -------
+        :class:`~relentless.math.Interpolator`
+            The interpolated spline potential.
+
+        """
+        dx = numpy.zeros(self.num_knots)
+        u = numpy.zeros(self.num_knots)
+        for i in range(self.num_knots):
+            dxi, ki = self.knot_params(i)
+            dx[i] = params[dxi]
+            u[i] = params[ki]
+
+        if numpy.any(dx[1:] <= 0):
+            raise ValueError("Knot differences must be positive")
+
+        # reconstruct r from differences
+        x = numpy.cumsum(dx)
+
+        # reconstruct the energies from differences, starting from the end
+        if self.mode == "diff":
+            u = numpy.flip(numpy.cumsum(numpy.flip(u)))
+
+        return math.AkimaSpline(x=x, y=u)
+
+    @property
+    def num_knots(self):
+        """int: Number of knots."""
+        return self._num_knots
+
+    @property
+    def mode(self):
+        """str: Spline construction mode."""
+        return self._mode
+
+    def knots(self, types):
+        r"""Generator for knot points.
+
+        Parameters
+        ----------
+        types : tuple[str]
+            The types for which to retrieve the knot points.
+
+        Yields
+        ------
+        :class:`~relentless.variable.Variable`
+            The next :math:`dr` variable in the parameters.
+        :class:`~relentless.variable.Variable`
+            The next knot variable in the parameters.
+
+        """
+        for i in range(self.num_knots):
+            dxi, ki = self.knot_params(i)
+            yield self.coeff[types][dxi], self.coeff[types][ki]
+
+    @property
+    def design_variables(self):
+        """tuple: Designable variables of the spline."""
+        dvars = []
+        for types in self.coeff:
+            for i, (dx, k) in enumerate(self.knots(types)):
+                if i != self.num_knots - 1:
+                    dvars.append(k)
+        return tuple(dvars)
diff --git a/src/relentless/simulate/__init__.py b/src/relentless/simulate/__init__.py
index b89b4f50..584ba405 100644
--- a/src/relentless/simulate/__init__.py
+++ b/src/relentless/simulate/__init__.py
@@ -122,6 +122,7 @@
 
     Potentials
     PotentialTabulator
+    AnglePotentialTabulator
     BondPotentialTabulator
     PairPotentialTabulator
 
@@ -158,6 +159,7 @@
 )
 from .simulate import (
     AnalysisOperation,
+    AnglePotentialTabulator,
     BondPotentialTabulator,
     PairPotentialTabulator,
     Potentials,
diff --git a/src/relentless/simulate/hoomd.py b/src/relentless/simulate/hoomd.py
index 018fbf12..09c9bd21 100644
--- a/src/relentless/simulate/hoomd.py
+++ b/src/relentless/simulate/hoomd.py
@@ -104,6 +104,21 @@ def _call_v3(self, sim):
                 bond_potential.params[i] = dict(r_min=r[0], r_max=r[-1], U=u[:], F=f[:])
             sim[self]["_potentials"].append(bond_potential)
 
+        sim[self]["_angles"] = self._get_angles_from_snapshot(sim, snap)
+        if sim.potentials.angle is not None:
+            sim.angle_types = sim["engine"]["_hoomd"].state.angle_types
+            angle_potential = hoomd.md.angle.Table(width=sim.potentials.angle.num)
+
+            for i in sim.angle_types:
+                u, tau = (
+                    sim.potentials.angle.energy(key=i),
+                    sim.potentials.angle.force(key=i),
+                )
+                if numpy.any(numpy.isinf(u)):
+                    raise ValueError("Angle potential/force is infinite at evaluated r")
+                angle_potential.params[i] = dict(U=u[:], tau=tau[:])
+            sim[self]["_potentials"].append(angle_potential)
+
     def _call_v2(self, sim):
         # initialize
         sim[self]["_system"] = self._initialize_v2(sim)
@@ -180,6 +195,12 @@ def _get_bonds_from_snapshot(self, sim, snap):
         else:
             return None
 
+    def _get_angles_from_snapshot(self, sim, snap):
+        if mpi.world.rank_is_root:
+            return snap.angles.group
+        else:
+            return None
+
     def _assert_dimension_safe(self, sim, snap):
         if sim.dimension == 3:
             dim_safe = True
@@ -1017,6 +1038,7 @@ def _pre_run_v3(self, sim, sim_op):
             constraints=self._get_constrained_quantities(sim, sim_op),
             exclusions=sim.potentials.pair.exclusions,
             bonds=sim[sim.initializer]["_bonds"],
+            angles=sim[sim.initializer]["_angles"],
         )
         sim[self]["_hoomd_thermo_callback"] = hoomd.write.CustomWriter(
             trigger=self.every,
@@ -1210,6 +1232,7 @@ def __init__(
             constraints=None,
             exclusions=None,
             bonds=None,
+            angles=None,
         ):
             if dimension not in (2, 3):
                 raise ValueError("Only 2 or 3 dimensions supported")
@@ -1222,6 +1245,7 @@ def __init__(
             self.constraints = constraints if constraints is not None else {}
             self.exclusion = exclusions
             self.bonds = bonds
+            self.angles = angles
 
             # this method handles all the initialization
             self.reset()
@@ -1333,6 +1357,21 @@ def act(self, timestep):
                             )
 
                             neighbors.filter(bond_exclusion_filter)
+                        # Similarly filter angles from the neighbor list
+                        if (
+                            self.rdf_params["exclude"]
+                            and snap.angles.N != 0
+                            and len(neighbors[:]) > 0
+                        ):
+                            angles = numpy.vstack(
+                                [self.angles, numpy.flip(self.angles, axis=1)],
+                            )
+
+                            angle_exclusion_filter = EnsembleAverage._cantor_pairing(
+                                self, angles, neighbors
+                            )
+
+                            neighbors.filter(angle_exclusion_filter)
 
                         for i in self.types:
                             for j in self.types:
diff --git a/src/relentless/simulate/lammps.py b/src/relentless/simulate/lammps.py
index 8253758f..ac9b87e1 100644
--- a/src/relentless/simulate/lammps.py
+++ b/src/relentless/simulate/lammps.py
@@ -186,12 +186,15 @@ def _call_commands(self, sim):
         has_dihedrals = None
         has_impropers = None
         bond_type_label = None
+        angle_type_label = None
         if mpi.world.rank_is_root:
             snap = lammpsio.DataFile(sim[self]["_datafile"]).read()
             has_bonds = snap.has_bonds()
             if has_bonds:
                 bond_type_label = sim[self]["_bonds"].type_label
             has_angles = snap.has_angles()
+            if has_angles:
+                angle_type_label = sim[self]["_angles"].type_label
             has_dihedrals = snap.has_dihedrals()
             has_impropers = snap.has_impropers()
             for i in sim.types:
@@ -221,6 +224,7 @@ def _call_commands(self, sim):
         has_dihedrals = mpi.world.bcast(has_dihedrals)
         has_impropers = mpi.world.bcast(has_impropers)
         bond_type_label = mpi.world.bcast(bond_type_label)
+        angle_type_label = mpi.world.bcast(angle_type_label)
 
         # attach the potentials
         if sim.potentials.pair.start == 0:
@@ -243,6 +247,19 @@ def _call_commands(self, sim):
                     raise ValueError("Bond potential/force is infinite at evaluated r")
             NrB = len(rB)
 
+        if has_angles:
+            uA = collections.FixedKeyDict(angle_type_label.types)
+            fA = collections.FixedKeyDict(angle_type_label.types)
+            for i in angle_type_label.types:
+                thetaA, uA[i], fA[i] = (
+                    sim.potentials.angle.linear_space,
+                    sim.potentials.angle.energy(key=i),
+                    sim.potentials.angle.force(key=i),
+                )
+                if numpy.any(numpy.isinf(uA[i])):
+                    raise ValueError("Angle potential/force is infinite at evaluated r")
+            NthetaA = len(thetaA)
+
         def pair_map(sim, pair):
             # Map lammps type indexes as a pair, lowest type first
             i, j = pair
@@ -307,6 +324,25 @@ def pair_map(sim, pair):
                                 )
                             )
 
+                if has_angles:
+                    fw.write("# LAMMPS tabulated angle potentials\n")
+                    for i in angle_type_label.types:
+                        fw.write("# angle ANGLE_TABLE_{}\n\n".format(i))
+                        fw.write("ANGLE_TABLE_{}\n".format(i))
+                        fw.write(
+                            "N {N} FP {f_low} {f_high} \n\n".format(
+                                N=NthetaA, f_low=fA[i][0], f_high=fA[i][-1]
+                            )
+                        )
+                        for idx, (thetaA_, uA_, fA_) in enumerate(
+                            zip(thetaA * 180 / numpy.pi, uA[i], fA[i]), start=1
+                        ):
+                            fw.write(
+                                "{id} {thetaA} {uA} {fA}\n".format(
+                                    id=idx, thetaA=thetaA_, uA=uA_, fA=fA_
+                                )
+                            )
+
         else:
             file_ = None
         file_ = mpi.world.bcast(file_)
@@ -324,11 +360,16 @@ def pair_map(sim, pair):
             if sim.potentials.pair.exclusions is not None:
                 if "1-2" in sim.potentials.pair.exclusions:
                     excl_12 = 0.0
+                if "1-3" in sim.potentials.pair.exclusions:
+                    excl_13 = 0.0
             cmds += [f"special_bonds lj/coul {excl_12} {excl_13} {excl_14}"]
 
         if has_bonds:
             cmds += ["bond_style table linear {Nb}".format(Nb=NrB)]
 
+        if has_angles:
+            cmds += ["angle_style table linear {Na}".format(Na=NthetaA)]
+
         for i, j in sim.pairs:
             # get lammps type indexes, lowest type first
             id_i, id_j = pair_map(sim, (i, j))
@@ -347,6 +388,16 @@ def pair_map(sim, pair):
                         id=bond_type_label.__getitem__(id),
                     )
                 ]
+
+        if has_angles:
+            for id in angle_type_label.typeid:
+                cmds += [
+                    ("angle_coeff {typeid} {filename}" " ANGLE_TABLE_{id}").format(
+                        typeid=id,
+                        filename=file_,
+                        id=angle_type_label.__getitem__(id),
+                    )
+                ]
         return cmds
 
     @abc.abstractmethod
@@ -393,8 +444,12 @@ def _convert_to_data_file(self, sim):
                         "Atomic atom style is not compatible with topology data."
                     )
                 # store topology data
+                sim[self]["_bonds"] = None
                 if snap.has_bonds():
                     sim[self]["_bonds"] = snap.bonds
+                sim[self]["_angles"] = None
+                if snap.has_angles():
+                    sim[self]["_angles"] = snap.angles
                 # figure out dimensions
                 dimension = self.dimension
                 if dimension is None:
@@ -414,6 +469,7 @@ def _convert_to_data_file(self, sim):
                 )
             else:
                 sim[self]["_bonds"] = None
+                sim[self]["_angles"] = None
                 data_filename = None
                 dimension = None
                 type_map = None
@@ -433,6 +489,7 @@ def _convert_to_data_file(self, sim):
                     typeids = numpy.unique(snap.typeid)
                     type_map = {str(typeid): typeid for typeid in typeids}
                     # store topology data
+                    sim[self]["_bonds"] = None
                     if snap.has_bonds():
                         sim[self]["_bonds"] = snap.bonds
                         unique_bond_types = set()
@@ -446,6 +503,20 @@ def _convert_to_data_file(self, sim):
                         sim[self]["_bonds"].type_label = lammpsio.LabelMap(
                             bond_type_map
                         )
+                    sim[self]["_angles"] = None
+                    if snap.has_angles():
+                        sim[self]["_angles"] = snap.angles
+                        unique_angle_types = set()
+                        for pot in sim.potentials.angle.potentials:
+                            for i in pot.coeff.types:
+                                unique_angle_types.add(i)
+
+                        angle_type_map = {
+                            i + 1: t for i, t in enumerate(unique_angle_types)
+                        }
+                        sim[self]["_angles"].type_label = lammpsio.LabelMap(
+                            angle_type_map
+                        )
                 else:
                     type_map = None
                 type_map = mpi.world.bcast(type_map)
@@ -1279,7 +1350,7 @@ def process(self, sim, sim_op):
                     # consider both (i,j) and (j,i) permutations
                     if (
                         sim[self]["_rdf_params"]["exclude"]
-                        and sim[sim.initializer]["_bonds"].N != 0
+                        and sim[sim.initializer]["_bonds"] is not None
                         and len(neighbors[:]) > 0
                     ):
                         bonds = numpy.vstack(
@@ -1299,6 +1370,29 @@ def process(self, sim, sim_op):
 
                         neighbors.filter(bond_exclusion_filter)
 
+                    # filter angles from the neighbor list if they are present
+                    if (
+                        sim[self]["_rdf_params"]["exclude"]
+                        and sim[sim.initializer]["_angles"] is not None
+                        and len(neighbors[:]) > 0
+                    ):
+                        angles = numpy.vstack(
+                            [
+                                sim[sim.initializer]["_angles"].members,
+                                numpy.flip(
+                                    sim[sim.initializer]["_angles"].members, axis=1
+                                ),
+                            ],
+                        )
+                        # Zero index angles
+                        angles -= 1
+
+                        angle_exclusion_filter = EnsembleAverage._cantor_pairing(
+                            self, angles, neighbors
+                        )
+
+                        neighbors.filter(angle_exclusion_filter)
+
                     for i in sim.types:
                         _rdf_density[i] += N[i] / box.volume
                         _rdf_num_origins[i] += N[i]
diff --git a/src/relentless/simulate/simulate.py b/src/relentless/simulate/simulate.py
index a28d558c..d48d6172 100644
--- a/src/relentless/simulate/simulate.py
+++ b/src/relentless/simulate/simulate.py
@@ -380,6 +380,7 @@ def __init__(self, kB=1.0):
         self.kB = kB
         self.pair = None
         self.bond = None
+        self.angle = None
 
     @property
     def pair(self):
@@ -403,6 +404,17 @@ def bond(self, val):
             raise TypeError("Bond potential must be tabulated")
         self._bond = val
 
+    @property
+    def angle(self):
+        """:class:`AnglePotentialTabulator`: Angle potentials."""
+        return self._angle
+
+    @angle.setter
+    def angle(self, val):
+        if val is not None and not isinstance(val, AnglePotentialTabulator):
+            raise TypeError("Angle potential must be tabulated")
+        self._angle = val
+
 
 class PotentialTabulator:
     r"""Tabulator for an interaction potential.
@@ -837,3 +849,24 @@ class BondPotentialTabulator(PotentialTabulator):
     """
 
     pass
+
+
+class AnglePotentialTabulator(PotentialTabulator):
+    r"""Tabulate one or more angle potentials.
+
+    Enables evaluation of energy, force, and derivative at different
+    angle values (i.e. :math:`\theta`) on a range :math:`\left[ 0, \pi \right]`.
+
+    Parameters
+    ----------
+    potentials : :class:`~relentless.potential.angle.AnglePotential` or array_like
+        The angle potential(s) to be tabulated. If array_like, all elements must
+        be :class:`~relentless.potential.angle.AnglePotential`\s.
+    num : int
+        The number of points (value of :math:`\theta`) at which to tabulate and
+        evaluate the potential.
+
+    """
+
+    def __init__(self, potentials, num):
+        super().__init__(potentials=potentials, start=0.0, stop=numpy.pi, num=num)
diff --git a/tests/model/potential/test_angle.py b/tests/model/potential/test_angle.py
new file mode 100644
index 00000000..576911b4
--- /dev/null
+++ b/tests/model/potential/test_angle.py
@@ -0,0 +1,561 @@
+"""Unit tests for angle module."""
+
+import tempfile
+import unittest
+
+import numpy
+
+import relentless
+
+
+class LinPotAngle(relentless.model.potential.angle.AnglePotential):
+    """Linear angle potential function"""
+
+    def __init__(self, types, params):
+        super().__init__(types, params)
+
+    def to_json(self):
+        data = super().to_json()
+        data["params"] = self.coeff.params
+        return data
+
+    def energy(self, types, theta):
+        m = self.coeff[types]["m"]
+
+        theta, u, s = self._zeros(theta)
+        u[:] = m * theta
+        if s:
+            u = u.item()
+        return u
+
+    def force(self, types, theta):
+        m = self.coeff[types]["m"]
+
+        theta, f, s = self._zeros(theta)
+        f[:] = -m
+        if s:
+            f = f.item()
+        return f
+
+    def _derivative(self, param, theta, **params):
+        theta, d, s = self._zeros(theta)
+        if param == "m":
+            d[:] = theta
+        if s:
+            d = d.item()
+        return d
+
+
+class TwoVarPot(relentless.model.potential.angle.AnglePotential):
+    """Mock potential function"""
+
+    def __init__(self, types, params):
+        super().__init__(types, params)
+
+    @classmethod
+    def from_json(cls, data):
+        raise NotImplementedError()
+
+    def energy(self, theta, x, y, **params):
+        pass
+
+    def force(self, theta, x, y, **params):
+        pass
+
+    def _derivative(self, param, theta, **params):
+        # not real derivative, just used to test functionality
+        theta, d, s = self._zeros(theta)
+        if param == "x":
+            d[:] = 2 * theta
+        elif param == "y":
+            d[:] = 3 * theta
+        if s:
+            d = d.item()
+        return d
+
+
+class test_AnglePotential(unittest.TestCase):
+    """Unit tests for relentless.model.angle.AnglePotential"""
+
+    def test_init(self):
+        """Test creation from data"""
+        # test creation with only m
+        p = LinPotAngle(types=("1",), params=("m",))
+        p.coeff["1"]["m"] = 3.5
+
+        coeff = relentless.model.potential.angle.AngleParameters(
+            types=("1",), params=("m")
+        )
+        coeff["1"]["m"] = 3.5
+
+        self.assertCountEqual(p.coeff.types, coeff.types)
+        self.assertCountEqual(p.coeff.params, coeff.params)
+        self.assertEqual(p.coeff.evaluate(("1")), coeff.evaluate(("1")))
+
+    def test_energy(self):
+        """Test energy method"""
+        p = LinPotAngle(types=("1",), params=("m",))
+        p.coeff["1"]["m"] = 2.0
+
+        # test with no cutoffs
+        u = p.energy(types=("1"), theta=0.5)
+        self.assertAlmostEqual(u, 1.0)
+        u = p.energy(types=("1"), theta=[0.25, 0.75])
+        numpy.testing.assert_allclose(u, [0.5, 1.5])
+
+    def test_force(self):
+        """Test force method"""
+        p = LinPotAngle(types=("1",), params=("m",))
+        p.coeff["1"]["m"] = 2.0
+
+        # test with no cutoffs
+        f = p.force(types=("1"), theta=0.5)
+        self.assertAlmostEqual(f, -2.0)
+        f = p.force(types=("1"), theta=[0.25, 0.75])
+        numpy.testing.assert_allclose(f, [-2.0, -2.0])
+
+    def test_derivative_values(self):
+        """Test derivative method with different param values"""
+        p = LinPotAngle(types=("1",), params=("m",))
+        x = relentless.model.IndependentVariable(value=2.0)
+        p.coeff["1"]["m"] = x
+
+        # test with no cutoffs
+        d = p.derivative(type_=("1"), var=x, theta=0.5)
+        self.assertAlmostEqual(d, 0.5)
+        d = p.derivative(type_=("1"), var=x, theta=[0.25, 0.75])
+        numpy.testing.assert_allclose(d, [0.25, 0.75])
+
+    def test_derivative_types(self):
+        """Test derivative method with different param types."""
+        q = LinPotAngle(types=("1",), params=("m",))
+        x = relentless.model.IndependentVariable(value=4.0)
+        y = relentless.model.IndependentVariable(value=64.0)
+        z = relentless.model.GeometricMean(x, y)
+        q.coeff["1"]["m"] = z
+
+        # test with respect to dependent variable parameter
+        d = q.derivative(type_=("1"), var=z, theta=2.0)
+        self.assertAlmostEqual(d, 2.0)
+
+        # test with respect to independent variable on which parameter is dependent
+        d = q.derivative(type_=("1"), var=x, theta=1.5)
+        self.assertAlmostEqual(d, 3.0)
+        d = q.derivative(type_=("1"), var=y, theta=4.0)
+        self.assertAlmostEqual(d, 0.5)
+
+        # test invalid derivative w.r.t. scalar
+        a = 2.5
+        q.coeff["1"]["m"] = a
+        with self.assertRaises(TypeError):
+            d = q.derivative(type_=("1"), var=a, theta=2.0)
+
+        # test with respect to independent variable which is
+        # related to a SameAs variable
+        r = TwoVarPot(types=("1",), params=("x", "y"))
+
+        r.coeff["1"]["x"] = x
+        r.coeff["1"]["y"] = relentless.model.variable.SameAs(x)
+        d = r.derivative(type_=("1"), var=x, theta=4.0)
+        self.assertAlmostEqual(d, 20.0)
+
+        r.coeff["1"]["y"] = x
+        r.coeff["1"]["x"] = relentless.model.variable.SameAs(x)
+        d = r.derivative(type_=("1"), var=x, theta=4.0)
+        self.assertAlmostEqual(d, 20.0)
+
+    def test_iteration(self):
+        """Test iteration on typesPotential object"""
+        p = LinPotAngle(types=("1",), params=("m",))
+        for types in p.coeff:
+            p.coeff[types]["m"] = 2.0
+
+        self.assertEqual(dict(p.coeff["1"]), {"m": 2.0})
+
+    def test_json(self):
+        """Test saving to file"""
+        p = LinPotAngle(types=("1",), params=("m"))
+        p.coeff["1"]["m"] = 2.0
+
+        data = p.to_json()
+        self.assertEqual(data["id"], p.id)
+        self.assertEqual(data["name"], p.name)
+
+        p2 = LinPotAngle.from_json(data)
+        self.assertEqual(p2.coeff["1"]["m"], 2.0)
+
+    def test_save(self):
+        """Test saving to file"""
+        temp = tempfile.NamedTemporaryFile()
+        p = LinPotAngle(types=("1",), params=("m"))
+        p.coeff["1"]["m"] = 2.0
+        p.save(temp.name)
+
+        p2 = LinPotAngle.from_file(temp.name)
+        self.assertEqual(p2.coeff["1"]["m"], 2.0)
+
+        temp.close()
+
+
+class test_HarmonicAngle(unittest.TestCase):
+    """Unit tests for relentless.model.potential.HarmonicAngle"""
+
+    def test_init(self):
+        """Test creation from data"""
+        harmonic_angle = relentless.model.potential.HarmonicAngle(types=("1",))
+        coeff = relentless.model.potential.AngleParameters(
+            types=("1",),
+            params=(
+                "k",
+                "theta0",
+            ),
+        )
+        self.assertCountEqual(harmonic_angle.coeff.types, coeff.types)
+        self.assertCountEqual(harmonic_angle.coeff.params, coeff.params)
+
+    def test_energy(self):
+        """Test _energy method"""
+        harmonic_angle = relentless.model.potential.HarmonicAngle(types=("1",))
+        harmonic_angle.coeff["1"].update(k=1000, theta0=1.0)
+        # test scalar r
+        theta_input = 0.5
+        u_actual = 125.0
+        u = harmonic_angle.energy(type_=("1"), theta=theta_input)
+        self.assertAlmostEqual(u, u_actual)
+
+        # test array r
+        theta_input = numpy.array([0.0, 0.5, 1.0])
+        u_actual = numpy.array([500.0, 125.0, 0.0])
+        u = harmonic_angle.energy(type_=("1"), theta=theta_input)
+        numpy.testing.assert_allclose(u, u_actual)
+
+    def test_force(self):
+        """Test _force method"""
+        harmonic_angle = relentless.model.potential.HarmonicAngle(types=("1",))
+        harmonic_angle.coeff["1"].update(k=1000, theta0=1.0)
+
+        # test scalar r
+        theta_input = 0.5
+        f_actual = 500
+        f = harmonic_angle.force(type_=("1"), theta=theta_input)
+        self.assertAlmostEqual(f, f_actual)
+
+        # test array r
+        theta_input = numpy.array([0.0, 0.5, 1.0])
+        f_actual = numpy.array([1000, 500, 0])
+        f = harmonic_angle.force(type_=("1"), theta=theta_input)
+        numpy.testing.assert_allclose(f, f_actual)
+
+    def test_derivative(self):
+        """Test _derivative method"""
+        harmonic_angle = relentless.model.potential.HarmonicAngle(types=("1",))
+
+        # w.r.t. k
+        # test scalar r
+        theta_input = 0.5
+        d_actual = 0.125
+        d = harmonic_angle._derivative(param="k", theta=theta_input, k=1000, theta0=1.0)
+        self.assertAlmostEqual(d, d_actual)
+
+        # test array r
+        theta_input = numpy.array([0, 1, 1.5])
+        d_actual = numpy.array([0.50, 0.0, 0.125])
+        d = harmonic_angle._derivative(param="k", theta=theta_input, k=1000, theta0=1.0)
+        numpy.testing.assert_allclose(d, d_actual)
+
+        # w.r.t. theta0
+        # test scalar r
+        theta_input = 0.5
+        d_actual = 500
+        d = harmonic_angle._derivative(
+            param="theta0", theta=theta_input, k=1000.0, theta0=1.0
+        )
+        self.assertAlmostEqual(d, d_actual)
+
+        # test array r
+        theta_input = numpy.array([0, 1, 1.5])
+        d_actual = numpy.array([1000, 0, -500])
+        d = harmonic_angle._derivative(
+            param="theta0", theta=theta_input, k=1000.0, theta0=1.0
+        )
+        numpy.testing.assert_allclose(d, d_actual)
+
+        # test invalid param
+        with self.assertRaises(ValueError):
+            harmonic_angle._derivative(
+                param="ro", theta=theta_input, k=1000.0, theta0=1.0
+            )
+
+    def test_json(self):
+        harmonic_angle = relentless.model.potential.HarmonicAngle(types=("1",))
+        harmonic_angle.coeff["1"].update(k=1000, theta0=1.0)
+        data = harmonic_angle.to_json()
+
+        harmonic_angle2 = relentless.model.potential.HarmonicAngle.from_json(data)
+        self.assertEqual(harmonic_angle2.coeff["1"]["k"], 1000)
+        self.assertEqual(harmonic_angle2.coeff["1"]["theta0"], 1.0)
+
+
+class test_CosineAngle(unittest.TestCase):
+    """Unit tests for relentless.model.potential.CosineAngle"""
+
+    def test_init(self):
+        """Test creation from data"""
+        cosine_angle = relentless.model.potential.CosineAngle(types=("1",))
+        coeff = relentless.model.potential.AngleParameters(
+            types=("1",),
+            params=("k",),
+        )
+        self.assertCountEqual(cosine_angle.coeff.types, coeff.types)
+        self.assertCountEqual(cosine_angle.coeff.params, coeff.params)
+
+    def test_energy(self):
+        """Test _energy method"""
+        cosine_angle = relentless.model.potential.CosineAngle(types=("1",))
+        cosine_angle.coeff["1"].update(
+            k=1000,
+        )
+        # test scalar r
+        theta_input = numpy.pi / 2
+        u_actual = 1000
+        u = cosine_angle.energy(type_=("1"), theta=theta_input)
+        self.assertAlmostEqual(u, u_actual)
+
+        # test array r
+        theta_input = numpy.array([0.0, numpy.pi / 4, numpy.pi])
+        u_actual = numpy.array([2000.0, 1707.10678119, 0])
+        u = cosine_angle.energy(type_=("1"), theta=theta_input)
+        numpy.testing.assert_allclose(u, u_actual)
+
+    def test_force(self):
+        """Test _force method"""
+        cosine_angle = relentless.model.potential.CosineAngle(types=("1",))
+        cosine_angle.coeff["1"].update(k=1000)
+
+        # test scalar r
+        theta_input = numpy.pi / 2
+        f_actual = 1000
+        f = cosine_angle.force(type_=("1"), theta=theta_input)
+        self.assertAlmostEqual(f, f_actual)
+
+        # test array r
+        theta_input = numpy.array([0.0, numpy.pi / 4, numpy.pi])
+        f_actual = numpy.array([0.0, 707.106781187, 0])
+        f = cosine_angle.force(type_=("1"), theta=theta_input)
+        numpy.testing.assert_allclose(f, f_actual, atol=1e-10)
+
+    def test_derivative(self):
+        """Test _derivative method"""
+        cosine_angle = relentless.model.potential.CosineAngle(types=("1",))
+
+        # w.r.t. k
+        # test scalar r
+        theta_input = numpy.pi
+        d_actual = 0.0
+        d = cosine_angle._derivative(param="k", theta=theta_input, k=1000)
+        self.assertAlmostEqual(d, d_actual)
+
+        # test array r
+        theta_input = numpy.array([0.0, numpy.pi / 4, numpy.pi])
+        d_actual = numpy.array([2.0, 1.70710678119, 0])
+        d = cosine_angle._derivative(param="k", theta=theta_input, k=1000)
+        numpy.testing.assert_allclose(d, d_actual)
+
+        # test invalid param
+        with self.assertRaises(ValueError):
+            cosine_angle._derivative(
+                param="thetao",
+                theta=theta_input,
+                k=1000.0,
+            )
+
+    def test_json(self):
+        cosine_angle = relentless.model.potential.CosineAngle(types=("1",))
+        cosine_angle.coeff["1"].update(k=1000)
+        data = cosine_angle.to_json()
+
+        cosine_angle2 = relentless.model.potential.CosineAngle.from_json(data)
+        self.assertEqual(cosine_angle2.coeff["1"]["k"], 1000)
+
+
+class test_HarmonicCosineAngle(unittest.TestCase):
+    """Unit tests for relentless.model.potential.HarmonicCosineAngle"""
+
+    def test_init(self):
+        """Test creation from data"""
+        cosine_squred_angle = relentless.model.potential.HarmonicCosineAngle(
+            types=("1",)
+        )
+        coeff = relentless.model.potential.AngleParameters(
+            types=("1",),
+            params=(
+                "k",
+                "theta0",
+            ),
+        )
+        self.assertCountEqual(cosine_squred_angle.coeff.types, coeff.types)
+        self.assertCountEqual(cosine_squred_angle.coeff.params, coeff.params)
+
+    def test_energy(self):
+        """Test _energy method"""
+        cosine_squred_angle = relentless.model.potential.HarmonicCosineAngle(
+            types=("1",)
+        )
+        cosine_squred_angle.coeff["1"].update(k=1000, theta0=numpy.pi / 2)
+        # test scalar r
+        theta_input = numpy.pi
+        u_actual = 500
+        u = cosine_squred_angle.energy(type_=("1"), theta=theta_input)
+        self.assertAlmostEqual(u, u_actual)
+
+        # test array r
+        theta_input = numpy.array([0.0, numpy.pi / 2, numpy.pi])
+        u_actual = numpy.array([500.0, 0, 500.0])
+        u = cosine_squred_angle.energy(type_=("1"), theta=theta_input)
+        numpy.testing.assert_allclose(u, u_actual)
+
+    def test_force(self):
+        """Test _force method"""
+        cosine_squred_angle = relentless.model.potential.HarmonicCosineAngle(
+            types=("1",)
+        )
+        cosine_squred_angle.coeff["1"].update(k=1000, theta0=numpy.pi / 2)
+
+        # test scalar r
+        theta_input = numpy.pi
+        f_actual = 0
+        f = cosine_squred_angle.force(type_=("1"), theta=theta_input)
+        self.assertAlmostEqual(f, f_actual)
+
+        # test array r
+        theta_input = numpy.array([numpy.pi / 4, numpy.pi / 2, 3 * numpy.pi / 4])
+        f_actual = numpy.array([500, 0, -500])
+        f = cosine_squred_angle.force(type_=("1"), theta=theta_input)
+        numpy.testing.assert_allclose(f, f_actual)
+
+    def test_derivative(self):
+        """Test _derivative method"""
+        harmonic_cosine_angle = relentless.model.potential.HarmonicCosineAngle(
+            types=("1",)
+        )
+
+        # w.r.t. k
+        # test scalar r
+        theta_input = numpy.pi
+        d_actual = 0.5
+        d = harmonic_cosine_angle._derivative(
+            param="k", theta=theta_input, k=1000, theta0=numpy.pi / 2
+        )
+        self.assertAlmostEqual(d, d_actual)
+
+        # test array r
+        theta_input = numpy.array([numpy.pi / 4, numpy.pi / 2, numpy.pi])
+        d_actual = numpy.array([0.25, 0.0, 0.5])
+        d = harmonic_cosine_angle._derivative(
+            param="k", theta=theta_input, k=1000, theta0=numpy.pi / 2
+        )
+        numpy.testing.assert_allclose(d, d_actual)
+
+        # w.r.t. theta0
+        # test scalar r
+        theta_input = numpy.pi
+        d_actual = -1000
+        d = harmonic_cosine_angle._derivative(
+            param="theta0", theta=theta_input, k=1000.0, theta0=numpy.pi / 2
+        )
+        self.assertAlmostEqual(d, d_actual)
+
+        # test array r
+        theta_input = numpy.array([numpy.pi / 4, numpy.pi / 2, numpy.pi])
+        d_actual = numpy.array([707.106781185, 0, -1000])
+        d = harmonic_cosine_angle._derivative(
+            param="theta0", theta=theta_input, k=1000.0, theta0=numpy.pi / 2
+        )
+        numpy.testing.assert_allclose(d, d_actual)
+
+        # test invalid param
+        with self.assertRaises(ValueError):
+            harmonic_cosine_angle._derivative(
+                param="thetao", theta=theta_input, k=1000.0, theta0=1.0
+            )
+
+    def test_json(self):
+        cosine_squred_angle = relentless.model.potential.HarmonicCosineAngle(
+            types=("1",)
+        )
+        cosine_squred_angle.coeff["1"].update(k=1000, theta0=1.0)
+        data = cosine_squred_angle.to_json()
+
+        cosine_squred_angle2 = relentless.model.potential.HarmonicCosineAngle.from_json(
+            data
+        )
+        self.assertEqual(cosine_squred_angle2.coeff["1"]["k"], 1000)
+        self.assertEqual(cosine_squred_angle2.coeff["1"]["theta0"], 1.0)
+
+
+class test_AngleSpline(unittest.TestCase):
+    """Unit tests for relentless.model.potential.AngleSpline"""
+
+    def test_init(self):
+        """Test creation from data"""
+        # test diff mode
+        s = relentless.model.potential.AngleSpline(types=("1",), num_knots=3)
+        self.assertEqual(s.num_knots, 3)
+        self.assertEqual(s.mode, "diff")
+        coeff = relentless.model.potential.AngleParameters(
+            types=("1",),
+            params=(
+                "dtheta-0",
+                "dtheta-1",
+                "dtheta-2",
+                "diff-0",
+                "diff-1",
+                "diff-2",
+            ),
+        )
+        self.assertCountEqual(s.coeff.types, coeff.types)
+        self.assertCountEqual(s.coeff.params, coeff.params)
+
+        # test value mode
+        s = relentless.model.potential.AngleSpline(
+            types=("1",), num_knots=3, mode="value"
+        )
+        self.assertEqual(s.num_knots, 3)
+        self.assertEqual(s.mode, "value")
+        coeff = relentless.model.potential.AngleParameters(
+            types=("1",),
+            params=(
+                "dtheta-0",
+                "dtheta-1",
+                "dtheta-2",
+                "value-0",
+                "value-1",
+                "value-2",
+            ),
+        )
+        self.assertCountEqual(s.coeff.types, coeff.types)
+        self.assertCountEqual(s.coeff.params, coeff.params)
+
+        # test invalid number of knots
+        with self.assertRaises(ValueError):
+            s = relentless.model.potential.AngleSpline(types=("1",), num_knots=1)
+
+        # test invalid mode
+        with self.assertRaises(ValueError):
+            s = relentless.model.potential.AngleSpline(
+                types=("1",), num_knots=3, mode="val"
+            )
+
+    def test_from_array(self):
+        """Test from_array method and knots generator"""
+        theta_arr = [1, 2, 3]
+        u_arr = [1, 2, 3]
+
+        # test that bounds are enforced
+        s = relentless.model.potential.AngleSpline(types=("1",), num_knots=3)
+        with self.assertRaises(ValueError):
+            s.from_array(types=("1"), theta=theta_arr, u=u_arr)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/simulate/test_hoomd.py b/tests/simulate/test_hoomd.py
index 2dad12a1..d77f6acb 100644
--- a/tests/simulate/test_hoomd.py
+++ b/tests/simulate/test_hoomd.py
@@ -14,6 +14,7 @@
 from parameterized import parameterized_class
 
 import relentless
+from tests.model.potential.test_angle import LinPotAngle
 from tests.model.potential.test_bond import LinPotBond
 from tests.model.potential.test_pair import LinPot
 
@@ -110,6 +111,36 @@ def ens_pot_bonds(self):
 
         return (ens, pots)
 
+    def ens_pot_angles(self):
+        if self.dim == 3:
+            ens = relentless.model.Ensemble(
+                T=2.0, V=relentless.model.Cube(L=20.0), N={"A": 3, "B": 3}
+            )
+        elif self.dim == 2:
+            ens = relentless.model.Ensemble(
+                T=2.0, V=relentless.model.Square(L=20.0), N={"A": 3, "B": 3}
+            )
+        else:
+            raise ValueError("HOOMD supports 2d and 3d simulations")
+
+        # setup potentials
+        pot = LinPot(ens.types, params=("m",))
+        for pair in pot.coeff:
+            pot.coeff[pair].update({"m": -2.0, "rmax": 1.0})
+        pots = relentless.simulate.Potentials()
+        pots.pair = relentless.simulate.PairPotentialTabulator(
+            pot, start=0.0, stop=2.0, num=3, neighbor_buffer=0.4
+        )
+        angle = LinPotAngle(("angleA", "angleB"), params=("m",))
+        for type_ in angle.coeff:
+            angle.coeff[type_].update({"m": -2.0})
+        pots.angle = relentless.simulate.AnglePotentialTabulator(
+            angle,
+            num=3,
+        )
+
+        return (ens, pots)
+
     # mock gsd file for testing
     def create_gsd_file(self):
         filename = self.directory.file("test.gsd")
@@ -171,6 +202,45 @@ def create_gsd_file_bonds(self):
         relentless.mpi.world.barrier()
         return filename
 
+    # mock gsd file with topology for testing
+    def create_gsd_file_angles(self):
+        filename = self.directory.file("test.gsd")
+        if relentless.mpi.world.rank_is_root:
+            with gsd.hoomd.open(name=filename, mode=gsd_write_mode) as f:
+                s = HOOMDFrame()
+                s.particles.N = 6
+                s.particles.types = ["A", "B"]
+                s.particles.typeid = [0, 0, 0, 1, 1, 1]
+                s.angles.N = 2
+                s.angles.types = ["angleA", "angleB"]
+                s.angles.group = [(0, 1, 2), (3, 4, 5)]
+                s.angles.typeid = [0, 1]
+                if self.dim == 3:
+                    s.particles.position = [
+                        [0, 0, 0],
+                        [1, 0, 0],
+                        [0, 1, 0],
+                        [0, 2, 1],
+                        [1, 2, 1],
+                        [0, 3, 1],
+                    ]
+                    s.configuration.box = [20, 20, 20, 0, 0, 0]
+                elif self.dim == 2:
+                    s.particles.position = [
+                        [0, 0, 0],
+                        [1, 0, 0],
+                        [0, 1, 0],
+                        [0, 2, 0],
+                        [1, 2, 0],
+                        [0, 3, 0],
+                    ]
+                    s.configuration.box = [20, 20, 0, 0, 0, 0]
+                else:
+                    raise ValueError("HOOMD supports 2d and 3d simulations")
+                f.append(s)
+        relentless.mpi.world.barrier()
+        return filename
+
     def create_lammps_file(self):
         file_ = self.directory.file("test.data")
 
@@ -226,6 +296,52 @@ def create_lammps_file_bonds(self):
 
         return file_
 
+    def create_lammps_file_angles(self):
+        file_ = self.directory.file("test.data")
+
+        if relentless.mpi.world.rank_is_root:
+            low = [-5, -5, -5 if self.dim == 3 else -0.1]
+            high = [5, 5, 5 if self.dim == 3 else 0.1]
+            snap = lammpsio.Snapshot(N=6, box=lammpsio.Box(low, high))
+            snap.position = [
+                [0, 0, 0],
+                [1, 0, 0],
+                [0, 1, 0],
+                [0, 2, 0],
+                [1, 2, 0],
+                [0, 3, 0],
+            ]
+            if self.dim == 3:
+                snap.position = [
+                    [0, 0, 0],
+                    [1, 0, 0],
+                    [0, 1, 0],
+                    [0, 2, 1],
+                    [1, 2, 1],
+                    [0, 3, 1],
+                ]
+            snap.typeid = [
+                1,
+                1,
+                1,
+                2,
+                2,
+                2,
+            ]
+            snap.mass = [0.3, 0.3, 0.3, 0.1, 0.1, 0.1]
+            snap.angles = lammpsio.topology.Angles(N=2, num_types=2)
+            snap.angles.id = [1, 2]
+            snap.angles.typeid = [1, 2]
+            snap.angles.members = [(1, 2, 3), (4, 5, 6)]
+            snap.angles.type_label = lammpsio.topology.LabelMap(
+                {1: "angleA", 2: "angleB"}
+            )
+
+            lammpsio.DataFile.create(file_, snap)
+        relentless.mpi.world.barrier()
+
+        return file_
+
     def test_initialize_from_gsd_file(self):
         ens, pot = self.ens_pot()
         f = self.create_gsd_file()
@@ -243,7 +359,7 @@ def test_initialize_from_gsd_file(self):
             h.initializer = relentless.simulate.InitializeFromFile(pathlib.Path(f).name)
             h.run(pot, d)
 
-    def test_initialize_from_gsd_file_topology(self):
+    def test_initialize_from_gsd_file_bonds(self):
         ens, pot = self.ens_pot_bonds()
         f = self.create_gsd_file_bonds()
         op = relentless.simulate.InitializeFromFile(filename=f)
@@ -260,6 +376,23 @@ def test_initialize_from_gsd_file_topology(self):
             h.initializer = relentless.simulate.InitializeFromFile(pathlib.Path(f).name)
             h.run(pot, d)
 
+    def test_initialize_from_gsd_file_angles(self):
+        ens, pot = self.ens_pot_angles()
+        f = self.create_gsd_file_angles()
+        op = relentless.simulate.InitializeFromFile(filename=f)
+        h = relentless.simulate.HOOMD(op)
+        h.run(pot, self.directory)
+
+        # Run in a different directory
+        with self.directory:
+            d = self.directory.directory(
+                "run", create=relentless.mpi.world.rank_is_root
+            )
+            relentless.mpi.world.barrier()
+            op.filename = pathlib.Path(f).name
+            h.initializer = relentless.simulate.InitializeFromFile(pathlib.Path(f).name)
+            h.run(pot, d)
+
     def test_initialize_from_lammps_file(self):
         """Test running initialization simulation operations."""
         pot = LinPot(("1", "2"), params=("m",))
@@ -275,7 +408,7 @@ def test_initialize_from_lammps_file(self):
         h = relentless.simulate.HOOMD(op)
         h.run(pots, self.directory)
 
-    def test_initialize_from_lammps_file_topology(self):
+    def test_initialize_from_lammps_file_bonds(self):
         """Test running initialization simulation operations."""
         pot = LinPot(("1", "2"), params=("m",))
         for pair in pot.coeff:
@@ -290,6 +423,21 @@ def test_initialize_from_lammps_file_topology(self):
         h = relentless.simulate.HOOMD(op)
         h.run(pots, self.directory)
 
+    def test_initialize_from_lammps_file_angles(self):
+        """Test running initialization simulation operations."""
+        pot = LinPot(("1", "2"), params=("m",))
+        for pair in pot.coeff:
+            pot.coeff[pair].update({"m": -2.0, "rmax": 1.0})
+        pots = relentless.simulate.Potentials()
+        pots.pair = relentless.simulate.PairPotentialTabulator(
+            pot, start=1e-6, stop=2.0, num=10, neighbor_buffer=0.1
+        )
+
+        f = self.create_lammps_file_angles()
+        op = relentless.simulate.InitializeFromFile(filename=f)
+        h = relentless.simulate.HOOMD(op)
+        h.run(pots, self.directory)
+
     def test_initialize_randomly(self):
         ens, pot = self.ens_pot()
         op = relentless.simulate.InitializeRandomly(seed=1, N=ens.N, V=ens.V, T=ens.T)
@@ -451,7 +599,7 @@ def test_molecular_dynamics(self):
                 vrl.barostat = relentless.simulate.MTKBarostat(P=1, tau=0.5)
                 h.run(pot, self.directory)
 
-    def test_topology_run(self):
+    def test_bonds_run(self):
         ens, pot = self.ens_pot_bonds()
         f = self.create_gsd_file_bonds()
         init = relentless.simulate.InitializeFromFile(filename=f)
@@ -464,6 +612,19 @@ def test_topology_run(self):
         lgv.friction = {"A": 1.5, "B": 2.5}
         h.run(pot, self.directory)
 
+    def test_angles_run(self):
+        ens, pot = self.ens_pot_angles()
+        f = self.create_gsd_file_angles()
+        init = relentless.simulate.InitializeFromFile(filename=f)
+        lgv = relentless.simulate.RunLangevinDynamics(
+            steps=1, timestep=1.0e-3, T=ens.T, friction=1.0, seed=2
+        )
+        h = relentless.simulate.HOOMD(init, lgv)
+        h.run(pot, self.directory)
+
+        lgv.friction = {"A": 1.5, "B": 2.5}
+        h.run(pot, self.directory)
+
     def test_temperature_ramp(self):
         if self.dim == 3:
             V = relentless.model.Cube(100.0)
diff --git a/tests/simulate/test_lammps.py b/tests/simulate/test_lammps.py
index 7a46d043..771e757b 100644
--- a/tests/simulate/test_lammps.py
+++ b/tests/simulate/test_lammps.py
@@ -34,6 +34,7 @@
 from packaging import version
 
 import relentless
+from tests.model.potential.test_angle import LinPotAngle
 from tests.model.potential.test_bond import LinPotBond
 from tests.model.potential.test_pair import LinPot
 
@@ -139,6 +140,36 @@ def ens_pot_bonds(self):
 
         return (ens, pots)
 
+    def ens_pot_angles(self):
+        if self.dim == 3:
+            ens = relentless.model.Ensemble(
+                T=2.0, V=relentless.model.Cube(L=20.0), N={"1": 3, "2": 3}
+            )
+        elif self.dim == 2:
+            ens = relentless.model.Ensemble(
+                T=2.0, V=relentless.model.Square(L=20.0), N={"1": 3, "2": 3}
+            )
+        else:
+            raise ValueError("HOOMD supports 2d and 3d simulations")
+
+        # setup potentials
+        pot = LinPot(ens.types, params=("m",))
+        for pair in pot.coeff:
+            pot.coeff[pair].update({"m": -2.0, "rmax": 1.0})
+        pots = relentless.simulate.Potentials()
+        pots.pair = relentless.simulate.PairPotentialTabulator(
+            pot, start=1e-6, stop=2.0, num=3, neighbor_buffer=0.4
+        )
+        angle = LinPotAngle(("angle1", "angle2"), params=("m",))
+        for type_ in angle.coeff:
+            angle.coeff[type_].update({"m": -1.0})
+        pots.angle = relentless.simulate.AnglePotentialTabulator(
+            angle,
+            num=10,
+        )
+
+        return (ens, pots)
+
     def create_gsd_file(self):
         filename = self.directory.file("test.gsd")
         if relentless.mpi.world.rank_is_root:
@@ -199,6 +230,44 @@ def create_gsd_file_bonds(self):
         relentless.mpi.world.barrier()
         return filename
 
+    def create_gsd_file_angles(self):
+        filename = self.directory.file("test.gsd")
+        if relentless.mpi.world.rank_is_root:
+            with gsd.hoomd.open(name=filename, mode=gsd_write_mode) as f:
+                s = HOOMDFrame()
+                s.particles.N = 6
+                s.particles.types = ["A", "B"]
+                s.particles.typeid = [0, 0, 0, 1, 1, 1]
+                s.angles.N = 2
+                s.angles.types = ["angleA", "angleB"]
+                s.angles.group = [(0, 1, 2), (3, 4, 5)]
+                s.angles.typeid = [0, 1]
+                if self.dim == 3:
+                    s.particles.position = [
+                        [0, 0, 0],
+                        [1, 0, 0],
+                        [0, 1, 0],
+                        [0, 2, 1],
+                        [1, 2, 1],
+                        [0, 3, 1],
+                    ]
+                    s.configuration.box = [20, 20, 20, 0, 0, 0]
+                elif self.dim == 2:
+                    s.particles.position = [
+                        [0, 0, 0],
+                        [1, 0, 0],
+                        [0, 1, 0],
+                        [0, 2, 0],
+                        [1, 2, 0],
+                        [0, 3, 0],
+                    ]
+                    s.configuration.box = [20, 20, 0, 0, 0, 0]
+                else:
+                    raise ValueError("HOOMD supports 2d and 3d simulations")
+                f.append(s)
+        relentless.mpi.world.barrier()
+        return filename
+
     def create_lammps_file(self):
         file_ = self.directory.file("test.data")
 
@@ -238,6 +307,52 @@ def create_lammps_file_bonds(self):
 
         return file_
 
+    def create_lammps_file_angles(self):
+        file_ = self.directory.file("test.data")
+
+        if relentless.mpi.world.rank_is_root:
+            low = [-5, -5, -5 if self.dim == 3 else -0.1]
+            high = [5, 5, 5 if self.dim == 3 else 0.1]
+            snap = lammpsio.Snapshot(N=6, box=lammpsio.Box(low, high))
+            snap.position = [
+                [0, 0, 0],
+                [1, 0, 0],
+                [0, 1, 0],
+                [1, 1, 0],
+                [2, 1, 0],
+                [1, 2, 0],
+            ]
+            if self.dim == 3:
+                snap.position = [
+                    [0, 0, 0],
+                    [1, 0, 0],
+                    [0, 1, 0],
+                    [0, 2, 1],
+                    [1, 2, 1],
+                    [0, 3, 1],
+                ]
+            snap.typeid = [
+                1,
+                1,
+                1,
+                2,
+                2,
+                2,
+            ]
+            snap.mass = [0.3, 0.3, 0.3, 0.1, 0.1, 0.1]
+            snap.angles = lammpsio.topology.Angles(N=2, num_types=2)
+            snap.angles.id = [1, 2]
+            snap.angles.typeid = [1, 2]
+            snap.angles.members = [(1, 2, 3), (4, 5, 6)]
+            snap.angles.type_label = lammpsio.topology.LabelMap(
+                {1: "angleA", 2: "angleB"}
+            )
+
+            lammpsio.DataFile.create(file_, snap)
+        relentless.mpi.world.barrier()
+
+        return file_
+
     def test_initialize_from_lammps_file(self):
         """Test running initialization simulation operations."""
         # InitializeFromFile
@@ -269,7 +384,7 @@ def test_initialize_from_lammps_file(self):
         lmp = relentless.simulate.LAMMPS(op, executable=self.executable)
         lmp.run(potentials=pot, directory=self.directory)
 
-    def test_initialize_from_lammps_file_topology(self):
+    def test_initialize_from_lammps_file_bonds(self):
         """Test running initialization simulation operations."""
         # InitializeFromFile
         ens, pot = self.ens_pot_bonds()
@@ -296,6 +411,33 @@ def test_initialize_from_lammps_file_topology(self):
 
             lmp.run(potentials=pot, directory=d)
 
+    def test_initialize_from_lammps_file_angles(self):
+        """Test running initialization simulation operations."""
+        # InitializeFromFile
+        ens, pot = self.ens_pot_angles()
+        f = self.create_lammps_file_angles()
+        op = relentless.simulate.InitializeFromFile(filename=f)
+        h = relentless.simulate.LAMMPS(op, atom_style="molecular")
+        h.run(pot, self.directory)
+
+        lmp = relentless.simulate.LAMMPS(
+            op,
+            executable=self.executable,
+            atom_style="molecular",
+        )
+        # Run in a different directory
+        with self.directory:
+            d = self.directory.directory(
+                "run", create=relentless.mpi.world.rank_is_root
+            )
+            relentless.mpi.world.barrier()
+            op.filename = pathlib.Path(f).name
+            lmp.initializer = relentless.simulate.InitializeFromFile(
+                pathlib.Path(f).name
+            )
+
+            lmp.run(potentials=pot, directory=d)
+
     def test_initialize_from_gsd_file(self):
         """Test running initialization simulation operations."""
         _, pot = self.ens_pot()
@@ -308,7 +450,7 @@ def test_initialize_from_gsd_file(self):
         )
         lmp.run(potentials=pot, directory=self.directory)
 
-    def test_initialize_from_gsd_file_topology(self):
+    def test_initialize_from_gsd_file_bonds(self):
         """Test running initialization simulation operations."""
         # setup potentials
         pot = LinPot(("A", "B"), params=("m",))
@@ -333,6 +475,29 @@ def test_initialize_from_gsd_file_topology(self):
         h = relentless.simulate.LAMMPS(op, atom_style="molecular")
         h.run(pots, self.directory)
 
+    def test_initialize_from_gsd_file_angles(self):
+        """Test running initialization simulation operations."""
+        # setup potentials
+        pot = LinPot(("A", "B"), params=("m",))
+        for pair in pot.coeff:
+            pot.coeff[pair].update({"m": -2.0, "rmax": 1.0})
+        pots = relentless.simulate.Potentials()
+        pots.pair = relentless.simulate.PairPotentialTabulator(
+            pot, start=1e-6, stop=2.0, num=3, neighbor_buffer=0.4
+        )
+        angle = LinPotAngle(("angleA", "angleB"), params=("m",))
+        for type_ in angle.coeff:
+            angle.coeff[type_].update({"m": -2.0})
+        pots.angle = relentless.simulate.AnglePotentialTabulator(
+            angle,
+            num=3,
+        )
+
+        f = self.create_gsd_file_angles()
+        op = relentless.simulate.InitializeFromFile(filename=f)
+        h = relentless.simulate.LAMMPS(op, atom_style="molecular")
+        h.run(pots, self.directory)
+
     def test_random_initialize_options(self):
         # no T
         ens, pot = self.ens_pot()
@@ -560,7 +725,7 @@ def test_molecular_dynamics(self):
         vrl.thermostat = None
         lmp.run(pot, self.directory)
 
-    def test_topology_run(self):
+    def test_bond_run(self):
         ens, pot = self.ens_pot_bonds()
         f = self.create_lammps_file_bonds()
         init = relentless.simulate.InitializeFromFile(filename=f)
@@ -571,6 +736,17 @@ def test_topology_run(self):
         lmp.operations = vrl
         lmp.run(pot, self.directory)
 
+    def test_angle_run(self):
+        ens, pot = self.ens_pot_angles()
+        f = self.create_lammps_file_angles()
+        init = relentless.simulate.InitializeFromFile(filename=f)
+        lmp = relentless.simulate.LAMMPS(init, atom_style="molecular")
+
+        # VerletIntegrator
+        vrl = relentless.simulate.RunMolecularDynamics(steps=1, timestep=1e-4)
+        lmp.operations = vrl
+        lmp.run(pot, self.directory)
+
     def test_temperature_ramp(self):
         if self.dim == 3:
             V = relentless.model.Cube(100.0)