diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index 1f16d9cafa..865f533ceb 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -129,6 +129,7 @@ class Chain(object):
         self.nuclides = []
         self.reactions = []
         self.nuclide_dict = OrderedDict()
+        self._default_fsn_yields = None
 
     def __contains__(self, nuclide):
         return nuclide in self.nuclide_dict
@@ -374,13 +375,27 @@ class Chain(object):
             clean_indentation(root_elem)
             tree.write(str(filename), encoding='utf-8')
 
-    def form_matrix(self, rates):
+    def _build_default_yields(self):
+        """Return dictionary {str: {str: float}}"""
+        # Take lowest energy for back compatability
+        # Should be removed by end of this feature
+        out = {}
+        for nuc in self.nuclides:
+            if len(nuc.yield_data) == 0:
+                continue
+            _energy, yield_data = sorted(nuc.yield_data.items())[0]
+            out[nuc.name] = {prod: frac for prod, frac in yield_data}
+        return out
+
+    def form_matrix(self, rates, fission_yields=None):
         """Forms depletion matrix.
 
         Parameters
         ----------
         rates : numpy.ndarray
             2D array indexed by (nuclide, reaction)
+        fission_yields : dictionary of tuple to float, optional
+            Option to use a custom set of fission yields
 
         Returns
         -------
@@ -391,6 +406,11 @@ class Chain(object):
         matrix = defaultdict(float)
         reactions = set()
 
+        if fission_yields is None:
+            if self._default_fsn_yields is None:
+                self._default_fsn_yields = self._build_default_yields()
+            fission_yields = self._default_fsn_yields
+
         for i, nuc in enumerate(self.nuclides):
 
             if nuc.n_decay_modes != 0:
@@ -438,8 +458,7 @@ class Chain(object):
                             # Assume that we should always use thermal fission
                             # yields. At some point it would be nice to account
                             # for the energy-dependence..
-                            energy, data = sorted(nuc.yield_data.items())[0]
-                            for product, y in data:
+                            for product, y in fission_yields[nuc.name].items():
                                 yield_val = y * path_rate
                                 if yield_val != 0.0:
                                     k = self.nuclide_dict[product]
diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py
index 33b13b596d..df82eab1b0 100644
--- a/tests/unit_tests/test_deplete_chain.py
+++ b/tests/unit_tests/test_deplete_chain.py
@@ -3,6 +3,7 @@
 from collections.abc import Mapping
 import os
 from pathlib import Path
+from itertools import product
 
 import numpy as np
 from openmc.data import zam, ATOMIC_SYMBOL
@@ -220,6 +221,12 @@ def test_form_matrix(simple_chain):
     assert mat[1, 2] == mat12
     assert mat[2, 2] == mat22
 
+    # Pass equivalent fission yields directly
+    # Ensure identical matrix is formed
+    f_yields = {"C": {"A": 0.0292737, "B": 0.002566345}}
+    new_mat = chain.form_matrix(react[0], f_yields)
+    for r, c in product(range(3), range(3)):
+        assert new_mat[r, c] == mat[r, c]
 
 def test_getitem():
     """ Test nuc_by_ind converter function. """