diff --git a/Dockerfile b/Dockerfile
index 755d9f2a06..341390c9ae 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -15,7 +15,7 @@
 
 # sudo docker run image_name:tag_name or ID with no tag sudo docker run ID number
 
-FROM debian:bullseye-slim
+FROM debian:bullseye-slim AS dependencies
 
 # By default this Dockerfile builds OpenMC without DAGMC and LIBMESH support
 ARG build_dagmc=off
@@ -172,6 +172,8 @@ RUN if [ "$build_libmesh" = "on" ]; then \
         && rm -rf ${LIBMESH_INSTALL_DIR}/build ${LIBMESH_INSTALL_DIR}/libmesh ; \
     fi
 
+FROM dependencies AS build
+
 # clone and install openmc
 RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
     && git clone --shallow-submodules --recurse-submodules --single-branch -b ${openmc_branch} --depth=1 ${OPENMC_REPO} \
@@ -207,5 +209,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
     && cd ../openmc && pip install .[test,depletion-mpi] \
     && python -c "import openmc"
 
+FROM build AS release
+
 # Download cross sections (NNDC and WMP) and ENDF data needed by test suite
 RUN ${HOME}/OpenMC/openmc/tools/ci/download-xs.sh
diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst
index 6f5c006c24..287738774c 100644
--- a/docs/source/pythonapi/data.rst
+++ b/docs/source/pythonapi/data.rst
@@ -61,6 +61,7 @@ Core Functions
 
     atomic_mass
     atomic_weight
+    decay_constant
     dose_coefficients
     gnd_name
     half_life
diff --git a/include/openmc/photon.h b/include/openmc/photon.h
index e8ec82d107..09fb3ba01b 100644
--- a/include/openmc/photon.h
+++ b/include/openmc/photon.h
@@ -36,7 +36,6 @@ public:
   int threshold;
   double n_electrons;
   double binding_energy;
-  xt::xtensor<double, 1> cross_section;
   vector<Transition> transitions;
 };
 
@@ -80,8 +79,11 @@ public:
   Tabulated1D coherent_anomalous_real_;
   Tabulated1D coherent_anomalous_imag_;
 
-  // Photoionization and atomic relaxation data
+  // Photoionization and atomic relaxation data. Subshell cross sections are
+  // stored separately to improve memory access pattern when calculating the
+  // total cross section
   vector<ElectronSubshell> shells_;
+  xt::xtensor<double, 2> cross_sections_;
 
   // Compton profile data
   xt::xtensor<double, 2> profile_pdf_;
diff --git a/openmc/cell.py b/openmc/cell.py
index 103a34fe63..4b419e1c6e 100644
--- a/openmc/cell.py
+++ b/openmc/cell.py
@@ -214,8 +214,8 @@ class Cell(IDManagerMixin):
                 self._atoms = self._fill.get_nuclide_atom_densities()
 
                 # Convert to total number of atoms
-                for key, nuclide in self._atoms.items():
-                    atom = nuclide[1] * self._volume * 1.0e+24
+                for key, atom_per_bcm in self._atoms.items():
+                    atom = atom_per_bcm * self._volume * 1.0e+24
                     self._atoms[key] = atom
 
             elif self.fill_type == 'distribmat':
@@ -224,12 +224,12 @@ class Cell(IDManagerMixin):
                 partial_volume = self.volume / len(self.fill)
                 self._atoms = OrderedDict()
                 for mat in self.fill:
-                    for key, nuclide in mat.get_nuclide_atom_densities().items():
+                    for key, atom_per_bcm in mat.get_nuclide_atom_densities().items():
                         # To account for overlap of nuclides between distribmat
                         # we need to append new atoms to any existing value
                         # hence it is necessary to ask for default.
                         atom = self._atoms.setdefault(key, 0)
-                        atom += nuclide[1] * partial_volume * 1.0e+24
+                        atom += atom_per_bcm * partial_volume * 1.0e+24
                         self._atoms[key] = atom
 
             else:
diff --git a/openmc/data/data.py b/openmc/data/data.py
index 3a3cf58268..71b93293e7 100644
--- a/openmc/data/data.py
+++ b/openmc/data/data.py
@@ -3,7 +3,7 @@ import json
 import os
 import re
 from pathlib import Path
-from math import sqrt
+from math import sqrt, log
 from warnings import warn
 
 # Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions
@@ -202,7 +202,7 @@ _GND_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)')
 
 # Used in half_life function as a cache
 _HALF_LIFE = {}
