diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index 4f92a64ee..3a079a932 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -15,7 +15,7 @@ void Mgxs::init(const std::string& in_name, const double in_awr,
      const double_1dvec& in_kTs, const bool in_fissionable,
      const int in_scatter_format, const int in_num_groups,
      const int in_num_delayed_groups, const double_1dvec& in_polar,
-     const double_1dvec& in_azimuthal)
+     const double_1dvec& in_azimuthal, const int n_threads)
 {
   name = in_name;
   awr = in_awr;
@@ -29,20 +29,23 @@ void Mgxs::init(const std::string& in_name, const double in_awr,
   azimuthal = in_azimuthal;
   n_pol = polar.size();
   n_azi = azimuthal.size();
-  index_temp = 0;
-  last_sqrtkT = 0.;
-  index_pol = 0;
-  index_azi = 0;
-  last_uvw[0] = 1.;
-  last_uvw[1] = 0.;
-  last_uvw[2] = 0.;
+  cache.resize(n_threads);
+  for (int thread = 0; thread < n_threads; thread++) {
+    cache[thread].sqrtkT = 0.;
+    cache[thread].t = 0;
+    cache[thread].p = 0;
+    cache[thread].a = 0;
+    cache[thread].uvw[0] = 1.;
+    cache[thread].uvw[1] = 0.;
+    cache[thread].uvw[2] = 0.;
+  }
 }
 
 
 void Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
      const int in_num_delayed_groups, double_1dvec& temperature, int& method,
      const double tolerance, int_1dvec& temps_to_read, int& order_dim,
-     bool& is_isotropic)
+     bool& is_isotropic, const int n_threads)
 {
   // get name
   char char_name[MAX_WORD_LEN];
@@ -242,21 +245,23 @@ void Mgxs::_metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
 
   // Finally use this data to initialize the MGXS Object
   init(in_name, in_awr, in_kTs, in_fissionable, in_scatter_format,
-       in_num_groups, in_num_delayed_groups, in_polar, in_azimuthal);
+       in_num_groups, in_num_delayed_groups, in_polar, in_azimuthal,
+       n_threads);
 }
 
 
-void Mgxs::from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups,
-                     double_1dvec& temperature, int& method, double tolerance,
-                     int max_order, bool legendre_to_tabular,
-                     int legendre_to_tabular_points)
+void Mgxs::from_hdf5(hid_t xs_id, const int energy_groups,
+     const int delayed_groups, double_1dvec& temperature, int& method,
+     const double tolerance, const int max_order,
+     const  bool legendre_to_tabular, const int legendre_to_tabular_points,
+     const int n_threads)
 {
   // Call generic data gathering routine (will populate the metadata)
   int order_data;
   int_1dvec temps_to_read;
   bool is_isotropic;
   _metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
-       method, tolerance, temps_to_read, order_data, is_isotropic);
+       method, tolerance, temps_to_read, order_data, is_isotropic, n_threads);
 
   // Set number of energy and delayed groups
   int final_scatter_format = scatter_format;
@@ -285,8 +290,8 @@ void Mgxs::from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups,
 
 
 void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs,
-                       std::vector<Mgxs*>& micros, double_1dvec& atom_densities,
-                       int& method, double tolerance)
+     std::vector<Mgxs*>& micros, double_1dvec& atom_densities, int& method,
+     const double tolerance, const int n_threads)
 {
   // Get the minimum data needed to initialize:
   // Dont need awr, but lets just initialize it anyways
@@ -305,7 +310,8 @@ void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs,
   double_1dvec in_azimuthal = micros[0]->azimuthal;
 
   init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format,
-       in_num_groups, in_num_delayed_groups, in_polar, in_azimuthal);
+       in_num_groups, in_num_delayed_groups, in_polar, in_azimuthal,
+       n_threads);
 
   // Create the xs data for each temperature
   for (int t = 0; t < mat_kTs.size(); t++) {
@@ -395,52 +401,52 @@ void Mgxs::combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
 }
 
 
