diff --git a/CMakeLists.txt b/CMakeLists.txt
index c5150745b..51bdcc4d9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -116,9 +116,9 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
   endif()
 
   # GCC compiler options
-  list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2)
+  list(APPEND f90flags -cpp -std=f2008ts -fbacktrace -O2 -fstack-arrays)
   if(debug)
-    list(REMOVE_ITEM f90flags -O2)
+    list(REMOVE_ITEM f90flags -O2 -fstack-arrays)
     list(APPEND f90flags -g -Wall -Wno-unused-dummy-argument -pedantic
       -fbounds-check -ffpe-trap=invalid,overflow,underflow)
     list(APPEND ldflags -g)
diff --git a/src/cross_section.F90 b/src/cross_section.F90
index 69e577966..47e1c4fc5 100644
--- a/src/cross_section.F90
+++ b/src/cross_section.F90
@@ -42,7 +42,6 @@ contains
 
     ! Set all material macroscopic cross sections to zero
     material_xs % total          = ZERO
-    material_xs % elastic        = ZERO
     material_xs % absorption     = ZERO
     material_xs % fission        = ZERO
     material_xs % nu_fission     = ZERO
@@ -115,10 +114,6 @@ contains
         material_xs % total = material_xs % total + &
              atom_density * micro_xs(i_nuclide) % total
 
-        ! Add contributions to material macroscopic scattering cross section
-        material_xs % elastic = material_xs % elastic + &
-             atom_density * micro_xs(i_nuclide) % elastic
-
         ! Add contributions to material macroscopic absorption cross section
         material_xs % absorption = material_xs % absorption + &
              atom_density * micro_xs(i_nuclide) % absorption
@@ -162,6 +157,7 @@ contains
     real(8) :: sig_t, sig_a, sig_f ! Intermediate multipole variables
 
     ! Initialize cached cross sections to zero
+    micro_xs(i_nuclide) % elastic         = CACHE_INVALID
     micro_xs(i_nuclide) % thermal         = ZERO
     micro_xs(i_nuclide) % thermal_elastic = ZERO
 
@@ -182,13 +178,11 @@ contains
 
         micro_xs(i_nuclide) % total = sig_t
         micro_xs(i_nuclide) % absorption = sig_a
-        micro_xs(i_nuclide) % elastic = sig_t - sig_a
+        micro_xs(i_nuclide) % fission = sig_f
 
         if (nuc % fissionable) then
-          micro_xs(i_nuclide) % fission = sig_f
           micro_xs(i_nuclide) % nu_fission = sig_f * nuc % nu(E, EMISSION_TOTAL)
         else
-          micro_xs(i_nuclide) % fission    = ZERO
           micro_xs(i_nuclide) % nu_fission = ZERO
         end if
 
@@ -233,7 +227,7 @@ contains
           if (f > prn()) i_temp = i_temp + 1
         end select
 
-        associate (grid => nuc % grid(i_temp), xs => nuc % sum_xs(i_temp))
+        associate (grid => nuc % grid(i_temp), xs => nuc % xs(i_temp))
           ! Determine the energy grid index using a logarithmic mapping to
           ! reduce the energy range over which a binary search needs to be
           ! performed
@@ -266,25 +260,21 @@ contains
           micro_xs(i_nuclide) % interp_factor = f
 
           ! Calculate microscopic nuclide total cross section
-          micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) &
-               + f * xs % total(i_grid + 1)
-
-          ! Calculate microscopic nuclide elastic cross section
-          micro_xs(i_nuclide) % elastic = (ONE - f) * xs % elastic(i_grid) &
-               + f * xs % elastic(i_grid + 1)
+          micro_xs(i_nuclide) % total = (ONE - f) * xs % value(XS_TOTAL,i_grid) &
+               + f * xs % value(XS_TOTAL,i_grid + 1)
 
           ! Calculate microscopic nuclide absorption cross section
-          micro_xs(i_nuclide) % absorption = (ONE - f) * xs % absorption( &
-               i_grid) + f * xs % absorption(i_grid + 1)
+          micro_xs(i_nuclide) % absorption = (ONE - f) * xs % value(XS_ABSORPTION, &
+               i_grid) + f * xs % value(XS_ABSORPTION,i_grid + 1)
 
           if (nuc % fissionable) then
             ! Calculate microscopic nuclide total cross section