-
+_LOG_TWO = log(2.0)
 
 def atomic_mass(isotope):
     """Return atomic mass of isotope in atomic mass units.
@@ -283,6 +283,8 @@ def half_life(isotope):
     Half-life values are from the `ENDF/B-VIII.0 decay sublibrary
     <https://www.nndc.bnl.gov/endf-b8.0/download.html>`_.
 
+    .. versionadded:: 0.13.1
+
     Parameters
     ----------
     isotope : str
@@ -302,6 +304,35 @@ def half_life(isotope):
 
     return _HALF_LIFE.get(isotope.lower())
 
+
+def decay_constant(isotope):
+    """Return decay constant of isotope in [s^-1]
+
+    Decay constants are based on half-life values from the
+    :func:`~openmc.data.half_life` function. When the isotope is stable, a decay
+    constant of zero is returned.
+
+    .. versionadded:: 0.13.1
+
+    Parameters
+    ----------
+    isotope : str
+        Name of isotope, e.g., 'Pu239'
+
+    Returns
+    -------
+    float
+        Decay constant of isotope in [s^-1]
+
+    See also
+    --------
+    openmc.data.half_life
+
+    """
+    t = half_life(isotope)
+    return _LOG_TWO / t if t else 0.0
+
+
 def water_density(temperature, pressure=0.1013):
     """Return the density of liquid water at a given temperature and pressure.
 
diff --git a/openmc/data/decay.py b/openmc/data/decay.py
index 5ca2b235a5..51a594bc80 100644
--- a/openmc/data/decay.py
+++ b/openmc/data/decay.py
@@ -465,11 +465,10 @@ class Decay(EqualityMixin):
 
     @property
     def decay_constant(self):
-        if hasattr(self.half_life, 'n'):
-            return log(2.)/self.half_life
-        else:
-            mu, sigma = self.half_life
-            return ufloat(log(2.)/mu, log(2.)/mu**2*sigma)
+        if self.half_life.n == 0.0:
+            name = self.nuclide['name']
+            raise ValueError(f"{name} is listed as unstable but has a zero half-life.")
+        return log(2.)/self.half_life
 
     @property
     def decay_energy(self):
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index 327c762e08..7ef874e87f 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -528,9 +528,9 @@ class Operator(TransportOperator):
         """
         mat_id = str(mat.id)
 
-        for nuclide, density in mat.get_nuclide_atom_densities().values():
-            number = density * 1.0e24
-            self.number.set_atom_density(mat_id, nuclide, number)
+        for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items():
+            atom_per_cc = atom_per_bcm * 1.0e24
+            self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
 
     def _set_number_from_results(self, mat, prev_res):
         """Extracts material nuclides and number densities.
@@ -556,16 +556,15 @@ class Operator(TransportOperator):
         # Merge lists of nuclides, with the same order for every calculation
         geom_nuc_densities.update(depl_nuc)
 
-        for nuclide in geom_nuc_densities.keys():
+        for nuclide, atom_per_bcm in geom_nuc_densities.items():
             if nuclide in depl_nuc:
                 concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1]
                 volume = prev_res[-1].volume[mat_id]
-                number = concentration / volume
+                atom_per_cc = concentration / volume
             else:
-                density = geom_nuc_densities[nuclide][1]
-                number = density * 1.0e24
+                atom_per_cc = atom_per_bcm * 1.0e24
 
-            self.number.set_atom_density(mat_id, nuclide, number)
+            self.number.set_atom_density(mat_id, nuclide, atom_per_cc)
 
     def initial_condition(self):
         """Performs final setup and returns initial condition.
diff --git a/openmc/material.py b/openmc/material.py
index dcb33ff29f..a837b1f17d 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -4,7 +4,6 @@ from copy import deepcopy
 from numbers import Real
 from pathlib import Path
 import os
-import math
 import re
 import typing  # imported separately as py3.8 requires typing.Iterable
 import warnings
@@ -90,6 +89,9 @@ class Material(IDManagerMixin):
     fissionable_mass : float
         Mass of fissionable nuclides in the material in [g]. Requires that the
         :attr:`volume` attribute is set.
+    activity : float
+        Activity of the material in [Bq]. Requires that the :attr:`volume`
+        attribute is set.
 
     """
 
@@ -146,20 +148,10 @@ class Material(IDManagerMixin):
 
         return string
 
-
     @property
     def activity(self):
         """Returns the total activity of the material in Becquerels."""
-
-        atoms_per_barn_cm = self.get_nuclide_atom_densities()
-        total_activity = 0
-        for key, value in atoms_per_barn_cm.items():
-            half_life = openmc.data.half_life(key)
-            if half_life:
-                total_activity += value[1] / half_life
-        total_activity *= math.log(2) * 1e24 * self.volume
-
-        return total_activity
+        return sum(self.get_nuclide_activity().values())
 
     @property
     def name(self):
@@ -263,10 +255,10 @@ class Material(IDManagerMixin):
         if self.volume is None:
             raise ValueError("Volume must be set in order to determine mass.")
         density = 0.0
-        for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
+        for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items():
             Z = openmc.data.zam(nuc)[0]
             if Z >= 90:
-                density += 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
+                density += 1e24 * atoms_per_bcm * openmc.data.atomic_mass(nuc) \
                            / openmc.data.AVOGADRO
         return density*self.volume
 