-double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
-                    int* dg)
+double Mgxs::get_xs(const int tid, const int xstype, const int gin,
+                    int* gout, double* mu, int* dg)
 {
   // This method assumes that the temperature and angle indices are set
   double val;
   switch(xstype) {
   case MG_GET_XS_TOTAL:
-    val = xs[index_temp].total[index_pol][index_azi][gin];
+    val = xs[cache[tid].t].total[cache[tid].p][cache[tid].a][gin];
     break;
   case MG_GET_XS_ABSORPTION:
-    val = xs[index_temp].absorption[index_pol][index_azi][gin];
+    val = xs[cache[tid].t].absorption[cache[tid].p][cache[tid].a][gin];
     break;
   case MG_GET_XS_INVERSE_VELOCITY:
-    val = xs[index_temp].inverse_velocity[index_pol][index_azi][gin];
+    val = xs[cache[tid].t].inverse_velocity[cache[tid].p][cache[tid].a][gin];
     break;
   case MG_GET_XS_DECAY_RATE:
     if (dg != nullptr) {
-      val = xs[index_temp].decay_rate[index_pol][index_azi][*dg + 1];
+      val = xs[cache[tid].t].decay_rate[cache[tid].p][cache[tid].a][*dg + 1];
     } else {
-      val = xs[index_temp].decay_rate[index_pol][index_azi][0];
+      val = xs[cache[tid].t].decay_rate[cache[tid].p][cache[tid].a][0];
     }
     break;
   case MG_GET_XS_SCATTER:
   case MG_GET_XS_SCATTER_MULT:
   case MG_GET_XS_SCATTER_FMU_MULT:
   case MG_GET_XS_SCATTER_FMU:
-    val = xs[index_temp].scatter[index_pol]
-                    [index_azi]->get_xs(xstype, gin, gout, mu);
+    val = xs[cache[tid].t].scatter[cache[tid].p]
+                    [cache[tid].a]->get_xs(xstype, gin, gout, mu);
     break;
   case MG_GET_XS_FISSION:
     if (fissionable) {
-      val = xs[index_temp].fission[index_pol][index_azi][gin];
+      val = xs[cache[tid].t].fission[cache[tid].p][cache[tid].a][gin];
     } else {
       val = 0.;
     }
     break;
   case MG_GET_XS_KAPPA_FISSION:
     if (fissionable) {
-      val = xs[index_temp].kappa_fission[index_pol][index_azi][gin];
+      val = xs[cache[tid].t].kappa_fission[cache[tid].p][cache[tid].a][gin];
     } else {
       val = 0.;
     }
     break;
   case MG_GET_XS_PROMPT_NU_FISSION:
     if (fissionable) {
-      val = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
+      val = xs[cache[tid].t].prompt_nu_fission[cache[tid].p][cache[tid].a][gin];
     } else {
       val = 0.;
     }
@@ -448,11 +454,11 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
   case MG_GET_XS_DELAYED_NU_FISSION:
     if (fissionable) {
       if (dg != nullptr) {
-        val = xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][*dg];
+        val = xs[cache[tid].t].delayed_nu_fission[cache[tid].p][cache[tid].a][gin][*dg];
       } else {
         val = 0.;
-        for (auto& num : xs[index_temp].delayed_nu_fission[index_pol]
-                                              [index_azi][gin]) {
+        for (auto& num : xs[cache[tid].t].delayed_nu_fission[cache[tid].p]
+             [cache[tid].a][gin]) {
           val += num;
         }
       }
@@ -462,7 +468,7 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
     break;
   case MG_GET_XS_NU_FISSION:
     if (fissionable) {
-      val = xs[index_temp].nu_fission[index_pol][index_azi][gin];
+      val = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
     } else {
       val = 0.;
     }
@@ -470,11 +476,11 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
   case MG_GET_XS_CHI_PROMPT:
     if (fissionable) {
       if (gout != nullptr) {
-        val = xs[index_temp].chi_prompt[index_pol][index_azi][gin][*gout];
+        val = xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][*gout];
       } else {
         // provide an outgoing group-wise sum
         val = 0.;
-        for (auto& num : xs[index_temp].chi_prompt[index_pol][index_azi][gin]) {
+        for (auto& num : xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin]) {
           val += num;
         }
       }
@@ -486,23 +492,23 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
     if (fissionable) {
       if (gout != nullptr) {
         if (dg != nullptr) {
-          val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][*dg];
+          val = xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][*gout][*dg];
         } else {
-          val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][0];
+          val = xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][*gout][0];
         }
       } else {
         if (dg != nullptr) {
           val = 0.;
-          for (int i = 0; i < xs[index_temp].chi_delayed[index_pol]
-                                            [index_azi][gin].size(); i++) {
-            val += xs[index_temp].chi_delayed[index_pol][index_azi][gin][i][*dg];
+          for (int i = 0; i < xs[cache[tid].t].chi_delayed[cache[tid].p]
+               [cache[tid].a][gin].size(); i++) {
+            val += xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][i][*dg];
           }
         } else {
           val = 0.;
-          for (int i = 0; i < xs[index_temp].chi_delayed[index_pol]
-                                            [index_azi][gin].size(); i++) {
-            for (auto& num : xs[index_temp].chi_delayed[index_pol]
-                                           [index_azi][gin][i]) {
+          for (int i = 0; i < xs[cache[tid].t].chi_delayed[cache[tid].p]
+                                            [cache[tid].a][gin].size(); i++) {
+            for (auto& num : xs[cache[tid].t].chi_delayed[cache[tid].p]
+                                           [cache[tid].a][gin][i]) {
               val += num;
             }
           }
@@ -519,14 +525,15 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
 }
 
 
-void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
+void Mgxs::sample_fission_energy(const int tid, const int gin, int& dg,
+     int& gout)
 {
   // This method assumes that the temperature and angle indices are set
-  double nu_fission = xs[index_temp].nu_fission[index_pol][index_azi][gin];
+  double nu_fission = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
 
   // Find the probability of having a prompt neutron
   double prob_prompt =
-       xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
+       xs[cache[tid].t].prompt_nu_fission[cache[tid].p][cache[tid].a][gin];
 
   // sample random numbers
   double xi_pd = prn() * nu_fission;
@@ -542,10 +549,10 @@ void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
     // sample the outgoing energy group
     gout = 0;
     double prob_gout =
-         xs[index_temp].chi_prompt[index_pol][index_azi][gin][gout];
+         xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][gout];
     while (prob_gout < xi_gout) {
       gout++;
-      prob_gout += xs[index_temp].chi_prompt[index_pol][index_azi][gin][gout];
+      prob_gout += xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][gout];
     }
 
   } else {
@@ -556,7 +563,7 @@ void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
     while (xi_pd >= prob_prompt) {
       dg++;
       prob_prompt +=
-           xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][dg];
+           xs[cache[tid].t].delayed_nu_fission[cache[tid].p][cache[tid].a][gin][dg];
     }
 
     // adjust dg in case of round-off error
@@ -565,35 +572,36 @@ void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
     // sample the outgoing energy group
     gout = 0;
     double prob_gout =
-         xs[index_temp].chi_delayed[index_pol][index_azi][gin][gout][dg];
+         xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][gout][dg];
     while (prob_gout < xi_gout) {
       gout++;
       prob_gout +=
-           xs[index_temp].chi_delayed[index_pol][index_azi][gin][gout][dg];
+           xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][gout][dg];
     }
   }
 }
 
 
-void Mgxs::sample_scatter(const int gin, int& gout, double& mu, double& wgt)
+void Mgxs::sample_scatter(const int tid, const int gin, int& gout, double& mu,
+     double& wgt)
 {
   // This method assumes that the temperature and angle indices are set
   // Sample the data
-  xs[index_temp].scatter[index_pol][index_azi]->sample(gin, gout, mu, wgt);
+  xs[cache[tid].t].scatter[cache[tid].p][cache[tid].a]->sample(gin, gout, mu, wgt);
 }
 
 
-void Mgxs::calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
-                        double& total_xs, double& abs_xs, double& nu_fiss_xs)
+void Mgxs::calculate_xs(const int tid, const int gin, const double sqrtkT,
+     const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs)
 {
   // Set our indices
-  set_temperature_index(sqrtkT);
-  set_angle_index(uvw);
-  total_xs = xs[index_temp].total[index_pol][index_azi][gin];
-  abs_xs = xs[index_temp].absorption[index_pol][index_azi][gin];
+  set_temperature_index(tid, sqrtkT);
+  set_angle_index(tid, uvw);
+  total_xs = xs[cache[tid].t].total[cache[tid].p][cache[tid].a][gin];
+  abs_xs = xs[cache[tid].t].absorption[cache[tid].p][cache[tid].a][gin];
 
   if (fissionable) {
-    nu_fiss_xs = xs[index_temp].nu_fission[index_pol][index_azi][gin];
+    nu_fiss_xs = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
   } else {
     nu_fiss_xs = 0.;
   }
@@ -617,10 +625,10 @@ bool Mgxs::equiv(const Mgxs& that)
 }
 
 