-            micro_xs(i_nuclide) % fission = (ONE - f) * xs % fission(i_grid) &
-                 + f * xs % fission(i_grid + 1)
+            micro_xs(i_nuclide) % fission = (ONE - f) * xs % value(XS_FISSION,i_grid) &
+                 + f * xs % value(XS_FISSION,i_grid + 1)
 
             ! Calculate microscopic nuclide nu-fission cross section
-            micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % nu_fission( &
-                 i_grid) + f * xs % nu_fission(i_grid + 1)
+            micro_xs(i_nuclide) % nu_fission = (ONE - f) * xs % value(XS_NU_FISSION, &
+                 i_grid) + f * xs % value(XS_NU_FISSION,i_grid + 1)
           else
             micro_xs(i_nuclide) % fission         = ZERO
             micro_xs(i_nuclide) % nu_fission      = ZERO
@@ -451,6 +441,9 @@ contains
     micro_xs(i_nuclide) % thermal = sab_frac * (elastic + inelastic)
     micro_xs(i_nuclide) % thermal_elastic = sab_frac * elastic
 
+    ! Calculate free atom elastic cross section
+    call calculate_elastic_xs(i_nuclide)
+
     ! Correct total and elastic cross sections
     micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total &
          + micro_xs(i_nuclide) % thermal &
@@ -579,6 +572,7 @@ contains
 
       ! Multiply by smooth cross-section if needed
       if (urr % multiply_smooth) then
+        call calculate_elastic_xs(i_nuclide)
         elastic = elastic * micro_xs(i_nuclide) % elastic
         capture = capture * (micro_xs(i_nuclide) % absorption - &
              micro_xs(i_nuclide) % fission)
@@ -607,6 +601,36 @@ contains
 
   end subroutine calculate_urr_xs
 
+!===============================================================================
+! CALCULATE_ELASTIC_XS precalculates the free atom elastic scattering cross
+! section. Normally it is not needed until a collision actually occurs in a
+! material. However, in the thermal and unresolved resonance regions, we have to
+! calculate it early to adjust the total cross section correctly.
+!===============================================================================
+
+  subroutine calculate_elastic_xs(i_nuclide)
+    integer, intent(in) :: i_nuclide
+
+    integer :: i_temp
+    integer :: i_grid
+    real(8) :: f
+
+    ! Get temperature index, grid index, and interpolation factor
+    i_temp =  micro_xs(i_nuclide) % index_temp
+    i_grid =  micro_xs(i_nuclide) % index_grid
+    f      =  micro_xs(i_nuclide) % interp_factor
+
+    if (i_temp > 0) then
+      associate (xs => nuclides(i_nuclide) % reactions(1) % xs(i_temp) % value)
+        micro_xs(i_nuclide) % elastic = (ONE - f)*xs(i_grid) + f*xs(i_grid + 1)
+      end associate
+    else
+      ! For multipole, elastic is total - absorption
+      micro_xs(i_nuclide) % elastic = micro_xs(i_nuclide) % total - &
+           micro_xs(i_nuclide) % absorption
+    end if
+  end subroutine calculate_elastic_xs
+
 !===============================================================================
 ! MULTIPOLE_EVAL evaluates the windowed multipole equations for cross
 ! sections in the resolved resonance regions
diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90
index 38e89b648..1ca584914 100644
--- a/src/nuclide_header.F90
+++ b/src/nuclide_header.F90
@@ -36,13 +36,18 @@ module nuclide_header
     real(8), allocatable :: energy(:)     ! energy values corresponding to xs
   end type EnergyGrid
 
+  ! Positions for first dimension of Nuclide % xs
+  integer, parameter :: &
+       XS_TOTAL      = 1, &
+       XS_ABSORPTION = 2, &
+       XS_FISSION    = 3, &
+       XS_NU_FISSION = 4
+
+  ! The array within SumXS is of shape (4, n_energy) where the first dimension
+  ! corresponds to the following values: 1) total, 2) absorption (MT > 100), 3)
+  ! fission, 4) neutron production
   type SumXS
