diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index bba260c3a4b..8c5e24a1308 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -285,6 +285,7 @@ enum class MgxsType {
   PROMPT_NU_FISSION,
   DELAYED_NU_FISSION,
   NU_FISSION,
+  CHI,
   CHI_PROMPT,
   CHI_DELAYED
 };
diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h
index feafde68dd3..134f738fa08 100644
--- a/include/openmc/xsdata.h
+++ b/include/openmc/xsdata.h
@@ -83,6 +83,9 @@ class XsData {
   // delayed_nu_fission has the following dimensions:
   // [angle][delayed group][incoming group]
   xt::xtensor<double, 3> delayed_nu_fission;
+  // chi has the following dimensions:
+  // [angle][incoming group][outgoing group]
+  xt::xtensor<double, 3> chi;
   // chi_prompt has the following dimensions:
   // [angle][incoming group][outgoing group]
   xt::xtensor<double, 3> chi_prompt;
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index a2c479f2156..1a3fa141310 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -434,31 +434,39 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
 {
   XsData* xs_t = &xs[t];
   double val;
+  std::string mgxs_type;
   switch (xstype) {
   case MgxsType::TOTAL:
     val = xs_t->total(a, gin);
+    mgxs_type = "total";
     break;
   case MgxsType::NU_FISSION:
     val = fissionable ? xs_t->nu_fission(a, gin) : 0.;
+    mgxs_type = "nu-fission";
     break;
   case MgxsType::ABSORPTION:
     val = xs_t->absorption(a, gin);
+    mgxs_type = "absorption";
     ;
     break;
   case MgxsType::FISSION:
     val = fissionable ? xs_t->fission(a, gin) : 0.;
+    mgxs_type = "fission";
     break;
   case MgxsType::KAPPA_FISSION:
     val = fissionable ? xs_t->kappa_fission(a, gin) : 0.;
+    mgxs_type = "kappa-fission";
     break;
   case MgxsType::NU_SCATTER:
   case MgxsType::SCATTER:
   case MgxsType::NU_SCATTER_FMU:
   case MgxsType::SCATTER_FMU:
     val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
+    mgxs_type = "scatter_data";
     break;
   case MgxsType::PROMPT_NU_FISSION:
     val = fissionable ? xs_t->prompt_nu_fission(a, gin) : 0.;
+    mgxs_type = "prompt-nu-fission";
     break;
   case MgxsType::DELAYED_NU_FISSION:
     if (fissionable) {
@@ -473,6 +481,23 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
     } else {
       val = 0.;
     }
+    mgxs_type = "delayed-nu-fission";
+    break;
+  case MgxsType::CHI:
+    if (fissionable) {
+      if (gout != nullptr) {
+        val = xs_t->chi(a, gin, *gout);
+      } else {
+        // provide an outgoing group-wise sum
+        val = 0.;
+        for (int g = 0; g < xs_t->chi.shape()[2]; g++) {
+          val += xs_t->chi(a, gin, g);
+        }
+      }
+    } else {
+      val = 0.;
+    }
+    mgxs_type = "chi";
     break;
   case MgxsType::CHI_PROMPT:
     if (fissionable) {
@@ -488,6 +513,7 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
     } else {
       val = 0.;
     }
+    mgxs_type = "chi-prompt";
     break;
   case MgxsType::CHI_DELAYED:
     if (fissionable) {
@@ -515,9 +541,11 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
     } else {
       val = 0.;
     }
+    mgxs_type = "chi-delayed";
     break;
   case MgxsType::INVERSE_VELOCITY:
     val = xs_t->inverse_velocity(a, gin);
+    mgxs_type = "inverse-velocity";
     break;
   case MgxsType::DECAY_RATE:
     if (dg != nullptr) {
@@ -525,10 +553,20 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
     } else {
       val = xs_t->decay_rate(a, 0);
     }
+    mgxs_type = "decay-rate";
     break;
   default:
     val = 0.;
   }
+  // TODO: need to add more robust handling of zero-values cross sections that
+  // produce NANs on normalization
+  if (val != val) {
+    warning(
+      fmt::format("The {} cross section for this material has a nan present "
+                  "(possibly due to a division by zero). Setting to zero...",
+        mgxs_type));
+    val = 0;
+  }
   return val;
 }
 
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index 2f5007fc927..1ed09a0a31e 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@@ -1156,8 +1156,7 @@ void FlatSourceDomain::flatten_xs()
           m.get_xs(MgxsType::FISSION, g_out, NULL, NULL, NULL, t, a);
         sigma_f_.push_back(sigma_f);
 
-        double chi =
-          m.get_xs(MgxsType::CHI_PROMPT, g_out, &g_out, NULL, NULL, t, a);
+        double chi = m.get_xs(MgxsType::CHI, g_out, &g_out, NULL, NULL, t, a);
         if (!std::isfinite(chi)) {
           // MGXS interface may return NaN in some cases, such as when material
           // is fissionable but has very small sigma_f.
diff --git a/src/xsdata.cpp b/src/xsdata.cpp
index 1929f51e6fa..2158765eff1 100644
--- a/src/xsdata.cpp
+++ b/src/xsdata.cpp
@@ -56,8 +56,9 @@ XsData::XsData(bool fissionable, AngleDistributionType scatter_format,
     // allocate delayed_nu_fission; [temperature][angle][delay group][in group]
     delayed_nu_fission = xt::zeros<double>(shape);
 
-    // chi_prompt; [temperature][angle][in group][out group]
+    // chi and chi_prompt; [temperature][angle][in group][out group]
     shape = {n_ang, n_g_, n_g_};
+    chi = xt::zeros<double>(shape);
     chi_prompt = xt::zeros<double>(shape);
 
     // chi_delayed; [temperature][angle][delay group][in group][out group]
@@ -133,8 +134,13 @@ void XsData::fission_vector_beta_from_hdf5(
   // Normalize chi by summing over the outgoing groups for each incoming angle
   temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis());
 
+  // Store chi in chi
+  chi = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
+
   // Now every incoming group in prompt_chi and delayed_chi is the normalized
   // chi we just made
+  // TODO: This is incorrect. This makes chi_prompt and chi_delayed identical
+  // to chi.
   chi_prompt = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
   chi_delayed =
     xt::view(temp_chi, xt::all(), xt::newaxis(), xt::newaxis(), xt::all());
@@ -176,20 +182,31 @@ void XsData::fission_vector_beta_from_hdf5(
 
 void XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
 {
-  // Data is provided separately as prompt + delayed nu-fission and chi
+  // If chi is included in the dataset, we should store it!
+  if (object_exists(xsdata_grp, "chi")) {
+    xt::xtensor<double, 2> temp_chi({n_ang, n_g_}, 0.);
+    read_nd_vector(xsdata_grp, "chi", temp_chi, true);
+    // Normalize chi by summing over the outgoing groups for each incoming angle
+    temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis());
+    // Now every incoming group in self.chi is the normalized chi we just made
+    chi = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
+  }
 
+  // Data is provided separately as prompt + delayed nu-fission and chi
   // Get chi-prompt
   xt::xtensor<double, 2> temp_chi_p({n_ang, n_g_}, 0.);
   read_nd_vector(xsdata_grp, "chi-prompt", temp_chi_p, true);
 
-  // Normalize chi by summing over the outgoing groups for each incoming angle
+  // Normalize chi-prompt by summing over the outgoing groups for each incoming
+  // angle
   temp_chi_p /= xt::view(xt::sum(temp_chi_p, {1}), xt::all(), xt::newaxis());
 
   // Get chi-delayed
   xt::xtensor<double, 3> temp_chi_d({n_ang, n_dg_, n_g_}, 0.);
   read_nd_vector(xsdata_grp, "chi-delayed", temp_chi_d, true);
 
-  // Normalize chi by summing over the outgoing groups for each incoming angle
+  // Normalize chi-delayed by summing over the outgoing groups for each incoming
+  // angle
   temp_chi_d /=
     xt::view(xt::sum(temp_chi_d, {2}), xt::all(), xt::all(), xt::newaxis());
 
@@ -217,6 +234,7 @@ void XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang)
   temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis());
 
   // Now every incoming group in self.chi is the normalized chi we just made
+  chi = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
   chi_prompt = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
 
   // Get nu-fission directly
@@ -268,6 +286,10 @@ void XsData::fission_matrix_beta_from_hdf5(
       xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis(), xt::newaxis()) *
       xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all());
 
+    // Store chi
+    chi =
+      chi_prompt * (1. - temp_beta_sum) + xt::sum(temp_beta * chi_delayed, {1});
+
   } else if (beta_ndims == ndim_target + 1) {
     xt::xtensor<double, 3> temp_beta({n_ang, n_dg_, n_g_}, 0.);
     read_nd_vector(xsdata_grp, "beta", temp_beta, true);
@@ -293,14 +315,20 @@ void XsData::fission_matrix_beta_from_hdf5(
     chi_delayed =
       xt::view(temp_beta, xt::all(), xt::all(), xt::all(), xt::newaxis()) *
       xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all());
+
+    // Store chi
+    chi =
+      chi_prompt * (1. - temp_beta_sum) + xt::sum(temp_beta * chi_delayed, {1});
   }
 