-void Mgxs::set_temperature_index(const double sqrtkT)
+void Mgxs::set_temperature_index(const int tid, const double sqrtkT)
 {
   // See if we need to find the new index
-  if (sqrtkT != last_sqrtkT) {
+  if (sqrtkT != cache[tid].sqrtkT) {
     double kT = sqrtkT * sqrtkT;
 
     // initialize vector for storage of the differences
@@ -628,34 +636,34 @@ void Mgxs::set_temperature_index(const double sqrtkT)
 
     // Find the minimum difference of kT and kTs
     temp_diff = std::abs(temp_diff - kT);
-    index_temp = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
+    cache[tid].t = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
          std::begin(temp_diff);
 
     // store this temperature as the last one used
-    last_sqrtkT = sqrtkT;
+    cache[tid].sqrtkT = sqrtkT;
   }
 }
 
 
-void Mgxs::set_angle_index(const double uvw[3])
+void Mgxs::set_angle_index(const int tid, const double uvw[3])
 {
   // See if we need to find the new index
-  if ((uvw[0] != last_uvw[0]) || (uvw[1] != last_uvw[1]) ||
-      (uvw[2] != last_uvw[2])) {
+  if ((uvw[0] != cache[tid].uvw[0]) || (uvw[1] != cache[tid].uvw[1]) ||
+      (uvw[2] != cache[tid].uvw[2])) {
     // convert uvw to polar and azimuthal angles
     double my_pol = std::acos(uvw[2]);
     double my_azi = std::atan2(uvw[1], uvw[0]);
 
     // Find the location, assuming equal-bin angles
     double delta_angle = PI / n_pol;
-    index_pol = std::floor(my_pol / delta_angle);
+    cache[tid].p = std::floor(my_pol / delta_angle);
     delta_angle = 2. * PI / n_azi;
-    index_azi = std::floor((my_azi + PI) / delta_angle);
+    cache[tid].a = std::floor((my_azi + PI) / delta_angle);
 
     // store this direction as the last one used
-    last_uvw[0] = uvw[0];
-    last_uvw[1] = uvw[1];
-    last_uvw[2] = uvw[2];
+    cache[tid].uvw[0] = uvw[0];
+    cache[tid].uvw[1] = uvw[1];
+    cache[tid].uvw[2] = uvw[2];
   }
 }
 
diff --git a/src/mgxs.h b/src/mgxs.h
index 2fc3e2bec..1c98417d6 100644
--- a/src/mgxs.h
+++ b/src/mgxs.h
@@ -23,6 +23,18 @@
 
 namespace openmc {
 
+//==============================================================================
+// Cache contains the cached data for an MGXS object
+//==============================================================================
+
+struct CacheData {
+  double sqrtkT; // last temperature corresponding to t
+  int t; // temperature index
+  int p; // polar angle index
+  int a; // azimuthal angle index
+  double uvw[3]; // last angle that corresponds to p and a
+};
+
 //==============================================================================
 // MGXS contains the mgxs data for a nuclide/material
 //==============================================================================
@@ -33,7 +45,6 @@ class Mgxs {
     int scatter_format; // flag for if this is legendre, histogram, or tabular
     int num_delayed_groups; // number of delayed neutron groups
     int num_groups;     // number of energy groups
-    double last_sqrtkT; // cache of the temperature corresponding to index_temp
     std::vector<XsData> xs; // Cross section data
     int n_pol;
     int n_azi;
@@ -42,7 +53,7 @@ class Mgxs {
     void _metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
          const int in_num_delayed_groups, double_1dvec& temperature,
          int& method, const double tolerance, int_1dvec& temps_to_read,
-         int& order_dim, bool& is_isotropic);
+         int& order_dim, bool& is_isotropic, const int n_threads);
     bool equiv(const Mgxs& that);
 
   public:
@@ -50,32 +61,31 @@ class Mgxs {
     double awr;         // atomic weight ratio
     bool fissionable;   // Is this fissionable
     // TODO: The following attributes be private when Fortran is fully replaced
-    int index_pol; // cache for the angle indices
-    int index_azi;
-    double last_uvw[3]; // cache of the angle corresponding to the above indices
-    int index_temp;     // cache of temperature index
+    std::vector<CacheData> cache; // index and data cache
     void init(const std::string& in_name, const double in_awr,
          const double_1dvec& in_kTs, const bool in_fissionable,
          const int in_scatter_format, const int in_num_groups,
          const int in_num_delayed_groups, const double_1dvec& in_polar,
-         const double_1dvec& in_azimuthal);
+         const double_1dvec& in_azimuthal, const int n_threads);
     void build_macro(const std::string& in_name, double_1dvec& mat_kTs,
                      std::vector<Mgxs*>& micros, double_1dvec& atom_densities,
-                     int& method, double tolerance);
+                     int& method, const double tolerance, const int n_threads);
     void combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
                  int_1dvec& micro_ts, int this_t);
-    void from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups,
-         double_1dvec& temperature, int& method, double tolerance,
-         int max_order, bool legendre_to_tabular,
-         int legendre_to_tabular_points);
-    double get_xs(const int xstype, const int gin, int* gout, double* mu,
-                  int* dg);
-    void sample_fission_energy(const int gin, int& dg, int& gout);
-    void sample_scatter(const int gin, int& gout, double& mu, double& wgt);
-    void calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
-         double& total_xs, double& abs_xs, double& nu_fiss_xs);
-    void set_temperature_index(const double sqrtkT);
-    void set_angle_index(const double uvw[3]);
+    void from_hdf5(hid_t xs_id, const int energy_groups,
+         const int delayed_groups, double_1dvec& temperature, int& method,
+         const double tolerance, const int max_order,
+         const bool legendre_to_tabular, const int legendre_to_tabular_points,
+         const int n_threads);
+    double get_xs(const int tid, const int xstype, const int gin, int* gout,
+         double* mu, int* dg);
+    void sample_fission_energy(const int tid, const int gin, int& dg, int& gout);
+    void sample_scatter(const int tid, const int gin, int& gout, double& mu,
+         double& wgt);
+    void calculate_xs(const int tid, const int gin, const double sqrtkT,
+         const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs);
+    void set_temperature_index(const int tid, const double sqrtkT);
+    void set_angle_index(const int tid, const double uvw[3]);
 };
 
 } // namespace openmc
diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90
index 8c48d5c98..ef416125e 100644
--- a/src/mgxs_data.F90
+++ b/src/mgxs_data.F90
@@ -2,6 +2,10 @@ module mgxs_data
 
   use, intrinsic :: ISO_C_BINDING
 
+#ifdef _OPENMP
+  use omp_lib
+#endif
+
   use constants
   use algorithm,       only: find
   use dict_header,     only: DictCharInt
@@ -36,6 +40,13 @@ contains
     type(VectorReal), allocatable, target :: temps(:)
     character(MAX_WORD_LEN) :: word
     integer, allocatable    :: array(:)
+    integer(C_INT) :: n_threads
+
+#ifdef _OPENMP
+    n_threads = OMP_GET_MAX_THREADS()
+#else
+    n_threads = 1
+#endif
 
     ! Check if MGXS Library exists
     inquire(FILE=path_cross_sections, EXIST=file_exists)
@@ -84,7 +95,8 @@ contains
                num_energy_groups, num_delayed_groups, &
                temps(i_nuclide) % size(), temps(i_nuclide) % data, &
                temperature_method, temperature_tolerance, max_order, &
-               logical(legendre_to_tabular, C_BOOL), legendre_to_tabular_points)
+               logical(legendre_to_tabular, C_BOOL), &
+               legendre_to_tabular_points, n_threads)
 
           call already_read % add(name)
         end if
@@ -110,6 +122,13 @@ contains
     type(Material), pointer       :: mat   ! current material
     type(VectorReal), allocatable :: kTs(:)
     character(MAX_WORD_LEN)       :: name  ! name of material
+    integer(C_INT) :: n_threads
+
+#ifdef _OPENMP
+    n_threads = OMP_GET_MAX_THREADS()
+#else
+    n_threads = 1
+#endif
 
     ! Get temperatures to read for each material
     call get_mat_kTs(kTs)
@@ -130,7 +149,7 @@ contains
       if (allocated(kTs(i_mat) % data)) then
         call create_macro_xs_c(name, mat % n_nuclides, mat % nuclide, &
              kTs(i_mat) % size(), kTs(i_mat) % data, mat % atom_density, &
-             temperature_method, temperature_tolerance)
+             temperature_method, temperature_tolerance, n_threads)
       end if
     end do
 
diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90
index 7404af686..34b0fb611 100644
--- a/src/mgxs_interface.F90
+++ b/src/mgxs_interface.F90
@@ -10,7 +10,7 @@ module mgxs_interface
 
     subroutine add_mgxs_c(file_id, name, energy_groups, delayed_groups, &
          n_temps, temps, method, tolerance, max_order, legendre_to_tabular, &
-         legendre_to_tabular_points) bind(C)
+         legendre_to_tabular_points, n_threads) bind(C)
       use ISO_C_BINDING
       import HID_T
       implicit none
@@ -25,6 +25,7 @@ module mgxs_interface
       integer(C_INT), value, intent(in) :: max_order
       logical(C_BOOL),value, intent(in) :: legendre_to_tabular
       integer(C_INT), value, intent(in) :: legendre_to_tabular_points