-    real(8), allocatable :: total(:)      ! total cross section
-    real(8), allocatable :: elastic(:)    ! elastic scattering
-    real(8), allocatable :: fission(:)    ! fission
-    real(8), allocatable :: nu_fission(:) ! neutron production
-    real(8), allocatable :: absorption(:) ! absorption (MT > 100)
-    real(8), allocatable :: heating(:)    ! heating
+    real(8), allocatable :: value(:,:)
   end type SumXS
 
   type :: Nuclide
@@ -61,7 +66,7 @@ module nuclide_header
     type(EnergyGrid), allocatable :: grid(:)
 
     ! Microscopic cross sections
-    type(SumXS), allocatable :: sum_xs(:)
+    type(SumXS), allocatable :: xs(:)
 
     ! Resonance scattering info
     logical              :: resonant = .false. ! resonant scatterer?
@@ -110,14 +115,19 @@ module nuclide_header
 ! nuclide at the current energy
 !===============================================================================
 
+  ! Arbitrary value to indicate invalid cache state for elastic scattering
+  ! (NuclideMicroXS % elastic)
+  real(8), parameter :: CACHE_INVALID = dble(Z"FFE0000000000000")
+
   type NuclideMicroXS
     ! Microscopic cross sections in barns
     real(8) :: total
-    real(8) :: elastic          ! If sab_frac is not 1 or 0, then this value is
-                                !   averaged over bound and non-bound nuclei
     real(8) :: absorption       ! absorption (disappearance)
     real(8) :: fission          ! fission
     real(8) :: nu_fission       ! neutron production from fission
+
+    real(8) :: elastic          ! If sab_frac is not 1 or 0, then this value is
+                                !   averaged over bound and non-bound nuclei
     real(8) :: thermal          ! Bound thermal elastic & inelastic scattering
     real(8) :: thermal_elastic  ! Bound thermal elastic scattering
 
@@ -148,7 +158,6 @@ module nuclide_header
 
   type MaterialMacroXS
     real(8) :: total         ! macroscopic total xs
-    real(8) :: elastic       ! macroscopic elastic scattering xs
     real(8) :: absorption    ! macroscopic absorption xs
     real(8) :: fission       ! macroscopic fission xs
     real(8) :: nu_fission    ! macroscopic production xs
@@ -566,21 +575,13 @@ contains
     type(VectorInt) :: MTs
 
     n_temperature = size(this % kTs)
-    allocate(this % sum_xs(n_temperature))
+    allocate(this % xs(n_temperature))
     this % reaction_index(:) = 0
     do i = 1, n_temperature
       ! Allocate and initialize derived cross sections
       n_grid = size(this % grid(i) % energy)
-      allocate(this % sum_xs(i) % total(n_grid))
-      allocate(this % sum_xs(i) % elastic(n_grid))
-      allocate(this % sum_xs(i) % fission(n_grid))
-      allocate(this % sum_xs(i) % nu_fission(n_grid))
-      allocate(this % sum_xs(i) % absorption(n_grid))
-      this % sum_xs(i) % total(:) = ZERO
-      this % sum_xs(i) % elastic(:) = ZERO
-      this % sum_xs(i) % fission(:) = ZERO
-      this % sum_xs(i) % nu_fission(:) = ZERO
-      this % sum_xs(i) % absorption(:) = ZERO
+      allocate(this % xs(i) % value(4,n_grid))
+      this % xs(i) % value(:,:) = ZERO
     end do
 
     i_fission = 0
@@ -607,17 +608,14 @@ contains
           j = rx % xs(t) % threshold
           n = size(rx % xs(t) % value)
 
-          ! Copy elastic
-          if (rx % MT == ELASTIC) this % sum_xs(t) % elastic(:) = rx % xs(t) % value
-
           ! Add contribution to total cross section
-          this % sum_xs(t) % total(j:j+n-1) = this % sum_xs(t) % total(j:j+n-1) + &
-               rx % xs(t) % value
+          this % xs(t) % value(XS_TOTAL,j:j+n-1) = this % xs(t) % &
+               value(XS_TOTAL,j:j+n-1) + rx % xs(t) % value
 
           ! Add contribution to absorption cross section
           if (is_disappearance(rx % MT)) then
