diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py
index 3c7a13687..0d1e8ad41 100644
--- a/openmc/deplete/abc.py
+++ b/openmc/deplete/abc.py
@@ -851,7 +851,7 @@ class Integrator(ABC):
         return (self.operator.prev_res[-1].time[-1],
                 len(self.operator.prev_res) - 1)
 
-    def integrate(self, final_step=True):
+    def integrate(self, final_step=True, output=True):
         """Perform the entire depletion process across all steps
 
         Parameters
@@ -861,12 +861,18 @@ class Integrator(ABC):
             of the last timestep.
 
             .. versionadded:: 0.12.1
+        output : bool, optional
+            Indicate whether to display information about progress
 
+            .. versionadded:: 0.13.1
         """
         with self.operator as conc:
             t, self._i_res = self._get_start_data()
 
             for i, (dt, source_rate) in enumerate(self):
+                if output:
+                    print(f"[openmc.deplete] t={t} s, dt={dt} s, source={source_rate}")
+
                 # Solve transport equation (or obtain result from restart)
                 if i > 0 or self.operator.prev_res is None:
                     conc, res = self._get_bos_data_from_operator(i, source_rate, conc)
@@ -892,6 +898,8 @@ class Integrator(ABC):
             # source rate is passed to the transport operator (which knows to
             # just return zero reaction rates without actually doing a transport
             # solve)
+            if output and final_step:
+                print(f"[openmc.deplete] t={t} (final operator evaluation)")
             res_list = [self.operator(conc, source_rate if final_step else 0.0)]
             Results.save(self.operator, [conc], res_list, [t, t],
                          source_rate, self._i_res + len(self), proc_time)
@@ -1004,12 +1012,23 @@ class SIIntegrator(Integrator):
             self.operator.settings.particles //= self.n_steps
         return inherited
 
-    def integrate(self):
-        """Perform the entire depletion process across all steps"""
+    def integrate(self, output=True):
+        """Perform the entire depletion process across all steps
+
+        Parameters
+        ----------
+        output : bool, optional
+            Indicate whether to display information about progress
+
+            .. versionadded:: 0.13.1
+        """
         with self.operator as conc:
             t, self._i_res = self._get_start_data()
 
             for i, (dt, p) in enumerate(self):
+                if output:
+                    print(f"[openmc.deplete] t={t} s, dt={dt} s, source={p}")
+
                 if i == 0:
                     if self.operator.prev_res is None:
                         conc, res = self._get_bos_data_from_operator(i, p, conc)
diff --git a/openmc/deplete/atom_number.py b/openmc/deplete/atom_number.py
index e4982f872..5430ddc35 100644
--- a/openmc/deplete/atom_number.py
+++ b/openmc/deplete/atom_number.py
@@ -6,6 +6,8 @@ from collections import OrderedDict
 
 import numpy as np
 
+from openmc import Material
+
 
 class AtomNumber:
     """Stores local material compositions (atoms of each nuclide).
@@ -58,6 +60,12 @@ class AtomNumber:
 
         self.number = np.zeros((len(local_mats), len(nuclides)))
 
+    def _get_mat_index(self, mat):
+        """Helper method for getting material index"""
+        if isinstance(mat, Material):
+            mat = str(mat.id)
+        return self.index_mat[mat] if isinstance(mat, str) else mat
+
     def __getitem__(self, pos):
         """Retrieves total atom number from AtomNumber.
 
@@ -75,8 +83,7 @@ class AtomNumber:
         """
 
         mat, nuc = pos
-        if isinstance(mat, str):
-            mat = self.index_mat[mat]
+        mat = self._get_mat_index(mat)
         if isinstance(nuc, str):
             nuc = self.index_nuc[nuc]
 
@@ -96,8 +103,7 @@ class AtomNumber:
 
         """
         mat, nuc = pos
-        if isinstance(mat, str):
-            mat = self.index_mat[mat]
+        mat = self._get_mat_index(mat)
         if isinstance(nuc, str):
             nuc = self.index_nuc[nuc]
 
@@ -121,11 +127,11 @@ class AtomNumber:
                 if ind < self.n_nuc_burn]
 
     def get_atom_density(self, mat, nuc):
-        """Accesses atom density instead of total number.
+        """Return atom density of given material and nuclide
 
         Parameters
         ----------
-        mat : str, int or slice
+        mat : str, int, openmc.Material or slice
             Material index.
         nuc : str, int or slice
             Nuclide index.
