diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 31479f9952..ddc48bc4de 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -190,7 +190,7 @@ class TransportOperator(ABC):
 
         Returns
         -------
-        volume : list of float
+        volume : dict of str to float
             Volumes corresponding to materials in burn_list
         nuc_list : list of str
             A list of all nuclide names. Used for sorting the simulation.
diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index f018aaac7b..960cf06456 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -10,7 +10,7 @@ import math
 import re
 from collections import OrderedDict, defaultdict
 from collections.abc import Mapping, Iterable
-from numbers import Real
+from numbers import Real, Integral
 from warnings import warn
 
 from openmc.checkvalue import check_type, check_greater_than
@@ -500,7 +500,7 @@ class Chain:
                 # Gain
                 for _, target, branching_ratio in nuc.decay_modes:
                     # Allow for total annihilation for debug purposes
-                    if target != 'Nothing':
+                    if target is not None:
                         branch_val = branching_ratio * decay_constant
 
                         if branch_val != 0.0:
@@ -525,17 +525,16 @@ class Chain:
                             matrix[i, i] -= path_rate
 
                     # Gain term; allow for total annihilation for debug purposes
-                    if target != 'Nothing':
-                        if r_type != 'fission':
-                            if path_rate != 0.0:
-                                k = self.nuclide_dict[target]
-                                matrix[k, i] += path_rate * br
-                        else:
-                            for product, y in fission_yields[nuc.name].items():
-                                yield_val = y * path_rate
-                                if yield_val != 0.0:
-                                    k = self.nuclide_dict[product]
-                                    matrix[k, i] += yield_val
+                    if r_type != 'fission':
+                        if target is not None and path_rate != 0.0:
+                            k = self.nuclide_dict[target]
+                            matrix[k, i] += path_rate * br
+                    else:
+                        for product, y in fission_yields[nuc.name].items():
+                            yield_val = y * path_rate
+                            if yield_val != 0.0:
+                                k = self.nuclide_dict[product]
+                                matrix[k, i] += yield_val
 
                 # Clear set of reactions
                 reactions.clear()
@@ -821,3 +820,168 @@ class Chain:
                 return stat
             valid = valid and stat
         return valid