@@ -794,11 +786,15 @@ class Material(IDManagerMixin):
         """Returns all nuclides in the material and their atomic densities in
         units of atom/b-cm
 
+        .. versionchanged:: 0.13.1
+            The values in the dictionary were changed from a tuple containing
+            the nuclide name and the density to just the density.
+
         Returns
         -------
         nuclides : dict
-            Dictionary whose keys are nuclide names and values are tuples of
-            (nuclide, density in atom/b-cm)
+            Dictionary whose keys are nuclide names and values are densities in
+            [atom/b-cm]
 
         """
 
@@ -858,10 +854,46 @@ class Material(IDManagerMixin):
 
         nuclides = OrderedDict()
         for n, nuc in enumerate(nucs):
-            nuclides[nuc] = (nuc, nuc_densities[n])
+            nuclides[nuc] = nuc_densities[n]
 
         return nuclides
 
+    def get_nuclide_activity(self):
+        """Return activity in [Bq] for each nuclide in the material
+
+        .. versionadded:: 0.13.1
+
+        Returns
+        -------
+        dict
+            Dictionary whose keys are nuclide names and values are activity in
+            [Bq].
+        """
+        activity = {}
+        for nuclide, atoms in self.get_nuclide_atoms().items():
+            inv_seconds = openmc.data.decay_constant(nuclide)
+            activity[nuclide] = inv_seconds * atoms
+        return activity
+
+    def get_nuclide_atoms(self):
+        """Return number of atoms of each nuclide in the material
+
+        .. versionadded:: 0.13.1
+
+        Returns
+        -------
+        dict
+            Dictionary whose keys are nuclide names and values are number of
+            atoms present in the material.
+
+        """
+        if self.volume is None:
+            raise ValueError("Volume must be set in order to determine atoms.")
+        atoms = {}
+        for nuclide, atom_per_bcm in self.get_nuclide_atom_densities().items():
+            atoms[nuclide] = 1.0e24 * atom_per_bcm * self.volume
+        return atoms
+
     def get_mass_density(self, nuclide: Optional[str] = None):
         """Return mass density of one or all nuclides
 
@@ -878,9 +910,9 @@ class Material(IDManagerMixin):
 
         """
         mass_density = 0.0
-        for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
+        for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items():
             if nuclide is None or nuclide == nuc:
-                density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
+                density_i = 1e24 * atoms_per_bcm * openmc.data.atomic_mass(nuc) \
                             / openmc.data.AVOGADRO
                 mass_density += density_i
         return mass_density
@@ -1101,7 +1133,7 @@ class Material(IDManagerMixin):
         nuclides_per_cc = defaultdict(float)
         mass_per_cc = defaultdict(float)
         for mat, wgt in zip(materials, wgts):
-            for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().values():
+            for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().items():
                 nuc_per_cc = wgt*1.e24*atoms_per_bcm
                 nuclides_per_cc[nuc] += nuc_per_cc
                 mass_per_cc[nuc] += nuc_per_cc*openmc.data.atomic_mass(nuc) / \
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 7f7b078795..d5a460baa8 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -193,6 +193,81 @@ class StructuredMesh(MeshBase):
         s1 = (slice(1, None),)*ndim + (slice(None),)
         return (vertices[s0] + vertices[s1]) / 2
 
+    @property
+    def num_mesh_cells(self):
+        return np.prod(self.dimension)
+
+    def write_data_to_vtk(self, points, filename, datasets, volume_normalization=True):
+        """Creates a VTK object of the mesh
+
+        Parameters
+        ----------
+        points : list or np.array
+            List of (X,Y,Z) tuples.
+        filename : str
+            Name of the VTK file to write.
+        datasets : dict
+            Dictionary whose keys are the data labels
+            and values are the data sets.
+        volume_normalization : bool, optional
+            Whether or not to normalize the data by
+            the volume of the mesh elements.
+
+        Raises
+        ------
+        RuntimeError
+            When the size of a dataset doesn't match the number of cells
+
+        Returns
+        -------
+        vtk.vtkStructuredGrid
+            the VTK object
+        """
+
+        import vtk
+        from vtk.util import numpy_support as nps
+
+        # check that the data sets are appropriately sized
+        errmsg = "The size of the dataset {} should be equal to the number of cells"
+        for label, dataset in datasets.items():
+            if isinstance(dataset, np.ndarray):
+                if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]:
+                    raise RuntimeError(errmsg.format(label))
+            else:
+                if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]:
+                    raise RuntimeError(errmsg.format(label))
+            cv.check_type('label', label, str)
+
+        vtk_grid = vtk.vtkStructuredGrid()
+
+        vtk_grid.SetDimensions(*self.dimension)
+
+        vtkPts = vtk.vtkPoints()
+        vtkPts.SetData(nps.numpy_to_vtk(points, deep=True))
+        vtk_grid.SetPoints(vtkPts)
+
+        # create VTK arrays for each of
+        # the data sets
+        for label, dataset in datasets.items():
+            dataset = np.asarray(dataset).flatten()
+
+            if volume_normalization:
+                dataset /= self.volumes.flatten()
+
+            dataset_array = vtk.vtkDoubleArray()
+            dataset_array.SetName(label)
+            dataset_array.SetArray(nps.numpy_to_vtk(dataset),
+                           dataset.size,
+                           True)
+            vtk_grid.GetCellData().AddArray(dataset_array)
+
+        # write the .vtk file
+        writer = vtk.vtkStructuredGridWriter()
+        writer.SetFileName(str(filename))
+        writer.SetInputData(vtk_grid)
+        writer.Write()
+
+        return vtk_grid
 
 class RegularMesh(StructuredMesh):
     """A regular Cartesian mesh in one, two, or three dimensions
@@ -273,10 +348,6 @@ class RegularMesh(StructuredMesh):
                 dims =  self._dimension
                 return [(u - l) / d for u, l, d in zip(us, ls, dims)]
 
-    @property
-    def num_mesh_cells(self):
-        return np.prod(self._dimension)
-
     @property
     def volumes(self):
         """Return Volumes for every mesh cell
