diff --git a/openmc/deplete/integrator/abc.py b/openmc/deplete/integrator/abc.py
index cb4ac13a25..94ed2a86c3 100644
--- a/openmc/deplete/integrator/abc.py
+++ b/openmc/deplete/integrator/abc.py
@@ -4,6 +4,7 @@ from collections.abc import Iterable
 
 from uncertainties import ufloat
 
+from openmc.capi import statepoint_write
 from openmc.deplete import Results, OperatorResult
 
 
@@ -111,6 +112,7 @@ class Integrator(ABC):
     def _get_bos_data_from_openmc(self, step_index, step_power, bos_conc):
         x = deepcopy(bos_conc)
         res = self.operator(x, step_power)
+        self._write_statepoint(step_index)
         return x, res
 
     def _get_bos_data_from_restart(self, step_index, step_power, bos_conc):
@@ -160,17 +162,21 @@ class Integrator(ABC):
             res_list = [self.operator(conc, p)]
             self._save_results(
                 [conc], res_list, [t, t], p, self._ires + len(self))
+            self._write_statepoint(len(self))
 
     def _save_results(self, conc_list, results_list, time_list, power,
                       index, proc_time=None):
-        """Save the results at the end of of one step
-
-        Abstracted to support the predictor's unique save location
-        """
+        """Save the results at the end of of one step"""
         Results.save(
             self.operator, conc_list, results_list, time_list,
             power, index, proc_time)
 
+    def _write_statepoint(self, step_index):
+        """Use capi to write a statepoint for this index"""
+        statepoint_write(
+            "openmc_simulation_n{}.h5".format(step_index + self._ires),
+            write_source=False)
+
 
 class SI_Integrator(Integrator):
     """Abstract for the Stochastic Implicit Euler integrators
@@ -252,3 +258,4 @@ class SI_Integrator(Integrator):
             # No final simulation for SIE, use last iteration results
             self._save_results(
                 [conc], [res_list[-1]], [t, t], p, self._ires + len(self))
+            self._write_statepoint(len(self))
diff --git a/tests/regression_tests/deplete/test.py b/tests/regression_tests/deplete/test.py
index 5f2d8cc454..59eca4a545 100644
--- a/tests/regression_tests/deplete/test.py
+++ b/tests/regression_tests/deplete/test.py
@@ -39,7 +39,7 @@ def test_full(run_in_tmpdir):
     space = openmc.stats.Box(lower_left, upper_right)
     settings.source = openmc.Source(space=space)
     settings.seed = 1
-    settings.verbosity = 3
+    settings.verbosity = 1
 
     # Create operator
     chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
@@ -101,3 +101,18 @@ def test_full(run_in_tmpdir):
 
             assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format(
                 mat, nuc, y_old, y_test)
+
+    # Compare statepoint files with depletion results
+
+    t_test, k_test = res_test.get_eigenvalue()
+    t_ref, k_ref = res_ref.get_eigenvalue()
+    k_state = np.empty_like(k_ref)
+
+    # Get statepoint files for all BOS points and EOL
+    for n in range(N + 1):
+        statepoint = openmc.StatePoint("openmc_simulation_n{}.h5".format(n))
+        k_n = statepoint.k_combined
+        k_state[n] = [k_n.nominal_value, k_n.std_dev]
+    # Look for exact match pulling from statepoint and depletion_results
+    assert np.all(k_state == k_test)
+    assert np.allclose(k_test, k_ref)