+
+    def reduce(self, initial_isotopes, level=None):
+        """Reduce the size of the chain by following transmutation paths
+
+        As an example, consider a simple chain with the following
+        isotopes and transmutation paths::
+
+            U235 (n,gamma) U236
+                 (n,fission) (Xe135, I135, Cs135)
+            I135 (beta decay) Xe135 (beta decay) Cs135
+            Xe135 (n,gamma) Xe136
+
+        Calling ``chain.reduce(["I135"])`` will produce a depletion
+        chain that contains only isotopes that would originate from
+        I135: I135, Xe135, Cs135, and Xe136. U235 and U236 will not
+        be included, but multiple isotopes can be used to start
+        the search.
+
+        The ``level`` value controls the depth of the search.
+        ``chain.reduce(["U235"], level=1)`` would return a chain
+        with all isotopes except Xe136, since it is two transmutations
+        removed from U235 in this case.
+
+        While targets will not be included in the new chain, the
+        total destruction rate and decay rate of included isotopes
+        will be preserved.
+
+        Parameters
+        ----------
+        initial_isotopes : iterable of str
+            Start the search based on the contents of these isotopes
+        level : int, optional
+            Depth of transmuation path to follow. Must be greater than
+            or equal to zero. A value of zero returns a chain with
+            ``initial_isotopes``. The default value of None implies
+            that all isotopes that appear in the transmutation paths
+            of the initial isotopes and their progeny should be
+            explored
+
+        Returns
+        -------
+        Chain
+            Depletion chain containing isotopes that would appear
+            after following up to ``level`` reactions and decay paths
+
+        """
+        check_type("initial_isotopes", initial_isotopes, Iterable, str)
+        if level is None:
+            level = math.inf
+        else:
+            check_type("level", level, Integral)
+            check_greater_than("level", level, 0, equality=True)
+
+        all_isotopes = self._follow(set(initial_isotopes), level)
+
+        # Avoid re-sorting for fission yields
+        name_sort = sorted(all_isotopes)
+
+        nuclides = []
+        nuclide_dict = {}
+        reactions = set()
+
+        for idx, iso in enumerate(sorted(all_isotopes, key=openmc.data.zam)):
+            previous = self[iso]
+            new_nuclide = Nuclide(previous.name)
+            new_nuclide.half_life = previous.half_life
+            new_nuclide.decay_energy = new_nuclide.decay_energy
+
+            new_decay = []
+            for mode in previous.decay_modes:
+                if mode.target in all_isotopes:
+                    new_decay.append(mode)
+                else:
+                    new_decay.append(DecayTuple(
+                        mode.type, None, mode.branching_ratio))
+            new_nuclide.decay_modes = new_decay
+
+            new_reactions = []
+            for rxn in previous.reactions:
+                if rxn.target in all_isotopes:
+                    new_reactions.append(rxn)
+                    reactions.add(rxn.type)
+                elif rxn.type == "fission":
+                    new_yields = new_nuclide.yield_data = (
+                        previous.yield_data.restrict_products(name_sort))
+                    if new_yields is not None:
+                        new_reactions.append(rxn)
+                        reactions.add("fission")
+                # Maintain total destruction rates but set no target
+                else:
+                    new_reactions.append(ReactionTuple(
+                        rxn.type, None, rxn.Q, rxn.branching_ratio))
+                    reactions.add(rxn.type)
+
+            new_nuclide.reactions = new_reactions
+
+            nuclides.append(new_nuclide)
+            nuclide_dict[iso] = idx
+
+        new_chain = type(self)()
+        new_chain.nuclides = nuclides
+        new_chain.nuclide_dict = nuclide_dict
+
+        # Doesn't appear that the ordering matters for the reactions,
+        # just the contents
+        new_chain.reactions = sorted(reactions)
+
+        return new_chain
+
+    def _follow(self, isotopes, level):
+        """Return all isotopes present up to depth level"""
+        found = isotopes.copy()
+        remaining = set(self.nuclide_dict)
+        if not found.issubset(remaining):
+            raise IndexError(
+                "The following isotopes were not found in the chain: "
+                "{}".format(", ".join(found - remaining)))
+
+        if level == 0:
+            return found
+
+        remaining -= found
+
+        depth = 0
+        next_iso = set()
+
+        while depth < level and remaining:
+            # Exhaust all isotopes at this level
+            while isotopes:
+                iso = isotopes.pop()
+                found.add(iso)
+                nuclide = self[iso]
+
+                # Follow all transmutation paths for this nuclide
+                for rxn in nuclide.reactions + nuclide.decay_modes:
+                    if rxn.type == "fission" or rxn.target is None:
+                        continue
+                    # Skip if we've already come across this isotope
+                    elif (rxn.target in next_iso
+                          or rxn.target in found or rxn.target in isotopes):
+                        continue
+                    next_iso.add(rxn.target)
+
+                if nuclide.yield_data is not None:
+                    for product in nuclide.yield_data.products:
+                        if (product in next_iso
+                                or product in found or product in isotopes):
+                            continue
+                        next_iso.add(product)
+
+            if not next_iso:
+                # No additional isotopes to process, nor to update the
+                # current set of discovered isotopes
+                return found
+
+            # Prepare for next dig
+            depth += 1
+            isotopes |= next_iso
+            remaining -= next_iso
+            next_iso.clear()
+
+        # Process isotope that would have started next depth
+        found.update(isotopes)
+
+        return found
diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
index 52dd2bf96d..854a2abb07 100644
--- a/openmc/deplete/nuclide.py
+++ b/openmc/deplete/nuclide.py
@@ -13,7 +13,7 @@ try:
 except ImportError:
     import xml.etree.ElementTree as ET
 
-from numpy import empty
+from numpy import empty, searchsorted
 
 from openmc.checkvalue import check_type
 
@@ -30,8 +30,10 @@ Parameters
 ----------
 type : str
     Type of the decay mode, e.g., 'beta-'
-target : str
-    Nuclide resulting from decay
+target : str or None
+    Nuclide resulting from decay. A value of ``None`` implies the
+    target does not exist in the currently configured depletion
+    chain
 branching_ratio : float
     Branching ratio of the decay mode
 
@@ -53,8 +55,11 @@ Parameters
 ----------
 type : str
     Type of the reaction, e.g., 'fission'
-target : str
-    nuclide resulting from reaction
+target : str or None
+    Nuclide resulting from reaction. A value of ``None``
+    implies either no single target, e.g. from fission,
+    or that the target nuclide is not considered
+    in the current depletion chain
 Q : float
     Q value of the reaction in [eV]
 branching_ratio : float
@@ -179,6 +184,8 @@ class Nuclide:
         for decay_elem in element.iter('decay'):
             d_type = decay_elem.get('type')
             target = decay_elem.get('target')