@@ -621,7 +692,42 @@ class RegularMesh(StructuredMesh):
         root_cell.fill = lattice
 
         return root_cell, cells
+    
+    def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
+        """Creates a VTK object of the mesh
 
+        Parameters
+        ----------
+        filename : str or pathlib.Path
+            Name of the VTK file to write.
+        datasets : dict
+            Dictionary whose keys are the data labels
+            and values are the data sets.
+        volume_normalization : bool, optional
+            Whether or not to normalize the data by
+            the volume of the mesh elements.
+            Defaults to True.
+
+        Returns
+        -------
+        vtk.vtkStructuredGrid
+            the VTK object
+        """
+
+        ll, ur = self.lower_left, self.upper_right
+        x_vals = np.linspace(ll[0], ur[0], num=self.dimension[0] + 1)
+        y_vals = np.linspace(ll[1], ur[1], num=self.dimension[1] + 1)
+        z_vals = np.linspace(ll[2], ur[2], num=self.dimension[2] + 1)
+
+        # create points
+        pts_cartesian = np.array([[x, y, z] for z in z_vals for y in y_vals for x in x_vals])
+
+        return super().write_data_to_vtk(
+            points=pts_cartesian,
+            filename=filename,
+            datasets=datasets,
+            volume_normalization=volume_normalization
+        )
 
 def Mesh(*args, **kwargs):
     warnings.warn("Mesh has been renamed RegularMesh. Future versions of "
@@ -821,6 +927,36 @@ class RectilinearMesh(StructuredMesh):
 
         return element
 
+    def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
+        """Creates a VTK object of the mesh
+
+        Parameters
+        ----------
+        filename : str or pathlib.Path
+            Name of the VTK file to write.
+        datasets : dict
+            Dictionary whose keys are the data labels
+            and values are the data sets.
+        volume_normalization : bool, optional
+            Whether or not to normalize the data by
+            the volume of the mesh elements.
+            Defaults to True.
+
+        Returns
+        -------
+        vtk.vtkStructuredGrid
+            the VTK object
+        """
+        # create points
+        pts_cartesian = np.array([[x, y, z] for z in self.z_grid for y in self.y_grid for x in self.x_grid])
+
+        return super().write_data_to_vtk(
+            points=pts_cartesian,
+            filename=filename,
+            datasets=datasets,
+            volume_normalization=volume_normalization
+        )
+
 
 class CylindricalMesh(StructuredMesh):
     """A 3D cylindrical mesh
@@ -1012,6 +1148,46 @@ class CylindricalMesh(StructuredMesh):
 
         return np.multiply.outer(np.outer(V_r, V_p), V_z)
 
+    def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
+        """Creates a VTK object of the mesh
+
+        Parameters
+        ----------
+        filename : str or pathlib.Path
+            Name of the VTK file to write.
+        datasets : dict
+            Dictionary whose keys are the data labels
+            and values are the data sets.
+        volume_normalization : bool, optional
+            Whether or not to normalize the data by
+            the volume of the mesh elements.
+            Defaults to True.
+
+        Returns
+        -------
+        vtk.vtkStructuredGrid
+            the VTK object
+        """
+        # create points
+        pts_cylindrical = np.array(
+            [
+                [r, phi, z]
+                for z in self.z_grid
+                for phi in self.phi_grid
+                for r in self.r_grid
+            ]
+        )
+        pts_cartesian = np.copy(pts_cylindrical)
+        r, phi = pts_cylindrical[:, 0], pts_cylindrical[:, 1]
+        pts_cartesian[:, 0] = r * np.cos(phi)
+        pts_cartesian[:, 1] = r * np.sin(phi)
+
+        return super().write_data_to_vtk(
+            points=pts_cartesian,
+            filename=filename,
+            datasets=datasets,
+            volume_normalization=volume_normalization
+        )
 
 class SphericalMesh(StructuredMesh):
     """A 3D spherical mesh
@@ -1204,6 +1380,49 @@ class SphericalMesh(StructuredMesh):
 
         return np.multiply.outer(np.outer(V_r, V_t), V_p)
 
+    def write_data_to_vtk(self, filename, datasets, volume_normalization=True):
+        """Creates a VTK object of the mesh
+
+        Parameters
+        ----------
+        filename : str or pathlib.Path
+            Name of the VTK file to write.
+        datasets : dict
+            Dictionary whose keys are the data labels
+            and values are the data sets.
+        volume_normalization : bool, optional
+            Whether or not to normalize the data by
+            the volume of the mesh elements.
+            Defaults to True.
+
+        Returns
+        -------
+        vtk.vtkStructuredGrid
+            the VTK object
+        """
+
+        # create points
+        pts_spherical = np.array(
+            [
+                [r, theta, phi]
+                for phi in self.phi_grid
+                for theta in self.theta_grid
+                for r in self.r_grid
+            ]
+        )
+        pts_cartesian = np.copy(pts_spherical)
+        r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2]
+        pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta)
+        pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta)
+        pts_cartesian[:, 2] = r * np.cos(phi)
+
+        return super().write_data_to_vtk(
+            points=pts_cartesian,
+            filename=filename,
+            datasets=datasets,
+            volume_normalization=volume_normalization
+        )
+
 
 class UnstructuredMesh(MeshBase):
     """A 3D unstructured mesh
@@ -1328,6 +1547,7 @@ class UnstructuredMesh(MeshBase):
                                "been loaded from a statepoint file.")
         return len(self._centroids)
 
+
     @centroids.setter
     def centroids(self, centroids):
         cv.check_type("Unstructured mesh centroids", centroids,
@@ -1373,7 +1593,7 @@ class UnstructuredMesh(MeshBase):
 
         Parameters
         ----------
-        filename : str
+        filename : str or pathlib.Path
             Name of the VTK file to write.
         datasets : dict
             Dictionary whose keys are the data labels
@@ -1381,6 +1601,11 @@ class UnstructuredMesh(MeshBase):
         volume_normalization : bool
             Whether or not to normalize the data by the
             volume of the mesh elements
+
+        Raises
+        ------
+            RuntimeError
+                when the size of a dataset doesn't match the number of cells 
         """
 
         import vtk