+      integer(C_INT), value, intent(in) :: n_threads
     end subroutine add_mgxs_c
 
     function query_fissionable_c(n_nuclides, i_nuclides) result(result) bind(C)
@@ -36,7 +37,7 @@ module mgxs_interface
     end function query_fissionable_c
 
     subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, &
-         atom_densities, method, tolerance) bind(C)
+         atom_densities, method, tolerance, n_threads) bind(C)
       use ISO_C_BINDING
       implicit none
       character(kind=C_CHAR),intent(in) :: name(*)
@@ -47,13 +48,15 @@ module mgxs_interface
       real(C_DOUBLE),        intent(in) :: atom_densities(1:n_nuclides)
       integer(C_INT),     intent(inout) :: method
       real(C_DOUBLE), value, intent(in) :: tolerance
+      integer(C_INT), value, intent(in) :: n_threads
     end subroutine create_macro_xs_c
 
-    subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, &
+    subroutine calculate_xs_c(i_mat, tid, gin, sqrtkT, uvw, total_xs, abs_xs, &
          nu_fiss_xs) bind(C)
       use ISO_C_BINDING
       implicit none
       integer(C_INT), value, intent(in) :: i_mat
+      integer(C_INT), value, intent(in) :: tid
       integer(C_INT), value, intent(in) :: gin
       real(C_DOUBLE), value, intent(in) :: sqrtkT
       real(C_DOUBLE),        intent(in) :: uvw(1:3)
@@ -62,10 +65,11 @@ module mgxs_interface
       real(C_DOUBLE),     intent(inout) :: nu_fiss_xs
     end subroutine calculate_xs_c
 
-    subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
+    subroutine sample_scatter_c(i_mat, tid, gin, gout, mu, wgt, uvw) bind(C)
       use ISO_C_BINDING
       implicit none
       integer(C_INT), value, intent(in) :: i_mat
+      integer(C_INT), value, intent(in) :: tid
       integer(C_INT), value, intent(in) :: gin
       integer(C_INT),     intent(inout) :: gout
       real(C_DOUBLE),     intent(inout) :: mu
@@ -73,10 +77,11 @@ module mgxs_interface
       real(C_DOUBLE),     intent(inout) :: uvw(1:3)
     end subroutine sample_scatter_c
 
-    subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
+    subroutine sample_fission_energy_c(i_mat, tid, gin, dg, gout) bind(C)
       use ISO_C_BINDING
       implicit none
       integer(C_INT), value, intent(in) :: i_mat
+      integer(C_INT), value, intent(in) :: tid
       integer(C_INT), value, intent(in) :: gin
       integer(C_INT),     intent(inout) :: dg
       integer(C_INT),     intent(inout) :: gout
@@ -97,11 +102,12 @@ module mgxs_interface
       real(C_DOUBLE)                    :: awr
     end function get_awr_c
 
-    function get_nuclide_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
+    function get_nuclide_xs_c(index, tid, xstype, gin, gout, mu, dg) result(val) &
          bind(C)
       use ISO_C_BINDING
       implicit none
       integer(C_INT),  value,   intent(in) :: index
+      integer(C_INT),  value,   intent(in) :: tid
       integer(C_INT),  value,   intent(in) :: xstype
       integer(C_INT),  value,   intent(in) :: gin
       integer(C_INT), optional, intent(in) :: gout
@@ -110,11 +116,12 @@ module mgxs_interface
       real(C_DOUBLE)                       :: val
     end function get_nuclide_xs_c
 
-    function get_macro_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
+    function get_macro_xs_c(index, tid, xstype, gin, gout, mu, dg) result(val) &
          bind(C)
       use ISO_C_BINDING
       implicit none
       integer(C_INT),  value,   intent(in) :: index
+      integer(C_INT),  value,   intent(in) :: tid
       integer(C_INT),  value,   intent(in) :: xstype
       integer(C_INT),  value,   intent(in) :: gin
       integer(C_INT), optional, intent(in) :: gout
@@ -123,63 +130,29 @@ module mgxs_interface
       real(C_DOUBLE)                       :: val
     end function get_macro_xs_c
 
-    subroutine set_nuclide_angle_index_c(index, uvw, last_pol, last_azi, &
-         last_uvw) bind(C)
+    subroutine set_nuclide_angle_index_c(index, tid, uvw) bind(C)
       use ISO_C_BINDING
       implicit none
       integer(C_INT), value, intent(in) :: index
+      integer(C_INT), value, intent(in) :: tid
       real(C_DOUBLE),        intent(in) :: uvw(1:3)
-      integer(C_INT),     intent(inout) :: last_pol
-      integer(C_INT),     intent(inout) :: last_azi
-      real(C_DOUBLE),     intent(inout) :: last_uvw(1:3)
     end subroutine set_nuclide_angle_index_c
 
-    subroutine reset_nuclide_angle_index_c(index, last_pol, last_azi, &
-         last_uvw) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: index
-      integer(C_INT), value, intent(in) :: last_pol
-      integer(C_INT), value, intent(in) :: last_azi
-      real(C_DOUBLE),        intent(in) :: last_uvw(1:3)
-    end subroutine reset_nuclide_angle_index_c
-
-    subroutine set_macro_angle_index_c(index, uvw, last_pol, last_azi, &
-         last_uvw) bind(C)
+    subroutine set_macro_angle_index_c(index, tid, uvw) bind(C)
       use ISO_C_BINDING
       implicit none
       integer(C_INT), value, intent(in) :: index
+      integer(C_INT), value, intent(in) :: tid
       real(C_DOUBLE),        intent(in) :: uvw(1:3)
-      integer(C_INT),     intent(inout) :: last_pol
-      integer(C_INT),     intent(inout) :: last_azi
-      real(C_DOUBLE),     intent(inout) :: last_uvw(1:3)
     end subroutine set_macro_angle_index_c
 
-    subroutine reset_macro_angle_index_c(index, last_pol, last_azi, &
-         last_uvw) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: index
-      integer(C_INT), value, intent(in) :: last_pol
-      integer(C_INT), value, intent(in) :: last_azi
-      real(C_DOUBLE),        intent(in) :: last_uvw(1:3)
-    end subroutine reset_macro_angle_index_c
-
-    function set_nuclide_temperature_index_c(index, sqrtkT) result(last_temp) &
-         bind(C)
+    subroutine set_nuclide_temperature_index_c(index, tid, sqrtkT) bind(C)
       use ISO_C_BINDING
       implicit none
       integer(C_INT), value, intent(in) :: index
+      integer(C_INT), value, intent(in) :: tid
       real(C_DOUBLE), value, intent(in) :: sqrtkT