+            if target is not None and target.lower() == "nothing":
+                target = None
             branching_ratio = float(decay_elem.get('branching_ratio'))
             nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))
 
@@ -192,6 +199,8 @@ class Nuclide:
             # just set null values
             if r_type != 'fission':
                 target = reaction_elem.get('target')
+                if target is not None and target.lower() == "nothing":
+                    target = None
             else:
                 target = None
                 if fission_q is not None:
@@ -226,7 +235,7 @@ class Nuclide:
             for mode, daughter, br in self.decay_modes:
                 mode_elem = ET.SubElement(elem, 'decay')
                 mode_elem.set('type', mode)
-                mode_elem.set('target', daughter)
+                mode_elem.set('target', daughter or "Nothing")
                 mode_elem.set('branching_ratio', str(br))
 
         elem.set('reactions', str(len(self.reactions)))
@@ -234,7 +243,7 @@ class Nuclide:
             rx_elem = ET.SubElement(elem, 'reaction')
             rx_elem.set('type', rx)
             rx_elem.set('Q', str(Q))
-            if rx != 'fission':
+            if rx != 'fission' or daughter is not None:
                 rx_elem.set('target', daughter)
             if br != 1.0:
                 rx_elem.set('branching_ratio', str(br))
@@ -397,9 +406,7 @@ class FissionYieldDistribution(Mapping):
         for g_index, energy in enumerate(energies):
             prod_map = fission_yields[energy]
             for prod_ix, product in enumerate(ordered_prod):
-                yield_val = prod_map.get(product)
-                yield_matrix[g_index, prod_ix] = (
-                    0.0 if yield_val is None else yield_val)
+                yield_matrix[g_index, prod_ix] = prod_map.get(product, 0.0)
         self.energies = tuple(energies)
         self.products = tuple(ordered_prod)
         self.yield_matrix = yield_matrix
@@ -435,7 +442,7 @@ class FissionYieldDistribution(Mapping):
         FissionYieldDistribution
         """
         all_yields = {}
-        for elem_index, yield_elem in enumerate(element.iter("fission_yields")):
+        for yield_elem in element.iter("fission_yields"):
             energy = float(yield_elem.get("energy"))
             products = yield_elem.find("products").text.split()
             yields = map(float, yield_elem.find("data").text.split())
@@ -460,6 +467,37 @@ class FissionYieldDistribution(Mapping):
             data_elem = ET.SubElement(yield_element, "data")
             data_elem.text = " ".join(map(str, yield_obj.yields))
 
+    def restrict_products(self, possible_products):
+        """Return a new distribution with select products
+
+        Parameters
+        ----------
+        possible_products : iterable of str
+            Candidate pool of fission products. Existing products
+            not contained here will not exist in the new instance
+
+        Returns
+        -------
+        FissionYieldDistribution or None
+            A value of None indicates no values in
+            ``possible_products`` exist in :attr:`products`
+
+        """
+
+        overlap = set(self.products).intersection(possible_products)
+        if not overlap:
+            return None
+
+        products = sorted(overlap)
+        indices = searchsorted(self.products, products)
+
+        # coerce back to dictionary to pass back to __init__
+        new_yields = {}
+        for ene, yields in zip(self.energies, self.yield_matrix.copy()):
+            new_yields[ene] = dict(zip(products, yields[indices]))
+
+        return type(self)(new_yields)
+
 
 class FissionYield(Mapping):
     """Mapping for fission yields of a parent at a specific energy
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 8eb033783f..ce7e709d9c 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -10,13 +10,10 @@ densities is all done in-memory instead of through the filesystem.
 import sys
 import copy
 from collections import OrderedDict
-from itertools import chain
 import os
-import time
 import xml.etree.ElementTree as ET
 from warnings import warn
 
-import h5py
 import numpy as np
 from uncertainties import ufloat
 
@@ -113,6 +110,13 @@ class Operator(TransportOperator):
         ``fission_yield_mode``. Will be passed directly on to the
         helper. Passing a value of None will use the defaults for
         the associated helper.
+    reduce_chain : bool, optional
+        If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
+        depletion chain up to ``reduce_chain_level``. Default is False.
+    reduce_chain_level : int, optional
+        Depth of the search when reducing the depletion chain. Only used
+        if ``reduce_chain`` evaluates to true. The default value of
+        ``None`` implies no limit on the depth.
 
     Attributes
     ----------