@@ -1397,11 +1622,14 @@ class UnstructuredMesh(MeshBase):
                                " mesh information from a statepoint file.")
 
         # check that the data sets are appropriately sized
+        errmsg = "The size of the dataset {} should be equal to the number of cells"
         for label, dataset in datasets.items():
             if isinstance(dataset, np.ndarray):
-                assert dataset.size == self.n_elements
+                if not dataset.size == self.dimension[0] * self.dimension[1]* self.dimension[2]:
+                    raise RuntimeError(errmsg.format(label))
             else:
-                assert len(dataset) == self.n_elements
+                if len(dataset) == self.dimension[0] * self.dimension[1]* self.dimension[2]:
+                    raise RuntimeError(errmsg.format(label))
             cv.check_type('label', label, str)
 
         # create data arrays for the cells/points
diff --git a/openmc/plotter.py b/openmc/plotter.py
index 92acaf5bb5..16760868e3 100644
--- a/openmc/plotter.py
+++ b/openmc/plotter.py
@@ -542,15 +542,11 @@ def _calculate_cexs_elem_mat(this, types, temperature=294.,
 
     if isinstance(this, openmc.Material):
         # Expand elements in to nuclides with atomic densities
-        nuclides = this.get_nuclide_atom_densities()
-        # For ease of processing split out the nuclide and its fraction
-        nuc_fractions = {nuclide[1][0]: nuclide[1][1]
-                         for nuclide in nuclides.items()}
+        nuc_fractions = this.get_nuclide_atom_densities()
         # Create a dict of [nuclide name] = nuclide object to carry forward
         # with a common nuclides format between openmc.Material and
         # openmc.Element objects
-        nuclides = {nuclide[1][0]: nuclide[1][0]
-                    for nuclide in nuclides.items()}
+        nuclides = {nuclide: nuclide for nuclide in nuc_fractions}
     else:
         # Expand elements in to nuclides with atomic densities
         nuclides = this.expand(1., 'ao', enrichment=enrichment,
@@ -885,13 +881,13 @@ def _calculate_mgxs_elem_mat(this, types, library, orders=None,
         # Check to see if we have nuclides/elements or a macroscopic object
         if this._macroscopic is not None:
             # We have macroscopics
-            nuclides = {this._macroscopic: (this._macroscopic, this.density)}
+            nuclides = {this._macroscopic: this.density}
         else:
             # Expand elements in to nuclides with atomic densities
             nuclides = this.get_nuclide_atom_densities()
 
         # For ease of processing split out nuc and nuc_density
-        nuc_fraction = [nuclide[1][1] for nuclide in nuclides.items()]
+        nuc_fraction = list(nuclides.values())
     else:
         T = temperature
         # Expand elements in to nuclides with atomic densities
diff --git a/src/photon.cpp b/src/photon.cpp
index b0cd8aa87d..2f9bc5e3dd 100644
--- a/src/photon.cpp
+++ b/src/photon.cpp
@@ -14,7 +14,9 @@
 #include "openmc/settings.h"
 
 #include "xtensor/xbuilder.hpp"
+#include "xtensor/xmath.hpp"
 #include "xtensor/xoperation.hpp"
+#include "xtensor/xslice.hpp"
 #include "xtensor/xview.hpp"
 
 #include <cmath>
@@ -44,6 +46,8 @@ vector<unique_ptr<PhotonInteraction>> elements;
 
 PhotonInteraction::PhotonInteraction(hid_t group)
 {
+  using namespace xt::placeholders;
+
   // Set index of element in global vector
   index_ = data::elements.size();
 
@@ -125,6 +129,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   }
 
   shells_.resize(n_shell);
+  cross_sections_ = xt::zeros<double>({energy_.size(), n_shell});
 
   // Create mapping from designator to index
   std::unordered_map<int, int> shell_map;
@@ -155,13 +160,15 @@ PhotonInteraction::PhotonInteraction(hid_t group)
     read_attribute(tgroup, "num_electrons", shell.n_electrons);
 
     // Read subshell cross section
+    xt::xtensor<double, 1> xs;
     dset = open_dataset(tgroup, "xs");
     read_attribute(dset, "threshold_idx", shell.threshold);
     close_dataset(dset);
-    read_dataset(tgroup, "xs", shell.cross_section);
+    read_dataset(tgroup, "xs", xs);
 
-    auto& xs = shell.cross_section;
-    xs = xt::where(xs > 0.0, xt::log(xs), -500.0);
+    auto cross_section =
+      xt::view(cross_sections_, xt::range(shell.threshold, _), i);
+    cross_section = xt::where(xs > 0, xt::log(xs), 0);
 
     if (object_exists(tgroup, "transitions")) {
       // Determine dimensions of transitions
@@ -566,18 +573,13 @@ void PhotonInteraction::calculate_xs(Particle& p) const
 
   // Calculate microscopic photoelectric cross section
   xs.photoelectric = 0.0;
-  for (const auto& shell : shells_) {
-    // Check threshold of reaction
-    int i_start = shell.threshold;
-    if (i_grid < i_start)
-      continue;
+  const auto& xs_lower = xt::row(cross_sections_, i_grid);
+  const auto& xs_upper = xt::row(cross_sections_, i_grid + 1);
 
-    // Evaluation subshell photoionization cross section
-    xs.photoelectric +=
-      std::exp(shell.cross_section(i_grid - i_start) +
-               f * (shell.cross_section(i_grid + 1 - i_start) -
-                     shell.cross_section(i_grid - i_start)));
-  }
+  for (int i = 0; i < xs_upper.size(); ++i)
+    if (xs_lower(i) != 0)
+      xs.photoelectric +=
+        std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i)));
 
   // Calculate microscopic pair production cross section
   xs.pair_production = std::exp(
diff --git a/src/physics.cpp b/src/physics.cpp
index 13dd063de5..bcdcf7ecf9 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -29,6 +29,7 @@
 
 #include <algorithm> // for max, min, max_element
 #include <cmath>     // for sqrt, exp, log, abs, copysign
+#include <xtensor/xview.hpp>
 
 namespace openmc {
 
@@ -344,25 +345,26 @@ void sample_photon_reaction(Particle& p)
 
   // Photoelectric effect
   double prob_after = prob + micro.photoelectric;
+
   if (prob_after > cutoff) {
+    // Get grid index, interpolation factor, and bounding subshell
+    // cross sections
+    int i_grid = micro.index_grid;
+    double f = micro.interp_factor;
+    const auto& xs_lower = xt::row(element.cross_sections_, i_grid);
+    const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1);
+
     for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) {
       const auto& shell {element.shells_[i_shell]};
 
-      // Get grid index and interpolation factor
-      int i_grid = micro.index_grid;
-      double f = micro.interp_factor;
-
       // Check threshold of reaction
-      int i_start = shell.threshold;
-      if (i_grid < i_start)
+      if (xs_lower(i_shell) == 0)
         continue;
 
-      // Evaluation subshell photoionization cross section
-      double xs = std::exp(shell.cross_section(i_grid - i_start) +
-                           f * (shell.cross_section(i_grid + 1 - i_start) -
-                                 shell.cross_section(i_grid - i_start)));
+      //  Evaluation subshell photoionization cross section
+      prob += std::exp(
+        xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell)));
 
-      prob += xs;
       if (prob > cutoff) {
         double E_electron = p.E() - shell.binding_energy;
 
diff --git a/tests/unit_tests/test_data_misc.py b/tests/unit_tests/test_data_misc.py
index 6a81fb1e08..c3f7e3cfff 100644
--- a/tests/unit_tests/test_data_misc.py
+++ b/tests/unit_tests/test_data_misc.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 
-from collections.abc import Mapping
+from math import log
 import os
 from pathlib import Path
 
@@ -126,3 +126,14 @@ def test_zam():
     assert openmc.data.zam('Am242_m10') == (95, 242, 10)
     with pytest.raises(ValueError):
         openmc.data.zam('garbage')
+
+
+def test_half_life():
+    assert openmc.data.half_life('H2') is None
+    assert openmc.data.half_life('U235') == pytest.approx(2.22102e16)
+    assert openmc.data.half_life('Am242') == pytest.approx(57672.0)
+    assert openmc.data.half_life('Am242_m1') == pytest.approx(4449622000.0)
+    assert openmc.data.decay_constant('H2') == 0.0
+    assert openmc.data.decay_constant('U235') == pytest.approx(log(2.0)/2.22102e16)
+    assert openmc.data.decay_constant('Am242') == pytest.approx(log(2.0)/57672.0)
+    assert openmc.data.decay_constant('Am242_m1') == pytest.approx(log(2.0)/4449622000.0)
diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py
index 0b82a5fbcf..cb3b86b9c8 100644
--- a/tests/unit_tests/test_deplete_activation.py
+++ b/tests/unit_tests/test_deplete_activation.py
@@ -92,7 +92,7 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts
 
     w = model.geometry.get_materials_by_name('tungsten')[0]
     atom_densities = w.get_nuclide_atom_densities()
-    atom_per_cc = 1e24 * atom_densities['W186'][1]  # Density in atom/cm^3
+    atom_per_cc = 1e24 * atom_densities['W186']  # Density in atom/cm^3
     n0 = atom_per_cc * w.volume  # Absolute number of atoms
 
     # Pick a random irradiation time and then determine necessary source rate to
diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py
index c15a4b9e55..fd0e4f6247 100644
--- a/tests/unit_tests/test_material.py
+++ b/tests/unit_tests/test_material.py
@@ -100,7 +100,7 @@ def test_add_elements_by_formula():
     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():
+    for nuclide, adens in m.get_nuclide_atom_densities().items():
         if nuclide.startswith("Li"):
             elem["Li"] += adens
         if nuclide.startswith("Si"):
@@ -117,7 +117,7 @@ def test_add_elements_by_formula():
                 '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)
+        assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2)
 
     # testing the correct nuclides are added to the Material when enriched
     m = openmc.Material()
@@ -129,7 +129,7 @@ def test_add_elements_by_formula():
                 '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)
+        assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2)
 
     # testing the use of brackets
     m = openmc.Material()
