diff --git a/src/global.F90 b/src/global.F90
index 700f43ddc0..a4398c0f5e 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -119,8 +119,13 @@ module global
 
   ! Energy group structure
   real(8), allocatable :: energy_bins(:)
+
+  ! Midpoint of the energy group structure
   real(8), allocatable :: energy_bin_avg(:)
 
+  ! Inverse velocities of the energy groups (provided or estimated)
+  real(8), allocatable :: inverse_velocities(:)
+
   ! Maximum Data Order
   integer :: max_order
 
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index c1aaaf53c5..47370e39dc 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -4390,6 +4390,19 @@ contains
       energy_bin_avg(i) = 0.5_8 * (energy_bins(i) + energy_bins(i + 1))
     end do
 
+    allocate(inverse_velocities(energy_groups))
+    if (check_for_node(doc, "inverse_velocities")) then
+      ! Get inverse velocities
+      call get_node_array(doc, "inverse_velocities", inverse_velocities)
+    else
+      ! If not given, estimate them by using average energy in group which is
+      ! assumed to be the midpoint
+      do i = 1, energy_groups
+        inverse_velocities(i) = &
+             (sqrt(TWO * energy_bin_avg(i) / (MASS_NEUTRON_MEV)) * C_LIGHT)
+      end do
+    end if
+
     ! Get node list of all <xsdata>
     call get_node_list(doc, "xsdata", node_xsdata_list)
     n_listings = get_list_size(node_xsdata_list)
diff --git a/src/tally.F90 b/src/tally.F90
index 848a6cb7eb..35524c1570 100644
--- a/src/tally.F90
+++ b/src/tally.F90
@@ -872,12 +872,11 @@ contains
           else
             score = p % last_wgt
           end if
-          score = score / material_xs % total &
-               / (sqrt(TWO * p % E / (MASS_NEUTRON_MEV)) * C_LIGHT)
+          score = score * inverse_velocities(p % last_g)
 
         else
           ! For inverse velocity, we need no cross section
-          score = flux / (sqrt(TWO * p % E / (MASS_NEUTRON_MEV)) * C_LIGHT)
+          score = score * inverse_velocities(p % g)
         end if