@@ -158,7 +162,8 @@ class Operator(TransportOperator):
     def __init__(self, geometry, settings, chain_file=None, prev_results=None,
                  diff_burnable_mats=False, energy_mode="fission-q",
                  fission_q=None, dilute_initial=1.0e3,
-                 fission_yield_mode="constant", fission_yield_opts=None):
+                 fission_yield_mode="constant", fission_yield_opts=None,
+                 reduce_chain=False, reduce_chain_level=None):
         if fission_yield_mode not in self._fission_helpers:
             raise KeyError(
                 "fission_yield_mode must be one of {}, not {}".format(
@@ -179,6 +184,16 @@ class Operator(TransportOperator):
         self.geometry = geometry
         self.diff_burnable_mats = diff_burnable_mats
 
+        # Reduce the chain before we create more materials
+        if reduce_chain:
+            all_isotopes = set()
+            for material in geometry.get_all_materials().values():
+                if not material.depletable:
+                    continue
+                for name, _dens_percent, _dens_type in material.nuclides:
+                    all_isotopes.add(name)
+            self.chain = self.chain.reduce(all_isotopes, reduce_chain_level)
+
         # Differentiate burnable materials with multiple instances
         if self.diff_burnable_mats:
             self._differentiate_burnable_mats()
diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py
index 2f102913de..28b489d1f1 100644
--- a/openmc/deplete/results.py
+++ b/openmc/deplete/results.py
@@ -36,7 +36,7 @@ class Results:
         Number of nuclides.
     rates : list of ReactionRates
         The reaction rates for each substep.
-    volume : OrderedDict of int to float
+    volume : OrderedDict of str to float
         Dictionary mapping mat id to volume.
     mat_to_ind : OrderedDict of str to int
         A dictionary mapping mat ID as string to index.
diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py
index 1b5b81ac4c..db1f0ff9db 100644
--- a/openmc/deplete/results_list.py
+++ b/openmc/deplete/results_list.py
@@ -2,7 +2,7 @@ import h5py
 import numpy as np
 
 from .results import Results, _VERSION_RESULTS
-from openmc.checkvalue import check_filetype_version
+from openmc.checkvalue import check_filetype_version, check_value
 
 
 __all__ = ["ResultsList"]
@@ -40,7 +40,7 @@ class ResultsList(list):
                 new.append(Results.from_hdf5(fh, i))
         return new
 
-    def get_atoms(self, mat, nuc):
+    def get_atoms(self, mat, nuc, nuc_units="atoms", time_units="s"):
         """Get number of nuclides over time from a single material
 
         .. note::
@@ -57,15 +57,24 @@ class ResultsList(list):
             Material name to evaluate
         nuc : str
             Nuclide name to evaluate
+        nuc_units : {"atoms", "atom/b-cm", "atom/cm3"}, optional
+            Units for the returned concentration. Default is ``"atoms"``
+        time_units : {"s", "min", "h", "d"}, optional
+            Units for the returned time array. Default is ``"s"`` to
+            return the value in seconds.
 
         Returns
         -------
         time : numpy.ndarray
-            Array of times in [s]
+            Array of times in units of ``time_units``
         concentration : numpy.ndarray
-            Total number of atoms for specified nuclide
+            Concentration of specified nuclide in units of ``nuc_units``
 
         """
+        check_value("time_units", time_units, {"s", "d", "min", "h"})
+        check_value("nuc_units", nuc_units,
+                    {"atoms", "atom/b-cm", "atom/cm3"})
+
         time = np.empty_like(self, dtype=float)
         concentration = np.empty_like(self, dtype=float)
 
@@ -74,6 +83,21 @@ class ResultsList(list):
             time[i] = result.time[0]
             concentration[i] = result[0, mat, nuc]
 
+        # Unit conversions
+        if time_units == "d":
+            time /= (60 * 60 * 24)
+        elif time_units == "h":
+            time /= (60 * 60)
+        elif time_units == "min":
+            time /= 60
+
+        if nuc_units != "atoms":
+            # Divide by volume to get density
+            concentration /= self[0].volume[mat]
+            if nuc_units == "atom/b-cm":
+                # 1 barn = 1e-24 cm^2
+                concentration *= 1e-24
+
         return time, concentration
 
     def get_reaction_rate(self, mat, nuc, rx):
diff --git a/openmc/material.py b/openmc/material.py
index cbd002c0bd..aabbf2e630 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -1,9 +1,10 @@
-from collections import OrderedDict, defaultdict, namedtuple
+from collections import OrderedDict, defaultdict, namedtuple, Counter
 from collections.abc import Iterable
 from copy import deepcopy
 from numbers import Real, Integral
 from pathlib import Path
 import warnings
+import re
 from xml.etree import ElementTree as ET
 
 import numpy as np
@@ -556,6 +557,100 @@ class Material(IDManagerMixin):
                                       enrichment_type):
             self.add_nuclide(*nuclide)
 
+    def add_elements_from_formula(self, formula, percent_type='ao', enrichment=None,
+                                  enrichment_target=None, enrichment_type=None):
+        """Add a elements from a chemical formula to the material.
+
+        Parameters
+        ----------
+        formula : str
+            Formula to add, e.g., 'C2O', 'C6H12O6', or (NH4)2SO4.
+            Note this is case sensitive, elements must start with an uppercase
+            character. Multiplier numbers must be integers.
+        percent_type : {'ao', 'wo'}, optional
+            'ao' for atom percent and 'wo' for weight percent. Defaults to atom
+            percent.
+        enrichment : float, optional
+            Enrichment of an enrichment_target nuclide in percent (ao or wo).
+            If enrichment_target is not supplied then it is enrichment for U235
+            in weight percent. For example, input 4.95 for 4.95 weight percent
+            enriched U. Default is None (natural composition).
+        enrichment_target : str, optional
+            Single nuclide name to enrich from a natural composition (e.g., 'O16')
+        enrichment_type : {'ao', 'wo'}, optional
+            'ao' for enrichment as atom percent and 'wo' for weight percent.
+            Default is: 'ao' for two-isotope enrichment; 'wo' for U enrichment
+
+        Notes
+        -----
+        General enrichment procedure is allowed only for elements composed of
+        two isotopes. If `enrichment_target` is given without `enrichment`
+        natural composition is added to the material.
+
+        """
+        cv.check_type('formula', formula, str)
+
+        if '.' in formula:
+            msg = 'Non-integer multiplier values are not accepted. The ' \
+                  'input formula {} contains a "." character.'.format(formula)
+            raise ValueError(msg)
+
+        # Tokenizes the formula and check validity of tokens
+        tokens = re.findall(r"([A-Z][a-z]*)(\d*)|(\()|(\))(\d*)", formula)
+        for row in tokens:
+            for token in row:
+                if token.isalpha():
+                    if token == "n" or token not in openmc.data.ATOMIC_NUMBER:
+                        msg = 'Formula entry {} not an element symbol.' \
+                              .format(token)
+                        raise ValueError(msg)
+                elif token not in ['(', ')', ''] and not token.isdigit():
+                        msg = 'Formula must be made from a sequence of ' \
+                              'element symbols, integers, and backets. ' \
+                              '{} is not an allowable entry.'.format(token)
+                        raise ValueError(msg)
+
+        # Checks that the number of opening and closing brackets are equal
+        if formula.count('(') != formula.count(')'):
+            msg = 'Number of opening and closing brackets is not equal ' \
+                  'in the input formula {}.'.format(formula)
+            raise ValueError(msg)
+
+        # Checks that every part of the original formula has been tokenized
+        for row in tokens:
+            for token in row:
+                formula = formula.replace(token, '', 1)
+        if len(formula) != 0:
+            msg = 'Part of formula was not successfully parsed as an ' \
+                  'element symbol, bracket or integer. {} was not parsed.' \
+                  .format(formula)
+            raise ValueError(msg)
+
+        # Works through the tokens building a stack
+        mat_stack = [Counter()]
+        for symbol, multi1, opening_bracket, closing_bracket, multi2 in tokens:
+            if symbol:
+                mat_stack[-1][symbol] += int(multi1 or 1)
+            if opening_bracket:
+                mat_stack.append(Counter())
+            if closing_bracket:
+                stack_top = mat_stack.pop()
+                for symbol, value in stack_top.items():
+                    mat_stack[-1][symbol] += int(multi2 or 1) * value
+
+        # Normalizing percentages
+        percents = mat_stack[0].values()
+        norm_percents = [float(i) / sum(percents) for i in percents]
+        elements = mat_stack[0].keys()
+
+        # Adds each element and percent to the material
+        for element, percent in zip(elements, norm_percents):
+            if enrichment_target is not None and element == re.sub(r'\d+$', '', enrichment_target):
+                self.add_element(element, percent, percent_type, enrichment,
+                                 enrichment_target, enrichment_type)
+            else:
+                self.add_element(element, percent, percent_type)
+
     def add_s_alpha_beta(self, name, fraction=1.0):
         r"""Add an :math:`S(\alpha,\beta)` table to the material
 
diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py
index 08c60d3ed1..ba3be4e37a 100644
--- a/tests/unit_tests/test_deplete_chain.py
+++ b/tests/unit_tests/test_deplete_chain.py
@@ -6,7 +6,6 @@ from pathlib import Path
 from itertools import product
 
 import numpy as np
-from openmc.data import zam, ATOMIC_SYMBOL
 from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram
 import pytest
 
@@ -405,7 +404,8 @@ def test_fission_yield_attribute(simple_chain):
     dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1)
     with pytest.raises(
             ValueError, match="fission yield.*not equal.*compositions"):
-       cram.deplete(empty_chain, dummy_conc, None, 0.5)
+        cram.deplete(empty_chain, dummy_conc, None, 0.5)
+
 
 def test_validate(simple_chain):
     """Test the validate method"""
@@ -457,3 +457,95 @@ def test_validate_inputs():
 
     with pytest.raises(ValueError, match="tolerance"):
         c.validate(tolerance=-1)
+
+
+@pytest.fixture
+def gnd_simple_chain():
+    chainfile = Path(__file__).parents[1] / "chain_simple.xml"
+    return Chain.from_xml(chainfile)
+
+
+def test_reduce(gnd_simple_chain):
+    ref_U5 = gnd_simple_chain["U235"]
+    ref_iodine = gnd_simple_chain["I135"]
+    ref_U5_yields = ref_U5.yield_data
+
+    no_depth = gnd_simple_chain.reduce(["U235", "I135"], 0)
+    # We should get a chain just containing U235 and I135
+    assert len(no_depth) == 2
+    assert set(no_depth.reactions) == set(gnd_simple_chain.reactions)
+
+    u5_round0 = no_depth["U235"]
+    assert u5_round0.n_decay_modes == ref_U5.n_decay_modes
+    for newmode, refmode in zip(u5_round0.decay_modes, ref_U5.decay_modes):
+        assert newmode.target is None
+        assert newmode.type == refmode.type
+        assert newmode.branching_ratio == refmode.branching_ratio
+
+    assert u5_round0.n_reaction_paths == ref_U5.n_reaction_paths
+    for newrxn, refrxn in zip(u5_round0.reactions, ref_U5.reactions):
+        assert newrxn.target is None
+        assert newrxn.type == refrxn.type
+        assert newrxn.Q == refrxn.Q
+        assert newrxn.branching_ratio == refrxn.branching_ratio
+
+    assert u5_round0.yield_data is not None
+    assert u5_round0.yield_data.products == ("I135",)
+    assert u5_round0.yield_data.yield_matrix == (
+            ref_U5_yields.yield_matrix[:, ref_U5_yields.products.index("I135")]
+    )
+
+    bareI5 = no_depth["I135"]
+    assert bareI5.n_decay_modes == ref_iodine.n_decay_modes
+    for newmode, refmode in zip(bareI5.decay_modes, ref_iodine.decay_modes):
+        assert newmode.target is None
+        assert newmode.type == refmode.type
+        assert newmode.branching_ratio == refmode.branching_ratio
+
+    assert bareI5.n_reaction_paths == ref_iodine.n_reaction_paths
+    for newrxn, refrxn in zip(bareI5.reactions, ref_iodine.reactions):
+        assert newrxn.target is None
+        assert newrxn.type == refrxn.type
+        assert newrxn.Q == refrxn.Q
+        assert newrxn.branching_ratio == refrxn.branching_ratio
+
+    follow_u5 = gnd_simple_chain.reduce(["U235"], 1)
+    u5_round1 = follow_u5["U235"]
+    assert u5_round1.decay_modes == ref_U5.decay_modes
+    assert u5_round1.reactions == ref_U5.reactions
+    assert u5_round1.yield_data is not None
+    assert (
+        u5_round1.yield_data.yield_matrix == ref_U5_yields.yield_matrix
+    ).all()
+
+    # Per the chain_simple.xml
+    # I135 -> Xe135 -> Cs135
+    # I135 -> Xe136
+    # No limit on depth
+    iodine_chain = gnd_simple_chain.reduce(["I135"])
+    truncated_iodine = gnd_simple_chain.reduce(["I135"], 1)
+    assert len(iodine_chain) == 4
+    assert len(truncated_iodine) == 3
+    assert set(iodine_chain.nuclide_dict) == {
+        "I135", "Xe135", "Xe136", "Cs135"}
+    assert set(truncated_iodine.nuclide_dict) == {"I135", "Xe135", "Xe136"}
+    assert iodine_chain.reactions == ["(n,gamma)"]
+    assert iodine_chain["I135"].decay_modes == ref_iodine.decay_modes
+    assert iodine_chain["I135"].reactions == ref_iodine.reactions
+    for mode in truncated_iodine["Xe135"].decay_modes:
+        assert mode.target is None
+
+    # Test that no FissionYieldDistribution is made if there are no
+    # fission products
+    u5_noyields = gnd_simple_chain.reduce(["U235"], 0)["U235"]
+    assert u5_noyields.yield_data is None
+
+    # Check early termination if the eventual full chain
+    # is specified by using the iodine isotopes
+    new_iodine = gnd_simple_chain.reduce(set(iodine_chain.nuclide_dict))
+    assert set(iodine_chain.nuclide_dict) == set(new_iodine.nuclide_dict)
+
+    # Failure if some requested isotopes not in chain
+
+    with pytest.raises(IndexError, match=".*not found.*Xx999"):
+        gnd_simple_chain.reduce(["U235", "Xx999"])
diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py
index b8a258df00..126206cc70 100644
--- a/tests/unit_tests/test_deplete_nuclide.py
+++ b/tests/unit_tests/test_deplete_nuclide.py
@@ -176,6 +176,21 @@ def test_fission_yield_distribution():
     with pytest.raises(TypeError):
         orig_yields *= similar
 
+    # Test restriction of fission products
+    strict_restrict = yield_dist.restrict_products(["Xe135", "Sm149"])
+    with_extras = yield_dist.restrict_products(
+        ["Xe135", "Sm149", "H1", "U235"])
+
+    assert strict_restrict.products == ("Sm149", "Xe135")
+    assert strict_restrict.energies == yield_dist.energies
+    assert with_extras.products == ("Sm149", "Xe135")
+    assert with_extras.energies == yield_dist.energies
+    for ene, new_yields in strict_restrict.items():
+        for product in strict_restrict.products:
+            assert new_yields[product] == yield_dist[ene][product]
+            assert with_extras[ene][product] == yield_dist[ene][product]
+
+    assert yield_dist.restrict_products(["U235"]) is None
 
 def test_validate():
 
diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py
index c3ceda5d56..9b40a1ac91 100644
--- a/tests/unit_tests/test_deplete_resultslist.py
+++ b/tests/unit_tests/test_deplete_resultslist.py
@@ -19,12 +19,28 @@ def test_get_atoms(res):
     """Tests evaluating single nuclide concentration."""
     t, n = res.get_atoms("1", "Xe135")
 
-    t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
-    n_ref = [6.67473282e+08, 3.76986925e+14, 3.68587383e+14, 3.91338675e+14]
+    t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0])
+    n_ref = np.array(
+        [6.67473282e+08, 3.76986925e+14, 3.68587383e+14, 3.91338675e+14])
 
     np.testing.assert_allclose(t, t_ref)
     np.testing.assert_allclose(n, n_ref)
 
+    # Check alternate units
+    volume = res[0].volume["1"]
+
+    t_days, n_cm3 = res.get_atoms("1", "Xe135", nuc_units="atom/cm3", time_units="d")
+
+    assert t_days == pytest.approx(t_ref / (60 * 60 * 24))
+    assert n_cm3 == pytest.approx(n_ref / volume)
+
+    t_min, n_bcm = res.get_atoms("1", "Xe135", nuc_units="atom/b-cm", time_units="min")
+    assert n_bcm == pytest.approx(n_cm3 * 1e-24)
+    assert t_min == pytest.approx(t_ref / 60)
+
+    t_hour, _n = res.get_atoms("1", "Xe135", time_units="h")
+    assert t_hour == pytest.approx(t_ref / (60 * 60))
+
 
 def test_get_reaction_rate(res):
     """Tests evaluating reaction rate."""
diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py
index 65fdb26239..d62691bb1e 100644
--- a/tests/unit_tests/test_material.py
+++ b/tests/unit_tests/test_material.py
@@ -1,8 +1,11 @@
+from collections import defaultdict
+
+import pytest
+
 import openmc
+import openmc.examples
 import openmc.model
 import openmc.stats
-import openmc.examples
-import pytest
 
 
 def test_attributes(uo2):
@@ -69,6 +72,97 @@ def test_elements_by_name():
     assert a._nuclides == b._nuclides
     assert b._nuclides == c._nuclides
 
+
+def test_add_elements_by_formula():
+    """Test adding elements from a formula"""
+    # testing the correct nuclides and elements are added to a material
+    m = openmc.Material()
+    m.add_elements_from_formula('Li4SiO4')
+    # checking the ratio of elements is 4:1:4 for Li:Si:O
+    elem = defaultdict(float)
+    for nuclide, adens in m.get_nuclide_atom_densities().values():
+        if nuclide.startswith("Li"):
+            elem["Li"] += adens
+        if nuclide.startswith("Si"):
+            elem["Si"] += adens
+        if nuclide.startswith("O"):
+            elem["O"] += adens
+    total_number_of_atoms = 9
+    assert elem["Li"] == pytest.approx(4./total_number_of_atoms)
+    assert elem["Si"] == pytest.approx(1./total_number_of_atoms)
+    assert elem["O"] == pytest.approx(4/total_number_of_atoms)
+    # testing the correct nuclides are added to the Material
+    ref_dens = {'Li6': 0.033728, 'Li7': 0.410715,
+                'Si28': 0.102477, 'Si29': 0.0052035, 'Si30': 0.0034301,
+                'O16': 0.443386, 'O17': 0.000168}
+    nuc_dens = m.get_nuclide_atom_densities()
+    for nuclide in ref_dens:
+        assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
+
+    # testing the correct nuclides are added to the Material when enriched
+    m = openmc.Material()
+    m.add_elements_from_formula('Li4SiO4',
+                                enrichment=60.,
+                                enrichment_target='Li6')
+    ref_dens = {'Li6': 0.2666, 'Li7': 0.1777,
+                'Si28': 0.102477, 'Si29': 0.0052035, 'Si30': 0.0034301,
+                'O16': 0.443386, 'O17': 0.000168}
+    nuc_dens = m.get_nuclide_atom_densities()
+    for nuclide in ref_dens:
+        assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
+
+    # testing the use of brackets
+    m = openmc.Material()
+    m.add_elements_from_formula('Mg2(NO3)2')
+
+    # checking the ratio of elements is 2:2:6 for Mg:N:O
+    elem = defaultdict(float)
+    for nuclide, adens in m.get_nuclide_atom_densities().values():
+        if nuclide.startswith("Mg"):
+            elem["Mg"] += adens
+        if nuclide.startswith("N"):
+            elem["N"] += adens
+        if nuclide.startswith("O"):
+            elem["O"] += adens
+    total_number_of_atoms = 10
+    assert elem["Mg"] == pytest.approx(2./total_number_of_atoms)
+    assert elem["N"] == pytest.approx(2./total_number_of_atoms)
+    assert elem["O"] == pytest.approx(6/total_number_of_atoms)
+
+    # testing the correct nuclides are added when brackets are used
+    ref_dens = {'Mg24': 0.157902, 'Mg25': 0.02004, 'Mg26': 0.022058,
+                'N14': 0.199267, 'N15': 0.000732,
+                'O16': 0.599772, 'O17': 0.000227}
+    nuc_dens = m.get_nuclide_atom_densities()
+    for nuclide in ref_dens:
+        assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
+
+    # testing non integer multiplier results in a value error
+    m = openmc.Material()
+    with pytest.raises(ValueError):
+        m.add_elements_from_formula('Li4.2SiO4')
+
+    # testing lowercase elements results in a value error
+    m = openmc.Material()
+    with pytest.raises(ValueError):
+        m.add_elements_from_formula('li4SiO4')
+
+    # testing lowercase elements results in a value error
+    m = openmc.Material()
+    with pytest.raises(ValueError):
+        m.add_elements_from_formula('Li4Sio4')
+
+    # testing incorrect character in formula results in a value error
+    m = openmc.Material()
+    with pytest.raises(ValueError):
+        m.add_elements_from_formula('Li4$SiO4')
+
+    # testing unequal opening and closing brackets
+    m = openmc.Material()
+    with pytest.raises(ValueError):
+        m.add_elements_from_formula('Fe(H2O)4(OH)2)')
+
+
 def test_density():
     m = openmc.Material()
     for unit in ['g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3']:
@@ -275,4 +369,3 @@ def test_mix_materials():
     assert m3.density == pytest.approx(dens3)
     assert m4.density == pytest.approx(dens4)
     assert m5.density == pytest.approx(dens5)
-