-      integer(C_INT)                    :: last_temp
-    end function set_nuclide_temperature_index_c
-
-    subroutine reset_nuclide_temperature_index_c(index, last_temp) bind(C)
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: index
-      integer(C_INT), value, intent(in) :: last_temp
-    end subroutine reset_nuclide_temperature_index_c
+    end subroutine set_nuclide_temperature_index_c
 
   end interface
 
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index 60afc2848..6c9dc9f42 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -6,10 +6,11 @@ namespace openmc {
 // Mgxs data loading interface methods
 //==============================================================================
 
-void add_mgxs_c(hid_t file_id, char* name, int energy_groups,
-     int delayed_groups, int n_temps, double temps[], int& method,
-     double tolerance, int max_order, bool legendre_to_tabular,
-     int legendre_to_tabular_points)
+void add_mgxs_c(hid_t file_id, char* name, const int energy_groups,
+     const int delayed_groups, const int n_temps, double temps[], int& method,
+     const double tolerance, const int max_order,
+     const bool legendre_to_tabular, const int legendre_to_tabular_points,
+     const int n_threads)
 {
   //!! mgxs_data.F90 will be modified to just create the list of names
   //!! in the order needed
@@ -32,7 +33,7 @@ void add_mgxs_c(hid_t file_id, char* name, int energy_groups,
   Mgxs mg;
   mg.from_hdf5(xs_grp, energy_groups, delayed_groups,
        temperature, method, tolerance, max_order, legendre_to_tabular,
-       legendre_to_tabular_points);
+       legendre_to_tabular_points, n_threads);
 
   nuclides_MG.push_back(mg);
 }
@@ -50,7 +51,8 @@ bool query_fissionable_c(const int n_nuclides, const int i_nuclides[])
 
 void create_macro_xs_c(char* mat_name, const int n_nuclides,
      const int i_nuclides[], const int n_temps, const double temps[],
-     const double atom_densities[], int& method, const double tolerance)
+     const double atom_densities[], int& method, const double tolerance,
+     const int n_threads)
 {
   Mgxs macro;
   if (n_temps > 0) {
@@ -70,7 +72,7 @@ void create_macro_xs_c(char* mat_name, const int n_nuclides,
     }
 
     macro.build_macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
-         method, tolerance);
+         method, tolerance, n_threads);
   }
   macro_xs.push_back(macro);
 }
@@ -79,19 +81,20 @@ void create_macro_xs_c(char* mat_name, const int n_nuclides,
 // Mgxs tracking/transport/tallying interface methods
 //==============================================================================
 
-void calculate_xs_c(const int i_mat, const int gin, const double sqrtkT,
-     const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs)
+void calculate_xs_c(const int i_mat, const int tid, const int gin,
+     const double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs,
+     double& nu_fiss_xs)
 {
-  macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
+  macro_xs[i_mat - 1].calculate_xs(tid, gin - 1, sqrtkT, uvw, total_xs, abs_xs,
        nu_fiss_xs);
 }
 
 
-void sample_scatter_c(const int i_mat, const int gin, int& gout, double& mu,
-     double& wgt, double uvw[3])
+void sample_scatter_c(const int i_mat, const int tid, const int gin, int& gout,
+     double& mu, double& wgt, double uvw[3])
 {
   int gout_c = gout - 1;
-  macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
+  macro_xs[i_mat - 1].sample_scatter(tid, gin - 1, gout_c, mu, wgt);
 
   // adjust return value for fortran indexing
   gout = gout_c + 1;
@@ -101,11 +104,12 @@ void sample_scatter_c(const int i_mat, const int gin, int& gout, double& mu,
 }
 
 
-void sample_fission_energy_c(const int i_mat, const int gin, int& dg, int& gout)
+void sample_fission_energy_c(const int i_mat, const int tid, const int gin,
+     int& dg, int& gout)
 {
   int dg_c = 0;
   int gout_c = 0;
-  macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
+  macro_xs[i_mat - 1].sample_fission_energy(tid, gin - 1, dg_c, gout_c);
 
   // adjust return values for fortran indexing
   dg = dg_c + 1;
@@ -113,6 +117,82 @@ void sample_fission_energy_c(const int i_mat, const int gin, int& dg, int& gout)
 }
 
 
+double get_nuclide_xs_c(const int index, const int tid, const int xstype,
+     const int gin, int* gout, double* mu, int* dg)
+{
+  int gout_c;
+  int* gout_c_p;
+  int dg_c;
+  int* dg_c_p;
+  if (gout != nullptr) {
+    gout_c = *gout - 1;
+    gout_c_p = &gout_c;
+  } else {
+    gout_c_p = gout;
+  }
+  if (dg != nullptr) {
+    dg_c = *dg - 1;
+    dg_c_p = &dg_c;
+  } else {
+    dg_c_p = dg;
+  }
+  return nuclides_MG[index - 1].get_xs(tid, xstype, gin - 1, gout_c_p, mu,
+                                       dg_c_p);
+}
+
+
+double get_macro_xs_c(const int index, const int tid, const int xstype,
+     const int gin, int* gout, double* mu, int* dg)
+{
+  int gout_c;
+  int* gout_c_p;
+  int dg_c;
+  int* dg_c_p;
+  if (gout != nullptr) {
+    gout_c = *gout - 1;
+    gout_c_p = &gout_c;
+  } else {
+    gout_c_p = gout;
+  }
+  if (dg != nullptr) {
+    dg_c = *dg - 1;
+    dg_c_p = &dg_c;
+  } else {
+    dg_c_p = dg;
+  }
+  return macro_xs[index - 1].get_xs(tid, xstype, gin - 1, gout_c_p, mu,
+                                    dg_c_p);
+}
+
+
+void set_nuclide_angle_index_c(const int index, const int tid,
+     const double uvw[3])
+{
+  // Update the values
+  nuclides_MG[index - 1].set_angle_index(tid, uvw);
+}
+
+
+void set_macro_angle_index_c(const int index, const int tid,
+     const double uvw[3])
+{
+  // Update the values
+  macro_xs[index - 1].set_angle_index(tid, uvw);
+}
+
+
+void set_nuclide_temperature_index_c(const int index, const int tid,
+     const double sqrtkT)
+{
+  // Update the values
+  nuclides_MG[index - 1].set_temperature_index(tid, sqrtkT);
+}
+
+
+//==============================================================================
+// Mgxs general methods
+//==============================================================================
+
 void get_name_c(const int index, int name_len, char* name)
 {
   // First blank out our input string
@@ -133,115 +213,4 @@ double get_awr_c(const int index)
   return nuclides_MG[index - 1].awr;
 }
 
-
-double get_nuclide_xs_c(const int index, const int xstype, const int gin,
-     int* gout, double* mu, int* dg)
-{
-  int gout_c;
-  int* gout_c_p;
-  int dg_c;
-  int* dg_c_p;
-  if (gout != nullptr) {
-    gout_c = *gout - 1;
-    gout_c_p = &gout_c;
-  } else {
-    gout_c_p = gout;
-  }
-  if (dg != nullptr) {
-    dg_c = *dg - 1;
-    dg_c_p = &dg_c;
-  } else {
-    dg_c_p = dg;
-  }
-  return nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
-}
-
-
-double get_macro_xs_c(const int index, const int xstype, const int gin,
-     int* gout, double* mu, int* dg)
-{
-  int gout_c;
-  int* gout_c_p;
-  int dg_c;
-  int* dg_c_p;
-  if (gout != nullptr) {
-    gout_c = *gout - 1;
-    gout_c_p = &gout_c;
-  } else {
-    gout_c_p = gout;
-  }
-  if (dg != nullptr) {
-    dg_c = *dg - 1;
-    dg_c_p = &dg_c;
-  } else {
-    dg_c_p = dg;
-  }
-  return macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
-}
-
-
-void set_nuclide_angle_index_c(const int index, const double uvw[3],
-     int& last_pol, int& last_azi, double last_uvw[3])
-{
-  // Store the old
-  last_pol = nuclides_MG[index - 1].index_pol;
-  last_azi = nuclides_MG[index - 1].index_azi;
-  last_uvw[0] = nuclides_MG[index - 1].last_uvw[0];
-  last_uvw[1] = nuclides_MG[index - 1].last_uvw[1];
-  last_uvw[2] = nuclides_MG[index - 1].last_uvw[2];
-
-  // Update the values
-  nuclides_MG[index - 1].set_angle_index(uvw);
-}
-
-
-void reset_nuclide_angle_index_c(const int index, const int last_pol,
-     const int last_azi, const double last_uvw[3])
-{
-  nuclides_MG[index - 1].index_pol = last_pol;
-  nuclides_MG[index - 1].index_azi = last_azi;
-  nuclides_MG[index - 1].last_uvw[0] = last_uvw[0];
-  nuclides_MG[index - 1].last_uvw[1] = last_uvw[1];
-  nuclides_MG[index - 1].last_uvw[2] = last_uvw[2];
-}
-
-
-void set_macro_angle_index_c(const int index, const double uvw[3],
-     int& last_pol, int& last_azi, double last_uvw[3])
-{
-  // Store the old
-  last_pol = macro_xs[index - 1].index_pol;
-  last_azi = macro_xs[index - 1].index_azi;
-  last_uvw[0] = macro_xs[index - 1].last_uvw[0];
-  last_uvw[1] = macro_xs[index - 1].last_uvw[1];
-  last_uvw[2] = macro_xs[index - 1].last_uvw[2];
-
-  // Update the values
-  macro_xs[index - 1].set_angle_index(uvw);
-}
-
-
-void reset_macro_angle_index_c(const int index, const int last_pol,
-     const int last_azi, const double last_uvw[3])
-{
-  macro_xs[index - 1].index_pol = last_pol;
-  macro_xs[index - 1].index_azi = last_azi;
-  macro_xs[index - 1].last_uvw[0] = last_uvw[0];
-  macro_xs[index - 1].last_uvw[1] = last_uvw[1];
-  macro_xs[index - 1].last_uvw[2] = last_uvw[2];
-}
-
-
-int set_nuclide_temperature_index_c(const int index, const double sqrtkT)
-{
-  int old = nuclides_MG[index - 1].index_temp;
-  nuclides_MG[index - 1].set_temperature_index(sqrtkT);
-  return old;
-}
-
-void reset_nuclide_temperature_index_c(const int index, const int last_temp)
-{
-  nuclides_MG[index - 1].index_temp = last_temp;
-}
-
 } // namespace openmc
\ No newline at end of file
diff --git a/src/mgxs_interface.h b/src/mgxs_interface.h
index 197ecb0bf..4e1068266 100644
--- a/src/mgxs_interface.h
+++ b/src/mgxs_interface.h
@@ -13,54 +13,47 @@ extern std::vector<Mgxs> nuclides_MG;
 extern std::vector<Mgxs> macro_xs;
 
 
-extern "C" void add_mgxs_c(hid_t file_id, char* name, int energy_groups,
-     int delayed_groups, int n_temps, double temps[], int& method,
-     double tolerance, int max_order, bool legendre_to_tabular,
-     int legendre_to_tabular_points);
+extern "C" void add_mgxs_c(hid_t file_id, char* name, const int energy_groups,
+     const int delayed_groups, const int n_temps, double temps[], int& method,
+     const double tolerance, const int max_order,
+     const bool legendre_to_tabular, const int legendre_to_tabular_points,
+     const int n_threads);
 
 extern "C" bool query_fissionable_c(const int n_nuclides, const int i_nuclides[]);
 
 extern "C" void create_macro_xs_c(char* mat_name, const int n_nuclides,
      const int i_nuclides[], const int n_temps, const double temps[],
-     const double atom_densities[], int& method, const double tolerance);
+     const double atom_densities[], int& method, const double tolerance,
+     const int n_threads);
 