-  // Normalize both chis
+  // Normalize chis
   chi_prompt /=
     xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::all(), xt::newaxis());
 
   chi_delayed /= xt::view(
     xt::sum(chi_delayed, {3}), xt::all(), xt::all(), xt::all(), xt::newaxis());
+
+  chi /= xt::view(xt::sum(chi, {2}), xt::all(), xt::all(), xt::newaxis());
 }
 
 void XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
@@ -326,8 +354,8 @@ void XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
   // delayed_nu_fission is the sum over outgoing groups
   delayed_nu_fission = xt::sum(temp_matrix_d, {3});
 
-  // chi_prompt is this matrix but normalized over outgoing groups, which we
-  // have already stored in prompt_nu_fission
+  // chi_delayee is this matrix but normalized over outgoing groups, which we
+  // have already stored in delayed_nu_fission
   chi_delayed = temp_matrix_d / xt::view(delayed_nu_fission, xt::all(),
                                   xt::all(), xt::all(), xt::newaxis());
 }
@@ -346,6 +374,8 @@ void XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang)
 
   // chi_prompt is this matrix but normalized over outgoing groups, which we
   // have already stored in prompt_nu_fission
+  chi = temp_matrix /
+        xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis());
   chi_prompt = temp_matrix /
                xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis());
 }
@@ -525,6 +555,7 @@ void XsData::combine(
       kappa_fission += scalar * that->kappa_fission;
       fission += scalar * that->fission;
       delayed_nu_fission += scalar * that->delayed_nu_fission;
+      chi += scalar * that->chi;
       chi_prompt += scalar *
                     xt::view(xt::sum(that->prompt_nu_fission, {1}), xt::all(),
                       xt::newaxis(), xt::newaxis()) *
@@ -537,8 +568,9 @@ void XsData::combine(
     decay_rate += scalar * that->decay_rate;
   }
 
-  // Ensure the chi_prompt and chi_delayed are normalized to 1 for each
+  // Ensure chi, chi_prompt and chi_delayed are normalized to 1 for each
   // azimuthal angle and delayed group (for chi_delayed)
+  chi /= xt::view(xt::sum(chi, {2}), xt::all(), xt::all(), xt::newaxis());
   chi_prompt /=
     xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::all(), xt::newaxis());
   chi_delayed /= xt::view(