diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 363c08867..91cd57ae2 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -2450,17 +2450,6 @@ contains
     n_nuclides = index_nuclide
     n_sab_tables = index_sab
 
-    ! Create direct address tables for tally optimization purposes
-    do i = 1, n_materials
-      mat => materials(i)
-      allocate(mat % mat_nuclide_index(n_nuclides))
-      mat % mat_nuclide_index(:) = 0
-
-      do j = 1, mat % n_nuclides
-        mat % mat_nuclide_index(mat % nuclide(j)) = j
-      end do
-    end do
-
     ! Close materials XML file
     call doc % clear()
 
diff --git a/src/material_header.F90 b/src/material_header.F90
index 3287e796f..d77220415 100644
--- a/src/material_header.F90
+++ b/src/material_header.F90
@@ -69,6 +69,7 @@ module material_header
 
   contains
     procedure :: set_density => material_set_density
+    procedure :: init_nuclide_index => material_init_nuclide_index
     procedure :: assign_sab_tables => material_assign_sab_tables
   end type Material
 
@@ -85,8 +86,8 @@ contains
 ! MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm.
 !===============================================================================
 
-  function material_set_density(m, density) result(err)
-    class(Material), intent(inout) :: m
+  function material_set_density(this, density) result(err)
+    class(Material), intent(inout) :: this
     real(8), intent(in) :: density
     integer :: err
 
@@ -94,23 +95,23 @@ contains
     real(8) :: sum_percent
     real(8) :: awr
 
-    if (allocated(m % atom_density)) then
+    if (allocated(this % atom_density)) then
       ! Set total density based on value provided
-      m % density = density
+      this % density = density
 
       ! Determine normalized atom percents
-      sum_percent = sum(m % atom_density)
-      m % atom_density(:) = m % atom_density / sum_percent
+      sum_percent = sum(this % atom_density)
+      this % atom_density(:) = this % atom_density / sum_percent
 
       ! Recalculate nuclide atom densities based on given density
-      m % atom_density(:) = density * m % atom_density
+      this % atom_density(:) = density * this % atom_density
 
       ! Calculate density in g/cm^3.
-      m % density_gpcc = ZERO
-      do i = 1, m % n_nuclides
-        awr = nuclides(m % nuclide(i)) % awr
-        m % density_gpcc = m % density_gpcc &
-             + m % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO
+      this % density_gpcc = ZERO
+      do i = 1, this % n_nuclides
+        awr = nuclides(this % nuclide(i)) % awr
+        this % density_gpcc = this % density_gpcc &
+             + this % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO
       end do
       err = 0
     else
@@ -119,6 +120,28 @@ contains
     end if
   end function material_set_density
 
+!===============================================================================
+! INIT_NUCLIDE_INDEX creates a mapping from indices in the global nuclides(:)
+! array to the Material % nuclides array
+!===============================================================================
+
+  subroutine material_init_nuclide_index(this)
+    class(Material), intent(inout) :: this
+
+    integer :: i
+
+    ! Allocate nuclide index array and set to zeros
+    if (allocated(this % mat_nuclide_index)) &
+         deallocate(this % mat_nuclide_index)
+    allocate(this % mat_nuclide_index(n_nuclides))
+    this % mat_nuclide_index(:) = 0
+
+    ! Assign entries in the index array
+    do i = 1, this % n_nuclides
+      this % mat_nuclide_index(this % nuclide(i)) = i
+    end do
+  end subroutine material_init_nuclide_index
+
 !===============================================================================
 ! ASSIGN_SAB_TABLES assigns S(alpha,beta) tables to specific nuclides within
 ! materials so the code knows when to apply bound thermal scattering data
diff --git a/src/simulation.F90 b/src/simulation.F90
index 5eacbd40b..ade209e36 100644
--- a/src/simulation.F90
+++ b/src/simulation.F90
@@ -18,6 +18,7 @@ module simulation
 #endif
   use error,           only: fatal_error
   use geometry_header, only: n_cells
+  use material_header, only: n_materials, materials
   use message_passing
   use mgxs_header,     only: energy_bins, energy_bin_avg
   use nuclide_header,  only: micro_xs, n_nuclides
@@ -401,6 +402,7 @@ contains
 !===============================================================================
 
   subroutine openmc_simulation_init() bind(C)
+    integer :: i
 
     ! Skip if simulation has already been initialized
     if (simulation_initialized) return
@@ -418,6 +420,11 @@ contains
     ! Allocate tally results arrays if they're not allocated yet
     call configure_tallies()
 
+    ! Set up material nuclide index mapping
+    do i = 1, n_materials
+      call materials(i) % init_nuclide_index()
+    end do
+
 !$omp parallel
     ! Allocate array for microscopic cross section cache
     allocate(micro_xs(n_nuclides))
@@ -459,8 +466,8 @@ contains
 
   subroutine openmc_simulation_finalize() bind(C)
 
+    integer    :: i       ! loop index
 #ifdef MPI
-    integer    :: i       ! loop index for tallies
     integer    :: n       ! size of arrays
     integer    :: mpi_err  ! MPI error code
     integer(8) :: temp
@@ -470,9 +477,14 @@ contains
     ! Skip if simulation was never run
     if (.not. simulation_initialized) return
 
-    ! Stop active batch timer
+    ! Stop active batch timer and start finalization timer
     call time_active % stop()
+    call time_finalize % start()
 
+    ! Free up simulation-specific memory
+    do i = 1, n_materials
+      deallocate(materials(i) % mat_nuclide_index)
+    end do
 !$omp parallel
     deallocate(micro_xs, filter_matches)
 !$omp end parallel
@@ -480,9 +492,6 @@ contains
     ! Increment total number of generations
     total_gen = total_gen + current_batch*gen_per_batch
 
-    ! Start finalization timer
-    call time_finalize % start()
-
 #ifdef MPI
     ! Broadcast tally results so that each process has access to results
     if (allocated(tallies)) then