-extern "C" void calculate_xs_c(const int i_mat, const int gin,
+extern "C" void calculate_xs_c(const int i_mat, const int tid, const int gin,
      const double sqrtkT, const double uvw[3], double& total_xs,
      double& abs_xs, double& nu_fiss_xs);
 
-extern "C" void sample_scatter_c(const int i_mat, const int gin, int& gout,
-     double& mu, double& wgt, double uvw[3]);
+extern "C" void sample_scatter_c(const int i_mat, const int tid, const int gin,
+     int& gout, double& mu, double& wgt, double uvw[3]);
 
-extern "C" void sample_fission_energy_c(const int i_mat, const int gin,
-     int& dg, int& gout);
+extern "C" void sample_fission_energy_c(const int i_mat, const int tid,
+     const int gin, int& dg, int& gout);
 
 extern "C" void get_name_c(const int index, int name_len, char* name);
 
 extern "C" double get_awr_c(const int index);
 
-extern "C" double get_nuclide_xs_c(const int index, const int xstype,
-     const int gin, int* gout, double* mu, int* dg);
+extern "C" double get_nuclide_xs_c(const int index, const int tid,
+     const int xstype, const int gin, int* gout, double* mu, int* dg);
 
-extern "C" double get_macro_xs_c(const int index, const int xstype,
-     const int gin, int* gout, double* mu, int* dg);
+extern "C" double get_macro_xs_c(const int index, const int tid,
+     const int xstype, const int gin, int* gout, double* mu, int* dg);
 
-extern "C" void set_nuclide_angle_index_c(const int index, const double uvw[3],
-     int& last_pol, int& last_azi, double last_uvw[3]);
+extern "C" void set_nuclide_angle_index_c(const int index, const int tid,
+     const double uvw[3]);
 
-extern "C" void reset_nuclide_angle_index_c(const int index, const int last_pol,
-     const int last_azi, const double last_uvw[3]);
+extern "C" void set_macro_angle_index_c(const int index, const int tid,
+     const double uvw[3]);
 
-extern "C" void set_macro_angle_index_c(const int index, const double uvw[3],
-     int& last_pol, int& last_azi, double last_uvw[3]);
-
-extern "C" void reset_macro_angle_index_c(const int index, const int last_pol,
-     const int last_azi, const double last_uvw[3]);
-
-extern "C" int set_nuclide_temperature_index_c(const int index,
+extern "C" void set_nuclide_temperature_index_c(const int index, const int tid,
      const double sqrtkT);
 
-extern "C" void reset_nuclide_temperature_index_c(const int index,
-     const int last_temp);
-
 } // namespace openmc
 #endif // MGXS_INTERFACE_H
\ No newline at end of file
diff --git a/src/physics_mg.F90 b/src/physics_mg.F90
index 797514e25..d04d4cb85 100644
--- a/src/physics_mg.F90
+++ b/src/physics_mg.F90
@@ -2,6 +2,10 @@ module physics_mg
   ! This module contains the multi-group specific physics routines so as to not
   ! hinder performance of the CE versions with multiple if-thens.
 
+#ifdef _OPENMP
+  use omp_lib
+#endif
+
   use bank_header
   use constants
   use error,                  only: fatal_error, warning, write_message
@@ -141,9 +145,15 @@ contains
   subroutine scatter(p)
 
     type(Particle), intent(inout)  :: p
+    integer(C_INT) :: tid
+#ifdef _OPENMP
+    tid = OMP_GET_THREAD_NUM()
+#else
+    tid = 0
+#endif
 
-    call sample_scatter_c(p % material, p % last_g, p % g, p % mu, p % wgt, &
-         p % coord(1) % uvw)
+    call sample_scatter_c(p % material, tid, p % last_g, p % g, p % mu, &
+         p % wgt, p % coord(1) % uvw)
 
     ! Update energy value for downstream compatability (in tallying)
     p % E = energy_bin_avg(p % g)
@@ -173,6 +183,12 @@ contains
     real(8) :: mu                       ! fission neutron angular cosine
     real(8) :: phi                      ! fission neutron azimuthal angle
     real(8) :: weight                   ! weight adjustment for ufs method
+    integer(C_INT) :: tid
+#ifdef _OPENMP
+    tid = OMP_GET_THREAD_NUM()
+#else
+    tid = 0
+#endif
 
     ! TODO: Heat generation from fission
 
@@ -252,7 +268,7 @@ contains
 
       ! Sample secondary energy distribution for fission reaction and set energy
       ! in fission bank
-      call sample_fission_energy_c(p % material, p % g, dg, gout)
+      call sample_fission_energy_c(p % material, tid, p % g, dg, gout)
 
       bank_array(i) % E = real(gout, 8)
       bank_array(i) % delayed_group = dg
diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90
index 4c24717b4..f6839081d 100644
--- a/src/tallies/tally.F90
+++ b/src/tallies/tally.F90
@@ -2,6 +2,10 @@ module tally
 
   use, intrinsic :: ISO_C_BINDING
 
+#ifdef _OPENMP
+  use omp_lib
+#endif
+
   use algorithm,        only: binary_search
   use constants
   use dict_header,      only: EMPTY
@@ -1229,14 +1233,12 @@ contains
     real(8) :: p_uvw(3)             ! Particle's current uvw
     integer :: p_g                  ! Particle group to use for getting info
                                     ! to tally with.
-    ! Storage of the indices the Mgxs object arrived with for resetting later
-    integer(C_INT) :: last_nuc_azi
-    integer(C_INT) :: last_nuc_pol
-    integer(C_INT) :: last_mat_azi
-    integer(C_INT) :: last_mat_pol
-    integer(C_INT) :: last_nuc_temp
-    real(C_DOUBLE) :: last_mat_uvw(3)
-    real(C_DOUBLE) :: last_nuc_uvw(3)
+    integer(C_INT) :: tid
+#ifdef _OPENMP
+    tid = OMP_GET_THREAD_NUM()
+#else
+    tid = 0
+#endif
 
     ! Set the direction and group to use with get_xs
     if (t % estimator == ESTIMATOR_ANALOG .or. &
@@ -1274,15 +1276,13 @@ contains
 
     ! To significantly reduce de-referencing, point matxs to the
     ! macroscopic Mgxs for the material of interest
-    call set_macro_angle_index_c(p % material, p_uvw, last_mat_pol, &
-                                 last_mat_azi, last_mat_uvw)
+    call set_macro_angle_index_c(p % material, tid, p_uvw)
 
     ! Do same for nucxs, point it to the microscopic nuclide data of interest
     if (i_nuclide > 0) then
       ! And since we haven't calculated this temperature index yet, do so now
-      last_nuc_temp = set_nuclide_temperature_index_c(i_nuclide, p % sqrtkT)
-      call set_nuclide_angle_index_c(i_nuclide, p_uvw, last_nuc_pol, &
-                                     last_nuc_azi, last_nuc_uvw)
+      call set_nuclide_temperature_index_c(i_nuclide, tid, p % sqrtkT)
+      call set_nuclide_angle_index_c(i_nuclide, tid, p_uvw)
     end if
 
     i = 0
@@ -1336,14 +1336,14 @@ contains
           end if
 
           if (i_nuclide > 0) then
-            score = score * atom_density * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_TOTAL, p_g) * flux
+            score = score * flux * atom_density * &
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_TOTAL, p_g) / &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_TOTAL, p_g)
           end if
 
         else
           if (i_nuclide > 0) then