@@ -137,7 +137,7 @@ def test_add_elements_by_formula():
 
     # 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():
+    for nuclide, adens in m.get_nuclide_atom_densities().items():
         if nuclide.startswith("Mg"):
             elem["Mg"] += adens
         if nuclide.startswith("N"):
@@ -155,7 +155,7 @@ def test_add_elements_by_formula():
                 '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)
+        assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2)
 
     # testing non integer multiplier results in a value error
     m = openmc.Material()
@@ -289,12 +289,21 @@ def test_get_nuclide_densities(uo2):
 
 
 def test_get_nuclide_atom_densities(uo2):
-    nucs = uo2.get_nuclide_atom_densities()
-    for nuc, density in nucs.values():
+    for nuc, density in uo2.get_nuclide_atom_densities().items():
         assert nuc in ('U235', 'O16')
         assert density > 0
 
 
+def test_get_nuclide_atoms():
+    mat = openmc.Material()
+    mat.add_nuclide('Li6', 1.0)
+    mat.set_density('atom/cm3', 3.26e20)
+    mat.volume = 100.0
+
+    atoms = mat.get_nuclide_atoms()
+    assert atoms['Li6'] == pytest.approx(mat.density * mat.volume)
+
+
 def test_mass():
     m = openmc.Material()
     m.add_nuclide('Zr90', 1.0, 'wo')