@@ -136,19 +142,43 @@ class AtomNumber:
             Density in [atom/cm^3]
 
         """
-        if isinstance(mat, str):
-            mat = self.index_mat[mat]
+        mat = self._get_mat_index(mat)
         if isinstance(nuc, str):
             nuc = self.index_nuc[nuc]
 
         return self[mat, nuc] / self.volume[mat]
 
+    def get_atom_densities(self, mat, units='atom/b-cm'):
+        """Return atom densities for a given material
+
+        Parameters
+        ----------
+        mat : str, int, openmc.Material or slice
+            Material index.
+        units : {"atom/b-cm", "atom/cm3"}, optional
+            Units for the returned concentration. Default is ``"atom/b-cm"``
+
+            .. versionadded:: 0.13.1
+
+        Returns
+        -------
+        dict
+            Dictionary mapping nuclides to atom densities
+
+        """
+        mat = self._get_mat_index(mat)
+        normalization = (1.0e-24 if units == 'atom/b-cm' else 1.0) / self.volume[mat]
+        return {
+            name: normalization * self[mat, nuc]
+            for name, nuc in self.index_nuc.items()
+        }
+
     def set_atom_density(self, mat, nuc, val):
         """Sets atom density instead of total number.
 
         Parameters
         ----------
-        mat : str, int or slice
+        mat : str, int, openmc.Material or slice
             Material index.
         nuc : str, int or slice
             Nuclide index.
@@ -156,8 +186,7 @@ class AtomNumber:
             Array of densities to set in [atom/cm^3]
 
         """
-        if isinstance(mat, str):
-            mat = self.index_mat[mat]
+        mat = self._get_mat_index(mat)
         if isinstance(nuc, str):
             nuc = self.index_nuc[nuc]
 
@@ -168,7 +197,7 @@ class AtomNumber:
 
         Parameters
         ----------
-        mat : str, int or slice
+        mat : str, int, openmc.Material or slice
             Material index.
 
         Returns
@@ -177,9 +206,7 @@ class AtomNumber:
             The slice requested in [atom].
 
         """
-        if isinstance(mat, str):
-            mat = self.index_mat[mat]
-
+        mat = self._get_mat_index(mat)
         return self[mat, :self.n_nuc_burn]
 
     def set_mat_slice(self, mat, val):
@@ -187,15 +214,13 @@ class AtomNumber:
 
         Parameters
         ----------