-            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) * &
+            score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_TOTAL, p_g) * &
                  atom_density * flux
           else
             score = material_xs % total * flux
@@ -1366,22 +1366,22 @@ contains
           end if
 
           if (i_nuclide > 0) then
-            score = score * get_nuclide_xs_c(i_nuclide, &
+            score = score * flux * get_nuclide_xs_c(i_nuclide, tid, &
                  MG_GET_XS_INVERSE_VELOCITY, p_g) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) * flux
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
           else
-            score = score * get_macro_xs_c(p % material, &
+            score = score * flux * get_macro_xs_c(p % material, tid, &
                  MG_GET_XS_INVERSE_VELOCITY, p_g) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) * flux
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
           end if
 
         else
 
           if (i_nuclide > 0) then
-            score = flux * get_nuclide_xs_c(i_nuclide, &
+            score = flux * get_nuclide_xs_c(i_nuclide, tid, &
                  MG_GET_XS_INVERSE_VELOCITY, p_g)
           else
-            score = flux * get_macro_xs_c(p % material, &
+            score = flux * get_macro_xs_c(p % material, tid, &
                  MG_GET_XS_INVERSE_VELOCITY, p_g)
           end if
         end if
@@ -1404,22 +1404,22 @@ contains
           ! adjust the score by the actual probability for that nuclide.
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT, &
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER_FMU_MULT, &
                                   p % last_g, p % g, MU=p % mu) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU_MULT, &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER_FMU_MULT, &
                                 p % last_g, p % g, MU=p % mu)
           end if
 
         else
           if (i_nuclide > 0) then
             score = atom_density * flux * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_MULT, &
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER_MULT, &
                                   p_g, MU=p % mu)
           else
             ! Get the scattering x/s and take away
             ! the multiplication baked in to sigS
             score = flux * &
-                 get_macro_xs_c(p % material, MG_GET_XS_SCATTER_MULT, &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER_MULT, &
                                 p_g, MU=p % mu)
           end if
         end if
@@ -1442,20 +1442,20 @@ contains
           ! adjust the score by the actual probability for that nuclide.
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU, &
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER_FMU, &
                                 p % last_g, p % g, MU=p % mu) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU, &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER_FMU, &
                                 p % last_g, p % g, MU=p % mu)
           end if
 
         else
           if (i_nuclide > 0) then
             score = atom_density * flux * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER, p_g)
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER, p_g)
           else
             ! Get the scattering x/s, which includes multiplication
             score = flux * &
-                   get_macro_xs_c(p % material, MG_GET_XS_SCATTER, p_g)
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER, p_g)
           end if
         end if
 
@@ -1475,13 +1475,13 @@ contains
           end if
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_ABSORPTION, p_g) / &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
           end if
         else
           if (i_nuclide > 0) then
             score = atom_density * flux * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g)
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_ABSORPTION, p_g)
           else
             score = material_xs % absorption * flux
           end if
@@ -1506,19 +1506,19 @@ contains
           end if
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
           else
             score = score * &
-                 get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g) / &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
           end if
         else
           if (i_nuclide > 0) then
-            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) * &
+            score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) * &
                  atom_density * flux
           else
-            score = get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
+            score = get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g) * flux
           end if
         end if
 
@@ -1544,12 +1544,12 @@ contains
             score = p % absorb_wgt * flux
             if (i_nuclide > 0) then
               score = score * atom_density * &
-                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) / &
-                   get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                   get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_NU_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
             else
               score = score * &
-                   get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) / &
-                   get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_NU_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
             end if
           else
             ! Skip any non-fission events
@@ -1562,17 +1562,17 @@ contains
             score = keff * p % wgt_bank * flux
             if (i_nuclide > 0) then
               score = score * atom_density * &
-                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
-                   get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
+                   get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g)
             end if
           end if
 
         else
           if (i_nuclide > 0) then
-            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * &
+            score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_NU_FISSION, p_g) * &
                  atom_density * flux
           else
-            score = get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) * flux
+            score = get_macro_xs_c(p % material, tid, MG_GET_XS_NU_FISSION, p_g) * flux
           end if
         end if
 