@@ -341,7 +350,7 @@ def test_borated_water():
                 'O16':2.4672e-02}
     nuc_dens = m.get_nuclide_atom_densities()
     for nuclide in ref_dens:
-        assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
+        assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2)
     assert m.id == 50
 
     # Test the Celsius conversion.
diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py
new file mode 100644
index 0000000000..26a393b411
--- /dev/null
+++ b/tests/unit_tests/test_mesh_to_vtk.py
@@ -0,0 +1,76 @@
+import numpy as np
+from pathlib import Path
+import pytest
+
+vtk = pytest.importorskip("vtk")
+from vtk.util import numpy_support as nps
+
+import openmc
+
+
+regular_mesh = openmc.RegularMesh()
+regular_mesh.lower_left = (0, 0, 0)
+regular_mesh.upper_right = (1, 1, 1)
+regular_mesh.dimension = [30, 20, 10]
+
+rectilinear_mesh = openmc.RectilinearMesh()
+rectilinear_mesh.x_grid = np.linspace(1, 2, num=30)
+rectilinear_mesh.y_grid = np.linspace(1, 2, num=30)
+rectilinear_mesh.z_grid = np.linspace(1, 2, num=30)
+
+cylinder_mesh = openmc.CylindricalMesh()
+cylinder_mesh.r_grid = np.linspace(1, 2, num=30)
+cylinder_mesh.phi_grid = np.linspace(0, np.pi, num=50)
+cylinder_mesh.z_grid = np.linspace(0, 1, num=30)
+
+spherical_mesh = openmc.SphericalMesh()
+spherical_mesh.r_grid = np.linspace(1, 2, num=30)
+spherical_mesh.phi_grid = np.linspace(0, np.pi, num=50)
+spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30)
+
+@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh])
+def test_write_data_to_vtk(mesh, tmpdir):
+    # BUILD
+    filename = Path(tmpdir) / "out.vtk"
+
+    data = np.random.random(mesh.num_mesh_cells)
+
+    # RUN
+    mesh.write_data_to_vtk(filename=filename, datasets={"label1": data, "label2": data})
+
+    # TEST
+    assert filename.is_file()
+
+    # read file
+    reader = vtk.vtkStructuredGridReader()
+    reader.SetFileName(filename)
+    reader.Update()
+
+    # check name of datasets
+    vtk_grid = reader.GetOutput()
+    array1 = vtk_grid.GetCellData().GetArray(0)
+    array2 = vtk_grid.GetCellData().GetArray(1)
+
+    assert array1.GetName() == "label1"
+    assert array2.GetName() == "label2"
+
+    # check size of datasets
+    assert nps.vtk_to_numpy(array1).size == data.size
+    assert nps.vtk_to_numpy(array2).size == data.size
+
+@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh])
+def test_write_data_to_vtk_size_mismatch(mesh):
+    """Checks that an error is raised when the size of the dataset
+    doesn't match the mesh number of cells
+
+    Parameters
+    ----------
+    mesh : openmc.StructuredMesh
+        The mesh to test
+    """
+    right_size = mesh.num_mesh_cells
+    data = np.random.random(right_size + 1)
+
+    expected_error_msg = "The size of the dataset label should be equal to the number of cells"
+    with pytest.raises(RuntimeError, match=expected_error_msg):
+        mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data})
\ No newline at end of file
diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py
index 10417b20ed..9b05ed2a5c 100644
--- a/tests/unit_tests/test_model.py
+++ b/tests/unit_tests/test_model.py
@@ -391,16 +391,14 @@ def test_py_lib_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     assert openmc.lib.materials[1].get_density('atom/b-cm') == \
         pytest.approx(0.06891296988603757, abs=1e-13)
     mat_a_dens = np.sum(
-        [v[1] for v in test_model.materials[0].
-            get_nuclide_atom_densities().values()])
+        list(test_model.materials[0].get_nuclide_atom_densities().values()))
     assert mat_a_dens == pytest.approx(0.06891296988603757, abs=1e-8)
     # Change the density
     test_model.update_densities(['UO2'], 2.)
     assert openmc.lib.materials[1].get_density('atom/b-cm') == \
         pytest.approx(2., abs=1e-13)
     mat_a_dens = np.sum(
-        [v[1] for v in test_model.materials[0].
-            get_nuclide_atom_densities().values()])
+        list(test_model.materials[0].get_nuclide_atom_densities().values()))
     assert mat_a_dens == pytest.approx(2., abs=1e-8)
 
     # Now lets do the cell temperature updates.