-        mat : str, int or slice
+        mat : str, int, openmc.Material, or slice
             Material index.
         val : numpy.ndarray
             The slice to set in [atom]
 
         """
-        if isinstance(mat, str):
-            mat = self.index_mat[mat]
-
+        mat = self._get_mat_index(mat)
         self[mat, :self.n_nuc_burn] = val
 
     def set_density(self, total_density):
diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py
index 49d1d62d7..c19ef076a 100644
--- a/openmc/deplete/integrators.py
+++ b/openmc/deplete/integrators.py
@@ -124,13 +124,13 @@ class CF4Integrator(Integrator):
 
     .. math::
         \begin{aligned}
-        \mathbf{A}_1 &= h\mathbf{A}(\mathbf{n}_0) \\
-        \hat{\mathbf{n}}_1 &= \exp \left ( \frac{\mathbf{A}_1}{2} \right ) \\
+        \mathbf{A}_1 &= h\mathbf{A}(\mathbf{n}_i) \\
+        \hat{\mathbf{n}}_1 &= \exp \left ( \frac{\mathbf{A}_1}{2} \right ) \mathbf{n}_i \\
         \mathbf{A}_2 &= h\mathbf{A}(\hat{\mathbf{n}}_1) \\
-        \hat{\mathbf{n}}_2 &= \exp \left ( \frac{\mathbf{A}_2}{2} \right ) \\
+        \hat{\mathbf{n}}_2 &= \exp \left ( \frac{\mathbf{A}_2}{2} \right ) \mathbf{n}_i \\
         \mathbf{A}_3 &= h \mathbf{A}(\hat{\mathbf{n}}_2) \\
         \hat{\mathbf{n}}_3 &= \exp \left ( -\frac{\mathbf{A}_1}{2} + \mathbf{A}_3
-            \right ) \\
+            \right ) \hat{\mathbf{n}}_1 \\
         \mathbf{A}_4 &= h\mathbf{A}(\hat{\mathbf{n}}_3) \\
         \mathbf{n}_{i+1} &= \exp \left ( \frac{\mathbf{A}_1}{4} + \frac{\mathbf{A}_2}{6}
             + \frac{\mathbf{A}_3}{6} - \frac{\mathbf{A}_4}{12} \right )
@@ -207,7 +207,8 @@ class CELIIntegrator(Integrator):
 
     .. math::
         \begin{aligned}
-        \mathbf{n}_{i+1}^p &= \exp \left ( h \mathbf{A}(\mathbf{n}_i ) \right ) \\
+        \mathbf{n}_{i+1}^p &= \exp \left ( h \mathbf{A}(\mathbf{n}_i ) \right )
+        \mathbf{n}_i \\
         \mathbf{n}_{i+1} &= \exp \left( \frac{h}{12} \mathbf{A}(\mathbf{n}_i) +
             \frac{5h}{12} \mathbf{A}(\mathbf{n}_{i+1}^p) \right)
         \exp \left( \frac{5h}{12} \mathbf{A}(\mathbf{n}_i) +
@@ -268,12 +269,12 @@ class EPCRK4Integrator(Integrator):
 
     .. math::
         \begin{aligned}
-        \mathbf{A}_1 &= h\mathbf{A}(\mathbf{n}_0) \\
-        \hat{\mathbf{n}}_1 &= \exp \left ( \frac{\mathbf{A}_1}{2} \right ) \\
+        \mathbf{A}_1 &= h\mathbf{A}(\mathbf{n}_i) \\
+        \hat{\mathbf{n}}_1 &= \exp \left ( \frac{\mathbf{A}_1}{2} \right ) \mathbf{n}_i \\
         \mathbf{A}_2 &= h\mathbf{A}(\hat{\mathbf{n}}_1) \\
-        \hat{\mathbf{n}}_2 &= \exp \left ( \frac{\mathbf{A}_2}{2} \right ) \\
+        \hat{\mathbf{n}}_2 &= \exp \left ( \frac{\mathbf{A}_2}{2} \right ) \mathbf{n}_i \\
         \mathbf{A}_3 &= h \mathbf{A}(\hat{\mathbf{n}}_2) \\
-        \hat{\mathbf{n}}_3 &= \exp \left ( \mathbf{A}_3 \right ) \\
+        \hat{\mathbf{n}}_3 &= \exp \left ( \mathbf{A}_3 \right ) \mathbf{n}_i \\
         \mathbf{A}_4 &= h\mathbf{A}(\hat{\mathbf{n}}_3) \\
         \mathbf{n}_{i+1} &= \exp \left ( \frac{\mathbf{A}_1}{6} + \frac{\mathbf{A}_2}{3}
         + \frac{\mathbf{A}_3}{3} + \frac{\mathbf{A}_4}{6} \right ) \mathbf{n}_i.
diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py
index 992248246..7b6ddd788 100644
--- a/openmc/deplete/nuclide.py
+++ b/openmc/deplete/nuclide.py
@@ -123,6 +123,10 @@ class Nuclide:
         # Neutron fission yields, if present
         self._yield_data = None
 
+    def __repr__(self):
+        n_modes, n_rx = self.n_decay_modes, self.n_reaction_paths
+        return f"<Nuclide: {self.name} ({n_modes} modes, {n_rx} reactions)>"
+
     @property
     def n_decay_modes(self):
         return len(self.decay_modes)
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index b5e377bed..337e8f8c0 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -84,11 +84,11 @@ class Operator(TransportOperator):
 
     .. versionchanged:: 0.13.0
         The geometry and settings parameters have been replaced with a
-        model parameter that takes an openmc.Model object
+        model parameter that takes a :class:`~openmc.model.Model` object
 
     Parameters
     ----------
-    model : openmc.Model
+    model : openmc.model.Model
         OpenMC model object
     chain_file : str, optional
         Path to the depletion chain XML file.  Defaults to the file
@@ -163,7 +163,7 @@ class Operator(TransportOperator):
 
     Attributes
     ----------
-    model : openmc.Model
+    model : openmc.model.Model
         OpenMC model object
     geometry : openmc.Geometry
         OpenMC geometry object
@@ -195,8 +195,6 @@ class Operator(TransportOperator):
     prev_res : ResultsList or None
         Results from a previous depletion calculation. ``None`` if no
         results are to be used.
-    diff_burnable_mats : bool
-        Whether to differentiate burnable materials with multiple instances
     cleanup_when_done : bool
         Whether to finalize and clear the shared library memory when the
         depletion operation is complete. Defaults to clearing the library.
@@ -248,7 +246,6 @@ class Operator(TransportOperator):
             )
         self.materials = model.materials
 
-        self.diff_burnable_mats = diff_burnable_mats
         self.cleanup_when_done = True
 
         # Reduce the chain before we create more materials
@@ -262,8 +259,11 @@ class Operator(TransportOperator):
             self.chain = self.chain.reduce(all_isotopes, reduce_chain_level)
 
         # Differentiate burnable materials with multiple instances
-        if self.diff_burnable_mats:
+        if diff_burnable_mats:
             self._differentiate_burnable_mats()
+            self.materials = openmc.Materials(
+                model.geometry.get_all_materials().values()
+            )
 
         # Clear out OpenMC, create task lists, distribute
         openmc.reset_auto_ids()
diff --git a/src/output.cpp b/src/output.cpp
index 58f081f67..6565a9e67 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -545,7 +545,7 @@ void print_results()
     fmt::print(" Leakage Fraction           = {:.5f}\n",
       gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n);
   }
-  fmt::print("\n");
+  std::cout << std::endl;
 }
 
 //==============================================================================