@@ -1598,12 +1598,12 @@ contains
             score = p % absorb_wgt * flux
             if (i_nuclide > 0) then
               score = score * atom_density * &
-                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
-                   get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                   get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
             else
               score = score * &
-                   get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
-                   get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
             end if
           else
             ! Skip any non-fission events
@@ -1617,17 +1617,17 @@ contains
                  / real(p % n_bank, 8)) * flux
             if (i_nuclide > 0) then
               score = score * atom_density * &
-                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
-                   get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
+                   get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g)
             end if
           end if
 
         else
           if (i_nuclide > 0) then
-            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
+            score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
                  atom_density * flux
           else
-            score = get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
+            score = get_macro_xs_c(p % material, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
           end if
         end if
 
@@ -1653,7 +1653,7 @@ contains
             ! No fission events occur if survival biasing is on -- need to
             ! calculate fraction of absorptions that would have resulted in
             ! nu-fission
-            if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
+            if (get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
 
               if (dg_filter > 0) then
                 select type(filt => filters(t % filter(dg_filter)) % obj)
@@ -1669,12 +1669,12 @@ contains
                     score = p % absorb_wgt * flux
                     if (i_nuclide > 0) then
                       score = score * &
-                           get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
-                           get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                           get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                           get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
                     else
                       score = score * &
-                           get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
-                           get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                           get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                           get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
                     end if
 
                     call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1685,12 +1685,12 @@ contains
                 score = p % absorb_wgt * flux
                 if (i_nuclide > 0) then
                   score = score * &
-                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
-                       get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                       get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
+                       get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
                 else
                   score = score * &
-                       get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
-                       get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                       get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
+                       get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
                 end if
               end if
             end if
@@ -1720,8 +1720,8 @@ contains
 
                   if (i_nuclide > 0) then
                     score = score * atom_density * &
-                         get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
-                         get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
+                         get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
+                         get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g)
                   end if
 
                   call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1732,8 +1732,8 @@ contains
               score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank) * flux
               if (i_nuclide > 0) then
                 score = score * atom_density * &
-                     get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
-                     get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
+                     get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
+                     get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g)
               end if
             end if
           end if
@@ -1753,10 +1753,10 @@ contains
 
                 if (i_nuclide > 0) then
                   score = atom_density * flux * &
-                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
+                       get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
                 else
                   score = flux * &
-                       get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
+                       get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
                 end if
 
                 call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1766,11 +1766,11 @@ contains
           else
             if (i_nuclide > 0) then
               score = atom_density * flux * &
-                   get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g)
+                   get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g)
 
             else
               score = flux * &
-                   get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g)
+                   get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g)
             end if
           end if
         end if
@@ -1785,7 +1785,7 @@ contains
             ! No fission events occur if survival biasing is on -- need to
             ! calculate fraction of absorptions that would have resulted in
             ! nu-fission
-            if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
+            if (get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
 
               if (dg_filter > 0) then
                 select type(filt => filters(t % filter(dg_filter)) % obj)
@@ -1801,14 +1801,14 @@ contains
                     score = p % absorb_wgt * flux
                     if (i_nuclide > 0) then
                       score = score * &
-                           get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                           get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
-                           get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                           get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                           get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                           get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
                     else
                       score = score * &
-                           get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                           get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
-                           get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                           get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                           get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                           get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
                     end if
 
                     call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1826,14 +1826,14 @@ contains
                 do d = 1, num_delayed_groups
                   if (i_nuclide > 0) then
                     score = score + p % absorb_wgt * flux * &
-                         get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                         get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
-                         get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                         get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                         get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                         get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
                   else
                     score = score + p % absorb_wgt * flux * &
-                         get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                         get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
-                         get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                         get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                         get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
+                         get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
                   end if
                 end do
               end if
@@ -1862,13 +1862,13 @@ contains
                 if (i_nuclide > 0) then
                   score = score + keff * atom_density * &
                        fission_bank(n_bank - p % n_bank + k) % wgt * &
-                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
-                       get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
+                       get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                       get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
+                       get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g) * flux
                 else
                   score = score + keff * &
                        fission_bank(n_bank - p % n_bank + k) % wgt * &
-                       get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
+                       get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
                 end if
 
                 ! if the delayed group filter is present, tally to corresponding
@@ -1921,12 +1921,12 @@ contains
 
                 if (i_nuclide > 0) then
                   score = atom_density * flux * &
-                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                       get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
+                       get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                       get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
                 else
                   score = flux * &
-                       get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                       get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
+                       get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                       get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
                 end if
 
                 call score_fission_delayed_dg(t, d_bin, score, score_index)
@@ -1943,12 +1943,12 @@ contains
             do d = 1, num_delayed_groups
               if (i_nuclide > 0) then
                 score = score + atom_density * flux * &
-                     get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                     get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
+                     get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                     get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
               else
                 score = score + flux * &
-                     get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
-                     get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
+                     get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
+                     get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
               end if
             end do
           end if
@@ -1973,20 +1973,20 @@ contains
           end if
           if (i_nuclide > 0) then
             score = score * atom_density * &
-                 get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                 get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_KAPPA_FISSION, p_g) / &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
           else
             score = score * &
-                 get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g) / &
-                 get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_KAPPA_FISSION, p_g) / &
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
           end if
         else
           if (i_nuclide > 0) then
-            score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * &
+            score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_KAPPA_FISSION, p_g) * &
                  atom_density * flux
           else
             score = flux * &
-                 get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g)
+                 get_macro_xs_c(p % material, tid, MG_GET_XS_KAPPA_FISSION, p_g)
 
           end if
         end if
@@ -2005,16 +2005,6 @@ contains
            t % results(RESULT_VALUE, score_index, filter_index) + score
 
     end do SCORE_LOOP
-
-    ! Reset temporary Mgxs indices
-    call reset_macro_angle_index_c(p % material, last_mat_pol, last_mat_azi, &
-                                   last_mat_uvw);
-
-    if (i_nuclide > 0) then
-      call reset_nuclide_temperature_index_c(i_nuclide, last_nuc_temp)
-      call reset_nuclide_angle_index_c(i_nuclide, last_nuc_pol, last_nuc_azi, &
-                                       last_nuc_uvw)
-    end if
   end subroutine score_general_mg
 
 !===============================================================================
diff --git a/src/tracking.F90 b/src/tracking.F90
index c832ca9cc..80fb55a80 100644
--- a/src/tracking.F90
+++ b/src/tracking.F90
@@ -2,6 +2,10 @@ module tracking
 
   use, intrinsic :: ISO_C_BINDING
 
+#ifdef _OPENMP
+  use omp_lib
+#endif
+
   use constants
   use error,              only: warning, write_message
   use geometry_header,    only: cells
@@ -48,6 +52,12 @@ contains
     real(8) :: d_collision            ! sampled distance to collision
     real(8) :: distance               ! distance particle travels
     logical :: found_cell             ! found cell which particle is in?
+    integer(C_INT) :: tid
+#ifdef _OPENMP
+    tid = OMP_GET_THREAD_NUM()
+#else
+    tid = 0
+#endif
 
     ! Display message if high verbosity or trace is on
     if (verbosity >= 9 .or. trace) then
@@ -114,7 +124,7 @@ contains
           end if
         else
           ! Get the MG data
-          call calculate_xs_c(p % material, p % g, p % sqrtkT, &
+          call calculate_xs_c(p % material, tid, p % g, p % sqrtkT, &
                p % coord(p % n_coord) % uvw, material_xs % total, &
                material_xs % absorption, material_xs % nu_fission)