@@ -441,7 +439,7 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     test_model = openmc.Model(geom, mats, settings, tals, plots)
 
     initial_mat = mats[0].clone()
-    initial_u = initial_mat.get_nuclide_atom_densities()['U235'][1]
+    initial_u = initial_mat.get_nuclide_atom_densities()['U235']
 
     # Note that the chain file includes only U-235 fission to a stable Xe136 w/
     # a yield of 100%. Thus all the U235 we lose becomes Xe136
@@ -453,8 +451,8 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
                        operator_kwargs=op_kwargs,
                        power=1., output=False)
     # Get the new Xe136 and U235 atom densities
-    after_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1]
-    after_u = mats[0].get_nuclide_atom_densities()['U235'][1]
+    after_xe = mats[0].get_nuclide_atom_densities()['Xe136']
+    after_u = mats[0].get_nuclide_atom_densities()['U235']
     assert after_xe + after_u == pytest.approx(initial_u, abs=1e-15)
     assert test_model.is_initialized is False
 
@@ -462,7 +460,7 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
     mats[0].nuclides.clear()
     densities = initial_mat.get_nuclide_atom_densities()
     tot_density = 0.
-    for nuc, density in densities.values():
+    for nuc, density in densities.items():
         mats[0].add_nuclide(nuc, density)
         tot_density += density
     mats[0].set_density('atom/b-cm', tot_density)
@@ -473,8 +471,8 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
                        operator_kwargs=op_kwargs,
                        power=1., output=False)
     # Get the new Xe136 and U235 atom densities
-    after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1]
-    after_lib_u = mats[0].get_nuclide_atom_densities()['U235'][1]
+    after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136']
+    after_lib_u = mats[0].get_nuclide_atom_densities()['U235']
     assert after_lib_xe + after_lib_u == pytest.approx(initial_u, abs=1e-15)
     assert test_model.is_initialized is True
 
diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py
index aeeb0612d5..bfff741a97 100644
--- a/tests/unit_tests/test_tallies.py
+++ b/tests/unit_tests/test_tallies.py
@@ -1,4 +1,5 @@
 import numpy as np
+
 import openmc