-            this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
-                 absorption(j:j+n-1) + rx % xs(t) % value
+            this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % &
+                 value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
           end if
 
           ! Information about fission reactions
@@ -633,12 +631,12 @@ contains
           ! Add contribution to fission cross section
           if (is_fission(rx % MT)) then
             this % fissionable = .true.
-            this % sum_xs(t) % fission(j:j+n-1) = this % sum_xs(t) % &
-                 fission(j:j+n-1) + rx % xs(t) % value
+            this % xs(t) % value(XS_FISSION,j:j+n-1) = this % xs(t) % &
+                 value(XS_FISSION,j:j+n-1) + rx % xs(t) % value
 
             ! Also need to add fission cross sections to absorption
-            this % sum_xs(t) % absorption(j:j+n-1) = this % sum_xs(t) % &
-                 absorption(j:j+n-1) + rx % xs(t) % value
+            this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % &
+                 value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
 
             ! If total fission reaction is present, there's no need to store the
             ! reaction cross-section since it was copied to this % fission
@@ -689,12 +687,11 @@ contains
     ! Calculate nu-fission cross section
     do t = 1, n_temperature
       if (this % fissionable) then
-        do i = 1, size(this % sum_xs(t) % fission)
-          this % sum_xs(t) % nu_fission(i) = this % nu(this % grid(t) % energy(i), &
-               EMISSION_TOTAL) * this % sum_xs(t) % fission(i)
+        do i = 1, n_grid
+          this % xs(t) % value(XS_NU_FISSION,i) = &
+               this % nu(this % grid(t) % energy(i), EMISSION_TOTAL) * &
+               this % xs(t) % value(XS_FISSION,i)
         end do
-      else
-        this % sum_xs(t) % nu_fission(:) = ZERO
       end if
     end do
   end subroutine nuclide_create_derived
diff --git a/src/physics.F90 b/src/physics.F90
index d13313e05..fe242dbb7 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -2,7 +2,7 @@ module physics
 
   use algorithm,              only: binary_search
   use constants
-  use cross_section,          only: elastic_xs_0K
+  use cross_section,          only: elastic_xs_0K, calculate_elastic_xs
   use endf,                   only: reaction_name
   use error,                  only: fatal_error, warning, write_message
   use material_header,        only: Material, materials
@@ -338,6 +338,11 @@ contains
          micro_xs(i_nuclide) % absorption)
     sampled = .false.
 
+    ! Calculate elastic cross section if it wasn't precalculated
+    if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then
+      call calculate_elastic_xs(i_nuclide)
+    end if
+
     prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal
     if (prob > cutoff) then
       ! =======================================================================
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index 7684a36aa..027e335e9 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -4,7 +4,7 @@ module tally
 
   use algorithm,        only: binary_search
   use constants
-  use cross_section,    only: multipole_deriv_eval
+  use cross_section,    only: multipole_deriv_eval, calculate_elastic_xs
   use dict_header,      only: EMPTY
   use error,            only: fatal_error
   use geometry_header
@@ -1017,9 +1017,26 @@ contains
 
         else
           if (i_nuclide > 0) then
+            if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then
+              call calculate_elastic_xs(i_nuclide)
+            end if
             score = micro_xs(i_nuclide) % elastic * atom_density * flux
           else
-            score = material_xs % elastic * flux
+            score = ZERO
+            if (p % material /= MATERIAL_VOID) then
+              do l = 1, materials(p % material) % n_nuclides
+                ! Get atom density
+                atom_density_ = materials(p % material) % atom_density(l)
+
+                ! Get index in nuclides array
+                i_nuc = materials(p % material) % nuclide(l)
+                if (micro_xs(i_nuc) % elastic == CACHE_INVALID) then
+                  call calculate_elastic_xs(i_nuc)
+                end if
+
+                score = score + micro_xs(i_nuc) % elastic * atom_density_ * flux
+              end do
+            end if
           end if
         end if