diff --git a/include/openmc/timer.h b/include/openmc/timer.h
index 8fbcee46b..ac8c81e0c 100644
--- a/include/openmc/timer.h
+++ b/include/openmc/timer.h
@@ -37,7 +37,7 @@ public:
   Timer() {};
 
   //! Start running the timer
-  void start ();
+  void start();
 
   //! Get total elapsed time in seconds
   //! \return Elapsed time in [s]
@@ -52,7 +52,7 @@ public:
 private:
   bool running_ {false}; //!< is timer running?
   std::chrono::time_point<clock> start_; //!< starting point for clock
-  double elapsed_ {0.0}; //!< elasped time in [s]
+  double elapsed_ {0.0}; //!< elapsed time in [s]
 };
 
 //==============================================================================
diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py
index 0b17b49b5..ac519cf87 100644
--- a/openmc/deplete/integrator/cecm.py
+++ b/openmc/deplete/integrator/cecm.py
@@ -7,7 +7,7 @@ from .cram import deplete
 from ..results import Results
 
 
-def cecm(operator, timesteps, power, print_out=True):
+def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
     r"""Deplete using the CE/CM algorithm.
 
     Implements the second order `CE/CM predictor-corrector algorithm
@@ -31,16 +31,29 @@ def cecm(operator, timesteps, power, print_out=True):
         The operator object to simulate on.
     timesteps : iterable of float
         Array of timesteps in units of [s]. Note that values are not cumulative.
-    power : float or iterable of float
+    power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that the power is
         constant over all timesteps. An iterable indicates potentially different
         power levels for each timestep. For a 2D problem, the power can be given
         in [W/cm] as long as the "volume" assigned to a depletion material is
-        actually an area in [cm^2].
+        actually an area in [cm^2]. Either `power` or `power_density` must be
+        specified.
+    power_density : float or iterable of float, optional
+        Power density of the reactor in [W/gHM]. It is multiplied by initial
+        heavy metal inventory to get total power if `power` is not speficied.
     print_out : bool, optional
         Whether or not to print out time.
 
     """
+    if power is None:
+        if power_density is None:
+            raise ValueError(
+                "Neither power nor power density was specified.")
+        if not isinstance(power_density, Iterable):
+            power = power_density*operator.heavy_metal
+        else:
+            power = [i*operator.heavy_metal for i in power_density]
+
     if not isinstance(power, Iterable):
         power = [power]*len(timesteps)
 
diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py
index a45017813..969144f68 100644
--- a/openmc/deplete/integrator/predictor.py
+++ b/openmc/deplete/integrator/predictor.py
@@ -7,7 +7,8 @@ from .cram import deplete
 from ..results import Results
 
 
-def predictor(operator, timesteps, power, print_out=True):
+def predictor(operator, timesteps, power=None, power_density=None,
+              print_out=True):
     r"""Deplete using a first-order predictor algorithm.
 
     Implements the first-order predictor algorithm. This algorithm is
@@ -26,16 +27,29 @@ def predictor(operator, timesteps, power, print_out=True):
         The operator object to simulate on.
     timesteps : iterable of float
         Array of timesteps in units of [s]. Note that values are not cumulative.
-    power : float or iterable of float
+    power : float or iterable of float, optional
         Power of the reactor in [W]. A single value indicates that the power is
         constant over all timesteps. An iterable indicates potentially different
         power levels for each timestep. For a 2D problem, the power can be given
         in [W/cm] as long as the "volume" assigned to a depletion material is
-        actually an area in [cm^2].
+        actually an area in [cm^2]. Either `power` or `power_density` must be
+        specified.
+    power_density : float or iterable of float, optional
+        Power density of the reactor in [W/gHM]. It is multiplied by initial
+        heavy metal inventory to get total power if `power` is not speficied.
     print_out : bool, optional
         Whether or not to print out time.
 
     """
+    if power is None:
+        if power_density is None:
+            raise ValueError(
+                "Neither power nor power density was specified.")
+        if not isinstance(power_density, Iterable):
+            power = power_density*operator.heavy_metal
+        else:
+            power = [i*operator.heavy_metal for i in power_density]
+
     if not isinstance(power, Iterable):
         power = [power]*len(timesteps)
 
diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py
index a0b2417c8..754e22f7d 100644
--- a/openmc/deplete/operator.py
+++ b/openmc/deplete/operator.py
@@ -70,6 +70,8 @@ class Operator(TransportOperator):
         Results from a previous depletion calculation. If this argument is
         specified, the depletion calculation will start from the latest state
         in the previous results.
+    diff_burnable_mats : bool, optional
+        Whether to differentiate burnable materials with multiple instances
 
     Attributes
     ----------
@@ -96,17 +98,23 @@ class Operator(TransportOperator):
         Reaction rates from the last operator step.
     burnable_mats : list of str
         All burnable material IDs
+    heavy_metal : float
+        Initial heavy metal inventory
     local_mats : list of str
         All burnable material IDs being managed by a single process
     prev_res : ResultsList
         Results from a previous depletion calculation
+    diff_burnable_mats : bool
+        Whether to differentiate burnable materials with multiple instances
 
     """
-    def __init__(self, geometry, settings, chain_file=None, prev_results=None):
+    def __init__(self, geometry, settings, chain_file=None, prev_results=None,
+                 diff_burnable_mats=False):
         super().__init__(chain_file)
         self.round_number = False
         self.settings = settings
         self.geometry = geometry
+        self.diff_burnable_mats = diff_burnable_mats
 
         if prev_results != None:
             # Reload volumes into geometry
@@ -185,8 +193,28 @@ class Operator(TransportOperator):
 
         return copy.deepcopy(op_result)
 
+    def _differentiate_burnable_mats(self):
+        """Assign distribmats for each burnable material
+
+        """
+
+        # Count the number of instances for each cell and material
+        self.geometry.determine_paths(instances_only=True)
+
+        # Extract all burnable materials which have multiple instances
+        distribmats = set(
+            [mat for mat in self.geometry.get_all_materials().values()
+             if mat.depletable and mat.num_instances > 1])
+
+        if distribmats:
+            # Assign distribmats to cells
+            for cell in self.geometry.get_all_material_cells().values():
+                if cell.fill in distribmats and cell.num_instances > 1:
+                    cell.fill = [cell.fill.clone()
+                                 for i in range(cell.num_instances)]
+
     def _get_burnable_mats(self):
-        """Determine depletable materials, volumes, and nuclids
+        """Determine depletable materials, volumes, and nuclides
 
         Returns
         -------
@@ -198,10 +226,17 @@ class Operator(TransportOperator):
             Nuclides in order of how they'll appear in the simulation.
 
         """
+
+        if self.diff_burnable_mats:
+            # Automatically distribute burnable materials
+            self._differentiate_burnable_mats()
+
         burnable_mats = set()
         model_nuclides = set()
         volume = OrderedDict()
 
+        self.heavy_metal = 0.0
+
         # Iterate once through the geometry to get dictionaries
         for mat in self.geometry.get_all_materials().values():
             for nuclide in mat.get_nuclides():
@@ -212,6 +247,7 @@ class Operator(TransportOperator):
                     raise RuntimeError("Volume not specified for depletable "
                                        "material with ID={}.".format(mat.id))
                 volume[str(mat.id)] = mat.volume
+                self.heavy_metal += mat.fissionable_mass
 
         # Make sure there are burnable materials
         if not burnable_mats:
diff --git a/openmc/material.py b/openmc/material.py
index 612a7dfb6..d7c93ebae 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -719,9 +719,9 @@ class Material(IDManagerMixin):
         """
         mass_density = 0.0
         for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
-            density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
-                        / openmc.data.AVOGADRO
             if nuclide is None or nuclide == nuc:
+                density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
+                            / openmc.data.AVOGADRO
                 mass_density += density_i
         return mass_density
 
diff --git a/src/finalize.cpp b/src/finalize.cpp
index 5885470c1..84b68e82e 100644
--- a/src/finalize.cpp
+++ b/src/finalize.cpp
@@ -23,6 +23,10 @@ int openmc_finalize()
   // Clear results
   openmc_reset();
 
+  // Reset timers
+  reset_timers();
+  reset_timers_f();
+
   // Reset global variables
   settings::assume_separate = false;
   settings::check_overlaps = false;
@@ -115,10 +119,6 @@ int openmc_reset()
   simulation::k_abs_tra = 0.0;
   simulation::k_sum = {0.0, 0.0};
 
-  // Reset timers
-  reset_timers();
-  reset_timers_f();
-
   return 0;
 }
 
@@ -126,6 +126,8 @@ int openmc_hard_reset()
 {
   // Reset all tallies and timers
   openmc_reset();
+  reset_timers();
+  reset_timers_f();
 
   // Reset total generations and keff guess
   simulation::keff = 1.0;
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index d723aae99..41bac0f88 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -2208,6 +2208,8 @@ contains
           call write_message("Maximum neutron transport energy: " // &
                trim(to_str(energy_max(NEUTRON))) // " eV for " // &
                trim(adjustl(nuclides(i) % name)), 7)
+          if (master .and. energy_max(NEUTRON) < 20.e6) call warning("Maximum &
+               &neutron energy is below 20 MeV. This may bias the results.")
           exit
         end if
       end if
@@ -2222,9 +2224,10 @@ contains
           exit
         end if
       end do
-      if (.not. mp_found) call warning("Windowed multipole functionality is &
-           &turned on, but no multipole libraries were found. Make sure that &
-           &windowed multipole data is present in your cross_sections.xml file.")
+      if (master .and. .not. mp_found) call warning("Windowed multipole &
+           &functionality is turned on, but no multipole libraries were found. &
+           &Make sure that windowed multipole data is present in your &
+           &cross_sections.xml file.")
     end if
 
     call already_read % clear()
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 296a4a6a2..fe35f7293 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -685,8 +685,8 @@ contains
 
   subroutine nuclide_init_grid(this, E_min, E_max, M)
     class(Nuclide), intent(inout) :: this
-    real(8), intent(in) :: E_min           ! Minimum energy in MeV
-    real(8), intent(in) :: E_max           ! Maximum energy in MeV
+    real(8), intent(in) :: E_min           ! Minimum energy in eV
+    real(8), intent(in) :: E_max           ! Maximum energy in eV
     integer, intent(in) :: M               ! Number of equally log-spaced bins
 
     integer :: i, j, k               ! Loop indices